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 9.2.p4)


  1 //                                             <<   1 //
  2 // ******************************************* <<   2 // ********************************************************************
  3 // * License and Disclaimer                    <<   3 // * License and Disclaimer                                           *
  4 // *                                           <<   4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th <<   5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided <<   6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License <<   7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/ <<   8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.      <<   9 // * include a list of copyright holders.                             *
 10 // *                                           <<  10 // *                                                                  *
 11 // * Neither the authors of this software syst <<  11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin <<  12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran <<  13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum <<  14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  <<  15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio <<  16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                           <<  17 // *                                                                  *
 18 // * This  code  implementation is the result  <<  18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio <<  19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri <<  20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag <<  21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati <<  22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof <<  23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ******************************************* <<  24 // ********************************************************************
 25 //                                             <<  25 //
 26 // ------------------------------------------- <<  26 // $Id: G4VEnergyLossProcess.cc,v 1.143.2.1 2010/01/26 14:33:54 gcosmo Exp $
 27 //                                             <<  27 // GEANT4 tag $Name: geant4-09-02-patch-04 $
 28 // GEANT4 Class file                           <<  28 //
 29 //                                             <<  29 // -------------------------------------------------------------------
 30 //                                             <<  30 //
 31 // File name:     G4VEnergyLossProcess         <<  31 // GEANT4 Class file
 32 //                                             <<  32 //
 33 // Author:        Vladimir Ivanchenko          <<  33 //
 34 //                                             <<  34 // File name:     G4VEnergyLossProcess
 35 // Creation date: 03.01.2002                   <<  35 //
 36 //                                             <<  36 // Author:        Vladimir Ivanchenko
 37 // Modifications: Vladimir Ivanchenko          <<  37 //
 38 //                                             <<  38 // Creation date: 03.01.2002
 39 //                                             <<  39 //
 40 // Class Description:                          <<  40 // Modifications:
 41 //                                             <<  41 //
 42 // It is the unified energy loss process it ca <<  42 // 13-11-02 Minor fix - use normalised direction (V.Ivanchenko)
 43 // energy loss for charged particles using a s <<  43 // 04-12-02 Minor change in PostStepDoIt (V.Ivanchenko)
 44 // models valid for different energy regions.  <<  44 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
 45 // to create and access to dE/dx and range tab <<  45 // 26-12-02 Secondary production moved to derived classes (V.Ivanchenko)
 46 // that information on fly.                    <<  46 // 04-01-03 Fix problem of very small steps for ions (V.Ivanchenko)
 47 // ------------------------------------------- <<  47 // 20-01-03 Migrade to cut per region (V.Ivanchenko)
 48 //                                             <<  48 // 24-01-03 Temporarily close a control on usage of couples (V.Ivanchenko)
 49 //....oooOO0OOooo........oooOO0OOooo........oo <<  49 // 24-01-03 Make models region aware (V.Ivanchenko)
 50 //....oooOO0OOooo........oooOO0OOooo........oo <<  50 // 05-02-03 Fix compilation warnings (V.Ivanchenko)
 51                                                <<  51 // 06-02-03 Add control on tmax in PostStepDoIt (V.Ivanchenko)
 52 #include "G4VEnergyLossProcess.hh"             <<  52 // 13-02-03 SubCutoffProcessors defined for regions (V.Ivanchenko)
 53 #include "G4PhysicalConstants.hh"              <<  53 // 15-02-03 Lambda table can be scaled (V.Ivanchenko)
 54 #include "G4SystemOfUnits.hh"                  <<  54 // 17-02-03 Fix problem of store/restore tables (V.Ivanchenko)
 55 #include "G4ProcessManager.hh"                 <<  55 // 18-02-03 Add control on CutCouple usage (V.Ivanchenko)
 56 #include "G4LossTableManager.hh"               <<  56 // 26-02-03 Simplify control on GenericIons (V.Ivanchenko)
 57 #include "G4LossTableBuilder.hh"               <<  57 // 06-03-03 Control on GenericIons using SubType + update verbose (V.Ivanchenko)
 58 #include "G4Step.hh"                           <<  58 // 10-03-03 Add Ion registration (V.Ivanchenko)
 59 #include "G4ParticleDefinition.hh"             <<  59 // 22-03-03 Add Initialisation of cash (V.Ivanchenko)
 60 #include "G4ParticleTable.hh"                  <<  60 // 26-03-03 Remove finalRange modification (V.Ivanchenko)
 61 #include "G4EmParameters.hh"                   <<  61 // 09-04-03 Fix problem of negative range limit for non integral (V.Ivanchenko)
 62 #include "G4EmUtility.hh"                      <<  62 // 26-04-03 Fix retrieve tables (V.Ivanchenko)
 63 #include "G4EmTableUtil.hh"                    <<  63 // 06-05-03 Set defalt finalRange = 1 mm (V.Ivanchenko)
 64 #include "G4VEmModel.hh"                       <<  64 // 12-05-03 Update range calculations + lowKinEnergy (V.Ivanchenko)
 65 #include "G4VEmFluctuationModel.hh"            <<  65 // 13-05-03 Add calculation of precise range (V.Ivanchenko)
 66 #include "G4DataVector.hh"                     <<  66 // 23-05-03 Remove tracking cuts (V.Ivanchenko)
 67 #include "G4PhysicsLogVector.hh"               <<  67 // 03-06-03 Fix initialisation problem for STD ionisation (V.Ivanchenko)
 68 #include "G4VParticleChange.hh"                <<  68 // 21-07-03 Add UpdateEmModel method (V.Ivanchenko)
 69 #include "G4Electron.hh"                       <<  69 // 03-11-03 Fix initialisation problem in RetrievePhysicsTable (V.Ivanchenko)
 70 #include "G4ProcessManager.hh"                 <<  70 // 04-11-03 Add checks in RetrievePhysicsTable (V.Ivanchenko)
 71 #include "G4UnitsTable.hh"                     <<  71 // 12-11-03 G4EnergyLossSTD -> G4EnergyLossProcess (V.Ivanchenko)
 72 #include "G4Region.hh"                         <<  72 // 21-01-04 Migrade to G4ParticleChangeForLoss (V.Ivanchenko)
 73 #include "G4RegionStore.hh"                    <<  73 // 27-02-04 Fix problem of loss in low presure gases, cleanup precise range
 74 #include "G4PhysicsTableHelper.hh"             <<  74 //          calculation, use functions ForLoss in AlongStepDoIt (V.Ivanchenko)
 75 #include "G4SafetyHelper.hh"                   <<  75 // 10-03-04 Fix a problem of Precise Range table (V.Ivanchenko)
 76 #include "G4EmDataHandler.hh"                  <<  76 // 19-03-04 Fix a problem energy below lowestKinEnergy (V.Ivanchenko)
 77 #include "G4TransportationManager.hh"          <<  77 // 31-03-04 Fix a problem of retrieve tables (V.Ivanchenko)
 78 #include "G4VAtomDeexcitation.hh"              <<  78 // 21-07-04 Check weather AtRest are active or not (V.Ivanchenko)
 79 #include "G4VSubCutProducer.hh"                <<  79 // 03-08-04 Add pointer of DEDX table to all processes (V.Ivanchenko)
 80 #include "G4EmBiasingManager.hh"               <<  80 // 06-08-04 Clear up names of member functions (V.Ivanchenko)
 81 #include "G4Log.hh"                            <<  81 // 06-08-04 Clear up names of member functions (V.Ivanchenko)
 82 #include <iostream>                            <<  82 // 27-08-04 Add NeedBuildTables method (V.Ivanchneko)
 83                                                <<  83 // 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivantchenko)
 84 //....oooOO0OOooo........oooOO0OOooo........oo <<  84 // 11-03-05 Shift verbose level by 1 (V.Ivantchenko)
 85                                                <<  85 // 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
 86 namespace                                      <<  86 // 11-04-05 Use MaxSecondaryEnergy from a model (V.Ivanchenko)
 87 {                                              <<  87 // 25-07-05 Add extra protection PostStep for non-integral mode (V.Ivanchenko)
 88   G4String tnames[7] =                         <<  88 // 12-08-05 Integral=false; SetStepFunction(0.2, 0.1*mm) (mma)
 89     {"DEDX","Ionisation","DEDXnr","CSDARange", <<  89 // 18-08-05 Return back both AlongStep and PostStep from 7.0 (V.Ivanchenko)
 90 }                                              <<  90 // 02-09-05 Default StepFunction 0.2 1 mm + integral (V.Ivanchenko)
 91                                                <<  91 // 04-09-05 default lambdaFactor 0.8 (V.Ivanchenko)
 92                                                <<  92 // 05-10-05 protection against 0 energy loss added (L.Urban)
 93 G4VEnergyLossProcess::G4VEnergyLossProcess(con <<  93 // 17-10-05 protection above has been removed (L.Urban)
 94                                            G4P <<  94 // 06-01-06 reset currentCouple when StepFunction is changed (V.Ivanchenko)
 95   G4VContinuousDiscreteProcess(name, type)     <<  95 // 10-01-06 PreciseRange -> CSDARange (V.Ivantchenko)
 96 {                                              <<  96 // 18-01-06 Clean up subcutoff including recalculation of presafety (VI)
 97   theParameters = G4EmParameters::Instance();  <<  97 // 20-01-06 Introduce G4EmTableType and reducing number of methods (VI)
 98   SetVerboseLevel(1);                          <<  98 // 22-03-06 Add control on warning printout AlongStep (VI)
 99                                                <<  99 // 23-03-06 Use isIonisation flag (V.Ivanchenko)
100   // low energy limit                          << 100 // 07-06-06 Do not reflect AlongStep in subcutoff regime (V.Ivanchenko)
101   lowestKinEnergy = theParameters->LowestElect << 101 // 14-01-07 add SetEmModel(index) and SetFluctModel() (mma)
102                                                << 102 // 16-01-07 add IonisationTable and IonisationSubTable (V.Ivanchenko)
103   // Size of tables                            << 103 // 16-02-07 set linLossLimit=1.e-6 (V.Ivanchenko)
104   minKinEnergy     = 0.1*CLHEP::keV;           << 104 // 13-03-07 use SafetyHelper instead of navigator (V.Ivanchenko)
105   maxKinEnergy     = 100.0*CLHEP::TeV;         << 105 // 10-04-07 use unique SafetyHelper (V.Ivanchenko)
106   maxKinEnergyCSDA = 1.0*CLHEP::GeV;           << 106 // 12-04-07 Add verbosity at destruction (V.Ivanchenko)
107   nBins            = 84;                       << 107 // 25-04-07 move initialisation of safety helper to BuildPhysicsTable (VI)
108   nBinsCSDA        = 35;                       << 108 // 27-10-07 Virtual functions moved to source (V.Ivanchenko)
109                                                << 109 //
110   invLambdaFactor = 1.0/lambdaFactor;          << 110 // Class Description:
111                                                << 111 //
112   // default linear loss limit                 << 112 // It is the unified energy loss process it calculates the continuous
113   finalRange = 1.*CLHEP::mm;                   << 113 // energy loss for charged particles using a set of Energy Loss
114                                                << 114 // models valid for different energy regions. There are a possibility
115   // run time objects                          << 115 // to create and access to dE/dx and range tables, or to calculate
116   pParticleChange = &fParticleChange;          << 116 // that information on fly.
117   fParticleChange.SetSecondaryWeightByProcess( << 117 // -------------------------------------------------------------------
118   modelManager = new G4EmModelManager();       << 118 //
119   safetyHelper = G4TransportationManager::GetT << 119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
120     ->GetSafetyHelper();                       << 120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
121   aGPILSelection = CandidateForSelection;      << 121 
122                                                << 122 #include "G4VEnergyLossProcess.hh"
123   // initialise model                          << 123 #include "G4LossTableManager.hh"
124   lManager = G4LossTableManager::Instance();   << 124 #include "G4Step.hh"
125   lManager->Register(this);                    << 125 #include "G4ParticleDefinition.hh"
126   isMaster = lManager->IsMaster();             << 126 #include "G4VEmModel.hh"
127                                                << 127 #include "G4VEmFluctuationModel.hh"
128   G4LossTableBuilder* bld = lManager->GetTable << 128 #include "G4DataVector.hh"
129   theDensityFactor = bld->GetDensityFactors(); << 129 #include "G4PhysicsLogVector.hh"
130   theDensityIdx = bld->GetCoupleIndexes();     << 130 #include "G4VParticleChange.hh"
131                                                << 131 #include "G4Gamma.hh"
132   scTracks.reserve(10);                        << 132 #include "G4Electron.hh"
133   secParticles.reserve(12);                    << 133 #include "G4Positron.hh"
134   emModels = new std::vector<G4VEmModel*>;     << 134 #include "G4Proton.hh"
135 }                                              << 135 #include "G4ProcessManager.hh"
136                                                << 136 #include "G4UnitsTable.hh"
137 //....oooOO0OOooo........oooOO0OOooo........oo << 137 #include "G4GenericIon.hh"
138                                                << 138 #include "G4ProductionCutsTable.hh"
139 G4VEnergyLossProcess::~G4VEnergyLossProcess()  << 139 #include "G4Region.hh"
140 {                                              << 140 #include "G4RegionStore.hh"
141   if (isMaster) {                              << 141 #include "G4PhysicsTableHelper.hh"
142     if(nullptr == baseParticle) { delete theDa << 142 #include "G4SafetyHelper.hh"
143     delete theEnergyOfCrossSectionMax;         << 143 #include "G4TransportationManager.hh"
144     if(nullptr != fXSpeaks) {                  << 144 
145       for(auto const & v : *fXSpeaks) { delete << 145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
146       delete fXSpeaks;                         << 146 
147     }                                          << 147 G4VEnergyLossProcess::G4VEnergyLossProcess(const G4String& name, 
148   }                                            << 148              G4ProcessType type): 
149   delete modelManager;                         << 149   G4VContinuousDiscreteProcess(name, type),
150   delete biasManager;                          << 150   secondaryParticle(0),
151   delete scoffRegions;                         << 151   nSCoffRegions(0),
152   delete emModels;                             << 152   idxSCoffRegions(0),
153   lManager->DeRegister(this);                  << 153   nProcesses(0),
154 }                                              << 154   theDEDXTable(0),
155                                                << 155   theDEDXSubTable(0),
156 //....oooOO0OOooo........oooOO0OOooo........oo << 156   theDEDXunRestrictedTable(0),
157                                                << 157   theIonisationTable(0),
158 G4double G4VEnergyLossProcess::MinPrimaryEnerg << 158   theIonisationSubTable(0),
159                                                << 159   theRangeTableForLoss(0),
160                                                << 160   theCSDARangeTable(0),
161 {                                              << 161   theSecondaryRangeTable(0),
162   return cut;                                  << 162   theInverseRangeTable(0),
163 }                                              << 163   theLambdaTable(0),
164                                                << 164   theSubLambdaTable(0),
165 //....oooOO0OOooo........oooOO0OOooo........oo << 165   theDEDXAtMaxEnergy(0),
166                                                << 166   theRangeAtMaxEnergy(0),
167 void G4VEnergyLossProcess::AddEmModel(G4int or << 167   theEnergyOfCrossSectionMax(0),
168                                       G4VEmFlu << 168   theCrossSectionMax(0),
169                                       const G4 << 169   baseParticle(0),
170 {                                              << 170   minSubRange(0.1),
171   if(nullptr == ptr) { return; }               << 171   lossFluctuationFlag(true),
172   G4VEmFluctuationModel* afluc = (nullptr == f << 172   rndmStepFlag(false),
173   modelManager->AddEmModel(order, ptr, afluc,  << 173   tablesAreBuilt(false),
174   ptr->SetParticleChange(pParticleChange, aflu << 174   integral(true),
175 }                                              << 175   isIon(false),
176                                                << 176   isIonisation(true),
177 //....oooOO0OOooo........oooOO0OOooo........oo << 177   useSubCutoff(false),
178                                                << 178   particle(0),
179 void G4VEnergyLossProcess::SetEmModel(G4VEmMod << 179   currentCouple(0),
180 {                                              << 180   nWarnings(0),
181   if(nullptr == ptr) { return; }               << 181   mfpKinEnergy(0.0)
182   if(!emModels->empty()) {                     << 182 {
183     for(auto & em : *emModels) { if(em == ptr) << 183   SetVerboseLevel(1);
184   }                                            << 184 
185   emModels->push_back(ptr);                    << 185   // low energy limit
186 }                                              << 186   lowestKinEnergy  = 1.*eV;
187                                                << 187 
188 //....oooOO0OOooo........oooOO0OOooo........oo << 188   // Size of tables assuming spline
189                                                << 189   minKinEnergy     = 0.1*keV;
190 void G4VEnergyLossProcess::SetDynamicMassCharg << 190   maxKinEnergy     = 100.0*TeV;
191                                                << 191   nBins            = 84;
192 {                                              << 192   maxKinEnergyCSDA = 1.0*GeV;
193   massRatio = massratio;                       << 193   nBinsCSDA        = 35;
194   logMassRatio = G4Log(massRatio);             << 194 
195   fFactor = charge2ratio*biasFactor;           << 195   // default linear loss limit for spline
196   if(baseMat) { fFactor *= (*theDensityFactor) << 196   linLossLimit  = 0.01;
197   chargeSqRatio = charge2ratio;                << 197 
198   reduceFactor  = 1.0/(fFactor*massRatio);     << 198   // default dRoverRange and finalRange
199 }                                              << 199   SetStepFunction(0.2, 1.0*mm);
200                                                << 200 
201 //....oooOO0OOooo........oooOO0OOooo........oo << 201   // default lambda factor
202                                                << 202   lambdaFactor  = 0.8;
203 void                                           << 203 
204 G4VEnergyLossProcess::PreparePhysicsTable(cons << 204   // particle types
205 {                                              << 205   theElectron   = G4Electron::Electron();
206   particle = G4EmTableUtil::CheckIon(this, &pa << 206   thePositron   = G4Positron::Positron();
207                                      verboseLe << 207   theGenericIon = 0;
208                                                << 208 
209   if( particle != &part ) {                    << 209   // run time objects
210     if(!isIon) { lManager->RegisterExtraPartic << 210   pParticleChange = &fParticleChange;
211     if(1 < verboseLevel) {                     << 211   modelManager = new G4EmModelManager();
212       G4cout << "### G4VEnergyLossProcess::Pre << 212   safetyHelper = G4TransportationManager::GetTransportationManager()
213              << " interrupted for " << GetProc << 213     ->GetSafetyHelper();
214              << part.GetParticleName() << " is << 214   aGPILSelection = CandidateForSelection;
215              << " spline=" << spline << G4endl << 215 
216     }                                          << 216   // initialise model
217     return;                                    << 217   (G4LossTableManager::Instance())->Register(this);
218   }                                            << 218   fluctModel = 0;
219                                                << 219 
220   tablesAreBuilt = false;                      << 220   scTracks.reserve(5);
221   if (GetProcessSubType() == fIonisation) { Se << 221   secParticles.reserve(5);
222                                                << 222 
223   G4LossTableBuilder* bld = lManager->GetTable << 223   // Data for stragling of ranges from ICRU'37 report
224   lManager->PreparePhysicsTable(&part, this);  << 224   const G4int nrbins = 7;
225                                                << 225   vstrag = new G4PhysicsLogVector(keV, GeV, nrbins-1);
226   // Base particle and set of models can be de << 226   vstrag->SetSpline(true);
227   InitialiseEnergyLossProcess(particle, basePa << 227   G4double s[nrbins] = {-0.2, -0.85, -1.3, -1.578, -1.76, -1.85, -1.9};
228                                                << 228   for(G4int i=0; i<nrbins; i++) {vstrag->PutValue(i, s[i]);}
229   // parameters of the process                 << 229 }
230   if(!actLossFluc) { lossFluctuationFlag = the << 230 
231   useCutAsFinalRange = theParameters->UseCutAs << 231 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
232   if(!actMinKinEnergy) { minKinEnergy = thePar << 232 
233   if(!actMaxKinEnergy) { maxKinEnergy = thePar << 233 G4VEnergyLossProcess::~G4VEnergyLossProcess()
234   if(!actBinning) { nBins = theParameters->Num << 234 {
235   maxKinEnergyCSDA = theParameters->MaxEnergyF << 235   if(1 < verboseLevel) 
236   nBinsCSDA = theParameters->NumberOfBinsPerDe << 236     G4cout << "G4VEnergyLossProcess destruct " << GetProcessName() 
237     *G4lrint(std::log10(maxKinEnergyCSDA/minKi << 237      << G4endl;
238   if(!actLinLossLimit) { linLossLimit = thePar << 238   delete vstrag;
239   lambdaFactor = theParameters->LambdaFactor() << 239   Clear();
240   invLambdaFactor = 1.0/lambdaFactor;          << 240 
241   if(isMaster) { SetVerboseLevel(theParameters << 241   if ( !baseParticle ) {
242   else { SetVerboseLevel(theParameters->Worker << 242     if(theDEDXTable && theRangeTableForLoss) {
243   // integral option may be disabled           << 243       if(theIonisationTable == theDEDXTable) theIonisationTable = 0;
244   if(!theParameters->Integral()) { fXSType = f << 244       theDEDXTable->clearAndDestroy();
245                                                << 245       delete theDEDXTable;
246   theParameters->DefineRegParamForLoss(this);  << 246       if(theDEDXSubTable) {
247                                                << 247   if(theIonisationSubTable == theDEDXSubTable) 
248   fRangeEnergy = 0.0;                          << 248     theIonisationSubTable = 0;
249                                                << 249   theDEDXSubTable->clearAndDestroy();
250   G4double initialCharge = particle->GetPDGCha << 250         delete theDEDXSubTable;
251   G4double initialMass   = particle->GetPDGMas << 251       }
252                                                << 252     }
253   theParameters->FillStepFunction(particle, th << 253     if(theIonisationTable) {
254                                                << 254       theIonisationTable->clearAndDestroy(); 
255   // parameters for scaling from the base part << 255       delete theIonisationTable;
256   if (nullptr != baseParticle) {               << 256     }
257     massRatio    = (baseParticle->GetPDGMass() << 257     if(theIonisationSubTable) {
258     logMassRatio = G4Log(massRatio);           << 258       theIonisationSubTable->clearAndDestroy(); 
259     G4double q = initialCharge/baseParticle->G << 259       delete theIonisationSubTable;
260     chargeSqRatio = q*q;                       << 260     }
261     if(chargeSqRatio > 0.0) { reduceFactor = 1 << 261     if(theDEDXunRestrictedTable && theCSDARangeTable) {
262   }                                            << 262        theDEDXunRestrictedTable->clearAndDestroy();
263   lowestKinEnergy = (initialMass < CLHEP::MeV) << 263        delete theDEDXunRestrictedTable;
264     ? theParameters->LowestElectronEnergy()    << 264     }
265     : theParameters->LowestMuHadEnergy();      << 265     if(theCSDARangeTable) {
266                                                << 266       theCSDARangeTable->clearAndDestroy();
267   // Tables preparation                        << 267       delete theCSDARangeTable;
268   if (isMaster && nullptr == baseParticle) {   << 268     }
269     if(nullptr == theData) { theData = new G4E << 269     if(theRangeTableForLoss) {
270                                                << 270       theRangeTableForLoss->clearAndDestroy();
271     if(nullptr != theDEDXTable && isIonisation << 271       delete theRangeTableForLoss;
272       if(nullptr != theIonisationTable && theD << 272     }
273   theData->CleanTable(0);                      << 273     if(theInverseRangeTable) {
274   theDEDXTable = theIonisationTable;           << 274       theInverseRangeTable->clearAndDestroy();
275   theIonisationTable = nullptr;                << 275       delete theInverseRangeTable;
276       }                                        << 276     }
277     }                                          << 277     if(theLambdaTable) {
278                                                << 278       theLambdaTable->clearAndDestroy();
279     theDEDXTable = theData->MakeTable(theDEDXT << 279       delete theLambdaTable;
280     bld->InitialiseBaseMaterials(theDEDXTable) << 280     }
281     theData->UpdateTable(theIonisationTable, 1 << 281     if(theSubLambdaTable) {
282                                                << 282       theSubLambdaTable->clearAndDestroy();
283     if (theParameters->BuildCSDARange()) {     << 283       delete theSubLambdaTable;
284       theDEDXunRestrictedTable = theData->Make << 284     }
285       if(isIonisation) { theCSDARangeTable = t << 285   }
286     }                                          << 286 
287                                                << 287   delete modelManager;
288     theLambdaTable = theData->MakeTable(4);    << 288   (G4LossTableManager::Instance())->DeRegister(this);
289     if(isIonisation) {                         << 289 }
290       theRangeTableForLoss = theData->MakeTabl << 290 
291       theInverseRangeTable = theData->MakeTabl << 291 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
292     }                                          << 292 
293   }                                            << 293 void G4VEnergyLossProcess::Clear()
294                                                << 294 {
295   // forced biasing                            << 295   if(1 < verboseLevel) { 
296   if(nullptr != biasManager) {                 << 296     G4cout << "G4VEnergyLossProcess::Clear() for " << GetProcessName() 
297     biasManager->Initialise(part,GetProcessNam << 297      << G4endl;
298     biasFlag = false;                          << 298   }
299   }                                            << 299   delete [] theDEDXAtMaxEnergy;
300   baseMat = bld->GetBaseMaterialFlag();        << 300   delete [] theRangeAtMaxEnergy;
301   numberOfModels = modelManager->NumberOfModel << 301   delete [] theEnergyOfCrossSectionMax;
302   currentModel = modelManager->GetModel(0);    << 302   delete [] theCrossSectionMax;
303   G4EmTableUtil::UpdateModels(this, modelManag << 303   delete [] idxSCoffRegions;
304                               numberOfModels,  << 304 
305                               mainSecondaries, << 305   theDEDXAtMaxEnergy = 0;
306                               theParameters->U << 306   theRangeAtMaxEnergy = 0;
307   theCuts = modelManager->Initialise(particle, << 307   theEnergyOfCrossSectionMax = 0;
308                                      verboseLe << 308   theCrossSectionMax = 0;
309   // subcut processor                          << 309   tablesAreBuilt = false;
310   if(isIonisation) {                           << 310 
311     subcutProducer = lManager->SubCutProducer( << 311   //scTracks.clear();
312   }                                            << 312   scProcesses.clear();
313   if(1 == nSCoffRegions) {                     << 313   nProcesses = 0;
314     if((*scoffRegions)[0]->GetName() == "Defau << 314 }
315       delete scoffRegions;                     << 315 
316       scoffRegions = nullptr;                  << 316 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
317       nSCoffRegions = 0;                       << 317 
318     }                                          << 318 void 
319   }                                            << 319 G4VEnergyLossProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
320                                                << 320 {
321   if(1 < verboseLevel) {                       << 321   if(1 < verboseLevel) {
322     G4cout << "G4VEnergyLossProcess::PrepearPh << 322     G4cout << "G4VEnergyLossProcess::PreparePhysicsTable for "
323            << " for " << GetProcessName() << " << 323            << GetProcessName()
324            << " isIon= " << isIon << " spline= << 324            << " for " << part.GetParticleName()
325     if(baseParticle) {                         << 325            << G4endl;
326       G4cout << "; base: " << baseParticle->Ge << 326   }
327     }                                          << 327 
328     G4cout << G4endl;                          << 328   currentCouple = 0;
329     G4cout << " chargeSqRatio= " << chargeSqRa << 329   preStepLambda = 0.0;
330            << " massRatio= " << massRatio      << 330   mfpKinEnergy  = DBL_MAX;
331            << " reduceFactor= " << reduceFacto << 331   fRange        = DBL_MAX;
332     if (nSCoffRegions > 0) {                   << 332   preStepKinEnergy = 0.0;
333       G4cout << " SubCut secondary production  << 333   chargeSqRatio = 1.0;
334       for (G4int i=0; i<nSCoffRegions; ++i) {  << 334   massRatio = 1.0;
335         const G4Region* r = (*scoffRegions)[i] << 335   reduceFactor = 1.0;
336         G4cout << "           " << r->GetName( << 336 
337       }                                        << 337   G4LossTableManager* lManager = G4LossTableManager::Instance();
338     } else if(nullptr != subcutProducer) {     << 338 
339       G4cout << " SubCut secondary production  << 339   // Are particle defined?
340     }                                          << 340   if( !particle ) {
341   }                                            << 341     particle = &part;
342 }                                              << 342     if(part.GetParticleType() == "nucleus") {
343                                                << 343       if(!theGenericIon) theGenericIon = G4GenericIon::GenericIon();
344 //....oooOO0OOooo........oooOO0OOooo........oo << 344       if(particle == theGenericIon) { isIon = true; }
345                                                << 345       else if(part.GetPDGCharge() > eplus) {
346 void G4VEnergyLossProcess::BuildPhysicsTable(c << 346   isIon = true; 
347 {                                              << 347 
348   if(1 < verboseLevel) {                       << 348   // generic ions created on-fly
349     G4cout << "### G4VEnergyLossProcess::Build << 349   if(part.GetPDGCharge() > 2.5*eplus) {
350            << GetProcessName()                 << 350     particle = theGenericIon;
351            << " and particle " << part.GetPart << 351   }
352            << "; the first particle " << parti << 352       }
353     if(baseParticle) {                         << 353     }
354       G4cout << "; base: " << baseParticle->Ge << 354   }
355     }                                          << 355 
356     G4cout << G4endl;                          << 356   if( particle != &part) {
357     G4cout << "    TablesAreBuilt= " << tables << 357     if(part.GetParticleType() == "nucleus") {
358            << " spline=" << spline << " ptr: " << 358       isIon = true;
359   }                                            << 359       lManager->RegisterIon(&part, this);
360                                                << 360     } else { 
361   if(&part == particle) {                      << 361       lManager->RegisterExtraParticle(&part, this);
362     if(isMaster) {                             << 362     }
363       lManager->BuildPhysicsTable(particle, th << 363     return;
364                                                << 364   }
365     } else {                                   << 365 
366       const auto masterProcess =               << 366   Clear();
367         static_cast<const G4VEnergyLossProcess << 367 
368                                                << 368   // Base particle and set of models can be defined here
369       numberOfModels = modelManager->NumberOfM << 369   InitialiseEnergyLossProcess(particle, baseParticle);
370       G4EmTableUtil::BuildLocalElossProcess(th << 370 
371                                             pa << 371   // Tables preparation
372       tablesAreBuilt = true;                   << 372   if (!baseParticle) {
373       baseMat = masterProcess->UseBaseMaterial << 373     
374       lManager->LocalPhysicsTables(particle, t << 374     theDEDXTable = G4PhysicsTableHelper::PreparePhysicsTable(theDEDXTable);
375     }                                          << 375     if (lManager->BuildCSDARange()) {
376                                                << 376       theDEDXunRestrictedTable = 
377     // needs to be done only once              << 377   G4PhysicsTableHelper::PreparePhysicsTable(theDEDXunRestrictedTable);
378     safetyHelper->InitialiseHelper();          << 378       theCSDARangeTable = 
379   }                                            << 379   G4PhysicsTableHelper::PreparePhysicsTable(theCSDARangeTable);
380   // Added tracking cut to avoid tracking arti << 380     }
381   // and identified deexcitation flag          << 381 
382   if(isIonisation) {                           << 382     theRangeTableForLoss = 
383     atomDeexcitation = lManager->AtomDeexcitat << 383       G4PhysicsTableHelper::PreparePhysicsTable(theRangeTableForLoss);
384     if(nullptr != atomDeexcitation) {          << 384     theInverseRangeTable = 
385       if(atomDeexcitation->IsPIXEActive()) { u << 385       G4PhysicsTableHelper::PreparePhysicsTable(theInverseRangeTable);
386     }                                          << 386   
387   }                                            << 387     theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable);
388                                                << 388     if (nSCoffRegions) {
389   // protection against double printout        << 389       theDEDXSubTable = 
390   if(theParameters->IsPrintLocked()) { return; << 390   G4PhysicsTableHelper::PreparePhysicsTable(theDEDXSubTable);
391                                                << 391       theSubLambdaTable = 
392   // explicitly defined printout by particle n << 392   G4PhysicsTableHelper::PreparePhysicsTable(theSubLambdaTable);
393   G4String num = part.GetParticleName();       << 393     }
394   if(1 < verboseLevel ||                       << 394   }
395      (0 < verboseLevel && (num == "e-" ||      << 395 
396                            num == "e+"    || n << 396   G4double initialCharge = particle->GetPDGCharge();
397                            num == "mu-"   || n << 397   G4double initialMass   = particle->GetPDGMass();
398                            num == "pi+"   || n << 398 
399                            num == "kaon+" || n << 399   if (baseParticle) {
400                            num == "alpha" || n << 400     massRatio = (baseParticle->GetPDGMass())/initialMass;
401                            num == "GenericIon" << 401     G4double q = initialCharge/baseParticle->GetPDGCharge();
402     StreamInfo(G4cout, part);                  << 402     chargeSqRatio = q*q;
403   }                                            << 403     if(chargeSqRatio > 0.0) reduceFactor = 1.0/(chargeSqRatio*massRatio);
404   if(1 < verboseLevel) {                       << 404   }
405     G4cout << "### G4VEnergyLossProcess::Build << 405 
406            << GetProcessName()                 << 406   theCuts = modelManager->Initialise(particle, secondaryParticle, 
407            << " and particle " << part.GetPart << 407              minSubRange, verboseLevel);
408     if(isIonisation) { G4cout << "  isIonisati << 408 
409     G4cout << " baseMat=" << baseMat << G4endl << 409   // Sub Cutoff Regime
410   }                                            << 410   if (nSCoffRegions>0) {
411 }                                              << 411     theSubCuts = modelManager->SubCutoff();
412                                                << 412 
413 //....oooOO0OOooo........oooOO0OOooo........oo << 413     const G4ProductionCutsTable* theCoupleTable=
414                                                << 414           G4ProductionCutsTable::GetProductionCutsTable();
415 G4PhysicsTable* G4VEnergyLossProcess::BuildDED << 415     size_t numOfCouples = theCoupleTable->GetTableSize();
416 {                                              << 416     idxSCoffRegions = new G4int[numOfCouples];
417   G4PhysicsTable* table = nullptr;             << 417   
418   G4double emax = maxKinEnergy;                << 418     for (size_t j=0; j<numOfCouples; j++) {
419   G4int bin = nBins;                           << 419 
420                                                << 420       const G4MaterialCutsCouple* couple = 
421   if(fTotal == tType) {                        << 421   theCoupleTable->GetMaterialCutsCouple(j);
422     emax  = maxKinEnergyCSDA;                  << 422       const G4ProductionCuts* pcuts = couple->GetProductionCuts();
423     bin   = nBinsCSDA;                         << 423       G4int reg = 0;
424     table = theDEDXunRestrictedTable;          << 424       for(G4int i=0; i<nSCoffRegions; i++) {
425   } else if(fRestricted == tType) {            << 425         if( pcuts == scoffRegions[i]->GetProductionCuts()) reg = 1;
426     table = theDEDXTable;                      << 426       }
427   } else {                                     << 427       idxSCoffRegions[j] = reg;
428     G4cout << "G4VEnergyLossProcess::BuildDEDX << 428     }
429            << tType << G4endl;                 << 429   }
430   }                                            << 430 
431   if(1 < verboseLevel) {                       << 431   lManager->EnergyLossProcessIsInitialised(particle, this);
432     G4cout << "G4VEnergyLossProcess::BuildDEDX << 432 
433            << " for " << GetProcessName()      << 433   if (1 < verboseLevel) {
434            << " and " << particle->GetParticle << 434     G4cout << "G4VEnergyLossProcess::Initialise() is done "
435      << "spline=" << spline << G4endl;         << 435            << " for local " << particle->GetParticleName()
436   }                                            << 436      << " isIon= " << isIon
437   if(nullptr == table) { return table; }       << 437            << " chargeSqRatio= " << chargeSqRatio
438                                                << 438            << " massRatio= " << massRatio
439   G4LossTableBuilder* bld = lManager->GetTable << 439            << " reduceFactor= " << reduceFactor << G4endl;
440   G4EmTableUtil::BuildDEDXTable(this, particle << 440     if (nSCoffRegions) {
441                                 table, minKinE << 441       G4cout << " SubCutoff Regime is ON for regions: " << G4endl;
442                                 verboseLevel,  << 442       for (G4int i=0; i<nSCoffRegions; i++) {
443   return table;                                << 443         const G4Region* r = scoffRegions[i];
444 }                                              << 444   G4cout << "           " << r->GetName() << G4endl;
445                                                << 445       }
446 //....oooOO0OOooo........oooOO0OOooo........oo << 446     }
447                                                << 447   }
448 G4PhysicsTable* G4VEnergyLossProcess::BuildLam << 448 }
449 {                                              << 449 
450   if(nullptr == theLambdaTable) { return theLa << 450 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
451                                                << 451 
452   G4double scale = theParameters->MaxKinEnergy << 452 void G4VEnergyLossProcess::BuildPhysicsTable(const G4ParticleDefinition& part)
453   G4int nbin =                                 << 453 {
454     theParameters->NumberOfBinsPerDecade()*G4l << 454   if(1 < verboseLevel) {
455   scale = nbin/G4Log(scale);                   << 455     G4cout << "### G4VEnergyLossProcess::BuildPhysicsTable() for "
456                                                << 456            << GetProcessName()
457   G4LossTableBuilder* bld = lManager->GetTable << 457            << " and particle " << part.GetParticleName()
458   G4EmTableUtil::BuildLambdaTable(this, partic << 458            << "; local: " << particle->GetParticleName();
459                                   bld, theLamb << 459     if(baseParticle) G4cout << "; base: " << baseParticle->GetParticleName();
460                                   minKinEnergy << 460     G4cout << G4endl;
461                                   verboseLevel << 461   }
462   return theLambdaTable;                       << 462 
463 }                                              << 463   if(&part == particle) {
464                                                << 464     if(!tablesAreBuilt) {
465 //....oooOO0OOooo........oooOO0OOooo........oo << 465       G4LossTableManager::Instance()->BuildPhysicsTable(particle, this);
466                                                << 466     }
467 void G4VEnergyLossProcess::StreamInfo(std::ost << 467     if(!baseParticle) {
468                 const G4ParticleDefinition& pa << 468       if(0 < verboseLevel) PrintInfoDefinition();
469 {                                              << 469     
470   G4String indent = (rst ? "  " : "");         << 470       // needs to be done only once
471   out << std::setprecision(6);                 << 471       safetyHelper->InitialiseHelper();
472   out << G4endl << indent << GetProcessName()  << 472     }
473   if (!rst) out << " for " << part.GetParticle << 473   }
474   out << "  XStype:" << fXSType                << 474 
475       << "  SubType=" << GetProcessSubType() < << 475   if(1 < verboseLevel) {
476       << "      dE/dx and range tables from "  << 476     G4cout << "### G4VEnergyLossProcess::BuildPhysicsTable() done for "
477       << G4BestUnit(minKinEnergy,"Energy")     << 477            << GetProcessName()
478       << " to " << G4BestUnit(maxKinEnergy,"En << 478            << " and particle " << part.GetParticleName();
479       << " in " << nBins << " bins" << G4endl  << 479     if(isIonisation) G4cout << "  isIonisation  flag = 1";
480       << "      Lambda tables from threshold t << 480     G4cout << G4endl;
481       << G4BestUnit(maxKinEnergy,"Energy")     << 481   }
482       << ", " << theParameters->NumberOfBinsPe << 482 }
483       << " bins/decade, spline: " << spline    << 483 
484       << G4endl;                               << 484 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
485   if(nullptr != theRangeTableForLoss && isIoni << 485 
486     out << "      StepFunction=(" << dRoverRan << 486 void G4VEnergyLossProcess::ActivateSubCutoff(G4bool val, const G4Region* r)
487         << finalRange/mm << " mm)"             << 487 {
488         << ", integ: " << fXSType              << 488   G4RegionStore* regionStore = G4RegionStore::GetInstance();
489         << ", fluct: " << lossFluctuationFlag  << 489   if(val) {
490         << ", linLossLim= " << linLossLimit    << 490     useSubCutoff = true;
491         << G4endl;                             << 491     if (!r) {r = regionStore->GetRegion("DefaultRegionForTheWorld", false);}
492   }                                            << 492     if (nSCoffRegions) {
493   StreamProcessInfo(out);                      << 493       for (G4int i=0; i<nSCoffRegions; i++) {
494   modelManager->DumpModelList(out, verboseLeve << 494   if (r == scoffRegions[i]) return;
495   if(nullptr != theCSDARangeTable && isIonisat << 495       }
496     out << "      CSDA range table up"         << 496     }
497         << " to " << G4BestUnit(maxKinEnergyCS << 497     scoffRegions.push_back(r);
498         << " in " << nBinsCSDA << " bins" << G << 498     nSCoffRegions++;
499   }                                            << 499   } else {
500   if(nSCoffRegions>0 && isIonisation) {        << 500     useSubCutoff = false;
501     out << "      Subcutoff sampling in " << n << 501   }
502         << " regions" << G4endl;               << 502 }
503   }                                            << 503 
504   if(2 < verboseLevel) {                       << 504 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
505     for(std::size_t i=0; i<7; ++i) {           << 505 
506       auto ta = theData->Table(i);             << 506 G4PhysicsTable* G4VEnergyLossProcess::BuildDEDXTable(G4EmTableType tType)
507       out << "      " << tnames[i] << " addres << 507 {
508       if(nullptr != ta) { out << *ta << G4endl << 508   if(1 < verboseLevel) {
509     }                                          << 509     G4cout << "G4VEnergyLossProcess::BuildDEDXTable() of type " << tType
510   }                                            << 510      << " for " << GetProcessName()
511 }                                              << 511            << " and particle " << particle->GetParticleName()
512                                                << 512            << G4endl;
513 //....oooOO0OOooo........oooOO0OOooo........oo << 513   }
514                                                << 514   G4PhysicsTable* table = 0;
515 void G4VEnergyLossProcess::ActivateSubCutoff(c << 515   G4double emin = minKinEnergy;
516 {                                              << 516   G4double emax = maxKinEnergy;
517   if(nullptr == scoffRegions) {                << 517   G4int bin = nBins;
518     scoffRegions = new std::vector<const G4Reg << 518 
519   }                                            << 519   if(fTotal == tType) {
520   // the region is in the list                 << 520     emax  = maxKinEnergyCSDA;
521   if(!scoffRegions->empty()) {                 << 521     bin   = nBinsCSDA;
522     for (auto & reg : *scoffRegions) {         << 522     table = theDEDXunRestrictedTable;
523       if (reg == r) { return; }                << 523   } else if(fRestricted == tType) {
524     }                                          << 524     table = theDEDXTable;
525   }                                            << 525     if(theIonisationTable) 
526   // new region                                << 526       table = G4PhysicsTableHelper::PreparePhysicsTable(theIonisationTable); 
527   scoffRegions->push_back(r);                  << 527   } else if(fSubRestricted == tType) {    
528   ++nSCoffRegions;                             << 528     table = theDEDXSubTable;
529 }                                              << 529     if(theIonisationSubTable) 
530                                                << 530       table = G4PhysicsTableHelper::PreparePhysicsTable(theIonisationSubTable); 
531 //....oooOO0OOooo........oooOO0OOooo........oo << 531   } else {
532                                                << 532     G4cout << "G4VEnergyLossProcess::BuildDEDXTable WARNING: wrong type "
533 G4bool G4VEnergyLossProcess::IsRegionForCubcut << 533      << tType << G4endl;
534 {                                              << 534   }
535   if(0 == nSCoffRegions) { return true; }      << 535 
536   const G4Region* r = aTrack.GetVolume()->GetL << 536   // Access to materials
537   for(auto & reg : *scoffRegions) {            << 537   const G4ProductionCutsTable* theCoupleTable=
538     if(r == reg) { return true; }              << 538         G4ProductionCutsTable::GetProductionCutsTable();
539   }                                            << 539   size_t numOfCouples = theCoupleTable->GetTableSize();
540   return false;                                << 540 
541 }                                              << 541   if(1 < verboseLevel) {
542                                                << 542     G4cout << numOfCouples << " materials"
543 //....oooOO0OOooo........oooOO0OOooo........oo << 543            << " minKinEnergy= " << minKinEnergy
544                                                << 544            << " maxKinEnergy= " << maxKinEnergy
545 void G4VEnergyLossProcess::StartTracking(G4Tra << 545            << " EmTableType= " << tType
546 {                                              << 546            << " table= " << table
547   // reset parameters for the new track        << 547            << G4endl;
548   theNumberOfInteractionLengthLeft = -1.0;     << 548   }
549   mfpKinEnergy = DBL_MAX;                      << 549   if(!table) return table;
550   preStepLambda = 0.0;                         << 550 
551   currentCouple = nullptr;                     << 551   for(size_t i=0; i<numOfCouples; i++) {
552                                                << 552 
553   // reset ion                                 << 553     if(1 < verboseLevel) {
554   if(isIon) {                                  << 554       G4cout << "G4VEnergyLossProcess::BuildDEDXVector flag=  " 
555     const G4double newmass = track->GetDefinit << 555        << table->GetFlag(i) << G4endl;
556     massRatio = (nullptr == baseParticle) ? CL << 556     }
557       : baseParticle->GetPDGMass()/newmass;    << 557     if (table->GetFlag(i)) {
558     logMassRatio = G4Log(massRatio);           << 558 
559   }                                            << 559       // create physics vector and fill it
560   // forced biasing only for primary particles << 560       const G4MaterialCutsCouple* couple = 
561   if(nullptr != biasManager) {                 << 561   theCoupleTable->GetMaterialCutsCouple(i);
562     if(0 == track->GetParentID()) {            << 562       G4PhysicsVector* aVector = new G4PhysicsLogVector(emin, emax, bin);
563       biasFlag = true;                         << 563       aVector->SetSpline((G4LossTableManager::Instance())->SplineFlag());
564       biasManager->ResetForcedInteraction();   << 564 
565     }                                          << 565       modelManager->FillDEDXVector(aVector, couple, tType);
566   }                                            << 566 
567 }                                              << 567       // Insert vector for this material into the table
568                                                << 568       G4PhysicsTableHelper::SetPhysicsVector(table, i, aVector);
569 //....oooOO0OOooo........oooOO0OOooo........oo << 569     }
570                                                << 570   }
571 G4double G4VEnergyLossProcess::AlongStepGetPhy << 571 
572                              const G4Track& tr << 572   if(1 < verboseLevel) {
573                              G4GPILSelection*  << 573     G4cout << "G4VEnergyLossProcess::BuildDEDXTable(): table is built for "
574 {                                              << 574            << particle->GetParticleName()
575   G4double x = DBL_MAX;                        << 575            << " and process " << GetProcessName()
576   *selection = aGPILSelection;                 << 576            << G4endl;
577   if(isIonisation && currentModel->IsActive(pr << 577     //    if(2 < verboseLevel) G4cout << (*table) << G4endl;
578     GetScaledRangeForScaledEnergy(preStepScale << 578   }
579     x = (useCutAsFinalRange) ? std::min(finalR << 579 
580       currentCouple->GetProductionCuts()->GetP << 580   return table;
581     x = (fRange > x) ? fRange*dRoverRange + x* << 581 }
582       : fRange;                                << 582 
583     /*                                         << 583 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
584       G4cout<<"AlongStepGPIL: " << GetProcessN << 584 
585   << " fRange=" << fRange << " finR=" << finR  << 585 G4PhysicsTable* G4VEnergyLossProcess::BuildLambdaTable(G4EmTableType tType)
586     */                                         << 586 {
587   }                                            << 587   G4PhysicsTable* table = 0;
588   return x;                                    << 588 
589 }                                              << 589   if(fRestricted == tType) {
590                                                << 590     table = theLambdaTable;
591 //....oooOO0OOooo........oooOO0OOooo........oo << 591   } else if(fSubRestricted == tType) {    
592                                                << 592     table = theSubLambdaTable;
593 G4double G4VEnergyLossProcess::PostStepGetPhys << 593   } else {
594                              const G4Track& tr << 594     G4cout << "G4VEnergyLossProcess::BuildLambdaTable WARNING: wrong type "
595                              G4double   previo << 595      << tType << G4endl;
596                              G4ForceCondition* << 596   }
597 {                                              << 597 
598   // condition is set to "Not Forced"          << 598   if(1 < verboseLevel) {
599   *condition = NotForced;                      << 599     G4cout << "G4VEnergyLossProcess::BuildLambdaTable() of type "
600   G4double x = DBL_MAX;                        << 600      << tType << " for process "
601                                                << 601            << GetProcessName() << " and particle "
602   // initialisation of material, mass, charge, << 602            << particle->GetParticleName()
603   // at the beginning of the step              << 603            << " EmTableType= " << tType
604   DefineMaterial(track.GetMaterialCutsCouple() << 604            << " table= " << table
605   preStepKinEnergy       = track.GetKineticEne << 605            << G4endl;
606   preStepScaledEnergy    = preStepKinEnergy*ma << 606   }
607   SelectModel(preStepScaledEnergy);            << 607   if(!table) {return table;}
608                                                << 608 
609   if(!currentModel->IsActive(preStepScaledEner << 609   // Access to materials
610     theNumberOfInteractionLengthLeft = -1.0;   << 610   const G4ProductionCutsTable* theCoupleTable=
611     mfpKinEnergy = DBL_MAX;                    << 611         G4ProductionCutsTable::GetProductionCutsTable();
612     preStepLambda = 0.0;                       << 612   size_t numOfCouples = theCoupleTable->GetTableSize();
613     currentInteractionLength = DBL_MAX;        << 613 
614     return x;                                  << 614   for(size_t i=0; i<numOfCouples; i++) {
615   }                                            << 615 
616                                                << 616     if (table->GetFlag(i)) {
617   // change effective charge of a charged part << 617 
618   if(isIon) {                                  << 618       // create physics vector and fill it
619     const G4double q2 = currentModel->ChargeSq << 619       const G4MaterialCutsCouple* couple = 
620     fFactor = q2*biasFactor;                   << 620   theCoupleTable->GetMaterialCutsCouple(i);
621     if(baseMat) { fFactor *= (*theDensityFacto << 621       G4double cut = (*theCuts)[i];
622     reduceFactor = 1.0/(fFactor*massRatio);    << 622       if(fSubRestricted == tType) cut = (*theSubCuts)[i]; 
623     if (lossFluctuationFlag) {                 << 623       G4PhysicsVector* aVector = LambdaPhysicsVector(couple, cut);
624       auto fluc = currentModel->GetModelOfFluc << 624       modelManager->FillLambdaVector(aVector, couple, true, tType);
625       fluc->SetParticleAndCharge(track.GetDefi << 625 
626     }                                          << 626       // Insert vector for this material into the table
627   }                                            << 627       G4PhysicsTableHelper::SetPhysicsVector(table, i, aVector);
628                                                << 628     }
629   // forced biasing only for primary particles << 629   }
630   if(biasManager) {                            << 630 
631     if(0 == track.GetParentID() && biasFlag && << 631   if(1 < verboseLevel) {
632        biasManager->ForcedInteractionRegion((G << 632     G4cout << "Lambda table is built for "
633       return biasManager->GetStepLimit((G4int) << 633            << particle->GetParticleName()
634     }                                          << 634            << G4endl;
635   }                                            << 635   }
636                                                << 636 
637   ComputeLambdaForScaledEnergy(preStepScaledEn << 637   return table;
638                                                << 638 }
639   // zero cross section                        << 639 
640   if(preStepLambda <= 0.0) {                   << 640 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
641     theNumberOfInteractionLengthLeft = -1.0;   << 641 
642     currentInteractionLength = DBL_MAX;        << 642 G4double G4VEnergyLossProcess::GetContinuousStepLimit(
643   } else {                                     << 643     const G4Track&,
644                                                << 644                 G4double, G4double, G4double&)
645     // non-zero cross section                  << 645 {
646     if (theNumberOfInteractionLengthLeft < 0.0 << 646   return DBL_MAX;
647                                                << 647 }
648       // beggining of tracking (or just after  << 648 
649       theNumberOfInteractionLengthLeft = -G4Lo << 649 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
650       theInitialNumberOfInteractionLength = th << 650 
651                                                << 651 G4double G4VEnergyLossProcess::AlongStepGetPhysicalInteractionLength(
652     } else if(currentInteractionLength < DBL_M << 652                              const G4Track&,
653                                                << 653                              G4double,
654       // subtract NumberOfInteractionLengthLef << 654                              G4double  currentMinStep,
655       theNumberOfInteractionLengthLeft -=      << 655                              G4double&,
656         previousStepSize/currentInteractionLen << 656                              G4GPILSelection* selection)
657                                                << 657 {
658       theNumberOfInteractionLengthLeft =       << 658   G4double x = DBL_MAX;
659         std::max(theNumberOfInteractionLengthL << 659   *selection = aGPILSelection;
660     }                                          << 660   if(isIonisation) {
661                                                << 661     fRange = GetScaledRangeForScaledEnergy(preStepScaledEnergy)*reduceFactor;
662     // new mean free path and step limit       << 662 
663     currentInteractionLength = 1.0/preStepLamb << 663     x = fRange;
664     x = theNumberOfInteractionLengthLeft * cur << 664     G4double y = x*dRoverRange;
665   }                                            << 665 
666 #ifdef G4VERBOSE                               << 666     if(x > finalRange && y < currentMinStep) { 
667   if (verboseLevel>2) {                        << 667       x = y + finalRange*(1.0 - dRoverRange)*(2.0 - finalRange/fRange);
668     G4cout << "G4VEnergyLossProcess::PostStepG << 668     } else if (rndmStepFlag) {x = SampleRange();}
669     G4cout << "[ " << GetProcessName() << "]"  << 669     //G4cout<<GetProcessName()<<": e= "<<preStepKinEnergy
670     G4cout << " for " << track.GetDefinition() << 670     //  <<" range= "<<fRange <<" cMinSt="<<currentMinStep
671            << " in Material  " <<  currentMate << 671     //  << " limit= " << x <<G4endl;
672            << " Ekin(MeV)= " << preStepKinEner << 672   }
673            << " track material: " << track.Get << 673   //G4cout<<GetProcessName()<<": e= "<<preStepKinEnergy
674            <<G4endl;                           << 674   //  <<" stepLimit= "<<x<<G4endl;
675     G4cout << "MeanFreePath = " << currentInte << 675   return x;
676            << "InteractionLength= " << x/cm << << 676 }
677   }                                            << 677 
678 #endif                                         << 678 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
679   return x;                                    << 679 
680 }                                              << 680 G4double G4VEnergyLossProcess::GetMeanFreePath(
681                                                << 681                              const G4Track& track,
682 //....oooOO0OOooo........oooOO0OOooo........oo << 682                              G4double,
683                                                << 683                              G4ForceCondition* condition)
684 void                                           << 684 
685 G4VEnergyLossProcess::ComputeLambdaForScaledEn << 685 {
686 {                                              << 686   *condition = NotForced;
687   // cross section increased with energy       << 687   return MeanFreePath(track);
688   if(fXSType == fEmIncreasing) {               << 688 }
689     if(e*invLambdaFactor < mfpKinEnergy) {     << 689 
690       preStepLambda = GetLambdaForScaledEnergy << 690 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
691       mfpKinEnergy = (preStepLambda > 0.0) ? e << 691 
692     }                                          << 692 G4double G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength(
693                                                << 693                              const G4Track& track,
694     // cross section has one peak              << 694                              G4double   previousStepSize,
695   } else if(fXSType == fEmOnePeak) {           << 695                              G4ForceCondition* condition)
696     const G4double epeak = (*theEnergyOfCrossS << 696 {
697     if(e <= epeak) {                           << 697   // condition is set to "Not Forced"
698       if(e*invLambdaFactor < mfpKinEnergy) {   << 698   *condition = NotForced;
699         preStepLambda = GetLambdaForScaledEner << 699   G4double x = DBL_MAX;
700         mfpKinEnergy = (preStepLambda > 0.0) ? << 700 
701       }                                        << 701   // initialisation of material, mass, charge, model at the beginning of the step
702     } else if(e < mfpKinEnergy) {              << 702   DefineMaterial(track.GetMaterialCutsCouple());
703       const G4double e1 = std::max(epeak, e*la << 703 
704       mfpKinEnergy = e1;                       << 704   const G4ParticleDefinition* currPart = track.GetDefinition();
705       preStepLambda = GetLambdaForScaledEnergy << 705   if(theGenericIon == particle) {
706     }                                          << 706     massRatio = proton_mass_c2/currPart->GetPDGMass();
707                                                << 707   }  
708     // cross section has more than one peaks   << 708   preStepKinEnergy    = track.GetKineticEnergy();
709   } else if(fXSType == fEmTwoPeaks) {          << 709   preStepScaledEnergy = preStepKinEnergy*massRatio;
710     G4TwoPeaksXS* xs = (*fXSpeaks)[basedCouple << 710   SelectModel(preStepScaledEnergy);
711     const G4double e1peak = xs->e1peak;        << 711 
712                                                << 712   if(isIon) {
713     // below the 1st peak                      << 713     chargeSqRatio = 
714     if(e <= e1peak) {                          << 714       currentModel->GetChargeSquareRatio(currPart,currentMaterial,preStepKinEnergy);
715       if(e*invLambdaFactor < mfpKinEnergy) {   << 715     reduceFactor  = 1.0/(chargeSqRatio*massRatio);
716         preStepLambda = GetLambdaForScaledEner << 716   }
717         mfpKinEnergy = (preStepLambda > 0.0) ? << 717   //G4cout << "q2= " << chargeSqRatio << " massRatio= " << massRatio << G4endl; 
718       }                                        << 718   // initialisation for sampling of the interaction length 
719       return;                                  << 719   if(previousStepSize <= DBL_MIN) theNumberOfInteractionLengthLeft = -1.0;
720     }                                          << 720   if(theNumberOfInteractionLengthLeft < 0.0) mfpKinEnergy = DBL_MAX;
721     const G4double e1deep = xs->e1deep;        << 721 
722     // above the 1st peak, below the deep      << 722   // compute mean free path
723     if(e <= e1deep) {                          << 723   if(preStepScaledEnergy < mfpKinEnergy) {
724       if(mfpKinEnergy >= e1deep || e <= mfpKin << 724     if (integral) ComputeLambdaForScaledEnergy(preStepScaledEnergy);
725         const G4double e1 = std::max(e1peak, e << 725     else  preStepLambda = GetLambdaForScaledEnergy(preStepScaledEnergy);
726         mfpKinEnergy = e1;                     << 726     if(preStepLambda <= DBL_MIN) mfpKinEnergy = 0.0;
727         preStepLambda = GetLambdaForScaledEner << 727   }
728       }                                        << 728 
729       return;                                  << 729   // non-zero cross section
730     }                                          << 730   if(preStepLambda > DBL_MIN) { 
731     const G4double e2peak = xs->e2peak;        << 731     if (theNumberOfInteractionLengthLeft < 0.0) {
732     // above the deep, below 2nd peak          << 732       // beggining of tracking (or just after DoIt of this process)
733     if(e <= e2peak) {                          << 733       //G4cout<<"G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength Reset"<<G4endl;
734       if(e*invLambdaFactor < mfpKinEnergy) {   << 734       ResetNumberOfInteractionLengthLeft();
735         mfpKinEnergy = e;                      << 735     } else if(currentInteractionLength < DBL_MAX) {
736         preStepLambda = GetLambdaForScaledEner << 736       // subtract NumberOfInteractionLengthLeft
737       }                                        << 737       SubtractNumberOfInteractionLengthLeft(previousStepSize);
738       return;                                  << 738       if(theNumberOfInteractionLengthLeft < 0.)
739     }                                          << 739   theNumberOfInteractionLengthLeft = perMillion;
740     const G4double e2deep = xs->e2deep;        << 740     }
741     // above the 2nd peak, below the deep      << 741 
742     if(e <= e2deep) {                          << 742     // get mean free path and step limit
743       if(mfpKinEnergy >= e2deep || e <= mfpKin << 743     currentInteractionLength = 1.0/preStepLambda;
744         const G4double e1 = std::max(e2peak, e << 744     x = theNumberOfInteractionLengthLeft * currentInteractionLength;
745         mfpKinEnergy = e1;                     << 745 
746         preStepLambda = GetLambdaForScaledEner << 746 #ifdef G4VERBOSE
747       }                                        << 747     if (verboseLevel>2){
748       return;                                  << 748       G4cout << "G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength ";
749     }                                          << 749       G4cout << "[ " << GetProcessName() << "]" << G4endl; 
750     const G4double e3peak = xs->e3peak;        << 750       G4cout << " for " << currPart->GetParticleName() 
751     // above the deep, below 3d peak           << 751              << " in Material  " <<  currentMaterial->GetName()
752     if(e <= e3peak) {                          << 752        << " Ekin(MeV)= " << preStepKinEnergy/MeV 
753       if(e*invLambdaFactor < mfpKinEnergy) {   << 753        <<G4endl;
754         mfpKinEnergy = e;                      << 754       G4cout << "MeanFreePath = " << currentInteractionLength/cm << "[cm]" 
755         preStepLambda = GetLambdaForScaledEner << 755        << "InteractionLength= " << x/cm <<"[cm] " <<G4endl;
756       }                                        << 756     }
757       return;                                  << 757 #endif
758     }                                          << 758     // zero cross section case
759     // above 3d peak                           << 759   } else {
760     if(e <= mfpKinEnergy) {                    << 760     if(theNumberOfInteractionLengthLeft > DBL_MIN && 
761       const G4double e1 = std::max(e3peak, e*l << 761        currentInteractionLength < DBL_MAX) {
762       mfpKinEnergy = e1;                       << 762       // subtract NumberOfInteractionLengthLeft
763       preStepLambda = GetLambdaForScaledEnergy << 763       SubtractNumberOfInteractionLengthLeft(previousStepSize);
764     }                                          << 764       if(theNumberOfInteractionLengthLeft < 0.)
765     // integral method is not used             << 765   theNumberOfInteractionLengthLeft = perMillion;
766   } else {                                     << 766     }
767     preStepLambda = GetLambdaForScaledEnergy(e << 767     currentInteractionLength = DBL_MAX;
768   }                                            << 768   }
769 }                                              << 769   return x;
770                                                << 770 }
771 //....oooOO0OOooo........oooOO0OOooo........oo << 771 
772                                                << 772 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
773 G4VParticleChange* G4VEnergyLossProcess::Along << 773 
774                                                << 774 G4VParticleChange* G4VEnergyLossProcess::AlongStepDoIt(const G4Track& track,
775 {                                              << 775                                                        const G4Step& step)
776   fParticleChange.InitializeForAlongStep(track << 776 {
777   // The process has range table - calculate e << 777   fParticleChange.InitializeForAlongStep(track);
778   if(!isIonisation || !currentModel->IsActive( << 778   // The process has range table - calculate energy loss
779     return &fParticleChange;                   << 779   if(!isIonisation) return &fParticleChange;
780   }                                            << 780 
781                                                << 781   // Get the actual (true) Step length
782   G4double length = step.GetStepLength();      << 782   G4double length = step.GetStepLength();
783   G4double eloss  = 0.0;                       << 783   if(length <= DBL_MIN) return &fParticleChange;
784                                                << 784   G4double eloss  = 0.0;
785   /*                                           << 785   G4double esecdep = 0.0;
786   if(-1 < verboseLevel) {                      << 786  
787     const G4ParticleDefinition* d = track.GetP << 787   /*  
788     G4cout << "AlongStepDoIt for "             << 788   if(-1 < verboseLevel) {
789            << GetProcessName() << " and partic << 789     const G4ParticleDefinition* d = track.GetDefinition();
790            << "  eScaled(MeV)=" << preStepScal << 790     G4cout << "AlongStepDoIt for "
791            << "  range(mm)=" << fRange/mm << " << 791            << GetProcessName() << " and particle "
792            << "  rf=" << reduceFactor << "  q^ << 792            << d->GetParticleName()
793            << " md=" << d->GetPDGMass() << "   << 793            << "  eScaled(MeV)= " << preStepScaledEnergy/MeV
794            << "  " << track.GetMaterial()->Get << 794            << "  range(mm)= " << fRange/mm
795   }                                            << 795            << "  s(mm)= " << length/mm
796   */                                           << 796            << "  q^2= " << chargeSqRatio
797   const G4DynamicParticle* dynParticle = track << 797            << " md= " << d->GetPDGMass()
798                                                << 798            << "  status= " << track.GetTrackStatus()
799   // define new weight for primary and seconda << 799            << G4endl;
800   G4double weight = fParticleChange.GetParentW << 800   }
801   if(weightFlag) {                             << 801   */
802     weight /= biasFactor;                      << 802 
803     fParticleChange.ProposeWeight(weight);     << 803   const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
804   }                                            << 804 
805                                                << 805   // stopping
806   // stopping, check actual range and kinetic  << 806   if (length >= fRange) {
807   if (length >= fRange || preStepKinEnergy <=  << 807     eloss = preStepKinEnergy;
808     eloss = preStepKinEnergy;                  << 808     currentModel->CorrectionsAlongStep(currentCouple, dynParticle, 
809     if (useDeexcitation) {                     << 809                eloss, esecdep, length);
810       atomDeexcitation->AlongStepDeexcitation( << 810     fParticleChange.SetProposedKineticEnergy(0.0);
811                                                << 811     fParticleChange.ProposeLocalEnergyDeposit(preStepKinEnergy);
812       if(scTracks.size() > 0) { FillSecondarie << 812     return &fParticleChange;
813       eloss = std::max(eloss, 0.0);            << 813   }
814     }                                          << 814 
815     fParticleChange.SetProposedKineticEnergy(0 << 815   // Short step
816     fParticleChange.ProposeLocalEnergyDeposit( << 816   eloss = GetDEDXForScaledEnergy(preStepScaledEnergy)*length;
817     return &fParticleChange;                   << 817 
818   }                                            << 818   // Long step
819   // zero step length with non-zero range      << 819   //} else {
820   if(length <= 0.0) { return &fParticleChange; << 820   if(eloss > preStepKinEnergy*linLossLimit) {
821                                                << 821 
822   // Short step                                << 822     G4double x = 
823   eloss = length*GetDEDXForScaledEnergy(preSte << 823       GetScaledRangeForScaledEnergy(preStepScaledEnergy) - length/reduceFactor;
824                                         LogSca << 824     eloss = preStepKinEnergy - ScaledKinEnergyForLoss(x)/massRatio;
825   /*                                           << 825    
826   G4cout << "##### Short STEP: eloss= " << elo << 826     /*
827    << " Escaled=" << preStepScaledEnergy       << 827     if(-1 < verboseLevel) 
828    << " R=" << fRange                          << 828       G4cout << "Long STEP: rPre(mm)= " 
829    << " L=" << length                          << 829              << GetScaledRangeForScaledEnergy(preStepScaledEnergy)/mm
830    << " fFactor=" << fFactor << " minE=" << mi << 830              << " rPost(mm)= " << x/mm
831    << " idxBase=" << basedCoupleIndex << G4end << 831              << " ePre(MeV)= " << preStepScaledEnergy/MeV
832   */                                           << 832              << " eloss(MeV)= " << eloss/MeV
833   // Long step                                 << 833              << " eloss0(MeV)= "
834   if(eloss > preStepKinEnergy*linLossLimit) {  << 834              << GetDEDXForScaledEnergy(preStepScaledEnergy)*length/MeV
835                                                << 835        << " lim(MeV)= " << preStepKinEnergy*linLossLimit/MeV
836     const G4double x = (fRange - length)/reduc << 836              << G4endl;
837     const G4double de = preStepKinEnergy - Sca << 837     */
838     if(de > 0.0) { eloss = de; }               << 838   }
839     /*                                         << 839 
840     if(-1 < verboseLevel)                      << 840   /*   
841       G4cout << "  Long STEP: rPre(mm)="       << 841   G4double eloss0 = eloss;
842              << GetScaledRangeForScaledEnergy( << 842   if(-1 < verboseLevel ) {
843              << " x(mm)=" << x/mm              << 843     G4cout << "Before fluct: eloss(MeV)= " << eloss/MeV
844              << " eloss(MeV)=" << eloss/MeV    << 844            << " e-eloss= " << preStepKinEnergy-eloss
845        << " rFactor=" << reduceFactor          << 845            << " step(mm)= " << length/mm
846        << " massRatio=" << massRatio           << 846            << " range(mm)= " << fRange/mm
847              << G4endl;                        << 847            << " fluct= " << lossFluctuationFlag
848     */                                         << 848            << G4endl;
849   }                                            << 849   }
850                                                << 850   */
851   /*                                           << 851 
852   if(-1 < verboseLevel ) {                     << 852   G4double cut  = (*theCuts)[currentMaterialIndex];
853     G4cout << "Before fluct: eloss(MeV)= " <<  << 853   G4double esec = 0.0;
854            << " e-eloss= " << preStepKinEnergy << 854 
855            << " step(mm)= " << length/mm << "  << 855   // SubCutOff 
856            << " fluct= " << lossFluctuationFla << 856   if(useSubCutoff) {
857   }                                            << 857     if(idxSCoffRegions[currentMaterialIndex]) {
858   */                                           << 858 
859                                                << 859       G4bool yes = false;
860   const G4double cut = (*theCuts)[currentCoupl << 860       G4StepPoint* prePoint  = step.GetPreStepPoint();
861   G4double esec = 0.0;                         << 861 
862                                                << 862       // Check boundary
863   // Corrections, which cannot be tabulated    << 863       if(prePoint->GetStepStatus() == fGeomBoundary) yes = true;
864   if(isIon) {                                  << 864 
865     currentModel->CorrectionsAlongStep(current << 865       // Check PrePoint
866                                        length, << 866       else {
867     eloss = std::max(eloss, 0.0);              << 867   G4double preSafety  = prePoint->GetSafety();
868   }                                            << 868   G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1);
869                                                << 869 
870   // Sample fluctuations if not full energy lo << 870   // recompute presafety
871   if(eloss >= preStepKinEnergy) {              << 871         if(preSafety < rcut) {
872     eloss = preStepKinEnergy;                  << 872     preSafety = safetyHelper->ComputeSafety(prePoint->GetPosition());
873                                                << 873   }
874   } else if (lossFluctuationFlag) {            << 874 
875     const G4double tmax = currentModel->MaxSec << 875         if(preSafety < rcut) yes = true;
876     const G4double tcut = std::min(cut, tmax); << 876 
877     G4VEmFluctuationModel* fluc = currentModel << 877   // Check PostPoint
878     eloss = fluc->SampleFluctuations(currentCo << 878   else {
879                                      tcut, tma << 879     G4double postSafety = preSafety - length; 
880     /*                                         << 880     if(postSafety < rcut) {
881     if(-1 < verboseLevel)                      << 881       postSafety = 
882       G4cout << "After fluct: eloss(MeV)= " << << 882         safetyHelper->ComputeSafety(step.GetPostStepPoint()->GetPosition());
883              << " fluc= " << (eloss-eloss0)/Me << 883       if(postSafety < rcut) yes = true;
884              << " ChargeSqRatio= " << chargeSq << 884     }
885              << " massRatio= " << massRatio << << 885   }
886     */                                         << 886       }
887   }                                            << 887   
888                                                << 888       // Decide to start subcut sampling
889   // deexcitation                              << 889       if(yes) {
890   if (useDeexcitation) {                       << 890 
891     G4double esecfluo = preStepKinEnergy;      << 891         cut = (*theSubCuts)[currentMaterialIndex];
892     G4double de = esecfluo;                    << 892   eloss -= GetSubDEDXForScaledEnergy(preStepScaledEnergy)*length;
893     atomDeexcitation->AlongStepDeexcitation(sc << 893   scTracks.clear();
894                                             de << 894   SampleSubCutSecondaries(scTracks, step, 
895                                                << 895         currentModel,currentMaterialIndex, 
896     // sum of de-excitation energies           << 896         esecdep);
897     esecfluo -= de;                            << 897   // add bremsstrahlung sampling
898                                                << 898   /*
899     // subtracted from energy loss             << 899   if(nProcesses > 0) {
900     if(eloss >= esecfluo) {                    << 900     for(G4int i=0; i<nProcesses; i++) {
901       esec  += esecfluo;                       << 901       (scProcesses[i])->SampleSubCutSecondaries(
902       eloss -= esecfluo;                       << 902     scTracks, step, (scProcesses[i])->
903     } else {                                   << 903     SelectModelForMaterial(preStepKinEnergy, currentMaterialIndex),
904       esec += esecfluo;                        << 904     currentMaterialIndex,esecdep);
905       eloss = 0.0;                             << 905     }
906     }                                          << 906   } 
907   }                                            << 907   */   
908   if(nullptr != subcutProducer && IsRegionForC << 908   G4int n = scTracks.size();
909     subcutProducer->SampleSecondaries(step, sc << 909   if(n>0) {
910   }                                            << 910     G4ThreeVector mom = dynParticle->GetMomentum();
911   // secondaries from atomic de-excitation and << 911     fParticleChange.SetNumberOfSecondaries(n);
912   if(!scTracks.empty()) { FillSecondariesAlong << 912     for(G4int i=0; i<n; i++) {
913                                                << 913       G4Track* t = scTracks[i];
914   // Energy balance                            << 914       G4double e = t->GetKineticEnergy();
915   G4double finalT = preStepKinEnergy - eloss - << 915       if (t->GetDefinition() == thePositron) e += 2.0*electron_mass_c2;
916   if (finalT <= lowestKinEnergy) {             << 916       esec += e;
917     eloss += finalT;                           << 917       pParticleChange->AddSecondary(t);
918     finalT = 0.0;                              << 918     }      
919   } else if(isIon) {                           << 919   }
920     fParticleChange.SetProposedCharge(         << 920       }
921       currentModel->GetParticleCharge(track.Ge << 921     }
922                                       currentM << 922   }
923   }                                            << 923 
924   eloss = std::max(eloss, 0.0);                << 924   // Corrections, which cannot be tabulated
925                                                << 925   currentModel->CorrectionsAlongStep(currentCouple, dynParticle, 
926   fParticleChange.SetProposedKineticEnergy(fin << 926              eloss, esecdep, length);
927   fParticleChange.ProposeLocalEnergyDeposit(el << 927 
928   /*                                           << 928   // Sample fluctuations
929   if(-1 < verboseLevel) {                      << 929   if (lossFluctuationFlag) {
930     G4double del = finalT + eloss + esec - pre << 930     G4VEmFluctuationModel* fluc = currentModel->GetModelOfFluctuations();
931     G4cout << "Final value eloss(MeV)= " << el << 931     if(fluc && 
932            << " preStepKinEnergy= " << preStep << 932       (eloss + esec + esecdep + lowestKinEnergy) < preStepKinEnergy) {
933            << " postStepKinEnergy= " << finalT << 933 
934            << " de(keV)= " << del/keV          << 934       G4double tmax = 
935            << " lossFlag= " << lossFluctuation << 935   std::min(currentModel->MaxSecondaryKinEnergy(dynParticle),cut);
936            << "  status= " << track.GetTrackSt << 936       eloss = fluc->SampleFluctuations(currentMaterial,dynParticle,
937            << G4endl;                          << 937                tmax,length,eloss);
938   }                                            << 938       /*                            
939   */                                           << 939       if(-1 < verboseLevel) 
940   return &fParticleChange;                     << 940       G4cout << "After fluct: eloss(MeV)= " << eloss/MeV
941 }                                              << 941              << " fluc= " << (eloss-eloss0)/MeV
942                                                << 942              << " ChargeSqRatio= " << chargeSqRatio
943 //....oooOO0OOooo........oooOO0OOooo........oo << 943              << " massRatio= " << massRatio
944                                                << 944              << " tmax= " << tmax
945 void G4VEnergyLossProcess::FillSecondariesAlon << 945              << G4endl;
946 {                                              << 946       */
947   const std::size_t n0 = scTracks.size();      << 947     }
948   G4double weight = wt;                        << 948   }
949   // weight may be changed by biasing manager  << 949   // add low-energy subcutoff particles
950   if(biasManager) {                            << 950   eloss += esecdep;
951     if(biasManager->SecondaryBiasingRegion((G4 << 951   if(eloss < 0.0) eloss = 0.0;
952       weight *=                                << 952 
953         biasManager->ApplySecondaryBiasing(scT << 953   // Energy balanse
954     }                                          << 954   G4double finalT = preStepKinEnergy - eloss - esec;
955   }                                            << 955   if (finalT <= lowestKinEnergy) {
956                                                << 956     eloss  = preStepKinEnergy - esec;
957   // fill secondaries                          << 957     finalT = 0.0;
958   const std::size_t n = scTracks.size();       << 958   } else if(isIon) {
959   fParticleChange.SetNumberOfSecondaries((G4in << 959     fParticleChange.SetProposedCharge(
960                                                << 960       currentModel->GetParticleCharge(track.GetDefinition(),currentMaterial,finalT));
961   for(std::size_t i=0; i<n; ++i) {             << 961   }
962     G4Track* t = scTracks[i];                  << 962 
963     if(nullptr != t) {                         << 963   fParticleChange.SetProposedKineticEnergy(finalT);
964       t->SetWeight(weight);                    << 964   fParticleChange.ProposeLocalEnergyDeposit(eloss);
965       pParticleChange->AddSecondary(t);        << 965 
966       G4int pdg = t->GetDefinition()->GetPDGEn << 966   /*  
967       if (i < n0) {                            << 967   if(-1 < verboseLevel) {
968         if (pdg == 22) {                       << 968     G4cout << "Final value eloss(MeV)= " << eloss/MeV
969     t->SetCreatorModelID(gpixeID);             << 969            << " preStepKinEnergy= " << preStepKinEnergy
970         } else if (pdg == 11) {                << 970            << " postStepKinEnergy= " << finalT
971           t->SetCreatorModelID(epixeID);       << 971            << " lossFlag= " << lossFluctuationFlag
972         } else {                               << 972            << "  status= " << track.GetTrackStatus()
973           t->SetCreatorModelID(biasID);        << 973            << G4endl;
974   }                                            << 974   }
975       } else {                                 << 975   */  
976   t->SetCreatorModelID(biasID);                << 976 
977       }                                        << 977   return &fParticleChange;
978     }                                          << 978 }
979   }                                            << 979 
980   scTracks.clear();                            << 980 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
981 }                                              << 981 
982                                                << 982 void G4VEnergyLossProcess::SampleSubCutSecondaries(
983 //....oooOO0OOooo........oooOO0OOooo........oo << 983        std::vector<G4Track*>& tracks, 
984                                                << 984        const G4Step& step, 
985 G4VParticleChange* G4VEnergyLossProcess::PostS << 985        G4VEmModel* model,
986                                                << 986        G4int idx,
987 {                                              << 987        G4double& /*extraEdep*/) 
988   // clear number of interaction lengths in an << 988 {
989   theNumberOfInteractionLengthLeft = -1.0;     << 989   // Fast check weather subcutoff can work
990   mfpKinEnergy = DBL_MAX;                      << 990   G4double subcut = (*theSubCuts)[idx];
991                                                << 991   G4double cut = (*theCuts)[idx];
992   fParticleChange.InitializeForPostStep(track) << 992   if(cut <= subcut) return;
993   const G4double finalT = track.GetKineticEner << 993 
994                                                << 994   const G4Track* track = step.GetTrack();
995   const G4double postStepScaledEnergy = finalT << 995   const G4DynamicParticle* dp = track->GetDynamicParticle();
996   SelectModel(postStepScaledEnergy);           << 996   G4double e = dp->GetKineticEnergy()*massRatio;
997                                                << 997   G4bool b;
998   if(!currentModel->IsActive(postStepScaledEne << 998   G4double cross = chargeSqRatio*(((*theSubLambdaTable)[idx])->GetValue(e,b));
999     return &fParticleChange;                   << 999   G4double length = step.GetStepLength();
1000   }                                           << 1000 
1001   /*                                          << 1001   // negligible probability to get any interaction
1002   if(1 < verboseLevel) {                      << 1002   if(length*cross < perMillion) return;
1003     G4cout<<GetProcessName()<<" PostStepDoIt: << 1003   /*      
1004   }                                           << 1004   if(-1 < verboseLevel) 
1005   */                                          << 1005     G4cout << "<<< Subcutoff for " << GetProcessName()
1006   // forced process - should happen only once << 1006      << " cross(1/mm)= " << cross*mm << ">>>"
1007   if(biasFlag) {                              << 1007      << " e(MeV)= " << preStepScaledEnergy
1008     if(biasManager->ForcedInteractionRegion(( << 1008      << " matIdx= " << currentMaterialIndex
1009       biasFlag = false;                       << 1009      << G4endl;
1010     }                                         << 1010   */
1011   }                                           << 1011 
1012   const G4DynamicParticle* dp = track.GetDyna << 1012   // Sample subcutoff secondaries
1013                                               << 1013   G4StepPoint* preStepPoint = step.GetPreStepPoint();
1014   // Integral approach                        << 1014   G4StepPoint* postStepPoint = step.GetPostStepPoint();
1015   if (fXSType != fEmNoIntegral) {             << 1015   G4ThreeVector prepoint = preStepPoint->GetPosition();
1016     const G4double logFinalT = dp->GetLogKine << 1016   G4ThreeVector dr = postStepPoint->GetPosition() - prepoint;
1017     G4double lx = GetLambdaForScaledEnergy(po << 1017   G4double pretime = preStepPoint->GetGlobalTime();
1018                                            lo << 1018   G4double dt = postStepPoint->GetGlobalTime() - pretime;
1019     lx = std::max(lx, 0.0);                   << 1019   //G4double dt = length/preStepPoint->GetVelocity();
1020                                               << 1020   G4double fragment = 0.0;
1021     // if both lg and lx are zero then no int << 1021 
1022     if(preStepLambda*G4UniformRand() >= lx) { << 1022   do {
1023       return &fParticleChange;                << 1023     G4double del = -std::log(G4UniformRand())/cross;
1024     }                                         << 1024     fragment += del/length;
1025   }                                           << 1025     if (fragment > 1.0) break;
1026                                               << 1026 
1027   // define new weight for primary and second << 1027     // sample secondaries
1028   G4double weight = fParticleChange.GetParent << 1028     secParticles.clear();
1029   if(weightFlag) {                            << 1029     model->SampleSecondaries(&secParticles,track->GetMaterialCutsCouple(),
1030     weight /= biasFactor;                     << 1030            dp,subcut,cut);
1031     fParticleChange.ProposeWeight(weight);    << 1031 
1032   }                                           << 1032     // position of subcutoff particles
1033                                               << 1033     G4ThreeVector r = prepoint + fragment*dr;
1034   const G4double tcut = (*theCuts)[currentCou << 1034     std::vector<G4DynamicParticle*>::iterator it;
1035                                               << 1035     for(it=secParticles.begin(); it!=secParticles.end(); it++) {
1036   // sample secondaries                       << 1036 
1037   secParticles.clear();                       << 1037       G4bool addSec = true;
1038   currentModel->SampleSecondaries(&secParticl << 1038       /*
1039                                               << 1039       // do not track very low-energy delta-electrons
1040   const G4int num0 = (G4int)secParticles.size << 1040       if(theSecondaryRangeTable && (*it)->GetDefinition() == theElectron) {
1041                                               << 1041   G4bool b;
1042   // bremsstrahlung splitting or Russian roul << 1042   G4double ekin = (*it)->GetKineticEnergy();
1043   if(biasManager) {                           << 1043   G4double rg = ((*theSecondaryRangeTable)[idx]->GetValue(ekin, b));
1044     if(biasManager->SecondaryBiasingRegion((G << 1044   //          if(rg < currentMinSafety) {
1045       G4double eloss = 0.0;                   << 1045   if(rg < safetyHelper->ComputeSafety(r)) {
1046       weight *= biasManager->ApplySecondaryBi << 1046     extraEdep += ekin;
1047                                       secPart << 1047     delete (*it);
1048                                       track,  << 1048     addSec = false;
1049                                       &fParti << 1049   }
1050                                       (G4int) << 1050       }
1051                                       step.Ge << 1051       */
1052       if(eloss > 0.0) {                       << 1052       if(addSec) {
1053         eloss += fParticleChange.GetLocalEner << 1053   G4Track* t = new G4Track((*it), pretime + fragment*dt, r);
1054         fParticleChange.ProposeLocalEnergyDep << 1054   //G4Track* t = new G4Track((*it), pretime, r);
1055       }                                       << 1055   t->SetTouchableHandle(track->GetTouchableHandle());
1056     }                                         << 1056   tracks.push_back(t);
1057   }                                           << 1057 
1058                                               << 1058   /*  
1059   // save secondaries                         << 1059   if(-1 < verboseLevel) 
1060   const G4int num = (G4int)secParticles.size( << 1060     G4cout << "New track " << t->GetDefinition()->GetParticleName()
1061   if(num > 0) {                               << 1061      << " e(keV)= " << t->GetKineticEnergy()/keV
1062                                               << 1062      << " fragment= " << fragment
1063     fParticleChange.SetNumberOfSecondaries(nu << 1063      << G4endl;
1064     G4double time = track.GetGlobalTime();    << 1064   */
1065                                               << 1065       }
1066     G4int n1(0), n2(0);                       << 1066     }
1067     if(num0 > mainSecondaries) {              << 1067   } while (fragment <= 1.0);
1068       currentModel->FillNumberOfSecondaries(n << 1068 } 
1069     }                                         << 1069 
1070                                               << 1070 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1071     for (G4int i=0; i<num; ++i) {             << 1071 
1072       if(nullptr != secParticles[i]) {        << 1072 G4VParticleChange* G4VEnergyLossProcess::PostStepDoIt(const G4Track& track,
1073         G4Track* t = new G4Track(secParticles << 1073                                                       const G4Step&)
1074         t->SetTouchableHandle(track.GetToucha << 1074 {
1075         if (biasManager) {                    << 1075   fParticleChange.InitializeForPostStep(track);
1076           t->SetWeight(weight * biasManager-> << 1076   G4double finalT = track.GetKineticEnergy();
1077         } else {                              << 1077   if(finalT <= lowestKinEnergy) return &fParticleChange;
1078           t->SetWeight(weight);               << 1078 
1079         }                                     << 1079   G4double postStepScaledEnergy = finalT*massRatio;
1080         if(i < num0) {                        << 1080   /*
1081           t->SetCreatorModelID(secID);        << 1081   if(-1 < verboseLevel) {
1082         } else if(i < num0 + n1) {            << 1082     G4cout << GetProcessName()
1083           t->SetCreatorModelID(tripletID);    << 1083            << "::PostStepDoIt: E(MeV)= " << finalT/MeV
1084         } else {                              << 1084      << G4endl;
1085           t->SetCreatorModelID(biasID);       << 1085   }
1086         }                                     << 1086   */
1087                                               << 1087   // Integral approach
1088         //G4cout << "Secondary(post step) has << 1088   if (integral) {
1089         //       << ", kenergy " << t->GetKin << 1089     G4double lx = GetLambdaForScaledEnergy(postStepScaledEnergy);
1090         //       << " time= " << time/ns << " << 1090     /*
1091         pParticleChange->AddSecondary(t);     << 1091     if(preStepLambda<lx && 1 < verboseLevel && nWarnings<200) {
1092       }                                       << 1092       G4cout << "WARNING: for " << particle->GetParticleName()
1093     }                                         << 1093              << " and " << GetProcessName()
1094   }                                           << 1094              << " E(MeV)= " << finalT/MeV
1095                                               << 1095              << " preLambda= " << preStepLambda 
1096   if(0.0 == fParticleChange.GetProposedKineti << 1096        << " < " << lx << " (postLambda) "
1097      fAlive == fParticleChange.GetTrackStatus << 1097        << G4endl;
1098     if(particle->GetProcessManager()->GetAtRe << 1098       nWarnings++;
1099          { fParticleChange.ProposeTrackStatus << 1099     }
1100     else { fParticleChange.ProposeTrackStatus << 1100     */
1101   }                                           << 1101     if(preStepLambda*G4UniformRand() > lx) {
1102                                               << 1102       ClearNumberOfInteractionLengthLeft();
1103   /*                                          << 1103       return &fParticleChange;
1104   if(-1 < verboseLevel) {                     << 1104     }
1105     G4cout << "::PostStepDoIt: Sample seconda << 1105   }
1106     << fParticleChange.GetProposedKineticEner << 1106 
1107            << " MeV; model= (" << currentMode << 1107   SelectModel(postStepScaledEnergy);
1108            << ", " <<  currentModel->HighEner << 1108   const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
1109            << "  preStepLambda= " << preStepL << 1109   G4double tcut = (*theCuts)[currentMaterialIndex];
1110            << "  dir= " << track.GetMomentumD << 1110 
1111            << "  status= " << track.GetTrackS << 1111   // sample secondaries
1112            << G4endl;                         << 1112   secParticles.clear();
1113   }                                           << 1113   currentModel->SampleSecondaries(&secParticles, currentCouple, dynParticle, tcut);
1114   */                                          << 1114 
1115   return &fParticleChange;                    << 1115   // save secondaries
1116 }                                             << 1116   G4int num = secParticles.size();
1117                                               << 1117   if(num > 0) {
1118 //....oooOO0OOooo........oooOO0OOooo........o << 1118     fParticleChange.SetNumberOfSecondaries(num);
1119                                               << 1119     for (G4int i=0; i<num; i++) {
1120 G4bool G4VEnergyLossProcess::StorePhysicsTabl << 1120       fParticleChange.AddSecondary(secParticles[i]);
1121        const G4ParticleDefinition* part, cons << 1121     }
1122 {                                             << 1122   }
1123   if (!isMaster || nullptr != baseParticle || << 1123 
1124   for(std::size_t i=0; i<7; ++i) {            << 1124   /*
1125     // ionisation table only for ionisation p << 1125   if(-1 < verboseLevel) {
1126     if (nullptr == theData->Table(i) || (!isI << 1126     G4cout << "::PostStepDoIt: Sample secondary; Efin= " 
1127       continue;                               << 1127     << fParticleChange.GetProposedKineticEnergy()/MeV
1128     }                                         << 1128            << " MeV; model= (" << currentModel->LowEnergyLimit()
1129     if (-1 < verboseLevel) {                  << 1129            << ", " <<  currentModel->HighEnergyLimit() << ")"
1130       G4cout << "G4VEnergyLossProcess::StoreP << 1130            << "  preStepLambda= " << preStepLambda
1131        << "  " << particle->GetParticleName() << 1131            << "  dir= " << track.GetMomentumDirection()
1132        << "  " << GetProcessName()            << 1132            << "  status= " << track.GetTrackStatus()
1133        << "  " << tnames[i] << "  " << theDat << 1133            << G4endl;
1134     }                                         << 1134   }
1135     if (!G4EmTableUtil::StoreTable(this, part << 1135   */
1136            dir, tnames[i], verboseLevel, asci << 1136   ClearNumberOfInteractionLengthLeft();
1137       return false;                           << 1137   return &fParticleChange;
1138     }                                         << 1138 }
1139   }                                           << 1139 
1140   return true;                                << 1140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1141 }                                             << 1141 
1142                                               << 1142 void G4VEnergyLossProcess::PrintInfoDefinition()
1143 //....oooOO0OOooo........oooOO0OOooo........o << 1143 {
1144                                               << 1144   if(0 < verboseLevel) {
1145 G4bool                                        << 1145     G4cout << G4endl << GetProcessName() << ":   for  "
1146 G4VEnergyLossProcess::RetrievePhysicsTable(co << 1146            << particle->GetParticleName()
1147                                            co << 1147      << "    SubType= " << GetProcessSubType() 
1148 {                                             << 1148            << G4endl
1149   if (!isMaster || nullptr != baseParticle || << 1149            << "      dE/dx and range tables from "
1150   for(std::size_t i=0; i<7; ++i) {            << 1150      << G4BestUnit(minKinEnergy,"Energy")
1151     // ionisation table only for ionisation p << 1151            << " to " << G4BestUnit(maxKinEnergy,"Energy")
1152     if (!isIonisation && 1 == i) { continue;  << 1152            << " in " << nBins << " bins" << G4endl
1153     if(!G4EmTableUtil::RetrieveTable(this, pa << 1153            << "      Lambda tables from threshold to "
1154                                      verboseL << 1154            << G4BestUnit(maxKinEnergy,"Energy")
1155       return false;                           << 1155            << " in " << nBins << " bins, spline: " 
1156     }                                         << 1156      << (G4LossTableManager::Instance())->SplineFlag()
1157   }                                           << 1157            << G4endl;
1158   return true;                                << 1158     if(theRangeTableForLoss && isIonisation) {
1159 }                                             << 1159       G4cout << "      finalRange(mm)= " << finalRange/mm
1160                                               << 1160              << ", dRoverRange= " << dRoverRange
1161 //....oooOO0OOooo........oooOO0OOooo........o << 1161              << ", integral: " << integral
1162                                               << 1162              << ", fluct: " << lossFluctuationFlag
1163 G4double G4VEnergyLossProcess::GetDEDXDispers << 1163        << ", linLossLimit= " << linLossLimit
1164                                   const G4Mat << 1164              << G4endl;
1165                                   const G4Dyn << 1165     }
1166                                         G4dou << 1166     PrintInfo();
1167 {                                             << 1167     modelManager->DumpModelList(verboseLevel);
1168   DefineMaterial(couple);                     << 1168     if(theCSDARangeTable && isIonisation) {
1169   G4double ekin = dp->GetKineticEnergy();     << 1169       G4cout << "      CSDA range table up"
1170   SelectModel(ekin*massRatio);                << 1170              << " to " << G4BestUnit(maxKinEnergyCSDA,"Energy")
1171   G4double tmax = currentModel->MaxSecondaryK << 1171              << " in " << nBinsCSDA << " bins" << G4endl;
1172   G4double tcut = std::min(tmax,(*theCuts)[cu << 1172     }
1173   G4double d = 0.0;                           << 1173     if(nSCoffRegions>0 && isIonisation) {
1174   G4VEmFluctuationModel* fm = currentModel->G << 1174       G4cout << "      Subcutoff sampling in " << nSCoffRegions 
1175   if(nullptr != fm) { d = fm->Dispersion(curr << 1175        << " regions" << G4endl;
1176   return d;                                   << 1176     }
1177 }                                             << 1177     if(2 < verboseLevel) {
1178                                               << 1178       G4cout << "      DEDXTable address= " << theDEDXTable << G4endl;
1179 //....oooOO0OOooo........oooOO0OOooo........o << 1179       if(theDEDXTable && isIonisation) G4cout << (*theDEDXTable) << G4endl;
1180                                               << 1180       G4cout << "non restricted DEDXTable address= " 
1181 G4double                                      << 1181        << theDEDXunRestrictedTable << G4endl;
1182 G4VEnergyLossProcess::CrossSectionPerVolume(G << 1182       if(theDEDXunRestrictedTable && isIonisation) {
1183                                             c << 1183            G4cout << (*theDEDXunRestrictedTable) << G4endl;
1184                                             G << 1184       }
1185 {                                             << 1185       if(theDEDXSubTable && isIonisation) {
1186   // Cross section per volume is calculated   << 1186   G4cout << (*theDEDXSubTable) << G4endl;
1187   DefineMaterial(couple);                     << 1187       }
1188   G4double cross = 0.0;                       << 1188       G4cout << "      CSDARangeTable address= " << theCSDARangeTable 
1189   if (nullptr != theLambdaTable) {            << 1189        << G4endl;
1190     cross = GetLambdaForScaledEnergy(kineticE << 1190       if(theCSDARangeTable && isIonisation) {
1191                                      logKinet << 1191   G4cout << (*theCSDARangeTable) << G4endl;
1192   } else {                                    << 1192       }
1193     SelectModel(kineticEnergy*massRatio);     << 1193       G4cout << "      RangeTableForLoss address= " << theRangeTableForLoss 
1194     cross = (!baseMat) ? biasFactor : biasFac << 1194        << G4endl;
1195     cross *= (currentModel->CrossSectionPerVo << 1195       if(theRangeTableForLoss && isIonisation) {
1196                                               << 1196              G4cout << (*theRangeTableForLoss) << G4endl;
1197   }                                           << 1197       }
1198   return std::max(cross, 0.0);                << 1198       G4cout << "      InverseRangeTable address= " << theInverseRangeTable 
1199 }                                             << 1199        << G4endl;
1200                                               << 1200       if(theInverseRangeTable && isIonisation) {
1201 //....oooOO0OOooo........oooOO0OOooo........o << 1201              G4cout << (*theInverseRangeTable) << G4endl;
1202                                               << 1202       }
1203 G4double G4VEnergyLossProcess::MeanFreePath(c << 1203       G4cout << "      LambdaTable address= " << theLambdaTable << G4endl;
1204 {                                             << 1204       if(theLambdaTable && isIonisation) {
1205   DefineMaterial(track.GetMaterialCutsCouple( << 1205   G4cout << (*theLambdaTable) << G4endl;
1206   const G4double kinEnergy    = track.GetKine << 1206       }
1207   const G4double logKinEnergy = track.GetDyna << 1207       G4cout << "      SubLambdaTable address= " << theSubLambdaTable << G4endl;
1208   const G4double cs = GetLambdaForScaledEnerg << 1208       if(theSubLambdaTable && isIonisation) {
1209                                               << 1209   G4cout << (*theSubLambdaTable) << G4endl;
1210   return (0.0 < cs) ? 1.0/cs : DBL_MAX;       << 1210       }
1211 }                                             << 1211     }
1212                                               << 1212   }
1213 //....oooOO0OOooo........oooOO0OOooo........o << 1213 }
1214                                               << 1214 
1215 G4double G4VEnergyLossProcess::ContinuousStep << 1215 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1216                                               << 1216 
1217                                               << 1217 void G4VEnergyLossProcess::SetDEDXTable(G4PhysicsTable* p, G4EmTableType tType)
1218 {                                             << 1218 {
1219   return AlongStepGetPhysicalInteractionLengt << 1219   if(fTotal == tType && theDEDXunRestrictedTable != p) {
1220 }                                             << 1220     if(theDEDXunRestrictedTable) theDEDXunRestrictedTable->clearAndDestroy();
1221                                               << 1221     theDEDXunRestrictedTable = p;
1222 //....oooOO0OOooo........oooOO0OOooo........o << 1222     if(p) {
1223                                               << 1223       size_t n = p->length();
1224 G4double G4VEnergyLossProcess::GetMeanFreePat << 1224       G4PhysicsVector* pv = (*p)[0];
1225                              const G4Track& t << 1225       G4double emax = maxKinEnergyCSDA;
1226                              G4double,        << 1226       G4bool b;
1227                              G4ForceCondition << 1227       theDEDXAtMaxEnergy = new G4double [n];
1228                                               << 1228 
1229 {                                             << 1229       for (size_t i=0; i<n; i++) {
1230   *condition = NotForced;                     << 1230   pv = (*p)[i];
1231   return MeanFreePath(track);                 << 1231   G4double dedx = pv->GetValue(emax, b);
1232 }                                             << 1232   theDEDXAtMaxEnergy[i] = dedx;
1233                                               << 1233   //G4cout << "i= " << i << " emax(MeV)= " << emax/MeV<< " dedx= " 
1234 //....oooOO0OOooo........oooOO0OOooo........o << 1234   //<< dedx << G4endl;
1235                                               << 1235       }
1236 G4double G4VEnergyLossProcess::GetContinuousS << 1236     }
1237                 const G4Track&,               << 1237 
1238                 G4double, G4double, G4double& << 1238   } else if(fRestricted == tType) {
1239 {                                             << 1239     theDEDXTable = p;
1240   return DBL_MAX;                             << 1240   } else if(fSubRestricted == tType) {    
1241 }                                             << 1241     theDEDXSubTable = p;
1242                                               << 1242   } else if(fIsIonisation == tType && theIonisationTable != p) {    
1243 //....oooOO0OOooo........oooOO0OOooo........o << 1243     if(theIonisationTable) theIonisationTable->clearAndDestroy();
1244                                               << 1244     theIonisationTable = p;
1245 G4PhysicsVector*                              << 1245   } else if(fIsSubIonisation == tType && theIonisationSubTable != p) {    
1246 G4VEnergyLossProcess::LambdaPhysicsVector(con << 1246     if(theIonisationSubTable) theIonisationSubTable->clearAndDestroy();
1247                                           G4d << 1247     theIonisationSubTable = p;
1248 {                                             << 1248   }
1249   DefineMaterial(couple);                     << 1249 }
1250   G4PhysicsVector* v = (*theLambdaTable)[base << 1250 
1251   return new G4PhysicsVector(*v);             << 1251 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1252 }                                             << 1252 
1253                                               << 1253 void G4VEnergyLossProcess::SetCSDARangeTable(G4PhysicsTable* p)
1254 //....oooOO0OOooo........oooOO0OOooo........o << 1254 {
1255                                               << 1255   if(theCSDARangeTable != p) theCSDARangeTable = p;
1256 void                                          << 1256 
1257 G4VEnergyLossProcess::SetDEDXTable(G4PhysicsT << 1257   if(p) {
1258 {                                             << 1258     size_t n = p->length();
1259   if(1 < verboseLevel) {                      << 1259     G4PhysicsVector* pv = (*p)[0];
1260     G4cout << "### Set DEDX table " << p << " << 1260     G4double emax = maxKinEnergyCSDA;
1261      << "  " <<  theDEDXunRestrictedTable <<  << 1261     G4bool b;
1262            << " for " << particle->GetParticl << 1262     theRangeAtMaxEnergy = new G4double [n];
1263            << " and process " << GetProcessNa << 1263 
1264      << " type=" << tType << " isIonisation:" << 1264     for (size_t i=0; i<n; i++) {
1265   }                                           << 1265       pv = (*p)[i];
1266   if(fTotal == tType) {                       << 1266       G4double r2 = pv->GetValue(emax, b);
1267     theDEDXunRestrictedTable = p;             << 1267       theRangeAtMaxEnergy[i] = r2;
1268   } else if(fRestricted == tType) {           << 1268       //G4cout << "i= " << i << " e2(MeV)= " << emax/MeV << " r2= " 
1269     theDEDXTable = p;                         << 1269       //<< r2<< G4endl;
1270     if(isMaster && nullptr == baseParticle) { << 1270     }
1271       theData->UpdateTable(theDEDXTable, 0);  << 1271   }
1272     }                                         << 1272 }
1273   } else if(fIsIonisation == tType) {         << 1273 
1274     theIonisationTable = p;                   << 1274 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1275     if(isMaster && nullptr == baseParticle) { << 1275 
1276       theData->UpdateTable(theIonisationTable << 1276 void G4VEnergyLossProcess::SetRangeTableForLoss(G4PhysicsTable* p)
1277     }                                         << 1277 {
1278   }                                           << 1278   if(theRangeTableForLoss != p) {
1279 }                                             << 1279     theRangeTableForLoss = p;
1280                                               << 1280     if(1 < verboseLevel) {
1281 //....oooOO0OOooo........oooOO0OOooo........o << 1281       G4cout << "### Set Range table " << p 
1282                                               << 1282        << " for " << particle->GetParticleName()
1283 void G4VEnergyLossProcess::SetCSDARangeTable( << 1283              << " and process " << GetProcessName() << G4endl;
1284 {                                             << 1284     }
1285   theCSDARangeTable = p;                      << 1285   }
1286 }                                             << 1286 }
1287                                               << 1287 
1288 //....oooOO0OOooo........oooOO0OOooo........o << 1288 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1289                                               << 1289 
1290 void G4VEnergyLossProcess::SetRangeTableForLo << 1290 void G4VEnergyLossProcess::SetSecondaryRangeTable(G4PhysicsTable* p)
1291 {                                             << 1291 {
1292   theRangeTableForLoss = p;                   << 1292   if(theSecondaryRangeTable != p) {
1293 }                                             << 1293     theSecondaryRangeTable = p;
1294                                               << 1294     if(1 < verboseLevel) {
1295 //....oooOO0OOooo........oooOO0OOooo........o << 1295       G4cout << "### Set SecondaryRange table " << p 
1296                                               << 1296        << " for " << particle->GetParticleName()
1297 void G4VEnergyLossProcess::SetInverseRangeTab << 1297              << " and process " << GetProcessName() << G4endl;
1298 {                                             << 1298     }
1299   theInverseRangeTable = p;                   << 1299   }
1300 }                                             << 1300 }
1301                                               << 1301 
1302 //....oooOO0OOooo........oooOO0OOooo........o << 1302 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1303                                               << 1303 
1304 void G4VEnergyLossProcess::SetLambdaTable(G4P << 1304 void G4VEnergyLossProcess::SetInverseRangeTable(G4PhysicsTable* p)
1305 {                                             << 1305 {
1306   if(1 < verboseLevel) {                      << 1306   if(theInverseRangeTable != p) {
1307     G4cout << "### Set Lambda table " << p << << 1307     theInverseRangeTable = p;
1308            << " for " << particle->GetParticl << 1308     if(1 < verboseLevel) {
1309            << " and process " << GetProcessNa << 1309       G4cout << "### Set InverseRange table " << p 
1310   }                                           << 1310        << " for " << particle->GetParticleName()
1311   theLambdaTable = p;                         << 1311              << " and process " << GetProcessName() << G4endl;
1312   tablesAreBuilt = true;                      << 1312     }
1313                                               << 1313   }
1314   if(isMaster && nullptr != p) {              << 1314 }
1315     delete theEnergyOfCrossSectionMax;        << 1315 
1316     theEnergyOfCrossSectionMax = nullptr;     << 1316 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1317     if(fEmTwoPeaks == fXSType) {              << 1317 
1318       if(nullptr != fXSpeaks) {               << 1318 void G4VEnergyLossProcess::SetLambdaTable(G4PhysicsTable* p)
1319   for(auto & ptr : *fXSpeaks) { delete ptr; } << 1319 {
1320   delete fXSpeaks;                            << 1320   if(1 < verboseLevel) {
1321       }                                       << 1321     G4cout << "### Set Lambda table " << p 
1322       G4LossTableBuilder* bld = lManager->Get << 1322      << " for " << particle->GetParticleName()
1323       fXSpeaks = G4EmUtility::FillPeaksStruct << 1323            << " and process " << GetProcessName() << G4endl;
1324       if(nullptr == fXSpeaks) { fXSType = fEm << 1324   }
1325     }                                         << 1325   if(theLambdaTable != p) theLambdaTable = p;
1326     if(fXSType == fEmOnePeak) {               << 1326   tablesAreBuilt = true;
1327       theEnergyOfCrossSectionMax = G4EmUtilit << 1327 
1328       if(nullptr == theEnergyOfCrossSectionMa << 1328   if(p) {
1329     }                                         << 1329     size_t n = p->length();
1330   }                                           << 1330     G4PhysicsVector* pv = (*p)[0];
1331 }                                             << 1331     G4double e, s, smax, emax;
1332                                               << 1332     theEnergyOfCrossSectionMax = new G4double [n];
1333 //....oooOO0OOooo........oooOO0OOooo........o << 1333     theCrossSectionMax = new G4double [n];
1334                                               << 1334     G4bool b;
1335 void G4VEnergyLossProcess::SetEnergyOfCrossSe << 1335 
1336 {                                             << 1336     for (size_t i=0; i<n; i++) {
1337   theEnergyOfCrossSectionMax = p;             << 1337       pv = (*p)[i];
1338 }                                             << 1338       emax = DBL_MAX;
1339                                               << 1339       smax = 0.0;
1340 //....oooOO0OOooo........oooOO0OOooo........o << 1340       if(pv) {
1341                                               << 1341         size_t nb = pv->GetVectorLength();
1342 void G4VEnergyLossProcess::SetTwoPeaksXS(std: << 1342         emax = pv->GetLowEdgeEnergy(nb);
1343 {                                             << 1343   for (size_t j=0; j<nb; j++) {
1344   fXSpeaks = ptr;                             << 1344     e = pv->GetLowEdgeEnergy(j);
1345 }                                             << 1345     s = pv->GetValue(e,b);
1346                                               << 1346     if(s > smax) {
1347 //....oooOO0OOooo........oooOO0OOooo........o << 1347       smax = s;
1348                                               << 1348       emax = e;
1349 const G4Element* G4VEnergyLossProcess::GetCur << 1349     }
1350 {                                             << 1350   }
1351   return (nullptr != currentModel)            << 1351       }
1352     ? currentModel->GetCurrentElement(current << 1352       theEnergyOfCrossSectionMax[i] = emax;
1353 }                                             << 1353       theCrossSectionMax[i] = smax;
1354                                               << 1354       if(1 < verboseLevel) {
1355 //....oooOO0OOooo........oooOO0OOooo........o << 1355         G4cout << "For " << particle->GetParticleName()
1356                                               << 1356                << " Max CS at i= " << i << " emax(MeV)= " << emax/MeV
1357 void G4VEnergyLossProcess::SetCrossSectionBia << 1357                << " lambda= " << smax << G4endl;
1358                                               << 1358       }
1359 {                                             << 1359     }
1360   if(f > 0.0) {                               << 1360   }
1361     biasFactor = f;                           << 1361 }
1362     weightFlag = flag;                        << 1362 
1363     if(1 < verboseLevel) {                    << 1363 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1364       G4cout << "### SetCrossSectionBiasingFa << 1364 
1365              << " process " << GetProcessName << 1365 void G4VEnergyLossProcess::SetSubLambdaTable(G4PhysicsTable* p)
1366              << " biasFactor= " << f << " wei << 1366 {
1367              << G4endl;                       << 1367   if(theSubLambdaTable != p) {
1368     }                                         << 1368     theSubLambdaTable = p;
1369   }                                           << 1369     if(1 < verboseLevel) {
1370 }                                             << 1370       G4cout << "### Set SebLambda table " << p 
1371                                               << 1371        << " for " << particle->GetParticleName()
1372 //....oooOO0OOooo........oooOO0OOooo........o << 1372              << " and process " << GetProcessName() << G4endl;
1373                                               << 1373     }
1374 void G4VEnergyLossProcess::ActivateForcedInte << 1374   }
1375                                               << 1375 }
1376                                               << 1376 
1377 {                                             << 1377 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1378   if(nullptr == biasManager) { biasManager =  << 1378 
1379   if(1 < verboseLevel) {                      << 1379 G4PhysicsVector* G4VEnergyLossProcess::LambdaPhysicsVector(
1380     G4cout << "### ActivateForcedInteraction: << 1380                  const G4MaterialCutsCouple* couple, G4double cut)
1381            << " process " << GetProcessName() << 1381 {
1382            << " length(mm)= " << length/mm    << 1382   //  G4double cut  = (*theCuts)[couple->GetIndex()];
1383            << " in G4Region <" << region      << 1383   //  G4int nbins = nLambdaBins;
1384            << "> weightFlag= " << flag        << 1384   G4double tmin = 
1385            << G4endl;                         << 1385     std::max(MinPrimaryEnergy(particle, couple->GetMaterial(), cut),
1386   }                                           << 1386        minKinEnergy);
1387   weightFlag = flag;                          << 1387   if(tmin >= maxKinEnergy) tmin = 0.5*maxKinEnergy;
1388   biasManager->ActivateForcedInteraction(leng << 1388   G4PhysicsVector* v = new G4PhysicsLogVector(tmin, maxKinEnergy, nBins);
1389 }                                             << 1389   v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
1390                                               << 1390 
1391 //....oooOO0OOooo........oooOO0OOooo........o << 1391   return v;
1392                                               << 1392 }
1393 void                                          << 1393 
1394 G4VEnergyLossProcess::ActivateSecondaryBiasin << 1394 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1395                                               << 1395 
1396                                               << 1396 G4double G4VEnergyLossProcess::CrossSectionPerVolume(
1397 {                                             << 1397    G4double kineticEnergy, const G4MaterialCutsCouple* couple)
1398   if (0.0 <= factor) {                        << 1398 {
1399     // Range cut can be applied only for e-   << 1399   // Cross section per volume is calculated
1400     if(0.0 == factor && secondaryParticle !=  << 1400   DefineMaterial(couple);
1401       { return; }                             << 1401   G4double cross = 0.0;
1402                                               << 1402   G4bool b;
1403     if(nullptr == biasManager) { biasManager  << 1403   if(theLambdaTable) {
1404     biasManager->ActivateSecondaryBiasing(reg << 1404     cross = 
1405     if(1 < verboseLevel) {                    << 1405       ((*theLambdaTable)[currentMaterialIndex])->GetValue(kineticEnergy, b);
1406       G4cout << "### ActivateSecondaryBiasing << 1406   } else {
1407              << " process " << GetProcessName << 1407     SelectModel(kineticEnergy);
1408              << " factor= " << factor         << 1408     cross = 
1409              << " in G4Region <" << region    << 1409       currentModel->CrossSectionPerVolume(currentMaterial,
1410              << "> energyLimit(MeV)= " << ene << 1410             particle, kineticEnergy,
1411              << G4endl;                       << 1411             (*theCuts)[currentMaterialIndex]);
1412     }                                         << 1412   }
1413   }                                           << 1413   if(cross < 0.0) { cross = 0.0; }
1414 }                                             << 1414 
1415                                               << 1415   return cross;
1416 //....oooOO0OOooo........oooOO0OOooo........o << 1416 }
1417                                               << 1417 
1418 void G4VEnergyLossProcess::SetIonisation(G4bo << 1418 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1419 {                                             << 1419 
1420   isIonisation = val;                         << 1420 G4bool G4VEnergyLossProcess::StorePhysicsTable(
1421   aGPILSelection = (val) ? CandidateForSelect << 1421        const G4ParticleDefinition* part, const G4String& directory, 
1422 }                                             << 1422        G4bool ascii)
1423                                               << 1423 {
1424 //....oooOO0OOooo........oooOO0OOooo........o << 1424   G4bool res = true;
1425                                               << 1425   if ( baseParticle || part != particle ) return res;
1426  void G4VEnergyLossProcess::SetLinearLossLimi << 1426 
1427 {                                             << 1427   if(!StoreTable(part,theDEDXTable,ascii,directory,"DEDX")) 
1428   if(0.0 < val && val < 1.0) {                << 1428     {res = false;}
1429     linLossLimit = val;                       << 1429 
1430     actLinLossLimit = true;                   << 1430   if(!StoreTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr")) 
1431   } else { PrintWarning("SetLinearLossLimit", << 1431     {res = false;}
1432 }                                             << 1432 
1433                                               << 1433   if(!StoreTable(part,theDEDXSubTable,ascii,directory,"SubDEDX")) 
1434 //....oooOO0OOooo........oooOO0OOooo........o << 1434     {res = false;}
1435                                               << 1435 
1436 void G4VEnergyLossProcess::SetStepFunction(G4 << 1436   if(!StoreTable(part,theIonisationTable,ascii,directory,"Ionisation")) 
1437 {                                             << 1437     {res = false;}
1438   if(0.0 < v1 && 0.0 < v2) {                  << 1438 
1439     dRoverRange = std::min(1.0, v1);          << 1439   if(!StoreTable(part,theIonisationSubTable,ascii,directory,"SubIonisation")) 
1440     finalRange = std::min(v2, 1.e+50);        << 1440     {res = false;}
1441   } else {                                    << 1441 
1442     PrintWarning("SetStepFunctionV1", v1);    << 1442   if(isIonisation &&
1443     PrintWarning("SetStepFunctionV2", v2);    << 1443      !StoreTable(part,theCSDARangeTable,ascii,directory,"CSDARange")) 
1444   }                                           << 1444     {res = false;}
1445 }                                             << 1445 
1446                                               << 1446   if(isIonisation &&
1447 //....oooOO0OOooo........oooOO0OOooo........o << 1447      !StoreTable(part,theRangeTableForLoss,ascii,directory,"Range")) 
1448                                               << 1448     {res = false;}
1449 void G4VEnergyLossProcess::SetLowestEnergyLim << 1449   
1450 {                                             << 1450   if(isIonisation &&
1451   if(1.e-18 < val && val < 1.e+50) { lowestKi << 1451      !StoreTable(part,theInverseRangeTable,ascii,directory,"InverseRange")) 
1452   else { PrintWarning("SetLowestEnergyLimit", << 1452     {res = false;}
1453 }                                             << 1453   
1454                                               << 1454   if(!StoreTable(part,theLambdaTable,ascii,directory,"Lambda")) 
1455 //....oooOO0OOooo........oooOO0OOooo........o << 1455     {res = false;}
1456                                               << 1456 
1457 void G4VEnergyLossProcess::SetDEDXBinning(G4i << 1457   if(!StoreTable(part,theSubLambdaTable,ascii,directory,"SubLambda")) 
1458 {                                             << 1458     {res = false;}
1459   if(2 < n && n < 1000000000) {               << 1459 
1460     nBins = n;                                << 1460   if ( res ) {
1461     actBinning = true;                        << 1461     if(0 < verboseLevel) {
1462   } else {                                    << 1462       G4cout << "Physics tables are stored for " << particle->GetParticleName()
1463     G4double e = (G4double)n;                 << 1463              << " and process " << GetProcessName()
1464     PrintWarning("SetDEDXBinning", e);        << 1464        << " in the directory <" << directory
1465   }                                           << 1465        << "> " << G4endl;
1466 }                                             << 1466     }
1467                                               << 1467   } else {
1468 //....oooOO0OOooo........oooOO0OOooo........o << 1468     G4cout << "Fail to store Physics Tables for " 
1469                                               << 1469      << particle->GetParticleName()
1470 void G4VEnergyLossProcess::SetMinKinEnergy(G4 << 1470            << " and process " << GetProcessName()
1471 {                                             << 1471      << " in the directory <" << directory
1472   if(1.e-18 < e && e < maxKinEnergy) {        << 1472      << "> " << G4endl;
1473     minKinEnergy = e;                         << 1473   }
1474     actMinKinEnergy = true;                   << 1474   return res;
1475   } else { PrintWarning("SetMinKinEnergy", e) << 1475 }
1476 }                                             << 1476 
1477                                               << 1477 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
1478 //....oooOO0OOooo........oooOO0OOooo........o << 1478 
1479                                               << 1479 G4bool G4VEnergyLossProcess::StoreTable(const G4ParticleDefinition* part, 
1480 void G4VEnergyLossProcess::SetMaxKinEnergy(G4 << 1480           G4PhysicsTable* aTable, G4bool ascii,
1481 {                                             << 1481           const G4String& directory,
1482   if(minKinEnergy < e && e < 1.e+50) {        << 1482           const G4String& tname)
1483     maxKinEnergy = e;                         << 1483 {
1484     actMaxKinEnergy = true;                   << 1484   G4bool res = true;
1485     if(e < maxKinEnergyCSDA) { maxKinEnergyCS << 1485   if ( aTable ) {
1486   } else { PrintWarning("SetMaxKinEnergy", e) << 1486     const G4String name = GetPhysicsTableFileName(part,directory,tname,ascii);
1487 }                                             << 1487     if( !aTable->StorePhysicsTable(name,ascii)) res = false;
1488                                               << 1488   }
1489 //....oooOO0OOooo........oooOO0OOooo........o << 1489   return res;
1490                                               << 1490 }
1491 void G4VEnergyLossProcess::PrintWarning(const << 1491 
1492 {                                             << 1492 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
1493   G4String ss = "G4VEnergyLossProcess::" + ti << 1493 
1494   G4ExceptionDescription ed;                  << 1494 G4bool G4VEnergyLossProcess::RetrievePhysicsTable(
1495   ed << "Parameter is out of range: " << val  << 1495        const G4ParticleDefinition* part, const G4String& directory,
1496      << " it will have no effect!\n" << "  Pr << 1496        G4bool ascii)
1497      << GetProcessName() << "  nbins= " << nB << 1497 {
1498      << " Emin(keV)= " << minKinEnergy/keV    << 1498   G4bool res = true;
1499      << " Emax(GeV)= " << maxKinEnergy/GeV;   << 1499   const G4String particleName = part->GetParticleName();
1500   G4Exception(ss, "em0044", JustWarning, ed); << 1500 
1501 }                                             << 1501   if(1 < verboseLevel) {
1502                                               << 1502     G4cout << "G4VEnergyLossProcess::RetrievePhysicsTable() for "
1503 //....oooOO0OOooo........oooOO0OOooo........o << 1503            << particleName << " and process " << GetProcessName()
1504                                               << 1504            << "; tables_are_built= " << tablesAreBuilt
1505 void G4VEnergyLossProcess::ProcessDescription << 1505            << G4endl;
1506 {                                             << 1506   }
1507   if(nullptr != particle) { StreamInfo(out, * << 1507   if(particle == part) {
1508 }                                             << 1508 
1509                                               << 1509     //    G4bool yes = true;
1510 //....oooOO0OOooo........oooOO0OOooo........o << 1510     if ( !baseParticle ) {
                                                   >> 1511 
                                                   >> 1512       G4bool fpi = true;
                                                   >> 1513       if(!RetrieveTable(part,theDEDXTable,ascii,directory,"DEDX",fpi)) 
                                                   >> 1514   {fpi = false;}
                                                   >> 1515 
                                                   >> 1516       if(!RetrieveTable(part,theIonisationTable,ascii,directory,"Ionisation",false)) 
                                                   >> 1517   {fpi = false;}
                                                   >> 1518 
                                                   >> 1519       if(!RetrieveTable(part,theRangeTableForLoss,ascii,directory,"Range",fpi)) 
                                                   >> 1520         {res = false;}
                                                   >> 1521 
                                                   >> 1522       if(!RetrieveTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr",false)) 
                                                   >> 1523   {res = false;}
                                                   >> 1524 
                                                   >> 1525       if(!RetrieveTable(part,theCSDARangeTable,ascii,directory,"CSDARange",false)) 
                                                   >> 1526   {res = false;}
                                                   >> 1527 
                                                   >> 1528       if(!RetrieveTable(part,theInverseRangeTable,ascii,directory,"InverseRange",fpi)) 
                                                   >> 1529         {res = false;}
                                                   >> 1530 
                                                   >> 1531       if(!RetrieveTable(part,theLambdaTable,ascii,directory,"Lambda",true)) 
                                                   >> 1532         {res = false;}
                                                   >> 1533 
                                                   >> 1534       G4bool yes = false;
                                                   >> 1535       if(nSCoffRegions > 0) {yes = true;}
                                                   >> 1536 
                                                   >> 1537       if(!RetrieveTable(part,theDEDXSubTable,ascii,directory,"SubDEDX",yes)) 
                                                   >> 1538         {res = false;}
                                                   >> 1539 
                                                   >> 1540       if(!RetrieveTable(part,theSubLambdaTable,ascii,directory,"SubLambda",yes)) 
                                                   >> 1541         {res = false;}
                                                   >> 1542 
                                                   >> 1543       if(!fpi) yes = false;
                                                   >> 1544       if(!RetrieveTable(part,theIonisationSubTable,ascii,directory,"SubIonisation",yes)) 
                                                   >> 1545         {res = false;}
                                                   >> 1546     }
                                                   >> 1547   }
                                                   >> 1548 
                                                   >> 1549   return res;
                                                   >> 1550 }
                                                   >> 1551 
                                                   >> 1552 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
                                                   >> 1553 
                                                   >> 1554 G4bool G4VEnergyLossProcess::RetrieveTable(const G4ParticleDefinition* part, 
                                                   >> 1555              G4PhysicsTable* aTable, G4bool ascii,
                                                   >> 1556              const G4String& directory,
                                                   >> 1557              const G4String& tname,
                                                   >> 1558              G4bool mandatory)
                                                   >> 1559 {
                                                   >> 1560   G4bool res = true;
                                                   >> 1561   if(!aTable) { 
                                                   >> 1562     if(mandatory) {res = false;} 
                                                   >> 1563     return res;
                                                   >> 1564   }
                                                   >> 1565   G4String filename = GetPhysicsTableFileName(part,directory,tname,ascii);
                                                   >> 1566   G4bool yes = aTable->ExistPhysicsTable(filename);
                                                   >> 1567   if(yes) {
                                                   >> 1568     yes = G4PhysicsTableHelper::RetrievePhysicsTable(aTable,filename,ascii);
                                                   >> 1569     if((G4LossTableManager::Instance())->SplineFlag()) {
                                                   >> 1570       size_t n = aTable->length();
                                                   >> 1571       for(size_t i=0; i<n; i++) {(*aTable)[i]->SetSpline(true);}
                                                   >> 1572     }
                                                   >> 1573   }
                                                   >> 1574   if(yes) {
                                                   >> 1575     if (0 < verboseLevel) {
                                                   >> 1576       G4cout << tname << " table for " << part->GetParticleName() 
                                                   >> 1577        << " is Retrieved from <" << filename << ">"
                                                   >> 1578        << G4endl;
                                                   >> 1579     }
                                                   >> 1580   } else {
                                                   >> 1581     if(mandatory) res = false;
                                                   >> 1582     if(mandatory || 1 < verboseLevel) {
                                                   >> 1583       G4cout << tname << " table for " << part->GetParticleName() 
                                                   >> 1584        << " from file <"
                                                   >> 1585        << filename << "> is not Retrieved"
                                                   >> 1586        << G4endl;
                                                   >> 1587     }
                                                   >> 1588   }
                                                   >> 1589   return res;
                                                   >> 1590 }
                                                   >> 1591 
                                                   >> 1592 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >> 1593   
                                                   >> 1594 void G4VEnergyLossProcess::AddCollaborativeProcess(
                                                   >> 1595             G4VEnergyLossProcess* p)
                                                   >> 1596 {
                                                   >> 1597   G4bool add = true;
                                                   >> 1598   if(p->GetProcessName() != "eBrem") add = false;
                                                   >> 1599   if(add && nProcesses > 0) {
                                                   >> 1600     for(G4int i=0; i<nProcesses; i++) {
                                                   >> 1601       if(p == scProcesses[i]) {
                                                   >> 1602         add = false;
                                                   >> 1603         break;
                                                   >> 1604       }
                                                   >> 1605     }
                                                   >> 1606   }
                                                   >> 1607   if(add) {
                                                   >> 1608     scProcesses.push_back(p);
                                                   >> 1609     nProcesses++;
                                                   >> 1610     if (1 < verboseLevel) { 
                                                   >> 1611       G4cout << "### The process " << p->GetProcessName() 
                                                   >> 1612        << " is added to the list of collaborative processes of "
                                                   >> 1613        << GetProcessName() << G4endl; 
                                                   >> 1614     }
                                                   >> 1615   }
                                                   >> 1616 }
                                                   >> 1617 
                                                   >> 1618 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >> 1619 
                                                   >> 1620 G4double G4VEnergyLossProcess::GetDEDXDispersion(
                                                   >> 1621                                   const G4MaterialCutsCouple *couple,
                                                   >> 1622                                   const G4DynamicParticle* dp,
                                                   >> 1623                                         G4double length)
                                                   >> 1624 {
                                                   >> 1625   DefineMaterial(couple);
                                                   >> 1626   G4double ekin = dp->GetKineticEnergy();
                                                   >> 1627   SelectModel(ekin*massRatio);
                                                   >> 1628   G4double tmax = currentModel->MaxSecondaryKinEnergy(dp);
                                                   >> 1629   tmax = std::min(tmax,(*theCuts)[currentMaterialIndex]);
                                                   >> 1630   G4double d = 0.0;
                                                   >> 1631   G4VEmFluctuationModel* fm = currentModel->GetModelOfFluctuations();
                                                   >> 1632   if(fm) d = fm->Dispersion(currentMaterial,dp,tmax,length);
                                                   >> 1633   return d;
                                                   >> 1634 }
                                                   >> 1635 
                                                   >> 1636 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >> 1637 
                                                   >> 1638 void G4VEnergyLossProcess::ActivateDeexcitation(G4bool, const G4Region*)
                                                   >> 1639 {}
                                                   >> 1640 
                                                   >> 1641 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >> 1642 
1511                                                  1643