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 7.0)


  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.11 2004/11/17 10:11:25 vnivanch Exp $
                                                   >>  24 // GEANT4 tag $Name: geant4-07-00-cand-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)
 38 //                                                 40 //
 39 // Class Description: V.Ivanchenko & M.Novak   <<  41 // Class Description:
 40 //                                                 42 //
 41 // -------------------------------------------     43 // -------------------------------------------------------------------
 42 //                                                 44 //
 43 //....oooOO0OOooo........oooOO0OOooo........oo     45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 44 //....oooOO0OOooo........oooOO0OOooo........oo     46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 45                                                    47 
 46 #include "G4EmCalculator.hh"                       48 #include "G4EmCalculator.hh"
 47 #include "G4SystemOfUnits.hh"                  << 
 48 #include "G4LossTableManager.hh"                   49 #include "G4LossTableManager.hh"
 49 #include "G4EmParameters.hh"                   << 
 50 #include "G4NistManager.hh"                    << 
 51 #include "G4DynamicParticle.hh"                << 
 52 #include "G4VEmProcess.hh"                         50 #include "G4VEmProcess.hh"
 53 #include "G4VEnergyLossProcess.hh"                 51 #include "G4VEnergyLossProcess.hh"
 54 #include "G4VMultipleScattering.hh"                52 #include "G4VMultipleScattering.hh"
 55 #include "G4Material.hh"                           53 #include "G4Material.hh"
 56 #include "G4MaterialCutsCouple.hh"                 54 #include "G4MaterialCutsCouple.hh"
 57 #include "G4ParticleDefinition.hh"                 55 #include "G4ParticleDefinition.hh"
 58 #include "G4ParticleTable.hh"                      56 #include "G4ParticleTable.hh"
 59 #include "G4IonTable.hh"                       << 
 60 #include "G4PhysicsTable.hh"                       57 #include "G4PhysicsTable.hh"
 61 #include "G4ProductionCutsTable.hh"                58 #include "G4ProductionCutsTable.hh"
 62 #include "G4ProcessManager.hh"                     59 #include "G4ProcessManager.hh"
 63 #include "G4ionEffectiveCharge.hh"                 60 #include "G4ionEffectiveCharge.hh"
 64 #include "G4RegionStore.hh"                        61 #include "G4RegionStore.hh"
 65 #include "G4Element.hh"                            62 #include "G4Element.hh"
 66 #include "G4EmCorrections.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                                                    63 
 74 //....oooOO0OOooo........oooOO0OOooo........oo     64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 75                                                    65 
 76 G4EmCalculator::G4EmCalculator()                   66 G4EmCalculator::G4EmCalculator()
 77 {                                                  67 {
 78   manager = G4LossTableManager::Instance();        68   manager = G4LossTableManager::Instance();
 79   nist    = G4NistManager::Instance();         <<  69   nLocalMaterials    = 0;
 80   theParameters = G4EmParameters::Instance();  <<  70   verbose            = 0;
 81   corr    = manager->EmCorrections();          <<  71   currentCoupleIndex = 0;
 82   cutenergy[0] = cutenergy[1] = cutenergy[2] = <<  72   currentCouple      = 0;
 83   theGenericIon = G4GenericIon::GenericIon();  <<  73   currentMaterial    = 0;
 84   ionEffCharge  = new G4ionEffectiveCharge();  <<  74   currentParticle    = 0;
 85   dynParticle   = new G4DynamicParticle();     <<  75   currentLambda      = 0;
 86   ionTable      = G4ParticleTable::GetParticle <<  76   chargeSquare       = 1.0;
                                                   >>  77   massRatio          = 1.0;
                                                   >>  78   currentParticleName= "";
                                                   >>  79   currentMaterialName= "";
                                                   >>  80   ionEffCharge       = new G4ionEffectiveCharge();
                                                   >>  81   isIon              = false;
                                                   >>  82   isApplicable       = false;
 87 }                                                  83 }
 88                                                    84 
 89 //....oooOO0OOooo........oooOO0OOooo........oo     85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 90                                                    86 
 91 G4EmCalculator::~G4EmCalculator()                  87 G4EmCalculator::~G4EmCalculator()
 92 {                                                  88 {
 93   delete ionEffCharge;                             89   delete ionEffCharge;
 94   delete dynParticle;                          <<  90   for (G4int i=0; i<nLocalMaterials; i++) {
 95   for (G4int i=0; i<nLocalMaterials; ++i) {    << 
 96     delete localCouples[i];                        91     delete localCouples[i];
 97   }                                                92   }
 98 }                                                  93 }
 99                                                    94 
100 //....oooOO0OOooo........oooOO0OOooo........oo     95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101                                                    96 
102 G4double G4EmCalculator::GetDEDX(G4double kinE <<  97 G4double G4EmCalculator::GetDEDX(G4double kinEnergy, const G4ParticleDefinition* p,
103                                  const G4Parti <<  98                                  const G4Material* mat, const G4Region* region)
104                                  const G4Mater << 
105                                  const G4Regio << 
106 {                                                  99 {
107   G4double res = 0.0;                             100   G4double res = 0.0;
108   const G4MaterialCutsCouple* couple = FindCou    101   const G4MaterialCutsCouple* couple = FindCouple(mat, region);
109   if(nullptr != couple && UpdateParticle(p, ki << 102   if(couple && UpdateParticle(p, kinEnergy) ) {
110     res = manager->GetDEDX(p, kinEnergy, coupl    103     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) {                               104     if(verbose>0) {
128       G4cout << "G4EmCalculator::GetDEDX: E(Me << 105       G4cout << "E(MeV)= " << kinEnergy/MeV
129              << " DEDX(MeV/mm)= " << res*mm/Me << 106        << " DEDX(MeV/mm)= " << res*mm/MeV
130              << " DEDX(MeV*cm^2/g)= " << res*g << 107        << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
131              << "  " <<  p->GetParticleName()  << 108        << "  " <<  p->GetParticleName()
132              << " in " <<  mat->GetName()      << 109        << " in " <<  mat->GetName()
133              << " isIon= " << isIon            << 110        << G4endl;
134              << G4endl;                        << 
135     }                                             111     }
136   }                                               112   }
137   return res;                                     113   return res;
138 }                                                 114 }
139                                                   115 
140 //....oooOO0OOooo........oooOO0OOooo........oo    116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
141                                                   117 
142 G4double G4EmCalculator::GetRangeFromRestricte << 118 G4double G4EmCalculator::GetDEDX(G4double kinEnergy, const G4String& particle,
143                                                << 119                                  const G4String& material, const G4String& reg)
144                                                << 
145                                                << 
146 {                                                 120 {
147   G4double res = 0.0;                          << 121   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 }                                                 122 }
162                                                   123 
163 //....oooOO0OOooo........oooOO0OOooo........oo    124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
164                                                   125 
165 G4double G4EmCalculator::GetCSDARange(G4double << 126 G4double G4EmCalculator::GetRange(G4double kinEnergy, const G4ParticleDefinition* p,
166                                       const G4 << 127                                   const G4Material* mat, const G4Region* region)
167                                       const G4 << 
168                                       const G4 << 
169 {                                                 128 {
170   G4double res = 0.0;                             129   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    130   const G4MaterialCutsCouple* couple = FindCouple(mat,region);
181   if(nullptr != couple && UpdateParticle(p, ki << 131   if(couple && UpdateParticle(p, kinEnergy)) {
182     res = manager->GetCSDARange(p, kinEnergy,  << 132     res = manager->GetRange(p, kinEnergy, couple);
183     if(verbose>1) {                            << 133     if(verbose>0) {
184       G4cout << " G4EmCalculator::GetCSDARange << 134       G4cout << "E(MeV)= " << kinEnergy/MeV
185              << " range(mm)= " << res/mm       << 135        << " range(mm)= " << res/mm
186              << "  " <<  p->GetParticleName()  << 136        << "  " <<  p->GetParticleName()
187              << " in " <<  mat->GetName()      << 137        << " in " <<  mat->GetName()
188              << G4endl;                        << 138        << G4endl;
189     }                                             139     }
190   }                                               140   }
191   return res;                                     141   return res;
192 }                                                 142 }
193                                                   143 
194 //....oooOO0OOooo........oooOO0OOooo........oo    144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
195                                                   145 
196 G4double G4EmCalculator::GetRange(G4double kin << 146 G4double G4EmCalculator::GetRange(G4double kinEnergy, const G4String& particle,
197                                   const G4Part << 147                                   const G4String& material, const G4String& reg)
198                                   const G4Mate << 
199                                   const G4Regi << 
200 {                                                 148 {
201   G4double res = 0.0;                          << 149   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 }                                                 150 }
209                                                   151 
210 //....oooOO0OOooo........oooOO0OOooo........oo    152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
211                                                   153 
212 G4double G4EmCalculator::GetKinEnergy(G4double << 154 G4double G4EmCalculator::GetKinEnergy(G4double range, const G4ParticleDefinition* p,
213                                       const G4 << 155                                       const G4Material* mat, const G4Region* region)
214                                       const G4 << 
215                                       const G4 << 
216 {                                                 156 {
217   G4double res = 0.0;                             157   G4double res = 0.0;
218   const G4MaterialCutsCouple* couple = FindCou    158   const G4MaterialCutsCouple* couple = FindCouple(mat,region);
219   if(nullptr != couple && UpdateParticle(p, 1. << 159   if(couple && UpdateParticle(p, 1.0*GeV)) {
220     res = manager->GetEnergy(p, range, couple)    160     res = manager->GetEnergy(p, range, couple);
221     if(verbose>0) {                               161     if(verbose>0) {
222       G4cout << "G4EmCalculator::GetKinEnergy: << 162       G4cout << "Range(mm)= " << range/mm
223              << " KinE(MeV)= " << res/MeV      << 163        << " KinE(MeV)= " << res/MeV
224              << "  " <<  p->GetParticleName()  << 164        << "  " <<  p->GetParticleName()
225              << " in " <<  mat->GetName()      << 165        << " in " <<  mat->GetName()
226              << G4endl;                        << 166        << G4endl;
227     }                                             167     }
228   }                                               168   }
229   return res;                                     169   return res;
230 }                                                 170 }
231                                                   171 
232 //....oooOO0OOooo........oooOO0OOooo........oo    172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
233                                                   173 
                                                   >> 174 G4double G4EmCalculator::GetKinEnergy(G4double range, const G4String& particle,
                                                   >> 175                                       const G4String& material, const G4String& reg)
                                                   >> 176 {
                                                   >> 177   return GetKinEnergy(range,FindParticle(particle),FindMaterial(material),FindRegion(reg));
                                                   >> 178 }
                                                   >> 179 
                                                   >> 180 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >> 181 
234 G4double G4EmCalculator::GetCrossSectionPerVol    182 G4double G4EmCalculator::GetCrossSectionPerVolume(G4double kinEnergy,
235                                             co    183                                             const G4ParticleDefinition* p,
236                                             co    184                                             const G4String& processName,
237                                             co << 185               const G4Material* mat,
238                                             co << 186               const G4Region* region)
239 {                                                 187 {
240   G4double res = 0.0;                             188   G4double res = 0.0;
241   const G4MaterialCutsCouple* couple = FindCou    189   const G4MaterialCutsCouple* couple = FindCouple(mat,region);
242                                                   190 
243   if(nullptr != couple && UpdateParticle(p, ki << 191   if(couple) {
244     if(FindEmModel(p, processName, kinEnergy)) << 192     G4int idx = couple->GetIndex();
245       G4int idx      = couple->GetIndex();     << 193     FindLambdaTable(p, processName);
246       G4int procType = -1;                     << 194     if(currentLambda && UpdateParticle(p, kinEnergy)) {
247       FindLambdaTable(p, processName, kinEnerg << 195       G4bool b;
248                                                << 196       G4double e = kinEnergy*massRatio;
249       G4VEmProcess* emproc = FindDiscreteProce << 197       res = (((*currentLambda)[idx])->GetValue(e,b))*chargeSquare;
250       if(nullptr != emproc) {                  << 198       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) {                             199       if(verbose>0) {
270         G4cout << "G4EmCalculator::GetXSPerVol << 200   G4cout << "E(MeV)= " << kinEnergy/MeV
271                << " cross(cm-1)= " << res*cm   << 201          << " cross(cm^2/g)= " << res*gram/cm2
272                << "  " <<  p->GetParticleName( << 202          << "  " <<  p->GetParticleName()
273                << " in " <<  mat->GetName();   << 203          << " in " <<  mat->GetName()
274         if(verbose>1)                          << 204          << G4endl;
275           G4cout << "  idx= " << idx << "  Esc << 205       }
276            << kinEnergy*massRatio              << 
277            << "  q2= " << chargeSquare;        << 
278         G4cout << G4endl;                      << 
279       }                                        << 
280     }                                             206     }
281   }                                               207   }
282   return res;                                     208   return res;
283 }                                                 209 }
284                                                   210 
285 //....oooOO0OOooo........oooOO0OOooo........oo    211 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
286                                                   212 
287 G4double G4EmCalculator::GetShellIonisationCro << 213 G4double G4EmCalculator::GetCrossSectionPerVolume(G4double kinEnergy,
288                                          const << 214                                             const G4String& particle,
289                                          G4int << 215               const G4String& processName,
290                                          G4Ato << 216                                             const G4String& material,
291                                          G4dou << 217               const G4String& reg)
292 {                                                 218 {
293   G4double res = 0.0;                          << 219   return GetCrossSectionPerVolume(kinEnergy,FindParticle(particle),processName,
294   const G4ParticleDefinition* p = FindParticle << 220                                   FindMaterial(material),FindRegion(reg));
295   G4VAtomDeexcitation* ad = manager->AtomDeexc << 221 }
296   if(nullptr != p && nullptr != ad) {          << 222 
297     res = ad->GetShellIonisationCrossSectionPe << 223 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
298   }                                            << 224 
                                                   >> 225 G4double G4EmCalculator::GetCrossSectionPerAtom(G4double kinEnergy,
                                                   >> 226                                           const G4ParticleDefinition* p,
                                                   >> 227                                           const G4String& processName,
                                                   >> 228             const G4Material* mat,
                                                   >> 229             const G4Region* region)
                                                   >> 230 {
                                                   >> 231   G4double res = GetCrossSectionPerVolume(kinEnergy,p,processName,mat,region);
                                                   >> 232   if(mat) res /= mat->GetTotNbOfAtomsPerVolume();
299   return res;                                     233   return res;
300 }                                                 234 }
301                                                   235 
302 //....oooOO0OOooo........oooOO0OOooo........oo    236 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
303                                                   237 
                                                   >> 238 G4double G4EmCalculator::GetCrossSectionPerAtom(  G4double kinEnergy,
                                                   >> 239                                             const G4String& particle,
                                                   >> 240               const G4String& processName,
                                                   >> 241                                             const G4String& material,
                                                   >> 242               const G4String& reg)
                                                   >> 243 {
                                                   >> 244   return GetCrossSectionPerAtom(kinEnergy,FindParticle(particle),processName,
                                                   >> 245                                 FindMaterial(material),FindRegion(reg));
                                                   >> 246 }
                                                   >> 247 
                                                   >> 248 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >> 249 
304 G4double G4EmCalculator::GetMeanFreePath(G4dou    250 G4double G4EmCalculator::GetMeanFreePath(G4double kinEnergy,
305                                          const    251                                          const G4ParticleDefinition* p,
306                                          const    252                                          const G4String& processName,
307                                          const << 253            const G4Material* mat,
308                                          const    254                                          const G4Region* region)
309 {                                                 255 {
310   G4double res = DBL_MAX;                         256   G4double res = DBL_MAX;
311   G4double x = GetCrossSectionPerVolume(kinEne    257   G4double x = GetCrossSectionPerVolume(kinEnergy,p, processName, mat,region);
312   if(x > 0.0) { res = 1.0/x; }                 << 258   if(x > 0.0) res = 1.0/(x*(mat->GetTotNbOfAtomsPerVolume()));
313   if(verbose>1) {                                 259   if(verbose>1) {
314     G4cout << "G4EmCalculator::GetMeanFreePath << 260     G4cout << "E(MeV)= " << kinEnergy/MeV
315            << " MFP(mm)= " << res/mm           << 261      << " MFP(mm)= " << res/mm
316            << "  " <<  p->GetParticleName()    << 262      << "  " <<  p->GetParticleName()
317            << " in " <<  mat->GetName()        << 263      << " in " <<  mat->GetName()
318            << G4endl;                          << 264      << G4endl;
319   }                                               265   }
320   return res;                                     266   return res;
321 }                                                 267 }
322                                                   268 
323 //....oooOO0OOooo........oooOO0OOooo........oo    269 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
324                                                   270 
                                                   >> 271 G4double G4EmCalculator::GetMeanFreePath(G4double kinEnergy,
                                                   >> 272                                          const G4String& particle,
                                                   >> 273            const G4String& processName,
                                                   >> 274                                          const G4String& material,
                                                   >> 275            const G4String& reg)
                                                   >> 276 {
                                                   >> 277   return GetMeanFreePath(kinEnergy,FindParticle(particle),processName,
                                                   >> 278                          FindMaterial(material),FindRegion(reg));
                                                   >> 279 }
                                                   >> 280 
                                                   >> 281 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >> 282 
325 void G4EmCalculator::PrintDEDXTable(const G4Pa    283 void G4EmCalculator::PrintDEDXTable(const G4ParticleDefinition* p)
326 {                                                 284 {
327   const G4VEnergyLossProcess* elp = manager->G << 285   const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p);
328   G4cout << "##### DEDX Table for " << p->GetP    286   G4cout << "##### DEDX Table for " << p->GetParticleName() << G4endl;
329   if(nullptr != elp) G4cout << *(elp->DEDXTabl << 287   if(elp) G4cout << *(elp->DEDXTable()) << G4endl;
330 }                                                 288 }
331                                                   289 
332 //....oooOO0OOooo........oooOO0OOooo........oo    290 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
333                                                   291 
334 void G4EmCalculator::PrintRangeTable(const G4P    292 void G4EmCalculator::PrintRangeTable(const G4ParticleDefinition* p)
335 {                                                 293 {
336   const G4VEnergyLossProcess* elp = manager->G << 294   const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p);
337   G4cout << "##### Range Table for " << p->Get    295   G4cout << "##### Range Table for " << p->GetParticleName() << G4endl;
338   if(nullptr != elp) G4cout << *(elp->RangeTab << 296   if(elp) G4cout << *(elp->RangeTableForLoss()) << G4endl;
339 }                                                 297 }
340                                                   298 
341 //....oooOO0OOooo........oooOO0OOooo........oo    299 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
342                                                   300 
343 void G4EmCalculator::PrintInverseRangeTable(co    301 void G4EmCalculator::PrintInverseRangeTable(const G4ParticleDefinition* p)
344 {                                                 302 {
345   const G4VEnergyLossProcess* elp = manager->G << 303   const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p);
346   G4cout << "### G4EmCalculator: Inverse Range << 304   G4cout << "### G4EmCalculator: Inverse Range Table for " << p->GetParticleName() << G4endl;
347          << p->GetParticleName() << G4endl;    << 305   if(elp) G4cout << *(elp->InverseRangeTable()) << G4endl;
348   if(nullptr != elp) G4cout << *(elp->InverseR << 
349 }                                                 306 }
350                                                   307 
351 //....oooOO0OOooo........oooOO0OOooo........oo    308 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
352                                                   309 
353 G4double G4EmCalculator::ComputeDEDX(G4double     310 G4double G4EmCalculator::ComputeDEDX(G4double kinEnergy,
354                                      const G4P    311                                      const G4ParticleDefinition* p,
355                                      const G4S    312                                      const G4String& processName,
356                                      const G4M << 313              const G4Material* mat,
357                                            G4d    314                                            G4double cut)
358 {                                                 315 {
359   SetupMaterial(mat);                          << 
360   G4double res = 0.0;                             316   G4double res = 0.0;
361   if(verbose > 1) {                               317   if(verbose > 1) {
362     G4cout << "### G4EmCalculator::ComputeDEDX << 318     G4cout << "ComputeDEDX: " << p->GetParticleName()
363            << " in " << currentMaterialName    << 319            << " in " << mat->GetName()
364            << " e(MeV)= " << kinEnergy/MeV <<  << 320            << " e(MeV)= " << kinEnergy/MeV << "  cut(MeV)= " << cut/MeV << G4endl;
365            << G4endl;                          << 321   }
366   }                                            << 322   if(FindEmModel(p, processName, kinEnergy)) {
367   if(UpdateParticle(p, kinEnergy)) {           << 323     //    G4cout << "currentModel= " << currentModel << G4endl;
368     if(FindEmModel(p, processName, kinEnergy)) << 324     if(UpdateCouple(mat, cut) && UpdateParticle(p, kinEnergy)) {
369       G4double escaled = kinEnergy*massRatio;  << 325       if(baseParticle) {
370       if(nullptr != baseParticle) {            << 326         G4double e = kinEnergy*massRatio;
371   res = currentModel->ComputeDEDXPerVolume(mat << 327         res = currentModel->ComputeDEDX(currentCouple, baseParticle, e, cut);
372                                                << 328         if(verbose > 1) {
373   if(verbose > 1) {                            << 329           G4cout <<  baseParticle->GetParticleName() << " E(MeV)= " << e
374     G4cout << "Particle: " << p->GetParticleNa << 330                  << " DEDX(MeV/mm)= " << res*mm/MeV
375      << " E(MeV)=" << kinEnergy                << 331                  << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
376      << " Base particle: " << baseParticle->Ge << 332                  << G4endl;
377      << " Escaled(MeV)= " << escaled           << 
378      << " q2=" << chargeSquare << G4endl;      << 
379   }                                            << 
380       } else {                                 << 
381   res = currentModel->ComputeDEDXPerVolume(mat << 
382   if(verbose > 1) {                            << 
383     G4cout << "Particle: " << p->GetParticleNa << 
384      << " E(MeV)=" << kinEnergy << G4endl;     << 
385   }                                            << 
386       }                                        << 
387       if(verbose > 1) {                        << 
388   G4cout << currentModel->GetName() << ": DEDX << 
389          << " DEDX(MeV*cm^2/g)= "              << 
390          << res*gram/(MeV*cm2*mat->GetDensity( << 
391          << G4endl;                            << 
392       }                                        << 
393       // emulate smoothing procedure           << 
394       if(applySmoothing && nullptr != loweMode << 
395   G4double eth = currentModel->LowEnergyLimit( << 
396   G4double res0 = 0.0;                         << 
397   G4double res1 = 0.0;                         << 
398   if(nullptr != baseParticle) {                << 
399     res1 = chargeSquare*                       << 
400       currentModel->ComputeDEDXPerVolume(mat,  << 
401     res0 = chargeSquare*                       << 
402       loweModel->ComputeDEDXPerVolume(mat, bas << 
403   } else {                                     << 
404     res1 = currentModel->ComputeDEDXPerVolume( << 
405     res0 = loweModel->ComputeDEDXPerVolume(mat << 
406   }                                            << 
407   if(res1 > 0.0 && escaled > 0.0) {            << 
408     res *= (1.0 + (res0/res1 - 1.0)*eth/escale << 
409   }                                            << 
410   if(verbose > 1) {                            << 
411     G4cout << "At boundary energy(MeV)= " << e << 
412      << " DEDX(MeV/mm)= " << res0*mm/MeV << "  << 
413      << " after correction DEDX(MeV/mm)=" << r << 
414   }                                            << 
415       }                                        << 
416       // correction for ions                   << 
417       if(isIon) {                              << 
418   const G4double length = CLHEP::nm;           << 
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   }                                               333   }
                                                   >> 334         res *= chargeSquare;
                                                   >> 335       } else {
                                                   >> 336         res = currentModel->ComputeDEDX(currentCouple, p, kinEnergy, cut);
432       }                                           337       }
433       if(verbose > 0) {                           338       if(verbose > 0) {
434   G4cout << "## E(MeV)= " << kinEnergy/MeV     << 339         G4cout << "E(MeV)= " << kinEnergy/MeV
435          << " DEDX(MeV/mm)= " << res*mm/MeV    << 340                << " DEDX(MeV/mm)= " << res*mm/MeV
436          << " DEDX(MeV*cm^2/g)= " << res*gram/ << 341                << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
437          << " cut(MeV)= " << cut/MeV           << 342                << " cut(MeV)= " << cut/MeV
438          << "  " <<  p->GetParticleName()      << 343                << "  " <<  p->GetParticleName()
439          << " in " <<  currentMaterialName     << 344                << " in " <<  mat->GetName()
440          << " Zi^2= " << chargeSquare          << 345                << " Zi^2= " << chargeSquare
441          << " isIon=" << isIon                 << 346                << G4endl;
442          << G4endl;                            << 
443       }                                           347       }
444     }                                             348     }
445   }                                               349   }
446   return res;                                     350   return res;
447 }                                                 351 }
448                                                   352 
449 //....oooOO0OOooo........oooOO0OOooo........oo    353 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
450                                                   354 
451 G4double G4EmCalculator::ComputeElectronicDEDX << 355 G4double G4EmCalculator::ComputeDEDX(G4double kinEnergy,
452                                                << 356                                      const G4String& particle,
453                                                << 357              const G4String& processName,
454                                                << 358                                      const G4String& material,
455 {                                              << 359                                            G4double cut)
456   SetupMaterial(mat);                          << 
457   G4double dedx = 0.0;                         << 
458   if(UpdateParticle(part, kinEnergy)) {        << 
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 }                                              << 
481                                                << 
482 //....oooOO0OOooo........oooOO0OOooo........oo << 
483                                                << 
484 G4double                                       << 
485 G4EmCalculator::ComputeDEDXForCutInRange(G4dou << 
486                                          const << 
487                                          const << 
488                                          G4dou << 
489 {                                              << 
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                                                << 
509     //G4cout << "ComputeElectronicDEDX for " < << 
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                                                << 
523           dedx += ComputeDEDX(kinEnergy,part,( << 
524                               mat,cutenergy[id << 
525         }                                      << 
526       }                                        << 
527     }                                          << 
528   }                                            << 
529   return dedx;                                 << 
530 }                                              << 
531                                                << 
532 //....oooOO0OOooo........oooOO0OOooo........oo << 
533                                                << 
534 G4double G4EmCalculator::ComputeTotalDEDX(G4do << 
535                                           cons << 
536                                           cons << 
537                                           G4do << 
538 {                                                 360 {
539   G4double dedx = ComputeElectronicDEDX(kinEne << 361   return ComputeDEDX(kinEnergy,FindParticle(particle),processName,
540   if(mass > 700.*MeV) { dedx += ComputeNuclear << 362                      FindMaterial(material),cut);
541   return dedx;                                 << 
542 }                                              << 
543                                                << 
544 //....oooOO0OOooo........oooOO0OOooo........oo << 
545                                                << 
546 G4double G4EmCalculator::ComputeNuclearDEDX(G4 << 
547                                       const G4 << 
548                                       const G4 << 
549 {                                              << 
550   G4double res = 0.0;                          << 
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 }                                                 363 }
569                                                   364 
570 //....oooOO0OOooo........oooOO0OOooo........oo    365 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
571                                                   366 
572 G4double G4EmCalculator::ComputeCrossSectionPe    367 G4double G4EmCalculator::ComputeCrossSectionPerVolume(
573                                                   368                                                    G4double kinEnergy,
574                                              c    369                                              const G4ParticleDefinition* p,
575                                              c    370                                              const G4String& processName,
576                                              c << 371                const G4Material* mat,
577                                                   372                                                    G4double cut)
578 {                                                 373 {
579   SetupMaterial(mat);                          << 
580   G4double res = 0.0;                             374   G4double res = 0.0;
581   if(UpdateParticle(p, kinEnergy)) {           << 375   if(FindEmModel(p, processName, kinEnergy)) {
582     if(FindEmModel(p, processName, kinEnergy)) << 376     if(UpdateCouple(mat, cut) && UpdateParticle(p, kinEnergy)) {
583       G4double e = kinEnergy;                     377       G4double e = kinEnergy;
584       G4double aCut = std::max(cut, theParamet << 
585       if(baseParticle) {                          378       if(baseParticle) {
586         e *= kinEnergy*massRatio;                 379         e *= kinEnergy*massRatio;
587         res = currentModel->CrossSectionPerVol << 380         res = currentModel->CrossSection(currentCouple, baseParticle, e, cut, e)*chargeSquare;
588               mat, baseParticle, e, aCut, e) * << 
589       } else {                                    381       } else {
590         res = currentModel->CrossSectionPerVol << 382         res = currentModel->CrossSection(currentCouple, p, e, cut, e);
591       }                                        << 
592       if(verbose>0) {                          << 
593         G4cout << "G4EmCalculator::ComputeXSPe << 
594                << kinEnergy/MeV                << 
595                << " cross(cm-1)= " << res*cm   << 
596                << " cut(keV)= " << aCut/keV    << 
597                << "  " <<  p->GetParticleName( << 
598                << " in " <<  mat->GetName()    << 
599                << G4endl;                      << 
600       }                                           383       }
601     }                                             384     }
                                                   >> 385     if(verbose>0) {
                                                   >> 386       G4cout << "E(MeV)= " << kinEnergy/MeV
                                                   >> 387        << " cross(cm^2/g)= " << res*gram/cm2
                                                   >> 388        << "  " <<  p->GetParticleName()
                                                   >> 389        << " in " <<  mat->GetName()
                                                   >> 390        << G4endl;
                                                   >> 391     }
602   }                                               392   }
603   return res;                                     393   return res;
604 }                                                 394 }
605                                                   395 
606 //....oooOO0OOooo........oooOO0OOooo........oo    396 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
607                                                   397 
608 G4double                                       << 398 G4double G4EmCalculator::ComputeCrossSectionPerVolume(
609 G4EmCalculator::ComputeCrossSectionPerAtom(G4d << 399                                                    G4double kinEnergy,
610                                            con << 400                                              const G4String& particle,
611                                            con << 401                const G4String& processName,
612                                            G4d << 402                                              const G4String& material,
613                                            G4d << 403                                                    G4double cut)
614 {                                                 404 {
615   G4double res = 0.0;                          << 405   return ComputeCrossSectionPerVolume(kinEnergy,FindParticle(particle),processName,
616   if(UpdateParticle(p, kinEnergy)) {           << 406                                       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 }                                                 407 }
643                                                   408 
644 //....oooOO0OOooo........oooOO0OOooo........oo    409 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
645                                                   410 
646 G4double                                       << 411 G4double G4EmCalculator::ComputeCrossSectionPerAtom(
647 G4EmCalculator::ComputeCrossSectionPerShell(G4 << 412                                                    G4double kinEnergy,
648                                             co << 413                const G4ParticleDefinition* p,
649                                             co << 414                                              const G4String& processName,
650                                             G4 << 415                      G4double Z, G4double A,
651                                             G4 << 416                                        G4double cut)
652 {                                                 417 {
653   G4double res = 0.0;                             418   G4double res = 0.0;
654   if(UpdateParticle(p, kinEnergy)) {           << 419   if(FindEmModel(p, processName, kinEnergy)) {
655     CheckMaterial(Z);                          << 420     if(UpdateParticle(p, kinEnergy)) {
656     if(FindEmModel(p, processName, kinEnergy)) << 
657       G4double e = kinEnergy;                     421       G4double e = kinEnergy;
658       G4double aCut = std::max(cut, theParamet << 422       if(baseParticle) {
659       if(nullptr != baseParticle) {            << 
660         e *= kinEnergy*massRatio;                 423         e *= kinEnergy*massRatio;
661         currentModel->InitialiseForElement(bas << 424         res = currentModel->ComputeCrossSectionPerAtom(baseParticle, e, Z, A, cut)*chargeSquare;
662         res =                                  << 
663           currentModel->ComputeCrossSectionPer << 
664                                                << 
665       } else {                                    425       } else {
666         currentModel->InitialiseForElement(p,  << 426         res = currentModel->ComputeCrossSectionPerAtom(p, e, Z, A, cut);
667         res = currentModel->ComputeCrossSectio << 
668       }                                        << 
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       }                                           427       }
677     }                                             428     }
                                                   >> 429     if(verbose>0) {
                                                   >> 430       G4cout << "E(MeV)= " << kinEnergy/MeV
                                                   >> 431        << " cross(barn)= " << res/barn
                                                   >> 432        << "  " <<  p->GetParticleName()
                                                   >> 433        << " Z= " <<  Z << " A= " << A
                                                   >> 434        << G4endl;
                                                   >> 435     }
678   }                                               436   }
679   return res;                                     437   return res;
680 }                                                 438 }
681                                                   439 
682 //....oooOO0OOooo........oooOO0OOooo........oo    440 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
683                                                   441 
684 G4double                                       << 442 G4double G4EmCalculator::ComputeCrossSectionPerAtom(G4double kinEnergy,
685 G4EmCalculator::ComputeGammaAttenuationLength( << 443                                               const G4String& particle,
686                                                << 444                                               const G4String& processName,
687 {                                              << 445                 const G4Element* elm,
688   G4double res = 0.0;                          << 446                                         G4double cut)
689   const G4ParticleDefinition* gamma = G4Gamma: << 
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 {                                                 447 {
707   G4double res = 0.0;                          << 448   return ComputeCrossSectionPerAtom(kinEnergy,FindParticle(particle),processName,
708   const G4ParticleDefinition* p = FindParticle << 449                                     elm->GetZ(),elm->GetA(),cut);
709   G4VAtomDeexcitation* ad = manager->AtomDeexc << 
710   if(p && ad) {                                << 
711     res = ad->ComputeShellIonisationCrossSecti << 
712                                                << 
713   }                                            << 
714   return res;                                  << 
715 }                                                 450 }
716                                                   451 
717 //....oooOO0OOooo........oooOO0OOooo........oo    452 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
718                                                   453 
719 G4double G4EmCalculator::ComputeMeanFreePath(G    454 G4double G4EmCalculator::ComputeMeanFreePath(G4double kinEnergy,
720                                              c    455                                              const G4ParticleDefinition* p,
721                                              c    456                                              const G4String& processName,
722                                              c << 457                const G4Material* mat,
723                                              G << 458                                                    G4double cut)
724 {                                                 459 {
725   G4double mfp = DBL_MAX;                         460   G4double mfp = DBL_MAX;
726   G4double x =                                 << 461   G4double x = ComputeCrossSectionPerVolume(kinEnergy, p, processName, mat, cut);
727     ComputeCrossSectionPerVolume(kinEnergy, p, << 462   if(x > 0.0) mfp = 1.0/(x*(currentMaterial->GetTotNbOfAtomsPerVolume()));
728   if(x > 0.0) { mfp = 1.0/x; }                 << 
729   if(verbose>1) {                                 463   if(verbose>1) {
730     G4cout << "E(MeV)= " << kinEnergy/MeV         464     G4cout << "E(MeV)= " << kinEnergy/MeV
731            << " MFP(mm)= " << mfp/mm           << 465      << " MFP(mm)= " << mfp/mm
732            << "  " <<  p->GetParticleName()    << 466      << "  " <<  p->GetParticleName()
733            << " in " <<  mat->GetName()        << 467      << " in " <<  mat->GetName()
734            << G4endl;                          << 468      << G4endl;
735   }                                               469   }
736   return mfp;                                     470   return mfp;
737 }                                                 471 }
738                                                   472 
739 //....oooOO0OOooo........oooOO0OOooo........oo    473 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
740                                                   474 
741 G4double G4EmCalculator::ComputeEnergyCutFromR << 475 G4double G4EmCalculator::ComputeMeanFreePath(G4double kinEnergy,
742                          G4double range,       << 476                                              const G4String& particle,
743                          const G4ParticleDefin << 477                                              const G4String& processName,
744                          const G4Material* mat << 478                                              const G4String& material,
                                                   >> 479                                                    G4double cut)
745 {                                                 480 {
746   return G4ProductionCutsTable::GetProductionC << 481   return ComputeMeanFreePath(kinEnergy,FindParticle(particle),processName,
747     ConvertRangeToEnergy(part, mat, range);    << 482                              FindMaterial(material),cut);
748 }                                                 483 }
749                                                   484 
750 //....oooOO0OOooo........oooOO0OOooo........oo    485 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
751                                                   486 
752 G4bool G4EmCalculator::UpdateParticle(const G4 << 487 G4bool G4EmCalculator::UpdateParticle(const G4ParticleDefinition* p, G4double kinEnergy)
753                                       G4double << 
754 {                                                 488 {
755   if(p != currentParticle) {                      489   if(p != currentParticle) {
756                                                << 
757     // new particle                            << 
758     currentParticle = p;                          490     currentParticle = p;
759     dynParticle->SetDefinition(const_cast<G4Pa << 
760     dynParticle->SetKineticEnergy(kinEnergy);  << 
761     baseParticle    = nullptr;                 << 
762     currentParticleName = p->GetParticleName()    491     currentParticleName = p->GetParticleName();
763     massRatio       = 1.0;                     << 492   }
764     mass            = p->GetPDGMass();         << 493   isApplicable    = false;
765     chargeSquare    = 1.0;                     << 494   if(p->GetPDGCharge() != 0.0) {
766     currentProcess  = manager->GetEnergyLossPr << 495     const G4VEnergyLossProcess* proc = FindEnergyLossProcess(p);
767     currentProcessName = "";                   << 496     if(proc) {
768     isIon = false;                             << 497       isApplicable    = true;
769                                                << 498       baseParticle = proc->BaseParticle();
770     // ionisation process exist                << 499       if(p->GetParticleType() == "nucleus") {
771     if(nullptr != currentProcess) {            << 500   isIon = true;
772       currentProcessName = currentProcess->Get << 501         chargeSquare =
773       baseParticle = currentProcess->BaseParti << 502           ionEffCharge->EffectiveChargeSquareRatio(p, currentMaterial, kinEnergy);
774       if(currentProcessName == "ionIoni" && p- << 503       } else {
775         baseParticle = theGenericIon;          << 504   isIon = false;
776         isIon = true;                          << 505   G4double q = p->GetPDGCharge()/eplus;
                                                   >> 506   chargeSquare = q*q;
777       }                                           507       }
778                                                << 508       massRatio = 1.0;
779       // base particle is used                 << 509       if(baseParticle) {
780       if(nullptr != baseParticle) {            << 
781         massRatio = baseParticle->GetPDGMass()    510         massRatio = baseParticle->GetPDGMass()/p->GetPDGMass();
782         G4double q = p->GetPDGCharge()/basePar << 511         G4double q = baseParticle->GetPDGCharge()/eplus;
783         chargeSquare = q*q;                    << 512         chargeSquare /= (q*q);
784       }                                        << 513       }
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     }                                             514     }
796   }                                               515   }
797   return true;                                 << 516   return isApplicable;
798 }                                                 517 }
799                                                   518 
800 //....oooOO0OOooo........oooOO0OOooo........oo    519 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
801                                                   520 
802 const G4ParticleDefinition* G4EmCalculator::Fi    521 const G4ParticleDefinition* G4EmCalculator::FindParticle(const G4String& name)
803 {                                                 522 {
804   const G4ParticleDefinition* p = nullptr;     << 
805   if(name != currentParticleName) {               523   if(name != currentParticleName) {
806     p = G4ParticleTable::GetParticleTable()->F << 524     currentParticle = G4ParticleTable::GetParticleTable()->FindParticle(name);
807     if(nullptr == p) {                         << 525     if(!currentParticle) {
808       G4cout << "### WARNING: G4EmCalculator:: << 526       G4cout << "### WARNING: G4EmCalculator::FindParticle fails to find " << name << G4endl;
809              << name << G4endl;                << 
810     }                                             527     }
811   } else {                                     << 528     currentParticleName = name;
812     p = currentParticle;                       << 
813   }                                               529   }
814   return p;                                    << 530   return currentParticle;
815 }                                              << 
816                                                << 
817 //....oooOO0OOooo........oooOO0OOooo........oo << 
818                                                << 
819 const G4ParticleDefinition* G4EmCalculator::Fi << 
820 {                                              << 
821   const G4ParticleDefinition* p = ionTable->Ge << 
822   return p;                                    << 
823 }                                                 531 }
824                                                   532 
825 //....oooOO0OOooo........oooOO0OOooo........oo    533 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
826                                                   534 
827 const G4Material* G4EmCalculator::FindMaterial    535 const G4Material* G4EmCalculator::FindMaterial(const G4String& name)
828 {                                                 536 {
829   if(name != currentMaterialName) {               537   if(name != currentMaterialName) {
830     SetupMaterial(G4Material::GetMaterial(name << 538     currentMaterial = G4Material::GetMaterial(name);
831     if(nullptr == currentMaterial) {           << 539     currentMaterialName = name;
832       G4cout << "### WARNING: G4EmCalculator:: << 540     if(!currentMaterial)
833              << name << G4endl;                << 541       G4cout << "### WARNING: G4EmCalculator::FindMaterial fails to find " << name << G4endl;
834     }                                          << 
835   }                                               542   }
836   return currentMaterial;                         543   return currentMaterial;
837 }                                                 544 }
838                                                   545 
839 //....oooOO0OOooo........oooOO0OOooo........oo    546 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
840                                                   547 
841 const G4Region* G4EmCalculator::FindRegion(con    548 const G4Region* G4EmCalculator::FindRegion(const G4String& reg)
842 {                                                 549 {
843   return G4EmUtility::FindRegion(reg);         << 550   return G4RegionStore::GetInstance()->GetRegion(reg);
844 }                                                 551 }
845                                                   552 
846 //....oooOO0OOooo........oooOO0OOooo........oo    553 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
847                                                   554 
848 const G4MaterialCutsCouple* G4EmCalculator::Fi << 555 const G4MaterialCutsCouple* G4EmCalculator::FindCouple(const G4Material* material,
849                             const G4Material*  << 556                                                        const G4Region* region)
850                             const G4Region* re << 557 {
851 {                                              << 558   // Access to materials
852   const G4MaterialCutsCouple* couple = nullptr << 559   const G4ProductionCutsTable* theCoupleTable=
853   SetupMaterial(material);                     << 560         G4ProductionCutsTable::GetProductionCutsTable();
854   if(nullptr != currentMaterial) {             << 561   const G4Region* r = region;
855     // Access to materials                     << 562   if(!r) r = G4RegionStore::GetInstance()->GetRegion("DefaultRegionForTheWorld");
856     const G4ProductionCutsTable* theCoupleTabl << 563 
857       G4ProductionCutsTable::GetProductionCuts << 564   return theCoupleTable->GetMaterialCutsCouple(material,r->GetProductionCuts());
858     const G4Region* r = region;                << 565 
859     if(nullptr != r) {                         << 
860       couple = theCoupleTable->GetMaterialCuts << 
861                                                << 
862     } else {                                   << 
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 }                                                 566 }
884                                                   567 
885 //....oooOO0OOooo........oooOO0OOooo........oo    568 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
886                                                   569 
887 G4bool G4EmCalculator::UpdateCouple(const G4Ma    570 G4bool G4EmCalculator::UpdateCouple(const G4Material* material, G4double cut)
888 {                                                 571 {
889   SetupMaterial(material);                     << 572   if(!material) return false;
890   if(!currentMaterial) { return false; }       << 573   currentMaterial = material;
891   for (G4int i=0; i<nLocalMaterials; ++i) {    << 574   currentMaterialName = material->GetName();
                                                   >> 575   for (G4int i=0; i<nLocalMaterials; i++) {
892     if(material == localMaterials[i] && cut ==    576     if(material == localMaterials[i] && cut == localCuts[i]) {
893       currentCouple = localCouples[i];            577       currentCouple = localCouples[i];
894       currentCoupleIndex = currentCouple->GetI    578       currentCoupleIndex = currentCouple->GetIndex();
895       currentCut = cut;                           579       currentCut = cut;
896       return true;                                580       return true;
897     }                                             581     }
898   }                                               582   }
899   const G4MaterialCutsCouple* cc = new G4Mater    583   const G4MaterialCutsCouple* cc = new G4MaterialCutsCouple(material);
900   localMaterials.push_back(material);             584   localMaterials.push_back(material);
901   localCouples.push_back(cc);                     585   localCouples.push_back(cc);
902   localCuts.push_back(cut);                       586   localCuts.push_back(cut);
903   ++nLocalMaterials;                           << 587   nLocalMaterials++;
904   currentCouple = cc;                             588   currentCouple = cc;
905   currentCoupleIndex = currentCouple->GetIndex    589   currentCoupleIndex = currentCouple->GetIndex();
906   currentCut = cut;                               590   currentCut = cut;
907   return true;                                    591   return true;
908 }                                                 592 }
909                                                   593 
910 //....oooOO0OOooo........oooOO0OOooo........oo    594 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
911                                                   595 
912 void G4EmCalculator::FindLambdaTable(const G4P    596 void G4EmCalculator::FindLambdaTable(const G4ParticleDefinition* p,
913                                      const G4S << 597                                      const G4String& processName)
914                                      G4double  << 
915 {                                                 598 {
916   // Search for the process                       599   // Search for the process
917   if (!currentLambda || p != lambdaParticle || << 600   if (p != currentParticle || processName != currentName) {
918     lambdaName     = processName;              << 601     currentName     = processName;
919     currentLambda  = nullptr;                  << 602     currentLambda   = 0;
920     lambdaParticle = p;                        << 603     G4LossTableManager* manager = G4LossTableManager::Instance();
921     isApplicable   = false;                    << 604     const std::vector<G4VEnergyLossProcess*> vel = manager->GetEnergyLossProcessVector();
922                                                << 605     G4int n = vel.size();
923     const G4ParticleDefinition* part = (isIon) << 606     for(G4int i=0; i<n; i++) {
924                                                << 607       if((vel[i])->GetProcessName() == currentName && (vel[i])->Particle() == p) {
925     // Search for energy loss process          << 608         currentLambda = (vel[i])->LambdaTable();
926     currentName = processName;                 << 609   isApplicable    = true;
927     currentModel = nullptr;                    << 610   break;
928     loweModel = nullptr;                       << 611       }
929                                                << 612     }
930     G4VEnergyLossProcess* elproc = FindEnLossP << 613     if(!currentLambda) {
931     if(nullptr != elproc) {                    << 614       const std::vector<G4VEmProcess*> vem = manager->GetEmProcessVector();
932       currentLambda = elproc->LambdaTable();   << 615       G4int n = vem.size();
933       proctype      = 0;                       << 616       for(G4int i=0; i<n; i++) {
934       if(nullptr != currentLambda) {           << 617         if((vem[i])->GetProcessName() == currentName && (vem[i])->Particle() == p) {
935         isApplicable = true;                   << 618           currentLambda = (vem[i])->LambdaTable();
936         if(verbose>1) {                        << 619           isApplicable    = true;
937           G4cout << "G4VEnergyLossProcess is f << 620     break;
938                  << G4endl;                    << 
939         }                                         621         }
940       }                                           622       }
941       curProcess = elproc;                     << 
942       return;                                  << 
943     }                                             623     }
944                                                << 624     if(!currentLambda) {
945     // Search for discrete process             << 625       const std::vector<G4VMultipleScattering*> vmsc = manager->GetMultipleScatteringVector();
946     G4VEmProcess* proc = FindDiscreteProcess(p << 626       G4int n = vmsc.size();
947     if(nullptr != proc) {                      << 627       for(G4int i=0; i<n; i++) {
948       currentLambda = proc->LambdaTable();     << 628         if((vmsc[i])->GetProcessName() == currentName && (vmsc[i])->Particle() == p) {
949       proctype      = 1;                       << 629           currentLambda = (vmsc[i])->LambdaTable();
950       if(nullptr != currentLambda) {           << 630     isApplicable    = true;
951         isApplicable = true;                   << 631     break;
952         if(verbose>1) {                        << 
953           G4cout << "G4VEmProcess is found out << 
954         }                                      << 
955       }                                        << 
956       curProcess = proc;                       << 
957       return;                                  << 
958     }                                          << 
959                                                << 
960     // Search for msc process                  << 
961     G4VMultipleScattering* msc = FindMscProces << 
962     if(nullptr != msc) {                       << 
963       currentModel = msc->SelectModel(kinEnerg << 
964       proctype     = 2;                        << 
965       if(nullptr != currentModel) {            << 
966         currentLambda = currentModel->GetCross << 
967         if(nullptr != currentLambda) {         << 
968           isApplicable = true;                 << 
969           if(verbose>1) {                      << 
970             G4cout << "G4VMultipleScattering i << 
971                    << G4endl;                  << 
972           }                                    << 
973         }                                         632         }
974       }                                           633       }
975       curProcess = msc;                        << 
976     }                                             634     }
977   }                                               635   }
978 }                                                 636 }
979                                                   637 
980 //....oooOO0OOooo........oooOO0OOooo........oo    638 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
981                                                   639 
982 G4bool G4EmCalculator::FindEmModel(const G4Par    640 G4bool G4EmCalculator::FindEmModel(const G4ParticleDefinition* p,
983                                    const G4Str    641                                    const G4String& processName,
984                                             G4 << 642                    G4double kinEnergy)
985 {                                                 643 {
986   isApplicable = false;                        << 644   G4bool res = false;
987   if(nullptr == p || nullptr == currentMateria << 
988     G4cout << "G4EmCalculator::FindEmModel WAR << 
989            << " or materail defined; particle: << 
990     return isApplicable;                       << 
991   }                                            << 
992   G4String partname =  p->GetParticleName();   << 
993   G4double scaledEnergy = kinEnergy*massRatio; << 
994   const G4ParticleDefinition* part = (isIon) ? << 
995                                                << 
996   if(verbose > 1) {                               645   if(verbose > 1) {
997     G4cout << "## G4EmCalculator::FindEmModel  << 646     G4cout << "G4EmCalculator::FindEmModel for " << p->GetParticleName()
998            << " (type= " << p->GetParticleType << 647            << " and " << processName << " at e(MeV)= " << kinEnergy
999            << ") and " << processName << " at  << 
1000            << G4endl;                            648            << G4endl;
1001     if(p != part) { G4cout << "  GenericIon i << 
1002   }                                              649   }
1003                                               << 650   // Search for the process
1004   // Search for energy loss process           << 
1005   currentName = processName;                     651   currentName = processName;
1006   currentModel = nullptr;                     << 652   currentModel = 0;
1007   loweModel = nullptr;                        << 653   size_t idx   = 0;
1008   std::size_t idx = 0;                        << 654   G4LossTableManager* manager = G4LossTableManager::Instance();
1009                                               << 655   const std::vector<G4VEnergyLossProcess*> vel = manager->GetEnergyLossProcessVector();
1010   G4VEnergyLossProcess* elproc = FindEnLossPr << 656   G4int n = vel.size();
1011   if(nullptr != elproc) {                     << 657   for(G4int i=0; i<n; i++) {
1012     currentModel = elproc->SelectModelForMate << 658     if((vel[i])->GetProcessName() == currentName && (vel[i])->Particle() == p) {
1013     currentModel->InitialiseForMaterial(part, << 659       const G4ParticleDefinition* bp = (vel[i])->BaseParticle();
1014     currentModel->SetupForMaterial(part, curr << 660       //      G4cout << "i= " << i << " bp= " << bp << G4endl;
1015     G4double eth = currentModel->LowEnergyLim << 661       if(!bp) {
1016     if(eth > 0.0) {                           << 662         currentModel = (vel[i])->SelectModelForMaterial(kinEnergy, idx);
1017       loweModel = elproc->SelectModelForMater << 663       } else {
1018       if(loweModel == currentModel) { loweMod << 664         G4double e = kinEnergy*bp->GetPDGMass()/p->GetPDGMass();
1019       else {                                  << 665         for(G4int j=0; j<n; j++) {
1020         loweModel->InitialiseForMaterial(part << 666           if((vel[j])->Particle() == bp) {
1021         loweModel->SetupForMaterial(part, cur << 667             currentModel = (vel[j])->SelectModelForMaterial(e, idx);
1022       }                                       << 668             break;
1023     }                                         << 669     }
1024   }                                           << 670   }
1025                                               << 
1026   // Search for discrete process              << 
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       }                                          671       }
                                                   >> 672       break;
1094     }                                            673     }
1095   }                                              674   }
1096   return proc;                                << 675   if(!currentModel) {
1097 }                                             << 676     const std::vector<G4VEmProcess*> vem = manager->GetEmProcessVector();
1098                                               << 677     G4int n = vem.size();
1099 //....oooOO0OOooo........oooOO0OOooo........o << 678     for(G4int i=0; i<n; i++) {
1100                                               << 679       if((vem[i])->GetProcessName() == currentName && (vem[i])->Particle() == p) {
1101 G4VEmProcess*                                 << 680         currentModel = (vem[i])->SelectModelForMaterial(kinEnergy, idx);
1102 G4EmCalculator::FindDiscreteProcess(const G4P << 681   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;                                   682         break;
1118       }                                          683       }
1119     }                                            684     }
1120   }                                              685   }
1121   return proc;                                << 686   if(!currentModel) {
1122 }                                             << 687     const std::vector<G4VMultipleScattering*> vmsc = manager->GetMultipleScatteringVector();
1123                                               << 688     G4int n = vmsc.size();
1124 //....oooOO0OOooo........oooOO0OOooo........o << 689     for(G4int i=0; i<n; i++) {
1125                                               << 690       if((vmsc[i])->GetProcessName() == currentName && (vmsc[i])->Particle() == p) {
1126 G4VMultipleScattering*                        << 691         currentModel = (vmsc[i])->SelectModelForMaterial(kinEnergy, idx);
1127 G4EmCalculator::FindMscProcess(const G4Partic << 692   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;                                   693         break;
1140       }                                          694       }
1141     }                                            695     }
1142   }                                              696   }
1143   return proc;                                << 697   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;                                    698   return res;
1180 }                                                699 }
1181                                                  700 
1182 //....oooOO0OOooo........oooOO0OOooo........o    701 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1183                                                  702 
1184 void G4EmCalculator::SetupMaterial(const G4Ma << 703 const G4VEnergyLossProcess* G4EmCalculator::FindEnergyLossProcess(
                                                   >> 704                       const G4ParticleDefinition* p)
1185 {                                                705 {
1186   if(mat) {                                   << 706   const G4VEnergyLossProcess* elp = 0;
1187     currentMaterial = mat;                    << 707   
1188     currentMaterialName = mat->GetName();     << 708   G4LossTableManager* manager = G4LossTableManager::Instance();
1189   } else {                                    << 709   const std::vector<G4VEnergyLossProcess*> vel = manager->GetEnergyLossProcessVector();
1190     currentMaterial = nullptr;                << 710   G4int n = vel.size();
1191     currentMaterialName = "";                 << 711   for(G4int i=0; i<n; i++) {
1192   }                                           << 712     if((vel[i])->Particle() == p) {
1193 }                                             << 713       elp = vel[i];
1194                                               << 714       break;
1195 //....oooOO0OOooo........oooOO0OOooo........o << 
1196                                               << 
1197 void G4EmCalculator::SetupMaterial(const G4St << 
1198 {                                             << 
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     }                                            715     }
1215   }                                              716   }
1216   if(!isFound) {                              << 717   return elp;
1217     SetupMaterial(nist->FindOrBuildSimpleMate << 
1218   }                                           << 
1219 }                                                718 }
1220                                                  719 
1221 //....oooOO0OOooo........oooOO0OOooo........o    720 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1222                                                  721 
1223 void G4EmCalculator::SetVerbose(G4int verb)      722 void G4EmCalculator::SetVerbose(G4int verb)
1224 {                                                723 {
1225   verbose = verb;                                724   verbose = verb;
1226 }                                                725 }
1227                                                  726 
1228 //....oooOO0OOooo........oooOO0OOooo........o    727 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1229                                                  728 
1230                                                  729