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


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