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 8.3)


  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.3 2006/12/13 15:44:27 gunter Exp $
                                                   >>  27 // GEANT4 tag $Name: geant4-08-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 << 
 40 //            Small bug fixing and refactoring << 
 41 // 13.11.2007 Use low-energy asymptotic from [ << 
 42 //                                                 41 //
 43 //                                                 42 //
 44 // -------------------------------------------     43 // -------------------------------------------------------------------
 45 // References                                  <<  44 //
 46 // [1] Steven P. Ahlen: Energy loss of relativ << 
 47 //     S.P. Ahlen, Rev. Mod. Phys 52(1980), p1 << 
 48 // [2] K.A. Milton arXiv:hep-ex/0602040        << 
 49 // [3] S.P. Ahlen and K. Kinoshita, Phys. Rev. << 
 50                                                    45 
 51                                                    46 
 52 //....oooOO0OOooo........oooOO0OOooo........oo     47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 53 //....oooOO0OOooo........oooOO0OOooo........oo     48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 54                                                    49 
 55 #include "G4mplIonisationModel.hh"                 50 #include "G4mplIonisationModel.hh"
 56 #include "Randomize.hh"                            51 #include "Randomize.hh"
 57 #include "G4PhysicalConstants.hh"              <<  52 #include "G4LossTableManager.hh"
 58 #include "G4SystemOfUnits.hh"                  << 
 59 #include "G4ParticleChangeForLoss.hh"              53 #include "G4ParticleChangeForLoss.hh"
 60 #include "G4ProductionCutsTable.hh"            << 
 61 #include "G4MaterialCutsCouple.hh"             << 
 62 #include "G4Log.hh"                            << 
 63 #include "G4Pow.hh"                            << 
 64                                                    54 
 65 //....oooOO0OOooo........oooOO0OOooo........oo     55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 66                                                    56 
 67 std::vector<G4double>* G4mplIonisationModel::d <<  57 using namespace std;
 68                                                    58 
 69 G4mplIonisationModel::G4mplIonisationModel(G4d <<  59 G4mplIonisationModel::G4mplIonisationModel(G4double mCharge, 
                                                   >>  60              const G4String& nam)
 70   : G4VEmModel(nam),G4VEmFluctuationModel(nam)     61   : G4VEmModel(nam),G4VEmFluctuationModel(nam),
 71   magCharge(mCharge),                              62   magCharge(mCharge),
 72   twoln10(G4Log(100.0)),                       <<  63   twoln10(2.0*log(10.0)),
 73   betalow(0.01),                               <<  64   beta2low(0.0001),
 74   betalim(0.1),                                <<  65   beta2lim(0.01),
 75   beta2lim(betalim*betalim),                   << 
 76   bg2lim(beta2lim*(1.0 + beta2lim))                66   bg2lim(beta2lim*(1.0 + beta2lim))
 77 {                                                  67 {
 78   nmpl = G4int(std::abs(magCharge) * 2 * CLHEP <<  68   nmpl         = G4int(abs(magCharge)/68.0);
 79   if(nmpl > 6)      { nmpl = 6; }              <<  69   if(nmpl > 6)      nmpl = 6;
 80   else if(nmpl < 1) { nmpl = 1; }              <<  70   else if(nmpl < 1) nmpl = 1;
 81   pi_hbarc2_over_mc2 = CLHEP::pi*CLHEP::hbarc* <<  71   G4double x   = 45.0*GeV*G4double(nmpl)/cm;
 82   chargeSquare = magCharge * magCharge;        <<  72   factlow      = x*x;
 83   dedxlim = 45.*nmpl*nmpl*CLHEP::GeV*CLHEP::cm <<  73   chargeSquare = magCharge*magCharge;
 84 }                                                  74 }
 85                                                    75 
 86 //....oooOO0OOooo........oooOO0OOooo........oo     76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 87                                                    77 
 88 G4mplIonisationModel::~G4mplIonisationModel()      78 G4mplIonisationModel::~G4mplIonisationModel()
 89 {                                              <<  79 {}
 90   if(IsMaster()) { delete dedx0; }             << 
 91 }                                              << 
 92                                                    80 
 93 //....oooOO0OOooo........oooOO0OOooo........oo     81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 94                                                    82 
 95 void G4mplIonisationModel::SetParticle(const G <<  83 void G4mplIonisationModel::Initialise(const G4ParticleDefinition* p,
                                                   >>  84               const G4DataVector&)
 96 {                                                  85 {
 97   monopole = p;                                    86   monopole = p;
 98   mass     = monopole->GetPDGMass();               87   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                                                << 
107 //....oooOO0OOooo........oooOO0OOooo........oo << 
108                                                    88 
109 void G4mplIonisationModel::Initialise(const G4 <<  89   if(pParticleChange) 
110               const G4DataVector&)             <<  90     fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
111 {                                              <<  91   else 
112   if(nullptr == monopole) { SetParticle(p); }  <<  92     fParticleChange = new G4ParticleChangeForLoss();
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 }                                                  93 }
136                                                    94 
137 //....oooOO0OOooo........oooOO0OOooo........oo     95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
138                                                    96 
139 G4double G4mplIonisationModel::ComputeDEDXPerV     97 G4double G4mplIonisationModel::ComputeDEDXPerVolume(const G4Material* material,
140                 const G4ParticleDefinition* p, <<  98                 const G4ParticleDefinition*,
141                 G4double kineticEnergy,            99                 G4double kineticEnergy,
142                 G4double)                         100                 G4double)
143 {                                                 101 {
144   if(nullptr == monopole) { SetParticle(p); }  << 102   G4double tau   = kineticEnergy/mass;
145   G4double tau   = kineticEnergy / mass;       << 
146   G4double gam   = tau + 1.0;                     103   G4double gam   = tau + 1.0;
147   G4double bg2   = tau * (tau + 2.0);          << 104   G4double bg2   = tau * (tau+2.0);
148   G4double beta2 = bg2 / (gam * gam);          << 105   G4double beta2 = bg2/(gam*gam);
149   G4double beta  = std::sqrt(beta2);           << 106 
150                                                << 107   G4double dedx0 = factlow*abs(beta2);
151   // low-energy asymptotic formula             << 108 
152   //G4double dedx  = dedxlim*beta*material->Ge << 109   if(beta2 > beta2low) {
153   G4double dedx = (*dedx0)[CurrentCouple()->Ge << 110 
154                                                << 111     G4double b2 = beta2;
155   // above asymptotic                          << 112     if(beta2 < beta2lim) {
156   if(beta > betalow) {                         << 113       beta2= beta2lim;
157                                                << 114       bg2  = bg2lim;
158     // high energy                             << 
159     if(beta >= betalim) {                      << 
160       dedx = ComputeDEDXAhlen(material, bg2);  << 
161                                                << 
162     } else {                                   << 
163                                                << 
164       //G4double dedx1 = dedxlim*betalow*mater << 
165       G4double dedx1 = (*dedx0)[CurrentCouple( << 
166       G4double dedx2 = ComputeDEDXAhlen(materi << 
167                                                << 
168       // extrapolation between two formula     << 
169       G4double kapa2 = beta - betalow;         << 
170       G4double kapa1 = betalim - beta;         << 
171       dedx = (kapa1*dedx1 + kapa2*dedx2)/(kapa << 
172     }                                             115     }
173   }                                            << 
174   return dedx;                                 << 
175 }                                              << 
176                                                   116 
177 //....oooOO0OOooo........oooOO0OOooo........oo << 117     G4double eexc  = material->GetIonisation()->GetMeanExcitationEnergy();
                                                   >> 118     G4double cden  = material->GetIonisation()->GetCdensity();
                                                   >> 119     G4double mden  = material->GetIonisation()->GetMdensity();
                                                   >> 120     G4double aden  = material->GetIonisation()->GetAdensity();
                                                   >> 121     G4double x0den = material->GetIonisation()->GetX0density();
                                                   >> 122     G4double x1den = material->GetIonisation()->GetX1density();
                                                   >> 123 
                                                   >> 124     G4double eDensity = material->GetElectronDensity();
                                                   >> 125 
                                                   >> 126     G4double dedx = 2.0*log(2.0*electron_mass_c2*bg2/eexc) - 1.0;
                                                   >> 127 
                                                   >> 128     G4double  k = 0.406;
                                                   >> 129     if(nmpl > 1) k = 0.346;
                                                   >> 130     const G4double B[7] = { 0.0, 0.248, 0.672, 1.022,  1.243, 1.464,  1.685}; 
                                                   >> 131 
                                                   >> 132     dedx += k - B[nmpl];
                                                   >> 133 
                                                   >> 134     // density correction
                                                   >> 135     G4double x = log(bg2)/twoln10;
                                                   >> 136     if ( x >= x0den ) {
                                                   >> 137       dedx -= twoln10*x - cden ;
                                                   >> 138       if ( x < x1den ) dedx -= aden*pow((x1den-x),mden) ;
                                                   >> 139     }
178                                                   140 
179 G4double G4mplIonisationModel::ComputeDEDXAhle << 141     // now compute the total ionization loss
180             G4double bg2)                      << 
181 {                                              << 
182   G4double eDensity = material->GetElectronDen << 
183   G4double eexc  = material->GetIonisation()-> << 
184   G4double cden  = material->GetIonisation()-> << 
185   G4double mden  = material->GetIonisation()-> << 
186   G4double aden  = material->GetIonisation()-> << 
187   G4double x0den = material->GetIonisation()-> << 
188   G4double x1den = material->GetIonisation()-> << 
189                                                << 
190   // Ahlen's formula for nonconductors, [1]p15 << 
191   G4double dedx = std::log(2.0 * electron_mass << 
192                                                << 
193   // Kazama et al. cross-section correction    << 
194   G4double  k = 0.406;                         << 
195   if(nmpl > 1) k = 0.346;                      << 
196                                                << 
197   // Bloch correction                          << 
198   const G4double B[7] = { 0.0, 0.248, 0.672, 1 << 
199                                                << 
200   dedx += 0.5 * k - B[nmpl];                   << 
201                                                << 
202   // density effect correction                 << 
203   G4double deltam;                             << 
204   G4double x = std::log(bg2) / twoln10;        << 
205   if ( x >= x0den ) {                          << 
206     deltam = twoln10 * x - cden;               << 
207     if ( x < x1den ) deltam += aden * std::pow << 
208     dedx -= 0.5 * deltam;                      << 
209   }                                            << 
210                                                   142 
211   // now compute the total ionization loss     << 143     if (dedx < 0.0) dedx = 0.0 ;
212   dedx *=  pi_hbarc2_over_mc2 * eDensity * nmp << 
213                                                   144 
214   if (dedx < 0.0) dedx = 0.;                   << 145     dedx *= twopi_mc2_rcl2*chargeSquare*eDensity;
215   return dedx;                                 << 
216 }                                              << 
217                                                   146 
218 //....oooOO0OOooo........oooOO0OOooo........oo << 147     // extrapolate between two formula
                                                   >> 148     if(beta2 < beta2lim) {
                                                   >> 149       x = log(dedx0) + log(dedx/dedx0)*log(b2/beta2low)/log(beta2lim/beta2low);
                                                   >> 150       dedx = exp(x);
                                                   >> 151     }
                                                   >> 152     dedx0 = dedx;
                                                   >> 153   }
219                                                   154 
220 void G4mplIonisationModel::SampleSecondaries(s << 155   return dedx0;
221                const G4MaterialCutsCouple*,    << 156 }
222                const G4DynamicParticle*,       << 
223                G4double,                       << 
224                G4double)                       << 
225 {}                                             << 
226                                                   157 
227 //....oooOO0OOooo........oooOO0OOooo........oo    158 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
228                                                   159 
229 G4double G4mplIonisationModel::SampleFluctuati    160 G4double G4mplIonisationModel::SampleFluctuations(
230                const G4MaterialCutsCouple* cou << 161                const G4Material* material,
231                const G4DynamicParticle* dp,       162                const G4DynamicParticle* dp,
232                                        const G << 163                G4double& tmax,
233                                        const G << 164                G4double& length,
234                const G4double length,          << 165                G4double& meanLoss)
235                const G4double meanLoss)        << 
236 {                                                 166 {
237   G4double siga = Dispersion(couple->GetMateri << 167   G4double siga = Dispersion(material,dp,tmax,length);
238   G4double loss = meanLoss;                       168   G4double loss = meanLoss;
239   siga = std::sqrt(siga);                      << 169   siga = sqrt(siga);
240   G4double twomeanLoss = meanLoss + meanLoss;     170   G4double twomeanLoss = meanLoss + meanLoss;
241                                                   171 
242   if(twomeanLoss < siga) {                        172   if(twomeanLoss < siga) {
243     G4double x;                                   173     G4double x;
244     do {                                          174     do {
245       loss = twomeanLoss*G4UniformRand();         175       loss = twomeanLoss*G4UniformRand();
246       x = (loss - meanLoss)/siga;                 176       x = (loss - meanLoss)/siga;
247       // Loop checking, 07-Aug-2015, Vladimir  << 
248     } while (1.0 - 0.5*x*x < G4UniformRand());    177     } while (1.0 - 0.5*x*x < G4UniformRand());
249   } else {                                        178   } else {
250     do {                                          179     do {
251       loss = G4RandGauss::shoot(meanLoss,siga)    180       loss = G4RandGauss::shoot(meanLoss,siga);
252       // Loop checking, 07-Aug-2015, Vladimir  << 
253     } while (0.0 > loss || loss > twomeanLoss)    181     } while (0.0 > loss || loss > twomeanLoss);
254   }                                               182   }
                                                   >> 183   //G4cout << "G4mplIonisationModel::SampleFluctuations:  loss= " << loss 
                                                   >> 184   //<< "  siga= " << siga << G4endl;
255   return loss;                                    185   return loss;
256 }                                              << 
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 }                                                 186 }
275                                                   187 
276 //....oooOO0OOooo........oooOO0OOooo........oo    188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
277                                                   189