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  BetheBlochModel = new G4BetheBlochModel(); // initialize the model for use in ComputeDEDXPerVolume and ComputeCrossSectionPerElectron
93  G4cout << "### Monopole ionisation model with d-electron production, Gmag= "
94  << magCharge/CLHEP::eplus << G4endl;
95 }
96 
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
98 
99 G4mplAtlasIonisationWithDeltaModel::~G4mplAtlasIonisationWithDeltaModel()
100 {}
101 
102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
103 
104 void
105 G4mplAtlasIonisationWithDeltaModel::Initialise(const G4ParticleDefinition* p,
106  const G4DataVector& G4DataVec)
107 {
108  monopole = p;
109  mass = monopole->GetPDGMass();
110  // if(!fParticleChange) { fParticleChange = GetParticleChangeForLoss(); }
111 
112  if(pParticleChange)
113  fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
114  else
115  fParticleChange = new G4ParticleChangeForLoss();
116 
117  BetheBlochModel->Initialise(p, G4DataVec);
118 }
119 
120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
121 
122 G4double
123 G4mplAtlasIonisationWithDeltaModel::ComputeDEDXPerVolume(const G4Material* material,
124  const G4ParticleDefinition* p,
125  G4double kineticEnergy,
126  G4double maxEnergy)
127 {
128  const G4double tmax = MaxSecondaryEnergy(p,kineticEnergy);
129  const G4double cutEnergy = std::min(tmax, maxEnergy);
130  const G4double tau = kineticEnergy / mass;
131  const G4double gam = tau + 1.0;
132  const G4double bg2 = tau * (tau + 2.0);
133  const G4double beta2 = bg2 / (gam * gam);
134  const G4double beta = std::sqrt(beta2);
135 
136  // low-energy asymptotic formula
137  G4double dedx = dedxlim*beta*material->GetDensity();
138 
139  // above asymptotic
140  if(beta > betalow) {
141 
142  // high energy
143  if(beta >= betalim) {
144  dedx = ComputeDEDXAhlen(material, bg2, cutEnergy);
145 
146  } else {
147 
148  const G4double dedx1 = dedxlim*betalow*material->GetDensity();
149  const G4double dedx2 = ComputeDEDXAhlen(material, bg2lim, cutEnergy);
150 
151  // extrapolation between two formula
152  const G4double kapa2 = beta - betalow;
153  const G4double kapa1 = betalim - beta;
154  dedx = (kapa1*dedx1 + kapa2*dedx2)/(kapa1 + kapa2);
155  }
156  }
157  if(p->GetPDGCharge() != 0) {
158  G4double dedxBetheBloch = BetheBlochModel->ComputeDEDXPerVolume(material, p, kineticEnergy, maxEnergy);
159  dedx += dedxBetheBloch; // only dyons would be electrically charged here, then add Bethe-Bloch DEDX to total
160  }
161  return dedx;
162 }
163 
164 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
165 
166 G4double
167 G4mplAtlasIonisationWithDeltaModel::ComputeDEDXAhlen(const G4Material* material,
168  G4double bg2,
169  G4double cutEnergy)
170 {
171  const G4double eDensity = material->GetElectronDensity();
172  const G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
173 
174  // Ahlen's formula for nonconductors, [1]p157, f(5.7)
175  G4double dedx =
176  0.5*(std::log(2.0 * CLHEP::electron_mass_c2 * bg2*cutEnergy / (eexc*eexc)) - 1.0);
177 
178  // Kazama et al. cross-section correction
179  G4double k = 0.406;
180  if(nmpl > 1) { k = 0.346; }
181 
182  // Bloch correction
183  const G4double B[7] = { 0.0, 0.248, 0.672, 1.022, 1.243, 1.464, 1.685};
184 
185  if(nmpl > 6.0) { dedx += 0.5 * k - B[6]; }
186  else if(nmpl < 1.0) { dedx += 0.5 * k - B[1]; }
187  else {
188  dedx += 0.5 * k - B[G4int(std::abs(magCharge) * 2 * CLHEP::fine_structure_const + 0.5)];
189  }
190  // density effect correction
191  const G4double x = std::log(bg2)/twoln10;
192 
193  dedx -= material->GetIonisation()->DensityCorrection(x);
194 
195  // now compute the total ionization loss
196  dedx *= pi_hbarc2_over_mc2 * eDensity * nmpl * nmpl;
197 
198  if (dedx < 0.0) { dedx = 0; }
199  return dedx;
200 }
201 
202 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
203 
204 G4double
205 G4mplAtlasIonisationWithDeltaModel::ComputeCrossSectionPerElectron(
206  const G4ParticleDefinition* p,
207  G4double kineticEnergy,
208  G4double cutEnergy,
209  G4double maxKinEnergy)
210 {
211  G4double cross = 0.0;
212  const G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
213  const G4double maxEnergy = std::min(tmax,maxKinEnergy);
214  if(cutEnergy < maxEnergy) {
215  cross = (1.0/cutEnergy - 1.0/maxEnergy)*CLHEP::twopi_mc2_rcl2*chargeSquare;
216  }
217 
218  if(p->GetPDGCharge() != 0) { // if dyon, include the Bethe-Bloch Cross-section to total calculation
219  G4double crossBetheBloch = BetheBlochModel->ComputeCrossSectionPerElectron(p, kineticEnergy, cutEnergy, maxKinEnergy);
220  cross += crossBetheBloch; // additive since they are different processes
221  }
222 
223  return cross;
224 }
225 
226 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
227 
228 G4double
229 G4mplAtlasIonisationWithDeltaModel::ComputeCrossSectionPerAtom(
230  const G4ParticleDefinition* p,
231  G4double kineticEnergy,
232  G4double Z, G4double,
233  G4double cutEnergy,
234  G4double maxEnergy)
235 {
236  const G4double cross =
237  Z*ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
238  return cross;
239 }
240 
241 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
242 
243 void
244 G4mplAtlasIonisationWithDeltaModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
245  const G4MaterialCutsCouple*,
246  const G4DynamicParticle* dp,
247  G4double minKinEnergy,
248  G4double maxEnergy)
249 {
250  G4double kineticEnergy = dp->GetKineticEnergy();
251  const G4double tmax = MaxSecondaryEnergy(dp->GetDefinition(),kineticEnergy);
252 
253  G4double maxKinEnergy = std::min(maxEnergy,tmax);
254  if(minKinEnergy >= maxKinEnergy) { return; }
255 
256  //G4cout << "G4mplAtlasIonisationWithDeltaModel::SampleSecondaries: E(GeV)= "
257  // << kineticEnergy/GeV << " M(GeV)= " << mass/GeV
258  // << " tmin(MeV)= " << minKinEnergy/MeV << G4endl;
259 
260  const G4double totEnergy = kineticEnergy + mass;
261  const G4double etot2 = totEnergy*totEnergy;
262  const G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/etot2;
263 
264  // sampling without nuclear size effect
265  const G4double q = G4UniformRand();
266  const G4double deltaKinEnergy = minKinEnergy*maxKinEnergy
267  /(minKinEnergy*(1.0 - q) + maxKinEnergy*q);
268 
269  // delta-electron is produced
270  const G4double totMomentum = totEnergy*std::sqrt(beta2);
271  const G4double deltaMomentum =
272  std::sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*CLHEP::electron_mass_c2));
273  G4double cost = deltaKinEnergy * (totEnergy + CLHEP::electron_mass_c2) /
274  (deltaMomentum * totMomentum);
275  if(cost > 1.0) { cost = 1.0; }
276 
277  // to << deltaKinEnergy << " " << deltaMomentum << G4endl;
278 
279  const G4double sint = std::sqrt((1.0 - cost)*(1.0 + cost));
280 
281  const G4double phi = CLHEP::twopi * G4UniformRand() ;
282 
283  G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost);
284  const G4ThreeVector& direction = dp->GetMomentumDirection();
285  deltaDirection.rotateUz(direction);
286 
287  // create G4DynamicParticle object for delta ray
288  G4DynamicParticle* delta =
289  new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
290 
291  vdp->push_back(delta);
292 
293  // Change kinematics of primary particle
294  kineticEnergy -= deltaKinEnergy;
295  G4ThreeVector finalP = direction*totMomentum - deltaDirection*deltaMomentum;
296  finalP = finalP.unit();
297 
298  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
299  fParticleChange->SetProposedMomentumDirection(finalP);
300 }
301 
302 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
303 
304 G4double G4mplAtlasIonisationWithDeltaModel::SampleFluctuations(
305 #if G4VERSION_NUMBER < 1100
306  const G4MaterialCutsCouple* material,
307  const G4DynamicParticle* dp,
308  G4double tmax,
309  G4double length,
310  G4double meanLoss
311 #else
312  const G4MaterialCutsCouple* material,
313  const G4DynamicParticle* dp,
314  const G4double tcut,
315  const G4double tmax,
316  const G4double length,
317  const G4double meanLoss
318 #endif
319  )
320 {
321 #if G4VERSION_NUMBER < 1100
322  G4double siga = Dispersion(material->GetMaterial(),dp,tmax,length);
323 #else
324  G4double siga = Dispersion(material->GetMaterial(),dp,tcut,tmax,length);
325 #endif
326  G4double loss = meanLoss;
327  siga = std::sqrt(siga);
328  const G4double twomeanLoss = meanLoss + meanLoss;
329 
330  if(twomeanLoss < siga) {
331  G4double x;
332  do {
333  loss = twomeanLoss*G4UniformRand();
334  x = (loss - meanLoss)/siga;
335  } while (1.0 - 0.5*x*x < G4UniformRand());
336  } else {
337  do {
338  loss = G4RandGauss::shoot(meanLoss,siga);
339  } while (0.0 > loss || loss > twomeanLoss);
340  }
341  return loss;
342 }
343 
344 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
345 
346 G4double
347 G4mplAtlasIonisationWithDeltaModel::Dispersion(const G4Material* material,
348  const G4DynamicParticle* dp,
349 #if G4VERSION_NUMBER < 1100
350  G4double tmax,
351  G4double length
352 #else
353  const G4double /*tcut*/,
354  const G4double tmax,
355  const G4double length
356 #endif
357  )
358 {
359  G4double siga = 0.0;
360  const G4double tau = dp->GetKineticEnergy()/mass;
361  if(tau > 0.0) {
362  const G4double electronDensity = material->GetElectronDensity();
363  const G4double gam = tau + 1.0;
364  const G4double invbeta2 = (gam*gam)/(tau * (tau+2.0));
365  siga = (invbeta2 - 0.5) * CLHEP::twopi_mc2_rcl2 * tmax * length
366  * electronDensity * chargeSquare;
367  }
368  return siga;
369 }
370 
371 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
372 
373 G4double
374 G4mplAtlasIonisationWithDeltaModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
375  G4double kinEnergy)
376 {
377  const G4double tau = kinEnergy/mass;
378  return 2.0*CLHEP::electron_mass_c2*tau*(tau + 2.);
379 }
380 
381 //....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:842
GeV
#define GeV
Definition: PhysicsAnalysis/TauID/TauAnalysisTools/Root/HelperFunctions.cxx:18
xAOD::Electron
Electron_v1 Electron
Definition of the current "egamma version".
Definition: Event/xAOD/xAODEgamma/xAODEgamma/Electron.h:17
python.PhysicalConstants.fine_structure_const
float fine_structure_const
Definition: PhysicalConstants.py:113
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:67
Base_Fragment.mass
mass
Definition: Sherpa_i/share/common/Base_Fragment.py:59
Monitored::Z
@ Z
Definition: HistogramFillerUtils.h:24
min
constexpr double min()
Definition: ap_fixedTest.cxx:26
python.PhysicalConstants.twopi_mc2_rcl2
tuple twopi_mc2_rcl2
Definition: PhysicalConstants.py:120
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
x
#define x
pi
#define pi
Definition: TileMuonFitter.cxx:65
python.PhysicalConstants.electron_mass_c2
float electron_mass_c2
Definition: PhysicalConstants.py:96
python.SystemOfUnits.eplus
float eplus
Definition: SystemOfUnits.py:155
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:209
python.PhysicalConstants.hbarc
float hbarc
Definition: PhysicalConstants.py:83
python.CaloCondTools.g
g
Definition: CaloCondTools.py:15
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:97
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
MuonParameters::beta
@ beta
Definition: MuonParamDefs.h:144
python.SystemOfUnits.cm2
float cm2
Definition: SystemOfUnits.py:103
length
double length(const pvec &v)
Definition: FPGATrackSimLLPDoubletHoughTransformTool.cxx:26
fitman.k
k
Definition: fitman.py:528
Generated on Tue Sep 30 2025 21:10:12 for ATLAS Offline Software by   doxygen 1.8.18