Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/utils/src/G4VEmProcess.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/G4VEmProcess.cc (Version 11.3.0) and /processes/electromagnetic/utils/src/G4VEmProcess.cc (Version 10.0.p1)


  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: G4VEmProcess.cc 79268 2014-02-20 16:46:31Z gcosmo $
                                                   >>  27 //
 26 // -------------------------------------------     28 // -------------------------------------------------------------------
 27 //                                                 29 //
 28 // GEANT4 Class file                               30 // GEANT4 Class file
 29 //                                                 31 //
 30 //                                                 32 //
 31 // File name:     G4VEmProcess                     33 // File name:     G4VEmProcess
 32 //                                                 34 //
 33 // Author:        Vladimir Ivanchenko on base      35 // Author:        Vladimir Ivanchenko on base of Laszlo Urban code
 34 //                                                 36 //
 35 // Creation date: 01.10.2003                       37 // Creation date: 01.10.2003
 36 //                                                 38 //
 37 // Modifications: by V.Ivanchenko              <<  39 // Modifications:
                                                   >>  40 // 30-06-04 make it to be pure discrete process (V.Ivanchenko)
                                                   >>  41 // 30-09-08 optimise integral option (V.Ivanchenko)
                                                   >>  42 // 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivanchenko)
                                                   >>  43 // 11-03-05 Shift verbose level by 1, add applyCuts and killPrimary flags (VI)
                                                   >>  44 // 14-03-05 Update logic PostStepDoIt (V.Ivanchenko)
                                                   >>  45 // 08-04-05 Major optimisation of internal interfaces (V.Ivanchenko)
                                                   >>  46 // 18-04-05 Use G4ParticleChangeForGamma (V.Ivanchenko)
                                                   >>  47 // 25-07-05 Add protection: integral mode only for charged particles (VI)
                                                   >>  48 // 04-09-05 default lambdaFactor 0.8 (V.Ivanchenko)
                                                   >>  49 // 11-01-06 add A to parameters of ComputeCrossSectionPerAtom (VI)
                                                   >>  50 // 12-09-06 add SetModel() (mma)
                                                   >>  51 // 12-04-07 remove double call to Clear model manager (V.Ivanchenko)
                                                   >>  52 // 27-10-07 Virtual functions moved to source (V.Ivanchenko)
                                                   >>  53 // 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko)
                                                   >>  54 // 17-02-10 Added pointer currentParticle (VI)
                                                   >>  55 // 30-05-12 allow Russian roulette, brem splitting (D. Sawkey)
 38 //                                                 56 //
 39 // Class Description: based class for discrete <<  57 // Class Description:
 40 //                                                 58 //
 41                                                    59 
 42 // -------------------------------------------     60 // -------------------------------------------------------------------
 43 //                                                 61 //
 44 //....oooOO0OOooo........oooOO0OOooo........oo     62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 45 //....oooOO0OOooo........oooOO0OOooo........oo     63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 46                                                    64 
 47 #include "G4VEmProcess.hh"                         65 #include "G4VEmProcess.hh"
 48 #include "G4PhysicalConstants.hh"                  66 #include "G4PhysicalConstants.hh"
 49 #include "G4SystemOfUnits.hh"                      67 #include "G4SystemOfUnits.hh"
 50 #include "G4ProcessManager.hh"                     68 #include "G4ProcessManager.hh"
 51 #include "G4LossTableManager.hh"                   69 #include "G4LossTableManager.hh"
 52 #include "G4LossTableBuilder.hh"                   70 #include "G4LossTableBuilder.hh"
 53 #include "G4Step.hh"                               71 #include "G4Step.hh"
 54 #include "G4ParticleDefinition.hh"                 72 #include "G4ParticleDefinition.hh"
 55 #include "G4VEmModel.hh"                           73 #include "G4VEmModel.hh"
 56 #include "G4DataVector.hh"                         74 #include "G4DataVector.hh"
 57 #include "G4PhysicsTable.hh"                       75 #include "G4PhysicsTable.hh"
 58 #include "G4EmDataHandler.hh"                  << 
 59 #include "G4PhysicsLogVector.hh"                   76 #include "G4PhysicsLogVector.hh"
 60 #include "G4VParticleChange.hh"                    77 #include "G4VParticleChange.hh"
 61 #include "G4ProductionCutsTable.hh"                78 #include "G4ProductionCutsTable.hh"
 62 #include "G4Region.hh"                             79 #include "G4Region.hh"
 63 #include "G4Gamma.hh"                              80 #include "G4Gamma.hh"
 64 #include "G4Electron.hh"                           81 #include "G4Electron.hh"
 65 #include "G4Positron.hh"                           82 #include "G4Positron.hh"
 66 #include "G4PhysicsTableHelper.hh"                 83 #include "G4PhysicsTableHelper.hh"
 67 #include "G4EmBiasingManager.hh"                   84 #include "G4EmBiasingManager.hh"
 68 #include "G4EmParameters.hh"                   << 
 69 #include "G4EmProcessSubType.hh"               << 
 70 #include "G4EmTableUtil.hh"                    << 
 71 #include "G4EmUtility.hh"                      << 
 72 #include "G4DNAModelSubType.hh"                << 
 73 #include "G4GenericIon.hh"                         85 #include "G4GenericIon.hh"
 74 #include "G4Log.hh"                                86 #include "G4Log.hh"
 75 #include <iostream>                            << 
 76                                                    87 
 77 //....oooOO0OOooo........oooOO0OOooo........oo     88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 78                                                    89 
 79 G4VEmProcess::G4VEmProcess(const G4String& nam     90 G4VEmProcess::G4VEmProcess(const G4String& name, G4ProcessType type):
 80   G4VDiscreteProcess(name, type)               <<  91   G4VDiscreteProcess(name, type),
                                                   >>  92   secondaryParticle(0),
                                                   >>  93   buildLambdaTable(true),
                                                   >>  94   numberOfModels(0),
                                                   >>  95   theLambdaTable(0),
                                                   >>  96   theLambdaTablePrim(0),
                                                   >>  97   theDensityFactor(0),
                                                   >>  98   theDensityIdx(0),
                                                   >>  99   integral(false),
                                                   >> 100   applyCuts(false),
                                                   >> 101   startFromNull(false),
                                                   >> 102   splineFlag(true),
                                                   >> 103   currentModel(0),
                                                   >> 104   particle(0),
                                                   >> 105   currentParticle(0),
                                                   >> 106   currentCouple(0)
 81 {                                                 107 {
 82   theParameters = G4EmParameters::Instance();  << 
 83   SetVerboseLevel(1);                             108   SetVerboseLevel(1);
 84                                                   109 
 85   // Size of tables                            << 110   // Size of tables assuming spline
 86   minKinEnergy = 0.1*CLHEP::keV;               << 111   minKinEnergy = 0.1*keV;
 87   maxKinEnergy = 100.0*CLHEP::TeV;             << 112   maxKinEnergy = 10.0*TeV;
                                                   >> 113   nLambdaBins  = 77;
                                                   >> 114   minKinEnergyPrim = DBL_MAX;
 88                                                   115 
 89   // default lambda factor                        116   // default lambda factor
 90   invLambdaFactor = 1.0/lambdaFactor;          << 117   lambdaFactor  = 0.8;
                                                   >> 118 
                                                   >> 119   // default limit on polar angle
                                                   >> 120   polarAngleLimit = 0.0;
                                                   >> 121   biasFactor = 1.0;
 91                                                   122 
 92   // particle types                               123   // particle types
 93   theGamma = G4Gamma::Gamma();                 << 124   theGamma     = G4Gamma::Gamma();
 94   theElectron = G4Electron::Electron();        << 125   theElectron  = G4Electron::Electron();
 95   thePositron = G4Positron::Positron();        << 126   thePositron  = G4Positron::Positron();
 96                                                   127 
 97   pParticleChange = &fParticleChange;             128   pParticleChange = &fParticleChange;
 98   fParticleChange.SetSecondaryWeightByProcess(    129   fParticleChange.SetSecondaryWeightByProcess(true);
 99   secParticles.reserve(5);                        130   secParticles.reserve(5);
100                                                   131 
                                                   >> 132   preStepLambda = 0.0;
                                                   >> 133   mfpKinEnergy  = DBL_MAX;
                                                   >> 134 
                                                   >> 135   idxLambda = idxLambdaPrim = 0;
                                                   >> 136 
101   modelManager = new G4EmModelManager();          137   modelManager = new G4EmModelManager();
                                                   >> 138   biasManager  = 0;
                                                   >> 139   biasFlag     = false; 
                                                   >> 140   weightFlag   = false;
102   lManager = G4LossTableManager::Instance();      141   lManager = G4LossTableManager::Instance();
103   lManager->Register(this);                       142   lManager->Register(this);
104   isTheMaster = lManager->IsMaster();          << 143   secID = fluoID = augerID = biasID = -1;
105   G4LossTableBuilder* bld = lManager->GetTable << 144   mainSecondaries = 100;
106   theDensityFactor = bld->GetDensityFactors(); << 145   if("phot" == GetProcessName() || "compt" == GetProcessName()) { 
107   theDensityIdx = bld->GetCoupleIndexes();     << 146     mainSecondaries = 1; 
                                                   >> 147   }
108 }                                                 148 }
109                                                   149 
110 //....oooOO0OOooo........oooOO0OOooo........oo    150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111                                                   151 
112 G4VEmProcess::~G4VEmProcess()                     152 G4VEmProcess::~G4VEmProcess()
113 {                                                 153 {
114   if(isTheMaster) {                            << 154   if(1 < verboseLevel) {
115     delete theData;                            << 155     G4cout << "G4VEmProcess destruct " << GetProcessName() 
116     delete theEnergyOfCrossSectionMax;         << 156      << "  " << this << "  " <<  theLambdaTable <<G4endl;
                                                   >> 157   }
                                                   >> 158   if(lManager->IsMaster()) {
                                                   >> 159     delete theLambdaTable;
                                                   >> 160     delete theLambdaTablePrim;
117   }                                               161   }
118   delete modelManager;                            162   delete modelManager;
119   delete biasManager;                             163   delete biasManager;
120   lManager->DeRegister(this);                     164   lManager->DeRegister(this);
                                                   >> 165   //G4cout << "G4VEmProcess removed " << G4endl; 
121 }                                                 166 }
122                                                   167 
123 //....oooOO0OOooo........oooOO0OOooo........oo    168 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
124                                                   169 
125 void G4VEmProcess::AddEmModel(G4int order, G4V << 170 void G4VEmProcess::Clear()
126                               const G4Region*  << 
127 {                                                 171 {
128   if(nullptr == ptr) { return; }               << 172   currentCouple = 0;
129   G4VEmFluctuationModel* fm = nullptr;         << 173   preStepLambda = 0.0;
130   modelManager->AddEmModel(order, ptr, fm, reg << 174   mfpKinEnergy  = DBL_MAX;
131   ptr->SetParticleChange(pParticleChange);     << 175   idxLambda = idxLambdaPrim = 0;
132 }                                                 176 }
133                                                   177 
134 //....oooOO0OOooo........oooOO0OOooo........oo    178 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
135                                                   179 
136 void G4VEmProcess::SetEmModel(G4VEmModel* ptr, << 180 G4double G4VEmProcess::MinPrimaryEnergy(const G4ParticleDefinition*,
                                                   >> 181           const G4Material*)
137 {                                                 182 {
138   if(nullptr == ptr) { return; }               << 183   return 0.0;
139   if(!emModels.empty()) {                      << 184 }
140     for(auto & em : emModels) { if(em == ptr)  << 185 
141   }                                            << 186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
142   emModels.push_back(ptr);                     << 187 
                                                   >> 188 void G4VEmProcess::AddEmModel(G4int order, G4VEmModel* p, 
                                                   >> 189             const G4Region* region)
                                                   >> 190 {
                                                   >> 191   G4VEmFluctuationModel* fm = 0;
                                                   >> 192   modelManager->AddEmModel(order, p, fm, region);
                                                   >> 193   if(p) { p->SetParticleChange(pParticleChange); }
                                                   >> 194 }
                                                   >> 195 
                                                   >> 196 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >> 197 
                                                   >> 198 void G4VEmProcess::SetEmModel(G4VEmModel* p, G4int index) 
                                                   >> 199 {
                                                   >> 200   G4int n = emModels.size();
                                                   >> 201   if(index >= n) { for(G4int i=n; i<=index; ++i) {emModels.push_back(0);} }
                                                   >> 202   emModels[index] = p;
                                                   >> 203 }
                                                   >> 204 
                                                   >> 205 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >> 206 
                                                   >> 207 G4VEmModel* G4VEmProcess::EmModel(G4int index) const
                                                   >> 208 {
                                                   >> 209   G4VEmModel* p = 0;
                                                   >> 210   if(index >= 0 && index <  G4int(emModels.size())) { p = emModels[index]; }
                                                   >> 211   return p;
                                                   >> 212 }
                                                   >> 213 
                                                   >> 214 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >> 215 
                                                   >> 216 void G4VEmProcess::UpdateEmModel(const G4String& nam, 
                                                   >> 217          G4double emin, G4double emax)
                                                   >> 218 {
                                                   >> 219   modelManager->UpdateEmModel(nam, emin, emax);
                                                   >> 220 }
                                                   >> 221 
                                                   >> 222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >> 223 
                                                   >> 224 G4VEmModel* G4VEmProcess::GetModelByIndex(G4int idx, G4bool ver) const
                                                   >> 225 {
                                                   >> 226   return modelManager->GetModel(idx, ver);
143 }                                                 227 }
144                                                   228 
145 //....oooOO0OOooo........oooOO0OOooo........oo    229 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
146                                                   230 
147 void G4VEmProcess::PreparePhysicsTable(const G    231 void G4VEmProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
148 {                                                 232 {
149   if(nullptr == particle) { SetParticle(&part) << 233   G4bool isMaster = false;
                                                   >> 234   if(GetMasterProcess() == this) { isMaster = true; }
                                                   >> 235   if(!particle) { SetParticle(&part); }
150                                                   236 
151   if(part.GetParticleType() == "nucleus" &&       237   if(part.GetParticleType() == "nucleus" && 
152      part.GetParticleSubType() == "generic") {    238      part.GetParticleSubType() == "generic") {
153                                                   239 
154     G4String pname = part.GetParticleName();      240     G4String pname = part.GetParticleName();
155     if(pname != "deuteron" && pname != "triton    241     if(pname != "deuteron" && pname != "triton" &&
156        pname != "He3" && pname != "alpha" && p << 242        pname != "alpha" && pname != "He3" &&
157        pname != "helium" && pname != "hydrogen << 243        pname != "alpha+"   && pname != "helium" &&
                                                   >> 244        pname != "hydrogen") {
158                                                   245 
159       particle = G4GenericIon::GenericIon();      246       particle = G4GenericIon::GenericIon();
160       isIon = true;                            << 
161     }                                             247     }
162   }                                               248   }
                                                   >> 249 
                                                   >> 250   if(1 < verboseLevel) {
                                                   >> 251     G4cout << "G4VEmProcess::PreparePhysicsTable() for "
                                                   >> 252            << GetProcessName()
                                                   >> 253            << " and particle " << part.GetParticleName()
                                                   >> 254      << " local particle " << particle->GetParticleName() 
                                                   >> 255            << G4endl;
                                                   >> 256   }
                                                   >> 257 
163   if(particle != &part) { return; }               258   if(particle != &part) { return; }
164                                                   259 
165   lManager->PreparePhysicsTable(&part, this);  << 260   G4LossTableBuilder* bld = lManager->GetTableBuilder();
166                                                   261 
167   // for new run                               << 262   lManager->PreparePhysicsTable(&part, this, isMaster);
168   currentCouple = nullptr;                     << 
169   preStepLambda = 0.0;                         << 
170   fLambdaEnergy = 0.0;                         << 
171                                                   263 
                                                   >> 264   Clear();
172   InitialiseProcess(particle);                    265   InitialiseProcess(particle);
173                                                   266 
174   G4LossTableBuilder* bld = lManager->GetTable << 
175   const G4ProductionCutsTable* theCoupleTable=    267   const G4ProductionCutsTable* theCoupleTable=
176     G4ProductionCutsTable::GetProductionCutsTa    268     G4ProductionCutsTable::GetProductionCutsTable();
                                                   >> 269   size_t n = theCoupleTable->GetTableSize();
                                                   >> 270 
                                                   >> 271   theEnergyOfCrossSectionMax.resize(n, 0.0);
                                                   >> 272   theCrossSectionMax.resize(n, DBL_MAX);
                                                   >> 273 
                                                   >> 274   // initialisation of models
                                                   >> 275   numberOfModels = modelManager->NumberOfModels();
                                                   >> 276   for(G4int i=0; i<numberOfModels; ++i) {
                                                   >> 277     G4VEmModel* mod = modelManager->GetModel(i);
                                                   >> 278     if(0 == i) { currentModel = mod; }
                                                   >> 279     mod->SetPolarAngleLimit(polarAngleLimit);
                                                   >> 280     mod->SetMasterThread(isMaster);
                                                   >> 281     if(mod->HighEnergyLimit() > maxKinEnergy) {
                                                   >> 282       mod->SetHighEnergyLimit(maxKinEnergy);
                                                   >> 283     }
                                                   >> 284   }
                                                   >> 285 
                                                   >> 286   if(lManager->AtomDeexcitation()) { modelManager->SetFluoFlag(true); }
                                                   >> 287   theCuts = modelManager->Initialise(particle,secondaryParticle,
                                                   >> 288              2.,verboseLevel);
177   theCutsGamma    = theCoupleTable->GetEnergyC    289   theCutsGamma    = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut);
178   theCutsElectron = theCoupleTable->GetEnergyC    290   theCutsElectron = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut);
179   theCutsPositron = theCoupleTable->GetEnergyC    291   theCutsPositron = theCoupleTable->GetEnergyCutsVector(idxG4PositronCut);
180                                                   292 
181   // initialisation of the process             << 
182   if(!actMinKinEnergy) { minKinEnergy = thePar << 
183   if(!actMaxKinEnergy) { maxKinEnergy = thePar << 
184                                                << 
185   applyCuts       = theParameters->ApplyCuts() << 
186   lambdaFactor    = theParameters->LambdaFacto << 
187   invLambdaFactor = 1.0/lambdaFactor;          << 
188   theParameters->DefineRegParamForEM(this);    << 
189                                                << 
190   // integral option may be disabled           << 
191   if(!theParameters->Integral()) { fXSType = f << 
192                                                << 
193   // prepare tables                               293   // prepare tables
194   if(isTheMaster) {                            << 294   if(buildLambdaTable && isMaster){
195     if(nullptr == theData) { theData = new G4E << 295     theLambdaTable = 
196                                                << 296       G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable);
197     if(buildLambdaTable) {                     << 297     bld->InitialiseBaseMaterials(theLambdaTable);
198       theLambdaTable = theData->MakeTable(0);  << 298   }
199       bld->InitialiseBaseMaterials(theLambdaTa << 299   // high energy table
200     }                                          << 300   if(isMaster && minKinEnergyPrim < maxKinEnergy){
201     // high energy table                       << 301     theLambdaTablePrim = 
202     if(minKinEnergyPrim < maxKinEnergy) {      << 302       G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTablePrim);
203       theLambdaTablePrim = theData->MakeTable( << 303     bld->InitialiseBaseMaterials(theLambdaTablePrim);
204       bld->InitialiseBaseMaterials(theLambdaTa << 
205     }                                          << 
206   }                                            << 
207   // models                                    << 
208   baseMat = bld->GetBaseMaterialFlag();        << 
209   numberOfModels = modelManager->NumberOfModel << 
210   currentModel = modelManager->GetModel(0);    << 
211   if(nullptr != lManager->AtomDeexcitation())  << 
212     modelManager->SetFluoFlag(true);           << 
213   }                                               304   }
214   // forced biasing                               305   // forced biasing
215   if(nullptr != biasManager) {                 << 306   if(biasManager) { 
216     biasManager->Initialise(part, GetProcessNa << 307     biasManager->Initialise(part,GetProcessName(),verboseLevel); 
217     biasFlag = false;                          << 308     biasFlag = false; 
                                                   >> 309   }
                                                   >> 310   // defined ID of secondary particles
                                                   >> 311   if(isMaster) {
                                                   >> 312     G4String nam1 = GetProcessName();
                                                   >> 313     G4String nam2 = nam1 + "_fluo" ;
                                                   >> 314     G4String nam3 = nam1 + "_auger";
                                                   >> 315     G4String nam4 = nam1 + "_split";
                                                   >> 316     secID   = G4PhysicsModelCatalog::Register(nam1); 
                                                   >> 317     fluoID  = G4PhysicsModelCatalog::Register(nam2); 
                                                   >> 318     augerID = G4PhysicsModelCatalog::Register(nam3); 
                                                   >> 319     biasID  = G4PhysicsModelCatalog::Register(nam4); 
218   }                                               320   }
219                                                << 
220   theCuts =                                    << 
221     G4EmTableUtil::PrepareEmProcess(this, part << 
222                                     modelManag << 
223                                     secID, tri << 
224                                     verboseLev << 
225 }                                                 321 }
226                                                   322 
227 //....oooOO0OOooo........oooOO0OOooo........oo    323 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
228                                                   324 
229 void G4VEmProcess::BuildPhysicsTable(const G4P    325 void G4VEmProcess::BuildPhysicsTable(const G4ParticleDefinition& part)
230 {                                                 326 {
231   if(nullptr == masterProc) {                  << 327   const G4VEmProcess* masterProc = 0;
232     if(isTheMaster) { masterProc = this; }     << 328   if(GetMasterProcess() != this) {
233     else { masterProc = static_cast<const G4VE << 329     masterProc = static_cast<const G4VEmProcess*>(GetMasterProcess()); 
234   }                                            << 330   }
235   G4int nModels = modelManager->NumberOfModels << 331 
236   G4bool isLocked = theParameters->IsPrintLock << 332   G4String num = part.GetParticleName();
237   G4bool toBuild = (buildLambdaTable || minKin << 333   if(1 < verboseLevel) {
238                                                << 334     G4cout << "G4VEmProcess::BuildPhysicsTable() for "
239   G4EmTableUtil::BuildEmProcess(this, masterPr << 335            << GetProcessName()
240                                 nModels, verbo << 336            << " and particle " << num
241                                 isLocked, toBu << 337      << " buildLambdaTable= " << buildLambdaTable
                                                   >> 338            << G4endl;
                                                   >> 339   }
                                                   >> 340 
                                                   >> 341   if(particle == &part) { 
                                                   >> 342 
                                                   >> 343     G4LossTableBuilder* bld = lManager->GetTableBuilder();
                                                   >> 344 
                                                   >> 345     // worker initialisation
                                                   >> 346     if(masterProc) {
                                                   >> 347       theLambdaTable = masterProc->LambdaTable();
                                                   >> 348       theLambdaTablePrim = masterProc->LambdaTablePrim();
                                                   >> 349 
                                                   >> 350       if(theLambdaTable) {
                                                   >> 351   bld->InitialiseBaseMaterials(theLambdaTable);
                                                   >> 352       } else if(theLambdaTablePrim) {
                                                   >> 353   bld->InitialiseBaseMaterials(theLambdaTablePrim);
                                                   >> 354       }
                                                   >> 355       theDensityFactor = bld->GetDensityFactors();
                                                   >> 356       theDensityIdx = bld->GetCoupleIndexes();
                                                   >> 357       if(theLambdaTable) { FindLambdaMax(); }
                                                   >> 358 
                                                   >> 359       // local initialisation of models
                                                   >> 360       G4bool printing = true;
                                                   >> 361       numberOfModels = modelManager->NumberOfModels();
                                                   >> 362       for(G4int i=0; i<numberOfModels; ++i) {
                                                   >> 363   G4VEmModel* mod = GetModelByIndex(i, printing);
                                                   >> 364   G4VEmModel* mod0= masterProc->GetModelByIndex(i, printing);
                                                   >> 365   mod->InitialiseLocal(particle, mod0);
                                                   >> 366       }
                                                   >> 367     // master thread
                                                   >> 368     } else {
                                                   >> 369       theDensityFactor = bld->GetDensityFactors();
                                                   >> 370       theDensityIdx = bld->GetCoupleIndexes();
                                                   >> 371       if(buildLambdaTable || minKinEnergyPrim < maxKinEnergy) {
                                                   >> 372   BuildLambdaTable();
                                                   >> 373       }
                                                   >> 374     }
                                                   >> 375   }
                                                   >> 376 
                                                   >> 377   // explicitly defined printout by particle name
                                                   >> 378   if(1 < verboseLevel || 
                                                   >> 379      (0 < verboseLevel && (num == "gamma" || num == "e-" || 
                                                   >> 380          num == "e+"    || num == "mu+" || 
                                                   >> 381          num == "mu-"   || num == "proton"|| 
                                                   >> 382          num == "pi+"   || num == "pi-" || 
                                                   >> 383          num == "kaon+" || num == "kaon-" || 
                                                   >> 384          num == "alpha" || num == "anti_proton" || 
                                                   >> 385          num == "GenericIon")))
                                                   >> 386     { 
                                                   >> 387       PrintInfoProcess(part); 
                                                   >> 388     }
                                                   >> 389 
                                                   >> 390   if(1 < verboseLevel) {
                                                   >> 391     G4cout << "G4VEmProcess::BuildPhysicsTable() done for "
                                                   >> 392            << GetProcessName()
                                                   >> 393            << " and particle " << num
                                                   >> 394            << G4endl;
                                                   >> 395   }
242 }                                                 396 }
243                                                   397 
244 //....oooOO0OOooo........oooOO0OOooo........oo    398 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
245                                                   399 
246 void G4VEmProcess::BuildLambdaTable()             400 void G4VEmProcess::BuildLambdaTable()
247 {                                                 401 {
248   G4double scale = theParameters->MaxKinEnergy << 402   if(1 < verboseLevel) {
249   G4int nbin =                                 << 403     G4cout << "G4EmProcess::BuildLambdaTable() for process "
250     theParameters->NumberOfBinsPerDecade()*G4l << 404            << GetProcessName() << " and particle "
251   if(actBinning) { nbin = std::max(nbin, nLamb << 405            << particle->GetParticleName() << "  " << this
252   scale = nbin/G4Log(scale);                   << 406            << G4endl;
253                                                << 407   }
                                                   >> 408 
                                                   >> 409   // Access to materials
                                                   >> 410   const G4ProductionCutsTable* theCoupleTable=
                                                   >> 411         G4ProductionCutsTable::GetProductionCutsTable();
                                                   >> 412   size_t numOfCouples = theCoupleTable->GetTableSize();
                                                   >> 413 
254   G4LossTableBuilder* bld = lManager->GetTable    414   G4LossTableBuilder* bld = lManager->GetTableBuilder();
255   G4EmTableUtil::BuildLambdaTable(this, partic << 415 
256                                   bld, theLamb << 416   G4PhysicsLogVector* aVector = 0;
257                                   minKinEnergy << 417   G4PhysicsLogVector* aVectorPrim = 0;
258                                   maxKinEnergy << 418   G4PhysicsLogVector* bVectorPrim = 0;
259                                   startFromNul << 419 
                                                   >> 420   G4double scale = G4Log(maxKinEnergy/minKinEnergy); 
                                                   >> 421   G4double emax1 = maxKinEnergy;
                                                   >> 422   if(minKinEnergyPrim < maxKinEnergy) { emax1 = minKinEnergyPrim; }
                                                   >> 423     
                                                   >> 424   for(size_t i=0; i<numOfCouples; ++i) {
                                                   >> 425 
                                                   >> 426     if (bld->GetFlag(i)) {
                                                   >> 427 
                                                   >> 428       // create physics vector and fill it
                                                   >> 429       const G4MaterialCutsCouple* couple = 
                                                   >> 430   theCoupleTable->GetMaterialCutsCouple(i);
                                                   >> 431 
                                                   >> 432       // build main table
                                                   >> 433       if(buildLambdaTable) {
                                                   >> 434   delete (*theLambdaTable)[i];
                                                   >> 435 
                                                   >> 436   // if start from zero then change the scale
                                                   >> 437   G4double emin = minKinEnergy;
                                                   >> 438   G4bool startNull = false;
                                                   >> 439   if(startFromNull) {
                                                   >> 440     G4double e = MinPrimaryEnergy(particle,couple->GetMaterial());
                                                   >> 441           if(e >= emin) {
                                                   >> 442       emin = e;
                                                   >> 443       startNull = true;
                                                   >> 444     }
                                                   >> 445   }
                                                   >> 446   G4double emax = emax1;
                                                   >> 447   if(emax <= emin) { emax = 2*emin; }
                                                   >> 448   G4int bin = G4lrint(nLambdaBins*G4Log(emax/emin)/scale);
                                                   >> 449   if(bin < 3) { bin = 3; }
                                                   >> 450   aVector = new G4PhysicsLogVector(emin, emax, bin);
                                                   >> 451   aVector->SetSpline(splineFlag);
                                                   >> 452   modelManager->FillLambdaVector(aVector, couple, startNull);
                                                   >> 453   if(splineFlag) { aVector->FillSecondDerivatives(); }
                                                   >> 454   G4PhysicsTableHelper::SetPhysicsVector(theLambdaTable, i, aVector);
                                                   >> 455       }
                                                   >> 456       // build high energy table 
                                                   >> 457       if(minKinEnergyPrim < maxKinEnergy) { 
                                                   >> 458   delete (*theLambdaTablePrim)[i];
                                                   >> 459 
                                                   >> 460   // start not from zero
                                                   >> 461   if(!bVectorPrim) {
                                                   >> 462     G4int bin = 
                                                   >> 463       G4lrint(nLambdaBins*G4Log(maxKinEnergy/minKinEnergyPrim)/scale);
                                                   >> 464     if(bin < 3) { bin = 3; }
                                                   >> 465     aVectorPrim = 
                                                   >> 466       new G4PhysicsLogVector(minKinEnergyPrim, maxKinEnergy, bin);
                                                   >> 467     bVectorPrim = aVectorPrim;
                                                   >> 468   } else {
                                                   >> 469     aVectorPrim = new G4PhysicsLogVector(*bVectorPrim);
                                                   >> 470   }
                                                   >> 471   // always use spline
                                                   >> 472   aVectorPrim->SetSpline(true);
                                                   >> 473   modelManager->FillLambdaVector(aVectorPrim, couple, false, 
                                                   >> 474                fIsCrossSectionPrim);
                                                   >> 475   aVectorPrim->FillSecondDerivatives();
                                                   >> 476   G4PhysicsTableHelper::SetPhysicsVector(theLambdaTablePrim, i, 
                                                   >> 477                  aVectorPrim);
                                                   >> 478       }
                                                   >> 479     }
                                                   >> 480   }
                                                   >> 481 
                                                   >> 482   if(buildLambdaTable) { FindLambdaMax(); }
                                                   >> 483 
                                                   >> 484   if(1 < verboseLevel) {
                                                   >> 485     G4cout << "Lambda table is built for "
                                                   >> 486            << particle->GetParticleName()
                                                   >> 487            << G4endl;
                                                   >> 488   }
260 }                                                 489 }
261                                                   490 
262 //....oooOO0OOooo........oooOO0OOooo........oo    491 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
263                                                   492 
264 void G4VEmProcess::StreamInfo(std::ostream& ou << 493 void G4VEmProcess::PrintInfoProcess(const G4ParticleDefinition& part)
265                   const G4ParticleDefinition&  << 494 {
266 {                                              << 495   if(verboseLevel > 0) {
267   G4String indent = (rst ? "  " : "");         << 496     G4cout << std::setprecision(6);
268   out << std::setprecision(6);                 << 497     G4cout << G4endl << GetProcessName() << ":   for  "
269   out << G4endl << indent << GetProcessName()  << 498            << part.GetParticleName();
270   if (!rst) {                                  << 499     if(integral)  { G4cout << ", integral: 1 "; }
271     out << " for " << part.GetParticleName();  << 500     if(applyCuts) { G4cout << ", applyCuts: 1 "; }
272   }                                            << 501     G4cout << "    SubType= " << GetProcessSubType();;
273   if(fXSType != fEmNoIntegral)  { out << " XSt << 502     if(biasFactor != 1.0) { G4cout << "   BiasingFactor= " << biasFactor; }
274   if(applyCuts) { out << " applyCuts:1 "; }    << 503     G4cout << G4endl;
275   G4int subtype = GetProcessSubType();         << 504     if(buildLambdaTable) {
276   out << " SubType=" << subtype;               << 505       size_t length = theLambdaTable->length();
277   if (subtype == fAnnihilation) {              << 506       for(size_t i=0; i<length; ++i) {
278     G4int mod = theParameters->PositronAtRestM << 507         G4PhysicsVector* v = (*theLambdaTable)[i];
279     const G4String namp[2] = {"Simple", "Allis << 508         if(v) { 
280     out << " AtRestModel:" << namp[mod];       << 509     G4cout << "      Lambda table from "
281   }                                            << 510      << G4BestUnit(minKinEnergy,"Energy") 
282   if(biasFactor != 1.0) { out << "  BiasingFac << 511      << " to "
283   out << " BuildTable=" << buildLambdaTable << << 512      << G4BestUnit(v->GetMaxEnergy(),"Energy")
284   if(buildLambdaTable) {                       << 513      << " in " << v->GetVectorLength()-1
285     if(particle == &part) {                    << 514      << " bins, spline: " 
286       for(auto & v : *theLambdaTable) {        << 515      << splineFlag
287         if(nullptr != v) {                     << 516      << G4endl;
288           out << "      Lambda table from ";   << 517     break;
289           G4double emin = v->Energy(0);        << 518   }
290           G4double emax = v->GetMaxEnergy();   << 
291           G4int nbin = G4int(v->GetVectorLengt << 
292           if(emin > minKinEnergy) { out << "th << 
293           else { out << G4BestUnit(emin,"Energ << 
294           out << " to "                        << 
295               << G4BestUnit(emax,"Energy")     << 
296               << ", " << G4lrint(nbin/std::log << 
297               << " bins/decade, spline: "      << 
298               << splineFlag << G4endl;         << 
299           break;                               << 
300         }                                      << 
301       }                                           519       }
302     } else {                                   << 
303       out << "      Used Lambda table of "     << 
304       << particle->GetParticleName() << G4endl << 
305     }                                             520     }
306   }                                            << 521     if(minKinEnergyPrim < maxKinEnergy) {
307   if(minKinEnergyPrim < maxKinEnergy) {        << 522       size_t length = theLambdaTablePrim->length();
308     if(particle == &part) {                    << 523       for(size_t i=0; i<length; ++i) {
309       for(auto & v : *theLambdaTablePrim) {    << 524         G4PhysicsVector* v = (*theLambdaTablePrim)[i];
310         if(nullptr != v) {                     << 525         if(v) { 
311           out << "      LambdaPrime table from << 526     G4cout << "      LambdaPrime table from "
312               << G4BestUnit(v->Energy(0),"Ener << 527      << G4BestUnit(minKinEnergyPrim,"Energy") 
313               << " to "                        << 528      << " to "
314               << G4BestUnit(v->GetMaxEnergy(), << 529      << G4BestUnit(maxKinEnergy,"Energy")
315               << " in " << v->GetVectorLength( << 530      << " in " << v->GetVectorLength()-1
316               << " bins " << G4endl;           << 531      << " bins " 
317           break;                               << 532      << G4endl;
318         }                                      << 533     break;
                                                   >> 534   }
319       }                                           535       }
320     } else {                                   << 
321       out << "      Used LambdaPrime table of  << 
322                << particle->GetParticleName()  << 
323     }                                             536     }
                                                   >> 537     PrintInfo();
                                                   >> 538     modelManager->DumpModelList(verboseLevel);
324   }                                               539   }
325   StreamProcessInfo(out);                      << 
326   modelManager->DumpModelList(out, verboseLeve << 
327                                                   540 
328   if(verboseLevel > 2 && buildLambdaTable) {      541   if(verboseLevel > 2 && buildLambdaTable) {
329     out << "      LambdaTable address= " << th << 542     G4cout << "      LambdaTable address= " << theLambdaTable << G4endl;
330     if(theLambdaTable && particle == &part) {  << 543     if(theLambdaTable) { G4cout << (*theLambdaTable) << G4endl; }
331       out << (*theLambdaTable) << G4endl;      << 
332     }                                          << 
333   }                                               544   }
334 }                                                 545 }
335                                                   546 
336 //....oooOO0OOooo........oooOO0OOooo........oo    547 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
337                                                   548 
338 void G4VEmProcess::StartTracking(G4Track* trac    549 void G4VEmProcess::StartTracking(G4Track* track)
339 {                                                 550 {
340   // reset parameters for the new track           551   // reset parameters for the new track
341   currentParticle = track->GetParticleDefiniti    552   currentParticle = track->GetParticleDefinition();
342   theNumberOfInteractionLengthLeft = -1.0;        553   theNumberOfInteractionLengthLeft = -1.0;
343   mfpKinEnergy = DBL_MAX;                      << 554   //currentInteractionLength = -1.0;
344   preStepLambda = 0.0;                         << 555   //  theInitialNumberOfInteractionLength=-1.0;
345                                                << 556   mfpKinEnergy = DBL_MAX; 
346   if(isIon) { massRatio = proton_mass_c2/curre << 
347                                                   557 
348   // forced biasing only for primary particles    558   // forced biasing only for primary particles
349   if(biasManager) {                               559   if(biasManager) {
350     if(0 == track->GetParentID()) {               560     if(0 == track->GetParentID()) {
351       // primary particle                         561       // primary particle
352       biasFlag = true;                            562       biasFlag = true; 
353       biasManager->ResetForcedInteraction();      563       biasManager->ResetForcedInteraction(); 
354     }                                             564     }
355   }                                               565   }
356 }                                                 566 }
357                                                   567 
358 //....oooOO0OOooo........oooOO0OOooo........oo    568 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
359                                                   569 
360 G4double G4VEmProcess::PostStepGetPhysicalInte    570 G4double G4VEmProcess::PostStepGetPhysicalInteractionLength(
361                              const G4Track& tr    571                              const G4Track& track,
362                              G4double   previo    572                              G4double   previousStepSize,
363                              G4ForceCondition*    573                              G4ForceCondition* condition)
364 {                                                 574 {
365   *condition = NotForced;                         575   *condition = NotForced;
366   G4double x = DBL_MAX;                           576   G4double x = DBL_MAX;
367                                                   577 
368   DefineMaterial(track.GetMaterialCutsCouple() << 
369   preStepKinEnergy = track.GetKineticEnergy();    578   preStepKinEnergy = track.GetKineticEnergy();
370   const G4double scaledEnergy = preStepKinEner << 579   DefineMaterial(track.GetMaterialCutsCouple());
371   SelectModel(scaledEnergy, currentCoupleIndex << 580   SelectModel(preStepKinEnergy, currentCoupleIndex);
372   /*                                           << 581 
373   G4cout << "PostStepGetPhysicalInteractionLen << 582   if(!currentModel->IsActive(preStepKinEnergy)) { 
374          << "  couple: " << currentCouple << G << 
375   */                                           << 
376   if(!currentModel->IsActive(scaledEnergy)) {  << 
377     theNumberOfInteractionLengthLeft = -1.0;   << 
378     currentInteractionLength = DBL_MAX;           583     currentInteractionLength = DBL_MAX;
379     mfpKinEnergy = DBL_MAX;                    << 
380     preStepLambda = 0.0;                       << 
381     return x;                                     584     return x; 
382   }                                               585   }
383                                                   586  
384   // forced biasing only for primary particles    587   // forced biasing only for primary particles
385   if(biasManager) {                               588   if(biasManager) {
386     if(0 == track.GetParentID()) {                589     if(0 == track.GetParentID()) {
387       if(biasFlag &&                              590       if(biasFlag && 
388          biasManager->ForcedInteractionRegion( << 591    biasManager->ForcedInteractionRegion(currentCoupleIndex)) {
389         return biasManager->GetStepLimit((G4in << 592         return biasManager->GetStepLimit(currentCoupleIndex, previousStepSize);
390       }                                           593       }
391     }                                             594     }
392   }                                               595   }
393                                                   596 
394   // compute mean free path                       597   // compute mean free path
                                                   >> 598   if(preStepKinEnergy < mfpKinEnergy) {
                                                   >> 599     if (integral) { ComputeIntegralLambda(preStepKinEnergy); }
                                                   >> 600     else { preStepLambda = GetCurrentLambda(preStepKinEnergy); }
                                                   >> 601 
                                                   >> 602     // zero cross section
                                                   >> 603     if(preStepLambda <= 0.0) { 
                                                   >> 604       theNumberOfInteractionLengthLeft = -1.0;
                                                   >> 605       currentInteractionLength = DBL_MAX;
                                                   >> 606     }
                                                   >> 607   }
395                                                   608 
396   ComputeIntegralLambda(preStepKinEnergy, trac << 609   // non-zero cross section
397                                                << 610   if(preStepLambda > 0.0) { 
398   // zero cross section                        << 
399   if(preStepLambda <= 0.0) {                   << 
400     theNumberOfInteractionLengthLeft = -1.0;   << 
401     currentInteractionLength = DBL_MAX;        << 
402                                                << 
403   } else {                                     << 
404                                                   611 
405     // non-zero cross section                  << 
406     if (theNumberOfInteractionLengthLeft < 0.0    612     if (theNumberOfInteractionLengthLeft < 0.0) {
407                                                   613 
408       // beggining of tracking (or just after     614       // beggining of tracking (or just after DoIt of this process)
409       theNumberOfInteractionLengthLeft = -G4Lo << 615       // ResetNumberOfInteractionLengthLeft();
                                                   >> 616 
                                                   >> 617       theNumberOfInteractionLengthLeft =  -G4Log( G4UniformRand() );
410       theInitialNumberOfInteractionLength = th    618       theInitialNumberOfInteractionLength = theNumberOfInteractionLengthLeft; 
411                                                   619 
412     } else {                                   << 620     } else if(currentInteractionLength < DBL_MAX) {
413                                                   621 
                                                   >> 622       // subtract NumberOfInteractionLengthLeft using previous step
414       theNumberOfInteractionLengthLeft -=         623       theNumberOfInteractionLengthLeft -= 
415         previousStepSize/currentInteractionLen << 624   previousStepSize/currentInteractionLength;
416       theNumberOfInteractionLengthLeft =       << 625       //SubtractNumberOfInteractionLengthLeft(previousStepSize);
417         std::max(theNumberOfInteractionLengthL << 626       if(theNumberOfInteractionLengthLeft < 0.) {
                                                   >> 627   theNumberOfInteractionLengthLeft = 0.0;
                                                   >> 628       }
418     }                                             629     }
419                                                   630 
420     // new mean free path and step limit for t    631     // new mean free path and step limit for the next step
421     currentInteractionLength = 1.0/preStepLamb    632     currentInteractionLength = 1.0/preStepLambda;
422     x = theNumberOfInteractionLengthLeft * cur    633     x = theNumberOfInteractionLengthLeft * currentInteractionLength;
423   }                                            << 634 #ifdef G4VERBOSE
424   return x;                                    << 635     if (verboseLevel>2){
425 }                                              << 636       G4cout << "G4VEmProcess::PostStepGetPhysicalInteractionLength ";
426                                                << 637       G4cout << "[ " << GetProcessName() << "]" << G4endl; 
427 //....oooOO0OOooo........oooOO0OOooo........oo << 638       G4cout << " for " << currentParticle->GetParticleName() 
428                                                << 639        << " in Material  " <<  currentMaterial->GetName()
429 void G4VEmProcess::ComputeIntegralLambda(G4dou << 640        << " Ekin(MeV)= " << preStepKinEnergy/MeV 
430 {                                              << 641        <<G4endl;
431   if (fXSType == fEmNoIntegral) {              << 642       G4cout << " MeanFreePath = " << currentInteractionLength/cm << "[cm]" 
432     preStepLambda = GetCurrentLambda(e, LogEki << 643        << " InteractionLength= " << x/cm <<"[cm] " <<G4endl;
433                                                << 
434   } else if (fXSType == fEmIncreasing) {       << 
435     if(e*invLambdaFactor < mfpKinEnergy) {     << 
436       preStepLambda = GetCurrentLambda(e, LogE << 
437       mfpKinEnergy = (preStepLambda > 0.0) ? e << 
438     }                                          << 
439                                                << 
440   } else if(fXSType == fEmDecreasing) {        << 
441     if(e < mfpKinEnergy) {                     << 
442       const G4double e1 = e*lambdaFactor;      << 
443       preStepLambda = GetCurrentLambda(e1);    << 
444       mfpKinEnergy = e1;                       << 
445     }                                          << 
446                                                << 
447   } else if(fXSType == fEmOnePeak) {           << 
448     const G4double epeak = (*theEnergyOfCrossS << 
449     if(e <= epeak) {                           << 
450       if(e*invLambdaFactor < mfpKinEnergy) {   << 
451         preStepLambda = GetCurrentLambda(e, Lo << 
452         mfpKinEnergy = (preStepLambda > 0.0) ? << 
453       }                                        << 
454     } else if(e < mfpKinEnergy) {              << 
455       const G4double e1 = std::max(epeak, e*la << 
456       preStepLambda = GetCurrentLambda(e1);    << 
457       mfpKinEnergy = e1;                       << 
458     }                                             644     }
459   } else {                                     << 645 #endif
460     preStepLambda = GetCurrentLambda(e, LogEki << 
461   }                                               646   }
                                                   >> 647   return x;
462 }                                                 648 }
463                                                   649 
464 //....oooOO0OOooo........oooOO0OOooo........oo    650 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
465                                                   651 
466 G4VParticleChange* G4VEmProcess::PostStepDoIt(    652 G4VParticleChange* G4VEmProcess::PostStepDoIt(const G4Track& track,
467                                                   653                                               const G4Step& step)
468 {                                                 654 {
469   // clear number of interaction lengths in an << 655   // In all cases clear number of interaction lengths
470   theNumberOfInteractionLengthLeft = -1.0;        656   theNumberOfInteractionLengthLeft = -1.0;
471   mfpKinEnergy = DBL_MAX;                      << 657   mfpKinEnergy = DBL_MAX; 
472                                                   658 
473   fParticleChange.InitializeForPostStep(track)    659   fParticleChange.InitializeForPostStep(track);
474                                                   660 
475   // Do not make anything if particle is stopp    661   // Do not make anything if particle is stopped, the annihilation then
476   // should be performed by the AtRestDoIt!       662   // should be performed by the AtRestDoIt!
477   if (track.GetTrackStatus() == fStopButAlive)    663   if (track.GetTrackStatus() == fStopButAlive) { return &fParticleChange; }
478                                                   664 
479   const G4double finalT = track.GetKineticEner << 665   G4double finalT = track.GetKineticEnergy();
480                                                   666 
481   // forced process - should happen only once     667   // forced process - should happen only once per track
482   if(biasFlag) {                                  668   if(biasFlag) {
483     if(biasManager->ForcedInteractionRegion((G << 669     if(biasManager->ForcedInteractionRegion(currentCoupleIndex)) {
484       biasFlag = false;                           670       biasFlag = false;
485     }                                             671     }
486   }                                               672   }
487                                                   673 
488   // check active and select model             << 
489   const G4double scaledEnergy = finalT*massRat << 
490   SelectModel(scaledEnergy, currentCoupleIndex << 
491   if(!currentModel->IsActive(scaledEnergy)) {  << 
492                                                << 
493   // Integral approach                            674   // Integral approach
494   if (fXSType != fEmNoIntegral) {              << 675   if (integral) {
495     const G4double logFinalT =                 << 676     G4double lx = GetLambda(finalT, currentCouple);
496       track.GetDynamicParticle()->GetLogKineti << 677     if(preStepLambda<lx && 1 < verboseLevel) {
497     const G4double lx = std::max(GetCurrentLam << 
498 #ifdef G4VERBOSE                               << 
499     if(preStepLambda < lx && 1 < verboseLevel) << 
500       G4cout << "WARNING: for " << currentPart    678       G4cout << "WARNING: for " << currentParticle->GetParticleName() 
501              << " and " << GetProcessName() << << 679              << " and " << GetProcessName()
502              << " preLambda= " << preStepLambd << 680              << " E(MeV)= " << finalT/MeV
503              << " < " << lx << " (postLambda)  << 681              << " preLambda= " << preStepLambda << " < " 
                                                   >> 682        << lx << " (postLambda) "
                                                   >> 683        << G4endl;  
504     }                                             684     }
505 #endif                                         << 685 
506     // if false interaction then use new cross << 686     if(preStepLambda*G4UniformRand() > lx) {
507     // if both values are zero - no interactio << 687       ClearNumberOfInteractionLengthLeft();
508     if(preStepLambda*G4UniformRand() >= lx) {  << 
509       return &fParticleChange;                    688       return &fParticleChange;
510     }                                             689     }
511   }                                               690   }
512                                                   691 
                                                   >> 692   SelectModel(finalT, currentCoupleIndex);
                                                   >> 693   if(!currentModel->IsActive(finalT)) { return &fParticleChange; }
                                                   >> 694 
513   // define new weight for primary and seconda    695   // define new weight for primary and secondaries
514   G4double weight = fParticleChange.GetParentW    696   G4double weight = fParticleChange.GetParentWeight();
515   if(weightFlag) {                                697   if(weightFlag) { 
516     weight /= biasFactor;                         698     weight /= biasFactor; 
517     fParticleChange.ProposeWeight(weight);        699     fParticleChange.ProposeWeight(weight);
518   }                                               700   }
519                                                   701   
520 #ifdef G4VERBOSE                               << 702   /*  
521   if(1 < verboseLevel) {                       << 703   if(0 < verboseLevel) {
522     G4cout << "G4VEmProcess::PostStepDoIt: Sam    704     G4cout << "G4VEmProcess::PostStepDoIt: Sample secondary; E= "
523            << finalT/MeV                          705            << finalT/MeV
524            << " MeV; model= (" << currentModel    706            << " MeV; model= (" << currentModel->LowEnergyLimit()
525            << ", " <<  currentModel->HighEnerg    707            << ", " <<  currentModel->HighEnergyLimit() << ")"
526            << G4endl;                             708            << G4endl;
527   }                                               709   }
528 #endif                                         << 710   */
529                                                   711 
530   // sample secondaries                           712   // sample secondaries
531   secParticles.clear();                           713   secParticles.clear();
532   currentModel->SampleSecondaries(&secParticle    714   currentModel->SampleSecondaries(&secParticles, 
533                                   currentCoupl << 715           currentCouple, 
534                                   track.GetDyn << 716           track.GetDynamicParticle(),
535                                   (*theCuts)[c << 717           (*theCuts)[currentCoupleIndex]);
536                                                   718 
537   G4int num0 = (G4int)secParticles.size();     << 719   G4int num0 = secParticles.size();
538                                                   720 
539   // splitting or Russian roulette                721   // splitting or Russian roulette
540   if(biasManager) {                               722   if(biasManager) {
541     if(biasManager->SecondaryBiasingRegion((G4 << 723     if(biasManager->SecondaryBiasingRegion(currentCoupleIndex)) {
542       G4double eloss = 0.0;                       724       G4double eloss = 0.0;
543       weight *= biasManager->ApplySecondaryBia    725       weight *= biasManager->ApplySecondaryBiasing(
544         secParticles, track, currentModel, &fP    726         secParticles, track, currentModel, &fParticleChange, eloss, 
545         (G4int)currentCoupleIndex, (*theCuts)[ << 727         currentCoupleIndex, (*theCuts)[currentCoupleIndex],
546         step.GetPostStepPoint()->GetSafety()); << 728   step.GetPostStepPoint()->GetSafety());
547       if(eloss > 0.0) {                           729       if(eloss > 0.0) {
548         eloss += fParticleChange.GetLocalEnerg << 730   eloss += fParticleChange.GetLocalEnergyDeposit();
549         fParticleChange.ProposeLocalEnergyDepo    731         fParticleChange.ProposeLocalEnergyDeposit(eloss);
550       }                                           732       }
551     }                                             733     }
552   }                                               734   }
553                                                   735 
554   // save secondaries                             736   // save secondaries
555   G4int num = (G4int)secParticles.size();      << 737   G4int num = secParticles.size();
556   if(num > 0) {                                   738   if(num > 0) {
557                                                   739 
558     fParticleChange.SetNumberOfSecondaries(num    740     fParticleChange.SetNumberOfSecondaries(num);
559     G4double edep = fParticleChange.GetLocalEn    741     G4double edep = fParticleChange.GetLocalEnergyDeposit();
560     G4double time = track.GetGlobalTime();        742     G4double time = track.GetGlobalTime();
561                                                << 
562     G4int n1(0), n2(0);                        << 
563     if(num0 > mainSecondaries) {               << 
564       currentModel->FillNumberOfSecondaries(n1 << 
565     }                                          << 
566                                                   743      
567     for (G4int i=0; i<num; ++i) {                 744     for (G4int i=0; i<num; ++i) {
568       G4DynamicParticle* dp = secParticles[i]; << 745       if (secParticles[i]) {
569       if (nullptr != dp) {                     << 746         G4DynamicParticle* dp = secParticles[i];
570         const G4ParticleDefinition* p = dp->Ge    747         const G4ParticleDefinition* p = dp->GetParticleDefinition();
571         G4double e = dp->GetKineticEnergy();      748         G4double e = dp->GetKineticEnergy();
572         G4bool good = true;                       749         G4bool good = true;
573         if(applyCuts) {                           750         if(applyCuts) {
574           if (p == theGamma) {                 << 751     if (p == theGamma) {
575             if (e < (*theCutsGamma)[currentCou << 752       if (e < (*theCutsGamma)[currentCoupleIndex]) { good = false; }
576                                                   753 
577           } else if (p == theElectron) {       << 754     } else if (p == theElectron) {
578             if (e < (*theCutsElectron)[current << 755       if (e < (*theCutsElectron)[currentCoupleIndex]) { good = false; }
579                                                   756 
580           } else if (p == thePositron) {       << 757     } else if (p == thePositron) {
581             if (electron_mass_c2 < (*theCutsGa << 758       if (electron_mass_c2 < (*theCutsGamma)[currentCoupleIndex] &&
582                 e < (*theCutsPositron)[current << 759     e < (*theCutsPositron)[currentCoupleIndex]) {
583               good = false;                    << 760         good = false;
584               e += 2.0*electron_mass_c2;       << 761         e += 2.0*electron_mass_c2;
585             }                                  << 762       }
586           }                                    << 763     }
587           // added secondary if it is good     << 764     // added secondary if it is good
588         }                                         765         }
589         if (good) {                               766         if (good) { 
590           G4Track* t = new G4Track(dp, time, t    767           G4Track* t = new G4Track(dp, time, track.GetPosition());
591           t->SetTouchableHandle(track.GetTouch    768           t->SetTouchableHandle(track.GetTouchableHandle());
592           if (biasManager) {                   << 769           t->SetWeight(weight);
593             t->SetWeight(weight * biasManager- << 
594           } else {                             << 
595             t->SetWeight(weight);              << 
596           }                                    << 
597           pParticleChange->AddSecondary(t);       770           pParticleChange->AddSecondary(t);
598                                                   771 
599           // define type of secondary          << 772     // define type of secondary
600           if(i < mainSecondaries) {            << 773           if(i < mainSecondaries) { t->SetCreatorModelIndex(secID); }
601             t->SetCreatorModelID(secID);       << 774           else if(i < num0) {
602             if(GetProcessSubType() == fCompton << 775             if(p == theGamma) { 
603               t->SetCreatorModelID(_ComptonGam << 776         t->SetCreatorModelIndex(fluoID);
604             }                                  << 777       } else {
605           } else if(i < mainSecondaries + n1)  << 778         t->SetCreatorModelIndex(augerID);
606             t->SetCreatorModelID(tripletID);   << 779       }
607           } else if(i < mainSecondaries + n1 + << 780     } else {
608             t->SetCreatorModelID(_IonRecoil);  << 781       t->SetCreatorModelIndex(biasID);
609           } else {                             << 782     }
610             if(i < num0) {                     << 783  
611               if(p == theGamma) {              << 784           //G4cout << "Secondary(post step) has weight " << t->GetWeight() 
612                 t->SetCreatorModelID(fluoID);  << 785     // << ", Ekin= " << t->GetKineticEnergy()/MeV << " MeV" <<G4endl;
613               } else {                         << 
614                 t->SetCreatorModelID(augerID); << 
615               }                                << 
616             } else {                           << 
617               t->SetCreatorModelID(biasID);    << 
618             }                                  << 
619           }                                    << 
620           /*                                   << 
621           G4cout << "Secondary(post step) has  << 
622                  << ", Ekin= " << t->GetKineti << 
623                  << GetProcessName() << " fluo << 
624                  << " augerID= " << augerID << << 
625           */                                   << 
626         } else {                                  786         } else {
627           delete dp;                           << 787     delete dp;
628           edep += e;                           << 788     edep += e;
629         }                                      << 789   }
630       }                                           790       } 
631     }                                             791     }
632     fParticleChange.ProposeLocalEnergyDeposit(    792     fParticleChange.ProposeLocalEnergyDeposit(edep);
633   }                                               793   }
634                                                   794 
635   if(0.0 == fParticleChange.GetProposedKinetic    795   if(0.0 == fParticleChange.GetProposedKineticEnergy() &&
636      fAlive == fParticleChange.GetTrackStatus(    796      fAlive == fParticleChange.GetTrackStatus()) {
637     if(particle->GetProcessManager()->GetAtRes    797     if(particle->GetProcessManager()->GetAtRestProcessVector()->size() > 0)
638          { fParticleChange.ProposeTrackStatus(    798          { fParticleChange.ProposeTrackStatus(fStopButAlive); }
639     else { fParticleChange.ProposeTrackStatus(    799     else { fParticleChange.ProposeTrackStatus(fStopAndKill); }
640   }                                               800   }
641                                                   801 
                                                   >> 802   //  ClearNumberOfInteractionLengthLeft();
642   return &fParticleChange;                        803   return &fParticleChange;
643 }                                                 804 }
644                                                   805 
645 //....oooOO0OOooo........oooOO0OOooo........oo    806 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
646                                                   807 
647 G4bool G4VEmProcess::StorePhysicsTable(const G    808 G4bool G4VEmProcess::StorePhysicsTable(const G4ParticleDefinition* part,
648                                        const G << 809                const G4String& directory,
649                                        G4bool  << 810                      G4bool ascii)
650 {                                                 811 {
651   if(!isTheMaster || part != particle) { retur << 812   G4bool yes = true;
652   if(G4EmTableUtil::StoreTable(this, part, the << 813 
653              directory, "Lambda",              << 814   if ( theLambdaTable && part == particle) {
654                                verboseLevel, a << 815     const G4String name = 
655      G4EmTableUtil::StoreTable(this, part, the << 816       GetPhysicsTableFileName(part,directory,"Lambda",ascii);
656              directory, "LambdaPrim",          << 817     yes = theLambdaTable->StorePhysicsTable(name,ascii);
657                                verboseLevel, a << 818 
658      return true;                              << 819     if ( yes ) {
                                                   >> 820       G4cout << "Physics table is stored for " << particle->GetParticleName()
                                                   >> 821              << " and process " << GetProcessName()
                                                   >> 822        << " in the directory <" << directory
                                                   >> 823        << "> " << G4endl;
                                                   >> 824     } else {
                                                   >> 825       G4cout << "Fail to store Physics Table for " 
                                                   >> 826        << particle->GetParticleName()
                                                   >> 827              << " and process " << GetProcessName()
                                                   >> 828        << " in the directory <" << directory
                                                   >> 829        << "> " << G4endl;
                                                   >> 830     }
659   }                                               831   }
660   return false;                                << 832   if ( theLambdaTablePrim && part == particle) {
                                                   >> 833     const G4String name = 
                                                   >> 834       GetPhysicsTableFileName(part,directory,"LambdaPrim",ascii);
                                                   >> 835     yes = theLambdaTablePrim->StorePhysicsTable(name,ascii);
                                                   >> 836 
                                                   >> 837     if ( yes ) {
                                                   >> 838       G4cout << "Physics table prim is stored for " 
                                                   >> 839        << particle->GetParticleName()
                                                   >> 840              << " and process " << GetProcessName()
                                                   >> 841        << " in the directory <" << directory
                                                   >> 842        << "> " << G4endl;
                                                   >> 843     } else {
                                                   >> 844       G4cout << "Fail to store Physics Table Prim for " 
                                                   >> 845        << particle->GetParticleName()
                                                   >> 846              << " and process " << GetProcessName()
                                                   >> 847        << " in the directory <" << directory
                                                   >> 848        << "> " << G4endl;
                                                   >> 849     }
                                                   >> 850   }
                                                   >> 851   return yes;
661 }                                                 852 }
662                                                   853 
663 //....oooOO0OOooo........oooOO0OOooo........oo    854 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
664                                                   855 
665 G4bool G4VEmProcess::RetrievePhysicsTable(cons    856 G4bool G4VEmProcess::RetrievePhysicsTable(const G4ParticleDefinition* part,
666                                           cons << 857                     const G4String& directory,
667                                           G4bo << 858             G4bool ascii)
668 {                                                 859 {
669   if(!isTheMaster || part != particle) { retur << 860   if(1 < verboseLevel) {
                                                   >> 861     G4cout << "G4VEmProcess::RetrievePhysicsTable() for "
                                                   >> 862            << part->GetParticleName() << " and process "
                                                   >> 863      << GetProcessName() << G4endl;
                                                   >> 864   }
670   G4bool yes = true;                              865   G4bool yes = true;
                                                   >> 866 
                                                   >> 867   if((!buildLambdaTable && minKinEnergyPrim > maxKinEnergy) 
                                                   >> 868      || particle != part) { return yes; }
                                                   >> 869 
                                                   >> 870   const G4String particleName = part->GetParticleName();
                                                   >> 871   G4String filename;
                                                   >> 872 
671   if(buildLambdaTable) {                          873   if(buildLambdaTable) {
672     yes = G4EmTableUtil::RetrieveTable(this, p << 874     filename = GetPhysicsTableFileName(part,directory,"Lambda",ascii);
673                                        "Lambda << 875     yes = G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTable,
674                                        ascii,  << 876                  filename,ascii);
675   }                                            << 877     if ( yes ) {
676   if(yes && minKinEnergyPrim < maxKinEnergy) { << 878       if (0 < verboseLevel) {
677     yes = G4EmTableUtil::RetrieveTable(this, p << 879   G4cout << "Lambda table for " << particleName 
678                                        "Lambda << 880          << " is Retrieved from <"
679                                        ascii,  << 881          << filename << ">"
                                                   >> 882          << G4endl;
                                                   >> 883       }
                                                   >> 884       if(lManager->SplineFlag()) {
                                                   >> 885   size_t n = theLambdaTable->length();
                                                   >> 886   for(size_t i=0; i<n; ++i) {
                                                   >> 887     if((* theLambdaTable)[i]) {
                                                   >> 888       (* theLambdaTable)[i]->SetSpline(true);
                                                   >> 889     }
                                                   >> 890   }
                                                   >> 891       }
                                                   >> 892     } else {
                                                   >> 893       if (1 < verboseLevel) {
                                                   >> 894   G4cout << "Lambda table for " << particleName << " in file <"
                                                   >> 895          << filename << "> is not exist"
                                                   >> 896          << G4endl;
                                                   >> 897       }
                                                   >> 898     }
                                                   >> 899   }
                                                   >> 900   if(minKinEnergyPrim < maxKinEnergy) {
                                                   >> 901     filename = GetPhysicsTableFileName(part,directory,"LambdaPrim",ascii);
                                                   >> 902     yes = G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTablePrim,
                                                   >> 903                  filename,ascii);
                                                   >> 904     if ( yes ) {
                                                   >> 905       if (0 < verboseLevel) {
                                                   >> 906   G4cout << "Lambda table prim for " << particleName 
                                                   >> 907          << " is Retrieved from <"
                                                   >> 908          << filename << ">"
                                                   >> 909          << G4endl;
                                                   >> 910       }
                                                   >> 911       if(lManager->SplineFlag()) {
                                                   >> 912   size_t n = theLambdaTablePrim->length();
                                                   >> 913   for(size_t i=0; i<n; ++i) {
                                                   >> 914     if((* theLambdaTablePrim)[i]) {
                                                   >> 915       (* theLambdaTablePrim)[i]->SetSpline(true);
                                                   >> 916     }
                                                   >> 917   }
                                                   >> 918       }
                                                   >> 919     } else {
                                                   >> 920       if (1 < verboseLevel) {
                                                   >> 921   G4cout << "Lambda table prim for " << particleName << " in file <"
                                                   >> 922          << filename << "> is not exist"
                                                   >> 923          << G4endl;
                                                   >> 924       }
                                                   >> 925     }
680   }                                               926   }
                                                   >> 927 
681   return yes;                                     928   return yes;
682 }                                                 929 }
683                                                   930 
684 //....oooOO0OOooo........oooOO0OOooo........oo    931 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
685                                                   932 
686 G4double G4VEmProcess::GetCrossSection(G4doubl << 933 G4double 
687                                        const G << 934 G4VEmProcess::CrossSectionPerVolume(G4double kineticEnergy,
                                                   >> 935             const G4MaterialCutsCouple* couple)
688 {                                                 936 {
689   CurrentSetup(couple, kinEnergy);             << 937   // Cross section per atom is calculated
690   return GetCurrentLambda(kinEnergy, G4Log(kin << 938   DefineMaterial(couple);
                                                   >> 939   G4double cross = 0.0;
                                                   >> 940   if(buildLambdaTable && theLambdaTable) {
                                                   >> 941     cross = GetCurrentLambda(kineticEnergy);
                                                   >> 942   } else {
                                                   >> 943     SelectModel(kineticEnergy, currentCoupleIndex);
                                                   >> 944     cross = fFactor*currentModel->CrossSectionPerVolume(currentMaterial,
                                                   >> 945               currentParticle,
                                                   >> 946               kineticEnergy);
                                                   >> 947   }
                                                   >> 948 
                                                   >> 949   if(cross < 0.0) { cross = 0.0; }
                                                   >> 950   return cross;
691 }                                                 951 }
692                                                   952 
693 //....oooOO0OOooo........oooOO0OOooo........oo    953 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
694                                                   954 
695 G4double G4VEmProcess::GetMeanFreePath(const G    955 G4double G4VEmProcess::GetMeanFreePath(const G4Track& track,
696                                        G4doubl << 956                G4double,
697                                        G4Force << 957                G4ForceCondition* condition)
698 {                                                 958 {
699   *condition = NotForced;                         959   *condition = NotForced;
700   return G4VEmProcess::MeanFreePath(track);       960   return G4VEmProcess::MeanFreePath(track);
701 }                                                 961 }
702                                                   962 
703 //....oooOO0OOooo........oooOO0OOooo........oo    963 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
704                                                   964 
705 G4double                                       << 965 G4double G4VEmProcess::MeanFreePath(const G4Track& track)
706 G4VEmProcess::ComputeCrossSectionPerAtom(G4dou << 
707                                          G4dou << 
708 {                                                 966 {
709   SelectModel(kinEnergy, currentCoupleIndex);  << 967   DefineMaterial(track.GetMaterialCutsCouple());
710   return (currentModel) ?                      << 968   preStepLambda = GetCurrentLambda(track.GetKineticEnergy());
711     currentModel->ComputeCrossSectionPerAtom(c << 969   G4double x = DBL_MAX;
712                                              Z << 970   if(0.0 < preStepLambda) { x = 1.0/preStepLambda; }
                                                   >> 971   return x;
713 }                                                 972 }
714                                                   973 
715 //....oooOO0OOooo........oooOO0OOooo........oo    974 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
716                                                   975 
717 G4PhysicsVector*                               << 976 G4double 
718 G4VEmProcess::LambdaPhysicsVector(const G4Mate << 977 G4VEmProcess::ComputeCrossSectionPerAtom(G4double kineticEnergy, 
                                                   >> 978            G4double Z, G4double A, G4double cut)
719 {                                                 979 {
720   DefineMaterial(couple);                      << 980   SelectModel(kineticEnergy, currentCoupleIndex);
721   G4PhysicsVector* newv = new G4PhysicsLogVect << 981   G4double x = 0.0;
722                                                << 982   if(currentModel) {
723   return newv;                                 << 983     x = currentModel->ComputeCrossSectionPerAtom(currentParticle,kineticEnergy,
                                                   >> 984              Z,A,cut);
                                                   >> 985   }
                                                   >> 986   return x;
724 }                                                 987 }
725                                                   988 
726 //....oooOO0OOooo........oooOO0OOooo........oo    989 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
727                                                   990 
728 const G4Element* G4VEmProcess::GetCurrentEleme << 991 void G4VEmProcess::FindLambdaMax()
729 {                                                 992 {
730   return (nullptr != currentModel) ?           << 993   if(1 < verboseLevel) {
731     currentModel->GetCurrentElement(currentMat << 994     G4cout << "### G4VEmProcess::FindLambdaMax: " 
                                                   >> 995      << particle->GetParticleName() 
                                                   >> 996            << " and process " << GetProcessName() << "  " << G4endl; 
                                                   >> 997   }
                                                   >> 998   size_t n = theLambdaTable->length();
                                                   >> 999   G4PhysicsVector* pv;
                                                   >> 1000   G4double e, ss, emax, smax;
                                                   >> 1001 
                                                   >> 1002   size_t i;
                                                   >> 1003 
                                                   >> 1004   // first loop on existing vectors
                                                   >> 1005   for (i=0; i<n; ++i) {
                                                   >> 1006     pv = (*theLambdaTable)[i];
                                                   >> 1007     if(pv) {
                                                   >> 1008       size_t nb = pv->GetVectorLength();
                                                   >> 1009       emax = DBL_MAX;
                                                   >> 1010       smax = 0.0;
                                                   >> 1011       if(nb > 0) {
                                                   >> 1012   for (size_t j=0; j<nb; ++j) {
                                                   >> 1013     e = pv->Energy(j);
                                                   >> 1014     ss = (*pv)(j);
                                                   >> 1015     if(ss > smax) {
                                                   >> 1016       smax = ss;
                                                   >> 1017       emax = e;
                                                   >> 1018     }
                                                   >> 1019   }
                                                   >> 1020       }
                                                   >> 1021       theEnergyOfCrossSectionMax[i] = emax;
                                                   >> 1022       theCrossSectionMax[i] = smax;
                                                   >> 1023       if(1 < verboseLevel) {
                                                   >> 1024   G4cout << "For " << particle->GetParticleName() 
                                                   >> 1025          << " Max CS at i= " << i << " emax(MeV)= " << emax/MeV
                                                   >> 1026          << " lambda= " << smax << G4endl;
                                                   >> 1027       }
                                                   >> 1028     }
                                                   >> 1029   }
                                                   >> 1030   // second loop using base materials
                                                   >> 1031   for (i=0; i<n; ++i) {
                                                   >> 1032     pv = (*theLambdaTable)[i];
                                                   >> 1033     if(!pv){
                                                   >> 1034       G4int j = (*theDensityIdx)[i];
                                                   >> 1035       theEnergyOfCrossSectionMax[i] = theEnergyOfCrossSectionMax[j];
                                                   >> 1036       theCrossSectionMax[i] = (*theDensityFactor)[i]*theCrossSectionMax[j];
                                                   >> 1037     }
                                                   >> 1038   }
732 }                                                 1039 }
733                                                   1040 
734 //....oooOO0OOooo........oooOO0OOooo........oo    1041 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
735                                                   1042 
736 const G4Element* G4VEmProcess::GetTargetElemen << 1043 G4PhysicsVector* 
                                                   >> 1044 G4VEmProcess::LambdaPhysicsVector(const G4MaterialCutsCouple*)
737 {                                                 1045 {
738   return (nullptr != currentModel) ?           << 1046   G4PhysicsVector* v = 
739     currentModel->GetCurrentElement(currentMat << 1047     new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nLambdaBins);
                                                   >> 1048   v->SetSpline(lManager->SplineFlag());
                                                   >> 1049   return v;
740 }                                                 1050 }
741                                                   1051 
742 //....oooOO0OOooo........oooOO0OOooo........oo    1052 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
743                                                   1053 
744 const G4Isotope* G4VEmProcess::GetTargetIsotop << 1054 const G4Element* G4VEmProcess::GetCurrentElement() const
745 {                                                 1055 {
746   return (nullptr != currentModel) ?           << 1056   const G4Element* elm = 0;
747     currentModel->GetCurrentIsotope(GetCurrent << 1057   if(currentModel) {elm = currentModel->GetCurrentElement(); }
                                                   >> 1058   return elm;
748 }                                                 1059 }
749                                                   1060 
750 //....oooOO0OOooo........oooOO0OOooo........oo    1061 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
751                                                   1062 
752 void G4VEmProcess::SetCrossSectionBiasingFacto    1063 void G4VEmProcess::SetCrossSectionBiasingFactor(G4double f, G4bool flag)
753 {                                                 1064 {
754   if(f > 0.0) {                                   1065   if(f > 0.0) { 
755     biasFactor = f;                               1066     biasFactor = f; 
756     weightFlag = flag;                            1067     weightFlag = flag;
757     if(1 < verboseLevel) {                        1068     if(1 < verboseLevel) {
758       G4cout << "### SetCrossSectionBiasingFac    1069       G4cout << "### SetCrossSectionBiasingFactor: for " 
759              << particle->GetParticleName()    << 1070        << particle->GetParticleName() 
760              << " and process " << GetProcessN << 1071        << " and process " << GetProcessName()
761              << " biasFactor= " << f << " weig << 1072        << " biasFactor= " << f << " weightFlag= " << flag 
762              << G4endl;                        << 1073        << G4endl; 
763     }                                             1074     }
764   }                                               1075   }
765 }                                                 1076 }
766                                                   1077 
767 //....oooOO0OOooo........oooOO0OOooo........oo    1078 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
768                                                   1079 
769 void                                              1080 void 
770 G4VEmProcess::ActivateForcedInteraction(G4doub    1081 G4VEmProcess::ActivateForcedInteraction(G4double length, const G4String& r,
771                                         G4bool << 1082           G4bool flag)
772 {                                                 1083 {
773   if(nullptr == biasManager) { biasManager = n << 1084   if(!biasManager) { biasManager = new G4EmBiasingManager(); }
774   if(1 < verboseLevel) {                          1085   if(1 < verboseLevel) {
775     G4cout << "### ActivateForcedInteraction:     1086     G4cout << "### ActivateForcedInteraction: for " 
776            << particle->GetParticleName()      << 1087      << particle->GetParticleName() 
777            << " and process " << GetProcessNam << 1088      << " and process " << GetProcessName()
778            << " length(mm)= " << length/mm     << 1089      << " length(mm)= " << length/mm
779            << " in G4Region <" << r            << 1090      << " in G4Region <" << r 
780            << "> weightFlag= " << flag         << 1091      << "> weightFlag= " << flag 
781            << G4endl;                          << 1092      << G4endl; 
782   }                                               1093   }
783   weightFlag = flag;                              1094   weightFlag = flag;
784   biasManager->ActivateForcedInteraction(lengt    1095   biasManager->ActivateForcedInteraction(length, r);
785 }                                                 1096 }
786                                                   1097 
787 //....oooOO0OOooo........oooOO0OOooo........oo    1098 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
788                                                   1099 
789 void                                              1100 void
790 G4VEmProcess::ActivateSecondaryBiasing(const G    1101 G4VEmProcess::ActivateSecondaryBiasing(const G4String& region,
791                  G4double factor,                 1102                  G4double factor,
792                  G4double energyLimit)            1103                  G4double energyLimit)
793 {                                                 1104 {
794   if (0.0 <= factor) {                            1105   if (0.0 <= factor) {
795                                                   1106 
796     // Range cut can be applied only for e-       1107     // Range cut can be applied only for e-
797     if(0.0 == factor && secondaryParticle != G    1108     if(0.0 == factor && secondaryParticle != G4Electron::Electron())
798       { return; }                                 1109       { return; }
799                                                   1110 
800     if(!biasManager) { biasManager = new G4EmB    1111     if(!biasManager) { biasManager = new G4EmBiasingManager(); }
801     biasManager->ActivateSecondaryBiasing(regi    1112     biasManager->ActivateSecondaryBiasing(region, factor, energyLimit);
802     if(1 < verboseLevel) {                        1113     if(1 < verboseLevel) {
803       G4cout << "### ActivateSecondaryBiasing:    1114       G4cout << "### ActivateSecondaryBiasing: for "
804        << " process " << GetProcessName()         1115        << " process " << GetProcessName()
805        << " factor= " << factor                   1116        << " factor= " << factor
806        << " in G4Region <" << region              1117        << " in G4Region <" << region
807        << "> energyLimit(MeV)= " << energyLimi    1118        << "> energyLimit(MeV)= " << energyLimit/MeV
808        << G4endl;                                 1119        << G4endl;
809     }                                             1120     }
810   }                                               1121   }
811 }                                                 1122 }
812                                                   1123 
813 //....oooOO0OOooo........oooOO0OOooo........oo    1124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
814                                                   1125 
815 void G4VEmProcess::SetLambdaBinning(G4int n)   << 
816 {                                              << 
817   if(5 < n && n < 10000000) {                  << 
818     nLambdaBins = n;                           << 
819     actBinning = true;                         << 
820   } else {                                     << 
821     G4double e = (G4double)n;                  << 
822     PrintWarning("SetLambdaBinning", e);       << 
823   }                                            << 
824 }                                              << 
825                                                << 
826 //....oooOO0OOooo........oooOO0OOooo........oo << 
827                                                << 
828 void G4VEmProcess::SetMinKinEnergy(G4double e)    1126 void G4VEmProcess::SetMinKinEnergy(G4double e)
829 {                                                 1127 {
830   if(1.e-3*eV < e && e < maxKinEnergy) {       << 1128   nLambdaBins = G4lrint(nLambdaBins*std::log(maxKinEnergy/e)
831     nLambdaBins = G4lrint(nLambdaBins*G4Log(ma << 1129       /std::log(maxKinEnergy/minKinEnergy));
832                           /G4Log(maxKinEnergy/ << 1130   minKinEnergy = e;
833     minKinEnergy = e;                          << 
834     actMinKinEnergy = true;                    << 
835   } else { PrintWarning("SetMinKinEnergy", e); << 
836 }                                                 1131 }
837                                                   1132 
838 //....oooOO0OOooo........oooOO0OOooo........oo    1133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
839                                                   1134 
840 void G4VEmProcess::SetMaxKinEnergy(G4double e)    1135 void G4VEmProcess::SetMaxKinEnergy(G4double e)
841 {                                                 1136 {
842   if(minKinEnergy < e && e < 1.e+6*TeV) {      << 1137   nLambdaBins = G4lrint(nLambdaBins*std::log(e/minKinEnergy)
843     nLambdaBins = G4lrint(nLambdaBins*G4Log(e/ << 1138       /std::log(maxKinEnergy/minKinEnergy));
844                           /G4Log(maxKinEnergy/ << 1139   maxKinEnergy = e;
845     maxKinEnergy = e;                          << 
846     actMaxKinEnergy = true;                    << 
847   } else { PrintWarning("SetMaxKinEnergy", e); << 
848 }                                              << 
849                                                << 
850 //....oooOO0OOooo........oooOO0OOooo........oo << 
851                                                << 
852 void G4VEmProcess::SetMinKinEnergyPrim(G4doubl << 
853 {                                              << 
854   if(theParameters->MinKinEnergy() <= e &&     << 
855      e <= theParameters->MaxKinEnergy()) { min << 
856   else { PrintWarning("SetMinKinEnergyPrim", e << 
857 }                                              << 
858                                                << 
859 //....oooOO0OOooo........oooOO0OOooo........oo << 
860                                                << 
861 G4VEmProcess* G4VEmProcess::GetEmProcess(const << 
862 {                                              << 
863   return (nam == GetProcessName()) ? this : nu << 
864 }                                              << 
865                                                << 
866 //....oooOO0OOooo........oooOO0OOooo........oo << 
867                                                << 
868 G4double G4VEmProcess::PolarAngleLimit() const << 
869 {                                              << 
870   return theParameters->MscThetaLimit();       << 
871 }                                              << 
872                                                << 
873 //....oooOO0OOooo........oooOO0OOooo........oo << 
874                                                << 
875 void G4VEmProcess::PrintWarning(G4String tit,  << 
876 {                                              << 
877   G4String ss = "G4VEmProcess::" + tit;        << 
878   G4ExceptionDescription ed;                   << 
879   ed << "Parameter is out of range: " << val   << 
880      << " it will have no effect!\n" << "  Pro << 
881      << GetProcessName() << "  nbins= " << the << 
882      << " Emin(keV)= " << theParameters->MinKi << 
883      << " Emax(GeV)= " << theParameters->MaxKi << 
884   G4Exception(ss, "em0044", JustWarning, ed);  << 
885 }                                              << 
886                                                << 
887 //....oooOO0OOooo........oooOO0OOooo........oo << 
888                                                << 
889 void G4VEmProcess::ProcessDescription(std::ost << 
890 {                                              << 
891   if(nullptr != particle) {                    << 
892     StreamInfo(out, *particle, true);          << 
893   }                                            << 
894 }                                                 1140 }
895                                                   1141 
896 //....oooOO0OOooo........oooOO0OOooo........oo    1142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
897                                                   1143