Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/highenergy/src/G4mplIonisationModel.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/highenergy/src/G4mplIonisationModel.cc (Version 11.3.0) and /processes/electromagnetic/highenergy/src/G4mplIonisationModel.cc (Version 9.1.p2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
                                                   >>  26 // $Id: G4mplIonisationModel.cc,v 1.5 2007/11/13 18:36:29 vnivanch Exp $
                                                   >>  27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
 26 //                                                 28 //
 27 // -------------------------------------------     29 // -------------------------------------------------------------------
 28 //                                                 30 //
 29 // GEANT4 Class header file                        31 // GEANT4 Class header file
 30 //                                                 32 //
 31 //                                                 33 //
 32 // File name:     G4mplIonisationModel             34 // File name:     G4mplIonisationModel
 33 //                                                 35 //
 34 // Author:        Vladimir Ivanchenko              36 // Author:        Vladimir Ivanchenko 
 35 //                                                 37 //
 36 // Creation date: 06.09.2005                       38 // Creation date: 06.09.2005
 37 //                                                 39 //
 38 // Modifications:                                  40 // Modifications:
 39 // 12.08.2007 Changing low energy approximatio     41 // 12.08.2007 Changing low energy approximation and extrapolation. 
 40 //            Small bug fixing and refactoring     42 //            Small bug fixing and refactoring (M. Vladymyrov)
 41 // 13.11.2007 Use low-energy asymptotic from [     43 // 13.11.2007 Use low-energy asymptotic from [3] (V.Ivanchenko) 
 42 //                                                 44 //
 43 //                                                 45 //
 44 // -------------------------------------------     46 // -------------------------------------------------------------------
 45 // References                                      47 // References
 46 // [1] Steven P. Ahlen: Energy loss of relativ     48 // [1] Steven P. Ahlen: Energy loss of relativistic heavy ionizing particles, 
 47 //     S.P. Ahlen, Rev. Mod. Phys 52(1980), p1     49 //     S.P. Ahlen, Rev. Mod. Phys 52(1980), p121
 48 // [2] K.A. Milton arXiv:hep-ex/0602040            50 // [2] K.A. Milton arXiv:hep-ex/0602040
 49 // [3] S.P. Ahlen and K. Kinoshita, Phys. Rev.     51 // [3] S.P. Ahlen and K. Kinoshita, Phys. Rev. D26 (1982) 2347
 50                                                    52 
 51                                                    53 
 52 //....oooOO0OOooo........oooOO0OOooo........oo     54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 53 //....oooOO0OOooo........oooOO0OOooo........oo     55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 54                                                    56 
 55 #include "G4mplIonisationModel.hh"                 57 #include "G4mplIonisationModel.hh"
 56 #include "Randomize.hh"                            58 #include "Randomize.hh"
 57 #include "G4PhysicalConstants.hh"              <<  59 #include "G4LossTableManager.hh"
 58 #include "G4SystemOfUnits.hh"                  << 
 59 #include "G4ParticleChangeForLoss.hh"              60 #include "G4ParticleChangeForLoss.hh"
 60 #include "G4ProductionCutsTable.hh"            << 
 61 #include "G4MaterialCutsCouple.hh"             << 
 62 #include "G4Log.hh"                            << 
 63 #include "G4Pow.hh"                            << 
 64                                                    61 
 65 //....oooOO0OOooo........oooOO0OOooo........oo     62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 66                                                    63 
 67 std::vector<G4double>* G4mplIonisationModel::d <<  64 using namespace std;
 68                                                    65 
 69 G4mplIonisationModel::G4mplIonisationModel(G4d     66 G4mplIonisationModel::G4mplIonisationModel(G4double mCharge, const G4String& nam)
 70   : G4VEmModel(nam),G4VEmFluctuationModel(nam)     67   : G4VEmModel(nam),G4VEmFluctuationModel(nam),
 71   magCharge(mCharge),                              68   magCharge(mCharge),
 72   twoln10(G4Log(100.0)),                       <<  69   twoln10(log(100.0)),
 73   betalow(0.01),                                   70   betalow(0.01),
 74   betalim(0.1),                                    71   betalim(0.1),
 75   beta2lim(betalim*betalim),                       72   beta2lim(betalim*betalim),
 76   bg2lim(beta2lim*(1.0 + beta2lim))                73   bg2lim(beta2lim*(1.0 + beta2lim))
 77 {                                                  74 {
 78   nmpl = G4int(std::abs(magCharge) * 2 * CLHEP <<  75   nmpl         = G4int(abs(magCharge) * 2 * fine_structure_const + 0.5);
 79   if(nmpl > 6)      { nmpl = 6; }              <<  76   if(nmpl > 6)      nmpl = 6;
 80   else if(nmpl < 1) { nmpl = 1; }              <<  77   else if(nmpl < 1) nmpl = 1;
 81   pi_hbarc2_over_mc2 = CLHEP::pi*CLHEP::hbarc* <<  78   pi_hbarc2_over_mc2 = pi * hbarc * hbarc / electron_mass_c2;
 82   chargeSquare = magCharge * magCharge;            79   chargeSquare = magCharge * magCharge;
 83   dedxlim = 45.*nmpl*nmpl*CLHEP::GeV*CLHEP::cm <<  80   dedxlim = 45.*nmpl*nmpl*GeV*cm2/g;
 84 }                                                  81 }
 85                                                    82 
 86 //....oooOO0OOooo........oooOO0OOooo........oo     83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 87                                                    84 
 88 G4mplIonisationModel::~G4mplIonisationModel()      85 G4mplIonisationModel::~G4mplIonisationModel()
 89 {                                              <<  86 {}
 90   if(IsMaster()) { delete dedx0; }             << 
 91 }                                              << 
 92                                                    87 
 93 //....oooOO0OOooo........oooOO0OOooo........oo     88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 94                                                    89 
 95 void G4mplIonisationModel::SetParticle(const G <<  90 void G4mplIonisationModel::Initialise(const G4ParticleDefinition* p,
                                                   >>  91               const G4DataVector&)
 96 {                                                  92 {
 97   monopole = p;                                    93   monopole = p;
 98   mass     = monopole->GetPDGMass();               94   mass     = monopole->GetPDGMass();
 99   G4double emin =                              << 
100     std::min(LowEnergyLimit(),0.1*mass*(1./std << 
101   G4double emax =                              << 
102     std::max(HighEnergyLimit(),10.*mass*(1./st << 
103   SetLowEnergyLimit(emin);                     << 
104   SetHighEnergyLimit(emax);                    << 
105 }                                              << 
106                                                    95 
107 //....oooOO0OOooo........oooOO0OOooo........oo <<  96   if(pParticleChange) 
108                                                <<  97     fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
109 void G4mplIonisationModel::Initialise(const G4 <<  98   else 
110               const G4DataVector&)             <<  99     fParticleChange = new G4ParticleChangeForLoss();
111 {                                              << 
112   if(nullptr == monopole) { SetParticle(p); }  << 
113   if(nullptr == fParticleChange) { fParticleCh << 
114   if(IsMaster()) {                             << 
115     if(nullptr == dedx0) { dedx0 = new std::ve << 
116     G4ProductionCutsTable* theCoupleTable=     << 
117       G4ProductionCutsTable::GetProductionCuts << 
118     G4int numOfCouples = (G4int)theCoupleTable << 
119     G4int n = (G4int)dedx0->size();            << 
120     if(n < numOfCouples) { dedx0->resize(numOf << 
121                                                << 
122     G4Pow* g4calc = G4Pow::GetInstance();      << 
123                                                << 
124     // initialise vector assuming low conducti << 
125     for(G4int i=0; i<numOfCouples; ++i) {      << 
126                                                << 
127       const G4Material* material =             << 
128         theCoupleTable->GetMaterialCutsCouple( << 
129       G4double eDensity = material->GetElectro << 
130       G4double vF2 = 2*electron_Compton_length << 
131       (*dedx0)[i] = pi_hbarc2_over_mc2*eDensit << 
132         (G4Log(vF2/fine_structure_const) - 0.5 << 
133     }                                          << 
134   }                                            << 
135 }                                                 100 }
136                                                   101 
137 //....oooOO0OOooo........oooOO0OOooo........oo    102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
138                                                   103 
139 G4double G4mplIonisationModel::ComputeDEDXPerV    104 G4double G4mplIonisationModel::ComputeDEDXPerVolume(const G4Material* material,
140                 const G4ParticleDefinition* p, << 105                 const G4ParticleDefinition*,
141                 G4double kineticEnergy,           106                 G4double kineticEnergy,
142                 G4double)                         107                 G4double)
143 {                                                 108 {
144   if(nullptr == monopole) { SetParticle(p); }  << 
145   G4double tau   = kineticEnergy / mass;          109   G4double tau   = kineticEnergy / mass;
146   G4double gam   = tau + 1.0;                     110   G4double gam   = tau + 1.0;
147   G4double bg2   = tau * (tau + 2.0);             111   G4double bg2   = tau * (tau + 2.0);
148   G4double beta2 = bg2 / (gam * gam);             112   G4double beta2 = bg2 / (gam * gam);
149   G4double beta  = std::sqrt(beta2);           << 113   G4double beta  = sqrt(beta2);
150                                                   114 
151   // low-energy asymptotic formula                115   // low-energy asymptotic formula
152   //G4double dedx  = dedxlim*beta*material->Ge << 116   G4double dedx  = dedxlim*beta*material->GetDensity();
153   G4double dedx = (*dedx0)[CurrentCouple()->Ge << 
154                                                   117 
155   // above asymptotic                             118   // above asymptotic
156   if(beta > betalow) {                            119   if(beta > betalow) {
157                                                   120 
158     // high energy                                121     // high energy
159     if(beta >= betalim) {                         122     if(beta >= betalim) {
160       dedx = ComputeDEDXAhlen(material, bg2);     123       dedx = ComputeDEDXAhlen(material, bg2);
161                                                   124 
162     } else {                                      125     } else {
163                                                   126 
164       //G4double dedx1 = dedxlim*betalow*mater << 127       G4double dedx1 = dedxlim*betalow*material->GetDensity();
165       G4double dedx1 = (*dedx0)[CurrentCouple( << 
166       G4double dedx2 = ComputeDEDXAhlen(materi    128       G4double dedx2 = ComputeDEDXAhlen(material, bg2lim);
167                                                   129 
168       // extrapolation between two formula        130       // extrapolation between two formula 
169       G4double kapa2 = beta - betalow;            131       G4double kapa2 = beta - betalow;
170       G4double kapa1 = betalim - beta;            132       G4double kapa1 = betalim - beta;
171       dedx = (kapa1*dedx1 + kapa2*dedx2)/(kapa    133       dedx = (kapa1*dedx1 + kapa2*dedx2)/(kapa1 + kapa2);
172     }                                             134     }
173   }                                               135   }
174   return dedx;                                    136   return dedx;
175 }                                                 137 }
176                                                   138 
177 //....oooOO0OOooo........oooOO0OOooo........oo    139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
178                                                   140 
179 G4double G4mplIonisationModel::ComputeDEDXAhle << 141 G4double G4mplIonisationModel::ComputeDEDXAhlen(const G4Material* material, G4double bg2)
180             G4double bg2)                      << 
181 {                                                 142 {
182   G4double eDensity = material->GetElectronDen    143   G4double eDensity = material->GetElectronDensity();
183   G4double eexc  = material->GetIonisation()->    144   G4double eexc  = material->GetIonisation()->GetMeanExcitationEnergy();
184   G4double cden  = material->GetIonisation()->    145   G4double cden  = material->GetIonisation()->GetCdensity();
185   G4double mden  = material->GetIonisation()->    146   G4double mden  = material->GetIonisation()->GetMdensity();
186   G4double aden  = material->GetIonisation()->    147   G4double aden  = material->GetIonisation()->GetAdensity();
187   G4double x0den = material->GetIonisation()->    148   G4double x0den = material->GetIonisation()->GetX0density();
188   G4double x1den = material->GetIonisation()->    149   G4double x1den = material->GetIonisation()->GetX1density();
189                                                   150 
190   // Ahlen's formula for nonconductors, [1]p15    151   // Ahlen's formula for nonconductors, [1]p157, f(5.7)
191   G4double dedx = std::log(2.0 * electron_mass << 152   G4double dedx = log(2.0 * electron_mass_c2 * bg2 / eexc) - 0.5;
192                                                   153 
193   // Kazama et al. cross-section correction       154   // Kazama et al. cross-section correction
194   G4double  k = 0.406;                            155   G4double  k = 0.406;
195   if(nmpl > 1) k = 0.346;                         156   if(nmpl > 1) k = 0.346;
196                                                   157 
197   // Bloch correction                             158   // Bloch correction
198   const G4double B[7] = { 0.0, 0.248, 0.672, 1    159   const G4double B[7] = { 0.0, 0.248, 0.672, 1.022, 1.243, 1.464, 1.685}; 
199                                                   160 
200   dedx += 0.5 * k - B[nmpl];                      161   dedx += 0.5 * k - B[nmpl];
201                                                   162 
202   // density effect correction                    163   // density effect correction
203   G4double deltam;                                164   G4double deltam;
204   G4double x = std::log(bg2) / twoln10;        << 165   G4double x = log(bg2) / twoln10;
205   if ( x >= x0den ) {                             166   if ( x >= x0den ) {
206     deltam = twoln10 * x - cden;                  167     deltam = twoln10 * x - cden;
207     if ( x < x1den ) deltam += aden * std::pow << 168     if ( x < x1den ) deltam += aden * pow((x1den-x), mden);
208     dedx -= 0.5 * deltam;                         169     dedx -= 0.5 * deltam;
209   }                                               170   }
210                                                   171 
211   // now compute the total ionization loss        172   // now compute the total ionization loss
212   dedx *=  pi_hbarc2_over_mc2 * eDensity * nmp    173   dedx *=  pi_hbarc2_over_mc2 * eDensity * nmpl * nmpl;
213                                                   174 
214   if (dedx < 0.0) dedx = 0.;                   << 175   if (dedx < 0.0) dedx = 0;
215   return dedx;                                    176   return dedx;
216 }                                                 177 }
217                                                   178 
218 //....oooOO0OOooo........oooOO0OOooo........oo    179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
219                                                   180 
220 void G4mplIonisationModel::SampleSecondaries(s << 
221                const G4MaterialCutsCouple*,    << 
222                const G4DynamicParticle*,       << 
223                G4double,                       << 
224                G4double)                       << 
225 {}                                             << 
226                                                << 
227 //....oooOO0OOooo........oooOO0OOooo........oo << 
228                                                << 
229 G4double G4mplIonisationModel::SampleFluctuati    181 G4double G4mplIonisationModel::SampleFluctuations(
230                const G4MaterialCutsCouple* cou << 182                const G4Material* material,
231                const G4DynamicParticle* dp,       183                const G4DynamicParticle* dp,
232                                        const G << 184                G4double& tmax,
233                                        const G << 185                G4double& length,
234                const G4double length,          << 186                G4double& meanLoss)
235                const G4double meanLoss)        << 
236 {                                                 187 {
237   G4double siga = Dispersion(couple->GetMateri << 188   G4double siga = Dispersion(material,dp,tmax,length);
238   G4double loss = meanLoss;                       189   G4double loss = meanLoss;
239   siga = std::sqrt(siga);                      << 190   siga = sqrt(siga);
240   G4double twomeanLoss = meanLoss + meanLoss;     191   G4double twomeanLoss = meanLoss + meanLoss;
241                                                   192 
242   if(twomeanLoss < siga) {                        193   if(twomeanLoss < siga) {
243     G4double x;                                   194     G4double x;
244     do {                                          195     do {
245       loss = twomeanLoss*G4UniformRand();         196       loss = twomeanLoss*G4UniformRand();
246       x = (loss - meanLoss)/siga;                 197       x = (loss - meanLoss)/siga;
247       // Loop checking, 07-Aug-2015, Vladimir  << 
248     } while (1.0 - 0.5*x*x < G4UniformRand());    198     } while (1.0 - 0.5*x*x < G4UniformRand());
249   } else {                                        199   } else {
250     do {                                          200     do {
251       loss = G4RandGauss::shoot(meanLoss,siga)    201       loss = G4RandGauss::shoot(meanLoss,siga);
252       // Loop checking, 07-Aug-2015, Vladimir  << 
253     } while (0.0 > loss || loss > twomeanLoss)    202     } while (0.0 > loss || loss > twomeanLoss);
254   }                                               203   }
255   return loss;                                    204   return loss;
256 }                                                 205 }
257                                                << 
258 //....oooOO0OOooo........oooOO0OOooo........oo << 
259                                                << 
260 G4double G4mplIonisationModel::Dispersion(cons << 
261             const G4DynamicParticle* dp,       << 
262             const G4double tcut,               << 
263             const G4double tmax,               << 
264             const G4double length)             << 
265 {                                              << 
266   G4double siga = 0.0;                         << 
267   G4double tau   = dp->GetKineticEnergy()/mass << 
268   if(tau > 0.0) {                              << 
269     const G4double beta = dp->GetBeta();       << 
270     siga  = (tmax/(beta*beta) - 0.5*tcut) * tw << 
271       * material->GetElectronDensity() * charg << 
272   }                                            << 
273   return siga;                                 << 
274 }                                              << 
275                                                << 
276 //....oooOO0OOooo........oooOO0OOooo........oo << 
277                                                   206