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 ]

  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 "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........oooOO0OOooo........oooOO0OOooo....
 85 
 86 namespace 
 87 {
 88   G4String tnames[7] =
 89     {"DEDX","Ionisation","DEDXnr","CSDARange","Lambda","Range","InverseRange"};
 90 }
 91 
 92 
 93 G4VEnergyLossProcess::G4VEnergyLossProcess(const G4String& name, 
 94                                            G4ProcessType type): 
 95   G4VContinuousDiscreteProcess(name, type)
 96 {
 97   theParameters = G4EmParameters::Instance();
 98   SetVerboseLevel(1);
 99 
100   // low energy limit
101   lowestKinEnergy = theParameters->LowestElectronEnergy();
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(true);
118   modelManager = new G4EmModelManager();
119   safetyHelper = G4TransportationManager::GetTransportationManager()
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->GetTableBuilder();
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........oooOO0OOooo........oooOO0OOooo....
138 
139 G4VEnergyLossProcess::~G4VEnergyLossProcess()
140 {
141   if (isMaster) {
142     if(nullptr == baseParticle) { delete theData; }
143     delete theEnergyOfCrossSectionMax;
144     if(nullptr != fXSpeaks) {
145       for(auto const & v : *fXSpeaks) { delete v; }
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........oooOO0OOooo........oooOO0OOooo....
157 
158 G4double G4VEnergyLossProcess::MinPrimaryEnergy(const G4ParticleDefinition*, 
159                                                 const G4Material*, 
160                                                 G4double cut)
161 {
162   return cut;
163 }
164 
165 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
166 
167 void G4VEnergyLossProcess::AddEmModel(G4int order, G4VEmModel* ptr,
168                                       G4VEmFluctuationModel* fluc,
169                                       const G4Region* region)
170 {
171   if(nullptr == ptr) { return; }
172   G4VEmFluctuationModel* afluc = (nullptr == fluc) ? fluctModel : fluc;
173   modelManager->AddEmModel(order, ptr, afluc, region);
174   ptr->SetParticleChange(pParticleChange, afluc);
175 }
176 
177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
178 
179 void G4VEnergyLossProcess::SetEmModel(G4VEmModel* ptr, G4int)
180 {
181   if(nullptr == ptr) { return; }
182   if(!emModels->empty()) {
183     for(auto & em : *emModels) { if(em == ptr) { return; } }
184   }
185   emModels->push_back(ptr);  
186 }
187 
188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
189 
190 void G4VEnergyLossProcess::SetDynamicMassCharge(G4double massratio,
191                                                 G4double charge2ratio)
192 {
193   massRatio = massratio;
194   logMassRatio = G4Log(massRatio);
195   fFactor = charge2ratio*biasFactor;
196   if(baseMat) { fFactor *= (*theDensityFactor)[currentCoupleIndex]; }
197   chargeSqRatio = charge2ratio;
198   reduceFactor  = 1.0/(fFactor*massRatio);
199 }
200 
201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
202 
203 void 
204 G4VEnergyLossProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
205 {
206   particle = G4EmTableUtil::CheckIon(this, &part, particle,
207                                      verboseLevel, isIon);
208 
209   if( particle != &part ) {
210     if(!isIon) { lManager->RegisterExtraParticle(&part, this); }
211     if(1 < verboseLevel) {
212       G4cout << "### G4VEnergyLossProcess::PreparePhysicsTable()"
213              << " interrupted for " << GetProcessName() << " and "
214              << part.GetParticleName() << " isIon=" << isIon 
215              << " spline=" << spline << G4endl;
216     }
217     return;
218   }
219 
220   tablesAreBuilt = false;
221   if (GetProcessSubType() == fIonisation) { SetIonisation(true); }
222 
223   G4LossTableBuilder* bld = lManager->GetTableBuilder();
224   lManager->PreparePhysicsTable(&part, this);
225 
226   // Base particle and set of models can be defined here
227   InitialiseEnergyLossProcess(particle, baseParticle);
228 
229   // parameters of the process
230   if(!actLossFluc) { lossFluctuationFlag = theParameters->LossFluctuation(); }
231   useCutAsFinalRange = theParameters->UseCutAsFinalRange();
232   if(!actMinKinEnergy) { minKinEnergy = theParameters->MinKinEnergy(); }
233   if(!actMaxKinEnergy) { maxKinEnergy = theParameters->MaxKinEnergy(); }
234   if(!actBinning) { nBins = theParameters->NumberOfBins(); }
235   maxKinEnergyCSDA = theParameters->MaxEnergyForCSDARange();
236   nBinsCSDA = theParameters->NumberOfBinsPerDecade()
237     *G4lrint(std::log10(maxKinEnergyCSDA/minKinEnergy));
238   if(!actLinLossLimit) { linLossLimit = theParameters->LinearLossLimit(); }
239   lambdaFactor = theParameters->LambdaFactor();
240   invLambdaFactor = 1.0/lambdaFactor;
241   if(isMaster) { SetVerboseLevel(theParameters->Verbose()); }
242   else { SetVerboseLevel(theParameters->WorkerVerbose()); }
243   // integral option may be disabled
244   if(!theParameters->Integral()) { fXSType = fEmNoIntegral; }
245 
246   theParameters->DefineRegParamForLoss(this);
247 
248   fRangeEnergy = 0.0;
249 
250   G4double initialCharge = particle->GetPDGCharge();
251   G4double initialMass   = particle->GetPDGMass();
252 
253   theParameters->FillStepFunction(particle, this);
254 
255   // parameters for scaling from the base particle
256   if (nullptr != baseParticle) {
257     massRatio    = (baseParticle->GetPDGMass())/initialMass;
258     logMassRatio = G4Log(massRatio);
259     G4double q = initialCharge/baseParticle->GetPDGCharge();
260     chargeSqRatio = q*q;
261     if(chargeSqRatio > 0.0) { reduceFactor = 1.0/(chargeSqRatio*massRatio); }
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 G4EmDataHandler(7); }
270 
271     if(nullptr != theDEDXTable && isIonisation) {
272       if(nullptr != theIonisationTable && theDEDXTable != theIonisationTable) {
273   theData->CleanTable(0);
274   theDEDXTable = theIonisationTable;
275   theIonisationTable = nullptr;
276       }
277     }
278     
279     theDEDXTable = theData->MakeTable(theDEDXTable, 0);
280     bld->InitialiseBaseMaterials(theDEDXTable);
281     theData->UpdateTable(theIonisationTable, 1);
282 
283     if (theParameters->BuildCSDARange()) {
284       theDEDXunRestrictedTable = theData->MakeTable(2);
285       if(isIonisation) { theCSDARangeTable = theData->MakeTable(3); }
286     }
287 
288     theLambdaTable = theData->MakeTable(4);
289     if(isIonisation) {
290       theRangeTableForLoss = theData->MakeTable(5);
291       theInverseRangeTable = theData->MakeTable(6);
292     }
293   }
294 
295   // forced biasing
296   if(nullptr != biasManager) { 
297     biasManager->Initialise(part,GetProcessName(),verboseLevel); 
298     biasFlag = false; 
299   }
300   baseMat = bld->GetBaseMaterialFlag();
301   numberOfModels = modelManager->NumberOfModels();
302   currentModel = modelManager->GetModel(0);
303   G4EmTableUtil::UpdateModels(this, modelManager, maxKinEnergy,
304                               numberOfModels, secID, biasID,
305                               mainSecondaries, baseMat, isMaster,
306                               theParameters->UseAngularGeneratorForIonisation());
307   theCuts = modelManager->Initialise(particle, secondaryParticle,
308                                      verboseLevel);
309   // subcut processor
310   if(isIonisation) { 
311     subcutProducer = lManager->SubCutProducer();
312   }
313   if(1 == nSCoffRegions) {
314     if((*scoffRegions)[0]->GetName() == "DefaultRegionForTheWorld") {
315       delete scoffRegions;
316       scoffRegions = nullptr;
317       nSCoffRegions = 0;
318     }
319   }
320 
321   if(1 < verboseLevel) {
322     G4cout << "G4VEnergyLossProcess::PrepearPhysicsTable() is done "
323            << " for " << GetProcessName() << " and " << particle->GetParticleName()
324            << " isIon= " << isIon << " spline=" << spline;
325     if(baseParticle) { 
326       G4cout << "; base: " << baseParticle->GetParticleName(); 
327     }
328     G4cout << G4endl;
329     G4cout << " chargeSqRatio= " << chargeSqRatio
330            << " massRatio= " << massRatio
331            << " reduceFactor= " << reduceFactor << G4endl;
332     if (nSCoffRegions > 0) {
333       G4cout << " SubCut secondary production is ON for regions: " << G4endl;
334       for (G4int i=0; i<nSCoffRegions; ++i) {
335         const G4Region* r = (*scoffRegions)[i];
336         G4cout << "           " << r->GetName() << G4endl;
337       }
338     } else if(nullptr != subcutProducer) {
339       G4cout << " SubCut secondary production is ON for all regions" << G4endl;
340     }
341   }
342 }
343 
344 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
345 
346 void G4VEnergyLossProcess::BuildPhysicsTable(const G4ParticleDefinition& part)
347 {
348   if(1 < verboseLevel) {  
349     G4cout << "### G4VEnergyLossProcess::BuildPhysicsTable() for "
350            << GetProcessName()
351            << " and particle " << part.GetParticleName()
352            << "; the first particle " << particle->GetParticleName();
353     if(baseParticle) { 
354       G4cout << "; base: " << baseParticle->GetParticleName(); 
355     }
356     G4cout << G4endl;
357     G4cout << "    TablesAreBuilt= " << tablesAreBuilt << " isIon= " << isIon
358            << " spline=" << spline << " ptr: " << this << G4endl;
359   }
360 
361   if(&part == particle) {
362     if(isMaster) {
363       lManager->BuildPhysicsTable(particle, this);
364 
365     } else {
366       const auto masterProcess =
367         static_cast<const G4VEnergyLossProcess*>(GetMasterProcess());
368 
369       numberOfModels = modelManager->NumberOfModels();
370       G4EmTableUtil::BuildLocalElossProcess(this, masterProcess,
371                                             particle, numberOfModels);
372       tablesAreBuilt = true;  
373       baseMat = masterProcess->UseBaseMaterial();
374       lManager->LocalPhysicsTables(particle, this);
375     }
376    
377     // needs to be done only once
378     safetyHelper->InitialiseHelper();
379   }
380   // Added tracking cut to avoid tracking artifacts
381   // and identified deexcitation flag
382   if(isIonisation) { 
383     atomDeexcitation = lManager->AtomDeexcitation();
384     if(nullptr != atomDeexcitation) { 
385       if(atomDeexcitation->IsPIXEActive()) { useDeexcitation = true; } 
386     }
387   }
388 
389   // protection against double printout
390   if(theParameters->IsPrintLocked()) { return; }
391 
392   // explicitly defined printout by particle name
393   G4String num = part.GetParticleName();
394   if(1 < verboseLevel || 
395      (0 < verboseLevel && (num == "e-" || 
396                            num == "e+"    || num == "mu+" || 
397                            num == "mu-"   || num == "proton"|| 
398                            num == "pi+"   || num == "pi-" || 
399                            num == "kaon+" || num == "kaon-" || 
400                            num == "alpha" || num == "anti_proton" || 
401                            num == "GenericIon"|| num == "alpha+" ))) { 
402     StreamInfo(G4cout, part); 
403   }
404   if(1 < verboseLevel) {
405     G4cout << "### G4VEnergyLossProcess::BuildPhysicsTable() done for "
406            << GetProcessName()
407            << " and particle " << part.GetParticleName();
408     if(isIonisation) { G4cout << "  isIonisation flag=1"; }
409     G4cout << " baseMat=" << baseMat << G4endl;
410   }
411 }
412 
413 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
414 
415 G4PhysicsTable* G4VEnergyLossProcess::BuildDEDXTable(G4EmTableType tType)
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::BuildDEDXTable WARNING: wrong type "
429            << tType << G4endl;
430   }
431   if(1 < verboseLevel) {
432     G4cout << "G4VEnergyLossProcess::BuildDEDXTable() of type " << tType
433            << " for " << GetProcessName()
434            << " and " << particle->GetParticleName() 
435      << "spline=" << spline << G4endl;
436   }
437   if(nullptr == table) { return table; }
438 
439   G4LossTableBuilder* bld = lManager->GetTableBuilder();
440   G4EmTableUtil::BuildDEDXTable(this, particle, modelManager, bld,
441                                 table, minKinEnergy, emax, bin,
442                                 verboseLevel, tType, spline);
443   return table;
444 }
445 
446 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
447 
448 G4PhysicsTable* G4VEnergyLossProcess::BuildLambdaTable(G4EmTableType)
449 {
450   if(nullptr == theLambdaTable) { return theLambdaTable; }
451 
452   G4double scale = theParameters->MaxKinEnergy()/theParameters->MinKinEnergy();
453   G4int nbin = 
454     theParameters->NumberOfBinsPerDecade()*G4lrint(std::log10(scale));
455   scale = nbin/G4Log(scale);
456   
457   G4LossTableBuilder* bld = lManager->GetTableBuilder();
458   G4EmTableUtil::BuildLambdaTable(this, particle, modelManager,
459                                   bld, theLambdaTable, theCuts,
460                                   minKinEnergy, maxKinEnergy, scale,
461                                   verboseLevel, spline);
462   return theLambdaTable;
463 }
464 
465 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
466 
467 void G4VEnergyLossProcess::StreamInfo(std::ostream& out,
468                 const G4ParticleDefinition& part, G4bool rst) const
469 {
470   G4String indent = (rst ? "  " : "");
471   out << std::setprecision(6);
472   out << G4endl << indent << GetProcessName()  << ": ";
473   if (!rst) out << " for " << part.GetParticleName();
474   out << "  XStype:" << fXSType 
475       << "  SubType=" << GetProcessSubType() << G4endl
476       << "      dE/dx and range tables from "
477       << G4BestUnit(minKinEnergy,"Energy")
478       << " to " << G4BestUnit(maxKinEnergy,"Energy")
479       << " in " << nBins << " bins" << G4endl
480       << "      Lambda tables from threshold to "
481       << G4BestUnit(maxKinEnergy,"Energy")
482       << ", " << theParameters->NumberOfBinsPerDecade() 
483       << " bins/decade, spline: " << spline
484       << G4endl;
485   if(nullptr != theRangeTableForLoss && isIonisation) {
486     out << "      StepFunction=(" << dRoverRange << ", "
487         << finalRange/mm << " mm)"
488         << ", integ: " << fXSType
489         << ", fluct: " << lossFluctuationFlag
490         << ", linLossLim= " << linLossLimit
491         << G4endl;
492   }
493   StreamProcessInfo(out);
494   modelManager->DumpModelList(out, verboseLevel);
495   if(nullptr != theCSDARangeTable && isIonisation) {
496     out << "      CSDA range table up"
497         << " to " << G4BestUnit(maxKinEnergyCSDA,"Energy")
498         << " in " << nBinsCSDA << " bins" << G4endl;
499   }
500   if(nSCoffRegions>0 && isIonisation) {
501     out << "      Subcutoff sampling in " << nSCoffRegions 
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] << " address: " << ta << G4endl; 
508       if(nullptr != ta) { out << *ta << G4endl; }
509     }
510   }
511 }
512 
513 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
514 
515 void G4VEnergyLossProcess::ActivateSubCutoff(const G4Region* r)
516 {
517   if(nullptr == scoffRegions) {
518     scoffRegions = new std::vector<const G4Region*>;
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........oooOO0OOooo........oooOO0OOooo....
532 
533 G4bool G4VEnergyLossProcess::IsRegionForCubcutProcessor(const G4Track& aTrack)
534 {
535   if(0 == nSCoffRegions) { return true; }
536   const G4Region* r = aTrack.GetVolume()->GetLogicalVolume()->GetRegion();
537   for(auto & reg : *scoffRegions) {
538     if(r == reg) { return true; }
539   }
540   return false;
541 }
542 
543 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
544 
545 void G4VEnergyLossProcess::StartTracking(G4Track* track)
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->GetDefinition()->GetPDGMass();
556     massRatio = (nullptr == baseParticle) ? CLHEP::proton_mass_c2/newmass
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........oooOO0OOooo........oooOO0OOooo....
570 
571 G4double G4VEnergyLossProcess::AlongStepGetPhysicalInteractionLength(
572                              const G4Track& track, G4double, G4double, G4double&,
573                              G4GPILSelection* selection)
574 {
575   G4double x = DBL_MAX;
576   *selection = aGPILSelection;
577   if(isIonisation && currentModel->IsActive(preStepScaledEnergy)) {
578     GetScaledRangeForScaledEnergy(preStepScaledEnergy, LogScaledEkin(track));
579     x = (useCutAsFinalRange) ? std::min(finalRange,
580       currentCouple->GetProductionCuts()->GetProductionCut(1)) : finalRange;
581     x = (fRange > x) ? fRange*dRoverRange + x*(1.0 - dRoverRange)*(2.0 - x/fRange)
582       : fRange;
583     /*    
584       G4cout<<"AlongStepGPIL: " << GetProcessName()<<": e="<<preStepKinEnergy
585   << " fRange=" << fRange << " finR=" << finR <<" stepLimit="<<x<<G4endl;
586     */
587   }
588   return x;
589 }
590 
591 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
592 
593 G4double G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength(
594                              const G4Track& track,
595                              G4double   previousStepSize,
596                              G4ForceCondition* condition)
597 {
598   // condition is set to "Not Forced"
599   *condition = NotForced;
600   G4double x = DBL_MAX;
601 
602   // initialisation of material, mass, charge, model 
603   // at the beginning of the step
604   DefineMaterial(track.GetMaterialCutsCouple());
605   preStepKinEnergy       = track.GetKineticEnergy();
606   preStepScaledEnergy    = preStepKinEnergy*massRatio;
607   SelectModel(preStepScaledEnergy);
608 
609   if(!currentModel->IsActive(preStepScaledEnergy)) { 
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 particle on fly
618   if(isIon) {
619     const G4double q2 = currentModel->ChargeSquareRatio(track);
620     fFactor = q2*biasFactor;
621     if(baseMat) { fFactor *= (*theDensityFactor)[currentCoupleIndex]; }
622     reduceFactor = 1.0/(fFactor*massRatio);
623     if (lossFluctuationFlag) {
624       auto fluc = currentModel->GetModelOfFluctuations();
625       fluc->SetParticleAndCharge(track.GetDefinition(), q2);
626     }
627   }
628 
629   // forced biasing only for primary particles
630   if(biasManager) {
631     if(0 == track.GetParentID() && biasFlag && 
632        biasManager->ForcedInteractionRegion((G4int)currentCoupleIndex)) {
633       return biasManager->GetStepLimit((G4int)currentCoupleIndex, previousStepSize);
634     }
635   }
636 
637   ComputeLambdaForScaledEnergy(preStepScaledEnergy, track);
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 DoIt of this process)
649       theNumberOfInteractionLengthLeft = -G4Log( G4UniformRand() );
650       theInitialNumberOfInteractionLength = theNumberOfInteractionLengthLeft; 
651 
652     } else if(currentInteractionLength < DBL_MAX) {
653 
654       // subtract NumberOfInteractionLengthLeft using previous step
655       theNumberOfInteractionLengthLeft -= 
656         previousStepSize/currentInteractionLength;
657 
658       theNumberOfInteractionLengthLeft = 
659         std::max(theNumberOfInteractionLengthLeft, 0.0);
660     }
661 
662     // new mean free path and step limit
663     currentInteractionLength = 1.0/preStepLambda;
664     x = theNumberOfInteractionLengthLeft * currentInteractionLength;
665   }
666 #ifdef G4VERBOSE
667   if (verboseLevel>2) {
668     G4cout << "G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength ";
669     G4cout << "[ " << GetProcessName() << "]" << G4endl; 
670     G4cout << " for " << track.GetDefinition()->GetParticleName() 
671            << " in Material  " <<  currentMaterial->GetName()
672            << " Ekin(MeV)= " << preStepKinEnergy/MeV 
673            << " track material: " << track.GetMaterial()->GetName()
674            <<G4endl;
675     G4cout << "MeanFreePath = " << currentInteractionLength/cm << "[cm]" 
676            << "InteractionLength= " << x/cm <<"[cm] " <<G4endl;
677   }
678 #endif
679   return x;
680 }
681 
682 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
683 
684 void
685 G4VEnergyLossProcess::ComputeLambdaForScaledEnergy(G4double e, const G4Track& track)
686 {
687   // cross section increased with energy
688   if(fXSType == fEmIncreasing) {
689     if(e*invLambdaFactor < mfpKinEnergy) {
690       preStepLambda = GetLambdaForScaledEnergy(e, LogScaledEkin(track));
691       mfpKinEnergy = (preStepLambda > 0.0) ? e : 0.0;
692     }
693 
694     // cross section has one peak
695   } else if(fXSType == fEmOnePeak) {
696     const G4double epeak = (*theEnergyOfCrossSectionMax)[basedCoupleIndex];
697     if(e <= epeak) {
698       if(e*invLambdaFactor < mfpKinEnergy) {
699         preStepLambda = GetLambdaForScaledEnergy(e, LogScaledEkin(track));
700         mfpKinEnergy = (preStepLambda > 0.0) ? e : 0.0;
701       }
702     } else if(e < mfpKinEnergy) { 
703       const G4double e1 = std::max(epeak, e*lambdaFactor);
704       mfpKinEnergy = e1;
705       preStepLambda = GetLambdaForScaledEnergy(e1);
706     }
707 
708     // cross section has more than one peaks
709   } else if(fXSType == fEmTwoPeaks) {
710     G4TwoPeaksXS* xs = (*fXSpeaks)[basedCoupleIndex];
711     const G4double e1peak = xs->e1peak;
712 
713     // below the 1st peak
714     if(e <= e1peak) {
715       if(e*invLambdaFactor < mfpKinEnergy) {
716         preStepLambda = GetLambdaForScaledEnergy(e, LogScaledEkin(track));
717         mfpKinEnergy = (preStepLambda > 0.0) ? e : 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 <= mfpKinEnergy) { 
725         const G4double e1 = std::max(e1peak, e*lambdaFactor);
726         mfpKinEnergy = e1;
727         preStepLambda = GetLambdaForScaledEnergy(e1); 
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 = GetLambdaForScaledEnergy(e, LogScaledEkin(track));
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 <= mfpKinEnergy) { 
744         const G4double e1 = std::max(e2peak, e*lambdaFactor);
745         mfpKinEnergy = e1;
746         preStepLambda = GetLambdaForScaledEnergy(e1); 
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 = GetLambdaForScaledEnergy(e, LogScaledEkin(track));
756       }
757       return;
758     }
759     // above 3d peak
760     if(e <= mfpKinEnergy) { 
761       const G4double e1 = std::max(e3peak, e*lambdaFactor);
762       mfpKinEnergy = e1;
763       preStepLambda = GetLambdaForScaledEnergy(e1); 
764     }
765     // integral method is not used
766   } else {
767     preStepLambda = GetLambdaForScaledEnergy(e, LogScaledEkin(track));
768   }
769 }
770 
771 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
772 
773 G4VParticleChange* G4VEnergyLossProcess::AlongStepDoIt(const G4Track& track,
774                                                        const G4Step& step)
775 {
776   fParticleChange.InitializeForAlongStep(track);
777   // The process has range table - calculate energy loss
778   if(!isIonisation || !currentModel->IsActive(preStepScaledEnergy)) {
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.GetParticleDefinition();
788     G4cout << "AlongStepDoIt for "
789            << GetProcessName() << " and particle " << d->GetParticleName()
790            << "  eScaled(MeV)=" << preStepScaledEnergy/MeV
791            << "  range(mm)=" << fRange/mm << "  s(mm)=" << length/mm
792            << "  rf=" << reduceFactor << "  q^2=" << chargeSqRatio
793            << " md=" << d->GetPDGMass() << "  status=" << track.GetTrackStatus()
794            << "  " << track.GetMaterial()->GetName() << G4endl;
795   }
796   */
797   const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
798 
799   // define new weight for primary and secondaries
800   G4double weight = fParticleChange.GetParentWeight();
801   if(weightFlag) {
802     weight /= biasFactor;
803     fParticleChange.ProposeWeight(weight);
804   }
805 
806   // stopping, check actual range and kinetic energy
807   if (length >= fRange || preStepKinEnergy <= lowestKinEnergy) {
808     eloss = preStepKinEnergy;
809     if (useDeexcitation) {
810       atomDeexcitation->AlongStepDeexcitation(scTracks, step, 
811                                               eloss, (G4int)currentCoupleIndex);
812       if(scTracks.size() > 0) { FillSecondariesAlongStep(weight); }
813       eloss = std::max(eloss, 0.0);
814     }
815     fParticleChange.SetProposedKineticEnergy(0.0);
816     fParticleChange.ProposeLocalEnergyDeposit(eloss);
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(preStepScaledEnergy,
824                                         LogScaledEkin(track));
825   /*
826   G4cout << "##### Short STEP: eloss= " << eloss 
827    << " Escaled=" << preStepScaledEnergy
828    << " R=" << fRange
829    << " L=" << length 
830    << " fFactor=" << fFactor << " minE=" << minKinEnergy 
831    << " idxBase=" << basedCoupleIndex << G4endl;
832   */
833   // Long step
834   if(eloss > preStepKinEnergy*linLossLimit) {
835 
836     const G4double x = (fRange - length)/reduceFactor;
837     const G4double de = preStepKinEnergy - ScaledKinEnergyForLoss(x)/massRatio;
838     if(de > 0.0) { eloss = de; }    
839     /*
840     if(-1 < verboseLevel) 
841       G4cout << "  Long STEP: rPre(mm)=" 
842              << GetScaledRangeForScaledEnergy(preStepScaledEnergy)/mm
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)= " << eloss/MeV
854            << " e-eloss= " << preStepKinEnergy-eloss
855            << " step(mm)= " << length/mm << " range(mm)= " << fRange/mm
856            << " fluct= " << lossFluctuationFlag << G4endl;
857   }
858   */
859 
860   const G4double cut = (*theCuts)[currentCoupleIndex];
861   G4double esec = 0.0;
862 
863   // Corrections, which cannot be tabulated
864   if(isIon) {
865     currentModel->CorrectionsAlongStep(currentCouple, dynParticle, 
866                                        length, eloss);
867     eloss = std::max(eloss, 0.0);
868   }
869 
870   // Sample fluctuations if not full energy loss
871   if(eloss >= preStepKinEnergy) {
872     eloss = preStepKinEnergy;
873 
874   } else if (lossFluctuationFlag) {
875     const G4double tmax = currentModel->MaxSecondaryKinEnergy(dynParticle);
876     const G4double tcut = std::min(cut, tmax);
877     G4VEmFluctuationModel* fluc = currentModel->GetModelOfFluctuations();
878     eloss = fluc->SampleFluctuations(currentCouple,dynParticle,
879                                      tcut, tmax, length, eloss);
880     /*
881     if(-1 < verboseLevel) 
882       G4cout << "After fluct: eloss(MeV)= " << eloss/MeV
883              << " fluc= " << (eloss-eloss0)/MeV
884              << " ChargeSqRatio= " << chargeSqRatio
885              << " massRatio= " << massRatio << " tmax= " << tmax << G4endl;
886     */
887   }
888 
889   // deexcitation
890   if (useDeexcitation) {
891     G4double esecfluo = preStepKinEnergy;
892     G4double de = esecfluo;
893     atomDeexcitation->AlongStepDeexcitation(scTracks, step, 
894                                             de, (G4int)currentCoupleIndex);
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 && IsRegionForCubcutProcessor(track)) {
909     subcutProducer->SampleSecondaries(step, scTracks, eloss, cut);
910   }
911   // secondaries from atomic de-excitation and subcut
912   if(!scTracks.empty()) { FillSecondariesAlongStep(weight); }
913 
914   // Energy balance
915   G4double finalT = preStepKinEnergy - eloss - esec;
916   if (finalT <= lowestKinEnergy) {
917     eloss += finalT;
918     finalT = 0.0;
919   } else if(isIon) {
920     fParticleChange.SetProposedCharge(
921       currentModel->GetParticleCharge(track.GetParticleDefinition(),
922                                       currentMaterial,finalT));
923   }
924   eloss = std::max(eloss, 0.0);
925 
926   fParticleChange.SetProposedKineticEnergy(finalT);
927   fParticleChange.ProposeLocalEnergyDeposit(eloss);
928   /*
929   if(-1 < verboseLevel) {
930     G4double del = finalT + eloss + esec - preStepKinEnergy;
931     G4cout << "Final value eloss(MeV)= " << eloss/MeV
932            << " preStepKinEnergy= " << preStepKinEnergy
933            << " postStepKinEnergy= " << finalT
934            << " de(keV)= " << del/keV
935            << " lossFlag= " << lossFluctuationFlag
936            << "  status= " << track.GetTrackStatus()
937            << G4endl;
938   }
939   */
940   return &fParticleChange;
941 }
942 
943 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
944 
945 void G4VEnergyLossProcess::FillSecondariesAlongStep(G4double wt)
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((G4int)currentCoupleIndex)) {
952       weight *=
953         biasManager->ApplySecondaryBiasing(scTracks, (G4int)currentCoupleIndex);
954     }
955   } 
956 
957   // fill secondaries
958   const std::size_t n = scTracks.size();
959   fParticleChange.SetNumberOfSecondaries((G4int)n);
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()->GetPDGEncoding();
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........oooOO0OOooo........oooOO0OOooo....
984 
985 G4VParticleChange* G4VEnergyLossProcess::PostStepDoIt(const G4Track& track,
986                                                       const G4Step& step)
987 {
988   // clear number of interaction lengths in any case
989   theNumberOfInteractionLengthLeft = -1.0;
990   mfpKinEnergy = DBL_MAX;
991 
992   fParticleChange.InitializeForPostStep(track);
993   const G4double finalT = track.GetKineticEnergy();
994 
995   const G4double postStepScaledEnergy = finalT*massRatio;
996   SelectModel(postStepScaledEnergy);
997 
998   if(!currentModel->IsActive(postStepScaledEnergy)) { 
999     return &fParticleChange; 
1000   }
1001   /*
1002   if(1 < verboseLevel) {
1003     G4cout<<GetProcessName()<<" PostStepDoIt: E(MeV)= "<< finalT/MeV<< G4endl;
1004   }
1005   */
1006   // forced process - should happen only once per track
1007   if(biasFlag) {
1008     if(biasManager->ForcedInteractionRegion((G4int)currentCoupleIndex)) {
1009       biasFlag = false;
1010     }
1011   }
1012   const G4DynamicParticle* dp = track.GetDynamicParticle();
1013 
1014   // Integral approach
1015   if (fXSType != fEmNoIntegral) {
1016     const G4double logFinalT = dp->GetLogKineticEnergy();
1017     G4double lx = GetLambdaForScaledEnergy(postStepScaledEnergy,
1018                                            logFinalT + logMassRatio);
1019     lx = std::max(lx, 0.0);
1020 
1021     // if both lg and lx are zero then no interaction
1022     if(preStepLambda*G4UniformRand() >= lx) {
1023       return &fParticleChange;
1024     }
1025   }
1026 
1027   // define new weight for primary and secondaries
1028   G4double weight = fParticleChange.GetParentWeight();
1029   if(weightFlag) {
1030     weight /= biasFactor;
1031     fParticleChange.ProposeWeight(weight);
1032   }
1033 
1034   const G4double tcut = (*theCuts)[currentCoupleIndex];
1035 
1036   // sample secondaries
1037   secParticles.clear();
1038   currentModel->SampleSecondaries(&secParticles, currentCouple, dp, tcut);
1039 
1040   const G4int num0 = (G4int)secParticles.size();
1041 
1042   // bremsstrahlung splitting or Russian roulette  
1043   if(biasManager) {
1044     if(biasManager->SecondaryBiasingRegion((G4int)currentCoupleIndex)) {
1045       G4double eloss = 0.0;
1046       weight *= biasManager->ApplySecondaryBiasing(
1047                                       secParticles,
1048                                       track, currentModel, 
1049                                       &fParticleChange, eloss,
1050                                       (G4int)currentCoupleIndex, tcut, 
1051                                       step.GetPostStepPoint()->GetSafety());
1052       if(eloss > 0.0) {
1053         eloss += fParticleChange.GetLocalEnergyDeposit();
1054         fParticleChange.ProposeLocalEnergyDeposit(eloss);
1055       }
1056     }
1057   }
1058 
1059   // save secondaries
1060   const G4int num = (G4int)secParticles.size();
1061   if(num > 0) {
1062 
1063     fParticleChange.SetNumberOfSecondaries(num);
1064     G4double time = track.GetGlobalTime();
1065 
1066     G4int n1(0), n2(0);
1067     if(num0 > mainSecondaries) { 
1068       currentModel->FillNumberOfSecondaries(n1, n2);
1069     }
1070 
1071     for (G4int i=0; i<num; ++i) {
1072       if(nullptr != secParticles[i]) {
1073         G4Track* t = new G4Track(secParticles[i], time, track.GetPosition());
1074         t->SetTouchableHandle(track.GetTouchableHandle());
1075         if (biasManager) {
1076           t->SetWeight(weight * biasManager->GetWeight(i));
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 weight " << t->GetWeight() 
1089         //       << ", kenergy " << t->GetKineticEnergy()/MeV << " MeV" 
1090         //       << " time= " << time/ns << " ns " << G4endl;
1091         pParticleChange->AddSecondary(t);
1092       }
1093     }
1094   }
1095 
1096   if(0.0 == fParticleChange.GetProposedKineticEnergy() &&
1097      fAlive == fParticleChange.GetTrackStatus()) {
1098     if(particle->GetProcessManager()->GetAtRestProcessVector()->size() > 0)
1099          { fParticleChange.ProposeTrackStatus(fStopButAlive); }
1100     else { fParticleChange.ProposeTrackStatus(fStopAndKill); }
1101   }
1102 
1103   /*
1104   if(-1 < verboseLevel) {
1105     G4cout << "::PostStepDoIt: Sample secondary; Efin= " 
1106     << fParticleChange.GetProposedKineticEnergy()/MeV
1107            << " MeV; model= (" << currentModel->LowEnergyLimit()
1108            << ", " <<  currentModel->HighEnergyLimit() << ")"
1109            << "  preStepLambda= " << preStepLambda
1110            << "  dir= " << track.GetMomentumDirection()
1111            << "  status= " << track.GetTrackStatus()
1112            << G4endl;
1113   }
1114   */
1115   return &fParticleChange;
1116 }
1117 
1118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1119 
1120 G4bool G4VEnergyLossProcess::StorePhysicsTable(
1121        const G4ParticleDefinition* part, const G4String& dir, G4bool ascii)
1122 {
1123   if (!isMaster || nullptr != baseParticle || part != particle ) return true;
1124   for(std::size_t i=0; i<7; ++i) {
1125     // ionisation table only for ionisation process
1126     if (nullptr == theData->Table(i) || (!isIonisation && 1 == i)) {
1127       continue;
1128     }
1129     if (-1 < verboseLevel) {
1130       G4cout << "G4VEnergyLossProcess::StorePhysicsTable i=" << i
1131        << "  " << particle->GetParticleName()
1132        << "  " << GetProcessName()
1133        << "  " << tnames[i] << "  " << theData->Table(i) << G4endl;
1134     }
1135     if (!G4EmTableUtil::StoreTable(this, part, theData->Table(i),
1136            dir, tnames[i], verboseLevel, ascii)) { 
1137       return false;
1138     }
1139   }
1140   return true;
1141 }
1142 
1143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
1144 
1145 G4bool 
1146 G4VEnergyLossProcess::RetrievePhysicsTable(const G4ParticleDefinition* part, 
1147                                            const G4String& dir, G4bool ascii)
1148 {
1149   if (!isMaster || nullptr != baseParticle || part != particle ) return true;
1150   for(std::size_t i=0; i<7; ++i) {
1151     // ionisation table only for ionisation process
1152     if (!isIonisation && 1 == i) { continue; }
1153     if(!G4EmTableUtil::RetrieveTable(this, part, theData->Table(i), dir, tnames[i],
1154                                      verboseLevel, ascii, spline)) { 
1155       return false;
1156     }
1157   }
1158   return true;
1159 }
1160 
1161 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1162 
1163 G4double G4VEnergyLossProcess::GetDEDXDispersion(
1164                                   const G4MaterialCutsCouple *couple,
1165                                   const G4DynamicParticle* dp,
1166                                         G4double length)
1167 {
1168   DefineMaterial(couple);
1169   G4double ekin = dp->GetKineticEnergy();
1170   SelectModel(ekin*massRatio);
1171   G4double tmax = currentModel->MaxSecondaryKinEnergy(dp);
1172   G4double tcut = std::min(tmax,(*theCuts)[currentCoupleIndex]);
1173   G4double d = 0.0;
1174   G4VEmFluctuationModel* fm = currentModel->GetModelOfFluctuations();
1175   if(nullptr != fm) { d = fm->Dispersion(currentMaterial,dp,tcut,tmax,length); }
1176   return d;
1177 }
1178 
1179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1180 
1181 G4double
1182 G4VEnergyLossProcess::CrossSectionPerVolume(G4double kineticEnergy,
1183                                             const G4MaterialCutsCouple* couple,
1184                                             G4double logKineticEnergy)
1185 {
1186   // Cross section per volume is calculated
1187   DefineMaterial(couple);
1188   G4double cross = 0.0;
1189   if (nullptr != theLambdaTable) {
1190     cross = GetLambdaForScaledEnergy(kineticEnergy * massRatio,
1191                                      logKineticEnergy + logMassRatio);
1192   } else {
1193     SelectModel(kineticEnergy*massRatio);
1194     cross = (!baseMat) ? biasFactor : biasFactor*(*theDensityFactor)[currentCoupleIndex];
1195     cross *= (currentModel->CrossSectionPerVolume(currentMaterial, particle, kineticEnergy,
1196                                                   (*theCuts)[currentCoupleIndex]));
1197   }
1198   return std::max(cross, 0.0);
1199 }
1200 
1201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1202 
1203 G4double G4VEnergyLossProcess::MeanFreePath(const G4Track& track)
1204 {
1205   DefineMaterial(track.GetMaterialCutsCouple());
1206   const G4double kinEnergy    = track.GetKineticEnergy();
1207   const G4double logKinEnergy = track.GetDynamicParticle()->GetLogKineticEnergy();
1208   const G4double cs = GetLambdaForScaledEnergy(kinEnergy * massRatio, 
1209                                                logKinEnergy + logMassRatio);
1210   return (0.0 < cs) ? 1.0/cs : DBL_MAX;
1211 }
1212 
1213 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1214 
1215 G4double G4VEnergyLossProcess::ContinuousStepLimit(const G4Track& track, 
1216                                                    G4double x, G4double y, 
1217                                                    G4double& z)
1218 {
1219   return AlongStepGetPhysicalInteractionLength(track, x, y, z, &aGPILSelection);
1220 }
1221 
1222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1223 
1224 G4double G4VEnergyLossProcess::GetMeanFreePath(
1225                              const G4Track& track,
1226                              G4double,
1227                              G4ForceCondition* condition)
1228 
1229 {
1230   *condition = NotForced;
1231   return MeanFreePath(track);
1232 }
1233 
1234 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1235 
1236 G4double G4VEnergyLossProcess::GetContinuousStepLimit(
1237                 const G4Track&,
1238                 G4double, G4double, G4double&)
1239 {
1240   return DBL_MAX;
1241 }
1242 
1243 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1244 
1245 G4PhysicsVector* 
1246 G4VEnergyLossProcess::LambdaPhysicsVector(const G4MaterialCutsCouple* couple, 
1247                                           G4double)
1248 {
1249   DefineMaterial(couple);
1250   G4PhysicsVector* v = (*theLambdaTable)[basedCoupleIndex];
1251   return new G4PhysicsVector(*v);
1252 }
1253 
1254 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1255 
1256 void 
1257 G4VEnergyLossProcess::SetDEDXTable(G4PhysicsTable* p, G4EmTableType tType)
1258 {
1259   if(1 < verboseLevel) {
1260     G4cout << "### Set DEDX table " << p << "  " << theDEDXTable
1261      << "  " <<  theDEDXunRestrictedTable << "  " << theIonisationTable
1262            << " for " << particle->GetParticleName()
1263            << " and process " << GetProcessName() 
1264      << " type=" << tType << " isIonisation:" << isIonisation << G4endl;
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, 1);
1277     }
1278   }
1279 }
1280 
1281 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1282 
1283 void G4VEnergyLossProcess::SetCSDARangeTable(G4PhysicsTable* p)
1284 {
1285   theCSDARangeTable = p;
1286 }
1287 
1288 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1289 
1290 void G4VEnergyLossProcess::SetRangeTableForLoss(G4PhysicsTable* p)
1291 {
1292   theRangeTableForLoss = p;
1293 }
1294 
1295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1296 
1297 void G4VEnergyLossProcess::SetInverseRangeTable(G4PhysicsTable* p)
1298 {
1299   theInverseRangeTable = p;
1300 }
1301 
1302 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1303 
1304 void G4VEnergyLossProcess::SetLambdaTable(G4PhysicsTable* p)
1305 {
1306   if(1 < verboseLevel) {
1307     G4cout << "### Set Lambda table " << p << " " << theLambdaTable 
1308            << " for " << particle->GetParticleName()
1309            << " and process " << GetProcessName() << G4endl;
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->GetTableBuilder();
1323       fXSpeaks = G4EmUtility::FillPeaksStructure(p, bld);
1324       if(nullptr == fXSpeaks) { fXSType = fEmOnePeak; }
1325     }
1326     if(fXSType == fEmOnePeak) { 
1327       theEnergyOfCrossSectionMax = G4EmUtility::FindCrossSectionMax(p);
1328       if(nullptr == theEnergyOfCrossSectionMax) { fXSType = fEmIncreasing; }
1329     }
1330   }
1331 }
1332 
1333 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1334 
1335 void G4VEnergyLossProcess::SetEnergyOfCrossSectionMax(std::vector<G4double>* p)
1336 {
1337   theEnergyOfCrossSectionMax = p;
1338 }
1339 
1340 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1341 
1342 void G4VEnergyLossProcess::SetTwoPeaksXS(std::vector<G4TwoPeaksXS*>* ptr)
1343 {
1344   fXSpeaks = ptr;
1345 }
1346 
1347 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1348 
1349 const G4Element* G4VEnergyLossProcess::GetCurrentElement() const
1350 {
1351   return (nullptr != currentModel) 
1352     ? currentModel->GetCurrentElement(currentMaterial) : nullptr;
1353 }
1354 
1355 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1356 
1357 void G4VEnergyLossProcess::SetCrossSectionBiasingFactor(G4double f, 
1358                                                         G4bool flag)
1359 {
1360   if(f > 0.0) { 
1361     biasFactor = f; 
1362     weightFlag = flag;
1363     if(1 < verboseLevel) {
1364       G4cout << "### SetCrossSectionBiasingFactor: for " 
1365              << " process " << GetProcessName()
1366              << " biasFactor= " << f << " weightFlag= " << flag 
1367              << G4endl; 
1368     }
1369   }
1370 }
1371 
1372 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1373 
1374 void G4VEnergyLossProcess::ActivateForcedInteraction(G4double length, 
1375                                                      const G4String& region,
1376                                                      G4bool flag)
1377 {
1378   if(nullptr == biasManager) { biasManager = new G4EmBiasingManager(); }
1379   if(1 < verboseLevel) {
1380     G4cout << "### ActivateForcedInteraction: for " 
1381            << " process " << GetProcessName()
1382            << " length(mm)= " << length/mm
1383            << " in G4Region <" << region
1384            << "> weightFlag= " << flag 
1385            << G4endl; 
1386   }
1387   weightFlag = flag;
1388   biasManager->ActivateForcedInteraction(length, region);
1389 }
1390 
1391 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1392 
1393 void 
1394 G4VEnergyLossProcess::ActivateSecondaryBiasing(const G4String& region, 
1395                                                G4double factor, 
1396                                                G4double energyLimit)
1397 {
1398   if (0.0 <= factor) {
1399     // Range cut can be applied only for e-
1400     if(0.0 == factor && secondaryParticle != G4Electron::Electron())
1401       { return; }
1402 
1403     if(nullptr == biasManager) { biasManager = new G4EmBiasingManager(); }
1404     biasManager->ActivateSecondaryBiasing(region, factor, energyLimit);
1405     if(1 < verboseLevel) {
1406       G4cout << "### ActivateSecondaryBiasing: for " 
1407              << " process " << GetProcessName()
1408              << " factor= " << factor
1409              << " in G4Region <" << region 
1410              << "> energyLimit(MeV)= " << energyLimit/MeV
1411              << G4endl; 
1412     }
1413   }
1414 }
1415 
1416 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1417 
1418 void G4VEnergyLossProcess::SetIonisation(G4bool val)
1419 {
1420   isIonisation = val;
1421   aGPILSelection = (val) ? CandidateForSelection : NotCandidateForSelection;
1422 }
1423 
1424 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1425 
1426  void G4VEnergyLossProcess::SetLinearLossLimit(G4double val)
1427 {
1428   if(0.0 < val && val < 1.0) { 
1429     linLossLimit = val;
1430     actLinLossLimit = true; 
1431   } else { PrintWarning("SetLinearLossLimit", val); }
1432 }
1433 
1434 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1435 
1436 void G4VEnergyLossProcess::SetStepFunction(G4double v1, G4double v2)
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........oooOO0OOooo........oooOO0OOooo....
1448 
1449 void G4VEnergyLossProcess::SetLowestEnergyLimit(G4double val)
1450 {
1451   if(1.e-18 < val && val < 1.e+50) { lowestKinEnergy = val; }
1452   else { PrintWarning("SetLowestEnergyLimit", val); }
1453 }
1454 
1455 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1456 
1457 void G4VEnergyLossProcess::SetDEDXBinning(G4int n)
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........oooOO0OOooo........oooOO0OOooo....
1469 
1470 void G4VEnergyLossProcess::SetMinKinEnergy(G4double e)
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........oooOO0OOooo........oooOO0OOooo....
1479 
1480 void G4VEnergyLossProcess::SetMaxKinEnergy(G4double e)
1481 {
1482   if(minKinEnergy < e && e < 1.e+50) { 
1483     maxKinEnergy = e;
1484     actMaxKinEnergy = true;
1485     if(e < maxKinEnergyCSDA) { maxKinEnergyCSDA = e; }
1486   } else { PrintWarning("SetMaxKinEnergy", e); } 
1487 }
1488 
1489 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1490 
1491 void G4VEnergyLossProcess::PrintWarning(const G4String& tit, G4double val) const
1492 {
1493   G4String ss = "G4VEnergyLossProcess::" + tit; 
1494   G4ExceptionDescription ed;
1495   ed << "Parameter is out of range: " << val 
1496      << " it will have no effect!\n" << "  Process " 
1497      << GetProcessName() << "  nbins= " << nBins 
1498      << " Emin(keV)= " << minKinEnergy/keV 
1499      << " Emax(GeV)= " << maxKinEnergy/GeV;
1500   G4Exception(ss, "em0044", JustWarning, ed);
1501 }
1502 
1503 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1504 
1505 void G4VEnergyLossProcess::ProcessDescription(std::ostream& out) const
1506 {
1507   if(nullptr != particle) { StreamInfo(out, *particle, true); }
1508 }
1509 
1510 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1511