Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/standard/src/G4BetheBlochModel.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/standard/src/G4BetheBlochModel.cc (Version 11.3.0) and /processes/electromagnetic/standard/src/G4BetheBlochModel.cc (Version 9.4.p4)


  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: G4BetheBlochModel.cc,v 1.41 2010-11-12 18:37:47 vnivanch Exp $
                                                   >>  27 // GEANT4 tag $Name: not supported by cvs2svn $
                                                   >>  28 //
 26 // -------------------------------------------     29 // -------------------------------------------------------------------
 27 //                                                 30 //
 28 // GEANT4 Class header file                        31 // GEANT4 Class header file
 29 //                                                 32 //
 30 //                                                 33 //
 31 // File name:     G4BetheBlochModel                34 // File name:     G4BetheBlochModel
 32 //                                                 35 //
 33 // Author:        Vladimir Ivanchenko on base      36 // Author:        Vladimir Ivanchenko on base of Laszlo Urban code
 34 //                                                 37 //
 35 // Creation date: 03.01.2002                       38 // Creation date: 03.01.2002
 36 //                                                 39 //
 37 // Modifications:                                  40 // Modifications:
 38 //                                                 41 //
 39 // 04-12-02 Fix problem of G4DynamicParticle c     42 // 04-12-02 Fix problem of G4DynamicParticle constructor (V.Ivanchenko)
 40 // 23-12-02 Change interface in order to move      43 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
 41 // 27-01-03 Make models region aware (V.Ivanch     44 // 27-01-03 Make models region aware (V.Ivanchenko)
 42 // 13-02-03 Add name (V.Ivanchenko)                45 // 13-02-03 Add name (V.Ivanchenko)
 43 // 24-03-05 Add G4EmCorrections (V.Ivanchenko)     46 // 24-03-05 Add G4EmCorrections (V.Ivanchenko)
 44 // 11-04-05 Major optimisation of internal int     47 // 11-04-05 Major optimisation of internal interfaces (V.Ivanchenko)
 45 // 11-02-06 ComputeCrossSectionPerElectron, Co     48 // 11-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma)
 46 // 12-02-06 move G4LossTableManager::Instance(     49 // 12-02-06 move G4LossTableManager::Instance()->EmCorrections() 
 47 //          in constructor (mma)                   50 //          in constructor (mma)
 48 // 12-08-08 Added methods GetParticleCharge, G     51 // 12-08-08 Added methods GetParticleCharge, GetChargeSquareRatio, 
 49 //          CorrectionsAlongStep needed for io     52 //          CorrectionsAlongStep needed for ions(V.Ivanchenko)
 50 //                                                 53 //
 51 // -------------------------------------------     54 // -------------------------------------------------------------------
 52 //                                                 55 //
 53                                                    56 
                                                   >>  57 
 54 //....oooOO0OOooo........oooOO0OOooo........oo     58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 55 //....oooOO0OOooo........oooOO0OOooo........oo     59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 56                                                    60 
 57 #include "G4BetheBlochModel.hh"                    61 #include "G4BetheBlochModel.hh"
 58 #include "Randomize.hh"                            62 #include "Randomize.hh"
 59 #include "G4PhysicalConstants.hh"              << 
 60 #include "G4SystemOfUnits.hh"                  << 
 61 #include "G4NistManager.hh"                    << 
 62 #include "G4Electron.hh"                           63 #include "G4Electron.hh"
 63 #include "G4LossTableManager.hh"                   64 #include "G4LossTableManager.hh"
 64 #include "G4EmCorrections.hh"                      65 #include "G4EmCorrections.hh"
 65 #include "G4EmParameters.hh"                   << 
 66 #include "G4ParticleChangeForLoss.hh"              66 #include "G4ParticleChangeForLoss.hh"
 67 #include "G4ICRU90StoppingData.hh"             << 
 68 #include "G4Log.hh"                            << 
 69 #include "G4DeltaAngle.hh"                     << 
 70 #include <vector>                              << 
 71                                                    67 
 72 //....oooOO0OOooo........oooOO0OOooo........oo     68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 73                                                    69 
 74 G4BetheBlochModel::G4BetheBlochModel(const G4P <<  70 using namespace std;
                                                   >>  71 
                                                   >>  72 G4BetheBlochModel::G4BetheBlochModel(const G4ParticleDefinition* p, 
 75                                      const G4S     73                                      const G4String& nam)
 76   : G4VEmModel(nam),                               74   : G4VEmModel(nam),
 77     twoln10(2.0*G4Log(10.0)),                  <<  75     particle(0),
 78     fAlphaTlimit(1*CLHEP::GeV),                <<  76     tlimit(DBL_MAX),
 79     fProtonTlimit(10*CLHEP::GeV)               <<  77     twoln10(2.0*log(10.0)),
                                                   >>  78     bg2lim(0.0169),
                                                   >>  79     taulim(8.4146e-3),
                                                   >>  80     isIon(false),
                                                   >>  81     isInitialised(false)
 80 {                                                  82 {
                                                   >>  83   fParticleChange = 0;
 81   theElectron = G4Electron::Electron();            84   theElectron = G4Electron::Electron();
                                                   >>  85   if(p) {
                                                   >>  86     SetGenericIon(p);
                                                   >>  87     SetParticle(p);
                                                   >>  88   } else {
                                                   >>  89     SetParticle(theElectron);
                                                   >>  90   }
 82   corr = G4LossTableManager::Instance()->EmCor     91   corr = G4LossTableManager::Instance()->EmCorrections();  
 83   nist = G4NistManager::Instance();                92   nist = G4NistManager::Instance();
 84   SetLowEnergyLimit(2.0*CLHEP::MeV);           <<  93   SetLowEnergyLimit(2.0*MeV);
 85 }                                                  94 }
 86                                                    95 
 87 //....oooOO0OOooo........oooOO0OOooo........oo     96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 88                                                    97 
 89 G4BetheBlochModel::~G4BetheBlochModel() = defa <<  98 G4BetheBlochModel::~G4BetheBlochModel()
                                                   >>  99 {}
                                                   >> 100 
                                                   >> 101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 102 
                                                   >> 103 G4double G4BetheBlochModel::MinEnergyCut(const G4ParticleDefinition*,
                                                   >> 104                                          const G4MaterialCutsCouple* couple)
                                                   >> 105 {
                                                   >> 106   return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy();
                                                   >> 107 }
 90                                                   108 
 91 //....oooOO0OOooo........oooOO0OOooo........oo    109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 92                                                   110 
 93 void G4BetheBlochModel::Initialise(const G4Par    111 void G4BetheBlochModel::Initialise(const G4ParticleDefinition* p,
 94                                    const G4Dat    112                                    const G4DataVector&)
 95 {                                                 113 {
 96   if(p != particle) { SetupParameters(p); }    << 114   SetGenericIon(p);
                                                   >> 115   SetParticle(p);
                                                   >> 116 
                                                   >> 117   //G4cout << "G4BetheBlochModel::Initialise for " << p->GetParticleName()
                                                   >> 118   //   << "  isIon= " << isIon 
                                                   >> 119   //   << G4endl;
 97                                                   120 
 98   // always false before the run                  121   // always false before the run
 99   SetDeexcitationFlag(false);                     122   SetDeexcitationFlag(false);
100                                                   123 
101   // initialisation once                       << 124   if(!isInitialised) {
102   if(nullptr == fParticleChange) {             << 125     isInitialised = true;
103     const G4String& pname = particle->GetParti << 
104     if(G4EmParameters::Instance()->UseICRU90Da << 
105        (pname == "proton" || pname == "Generic << 
106       fICRU90 = nist->GetICRU90StoppingData(); << 
107     }                                          << 
108     if (pname == "GenericIon") {               << 
109       isIon = true;                            << 
110     } else if (pname == "alpha") {             << 
111       isAlpha = true;                          << 
112     } else if (particle->GetPDGCharge() > 1.1* << 
113       isIon = true;                            << 
114     }                                          << 
115                                                << 
116     fParticleChange = GetParticleChangeForLoss    126     fParticleChange = GetParticleChangeForLoss();
117     if(UseAngularGeneratorFlag() && nullptr == << 
118       SetAngularDistribution(new G4DeltaAngle( << 
119     }                                          << 
120   }                                            << 
121   // initialisation for each new run           << 
122   if(IsMaster() && nullptr != fICRU90) {       << 
123     fICRU90->Initialise();                     << 
124   }                                               127   }
125 }                                                 128 }
126                                                   129 
127 //....oooOO0OOooo........oooOO0OOooo........oo    130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
128                                                   131 
129 G4double G4BetheBlochModel::GetChargeSquareRat    132 G4double G4BetheBlochModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
130                                                << 133              const G4Material* mat,
131                                                << 134              G4double kineticEnergy)
132 {                                                 135 {
133   // this method is called only for ions, so n << 136   // this method is called only for ions
134   if(isAlpha) { return 1.0; }                  << 137   G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,kineticEnergy);
135   chargeSquare = corr->EffectiveChargeSquareRa << 138   corrFactor = q2*corr->EffectiveChargeCorrection(p,mat,kineticEnergy);
136   return chargeSquare;                         << 139   return corrFactor;
137 }                                                 140 }
138                                                   141 
139 //....oooOO0OOooo........oooOO0OOooo........oo    142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
140                                                   143 
141 G4double G4BetheBlochModel::GetParticleCharge(    144 G4double G4BetheBlochModel::GetParticleCharge(const G4ParticleDefinition* p,
142                                                << 145                 const G4Material* mat,
143                                                << 146                 G4double kineticEnergy)
144 {                                                 147 {
145   // this method is called only for ions, so n    148   // this method is called only for ions, so no check if it is an ion
146   return corr->GetParticleCharge(p, mat, kinet << 149   return corr->GetParticleCharge(p,mat,kineticEnergy);
147 }                                                 150 }
148                                                   151 
149 //....oooOO0OOooo........oooOO0OOooo........oo    152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
150                                                   153 
151 void G4BetheBlochModel::SetupParameters(const  << 154 void G4BetheBlochModel::SetupParameters()
152 {                                                 155 {
153   particle = p;                                << 
154   mass = particle->GetPDGMass();                  156   mass = particle->GetPDGMass();
155   spin = particle->GetPDGSpin();                  157   spin = particle->GetPDGSpin();
156   G4double q = particle->GetPDGCharge()*invepl << 158   G4double q = particle->GetPDGCharge()/eplus;
157   chargeSquare = q*q;                             159   chargeSquare = q*q;
                                                   >> 160   corrFactor = chargeSquare;
158   ratio = electron_mass_c2/mass;                  161   ratio = electron_mass_c2/mass;
159   constexpr G4double aMag = 1./(0.5*eplus*CLHE << 162   G4double magmom = 
160   G4double magmom = particle->GetPDGMagneticMo << 163     particle->GetPDGMagneticMoment()*mass/(0.5*eplus*hbar_Planck*c_squared);
161   magMoment2 = magmom*magmom - 1.0;               164   magMoment2 = magmom*magmom - 1.0;
162   formfact = 0.0;                                 165   formfact = 0.0;
163   tlimit = DBL_MAX;                            << 
164   if(particle->GetLeptonNumber() == 0) {          166   if(particle->GetLeptonNumber() == 0) {
165     G4double x = 0.8426*CLHEP::GeV;            << 167     G4double x = 0.8426*GeV;
166     if(spin == 0.0 && mass < CLHEP::GeV) { x = << 168     if(spin == 0.0 && mass < GeV) {x = 0.736*GeV;}
167     else if (mass > CLHEP::GeV) {              << 169     else if(mass > GeV) {
168       G4int iz = G4lrint(std::abs(q));         << 170       x /= nist->GetZ13(mass/proton_mass_c2);
169       if(iz > 1) { x /= nist->GetA27(iz); }    << 171       //  tlimit = 51.2*GeV*A13[iz]*A13[iz];
170     }                                             172     }
171     formfact = 2.0*CLHEP::electron_mass_c2/(x* << 173     formfact = 2.0*electron_mass_c2/(x*x);
172     tlimit = 2.0/formfact;                     << 174     tlimit   = 2.0/formfact;
173   }                                               175   }
174 }                                                 176 }
175                                                   177 
176 //....oooOO0OOooo........oooOO0OOooo........oo    178 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
177                                                   179 
178 G4double G4BetheBlochModel::MinEnergyCut(const << 
179                                          const << 
180 {                                              << 
181   return couple->GetMaterial()->GetIonisation( << 
182 }                                              << 
183                                                << 
184 //....oooOO0OOooo........oooOO0OOooo........oo << 
185                                                << 
186 G4double                                          180 G4double 
187 G4BetheBlochModel::ComputeCrossSectionPerElect    181 G4BetheBlochModel::ComputeCrossSectionPerElectron(const G4ParticleDefinition* p,
188                                                << 182               G4double kineticEnergy,
189                                                << 183               G4double cutEnergy,
190                                                << 184               G4double maxKinEnergy)  
191 {                                                 185 {
192   G4double cross = 0.0;                           186   G4double cross = 0.0;
193   const G4double tmax = MaxSecondaryEnergy(p,  << 187   G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
194   const G4double cutEnergy = std::min(std::min << 188   G4double maxEnergy = min(tmax,maxKinEnergy);
195   const G4double maxEnergy = std::min(tmax, ma << 
196   if(cutEnergy < maxEnergy) {                     189   if(cutEnergy < maxEnergy) {
197                                                   190 
198     G4double totEnergy = kineticEnergy + mass;    191     G4double totEnergy = kineticEnergy + mass;
199     G4double energy2   = totEnergy*totEnergy;     192     G4double energy2   = totEnergy*totEnergy;
200     G4double beta2     = kineticEnergy*(kineti    193     G4double beta2     = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
201                                                   194 
202     cross = (maxEnergy - cutEnergy)/(cutEnergy << 195     cross = 1.0/cutEnergy - 1.0/maxEnergy 
203       - beta2*G4Log(maxEnergy/cutEnergy)/tmax; << 196       - beta2*log(maxEnergy/cutEnergy)/tmax;
204                                                   197 
205     // +term for spin=1/2 particle                198     // +term for spin=1/2 particle
206     if( 0.0 < spin ) { cross += 0.5*(maxEnergy << 199     if( 0.5 == spin ) { cross += 0.5*(maxEnergy - cutEnergy)/energy2; }
                                                   >> 200 
                                                   >> 201     // High order correction different for hadrons and ions
                                                   >> 202     // nevetheless they are applied to reduce high energy transfers
                                                   >> 203     //    if(!isIon) 
                                                   >> 204     //cross += corr->FiniteSizeCorrectionXS(p,currentMaterial,
                                                   >> 205     //            kineticEnergy,cutEnergy);
207                                                   206 
208     cross *= CLHEP::twopi_mc2_rcl2*chargeSquar << 207     cross *= twopi_mc2_rcl2*chargeSquare/beta2;
209   }                                               208   }
210                                                   209   
211    // G4cout << "BB: e= " << kineticEnergy <<     210    // G4cout << "BB: e= " << kineticEnergy << " tmin= " << cutEnergy 
212    //        << " tmax= " << tmax << " cross=     211    //        << " tmax= " << tmax << " cross= " << cross << G4endl;
213                                                   212   
214   return cross;                                   213   return cross;
215 }                                                 214 }
216                                                   215 
217 //....oooOO0OOooo........oooOO0OOooo........oo    216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
218                                                   217 
219 G4double G4BetheBlochModel::ComputeCrossSectio    218 G4double G4BetheBlochModel::ComputeCrossSectionPerAtom(
220                                            con    219                                            const G4ParticleDefinition* p,
221                                                << 220                                                  G4double kineticEnergy,
222                                                << 221              G4double Z, G4double,
223                                                   222                                                  G4double cutEnergy,
224                                                   223                                                  G4double maxEnergy)
225 {                                                 224 {
226   return Z*ComputeCrossSectionPerElectron(p,ki << 225   G4double cross = Z*ComputeCrossSectionPerElectron
                                                   >> 226                                          (p,kineticEnergy,cutEnergy,maxEnergy);
                                                   >> 227   return cross;
227 }                                                 228 }
228                                                   229 
229 //....oooOO0OOooo........oooOO0OOooo........oo    230 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
230                                                   231 
231 G4double G4BetheBlochModel::CrossSectionPerVol    232 G4double G4BetheBlochModel::CrossSectionPerVolume(
232                                            con << 233              const G4Material* material,
233                                            con    234                                            const G4ParticleDefinition* p,
234                                                << 235                                                  G4double kineticEnergy,
235                                                   236                                                  G4double cutEnergy,
236                                                   237                                                  G4double maxEnergy)
237 {                                                 238 {
238   G4double sigma = mat->GetElectronDensity()   << 239   currentMaterial   = material;
239     *ComputeCrossSectionPerElectron(p,kinEnerg << 240   G4double eDensity = material->GetElectronDensity();
240   if(isAlpha) {                                << 241   G4double cross = eDensity*ComputeCrossSectionPerElectron
241     sigma *= corr->EffectiveChargeSquareRatio( << 242                                          (p,kineticEnergy,cutEnergy,maxEnergy);
242   }                                            << 243   return cross;
243   return sigma;                                << 
244 }                                                 244 }
245                                                   245 
246 //....oooOO0OOooo........oooOO0OOooo........oo    246 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
247                                                   247 
248 G4double G4BetheBlochModel::ComputeDEDXPerVolu    248 G4double G4BetheBlochModel::ComputeDEDXPerVolume(const G4Material* material,
249                                                << 249              const G4ParticleDefinition* p,
250                                                << 250              G4double kineticEnergy,
251                                                << 251              G4double cut)
252 {                                                 252 {
253   const G4double tmax = MaxSecondaryEnergy(p,  << 253   G4double tmax      = MaxSecondaryEnergy(p, kineticEnergy);
254   // projectile formfactor limit energy loss   << 254   G4double cutEnergy = std::min(cut,tmax);
255   const G4double cutEnergy = std::min(std::min << 
256                                                   255 
257   G4double tau   = kineticEnergy/mass;            256   G4double tau   = kineticEnergy/mass;
258   G4double gam   = tau + 1.0;                     257   G4double gam   = tau + 1.0;
259   G4double bg2   = tau * (tau+2.0);               258   G4double bg2   = tau * (tau+2.0);
260   G4double beta2 = bg2/(gam*gam);                 259   G4double beta2 = bg2/(gam*gam);
261   G4double xc    = cutEnergy/tmax;             << 
262                                                   260 
263   G4double eexc  = material->GetIonisation()->    261   G4double eexc  = material->GetIonisation()->GetMeanExcitationEnergy();
264   G4double eexc2 = eexc*eexc;                     262   G4double eexc2 = eexc*eexc;
265                                                   263 
266   G4double eDensity = material->GetElectronDen    264   G4double eDensity = material->GetElectronDensity();
267                                                   265 
268   // added ICRU90 stopping data for limited li << 266   G4double dedx = log(2.0*electron_mass_c2*bg2*cutEnergy/eexc2)
269   /*                                           << 267                 - (1.0 + cutEnergy/tmax)*beta2;
270   G4cout << "### DEDX ICRI90:" << (nullptr !=  << 
271    << " Ekin=" << kineticEnergy                << 
272    << "  " << p->GetParticleName()             << 
273    << " q2=" << chargeSquare                   << 
274    << " inside  " << material->GetName() << G4 << 
275   */                                           << 
276   if(nullptr != fICRU90 && kineticEnergy < fPr << 
277     if(material != currentMaterial) {          << 
278       currentMaterial = material;              << 
279       baseMaterial = material->GetBaseMaterial << 
280         ? material->GetBaseMaterial() : materi << 
281       iICRU90 = fICRU90->GetIndex(baseMaterial << 
282     }                                          << 
283     if(iICRU90 >= 0) {                         << 
284       G4double dedx = 0.0;                     << 
285       // only for alpha                        << 
286       if(isAlpha) {                            << 
287   if(kineticEnergy <= fAlphaTlimit) {          << 
288     dedx = fICRU90->GetElectronicDEDXforAlpha( << 
289   } else {                                     << 
290           const G4double e = kineticEnergy*CLH << 
291     dedx = fICRU90->GetElectronicDEDXforProton << 
292   }                                            << 
293       } else {                                 << 
294         dedx = fICRU90->GetElectronicDEDXforPr << 
295     *chargeSquare;                             << 
296       }                                        << 
297       dedx *= material->GetDensity();          << 
298       if(cutEnergy < tmax) {                   << 
299         dedx += (G4Log(xc) + (1.0 - xc)*beta2) << 
300           *(eDensity*chargeSquare/beta2);      << 
301       }                                        << 
302       //G4cout << "   iICRU90=" << iICRU90 <<  << 
303       if(dedx > 0.0) { return dedx; }          << 
304     }                                          << 
305   }                                            << 
306   // general Bethe-Bloch formula               << 
307   G4double dedx = G4Log(2.0*CLHEP::electron_ma << 
308                 - (1.0 + xc)*beta2;            << 
309                                                   268 
310   if(0.0 < spin) {                             << 269   if(0.5 == spin) {
311     G4double del = 0.5*cutEnergy/(kineticEnerg    270     G4double del = 0.5*cutEnergy/(kineticEnergy + mass);
312     dedx += del*del;                              271     dedx += del*del;
313   }                                               272   }
314                                                   273 
315   // density correction                           274   // density correction
316   G4double x = G4Log(bg2)/twoln10;             << 275   G4double x = log(bg2)/twoln10;
317   dedx -= material->GetIonisation()->DensityCo    276   dedx -= material->GetIonisation()->DensityCorrection(x);
318                                                   277 
319   // shell correction                             278   // shell correction
320   dedx -= 2.0*corr->ShellCorrection(p,material    279   dedx -= 2.0*corr->ShellCorrection(p,material,kineticEnergy);
321                                                   280 
322   // now compute the total ionization loss        281   // now compute the total ionization loss
323   dedx *= CLHEP::twopi_mc2_rcl2*chargeSquare*e << 282   dedx *= twopi_mc2_rcl2*chargeSquare*eDensity/beta2;
324                                                   283 
325   //High order correction different for hadron    284   //High order correction different for hadrons and ions
326   if(isIon) {                                     285   if(isIon) {
327     dedx += corr->IonBarkasCorrection(p,materi    286     dedx += corr->IonBarkasCorrection(p,material,kineticEnergy);
328   } else {                                        287   } else {      
329     dedx += corr->HighOrderCorrections(p,mater    288     dedx += corr->HighOrderCorrections(p,material,kineticEnergy,cutEnergy);
330   }                                               289   }
331                                                   290 
332   dedx = std::max(dedx, 0.0);                  << 291   if (dedx < 0.0) { dedx = 0.0; }
333   /*                                           << 
334   G4cout << "E(MeV)= " << kineticEnergy/CLHEP: << 
335            << "  " << material->GetName() << G << 
336   */                                           << 
337   return dedx;                                    292   return dedx;
338 }                                                 293 }
339                                                   294 
340 //....oooOO0OOooo........oooOO0OOooo........oo    295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
341                                                   296 
342 void G4BetheBlochModel::CorrectionsAlongStep(c    297 void G4BetheBlochModel::CorrectionsAlongStep(const G4MaterialCutsCouple* couple,
343                                              c << 298                const G4DynamicParticle* dp,
344                                              c << 299                G4double& eloss,
345                                              G << 300                G4double&,
346 {                                              << 301                G4double length)
347   // no correction for alpha                   << 302 {
348   if(isAlpha) { return; }                      << 303   if(isIon) {
349                                                << 304     const G4ParticleDefinition* p = dp->GetDefinition();
350   // no correction at the last step or at smal << 305     const G4Material* mat = couple->GetMaterial();
351   const G4double preKinEnergy = dp->GetKinetic << 306     G4double preKinEnergy = dp->GetKineticEnergy();
352   if(eloss >= preKinEnergy || eloss < preKinEn << 307     G4double e = preKinEnergy - eloss*0.5;
353                                                << 308     if(e < 0.0) e = preKinEnergy*0.5;
354   // corrections for all charged particles wit << 309 
355   const G4ParticleDefinition* p = dp->GetDefin << 310     G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,e);
356   if(p != particle) { SetupParameters(p); }    << 311     GetModelOfFluctuations()->SetParticleAndCharge(p, q2);
357   if(!isIon) { return; }                       << 312     G4double qfactor = q2*corr->EffectiveChargeCorrection(p,mat,e)/corrFactor;
358                                                << 313     G4double highOrder = length*corr->IonHighOrderCorrections(p,couple,e);
359   // effective energy and charge at a step     << 314     eloss *= qfactor; 
360   const G4double e = std::max(preKinEnergy - e << 315     eloss += highOrder;
361   const G4Material* mat = couple->GetMaterial( << 316     //G4cout << "G4BetheBlochModel::CorrectionsAlongStep: e= " << preKinEnergy
362   const G4double q20 = corr->EffectiveChargeSq << 317     //     << " qfactor= " << qfactor 
363   const G4double q2 = corr->EffectiveChargeSqu << 318     //     << " highOrder= " << highOrder << " (" << highOrder/eloss << ")" << G4endl;    
364   const G4double qfactor = q2/q20;             << 319   }
365                                                << 
366   /*                                           << 
367     G4cout << "G4BetheBlochModel::CorrectionsA << 
368     << preKinEnergy << " Eeff(MeV)=" << e      << 
369     << " eloss=" << eloss << " elossnew=" << e << 
370     << " qfactor=" << qfactor << " Qpre=" << q << 
371     << p->GetParticleName() <<G4endl;          << 
372   */                                           << 
373   eloss *= qfactor;                            << 
374 }                                                 320 }
375                                                   321 
376 //....oooOO0OOooo........oooOO0OOooo........oo    322 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
377                                                   323 
378 void G4BetheBlochModel::SampleSecondaries(std: << 324 void G4BetheBlochModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
379                                           cons << 325             const G4MaterialCutsCouple*,
380                                           cons << 326             const G4DynamicParticle* dp,
381                                           G4do << 327             G4double minKinEnergy,
382                                           G4do << 328             G4double maxEnergy)
383 {                                              << 329 {
384   G4double kinEnergy = dp->GetKineticEnergy(); << 330   G4double kineticEnergy = dp->GetKineticEnergy();
385   const G4double tmax = MaxSecondaryEnergy(dp- << 331   G4double tmax = MaxSecondaryEnergy(dp->GetDefinition(),kineticEnergy);
386   const G4double minKinEnergy = std::min(cut,  << 332 
387   const G4double maxKinEnergy = std::min(maxEn << 333   G4double maxKinEnergy = std::min(maxEnergy,tmax);
388   if(minKinEnergy >= maxKinEnergy) { return; } << 334   if(minKinEnergy >= maxKinEnergy) return;
389                                                << 335 
390   //G4cout << "G4BetheBlochModel::SampleSecond << 336   G4double totEnergy     = kineticEnergy + mass;
391   //         << " Emax= " << maxKinEnergy << G << 337   G4double etot2         = totEnergy*totEnergy;
392                                                << 338   G4double beta2         = kineticEnergy*(kineticEnergy + 2.0*mass)/etot2;
393   const G4double totEnergy = kinEnergy + mass; << 
394   const G4double etot2 = totEnergy*totEnergy;  << 
395   const G4double beta2 = kinEnergy*(kinEnergy  << 
396                                                   339 
397   G4double deltaKinEnergy, f;                     340   G4double deltaKinEnergy, f; 
398   G4double f1 = 0.0;                              341   G4double f1 = 0.0;
399   G4double fmax = 1.0;                            342   G4double fmax = 1.0;
400   if( 0.0 < spin ) { fmax += 0.5*maxKinEnergy* << 343   if( 0.5 == spin ) fmax += 0.5*maxKinEnergy*maxKinEnergy/etot2; 
401                                                << 
402   CLHEP::HepRandomEngine* rndmEngineMod = G4Ra << 
403   G4double rndm[2];                            << 
404                                                   344 
405   // sampling without nuclear size effect         345   // sampling without nuclear size effect
406   do {                                            346   do {
407     rndmEngineMod->flatArray(2, rndm);         << 347     G4double q = G4UniformRand();
408     deltaKinEnergy = minKinEnergy*maxKinEnergy    348     deltaKinEnergy = minKinEnergy*maxKinEnergy
409                     /(minKinEnergy*(1.0 - rndm << 349                     /(minKinEnergy*(1.0 - q) + maxKinEnergy*q);
410                                                   350 
411     f = 1.0 - beta2*deltaKinEnergy/tmax;          351     f = 1.0 - beta2*deltaKinEnergy/tmax;
412     if( 0.0 < spin ) {                         << 352     if( 0.5 == spin ) {
413       f1 = 0.5*deltaKinEnergy*deltaKinEnergy/e    353       f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
414       f += f1;                                    354       f += f1;
415     }                                             355     }
416                                                   356 
417     // Loop checking, 03-Aug-2015, Vladimir Iv << 357   } while( fmax*G4UniformRand() > f);
418   } while( fmax*rndm[1] > f);                  << 
419                                                   358 
420   // projectile formfactor - suppresion of hig    359   // projectile formfactor - suppresion of high energy
421   // delta-electron production at high energy     360   // delta-electron production at high energy
422                                                   361   
423   G4double x = formfact*deltaKinEnergy;           362   G4double x = formfact*deltaKinEnergy;
424   if(x > 1.e-6) {                                 363   if(x > 1.e-6) {
425                                                   364 
426     G4double x1 = 1.0 + x;                        365     G4double x1 = 1.0 + x;
427     G4double grej  = 1.0/(x1*x1);              << 366     G4double g  = 1.0/(x1*x1);
428     if( 0.0 < spin ) {                         << 367     if( 0.5 == spin ) {
429       G4double x2 = 0.5*electron_mass_c2*delta    368       G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
430       grej *= (1.0 + magMoment2*(x2 - f1/f)/(1 << 369       g *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
431     }                                             370     }
432     if(grej > 1.1) {                           << 371     if(g > 1.0) {
433       G4cout << "### G4BetheBlochModel WARNING << 372       G4cout << "### G4BetheBlochModel WARNING: g= " << g
434              << "  " << dp->GetDefinition()->G << 373        << dp->GetDefinition()->GetParticleName()
435              << " Ekin(MeV)= " <<  kinEnergy   << 374        << " Ekin(MeV)= " <<  kineticEnergy
436              << " delEkin(MeV)= " << deltaKinE << 375        << " delEkin(MeV)= " << deltaKinEnergy
437              << G4endl;                        << 376        << G4endl;
438     }                                             377     }
439     if(rndmEngineMod->flat() > grej) { return; << 378     if(G4UniformRand() > g) return;
440   }                                               379   }
441                                                   380 
442   G4ThreeVector deltaDirection;                << 381   // delta-electron is produced
443                                                << 382   G4double totMomentum = totEnergy*sqrt(beta2);
444   if(UseAngularGeneratorFlag()) {              << 383   G4double deltaMomentum =
445     const G4Material* mat = couple->GetMateria << 384            sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
446     deltaDirection =                           << 385   G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) /
447       GetAngularDistribution()->SampleDirectio << 386                                    (deltaMomentum * totMomentum);
448             SelectRandomAtomNumber(mat),       << 
449             mat);                              << 
450   } else {                                     << 
451                                                << 
452     G4double deltaMomentum =                   << 
453       std::sqrt(deltaKinEnergy * (deltaKinEner << 
454     G4double cost = deltaKinEnergy * (totEnerg << 
455       (deltaMomentum * dp->GetTotalMomentum()) << 
456     cost = std::min(cost, 1.0);                << 
457     const G4double sint = std::sqrt((1.0 - cos << 
458     const G4double phi = twopi*rndmEngineMod-> << 
459                                                << 
460     deltaDirection.set(sint*std::cos(phi),sint << 
461     deltaDirection.rotateUz(dp->GetMomentumDir << 
462   }                                            << 
463   /*                                              387   /*
464     G4cout << "### G4BetheBlochModel "         << 388   if(cost > 1.0) {
465            << dp->GetDefinition()->GetParticle << 389     G4cout << "### G4BetheBlochModel WARNING: cost= " 
466            << " Ekin(MeV)= " <<  kinEnergy     << 390      << cost << " > 1 for "
467            << " delEkin(MeV)= " << deltaKinEne << 391      << dp->GetDefinition()->GetParticleName()
468            << " tmin(MeV)= " << minKinEnergy   << 392      << " Ekin(MeV)= " <<  kineticEnergy
469            << " tmax(MeV)= " << maxKinEnergy   << 393      << " p(MeV/c)= " <<  totMomentum
                                                   >> 394      << " delEkin(MeV)= " << deltaKinEnergy
                                                   >> 395      << " delMom(MeV/c)= " << deltaMomentum
                                                   >> 396      << " tmin(MeV)= " << minKinEnergy
                                                   >> 397      << " tmax(MeV)= " << maxKinEnergy
470            << " dir= " << dp->GetMomentumDirec    398            << " dir= " << dp->GetMomentumDirection()
471            << " dirDelta= " << deltaDirection  << 399      << G4endl;
472            << G4endl;                          << 400     cost = 1.0;
                                                   >> 401   }
473   */                                              402   */
                                                   >> 403   G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
                                                   >> 404 
                                                   >> 405   G4double phi = twopi * G4UniformRand() ;
                                                   >> 406 
                                                   >> 407 
                                                   >> 408   G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost);
                                                   >> 409   G4ThreeVector direction = dp->GetMomentumDirection();
                                                   >> 410   deltaDirection.rotateUz(direction);
                                                   >> 411 
474   // create G4DynamicParticle object for delta    412   // create G4DynamicParticle object for delta ray
475   auto delta = new G4DynamicParticle(theElectr << 413   G4DynamicParticle* delta = new G4DynamicParticle(theElectron,
                                                   >> 414                deltaDirection,deltaKinEnergy);
476                                                   415 
477   vdp->push_back(delta);                          416   vdp->push_back(delta);
478                                                   417 
479   // Change kinematics of primary particle        418   // Change kinematics of primary particle
480   kinEnergy -= deltaKinEnergy;                 << 419   kineticEnergy       -= deltaKinEnergy;
481   G4ThreeVector finalP = dp->GetMomentum() - d << 420   G4ThreeVector finalP = direction*totMomentum - deltaDirection*deltaMomentum;
482   finalP = finalP.unit();                      << 421   finalP               = finalP.unit();
483                                                   422   
484   fParticleChange->SetProposedKineticEnergy(ki << 423   fParticleChange->SetProposedKineticEnergy(kineticEnergy);
485   fParticleChange->SetProposedMomentumDirectio    424   fParticleChange->SetProposedMomentumDirection(finalP);
486 }                                                 425 }
487                                                   426 
488 //....oooOO0OOooo........oooOO0OOooo........oo    427 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
489                                                   428 
490 G4double G4BetheBlochModel::MaxSecondaryEnergy    429 G4double G4BetheBlochModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd,
491                                                << 430                  G4double kinEnergy) 
492 {                                                 431 {
493   // here particle type is checked for the cas << 432   // here particle type is checked for any method
494   // when this model is shared between particl << 433   SetParticle(pd);
495   if(pd != particle) { SetupParameters(pd); }  << 
496   G4double tau  = kinEnergy/mass;                 434   G4double tau  = kinEnergy/mass;
497   return 2.0*CLHEP::electron_mass_c2*tau*(tau  << 435   G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
498     (1. + 2.0*(tau + 1.)*ratio + ratio*ratio); << 436                   (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
                                                   >> 437   return std::min(tmax,tlimit);
499 }                                                 438 }
500                                                   439 
501 //....oooOO0OOooo........oooOO0OOooo........oo    440 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
502                                                   441