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.1)


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