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


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