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.1)


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