ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4mplIonisationModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4mplIonisationModel.cc
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 using namespace std;
68 
69 std::vector<G4double>* G4mplIonisationModel::dedx0 = nullptr;
70 
71 G4mplIonisationModel::G4mplIonisationModel(G4double mCharge, const G4String& nam)
72  : G4VEmModel(nam),G4VEmFluctuationModel(nam),
73  magCharge(mCharge),
74  twoln10(log(100.0)),
75  betalow(0.01),
76  betalim(0.1),
77  beta2lim(betalim*betalim),
78  bg2lim(beta2lim*(1.0 + beta2lim))
79 {
80  nmpl = G4int(abs(magCharge) * 2 * fine_structure_const + 0.5);
81  if(nmpl > 6) { nmpl = 6; }
82  else if(nmpl < 1) { nmpl = 1; }
83  pi_hbarc2_over_mc2 = pi * hbarc * hbarc / electron_mass_c2;
84  chargeSquare = magCharge * magCharge;
85  dedxlim = 45.*nmpl*nmpl*GeV*cm2/g;
86  fParticleChange = nullptr;
87  monopole = nullptr;
88  mass = 0.0;
89 }
90 
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92 
93 G4mplIonisationModel::~G4mplIonisationModel()
94 {
95  if(IsMaster()) { delete dedx0; }
96 }
97 
98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
99 
100 void G4mplIonisationModel::SetParticle(const G4ParticleDefinition* p)
101 {
102  monopole = p;
103  mass = monopole->GetPDGMass();
104  G4double emin =
105  std::min(LowEnergyLimit(),0.1*mass*(1./sqrt(1. - betalow*betalow) - 1.));
106  G4double emax =
107  std::max(HighEnergyLimit(),10.*mass*(1./sqrt(1. - beta2lim) - 1.));
108  SetLowEnergyLimit(emin);
109  SetHighEnergyLimit(emax);
110 }
111 
112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
113 
114 void G4mplIonisationModel::Initialise(const G4ParticleDefinition* p,
115  const G4DataVector&)
116 {
117  if(!monopole) { SetParticle(p); }
118  if(!fParticleChange) { fParticleChange = GetParticleChangeForLoss(); }
119  if(IsMaster()) {
120  if(!dedx0) { dedx0 = new std::vector<G4double>; }
121  G4ProductionCutsTable* theCoupleTable=
122  G4ProductionCutsTable::GetProductionCutsTable();
123  G4int numOfCouples = theCoupleTable->GetTableSize();
124  G4int n = dedx0->size();
125  if(n < numOfCouples) { dedx0->resize(numOfCouples); }
126 
127  G4Pow* g4calc = G4Pow::GetInstance();
128 
129  // initialise vector assuming low conductivity
130  for(G4int i=0; i<numOfCouples; ++i) {
131 
132  const G4Material* material =
133  theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
134  G4double eDensity = material->GetElectronDensity();
135  G4double vF2 = 2*electron_Compton_length*g4calc->A13(3.*pi*pi*eDensity);
136  (*dedx0)[i] = pi_hbarc2_over_mc2*eDensity*nmpl*nmpl*
137  (G4Log(vF2/fine_structure_const) - 0.5)/vF2;
138  }
139  }
140 }
141 
142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
143 
144 G4double G4mplIonisationModel::ComputeDEDXPerVolume(const G4Material* material,
145  const G4ParticleDefinition* p,
146  G4double kineticEnergy,
147  G4double)
148 {
149  if(!monopole) { SetParticle(p); }
150  G4double tau = kineticEnergy / mass;
151  G4double gam = tau + 1.0;
152  G4double bg2 = tau * (tau + 2.0);
153  G4double beta2 = bg2 / (gam * gam);
154  G4double beta = sqrt(beta2);
155 
156  // low-energy asymptotic formula
157  //G4double dedx = dedxlim*beta*material->GetDensity();
158  G4double dedx = (*dedx0)[CurrentCouple()->GetIndex()]*beta;
159 
160  // above asymptotic
161  if(beta > betalow) {
162 
163  // high energy
164  if(beta >= betalim) {
165  dedx = ComputeDEDXAhlen(material, bg2);
166 
167  } else {
168 
169  //G4double dedx1 = dedxlim*betalow*material->GetDensity();
170  G4double dedx1 = (*dedx0)[CurrentCouple()->GetIndex()]*betalow;
171  G4double dedx2 = ComputeDEDXAhlen(material, bg2lim);
172 
173  // extrapolation between two formula
174  G4double kapa2 = beta - betalow;
175  G4double kapa1 = betalim - beta;
176  dedx = (kapa1*dedx1 + kapa2*dedx2)/(kapa1 + kapa2);
177  }
178  }
179  return dedx;
180 }
181 
182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
183 
184 G4double G4mplIonisationModel::ComputeDEDXAhlen(const G4Material* material,
185  G4double bg2)
186 {
187  G4double eDensity = material->GetElectronDensity();
188  G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
189  G4double cden = material->GetIonisation()->GetCdensity();
190  G4double mden = material->GetIonisation()->GetMdensity();
191  G4double aden = material->GetIonisation()->GetAdensity();
192  G4double x0den = material->GetIonisation()->GetX0density();
193  G4double x1den = material->GetIonisation()->GetX1density();
194 
195  // Ahlen's formula for nonconductors, [1]p157, f(5.7)
196  G4double dedx = log(2.0 * electron_mass_c2 * bg2 / eexc) - 0.5;
197 
198  // Kazama et al. cross-section correction
199  G4double k = 0.406;
200  if(nmpl > 1) k = 0.346;
201 
202  // Bloch correction
203  const G4double B[7] = { 0.0, 0.248, 0.672, 1.022, 1.243, 1.464, 1.685};
204 
205  dedx += 0.5 * k - B[nmpl];
206 
207  // density effect correction
208  G4double deltam;
209  G4double x = log(bg2) / twoln10;
210  if ( x >= x0den ) {
211  deltam = twoln10 * x - cden;
212  if ( x < x1den ) deltam += aden * pow((x1den-x), mden);
213  dedx -= 0.5 * deltam;
214  }
215 
216  // now compute the total ionization loss
217  dedx *= pi_hbarc2_over_mc2 * eDensity * nmpl * nmpl;
218 
219  if (dedx < 0.0) dedx = 0.;
220  return dedx;
221 }
222 
223 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
224 
225 void G4mplIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>*,
226  const G4MaterialCutsCouple*,
227  const G4DynamicParticle*,
228  G4double,
229  G4double)
230 {}
231 
232 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
233 
234 G4double G4mplIonisationModel::SampleFluctuations(
235  const G4MaterialCutsCouple* couple,
236  const G4DynamicParticle* dp,
237  G4double tmax,
238  G4double length,
239  G4double meanLoss)
240 {
241  G4double siga = Dispersion(couple->GetMaterial(),dp,tmax,length);
242  G4double loss = meanLoss;
243  siga = sqrt(siga);
244  G4double twomeanLoss = meanLoss + meanLoss;
245 
246  if(twomeanLoss < siga) {
247  G4double x;
248  do {
249  loss = twomeanLoss*G4UniformRand();
250  x = (loss - meanLoss)/siga;
251  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
252  } while (1.0 - 0.5*x*x < G4UniformRand());
253  } else {
254  do {
255  loss = G4RandGauss::shoot(meanLoss,siga);
256  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
257  } while (0.0 > loss || loss > twomeanLoss);
258  }
259  return loss;
260 }
261 
262 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
263 
264 G4double G4mplIonisationModel::Dispersion(const G4Material* material,
265  const G4DynamicParticle* dp,
266  G4double tmax,
267  G4double length)
268 {
269  G4double siga = 0.0;
270  G4double tau = dp->GetKineticEnergy()/mass;
271  if(tau > 0.0) {
272  G4double electronDensity = material->GetElectronDensity();
273  G4double gam = tau + 1.0;
274  G4double invbeta2 = (gam*gam)/(tau * (tau+2.0));
275  siga = (invbeta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
276  * electronDensity * chargeSquare;
277  }
278  return siga;
279 }
280 
281 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....