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