Belle II Software  release-05-01-25
G4mplIonisationWithDeltaModel.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2018 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Dmitrii Neverov *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 // References
12 // [1] Steven P. Ahlen: Energy loss of relativistic heavy ionizing particles,
13 // S.P. Ahlen, Rev. Mod. Phys 52(1980), p121
14 // [2] K.A. Milton arXiv:hep-ex/0602040
15 // [3] S.P. Ahlen and K. Kinoshita, Phys. Rev. D26 (1982) 2347
16 // [4] Y. Kazama et al., Phys. Rev. D15 (1977) 2287-2299
17 // [5] S.P. Ahlen, Phys. Rev. D17 (1978) 229-233
18 
19 // modified from GEANT4 exoticphysics/monopole/*
20 
21 #include <simulation/monopoles/G4mplIonisationWithDeltaModel.h>
22 
23 #include <Randomize.hh>
24 #include <CLHEP/Units/PhysicalConstants.h>
25 #include <CLHEP/Units/SystemOfUnits.h>
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 #include <framework/logging/Logger.h>
35 
36 using namespace std;
37 using namespace Belle2;
38 using namespace Belle2::Monopoles;
39 using namespace CLHEP;
40 
41 std::vector<G4double>* G4mplIonisationWithDeltaModel::dedx0 = nullptr;
42 
43 G4mplIonisationWithDeltaModel::G4mplIonisationWithDeltaModel(G4double mCharge,
44  const G4String& nam)
45  : G4VEmModel(nam), G4VEmFluctuationModel(nam),
46  magCharge(mCharge),
47  twoln10(log(100.0)),
48  betalow(0.01),
49  betalim(0.1),
50  beta2lim(betalim * betalim),
51  bg2lim(beta2lim * (1.0 + beta2lim)),
52  pi_hbarc2_over_mc2(pi * hbarc * hbarc / electron_mass_c2)
53 {
54  nmpl = magCharge * 2 * fine_structure_const;
55  chargeSquare = magCharge * magCharge * 4 * fine_structure_const *
56  fine_structure_const;
57  //NOTE Formulas below assume Dirac charge units for magnetic charge, g_D = 68.5e
58  dedxlim = 45. * chargeSquare * GeV * cm2 / g;
59  fParticleChange = nullptr;
60  theElectron = G4Electron::Electron();
61  B2INFO("Monopole ionisation model with d-electron production, Gmag= " << magCharge / eplus);
62  if (nmpl >= 6)
63  B2WARNING("Monopole charge Gmag= " << magCharge / eplus << "e not reasonable. Please choose a value smaller than 411e.");
64  monopole = nullptr;
65  mass = 0.0;
66 }
67 
69 {
70  if (IsMaster()) { delete dedx0; }
71 }
72 
73 void G4mplIonisationWithDeltaModel::SetParticle(const G4ParticleDefinition* p)
74 {
75  monopole = p;
76  mass = monopole->GetPDGMass();
77  G4double emin =
78  std::min(LowEnergyLimit(), 0.1 * mass * (1. / sqrt(1. - betalow * betalow) - 1.));
79  G4double emax =
80  std::max(HighEnergyLimit(), 10 * mass * (1. / sqrt(1. - beta2lim) - 1.));
81  SetLowEnergyLimit(emin);
82  SetHighEnergyLimit(emax);
83 }
84 
85 void
86 G4mplIonisationWithDeltaModel::Initialise(const G4ParticleDefinition* p,
87  const G4DataVector&)
88 {
89  if (!monopole) { SetParticle(p); }
90  if (!fParticleChange) { fParticleChange = GetParticleChangeForLoss(); }
91  if (IsMaster()) {
92  if (!dedx0) { dedx0 = new std::vector<G4double>; }
93  G4ProductionCutsTable* theCoupleTable =
94  G4ProductionCutsTable::GetProductionCutsTable();
95  G4int numOfCouples = theCoupleTable->GetTableSize();
96  G4int n = dedx0->size();
97  if (n < numOfCouples) { dedx0->resize(numOfCouples); }
98  G4Pow* g4calc = G4Pow::GetInstance();
99 
100  // initialise vector
101  for (G4int i = 0; i < numOfCouples; ++i) {
102 
103  const G4Material* material =
104  theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
105  G4double eDensity = material->GetElectronDensity();
106  G4double vF = electron_Compton_length * g4calc->A13(3.*pi * pi * eDensity);
107  (*dedx0)[i] = pi_hbarc2_over_mc2 * eDensity * chargeSquare *
108  (G4Log(2 * vF / fine_structure_const) - 0.5) / vF;
109  }
110  }
111 }
112 
113 G4double
115  const G4MaterialCutsCouple* couple)
116 {
117  return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy();
118 }
119 
120 G4double
122  const G4ParticleDefinition* p,
123  G4double kineticEnergy,
124  G4double maxEnergy)
125 {
126  if (!monopole) { SetParticle(p); }
127  G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
128  G4double cutEnergy = std::min(tmax, maxEnergy);
129  cutEnergy = std::max(LowEnergyLimit(), cutEnergy);
130  G4double tau = kineticEnergy / mass;
131  G4double gam = tau + 1.0;
132  G4double bg2 = tau * (tau + 2.0);
133  G4double beta2 = bg2 / (gam * gam);
134  G4double beta = sqrt(beta2);
135 
136  // low-energy asymptotic formula
137  //G4double dedx = dedxlim*beta*material->GetDensity();
138  G4double dedx = (*dedx0)[CurrentCouple()->GetIndex()] * beta;
139 
140  // above asymptotic
141  if (beta > betalow) {
142 
143  // high energy
144  if (beta >= betalim) {
145  dedx = ComputeDEDXAhlen(material, bg2, cutEnergy);
146 
147  } else {
148 
149  //G4double dedx1 = dedxlim*betalow*material->GetDensity();
150  G4double dedx1 = (*dedx0)[CurrentCouple()->GetIndex()] * betalow;
151  G4double dedx2 = ComputeDEDXAhlen(material, bg2lim, cutEnergy);
152 
153  // extrapolation between two formula
154  G4double kapa2 = beta - betalow;
155  G4double kapa1 = betalim - beta;
156  dedx = (kapa1 * dedx1 + kapa2 * dedx2) / (kapa1 + kapa2);
157  }
158  }
159  return dedx;
160 }
161 
162 G4double
164  G4double bg2,
165  G4double cutEnergy)
166 {
167  G4double eDensity = material->GetElectronDensity();
168  G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
169 
170  // Ahlen's formula for nonconductors, [1]p157, f(5.7)
171  G4double dedx =
172  0.5 * (G4Log(2.0 * electron_mass_c2 * bg2 * cutEnergy / (eexc * eexc)) - 1.0);//"Conventional" ionisation
173 // G4double dedx =
174 // 1.0 * (G4Log(2.0 * electron_mass_c2 * bg2 * cutEnergy / (eexc * eexc)));//Fryberger magneticon double ionisation
175 
176 
177  G4double k = 0; // Cross-section correction
178  if (nmpl >= 0.5) { k = 0.406; }
179  if (nmpl >= 1) { k = 0.346; }
180  if (nmpl >= 1.5) { k = 0.3; }
181  const G4double B[7] = { 0.0, 0.248, 0.672, 1.022, 1.243, 1.464, 1.685}; // Bloch correction
182  if (nmpl < 6)
183  dedx += 0.5 * k - B[int(floor(nmpl + 0.5))];
184 
185 
186  // density effect correction
187  G4double x = G4Log(bg2) / twoln10;
188  dedx -= material->GetIonisation()->DensityCorrection(x);
189 
190  // now compute the total ionization loss
191  dedx *= pi_hbarc2_over_mc2 * eDensity * chargeSquare;
192 
193  dedx = std::max(dedx, 0.0);
194  return dedx;
195 }
196 
197 G4double
199  const G4ParticleDefinition* p,
200  G4double kineticEnergy,
201  G4double cut,
202  G4double maxKinEnergy)
203 {
204  if (!monopole) { SetParticle(p); }
205  G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
206  G4double maxEnergy = std::min(tmax, maxKinEnergy);
207  G4double cutEnergy = std::max(LowEnergyLimit(), cut);
208  G4double cross = (cutEnergy < maxEnergy)
209  ? (1.0 / cutEnergy - 1.0 / maxEnergy) * twopi_mc2_rcl2 * chargeSquare : 0.0; //cross section used in ATLAS
210 // G4double cross = (cutEnergy < maxEnergy)
211 // ? (0.5 / cutEnergy - 0.5 / maxEnergy) * pi_hbarc2_over_mc2 * chargeSquare : 0.0; //one coming with GEANT4; causes 4 orders higher values
212  return cross;
213 }
214 
215 G4double
217  const G4ParticleDefinition* p,
218  G4double kineticEnergy,
219  G4double Z, G4double,
220  G4double cutEnergy,
221  G4double maxEnergy)
222 {
223  G4double cross =
224  Z * ComputeCrossSectionPerElectron(p, kineticEnergy, cutEnergy, maxEnergy);
225  return cross;
226 }
227 
228 void
229 G4mplIonisationWithDeltaModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
230  const G4MaterialCutsCouple*,
231  const G4DynamicParticle* dp,
232  G4double minKinEnergy,
233  G4double maxEnergy)
234 {
235  G4double kineticEnergy = dp->GetKineticEnergy();
236  G4double tmax = MaxSecondaryEnergy(dp->GetDefinition(), kineticEnergy);
237 
238  G4double maxKinEnergy = std::min(maxEnergy, tmax);
239  if (minKinEnergy >= maxKinEnergy) { return; }
240 
241  G4double totEnergy = kineticEnergy + mass;
242  G4double etot2 = totEnergy * totEnergy;
243  G4double beta2 = kineticEnergy * (kineticEnergy + 2.0 * mass) / etot2;
244 
245  // sampling without nuclear size effect
246  G4double q = G4UniformRand();
247  G4double deltaKinEnergy = minKinEnergy * maxKinEnergy
248  / (minKinEnergy * (1.0 - q) + maxKinEnergy * q);
249 
250  // delta-electron is produced
251  G4double totMomentum = totEnergy * sqrt(beta2);
252  G4double deltaMomentum =
253  sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0 * electron_mass_c2));
254  G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) /
255  (deltaMomentum * totMomentum);
256  cost = std::min(cost, 1.0);
257 
258  G4double sint = sqrt((1.0 - cost) * (1.0 + cost));
259 
260  G4double phi = twopi * G4UniformRand() ;
261 
262  G4ThreeVector deltaDirection(sint * cos(phi), sint * sin(phi), cost);
263  G4ThreeVector direction = dp->GetMomentumDirection();
264  deltaDirection.rotateUz(direction);
265 
266  // create G4DynamicParticle object for delta ray
267  G4DynamicParticle* delta =
268  new G4DynamicParticle(theElectron, deltaDirection, deltaKinEnergy);
269 
270  vdp->push_back(delta);
271 
272  // Change kinematics of primary particle
273  kineticEnergy -= deltaKinEnergy;
274  G4ThreeVector finalP = direction * totMomentum - deltaDirection * deltaMomentum;
275  finalP = finalP.unit();
276 
277  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
278  fParticleChange->SetProposedMomentumDirection(finalP);
279 }
280 
282  const G4MaterialCutsCouple* couple,
283  const G4DynamicParticle* dp,
284  G4double tmax,
285  G4double length,
286  G4double meanLoss)
287 {
288  G4double siga = Dispersion(couple->GetMaterial(), dp, tmax, length);
289  G4double loss = meanLoss;
290  siga = sqrt(siga);
291  G4double twomeanLoss = meanLoss + meanLoss;
292 
293  if (twomeanLoss < siga) {
294  G4double x;
295  do {
296  loss = twomeanLoss * G4UniformRand();
297  x = (loss - meanLoss) / siga;
298  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
299  } while (1.0 - 0.5 * x * x < G4UniformRand());
300  } else {
301  do {
302  loss = G4RandGauss::shoot(meanLoss, siga);
303  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
304  } while (0.0 > loss || loss > twomeanLoss);
305  }
306  return loss;
307 }
308 
309 G4double
310 G4mplIonisationWithDeltaModel::Dispersion(const G4Material* material,
311  const G4DynamicParticle* dp,
312  G4double tmax,
313  G4double length)
314 {
315  G4double siga = 0.0;
316  G4double tau = dp->GetKineticEnergy() / mass;
317  if (tau > 0.0) {
318  G4double electronDensity = material->GetElectronDensity();
319  G4double gam = tau + 1.0;
320  G4double invbeta2 = (gam * gam) / (tau * (tau + 2.0));
321  siga = (invbeta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
322  * electronDensity * chargeSquare;
323  }
324  return siga;
325 }
326 
327 G4double
329  G4double kinEnergy)
330 {
331  G4double tau = kinEnergy / mass;
332  return 2.0 * electron_mass_c2 * tau * (tau + 2.);
333 }
Belle2::Monopoles::G4mplIonisationWithDeltaModel::MinEnergyCut
virtual G4double MinEnergyCut(const G4ParticleDefinition *, const G4MaterialCutsCouple *couple) override
Threshold for zero value.
Definition: G4mplIonisationWithDeltaModel.cc:114
Belle2::Monopoles::G4mplIonisationWithDeltaModel::betalim
G4double betalim
Beta threshold for high energy (only Ahlen formula)
Definition: G4mplIonisationWithDeltaModel.h:219
Belle2::Monopoles::G4mplIonisationWithDeltaModel::beta2lim
G4double beta2lim
Square of betalim.
Definition: G4mplIonisationWithDeltaModel.h:220
Belle2::Monopoles::G4mplIonisationWithDeltaModel::ComputeDEDXAhlen
G4double ComputeDEDXAhlen(const G4Material *material, G4double bg2, G4double cut)
Caculate dedx based on extrapolated Ahlen formula.
Definition: G4mplIonisationWithDeltaModel.cc:163
Belle2::Monopoles::G4mplIonisationWithDeltaModel::SetParticle
void SetParticle(const G4ParticleDefinition *p)
Read definition of the monopole.
Definition: G4mplIonisationWithDeltaModel.cc:73
Belle2::Monopoles::G4mplIonisationWithDeltaModel::ComputeDEDXPerVolume
virtual G4double ComputeDEDXPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double kineticEnergy, G4double maxEnergy) override
G4VEmModel::ComputeDEDXPerVolume implementation.
Definition: G4mplIonisationWithDeltaModel.cc:121
Belle2::Monopoles::G4mplIonisationWithDeltaModel::magCharge
G4double magCharge
Monopole magnetic charge in e+ units.
Definition: G4mplIonisationWithDeltaModel.h:215
Belle2::Monopoles::G4mplIonisationWithDeltaModel::ComputeCrossSectionPerAtom
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *p, G4double kineticEnergy, G4double Z, G4double A, G4double cutEnergy, G4double maxEnergy) override
Compute cross section per atom for delta electrons emission.
Definition: G4mplIonisationWithDeltaModel.cc:216
Belle2::Monopoles::G4mplIonisationWithDeltaModel::SampleFluctuations
virtual G4double SampleFluctuations(const G4MaterialCutsCouple *couple, const G4DynamicParticle *dp, G4double tmax, G4double length, G4double meanLoss) override
Create fluctuations in the energies lost to a secondary delta electron.
Definition: G4mplIonisationWithDeltaModel.cc:281
Belle2::Monopoles::G4mplIonisationWithDeltaModel::dedxlim
G4double dedxlim
dedx limit in asymptotic formula, not used
Definition: G4mplIonisationWithDeltaModel.h:223
Belle2::Monopoles::G4mplIonisationWithDeltaModel::ComputeCrossSectionPerElectron
virtual G4double ComputeCrossSectionPerElectron(const G4ParticleDefinition *p, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
Compute cross section per electron for delta electrons emission.
Definition: G4mplIonisationWithDeltaModel.cc:198
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::Monopoles::G4mplIonisationWithDeltaModel::Initialise
virtual void Initialise(const G4ParticleDefinition *p, const G4DataVector &) override
G4VEmModel::Initialise implementation.
Definition: G4mplIonisationWithDeltaModel.cc:86
Belle2::Monopoles::G4mplIonisationWithDeltaModel::dedx0
static std::vector< G4double > * dedx0
Base dedx for each couple in current material.
Definition: G4mplIonisationWithDeltaModel.h:227
Belle2::Monopoles::G4mplIonisationWithDeltaModel::mass
G4double mass
Mass of the monopole.
Definition: G4mplIonisationWithDeltaModel.h:214
Belle2::Monopoles::G4mplIonisationWithDeltaModel::twoln10
G4double twoln10
log(100.0)
Definition: G4mplIonisationWithDeltaModel.h:217
Belle2::Monopoles::G4mplIonisationWithDeltaModel::SampleSecondaries
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *vdp, const G4MaterialCutsCouple *materialCutsCouple, const G4DynamicParticle *dp, G4double tmin, G4double maxEnergy) override
Create the sample of secondary delta electrons.
Definition: G4mplIonisationWithDeltaModel.cc:229
Belle2::Monopoles::G4mplIonisationWithDeltaModel::fParticleChange
G4ParticleChangeForLoss * fParticleChange
Pointer to ionising particle.
Definition: G4mplIonisationWithDeltaModel.h:212
Belle2::Monopoles::G4mplIonisationWithDeltaModel::theElectron
G4ParticleDefinition * theElectron
Electron definition.
Definition: G4mplIonisationWithDeltaModel.h:211
Belle2::Monopoles::G4mplIonisationWithDeltaModel::bg2lim
G4double bg2lim
(beta*gamma)^2 for betalim
Definition: G4mplIonisationWithDeltaModel.h:221
Belle2::Monopoles::G4mplIonisationWithDeltaModel::nmpl
G4double nmpl
Magn.
Definition: G4mplIonisationWithDeltaModel.h:225
Belle2::Monopoles::G4mplIonisationWithDeltaModel::Dispersion
virtual G4double Dispersion(const G4Material *material, const G4DynamicParticle *, G4double tmax, G4double length) override
Calculate dispersion.
Definition: G4mplIonisationWithDeltaModel.cc:310
Belle2::Monopoles::G4mplIonisationWithDeltaModel::MaxSecondaryEnergy
virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition *p, G4double kinEnergy) override
Calculate maximum energy available for secondary particle emission.
Definition: G4mplIonisationWithDeltaModel.cc:328
Belle2::Monopoles::G4mplIonisationWithDeltaModel::betalow
G4double betalow
Beta threshold for low asymptotic.
Definition: G4mplIonisationWithDeltaModel.h:218
Belle2::Monopoles::G4mplIonisationWithDeltaModel::~G4mplIonisationWithDeltaModel
virtual ~G4mplIonisationWithDeltaModel()
Destructor.
Definition: G4mplIonisationWithDeltaModel.cc:68
Belle2::Monopoles::G4mplIonisationWithDeltaModel::chargeSquare
G4double chargeSquare
Square of magnetic charge in units of Dirac charge.
Definition: G4mplIonisationWithDeltaModel.h:222
Belle2::Monopoles::G4mplIonisationWithDeltaModel::pi_hbarc2_over_mc2
G4double pi_hbarc2_over_mc2
Convenient constants combination with mass.
Definition: G4mplIonisationWithDeltaModel.h:224
Belle2::Monopoles::G4mplIonisationWithDeltaModel::monopole
const G4ParticleDefinition * monopole
Monopole definition.
Definition: G4mplIonisationWithDeltaModel.h:210