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 ]

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