Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/utils/src/G4VEnergyLossProcess.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/utils/src/G4VEnergyLossProcess.cc (Version 11.3.0) and /processes/electromagnetic/utils/src/G4VEnergyLossProcess.cc (Version 5.2.p1)


  1 //                                                  1 
  2 // *******************************************    
  3 // * License and Disclaimer                       
  4 // *                                              
  5 // * The  Geant4 software  is  copyright of th    
  6 // * the Geant4 Collaboration.  It is provided    
  7 // * conditions of the Geant4 Software License    
  8 // * LICENSE and available at  http://cern.ch/    
  9 // * include a list of copyright holders.         
 10 // *                                              
 11 // * Neither the authors of this software syst    
 12 // * institutes,nor the agencies providing fin    
 13 // * work  make  any representation or  warran    
 14 // * regarding  this  software system or assum    
 15 // * use.  Please see the license in the file     
 16 // * for the full disclaimer and the limitatio    
 17 // *                                              
 18 // * This  code  implementation is the result     
 19 // * technical work of the GEANT4 collaboratio    
 20 // * By using,  copying,  modifying or  distri    
 21 // * any work based  on the software)  you  ag    
 22 // * use  in  resulting  scientific  publicati    
 23 // * acceptance of all terms of the Geant4 Sof    
 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 ca    
 43 // energy loss for charged particles using a s    
 44 // models valid for different energy regions.     
 45 // to create and access to dE/dx and range tab    
 46 // that information on fly.                       
 47 // -------------------------------------------    
 48 //                                                
 49 //....oooOO0OOooo........oooOO0OOooo........oo    
 50 //....oooOO0OOooo........oooOO0OOooo........oo    
 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 "G4EmParameters.hh"                      
 62 #include "G4EmUtility.hh"                         
 63 #include "G4EmTableUtil.hh"                       
 64 #include "G4VEmModel.hh"                          
 65 #include "G4VEmFluctuationModel.hh"               
 66 #include "G4DataVector.hh"                        
 67 #include "G4PhysicsLogVector.hh"                  
 68 #include "G4VParticleChange.hh"                   
 69 #include "G4Electron.hh"                          
 70 #include "G4ProcessManager.hh"                    
 71 #include "G4UnitsTable.hh"                        
 72 #include "G4Region.hh"                            
 73 #include "G4RegionStore.hh"                       
 74 #include "G4PhysicsTableHelper.hh"                
 75 #include "G4SafetyHelper.hh"                      
 76 #include "G4EmDataHandler.hh"                     
 77 #include "G4TransportationManager.hh"             
 78 #include "G4VAtomDeexcitation.hh"                 
 79 #include "G4VSubCutProducer.hh"                   
 80 #include "G4EmBiasingManager.hh"                  
 81 #include "G4Log.hh"                               
 82 #include <iostream>                               
 83                                                   
 84 //....oooOO0OOooo........oooOO0OOooo........oo    
 85                                                   
 86 namespace                                         
 87 {                                                 
 88   G4String tnames[7] =                            
 89     {"DEDX","Ionisation","DEDXnr","CSDARange",    
 90 }                                                 
 91                                                   
 92                                                   
 93 G4VEnergyLossProcess::G4VEnergyLossProcess(con    
 94                                            G4P    
 95   G4VContinuousDiscreteProcess(name, type)        
 96 {                                                 
 97   theParameters = G4EmParameters::Instance();     
 98   SetVerboseLevel(1);                             
 99                                                   
100   // low energy limit                             
101   lowestKinEnergy = theParameters->LowestElect    
102                                                   
103   // Size of tables                               
104   minKinEnergy     = 0.1*CLHEP::keV;              
105   maxKinEnergy     = 100.0*CLHEP::TeV;            
106   maxKinEnergyCSDA = 1.0*CLHEP::GeV;              
107   nBins            = 84;                          
108   nBinsCSDA        = 35;                          
109                                                   
110   invLambdaFactor = 1.0/lambdaFactor;             
111                                                   
112   // default linear loss limit                    
113   finalRange = 1.*CLHEP::mm;                      
114                                                   
115   // run time objects                             
116   pParticleChange = &fParticleChange;             
117   fParticleChange.SetSecondaryWeightByProcess(    
118   modelManager = new G4EmModelManager();          
119   safetyHelper = G4TransportationManager::GetT    
120     ->GetSafetyHelper();                          
121   aGPILSelection = CandidateForSelection;         
122                                                   
123   // initialise model                             
124   lManager = G4LossTableManager::Instance();      
125   lManager->Register(this);                       
126   isMaster = lManager->IsMaster();                
127                                                   
128   G4LossTableBuilder* bld = lManager->GetTable    
129   theDensityFactor = bld->GetDensityFactors();    
130   theDensityIdx = bld->GetCoupleIndexes();        
131                                                   
132   scTracks.reserve(10);                           
133   secParticles.reserve(12);                       
134   emModels = new std::vector<G4VEmModel*>;        
135 }                                                 
136                                                   
137 //....oooOO0OOooo........oooOO0OOooo........oo    
138                                                   
139 G4VEnergyLossProcess::~G4VEnergyLossProcess()     
140 {                                                 
141   if (isMaster) {                                 
142     if(nullptr == baseParticle) { delete theDa    
143     delete theEnergyOfCrossSectionMax;            
144     if(nullptr != fXSpeaks) {                     
145       for(auto const & v : *fXSpeaks) { delete    
146       delete fXSpeaks;                            
147     }                                             
148   }                                               
149   delete modelManager;                            
150   delete biasManager;                             
151   delete scoffRegions;                            
152   delete emModels;                                
153   lManager->DeRegister(this);                     
154 }                                                 
155                                                   
156 //....oooOO0OOooo........oooOO0OOooo........oo    
157                                                   
158 G4double G4VEnergyLossProcess::MinPrimaryEnerg    
159                                                   
160                                                   
161 {                                                 
162   return cut;                                     
163 }                                                 
164                                                   
165 //....oooOO0OOooo........oooOO0OOooo........oo    
166                                                   
167 void G4VEnergyLossProcess::AddEmModel(G4int or    
168                                       G4VEmFlu    
169                                       const G4    
170 {                                                 
171   if(nullptr == ptr) { return; }                  
172   G4VEmFluctuationModel* afluc = (nullptr == f    
173   modelManager->AddEmModel(order, ptr, afluc,     
174   ptr->SetParticleChange(pParticleChange, aflu    
175 }                                                 
176                                                   
177 //....oooOO0OOooo........oooOO0OOooo........oo    
178                                                   
179 void G4VEnergyLossProcess::SetEmModel(G4VEmMod    
180 {                                                 
181   if(nullptr == ptr) { return; }                  
182   if(!emModels->empty()) {                        
183     for(auto & em : *emModels) { if(em == ptr)    
184   }                                               
185   emModels->push_back(ptr);                       
186 }                                                 
187                                                   
188 //....oooOO0OOooo........oooOO0OOooo........oo    
189                                                   
190 void G4VEnergyLossProcess::SetDynamicMassCharg    
191                                                   
192 {                                                 
193   massRatio = massratio;                          
194   logMassRatio = G4Log(massRatio);                
195   fFactor = charge2ratio*biasFactor;              
196   if(baseMat) { fFactor *= (*theDensityFactor)    
197   chargeSqRatio = charge2ratio;                   
198   reduceFactor  = 1.0/(fFactor*massRatio);        
199 }                                                 
200                                                   
201 //....oooOO0OOooo........oooOO0OOooo........oo    
202                                                   
203 void                                              
204 G4VEnergyLossProcess::PreparePhysicsTable(cons    
205 {                                                 
206   particle = G4EmTableUtil::CheckIon(this, &pa    
207                                      verboseLe    
208                                                   
209   if( particle != &part ) {                       
210     if(!isIon) { lManager->RegisterExtraPartic    
211     if(1 < verboseLevel) {                        
212       G4cout << "### G4VEnergyLossProcess::Pre    
213              << " interrupted for " << GetProc    
214              << part.GetParticleName() << " is    
215              << " spline=" << spline << G4endl    
216     }                                             
217     return;                                       
218   }                                               
219                                                   
220   tablesAreBuilt = false;                         
221   if (GetProcessSubType() == fIonisation) { Se    
222                                                   
223   G4LossTableBuilder* bld = lManager->GetTable    
224   lManager->PreparePhysicsTable(&part, this);     
225                                                   
226   // Base particle and set of models can be de    
227   InitialiseEnergyLossProcess(particle, basePa    
228                                                   
229   // parameters of the process                    
230   if(!actLossFluc) { lossFluctuationFlag = the    
231   useCutAsFinalRange = theParameters->UseCutAs    
232   if(!actMinKinEnergy) { minKinEnergy = thePar    
233   if(!actMaxKinEnergy) { maxKinEnergy = thePar    
234   if(!actBinning) { nBins = theParameters->Num    
235   maxKinEnergyCSDA = theParameters->MaxEnergyF    
236   nBinsCSDA = theParameters->NumberOfBinsPerDe    
237     *G4lrint(std::log10(maxKinEnergyCSDA/minKi    
238   if(!actLinLossLimit) { linLossLimit = thePar    
239   lambdaFactor = theParameters->LambdaFactor()    
240   invLambdaFactor = 1.0/lambdaFactor;             
241   if(isMaster) { SetVerboseLevel(theParameters    
242   else { SetVerboseLevel(theParameters->Worker    
243   // integral option may be disabled              
244   if(!theParameters->Integral()) { fXSType = f    
245                                                   
246   theParameters->DefineRegParamForLoss(this);     
247                                                   
248   fRangeEnergy = 0.0;                             
249                                                   
250   G4double initialCharge = particle->GetPDGCha    
251   G4double initialMass   = particle->GetPDGMas    
252                                                   
253   theParameters->FillStepFunction(particle, th    
254                                                   
255   // parameters for scaling from the base part    
256   if (nullptr != baseParticle) {                  
257     massRatio    = (baseParticle->GetPDGMass()    
258     logMassRatio = G4Log(massRatio);              
259     G4double q = initialCharge/baseParticle->G    
260     chargeSqRatio = q*q;                          
261     if(chargeSqRatio > 0.0) { reduceFactor = 1    
262   }                                               
263   lowestKinEnergy = (initialMass < CLHEP::MeV)    
264     ? theParameters->LowestElectronEnergy()       
265     : theParameters->LowestMuHadEnergy();         
266                                                   
267   // Tables preparation                           
268   if (isMaster && nullptr == baseParticle) {      
269     if(nullptr == theData) { theData = new G4E    
270                                                   
271     if(nullptr != theDEDXTable && isIonisation    
272       if(nullptr != theIonisationTable && theD    
273   theData->CleanTable(0);                         
274   theDEDXTable = theIonisationTable;              
275   theIonisationTable = nullptr;                   
276       }                                           
277     }                                             
278                                                   
279     theDEDXTable = theData->MakeTable(theDEDXT    
280     bld->InitialiseBaseMaterials(theDEDXTable)    
281     theData->UpdateTable(theIonisationTable, 1    
282                                                   
283     if (theParameters->BuildCSDARange()) {        
284       theDEDXunRestrictedTable = theData->Make    
285       if(isIonisation) { theCSDARangeTable = t    
286     }                                             
287                                                   
288     theLambdaTable = theData->MakeTable(4);       
289     if(isIonisation) {                            
290       theRangeTableForLoss = theData->MakeTabl    
291       theInverseRangeTable = theData->MakeTabl    
292     }                                             
293   }                                               
294                                                   
295   // forced biasing                               
296   if(nullptr != biasManager) {                    
297     biasManager->Initialise(part,GetProcessNam    
298     biasFlag = false;                             
299   }                                               
300   baseMat = bld->GetBaseMaterialFlag();           
301   numberOfModels = modelManager->NumberOfModel    
302   currentModel = modelManager->GetModel(0);       
303   G4EmTableUtil::UpdateModels(this, modelManag    
304                               numberOfModels,     
305                               mainSecondaries,    
306                               theParameters->U    
307   theCuts = modelManager->Initialise(particle,    
308                                      verboseLe    
309   // subcut processor                             
310   if(isIonisation) {                              
311     subcutProducer = lManager->SubCutProducer(    
312   }                                               
313   if(1 == nSCoffRegions) {                        
314     if((*scoffRegions)[0]->GetName() == "Defau    
315       delete scoffRegions;                        
316       scoffRegions = nullptr;                     
317       nSCoffRegions = 0;                          
318     }                                             
319   }                                               
320                                                   
321   if(1 < verboseLevel) {                          
322     G4cout << "G4VEnergyLossProcess::PrepearPh    
323            << " for " << GetProcessName() << "    
324            << " isIon= " << isIon << " spline=    
325     if(baseParticle) {                            
326       G4cout << "; base: " << baseParticle->Ge    
327     }                                             
328     G4cout << G4endl;                             
329     G4cout << " chargeSqRatio= " << chargeSqRa    
330            << " massRatio= " << massRatio         
331            << " reduceFactor= " << reduceFacto    
332     if (nSCoffRegions > 0) {                      
333       G4cout << " SubCut secondary production     
334       for (G4int i=0; i<nSCoffRegions; ++i) {     
335         const G4Region* r = (*scoffRegions)[i]    
336         G4cout << "           " << r->GetName(    
337       }                                           
338     } else if(nullptr != subcutProducer) {        
339       G4cout << " SubCut secondary production     
340     }                                             
341   }                                               
342 }                                                 
343                                                   
344 //....oooOO0OOooo........oooOO0OOooo........oo    
345                                                   
346 void G4VEnergyLossProcess::BuildPhysicsTable(c    
347 {                                                 
348   if(1 < verboseLevel) {                          
349     G4cout << "### G4VEnergyLossProcess::Build    
350            << GetProcessName()                    
351            << " and particle " << part.GetPart    
352            << "; the first particle " << parti    
353     if(baseParticle) {                            
354       G4cout << "; base: " << baseParticle->Ge    
355     }                                             
356     G4cout << G4endl;                             
357     G4cout << "    TablesAreBuilt= " << tables    
358            << " spline=" << spline << " ptr: "    
359   }                                               
360                                                   
361   if(&part == particle) {                         
362     if(isMaster) {                                
363       lManager->BuildPhysicsTable(particle, th    
364                                                   
365     } else {                                      
366       const auto masterProcess =                  
367         static_cast<const G4VEnergyLossProcess    
368                                                   
369       numberOfModels = modelManager->NumberOfM    
370       G4EmTableUtil::BuildLocalElossProcess(th    
371                                             pa    
372       tablesAreBuilt = true;                      
373       baseMat = masterProcess->UseBaseMaterial    
374       lManager->LocalPhysicsTables(particle, t    
375     }                                             
376                                                   
377     // needs to be done only once                 
378     safetyHelper->InitialiseHelper();             
379   }                                               
380   // Added tracking cut to avoid tracking arti    
381   // and identified deexcitation flag             
382   if(isIonisation) {                              
383     atomDeexcitation = lManager->AtomDeexcitat    
384     if(nullptr != atomDeexcitation) {             
385       if(atomDeexcitation->IsPIXEActive()) { u    
386     }                                             
387   }                                               
388                                                   
389   // protection against double printout           
390   if(theParameters->IsPrintLocked()) { return;    
391                                                   
392   // explicitly defined printout by particle n    
393   G4String num = part.GetParticleName();          
394   if(1 < verboseLevel ||                          
395      (0 < verboseLevel && (num == "e-" ||         
396                            num == "e+"    || n    
397                            num == "mu-"   || n    
398                            num == "pi+"   || n    
399                            num == "kaon+" || n    
400                            num == "alpha" || n    
401                            num == "GenericIon"    
402     StreamInfo(G4cout, part);                     
403   }                                               
404   if(1 < verboseLevel) {                          
405     G4cout << "### G4VEnergyLossProcess::Build    
406            << GetProcessName()                    
407            << " and particle " << part.GetPart    
408     if(isIonisation) { G4cout << "  isIonisati    
409     G4cout << " baseMat=" << baseMat << G4endl    
410   }                                               
411 }                                                 
412                                                   
413 //....oooOO0OOooo........oooOO0OOooo........oo    
414                                                   
415 G4PhysicsTable* G4VEnergyLossProcess::BuildDED    
416 {                                                 
417   G4PhysicsTable* table = nullptr;                
418   G4double emax = maxKinEnergy;                   
419   G4int bin = nBins;                              
420                                                   
421   if(fTotal == tType) {                           
422     emax  = maxKinEnergyCSDA;                     
423     bin   = nBinsCSDA;                            
424     table = theDEDXunRestrictedTable;             
425   } else if(fRestricted == tType) {               
426     table = theDEDXTable;                         
427   } else {                                        
428     G4cout << "G4VEnergyLossProcess::BuildDEDX    
429            << tType << G4endl;                    
430   }                                               
431   if(1 < verboseLevel) {                          
432     G4cout << "G4VEnergyLossProcess::BuildDEDX    
433            << " for " << GetProcessName()         
434            << " and " << particle->GetParticle    
435      << "spline=" << spline << G4endl;            
436   }                                               
437   if(nullptr == table) { return table; }          
438                                                   
439   G4LossTableBuilder* bld = lManager->GetTable    
440   G4EmTableUtil::BuildDEDXTable(this, particle    
441                                 table, minKinE    
442                                 verboseLevel,     
443   return table;                                   
444 }                                                 
445                                                   
446 //....oooOO0OOooo........oooOO0OOooo........oo    
447                                                   
448 G4PhysicsTable* G4VEnergyLossProcess::BuildLam    
449 {                                                 
450   if(nullptr == theLambdaTable) { return theLa    
451                                                   
452   G4double scale = theParameters->MaxKinEnergy    
453   G4int nbin =                                    
454     theParameters->NumberOfBinsPerDecade()*G4l    
455   scale = nbin/G4Log(scale);                      
456                                                   
457   G4LossTableBuilder* bld = lManager->GetTable    
458   G4EmTableUtil::BuildLambdaTable(this, partic    
459                                   bld, theLamb    
460                                   minKinEnergy    
461                                   verboseLevel    
462   return theLambdaTable;                          
463 }                                                 
464                                                   
465 //....oooOO0OOooo........oooOO0OOooo........oo    
466                                                   
467 void G4VEnergyLossProcess::StreamInfo(std::ost    
468                 const G4ParticleDefinition& pa    
469 {                                                 
470   G4String indent = (rst ? "  " : "");            
471   out << std::setprecision(6);                    
472   out << G4endl << indent << GetProcessName()     
473   if (!rst) out << " for " << part.GetParticle    
474   out << "  XStype:" << fXSType                   
475       << "  SubType=" << GetProcessSubType() <    
476       << "      dE/dx and range tables from "     
477       << G4BestUnit(minKinEnergy,"Energy")        
478       << " to " << G4BestUnit(maxKinEnergy,"En    
479       << " in " << nBins << " bins" << G4endl     
480       << "      Lambda tables from threshold t    
481       << G4BestUnit(maxKinEnergy,"Energy")        
482       << ", " << theParameters->NumberOfBinsPe    
483       << " bins/decade, spline: " << spline       
484       << G4endl;                                  
485   if(nullptr != theRangeTableForLoss && isIoni    
486     out << "      StepFunction=(" << dRoverRan    
487         << finalRange/mm << " mm)"                
488         << ", integ: " << fXSType                 
489         << ", fluct: " << lossFluctuationFlag     
490         << ", linLossLim= " << linLossLimit       
491         << G4endl;                                
492   }                                               
493   StreamProcessInfo(out);                         
494   modelManager->DumpModelList(out, verboseLeve    
495   if(nullptr != theCSDARangeTable && isIonisat    
496     out << "      CSDA range table up"            
497         << " to " << G4BestUnit(maxKinEnergyCS    
498         << " in " << nBinsCSDA << " bins" << G    
499   }                                               
500   if(nSCoffRegions>0 && isIonisation) {           
501     out << "      Subcutoff sampling in " << n    
502         << " regions" << G4endl;                  
503   }                                               
504   if(2 < verboseLevel) {                          
505     for(std::size_t i=0; i<7; ++i) {              
506       auto ta = theData->Table(i);                
507       out << "      " << tnames[i] << " addres    
508       if(nullptr != ta) { out << *ta << G4endl    
509     }                                             
510   }                                               
511 }                                                 
512                                                   
513 //....oooOO0OOooo........oooOO0OOooo........oo    
514                                                   
515 void G4VEnergyLossProcess::ActivateSubCutoff(c    
516 {                                                 
517   if(nullptr == scoffRegions) {                   
518     scoffRegions = new std::vector<const G4Reg    
519   }                                               
520   // the region is in the list                    
521   if(!scoffRegions->empty()) {                    
522     for (auto & reg : *scoffRegions) {            
523       if (reg == r) { return; }                   
524     }                                             
525   }                                               
526   // new region                                   
527   scoffRegions->push_back(r);                     
528   ++nSCoffRegions;                                
529 }                                                 
530                                                   
531 //....oooOO0OOooo........oooOO0OOooo........oo    
532                                                   
533 G4bool G4VEnergyLossProcess::IsRegionForCubcut    
534 {                                                 
535   if(0 == nSCoffRegions) { return true; }         
536   const G4Region* r = aTrack.GetVolume()->GetL    
537   for(auto & reg : *scoffRegions) {               
538     if(r == reg) { return true; }                 
539   }                                               
540   return false;                                   
541 }                                                 
542                                                   
543 //....oooOO0OOooo........oooOO0OOooo........oo    
544                                                   
545 void G4VEnergyLossProcess::StartTracking(G4Tra    
546 {                                                 
547   // reset parameters for the new track           
548   theNumberOfInteractionLengthLeft = -1.0;        
549   mfpKinEnergy = DBL_MAX;                         
550   preStepLambda = 0.0;                            
551   currentCouple = nullptr;                        
552                                                   
553   // reset ion                                    
554   if(isIon) {                                     
555     const G4double newmass = track->GetDefinit    
556     massRatio = (nullptr == baseParticle) ? CL    
557       : baseParticle->GetPDGMass()/newmass;       
558     logMassRatio = G4Log(massRatio);              
559   }                                               
560   // forced biasing only for primary particles    
561   if(nullptr != biasManager) {                    
562     if(0 == track->GetParentID()) {               
563       biasFlag = true;                            
564       biasManager->ResetForcedInteraction();      
565     }                                             
566   }                                               
567 }                                                 
568                                                   
569 //....oooOO0OOooo........oooOO0OOooo........oo    
570                                                   
571 G4double G4VEnergyLossProcess::AlongStepGetPhy    
572                              const G4Track& tr    
573                              G4GPILSelection*     
574 {                                                 
575   G4double x = DBL_MAX;                           
576   *selection = aGPILSelection;                    
577   if(isIonisation && currentModel->IsActive(pr    
578     GetScaledRangeForScaledEnergy(preStepScale    
579     x = (useCutAsFinalRange) ? std::min(finalR    
580       currentCouple->GetProductionCuts()->GetP    
581     x = (fRange > x) ? fRange*dRoverRange + x*    
582       : fRange;                                   
583     /*                                            
584       G4cout<<"AlongStepGPIL: " << GetProcessN    
585   << " fRange=" << fRange << " finR=" << finR     
586     */                                            
587   }                                               
588   return x;                                       
589 }                                                 
590                                                   
591 //....oooOO0OOooo........oooOO0OOooo........oo    
592                                                   
593 G4double G4VEnergyLossProcess::PostStepGetPhys    
594                              const G4Track& tr    
595                              G4double   previo    
596                              G4ForceCondition*    
597 {                                                 
598   // condition is set to "Not Forced"             
599   *condition = NotForced;                         
600   G4double x = DBL_MAX;                           
601                                                   
602   // initialisation of material, mass, charge,    
603   // at the beginning of the step                 
604   DefineMaterial(track.GetMaterialCutsCouple()    
605   preStepKinEnergy       = track.GetKineticEne    
606   preStepScaledEnergy    = preStepKinEnergy*ma    
607   SelectModel(preStepScaledEnergy);               
608                                                   
609   if(!currentModel->IsActive(preStepScaledEner    
610     theNumberOfInteractionLengthLeft = -1.0;      
611     mfpKinEnergy = DBL_MAX;                       
612     preStepLambda = 0.0;                          
613     currentInteractionLength = DBL_MAX;           
614     return x;                                     
615   }                                               
616                                                   
617   // change effective charge of a charged part    
618   if(isIon) {                                     
619     const G4double q2 = currentModel->ChargeSq    
620     fFactor = q2*biasFactor;                      
621     if(baseMat) { fFactor *= (*theDensityFacto    
622     reduceFactor = 1.0/(fFactor*massRatio);       
623     if (lossFluctuationFlag) {                    
624       auto fluc = currentModel->GetModelOfFluc    
625       fluc->SetParticleAndCharge(track.GetDefi    
626     }                                             
627   }                                               
628                                                   
629   // forced biasing only for primary particles    
630   if(biasManager) {                               
631     if(0 == track.GetParentID() && biasFlag &&    
632        biasManager->ForcedInteractionRegion((G    
633       return biasManager->GetStepLimit((G4int)    
634     }                                             
635   }                                               
636                                                   
637   ComputeLambdaForScaledEnergy(preStepScaledEn    
638                                                   
639   // zero cross section                           
640   if(preStepLambda <= 0.0) {                      
641     theNumberOfInteractionLengthLeft = -1.0;      
642     currentInteractionLength = DBL_MAX;           
643   } else {                                        
644                                                   
645     // non-zero cross section                     
646     if (theNumberOfInteractionLengthLeft < 0.0    
647                                                   
648       // beggining of tracking (or just after     
649       theNumberOfInteractionLengthLeft = -G4Lo    
650       theInitialNumberOfInteractionLength = th    
651                                                   
652     } else if(currentInteractionLength < DBL_M    
653                                                   
654       // subtract NumberOfInteractionLengthLef    
655       theNumberOfInteractionLengthLeft -=         
656         previousStepSize/currentInteractionLen    
657                                                   
658       theNumberOfInteractionLengthLeft =          
659         std::max(theNumberOfInteractionLengthL    
660     }                                             
661                                                   
662     // new mean free path and step limit          
663     currentInteractionLength = 1.0/preStepLamb    
664     x = theNumberOfInteractionLengthLeft * cur    
665   }                                               
666 #ifdef G4VERBOSE                                  
667   if (verboseLevel>2) {                           
668     G4cout << "G4VEnergyLossProcess::PostStepG    
669     G4cout << "[ " << GetProcessName() << "]"     
670     G4cout << " for " << track.GetDefinition()    
671            << " in Material  " <<  currentMate    
672            << " Ekin(MeV)= " << preStepKinEner    
673            << " track material: " << track.Get    
674            <<G4endl;                              
675     G4cout << "MeanFreePath = " << currentInte    
676            << "InteractionLength= " << x/cm <<    
677   }                                               
678 #endif                                            
679   return x;                                       
680 }                                                 
681                                                   
682 //....oooOO0OOooo........oooOO0OOooo........oo    
683                                                   
684 void                                              
685 G4VEnergyLossProcess::ComputeLambdaForScaledEn    
686 {                                                 
687   // cross section increased with energy          
688   if(fXSType == fEmIncreasing) {                  
689     if(e*invLambdaFactor < mfpKinEnergy) {        
690       preStepLambda = GetLambdaForScaledEnergy    
691       mfpKinEnergy = (preStepLambda > 0.0) ? e    
692     }                                             
693                                                   
694     // cross section has one peak                 
695   } else if(fXSType == fEmOnePeak) {              
696     const G4double epeak = (*theEnergyOfCrossS    
697     if(e <= epeak) {                              
698       if(e*invLambdaFactor < mfpKinEnergy) {      
699         preStepLambda = GetLambdaForScaledEner    
700         mfpKinEnergy = (preStepLambda > 0.0) ?    
701       }                                           
702     } else if(e < mfpKinEnergy) {                 
703       const G4double e1 = std::max(epeak, e*la    
704       mfpKinEnergy = e1;                          
705       preStepLambda = GetLambdaForScaledEnergy    
706     }                                             
707                                                   
708     // cross section has more than one peaks      
709   } else if(fXSType == fEmTwoPeaks) {             
710     G4TwoPeaksXS* xs = (*fXSpeaks)[basedCouple    
711     const G4double e1peak = xs->e1peak;           
712                                                   
713     // below the 1st peak                         
714     if(e <= e1peak) {                             
715       if(e*invLambdaFactor < mfpKinEnergy) {      
716         preStepLambda = GetLambdaForScaledEner    
717         mfpKinEnergy = (preStepLambda > 0.0) ?    
718       }                                           
719       return;                                     
720     }                                             
721     const G4double e1deep = xs->e1deep;           
722     // above the 1st peak, below the deep         
723     if(e <= e1deep) {                             
724       if(mfpKinEnergy >= e1deep || e <= mfpKin    
725         const G4double e1 = std::max(e1peak, e    
726         mfpKinEnergy = e1;                        
727         preStepLambda = GetLambdaForScaledEner    
728       }                                           
729       return;                                     
730     }                                             
731     const G4double e2peak = xs->e2peak;           
732     // above the deep, below 2nd peak             
733     if(e <= e2peak) {                             
734       if(e*invLambdaFactor < mfpKinEnergy) {      
735         mfpKinEnergy = e;                         
736         preStepLambda = GetLambdaForScaledEner    
737       }                                           
738       return;                                     
739     }                                             
740     const G4double e2deep = xs->e2deep;           
741     // above the 2nd peak, below the deep         
742     if(e <= e2deep) {                             
743       if(mfpKinEnergy >= e2deep || e <= mfpKin    
744         const G4double e1 = std::max(e2peak, e    
745         mfpKinEnergy = e1;                        
746         preStepLambda = GetLambdaForScaledEner    
747       }                                           
748       return;                                     
749     }                                             
750     const G4double e3peak = xs->e3peak;           
751     // above the deep, below 3d peak              
752     if(e <= e3peak) {                             
753       if(e*invLambdaFactor < mfpKinEnergy) {      
754         mfpKinEnergy = e;                         
755         preStepLambda = GetLambdaForScaledEner    
756       }                                           
757       return;                                     
758     }                                             
759     // above 3d peak                              
760     if(e <= mfpKinEnergy) {                       
761       const G4double e1 = std::max(e3peak, e*l    
762       mfpKinEnergy = e1;                          
763       preStepLambda = GetLambdaForScaledEnergy    
764     }                                             
765     // integral method is not used                
766   } else {                                        
767     preStepLambda = GetLambdaForScaledEnergy(e    
768   }                                               
769 }                                                 
770                                                   
771 //....oooOO0OOooo........oooOO0OOooo........oo    
772                                                   
773 G4VParticleChange* G4VEnergyLossProcess::Along    
774                                                   
775 {                                                 
776   fParticleChange.InitializeForAlongStep(track    
777   // The process has range table - calculate e    
778   if(!isIonisation || !currentModel->IsActive(    
779     return &fParticleChange;                      
780   }                                               
781                                                   
782   G4double length = step.GetStepLength();         
783   G4double eloss  = 0.0;                          
784                                                   
785   /*                                              
786   if(-1 < verboseLevel) {                         
787     const G4ParticleDefinition* d = track.GetP    
788     G4cout << "AlongStepDoIt for "                
789            << GetProcessName() << " and partic    
790            << "  eScaled(MeV)=" << preStepScal    
791            << "  range(mm)=" << fRange/mm << "    
792            << "  rf=" << reduceFactor << "  q^    
793            << " md=" << d->GetPDGMass() << "      
794            << "  " << track.GetMaterial()->Get    
795   }                                               
796   */                                              
797   const G4DynamicParticle* dynParticle = track    
798                                                   
799   // define new weight for primary and seconda    
800   G4double weight = fParticleChange.GetParentW    
801   if(weightFlag) {                                
802     weight /= biasFactor;                         
803     fParticleChange.ProposeWeight(weight);        
804   }                                               
805                                                   
806   // stopping, check actual range and kinetic     
807   if (length >= fRange || preStepKinEnergy <=     
808     eloss = preStepKinEnergy;                     
809     if (useDeexcitation) {                        
810       atomDeexcitation->AlongStepDeexcitation(    
811                                                   
812       if(scTracks.size() > 0) { FillSecondarie    
813       eloss = std::max(eloss, 0.0);               
814     }                                             
815     fParticleChange.SetProposedKineticEnergy(0    
816     fParticleChange.ProposeLocalEnergyDeposit(    
817     return &fParticleChange;                      
818   }                                               
819   // zero step length with non-zero range         
820   if(length <= 0.0) { return &fParticleChange;    
821                                                   
822   // Short step                                   
823   eloss = length*GetDEDXForScaledEnergy(preSte    
824                                         LogSca    
825   /*                                              
826   G4cout << "##### Short STEP: eloss= " << elo    
827    << " Escaled=" << preStepScaledEnergy          
828    << " R=" << fRange                             
829    << " L=" << length                             
830    << " fFactor=" << fFactor << " minE=" << mi    
831    << " idxBase=" << basedCoupleIndex << G4end    
832   */                                              
833   // Long step                                    
834   if(eloss > preStepKinEnergy*linLossLimit) {     
835                                                   
836     const G4double x = (fRange - length)/reduc    
837     const G4double de = preStepKinEnergy - Sca    
838     if(de > 0.0) { eloss = de; }                  
839     /*                                            
840     if(-1 < verboseLevel)                         
841       G4cout << "  Long STEP: rPre(mm)="          
842              << GetScaledRangeForScaledEnergy(    
843              << " x(mm)=" << x/mm                 
844              << " eloss(MeV)=" << eloss/MeV       
845        << " rFactor=" << reduceFactor             
846        << " massRatio=" << massRatio              
847              << G4endl;                           
848     */                                            
849   }                                               
850                                                   
851   /*                                              
852   if(-1 < verboseLevel ) {                        
853     G4cout << "Before fluct: eloss(MeV)= " <<     
854            << " e-eloss= " << preStepKinEnergy    
855            << " step(mm)= " << length/mm << "     
856            << " fluct= " << lossFluctuationFla    
857   }                                               
858   */                                              
859                                                   
860   const G4double cut = (*theCuts)[currentCoupl    
861   G4double esec = 0.0;                            
862                                                   
863   // Corrections, which cannot be tabulated       
864   if(isIon) {                                     
865     currentModel->CorrectionsAlongStep(current    
866                                        length,    
867     eloss = std::max(eloss, 0.0);                 
868   }                                               
869                                                   
870   // Sample fluctuations if not full energy lo    
871   if(eloss >= preStepKinEnergy) {                 
872     eloss = preStepKinEnergy;                     
873                                                   
874   } else if (lossFluctuationFlag) {               
875     const G4double tmax = currentModel->MaxSec    
876     const G4double tcut = std::min(cut, tmax);    
877     G4VEmFluctuationModel* fluc = currentModel    
878     eloss = fluc->SampleFluctuations(currentCo    
879                                      tcut, tma    
880     /*                                            
881     if(-1 < verboseLevel)                         
882       G4cout << "After fluct: eloss(MeV)= " <<    
883              << " fluc= " << (eloss-eloss0)/Me    
884              << " ChargeSqRatio= " << chargeSq    
885              << " massRatio= " << massRatio <<    
886     */                                            
887   }                                               
888                                                   
889   // deexcitation                                 
890   if (useDeexcitation) {                          
891     G4double esecfluo = preStepKinEnergy;         
892     G4double de = esecfluo;                       
893     atomDeexcitation->AlongStepDeexcitation(sc    
894                                             de    
895                                                   
896     // sum of de-excitation energies              
897     esecfluo -= de;                               
898                                                   
899     // subtracted from energy loss                
900     if(eloss >= esecfluo) {                       
901       esec  += esecfluo;                          
902       eloss -= esecfluo;                          
903     } else {                                      
904       esec += esecfluo;                           
905       eloss = 0.0;                                
906     }                                             
907   }                                               
908   if(nullptr != subcutProducer && IsRegionForC    
909     subcutProducer->SampleSecondaries(step, sc    
910   }                                               
911   // secondaries from atomic de-excitation and    
912   if(!scTracks.empty()) { FillSecondariesAlong    
913                                                   
914   // Energy balance                               
915   G4double finalT = preStepKinEnergy - eloss -    
916   if (finalT <= lowestKinEnergy) {                
917     eloss += finalT;                              
918     finalT = 0.0;                                 
919   } else if(isIon) {                              
920     fParticleChange.SetProposedCharge(            
921       currentModel->GetParticleCharge(track.Ge    
922                                       currentM    
923   }                                               
924   eloss = std::max(eloss, 0.0);                   
925                                                   
926   fParticleChange.SetProposedKineticEnergy(fin    
927   fParticleChange.ProposeLocalEnergyDeposit(el    
928   /*                                              
929   if(-1 < verboseLevel) {                         
930     G4double del = finalT + eloss + esec - pre    
931     G4cout << "Final value eloss(MeV)= " << el    
932            << " preStepKinEnergy= " << preStep    
933            << " postStepKinEnergy= " << finalT    
934            << " de(keV)= " << del/keV             
935            << " lossFlag= " << lossFluctuation    
936            << "  status= " << track.GetTrackSt    
937            << G4endl;                             
938   }                                               
939   */                                              
940   return &fParticleChange;                        
941 }                                                 
942                                                   
943 //....oooOO0OOooo........oooOO0OOooo........oo    
944                                                   
945 void G4VEnergyLossProcess::FillSecondariesAlon    
946 {                                                 
947   const std::size_t n0 = scTracks.size();         
948   G4double weight = wt;                           
949   // weight may be changed by biasing manager     
950   if(biasManager) {                               
951     if(biasManager->SecondaryBiasingRegion((G4    
952       weight *=                                   
953         biasManager->ApplySecondaryBiasing(scT    
954     }                                             
955   }                                               
956                                                   
957   // fill secondaries                             
958   const std::size_t n = scTracks.size();          
959   fParticleChange.SetNumberOfSecondaries((G4in    
960                                                   
961   for(std::size_t i=0; i<n; ++i) {                
962     G4Track* t = scTracks[i];                     
963     if(nullptr != t) {                            
964       t->SetWeight(weight);                       
965       pParticleChange->AddSecondary(t);           
966       G4int pdg = t->GetDefinition()->GetPDGEn    
967       if (i < n0) {                               
968         if (pdg == 22) {                          
969     t->SetCreatorModelID(gpixeID);                
970         } else if (pdg == 11) {                   
971           t->SetCreatorModelID(epixeID);          
972         } else {                                  
973           t->SetCreatorModelID(biasID);           
974   }                                               
975       } else {                                    
976   t->SetCreatorModelID(biasID);                   
977       }                                           
978     }                                             
979   }                                               
980   scTracks.clear();                               
981 }                                                 
982                                                   
983 //....oooOO0OOooo........oooOO0OOooo........oo    
984                                                   
985 G4VParticleChange* G4VEnergyLossProcess::PostS    
986                                                   
987 {                                                 
988   // clear number of interaction lengths in an    
989   theNumberOfInteractionLengthLeft = -1.0;        
990   mfpKinEnergy = DBL_MAX;                         
991                                                   
992   fParticleChange.InitializeForPostStep(track)    
993   const G4double finalT = track.GetKineticEner    
994                                                   
995   const G4double postStepScaledEnergy = finalT    
996   SelectModel(postStepScaledEnergy);              
997                                                   
998   if(!currentModel->IsActive(postStepScaledEne    
999     return &fParticleChange;                      
1000   }                                              
1001   /*                                             
1002   if(1 < verboseLevel) {                         
1003     G4cout<<GetProcessName()<<" PostStepDoIt:    
1004   }                                              
1005   */                                             
1006   // forced process - should happen only once    
1007   if(biasFlag) {                                 
1008     if(biasManager->ForcedInteractionRegion((    
1009       biasFlag = false;                          
1010     }                                            
1011   }                                              
1012   const G4DynamicParticle* dp = track.GetDyna    
1013                                                  
1014   // Integral approach                           
1015   if (fXSType != fEmNoIntegral) {                
1016     const G4double logFinalT = dp->GetLogKine    
1017     G4double lx = GetLambdaForScaledEnergy(po    
1018                                            lo    
1019     lx = std::max(lx, 0.0);                      
1020                                                  
1021     // if both lg and lx are zero then no int    
1022     if(preStepLambda*G4UniformRand() >= lx) {    
1023       return &fParticleChange;                   
1024     }                                            
1025   }                                              
1026                                                  
1027   // define new weight for primary and second    
1028   G4double weight = fParticleChange.GetParent    
1029   if(weightFlag) {                               
1030     weight /= biasFactor;                        
1031     fParticleChange.ProposeWeight(weight);       
1032   }                                              
1033                                                  
1034   const G4double tcut = (*theCuts)[currentCou    
1035                                                  
1036   // sample secondaries                          
1037   secParticles.clear();                          
1038   currentModel->SampleSecondaries(&secParticl    
1039                                                  
1040   const G4int num0 = (G4int)secParticles.size    
1041                                                  
1042   // bremsstrahlung splitting or Russian roul    
1043   if(biasManager) {                              
1044     if(biasManager->SecondaryBiasingRegion((G    
1045       G4double eloss = 0.0;                      
1046       weight *= biasManager->ApplySecondaryBi    
1047                                       secPart    
1048                                       track,     
1049                                       &fParti    
1050                                       (G4int)    
1051                                       step.Ge    
1052       if(eloss > 0.0) {                          
1053         eloss += fParticleChange.GetLocalEner    
1054         fParticleChange.ProposeLocalEnergyDep    
1055       }                                          
1056     }                                            
1057   }                                              
1058                                                  
1059   // save secondaries                            
1060   const G4int num = (G4int)secParticles.size(    
1061   if(num > 0) {                                  
1062                                                  
1063     fParticleChange.SetNumberOfSecondaries(nu    
1064     G4double time = track.GetGlobalTime();       
1065                                                  
1066     G4int n1(0), n2(0);                          
1067     if(num0 > mainSecondaries) {                 
1068       currentModel->FillNumberOfSecondaries(n    
1069     }                                            
1070                                                  
1071     for (G4int i=0; i<num; ++i) {                
1072       if(nullptr != secParticles[i]) {           
1073         G4Track* t = new G4Track(secParticles    
1074         t->SetTouchableHandle(track.GetToucha    
1075         if (biasManager) {                       
1076           t->SetWeight(weight * biasManager->    
1077         } else {                                 
1078           t->SetWeight(weight);                  
1079         }                                        
1080         if(i < num0) {                           
1081           t->SetCreatorModelID(secID);           
1082         } else if(i < num0 + n1) {               
1083           t->SetCreatorModelID(tripletID);       
1084         } else {                                 
1085           t->SetCreatorModelID(biasID);          
1086         }                                        
1087                                                  
1088         //G4cout << "Secondary(post step) has    
1089         //       << ", kenergy " << t->GetKin    
1090         //       << " time= " << time/ns << "    
1091         pParticleChange->AddSecondary(t);        
1092       }                                          
1093     }                                            
1094   }                                              
1095                                                  
1096   if(0.0 == fParticleChange.GetProposedKineti    
1097      fAlive == fParticleChange.GetTrackStatus    
1098     if(particle->GetProcessManager()->GetAtRe    
1099          { fParticleChange.ProposeTrackStatus    
1100     else { fParticleChange.ProposeTrackStatus    
1101   }                                              
1102                                                  
1103   /*                                             
1104   if(-1 < verboseLevel) {                        
1105     G4cout << "::PostStepDoIt: Sample seconda    
1106     << fParticleChange.GetProposedKineticEner    
1107            << " MeV; model= (" << currentMode    
1108            << ", " <<  currentModel->HighEner    
1109            << "  preStepLambda= " << preStepL    
1110            << "  dir= " << track.GetMomentumD    
1111            << "  status= " << track.GetTrackS    
1112            << G4endl;                            
1113   }                                              
1114   */                                             
1115   return &fParticleChange;                       
1116 }                                                
1117                                                  
1118 //....oooOO0OOooo........oooOO0OOooo........o    
1119                                                  
1120 G4bool G4VEnergyLossProcess::StorePhysicsTabl    
1121        const G4ParticleDefinition* part, cons    
1122 {                                                
1123   if (!isMaster || nullptr != baseParticle ||    
1124   for(std::size_t i=0; i<7; ++i) {               
1125     // ionisation table only for ionisation p    
1126     if (nullptr == theData->Table(i) || (!isI    
1127       continue;                                  
1128     }                                            
1129     if (-1 < verboseLevel) {                     
1130       G4cout << "G4VEnergyLossProcess::StoreP    
1131        << "  " << particle->GetParticleName()    
1132        << "  " << GetProcessName()               
1133        << "  " << tnames[i] << "  " << theDat    
1134     }                                            
1135     if (!G4EmTableUtil::StoreTable(this, part    
1136            dir, tnames[i], verboseLevel, asci    
1137       return false;                              
1138     }                                            
1139   }                                              
1140   return true;                                   
1141 }                                                
1142                                                  
1143 //....oooOO0OOooo........oooOO0OOooo........o    
1144                                                  
1145 G4bool                                           
1146 G4VEnergyLossProcess::RetrievePhysicsTable(co    
1147                                            co    
1148 {                                                
1149   if (!isMaster || nullptr != baseParticle ||    
1150   for(std::size_t i=0; i<7; ++i) {               
1151     // ionisation table only for ionisation p    
1152     if (!isIonisation && 1 == i) { continue;     
1153     if(!G4EmTableUtil::RetrieveTable(this, pa    
1154                                      verboseL    
1155       return false;                              
1156     }                                            
1157   }                                              
1158   return true;                                   
1159 }                                                
1160                                                  
1161 //....oooOO0OOooo........oooOO0OOooo........o    
1162                                                  
1163 G4double G4VEnergyLossProcess::GetDEDXDispers    
1164                                   const G4Mat    
1165                                   const G4Dyn    
1166                                         G4dou    
1167 {                                                
1168   DefineMaterial(couple);                        
1169   G4double ekin = dp->GetKineticEnergy();        
1170   SelectModel(ekin*massRatio);                   
1171   G4double tmax = currentModel->MaxSecondaryK    
1172   G4double tcut = std::min(tmax,(*theCuts)[cu    
1173   G4double d = 0.0;                              
1174   G4VEmFluctuationModel* fm = currentModel->G    
1175   if(nullptr != fm) { d = fm->Dispersion(curr    
1176   return d;                                      
1177 }                                                
1178                                                  
1179 //....oooOO0OOooo........oooOO0OOooo........o    
1180                                                  
1181 G4double                                         
1182 G4VEnergyLossProcess::CrossSectionPerVolume(G    
1183                                             c    
1184                                             G    
1185 {                                                
1186   // Cross section per volume is calculated      
1187   DefineMaterial(couple);                        
1188   G4double cross = 0.0;                          
1189   if (nullptr != theLambdaTable) {               
1190     cross = GetLambdaForScaledEnergy(kineticE    
1191                                      logKinet    
1192   } else {                                       
1193     SelectModel(kineticEnergy*massRatio);        
1194     cross = (!baseMat) ? biasFactor : biasFac    
1195     cross *= (currentModel->CrossSectionPerVo    
1196                                                  
1197   }                                              
1198   return std::max(cross, 0.0);                   
1199 }                                                
1200                                                  
1201 //....oooOO0OOooo........oooOO0OOooo........o    
1202                                                  
1203 G4double G4VEnergyLossProcess::MeanFreePath(c    
1204 {                                                
1205   DefineMaterial(track.GetMaterialCutsCouple(    
1206   const G4double kinEnergy    = track.GetKine    
1207   const G4double logKinEnergy = track.GetDyna    
1208   const G4double cs = GetLambdaForScaledEnerg    
1209                                                  
1210   return (0.0 < cs) ? 1.0/cs : DBL_MAX;          
1211 }                                                
1212                                                  
1213 //....oooOO0OOooo........oooOO0OOooo........o    
1214                                                  
1215 G4double G4VEnergyLossProcess::ContinuousStep    
1216                                                  
1217                                                  
1218 {                                                
1219   return AlongStepGetPhysicalInteractionLengt    
1220 }                                                
1221                                                  
1222 //....oooOO0OOooo........oooOO0OOooo........o    
1223                                                  
1224 G4double G4VEnergyLossProcess::GetMeanFreePat    
1225                              const G4Track& t    
1226                              G4double,           
1227                              G4ForceCondition    
1228                                                  
1229 {                                                
1230   *condition = NotForced;                        
1231   return MeanFreePath(track);                    
1232 }                                                
1233                                                  
1234 //....oooOO0OOooo........oooOO0OOooo........o    
1235                                                  
1236 G4double G4VEnergyLossProcess::GetContinuousS    
1237                 const G4Track&,                  
1238                 G4double, G4double, G4double&    
1239 {                                                
1240   return DBL_MAX;                                
1241 }                                                
1242                                                  
1243 //....oooOO0OOooo........oooOO0OOooo........o    
1244                                                  
1245 G4PhysicsVector*                                 
1246 G4VEnergyLossProcess::LambdaPhysicsVector(con    
1247                                           G4d    
1248 {                                                
1249   DefineMaterial(couple);                        
1250   G4PhysicsVector* v = (*theLambdaTable)[base    
1251   return new G4PhysicsVector(*v);                
1252 }                                                
1253                                                  
1254 //....oooOO0OOooo........oooOO0OOooo........o    
1255                                                  
1256 void                                             
1257 G4VEnergyLossProcess::SetDEDXTable(G4PhysicsT    
1258 {                                                
1259   if(1 < verboseLevel) {                         
1260     G4cout << "### Set DEDX table " << p << "    
1261      << "  " <<  theDEDXunRestrictedTable <<     
1262            << " for " << particle->GetParticl    
1263            << " and process " << GetProcessNa    
1264      << " type=" << tType << " isIonisation:"    
1265   }                                              
1266   if(fTotal == tType) {                          
1267     theDEDXunRestrictedTable = p;                
1268   } else if(fRestricted == tType) {              
1269     theDEDXTable = p;                            
1270     if(isMaster && nullptr == baseParticle) {    
1271       theData->UpdateTable(theDEDXTable, 0);     
1272     }                                            
1273   } else if(fIsIonisation == tType) {            
1274     theIonisationTable = p;                      
1275     if(isMaster && nullptr == baseParticle) {    
1276       theData->UpdateTable(theIonisationTable    
1277     }                                            
1278   }                                              
1279 }                                                
1280                                                  
1281 //....oooOO0OOooo........oooOO0OOooo........o    
1282                                                  
1283 void G4VEnergyLossProcess::SetCSDARangeTable(    
1284 {                                                
1285   theCSDARangeTable = p;                         
1286 }                                                
1287                                                  
1288 //....oooOO0OOooo........oooOO0OOooo........o    
1289                                                  
1290 void G4VEnergyLossProcess::SetRangeTableForLo    
1291 {                                                
1292   theRangeTableForLoss = p;                      
1293 }                                                
1294                                                  
1295 //....oooOO0OOooo........oooOO0OOooo........o    
1296                                                  
1297 void G4VEnergyLossProcess::SetInverseRangeTab    
1298 {                                                
1299   theInverseRangeTable = p;                      
1300 }                                                
1301                                                  
1302 //....oooOO0OOooo........oooOO0OOooo........o    
1303                                                  
1304 void G4VEnergyLossProcess::SetLambdaTable(G4P    
1305 {                                                
1306   if(1 < verboseLevel) {                         
1307     G4cout << "### Set Lambda table " << p <<    
1308            << " for " << particle->GetParticl    
1309            << " and process " << GetProcessNa    
1310   }                                              
1311   theLambdaTable = p;                            
1312   tablesAreBuilt = true;                         
1313                                                  
1314   if(isMaster && nullptr != p) {                 
1315     delete theEnergyOfCrossSectionMax;           
1316     theEnergyOfCrossSectionMax = nullptr;        
1317     if(fEmTwoPeaks == fXSType) {                 
1318       if(nullptr != fXSpeaks) {                  
1319   for(auto & ptr : *fXSpeaks) { delete ptr; }    
1320   delete fXSpeaks;                               
1321       }                                          
1322       G4LossTableBuilder* bld = lManager->Get    
1323       fXSpeaks = G4EmUtility::FillPeaksStruct    
1324       if(nullptr == fXSpeaks) { fXSType = fEm    
1325     }                                            
1326     if(fXSType == fEmOnePeak) {                  
1327       theEnergyOfCrossSectionMax = G4EmUtilit    
1328       if(nullptr == theEnergyOfCrossSectionMa    
1329     }                                            
1330   }                                              
1331 }                                                
1332                                                  
1333 //....oooOO0OOooo........oooOO0OOooo........o    
1334                                                  
1335 void G4VEnergyLossProcess::SetEnergyOfCrossSe    
1336 {                                                
1337   theEnergyOfCrossSectionMax = p;                
1338 }                                                
1339                                                  
1340 //....oooOO0OOooo........oooOO0OOooo........o    
1341                                                  
1342 void G4VEnergyLossProcess::SetTwoPeaksXS(std:    
1343 {                                                
1344   fXSpeaks = ptr;                                
1345 }                                                
1346                                                  
1347 //....oooOO0OOooo........oooOO0OOooo........o    
1348                                                  
1349 const G4Element* G4VEnergyLossProcess::GetCur    
1350 {                                                
1351   return (nullptr != currentModel)               
1352     ? currentModel->GetCurrentElement(current    
1353 }                                                
1354                                                  
1355 //....oooOO0OOooo........oooOO0OOooo........o    
1356                                                  
1357 void G4VEnergyLossProcess::SetCrossSectionBia    
1358                                                  
1359 {                                                
1360   if(f > 0.0) {                                  
1361     biasFactor = f;                              
1362     weightFlag = flag;                           
1363     if(1 < verboseLevel) {                       
1364       G4cout << "### SetCrossSectionBiasingFa    
1365              << " process " << GetProcessName    
1366              << " biasFactor= " << f << " wei    
1367              << G4endl;                          
1368     }                                            
1369   }                                              
1370 }                                                
1371                                                  
1372 //....oooOO0OOooo........oooOO0OOooo........o    
1373                                                  
1374 void G4VEnergyLossProcess::ActivateForcedInte    
1375                                                  
1376                                                  
1377 {                                                
1378   if(nullptr == biasManager) { biasManager =     
1379   if(1 < verboseLevel) {                         
1380     G4cout << "### ActivateForcedInteraction:    
1381            << " process " << GetProcessName()    
1382            << " length(mm)= " << length/mm       
1383            << " in G4Region <" << region         
1384            << "> weightFlag= " << flag           
1385            << G4endl;                            
1386   }                                              
1387   weightFlag = flag;                             
1388   biasManager->ActivateForcedInteraction(leng    
1389 }                                                
1390                                                  
1391 //....oooOO0OOooo........oooOO0OOooo........o    
1392                                                  
1393 void                                             
1394 G4VEnergyLossProcess::ActivateSecondaryBiasin    
1395                                                  
1396                                                  
1397 {                                                
1398   if (0.0 <= factor) {                           
1399     // Range cut can be applied only for e-      
1400     if(0.0 == factor && secondaryParticle !=     
1401       { return; }                                
1402                                                  
1403     if(nullptr == biasManager) { biasManager     
1404     biasManager->ActivateSecondaryBiasing(reg    
1405     if(1 < verboseLevel) {                       
1406       G4cout << "### ActivateSecondaryBiasing    
1407              << " process " << GetProcessName    
1408              << " factor= " << factor            
1409              << " in G4Region <" << region       
1410              << "> energyLimit(MeV)= " << ene    
1411              << G4endl;                          
1412     }                                            
1413   }                                              
1414 }                                                
1415                                                  
1416 //....oooOO0OOooo........oooOO0OOooo........o    
1417                                                  
1418 void G4VEnergyLossProcess::SetIonisation(G4bo    
1419 {                                                
1420   isIonisation = val;                            
1421   aGPILSelection = (val) ? CandidateForSelect    
1422 }                                                
1423                                                  
1424 //....oooOO0OOooo........oooOO0OOooo........o    
1425                                                  
1426  void G4VEnergyLossProcess::SetLinearLossLimi    
1427 {                                                
1428   if(0.0 < val && val < 1.0) {                   
1429     linLossLimit = val;                          
1430     actLinLossLimit = true;                      
1431   } else { PrintWarning("SetLinearLossLimit",    
1432 }                                                
1433                                                  
1434 //....oooOO0OOooo........oooOO0OOooo........o    
1435                                                  
1436 void G4VEnergyLossProcess::SetStepFunction(G4    
1437 {                                                
1438   if(0.0 < v1 && 0.0 < v2) {                     
1439     dRoverRange = std::min(1.0, v1);             
1440     finalRange = std::min(v2, 1.e+50);           
1441   } else {                                       
1442     PrintWarning("SetStepFunctionV1", v1);       
1443     PrintWarning("SetStepFunctionV2", v2);       
1444   }                                              
1445 }                                                
1446                                                  
1447 //....oooOO0OOooo........oooOO0OOooo........o    
1448                                                  
1449 void G4VEnergyLossProcess::SetLowestEnergyLim    
1450 {                                                
1451   if(1.e-18 < val && val < 1.e+50) { lowestKi    
1452   else { PrintWarning("SetLowestEnergyLimit",    
1453 }                                                
1454                                                  
1455 //....oooOO0OOooo........oooOO0OOooo........o    
1456                                                  
1457 void G4VEnergyLossProcess::SetDEDXBinning(G4i    
1458 {                                                
1459   if(2 < n && n < 1000000000) {                  
1460     nBins = n;                                   
1461     actBinning = true;                           
1462   } else {                                       
1463     G4double e = (G4double)n;                    
1464     PrintWarning("SetDEDXBinning", e);           
1465   }                                              
1466 }                                                
1467                                                  
1468 //....oooOO0OOooo........oooOO0OOooo........o    
1469                                                  
1470 void G4VEnergyLossProcess::SetMinKinEnergy(G4    
1471 {                                                
1472   if(1.e-18 < e && e < maxKinEnergy) {           
1473     minKinEnergy = e;                            
1474     actMinKinEnergy = true;                      
1475   } else { PrintWarning("SetMinKinEnergy", e)    
1476 }                                                
1477                                                  
1478 //....oooOO0OOooo........oooOO0OOooo........o    
1479                                                  
1480 void G4VEnergyLossProcess::SetMaxKinEnergy(G4    
1481 {                                                
1482   if(minKinEnergy < e && e < 1.e+50) {           
1483     maxKinEnergy = e;                            
1484     actMaxKinEnergy = true;                      
1485     if(e < maxKinEnergyCSDA) { maxKinEnergyCS    
1486   } else { PrintWarning("SetMaxKinEnergy", e)    
1487 }                                                
1488                                                  
1489 //....oooOO0OOooo........oooOO0OOooo........o    
1490                                                  
1491 void G4VEnergyLossProcess::PrintWarning(const    
1492 {                                                
1493   G4String ss = "G4VEnergyLossProcess::" + ti    
1494   G4ExceptionDescription ed;                     
1495   ed << "Parameter is out of range: " << val     
1496      << " it will have no effect!\n" << "  Pr    
1497      << GetProcessName() << "  nbins= " << nB    
1498      << " Emin(keV)= " << minKinEnergy/keV       
1499      << " Emax(GeV)= " << maxKinEnergy/GeV;      
1500   G4Exception(ss, "em0044", JustWarning, ed);    
1501 }                                                
1502                                                  
1503 //....oooOO0OOooo........oooOO0OOooo........o    
1504                                                  
1505 void G4VEnergyLossProcess::ProcessDescription    
1506 {                                                
1507   if(nullptr != particle) { StreamInfo(out, *    
1508 }                                                
1509                                                  
1510 //....oooOO0OOooo........oooOO0OOooo........o    
1511