Geant4  10.00.p03
G4VEnergyLossProcess.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4VEnergyLossProcess.cc 81864 2014-06-06 11:30:54Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4VEnergyLossProcess
34 //
35 // Author: Vladimir Ivanchenko
36 //
37 // Creation date: 03.01.2002
38 //
39 // Modifications:
40 //
41 // 13-11-02 Minor fix - use normalised direction (V.Ivanchenko)
42 // 04-12-02 Minor change in PostStepDoIt (V.Ivanchenko)
43 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
44 // 26-12-02 Secondary production moved to derived classes (V.Ivanchenko)
45 // 04-01-03 Fix problem of very small steps for ions (V.Ivanchenko)
46 // 20-01-03 Migrade to cut per region (V.Ivanchenko)
47 // 24-01-03 Temporarily close a control on usage of couples (V.Ivanchenko)
48 // 24-01-03 Make models region aware (V.Ivanchenko)
49 // 05-02-03 Fix compilation warnings (V.Ivanchenko)
50 // 06-02-03 Add control on tmax in PostStepDoIt (V.Ivanchenko)
51 // 13-02-03 SubCutoffProcessors defined for regions (V.Ivanchenko)
52 // 15-02-03 Lambda table can be scaled (V.Ivanchenko)
53 // 17-02-03 Fix problem of store/restore tables (V.Ivanchenko)
54 // 18-02-03 Add control on CutCouple usage (V.Ivanchenko)
55 // 26-02-03 Simplify control on GenericIons (V.Ivanchenko)
56 // 06-03-03 Control on GenericIons using SubType+ update verbose (V.Ivanchenko)
57 // 10-03-03 Add Ion registration (V.Ivanchenko)
58 // 22-03-03 Add Initialisation of cash (V.Ivanchenko)
59 // 26-03-03 Remove finalRange modification (V.Ivanchenko)
60 // 09-04-03 Fix problem of negative range limit for non integral (V.Ivanchenko)
61 // 26-04-03 Fix retrieve tables (V.Ivanchenko)
62 // 06-05-03 Set defalt finalRange = 1 mm (V.Ivanchenko)
63 // 12-05-03 Update range calculations + lowKinEnergy (V.Ivanchenko)
64 // 13-05-03 Add calculation of precise range (V.Ivanchenko)
65 // 23-05-03 Remove tracking cuts (V.Ivanchenko)
66 // 03-06-03 Fix initialisation problem for STD ionisation (V.Ivanchenko)
67 // 21-07-03 Add UpdateEmModel method (V.Ivanchenko)
68 // 03-11-03 Fix initialisation problem in RetrievePhysicsTable (V.Ivanchenko)
69 // 04-11-03 Add checks in RetrievePhysicsTable (V.Ivanchenko)
70 // 12-11-03 G4EnergyLossSTD -> G4EnergyLossProcess (V.Ivanchenko)
71 // 21-01-04 Migrade to G4ParticleChangeForLoss (V.Ivanchenko)
72 // 27-02-04 Fix problem of loss in low presure gases, cleanup precise range
73 // calculation, use functions ForLoss in AlongStepDoIt (V.Ivanchenko)
74 // 10-03-04 Fix a problem of Precise Range table (V.Ivanchenko)
75 // 19-03-04 Fix a problem energy below lowestKinEnergy (V.Ivanchenko)
76 // 31-03-04 Fix a problem of retrieve tables (V.Ivanchenko)
77 // 21-07-04 Check weather AtRest are active or not (V.Ivanchenko)
78 // 03-08-04 Add pointer of DEDX table to all processes (V.Ivanchenko)
79 // 06-08-04 Clear up names of member functions (V.Ivanchenko)
80 // 06-08-04 Clear up names of member functions (V.Ivanchenko)
81 // 27-08-04 Add NeedBuildTables method (V.Ivanchneko)
82 // 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivantchenko)
83 // 11-03-05 Shift verbose level by 1 (V.Ivantchenko)
84 // 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
85 // 11-04-05 Use MaxSecondaryEnergy from a model (V.Ivanchenko)
86 // 25-07-05 Add extra protection PostStep for non-integral mode (V.Ivanchenko)
87 // 12-08-05 Integral=false; SetStepFunction(0.2, 0.1*mm) (mma)
88 // 18-08-05 Return back both AlongStep and PostStep from 7.0 (V.Ivanchenko)
89 // 02-09-05 Default StepFunction 0.2 1 mm + integral (V.Ivanchenko)
90 // 04-09-05 default lambdaFactor 0.8 (V.Ivanchenko)
91 // 05-10-05 protection against 0 energy loss added (L.Urban)
92 // 17-10-05 protection above has been removed (L.Urban)
93 // 06-01-06 reset currentCouple when StepFunction is changed (V.Ivanchenko)
94 // 10-01-06 PreciseRange -> CSDARange (V.Ivantchenko)
95 // 18-01-06 Clean up subcutoff including recalculation of presafety (VI)
96 // 20-01-06 Introduce G4EmTableType and reducing number of methods (VI)
97 // 22-03-06 Add control on warning printout AlongStep (VI)
98 // 23-03-06 Use isIonisation flag (V.Ivanchenko)
99 // 07-06-06 Do not reflect AlongStep in subcutoff regime (V.Ivanchenko)
100 // 14-01-07 add SetEmModel(index) and SetFluctModel() (mma)
101 // 16-01-07 add IonisationTable and IonisationSubTable (V.Ivanchenko)
102 // 16-02-07 set linLossLimit=1.e-6 (V.Ivanchenko)
103 // 13-03-07 use SafetyHelper instead of navigator (V.Ivanchenko)
104 // 10-04-07 use unique SafetyHelper (V.Ivanchenko)
105 // 12-04-07 Add verbosity at destruction (V.Ivanchenko)
106 // 25-04-07 move initialisation of safety helper to BuildPhysicsTable (VI)
107 // 27-10-07 Virtual functions moved to source (V.Ivanchenko)
108 // 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko)
109 // 01-25-09 (Xin Dong) Phase II change for Geant4 multi-threading:
110 // New methods SlavePreparePhysicsTable, SlaveBuildPhysicsTable
111 // Worker threads share physics tables with the master thread for
112 // this kind of process. This member function is used by worker
113 // threads to achieve the partial effect of the master thread when
114 // it builds physcis tables.
115 // 15-10-10 Fixed 4-momentum balance if deexcitation is active (L.Pandola)
116 // 30-05-12 Call new ApplySecondaryBiasing so 2ries may be unique (D. Sawkey)
117 // 30-05-12 Fix bug in forced biasing: now called on first step (D. Sawkey)
118 // 04-06-13 Adoptation to MT mode, adding internal cache to GetRangeForLoss,
119 // more accurate initialisation for ions (V.Ivanchenko)
120 //
121 // Class Description:
122 //
123 // It is the unified energy loss process it calculates the continuous
124 // energy loss for charged particles using a set of Energy Loss
125 // models valid for different energy regions. There are a possibility
126 // to create and access to dE/dx and range tables, or to calculate
127 // that information on fly.
128 // -------------------------------------------------------------------
129 //
130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
132 
133 #include "G4VEnergyLossProcess.hh"
134 #include "G4PhysicalConstants.hh"
135 #include "G4SystemOfUnits.hh"
136 #include "G4ProcessManager.hh"
137 #include "G4LossTableManager.hh"
138 #include "G4LossTableBuilder.hh"
139 #include "G4Step.hh"
140 #include "G4ParticleDefinition.hh"
141 #include "G4ParticleTable.hh"
142 #include "G4VEmModel.hh"
143 #include "G4VEmFluctuationModel.hh"
144 #include "G4DataVector.hh"
145 #include "G4PhysicsLogVector.hh"
146 #include "G4VParticleChange.hh"
147 #include "G4Gamma.hh"
148 #include "G4Electron.hh"
149 #include "G4Positron.hh"
150 #include "G4ProcessManager.hh"
151 #include "G4UnitsTable.hh"
152 #include "G4ProductionCutsTable.hh"
153 #include "G4Region.hh"
154 #include "G4RegionStore.hh"
155 #include "G4PhysicsTableHelper.hh"
156 #include "G4SafetyHelper.hh"
158 #include "G4EmConfigurator.hh"
159 #include "G4VAtomDeexcitation.hh"
160 #include "G4EmBiasingManager.hh"
161 #include "G4Log.hh"
162 
163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
164 
166  G4ProcessType type):
167  G4VContinuousDiscreteProcess(name, type),
168  secondaryParticle(0),
169  nSCoffRegions(0),
170  idxSCoffRegions(0),
171  nProcesses(0),
172  theDEDXTable(0),
173  theDEDXSubTable(0),
174  theDEDXunRestrictedTable(0),
175  theIonisationTable(0),
176  theIonisationSubTable(0),
177  theRangeTableForLoss(0),
178  theCSDARangeTable(0),
179  theSecondaryRangeTable(0),
180  theInverseRangeTable(0),
181  theLambdaTable(0),
182  theSubLambdaTable(0),
183  theDensityFactor(0),
184  theDensityIdx(0),
185  baseParticle(0),
186  minSubRange(0.1),
187  lossFluctuationFlag(true),
188  rndmStepFlag(false),
189  tablesAreBuilt(false),
190  integral(true),
191  isIon(false),
192  isIonisation(true),
193  useSubCutoff(false),
194  useDeexcitation(false),
195  particle(0),
196  currentCouple(0),
197  nWarnings(0),
198  mfpKinEnergy(0.0)
199 {
200  SetVerboseLevel(1);
201 
202  // low energy limit
203  lowestKinEnergy = 1.*eV;
204  preStepKinEnergy = 0.0;
205  preStepRangeEnergy = 0.0;
207 
208  // Size of tables assuming spline
209  minKinEnergy = 0.1*keV;
210  maxKinEnergy = 10.0*TeV;
211  nBins = 77;
212  maxKinEnergyCSDA = 1.0*GeV;
213  nBinsCSDA = 35;
214 
215  // default linear loss limit for spline
216  linLossLimit = 0.01;
217 
218  // default dRoverRange and finalRange
219  SetStepFunction(0.2, 1.0*mm);
220 
221  // default lambda factor
222  lambdaFactor = 0.8;
223 
224  // cross section biasing
225  biasFactor = 1.0;
226 
227  // particle types
231  theGenericIon = 0;
232 
233  // run time objects
238  ->GetSafetyHelper();
240 
241  // initialise model
243  lManager->Register(this);
244  fluctModel = 0;
245  currentModel = 0;
246  atomDeexcitation = 0;
247 
248  biasManager = 0;
249  biasFlag = false;
250  weightFlag = false;
251  isMaster = true;
252  lastIdx = 0;
253 
257 
258  scTracks.reserve(5);
259  secParticles.reserve(5);
260 
261  theCuts = theSubCuts = 0;
262  currentMaterial = 0;
263  currentCoupleIndex = 0;
266 
267  secID = biasID = subsecID = -1;
268 }
269 
270 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
271 
273 {
274  /*
275  G4cout << "** G4VEnergyLossProcess::~G4VEnergyLossProcess() for "
276  << GetProcessName()
277  << " isMaster: " << isMaster << G4endl;
278  */
279  Clean();
280 
281  if ( !baseParticle && isMaster ) {
282  //G4cout << " isIonisation " << isIonisation << " "
283  // << theDEDXTable << G4endl;
284 
285  if(theDEDXTable) {
287  if(isIonisation) {
288  delete theDEDXTable;
289  theDEDXTable = 0;
290  if(theDEDXSubTable) {
292  { theIonisationSubTable = 0; }
293  delete theDEDXSubTable;
294  theDEDXSubTable = 0;
295  }
296  }
297  }
298  if(theIonisationTable) {
299  delete theIonisationTable;
300  theIonisationTable = 0;
301  }
303  delete theIonisationSubTable;
305  }
309  }
311  delete theCSDARangeTable;
312  theCSDARangeTable = 0;
313  }
315  delete theRangeTableForLoss;
317  }
319  delete theInverseRangeTable;
321  }
322  if(theLambdaTable) {
323  delete theLambdaTable;
324  theLambdaTable = 0;
325  }
326  if(theSubLambdaTable) {
327  delete theSubLambdaTable;
328  theSubLambdaTable = 0;
329  }
330  }
331 
332  delete modelManager;
333  delete biasManager;
334  lManager->DeRegister(this);
335  //G4cout << "** all removed" << G4endl;
336 }
337 
338 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
339 
341 {
342  /*
343  if(1 < verboseLevel) {
344  G4cout << "G4VEnergyLossProcess::Clear() for " << GetProcessName()
345  << G4endl;
346  }
347  */
348  delete [] idxSCoffRegions;
349 
350  tablesAreBuilt = false;
351 
352  scProcesses.clear();
353  nProcesses = 0;
354 
358 }
359 
360 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
361 
363  const G4Material*,
364  G4double cut)
365 {
366  return cut;
367 }
368 
369 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
370 
372  G4VEmFluctuationModel* fluc,
373  const G4Region* region)
374 {
375  modelManager->AddEmModel(order, p, fluc, region);
376  if(p) { p->SetParticleChange(pParticleChange, fluc); }
377 }
378 
379 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
380 
382  G4double emin, G4double emax)
383 {
384  modelManager->UpdateEmModel(nam, emin, emax);
385 }
386 
387 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
388 
390 {
391  G4int n = emModels.size();
392  if(index >= n) { for(G4int i=n; i<=index; ++i) {emModels.push_back(0);} }
393  emModels[index] = p;
394 }
395 
396 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
397 
399 {
400  G4VEmModel* p = 0;
401  if(index >= 0 && index < G4int(emModels.size())) { p = emModels[index]; }
402  return p;
403 }
404 
405 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
406 
408 {
409  return modelManager->GetModel(idx, ver);
410 }
411 
412 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
413 
415 {
416  return modelManager->NumberOfModels();
417 }
418 
419 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
420 
421 void
423 {
424  if(1 < verboseLevel) {
425  G4cout << "G4VEnergyLossProcess::PreparePhysicsTable for "
426  << GetProcessName() << " for " << part.GetParticleName()
427  << " " << this << G4endl;
428  }
429 
430  const G4VEnergyLossProcess* masterProcess =
431  static_cast<const G4VEnergyLossProcess*>(GetMasterProcess());
432  if(masterProcess && masterProcess != this) { isMaster = false; }
433 
434  currentCouple = 0;
435  preStepLambda = 0.0;
437  fRange = DBL_MAX;
438  preStepKinEnergy = 0.0;
439  preStepRangeEnergy = 0.0;
440  chargeSqRatio = 1.0;
441  massRatio = 1.0;
442  reduceFactor = 1.0;
443  fFactor = 1.0;
444  lastIdx = 0;
445 
446  // Are particle defined?
447  if( !particle ) { particle = &part; }
448 
449  if(part.GetParticleType() == "nucleus") {
450 
451  G4String pname = part.GetParticleName();
452  if(pname != "deuteron" && pname != "triton" &&
453  pname != "alpha+" && pname != "helium" &&
454  pname != "hydrogen") {
455 
456  if(!theGenericIon) {
457  theGenericIon =
459  }
460  isIon = true;
464  size_t n = v->size();
465  for(size_t j=0; j<n; ++j) {
466  if((*v)[j] == this) {
468  break;
469  }
470  }
471  }
472  }
473  }
474 
475  if( particle != &part ) {
476  if(!isIon) {
477  lManager->RegisterExtraParticle(&part, this);
478  }
479  if(1 < verboseLevel) {
480  G4cout << "### G4VEnergyLossProcess::PreparePhysicsTable()"
481  << " interrupted for "
482  << part.GetParticleName() << " isIon= " << isIon
483  << " particle " << particle << " GenericIon " << theGenericIon
484  << G4endl;
485  }
486  return;
487  }
488 
489  Clean();
490  lManager->PreparePhysicsTable(&part, this, isMaster);
492 
493  // Base particle and set of models can be defined here
495 
496  const G4ProductionCutsTable* theCoupleTable=
498  size_t n = theCoupleTable->GetTableSize();
499 
500  theDEDXAtMaxEnergy.resize(n, 0.0);
501  theRangeAtMaxEnergy.resize(n, 0.0);
502  theEnergyOfCrossSectionMax.resize(n, 0.0);
503  theCrossSectionMax.resize(n, DBL_MAX);
504 
505  // Tables preparation
506  if (isMaster && !baseParticle) {
507 
508  if(theDEDXTable && isIonisation) {
512  }
517  }
518  }
519 
522 
523  if(theDEDXSubTable) {
524  theDEDXSubTable =
526  }
527 
528  if (lManager->BuildCSDARange()) {
533  }
534 
536 
537  if(isIonisation) {
542  }
543 
544  if (nSCoffRegions) {
545  theDEDXSubTable =
549  }
550  }
551 
552  // forced biasing
553  if(biasManager) {
555  biasFlag = false;
556  }
557 
558  G4double initialCharge = particle->GetPDGCharge();
559  G4double initialMass = particle->GetPDGMass();
560 
561  if (baseParticle) {
562  massRatio = (baseParticle->GetPDGMass())/initialMass;
563  G4double q = initialCharge/baseParticle->GetPDGCharge();
564  chargeSqRatio = q*q;
565  if(chargeSqRatio > 0.0) { reduceFactor = 1.0/(chargeSqRatio*massRatio); }
566  }
567 
568  // defined ID of secondary particles
569  if(isMaster) {
570  G4String nam1 = GetProcessName();
571  G4String nam4 = nam1 + "_split";
572  G4String nam5 = nam1 + "_subcut";
576  }
577 
578  // initialisation of models
580  for(G4int i=0; i<nmod; ++i) {
581  G4VEmModel* mod = modelManager->GetModel(i);
583  if(mod->HighEnergyLimit() > maxKinEnergy) {
585  }
586  }
587 
590 
591  // Sub Cutoff
592  if (nSCoffRegions>0) {
594 
595  if(nSCoffRegions>0) { idxSCoffRegions = new G4bool[n]; }
596  for (size_t j=0; j<n; ++j) {
597 
598  const G4MaterialCutsCouple* couple =
599  theCoupleTable->GetMaterialCutsCouple(j);
600  const G4ProductionCuts* pcuts = couple->GetProductionCuts();
601 
602  if(nSCoffRegions>0) {
603  G4bool reg = false;
604  for(G4int i=0; i<nSCoffRegions; ++i) {
605  if( pcuts == scoffRegions[i]->GetProductionCuts()) { reg = true; }
606  }
607  idxSCoffRegions[j] = reg;
608  }
609  }
610  }
611 
612  if(1 < verboseLevel) {
613  G4cout << "G4VEnergyLossProcess::PrepearPhysicsTable() is done "
614  << " for local " << particle->GetParticleName()
615  << " isIon= " << isIon;
616  if(baseParticle) { G4cout << "; base: " << baseParticle->GetParticleName(); }
617  G4cout << " chargeSqRatio= " << chargeSqRatio
618  << " massRatio= " << massRatio
619  << " reduceFactor= " << reduceFactor << G4endl;
620  if (nSCoffRegions) {
621  G4cout << " SubCutoff Regime is ON for regions: " << G4endl;
622  for (G4int i=0; i<nSCoffRegions; ++i) {
623  const G4Region* r = scoffRegions[i];
624  G4cout << " " << r->GetName() << G4endl;
625  }
626  }
627  }
628 }
629 
630 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
631 
633 {
634  G4bool verb = false;
635  if(1 < verboseLevel || verb) {
636 
637  //if(1 < verboseLevel) {
638  G4cout << "### G4VEnergyLossProcess::BuildPhysicsTable() for "
639  << GetProcessName()
640  << " and particle " << part.GetParticleName()
641  << "; local: " << particle->GetParticleName();
642  if(baseParticle) {
643  G4cout << "; base: " << baseParticle->GetParticleName();
644  }
645  G4cout << " TablesAreBuilt= " << tablesAreBuilt
646  << " isIon= " << isIon << " " << this << G4endl;
647  }
648 
649  if(&part == particle) {
650 
652  if(isMaster) {
656 
657  } else {
658 
659  const G4VEnergyLossProcess* masterProcess =
660  static_cast<const G4VEnergyLossProcess*>(GetMasterProcess());
661 
662  // define density factors for worker thread
663  bld->InitialiseBaseMaterials(masterProcess->DEDXTable());
666 
667  // copy table pointers from master thread
668  SetDEDXTable(masterProcess->DEDXTable(),fRestricted);
670  SetDEDXTable(masterProcess->DEDXunRestrictedTable(),fTotal);
671  SetDEDXTable(masterProcess->IonisationTable(),fIsIonisation);
673  SetRangeTableForLoss(masterProcess->RangeTableForLoss());
674  SetCSDARangeTable(masterProcess->CSDARangeTable());
676  SetInverseRangeTable(masterProcess->InverseRangeTable());
677  SetLambdaTable(masterProcess->LambdaTable());
678  SetSubLambdaTable(masterProcess->SubLambdaTable());
679  isIonisation = masterProcess->IsIonisationProcess();
680 
681  tablesAreBuilt = true;
682  // local initialisation of models
683  G4bool printing = true;
684  G4int numberOfModels = modelManager->NumberOfModels();
685  for(G4int i=0; i<numberOfModels; ++i) {
686  G4VEmModel* mod = GetModelByIndex(i, printing);
687  G4VEmModel* mod0= masterProcess->GetModelByIndex(i,printing);
688  mod->InitialiseLocal(particle, mod0);
689  }
690 
692  }
693 
694  // needs to be done only once
696  }
697  // explicitly defined printout by particle name
698  G4String num = part.GetParticleName();
699  if(1 < verboseLevel ||
700  (0 < verboseLevel && (num == "e-" ||
701  num == "e+" || num == "mu+" ||
702  num == "mu-" || num == "proton"||
703  num == "pi+" || num == "pi-" ||
704  num == "kaon+" || num == "kaon-" ||
705  num == "alpha" || num == "anti_proton" ||
706  num == "GenericIon")))
707  {
708  PrintInfoDefinition(part);
709  }
710 
711  // Added tracking cut to avoid tracking artifacts
712  // identify deexcitation flag
713  if(isIonisation) {
716  if(atomDeexcitation) {
717  if(atomDeexcitation->IsPIXEActive()) { useDeexcitation = true; }
718  }
719  }
720 
721  //if(1 < verboseLevel || verb) {
722  if(1 < verboseLevel) {
723  G4cout << "### G4VEnergyLossProcess::BuildPhysicsTable() done for "
724  << GetProcessName()
725  << " and particle " << part.GetParticleName();
726  if(isIonisation) { G4cout << " isIonisation flag = 1"; }
727  G4cout << G4endl;
728  }
729 }
730 
731 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
732 
734 {
735  G4bool verb = false;
736  if(1 < verboseLevel || verb) {
737  G4cout << "G4VEnergyLossProcess::BuildDEDXTable() of type " << tType
738  << " for " << GetProcessName()
739  << " and particle " << particle->GetParticleName()
740  << G4endl;
741  }
742  G4PhysicsTable* table = 0;
743  G4double emax = maxKinEnergy;
744  G4int bin = nBins;
745 
746  if(fTotal == tType) {
747  emax = maxKinEnergyCSDA;
748  bin = nBinsCSDA;
749  table = theDEDXunRestrictedTable;
750  } else if(fRestricted == tType) {
751  table = theDEDXTable;
752  } else if(fSubRestricted == tType) {
753  table = theDEDXSubTable;
754  } else {
755  G4cout << "G4VEnergyLossProcess::BuildDEDXTable WARNING: wrong type "
756  << tType << G4endl;
757  }
758 
759  // Access to materials
760  const G4ProductionCutsTable* theCoupleTable=
762  size_t numOfCouples = theCoupleTable->GetTableSize();
763 
764  if(1 < verboseLevel || verb) {
765  G4cout << numOfCouples << " materials"
766  << " minKinEnergy= " << minKinEnergy
767  << " maxKinEnergy= " << emax
768  << " nbin= " << bin
769  << " EmTableType= " << tType
770  << " table= " << table << " " << this
771  << G4endl;
772  }
773  if(!table) { return table; }
774 
776  G4bool splineFlag = lManager->SplineFlag();
777  G4PhysicsLogVector* aVector = 0;
778  G4PhysicsLogVector* bVector = 0;
779 
780  for(size_t i=0; i<numOfCouples; ++i) {
781 
782  if(1 < verboseLevel || verb) {
783  G4cout << "G4VEnergyLossProcess::BuildDEDXVector Idx= " << i
784  << " flagTable= " << table->GetFlag(i)
785  << " Flag= " << bld->GetFlag(i) << G4endl;
786  }
787  if(bld->GetFlag(i)) {
788 
789  // create physics vector and fill it
790  const G4MaterialCutsCouple* couple =
791  theCoupleTable->GetMaterialCutsCouple(i);
792  delete (*table)[i];
793  if(!bVector) {
794  aVector = new G4PhysicsLogVector(minKinEnergy, emax, bin);
795  bVector = aVector;
796  } else {
797  aVector = new G4PhysicsLogVector(*bVector);
798  }
799  aVector->SetSpline(splineFlag);
800 
801  modelManager->FillDEDXVector(aVector, couple, tType);
802  if(splineFlag) { aVector->FillSecondDerivatives(); }
803 
804  // Insert vector for this material into the table
805  G4PhysicsTableHelper::SetPhysicsVector(table, i, aVector);
806  }
807  }
808 
809  if(1 < verboseLevel || verb) {
810  G4cout << "G4VEnergyLossProcess::BuildDEDXTable(): table is built for "
812  << " and process " << GetProcessName()
813  << G4endl;
814  //if(2 < verboseLevel) G4cout << (*table) << G4endl;
815  }
816 
817  return table;
818 }
819 
820 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
821 
823 {
824  G4PhysicsTable* table = 0;
825 
826  if(fRestricted == tType) {
827  table = theLambdaTable;
828  } else if(fSubRestricted == tType) {
829  table = theSubLambdaTable;
830  } else {
831  G4cout << "G4VEnergyLossProcess::BuildLambdaTable WARNING: wrong type "
832  << tType << G4endl;
833  }
834 
835  if(1 < verboseLevel) {
836  G4cout << "G4VEnergyLossProcess::BuildLambdaTable() of type "
837  << tType << " for process "
838  << GetProcessName() << " and particle "
840  << " EmTableType= " << tType
841  << " table= " << table
842  << G4endl;
843  }
844  if(!table) {return table;}
845 
846  // Access to materials
847  const G4ProductionCutsTable* theCoupleTable=
849  size_t numOfCouples = theCoupleTable->GetTableSize();
850 
854 
855  G4bool splineFlag = lManager->SplineFlag();
856  G4PhysicsLogVector* aVector = 0;
858 
859  for(size_t i=0; i<numOfCouples; ++i) {
860 
861  if (bld->GetFlag(i)) {
862 
863  // create physics vector and fill it
864  const G4MaterialCutsCouple* couple =
865  theCoupleTable->GetMaterialCutsCouple(i);
866  delete (*table)[i];
867 
868  G4bool startNull = true;
869  G4double emin =
870  MinPrimaryEnergy(particle,couple->GetMaterial(),(*theCuts)[i]);
871  if(minKinEnergy > emin) {
872  emin = minKinEnergy;
873  startNull = false;
874  }
875 
876  G4double emax = maxKinEnergy;
877  if(emax <= emin) { emax = 2*emin; }
878  G4int bin = G4lrint(nBins*G4Log(emax/emin)/scale);
879  if(bin < 3) { bin = 3; }
880  aVector = new G4PhysicsLogVector(emin, emax, bin);
881  aVector->SetSpline(splineFlag);
882 
883  modelManager->FillLambdaVector(aVector, couple, startNull, tType);
884  if(splineFlag) { aVector->FillSecondDerivatives(); }
885 
886  // Insert vector for this material into the table
887  G4PhysicsTableHelper::SetPhysicsVector(table, i, aVector);
888  }
889  }
890 
891  if(1 < verboseLevel) {
892  G4cout << "Lambda table is built for "
894  << G4endl;
895  }
896 
897  return table;
898 }
899 
900 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
901 
902 void
904 {
905  if(0 < verboseLevel) {
906  G4cout << std::setprecision(6);
907  G4cout << G4endl << GetProcessName() << ": for "
908  << part.GetParticleName()
909  << " SubType= " << GetProcessSubType()
910  << G4endl;
911  G4cout << " dE/dx and range tables from "
912  << G4BestUnit(minKinEnergy,"Energy")
913  << " to " << G4BestUnit(maxKinEnergy,"Energy")
914  << " in " << nBins << " bins" << G4endl
915  << " Lambda tables from threshold to "
916  << G4BestUnit(maxKinEnergy,"Energy")
917  << " in " << nBins << " bins, spline: "
918  << lManager->SplineFlag()
919  << G4endl;
921  G4cout << " finalRange(mm)= " << finalRange/mm
922  << ", dRoverRange= " << dRoverRange
923  << ", integral: " << integral
924  << ", fluct: " << lossFluctuationFlag
925  << ", linLossLimit= " << linLossLimit
926  << G4endl;
927  }
928  PrintInfo();
931  G4cout << " CSDA range table up"
932  << " to " << G4BestUnit(maxKinEnergyCSDA,"Energy")
933  << " in " << nBinsCSDA << " bins" << G4endl;
934  }
935  if(nSCoffRegions>0 && isIonisation) {
936  G4cout << " Subcutoff sampling in " << nSCoffRegions
937  << " regions" << G4endl;
938  }
939  if(2 < verboseLevel) {
940  G4cout << " DEDXTable address= " << theDEDXTable << G4endl;
941  if(theDEDXTable && isIonisation) G4cout << (*theDEDXTable) << G4endl;
942  G4cout << "non restricted DEDXTable address= "
945  G4cout << (*theDEDXunRestrictedTable) << G4endl;
946  }
948  G4cout << (*theDEDXSubTable) << G4endl;
949  }
950  G4cout << " CSDARangeTable address= " << theCSDARangeTable
951  << G4endl;
953  G4cout << (*theCSDARangeTable) << G4endl;
954  }
955  G4cout << " RangeTableForLoss address= " << theRangeTableForLoss
956  << G4endl;
958  G4cout << (*theRangeTableForLoss) << G4endl;
959  }
960  G4cout << " InverseRangeTable address= " << theInverseRangeTable
961  << G4endl;
963  G4cout << (*theInverseRangeTable) << G4endl;
964  }
965  G4cout << " LambdaTable address= " << theLambdaTable << G4endl;
967  G4cout << (*theLambdaTable) << G4endl;
968  }
969  G4cout << " SubLambdaTable address= " << theSubLambdaTable
970  << G4endl;
972  G4cout << (*theSubLambdaTable) << G4endl;
973  }
974  }
975  }
976 }
977 
978 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
979 
981 {
982  G4RegionStore* regionStore = G4RegionStore::GetInstance();
983  const G4Region* reg = r;
984  if (!reg) {
985  reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);
986  }
987 
988  // the region is in the list
989  if (nSCoffRegions) {
990  for (G4int i=0; i<nSCoffRegions; ++i) {
991  if (reg == scoffRegions[i]) {
992  return;
993  }
994  }
995  }
996 
997  // new region
998  if(val) {
999  useSubCutoff = true;
1000  scoffRegions.push_back(reg);
1001  ++nSCoffRegions;
1002  } else {
1003  useSubCutoff = false;
1004  }
1005 }
1006 
1007 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1008 
1010 {
1011  /*
1012  G4cout << track->GetDefinition()->GetParticleName()
1013  << " e(MeV)= " << track->GetKineticEnergy()
1014  << " baseParticle " << baseParticle << " proc " << this;
1015  if(particle) G4cout << " " << particle->GetParticleName();
1016  G4cout << " isIon= " << isIon << " dedx " << theDEDXTable <<G4endl;
1017  */
1018  // reset parameters for the new track
1021  preStepRangeEnergy = 0.0;
1022 
1023  // reset ion
1024  if(isIon) {
1025  chargeSqRatio = 0.5;
1026 
1027  G4double newmass = track->GetDefinition()->GetPDGMass();
1028  if(baseParticle) {
1029  massRatio = baseParticle->GetPDGMass()/newmass;
1030  } else if(theGenericIon) {
1031  massRatio = proton_mass_c2/newmass;
1032  } else {
1033  massRatio = 1.0;
1034  }
1035  }
1036  // forced biasing only for primary particles
1037  if(biasManager) {
1038  if(0 == track->GetParentID()) {
1039  // primary particle
1040  biasFlag = true;
1042  }
1043  }
1044 }
1045 
1046 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1047 
1050  G4GPILSelection* selection)
1051 {
1052  G4double x = DBL_MAX;
1053  *selection = aGPILSelection;
1054  if(isIonisation) {
1056  x = fRange;
1057  G4double finR = finalRange;
1058  if(rndmStepFlag) {
1059  finR = std::min(finR,
1061  }
1062  if(fRange > finR) {
1063  x = fRange*dRoverRange + finR*(1.0 - dRoverRange)*(2.0 - finR/fRange);
1064  }
1065  /*
1066  if(particle->GetPDGMass() > 0.9*GeV)
1067  G4cout<<GetProcessName()<<": e= "<<preStepKinEnergy
1068  <<" range= "<<fRange << " idx= " << basedCoupleIndex
1069  << " finR= " << finR
1070  << " limit= " << x <<G4endl;
1071  G4cout << "massRatio= " << massRatio << " Q^2= " << chargeSqRatio
1072  << " finR= " << finR << " dRoverRange= " << dRoverRange
1073  << " finalRange= " << finalRange << G4endl;
1074  */
1075  }
1076  //G4cout<<GetProcessName()<<": e= "<<preStepKinEnergy
1077  //<<" stepLimit= "<<x<<G4endl;
1078  return x;
1079 }
1080 
1081 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1082 
1084  const G4Track& track,
1085  G4double previousStepSize,
1087 {
1088  // condition is set to "Not Forced"
1089  *condition = NotForced;
1090  G4double x = DBL_MAX;
1091 
1092  // initialisation of material, mass, charge, model
1093  // at the beginning of the step
1096  preStepScaledEnergy = preStepKinEnergy*massRatio;
1098 
1101  return x;
1102  }
1103 
1104  // change effective charge of an ion on fly
1105  if(isIon) {
1107  if(q2 != chargeSqRatio && q2 > 0.0) {
1108  chargeSqRatio = q2;
1109  fFactor = q2*biasFactor*(*theDensityFactor)[currentCoupleIndex];
1110  reduceFactor = 1.0/(fFactor*massRatio);
1111  }
1112  }
1113  // if(particle->GetPDGMass() > 0.9*GeV)
1114  //G4cout << "q2= "<<chargeSqRatio << " massRatio= " << massRatio << G4endl;
1115  // initialisation for sampling of the interaction length
1116  //if(previousStepSize <= 0.0) { theNumberOfInteractionLengthLeft = -1.0; }
1117  //if(theNumberOfInteractionLengthLeft < 0.0) { mfpKinEnergy = DBL_MAX; }
1118 
1119  // forced biasing only for primary particles
1120  if(biasManager) {
1121  if(0 == track.GetParentID()) {
1122  if(biasFlag &&
1124  return biasManager->GetStepLimit(currentCoupleIndex, previousStepSize);
1125  }
1126  }
1127  }
1128 
1129  // compute mean free path
1133 
1134  // zero cross section
1135  if(preStepLambda <= 0.0) {
1138  }
1139  }
1140 
1141  // non-zero cross section
1142  if(preStepLambda > 0.0) {
1144 
1145  // beggining of tracking (or just after DoIt of this process)
1146  // ResetNumberOfInteractionLengthLeft();
1149 
1150  } else if(currentInteractionLength < DBL_MAX) {
1151 
1152  // subtract NumberOfInteractionLengthLeft using previous step
1154  previousStepSize/currentInteractionLength;
1155  // SubtractNumberOfInteractionLengthLeft(previousStepSize);
1158  //theNumberOfInteractionLengthLeft = perMillion;
1159  }
1160  }
1161 
1162  // new mean free path and step limit
1165  }
1166 #ifdef G4VERBOSE
1167  if (verboseLevel>2){
1168  // if(particle->GetPDGMass() > 0.9*GeV){
1169  G4cout << "G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength ";
1170  G4cout << "[ " << GetProcessName() << "]" << G4endl;
1171  G4cout << " for " << track.GetDefinition()->GetParticleName()
1172  << " in Material " << currentMaterial->GetName()
1173  << " Ekin(MeV)= " << preStepKinEnergy/MeV
1174  << " " << track.GetMaterial()->GetName()
1175  <<G4endl;
1176  G4cout << "MeanFreePath = " << currentInteractionLength/cm << "[cm]"
1177  << "InteractionLength= " << x/cm <<"[cm] " <<G4endl;
1178  }
1179 #endif
1180  return x;
1181 }
1182 
1183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1184 
1186  const G4Step& step)
1187 {
1189  // The process has range table - calculate energy loss
1191  return &fParticleChange;
1192  }
1193 
1194  // Get the actual (true) Step length
1195  G4double length = step.GetStepLength();
1196  if(length <= 0.0) { return &fParticleChange; }
1197  G4double eloss = 0.0;
1198 
1199  /*
1200  if(-1 < verboseLevel) {
1201  const G4ParticleDefinition* d = track.GetParticleDefinition();
1202  G4cout << "AlongStepDoIt for "
1203  << GetProcessName() << " and particle "
1204  << d->GetParticleName()
1205  << " eScaled(MeV)= " << preStepScaledEnergy/MeV
1206  << " range(mm)= " << fRange/mm
1207  << " s(mm)= " << length/mm
1208  << " rf= " << reduceFactor
1209  << " q^2= " << chargeSqRatio
1210  << " md= " << d->GetPDGMass()
1211  << " status= " << track.GetTrackStatus()
1212  << " " << track.GetMaterial()->GetName()
1213  << G4endl;
1214  }
1215  */
1216 
1217  const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
1218 
1219  // define new weight for primary and secondaries
1221  if(weightFlag) {
1222  weight /= biasFactor;
1224  }
1225 
1226  // stopping
1227  if (length >= fRange) {
1228  eloss = preStepKinEnergy;
1229  if (useDeexcitation) {
1231  eloss, currentCoupleIndex);
1232  if(scTracks.size() > 0) { FillSecondariesAlongStep(eloss, weight); }
1233  if(eloss < 0.0) { eloss = 0.0; }
1234  }
1237  return &fParticleChange;
1238  }
1239  //G4cout << theDEDXTable << " idx= " << basedCoupleIndex
1240  // << " " << GetProcessName() << " "<< currentMaterial->GetName()<<G4endl;
1241  //if(particle->GetParticleName() == "e-")G4cout << (*theDEDXTable) <<G4endl;
1242  // Short step
1244 
1245  //G4cout << "eloss= " << eloss << G4endl;
1246 
1247  // Long step
1248  if(eloss > preStepKinEnergy*linLossLimit) {
1249 
1250  G4double x = (fRange - length)/reduceFactor;
1251  //G4cout << "x= " << x << " " << theInverseRangeTable << G4endl;
1253 
1254  /*
1255  if(-1 < verboseLevel)
1256  G4cout << "Long STEP: rPre(mm)= "
1257  << GetScaledRangeForScaledEnergy(preStepScaledEnergy)/mm
1258  << " rPost(mm)= " << x/mm
1259  << " ePre(MeV)= " << preStepScaledEnergy/MeV
1260  << " eloss(MeV)= " << eloss/MeV
1261  << " eloss0(MeV)= "
1262  << GetDEDXForScaledEnergy(preStepScaledEnergy)*length/MeV
1263  << " lim(MeV)= " << preStepKinEnergy*linLossLimit/MeV
1264  << G4endl;
1265  */
1266  }
1267 
1268  /*
1269  G4double eloss0 = eloss;
1270  if(-1 < verboseLevel ) {
1271  G4cout << "Before fluct: eloss(MeV)= " << eloss/MeV
1272  << " e-eloss= " << preStepKinEnergy-eloss
1273  << " step(mm)= " << length/mm
1274  << " range(mm)= " << fRange/mm
1275  << " fluct= " << lossFluctuationFlag
1276  << G4endl;
1277  }
1278  */
1279 
1280  G4double cut = (*theCuts)[currentCoupleIndex];
1281  G4double esec = 0.0;
1282 
1283  // SubCutOff
1284  if(useSubCutoff) {
1286 
1287  G4bool yes = false;
1288  G4StepPoint* prePoint = step.GetPreStepPoint();
1289 
1290  // Check boundary
1291  if(prePoint->GetStepStatus() == fGeomBoundary) { yes = true; }
1292 
1293  // Check PrePoint
1294  else {
1295  G4double preSafety = prePoint->GetSafety();
1296  G4double rcut =
1298 
1299  // recompute presafety
1300  if(preSafety < rcut) {
1301  preSafety = safetyHelper->ComputeSafety(prePoint->GetPosition());
1302  }
1303 
1304  if(preSafety < rcut) { yes = true; }
1305 
1306  // Check PostPoint
1307  else {
1308  G4double postSafety = preSafety - length;
1309  if(postSafety < rcut) {
1310  postSafety = safetyHelper->ComputeSafety(
1311  step.GetPostStepPoint()->GetPosition());
1312  if(postSafety < rcut) { yes = true; }
1313  }
1314  }
1315  }
1316 
1317  // Decided to start subcut sampling
1318  if(yes) {
1319 
1320  cut = (*theSubCuts)[currentCoupleIndex];
1322  esec = SampleSubCutSecondaries(scTracks, step,
1323  currentModel,currentCoupleIndex);
1324  // add bremsstrahlung sampling
1325  /*
1326  if(nProcesses > 0) {
1327  for(G4int i=0; i<nProcesses; ++i) {
1328  (scProcesses[i])->SampleSubCutSecondaries(
1329  scTracks, step, (scProcesses[i])->
1330  SelectModelForMaterial(preStepKinEnergy, currentCoupleIndex),
1331  currentCoupleIndex);
1332  }
1333  }
1334  */
1335  }
1336  }
1337  }
1338 
1339  // Corrections, which cannot be tabulated
1340  if(isIon) {
1341  G4double eadd = 0.0;
1342  G4double eloss_before = eloss;
1344  eloss, eadd, length);
1345  if(eloss < 0.0) { eloss = 0.5*eloss_before; }
1346  }
1347 
1348  // Sample fluctuations
1349  if (lossFluctuationFlag) {
1351  if(fluc &&
1352  (eloss + esec + lowestKinEnergy) < preStepKinEnergy) {
1353 
1354  G4double tmax =
1355  std::min(currentModel->MaxSecondaryKinEnergy(dynParticle),cut);
1356  eloss = fluc->SampleFluctuations(currentCouple,dynParticle,
1357  tmax,length,eloss);
1358  /*
1359  if(-1 < verboseLevel)
1360  G4cout << "After fluct: eloss(MeV)= " << eloss/MeV
1361  << " fluc= " << (eloss-eloss0)/MeV
1362  << " ChargeSqRatio= " << chargeSqRatio
1363  << " massRatio= " << massRatio
1364  << " tmax= " << tmax
1365  << G4endl;
1366  */
1367  }
1368  }
1369 
1370  // deexcitation
1371  if (useDeexcitation) {
1372  G4double esecfluo = preStepKinEnergy - esec;
1373  G4double de = esecfluo;
1374  //G4double eloss0 = eloss;
1375  /*
1376  G4cout << "### 1: E(keV)= " << preStepKinEnergy/keV
1377  << " Efluomax(keV)= " << de/keV
1378  << " Eloss(keV)= " << eloss/keV << G4endl;
1379  */
1381  de, currentCoupleIndex);
1382 
1383  // sum of de-excitation energies
1384  esecfluo -= de;
1385 
1386  // subtracted from energy loss
1387  if(eloss >= esecfluo) {
1388  esec += esecfluo;
1389  eloss -= esecfluo;
1390  } else {
1391  esec += esecfluo;
1392  eloss = 0.0;
1393  }
1394  /*
1395  if(esecfluo > 0.0) {
1396  G4cout << "### 2: E(keV)= " << preStepKinEnergy/keV
1397  << " Esec(keV)= " << esec/keV
1398  << " Esecf(kV)= " << esecfluo/keV
1399  << " Eloss0(kV)= " << eloss0/keV
1400  << " Eloss(keV)= " << eloss/keV
1401  << G4endl;
1402  }
1403  */
1404  }
1405  if(scTracks.size() > 0) { FillSecondariesAlongStep(eloss, weight); }
1406 
1407  // Energy balanse
1408  G4double finalT = preStepKinEnergy - eloss - esec;
1409  if (finalT <= lowestKinEnergy) {
1410  eloss += finalT;
1411  finalT = 0.0;
1412  } else if(isIon) {
1415  currentMaterial,finalT));
1416  }
1417 
1418  if(eloss < 0.0) { eloss = 0.0; }
1421  /*
1422  if(-1 < verboseLevel) {
1423  G4double del = finalT + eloss + esec - preStepKinEnergy;
1424  G4cout << "Final value eloss(MeV)= " << eloss/MeV
1425  << " preStepKinEnergy= " << preStepKinEnergy
1426  << " postStepKinEnergy= " << finalT
1427  << " de(keV)= " << del/keV
1428  << " lossFlag= " << lossFluctuationFlag
1429  << " status= " << track.GetTrackStatus()
1430  << G4endl;
1431  }
1432  */
1433  return &fParticleChange;
1434 }
1435 
1436 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1437 
1438 void
1440 {
1441  G4int n0 = scTracks.size();
1442 
1443  // weight may be changed by biasing manager
1444  if(biasManager) {
1446  weight *=
1448  }
1449  }
1450 
1451  // fill secondaries
1452  G4int n = scTracks.size();
1454 
1455  for(G4int i=0; i<n; ++i) {
1456  G4Track* t = scTracks[i];
1457  if(t) {
1458  t->SetWeight(weight);
1460  if(i < n0) { t->SetCreatorModelIndex(secID); }
1461  else { t->SetCreatorModelIndex(biasID); }
1462  //G4cout << "Secondary(along step) has weight " << t->GetWeight()
1463  //<< ", kenergy " << t->GetKineticEnergy()/MeV << " MeV" <<G4endl;
1464  }
1465  }
1466  scTracks.clear();
1467 }
1468 
1469 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1470 
1471 G4double
1472 G4VEnergyLossProcess::SampleSubCutSecondaries(std::vector<G4Track*>& tracks,
1473  const G4Step& step,
1474  G4VEmModel* model,
1475  G4int idx)
1476 {
1477  // Fast check weather subcutoff can work
1478  G4double esec = 0.0;
1479  G4double subcut = (*theSubCuts)[idx];
1480  G4double cut = (*theCuts)[idx];
1481  if(cut <= subcut) { return esec; }
1482 
1483  const G4Track* track = step.GetTrack();
1484  const G4DynamicParticle* dp = track->GetDynamicParticle();
1486  G4double cross = (*theDensityFactor)[idx]*chargeSqRatio
1487  *(((*theSubLambdaTable)[(*theDensityIdx)[idx]])->Value(e, idxSubLambda));
1488  G4double length = step.GetStepLength();
1489 
1490  // negligible probability to get any interaction
1491  if(length*cross < perMillion) { return esec; }
1492  /*
1493  if(-1 < verboseLevel)
1494  G4cout << "<<< Subcutoff for " << GetProcessName()
1495  << " cross(1/mm)= " << cross*mm << ">>>"
1496  << " e(MeV)= " << preStepScaledEnergy
1497  << " matIdx= " << currentCoupleIndex
1498  << G4endl;
1499  */
1500 
1501  // Sample subcutoff secondaries
1502  G4StepPoint* preStepPoint = step.GetPreStepPoint();
1503  G4StepPoint* postStepPoint = step.GetPostStepPoint();
1504  G4ThreeVector prepoint = preStepPoint->GetPosition();
1505  G4ThreeVector dr = postStepPoint->GetPosition() - prepoint;
1506  G4double pretime = preStepPoint->GetGlobalTime();
1507  G4double dt = postStepPoint->GetGlobalTime() - pretime;
1508  //G4double dt = length/preStepPoint->GetVelocity();
1509  G4double fragment = 0.0;
1510 
1511  do {
1512  G4double del = -G4Log(G4UniformRand())/cross;
1513  fragment += del/length;
1514  if (fragment > 1.0) break;
1515 
1516  // sample secondaries
1517  secParticles.clear();
1519  dp,subcut,cut);
1520 
1521  // position of subcutoff particles
1522  G4ThreeVector r = prepoint + fragment*dr;
1523  std::vector<G4DynamicParticle*>::iterator it;
1524  for(it=secParticles.begin(); it!=secParticles.end(); ++it) {
1525 
1526  G4Track* t = new G4Track((*it), pretime + fragment*dt, r);
1529  tracks.push_back(t);
1530  esec += t->GetKineticEnergy();
1531  if (t->GetParticleDefinition() == thePositron) {
1532  esec += 2.0*electron_mass_c2;
1533  }
1534 
1535  /*
1536  if(-1 < verboseLevel)
1537  G4cout << "New track "
1538  << t->GetParticleDefinition()->GetParticleName()
1539  << " e(keV)= " << t->GetKineticEnergy()/keV
1540  << " fragment= " << fragment
1541  << G4endl;
1542  */
1543  }
1544  } while (fragment <= 1.0);
1545  return esec;
1546 }
1547 
1548 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1549 
1551  const G4Step& step)
1552 {
1553  // In all cases clear number of interaction lengths
1556 
1558  G4double finalT = track.GetKineticEnergy();
1559  if(finalT <= lowestKinEnergy) { return &fParticleChange; }
1560 
1561  G4double postStepScaledEnergy = finalT*massRatio;
1562 
1563  if(!currentModel->IsActive(postStepScaledEnergy)) {
1564  return &fParticleChange;
1565  }
1566  /*
1567  if(-1 < verboseLevel) {
1568  G4cout << GetProcessName()
1569  << "::PostStepDoIt: E(MeV)= " << finalT/MeV
1570  << G4endl;
1571  }
1572  */
1573 
1574  // forced process - should happen only once per track
1575  if(biasFlag) {
1577  biasFlag = false;
1578  }
1579  }
1580 
1581  // Integral approach
1582  if (integral) {
1583  G4double lx = GetLambdaForScaledEnergy(postStepScaledEnergy);
1584  /*
1585  if(preStepLambda<lx && 1 < verboseLevel && nWarnings<200) {
1586  G4cout << "WARNING: for " << particle->GetParticleName()
1587  << " and " << GetProcessName()
1588  << " E(MeV)= " << finalT/MeV
1589  << " preLambda= " << preStepLambda
1590  << " < " << lx << " (postLambda) "
1591  << G4endl;
1592  ++nWarnings;
1593  }
1594  */
1595  if(lx <= 0.0 || preStepLambda*G4UniformRand() > lx) {
1596  return &fParticleChange;
1597  }
1598  }
1599 
1600  SelectModel(postStepScaledEnergy);
1601 
1602  // define new weight for primary and secondaries
1604  if(weightFlag) {
1605  weight /= biasFactor;
1607  }
1608 
1609  const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
1610  G4double tcut = (*theCuts)[currentCoupleIndex];
1611 
1612  // sample secondaries
1613  secParticles.clear();
1614  //G4cout<< "Eprimary: "<<dynParticle->GetKineticEnergy()/MeV<<G4endl;
1616  dynParticle, tcut);
1617 
1618  G4int num0 = secParticles.size();
1619 
1620  // bremsstrahlung splitting or Russian roulette
1621  if(biasManager) {
1623  G4double eloss = 0.0;
1625  secParticles,
1626  track, currentModel,
1627  &fParticleChange, eloss,
1628  currentCoupleIndex, tcut,
1629  step.GetPostStepPoint()->GetSafety());
1630  if(eloss > 0.0) {
1633  }
1634  }
1635  }
1636 
1637  // save secondaries
1638  G4int num = secParticles.size();
1639  if(num > 0) {
1640 
1642  G4double time = track.GetGlobalTime();
1643 
1644  for (G4int i=0; i<num; ++i) {
1645  if(secParticles[i]) {
1646  G4Track* t = new G4Track(secParticles[i], time, track.GetPosition());
1648  t->SetWeight(weight);
1649  if(i < num0) { t->SetCreatorModelIndex(secID); }
1650  else { t->SetCreatorModelIndex(biasID); }
1651 
1652  //G4cout << "Secondary(post step) has weight " << t->GetWeight()
1653  // << ", kenergy " << t->GetKineticEnergy()/MeV << " MeV"
1654  // << " time= " << time/ns << " ns " << G4endl;
1656  }
1657  }
1658  }
1659 
1665  }
1666 
1667  /*
1668  if(-1 < verboseLevel) {
1669  G4cout << "::PostStepDoIt: Sample secondary; Efin= "
1670  << fParticleChange.GetProposedKineticEnergy()/MeV
1671  << " MeV; model= (" << currentModel->LowEnergyLimit()
1672  << ", " << currentModel->HighEnergyLimit() << ")"
1673  << " preStepLambda= " << preStepLambda
1674  << " dir= " << track.GetMomentumDirection()
1675  << " status= " << track.GetTrackStatus()
1676  << G4endl;
1677  }
1678  */
1679  return &fParticleChange;
1680 }
1681 
1682 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1683 
1685  const G4ParticleDefinition* part, const G4String& directory,
1686  G4bool ascii)
1687 {
1688  G4bool res = true;
1689  if ( baseParticle || part != particle ) return res;
1690 
1691  if(!StoreTable(part,theDEDXTable,ascii,directory,"DEDX"))
1692  {res = false;}
1693 
1694  if(!StoreTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr"))
1695  {res = false;}
1696 
1697  if(!StoreTable(part,theDEDXSubTable,ascii,directory,"SubDEDX"))
1698  {res = false;}
1699 
1700  if(!StoreTable(part,theIonisationTable,ascii,directory,"Ionisation"))
1701  {res = false;}
1702 
1703  if(!StoreTable(part,theIonisationSubTable,ascii,directory,"SubIonisation"))
1704  {res = false;}
1705 
1706  if(isIonisation &&
1707  !StoreTable(part,theCSDARangeTable,ascii,directory,"CSDARange"))
1708  {res = false;}
1709 
1710  if(isIonisation &&
1711  !StoreTable(part,theRangeTableForLoss,ascii,directory,"Range"))
1712  {res = false;}
1713 
1714  if(isIonisation &&
1715  !StoreTable(part,theInverseRangeTable,ascii,directory,"InverseRange"))
1716  {res = false;}
1717 
1718  if(!StoreTable(part,theLambdaTable,ascii,directory,"Lambda"))
1719  {res = false;}
1720 
1721  if(!StoreTable(part,theSubLambdaTable,ascii,directory,"SubLambda"))
1722  {res = false;}
1723 
1724  if ( res ) {
1725  if(0 < verboseLevel) {
1726  G4cout << "Physics tables are stored for "
1728  << " and process " << GetProcessName()
1729  << " in the directory <" << directory
1730  << "> " << G4endl;
1731  }
1732  } else {
1733  G4cout << "Fail to store Physics Tables for "
1735  << " and process " << GetProcessName()
1736  << " in the directory <" << directory
1737  << "> " << G4endl;
1738  }
1739  return res;
1740 }
1741 
1742 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
1743 
1744 G4bool
1746  const G4String& directory,
1747  G4bool ascii)
1748 {
1749  G4bool res = true;
1750  const G4String particleName = part->GetParticleName();
1751 
1752  if(1 < verboseLevel) {
1753  G4cout << "G4VEnergyLossProcess::RetrievePhysicsTable() for "
1754  << particleName << " and process " << GetProcessName()
1755  << "; tables_are_built= " << tablesAreBuilt
1756  << G4endl;
1757  }
1758  if(particle == part) {
1759 
1760  if ( !baseParticle ) {
1761 
1762  G4bool fpi = true;
1763  if(!RetrieveTable(part,theDEDXTable,ascii,directory,"DEDX",fpi))
1764  {fpi = false;}
1765 
1766  // ionisation table keeps individual dEdx and not sum of sub-processes
1767  if(!RetrieveTable(part,theDEDXTable,ascii,directory,"Ionisation",false))
1768  {fpi = false;}
1769 
1770  if(!RetrieveTable(part,theRangeTableForLoss,ascii,directory,"Range",fpi))
1771  {res = false;}
1772 
1773  if(!RetrieveTable(part,theDEDXunRestrictedTable,ascii,directory,
1774  "DEDXnr",false))
1775  {res = false;}
1776 
1777  if(!RetrieveTable(part,theCSDARangeTable,ascii,directory,
1778  "CSDARange",false))
1779  {res = false;}
1780 
1781  if(!RetrieveTable(part,theInverseRangeTable,ascii,directory,
1782  "InverseRange",fpi))
1783  {res = false;}
1784 
1785  if(!RetrieveTable(part,theLambdaTable,ascii,directory,"Lambda",true))
1786  {res = false;}
1787 
1788  G4bool yes = false;
1789  if(nSCoffRegions > 0) {yes = true;}
1790 
1791  if(!RetrieveTable(part,theDEDXSubTable,ascii,directory,"SubDEDX",yes))
1792  {res = false;}
1793 
1794  if(!RetrieveTable(part,theSubLambdaTable,ascii,directory,
1795  "SubLambda",yes))
1796  {res = false;}
1797 
1798  if(!fpi) yes = false;
1799  if(!RetrieveTable(part,theIonisationSubTable,ascii,directory,
1800  "SubIonisation",yes))
1801  {res = false;}
1802  }
1803  }
1804 
1805  return res;
1806 }
1807 
1808 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
1809 
1811  G4PhysicsTable* aTable, G4bool ascii,
1812  const G4String& directory,
1813  const G4String& tname)
1814 {
1815  G4bool res = true;
1816  if ( aTable ) {
1817  const G4String name = GetPhysicsTableFileName(part,directory,tname,ascii);
1818  if( !aTable->StorePhysicsTable(name,ascii)) res = false;
1819  }
1820  return res;
1821 }
1822 
1823 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
1824 
1825 G4bool
1827  G4PhysicsTable* aTable,
1828  G4bool ascii,
1829  const G4String& directory,
1830  const G4String& tname,
1831  G4bool mandatory)
1832 {
1833  G4bool isRetrieved = false;
1834  G4String filename = GetPhysicsTableFileName(part,directory,tname,ascii);
1835  if(aTable) {
1836  if(aTable->ExistPhysicsTable(filename)) {
1837  if(G4PhysicsTableHelper::RetrievePhysicsTable(aTable,filename,ascii)) {
1838  isRetrieved = true;
1839  if(lManager->SplineFlag()) {
1840  size_t n = aTable->length();
1841  for(size_t i=0; i<n; ++i) {
1842  if((*aTable)[i]) { (*aTable)[i]->SetSpline(true); }
1843  }
1844  }
1845  if (0 < verboseLevel) {
1846  G4cout << tname << " table for " << part->GetParticleName()
1847  << " is Retrieved from <" << filename << ">"
1848  << G4endl;
1849  }
1850  }
1851  }
1852  }
1853  if(mandatory && !isRetrieved) {
1854  if(0 < verboseLevel) {
1855  G4cout << tname << " table for " << part->GetParticleName()
1856  << " from file <"
1857  << filename << "> is not Retrieved"
1858  << G4endl;
1859  }
1860  return false;
1861  }
1862  return true;
1863 }
1864 
1865 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1866 
1868  const G4MaterialCutsCouple *couple,
1869  const G4DynamicParticle* dp,
1870  G4double length)
1871 {
1872  DefineMaterial(couple);
1873  G4double ekin = dp->GetKineticEnergy();
1874  SelectModel(ekin*massRatio);
1876  tmax = std::min(tmax,(*theCuts)[currentCoupleIndex]);
1877  G4double d = 0.0;
1879  if(fm) { d = fm->Dispersion(currentMaterial,dp,tmax,length); }
1880  return d;
1881 }
1882 
1883 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1884 
1886  G4double kineticEnergy, const G4MaterialCutsCouple* couple)
1887 {
1888  // Cross section per volume is calculated
1889  DefineMaterial(couple);
1890  G4double cross = 0.0;
1891  if(theLambdaTable) {
1892  cross = GetLambdaForScaledEnergy(kineticEnergy*massRatio);
1893  } else {
1894  SelectModel(kineticEnergy*massRatio);
1895  cross = biasFactor*(*theDensityFactor)[currentCoupleIndex]
1897  particle, kineticEnergy,
1899  }
1900  if(cross < 0.0) { cross = 0.0; }
1901  return cross;
1902 }
1903 
1904 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1905 
1907 {
1910  G4double x = DBL_MAX;
1911  if(0.0 < preStepLambda) { x = 1.0/preStepLambda; }
1912  return x;
1913 }
1914 
1915 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1916 
1918  G4double x, G4double y,
1919  G4double& z)
1920 {
1921  G4GPILSelection sel;
1922  return AlongStepGetPhysicalInteractionLength(track, x, y, z, &sel);
1923 }
1924 
1925 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1926 
1928  const G4Track& track,
1929  G4double,
1931 
1932 {
1933  *condition = NotForced;
1934  return MeanFreePath(track);
1935 }
1936 
1937 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1938 
1940  const G4Track&,
1942 {
1943  return DBL_MAX;
1944 }
1945 
1946 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1947 
1950  G4double)
1951 {
1952  G4PhysicsVector* v =
1954  v->SetSpline(lManager->SplineFlag());
1955  return v;
1956 }
1957 
1958 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1959 
1962 {
1963  G4bool add = true;
1964  if(p->GetProcessName() != "eBrem") { add = false; }
1965  if(add && nProcesses > 0) {
1966  for(G4int i=0; i<nProcesses; ++i) {
1967  if(p == scProcesses[i]) {
1968  add = false;
1969  break;
1970  }
1971  }
1972  }
1973  if(add) {
1974  scProcesses.push_back(p);
1975  ++nProcesses;
1976  if (1 < verboseLevel) {
1977  G4cout << "### The process " << p->GetProcessName()
1978  << " is added to the list of collaborative processes of "
1979  << GetProcessName() << G4endl;
1980  }
1981  }
1982 }
1983 
1984 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1985 
1986 void
1988 {
1989  G4bool verb = false;
1990  if(fTotal == tType) {
1992  if(p) {
1993  size_t n = p->length();
1994  G4PhysicsVector* pv = (*p)[0];
1995  G4double emax = maxKinEnergyCSDA;
1996 
1997  for (size_t i=0; i<n; ++i) {
1998  G4double dedx = 0.0;
1999  pv = (*p)[i];
2000  if(pv) {
2001  dedx = pv->Value(emax, idxDEDXunRestricted);
2002  } else {
2003  pv = (*p)[(*theDensityIdx)[i]];
2004  if(pv) {
2005  dedx =
2006  pv->Value(emax, idxDEDXunRestricted)*(*theDensityFactor)[i];
2007  }
2008  }
2009  theDEDXAtMaxEnergy[i] = dedx;
2010  //G4cout << "i= " << i << " emax(MeV)= " << emax/MeV<< " dedx= "
2011  //<< dedx << G4endl;
2012  }
2013  }
2014 
2015  } else if(fRestricted == tType) {
2016  if(verb) {
2017  G4cout<< "G4VEnergyLossProcess::SetDEDXTable "
2019  << " oldTable " << theDEDXTable << " newTable " << p
2020  << " ion " << theIonisationTable
2021  << " IsMaster " << isMaster
2022  << " " << GetProcessName() << G4endl;
2023  G4cout << (*p) << G4endl;
2024  }
2025  theDEDXTable = p;
2026  } else if(fSubRestricted == tType) {
2027  theDEDXSubTable = p;
2028  } else if(fIsIonisation == tType) {
2029  if(verb) {
2030  G4cout<< "G4VEnergyLossProcess::SetIonisationTable "
2032  << " oldTable " << theDEDXTable << " newTable " << p
2033  << " ion " << theIonisationTable
2034  << " IsMaster " << isMaster
2035  << " " << GetProcessName() << G4endl;
2036  }
2037  theIonisationTable = p;
2038  } else if(fIsSubIonisation == tType) {
2040  }
2041 }
2042 
2043 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2044 
2046 {
2047  theCSDARangeTable = p;
2048 
2049  if(p) {
2050  size_t n = p->length();
2051  G4PhysicsVector* pv;
2052  G4double emax = maxKinEnergyCSDA;
2053 
2054  for (size_t i=0; i<n; ++i) {
2055  pv = (*p)[i];
2056  G4double rmax = 0.0;
2057  if(pv) { rmax = pv->Value(emax, idxCSDA); }
2058  else {
2059  pv = (*p)[(*theDensityIdx)[i]];
2060  if(pv) { rmax = pv->Value(emax, idxCSDA)/(*theDensityFactor)[i]; }
2061  }
2062  theRangeAtMaxEnergy[i] = rmax;
2063  //G4cout << "i= " << i << " Emax(MeV)= " << emax/MeV << " Rmax= "
2064  //<< rmax<< G4endl;
2065  }
2066  }
2067 }
2068 
2069 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2070 
2072 {
2074  if(1 < verboseLevel) {
2075  G4cout << "### Set Range table " << p
2076  << " for " << particle->GetParticleName()
2077  << " and process " << GetProcessName() << G4endl;
2078  }
2079 }
2080 
2081 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2082 
2084 {
2086  if(1 < verboseLevel) {
2087  G4cout << "### Set SecondaryRange table " << p
2088  << " for " << particle->GetParticleName()
2089  << " and process " << GetProcessName() << G4endl;
2090  }
2091 }
2092 
2093 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2094 
2096 {
2098  if(1 < verboseLevel) {
2099  G4cout << "### Set InverseRange table " << p
2100  << " for " << particle->GetParticleName()
2101  << " and process " << GetProcessName() << G4endl;
2102  }
2103 }
2104 
2105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2106 
2108 {
2109  if(1 < verboseLevel) {
2110  G4cout << "### Set Lambda table " << p
2111  << " for " << particle->GetParticleName()
2112  << " and process " << GetProcessName() << G4endl;
2113  }
2114  theLambdaTable = p;
2115  tablesAreBuilt = true;
2116 
2120 
2121  if(theLambdaTable) {
2122  size_t n = theLambdaTable->length();
2123  G4PhysicsVector* pv = (*theLambdaTable)[0];
2124  G4double e, ss, smax, emax;
2125 
2126  size_t i;
2127 
2128  // first loop on existing vectors
2129  for (i=0; i<n; ++i) {
2130  pv = (*theLambdaTable)[i];
2131  if(pv) {
2132  size_t nb = pv->GetVectorLength();
2133  emax = DBL_MAX;
2134  smax = 0.0;
2135  if(nb > 0) {
2136  for (size_t j=0; j<nb; ++j) {
2137  e = pv->Energy(j);
2138  ss = (*pv)(j);
2139  if(ss > smax) {
2140  smax = ss;
2141  emax = e;
2142  }
2143  }
2144  }
2145  theEnergyOfCrossSectionMax[i] = emax;
2146  theCrossSectionMax[i] = smax;
2147  if(1 < verboseLevel) {
2148  G4cout << "For " << particle->GetParticleName()
2149  << " Max CS at i= " << i << " emax(MeV)= " << emax/MeV
2150  << " lambda= " << smax << G4endl;
2151  }
2152  }
2153  }
2154  // second loop using base materials
2155  for (i=0; i<n; ++i) {
2156  pv = (*theLambdaTable)[i];
2157  if(!pv){
2158  G4int j = (*theDensityIdx)[i];
2160  theCrossSectionMax[i] = (*theDensityFactor)[i]*theCrossSectionMax[j];
2161  }
2162  }
2163  }
2164 }
2165 
2166 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2167 
2169 {
2170  theSubLambdaTable = p;
2171  if(1 < verboseLevel) {
2172  G4cout << "### Set SebLambda table " << p
2173  << " for " << particle->GetParticleName()
2174  << " and process " << GetProcessName() << G4endl;
2175  }
2176 }
2177 
2178 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2179 
2181 {
2182  const G4Element* elm = 0;
2183  if(currentModel) { elm = currentModel->GetCurrentElement(); }
2184  return elm;
2185 }
2186 
2187 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2188 
2190  G4bool flag)
2191 {
2192  if(f > 0.0) {
2193  biasFactor = f;
2194  weightFlag = flag;
2195  if(1 < verboseLevel) {
2196  G4cout << "### SetCrossSectionBiasingFactor: for "
2197  << " process " << GetProcessName()
2198  << " biasFactor= " << f << " weightFlag= " << flag
2199  << G4endl;
2200  }
2201  }
2202 }
2203 
2204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2205 
2206 void
2208  const G4String& region,
2209  G4bool flag)
2210 {
2211  if(!biasManager) { biasManager = new G4EmBiasingManager(); }
2212  if(1 < verboseLevel) {
2213  G4cout << "### ActivateForcedInteraction: for "
2214  << " process " << GetProcessName()
2215  << " length(mm)= " << length/mm
2216  << " in G4Region <" << region
2217  << "> weightFlag= " << flag
2218  << G4endl;
2219  }
2220  weightFlag = flag;
2221  biasManager->ActivateForcedInteraction(length, region);
2222 }
2223 
2224 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2225 
2226 void
2228  G4double factor,
2229  G4double energyLimit)
2230 {
2231  if (0.0 <= factor) {
2232 
2233  // Range cut can be applied only for e-
2234  if(0.0 == factor && secondaryParticle != G4Electron::Electron())
2235  { return; }
2236 
2237  if(!biasManager) { biasManager = new G4EmBiasingManager(); }
2238  biasManager->ActivateSecondaryBiasing(region, factor, energyLimit);
2239  if(1 < verboseLevel) {
2240  G4cout << "### ActivateSecondaryBiasing: for "
2241  << " process " << GetProcessName()
2242  << " factor= " << factor
2243  << " in G4Region <" << region
2244  << "> energyLimit(MeV)= " << energyLimit/MeV
2245  << G4endl;
2246  }
2247  }
2248 }
2249 
2250 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2251 
G4VEnergyLossProcess(const G4String &name="EnergyLoss", G4ProcessType type=fElectromagnetic)
static const double cm
Definition: G4SIunits.hh:106
G4double condition(const G4ErrorSymMatrix &m)
G4double MaxSecondaryKinEnergy(const G4DynamicParticle *dynParticle)
Definition: G4VEmModel.hh:448
G4ParticleDefinition * GetDefinition() const
const G4VProcess * GetMasterProcess() const
Definition: G4VProcess.hh:538
G4LossTableManager * lManager
std::vector< G4double > theEnergyOfCrossSectionMax
void SetDEDXTable(G4PhysicsTable *p, G4EmTableType tType)
G4int GetParentID() const
G4bool SplineFlag() const
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:245
static const double MeV
Definition: G4SIunits.hh:193
G4Region * GetRegion(const G4String &name, G4bool verbose=true) const
void InitializeForPostStep(const G4Track &)
const G4String & GetName() const
G4SafetyHelper * GetSafetyHelper() const
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
const std::vector< G4double > * GetDensityFactors()
static void SetPhysicsVector(G4PhysicsTable *physTable, size_t idx, G4PhysicsVector *vec)
void InitialiseHelper()
static G4LossTableManager * Instance()
void ActivateForcedInteraction(G4double length=0.0, const G4String &r="")
G4GPILSelection aGPILSelection
virtual G4double SampleFluctuations(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmax, G4double length, G4double meanLoss)=0
void AlongStepDeexcitation(std::vector< G4Track * > &tracks, const G4Step &step, G4double &eLoss, G4int coupleIndex)
G4double ApplySecondaryBiasing(std::vector< G4DynamicParticle * > &, const G4Track &track, G4VEmModel *currentModel, G4ParticleChangeForGamma *pParticleChange, G4double &eloss, G4int coupleIdx, G4double tcut, G4double safety=0.0)
G4int verboseLevel
Definition: G4VProcess.hh:368
virtual void CorrectionsAlongStep(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &eloss, G4double &niel, G4double length)
Definition: G4VEmModel.cc:339
void UpdateEmModel(const G4String &, G4double, G4double)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin=0.0, G4double tmax=DBL_MAX)=0
G4double GetKineticEnergy() const
void PrintInfoDefinition(const G4ParticleDefinition &part)
CLHEP::Hep3Vector G4ThreeVector
G4double GetStepLength() const
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:592
void BuildPhysicsTable(const G4ParticleDefinition &)
const G4DynamicParticle * GetDynamicParticle() const
G4bool SecondaryBiasingRegion(G4int coupleIdx)
G4bool IsPIXEActive() const
G4PhysicsTable * SubLambdaTable() const
G4double z
Definition: TRTMaterials.hh:39
void SetLowEnergyLimit(G4double elimit)
G4VEmModel * GetModel(G4int, G4bool ver=false)
G4double GetProductionCut(G4int index) const
void DeRegister(G4VEnergyLossProcess *p)
G4PhysicsTable * RangeTableForLoss() const
G4double GetContinuousStepLimit(const G4Track &track, G4double previousStepSize, G4double currentMinimumStep, G4double &currentSafety)
G4String name
Definition: TRTMaterials.hh:40
const std::vector< G4double > * theDensityFactor
G4StepStatus GetStepStatus() const
G4bool BuildCSDARange() const
const G4String & GetName() const
Definition: G4Material.hh:176
G4PhysicsTable * SecondaryRangeTable() const
G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
void AddCollaborativeProcess(G4VEnergyLossProcess *)
const G4ThreeVector & GetPosition() const
const std::vector< G4int > * theDensityIdx
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel)
Definition: G4VEmModel.cc:206
void ActivateForcedInteraction(G4double length=0.0, const G4String &region="", G4bool flag=true)
const G4ParticleDefinition * secondaryParticle
G4double GetMeanFreePath(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
void SetStepFunction(G4double v1, G4double v2)
G4bool ForcedInteractionRegion(G4int coupleIdx)
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
G4PhysicsTable * CSDARangeTable() const
G4PhysicsTable * theSecondaryRangeTable
G4bool StorePhysicsTable(const G4ParticleDefinition *, const G4String &directory, G4bool ascii=false)
G4VEmFluctuationModel * fluctModel
G4ParticleChangeForLoss fParticleChange
G4double theNumberOfInteractionLengthLeft
Definition: G4VProcess.hh:293
void SelectModel(G4double kinEnergy)
G4bool GetFlag(size_t idx) const
G4PhysicsTable * theIonisationTable
void SetTouchableHandle(const G4TouchableHandle &apValue)
void AddEmModel(G4int, G4VEmModel *, G4VEmFluctuationModel *, const G4Region *)
G4VEmModel * GetModelByIndex(G4int idx=0, G4bool ver=false) const
G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
G4PhysicsTable * IonisationTableForSubsec() const
G4double GetParentWeight() const
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
void ActivateSecondaryBiasing(const G4String &region, G4double factor, G4double energyLimit)
G4VEmFluctuationModel * GetModelOfFluctuations()
Definition: G4VEmModel.hh:571
const G4DataVector * SubCutoff() const
size_t GetVectorLength() const
G4PhysicsVector * LambdaPhysicsVector(const G4MaterialCutsCouple *, G4double cut)
void InitializeForAlongStep(const G4Track &)
G4double GetProposedKineticEnergy() const
G4PhysicsTable * IonisationTable() const
void FillLambdaVector(G4PhysicsVector *, const G4MaterialCutsCouple *, G4bool startFromNull=true, G4EmTableType t=fRestricted)
G4ProcessManager * GetProcessManager() const
G4VAtomDeexcitation * atomDeexcitation
int G4int
Definition: G4Types.hh:78
void FillSecondDerivatives()
G4PhysicsTable * BuildLambdaTable(G4EmTableType tType=fRestricted)
const G4String & GetPhysicsTableFileName(const G4ParticleDefinition *, const G4String &directory, const G4String &tableName, G4bool ascii=false)
Definition: G4VProcess.cc:186
static G4PhysicsTable * PreparePhysicsTable(G4PhysicsTable *physTable)
const G4int smax
const G4String & GetParticleName() const
G4LossTableBuilder * GetTableBuilder()
void SetWeight(G4double aValue)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetCreatorModelIndex(G4int idx)
Definition: G4Track.hh:245
G4bool ExistPhysicsTable(const G4String &fileName) const
const G4ParticleDefinition * theGamma
static const G4double reg
const G4MaterialCutsCouple * currentCouple
void SetSpline(G4bool)
void Initialise(const G4ParticleDefinition &part, const G4String &procName, G4int verbose)
const G4DataVector * theCuts
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:683
virtual G4double Dispersion(const G4Material *, const G4DynamicParticle *, G4double tmax, G4double length)=0
static G4RegionStore * GetInstance()
G4PhysicsTable * theSubLambdaTable
G4StepPoint * GetPreStepPoint() const
G4PhysicsTable * LambdaTable() const
void ProposeWeight(G4double finalWeight)
void SetSecondaryWeightByProcess(G4bool)
void SetInverseRangeTable(G4PhysicsTable *p)
void SetCrossSectionBiasingFactor(G4double f, G4bool flag=true)
G4ProcessVector * GetAtRestProcessVector(G4ProcessVectorTypeIndex typ=typeGPIL) const
G4double GetKineticEnergy() const
G4double SampleSubCutSecondaries(std::vector< G4Track * > &, const G4Step &, G4VEmModel *model, G4int matIdx)
const G4ParticleDefinition * baseParticle
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
void FillDEDXVector(G4PhysicsVector *, const G4MaterialCutsCouple *, G4EmTableType t=fRestricted)
void FillSecondariesAlongStep(G4double &eloss, G4double &weight)
const G4DataVector * theSubCuts
G4PhysicsTable * DEDXTable() const
G4double ScaledKinEnergyForLoss(G4double range)
G4double GetDEDXForScaledEnergy(G4double scaledKinEnergy)
void SetParticleChange(G4VParticleChange *, G4VEmFluctuationModel *f=0)
Definition: G4VEmModel.cc:387
const G4ThreeVector & GetPosition() const
const G4Element * GetCurrentElement() const
const G4ParticleDefinition * theGenericIon
const G4ParticleDefinition * thePositron
void ComputeLambdaForScaledEnergy(G4double scaledKinEnergy)
bool G4bool
Definition: G4Types.hh:79
std::vector< G4Track * > scTracks
G4double GetStepLimit(G4int coupleIdx, G4double previousStep)
G4double currentInteractionLength
Definition: G4VProcess.hh:297
void SetProposedKineticEnergy(G4double proposedKinEnergy)
G4double GetSubDEDXForScaledEnergy(G4double scaledKinEnergy)
const G4ParticleDefinition * GetParticleDefinition() const
G4PhysicsTable * theCSDARangeTable
static const double GeV
Definition: G4SIunits.hh:196
const G4String & GetParticleType() const
void Register(G4VEnergyLossProcess *p)
std::vector< G4double > theDEDXAtMaxEnergy
G4VParticleChange * AlongStepDoIt(const G4Track &, const G4Step &)
Definition: G4Step.hh:76
void AddEmModel(G4int, G4VEmModel *, G4VEmFluctuationModel *fluc=0, const G4Region *region=0)
const G4ParticleDefinition * particle
const G4int n
G4double GetGlobalTime() const
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4PhysicsTable * DEDXTableForSubsec() const
G4double Energy(size_t index) const
G4double Value(G4double theEnergy, size_t &lastidx) const
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
G4double AlongStepGetPhysicalInteractionLength(const G4Track &, G4double previousStepSize, G4double currentMinimumStep, G4double &currentSafety, G4GPILSelection *selection)
void AddSecondary(G4Track *aSecondary)
virtual G4double ChargeSquareRatio(const G4Track &)
Definition: G4VEmModel.cc:314
const G4TouchableHandle & GetTouchableHandle() const
G4Material * GetMaterial() const
void RegisterExtraParticle(const G4ParticleDefinition *aParticle, G4VEnergyLossProcess *p)
G4bool IsActive(G4double kinEnergy)
Definition: G4VEmModel.hh:711
void DumpModelList(G4int verb)
std::vector< const G4Region * > scoffRegions
const G4DataVector * Initialise(const G4ParticleDefinition *part, const G4ParticleDefinition *secPart, G4double minSubRange, G4int verb)
G4double ComputeSafety(const G4ThreeVector &pGlobalPoint, G4double maxRadius=DBL_MAX)
G4PhysicsTable * theLambdaTable
static G4TransportationManager * GetTransportationManager()
G4bool RetrievePhysicsTable(const G4ParticleDefinition *, const G4String &directory, G4bool ascii)
void SetMasterThread(G4bool val)
Definition: G4VEmModel.hh:669
void PreparePhysicsTable(const G4ParticleDefinition &)
size_t length() const
G4double G4Log(G4double x)
Definition: G4Log.hh:230
const G4ParticleDefinition * theElectron
G4int size() const
void ActivateSecondaryBiasing(const G4String &region, G4double factor, G4double energyLimit)
static const double perMillion
Definition: G4SIunits.hh:298
virtual G4double GetParticleCharge(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4VEmModel.cc:331
static G4ProductionCutsTable * GetProductionCutsTable()
const G4Material * currentMaterial
void SetLambdaTable(G4PhysicsTable *p)
static const double eV
Definition: G4SIunits.hh:194
#define fm
G4PhysicsTable * theInverseRangeTable
static G4Positron * Positron()
Definition: G4Positron.cc:94
G4int NumberOfModels() const
G4SafetyHelper * safetyHelper
std::vector< G4double > theCrossSectionMax
virtual void PrintInfo()=0
G4PhysicsTable * InverseRangeTable() const
void BuildPhysicsTable(const G4ParticleDefinition *aParticle)
G4double GetPDGMass() const
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
static G4ParticleTable * GetParticleTable()
int G4lrint(double ad)
Definition: templates.hh:163
void InitialiseBaseMaterials(G4PhysicsTable *table)
void SetNumberOfSecondaries(G4int totSecondaries)
G4double GetLocalEnergyDeposit() const
G4double MeanFreePath(const G4Track &track)
G4PhysicsTable * theDEDXunRestrictedTable
G4StepPoint * GetPostStepPoint() const
void UpdateEmModel(const G4String &, G4double, G4double)
G4bool StoreTable(const G4ParticleDefinition *p, G4PhysicsTable *, G4bool ascii, const G4String &directory, const G4String &tname)
G4VParticleChange * pParticleChange
Definition: G4VProcess.hh:283
G4PhysicsTable * theRangeTableForLoss
void SetEmModel(G4VEmModel *, G4int index=1)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double ContinuousStepLimit(const G4Track &track, G4double previousStepSize, G4double currentMinimumStep, G4double &currentSafety)
std::vector< G4VEmModel * > emModels
G4double GetSafety() const
void SetProposedCharge(G4double theCharge)
void PreparePhysicsTable(const G4ParticleDefinition *aParticle, G4VEnergyLossProcess *p, G4bool theMaster)
void SetCSDARangeTable(G4PhysicsTable *pRange)
G4double GetDEDXDispersion(const G4MaterialCutsCouple *couple, const G4DynamicParticle *dp, G4double length)
void LocalPhysicsTables(const G4ParticleDefinition *aParticle, G4VEnergyLossProcess *p)
G4EmTableType
G4double GetLambdaForScaledEnergy(G4double scaledKinEnergy)
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
G4double GetGlobalTime() const
void SetSecondaryRangeTable(G4PhysicsTable *p)
static const double TeV
Definition: G4SIunits.hh:197
G4bool StorePhysicsTable(const G4String &filename, G4bool ascii=false)
G4TrackStatus GetTrackStatus() const
G4VAtomDeexcitation * AtomDeexcitation()
static const double keV
Definition: G4SIunits.hh:195
G4PhysicsTable * DEDXunRestrictedTable() const
std::vector< G4double > theRangeAtMaxEnergy
G4bool RetrieveTable(const G4ParticleDefinition *p, G4PhysicsTable *, G4bool ascii, const G4String &directory, const G4String &tname, G4bool mandatory)
G4PhysicsTable * theDEDXTable
G4double theInitialNumberOfInteractionLength
Definition: G4VProcess.hh:300
G4bool GetFlag(size_t i) const
static G4int Register(G4String &)
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
void SetSubLambdaTable(G4PhysicsTable *p)
virtual G4double MinPrimaryEnergy(const G4ParticleDefinition *, const G4Material *, G4double cut)
G4ProcessVector * GetAlongStepProcessVector(G4ProcessVectorTypeIndex typ=typeGPIL) const
G4ForceCondition
G4Track * GetTrack() const
G4double GetPDGCharge() const
G4EmBiasingManager * biasManager
void ActivateSubCutoff(G4bool val, const G4Region *region=0)
void SetRangeTableForLoss(G4PhysicsTable *p)
static G4bool RetrievePhysicsTable(G4PhysicsTable *physTable, const G4String &fileName, G4bool ascii)
G4ProductionCuts * GetProductionCuts() const
G4PhysicsTable * theDEDXSubTable
#define DBL_MAX
Definition: templates.hh:83
static const double mm
Definition: G4SIunits.hh:102
G4int GetProcessSubType() const
Definition: G4VProcess.hh:426
G4double GetScaledRangeForScaledEnergy(G4double scaledKinEnergy)
std::vector< G4VEnergyLossProcess * > scProcesses
G4PhysicsTable * BuildDEDXTable(G4EmTableType tType=fRestricted)
G4VEmModel * EmModel(G4int index=1) const
void DefineMaterial(const G4MaterialCutsCouple *couple)
void clearAndDestroy()
void SetVerboseLevel(G4int value)
Definition: G4VProcess.hh:437
virtual void InitialiseEnergyLossProcess(const G4ParticleDefinition *, const G4ParticleDefinition *)=0
G4GPILSelection
const G4Material * GetMaterial() const
G4PhysicsTable * theIonisationSubTable
G4double CrossSectionPerVolume(G4double kineticEnergy, const G4MaterialCutsCouple *couple)
std::vector< G4DynamicParticle * > secParticles
const std::vector< G4int > * GetCoupleIndexes()
G4ProcessType
const G4Element * GetCurrentElement() const
Definition: G4VEmModel.hh:440
G4bool IsIonisationProcess() const
G4EmModelManager * modelManager