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