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