Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4VEnergyLossProcess.cc
이 파일의 문서화 페이지로 가기
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 107959 2017-12-14 13:05:59Z 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 "G4VSubCutProducer.hh"
161 #include "G4EmBiasingManager.hh"
162 #include "G4Log.hh"
163 #include <iostream>
164 
165 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
166 
168  G4ProcessType type):
169  G4VContinuousDiscreteProcess(name, type),
170  secondaryParticle(nullptr),
171  nSCoffRegions(0),
172  idxSCoffRegions(nullptr),
173  nProcesses(0),
174  theDEDXTable(nullptr),
175  theDEDXSubTable(nullptr),
176  theDEDXunRestrictedTable(nullptr),
177  theIonisationTable(nullptr),
178  theIonisationSubTable(nullptr),
179  theRangeTableForLoss(nullptr),
180  theCSDARangeTable(nullptr),
181  theSecondaryRangeTable(nullptr),
182  theInverseRangeTable(nullptr),
183  theLambdaTable(nullptr),
184  theSubLambdaTable(nullptr),
185  theDensityFactor(nullptr),
186  theDensityIdx(nullptr),
187  baseParticle(nullptr),
188  lossFluctuationFlag(true),
189  rndmStepFlag(false),
190  tablesAreBuilt(false),
191  integral(true),
192  isIon(false),
193  isIonisation(true),
194  useSubCutoff(false),
195  useDeexcitation(false),
196  particle(nullptr),
197  currentCouple(nullptr),
198  mfpKinEnergy(0.0)
199 {
201  SetVerboseLevel(1);
202 
203  // low energy limit
205  preStepKinEnergy = 0.0;
206  preStepRangeEnergy = 0.0;
208 
209  // Size of tables assuming spline
210  minKinEnergy = 0.1*keV;
211  maxKinEnergy = 100.0*TeV;
212  nBins = 84;
213  maxKinEnergyCSDA = 1.0*GeV;
214  nBinsCSDA = 35;
216  = actLossFluc = actIntegral = actStepFunc = false;
217 
218  // default linear loss limit for spline
219  linLossLimit = 0.01;
220  dRoverRange = 0.2;
222 
223  // default lambda factor
224  lambdaFactor = 0.8;
225 
226  // cross section biasing
227  biasFactor = 1.0;
228 
229  // particle types
233  theGenericIon = nullptr;
234 
235  // run time objects
240  ->GetSafetyHelper();
242 
243  // initialise model
245  lManager->Register(this);
246  fluctModel = nullptr;
247  currentModel = nullptr;
248  atomDeexcitation = nullptr;
249  subcutProducer = nullptr;
250 
251  biasManager = nullptr;
252  biasFlag = false;
253  weightFlag = false;
254  isMaster = true;
255  lastIdx = 0;
256 
260 
261  scTracks.reserve(5);
262  secParticles.reserve(5);
263 
264  theCuts = theSubCuts = nullptr;
265  currentMaterial = nullptr;
269 
270  secID = biasID = subsecID = -1;
271 }
272 
273 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
274 
276 {
277  /*
278  G4cout << "** G4VEnergyLossProcess::~G4VEnergyLossProcess() for "
279  << GetProcessName() << " isMaster: " << isMaster
280  << " basePart: " << baseParticle
281  << G4endl;
282  */
283  Clean();
284 
285  // G4cout << " isIonisation " << isIonisation << " "
286  // << theDEDXTable << " " << theIonisationTable << G4endl;
287 
288  if (isMaster && !baseParticle) {
289  if(theDEDXTable) {
290 
291  //G4cout << " theIonisationTable " << theIonisationTable << G4endl;
293  //G4cout << " delete theDEDXTable " << theDEDXTable << G4endl;
295  delete theDEDXTable;
296  theDEDXTable = nullptr;
297  if(theDEDXSubTable) {
299  { theIonisationSubTable = nullptr; }
301  delete theDEDXSubTable;
302  theDEDXSubTable = nullptr;
303  }
304  }
305  //G4cout << " theIonisationTable " << theIonisationTable << G4endl;
306  if(theIonisationTable) {
307  //G4cout << " delete theIonisationTable " << theIonisationTable << G4endl;
309  delete theIonisationTable;
310  theIonisationTable = nullptr;
311  }
314  delete theIonisationSubTable;
315  theIonisationSubTable = nullptr;
316  }
320  theDEDXunRestrictedTable = nullptr;
321  }
324  delete theCSDARangeTable;
325  theCSDARangeTable = nullptr;
326  }
327  //G4cout << "delete RangeTable: " << theRangeTableForLoss << G4endl;
330  delete theRangeTableForLoss;
331  theRangeTableForLoss = nullptr;
332  }
333  //G4cout << "delete InvRangeTable: " << theInverseRangeTable << G4endl;
334  if(theInverseRangeTable && isIonisation /*&& !isIon*/) {
336  delete theInverseRangeTable;
337  theInverseRangeTable = nullptr;
338  }
339  //G4cout << "delete LambdaTable: " << theLambdaTable << G4endl;
340  if(theLambdaTable) {
342  delete theLambdaTable;
343  theLambdaTable = nullptr;
344  }
345  if(theSubLambdaTable) {
347  delete theSubLambdaTable;
348  theSubLambdaTable = nullptr;
349  }
350  }
351 
352  delete modelManager;
353  delete biasManager;
354  lManager->DeRegister(this);
355  //G4cout << "** all removed" << G4endl;
356 }
357 
358 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
359 
361 {
362  /*
363  if(1 < verboseLevel) {
364  G4cout << "G4VEnergyLossProcess::Clear() for " << GetProcessName()
365  << G4endl;
366  }
367  */
368  delete [] idxSCoffRegions;
369 
370  tablesAreBuilt = false;
371 
372  scProcesses.clear();
373  nProcesses = 0;
374 
378 }
379 
380 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
381 
383  const G4Material*,
384  G4double cut)
385 {
386  return cut;
387 }
388 
389 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
390 
392  G4VEmFluctuationModel* fluc,
393  const G4Region* region)
394 {
395  modelManager->AddEmModel(order, p, fluc, region);
396  if(p) { p->SetParticleChange(pParticleChange, fluc); }
397 }
398 
399 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
400 
402  G4double emin, G4double emax)
403 {
404  modelManager->UpdateEmModel(nam, emin, emax);
405 }
406 
407 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
408 
410 {
411  for(auto & em : emModels) { if(em == ptr) { return; } }
412  emModels.push_back(ptr);
413 }
414 
415 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
416 
418 {
419  return (index < emModels.size()) ? emModels[index] : nullptr;
420 }
421 
422 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
423 
425 {
426  return modelManager->GetModel(idx, ver);
427 }
428 
429 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
430 
432 {
433  return modelManager->NumberOfModels();
434 }
435 
436 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
437 
438 void
440 {
441  if(1 < verboseLevel) {
442  G4cout << "G4VEnergyLossProcess::PreparePhysicsTable for "
443  << GetProcessName() << " for " << part.GetParticleName()
444  << " " << this << G4endl;
445  }
447 
448  currentCouple = nullptr;
449  preStepLambda = 0.0;
451  fRange = DBL_MAX;
452  preStepKinEnergy = 0.0;
453  preStepRangeEnergy = 0.0;
454  chargeSqRatio = 1.0;
455  massRatio = 1.0;
456  reduceFactor = 1.0;
457  fFactor = 1.0;
458  lastIdx = 0;
459 
460  // Are particle defined?
461  if( !particle ) { particle = &part; }
462 
463  if(part.GetParticleType() == "nucleus") {
464 
465  G4String pname = part.GetParticleName();
466  if(pname != "deuteron" && pname != "triton" &&
467  pname != "alpha+" && pname != "helium" &&
468  pname != "hydrogen") {
469 
470  if(!theGenericIon) {
471  theGenericIon =
473  }
474  isIon = true;
478  size_t n = v->size();
479  for(size_t j=0; j<n; ++j) {
480  if((*v)[j] == this) {
482  break;
483  }
484  }
485  }
486  }
487  }
488 
489  if( particle != &part ) {
490  if(!isIon) {
491  lManager->RegisterExtraParticle(&part, this);
492  }
493  if(1 < verboseLevel) {
494  G4cout << "### G4VEnergyLossProcess::PreparePhysicsTable()"
495  << " interrupted for "
496  << part.GetParticleName() << " isIon= " << isIon
497  << " particle " << particle << " GenericIon " << theGenericIon
498  << G4endl;
499  }
500  return;
501  }
502 
503  Clean();
504  lManager->PreparePhysicsTable(&part, this, isMaster);
506 
507  // Base particle and set of models can be defined here
509 
510  const G4ProductionCutsTable* theCoupleTable=
512  size_t n = theCoupleTable->GetTableSize();
513 
514  theDEDXAtMaxEnergy.resize(n, 0.0);
515  theRangeAtMaxEnergy.resize(n, 0.0);
516  theEnergyOfCrossSectionMax.resize(n, 0.0);
517  theCrossSectionMax.resize(n, DBL_MAX);
518 
519  // parameters of the process
525  if(!actBinning) {
527  *G4lrint(std::log10(maxKinEnergy/minKinEnergy));
528  }
531  *G4lrint(std::log10(maxKinEnergyCSDA/minKinEnergy));
536 
537  G4bool isElec = true;
538  if(particle->GetPDGMass() > CLHEP::MeV) { isElec = false; }
539  theParameters->DefineRegParamForLoss(this, isElec);
540 
541  G4double initialCharge = particle->GetPDGCharge();
542  G4double initialMass = particle->GetPDGMass();
543 
544  if (baseParticle) {
545  massRatio = (baseParticle->GetPDGMass())/initialMass;
546  G4double q = initialCharge/baseParticle->GetPDGCharge();
547  chargeSqRatio = q*q;
548  if(chargeSqRatio > 0.0) { reduceFactor = 1.0/(chargeSqRatio*massRatio); }
549  }
550  if(initialMass < MeV) {
552  } else {
554  }
555 
556  // Tables preparation
557  if (isMaster && !baseParticle) {
558 
559  if(theDEDXTable && isIonisation) {
562  delete theDEDXTable;
564  }
568  delete theDEDXSubTable;
570  }
571  }
572 
575 
576  if(theDEDXSubTable) {
577  theDEDXSubTable =
579  }
580 
581  if (theParameters->BuildCSDARange()) {
586  }
587 
589 
590  if(isIonisation) {
595  }
596 
598  theDEDXSubTable =
602  }
603  }
604  /*
605  G4cout << "** G4VEnergyLossProcess::PreparePhysicsTable() for "
606  << GetProcessName() << " and " << particle->GetParticleName()
607  << " isMaster: " << isMaster << " isIonisation: "
608  << isIonisation << G4endl;
609  G4cout << " theDEDX: " << theDEDXTable
610  << " theRange: " << theRangeTableForLoss
611  << " theInverse: " << theInverseRangeTable
612  << " theLambda: " << theLambdaTable << G4endl;
613  */
614  // forced biasing
615  if(biasManager) {
617  biasFlag = false;
618  }
619 
620  // defined ID of secondary particles
621  if(isMaster) {
622  G4String nam1 = GetProcessName();
623  G4String nam4 = nam1 + "_split";
624  G4String nam5 = nam1 + "_subcut";
628  }
629 
630  // initialisation of models
632  for(G4int i=0; i<nmod; ++i) {
633  G4VEmModel* mod = modelManager->GetModel(i);
637  if(mod->HighEnergyLimit() > maxKinEnergy) {
639  }
640  }
643  verboseLevel);
644 
645  // Sub Cutoff
646  if(nSCoffRegions > 0) {
647  if(theParameters->MinSubRange() < 1.0) { useSubCutoff = true; }
648 
650 
651  idxSCoffRegions = new G4bool[n];
652  for (size_t j=0; j<n; ++j) {
653 
654  const G4MaterialCutsCouple* couple =
655  theCoupleTable->GetMaterialCutsCouple(j);
656  const G4ProductionCuts* pcuts = couple->GetProductionCuts();
657 
658  G4bool reg = false;
659  for(G4int i=0; i<nSCoffRegions; ++i) {
660  if( pcuts == scoffRegions[i]->GetProductionCuts()) {
661  reg = true;
662  break;
663  }
664  }
665  idxSCoffRegions[j] = reg;
666  }
667  }
668 
669  if(1 < verboseLevel) {
670  G4cout << "G4VEnergyLossProcess::PrepearPhysicsTable() is done "
671  << " for local " << particle->GetParticleName()
672  << " isIon= " << isIon;
673  if(baseParticle) {
674  G4cout << "; base: " << baseParticle->GetParticleName();
675  }
676  G4cout << " chargeSqRatio= " << chargeSqRatio
677  << " massRatio= " << massRatio
678  << " reduceFactor= " << reduceFactor << G4endl;
679  if (nSCoffRegions) {
680  G4cout << " SubCutoff Regime is ON for regions: " << G4endl;
681  for (G4int i=0; i<nSCoffRegions; ++i) {
682  const G4Region* r = scoffRegions[i];
683  G4cout << " " << r->GetName() << G4endl;
684  }
685  }
686  }
687 }
688 
689 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
690 
692 {
693  if(1 < verboseLevel) {
694  G4cout << "### G4VEnergyLossProcess::BuildPhysicsTable() for "
695  << GetProcessName()
696  << " and particle " << part.GetParticleName()
697  << "; local: " << particle->GetParticleName();
698  if(baseParticle) {
699  G4cout << "; base: " << baseParticle->GetParticleName();
700  }
701  G4cout << " TablesAreBuilt= " << tablesAreBuilt
702  << " isIon= " << isIon << " " << this << G4endl;
703  }
704 
705  if(&part == particle) {
706 
708  if(isMaster) {
712 
713  } else {
714 
715  const G4VEnergyLossProcess* masterProcess =
716  static_cast<const G4VEnergyLossProcess*>(GetMasterProcess());
717 
718  // define density factors for worker thread
719  bld->InitialiseBaseMaterials(masterProcess->DEDXTable());
722 
723  // copy table pointers from master thread
724  SetDEDXTable(masterProcess->DEDXTable(),fRestricted);
726  SetDEDXTable(masterProcess->DEDXunRestrictedTable(),fTotal);
727  SetDEDXTable(masterProcess->IonisationTable(),fIsIonisation);
729  SetRangeTableForLoss(masterProcess->RangeTableForLoss());
730  SetCSDARangeTable(masterProcess->CSDARangeTable());
732  SetInverseRangeTable(masterProcess->InverseRangeTable());
733  SetLambdaTable(masterProcess->LambdaTable());
734  SetSubLambdaTable(masterProcess->SubLambdaTable());
735  isIonisation = masterProcess->IsIonisationProcess();
736 
737  tablesAreBuilt = true;
738  // local initialisation of models
739  G4bool printing = true;
740  G4int numberOfModels = modelManager->NumberOfModels();
741  for(G4int i=0; i<numberOfModels; ++i) {
742  G4VEmModel* mod = GetModelByIndex(i, printing);
743  G4VEmModel* mod0= masterProcess->GetModelByIndex(i,printing);
744  mod->InitialiseLocal(particle, mod0);
745  }
746 
748  }
749 
750  // needs to be done only once
752  }
753  // explicitly defined printout by particle name
754  G4String num = part.GetParticleName();
755  if(1 < verboseLevel ||
756  (0 < verboseLevel && (num == "e-" ||
757  num == "e+" || num == "mu+" ||
758  num == "mu-" || num == "proton"||
759  num == "pi+" || num == "pi-" ||
760  num == "kaon+" || num == "kaon-" ||
761  num == "alpha" || num == "anti_proton" ||
762  num == "GenericIon"|| num == "alpha++" ||
763  num == "alpha+" )))
764  {
765  StreamInfo(G4cout, part);
766  }
767 
768  // Added tracking cut to avoid tracking artifacts
769  // identify deexcitation flag
770  if(isIonisation) {
773  if(atomDeexcitation) {
774  if(atomDeexcitation->IsPIXEActive()) { useDeexcitation = true; }
775  }
776  }
777  /*
778  G4cout << "** G4VEnergyLossProcess::BuildPhysicsTable() for "
779  << GetProcessName() << " and " << particle->GetParticleName()
780  << " isMaster: " << isMaster << " isIonisation: "
781  << isIonisation << G4endl;
782  G4cout << " theDEDX: " << theDEDXTable
783  << " theRange: " << theRangeTableForLoss
784  << " theInverse: " << theInverseRangeTable
785  << " theLambda: " << theLambdaTable << G4endl;
786  */
787  //if(1 < verboseLevel || verb) {
788  if(1 < verboseLevel) {
789  G4cout << "### G4VEnergyLossProcess::BuildPhysicsTable() done for "
790  << GetProcessName()
791  << " and particle " << part.GetParticleName();
792  if(isIonisation) { G4cout << " isIonisation flag = 1"; }
793  G4cout << G4endl;
794  }
795 }
796 
797 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
798 
800 {
801  if(1 < verboseLevel ) {
802  G4cout << "G4VEnergyLossProcess::BuildDEDXTable() of type " << tType
803  << " for " << GetProcessName()
804  << " and particle " << particle->GetParticleName()
805  << G4endl;
806  }
807  G4PhysicsTable* table = nullptr;
809  G4int bin = nBins;
810 
811  if(fTotal == tType) {
812  emax = maxKinEnergyCSDA;
813  bin = nBinsCSDA;
814  table = theDEDXunRestrictedTable;
815  } else if(fRestricted == tType) {
816  table = theDEDXTable;
817  } else if(fSubRestricted == tType) {
818  table = theDEDXSubTable;
819  } else {
820  G4cout << "G4VEnergyLossProcess::BuildDEDXTable WARNING: wrong type "
821  << tType << G4endl;
822  }
823 
824  // Access to materials
825  const G4ProductionCutsTable* theCoupleTable=
827  size_t numOfCouples = theCoupleTable->GetTableSize();
828 
829  if(1 < verboseLevel) {
830  G4cout << numOfCouples << " materials"
831  << " minKinEnergy= " << minKinEnergy
832  << " maxKinEnergy= " << emax
833  << " nbin= " << bin
834  << " EmTableType= " << tType
835  << " table= " << table << " " << this
836  << G4endl;
837  }
838  if(!table) { return table; }
839 
841  G4bool splineFlag = theParameters->Spline();
842  G4PhysicsLogVector* aVector = nullptr;
843  G4PhysicsLogVector* bVector = nullptr;
844 
845  for(size_t i=0; i<numOfCouples; ++i) {
846 
847  if(1 < verboseLevel) {
848  G4cout << "G4VEnergyLossProcess::BuildDEDXVector Idx= " << i
849  << " flagTable= " << table->GetFlag(i)
850  << " Flag= " << bld->GetFlag(i) << G4endl;
851  }
852  if(bld->GetFlag(i)) {
853 
854  // create physics vector and fill it
855  const G4MaterialCutsCouple* couple =
856  theCoupleTable->GetMaterialCutsCouple(i);
857  if((*table)[i]) { delete (*table)[i]; }
858  if(bVector) {
859  aVector = new G4PhysicsLogVector(*bVector);
860  } else {
861  bVector = new G4PhysicsLogVector(minKinEnergy, emax, bin);
862  aVector = bVector;
863  }
864  aVector->SetSpline(splineFlag);
865 
866  modelManager->FillDEDXVector(aVector, couple, tType);
867  if(splineFlag) { aVector->FillSecondDerivatives(); }
868 
869  // Insert vector for this material into the table
870  G4PhysicsTableHelper::SetPhysicsVector(table, i, aVector);
871  }
872  }
873 
874  if(1 < verboseLevel) {
875  G4cout << "G4VEnergyLossProcess::BuildDEDXTable(): table is built for "
877  << " and process " << GetProcessName()
878  << G4endl;
879  if(2 < verboseLevel) G4cout << (*table) << G4endl;
880  }
881 
882  return table;
883 }
884 
885 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
886 
888 {
889  G4PhysicsTable* table = nullptr;
890 
891  if(fRestricted == tType) {
892  table = theLambdaTable;
893  } else if(fSubRestricted == tType) {
894  table = theSubLambdaTable;
895  } else {
896  G4cout << "G4VEnergyLossProcess::BuildLambdaTable WARNING: wrong type "
897  << tType << G4endl;
898  }
899 
900  if(1 < verboseLevel) {
901  G4cout << "G4VEnergyLossProcess::BuildLambdaTable() of type "
902  << tType << " for process "
903  << GetProcessName() << " and particle "
905  << " EmTableType= " << tType
906  << " table= " << table
907  << G4endl;
908  }
909  if(!table) {return table;}
910 
911  // Access to materials
912  const G4ProductionCutsTable* theCoupleTable=
914  size_t numOfCouples = theCoupleTable->GetTableSize();
915 
919 
920  G4bool splineFlag = theParameters->Spline();
921  G4PhysicsLogVector* aVector = nullptr;
923 
924  for(size_t i=0; i<numOfCouples; ++i) {
925 
926  if (bld->GetFlag(i)) {
927 
928  // create physics vector and fill it
929  const G4MaterialCutsCouple* couple =
930  theCoupleTable->GetMaterialCutsCouple(i);
931  delete (*table)[i];
932 
933  G4bool startNull = true;
934  G4double emin =
935  MinPrimaryEnergy(particle,couple->GetMaterial(),(*theCuts)[i]);
936  if(minKinEnergy > emin) {
937  emin = minKinEnergy;
938  startNull = false;
939  }
940 
942  if(emax <= emin) { emax = 2*emin; }
943  G4int bin = G4lrint(nBins*G4Log(emax/emin)/scale);
944  bin = std::max(bin, 3);
945  aVector = new G4PhysicsLogVector(emin, emax, bin);
946  aVector->SetSpline(splineFlag);
947 
948  modelManager->FillLambdaVector(aVector, couple, startNull, tType);
949  if(splineFlag) { aVector->FillSecondDerivatives(); }
950 
951  // Insert vector for this material into the table
952  G4PhysicsTableHelper::SetPhysicsVector(table, i, aVector);
953  }
954  }
955 
956  if(1 < verboseLevel) {
957  G4cout << "Lambda table is built for "
959  << G4endl;
960  }
961 
962  return table;
963 }
964 
965 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
966 
967 void G4VEnergyLossProcess::StreamInfo(std::ostream& out,
968  const G4ParticleDefinition& part, G4String endOfLine) const
969 {
970  out << std::setprecision(6);
971  out << endOfLine << GetProcessName() << ": ";
972  if (endOfLine != G4String("<br>\n")) {
973  out << " for " << part.GetParticleName();
974  }
975  out << " SubType= " << GetProcessSubType() << endOfLine
976  << " dE/dx and range tables from "
977  << G4BestUnit(minKinEnergy,"Energy")
978  << " to " << G4BestUnit(maxKinEnergy,"Energy")
979  << " in " << nBins << " bins" << endOfLine
980  << " Lambda tables from threshold to "
981  << G4BestUnit(maxKinEnergy,"Energy")
982  << ", " << theParameters->NumberOfBinsPerDecade()
983  << " bins per decade, spline: "
984  << theParameters->Spline()
985  << endOfLine;
987  out << " finalRange(mm)= " << finalRange/mm
988  << ", dRoverRange= " << dRoverRange
989  << ", integral: " << integral
990  << ", fluct: " << lossFluctuationFlag
991  << ", linLossLimit= " << linLossLimit
992  << endOfLine;
993  }
994  StreamProcessInfo(out, endOfLine);
995  modelManager->DumpModelList(out, verboseLevel, endOfLine);
997  out << " CSDA range table up"
998  << " to " << G4BestUnit(maxKinEnergyCSDA,"Energy")
999  << " in " << nBinsCSDA << " bins" << endOfLine;
1000  }
1001  if(nSCoffRegions>0 && isIonisation) {
1002  out << " Subcutoff sampling in " << nSCoffRegions
1003  << " regions" << endOfLine;
1004  }
1005  if(2 < verboseLevel) {
1006  out << " DEDXTable address= " << theDEDXTable << endOfLine;
1007  if(theDEDXTable && isIonisation) out << (*theDEDXTable) << endOfLine;
1008  out << "non restricted DEDXTable address= "
1009  << theDEDXunRestrictedTable << endOfLine;
1011  out << (*theDEDXunRestrictedTable) << endOfLine;
1012  }
1013  if(theDEDXSubTable && isIonisation) {
1014  out << (*theDEDXSubTable) << endOfLine;
1015  }
1016  out << " CSDARangeTable address= " << theCSDARangeTable
1017  << endOfLine;
1019  out << (*theCSDARangeTable) << endOfLine;
1020  }
1021  out << " RangeTableForLoss address= " << theRangeTableForLoss
1022  << endOfLine;
1024  out << (*theRangeTableForLoss) << endOfLine;
1025  }
1026  out << " InverseRangeTable address= " << theInverseRangeTable
1027  << endOfLine;
1029  out << (*theInverseRangeTable) << endOfLine;
1030  }
1031  out << " LambdaTable address= " << theLambdaTable << endOfLine;
1032  if(theLambdaTable && isIonisation) {
1033  out << (*theLambdaTable) << endOfLine;
1034  }
1035  out << " SubLambdaTable address= " << theSubLambdaTable
1036  << endOfLine;
1038  out << (*theSubLambdaTable) << endOfLine;
1039  }
1040  }
1041 }
1042 
1043 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1044 
1046 {
1047  G4RegionStore* regionStore = G4RegionStore::GetInstance();
1048  const G4Region* reg = r;
1049  if (!reg) {
1050  reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);
1051  }
1052 
1053  // the region is in the list
1054  if (nSCoffRegions > 0) {
1055  for (G4int i=0; i<nSCoffRegions; ++i) {
1056  if (reg == scoffRegions[i]) {
1057  return;
1058  }
1059  }
1060  }
1061  // new region
1062  if(val) {
1063  scoffRegions.push_back(reg);
1064  ++nSCoffRegions;
1065  }
1066 }
1067 
1068 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1069 
1071 {
1072  /*
1073  G4cout << track->GetDefinition()->GetParticleName()
1074  << " e(MeV)= " << track->GetKineticEnergy()
1075  << " baseParticle " << baseParticle << " proc " << this;
1076  if(particle) G4cout << " " << particle->GetParticleName();
1077  G4cout << " isIon= " << isIon << " dedx " << theDEDXTable <<G4endl;
1078  */
1079  // reset parameters for the new track
1082  preStepRangeEnergy = 0.0;
1083 
1084  // reset ion
1085  if(isIon) {
1086  chargeSqRatio = 0.5;
1087 
1088  G4double newmass = track->GetDefinition()->GetPDGMass();
1089  if(baseParticle) {
1090  massRatio = baseParticle->GetPDGMass()/newmass;
1091  } else if(theGenericIon) {
1092  massRatio = proton_mass_c2/newmass;
1093  } else {
1094  massRatio = 1.0;
1095  }
1096  }
1097  // forced biasing only for primary particles
1098  if(biasManager) {
1099  if(0 == track->GetParentID()) {
1100  // primary particle
1101  biasFlag = true;
1103  }
1104  }
1105 }
1106 
1107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1108 
1111  G4GPILSelection* selection)
1112 {
1113  G4double x = DBL_MAX;
1114  *selection = aGPILSelection;
1119  x = (fRange > finR) ?
1120  fRange*dRoverRange + finR*(1.0 - dRoverRange)*(2.0 - finR/fRange) : fRange;
1121  // if(particle->GetPDGMass() > 0.9*GeV)
1122  /*
1123  G4cout<<GetProcessName()<<": e= "<<preStepKinEnergy
1124  <<" range= "<<fRange << " idx= " << basedCoupleIndex
1125  << " finR= " << finR
1126  << " limit= " << x <<G4endl;
1127  G4cout << "massRatio= " << massRatio << " Q^2= " << chargeSqRatio
1128  << " finR= " << finR << " dRoverRange= " << dRoverRange
1129  << " finalRange= " << finalRange << G4endl;
1130  */
1131  }
1132  //G4cout<<GetProcessName()<<": e= "<<preStepKinEnergy
1133  //<<" stepLimit= "<<x<<G4endl;
1134  return x;
1135 }
1136 
1137 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1138 
1140  const G4Track& track,
1141  G4double previousStepSize,
1143 {
1144  // condition is set to "Not Forced"
1145  *condition = NotForced;
1146  G4double x = DBL_MAX;
1147 
1148  // initialisation of material, mass, charge, model
1149  // at the beginning of the step
1152  preStepScaledEnergy = preStepKinEnergy*massRatio;
1154 
1158  return x;
1159  }
1160 
1161  // change effective charge of an ion on fly
1162  if(isIon) {
1164  if(q2 != chargeSqRatio && q2 > 0.0) {
1165  chargeSqRatio = q2;
1166  fFactor = q2*biasFactor*(*theDensityFactor)[currentCoupleIndex];
1167  reduceFactor = 1.0/(fFactor*massRatio);
1168  }
1169  }
1170  // if(particle->GetPDGMass() > 0.9*GeV)
1171  //G4cout << "q2= "<<chargeSqRatio << " massRatio= " << massRatio << G4endl;
1172 
1173  // forced biasing only for primary particles
1174  if(biasManager) {
1175  if(0 == track.GetParentID() && biasFlag &&
1177  return biasManager->GetStepLimit(currentCoupleIndex, previousStepSize);
1178  }
1179  }
1180 
1181  // compute mean free path
1185 
1186  // zero cross section
1187  if(preStepLambda <= 0.0) {
1190  }
1191  }
1192 
1193  // non-zero cross section
1194  if(preStepLambda > 0.0) {
1196 
1197  // beggining of tracking (or just after DoIt of this process)
1200 
1201  } else if(currentInteractionLength < DBL_MAX) {
1202 
1203  // subtract NumberOfInteractionLengthLeft using previous step
1205  previousStepSize/currentInteractionLength;
1206 
1209  }
1210 
1211  // new mean free path and step limit
1214  }
1215 #ifdef G4VERBOSE
1216  if (verboseLevel>2){
1217  // if(particle->GetPDGMass() > 0.9*GeV){
1218  G4cout << "G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength ";
1219  G4cout << "[ " << GetProcessName() << "]" << G4endl;
1220  G4cout << " for " << track.GetDefinition()->GetParticleName()
1221  << " in Material " << currentMaterial->GetName()
1222  << " Ekin(MeV)= " << preStepKinEnergy/MeV
1223  << " " << track.GetMaterial()->GetName()
1224  <<G4endl;
1225  G4cout << "MeanFreePath = " << currentInteractionLength/cm << "[cm]"
1226  << "InteractionLength= " << x/cm <<"[cm] " <<G4endl;
1227  }
1228 #endif
1229  return x;
1230 }
1231 
1232 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1233 
1235  const G4Step& step)
1236 {
1238  // The process has range table - calculate energy loss
1240  return &fParticleChange;
1241  }
1242 
1243  // Get the actual (true) Step length
1244  G4double length = step.GetStepLength();
1245  if(length <= 0.0) { return &fParticleChange; }
1246  G4double eloss = 0.0;
1247 
1248  /*
1249  if(-1 < verboseLevel) {
1250  const G4ParticleDefinition* d = track.GetParticleDefinition();
1251  G4cout << "AlongStepDoIt for "
1252  << GetProcessName() << " and particle "
1253  << d->GetParticleName()
1254  << " eScaled(MeV)= " << preStepScaledEnergy/MeV
1255  << " range(mm)= " << fRange/mm
1256  << " s(mm)= " << length/mm
1257  << " rf= " << reduceFactor
1258  << " q^2= " << chargeSqRatio
1259  << " md= " << d->GetPDGMass()
1260  << " status= " << track.GetTrackStatus()
1261  << " " << track.GetMaterial()->GetName()
1262  << G4endl;
1263  }
1264  */
1265 
1266  const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
1267 
1268  // define new weight for primary and secondaries
1270  if(weightFlag) {
1271  weight /= biasFactor;
1273  }
1274 
1275  // stopping
1276  if (length >= fRange || preStepKinEnergy <= lowestKinEnergy) {
1277  eloss = preStepKinEnergy;
1278  if (useDeexcitation) {
1280  eloss, currentCoupleIndex);
1281  if(scTracks.size() > 0) { FillSecondariesAlongStep(eloss, weight); }
1282  eloss = std::max(eloss, 0.0);
1283  }
1286  return &fParticleChange;
1287  }
1288  //G4cout << theDEDXTable << " idx= " << basedCoupleIndex
1289  // << " " << GetProcessName() << " "<< currentMaterial->GetName()<<G4endl;
1290  //if(particle->GetParticleName() == "e-")G4cout << (*theDEDXTable) <<G4endl;
1291  // Short step
1293 
1294  //G4cout << "eloss= " << eloss << G4endl;
1295 
1296  // Long step
1297  if(eloss > preStepKinEnergy*linLossLimit) {
1298 
1299  G4double x = (fRange - length)/reduceFactor;
1300  //G4cout << "x= " << x << " " << theInverseRangeTable << G4endl;
1302 
1303  /*
1304  if(-1 < verboseLevel)
1305  G4cout << "Long STEP: rPre(mm)= "
1306  << GetScaledRangeForScaledEnergy(preStepScaledEnergy)/mm
1307  << " rPost(mm)= " << x/mm
1308  << " ePre(MeV)= " << preStepScaledEnergy/MeV
1309  << " eloss(MeV)= " << eloss/MeV
1310  << " eloss0(MeV)= "
1311  << GetDEDXForScaledEnergy(preStepScaledEnergy)*length/MeV
1312  << " lim(MeV)= " << preStepKinEnergy*linLossLimit/MeV
1313  << G4endl;
1314  */
1315  }
1316 
1317  /*
1318  G4double eloss0 = eloss;
1319  if(-1 < verboseLevel ) {
1320  G4cout << "Before fluct: eloss(MeV)= " << eloss/MeV
1321  << " e-eloss= " << preStepKinEnergy-eloss
1322  << " step(mm)= " << length/mm
1323  << " range(mm)= " << fRange/mm
1324  << " fluct= " << lossFluctuationFlag
1325  << G4endl;
1326  }
1327  */
1328 
1329  G4double cut = (*theCuts)[currentCoupleIndex];
1330  G4double esec = 0.0;
1331 
1332  //G4cout << "cut= " << cut << " useSubCut= " << useSubCutoff << G4endl;
1333 
1334  // SubCutOff
1335  if(useSubCutoff && !subcutProducer) {
1337 
1338  G4bool yes = false;
1339  const G4StepPoint* prePoint = step.GetPreStepPoint();
1340 
1341  // Check boundary
1342  if(prePoint->GetStepStatus() == fGeomBoundary) { yes = true; }
1343 
1344  // Check PrePoint
1345  else {
1346  G4double preSafety = prePoint->GetSafety();
1347  G4double rcut =
1349 
1350  // recompute presafety
1351  if(preSafety < rcut) {
1352  preSafety = safetyHelper->ComputeSafety(prePoint->GetPosition(),
1353  rcut);
1354  }
1355 
1356  if(preSafety < rcut) { yes = true; }
1357 
1358  // Check PostPoint
1359  else {
1360  G4double postSafety = preSafety - length;
1361  if(postSafety < rcut) {
1362  postSafety = safetyHelper->ComputeSafety(
1363  step.GetPostStepPoint()->GetPosition(), rcut);
1364  if(postSafety < rcut) { yes = true; }
1365  }
1366  }
1367  }
1368 
1369  // Decided to start subcut sampling
1370  if(yes) {
1371 
1372  cut = (*theSubCuts)[currentCoupleIndex];
1374  esec = SampleSubCutSecondaries(scTracks, step,
1375  currentModel,currentCoupleIndex);
1376  // add bremsstrahlung sampling
1377  /*
1378  if(nProcesses > 0) {
1379  for(G4int i=0; i<nProcesses; ++i) {
1380  (scProcesses[i])->SampleSubCutSecondaries(
1381  scTracks, step, (scProcesses[i])->
1382  SelectModelForMaterial(preStepKinEnergy, currentCoupleIndex),
1383  currentCoupleIndex);
1384  }
1385  }
1386  */
1387  }
1388  }
1389  }
1390 
1391  // Corrections, which cannot be tabulated
1392  if(isIon) {
1393  G4double eadd = 0.0;
1394  G4double eloss_before = eloss;
1396  eloss, eadd, length);
1397  if(eloss < 0.0) { eloss = 0.5*eloss_before; }
1398  }
1399 
1400  // Sample fluctuations
1401  if (lossFluctuationFlag) {
1403  if(eloss + esec < preStepKinEnergy) {
1404 
1405  G4double tmax =
1406  std::min(currentModel->MaxSecondaryKinEnergy(dynParticle),cut);
1407  eloss = fluc->SampleFluctuations(currentCouple,dynParticle,
1408  tmax,length,eloss);
1409  /*
1410  if(-1 < verboseLevel)
1411  G4cout << "After fluct: eloss(MeV)= " << eloss/MeV
1412  << " fluc= " << (eloss-eloss0)/MeV
1413  << " ChargeSqRatio= " << chargeSqRatio
1414  << " massRatio= " << massRatio
1415  << " tmax= " << tmax
1416  << G4endl;
1417  */
1418  }
1419  }
1420 
1421  // deexcitation
1422  if (useDeexcitation) {
1423  G4double esecfluo = preStepKinEnergy - esec;
1424  G4double de = esecfluo;
1425  //G4double eloss0 = eloss;
1426  /*
1427  G4cout << "### 1: E(keV)= " << preStepKinEnergy/keV
1428  << " Efluomax(keV)= " << de/keV
1429  << " Eloss(keV)= " << eloss/keV << G4endl;
1430  */
1432  de, currentCoupleIndex);
1433 
1434  // sum of de-excitation energies
1435  esecfluo -= de;
1436 
1437  // subtracted from energy loss
1438  if(eloss >= esecfluo) {
1439  esec += esecfluo;
1440  eloss -= esecfluo;
1441  } else {
1442  esec += esecfluo;
1443  eloss = 0.0;
1444  }
1445  /*
1446  if(esecfluo > 0.0) {
1447  G4cout << "### 2: E(keV)= " << preStepKinEnergy/keV
1448  << " Esec(keV)= " << esec/keV
1449  << " Esecf(kV)= " << esecfluo/keV
1450  << " Eloss0(kV)= " << eloss0/keV
1451  << " Eloss(keV)= " << eloss/keV
1452  << G4endl;
1453  }
1454  */
1455  }
1457  subcutProducer->SampleSecondaries(step, scTracks, eloss, cut);
1458  }
1459  if(scTracks.size() > 0) { FillSecondariesAlongStep(eloss, weight); }
1460 
1461  // Energy balance
1462  G4double finalT = preStepKinEnergy - eloss - esec;
1463  if (finalT <= lowestKinEnergy) {
1464  eloss += finalT;
1465  finalT = 0.0;
1466  } else if(isIon) {
1469  currentMaterial,finalT));
1470  }
1471 
1472  eloss = std::max(eloss, 0.0);
1473 
1476  /*
1477  if(-1 < verboseLevel) {
1478  G4double del = finalT + eloss + esec - preStepKinEnergy;
1479  G4cout << "Final value eloss(MeV)= " << eloss/MeV
1480  << " preStepKinEnergy= " << preStepKinEnergy
1481  << " postStepKinEnergy= " << finalT
1482  << " de(keV)= " << del/keV
1483  << " lossFlag= " << lossFluctuationFlag
1484  << " status= " << track.GetTrackStatus()
1485  << G4endl;
1486  }
1487  */
1488  return &fParticleChange;
1489 }
1490 
1491 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1492 
1493 void
1495 {
1496  G4int n0 = scTracks.size();
1497 
1498  // weight may be changed by biasing manager
1499  if(biasManager) {
1501  weight *=
1503  }
1504  }
1505 
1506  // fill secondaries
1507  G4int n = scTracks.size();
1509 
1510  for(G4int i=0; i<n; ++i) {
1511  G4Track* t = scTracks[i];
1512  if(t) {
1513  t->SetWeight(weight);
1515  if(i >= n0) { t->SetCreatorModelIndex(biasID); }
1516  //G4cout << "Secondary(along step) has weight " << t->GetWeight()
1517  //<< ", kenergy " << t->GetKineticEnergy()/MeV << " MeV" <<G4endl;
1518  }
1519  }
1520  scTracks.clear();
1521 }
1522 
1523 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1524 
1525 G4double
1526 G4VEnergyLossProcess::SampleSubCutSecondaries(std::vector<G4Track*>& tracks,
1527  const G4Step& step,
1528  G4VEmModel* model,
1529  G4int idx)
1530 {
1531  // Fast check weather subcutoff can work
1532  G4double esec = 0.0;
1533  G4double subcut = (*theSubCuts)[idx];
1534  G4double cut = (*theCuts)[idx];
1535  if(cut <= subcut) { return esec; }
1536 
1537  const G4Track* track = step.GetTrack();
1538  const G4DynamicParticle* dp = track->GetDynamicParticle();
1540  G4double cross = (*theDensityFactor)[idx]*chargeSqRatio
1541  *(((*theSubLambdaTable)[(*theDensityIdx)[idx]])->Value(e, idxSubLambda));
1542  G4double length = step.GetStepLength();
1543 
1544  // negligible probability to get any interaction
1545  if(length*cross < perMillion) { return esec; }
1546  /*
1547  if(-1 < verboseLevel)
1548  G4cout << "<<< Subcutoff for " << GetProcessName()
1549  << " cross(1/mm)= " << cross*mm << ">>>"
1550  << " e(MeV)= " << preStepScaledEnergy
1551  << " matIdx= " << currentCoupleIndex
1552  << G4endl;
1553  */
1554 
1555  // Sample subcutoff secondaries
1556  G4StepPoint* preStepPoint = step.GetPreStepPoint();
1557  G4StepPoint* postStepPoint = step.GetPostStepPoint();
1558  G4ThreeVector prepoint = preStepPoint->GetPosition();
1559  G4ThreeVector dr = postStepPoint->GetPosition() - prepoint;
1560  G4double pretime = preStepPoint->GetGlobalTime();
1561  G4double dt = postStepPoint->GetGlobalTime() - pretime;
1562  G4double fragment = 0.0;
1563 
1564  do {
1565  G4double del = -G4Log(G4UniformRand())/cross;
1566  fragment += del/length;
1567  if (fragment > 1.0) { break; }
1568 
1569  // sample secondaries
1570  secParticles.clear();
1572  dp,subcut,cut);
1573 
1574  // position of subcutoff particles
1575  G4ThreeVector r = prepoint + fragment*dr;
1576  std::vector<G4DynamicParticle*>::iterator it;
1577  for(it=secParticles.begin(); it!=secParticles.end(); ++it) {
1578 
1579  G4Track* t = new G4Track((*it), pretime + fragment*dt, r);
1582  tracks.push_back(t);
1583  esec += t->GetKineticEnergy();
1584  if (t->GetParticleDefinition() == thePositron) {
1585  esec += 2.0*electron_mass_c2;
1586  }
1587 
1588  /*
1589  if(-1 < verboseLevel)
1590  G4cout << "New track "
1591  << t->GetParticleDefinition()->GetParticleName()
1592  << " e(keV)= " << t->GetKineticEnergy()/keV
1593  << " fragment= " << fragment
1594  << G4endl;
1595  */
1596  }
1597  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
1598  } while (fragment <= 1.0);
1599  return esec;
1600 }
1601 
1602 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1603 
1605  const G4Step& step)
1606 {
1607  // In all cases clear number of interaction lengths
1610 
1612  G4double finalT = track.GetKineticEnergy();
1613  if(finalT <= lowestKinEnergy) { return &fParticleChange; }
1614 
1615  G4double postStepScaledEnergy = finalT*massRatio;
1616  SelectModel(postStepScaledEnergy);
1617 
1618  if(!currentModel->IsActive(postStepScaledEnergy)) {
1619  return &fParticleChange;
1620  }
1621  /*
1622  if(-1 < verboseLevel) {
1623  G4cout << GetProcessName()
1624  << "::PostStepDoIt: E(MeV)= " << finalT/MeV
1625  << G4endl;
1626  }
1627  */
1628 
1629  // forced process - should happen only once per track
1630  if(biasFlag) {
1632  biasFlag = false;
1633  }
1634  }
1635 
1636  // Integral approach
1637  if (integral) {
1638  G4double lx = GetLambdaForScaledEnergy(postStepScaledEnergy);
1639  /*
1640  if(preStepLambda<lx && 1 < verboseLevel) {
1641  G4cout << "WARNING: for " << particle->GetParticleName()
1642  << " and " << GetProcessName()
1643  << " E(MeV)= " << finalT/MeV
1644  << " preLambda= " << preStepLambda
1645  << " < " << lx << " (postLambda) "
1646  << G4endl;
1647  }
1648  */
1649  if(lx <= 0.0 || preStepLambda*G4UniformRand() > lx) {
1650  return &fParticleChange;
1651  }
1652  }
1653 
1654  SelectModel(postStepScaledEnergy);
1655 
1656  // define new weight for primary and secondaries
1658  if(weightFlag) {
1659  weight /= biasFactor;
1661  }
1662 
1663  const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
1664  G4double tcut = (*theCuts)[currentCoupleIndex];
1665 
1666  // sample secondaries
1667  secParticles.clear();
1668  //G4cout<< "@@@ Eprimary= "<<dynParticle->GetKineticEnergy()/MeV
1669  // << " cut= " << tcut/MeV << G4endl;
1671  dynParticle, tcut);
1672 
1673  G4int num0 = secParticles.size();
1674 
1675  // bremsstrahlung splitting or Russian roulette
1676  if(biasManager) {
1678  G4double eloss = 0.0;
1680  secParticles,
1681  track, currentModel,
1682  &fParticleChange, eloss,
1683  currentCoupleIndex, tcut,
1684  step.GetPostStepPoint()->GetSafety());
1685  if(eloss > 0.0) {
1688  }
1689  }
1690  }
1691 
1692  // save secondaries
1693  G4int num = secParticles.size();
1694  if(num > 0) {
1695 
1697  G4double time = track.GetGlobalTime();
1698 
1699  for (G4int i=0; i<num; ++i) {
1700  if(secParticles[i]) {
1701  G4Track* t = new G4Track(secParticles[i], time, track.GetPosition());
1703  t->SetWeight(weight);
1704  if(i < num0) { t->SetCreatorModelIndex(secID); }
1705  else { t->SetCreatorModelIndex(biasID); }
1706 
1707  //G4cout << "Secondary(post step) has weight " << t->GetWeight()
1708  // << ", kenergy " << t->GetKineticEnergy()/MeV << " MeV"
1709  // << " time= " << time/ns << " ns " << G4endl;
1711  }
1712  }
1713  }
1714 
1720  }
1721 
1722  /*
1723  if(-1 < verboseLevel) {
1724  G4cout << "::PostStepDoIt: Sample secondary; Efin= "
1725  << fParticleChange.GetProposedKineticEnergy()/MeV
1726  << " MeV; model= (" << currentModel->LowEnergyLimit()
1727  << ", " << currentModel->HighEnergyLimit() << ")"
1728  << " preStepLambda= " << preStepLambda
1729  << " dir= " << track.GetMomentumDirection()
1730  << " status= " << track.GetTrackStatus()
1731  << G4endl;
1732  }
1733  */
1734  return &fParticleChange;
1735 }
1736 
1737 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1738 
1740  const G4ParticleDefinition* part, const G4String& directory,
1741  G4bool ascii)
1742 {
1743  G4bool res = true;
1744  //G4cout << "G4VEnergyLossProcess::StorePhysicsTable: " << part->GetParticleName()
1745  // << " " << directory << " " << ascii << G4endl;
1746  if (!isMaster || baseParticle || part != particle ) return res;
1747 
1748  if(!StoreTable(part,theDEDXTable,ascii,directory,"DEDX"))
1749  {res = false;}
1750 
1751  if(!StoreTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr"))
1752  {res = false;}
1753 
1754  if(!StoreTable(part,theDEDXSubTable,ascii,directory,"SubDEDX"))
1755  {res = false;}
1756 
1757  if(!StoreTable(part,theIonisationTable,ascii,directory,"Ionisation"))
1758  {res = false;}
1759 
1760  if(!StoreTable(part,theIonisationSubTable,ascii,directory,"SubIonisation"))
1761  {res = false;}
1762 
1763  if(isIonisation &&
1764  !StoreTable(part,theCSDARangeTable,ascii,directory,"CSDARange"))
1765  {res = false;}
1766 
1767  if(isIonisation &&
1768  !StoreTable(part,theRangeTableForLoss,ascii,directory,"Range"))
1769  {res = false;}
1770 
1771  if(isIonisation &&
1772  !StoreTable(part,theInverseRangeTable,ascii,directory,"InverseRange"))
1773  {res = false;}
1774 
1775  if(!StoreTable(part,theLambdaTable,ascii,directory,"Lambda"))
1776  {res = false;}
1777 
1778  if(!StoreTable(part,theSubLambdaTable,ascii,directory,"SubLambda"))
1779  {res = false;}
1780 
1781  if ( !res ) {
1782  if(1 < verboseLevel) {
1783  G4cout << "Physics tables are stored for "
1785  << " and process " << GetProcessName()
1786  << " in the directory <" << directory
1787  << "> " << G4endl;
1788  }
1789  } else {
1790  G4cout << "Fail to store Physics Tables for "
1792  << " and process " << GetProcessName()
1793  << " in the directory <" << directory
1794  << "> " << G4endl;
1795  }
1796  return res;
1797 }
1798 
1799 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
1800 
1801 G4bool
1803  const G4String& directory,
1804  G4bool ascii)
1805 {
1806  G4bool res = true;
1807  if (!isMaster) return res;
1808  const G4String particleName = part->GetParticleName();
1809 
1810  if(1 < verboseLevel) {
1811  G4cout << "G4VEnergyLossProcess::RetrievePhysicsTable() for "
1812  << particleName << " and process " << GetProcessName()
1813  << "; tables_are_built= " << tablesAreBuilt
1814  << G4endl;
1815  }
1816  if(particle == part) {
1817 
1818  if ( !baseParticle ) {
1819 
1820  G4bool fpi = true;
1821  if(!RetrieveTable(part,theDEDXTable,ascii,directory,"DEDX",fpi))
1822  {fpi = false;}
1823 
1824  // ionisation table keeps individual dEdx and not sum of sub-processes
1825  if(!RetrieveTable(part,theDEDXTable,ascii,directory,"Ionisation",false))
1826  {fpi = false;}
1827 
1828  if(!RetrieveTable(part,theRangeTableForLoss,ascii,directory,"Range",fpi))
1829  {res = false;}
1830 
1831  if(!RetrieveTable(part,theDEDXunRestrictedTable,ascii,directory,
1832  "DEDXnr",false))
1833  {res = false;}
1834 
1835  if(!RetrieveTable(part,theCSDARangeTable,ascii,directory,
1836  "CSDARange",false))
1837  {res = false;}
1838 
1839  if(!RetrieveTable(part,theInverseRangeTable,ascii,directory,
1840  "InverseRange",fpi))
1841  {res = false;}
1842 
1843  if(!RetrieveTable(part,theLambdaTable,ascii,directory,"Lambda",true))
1844  {res = false;}
1845 
1846  G4bool yes = false;
1847  if(nSCoffRegions > 0) {yes = true;}
1848 
1849  if(!RetrieveTable(part,theDEDXSubTable,ascii,directory,"SubDEDX",yes))
1850  {res = false;}
1851 
1852  if(!RetrieveTable(part,theSubLambdaTable,ascii,directory,
1853  "SubLambda",yes))
1854  {res = false;}
1855 
1856  if(!fpi) yes = false;
1857  if(!RetrieveTable(part,theIonisationSubTable,ascii,directory,
1858  "SubIonisation",yes))
1859  {res = false;}
1860  }
1861  }
1862 
1863  return res;
1864 }
1865 
1866 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
1867 
1869  G4PhysicsTable* aTable, G4bool ascii,
1870  const G4String& directory,
1871  const G4String& tname)
1872 {
1873  //G4cout << "G4VEnergyLossProcess::StoreTable: " << aTable
1874  // << " " << directory << " " << tname << G4endl;
1875  G4bool res = true;
1876  if ( aTable ) {
1877  const G4String name = GetPhysicsTableFileName(part,directory,tname,ascii);
1878  G4cout << name << G4endl;
1879  //G4cout << *aTable << G4endl;
1880  if( !aTable->StorePhysicsTable(name,ascii)) res = false;
1881  }
1882  return res;
1883 }
1884 
1885 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
1886 
1887 G4bool
1889  G4PhysicsTable* aTable,
1890  G4bool ascii,
1891  const G4String& directory,
1892  const G4String& tname,
1893  G4bool mandatory)
1894 {
1895  G4bool isRetrieved = false;
1896  G4String filename = GetPhysicsTableFileName(part,directory,tname,ascii);
1897  if(aTable) {
1898  if(aTable->ExistPhysicsTable(filename)) {
1899  if(G4PhysicsTableHelper::RetrievePhysicsTable(aTable,filename,ascii)) {
1900  isRetrieved = true;
1901  if(theParameters->Spline()) {
1902  size_t n = aTable->length();
1903  for(size_t i=0; i<n; ++i) {
1904  if((*aTable)[i]) { (*aTable)[i]->SetSpline(true); }
1905  }
1906  }
1907  if (0 < verboseLevel) {
1908  G4cout << tname << " table for " << part->GetParticleName()
1909  << " is Retrieved from <" << filename << ">"
1910  << G4endl;
1911  }
1912  }
1913  }
1914  }
1915  if(mandatory && !isRetrieved) {
1916  if(0 < verboseLevel) {
1917  G4cout << tname << " table for " << part->GetParticleName()
1918  << " from file <"
1919  << filename << "> is not Retrieved"
1920  << G4endl;
1921  }
1922  return false;
1923  }
1924  return true;
1925 }
1926 
1927 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1928 
1930  const G4MaterialCutsCouple *couple,
1931  const G4DynamicParticle* dp,
1932  G4double length)
1933 {
1934  DefineMaterial(couple);
1935  G4double ekin = dp->GetKineticEnergy();
1936  SelectModel(ekin*massRatio);
1938  tmax = std::min(tmax,(*theCuts)[currentCoupleIndex]);
1939  G4double d = 0.0;
1941  if(fm) { d = fm->Dispersion(currentMaterial,dp,tmax,length); }
1942  return d;
1943 }
1944 
1945 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1946 
1949 {
1950  // Cross section per volume is calculated
1951  DefineMaterial(couple);
1952  G4double cross = 0.0;
1953  if(theLambdaTable) {
1954  cross = GetLambdaForScaledEnergy(kineticEnergy*massRatio);
1955  } else {
1956  SelectModel(kineticEnergy*massRatio);
1957  cross = biasFactor*(*theDensityFactor)[currentCoupleIndex]
1959  particle, kineticEnergy,
1961  }
1962  if(cross < 0.0) { cross = 0.0; }
1963  return cross;
1964 }
1965 
1966 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1967 
1969 {
1972  G4double x = DBL_MAX;
1973  if(0.0 < preStepLambda) { x = 1.0/preStepLambda; }
1974  return x;
1975 }
1976 
1977 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1978 
1980  G4double x, G4double y,
1981  G4double& z)
1982 {
1983  G4GPILSelection sel;
1984  return AlongStepGetPhysicalInteractionLength(track, x, y, z, &sel);
1985 }
1986 
1987 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1988 
1990  const G4Track& track,
1991  G4double,
1993 
1994 {
1995  *condition = NotForced;
1996  return MeanFreePath(track);
1997 }
1998 
1999 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2000 
2002  const G4Track&,
2004 {
2005  return DBL_MAX;
2006 }
2007 
2008 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2009 
2012  G4double)
2013 {
2014  G4PhysicsVector* v =
2017  return v;
2018 }
2019 
2020 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2021 
2024 {
2025  G4bool add = true;
2026  if(p->GetProcessName() != "eBrem") { add = false; }
2027  if(add && nProcesses > 0) {
2028  for(G4int i=0; i<nProcesses; ++i) {
2029  if(p == scProcesses[i]) {
2030  add = false;
2031  break;
2032  }
2033  }
2034  }
2035  if(add) {
2036  scProcesses.push_back(p);
2037  ++nProcesses;
2038  if (1 < verboseLevel) {
2039  G4cout << "### The process " << p->GetProcessName()
2040  << " is added to the list of collaborative processes of "
2041  << GetProcessName() << G4endl;
2042  }
2043  }
2044 }
2045 
2046 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2047 
2048 void
2050 {
2051  if(fTotal == tType) {
2053  if(p) {
2054  size_t n = p->length();
2055  G4PhysicsVector* pv = (*p)[0];
2057 
2061 
2062  for (size_t i=0; i<n; ++i) {
2063  G4double dedx = 0.0;
2064  pv = (*p)[i];
2065  if(pv) {
2066  dedx = pv->Value(emax, idxDEDXunRestricted);
2067  } else {
2068  pv = (*p)[(*theDensityIdx)[i]];
2069  if(pv) {
2070  dedx =
2071  pv->Value(emax, idxDEDXunRestricted)*(*theDensityFactor)[i];
2072  }
2073  }
2074  theDEDXAtMaxEnergy[i] = dedx;
2075  //G4cout << "i= " << i << " emax(MeV)= " << emax/MeV<< " dedx= "
2076  // << dedx << G4endl;
2077  }
2078  }
2079 
2080  } else if(fRestricted == tType) {
2081  /*
2082  G4cout<< "G4VEnergyLossProcess::SetDEDXTable "
2083  << particle->GetParticleName()
2084  << " oldTable " << theDEDXTable << " newTable " << p
2085  << " ion " << theIonisationTable
2086  << " IsMaster " << isMaster
2087  << " " << GetProcessName() << G4endl;
2088  G4cout << (*p) << G4endl;
2089  */
2090  theDEDXTable = p;
2091  } else if(fSubRestricted == tType) {
2092  theDEDXSubTable = p;
2093  } else if(fIsIonisation == tType) {
2094  /*
2095  G4cout<< "G4VEnergyLossProcess::SetIonisationTable "
2096  << particle->GetParticleName()
2097  << " oldTable " << theDEDXTable << " newTable " << p
2098  << " ion " << theIonisationTable
2099  << " IsMaster " << isMaster
2100  << " " << GetProcessName() << G4endl;
2101  */
2103  } else if(fIsSubIonisation == tType) {
2105  }
2106 }
2107 
2108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2109 
2111 {
2112  theCSDARangeTable = p;
2113 
2114  if(p) {
2115  size_t n = p->length();
2116  G4PhysicsVector* pv;
2118 
2119  for (size_t i=0; i<n; ++i) {
2120  pv = (*p)[i];
2121  G4double rmax = 0.0;
2122  if(pv) { rmax = pv->Value(emax, idxCSDA); }
2123  else {
2124  pv = (*p)[(*theDensityIdx)[i]];
2125  if(pv) { rmax = pv->Value(emax, idxCSDA)/(*theDensityFactor)[i]; }
2126  }
2127  theRangeAtMaxEnergy[i] = rmax;
2128  //G4cout << "i= " << i << " Emax(MeV)= " << emax/MeV << " Rmax= "
2129  //<< rmax<< G4endl;
2130  }
2131  }
2132 }
2133 
2134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2135 
2137 {
2139  if(1 < verboseLevel) {
2140  G4cout << "### Set Range table " << p
2141  << " for " << particle->GetParticleName()
2142  << " and process " << GetProcessName() << G4endl;
2143  }
2144 }
2145 
2146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2147 
2149 {
2151  if(1 < verboseLevel) {
2152  G4cout << "### Set SecondaryRange table " << p
2153  << " for " << particle->GetParticleName()
2154  << " and process " << GetProcessName() << G4endl;
2155  }
2156 }
2157 
2158 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2159 
2161 {
2163  if(1 < verboseLevel) {
2164  G4cout << "### Set InverseRange table " << p
2165  << " for " << particle->GetParticleName()
2166  << " and process " << GetProcessName() << G4endl;
2167  }
2168 }
2169 
2170 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2171 
2173 {
2174  if(1 < verboseLevel) {
2175  G4cout << "### Set Lambda table " << p
2176  << " for " << particle->GetParticleName()
2177  << " and process " << GetProcessName() << G4endl;
2178  //G4cout << *p << G4endl;
2179  }
2180  theLambdaTable = p;
2181  tablesAreBuilt = true;
2182 
2186 
2187  if(theLambdaTable) {
2188  size_t n = theLambdaTable->length();
2189  G4PhysicsVector* pv = (*theLambdaTable)[0];
2190  G4double e, ss, smax, emax;
2191 
2192  size_t i;
2193 
2194  // first loop on existing vectors
2195  for (i=0; i<n; ++i) {
2196  pv = (*theLambdaTable)[i];
2197  if(pv) {
2198  size_t nb = pv->GetVectorLength();
2199  emax = DBL_MAX;
2200  smax = 0.0;
2201  if(nb > 0) {
2202  for (size_t j=0; j<nb; ++j) {
2203  e = pv->Energy(j);
2204  ss = (*pv)(j);
2205  if(ss > smax) {
2206  smax = ss;
2207  emax = e;
2208  }
2209  }
2210  }
2212  theCrossSectionMax[i] = smax;
2213  if(1 < verboseLevel) {
2214  G4cout << "For " << particle->GetParticleName()
2215  << " Max CS at i= " << i << " emax(MeV)= " << emax/MeV
2216  << " lambda= " << smax << G4endl;
2217  }
2218  }
2219  }
2220  // second loop using base materials
2221  for (i=0; i<n; ++i) {
2222  pv = (*theLambdaTable)[i];
2223  if(!pv){
2224  G4int j = (*theDensityIdx)[i];
2226  theCrossSectionMax[i] = (*theDensityFactor)[i]*theCrossSectionMax[j];
2227  }
2228  }
2229  }
2230 }
2231 
2232 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2233 
2235 {
2236  theSubLambdaTable = p;
2237  if(1 < verboseLevel) {
2238  G4cout << "### Set SebLambda table " << p
2239  << " for " << particle->GetParticleName()
2240  << " and process " << GetProcessName() << G4endl;
2241  }
2242 }
2243 
2244 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2245 
2247 {
2248  const G4Element* elm = nullptr;
2249  if(currentModel) { elm = currentModel->GetCurrentElement(); }
2250  return elm;
2251 }
2252 
2253 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2254 
2256  G4bool flag)
2257 {
2258  if(f > 0.0) {
2259  biasFactor = f;
2260  weightFlag = flag;
2261  if(1 < verboseLevel) {
2262  G4cout << "### SetCrossSectionBiasingFactor: for "
2263  << " process " << GetProcessName()
2264  << " biasFactor= " << f << " weightFlag= " << flag
2265  << G4endl;
2266  }
2267  }
2268 }
2269 
2270 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2271 
2272 void
2274  const G4String& region,
2275  G4bool flag)
2276 {
2277  if(!biasManager) { biasManager = new G4EmBiasingManager(); }
2278  if(1 < verboseLevel) {
2279  G4cout << "### ActivateForcedInteraction: for "
2280  << " process " << GetProcessName()
2281  << " length(mm)= " << length/mm
2282  << " in G4Region <" << region
2283  << "> weightFlag= " << flag
2284  << G4endl;
2285  }
2286  weightFlag = flag;
2287  biasManager->ActivateForcedInteraction(length, region);
2288 }
2289 
2290 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2291 
2292 void
2294  G4double factor,
2295  G4double energyLimit)
2296 {
2297  if (0.0 <= factor) {
2298 
2299  // Range cut can be applied only for e-
2300  if(0.0 == factor && secondaryParticle != G4Electron::Electron())
2301  { return; }
2302 
2303  if(!biasManager) { biasManager = new G4EmBiasingManager(); }
2304  biasManager->ActivateSecondaryBiasing(region, factor, energyLimit);
2305  if(1 < verboseLevel) {
2306  G4cout << "### ActivateSecondaryBiasing: for "
2307  << " process " << GetProcessName()
2308  << " factor= " << factor
2309  << " in G4Region <" << region
2310  << "> energyLimit(MeV)= " << energyLimit/MeV
2311  << G4endl;
2312  }
2313  }
2314 }
2315 
2316 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2317 
2319 {
2320  isIonisation = val;
2321  if(val) { aGPILSelection = CandidateForSelection; }
2323 }
2324 
2325 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2326 
2328 {
2329  if(0.0 < val && val < 1.0) {
2330  linLossLimit = val;
2331  actLinLossLimit = true;
2332  } else { PrintWarning("SetLinearLossLimit", val); }
2333 }
2334 
2335 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2336 
2337 void
2339 {
2340  if(actStepFunc) { return; }
2341  actStepFunc = lock;
2342  if(0.0 < v1 && 0.0 < v2 && v2 < 1.e+50) {
2343  dRoverRange = std::min(1.0, v1);
2344  finalRange = v2;
2345  } else if(v1 <= 0.0) {
2346  PrintWarning("SetStepFunction", v1);
2347  } else {
2348  PrintWarning("SetStepFunction", v2);
2349  }
2350 }
2351 
2352 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2353 
2355 {
2356  if(1.e-18 < val && val < 1.e+50) { lowestKinEnergy = val; }
2357  else { PrintWarning("SetLowestEnergyLimit", val); }
2358 }
2359 
2360 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2361 
2363 {
2364  if(2 < n && n < 1000000000) {
2365  nBins = n;
2366  actBinning = true;
2367  } else {
2368  G4double e = (G4double)n;
2369  PrintWarning("SetDEDXBinning", e);
2370  }
2371 }
2372 
2373 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2374 
2376 {
2377  if(1.e-18 < e && e < maxKinEnergy) {
2378  minKinEnergy = e;
2379  actMinKinEnergy = true;
2380  } else { PrintWarning("SetMinKinEnergy", e); }
2381 }
2382 
2383 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2384 
2386 {
2387  if(minKinEnergy < e && e < 1.e+50) {
2388  maxKinEnergy = e;
2389  actMaxKinEnergy = true;
2390  if(e < maxKinEnergyCSDA) { maxKinEnergyCSDA = e; }
2391  } else { PrintWarning("SetMaxKinEnergy", e); }
2392 }
2393 
2394 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2395 
2397 {
2398  G4String ss = "G4VEnergyLossProcess::" + tit;
2400  ed << "Parameter is out of range: " << val
2401  << " it will have no effect!\n" << " Process "
2402  << GetProcessName() << " nbins= " << nBins
2403  << " Emin(keV)= " << minKinEnergy/keV
2404  << " Emax(GeV)= " << maxKinEnergy/GeV;
2405  G4Exception(ss, "em0044", JustWarning, ed);
2406 }
2407 
2408 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2409 
2410 void G4VEnergyLossProcess::ProcessDescription(std::ostream& out) const
2411 {
2412  if(particle) { StreamInfo(out, *particle, G4String("<br>\n")); }
2413 }
2414 
2415 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
virtual G4bool RetrievePhysicsTable(const G4ParticleDefinition *, const G4String &directory, G4bool ascii) override
Float_t x
Definition: compare.C:6
void SetDEDXTable(G4PhysicsTable *p, G4EmTableType tType)
const std::vector< G4double > * theDensityFactor
G4double Energy(size_t index) const
void SetAngularGeneratorFlag(G4bool)
Definition: G4VEmModel.hh:679
static void SetPhysicsVector(G4PhysicsTable *physTable, size_t idx, G4PhysicsVector *vec)
G4double GetKineticEnergy() const
const G4Element * GetCurrentElement() const
void AddSecondary(G4Track *aSecondary)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
const XML_Char * name
Definition: expat.h:151
G4PhysicsTable * theCSDARangeTable
G4VEmFluctuationModel * fluctModel
std::vector< G4double > theCrossSectionMax
G4bool IsActive(G4double kinEnergy)
Definition: G4VEmModel.hh:735
void FillDEDXVector(G4PhysicsVector *, const G4MaterialCutsCouple *, G4EmTableType t=fRestricted)
std::vector< G4DynamicParticle * > secParticles
void SetMinKinEnergy(G4double e)
G4PhysicsTable * DEDXTable() const
const G4ParticleDefinition * theElectron
void AddCollaborativeProcess(G4VEnergyLossProcess *)
void SetStepFunction(G4double v1, G4double v2, G4bool lock=true)
G4SafetyHelper * GetSafetyHelper() const
G4PhysicsTable * BuildDEDXTable(G4EmTableType tType=fRestricted)
G4int NumberOfBinsPerDecade() const
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:45
G4PhysicsTable * theLambdaTable
void SetEmModel(G4VEmModel *, G4int index=0)
std::vector< G4double > theRangeAtMaxEnergy
static G4PhysicsTable * PreparePhysicsTable(G4PhysicsTable *physTable)
static G4ParticleTable * GetParticleTable()
G4VAtomDeexcitation * AtomDeexcitation()
void ActivateSecondaryBiasing(const G4String &region, G4double factor, G4double energyLimit)
static constexpr double MeV
Definition: G4SIunits.hh:214
void SetWeight(G4double aValue)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:707
G4bool UseAngularGeneratorForIonisation() const
G4PhysicsTable * IonisationTableForSubsec() const
G4StepPoint * GetPreStepPoint() const
const G4DataVector * Initialise(const G4ParticleDefinition *part, const G4ParticleDefinition *secPart, G4double minSubRange, G4int verb)
static constexpr double keV
Definition: G4SIunits.hh:216
void AlongStepDeexcitation(std::vector< G4Track * > &tracks, const G4Step &step, G4double &eLoss, G4int coupleIndex)
static constexpr double mm
Definition: G4SIunits.hh:115
G4double MaxSecondaryKinEnergy(const G4DynamicParticle *dynParticle)
Definition: G4VEmModel.hh:492
void ProposeWeight(G4double finalWeight)
G4PhysicsTable * theIonisationTable
const G4String & GetName() const
void SetCreatorModelIndex(G4int idx)
const std::vector< G4int > * GetCoupleIndexes()
void InitializeForAlongStep(const G4Track &)
#define G4endl
Definition: G4ios.hh:61
G4PhysicsTable * theDEDXTable
G4PhysicsTable * theRangeTableForLoss
Float_t y
Definition: compare.C:6
const G4VProcess * GetMasterProcess() const
Definition: G4VProcess.hh:541
const char * p
Definition: xmltok.h:285
void DefineRegParamForLoss(G4VEnergyLossProcess *, G4bool isElectron) const
void BuildPhysicsTable(const G4ParticleDefinition *aParticle)
Double_t z
G4double SampleSubCutSecondaries(std::vector< G4Track * > &, const G4Step &, G4VEmModel *model, G4int matIdx)
G4double MaxKinEnergy() const
G4bool BuildCSDARange() const
const G4DataVector * SubCutoff() const
std::vector< G4double > theEnergyOfCrossSectionMax
G4PhysicsTable * RangeTableForLoss() const
void PreparePhysicsTable(const G4ParticleDefinition *aParticle, G4VEnergyLossProcess *p, G4bool theMaster)
const G4String & GetParticleName() const
void ActivateForcedInteraction(G4double length=0.0, const G4String &r="")
std::vector< G4double > theDEDXAtMaxEnergy
G4double GetScaledRangeForScaledEnergy(G4double scaledKinEnergy)
const G4String & GetParticleType() const
const G4ParticleDefinition * thePositron
G4bool IsMaster() const
const G4TouchableHandle & GetTouchableHandle() const
virtual void StartTracking(G4Track *) override
G4double GetPDGCharge() const
static constexpr double perMillion
Definition: G4SIunits.hh:334
const G4DataVector * theSubCuts
void SetSecondaryWeightByProcess(G4bool)
G4bool GetFlag(size_t idx) const
G4bool UseCutAsFinalRange() const
virtual void PreparePhysicsTable(const G4ParticleDefinition &) override
G4bool StorePhysicsTable(const G4String &filename, G4bool ascii=false)
G4double ApplySecondaryBiasing(std::vector< G4DynamicParticle * > &, const G4Track &track, G4VEmModel *currentModel, G4ParticleChangeForGamma *pParticleChange, G4double &eloss, G4int coupleIdx, G4double tcut, G4double safety=0.0)
void ActivateForcedInteraction(G4double length, const G4String &region, G4bool flag=true)
const G4ParticleDefinition * theGenericIon
void SetSubLambdaTable(G4PhysicsTable *p)
static const G4double emax
const G4ParticleDefinition * GetParticleDefinition() const
G4PhysicsTable * SubLambdaTable() const
G4double Value(G4double theEnergy, size_t &lastidx) const
G4double condition(const G4ErrorSymMatrix &m)
G4bool ForcedInteractionRegion(G4int coupleIdx)
const std::vector< G4int > * theDensityIdx
G4PhysicsTable * theInverseRangeTable
void SetMasterThread(G4bool val)
Definition: G4VEmModel.hh:693
void SetDEDXBinning(G4int nbins)
void SetProposedCharge(G4double theCharge)
G4bool SecondaryBiasingRegion(G4int coupleIdx)
G4ParticleDefinition * GetDefinition() const
std::vector< G4VEmModel * > emModels
G4EmTableType
G4VEmModel * GetModel(G4int idx, G4bool ver=false)
G4double GetPDGMass() const
void DumpModelList(std::ostream &out, G4int verb)
G4ProcessVector * GetAlongStepProcessVector(G4ProcessVectorTypeIndex typ=typeGPIL) const
G4PhysicsTable * theDEDXSubTable
G4EmParameters * theParameters
G4double G4Log(G4double x)
Definition: G4Log.hh:230
void Initialise(const G4ParticleDefinition &part, const G4String &procName, G4int verbose)
G4Region * GetRegion(const G4String &name, G4bool verbose=true) const
G4double GetProductionCut(G4int index) const
virtual G4double Dispersion(const G4Material *, const G4DynamicParticle *, G4double tmax, G4double length)=0
std::vector< G4Track * > scTracks
const G4ParticleDefinition * theGamma
G4double currentInteractionLength
Definition: G4VProcess.hh:297
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
const G4MaterialCutsCouple * currentCouple
static constexpr double proton_mass_c2
G4double GetLocalEnergyDeposit() const
void SetLowestEnergyLimit(G4double)
const G4String & GetName() const
Definition: G4Material.hh:179
void SetVerboseLevel(G4int value)
Definition: G4VProcess.hh:440
void SetRangeTableForLoss(G4PhysicsTable *p)
static constexpr double TeV
Definition: G4SIunits.hh:218
G4int NumberOfModels() const
void Register(G4VEnergyLossProcess *p)
G4double LinearLossLimit() const
double G4double
Definition: G4Types.hh:76
G4double MaxEnergyForCSDARange() const
bool G4bool
Definition: G4Types.hh:79
Double_t scale
G4StepStatus GetStepStatus() const
TString part[npart]
Definition: Style.C:32
G4EmBiasingManager * biasManager
void AddEmModel(G4int, G4VEmModel *, G4VEmFluctuationModel *, const G4Region *)
G4double GetStepLimit(G4int coupleIdx, G4double previousStep)
G4VEmModel * GetModelByIndex(G4int idx=0, G4bool ver=false) const
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
void DefineMaterial(const G4MaterialCutsCouple *couple)
G4VParticleChange * pParticleChange
Definition: G4VProcess.hh:283
void SetIonisation(G4bool val)
virtual G4double SampleFluctuations(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmax, G4double length, G4double meanLoss)=0
virtual G4double GetContinuousStepLimit(const G4Track &track, G4double previousStepSize, G4double currentMinimumStep, G4double &currentSafety) override
G4PhysicsTable * InverseRangeTable() const
G4double MinSubRange() const
G4ProcessType
void SetCSDARangeTable(G4PhysicsTable *pRange)
G4double GetStepLength() const
G4double GetSubDEDXForScaledEnergy(G4double scaledKinEnergy)
void AddEmModel(G4int, G4VEmModel *, G4VEmFluctuationModel *fluc=0, const G4Region *region=nullptr)
G4VEmFluctuationModel * GetModelOfFluctuations()
Definition: G4VEmModel.hh:568
static constexpr double electron_mass_c2
G4double ScaledKinEnergyForLoss(G4double range)
const G4ThreeVector & GetPosition() const
G4double LowestMuHadEnergy() const
const G4DataVector * theCuts
virtual void InitialiseEnergyLossProcess(const G4ParticleDefinition *, const G4ParticleDefinition *)=0
void InitialiseBaseMaterials(G4PhysicsTable *table)
virtual G4VParticleChange * AlongStepDoIt(const G4Track &, const G4Step &) override
static constexpr double MeV
G4GPILSelection
#define G4UniformRand()
Definition: Randomize.hh:53
G4PhysicsTable * theIonisationSubTable
G4double GetGlobalTime() const
virtual void CorrectionsAlongStep(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &eloss, G4double &niel, G4double length)
Definition: G4VEmModel.cc:406
void SetSpline(G4bool)
const std::vector< G4double > * GetDensityFactors()
void FillSecondDerivatives()
G4int WorkerVerbose() const
G4double GetProposedKineticEnergy() const
G4PhysicsTable * theDEDXunRestrictedTable
void SetTouchableHandle(const G4TouchableHandle &apValue)
void SetLambdaTable(G4PhysicsTable *p)
G4PhysicsTable * BuildLambdaTable(G4EmTableType tType=fRestricted)
G4PhysicsVector * LambdaPhysicsVector(const G4MaterialCutsCouple *, G4double cut)
G4Track * GetTrack() const
Float_t d
const G4String & GetProcessName() const
Definition: G4VProcess.hh:411
const G4ThreeVector & GetPosition() const
G4GPILSelection aGPILSelection
G4double LambdaFactor() const
G4VSubCutProducer * SubCutProducer()
void DeRegister(G4VEnergyLossProcess *p)
void UpdateEmModel(const G4String &model_name, G4double emin, G4double emax)
void SetMaxKinEnergy(G4double e)
G4EmModelManager * modelManager
static G4Positron * Positron()
Definition: G4Positron.cc:94
G4double GetGlobalTime() const
void LocalPhysicsTables(const G4ParticleDefinition *aParticle, G4VEnergyLossProcess *p)
G4PhysicsTable * DEDXTableForSubsec() const
Definition: G4Step.hh:76
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &) override
G4StepPoint * GetPostStepPoint() const
G4ParticleChangeForLoss fParticleChange
G4Material * GetMaterial() const
const G4Material * currentMaterial
G4int Verbose() const
void InitialiseHelper()
static G4Electron * Electron()
Definition: G4Electron.cc:94
float bin[41]
Definition: plottest35.C:14
G4LossTableManager * lManager
const G4ParticleDefinition * secondaryParticle
G4LossTableBuilder * GetTableBuilder()
G4double MinKinEnergy() const
void ComputeLambdaForScaledEnergy(G4double scaledKinEnergy)
G4SafetyHelper * safetyHelper
virtual G4double AlongStepGetPhysicalInteractionLength(const G4Track &, G4double previousStepSize, G4double currentMinimumStep, G4double &currentSafety, G4GPILSelection *selection) override
std::vector< G4VEnergyLossProcess * > scProcesses
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
G4double theInitialNumberOfInteractionLength
Definition: G4VProcess.hh:300
static G4TransportationManager * GetTransportationManager()
G4VAtomDeexcitation * atomDeexcitation
G4bool StoreTable(const G4ParticleDefinition *p, G4PhysicsTable *, G4bool ascii, const G4String &directory, const G4String &tname)
static G4ProductionCutsTable * GetProductionCutsTable()
virtual G4double GetParticleCharge(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4VEmModel.cc:398
void SetSecondaryRangeTable(G4PhysicsTable *p)
G4PhysicsTable * theSubLambdaTable
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
void PrintWarning(G4String, G4double val)
G4double GetDEDXForScaledEnergy(G4double scaledKinEnergy)
G4bool ExistPhysicsTable(const G4String &fileName) const
const G4ParticleDefinition * particle
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:252
G4bool IsIonisationProcess() const
virtual G4bool StorePhysicsTable(const G4ParticleDefinition *, const G4String &directory, G4bool ascii=false) override
G4bool IsPIXEActive() const
virtual G4double MinPrimaryEnergy(const G4ParticleDefinition *, const G4Material *, G4double cut)
virtual G4double GetMeanFreePath(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition) override
int G4lrint(double ad)
Definition: templates.hh:151
G4double ContinuousStepLimit(const G4Track &track, G4double previousStepSize, G4double currentMinimumStep, G4double &currentSafety)
G4PhysicsTable * SecondaryRangeTable() const
double weight
Definition: plottest35.C:25
static G4RegionStore * GetInstance()
int G4int
Definition: G4Types.hh:78
G4VSubCutProducer * subcutProducer
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin=0.0, G4double tmax=DBL_MAX)=0
void SetParticleChange(G4VParticleChange *, G4VEmFluctuationModel *f=nullptr)
Definition: G4VEmModel.cc:454
G4ProcessVector * GetAtRestProcessVector(G4ProcessVectorTypeIndex typ=typeGPIL) const
G4ProcessManager * GetProcessManager() const
void ActivateSecondaryBiasing(const G4String &region, G4double factor, G4double energyLimit)
G4int verboseLevel
Definition: G4VProcess.hh:371
G4VEmModel * EmModel(size_t index=0) const
void SelectModel(G4double kinEnergy)
G4PhysicsTable * LambdaTable() const
void InitializeForPostStep(const G4Track &)
void SetProposedKineticEnergy(G4double proposedKinEnergy)
G4bool RetrieveTable(const G4ParticleDefinition *p, G4PhysicsTable *, G4bool ascii, const G4String &directory, const G4String &tname, G4bool mandatory)
G4double GetSafety() const
G4TrackStatus GetTrackStatus() const
void FillSecondariesAlongStep(G4double &eloss, G4double &weight)
G4VEnergyLossProcess(const G4String &name="EnergyLoss", G4ProcessType type=fElectromagnetic)
G4double CrossSectionPerVolume(G4double kineticEnergy, const G4MaterialCutsCouple *couple)
static constexpr double mm
Definition: SystemOfUnits.h:95
G4ForceCondition
virtual void SampleSecondaries(const G4Step &step, std::vector< G4Track * > &tracks, G4double &eloss, G4double cut) const =0
void FillLambdaVector(G4PhysicsVector *, const G4MaterialCutsCouple *, G4bool startFromNull=true, G4EmTableType t=fRestricted)
G4double GetKineticEnergy() const
virtual void BuildPhysicsTable(const G4ParticleDefinition &) override
G4int size() const
static constexpr double cm
Definition: G4SIunits.hh:119
G4GLOB_DLL std::ostream G4cout
std::vector< const G4Region * > scoffRegions
const G4int smax
static G4LossTableManager * Instance()
G4PhysicsTable * IonisationTable() const
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
G4ProductionCuts * GetProductionCuts() const
void clearAndDestroy()
G4PhysicsTable * CSDARangeTable() const
Char_t n[5]
const G4Element * GetCurrentElement() const
Definition: G4VEmModel.hh:477
G4double ComputeSafety(const G4ThreeVector &pGlobalPoint, G4double maxRadius=DBL_MAX)
const G4Material * GetMaterial() const
void SetInverseRangeTable(G4PhysicsTable *p)
static const G4double reg
G4PhysicsTable * DEDXunRestrictedTable() const
G4bool GetFlag(size_t i) const
G4double theNumberOfInteractionLengthLeft
Definition: G4VProcess.hh:293
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetCrossSectionBiasingFactor(G4double f, G4bool flag=true)
virtual G4double ChargeSquareRatio(const G4Track &)
Definition: G4VEmModel.cc:381
virtual G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition) override
G4double MeanFreePath(const G4Track &track)
G4bool LossFluctuation() const
void StreamInfo(std::ostream &out, const G4ParticleDefinition &part, G4String endOfLine=G4String("\n")) const
virtual void StreamProcessInfo(std::ostream &, G4String) const
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:609
#define DBL_MAX
Definition: templates.hh:83
const XML_Char XML_Content * model
Definition: expat.h:151
void ProposeTrackStatus(G4TrackStatus status)
G4double GetDEDXDispersion(const G4MaterialCutsCouple *couple, const G4DynamicParticle *dp, G4double length)
static constexpr double GeV
Definition: G4SIunits.hh:217
G4double GetLambdaForScaledEnergy(G4double scaledKinEnergy)
void UpdateEmModel(const G4String &, G4double, G4double)
const G4String & GetPhysicsTableFileName(const G4ParticleDefinition *, const G4String &directory, const G4String &tableName, G4bool ascii=false)
Definition: G4VProcess.cc:191
const G4DynamicParticle * GetDynamicParticle() const
G4double LowestElectronEnergy() const
G4bool Spline() const
G4PhysicsTable * theSecondaryRangeTable
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel)
Definition: G4VEmModel.cc:218
void RegisterExtraParticle(const G4ParticleDefinition *aParticle, G4VEnergyLossProcess *p)
void ActivateSubCutoff(G4bool val, const G4Region *region=nullptr)
static G4EmParameters * Instance()
size_t GetVectorLength() const
G4bool Integral() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
void SetLinearLossLimit(G4double val)
static G4bool RetrievePhysicsTable(G4PhysicsTable *physTable, const G4String &fileName, G4bool ascii)
G4double GetParentWeight() const
G4int GetParentID() const
virtual void ProcessDescription(std::ostream &outFile) const override
const G4ParticleDefinition * baseParticle
G4int GetProcessSubType() const
Definition: G4VProcess.hh:429
void SetNumberOfSecondaries(G4int totSecondaries)
static G4int Register(const G4String &)
size_t length() const