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 10.4)


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