Geant4 Cross Reference

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


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                    <<   3 // * DISCLAIMER                                                       *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th <<   5 // * The following disclaimer summarizes all the specific disclaimers *
  6 // * the Geant4 Collaboration.  It is provided <<   6 // * of contributors to this software. The specific disclaimers,which *
  7 // * conditions of the Geant4 Software License <<   7 // * govern, are listed with their locations in:                      *
  8 // * LICENSE and available at  http://cern.ch/ <<   8 // *   http://cern.ch/geant4/license                                  *
  9 // * include a list of copyright holders.      << 
 10 // *                                                9 // *                                                                  *
 11 // * Neither the authors of this software syst     10 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     11 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     12 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     13 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  <<  14 // * use.                                                             *
 16 // * for the full disclaimer and the limitatio << 
 17 // *                                               15 // *                                                                  *
 18 // * This  code  implementation is the result  <<  16 // * This  code  implementation is the  intellectual property  of the *
 19 // * technical work of the GEANT4 collaboratio <<  17 // * GEANT4 collaboration.                                            *
 20 // * By using,  copying,  modifying or  distri <<  18 // * By copying,  distributing  or modifying the Program (or any work *
 21 // * any work based  on the software)  you  ag <<  19 // * based  on  the Program)  you indicate  your  acceptance of  this *
 22 // * use  in  resulting  scientific  publicati <<  20 // * statement, and all its terms.                                    *
 23 // * acceptance of all terms of the Geant4 Sof << 
 24 // *******************************************     21 // ********************************************************************
 25 //                                                 22 //
                                                   >>  23 // $Id: G4EmCalculator.cc,v 1.24 2005/11/26 15:33:08 vnivanch Exp $
                                                   >>  24 // GEANT4 tag $Name: geant4-08-00-patch-01 $
 26 //                                                 25 //
 27 // -------------------------------------------     26 // -------------------------------------------------------------------
 28 //                                                 27 //
 29 // GEANT4 Class file                               28 // GEANT4 Class file
 30 //                                                 29 //
 31 //                                                 30 //
 32 // File name:     G4EmCalculator                   31 // File name:     G4EmCalculator
 33 //                                                 32 //
 34 // Author:        Vladimir Ivanchenko              33 // Author:        Vladimir Ivanchenko
 35 //                                                 34 //
 36 // Creation date: 28.06.2004                       35 // Creation date: 28.06.2004
 37 //                                                 36 //
                                                   >>  37 // Modifications:
                                                   >>  38 // 12.09.2004 Add verbosity (V.Ivanchenko)
                                                   >>  39 // 17.11.2004 Change signature of methods, add new methods (V.Ivanchenko)
                                                   >>  40 // 08.04.2005 Major optimisation of internal interfaces (V.Ivantchenko)
                                                   >>  41 // 08.05.2005 Use updated interfaces (V.Ivantchenko)
                                                   >>  42 // 23.10.2005 Fix computations for ions (V.Ivantchenko)
 38 //                                                 43 //
 39 // Class Description: V.Ivanchenko & M.Novak   <<  44 // Class Description:
 40 //                                                 45 //
 41 // -------------------------------------------     46 // -------------------------------------------------------------------
 42 //                                                 47 //
 43 //....oooOO0OOooo........oooOO0OOooo........oo     48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 44 //....oooOO0OOooo........oooOO0OOooo........oo     49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 45                                                    50 
 46 #include "G4EmCalculator.hh"                       51 #include "G4EmCalculator.hh"
 47 #include "G4SystemOfUnits.hh"                  << 
 48 #include "G4LossTableManager.hh"                   52 #include "G4LossTableManager.hh"
 49 #include "G4EmParameters.hh"                   << 
 50 #include "G4NistManager.hh"                    << 
 51 #include "G4DynamicParticle.hh"                << 
 52 #include "G4VEmProcess.hh"                         53 #include "G4VEmProcess.hh"
 53 #include "G4VEnergyLossProcess.hh"                 54 #include "G4VEnergyLossProcess.hh"
 54 #include "G4VMultipleScattering.hh"                55 #include "G4VMultipleScattering.hh"
 55 #include "G4Material.hh"                           56 #include "G4Material.hh"
 56 #include "G4MaterialCutsCouple.hh"                 57 #include "G4MaterialCutsCouple.hh"
 57 #include "G4ParticleDefinition.hh"                 58 #include "G4ParticleDefinition.hh"
 58 #include "G4ParticleTable.hh"                      59 #include "G4ParticleTable.hh"
 59 #include "G4IonTable.hh"                       << 
 60 #include "G4PhysicsTable.hh"                       60 #include "G4PhysicsTable.hh"
 61 #include "G4ProductionCutsTable.hh"                61 #include "G4ProductionCutsTable.hh"
 62 #include "G4ProcessManager.hh"                     62 #include "G4ProcessManager.hh"
 63 #include "G4ionEffectiveCharge.hh"                 63 #include "G4ionEffectiveCharge.hh"
 64 #include "G4RegionStore.hh"                        64 #include "G4RegionStore.hh"
 65 #include "G4Element.hh"                            65 #include "G4Element.hh"
 66 #include "G4EmCorrections.hh"                      66 #include "G4EmCorrections.hh"
 67 #include "G4GenericIon.hh"                         67 #include "G4GenericIon.hh"
 68 #include "G4ProcessVector.hh"                  << 
 69 #include "G4Gamma.hh"                          << 
 70 #include "G4Electron.hh"                       << 
 71 #include "G4Positron.hh"                       << 
 72 #include "G4EmUtility.hh"                      << 
 73                                                    68 
 74 //....oooOO0OOooo........oooOO0OOooo........oo     69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 75                                                    70 
 76 G4EmCalculator::G4EmCalculator()                   71 G4EmCalculator::G4EmCalculator()
 77 {                                                  72 {
 78   manager = G4LossTableManager::Instance();        73   manager = G4LossTableManager::Instance();
 79   nist    = G4NistManager::Instance();         << 
 80   theParameters = G4EmParameters::Instance();  << 
 81   corr    = manager->EmCorrections();              74   corr    = manager->EmCorrections();
 82   cutenergy[0] = cutenergy[1] = cutenergy[2] = <<  75   nLocalMaterials    = 0;
 83   theGenericIon = G4GenericIon::GenericIon();  <<  76   verbose            = 0;
 84   ionEffCharge  = new G4ionEffectiveCharge();  <<  77   currentCoupleIndex = 0;
 85   dynParticle   = new G4DynamicParticle();     <<  78   currentCouple      = 0;
 86   ionTable      = G4ParticleTable::GetParticle <<  79   currentMaterial    = 0;
                                                   >>  80   currentParticle    = 0;
                                                   >>  81   baseParticle       = 0;
                                                   >>  82   currentLambda      = 0;
                                                   >>  83   chargeSquare       = 1.0;
                                                   >>  84   massRatio          = 1.0;
                                                   >>  85   currentParticleName= "";
                                                   >>  86   currentMaterialName= "";
                                                   >>  87   ionEffCharge       = new G4ionEffectiveCharge();
                                                   >>  88   isIon              = false;
                                                   >>  89   isApplicable       = false;
 87 }                                                  90 }
 88                                                    91 
 89 //....oooOO0OOooo........oooOO0OOooo........oo     92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 90                                                    93 
 91 G4EmCalculator::~G4EmCalculator()                  94 G4EmCalculator::~G4EmCalculator()
 92 {                                                  95 {
 93   delete ionEffCharge;                             96   delete ionEffCharge;
 94   delete dynParticle;                          <<  97   for (G4int i=0; i<nLocalMaterials; i++) {
 95   for (G4int i=0; i<nLocalMaterials; ++i) {    << 
 96     delete localCouples[i];                        98     delete localCouples[i];
 97   }                                                99   }
 98 }                                                 100 }
 99                                                   101 
100 //....oooOO0OOooo........oooOO0OOooo........oo    102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101                                                   103 
102 G4double G4EmCalculator::GetDEDX(G4double kinE << 104 G4double G4EmCalculator::GetDEDX(G4double kinEnergy, const G4ParticleDefinition* p,
103                                  const G4Parti << 105                                  const G4Material* mat, const G4Region* region)
104                                  const G4Mater << 
105                                  const G4Regio << 
106 {                                                 106 {
107   G4double res = 0.0;                             107   G4double res = 0.0;
108   const G4MaterialCutsCouple* couple = FindCou    108   const G4MaterialCutsCouple* couple = FindCouple(mat, region);
109   if(nullptr != couple && UpdateParticle(p, ki << 109   if(couple && UpdateParticle(p, kinEnergy) ) {
110     res = manager->GetDEDX(p, kinEnergy, coupl    110     res = manager->GetDEDX(p, kinEnergy, couple);
111                                                << 
112     if(isIon) {                                << 
113       if(FindEmModel(p, currentProcessName, ki << 
114         G4double length = CLHEP::nm;           << 
115         G4double eloss = res*length;           << 
116         //G4cout << "### GetDEDX: E= " << kinE << 
117         //       << " de= " << eloss << G4endl << 
118         dynParticle->SetKineticEnergy(kinEnerg << 
119         currentModel->GetChargeSquareRatio(p,  << 
120         currentModel->CorrectionsAlongStep(cou << 
121         res = eloss/length;                    << 
122              //G4cout << " de1= " << eloss <<  << 
123         //       << " " << p->GetParticleName( << 
124       }                                        << 
125     }                                          << 
126                                                << 
127     if(verbose>0) {                               111     if(verbose>0) {
128       G4cout << "G4EmCalculator::GetDEDX: E(Me    112       G4cout << "G4EmCalculator::GetDEDX: E(MeV)= " << kinEnergy/MeV
129              << " DEDX(MeV/mm)= " << res*mm/Me << 113        << " DEDX(MeV/mm)= " << res*mm/MeV
130              << " DEDX(MeV*cm^2/g)= " << res*g << 114        << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
131              << "  " <<  p->GetParticleName()  << 115        << "  " <<  p->GetParticleName()
132              << " in " <<  mat->GetName()      << 116        << " in " <<  mat->GetName()
133              << " isIon= " << isIon            << 117        << G4endl;
134              << G4endl;                        << 
135     }                                             118     }
136   }                                               119   }
137   return res;                                     120   return res;
138 }                                                 121 }
139                                                   122 
140 //....oooOO0OOooo........oooOO0OOooo........oo    123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
141                                                   124 
142 G4double G4EmCalculator::GetRangeFromRestricte << 125 G4double G4EmCalculator::GetDEDX(G4double kinEnergy, const G4String& particle,
143                                                << 126                                  const G4String& material, const G4String& reg)
144                                                << 
145                                                << 
146 {                                                 127 {
147   G4double res = 0.0;                          << 128   return GetDEDX(kinEnergy,FindParticle(particle),FindMaterial(material),FindRegion(reg));
148   const G4MaterialCutsCouple* couple = FindCou << 
149   if(couple && UpdateParticle(p, kinEnergy)) { << 
150     res = manager->GetRangeFromRestricteDEDX(p << 
151     if(verbose>1) {                            << 
152       G4cout << " G4EmCalculator::GetRangeFrom << 
153        << kinEnergy/MeV                        << 
154              << " range(mm)= " << res/mm       << 
155              << "  " <<  p->GetParticleName()  << 
156              << " in " <<  mat->GetName()      << 
157              << G4endl;                        << 
158     }                                          << 
159   }                                            << 
160   return res;                                  << 
161 }                                                 129 }
162                                                   130 
163 //....oooOO0OOooo........oooOO0OOooo........oo    131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
164                                                   132 
165 G4double G4EmCalculator::GetCSDARange(G4double << 133 G4double G4EmCalculator::GetRange(G4double kinEnergy, const G4ParticleDefinition* p,
166                                       const G4 << 134                                   const G4Material* mat, const G4Region* region)
167                                       const G4 << 
168                                       const G4 << 
169 {                                                 135 {
170   G4double res = 0.0;                             136   G4double res = 0.0;
171   if(!theParameters->BuildCSDARange()) {       << 
172     G4ExceptionDescription ed;                 << 
173     ed << "G4EmCalculator::GetCSDARange: CSDA  << 
174        << " use UI command: /process/eLoss/CSD << 
175     G4Exception("G4EmCalculator::GetCSDARange" << 
176                 JustWarning, ed);              << 
177     return res;                                << 
178   }                                            << 
179                                                << 
180   const G4MaterialCutsCouple* couple = FindCou    137   const G4MaterialCutsCouple* couple = FindCouple(mat,region);
181   if(nullptr != couple && UpdateParticle(p, ki << 138   if(couple && UpdateParticle(p, kinEnergy)) {
182     res = manager->GetCSDARange(p, kinEnergy,  << 139     res = manager->GetRange(p, kinEnergy, couple);
183     if(verbose>1) {                            << 140     if(verbose>0) {
184       G4cout << " G4EmCalculator::GetCSDARange << 141       G4cout << "G4EmCalculator::GetRange: E(MeV)= " << kinEnergy/MeV
185              << " range(mm)= " << res/mm       << 142        << " range(mm)= " << res/mm
186              << "  " <<  p->GetParticleName()  << 143        << "  " <<  p->GetParticleName()
187              << " in " <<  mat->GetName()      << 144        << " in " <<  mat->GetName()
188              << G4endl;                        << 145        << G4endl;
189     }                                             146     }
190   }                                               147   }
191   return res;                                     148   return res;
192 }                                                 149 }
193                                                   150 
194 //....oooOO0OOooo........oooOO0OOooo........oo    151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
195                                                   152 
196 G4double G4EmCalculator::GetRange(G4double kin << 153 G4double G4EmCalculator::GetRange(G4double kinEnergy, const G4String& particle,
197                                   const G4Part << 154                                   const G4String& material, const G4String& reg)
198                                   const G4Mate << 
199                                   const G4Regi << 
200 {                                                 155 {
201   G4double res = 0.0;                          << 156   return GetRange(kinEnergy,FindParticle(particle),FindMaterial(material),FindRegion(reg));
202   if(theParameters->BuildCSDARange()) {        << 
203     res = GetCSDARange(kinEnergy, p, mat, regi << 
204   } else {                                     << 
205     res = GetRangeFromRestricteDEDX(kinEnergy, << 
206   }                                            << 
207   return res;                                  << 
208 }                                                 157 }
209                                                   158 
210 //....oooOO0OOooo........oooOO0OOooo........oo    159 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
211                                                   160 
212 G4double G4EmCalculator::GetKinEnergy(G4double << 161 G4double G4EmCalculator::GetKinEnergy(G4double range, const G4ParticleDefinition* p,
213                                       const G4 << 162                                       const G4Material* mat, const G4Region* region)
214                                       const G4 << 
215                                       const G4 << 
216 {                                                 163 {
217   G4double res = 0.0;                             164   G4double res = 0.0;
218   const G4MaterialCutsCouple* couple = FindCou    165   const G4MaterialCutsCouple* couple = FindCouple(mat,region);
219   if(nullptr != couple && UpdateParticle(p, 1. << 166   if(couple && UpdateParticle(p, 1.0*GeV)) {
220     res = manager->GetEnergy(p, range, couple)    167     res = manager->GetEnergy(p, range, couple);
221     if(verbose>0) {                               168     if(verbose>0) {
222       G4cout << "G4EmCalculator::GetKinEnergy:    169       G4cout << "G4EmCalculator::GetKinEnergy: Range(mm)= " << range/mm
223              << " KinE(MeV)= " << res/MeV      << 170        << " KinE(MeV)= " << res/MeV
224              << "  " <<  p->GetParticleName()  << 171        << "  " <<  p->GetParticleName()
225              << " in " <<  mat->GetName()      << 172        << " in " <<  mat->GetName()
226              << G4endl;                        << 173        << G4endl;
227     }                                             174     }
228   }                                               175   }
229   return res;                                     176   return res;
230 }                                                 177 }
231                                                   178 
232 //....oooOO0OOooo........oooOO0OOooo........oo    179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
233                                                   180 
                                                   >> 181 G4double G4EmCalculator::GetKinEnergy(G4double range, const G4String& particle,
                                                   >> 182                                       const G4String& material, const G4String& reg)
                                                   >> 183 {
                                                   >> 184   return GetKinEnergy(range,FindParticle(particle),FindMaterial(material),FindRegion(reg));
                                                   >> 185 }
                                                   >> 186 
                                                   >> 187 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >> 188 
234 G4double G4EmCalculator::GetCrossSectionPerVol    189 G4double G4EmCalculator::GetCrossSectionPerVolume(G4double kinEnergy,
235                                             co    190                                             const G4ParticleDefinition* p,
236                                             co    191                                             const G4String& processName,
237                                             co << 192               const G4Material* mat,
238                                             co << 193               const G4Region* region)
239 {                                                 194 {
240   G4double res = 0.0;                             195   G4double res = 0.0;
241   const G4MaterialCutsCouple* couple = FindCou    196   const G4MaterialCutsCouple* couple = FindCouple(mat,region);
242                                                   197 
243   if(nullptr != couple && UpdateParticle(p, ki << 198   if(couple) {
244     if(FindEmModel(p, processName, kinEnergy)) << 199     G4int idx = couple->GetIndex();
245       G4int idx      = couple->GetIndex();     << 200     FindLambdaTable(p, processName);
246       G4int procType = -1;                     << 201     if(currentLambda && UpdateParticle(p, kinEnergy)) {
247       FindLambdaTable(p, processName, kinEnerg << 202       G4bool b;
248                                                << 203       G4double e = kinEnergy*massRatio;
249       G4VEmProcess* emproc = FindDiscreteProce << 204       res = (((*currentLambda)[idx])->GetValue(e,b))*chargeSquare;
250       if(nullptr != emproc) {                  << 205       res /= mat->GetTotNbOfAtomsPerVolume();
251   res = emproc->GetCrossSection(kinEnergy, cou << 
252       } else if(currentLambda) {               << 
253         // special tables are built for Msc mo << 
254   // procType is set in FindLambdaTable        << 
255         if(procType==2) {                      << 
256           auto mscM = static_cast<G4VMscModel* << 
257           mscM->SetCurrentCouple(couple);      << 
258           G4double tr1Mfp = mscM->GetTransport << 
259           if (tr1Mfp<DBL_MAX) {                << 
260             res = 1./tr1Mfp;                   << 
261           }                                    << 
262         } else {                               << 
263           G4double e = kinEnergy*massRatio;    << 
264           res = (((*currentLambda)[idx])->Valu << 
265         }                                      << 
266       } else {                                 << 
267         res = ComputeCrossSectionPerVolume(kin << 
268       }                                        << 
269       if(verbose>0) {                             206       if(verbose>0) {
270         G4cout << "G4EmCalculator::GetXSPerVol << 207   G4cout << "E(MeV)= " << kinEnergy/MeV
271                << " cross(cm-1)= " << res*cm   << 208          << " cross(cm^2/g)= " << res*gram/cm2
272                << "  " <<  p->GetParticleName( << 209          << "  " <<  p->GetParticleName()
273                << " in " <<  mat->GetName();   << 210          << " in " <<  mat->GetName()
274         if(verbose>1)                          << 211          << G4endl;
275           G4cout << "  idx= " << idx << "  Esc << 212       }
276            << kinEnergy*massRatio              << 
277            << "  q2= " << chargeSquare;        << 
278         G4cout << G4endl;                      << 
279       }                                        << 
280     }                                             213     }
281   }                                               214   }
282   return res;                                     215   return res;
283 }                                                 216 }
284                                                   217 
285 //....oooOO0OOooo........oooOO0OOooo........oo    218 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
286                                                   219 
287 G4double G4EmCalculator::GetShellIonisationCro << 220 G4double G4EmCalculator::GetCrossSectionPerVolume(G4double kinEnergy,
288                                          const << 221                                             const G4String& particle,
289                                          G4int << 222               const G4String& processName,
290                                          G4Ato << 223                                             const G4String& material,
291                                          G4dou << 224               const G4String& reg)
292 {                                                 225 {
293   G4double res = 0.0;                          << 226   return GetCrossSectionPerVolume(kinEnergy,FindParticle(particle),processName,
294   const G4ParticleDefinition* p = FindParticle << 227                                   FindMaterial(material),FindRegion(reg));
295   G4VAtomDeexcitation* ad = manager->AtomDeexc << 228 }
296   if(nullptr != p && nullptr != ad) {          << 229 
297     res = ad->GetShellIonisationCrossSectionPe << 230 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
298   }                                            << 231 
                                                   >> 232 G4double G4EmCalculator::GetCrossSectionPerAtom(G4double kinEnergy,
                                                   >> 233                                           const G4ParticleDefinition* p,
                                                   >> 234                                           const G4String& processName,
                                                   >> 235             const G4Material* mat,
                                                   >> 236             const G4Region* region)
                                                   >> 237 {
                                                   >> 238   G4double res = GetCrossSectionPerVolume(kinEnergy,p,processName,mat,region);
                                                   >> 239   if(mat) res /= mat->GetTotNbOfAtomsPerVolume();
299   return res;                                     240   return res;
300 }                                                 241 }
301                                                   242 
302 //....oooOO0OOooo........oooOO0OOooo........oo    243 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
303                                                   244 
                                                   >> 245 G4double G4EmCalculator::GetCrossSectionPerAtom(  G4double kinEnergy,
                                                   >> 246                                             const G4String& particle,
                                                   >> 247               const G4String& processName,
                                                   >> 248                                             const G4String& material,
                                                   >> 249               const G4String& reg)
                                                   >> 250 {
                                                   >> 251   return GetCrossSectionPerAtom(kinEnergy,FindParticle(particle),processName,
                                                   >> 252                                 FindMaterial(material),FindRegion(reg));
                                                   >> 253 }
                                                   >> 254 
                                                   >> 255 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >> 256 
304 G4double G4EmCalculator::GetMeanFreePath(G4dou    257 G4double G4EmCalculator::GetMeanFreePath(G4double kinEnergy,
305                                          const    258                                          const G4ParticleDefinition* p,
306                                          const    259                                          const G4String& processName,
307                                          const << 260            const G4Material* mat,
308                                          const    261                                          const G4Region* region)
309 {                                                 262 {
310   G4double res = DBL_MAX;                         263   G4double res = DBL_MAX;
311   G4double x = GetCrossSectionPerVolume(kinEne    264   G4double x = GetCrossSectionPerVolume(kinEnergy,p, processName, mat,region);
312   if(x > 0.0) { res = 1.0/x; }                 << 265   if(x > 0.0) res = 1.0/x;
313   if(verbose>1) {                                 266   if(verbose>1) {
314     G4cout << "G4EmCalculator::GetMeanFreePath    267     G4cout << "G4EmCalculator::GetMeanFreePath: E(MeV)= " << kinEnergy/MeV
315            << " MFP(mm)= " << res/mm           << 268      << " MFP(mm)= " << res/mm
316            << "  " <<  p->GetParticleName()    << 269      << "  " <<  p->GetParticleName()
317            << " in " <<  mat->GetName()        << 270      << " in " <<  mat->GetName()
318            << G4endl;                          << 271      << G4endl;
319   }                                               272   }
320   return res;                                     273   return res;
321 }                                                 274 }
322                                                   275 
323 //....oooOO0OOooo........oooOO0OOooo........oo    276 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
324                                                   277 
                                                   >> 278 G4double G4EmCalculator::GetMeanFreePath(G4double kinEnergy,
                                                   >> 279                                          const G4String& particle,
                                                   >> 280            const G4String& processName,
                                                   >> 281                                          const G4String& material,
                                                   >> 282            const G4String& reg)
                                                   >> 283 {
                                                   >> 284   return GetMeanFreePath(kinEnergy,FindParticle(particle),processName,
                                                   >> 285                          FindMaterial(material),FindRegion(reg));
                                                   >> 286 }
                                                   >> 287 
                                                   >> 288 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >> 289 
325 void G4EmCalculator::PrintDEDXTable(const G4Pa    290 void G4EmCalculator::PrintDEDXTable(const G4ParticleDefinition* p)
326 {                                                 291 {
327   const G4VEnergyLossProcess* elp = manager->G << 292   const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p);
328   G4cout << "##### DEDX Table for " << p->GetP    293   G4cout << "##### DEDX Table for " << p->GetParticleName() << G4endl;
329   if(nullptr != elp) G4cout << *(elp->DEDXTabl << 294   if(elp) G4cout << *(elp->DEDXTable()) << G4endl;
330 }                                                 295 }
331                                                   296 
332 //....oooOO0OOooo........oooOO0OOooo........oo    297 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
333                                                   298 
334 void G4EmCalculator::PrintRangeTable(const G4P    299 void G4EmCalculator::PrintRangeTable(const G4ParticleDefinition* p)
335 {                                                 300 {
336   const G4VEnergyLossProcess* elp = manager->G << 301   const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p);
337   G4cout << "##### Range Table for " << p->Get    302   G4cout << "##### Range Table for " << p->GetParticleName() << G4endl;
338   if(nullptr != elp) G4cout << *(elp->RangeTab << 303   if(elp) G4cout << *(elp->RangeTableForLoss()) << G4endl;
339 }                                                 304 }
340                                                   305 
341 //....oooOO0OOooo........oooOO0OOooo........oo    306 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
342                                                   307 
343 void G4EmCalculator::PrintInverseRangeTable(co    308 void G4EmCalculator::PrintInverseRangeTable(const G4ParticleDefinition* p)
344 {                                                 309 {
345   const G4VEnergyLossProcess* elp = manager->G << 310   const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p);
346   G4cout << "### G4EmCalculator: Inverse Range << 311   G4cout << "### G4EmCalculator: Inverse Range Table for " << p->GetParticleName() << G4endl;
347          << p->GetParticleName() << G4endl;    << 312   if(elp) G4cout << *(elp->InverseRangeTable()) << G4endl;
348   if(nullptr != elp) G4cout << *(elp->InverseR << 
349 }                                                 313 }
350                                                   314 
351 //....oooOO0OOooo........oooOO0OOooo........oo    315 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
352                                                   316 
353 G4double G4EmCalculator::ComputeDEDX(G4double     317 G4double G4EmCalculator::ComputeDEDX(G4double kinEnergy,
354                                      const G4P    318                                      const G4ParticleDefinition* p,
355                                      const G4S    319                                      const G4String& processName,
356                                      const G4M << 320              const G4Material* mat,
357                                            G4d    321                                            G4double cut)
358 {                                                 322 {
359   SetupMaterial(mat);                          << 323   currentMaterial = mat;
                                                   >> 324   currentMaterialName = mat->GetName();
360   G4double res = 0.0;                             325   G4double res = 0.0;
361   if(verbose > 1) {                               326   if(verbose > 1) {
362     G4cout << "### G4EmCalculator::ComputeDEDX    327     G4cout << "### G4EmCalculator::ComputeDEDX: " << p->GetParticleName()
363            << " in " << currentMaterialName       328            << " in " << currentMaterialName
364            << " e(MeV)= " << kinEnergy/MeV <<  << 329            << " e(MeV)= " << kinEnergy/MeV << "  cut(MeV)= " << cut/MeV << G4endl;
365            << G4endl;                          << 
366   }                                               330   }
367   if(UpdateParticle(p, kinEnergy)) {           << 331   if(FindEmModel(p, processName, kinEnergy)) {
368     if(FindEmModel(p, processName, kinEnergy)) << 332     //    G4cout << "currentModel= " << currentModel << G4endl;
                                                   >> 333     if(UpdateParticle(p, kinEnergy)) {
369       G4double escaled = kinEnergy*massRatio;     334       G4double escaled = kinEnergy*massRatio;
370       if(nullptr != baseParticle) {            << 335       if(baseParticle) {
371   res = currentModel->ComputeDEDXPerVolume(mat << 336         res = currentModel->ComputeDEDXPerVolume(mat, baseParticle, escaled, cut)
372                                                << 337             * chargeSquare;
373   if(verbose > 1) {                            << 338         if(verbose > 1)
374     G4cout << "Particle: " << p->GetParticleNa << 339           G4cout <<  baseParticle->GetParticleName() << " Escaled(MeV)= " << escaled;
375      << " E(MeV)=" << kinEnergy                << 
376      << " Base particle: " << baseParticle->Ge << 
377      << " Escaled(MeV)= " << escaled           << 
378      << " q2=" << chargeSquare << G4endl;      << 
379   }                                            << 
380       } else {                                    340       } else {
381   res = currentModel->ComputeDEDXPerVolume(mat << 341         res = currentModel->ComputeDEDXPerVolume(mat, p, kinEnergy, cut);
382   if(verbose > 1) {                            << 342         if(verbose > 1) G4cout <<  " no basePart E(MeV)= " << kinEnergy;
383     G4cout << "Particle: " << p->GetParticleNa << 
384      << " E(MeV)=" << kinEnergy << G4endl;     << 
385   }                                            << 
386       }                                           343       }
387       if(verbose > 1) {                        << 344       if(verbose > 1)
388   G4cout << currentModel->GetName() << ": DEDX << 345   G4cout << " DEDX(MeV/mm)= " << res*mm/MeV
389          << " DEDX(MeV*cm^2/g)= "              << 346          << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
390          << res*gram/(MeV*cm2*mat->GetDensity( << 
391          << G4endl;                               347          << G4endl;
                                                   >> 348       if(isIon && currentModel->HighEnergyLimit() > 100.*MeV) {
                                                   >> 349         res += corr->HighOrderCorrections(p,mat,kinEnergy);
                                                   >> 350   if(verbose > 1)
                                                   >> 351     G4cout << "After Corrections: DEDX(MeV/mm)= " << res*mm/MeV
                                                   >> 352      << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
                                                   >> 353      << G4endl;
392       }                                           354       }
393       // emulate smoothing procedure           << 355       // emulate boundary region for different parameterisations
394       if(applySmoothing && nullptr != loweMode << 356       G4double eth = currentModel->LowEnergyLimit();
395   G4double eth = currentModel->LowEnergyLimit( << 357       if(eth > 0.05*MeV && eth < 10.*MeV && escaled > eth) {
396   G4double res0 = 0.0;                         << 358         G4double res1 = 0.0;
397   G4double res1 = 0.0;                         << 359         if(baseParticle) {
398   if(nullptr != baseParticle) {                << 360           res1 = currentModel->ComputeDEDXPerVolume(mat, baseParticle, eth, cut)
399     res1 = chargeSquare*                       << 361                * chargeSquare;
400       currentModel->ComputeDEDXPerVolume(mat,  << 
401     res0 = chargeSquare*                       << 
402       loweModel->ComputeDEDXPerVolume(mat, bas << 
403   } else {                                        362   } else {
404     res1 = currentModel->ComputeDEDXPerVolume(    363     res1 = currentModel->ComputeDEDXPerVolume(mat, p, eth, cut);
405     res0 = loweModel->ComputeDEDXPerVolume(mat << 
406   }                                               364   }
407   if(res1 > 0.0 && escaled > 0.0) {            << 365   if(verbose > 1)
408     res *= (1.0 + (res0/res1 - 1.0)*eth/escale << 366     G4cout << "At boundary energy(MeV)= " << eth/MeV << " DEDX(MeV/mm)= " << res1*mm/MeV
409   }                                            << 367      << G4endl;
410   if(verbose > 1) {                            << 368         if(isIon) res1 += corr->HighOrderCorrections(p,mat,eth/massRatio);
411     G4cout << "At boundary energy(MeV)= " << e << 369         G4double res0 = res1;
412      << " DEDX(MeV/mm)= " << res0*mm/MeV << "  << 370         if(FindEmModel(p, processName, eth-1.0*keV)) {
413      << " after correction DEDX(MeV/mm)=" << r << 371     if(baseParticle) {
414   }                                            << 372       res0 = currentModel->ComputeDEDXPerVolume(mat, baseParticle, eth, cut)
415       }                                        << 373            * chargeSquare;
416       // correction for ions                   << 374     } else {
417       if(isIon) {                              << 375       res0 = currentModel->ComputeDEDXPerVolume(mat, p, eth, cut);
418   const G4double length = CLHEP::nm;           << 376     }
419   if(UpdateCouple(mat, cut)) {                 << 
420     G4double eloss = res*length;               << 
421     dynParticle->SetKineticEnergy(kinEnergy);  << 
422     currentModel->CorrectionsAlongStep(current << 
423                                              l << 
424     res = eloss/length;                        << 
425                                                << 
426     if(verbose > 1) {                          << 
427       G4cout << "After Corrections: DEDX(MeV/m << 
428        << " DEDX(MeV*cm^2/g)= "                << 
429        << res*gram/(MeV*cm2*mat->GetDensity()) << 
430     }                                          << 
431   }                                               377   }
                                                   >> 378         //G4cout << "eth= " << eth << " escaled= " << escaled << " res0= " << res0 << " res1= "
                                                   >> 379         //       << res1 <<  "  q2= " << chargeSquare << G4endl;
                                                   >> 380         res *= (1.0 + (res0/res1 - 1.0)*eth/escaled);
432       }                                           381       }
433       if(verbose > 0) {                           382       if(verbose > 0) {
434   G4cout << "## E(MeV)= " << kinEnergy/MeV     << 383         G4cout << "E(MeV)= " << kinEnergy/MeV
435          << " DEDX(MeV/mm)= " << res*mm/MeV    << 384                << " DEDX(MeV/mm)= " << res*mm/MeV
436          << " DEDX(MeV*cm^2/g)= " << res*gram/ << 385                << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
437          << " cut(MeV)= " << cut/MeV           << 386                << " cut(MeV)= " << cut/MeV
438          << "  " <<  p->GetParticleName()      << 387                << "  " <<  p->GetParticleName()
439          << " in " <<  currentMaterialName     << 388                << " in " <<  currentMaterialName
440          << " Zi^2= " << chargeSquare          << 389                << " Zi^2= " << chargeSquare
441          << " isIon=" << isIon                 << 390                << G4endl;
442          << G4endl;                            << 
443       }                                           391       }
444     }                                             392     }
445   }                                               393   }
446   return res;                                     394   return res;
447 }                                                 395 }
448                                                   396 
449 //....oooOO0OOooo........oooOO0OOooo........oo    397 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
450                                                   398 
451 G4double G4EmCalculator::ComputeElectronicDEDX << 399 G4double G4EmCalculator::ComputeDEDX(G4double kinEnergy,
452                                                << 400                                      const G4String& particle,
453                                                << 401              const G4String& processName,
454                                                << 402                                      const G4String& material,
455 {                                              << 403                                            G4double cut)
456   SetupMaterial(mat);                          << 404 {
457   G4double dedx = 0.0;                         << 405   return ComputeDEDX(kinEnergy,FindParticle(particle),processName,
458   if(UpdateParticle(part, kinEnergy)) {        << 406                      FindMaterial(material),cut);
459                                                << 
460     G4LossTableManager* lManager = G4LossTable << 
461     const std::vector<G4VEnergyLossProcess*> v << 
462       lManager->GetEnergyLossProcessVector();  << 
463     std::size_t n = vel.size();                << 
464                                                << 
465     //G4cout << "ComputeElectronicDEDX for " < << 
466     //           << " n= " << n << G4endl;     << 
467                                                << 
468     for(std::size_t i=0; i<n; ++i) {           << 
469       if(vel[i]) {                             << 
470         auto p = static_cast<G4VProcess*>(vel[ << 
471         if(ActiveForParticle(part, p)) {       << 
472           //G4cout << "idx= " << i << " " << ( << 
473           //  << "  " << (vel[i])->Particle()- << 
474           dedx += ComputeDEDX(kinEnergy,part,( << 
475         }                                      << 
476       }                                        << 
477     }                                          << 
478   }                                            << 
479   return dedx;                                 << 
480 }                                                 407 }
481                                                   408 
482 //....oooOO0OOooo........oooOO0OOooo........oo    409 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
483                                                   410 
484 G4double                                       << 411 G4double G4EmCalculator::ComputeNuclearDEDX(G4double kinEnergy,
485 G4EmCalculator::ComputeDEDXForCutInRange(G4dou << 412                                       const G4ParticleDefinition* p,
486                                          const << 413               const G4Material* mat)
487                                          const << 
488                                          G4dou << 
489 {                                                 414 {
490   SetupMaterial(mat);                          << 
491   G4double dedx = 0.0;                         << 
492   if(UpdateParticle(part, kinEnergy)) {        << 
493                                                << 
494     G4LossTableManager* lManager = G4LossTable << 
495     const std::vector<G4VEnergyLossProcess*> v << 
496       lManager->GetEnergyLossProcessVector();  << 
497     std::size_t n = vel.size();                << 
498                                                << 
499     if(mat != cutMaterial) {                   << 
500       cutMaterial = mat;                       << 
501       cutenergy[0] =                           << 
502         ComputeEnergyCutFromRangeCut(rangecut, << 
503       cutenergy[1] =                           << 
504         ComputeEnergyCutFromRangeCut(rangecut, << 
505       cutenergy[2] =                           << 
506         ComputeEnergyCutFromRangeCut(rangecut, << 
507     }                                          << 
508                                                   415 
509     //G4cout << "ComputeElectronicDEDX for " < << 416   G4double res = corr->NuclearDEDX(p, mat, kinEnergy, false);
510     //           << " n= " << n << G4endl;     << 
511                                                << 
512     for(std::size_t i=0; i<n; ++i) {           << 
513       if(vel[i]) {                             << 
514         auto p = static_cast<G4VProcess*>(vel[ << 
515         if(ActiveForParticle(part, p)) {       << 
516           //G4cout << "idx= " << i << " " << ( << 
517           // << "  " << (vel[i])->Particle()-> << 
518           const G4ParticleDefinition* sec = (v << 
519           std::size_t idx = 0;                 << 
520           if(sec == G4Electron::Electron()) {  << 
521           else if(sec == G4Positron::Positron( << 
522                                                   417 
523           dedx += ComputeDEDX(kinEnergy,part,( << 418   if(verbose > 1) {
524                               mat,cutenergy[id << 419     G4cout <<  p->GetParticleName() << " E(MeV)= " << kinEnergy/MeV
525         }                                      << 420      << " NuclearDEDX(MeV/mm)= " << res*mm/MeV
526       }                                        << 421      << " NuclearDEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
527     }                                          << 422      << G4endl;
528   }                                               423   }
529   return dedx;                                 << 424   return res;
530 }                                              << 
531                                                << 
532 //....oooOO0OOooo........oooOO0OOooo........oo << 
533                                                << 
534 G4double G4EmCalculator::ComputeTotalDEDX(G4do << 
535                                           cons << 
536                                           cons << 
537                                           G4do << 
538 {                                              << 
539   G4double dedx = ComputeElectronicDEDX(kinEne << 
540   if(mass > 700.*MeV) { dedx += ComputeNuclear << 
541   return dedx;                                 << 
542 }                                                 425 }
543                                                   426 
544 //....oooOO0OOooo........oooOO0OOooo........oo    427 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
545                                                   428 
546 G4double G4EmCalculator::ComputeNuclearDEDX(G4    429 G4double G4EmCalculator::ComputeNuclearDEDX(G4double kinEnergy,
547                                       const G4 << 430                                       const G4String& particle,
548                                       const G4 << 431               const G4String& material)
549 {                                                 432 {
550   G4double res = 0.0;                          << 433   return ComputeNuclearDEDX(kinEnergy,FindParticle(particle),FindMaterial(material));
551   G4VEmProcess* nucst = FindDiscreteProcess(p, << 
552   if(nucst) {                                  << 
553     G4VEmModel* mod = nucst->EmModel();        << 
554     if(mod) {                                  << 
555       mod->SetFluctuationFlag(false);          << 
556       res = mod->ComputeDEDXPerVolume(mat, p,  << 
557     }                                          << 
558   }                                            << 
559                                                << 
560   if(verbose > 1) {                            << 
561     G4cout <<  p->GetParticleName() << " E(MeV << 
562            << " NuclearDEDX(MeV/mm)= " << res* << 
563            << " NuclearDEDX(MeV*cm^2/g)= "     << 
564            << res*gram/(MeV*cm2*mat->GetDensit << 
565            << G4endl;                          << 
566   }                                            << 
567   return res;                                  << 
568 }                                                 434 }
569                                                   435 
570 //....oooOO0OOooo........oooOO0OOooo........oo    436 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
571                                                   437 
572 G4double G4EmCalculator::ComputeCrossSectionPe    438 G4double G4EmCalculator::ComputeCrossSectionPerVolume(
573                                                   439                                                    G4double kinEnergy,
574                                              c    440                                              const G4ParticleDefinition* p,
575                                              c    441                                              const G4String& processName,
576                                              c << 442                const G4Material* mat,
577                                                   443                                                    G4double cut)
578 {                                                 444 {
579   SetupMaterial(mat);                          << 445   currentMaterial = mat;
                                                   >> 446   currentMaterialName = mat->GetName();
580   G4double res = 0.0;                             447   G4double res = 0.0;
581   if(UpdateParticle(p, kinEnergy)) {           << 448   if(FindEmModel(p, processName, kinEnergy)) {
582     if(FindEmModel(p, processName, kinEnergy)) << 449     UpdateParticle(p, kinEnergy);
583       G4double e = kinEnergy;                  << 450     G4double e = kinEnergy;
584       G4double aCut = std::max(cut, theParamet << 451     if(baseParticle) {
585       if(baseParticle) {                       << 452       e *= kinEnergy*massRatio;
586         e *= kinEnergy*massRatio;              << 453       res = currentModel->CrossSectionPerVolume(mat, baseParticle, e, cut, e)*chargeSquare;
587         res = currentModel->CrossSectionPerVol << 454     } else {
588               mat, baseParticle, e, aCut, e) * << 455       res = currentModel->CrossSectionPerVolume(mat, p, e, cut, e);
589       } else {                                 << 456     }
590         res = currentModel->CrossSectionPerVol << 457     if(verbose>0) {
591       }                                        << 458       G4cout << "E(MeV)= " << kinEnergy/MeV
592       if(verbose>0) {                          << 459        << " cross(cm^2/g)= " << res*gram/cm2
593         G4cout << "G4EmCalculator::ComputeXSPe << 460        << "  " <<  p->GetParticleName()
594                << kinEnergy/MeV                << 461        << " in " <<  mat->GetName()
595                << " cross(cm-1)= " << res*cm   << 462        << G4endl;
596                << " cut(keV)= " << aCut/keV    << 
597                << "  " <<  p->GetParticleName( << 
598                << " in " <<  mat->GetName()    << 
599                << G4endl;                      << 
600       }                                        << 
601     }                                             463     }
602   }                                               464   }
603   return res;                                     465   return res;
604 }                                                 466 }
605                                                   467 
606 //....oooOO0OOooo........oooOO0OOooo........oo    468 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
607                                                   469 
608 G4double                                       << 470 G4double G4EmCalculator::ComputeCrossSectionPerVolume(
609 G4EmCalculator::ComputeCrossSectionPerAtom(G4d << 471                                                    G4double kinEnergy,
610                                            con << 472                                              const G4String& particle,
611                                            con << 473                const G4String& processName,
612                                            G4d << 474                                              const G4String& material,
613                                            G4d << 475                                                    G4double cut)
614 {                                                 476 {
615   G4double res = 0.0;                          << 477   return ComputeCrossSectionPerVolume(kinEnergy,FindParticle(particle),processName,
616   if(UpdateParticle(p, kinEnergy)) {           << 478                                       FindMaterial(material),cut);
617     G4int iz = G4lrint(Z);                     << 
618     CheckMaterial(iz);                         << 
619     if(FindEmModel(p, processName, kinEnergy)) << 
620       G4double e = kinEnergy;                  << 
621       G4double aCut = std::max(cut, theParamet << 
622       if(baseParticle) {                       << 
623         e *= kinEnergy*massRatio;              << 
624         currentModel->InitialiseForElement(bas << 
625         res = currentModel->ComputeCrossSectio << 
626               baseParticle, e, Z, A, aCut) * c << 
627       } else {                                 << 
628         currentModel->InitialiseForElement(p,  << 
629         res = currentModel->ComputeCrossSectio << 
630       }                                        << 
631       if(verbose>0) {                          << 
632         G4cout << "E(MeV)= " << kinEnergy/MeV  << 
633                << " cross(barn)= " << res/barn << 
634                << "  " <<  p->GetParticleName( << 
635                << " Z= " <<  Z << " A= " << A/ << 
636                << " cut(keV)= " << aCut/keV    << 
637                << G4endl;                      << 
638       }                                        << 
639     }                                          << 
640   }                                            << 
641   return res;                                  << 
642 }                                                 479 }
643                                                   480 
644 //....oooOO0OOooo........oooOO0OOooo........oo    481 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
645                                                   482 
646 G4double                                       << 483 G4double G4EmCalculator::ComputeCrossSectionPerAtom(
647 G4EmCalculator::ComputeCrossSectionPerShell(G4 << 484                                                    G4double kinEnergy,
648                                             co << 485                const G4ParticleDefinition* p,
649                                             co << 486                                              const G4String& processName,
650                                             G4 << 487                      G4double Z, G4double A,
651                                             G4 << 488                                        G4double cut)
652 {                                                 489 {
653   G4double res = 0.0;                             490   G4double res = 0.0;
654   if(UpdateParticle(p, kinEnergy)) {           << 491   if(FindEmModel(p, processName, kinEnergy)) {
655     CheckMaterial(Z);                          << 492     UpdateParticle(p, kinEnergy);
656     if(FindEmModel(p, processName, kinEnergy)) << 493     G4double e = kinEnergy;
657       G4double e = kinEnergy;                  << 494     if(baseParticle) {
658       G4double aCut = std::max(cut, theParamet << 495       e *= kinEnergy*massRatio;
659       if(nullptr != baseParticle) {            << 496       res = currentModel->ComputeCrossSectionPerAtom(baseParticle, e, Z, A, cut)*chargeSquare;
660         e *= kinEnergy*massRatio;              << 497     } else {
661         currentModel->InitialiseForElement(bas << 498       res = currentModel->ComputeCrossSectionPerAtom(p, e, Z, A, cut);
662         res =                                  << 499     }
663           currentModel->ComputeCrossSectionPer << 500     if(verbose>0) {
664                                                << 501       G4cout << "E(MeV)= " << kinEnergy/MeV
665       } else {                                 << 502        << " cross(barn)= " << res/barn
666         currentModel->InitialiseForElement(p,  << 503        << "  " <<  p->GetParticleName()
667         res = currentModel->ComputeCrossSectio << 504        << " Z= " <<  Z << " A= " << A
668       }                                        << 505        << G4endl;
669       if(verbose>0) {                          << 
670         G4cout << "E(MeV)= " << kinEnergy/MeV  << 
671                << " cross(barn)= " << res/barn << 
672                << "  " <<  p->GetParticleName( << 
673                << " Z= " <<  Z << " shellIdx=  << 
674                << " cut(keV)= " << aCut/keV    << 
675          << G4endl;                            << 
676       }                                        << 
677     }                                             506     }
678   }                                               507   }
679   return res;                                     508   return res;
680 }                                                 509 }
681                                                   510 
682 //....oooOO0OOooo........oooOO0OOooo........oo    511 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
683                                                   512 
684 G4double                                       << 513 G4double G4EmCalculator::ComputeCrossSectionPerAtom(G4double kinEnergy,
685 G4EmCalculator::ComputeGammaAttenuationLength( << 514                                               const G4String& particle,
686                                                << 515                                               const G4String& processName,
                                                   >> 516                 const G4Element* elm,
                                                   >> 517                                         G4double cut)
687 {                                                 518 {
688   G4double res = 0.0;                          << 519   return ComputeCrossSectionPerAtom(kinEnergy,FindParticle(particle),processName,
689   const G4ParticleDefinition* gamma = G4Gamma: << 520                                     elm->GetZ(),elm->GetN(),cut);
690   res += ComputeCrossSectionPerVolume(kinEnerg << 
691   res += ComputeCrossSectionPerVolume(kinEnerg << 
692   res += ComputeCrossSectionPerVolume(kinEnerg << 
693   res += ComputeCrossSectionPerVolume(kinEnerg << 
694   if(res > 0.0) { res = 1.0/res; }             << 
695   return res;                                  << 
696 }                                              << 
697                                                << 
698 //....oooOO0OOooo........oooOO0OOooo........oo << 
699                                                << 
700 G4double G4EmCalculator::ComputeShellIonisatio << 
701                                          const << 
702                                          G4int << 
703                                          G4Ato << 
704                                          G4dou << 
705                                          const << 
706 {                                              << 
707   G4double res = 0.0;                          << 
708   const G4ParticleDefinition* p = FindParticle << 
709   G4VAtomDeexcitation* ad = manager->AtomDeexc << 
710   if(p && ad) {                                << 
711     res = ad->ComputeShellIonisationCrossSecti << 
712                                                << 
713   }                                            << 
714   return res;                                  << 
715 }                                                 521 }
716                                                   522 
717 //....oooOO0OOooo........oooOO0OOooo........oo    523 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
718                                                   524 
719 G4double G4EmCalculator::ComputeMeanFreePath(G    525 G4double G4EmCalculator::ComputeMeanFreePath(G4double kinEnergy,
720                                              c    526                                              const G4ParticleDefinition* p,
721                                              c    527                                              const G4String& processName,
722                                              c << 528                const G4Material* mat,
723                                              G << 529                                                    G4double cut)
724 {                                                 530 {
725   G4double mfp = DBL_MAX;                         531   G4double mfp = DBL_MAX;
726   G4double x =                                 << 532   G4double x = ComputeCrossSectionPerVolume(kinEnergy, p, processName, mat, cut);
727     ComputeCrossSectionPerVolume(kinEnergy, p, << 533   if(x > 0.0) mfp = 1.0/x;
728   if(x > 0.0) { mfp = 1.0/x; }                 << 
729   if(verbose>1) {                                 534   if(verbose>1) {
730     G4cout << "E(MeV)= " << kinEnergy/MeV         535     G4cout << "E(MeV)= " << kinEnergy/MeV
731            << " MFP(mm)= " << mfp/mm           << 536      << " MFP(mm)= " << mfp/mm
732            << "  " <<  p->GetParticleName()    << 537      << "  " <<  p->GetParticleName()
733            << " in " <<  mat->GetName()        << 538      << " in " <<  mat->GetName()
734            << G4endl;                          << 539      << G4endl;
735   }                                               540   }
736   return mfp;                                     541   return mfp;
737 }                                                 542 }
738                                                   543 
739 //....oooOO0OOooo........oooOO0OOooo........oo    544 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
740                                                   545 
741 G4double G4EmCalculator::ComputeEnergyCutFromR << 546 G4double G4EmCalculator::ComputeMeanFreePath(G4double kinEnergy,
742                          G4double range,       << 547                                              const G4String& particle,
743                          const G4ParticleDefin << 548                                              const G4String& processName,
744                          const G4Material* mat << 549                                              const G4String& material,
                                                   >> 550                                                    G4double cut)
745 {                                                 551 {
746   return G4ProductionCutsTable::GetProductionC << 552   return ComputeMeanFreePath(kinEnergy,FindParticle(particle),processName,
747     ConvertRangeToEnergy(part, mat, range);    << 553                              FindMaterial(material),cut);
748 }                                                 554 }
749                                                   555 
750 //....oooOO0OOooo........oooOO0OOooo........oo    556 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
751                                                   557 
752 G4bool G4EmCalculator::UpdateParticle(const G4 << 558 G4bool G4EmCalculator::UpdateParticle(const G4ParticleDefinition* p, G4double kinEnergy)
753                                       G4double << 
754 {                                                 559 {
755   if(p != currentParticle) {                      560   if(p != currentParticle) {
756                                                << 
757     // new particle                            << 
758     currentParticle = p;                          561     currentParticle = p;
759     dynParticle->SetDefinition(const_cast<G4Pa << 562     baseParticle    = 0;
760     dynParticle->SetKineticEnergy(kinEnergy);  << 
761     baseParticle    = nullptr;                 << 
762     currentParticleName = p->GetParticleName()    563     currentParticleName = p->GetParticleName();
763     massRatio       = 1.0;                        564     massRatio       = 1.0;
764     mass            = p->GetPDGMass();         << 
765     chargeSquare    = 1.0;                        565     chargeSquare    = 1.0;
766     currentProcess  = manager->GetEnergyLossPr << 566     if(p->GetParticleType() == "nucleus" &&
767     currentProcessName = "";                   << 567        currentParticleName != "deuteron" && currentParticleName != "triton") {
768     isIon = false;                             << 568       baseParticle = G4GenericIon::GenericIon();
769                                                << 569       massRatio = baseParticle->GetPDGMass()/p->GetPDGMass();
770     // ionisation process exist                << 570       isIon = true;
771     if(nullptr != currentProcess) {            << 571       //      G4cout << p->GetParticleName() << " in " << currentMaterial->GetName()
772       currentProcessName = currentProcess->Get << 572       //       << "  e= " << kinEnergy << G4endl;
773       baseParticle = currentProcess->BaseParti << 573       chargeSquare =
774       if(currentProcessName == "ionIoni" && p- << 574           ionEffCharge->EffectiveChargeSquareRatio(p, currentMaterial, kinEnergy);
775         baseParticle = theGenericIon;          << 575       //G4cout << "q2= " << chargeSquare << G4endl;
776         isIon = true;                          << 576     } else {
                                                   >> 577       isIon = false;
                                                   >> 578       const G4VEnergyLossProcess* proc = FindEnergyLossProcess(p);
                                                   >> 579       if(proc) {
                                                   >> 580         baseParticle    = proc->BaseParticle();
                                                   >> 581 
                                                   >> 582   if(baseParticle) {
                                                   >> 583           massRatio = baseParticle->GetPDGMass()/p->GetPDGMass();
                                                   >> 584           G4double q = p->GetPDGCharge()/baseParticle->GetPDGCharge();
                                                   >> 585     chargeSquare *= q*q;
                                                   >> 586         }
777       }                                           587       }
778                                                << 
779       // base particle is used                 << 
780       if(nullptr != baseParticle) {            << 
781         massRatio = baseParticle->GetPDGMass() << 
782         G4double q = p->GetPDGCharge()/basePar << 
783         chargeSquare = q*q;                    << 
784       }                                        << 
785     }                                          << 
786   }                                            << 
787   // Effective charge for ions                 << 
788   if(isIon && nullptr != currentProcess) {     << 
789     chargeSquare =                             << 
790       corr->EffectiveChargeSquareRatio(p, curr << 
791     currentProcess->SetDynamicMassCharge(massR << 
792     if(verbose>1) {                            << 
793       G4cout <<"\n NewIon: massR= "<< massRati << 
794        << chargeSquare << "  " << currentProce << 
795     }                                             588     }
796   }                                               589   }
                                                   >> 590   if(isIon) chargeSquare =
                                                   >> 591      ionEffCharge->EffectiveChargeSquareRatio(p, currentMaterial, kinEnergy);
797   return true;                                    592   return true;
798 }                                                 593 }
799                                                   594 
800 //....oooOO0OOooo........oooOO0OOooo........oo    595 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
801                                                   596 
802 const G4ParticleDefinition* G4EmCalculator::Fi    597 const G4ParticleDefinition* G4EmCalculator::FindParticle(const G4String& name)
803 {                                                 598 {
804   const G4ParticleDefinition* p = nullptr;     << 599   const G4ParticleDefinition* p = 0;
805   if(name != currentParticleName) {               600   if(name != currentParticleName) {
806     p = G4ParticleTable::GetParticleTable()->F    601     p = G4ParticleTable::GetParticleTable()->FindParticle(name);
807     if(nullptr == p) {                         << 602     if(!p) {
808       G4cout << "### WARNING: G4EmCalculator:: << 603       G4cout << "### WARNING: G4EmCalculator::FindParticle fails to find " << name << G4endl;
809              << name << G4endl;                << 
810     }                                             604     }
811   } else {                                        605   } else {
812     p = currentParticle;                          606     p = currentParticle;
813   }                                               607   }
814   return p;                                       608   return p;
815 }                                                 609 }
816                                                   610 
817 //....oooOO0OOooo........oooOO0OOooo........oo    611 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
818                                                   612 
819 const G4ParticleDefinition* G4EmCalculator::Fi << 
820 {                                              << 
821   const G4ParticleDefinition* p = ionTable->Ge << 
822   return p;                                    << 
823 }                                              << 
824                                                << 
825 //....oooOO0OOooo........oooOO0OOooo........oo << 
826                                                << 
827 const G4Material* G4EmCalculator::FindMaterial    613 const G4Material* G4EmCalculator::FindMaterial(const G4String& name)
828 {                                                 614 {
829   if(name != currentMaterialName) {               615   if(name != currentMaterialName) {
830     SetupMaterial(G4Material::GetMaterial(name << 616     currentMaterial = G4Material::GetMaterial(name);
831     if(nullptr == currentMaterial) {           << 617     currentMaterialName = name;
832       G4cout << "### WARNING: G4EmCalculator:: << 618     if(!currentMaterial)
833              << name << G4endl;                << 619       G4cout << "### WARNING: G4EmCalculator::FindMaterial fails to find " << name << G4endl;
834     }                                          << 
835   }                                               620   }
836   return currentMaterial;                         621   return currentMaterial;
837 }                                                 622 }
838                                                   623 
839 //....oooOO0OOooo........oooOO0OOooo........oo    624 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
840                                                   625 
841 const G4Region* G4EmCalculator::FindRegion(con    626 const G4Region* G4EmCalculator::FindRegion(const G4String& reg)
842 {                                                 627 {
843   return G4EmUtility::FindRegion(reg);         << 628   return G4RegionStore::GetInstance()->GetRegion(reg);
844 }                                                 629 }
845                                                   630 
846 //....oooOO0OOooo........oooOO0OOooo........oo    631 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
847                                                   632 
848 const G4MaterialCutsCouple* G4EmCalculator::Fi << 633 const G4MaterialCutsCouple* G4EmCalculator::FindCouple(const G4Material* material,
849                             const G4Material*  << 634                                                        const G4Region* region)
850                             const G4Region* re << 
851 {                                                 635 {
852   const G4MaterialCutsCouple* couple = nullptr << 636   if(!material) return 0;
853   SetupMaterial(material);                     << 637   currentMaterial = material;
854   if(nullptr != currentMaterial) {             << 638   currentMaterialName = material->GetName();
855     // Access to materials                     << 639   // Access to materials
856     const G4ProductionCutsTable* theCoupleTabl << 640   const G4ProductionCutsTable* theCoupleTable=
857       G4ProductionCutsTable::GetProductionCuts << 641         G4ProductionCutsTable::GetProductionCutsTable();
858     const G4Region* r = region;                << 642   const G4Region* r = region;
859     if(nullptr != r) {                         << 643   if(!r) r = G4RegionStore::GetInstance()->GetRegion("DefaultRegionForTheWorld");
860       couple = theCoupleTable->GetMaterialCuts << 644 
861                                                << 645   return theCoupleTable->GetMaterialCutsCouple(material,r->GetProductionCuts());
862     } else {                                   << 646 
863       G4RegionStore* store = G4RegionStore::Ge << 
864       std::size_t nr = store->size();          << 
865       if(0 < nr) {                             << 
866         for(std::size_t i=0; i<nr; ++i) {      << 
867           couple = theCoupleTable->GetMaterial << 
868             material, ((*store)[i])->GetProduc << 
869           if(nullptr != couple) { break; }     << 
870         }                                      << 
871       }                                        << 
872     }                                          << 
873   }                                            << 
874   if(nullptr == couple) {                      << 
875     G4ExceptionDescription ed;                 << 
876     ed << "G4EmCalculator::FindCouple: fail fo << 
877        << currentMaterialName << ">";          << 
878     if(region) { ed << " and region " << regio << 
879     G4Exception("G4EmCalculator::FindCouple",  << 
880                 FatalException, ed);           << 
881   }                                            << 
882   return couple;                               << 
883 }                                                 647 }
884                                                   648 
885 //....oooOO0OOooo........oooOO0OOooo........oo    649 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
886                                                   650 
887 G4bool G4EmCalculator::UpdateCouple(const G4Ma    651 G4bool G4EmCalculator::UpdateCouple(const G4Material* material, G4double cut)
888 {                                                 652 {
889   SetupMaterial(material);                     << 653   if(!material) return false;
890   if(!currentMaterial) { return false; }       << 654   currentMaterial = material;
891   for (G4int i=0; i<nLocalMaterials; ++i) {    << 655   currentMaterialName = material->GetName();
                                                   >> 656   for (G4int i=0; i<nLocalMaterials; i++) {
892     if(material == localMaterials[i] && cut ==    657     if(material == localMaterials[i] && cut == localCuts[i]) {
893       currentCouple = localCouples[i];            658       currentCouple = localCouples[i];
894       currentCoupleIndex = currentCouple->GetI    659       currentCoupleIndex = currentCouple->GetIndex();
895       currentCut = cut;                           660       currentCut = cut;
896       return true;                                661       return true;
897     }                                             662     }
898   }                                               663   }
899   const G4MaterialCutsCouple* cc = new G4Mater    664   const G4MaterialCutsCouple* cc = new G4MaterialCutsCouple(material);
900   localMaterials.push_back(material);             665   localMaterials.push_back(material);
901   localCouples.push_back(cc);                     666   localCouples.push_back(cc);
902   localCuts.push_back(cut);                       667   localCuts.push_back(cut);
903   ++nLocalMaterials;                           << 668   nLocalMaterials++;
904   currentCouple = cc;                             669   currentCouple = cc;
905   currentCoupleIndex = currentCouple->GetIndex    670   currentCoupleIndex = currentCouple->GetIndex();
906   currentCut = cut;                               671   currentCut = cut;
907   return true;                                    672   return true;
908 }                                                 673 }
909                                                   674 
910 //....oooOO0OOooo........oooOO0OOooo........oo    675 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
911                                                   676 
912 void G4EmCalculator::FindLambdaTable(const G4P    677 void G4EmCalculator::FindLambdaTable(const G4ParticleDefinition* p,
913                                      const G4S << 678                                      const G4String& processName)
914                                      G4double  << 
915 {                                                 679 {
916   // Search for the process                       680   // Search for the process
917   if (!currentLambda || p != lambdaParticle || << 681   if (p != currentParticle || processName != currentName) {
918     lambdaName     = processName;              << 682     currentName     = processName;
919     currentLambda  = nullptr;                  << 683     currentLambda   = 0;
920     lambdaParticle = p;                        << 684 
921     isApplicable   = false;                    << 685     G4String partname =  p->GetParticleName();
922                                                << 686     const G4ParticleDefinition* part = p;
923     const G4ParticleDefinition* part = (isIon) << 687     if(p->GetParticleType() == "nucleus" && partname != "deuteron" && partname != "triton")
924                                                << 688       part = G4GenericIon::GenericIon();
925     // Search for energy loss process          << 689 
926     currentName = processName;                 << 
927     currentModel = nullptr;                    << 
928     loweModel = nullptr;                       << 
929                                                << 
930     G4VEnergyLossProcess* elproc = FindEnLossP << 
931     if(nullptr != elproc) {                    << 
932       currentLambda = elproc->LambdaTable();   << 
933       proctype      = 0;                       << 
934       if(nullptr != currentLambda) {           << 
935         isApplicable = true;                   << 
936         if(verbose>1) {                        << 
937           G4cout << "G4VEnergyLossProcess is f << 
938                  << G4endl;                    << 
939         }                                      << 
940       }                                        << 
941       curProcess = elproc;                     << 
942       return;                                  << 
943     }                                          << 
944                                                   690 
945     // Search for discrete process             << 691     G4LossTableManager* lManager = G4LossTableManager::Instance();
946     G4VEmProcess* proc = FindDiscreteProcess(p << 692     const std::vector<G4VEnergyLossProcess*> vel = lManager->GetEnergyLossProcessVector();
947     if(nullptr != proc) {                      << 693     G4int n = vel.size();
948       currentLambda = proc->LambdaTable();     << 694     for(G4int i=0; i<n; i++) {
949       proctype      = 1;                       << 695       if((vel[i])->GetProcessName() == currentName && (vel[i])->Particle() == part) {
950       if(nullptr != currentLambda) {           << 696         currentLambda = (vel[i])->LambdaTable();
951         isApplicable = true;                   << 697   isApplicable    = true;
952         if(verbose>1) {                        << 698   break;
953           G4cout << "G4VEmProcess is found out << 699       }
                                                   >> 700     }
                                                   >> 701     if(!currentLambda) {
                                                   >> 702       const std::vector<G4VEmProcess*> vem = lManager->GetEmProcessVector();
                                                   >> 703       G4int n = vem.size();
                                                   >> 704       for(G4int i=0; i<n; i++) {
                                                   >> 705         if((vem[i])->GetProcessName() == currentName && (vem[i])->Particle() == part) {
                                                   >> 706           currentLambda = (vem[i])->LambdaTable();
                                                   >> 707           isApplicable    = true;
                                                   >> 708     break;
954         }                                         709         }
955       }                                           710       }
956       curProcess = proc;                       << 
957       return;                                  << 
958     }                                             711     }
959                                                << 712     if(!currentLambda) {
960     // Search for msc process                  << 713       const std::vector<G4VMultipleScattering*> vmsc = lManager->GetMultipleScatteringVector();
961     G4VMultipleScattering* msc = FindMscProces << 714       G4int n = vmsc.size();
962     if(nullptr != msc) {                       << 715       for(G4int i=0; i<n; i++) {
963       currentModel = msc->SelectModel(kinEnerg << 716         if((vmsc[i])->GetProcessName() == currentName && (vmsc[i])->Particle() == part) {
964       proctype     = 2;                        << 717           currentLambda = (vmsc[i])->LambdaTable();
965       if(nullptr != currentModel) {            << 718     isApplicable    = true;
966         currentLambda = currentModel->GetCross << 719     break;
967         if(nullptr != currentLambda) {         << 
968           isApplicable = true;                 << 
969           if(verbose>1) {                      << 
970             G4cout << "G4VMultipleScattering i << 
971                    << G4endl;                  << 
972           }                                    << 
973         }                                         720         }
974       }                                           721       }
975       curProcess = msc;                        << 
976     }                                             722     }
977   }                                               723   }
978 }                                                 724 }
979                                                   725 
980 //....oooOO0OOooo........oooOO0OOooo........oo    726 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
981                                                   727 
982 G4bool G4EmCalculator::FindEmModel(const G4Par    728 G4bool G4EmCalculator::FindEmModel(const G4ParticleDefinition* p,
983                                    const G4Str    729                                    const G4String& processName,
984                                             G4 << 730                    G4double kineticEnergy)
985 {                                                 731 {
986   isApplicable = false;                        << 732   G4bool res = false;
987   if(nullptr == p || nullptr == currentMateria << 733   if(!p) {
988     G4cout << "G4EmCalculator::FindEmModel WAR << 734     G4cout << "G4EmCalculator::FindEmModel WARNING: no particle defined" << G4endl;
989            << " or materail defined; particle: << 735     return res;
990     return isApplicable;                       << 
991   }                                               736   }
992   G4String partname =  p->GetParticleName();      737   G4String partname =  p->GetParticleName();
993   G4double scaledEnergy = kinEnergy*massRatio; << 738   const G4ParticleDefinition* part = p;
994   const G4ParticleDefinition* part = (isIon) ? << 739   G4double kinEnergy = kineticEnergy;
                                                   >> 740   if(p->GetParticleType() == "nucleus" && partname != "deuteron" && partname != "triton") {
                                                   >> 741     part = G4GenericIon::GenericIon(); 
                                                   >> 742     kinEnergy *= part->GetPDGMass()/p->GetPDGMass();
                                                   >> 743   }
995                                                   744 
996   if(verbose > 1) {                               745   if(verbose > 1) {
997     G4cout << "## G4EmCalculator::FindEmModel  << 746     G4cout << "G4EmCalculator::FindEmModel for " << partname
998            << " (type= " << p->GetParticleType    747            << " (type= " << p->GetParticleType()
999            << ") and " << processName << " at  << 748            << ") and " << processName << " at e(MeV)= " << kinEnergy;
1000            << G4endl;                         << 749     if(p != part) G4cout << "  GenericIon is the base particle";       
1001     if(p != part) { G4cout << "  GenericIon i << 750     G4cout << G4endl;
1002   }                                              751   }
1003                                                  752 
1004   // Search for energy loss process           << 753   // Search for the process
1005   currentName = processName;                     754   currentName = processName;
1006   currentModel = nullptr;                     << 755   currentModel = 0;
1007   loweModel = nullptr;                        << 756   size_t idx   = 0;
1008   std::size_t idx = 0;                        << 757   G4LossTableManager* lManager = G4LossTableManager::Instance();
1009                                               << 758   const std::vector<G4VEnergyLossProcess*> vel = lManager->GetEnergyLossProcessVector();
1010   G4VEnergyLossProcess* elproc = FindEnLossPr << 759   G4int n = vel.size();
1011   if(nullptr != elproc) {                     << 760   for(G4int i=0; i<n; i++) {
1012     currentModel = elproc->SelectModelForMate << 761     //    G4cout << "i= " << i << " part= " << (vel[i])->Particle()->GetParticleName()
1013     currentModel->InitialiseForMaterial(part, << 762     //     << "   proc= " << (vel[i])->GetProcessName()  << G4endl;
1014     currentModel->SetupForMaterial(part, curr << 763     if((vel[i])->GetProcessName() == currentName && (vel[i])->Particle() == part) {
1015     G4double eth = currentModel->LowEnergyLim << 764       const G4ParticleDefinition* bp = (vel[i])->BaseParticle();
1016     if(eth > 0.0) {                           << 765       //      G4cout << "i= " << i << " bp= " << bp << G4endl;
1017       loweModel = elproc->SelectModelForMater << 766       if(!bp) {
1018       if(loweModel == currentModel) { loweMod << 767         currentModel = (vel[i])->SelectModelForMaterial(kinEnergy, idx);
1019       else {                                  << 768       } else {
1020         loweModel->InitialiseForMaterial(part << 769         G4double e = kinEnergy*bp->GetPDGMass()/p->GetPDGMass();
1021         loweModel->SetupForMaterial(part, cur << 770         for(G4int j=0; j<n; j++) {
1022       }                                       << 771           if((vel[j])->Particle() == bp) {
1023     }                                         << 772             currentModel = (vel[j])->SelectModelForMaterial(e, idx);
1024   }                                           << 773             break;
1025                                               << 774     }
1026   // Search for discrete process              << 775   }
1027   if(nullptr == currentModel) {               << 
1028     G4VEmProcess* proc = FindDiscreteProcess( << 
1029     if(nullptr != proc) {                     << 
1030       currentModel = proc->SelectModelForMate << 
1031       currentModel->InitialiseForMaterial(par << 
1032       currentModel->SetupForMaterial(part, cu << 
1033       G4double eth = currentModel->LowEnergyL << 
1034       if(eth > 0.0) {                         << 
1035         loweModel = proc->SelectModelForMater << 
1036         if(loweModel == currentModel) { loweM << 
1037         else {                                << 
1038           loweModel->InitialiseForMaterial(pa << 
1039           loweModel->SetupForMaterial(part, c << 
1040         }                                     << 
1041       }                                       << 
1042     }                                         << 
1043   }                                           << 
1044                                               << 
1045   // Search for msc process                   << 
1046   if(nullptr == currentModel) {               << 
1047     G4VMultipleScattering* proc = FindMscProc << 
1048     if(nullptr != proc) {                     << 
1049       currentModel = proc->SelectModel(kinEne << 
1050       loweModel = nullptr;                    << 
1051     }                                         << 
1052   }                                           << 
1053   if(nullptr != currentModel) {               << 
1054     if(loweModel == currentModel) { loweModel << 
1055     isApplicable = true;                      << 
1056     currentModel->InitialiseForMaterial(part, << 
1057     if(loweModel) {                           << 
1058       loweModel->InitialiseForMaterial(part,  << 
1059     }                                         << 
1060     if(verbose > 1) {                         << 
1061       G4cout << "   Model <" << currentModel- << 
1062              << "> Emin(MeV)= " << currentMod << 
1063              << " for " << part->GetParticleN << 
1064       if(nullptr != elproc) {                 << 
1065         G4cout << " and " << elproc->GetProce << 
1066                << G4endl;                     << 
1067       }                                       << 
1068       if(nullptr != loweModel) {              << 
1069         G4cout << " LowEnergy model <" << low << 
1070       }                                       << 
1071       G4cout << G4endl;                       << 
1072     }                                         << 
1073   }                                           << 
1074   return isApplicable;                        << 
1075 }                                             << 
1076                                               << 
1077 //....oooOO0OOooo........oooOO0OOooo........o << 
1078                                               << 
1079 G4VEnergyLossProcess*                         << 
1080 G4EmCalculator::FindEnLossProcess(const G4Par << 
1081                                   const G4Str << 
1082 {                                             << 
1083   G4VEnergyLossProcess* proc = nullptr;       << 
1084   const std::vector<G4VEnergyLossProcess*> v  << 
1085     manager->GetEnergyLossProcessVector();    << 
1086   std::size_t n = v.size();                   << 
1087   for(std::size_t i=0; i<n; ++i) {            << 
1088     if((v[i])->GetProcessName() == processNam << 
1089       auto p = static_cast<G4VProcess*>(v[i]) << 
1090       if(ActiveForParticle(part, p)) {        << 
1091         proc = v[i];                          << 
1092         break;                                << 
1093       }                                          776       }
                                                   >> 777       //      G4cout << "model= " << currentModel  << G4endl;
                                                   >> 778       break;
1094     }                                            779     }
1095   }                                              780   }
1096   return proc;                                << 781   if(!currentModel) {
1097 }                                             << 782     const std::vector<G4VEmProcess*> vem = lManager->GetEmProcessVector();
1098                                               << 783     G4int n = vem.size();
1099 //....oooOO0OOooo........oooOO0OOooo........o << 784     for(G4int i=0; i<n; i++) {
1100                                               << 785       if((vem[i])->GetProcessName() == currentName && (vem[i])->Particle() == part) {
1101 G4VEmProcess*                                 << 786         currentModel = (vem[i])->SelectModelForMaterial(kinEnergy, idx);
1102 G4EmCalculator::FindDiscreteProcess(const G4P << 787   isApplicable    = true;
1103                                     const G4S << 
1104 {                                             << 
1105   G4VEmProcess* proc = nullptr;               << 
1106   auto v = manager->GetEmProcessVector();     << 
1107   std::size_t n = v.size();                   << 
1108   for(std::size_t i=0; i<n; ++i) {            << 
1109     const G4String& pName = v[i]->GetProcessN << 
1110     if(pName == "GammaGeneralProc") {         << 
1111       proc = v[i]->GetEmProcess(processName); << 
1112       break;                                  << 
1113     } else if(pName == processName) {         << 
1114       const auto p = static_cast<G4VProcess*> << 
1115       if(ActiveForParticle(part, p)) {        << 
1116         proc = v[i];                          << 
1117         break;                                   788         break;
1118       }                                          789       }
1119     }                                            790     }
1120   }                                              791   }
1121   return proc;                                << 792   if(!currentModel) {
1122 }                                             << 793     const std::vector<G4VMultipleScattering*> vmsc = lManager->GetMultipleScatteringVector();
1123                                               << 794     G4int n = vmsc.size();
1124 //....oooOO0OOooo........oooOO0OOooo........o << 795     for(G4int i=0; i<n; i++) {
1125                                               << 796       if((vmsc[i])->GetProcessName() == currentName && (vmsc[i])->Particle() == part) {
1126 G4VMultipleScattering*                        << 797         currentModel = (vmsc[i])->SelectModelForMaterial(kinEnergy, idx);
1127 G4EmCalculator::FindMscProcess(const G4Partic << 798   isApplicable    = true;
1128                                const G4String << 
1129 {                                             << 
1130   G4VMultipleScattering* proc = nullptr;      << 
1131   const std::vector<G4VMultipleScattering*> v << 
1132     manager->GetMultipleScatteringVector();   << 
1133   std::size_t n = v.size();                   << 
1134   for(std::size_t i=0; i<n; ++i) {            << 
1135     if((v[i])->GetProcessName() == processNam << 
1136       auto p = static_cast<G4VProcess*>(v[i]) << 
1137       if(ActiveForParticle(part, p)) {        << 
1138         proc = v[i];                          << 
1139         break;                                   799         break;
1140       }                                          800       }
1141     }                                            801     }
1142   }                                              802   }
1143   return proc;                                << 803   if(currentModel) res = true;
1144 }                                             << 
1145                                               << 
1146 //....oooOO0OOooo........oooOO0OOooo........o << 
1147                                               << 
1148 G4VProcess* G4EmCalculator::FindProcess(const << 
1149                                         const << 
1150 {                                             << 
1151   G4VProcess* proc = nullptr;                 << 
1152   const G4ProcessManager* procman = part->Get << 
1153   G4ProcessVector* pv = procman->GetProcessLi << 
1154   G4int nproc = (G4int)pv->size();            << 
1155   for(G4int i=0; i<nproc; ++i) {              << 
1156     if(processName == (*pv)[i]->GetProcessNam << 
1157       proc = (*pv)[i];                        << 
1158       break;                                  << 
1159     }                                         << 
1160   }                                           << 
1161   return proc;                                << 
1162 }                                             << 
1163                                               << 
1164 //....oooOO0OOooo........oooOO0OOooo........o << 
1165                                               << 
1166 G4bool G4EmCalculator::ActiveForParticle(cons << 
1167                                          G4VP << 
1168 {                                             << 
1169   G4ProcessManager* pm = part->GetProcessMana << 
1170   G4ProcessVector* pv = pm->GetProcessList(); << 
1171   G4int n = (G4int)pv->size();                << 
1172   G4bool res = false;                         << 
1173   for(G4int i=0; i<n; ++i) {                  << 
1174     if((*pv)[i] == proc) {                    << 
1175       if(pm->GetProcessActivation(i)) { res = << 
1176       break;                                  << 
1177     }                                         << 
1178   }                                           << 
1179   return res;                                    804   return res;
1180 }                                                805 }
1181                                                  806 
1182 //....oooOO0OOooo........oooOO0OOooo........o    807 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1183                                                  808 
1184 void G4EmCalculator::SetupMaterial(const G4Ma << 809 const G4VEnergyLossProcess* G4EmCalculator::FindEnergyLossProcess(
                                                   >> 810                       const G4ParticleDefinition* p)
1185 {                                                811 {
1186   if(mat) {                                   << 812   const G4VEnergyLossProcess* elp = 0;
1187     currentMaterial = mat;                    << 813   G4String partname =  p->GetParticleName();
1188     currentMaterialName = mat->GetName();     << 814   const G4ParticleDefinition* part = p;
1189   } else {                                    << 815   if(p->GetParticleType() == "nucleus" && partname != "deuteron" && partname != "triton")
1190     currentMaterial = nullptr;                << 816     part = G4GenericIon::GenericIon(); 
1191     currentMaterialName = "";                 << 817   
1192   }                                           << 818   G4LossTableManager* lManager = G4LossTableManager::Instance();
1193 }                                             << 819   const std::vector<G4VEnergyLossProcess*> vel = lManager->GetEnergyLossProcessVector();
1194                                               << 820   G4int n = vel.size();
1195 //....oooOO0OOooo........oooOO0OOooo........o << 821   for(G4int i=0; i<n; i++) {
1196                                               << 822     if((vel[i])->Particle() == part) {
1197 void G4EmCalculator::SetupMaterial(const G4St << 823       elp = vel[i];
1198 {                                             << 824       break;
1199   SetupMaterial(nist->FindOrBuildMaterial(mna << 
1200 }                                             << 
1201                                               << 
1202 //....oooOO0OOooo........oooOO0OOooo........o << 
1203                                               << 
1204 void G4EmCalculator::CheckMaterial(G4int Z)   << 
1205 {                                             << 
1206   G4bool isFound = false;                     << 
1207   if(nullptr != currentMaterial) {            << 
1208     G4int nn = (G4int)currentMaterial->GetNum << 
1209     for(G4int i=0; i<nn; ++i) {               << 
1210       if(Z == currentMaterial->GetElement(i)- << 
1211         isFound = true;                       << 
1212         break;                                << 
1213       }                                       << 
1214     }                                            825     }
1215   }                                              826   }
1216   if(!isFound) {                              << 827   return elp;
1217     SetupMaterial(nist->FindOrBuildSimpleMate << 
1218   }                                           << 
1219 }                                                828 }
1220                                                  829 
1221 //....oooOO0OOooo........oooOO0OOooo........o    830 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1222                                                  831 
1223 void G4EmCalculator::SetVerbose(G4int verb)      832 void G4EmCalculator::SetVerbose(G4int verb)
1224 {                                                833 {
1225   verbose = verb;                                834   verbose = verb;
1226 }                                                835 }
1227                                                  836 
1228 //....oooOO0OOooo........oooOO0OOooo........o    837 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1229                                                  838 
1230                                                  839