ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4mplIonisationWithDeltaModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4mplIonisationWithDeltaModel.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: G4mplIonisationWithDeltaModel
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 "G4mplIonisationWithDeltaModel.hh"
56 #include "Randomize.hh"
57 #include "G4PhysicalConstants.hh"
58 #include "G4SystemOfUnits.hh"
59 #include "G4ParticleChangeForLoss.hh"
60 #include "G4Electron.hh"
61 #include "G4DynamicParticle.hh"
62 #include "G4ProductionCutsTable.hh"
63 #include "G4MaterialCutsCouple.hh"
64 #include "G4Log.hh"
65 #include "G4Pow.hh"
66 
67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
68 
69 using namespace std;
70 
71 std::vector<G4double>* G4mplIonisationWithDeltaModel::dedx0 = nullptr;
72 
73 G4mplIonisationWithDeltaModel::G4mplIonisationWithDeltaModel(G4double mCharge,
74  const G4String& nam)
75  : G4VEmModel(nam),G4VEmFluctuationModel(nam),
76  magCharge(mCharge),
77  twoln10(std::log(100.0)),
78  betalow(0.01),
79  betalim(0.1),
80  beta2lim(betalim*betalim),
81  bg2lim(beta2lim*(1.0 + beta2lim))
82 {
83  nmpl = G4lrint(std::abs(magCharge) * 2 * fine_structure_const);
84  if(nmpl > 6) { nmpl = 6; }
85  else if(nmpl < 1) { nmpl = 1; }
86  pi_hbarc2_over_mc2 = pi * hbarc * hbarc / electron_mass_c2;
87  chargeSquare = magCharge * magCharge;
88  dedxlim = 45.*nmpl*nmpl*GeV*cm2/g;
89  fParticleChange = nullptr;
90  theElectron = G4Electron::Electron();
91  G4cout << "### Monopole ionisation model with d-electron production, Gmag= "
92  << magCharge/eplus << G4endl;
93  monopole = nullptr;
94  mass = 0.0;
95 }
96 
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
98 
99 G4mplIonisationWithDeltaModel::~G4mplIonisationWithDeltaModel()
100 {
101  if(IsMaster()) { delete dedx0; }
102 }
103 
104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
105 
106 void G4mplIonisationWithDeltaModel::SetParticle(const G4ParticleDefinition* p)
107 {
108  monopole = p;
109  mass = monopole->GetPDGMass();
110  G4double emin =
111  std::min(LowEnergyLimit(),0.1*mass*(1./sqrt(1. - betalow*betalow) - 1.));
112  G4double emax =
113  std::max(HighEnergyLimit(),10*mass*(1./sqrt(1. - beta2lim) - 1.));
114  SetLowEnergyLimit(emin);
115  SetHighEnergyLimit(emax);
116 }
117 
118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
119 
120 void
121 G4mplIonisationWithDeltaModel::Initialise(const G4ParticleDefinition* p,
122  const G4DataVector&)
123 {
124  if(!monopole) { SetParticle(p); }
125  if(!fParticleChange) { fParticleChange = GetParticleChangeForLoss(); }
126  if(IsMaster()) {
127  if(!dedx0) { dedx0 = new std::vector<G4double>; }
128  G4ProductionCutsTable* theCoupleTable=
129  G4ProductionCutsTable::GetProductionCutsTable();
130  G4int numOfCouples = theCoupleTable->GetTableSize();
131  G4int n = dedx0->size();
132  if(n < numOfCouples) { dedx0->resize(numOfCouples); }
133  G4Pow* g4calc = G4Pow::GetInstance();
134 
135  // initialise vector assuming low conductivity
136  for(G4int i=0; i<numOfCouples; ++i) {
137 
138  const G4Material* material =
139  theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
140  G4double eDensity = material->GetElectronDensity();
141  G4double vF2 = 2*electron_Compton_length*g4calc->A13(3.*pi*pi*eDensity);
142  (*dedx0)[i] = pi_hbarc2_over_mc2*eDensity*nmpl*nmpl*
143  (G4Log(vF2/fine_structure_const) - 0.5)/vF2;
144  }
145  }
146 }
147 
148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
149 
150 G4double
151 G4mplIonisationWithDeltaModel::MinEnergyCut(const G4ParticleDefinition*,
152  const G4MaterialCutsCouple* couple)
153 {
154  return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy();
155 }
156 
157 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
158 
159 G4double
160 G4mplIonisationWithDeltaModel::ComputeDEDXPerVolume(const G4Material* material,
161  const G4ParticleDefinition* p,
162  G4double kineticEnergy,
163  G4double maxEnergy)
164 {
165  if(!monopole) { SetParticle(p); }
166  G4double tmax = MaxSecondaryEnergy(p,kineticEnergy);
167  G4double cutEnergy = std::min(tmax, maxEnergy);
168  cutEnergy = std::max(LowEnergyLimit(), cutEnergy);
169  G4double tau = kineticEnergy / mass;
170  G4double gam = tau + 1.0;
171  G4double bg2 = tau * (tau + 2.0);
172  G4double beta2 = bg2 / (gam * gam);
173  G4double beta = sqrt(beta2);
174 
175  // low-energy asymptotic formula
176  G4double dedx = (*dedx0)[CurrentCouple()->GetIndex()]*beta;
177 
178  // above asymptotic
179  if(beta > betalow) {
180 
181  // high energy
182  if(beta >= betalim) {
183  dedx = ComputeDEDXAhlen(material, bg2, cutEnergy);
184 
185  } else {
186  G4double dedx1 = (*dedx0)[CurrentCouple()->GetIndex()]*betalow;
187  G4double dedx2 = ComputeDEDXAhlen(material, bg2lim, cutEnergy);
188 
189  // extrapolation between two formula
190  G4double kapa2 = beta - betalow;
191  G4double kapa1 = betalim - beta;
192  dedx = (kapa1*dedx1 + kapa2*dedx2)/(kapa1 + kapa2);
193  }
194  }
195  return dedx;
196 }
197 
198 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
199 
200 G4double
201 G4mplIonisationWithDeltaModel::ComputeDEDXAhlen(const G4Material* material,
202  G4double bg2,
203  G4double cutEnergy)
204 {
205  G4double eDensity = material->GetElectronDensity();
206  G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
207 
208  // Ahlen's formula for nonconductors, [1]p157, f(5.7)
209  G4double dedx =
210  0.5*(G4Log(2.0*electron_mass_c2*bg2*cutEnergy/(eexc*eexc)) -1.0);
211 
212  // Kazama et al. cross-section correction
213  G4double k = 0.406;
214  if(nmpl > 1) { k = 0.346; }
215 
216  // Bloch correction
217  const G4double B[7] = { 0.0, 0.248, 0.672, 1.022, 1.243, 1.464, 1.685};
218 
219  dedx += 0.5 * k - B[nmpl];
220 
221  // density effect correction
222  G4double x = G4Log(bg2)/twoln10;
223  dedx -= material->GetIonisation()->DensityCorrection(x);
224 
225  // now compute the total ionization loss
226  dedx *= pi_hbarc2_over_mc2 * eDensity * nmpl * nmpl;
227 
228  dedx = std::max(dedx, 0.0);
229  return dedx;
230 }
231 
232 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
233 
234 G4double
235 G4mplIonisationWithDeltaModel::ComputeCrossSectionPerElectron(
236  const G4ParticleDefinition* p,
237  G4double kineticEnergy,
238  G4double cut,
239  G4double maxKinEnergy)
240 {
241  if(!monopole) { SetParticle(p); }
242  G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
243  G4double maxEnergy = std::min(tmax, maxKinEnergy);
244  G4double cutEnergy = std::max(LowEnergyLimit(), cut);
245  G4double cross = (cutEnergy < maxEnergy)
246  ? (0.5/cutEnergy - 0.5/maxEnergy)*pi_hbarc2_over_mc2 * nmpl * nmpl : 0.0;
247  return cross;
248 }
249 
250 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
251 
252 G4double
253 G4mplIonisationWithDeltaModel::ComputeCrossSectionPerAtom(
254  const G4ParticleDefinition* p,
255  G4double kineticEnergy,
256  G4double Z, G4double,
257  G4double cutEnergy,
258  G4double maxEnergy)
259 {
260  G4double cross =
261  Z*ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
262  return cross;
263 }
264 
265 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
266 
267 void
268 G4mplIonisationWithDeltaModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
269  const G4MaterialCutsCouple*,
270  const G4DynamicParticle* dp,
271  G4double minKinEnergy,
272  G4double maxEnergy)
273 {
274  G4double kineticEnergy = dp->GetKineticEnergy();
275  G4double tmax = MaxSecondaryEnergy(dp->GetDefinition(),kineticEnergy);
276 
277  G4double maxKinEnergy = std::min(maxEnergy,tmax);
278  if(minKinEnergy >= maxKinEnergy) { return; }
279 
280  //G4cout << "G4mplIonisationWithDeltaModel::SampleSecondaries: E(GeV)= "
281  // << kineticEnergy/GeV << " M(GeV)= " << mass/GeV
282  // << " tmin(MeV)= " << minKinEnergy/MeV << G4endl;
283 
284  G4double totEnergy = kineticEnergy + mass;
285  G4double etot2 = totEnergy*totEnergy;
286  G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/etot2;
287 
288  // sampling without nuclear size effect
289  G4double q = G4UniformRand();
290  G4double deltaKinEnergy = minKinEnergy*maxKinEnergy
291  /(minKinEnergy*(1.0 - q) + maxKinEnergy*q);
292 
293  // delta-electron is produced
294  G4double totMomentum = totEnergy*sqrt(beta2);
295  G4double deltaMomentum =
296  sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
297  G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) /
298  (deltaMomentum * totMomentum);
299  cost = std::min(cost, 1.0);
300 
301  G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
302 
303  G4double phi = twopi * G4UniformRand() ;
304 
305  G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost);
306  G4ThreeVector direction = dp->GetMomentumDirection();
307  deltaDirection.rotateUz(direction);
308 
309  // create G4DynamicParticle object for delta ray
310  G4DynamicParticle* delta =
311  new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
312 
313  vdp->push_back(delta);
314 
315  // Change kinematics of primary particle
316  kineticEnergy -= deltaKinEnergy;
317  G4ThreeVector finalP = direction*totMomentum - deltaDirection*deltaMomentum;
318  finalP = finalP.unit();
319 
320  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
321  fParticleChange->SetProposedMomentumDirection(finalP);
322 }
323 
324 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
325 
326 G4double G4mplIonisationWithDeltaModel::SampleFluctuations(
327  const G4MaterialCutsCouple* couple,
328  const G4DynamicParticle* dp,
329  G4double tmax,
330  G4double length,
331  G4double meanLoss)
332 {
333  G4double siga = Dispersion(couple->GetMaterial(),dp,tmax,length);
334  G4double loss = meanLoss;
335  siga = sqrt(siga);
336  G4double twomeanLoss = meanLoss + meanLoss;
337 
338  if(twomeanLoss < siga) {
339  G4double x;
340  do {
341  loss = twomeanLoss*G4UniformRand();
342  x = (loss - meanLoss)/siga;
343  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
344  } while (1.0 - 0.5*x*x < G4UniformRand());
345  } else {
346  do {
347  loss = G4RandGauss::shoot(meanLoss,siga);
348  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
349  } while (0.0 > loss || loss > twomeanLoss);
350  }
351  return loss;
352 }
353 
354 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
355 
356 G4double
357 G4mplIonisationWithDeltaModel::Dispersion(const G4Material* material,
358  const G4DynamicParticle* dp,
359  G4double tmax,
360  G4double length)
361 {
362  G4double siga = 0.0;
363  G4double tau = dp->GetKineticEnergy()/mass;
364  if(tau > 0.0) {
365  G4double electronDensity = material->GetElectronDensity();
366  G4double gam = tau + 1.0;
367  G4double invbeta2 = (gam*gam)/(tau * (tau+2.0));
368  siga = (invbeta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
369  * electronDensity * chargeSquare;
370  }
371  return siga;
372 }
373 
374 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
375 
376 G4double
377 G4mplIonisationWithDeltaModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
378  G4double kinEnergy)
379 {
380  G4double tau = kinEnergy/mass;
381  return 2.0*electron_mass_c2*tau*(tau + 2.);
382 }
383 
384 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....