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