ATLAS Offline Software
G4mplAtlasIonisationWithDeltaModel.cxx
Go to the documentation of this file.
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 // $Id: G4mplAtlasIonisationWithDeltaModel.cxx 729653 2016-03-14 15:55:47Z jchapman $
27 // GEANT4 tag $Name: geant4-09-03-ref-00 $
28 //
29 // -------------------------------------------------------------------
30 //
31 // GEANT4 Class header file
32 //
33 //
34 // File name: G4mplAtlasIonisationWithDeltaModel
35 //
36 // Author: Vladimir Ivanchenko
37 //
38 // Creation date: 06.09.2005
39 //
40 // Modifications:
41 // 12.08.2007 Changing low energy approximation and extrapolation.
42 // Small bug fixing and refactoring (M. Vladymyrov)
43 // 13.11.2007 Use low-energy asymptotic from [3] (V.Ivanchenko)
44 // 2014-03-28 Allow for fractional magnetic charges (G. Palacino)
45 //
46 //
47 // -------------------------------------------------------------------
48 // References
49 // [1] Steven P. Ahlen: Energy loss of relativistic heavy ionizing particles,
50 // S.P. Ahlen, Rev. Mod. Phys 52(1980), p121
51 // [2] K.A. Milton arXiv:hep-ex/0602040
52 // [3] S.P. Ahlen and K. Kinoshita, Phys. Rev. D26 (1982) 2347
53 
54 
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
57 
58 // class header
59 #include "G4mplAtlasIonisationWithDeltaModel.hh"
60 // Geant4 headers
61 #include "Randomize.hh"
62 #include "G4LossTableManager.hh"
63 #include "G4ParticleChangeForLoss.hh"
64 #include "G4Electron.hh"
65 #include "G4DynamicParticle.hh"
66 #include "G4Version.hh"
67 // CLHEP headers
68 #include "CLHEP/Units/SystemOfUnits.h"
69 #include "CLHEP/Units/PhysicalConstants.h"
70 
71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
72 
73 G4mplAtlasIonisationWithDeltaModel::G4mplAtlasIonisationWithDeltaModel(G4double mCharge, const G4String& nam)
74  : G4VEmModel(nam),G4VEmFluctuationModel(nam),
75  monopole(0),
76  mass(0.0),
77  magCharge(mCharge),
78  twoln10(std::log(100.0)),
79  betalow(0.01),
80  betalim(0.1),
81  beta2lim(betalim*betalim),
82  bg2lim(beta2lim*(1.0 + beta2lim))
83 {
84  nmpl = std::abs(magCharge) * 2 * CLHEP::fine_structure_const;
85  //if(nmpl > 6) { nmpl = 6; }
86  //else if(nmpl < 1) { nmpl = 1; }
87  pi_hbarc2_over_mc2 = CLHEP::pi * CLHEP::hbarc * CLHEP::hbarc / CLHEP::electron_mass_c2;
88  chargeSquare = magCharge * magCharge;
89  dedxlim = 45.*nmpl*nmpl*CLHEP::GeV*CLHEP::cm2/CLHEP::g;
90  fParticleChange = 0;
91  theElectron = G4Electron::Electron();
92  G4cout << "### Monopole ionisation model with d-electron production, Gmag= "
93  << magCharge/CLHEP::eplus << G4endl;
94 }
95 
96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
97 
98 G4mplAtlasIonisationWithDeltaModel::~G4mplAtlasIonisationWithDeltaModel()
99 {}
100 
101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
102 
103 void
104 G4mplAtlasIonisationWithDeltaModel::Initialise(const G4ParticleDefinition* p,
105  const G4DataVector&)
106 {
107  monopole = p;
108  mass = monopole->GetPDGMass();
109  // if(!fParticleChange) { fParticleChange = GetParticleChangeForLoss(); }
110 
111  if(pParticleChange)
112  fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
113  else
114  fParticleChange = new G4ParticleChangeForLoss();
115 
116 }
117 
118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
119 
120 G4double
121 G4mplAtlasIonisationWithDeltaModel::ComputeDEDXPerVolume(const G4Material* material,
122  const G4ParticleDefinition* p,
123  G4double kineticEnergy,
124  G4double maxEnergy)
125 {
126  const G4double tmax = MaxSecondaryEnergy(p,kineticEnergy);
127  const G4double cutEnergy = std::min(tmax, maxEnergy);
128  const G4double tau = kineticEnergy / mass;
129  const G4double gam = tau + 1.0;
130  const G4double bg2 = tau * (tau + 2.0);
131  const G4double beta2 = bg2 / (gam * gam);
132  const G4double beta = std::sqrt(beta2);
133 
134  // low-energy asymptotic formula
135  G4double dedx = dedxlim*beta*material->GetDensity();
136 
137  // above asymptotic
138  if(beta > betalow) {
139 
140  // high energy
141  if(beta >= betalim) {
142  dedx = ComputeDEDXAhlen(material, bg2, cutEnergy);
143 
144  } else {
145 
146  const G4double dedx1 = dedxlim*betalow*material->GetDensity();
147  const G4double dedx2 = ComputeDEDXAhlen(material, bg2lim, cutEnergy);
148 
149  // extrapolation between two formula
150  const G4double kapa2 = beta - betalow;
151  const G4double kapa1 = betalim - beta;
152  dedx = (kapa1*dedx1 + kapa2*dedx2)/(kapa1 + kapa2);
153  }
154  }
155  return dedx;
156 }
157 
158 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
159 
160 G4double
161 G4mplAtlasIonisationWithDeltaModel::ComputeDEDXAhlen(const G4Material* material,
162  G4double bg2,
163  G4double cutEnergy)
164 {
165  const G4double eDensity = material->GetElectronDensity();
166  const G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
167 
168  // Ahlen's formula for nonconductors, [1]p157, f(5.7)
169  G4double dedx =
170  0.5*(std::log(2.0 * CLHEP::electron_mass_c2 * bg2*cutEnergy / (eexc*eexc)) - 1.0);
171 
172  // Kazama et al. cross-section correction
173  G4double k = 0.406;
174  if(nmpl > 1) { k = 0.346; }
175 
176  // Bloch correction
177  const G4double B[7] = { 0.0, 0.248, 0.672, 1.022, 1.243, 1.464, 1.685};
178 
179  if(nmpl > 6.0) { dedx += 0.5 * k - B[6]; }
180  else if(nmpl < 1.0) { dedx += 0.5 * k - B[1]; }
181  else {
182  dedx += 0.5 * k - B[G4int(std::abs(magCharge) * 2 * CLHEP::fine_structure_const + 0.5)];
183  }
184  // density effect correction
185  const G4double x = std::log(bg2)/twoln10;
186 
187  dedx -= material->GetIonisation()->DensityCorrection(x);
188 
189  // now compute the total ionization loss
190  dedx *= pi_hbarc2_over_mc2 * eDensity * nmpl * nmpl;
191 
192  if (dedx < 0.0) { dedx = 0; }
193  return dedx;
194 }
195 
196 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
197 
198 G4double
199 G4mplAtlasIonisationWithDeltaModel::ComputeCrossSectionPerElectron(
200  const G4ParticleDefinition* p,
201  G4double kineticEnergy,
202  G4double cutEnergy,
203  G4double maxKinEnergy)
204 {
205  G4double cross = 0.0;
206  const G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
207  const G4double maxEnergy = std::min(tmax,maxKinEnergy);
208  if(cutEnergy < maxEnergy) {
209  cross = (1.0/cutEnergy - 1.0/maxEnergy)*CLHEP::twopi_mc2_rcl2*chargeSquare;
210  }
211  return cross;
212 }
213 
214 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
215 
216 G4double
217 G4mplAtlasIonisationWithDeltaModel::ComputeCrossSectionPerAtom(
218  const G4ParticleDefinition* p,
219  G4double kineticEnergy,
220  G4double Z, G4double,
221  G4double cutEnergy,
222  G4double maxEnergy)
223 {
224  const G4double cross =
225  Z*ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
226  return cross;
227 }
228 
229 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
230 
231 void
232 G4mplAtlasIonisationWithDeltaModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
233  const G4MaterialCutsCouple*,
234  const G4DynamicParticle* dp,
235  G4double minKinEnergy,
236  G4double maxEnergy)
237 {
238  G4double kineticEnergy = dp->GetKineticEnergy();
239  const G4double tmax = MaxSecondaryEnergy(dp->GetDefinition(),kineticEnergy);
240 
241  G4double maxKinEnergy = std::min(maxEnergy,tmax);
242  if(minKinEnergy >= maxKinEnergy) { return; }
243 
244  //G4cout << "G4mplAtlasIonisationWithDeltaModel::SampleSecondaries: E(GeV)= "
245  // << kineticEnergy/GeV << " M(GeV)= " << mass/GeV
246  // << " tmin(MeV)= " << minKinEnergy/MeV << G4endl;
247 
248  const G4double totEnergy = kineticEnergy + mass;
249  const G4double etot2 = totEnergy*totEnergy;
250  const G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/etot2;
251 
252  // sampling without nuclear size effect
253  const G4double q = G4UniformRand();
254  const G4double deltaKinEnergy = minKinEnergy*maxKinEnergy
255  /(minKinEnergy*(1.0 - q) + maxKinEnergy*q);
256 
257  // delta-electron is produced
258  const G4double totMomentum = totEnergy*std::sqrt(beta2);
259  const G4double deltaMomentum =
260  std::sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*CLHEP::electron_mass_c2));
261  G4double cost = deltaKinEnergy * (totEnergy + CLHEP::electron_mass_c2) /
262  (deltaMomentum * totMomentum);
263  if(cost > 1.0) { cost = 1.0; }
264 
265  // to << deltaKinEnergy << " " << deltaMomentum << G4endl;
266 
267  const G4double sint = std::sqrt((1.0 - cost)*(1.0 + cost));
268 
269  const G4double phi = CLHEP::twopi * G4UniformRand() ;
270 
271  G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost);
272  const G4ThreeVector& direction = dp->GetMomentumDirection();
273  deltaDirection.rotateUz(direction);
274 
275  // create G4DynamicParticle object for delta ray
276  G4DynamicParticle* delta =
277  new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
278 
279  vdp->push_back(delta);
280 
281  // Change kinematics of primary particle
282  kineticEnergy -= deltaKinEnergy;
283  G4ThreeVector finalP = direction*totMomentum - deltaDirection*deltaMomentum;
284  finalP = finalP.unit();
285 
286  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
287  fParticleChange->SetProposedMomentumDirection(finalP);
288 }
289 
290 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
291 
292 G4double G4mplAtlasIonisationWithDeltaModel::SampleFluctuations(
293 #if G4VERSION_NUMBER < 1100
294  const G4MaterialCutsCouple* material,
295  const G4DynamicParticle* dp,
296  G4double tmax,
297  G4double length,
298  G4double meanLoss
299 #else
300  const G4MaterialCutsCouple* material,
301  const G4DynamicParticle* dp,
302  const G4double tcut,
303  const G4double tmax,
304  const G4double length,
305  const G4double meanLoss
306 #endif
307  )
308 {
309 #if G4VERSION_NUMBER < 1100
310  G4double siga = Dispersion(material->GetMaterial(),dp,tmax,length);
311 #else
312  G4double siga = Dispersion(material->GetMaterial(),dp,tcut,tmax,length);
313 #endif
314  G4double loss = meanLoss;
315  siga = std::sqrt(siga);
316  const G4double twomeanLoss = meanLoss + meanLoss;
317 
318  if(twomeanLoss < siga) {
319  G4double x;
320  do {
321  loss = twomeanLoss*G4UniformRand();
322  x = (loss - meanLoss)/siga;
323  } while (1.0 - 0.5*x*x < G4UniformRand());
324  } else {
325  do {
326  loss = G4RandGauss::shoot(meanLoss,siga);
327  } while (0.0 > loss || loss > twomeanLoss);
328  }
329  return loss;
330 }
331 
332 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
333 
334 G4double
335 G4mplAtlasIonisationWithDeltaModel::Dispersion(const G4Material* material,
336  const G4DynamicParticle* dp,
337 #if G4VERSION_NUMBER < 1100
338  G4double tmax,
339  G4double length
340 #else
341  const G4double /*tcut*/,
342  const G4double tmax,
343  const G4double length
344 #endif
345  )
346 {
347  G4double siga = 0.0;
348  const G4double tau = dp->GetKineticEnergy()/mass;
349  if(tau > 0.0) {
350  const G4double electronDensity = material->GetElectronDensity();
351  const G4double gam = tau + 1.0;
352  const G4double invbeta2 = (gam*gam)/(tau * (tau+2.0));
353  siga = (invbeta2 - 0.5) * CLHEP::twopi_mc2_rcl2 * tmax * length
354  * electronDensity * chargeSquare;
355  }
356  return siga;
357 }
358 
359 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
360 
361 G4double
362 G4mplAtlasIonisationWithDeltaModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
363  G4double kinEnergy)
364 {
365  const G4double tau = kinEnergy/mass;
366  return 2.0*CLHEP::electron_mass_c2*tau*(tau + 2.);
367 }
368 
369 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
cost
int cost(std::vector< std::string > &files, node &n, const std::string &directory="", bool deleteref=false, bool relocate=false)
Definition: hcg.cxx:921
TileDCSDataPlotter.dp
dp
Definition: TileDCSDataPlotter.py:840
xAOD::Electron
Electron_v1 Electron
Definition of the current "egamma version".
Definition: Event/xAOD/xAODEgamma/xAODEgamma/Electron.h:17
python.PerfMonSerializer.p
def p
Definition: PerfMonSerializer.py:743
python.PhysicalConstants.fine_structure_const
float fine_structure_const
Definition: PhysicalConstants.py:103
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:64
Monitored::Z
@ Z
Definition: HistogramFillerUtils.h:24
python.PhysicalConstants.twopi_mc2_rcl2
int twopi_mc2_rcl2
Definition: PhysicalConstants.py:110
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
x
#define x
dqt_zlumi_pandas.mass
mass
Definition: dqt_zlumi_pandas.py:170
pi
#define pi
Definition: TileMuonFitter.cxx:65
python.PhysicalConstants.electron_mass_c2
float electron_mass_c2
Definition: PhysicalConstants.py:86
python.PhysicalConstants.hbarc
float hbarc
Definition: PhysicalConstants.py:73
python.SystemOfUnits.cm2
int cm2
Definition: SystemOfUnits.py:88
python.CaloCondTools.g
g
Definition: CaloCondTools.py:15
min
#define min(a, b)
Definition: cfImp.cxx:40
twopi
constexpr double twopi
Definition: VertexPointEstimator.cxx:16
dqt_zlumi_alleff_HIST.B
B
Definition: dqt_zlumi_alleff_HIST.py:110
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
extractSporadic.q
list q
Definition: extractSporadic.py:98
python.SystemOfUnits.eplus
int eplus
Definition: SystemOfUnits.py:137
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
MuonParameters::beta
@ beta
Definition: MuonParamDefs.h:144
GeV
#define GeV
Definition: CaloTransverseBalanceVecMon.cxx:30
length
double length(const pvec &v)
Definition: FPGATrackSimLLPDoubletHoughTransformTool.cxx:26
fitman.k
k
Definition: fitman.py:528
Generated on Fri Jul 5 2024 21:16:12 for ATLAS Offline Software by   doxygen 1.8.18