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 9.4.p2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
                                                   >>  26 // $Id: G4VEmProcess.cc,v 1.88 2010-08-17 17:36:59 vnivanch Exp $
                                                   >>  27 // GEANT4 tag $Name: geant4-09-04-patch-02 $
                                                   >>  28 //
 26 // -------------------------------------------     29 // -------------------------------------------------------------------
 27 //                                                 30 //
 28 // GEANT4 Class file                               31 // GEANT4 Class file
 29 //                                                 32 //
 30 //                                                 33 //
 31 // File name:     G4VEmProcess                     34 // File name:     G4VEmProcess
 32 //                                                 35 //
 33 // Author:        Vladimir Ivanchenko on base      36 // Author:        Vladimir Ivanchenko on base of Laszlo Urban code
 34 //                                                 37 //
 35 // Creation date: 01.10.2003                       38 // Creation date: 01.10.2003
 36 //                                                 39 //
 37 // Modifications: by V.Ivanchenko              <<  40 // Modifications:
                                                   >>  41 // 30-06-04 make it to be pure discrete process (V.Ivanchenko)
                                                   >>  42 // 30-09-08 optimise integral option (V.Ivanchenko)
                                                   >>  43 // 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivanchenko)
                                                   >>  44 // 11-03-05 Shift verbose level by 1, add applyCuts and killPrimary flags (VI)
                                                   >>  45 // 14-03-05 Update logic PostStepDoIt (V.Ivanchenko)
                                                   >>  46 // 08-04-05 Major optimisation of internal interfaces (V.Ivanchenko)
                                                   >>  47 // 18-04-05 Use G4ParticleChangeForGamma (V.Ivanchenko)
                                                   >>  48 // 25-07-05 Add protection: integral mode only for charged particles (VI)
                                                   >>  49 // 04-09-05 default lambdaFactor 0.8 (V.Ivanchenko)
                                                   >>  50 // 11-01-06 add A to parameters of ComputeCrossSectionPerAtom (VI)
                                                   >>  51 // 12-09-06 add SetModel() (mma)
                                                   >>  52 // 12-04-07 remove double call to Clear model manager (V.Ivanchenko)
                                                   >>  53 // 27-10-07 Virtual functions moved to source (V.Ivanchenko)
                                                   >>  54 // 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko)
                                                   >>  55 // 17-02-10 Added pointer currentParticle (VI)
 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"              << 
 49 #include "G4SystemOfUnits.hh"                  << 
 50 #include "G4ProcessManager.hh"                 << 
 51 #include "G4LossTableManager.hh"                   66 #include "G4LossTableManager.hh"
 52 #include "G4LossTableBuilder.hh"               << 
 53 #include "G4Step.hh"                               67 #include "G4Step.hh"
 54 #include "G4ParticleDefinition.hh"                 68 #include "G4ParticleDefinition.hh"
 55 #include "G4VEmModel.hh"                           69 #include "G4VEmModel.hh"
 56 #include "G4DataVector.hh"                         70 #include "G4DataVector.hh"
 57 #include "G4PhysicsTable.hh"                       71 #include "G4PhysicsTable.hh"
 58 #include "G4EmDataHandler.hh"                  <<  72 #include "G4PhysicsVector.hh"
 59 #include "G4PhysicsLogVector.hh"                   73 #include "G4PhysicsLogVector.hh"
 60 #include "G4VParticleChange.hh"                    74 #include "G4VParticleChange.hh"
 61 #include "G4ProductionCutsTable.hh"                75 #include "G4ProductionCutsTable.hh"
 62 #include "G4Region.hh"                             76 #include "G4Region.hh"
                                                   >>  77 #include "G4RegionStore.hh"
 63 #include "G4Gamma.hh"                              78 #include "G4Gamma.hh"
 64 #include "G4Electron.hh"                           79 #include "G4Electron.hh"
 65 #include "G4Positron.hh"                           80 #include "G4Positron.hh"
 66 #include "G4PhysicsTableHelper.hh"                 81 #include "G4PhysicsTableHelper.hh"
 67 #include "G4EmBiasingManager.hh"               <<  82 #include "G4EmConfigurator.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"                     << 
 74 #include "G4Log.hh"                            << 
 75 #include <iostream>                            << 
 76                                                    83 
 77 //....oooOO0OOooo........oooOO0OOooo........oo     84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 78                                                    85 
 79 G4VEmProcess::G4VEmProcess(const G4String& nam     86 G4VEmProcess::G4VEmProcess(const G4String& name, G4ProcessType type):
 80   G4VDiscreteProcess(name, type)               <<  87   G4VDiscreteProcess(name, type),
                                                   >>  88   secondaryParticle(0),
                                                   >>  89   buildLambdaTable(true),
                                                   >>  90   theLambdaTable(0),
                                                   >>  91   theEnergyOfCrossSectionMax(0),
                                                   >>  92   theCrossSectionMax(0),
                                                   >>  93   integral(false),
                                                   >>  94   applyCuts(false),
                                                   >>  95   startFromNull(true),
                                                   >>  96   useDeexcitation(false),
                                                   >>  97   nDERegions(0),
                                                   >>  98   idxDERegions(0),
                                                   >>  99   currentModel(0),
                                                   >> 100   particle(0),
                                                   >> 101   currentParticle(0),
                                                   >> 102   currentCouple(0)
 81 {                                                 103 {
 82   theParameters = G4EmParameters::Instance();  << 
 83   SetVerboseLevel(1);                             104   SetVerboseLevel(1);
 84                                                   105 
 85   // Size of tables                            << 106   // Size of tables assuming spline
 86   minKinEnergy = 0.1*CLHEP::keV;               << 107   minKinEnergy = 0.1*keV;
 87   maxKinEnergy = 100.0*CLHEP::TeV;             << 108   maxKinEnergy = 10.0*TeV;
                                                   >> 109   nLambdaBins  = 77;
 88                                                   110 
 89   // default lambda factor                        111   // default lambda factor
 90   invLambdaFactor = 1.0/lambdaFactor;          << 112   lambdaFactor  = 0.8;
                                                   >> 113 
                                                   >> 114   // default limit on polar angle
                                                   >> 115   polarAngleLimit = 0.0;
 91                                                   116 
 92   // particle types                               117   // particle types
 93   theGamma = G4Gamma::Gamma();                 << 118   theGamma     = G4Gamma::Gamma();
 94   theElectron = G4Electron::Electron();        << 119   theElectron  = G4Electron::Electron();
 95   thePositron = G4Positron::Positron();        << 120   thePositron  = G4Positron::Positron();
 96                                                   121 
 97   pParticleChange = &fParticleChange;             122   pParticleChange = &fParticleChange;
 98   fParticleChange.SetSecondaryWeightByProcess( << 
 99   secParticles.reserve(5);                        123   secParticles.reserve(5);
100                                                   124 
101   modelManager = new G4EmModelManager();          125   modelManager = new G4EmModelManager();
102   lManager = G4LossTableManager::Instance();   << 126   (G4LossTableManager::Instance())->Register(this);
103   lManager->Register(this);                    << 
104   isTheMaster = lManager->IsMaster();          << 
105   G4LossTableBuilder* bld = lManager->GetTable << 
106   theDensityFactor = bld->GetDensityFactors(); << 
107   theDensityIdx = bld->GetCoupleIndexes();     << 
108 }                                                 127 }
109                                                   128 
110 //....oooOO0OOooo........oooOO0OOooo........oo    129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111                                                   130 
112 G4VEmProcess::~G4VEmProcess()                     131 G4VEmProcess::~G4VEmProcess()
113 {                                                 132 {
114   if(isTheMaster) {                            << 133   if(1 < verboseLevel) 
115     delete theData;                            << 134     G4cout << "G4VEmProcess destruct " << GetProcessName() 
116     delete theEnergyOfCrossSectionMax;         << 135      << G4endl;
                                                   >> 136   Clear();
                                                   >> 137   if(theLambdaTable) {
                                                   >> 138     theLambdaTable->clearAndDestroy();
                                                   >> 139     delete theLambdaTable;
117   }                                               140   }
118   delete modelManager;                            141   delete modelManager;
119   delete biasManager;                          << 142   (G4LossTableManager::Instance())->DeRegister(this);
120   lManager->DeRegister(this);                  << 
121 }                                                 143 }
122                                                   144 
123 //....oooOO0OOooo........oooOO0OOooo........oo    145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
124                                                   146 
125 void G4VEmProcess::AddEmModel(G4int order, G4V << 147 void G4VEmProcess::Clear()
126                               const G4Region*  << 
127 {                                                 148 {
128   if(nullptr == ptr) { return; }               << 149   delete [] theEnergyOfCrossSectionMax;
129   G4VEmFluctuationModel* fm = nullptr;         << 150   delete [] theCrossSectionMax;
130   modelManager->AddEmModel(order, ptr, fm, reg << 151   delete [] idxDERegions;
131   ptr->SetParticleChange(pParticleChange);     << 152   theEnergyOfCrossSectionMax = 0;
                                                   >> 153   theCrossSectionMax = 0;
                                                   >> 154   idxDERegions = 0;
                                                   >> 155   currentCouple = 0;
                                                   >> 156   preStepLambda = 0.0;
                                                   >> 157   mfpKinEnergy  = DBL_MAX;
                                                   >> 158   deRegions.clear();
                                                   >> 159   nDERegions = 0;
132 }                                                 160 }
133                                                   161 
134 //....oooOO0OOooo........oooOO0OOooo........oo    162 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
135                                                   163 
136 void G4VEmProcess::SetEmModel(G4VEmModel* ptr, << 164 void G4VEmProcess::AddEmModel(G4int order, G4VEmModel* p, 
                                                   >> 165             const G4Region* region)
137 {                                                 166 {
138   if(nullptr == ptr) { return; }               << 167   G4VEmFluctuationModel* fm = 0;
139   if(!emModels.empty()) {                      << 168   modelManager->AddEmModel(order, p, fm, region);
140     for(auto & em : emModels) { if(em == ptr)  << 169   if(p) { p->SetParticleChange(pParticleChange); }
141   }                                            << 
142   emModels.push_back(ptr);                     << 
143 }                                                 170 }
144                                                   171 
145 //....oooOO0OOooo........oooOO0OOooo........oo    172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
146                                                   173 
147 void G4VEmProcess::PreparePhysicsTable(const G << 174 void G4VEmProcess::SetModel(G4VEmModel* p, G4int index)
148 {                                                 175 {
149   if(nullptr == particle) { SetParticle(&part) << 176   G4int n = emModels.size();
                                                   >> 177   if(index >= n) { for(G4int i=n; i<=index; ++i) {emModels.push_back(0);} }
                                                   >> 178   emModels[index] = p;
                                                   >> 179 }
                                                   >> 180 
                                                   >> 181 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
150                                                   182 
151   if(part.GetParticleType() == "nucleus" &&    << 183 G4VEmModel* G4VEmProcess::Model(G4int index)
152      part.GetParticleSubType() == "generic") { << 184 {
                                                   >> 185   G4VEmModel* p = 0;
                                                   >> 186   if(index >= 0 && index <  G4int(emModels.size())) { p = emModels[index]; }
                                                   >> 187   return p;
                                                   >> 188 }
153                                                   189 
154     G4String pname = part.GetParticleName();   << 190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
155     if(pname != "deuteron" && pname != "triton << 
156        pname != "He3" && pname != "alpha" && p << 
157        pname != "helium" && pname != "hydrogen << 
158                                                   191 
159       particle = G4GenericIon::GenericIon();   << 192 void G4VEmProcess::UpdateEmModel(const G4String& nam, 
160       isIon = true;                            << 193          G4double emin, G4double emax)
161     }                                          << 194 {
162   }                                            << 195   modelManager->UpdateEmModel(nam, emin, emax);
163   if(particle != &part) { return; }            << 196 }
164                                                   197 
165   lManager->PreparePhysicsTable(&part, this);  << 198 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
166                                                   199 
167   // for new run                               << 200 G4VEmModel* G4VEmProcess::GetModelByIndex(G4int idx, G4bool ver)
168   currentCouple = nullptr;                     << 201 {
169   preStepLambda = 0.0;                         << 202   return modelManager->GetModel(idx, ver);
170   fLambdaEnergy = 0.0;                         << 203 }
171                                                   204 
172   InitialiseProcess(particle);                 << 205 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
173                                                   206 
174   G4LossTableBuilder* bld = lManager->GetTable << 207 void G4VEmProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
175   const G4ProductionCutsTable* theCoupleTable= << 208 {
176     G4ProductionCutsTable::GetProductionCutsTa << 209   if(!particle) { SetParticle(&part); }
177   theCutsGamma    = theCoupleTable->GetEnergyC << 210   if(1 < verboseLevel) {
178   theCutsElectron = theCoupleTable->GetEnergyC << 211     G4cout << "G4VEmProcess::PreparePhysicsTable() for "
179   theCutsPositron = theCoupleTable->GetEnergyC << 212            << GetProcessName()
180                                                << 213            << " and particle " << part.GetParticleName()
181   // initialisation of the process             << 214      << " local particle " << particle->GetParticleName() 
182   if(!actMinKinEnergy) { minKinEnergy = thePar << 215            << G4endl;
183   if(!actMaxKinEnergy) { maxKinEnergy = thePar << 216   }
184                                                << 
185   applyCuts       = theParameters->ApplyCuts() << 
186   lambdaFactor    = theParameters->LambdaFacto << 
187   invLambdaFactor = 1.0/lambdaFactor;          << 
188   theParameters->DefineRegParamForEM(this);    << 
189                                                << 
190   // integral option may be disabled           << 
191   if(!theParameters->Integral()) { fXSType = f << 
192                                                << 
193   // prepare tables                            << 
194   if(isTheMaster) {                            << 
195     if(nullptr == theData) { theData = new G4E << 
196                                                   217 
197     if(buildLambdaTable) {                     << 218   (G4LossTableManager::Instance())->PreparePhysicsTable(&part, this);
198       theLambdaTable = theData->MakeTable(0);  << 219 
199       bld->InitialiseBaseMaterials(theLambdaTa << 220   if(particle == &part) {
                                                   >> 221     Clear();
                                                   >> 222     InitialiseProcess(particle);
                                                   >> 223 
                                                   >> 224     // initialisation of models
                                                   >> 225     G4int nmod = modelManager->NumberOfModels();
                                                   >> 226     for(G4int i=0; i<nmod; ++i) {
                                                   >> 227       G4VEmModel* mod = modelManager->GetModel(i);
                                                   >> 228       mod->SetPolarAngleLimit(polarAngleLimit);
                                                   >> 229       if(mod->HighEnergyLimit() > maxKinEnergy) {
                                                   >> 230   mod->SetHighEnergyLimit(maxKinEnergy);
                                                   >> 231       }
                                                   >> 232     }
                                                   >> 233 
                                                   >> 234     theCuts = modelManager->Initialise(particle,secondaryParticle,
                                                   >> 235                2.,verboseLevel);
                                                   >> 236     const G4ProductionCutsTable* theCoupleTable=
                                                   >> 237           G4ProductionCutsTable::GetProductionCutsTable();
                                                   >> 238     theCutsGamma    = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut);
                                                   >> 239     theCutsElectron = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut);
                                                   >> 240     theCutsPositron = theCoupleTable->GetEnergyCutsVector(idxG4PositronCut);
                                                   >> 241 
                                                   >> 242     // prepare tables
                                                   >> 243     if(buildLambdaTable){
                                                   >> 244       theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable);
                                                   >> 245     }
                                                   >> 246   }
                                                   >> 247   // Deexcitation
                                                   >> 248   if (nDERegions>0) {
                                                   >> 249 
                                                   >> 250     const G4ProductionCutsTable* theCoupleTable=
                                                   >> 251           G4ProductionCutsTable::GetProductionCutsTable();
                                                   >> 252     size_t numOfCouples = theCoupleTable->GetTableSize();
                                                   >> 253 
                                                   >> 254     idxDERegions = new G4bool[numOfCouples];
                                                   >> 255 
                                                   >> 256     for (size_t j=0; j<numOfCouples; ++j) {
                                                   >> 257 
                                                   >> 258       const G4MaterialCutsCouple* couple =
                                                   >> 259         theCoupleTable->GetMaterialCutsCouple(j);
                                                   >> 260       const G4ProductionCuts* pcuts = couple->GetProductionCuts();
                                                   >> 261       G4bool reg = false;
                                                   >> 262       for(G4int i=0; i<nDERegions; ++i) {
                                                   >> 263   if(deRegions[i]) {
                                                   >> 264     if(pcuts == deRegions[i]->GetProductionCuts()) reg = true;
                                                   >> 265   }
                                                   >> 266       }
                                                   >> 267       idxDERegions[j] = reg;
200     }                                             268     }
201     // high energy table                       << 269   }
202     if(minKinEnergyPrim < maxKinEnergy) {      << 270   if (1 < verboseLevel && nDERegions>0) {
203       theLambdaTablePrim = theData->MakeTable( << 271     G4cout << " Deexcitation is activated for regions: " << G4endl;
204       bld->InitialiseBaseMaterials(theLambdaTa << 272     for (G4int i=0; i<nDERegions; ++i) {
205     }                                          << 273       const G4Region* r = deRegions[i];
206   }                                            << 274       G4cout << "           " << r->GetName() << G4endl;
207   // models                                    << 275     }
208   baseMat = bld->GetBaseMaterialFlag();        << 276   }
209   numberOfModels = modelManager->NumberOfModel << 
210   currentModel = modelManager->GetModel(0);    << 
211   if(nullptr != lManager->AtomDeexcitation())  << 
212     modelManager->SetFluoFlag(true);           << 
213   }                                            << 
214   // forced biasing                            << 
215   if(nullptr != biasManager) {                 << 
216     biasManager->Initialise(part, GetProcessNa << 
217     biasFlag = false;                          << 
218   }                                            << 
219                                                << 
220   theCuts =                                    << 
221     G4EmTableUtil::PrepareEmProcess(this, part << 
222                                     modelManag << 
223                                     secID, tri << 
224                                     verboseLev << 
225 }                                                 277 }
226                                                   278 
227 //....oooOO0OOooo........oooOO0OOooo........oo    279 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
228                                                   280 
229 void G4VEmProcess::BuildPhysicsTable(const G4P    281 void G4VEmProcess::BuildPhysicsTable(const G4ParticleDefinition& part)
230 {                                                 282 {
231   if(nullptr == masterProc) {                  << 283   G4String partname = part.GetParticleName();
232     if(isTheMaster) { masterProc = this; }     << 284   if(1 < verboseLevel) {
233     else { masterProc = static_cast<const G4VE << 285     G4cout << "G4VEmProcess::BuildPhysicsTable() for "
234   }                                            << 286            << GetProcessName()
235   G4int nModels = modelManager->NumberOfModels << 287            << " and particle " << partname
236   G4bool isLocked = theParameters->IsPrintLock << 288      << " buildLambdaTable= " << buildLambdaTable
237   G4bool toBuild = (buildLambdaTable || minKin << 289            << G4endl;
238                                                << 290   }
239   G4EmTableUtil::BuildEmProcess(this, masterPr << 291 
240                                 nModels, verbo << 292   (G4LossTableManager::Instance())->BuildPhysicsTable(particle);
241                                 isLocked, toBu << 293   if(buildLambdaTable) {
                                                   >> 294     BuildLambdaTable();
                                                   >> 295     FindLambdaMax();
                                                   >> 296   }
                                                   >> 297 
                                                   >> 298   // reduce printout for nuclear stopping
                                                   >> 299   G4bool gproc = true;
                                                   >> 300   G4int st = GetProcessSubType();
                                                   >> 301   if(st == fCoulombScattering && part.GetParticleType() == "nucleus" && 
                                                   >> 302      partname != "GenericIon" && partname != "alpha") { gproc = false; } 
                                                   >> 303 
                                                   >> 304   if(gproc && 0 < verboseLevel) { PrintInfoDefinition(); }
                                                   >> 305 
                                                   >> 306   if(1 < verboseLevel) {
                                                   >> 307     G4cout << "G4VEmProcess::BuildPhysicsTable() done for "
                                                   >> 308            << GetProcessName()
                                                   >> 309            << " and particle " << partname
                                                   >> 310            << G4endl;
                                                   >> 311   }
242 }                                                 312 }
243                                                   313 
244 //....oooOO0OOooo........oooOO0OOooo........oo    314 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
245                                                   315 
246 void G4VEmProcess::BuildLambdaTable()             316 void G4VEmProcess::BuildLambdaTable()
247 {                                                 317 {
248   G4double scale = theParameters->MaxKinEnergy << 318   if(1 < verboseLevel) {
249   G4int nbin =                                 << 319     G4cout << "G4EmProcess::BuildLambdaTable() for process "
250     theParameters->NumberOfBinsPerDecade()*G4l << 320            << GetProcessName() << " and particle "
251   if(actBinning) { nbin = std::max(nbin, nLamb << 321            << particle->GetParticleName()
252   scale = nbin/G4Log(scale);                   << 322            << G4endl;
253                                                << 323   }
254   G4LossTableBuilder* bld = lManager->GetTable << 
255   G4EmTableUtil::BuildLambdaTable(this, partic << 
256                                   bld, theLamb << 
257                                   minKinEnergy << 
258                                   maxKinEnergy << 
259                                   startFromNul << 
260 }                                              << 
261                                                   324 
262 //....oooOO0OOooo........oooOO0OOooo........oo << 325   // Access to materials
                                                   >> 326   const G4ProductionCutsTable* theCoupleTable=
                                                   >> 327         G4ProductionCutsTable::GetProductionCutsTable();
                                                   >> 328   size_t numOfCouples = theCoupleTable->GetTableSize();
263                                                   329 
264 void G4VEmProcess::StreamInfo(std::ostream& ou << 330   G4bool splineFlag = (G4LossTableManager::Instance())->SplineFlag();
265                   const G4ParticleDefinition&  << 331 
266 {                                              << 332   G4PhysicsLogVector* aVector = 0;
267   G4String indent = (rst ? "  " : "");         << 333   G4PhysicsLogVector* bVector = 0;
268   out << std::setprecision(6);                 << 334 
269   out << G4endl << indent << GetProcessName()  << 335   for(size_t i=0; i<numOfCouples; ++i) {
270   if (!rst) {                                  << 336 
271     out << " for " << part.GetParticleName();  << 337     if (theLambdaTable->GetFlag(i)) {
272   }                                            << 338 
273   if(fXSType != fEmNoIntegral)  { out << " XSt << 339       // create physics vector and fill it
274   if(applyCuts) { out << " applyCuts:1 "; }    << 340       const G4MaterialCutsCouple* couple = 
275   G4int subtype = GetProcessSubType();         << 341   theCoupleTable->GetMaterialCutsCouple(i);
276   out << " SubType=" << subtype;               << 342       if(!bVector) {
277   if (subtype == fAnnihilation) {              << 343   aVector = 
278     G4int mod = theParameters->PositronAtRestM << 344     static_cast<G4PhysicsLogVector*>(LambdaPhysicsVector(couple));
279     const G4String namp[2] = {"Simple", "Allis << 345         bVector = aVector;
280     out << " AtRestModel:" << namp[mod];       << 346       } else {
281   }                                            << 347         aVector = new G4PhysicsLogVector(*bVector);
282   if(biasFactor != 1.0) { out << "  BiasingFac << 
283   out << " BuildTable=" << buildLambdaTable << << 
284   if(buildLambdaTable) {                       << 
285     if(particle == &part) {                    << 
286       for(auto & v : *theLambdaTable) {        << 
287         if(nullptr != v) {                     << 
288           out << "      Lambda table from ";   << 
289           G4double emin = v->Energy(0);        << 
290           G4double emax = v->GetMaxEnergy();   << 
291           G4int nbin = G4int(v->GetVectorLengt << 
292           if(emin > minKinEnergy) { out << "th << 
293           else { out << G4BestUnit(emin,"Energ << 
294           out << " to "                        << 
295               << G4BestUnit(emax,"Energy")     << 
296               << ", " << G4lrint(nbin/std::log << 
297               << " bins/decade, spline: "      << 
298               << splineFlag << G4endl;         << 
299           break;                               << 
300         }                                      << 
301       }                                        << 
302     } else {                                   << 
303       out << "      Used Lambda table of "     << 
304       << particle->GetParticleName() << G4endl << 
305     }                                          << 
306   }                                            << 
307   if(minKinEnergyPrim < maxKinEnergy) {        << 
308     if(particle == &part) {                    << 
309       for(auto & v : *theLambdaTablePrim) {    << 
310         if(nullptr != v) {                     << 
311           out << "      LambdaPrime table from << 
312               << G4BestUnit(v->Energy(0),"Ener << 
313               << " to "                        << 
314               << G4BestUnit(v->GetMaxEnergy(), << 
315               << " in " << v->GetVectorLength( << 
316               << " bins " << G4endl;           << 
317           break;                               << 
318         }                                      << 
319       }                                           348       }
320     } else {                                   << 349   //      G4PhysicsVector* aVector = LambdaPhysicsVector(couple);
321       out << "      Used LambdaPrime table of  << 350       aVector->SetSpline(splineFlag);
322                << particle->GetParticleName()  << 351       modelManager->FillLambdaVector(aVector, couple, startFromNull);
                                                   >> 352       if(splineFlag) aVector->FillSecondDerivatives();
                                                   >> 353       G4PhysicsTableHelper::SetPhysicsVector(theLambdaTable, i, aVector);
323     }                                             354     }
324   }                                               355   }
325   StreamProcessInfo(out);                      << 
326   modelManager->DumpModelList(out, verboseLeve << 
327                                                   356 
328   if(verboseLevel > 2 && buildLambdaTable) {   << 357   if(1 < verboseLevel) {
329     out << "      LambdaTable address= " << th << 358     G4cout << "Lambda table is built for "
330     if(theLambdaTable && particle == &part) {  << 359            << particle->GetParticleName()
331       out << (*theLambdaTable) << G4endl;      << 360            << G4endl;
332     }                                          << 
333   }                                               361   }
334 }                                                 362 }
335                                                   363 
336 //....oooOO0OOooo........oooOO0OOooo........oo    364 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
337                                                   365 
338 void G4VEmProcess::StartTracking(G4Track* trac << 366 void G4VEmProcess::PrintInfoDefinition()
339 {                                                 367 {
340   // reset parameters for the new track        << 368   if(verboseLevel > 0) {
341   currentParticle = track->GetParticleDefiniti << 369     G4cout << G4endl << GetProcessName() << ":   for  "
342   theNumberOfInteractionLengthLeft = -1.0;     << 370            << particle->GetParticleName();
343   mfpKinEnergy = DBL_MAX;                      << 371     if(integral) G4cout << ", integral: 1 ";
344   preStepLambda = 0.0;                         << 372     if(applyCuts) G4cout << ", applyCuts: 1 ";
345                                                << 373     G4cout << "    SubType= " << GetProcessSubType() << G4endl;
346   if(isIon) { massRatio = proton_mass_c2/curre << 374     if(buildLambdaTable) {
347                                                << 375       G4cout << "      Lambda tables from "
348   // forced biasing only for primary particles << 376        << G4BestUnit(minKinEnergy,"Energy") 
349   if(biasManager) {                            << 377        << " to "
350     if(0 == track->GetParentID()) {            << 378        << G4BestUnit(maxKinEnergy,"Energy")
351       // primary particle                      << 379        << " in " << nLambdaBins << " bins, spline: " 
352       biasFlag = true;                         << 380        << (G4LossTableManager::Instance())->SplineFlag()
353       biasManager->ResetForcedInteraction();   << 381        << G4endl;
354     }                                             382     }
                                                   >> 383     PrintInfo();
                                                   >> 384     modelManager->DumpModelList(verboseLevel);
                                                   >> 385   }
                                                   >> 386 
                                                   >> 387   if(verboseLevel > 2 && buildLambdaTable) {
                                                   >> 388     G4cout << "      LambdaTable address= " << theLambdaTable << G4endl;
                                                   >> 389     if(theLambdaTable) { G4cout << (*theLambdaTable) << G4endl; }
355   }                                               390   }
356 }                                                 391 }
357                                                   392 
358 //....oooOO0OOooo........oooOO0OOooo........oo    393 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
359                                                   394 
360 G4double G4VEmProcess::PostStepGetPhysicalInte    395 G4double G4VEmProcess::PostStepGetPhysicalInteractionLength(
361                              const G4Track& tr    396                              const G4Track& track,
362                              G4double   previo    397                              G4double   previousStepSize,
363                              G4ForceCondition*    398                              G4ForceCondition* condition)
364 {                                                 399 {
                                                   >> 400   // condition is set to "Not Forced"
365   *condition = NotForced;                         401   *condition = NotForced;
366   G4double x = DBL_MAX;                           402   G4double x = DBL_MAX;
367                                                << 403   if(previousStepSize <= DBL_MIN) { theNumberOfInteractionLengthLeft = -1.0; }
368   DefineMaterial(track.GetMaterialCutsCouple() << 404   InitialiseStep(track);
369   preStepKinEnergy = track.GetKineticEnergy(); << 405   if(!currentModel->IsActive(preStepKinEnergy)) { return x; }
370   const G4double scaledEnergy = preStepKinEner << 406 
371   SelectModel(scaledEnergy, currentCoupleIndex << 407   if(preStepKinEnergy < mfpKinEnergy) {
372   /*                                           << 408     if (integral) ComputeIntegralLambda(preStepKinEnergy);
373   G4cout << "PostStepGetPhysicalInteractionLen << 409     else  preStepLambda = GetCurrentLambda(preStepKinEnergy);
374          << "  couple: " << currentCouple << G << 410     if(preStepLambda <= DBL_MIN) mfpKinEnergy = 0.0;
375   */                                           << 
376   if(!currentModel->IsActive(scaledEnergy)) {  << 
377     theNumberOfInteractionLengthLeft = -1.0;   << 
378     currentInteractionLength = DBL_MAX;        << 
379     mfpKinEnergy = DBL_MAX;                    << 
380     preStepLambda = 0.0;                       << 
381     return x;                                  << 
382   }                                            << 
383                                                << 
384   // forced biasing only for primary particles << 
385   if(biasManager) {                            << 
386     if(0 == track.GetParentID()) {             << 
387       if(biasFlag &&                           << 
388          biasManager->ForcedInteractionRegion( << 
389         return biasManager->GetStepLimit((G4in << 
390       }                                        << 
391     }                                          << 
392   }                                               411   }
393                                                   412 
394   // compute mean free path                    << 413   // non-zero cross section
395                                                << 414   if(preStepLambda > DBL_MIN) { 
396   ComputeIntegralLambda(preStepKinEnergy, trac << 
397                                                << 
398   // zero cross section                        << 
399   if(preStepLambda <= 0.0) {                   << 
400     theNumberOfInteractionLengthLeft = -1.0;   << 
401     currentInteractionLength = DBL_MAX;        << 
402                                                << 
403   } else {                                     << 
404                                                << 
405     // non-zero cross section                  << 
406     if (theNumberOfInteractionLengthLeft < 0.0    415     if (theNumberOfInteractionLengthLeft < 0.0) {
407                                                << 
408       // beggining of tracking (or just after     416       // beggining of tracking (or just after DoIt of this process)
409       theNumberOfInteractionLengthLeft = -G4Lo << 417       ResetNumberOfInteractionLengthLeft();
410       theInitialNumberOfInteractionLength = th << 418     } else if(currentInteractionLength < DBL_MAX) {
411                                                << 419       // subtract NumberOfInteractionLengthLeft
412     } else {                                   << 420       SubtractNumberOfInteractionLengthLeft(previousStepSize);
413                                                << 421       if(theNumberOfInteractionLengthLeft < 0.)
414       theNumberOfInteractionLengthLeft -=      << 422   theNumberOfInteractionLengthLeft = perMillion;
415         previousStepSize/currentInteractionLen << 
416       theNumberOfInteractionLengthLeft =       << 
417         std::max(theNumberOfInteractionLengthL << 
418     }                                             423     }
419                                                   424 
420     // new mean free path and step limit for t << 425     // get mean free path and step limit
421     currentInteractionLength = 1.0/preStepLamb    426     currentInteractionLength = 1.0/preStepLambda;
422     x = theNumberOfInteractionLengthLeft * cur    427     x = theNumberOfInteractionLengthLeft * currentInteractionLength;
423   }                                            << 428 #ifdef G4VERBOSE
424   return x;                                    << 429     if (verboseLevel>2){
425 }                                              << 430       G4cout << "G4VEmProcess::PostStepGetPhysicalInteractionLength ";
426                                                << 431       G4cout << "[ " << GetProcessName() << "]" << G4endl; 
427 //....oooOO0OOooo........oooOO0OOooo........oo << 432       G4cout << " for " << currentParticle->GetParticleName() 
                                                   >> 433        << " in Material  " <<  currentMaterial->GetName()
                                                   >> 434        << " Ekin(MeV)= " << preStepKinEnergy/MeV 
                                                   >> 435        <<G4endl;
                                                   >> 436       G4cout << "MeanFreePath = " << currentInteractionLength/cm << "[cm]" 
                                                   >> 437        << "InteractionLength= " << x/cm <<"[cm] " <<G4endl;
                                                   >> 438     }
                                                   >> 439 #endif
428                                                   440 
429 void G4VEmProcess::ComputeIntegralLambda(G4dou << 441     // zero cross section case
430 {                                              << 442   } else {
431   if (fXSType == fEmNoIntegral) {              << 443     if(theNumberOfInteractionLengthLeft > DBL_MIN && 
432     preStepLambda = GetCurrentLambda(e, LogEki << 444        currentInteractionLength < DBL_MAX) {
433                                                   445 
434   } else if (fXSType == fEmIncreasing) {       << 446       // subtract NumberOfInteractionLengthLeft
435     if(e*invLambdaFactor < mfpKinEnergy) {     << 447       SubtractNumberOfInteractionLengthLeft(previousStepSize);
436       preStepLambda = GetCurrentLambda(e, LogE << 448       if(theNumberOfInteractionLengthLeft < 0.)
437       mfpKinEnergy = (preStepLambda > 0.0) ? e << 449   theNumberOfInteractionLengthLeft = perMillion;
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     }                                             450     }
459   } else {                                     << 451     currentInteractionLength = DBL_MAX;
460     preStepLambda = GetCurrentLambda(e, LogEki << 
461   }                                               452   }
                                                   >> 453   return x;
462 }                                                 454 }
463                                                   455 
464 //....oooOO0OOooo........oooOO0OOooo........oo    456 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
465                                                   457 
466 G4VParticleChange* G4VEmProcess::PostStepDoIt(    458 G4VParticleChange* G4VEmProcess::PostStepDoIt(const G4Track& track,
467                                                << 459                                               const G4Step&)
468 {                                                 460 {
469   // clear number of interaction lengths in an << 
470   theNumberOfInteractionLengthLeft = -1.0;     << 
471   mfpKinEnergy = DBL_MAX;                      << 
472                                                << 
473   fParticleChange.InitializeForPostStep(track)    461   fParticleChange.InitializeForPostStep(track);
474                                                   462 
475   // Do not make anything if particle is stopp    463   // Do not make anything if particle is stopped, the annihilation then
476   // should be performed by the AtRestDoIt!       464   // should be performed by the AtRestDoIt!
477   if (track.GetTrackStatus() == fStopButAlive)    465   if (track.GetTrackStatus() == fStopButAlive) { return &fParticleChange; }
478                                                   466 
479   const G4double finalT = track.GetKineticEner << 467   G4double finalT = track.GetKineticEnergy();
480                                                << 
481   // forced process - should happen only once  << 
482   if(biasFlag) {                               << 
483     if(biasManager->ForcedInteractionRegion((G << 
484       biasFlag = false;                        << 
485     }                                          << 
486   }                                            << 
487                                                << 
488   // check active and select model             << 
489   const G4double scaledEnergy = finalT*massRat << 
490   SelectModel(scaledEnergy, currentCoupleIndex << 
491   if(!currentModel->IsActive(scaledEnergy)) {  << 
492                                                   468 
493   // Integral approach                            469   // Integral approach
494   if (fXSType != fEmNoIntegral) {              << 470   if (integral) {
495     const G4double logFinalT =                 << 471     G4double lx = GetLambda(finalT, currentCouple);
496       track.GetDynamicParticle()->GetLogKineti << 472     if(preStepLambda<lx && 1 < verboseLevel) {
497     const G4double lx = std::max(GetCurrentLam << 473       G4cout << "WARING: for " << currentParticle->GetParticleName() 
498 #ifdef G4VERBOSE                               << 474              << " and " << GetProcessName()
499     if(preStepLambda < lx && 1 < verboseLevel) << 475              << " E(MeV)= " << finalT/MeV
500       G4cout << "WARNING: for " << currentPart << 476              << " preLambda= " << preStepLambda << " < " << lx << " (postLambda) "
501              << " and " << GetProcessName() << << 477        << G4endl;  
502              << " preLambda= " << preStepLambd << 
503              << " < " << lx << " (postLambda)  << 
504     }                                             478     }
505 #endif                                         << 479 
506     // if false interaction then use new cross << 480     if(preStepLambda*G4UniformRand() > lx) {
507     // if both values are zero - no interactio << 481       ClearNumberOfInteractionLengthLeft();
508     if(preStepLambda*G4UniformRand() >= lx) {  << 
509       return &fParticleChange;                    482       return &fParticleChange;
510     }                                             483     }
511   }                                               484   }
512                                                   485 
513   // define new weight for primary and seconda << 486   SelectModel(finalT, currentCoupleIndex);
514   G4double weight = fParticleChange.GetParentW << 487   if(!currentModel->IsActive(finalT)) { return &fParticleChange; }
515   if(weightFlag) {                             << 488   if(useDeexcitation) {
516     weight /= biasFactor;                      << 489     currentModel->SetDeexcitationFlag(idxDERegions[currentCoupleIndex]);
517     fParticleChange.ProposeWeight(weight);     << 
518   }                                               490   }
519                                                << 491   /*  
520 #ifdef G4VERBOSE                               << 492   if(0 < verboseLevel) {
521   if(1 < verboseLevel) {                       << 
522     G4cout << "G4VEmProcess::PostStepDoIt: Sam    493     G4cout << "G4VEmProcess::PostStepDoIt: Sample secondary; E= "
523            << finalT/MeV                          494            << finalT/MeV
524            << " MeV; model= (" << currentModel    495            << " MeV; model= (" << currentModel->LowEnergyLimit()
525            << ", " <<  currentModel->HighEnerg    496            << ", " <<  currentModel->HighEnergyLimit() << ")"
526            << G4endl;                             497            << G4endl;
527   }                                               498   }
528 #endif                                         << 499   */
529                                                   500 
                                                   >> 501   
530   // sample secondaries                           502   // sample secondaries
531   secParticles.clear();                           503   secParticles.clear();
532   currentModel->SampleSecondaries(&secParticle    504   currentModel->SampleSecondaries(&secParticles, 
533                                   currentCoupl << 505           currentCouple, 
534                                   track.GetDyn << 506           track.GetDynamicParticle(),
535                                   (*theCuts)[c << 507           (*theCuts)[currentCoupleIndex]);
536                                                << 
537   G4int num0 = (G4int)secParticles.size();     << 
538                                                << 
539   // splitting or Russian roulette             << 
540   if(biasManager) {                            << 
541     if(biasManager->SecondaryBiasingRegion((G4 << 
542       G4double eloss = 0.0;                    << 
543       weight *= biasManager->ApplySecondaryBia << 
544         secParticles, track, currentModel, &fP << 
545         (G4int)currentCoupleIndex, (*theCuts)[ << 
546         step.GetPostStepPoint()->GetSafety()); << 
547       if(eloss > 0.0) {                        << 
548         eloss += fParticleChange.GetLocalEnerg << 
549         fParticleChange.ProposeLocalEnergyDepo << 
550       }                                        << 
551     }                                          << 
552   }                                            << 
553                                                   508 
554   // save secondaries                             509   // save secondaries
555   G4int num = (G4int)secParticles.size();      << 510   G4int num = secParticles.size();
556   if(num > 0) {                                   511   if(num > 0) {
557                                                   512 
558     fParticleChange.SetNumberOfSecondaries(num    513     fParticleChange.SetNumberOfSecondaries(num);
559     G4double edep = fParticleChange.GetLocalEn    514     G4double edep = fParticleChange.GetLocalEnergyDeposit();
560     G4double time = track.GetGlobalTime();     << 
561                                                << 
562     G4int n1(0), n2(0);                        << 
563     if(num0 > mainSecondaries) {               << 
564       currentModel->FillNumberOfSecondaries(n1 << 
565     }                                          << 
566                                                   515      
567     for (G4int i=0; i<num; ++i) {                 516     for (G4int i=0; i<num; ++i) {
568       G4DynamicParticle* dp = secParticles[i];    517       G4DynamicParticle* dp = secParticles[i];
569       if (nullptr != dp) {                     << 518       const G4ParticleDefinition* p = dp->GetParticleDefinition();
570         const G4ParticleDefinition* p = dp->Ge << 519       G4double e = dp->GetKineticEnergy();
571         G4double e = dp->GetKineticEnergy();   << 520       G4bool good = true;
572         G4bool good = true;                    << 521       if(applyCuts) {
573         if(applyCuts) {                        << 522   if (p == theGamma) {
574           if (p == theGamma) {                 << 523     if (e < (*theCutsGamma)[currentCoupleIndex]) good = false;
575             if (e < (*theCutsGamma)[currentCou << 524 
576                                                << 525   } else if (p == theElectron) {
577           } else if (p == theElectron) {       << 526     if (e < (*theCutsElectron)[currentCoupleIndex]) good = false;
578             if (e < (*theCutsElectron)[current << 527 
579                                                << 528   } else if (p == thePositron) {
580           } else if (p == thePositron) {       << 529     if (electron_mass_c2 < (*theCutsGamma)[currentCoupleIndex] &&
581             if (electron_mass_c2 < (*theCutsGa << 530         e < (*theCutsPositron)[currentCoupleIndex]) {
582                 e < (*theCutsPositron)[current << 531       good = false;
583               good = false;                    << 532       e += 2.0*electron_mass_c2;
584               e += 2.0*electron_mass_c2;       << 533     }
585             }                                  << 534   }
586           }                                    << 535         if(!good) {
587           // added secondary if it is good     << 536     delete dp;
588         }                                      << 537     edep += e;
589         if (good) {                            << 538   }
590           G4Track* t = new G4Track(dp, time, t << 539       }
591           t->SetTouchableHandle(track.GetTouch << 540       if (good) fParticleChange.AddSecondary(dp);
592           if (biasManager) {                   << 541     } 
593             t->SetWeight(weight * biasManager- << 
594           } else {                             << 
595             t->SetWeight(weight);              << 
596           }                                    << 
597           pParticleChange->AddSecondary(t);    << 
598                                                << 
599           // define type of secondary          << 
600           if(i < mainSecondaries) {            << 
601             t->SetCreatorModelID(secID);       << 
602             if(GetProcessSubType() == fCompton << 
603               t->SetCreatorModelID(_ComptonGam << 
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 {                           << 
617               t->SetCreatorModelID(biasID);    << 
618             }                                  << 
619           }                                    << 
620           /*                                   << 
621           G4cout << "Secondary(post step) has  << 
622                  << ", Ekin= " << t->GetKineti << 
623                  << GetProcessName() << " fluo << 
624                  << " augerID= " << augerID << << 
625           */                                   << 
626         } else {                               << 
627           delete dp;                           << 
628           edep += e;                           << 
629         }                                      << 
630       }                                        << 
631     }                                          << 
632     fParticleChange.ProposeLocalEnergyDeposit(    542     fParticleChange.ProposeLocalEnergyDeposit(edep);
633   }                                               543   }
634                                                   544 
635   if(0.0 == fParticleChange.GetProposedKinetic << 545   ClearNumberOfInteractionLengthLeft();
636      fAlive == fParticleChange.GetTrackStatus( << 
637     if(particle->GetProcessManager()->GetAtRes << 
638          { fParticleChange.ProposeTrackStatus( << 
639     else { fParticleChange.ProposeTrackStatus( << 
640   }                                            << 
641                                                << 
642   return &fParticleChange;                        546   return &fParticleChange;
643 }                                                 547 }
644                                                   548 
645 //....oooOO0OOooo........oooOO0OOooo........oo    549 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
646                                                   550 
647 G4bool G4VEmProcess::StorePhysicsTable(const G    551 G4bool G4VEmProcess::StorePhysicsTable(const G4ParticleDefinition* part,
648                                        const G << 552                const G4String& directory,
649                                        G4bool  << 553                      G4bool ascii)
650 {                                                 554 {
651   if(!isTheMaster || part != particle) { retur << 555   G4bool yes = true;
652   if(G4EmTableUtil::StoreTable(this, part, the << 556 
653              directory, "Lambda",              << 557   if ( theLambdaTable && part == particle) {
654                                verboseLevel, a << 558     const G4String name = 
655      G4EmTableUtil::StoreTable(this, part, the << 559       GetPhysicsTableFileName(part,directory,"Lambda",ascii);
656              directory, "LambdaPrim",          << 560     yes = theLambdaTable->StorePhysicsTable(name,ascii);
657                                verboseLevel, a << 561 
658      return true;                              << 562     if ( yes ) {
                                                   >> 563       G4cout << "Physics tables are stored for " << particle->GetParticleName()
                                                   >> 564              << " and process " << GetProcessName()
                                                   >> 565        << " in the directory <" << directory
                                                   >> 566        << "> " << G4endl;
                                                   >> 567     } else {
                                                   >> 568       G4cout << "Fail to store Physics Tables for " 
                                                   >> 569        << particle->GetParticleName()
                                                   >> 570              << " and process " << GetProcessName()
                                                   >> 571        << " in the directory <" << directory
                                                   >> 572        << "> " << G4endl;
                                                   >> 573     }
659   }                                               574   }
660   return false;                                << 575   return yes;
661 }                                                 576 }
662                                                   577 
663 //....oooOO0OOooo........oooOO0OOooo........oo    578 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
664                                                   579 
665 G4bool G4VEmProcess::RetrievePhysicsTable(cons    580 G4bool G4VEmProcess::RetrievePhysicsTable(const G4ParticleDefinition* part,
666                                           cons << 581                     const G4String& directory,
667                                           G4bo << 582             G4bool ascii)
668 {                                                 583 {
669   if(!isTheMaster || part != particle) { retur << 584   if(1 < verboseLevel) {
670   G4bool yes = true;                           << 585     G4cout << "G4VEmProcess::RetrievePhysicsTable() for "
671   if(buildLambdaTable) {                       << 586            << part->GetParticleName() << " and process "
672     yes = G4EmTableUtil::RetrieveTable(this, p << 587      << GetProcessName() << G4endl;
673                                        "Lambda << 
674                                        ascii,  << 
675   }                                            << 
676   if(yes && minKinEnergyPrim < maxKinEnergy) { << 
677     yes = G4EmTableUtil::RetrieveTable(this, p << 
678                                        "Lambda << 
679                                        ascii,  << 
680   }                                               588   }
681   return yes;                                  << 589   G4bool yes = true;
682 }                                              << 
683                                                   590 
684 //....oooOO0OOooo........oooOO0OOooo........oo << 591   if(!buildLambdaTable || particle != part) return yes;
685                                                   592 
686 G4double G4VEmProcess::GetCrossSection(G4doubl << 593   const G4String particleName = part->GetParticleName();
687                                        const G << 594   G4String filename;
688 {                                              << 
689   CurrentSetup(couple, kinEnergy);             << 
690   return GetCurrentLambda(kinEnergy, G4Log(kin << 
691 }                                              << 
692                                                   595 
693 //....oooOO0OOooo........oooOO0OOooo........oo << 596   filename = GetPhysicsTableFileName(part,directory,"Lambda",ascii);
                                                   >> 597   yes = G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTable,
                                                   >> 598                filename,ascii);
                                                   >> 599   if ( yes ) {
                                                   >> 600     if (0 < verboseLevel) {
                                                   >> 601       G4cout << "Lambda table for " << particleName 
                                                   >> 602        << " is Retrieved from <"
                                                   >> 603              << filename << ">"
                                                   >> 604              << G4endl;
                                                   >> 605     }
                                                   >> 606     if((G4LossTableManager::Instance())->SplineFlag()) {
                                                   >> 607       size_t n = theLambdaTable->length();
                                                   >> 608       for(size_t i=0; i<n; ++i) {
                                                   >> 609         if((* theLambdaTable)[i]) {
                                                   >> 610     (* theLambdaTable)[i]->SetSpline(true);
                                                   >> 611   }
                                                   >> 612       }
                                                   >> 613     }
                                                   >> 614   } else {
                                                   >> 615     if (1 < verboseLevel) {
                                                   >> 616       G4cout << "Lambda table for " << particleName << " in file <"
                                                   >> 617              << filename << "> is not exist"
                                                   >> 618              << G4endl;
                                                   >> 619     }
                                                   >> 620   }
694                                                   621 
695 G4double G4VEmProcess::GetMeanFreePath(const G << 622   return yes;
696                                        G4doubl << 
697                                        G4Force << 
698 {                                              << 
699   *condition = NotForced;                      << 
700   return G4VEmProcess::MeanFreePath(track);    << 
701 }                                                 623 }
702                                                   624 
703 //....oooOO0OOooo........oooOO0OOooo........oo    625 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
704                                                   626 
705 G4double                                       << 627 void G4VEmProcess::ActivateDeexcitation(G4bool val, const G4Region* r)
706 G4VEmProcess::ComputeCrossSectionPerAtom(G4dou << 
707                                          G4dou << 
708 {                                                 628 {
709   SelectModel(kinEnergy, currentCoupleIndex);  << 629   G4RegionStore* regionStore = G4RegionStore::GetInstance();
710   return (currentModel) ?                      << 630   const G4Region* reg = r;
711     currentModel->ComputeCrossSectionPerAtom(c << 631   if (!reg) {reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);}
712                                              Z << 632 
                                                   >> 633   // the region is in the list
                                                   >> 634   if (nDERegions) {
                                                   >> 635     for (G4int i=0; i<nDERegions; ++i) {
                                                   >> 636       if (reg == deRegions[i]) {
                                                   >> 637   if(!val) deRegions[i] = 0;
                                                   >> 638         return;
                                                   >> 639       }
                                                   >> 640     }
                                                   >> 641   }
                                                   >> 642 
                                                   >> 643   // new region 
                                                   >> 644   if(val) {
                                                   >> 645     useDeexcitation = true;
                                                   >> 646     deRegions.push_back(reg);
                                                   >> 647     nDERegions++;
                                                   >> 648   } else {
                                                   >> 649     useDeexcitation = false;
                                                   >> 650   }
713 }                                                 651 }
714                                                   652 
715 //....oooOO0OOooo........oooOO0OOooo........oo    653 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
716                                                   654 
717 G4PhysicsVector*                               << 655 G4double G4VEmProcess::CrossSectionPerVolume(G4double kineticEnergy,
718 G4VEmProcess::LambdaPhysicsVector(const G4Mate << 656                const G4MaterialCutsCouple* couple)
719 {                                                 657 {
                                                   >> 658   // Cross section per atom is calculated
720   DefineMaterial(couple);                         659   DefineMaterial(couple);
721   G4PhysicsVector* newv = new G4PhysicsLogVect << 660   G4double cross = 0.0;
722                                                << 661   if(theLambdaTable) {
723   return newv;                                 << 662     cross = (((*theLambdaTable)[currentCoupleIndex])->Value(kineticEnergy));
724 }                                              << 663   } else {
725                                                << 664     SelectModel(kineticEnergy, currentCoupleIndex);
726 //....oooOO0OOooo........oooOO0OOooo........oo << 665     cross = currentModel->CrossSectionPerVolume(currentMaterial,
                                                   >> 666             currentParticle,kineticEnergy);
                                                   >> 667   }
727                                                   668 
728 const G4Element* G4VEmProcess::GetCurrentEleme << 669   if(cross < 0.0) { cross = 0.0; }
729 {                                              << 670   return cross;
730   return (nullptr != currentModel) ?           << 
731     currentModel->GetCurrentElement(currentMat << 
732 }                                                 671 }
733                                                   672 
734 //....oooOO0OOooo........oooOO0OOooo........oo    673 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
735                                                   674 
736 const G4Element* G4VEmProcess::GetTargetElemen << 675 G4double G4VEmProcess::GetMeanFreePath(const G4Track& track,
                                                   >> 676                G4double,
                                                   >> 677                G4ForceCondition* condition)
737 {                                                 678 {
738   return (nullptr != currentModel) ?           << 679   *condition = NotForced;
739     currentModel->GetCurrentElement(currentMat << 680   return G4VEmProcess::MeanFreePath(track);
740 }                                                 681 }
741                                                   682 
742 //....oooOO0OOooo........oooOO0OOooo........oo    683 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
743                                                   684 
744 const G4Isotope* G4VEmProcess::GetTargetIsotop << 685 G4double G4VEmProcess::MeanFreePath(const G4Track& track)
745 {                                                 686 {
746   return (nullptr != currentModel) ?           << 687   DefineMaterial(track.GetMaterialCutsCouple());
747     currentModel->GetCurrentIsotope(GetCurrent << 688   preStepLambda = GetCurrentLambda(track.GetKineticEnergy());
                                                   >> 689   G4double x = DBL_MAX;
                                                   >> 690   if(DBL_MIN < preStepLambda) x = 1.0/preStepLambda;
                                                   >> 691   return x;
748 }                                                 692 }
749                                                   693 
750 //....oooOO0OOooo........oooOO0OOooo........oo    694 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
751                                                   695 
752 void G4VEmProcess::SetCrossSectionBiasingFacto << 696 G4double 
                                                   >> 697 G4VEmProcess::ComputeCrossSectionPerAtom(G4double kineticEnergy, 
                                                   >> 698            G4double Z, G4double A, G4double cut)
753 {                                                 699 {
754   if(f > 0.0) {                                << 700   SelectModel(kineticEnergy, currentCoupleIndex);
755     biasFactor = f;                            << 701   G4double x = 0.0;
756     weightFlag = flag;                         << 702   if(currentModel) {
757     if(1 < verboseLevel) {                     << 703    x = currentModel->ComputeCrossSectionPerAtom(currentParticle,kineticEnergy,
758       G4cout << "### SetCrossSectionBiasingFac << 704              Z,A,cut);
759              << particle->GetParticleName()    << 
760              << " and process " << GetProcessN << 
761              << " biasFactor= " << f << " weig << 
762              << G4endl;                        << 
763     }                                          << 
764   }                                               705   }
                                                   >> 706   return x;
765 }                                                 707 }
766                                                   708 
767 //....oooOO0OOooo........oooOO0OOooo........oo    709 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
768                                                   710 
769 void                                           << 711 void G4VEmProcess::FindLambdaMax()
770 G4VEmProcess::ActivateForcedInteraction(G4doub << 
771                                         G4bool << 
772 {                                                 712 {
773   if(nullptr == biasManager) { biasManager = n << 
774   if(1 < verboseLevel) {                          713   if(1 < verboseLevel) {
775     G4cout << "### ActivateForcedInteraction:  << 714     G4cout << "### G4VEmProcess::FindLambdaMax: " 
776            << particle->GetParticleName()      << 715      << particle->GetParticleName() 
777            << " and process " << GetProcessNam << 716            << " and process " << GetProcessName() << G4endl; 
778            << " length(mm)= " << length/mm     << 717   }
779            << " in G4Region <" << r            << 718   size_t n = theLambdaTable->length();
780            << "> weightFlag= " << flag         << 719   G4PhysicsVector* pv = (*theLambdaTable)[0];
781            << G4endl;                          << 720   G4double e, s, emax, smax;
782   }                                            << 721   theEnergyOfCrossSectionMax = new G4double [n];
783   weightFlag = flag;                           << 722   theCrossSectionMax = new G4double [n];
784   biasManager->ActivateForcedInteraction(lengt << 723 
785 }                                              << 724   for (size_t i=0; i<n; ++i) {
786                                                << 725     pv = (*theLambdaTable)[i];
787 //....oooOO0OOooo........oooOO0OOooo........oo << 726     emax = DBL_MAX;
788                                                << 727     smax = 0.0;
789 void                                           << 728     if(pv) {
790 G4VEmProcess::ActivateSecondaryBiasing(const G << 729       size_t nb = pv->GetVectorLength();
791                  G4double factor,              << 730       emax = DBL_MAX;
792                  G4double energyLimit)         << 731       smax = 0.0;
793 {                                              << 732       if(nb > 0) {
794   if (0.0 <= factor) {                         << 733   for (size_t j=0; j<nb; ++j) {
795                                                << 734     e = pv->Energy(j);
796     // Range cut can be applied only for e-    << 735     s = (*pv)(j);
797     if(0.0 == factor && secondaryParticle != G << 736     if(s > smax) {
798       { return; }                              << 737       smax = s;
799                                                << 738       emax = e;
800     if(!biasManager) { biasManager = new G4EmB << 739     }
801     biasManager->ActivateSecondaryBiasing(regi << 740   }
                                                   >> 741       }
                                                   >> 742     }
                                                   >> 743     theEnergyOfCrossSectionMax[i] = emax;
                                                   >> 744     theCrossSectionMax[i] = smax;
802     if(1 < verboseLevel) {                        745     if(1 < verboseLevel) {
803       G4cout << "### ActivateSecondaryBiasing: << 746       G4cout << "For " << particle->GetParticleName() 
804        << " process " << GetProcessName()      << 747        << " Max CS at i= " << i << " emax(MeV)= " << emax/MeV
805        << " factor= " << factor                << 748        << " lambda= " << smax << G4endl;
806        << " in G4Region <" << region           << 
807        << "> energyLimit(MeV)= " << energyLimi << 
808        << G4endl;                              << 
809     }                                             749     }
810   }                                               750   }
811 }                                                 751 }
812                                                   752 
813 //....oooOO0OOooo........oooOO0OOooo........oo    753 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
814                                                   754 
815 void G4VEmProcess::SetLambdaBinning(G4int n)   << 755 G4PhysicsVector* G4VEmProcess::LambdaPhysicsVector(const G4MaterialCutsCouple*)
816 {                                              << 
817   if(5 < n && n < 10000000) {                  << 
818     nLambdaBins = n;                           << 
819     actBinning = true;                         << 
820   } else {                                     << 
821     G4double e = (G4double)n;                  << 
822     PrintWarning("SetLambdaBinning", e);       << 
823   }                                            << 
824 }                                              << 
825                                                << 
826 //....oooOO0OOooo........oooOO0OOooo........oo << 
827                                                << 
828 void G4VEmProcess::SetMinKinEnergy(G4double e) << 
829 {                                                 756 {
830   if(1.e-3*eV < e && e < maxKinEnergy) {       << 757   G4PhysicsVector* v = 
831     nLambdaBins = G4lrint(nLambdaBins*G4Log(ma << 758     new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nLambdaBins);
832                           /G4Log(maxKinEnergy/ << 759   v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
833     minKinEnergy = e;                          << 760   return v;
834     actMinKinEnergy = true;                    << 
835   } else { PrintWarning("SetMinKinEnergy", e); << 
836 }                                                 761 }
837                                                   762 
838 //....oooOO0OOooo........oooOO0OOooo........oo    763 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
839                                                   764 
840 void G4VEmProcess::SetMaxKinEnergy(G4double e) << 765 const G4Element* G4VEmProcess::GetCurrentElement() const
841 {                                              << 
842   if(minKinEnergy < e && e < 1.e+6*TeV) {      << 
843     nLambdaBins = G4lrint(nLambdaBins*G4Log(e/ << 
844                           /G4Log(maxKinEnergy/ << 
845     maxKinEnergy = e;                          << 
846     actMaxKinEnergy = true;                    << 
847   } else { PrintWarning("SetMaxKinEnergy", e); << 
848 }                                              << 
849                                                << 
850 //....oooOO0OOooo........oooOO0OOooo........oo << 
851                                                << 
852 void G4VEmProcess::SetMinKinEnergyPrim(G4doubl << 
853 {                                              << 
854   if(theParameters->MinKinEnergy() <= e &&     << 
855      e <= theParameters->MaxKinEnergy()) { min << 
856   else { PrintWarning("SetMinKinEnergyPrim", e << 
857 }                                              << 
858                                                << 
859 //....oooOO0OOooo........oooOO0OOooo........oo << 
860                                                << 
861 G4VEmProcess* G4VEmProcess::GetEmProcess(const << 
862 {                                              << 
863   return (nam == GetProcessName()) ? this : nu << 
864 }                                              << 
865                                                << 
866 //....oooOO0OOooo........oooOO0OOooo........oo << 
867                                                << 
868 G4double G4VEmProcess::PolarAngleLimit() const << 
869 {                                              << 
870   return theParameters->MscThetaLimit();       << 
871 }                                              << 
872                                                << 
873 //....oooOO0OOooo........oooOO0OOooo........oo << 
874                                                << 
875 void G4VEmProcess::PrintWarning(G4String tit,  << 
876 {                                              << 
877   G4String ss = "G4VEmProcess::" + tit;        << 
878   G4ExceptionDescription ed;                   << 
879   ed << "Parameter is out of range: " << val   << 
880      << " it will have no effect!\n" << "  Pro << 
881      << GetProcessName() << "  nbins= " << the << 
882      << " Emin(keV)= " << theParameters->MinKi << 
883      << " Emax(GeV)= " << theParameters->MaxKi << 
884   G4Exception(ss, "em0044", JustWarning, ed);  << 
885 }                                              << 
886                                                << 
887 //....oooOO0OOooo........oooOO0OOooo........oo << 
888                                                << 
889 void G4VEmProcess::ProcessDescription(std::ost << 
890 {                                                 766 {
891   if(nullptr != particle) {                    << 767   const G4Element* elm = 0;
892     StreamInfo(out, *particle, true);          << 768   if(currentModel) {elm = currentModel->GetCurrentElement(); }
893   }                                            << 769   return elm;
894 }                                                 770 }
895                                                   771 
896 //....oooOO0OOooo........oooOO0OOooo........oo    772 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
897                                                   773