Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/utils/src/G4EmCalculator.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/utils/src/G4EmCalculator.cc (Version 11.3.0) and /processes/electromagnetic/utils/src/G4EmCalculator.cc (Version 8.2)


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