CMS 3D CMS Logo
CMSmplIonisationWithDeltaModel.cc
Go to the documentation of this file.
1 //
2 // -------------------------------------------------------------------
3 //
4 //
5 // File name: CMSmplIonisationWithDeltaModel
6 //
7 // Author: Vladimir Ivanchenko
8 //
9 // Creation date: 02.03.2019 copied from Geant4 10.5p01
10 //
11 //
12 // -------------------------------------------------------------------
13 // References
14 // [1] Steven P. Ahlen: Energy loss of relativistic heavy ionizing particles,
15 // S.P. Ahlen, Rev. Mod. Phys 52(1980), p121
16 // [2] K.A. Milton arXiv:hep-ex/0602040
17 // [3] S.P. Ahlen and K. Kinoshita, Phys. Rev. D26 (1982) 2347
18 
19 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
20 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
21 
22 #include "SimG4Core/PhysicsLists/interface/CMSmplIonisationWithDeltaModel.h"
23 #include "Randomize.hh"
24 #include "G4PhysicalConstants.hh"
25 #include "G4SystemOfUnits.hh"
26 #include "G4ParticleChangeForLoss.hh"
27 #include "G4Electron.hh"
28 #include "G4DynamicParticle.hh"
29 #include "G4ProductionCutsTable.hh"
30 #include "G4MaterialCutsCouple.hh"
31 #include "G4Log.hh"
32 #include "G4Pow.hh"
33 
34 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
35 
36 using namespace std;
37 
38 std::vector<G4double>* CMSmplIonisationWithDeltaModel::dedx0 = nullptr;
39 
40 CMSmplIonisationWithDeltaModel::CMSmplIonisationWithDeltaModel(G4double mCharge, const G4String& nam)
41  : G4VEmModel(nam),
42  G4VEmFluctuationModel(nam),
43  magCharge(mCharge),
44  twoln10(std::log(100.0)),
45  betalow(0.01),
46  betalim(0.1),
47  beta2lim(betalim * betalim),
48  bg2lim(beta2lim * (1.0 + beta2lim)) {
49  nmpl = G4lrint(std::abs(magCharge) * 2 * fine_structure_const);
50  if (nmpl > 6) {
51  nmpl = 6;
52  } else if (nmpl < 1) {
53  nmpl = 1;
54  }
55  pi_hbarc2_over_mc2 = pi * hbarc * hbarc / electron_mass_c2;
56  chargeSquare = magCharge * magCharge;
57  dedxlim = 45. * nmpl * nmpl * GeV * cm2 / g;
58  fParticleChange = nullptr;
59  theElectron = G4Electron::Electron();
60  G4cout << "### Monopole ionisation model with d-electron production, Gmag= " << magCharge / eplus << G4endl;
61  monopole = nullptr;
62  mass = 0.0;
63 }
64 
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66 
67 CMSmplIonisationWithDeltaModel::~CMSmplIonisationWithDeltaModel() {
68  if (IsMaster()) {
69  delete dedx0;
70  }
71 }
72 
73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
74 
75 void CMSmplIonisationWithDeltaModel::SetParticle(const G4ParticleDefinition* p) {
76  monopole = p;
77  mass = monopole->GetPDGMass();
78  G4double emin = std::min(LowEnergyLimit(), 0.1 * mass * (1. / sqrt(1. - betalow * betalow) - 1.));
79  G4double emax = std::max(HighEnergyLimit(), 10 * mass * (1. / sqrt(1. - beta2lim) - 1.));
80  SetLowEnergyLimit(emin);
81  SetHighEnergyLimit(emax);
82 }
83 
84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
85 
86 void CMSmplIonisationWithDeltaModel::Initialise(const G4ParticleDefinition* p, const G4DataVector&) {
87  if (!monopole) {
88  SetParticle(p);
89  }
90  if (!fParticleChange) {
91  fParticleChange = GetParticleChangeForLoss();
92  }
93  if (IsMaster()) {
94  if (!dedx0) {
95  dedx0 = new std::vector<G4double>;
96  }
97  G4ProductionCutsTable* theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
98  G4int numOfCouples = theCoupleTable->GetTableSize();
99  G4int n = dedx0->size();
100  if (n < numOfCouples) {
101  dedx0->resize(numOfCouples);
102  }
103  G4Pow* g4calc = G4Pow::GetInstance();
104 
105  // initialise vector
106  for (G4int i = 0; i < numOfCouples; ++i) {
107  const G4Material* material = theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
108  G4double eDensity = material->GetElectronDensity();
109  G4double vF2 = 2 * electron_Compton_length * g4calc->A13(3. * pi * pi * eDensity);
110  (*dedx0)[i] = pi_hbarc2_over_mc2 * eDensity * nmpl * nmpl * (G4Log(vF2 / fine_structure_const) - 0.5) / vF2;
111  }
112  }
113 }
114 
115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
116 
117 G4double CMSmplIonisationWithDeltaModel::MinEnergyCut(const G4ParticleDefinition*, const G4MaterialCutsCouple* couple) {
118  return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy();
119 }
120 
121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
122 
123 G4double CMSmplIonisationWithDeltaModel::ComputeDEDXPerVolume(const G4Material* material,
124  const G4ParticleDefinition* p,
125  G4double kineticEnergy,
126  G4double maxEnergy) {
127  if (!monopole) {
128  SetParticle(p);
129  }
130  G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
131  G4double cutEnergy = std::min(tmax, maxEnergy);
132  cutEnergy = std::max(LowEnergyLimit(), cutEnergy);
133  G4double tau = kineticEnergy / mass;
134  G4double gam = tau + 1.0;
135  G4double bg2 = tau * (tau + 2.0);
136  G4double beta2 = bg2 / (gam * gam);
137  G4double beta = sqrt(beta2);
138 
139  // low-energy asymptotic formula
140  G4double dedx = (*dedx0)[CurrentCouple()->GetIndex()] * beta;
141 
142  // above asymptotic
143  if (beta > betalow) {
144  // high energy
145  if (beta >= betalim) {
146  dedx = ComputeDEDXAhlen(material, bg2, cutEnergy);
147 
148  } else {
149  G4double dedx1 = (*dedx0)[CurrentCouple()->GetIndex()] * betalow;
150  G4double dedx2 = ComputeDEDXAhlen(material, bg2lim, cutEnergy);
151 
152  // extrapolation between two formula
153  G4double kapa2 = beta - betalow;
154  G4double kapa1 = betalim - beta;
155  dedx = (kapa1 * dedx1 + kapa2 * dedx2) / (kapa1 + kapa2);
156  }
157  }
158  return dedx;
159 }
160 
161 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
162 
163 G4double CMSmplIonisationWithDeltaModel::ComputeDEDXAhlen(const G4Material* material,
164  G4double bg2,
165  G4double cutEnergy) {
166  G4double eDensity = material->GetElectronDensity();
167  G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
168 
169  // Ahlen's formula for nonconductors, [1]p157, f(5.7)
170  G4double dedx = 0.5 * (G4Log(2.0 * 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) {
175  k = 0.346;
176  }
177 
178  // Bloch correction
179  const G4double B[7] = {0.0, 0.248, 0.672, 1.022, 1.243, 1.464, 1.685};
180 
181  dedx += 0.5 * k - B[nmpl];
182 
183  // density effect correction
184  G4double x = G4Log(bg2) / twoln10;
185  dedx -= material->GetIonisation()->DensityCorrection(x);
186 
187  // now compute the total ionization loss
188  dedx *= pi_hbarc2_over_mc2 * eDensity * nmpl * nmpl;
189 
190  dedx = std::max(dedx, 0.0);
191  return dedx;
192 }
193 
194 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
195 
196 G4double CMSmplIonisationWithDeltaModel::ComputeCrossSectionPerElectron(const G4ParticleDefinition* p,
197  G4double kineticEnergy,
198  G4double cut,
199  G4double maxKinEnergy) {
200  if (!monopole) {
201  SetParticle(p);
202  }
203  G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
204  G4double maxEnergy = std::min(tmax, maxKinEnergy);
205  G4double cutEnergy = std::max(LowEnergyLimit(), cut);
206  G4double cross =
207  (cutEnergy < maxEnergy) ? (0.5 / cutEnergy - 0.5 / maxEnergy) * pi_hbarc2_over_mc2 * nmpl * nmpl : 0.0;
208  return cross;
209 }
210 
211 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
212 
213 G4double CMSmplIonisationWithDeltaModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition* p,
214  G4double kineticEnergy,
215  G4double Z,
216  G4double,
217  G4double cutEnergy,
218  G4double maxEnergy) {
219  G4double cross = Z * ComputeCrossSectionPerElectron(p, kineticEnergy, cutEnergy, maxEnergy);
220  return cross;
221 }
222 
223 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
224 
225 void CMSmplIonisationWithDeltaModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
226  const G4MaterialCutsCouple*,
227  const G4DynamicParticle* dp,
228  G4double minKinEnergy,
229  G4double maxEnergy) {
230  G4double kineticEnergy = dp->GetKineticEnergy();
231  G4double tmax = MaxSecondaryEnergy(dp->GetDefinition(), kineticEnergy);
232 
233  G4double maxKinEnergy = std::min(maxEnergy, tmax);
234  if (minKinEnergy >= maxKinEnergy) {
235  return;
236  }
237 
238  //G4cout << "CMSmplIonisationWithDeltaModel::SampleSecondaries: E(GeV)= "
239  // << kineticEnergy/GeV << " M(GeV)= " << mass/GeV
240  // << " tmin(MeV)= " << minKinEnergy/MeV << G4endl;
241 
242  G4double totEnergy = kineticEnergy + mass;
243  G4double etot2 = totEnergy * totEnergy;
244  G4double beta2 = kineticEnergy * (kineticEnergy + 2.0 * mass) / etot2;
245 
246  // sampling without nuclear size effect
247  G4double q = G4UniformRand();
248  G4double deltaKinEnergy = minKinEnergy * maxKinEnergy / (minKinEnergy * (1.0 - q) + maxKinEnergy * q);
249 
250  // delta-electron is produced
251  G4double totMomentum = totEnergy * sqrt(beta2);
252  G4double deltaMomentum = sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0 * electron_mass_c2));
253  G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) / (deltaMomentum * totMomentum);
254  cost = std::min(cost, 1.0);
255 
256  G4double sint = sqrt((1.0 - cost) * (1.0 + cost));
257 
258  G4double phi = twopi * G4UniformRand();
259 
260  G4ThreeVector deltaDirection(sint * cos(phi), sint * sin(phi), cost);
261  const G4ThreeVector& direction = dp->GetMomentumDirection();
262  deltaDirection.rotateUz(direction);
263 
264  // create G4DynamicParticle object for delta ray
265  G4DynamicParticle* delta = new G4DynamicParticle(theElectron, deltaDirection, deltaKinEnergy);
266 
267  vdp->push_back(delta);
268 
269  // Change kinematics of primary particle
270  kineticEnergy -= deltaKinEnergy;
271  G4ThreeVector finalP = direction * totMomentum - deltaDirection * deltaMomentum;
272  finalP = finalP.unit();
273 
274  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
275  fParticleChange->SetProposedMomentumDirection(finalP);
276 }
277 
278 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
279 
280 G4double CMSmplIonisationWithDeltaModel::SampleFluctuations(const G4MaterialCutsCouple* couple,
281  const G4DynamicParticle* dp,
282  G4double tmax,
283  G4double length,
284  G4double meanLoss) {
285  G4double siga = Dispersion(couple->GetMaterial(), dp, tmax, length);
286  G4double loss = meanLoss;
287  siga = sqrt(siga);
288  G4double twomeanLoss = meanLoss + meanLoss;
289 
290  if (twomeanLoss < siga) {
291  G4double x;
292  do {
293  loss = twomeanLoss * G4UniformRand();
294  x = (loss - meanLoss) / siga;
295  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
296  } while (1.0 - 0.5 * x * x < G4UniformRand());
297  } else {
298  do {
299  loss = G4RandGauss::shoot(meanLoss, siga);
300  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
301  } while (0.0 > loss || loss > twomeanLoss);
302  }
303  return loss;
304 }
305 
306 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
307 
308 G4double CMSmplIonisationWithDeltaModel::Dispersion(const G4Material* material,
309  const G4DynamicParticle* dp,
310  G4double tmax,
311  G4double length) {
312  G4double siga = 0.0;
313  G4double tau = dp->GetKineticEnergy() / mass;
314  if (tau > 0.0) {
315  G4double electronDensity = material->GetElectronDensity();
316  G4double gam = tau + 1.0;
317  G4double invbeta2 = (gam * gam) / (tau * (tau + 2.0));
318  siga = (invbeta2 - 0.5) * twopi_mc2_rcl2 * tmax * length * electronDensity * chargeSquare;
319  }
320  return siga;
321 }
322 
323 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
324 
325 G4double CMSmplIonisationWithDeltaModel::MaxSecondaryEnergy(const G4ParticleDefinition*, G4double kinEnergy) {
326  G4double tau = kinEnergy / mass;
327  return 2.0 * electron_mass_c2 * tau * (tau + 2.);
328 }
329 
330 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
CMSmplIonisationWithDeltaModel::magCharge
G4double magCharge
Definition: CMSmplIonisationWithDeltaModel.h:83
mps_fire.i
i
Definition: mps_fire.py:428
dqmiodumpmetadata.n
n
Definition: dqmiodumpmetadata.py:28
CMSmplIonisationWithDeltaModel::fParticleChange
G4ParticleChangeForLoss * fParticleChange
Definition: CMSmplIonisationWithDeltaModel.h:80
TkAlMuonSelectors_cfi.cut
cut
Definition: TkAlMuonSelectors_cfi.py:5
metsig::tau
Definition: SignAlgoResolutions.h:49
CMSmplIonisationWithDeltaModel::betalow
G4double betalow
Definition: CMSmplIonisationWithDeltaModel.h:85
CMSmplIonisationWithDeltaModel::ComputeDEDXPerVolume
G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy) override
Definition: CMSmplIonisationWithDeltaModel.cc:123
min
T min(T a, T b)
Definition: MathUtil.h:58
HLT_FULL_cff.beta
beta
Definition: HLT_FULL_cff.py:8686
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
CMSmplIonisationWithDeltaModel::theElectron
G4ParticleDefinition * theElectron
Definition: CMSmplIonisationWithDeltaModel.h:79
CMSmplIonisationWithDeltaModel::pi_hbarc2_over_mc2
G4double pi_hbarc2_over_mc2
Definition: CMSmplIonisationWithDeltaModel.h:92
CMSmplIonisationWithDeltaModel::chargeSquare
G4double chargeSquare
Definition: CMSmplIonisationWithDeltaModel.h:89
DDAxes::x
CMSmplIonisationWithDeltaModel::dedxlim
G4double dedxlim
Definition: CMSmplIonisationWithDeltaModel.h:90
CMSmplIonisationWithDeltaModel::ComputeDEDXAhlen
G4double ComputeDEDXAhlen(const G4Material *material, G4double bg2, G4double cut)
Definition: CMSmplIonisationWithDeltaModel.cc:163
CMSmplIonisationWithDeltaModel::MinEnergyCut
G4double MinEnergyCut(const G4ParticleDefinition *, const G4MaterialCutsCouple *couple) override
Definition: CMSmplIonisationWithDeltaModel.cc:117
tmax
static const double tmax[3]
Definition: CastorTimeSlew.cc:7
CMSmplIonisationWithDeltaModel::ComputeCrossSectionPerAtom
G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kineticEnergy, G4double Z, G4double A, G4double cutEnergy, G4double maxEnergy) override
Definition: CMSmplIonisationWithDeltaModel.cc:213
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
CMSmplIonisationWithDeltaModel::Initialise
void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
Definition: CMSmplIonisationWithDeltaModel.cc:86
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
CMSmplIonisationWithDeltaModel::bg2lim
G4double bg2lim
Definition: CMSmplIonisationWithDeltaModel.h:88
CMSmplIonisationWithDeltaModel::mass
G4double mass
Definition: CMSmplIonisationWithDeltaModel.h:82
Calorimetry_cff.dp
dp
Definition: Calorimetry_cff.py:158
CMSmplIonisationWithDeltaModel::Dispersion
G4double Dispersion(const G4Material *, const G4DynamicParticle *, G4double tmax, G4double length) override
Definition: CMSmplIonisationWithDeltaModel.cc:308
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
CMSmplIonisationWithDeltaModel::CMSmplIonisationWithDeltaModel
CMSmplIonisationWithDeltaModel(G4double mCharge, const G4String &nam="mplIonisationWithDelta")
Definition: CMSmplIonisationWithDeltaModel.cc:40
CMSmplIonisationWithDeltaModel::SampleFluctuations
G4double SampleFluctuations(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmax, G4double length, G4double meanLoss) override
Definition: CMSmplIonisationWithDeltaModel.cc:280
dqmdumpme.k
k
Definition: dqmdumpme.py:60
CMSmplIonisationWithDeltaModel::nmpl
G4int nmpl
Definition: CMSmplIonisationWithDeltaModel.h:91
CMSmplIonisationWithDeltaModel::ComputeCrossSectionPerElectron
virtual G4double ComputeCrossSectionPerElectron(const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
Definition: CMSmplIonisationWithDeltaModel.cc:196
CMSmplIonisationWithDeltaModel::~CMSmplIonisationWithDeltaModel
~CMSmplIonisationWithDeltaModel() override
Definition: CMSmplIonisationWithDeltaModel.cc:67
cross
Basic3DVector cross(const Basic3DVector &v) const
Vector product, or "cross" product, with a vector of same type.
Definition: Basic3DVectorLD.h:225
SiStripPI::max
Definition: SiStripPayloadInspectorHelper.h:169
GeV
const double GeV
Definition: MathUtil.h:16
dumpMFGeometry_cfg.delta
delta
Definition: dumpMFGeometry_cfg.py:25
CMSmplIonisationWithDeltaModel.h
hbarc
const double hbarc
Definition: MathUtil.h:18
submitPVResolutionJobs.q
q
Definition: submitPVResolutionJobs.py:84
particleFlowClusterECALTimeSelected_cfi.maxEnergy
maxEnergy
Definition: particleFlowClusterECALTimeSelected_cfi.py:10
nanoDQM_cff.Electron
Electron
Definition: nanoDQM_cff.py:62
CMSmplIonisationWithDeltaModel::twoln10
G4double twoln10
Definition: CMSmplIonisationWithDeltaModel.h:84
DDAxes::phi
TtFullHadDaughter::B
static const std::string B
Definition: TtFullHadronicEvent.h:9
std
Definition: JetResolutionObject.h:76
CMSmplIonisationWithDeltaModel::beta2lim
G4double beta2lim
Definition: CMSmplIonisationWithDeltaModel.h:87
CMSmplIonisationWithDeltaModel::dedx0
static std::vector< G4double > * dedx0
Definition: CMSmplIonisationWithDeltaModel.h:94
CMSmplIonisationWithDeltaModel::betalim
G4double betalim
Definition: CMSmplIonisationWithDeltaModel.h:86
dqm-mbProfile.log
log
Definition: dqm-mbProfile.py:17
CMSmplIonisationWithDeltaModel::SetParticle
void SetParticle(const G4ParticleDefinition *p)
Definition: CMSmplIonisationWithDeltaModel.cc:75
CMSmplIonisationWithDeltaModel::monopole
const G4ParticleDefinition * monopole
Definition: CMSmplIonisationWithDeltaModel.h:78
pi
const Double_t pi
Definition: trackSplitPlot.h:36
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
BeamSpotPI::Z
Definition: BeamSpotPayloadInspectorHelper.h:33
CMSmplIonisationWithDeltaModel::MaxSecondaryEnergy
G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kinEnergy) override
Definition: CMSmplIonisationWithDeltaModel.cc:325
ecalTB2006H4_GenSimDigiReco_cfg.G4cout
G4cout
Definition: ecalTB2006H4_GenSimDigiReco_cfg.py:285
g
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g
Definition: Activities.doc:4
CMSmplIonisationWithDeltaModel::SampleSecondaries
void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
Definition: CMSmplIonisationWithDeltaModel.cc:225
Generated for CMSSW Reference Manual by doxygen 1.8.16