ATLAS Offline Software
Loading...
Searching...
No Matches
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
73G4mplAtlasIonisationWithDeltaModel::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
99G4mplAtlasIonisationWithDeltaModel::~G4mplAtlasIonisationWithDeltaModel()
100{}
101
102//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
103
104void
105G4mplAtlasIonisationWithDeltaModel::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
122G4double
123G4mplAtlasIonisationWithDeltaModel::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
166G4double
167G4mplAtlasIonisationWithDeltaModel::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
204G4double
205G4mplAtlasIonisationWithDeltaModel::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
228G4double
229G4mplAtlasIonisationWithDeltaModel::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
243void
244G4mplAtlasIonisationWithDeltaModel::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
304G4double 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
346G4double
347G4mplAtlasIonisationWithDeltaModel::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
373G4double
374G4mplAtlasIonisationWithDeltaModel::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....
phi
Scalar phi() const
phi method
Definition AmgMatrixBasePlugin.h:67
length
double length(const pvec &v)
Definition FPGATrackSimLLPDoubletHoughTransformTool.cxx:26
x
#define x
cost
int cost(std::vector< std::string > &files, node &n, const std::string &directory="", bool deleteref=false, bool relocate=false)
Definition hcg.cxx:922
Base_Fragment.mass
mass
Definition Sherpa_i/share/common/Base_Fragment.py:59
dqt_zlumi_alleff_HIST.B
B
Definition dqt_zlumi_alleff_HIST.py:110
std
STL namespace.
Generated on for ATLAS Offline Software by doxygen 1.14.0