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.3.p3)


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