21 #include <simulation/monopoles/G4mplIonisationWithDeltaModel.h>
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>
34 #include <framework/logging/Logger.h>
38 using namespace Belle2::Monopoles;
39 using namespace CLHEP;
41 std::vector<G4double>* G4mplIonisationWithDeltaModel::dedx0 =
nullptr;
43 G4mplIonisationWithDeltaModel::G4mplIonisationWithDeltaModel(G4double mCharge,
45 : G4VEmModel(nam), G4VEmFluctuationModel(nam),
50 beta2lim(betalim * betalim),
51 bg2lim(beta2lim * (1.0 + beta2lim)),
52 pi_hbarc2_over_mc2(pi * hbarc * hbarc / electron_mass_c2)
61 B2INFO(
"Monopole ionisation model with d-electron production, Gmag= " <<
magCharge / eplus);
63 B2WARNING(
"Monopole charge Gmag= " <<
magCharge / eplus <<
"e not reasonable. Please choose a value smaller than 411e.");
70 if (IsMaster()) {
delete dedx0; }
80 std::max(HighEnergyLimit(), 10 *
mass * (1. / sqrt(1. -
beta2lim) - 1.));
81 SetLowEnergyLimit(emin);
82 SetHighEnergyLimit(emax);
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();
101 for (G4int i = 0; i < numOfCouples; ++i) {
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);
108 (G4Log(2 * vF / fine_structure_const) - 0.5) / vF;
115 const G4MaterialCutsCouple* couple)
117 return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy();
122 const G4ParticleDefinition* p,
123 G4double 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);
138 G4double dedx = (*dedx0)[CurrentCouple()->GetIndex()] * beta;
150 G4double dedx1 = (*dedx0)[CurrentCouple()->GetIndex()] *
betalow;
154 G4double kapa2 = beta -
betalow;
155 G4double kapa1 =
betalim - beta;
156 dedx = (kapa1 * dedx1 + kapa2 * dedx2) / (kapa1 + kapa2);
167 G4double eDensity = material->GetElectronDensity();
168 G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
172 0.5 * (G4Log(2.0 * electron_mass_c2 * bg2 * cutEnergy / (eexc * eexc)) - 1.0);
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};
183 dedx += 0.5 * k - B[int(floor(
nmpl + 0.5))];
187 G4double x = G4Log(bg2) /
twoln10;
188 dedx -= material->GetIonisation()->DensityCorrection(x);
193 dedx = std::max(dedx, 0.0);
199 const G4ParticleDefinition* p,
200 G4double kineticEnergy,
202 G4double maxKinEnergy)
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;
217 const G4ParticleDefinition* p,
218 G4double kineticEnergy,
219 G4double Z, G4double,
230 const G4MaterialCutsCouple*,
231 const G4DynamicParticle* dp,
232 G4double minKinEnergy,
235 G4double kineticEnergy = dp->GetKineticEnergy();
238 G4double maxKinEnergy = std::min(maxEnergy, tmax);
239 if (minKinEnergy >= maxKinEnergy) {
return; }
241 G4double totEnergy = kineticEnergy +
mass;
242 G4double etot2 = totEnergy * totEnergy;
243 G4double beta2 = kineticEnergy * (kineticEnergy + 2.0 *
mass) / etot2;
246 G4double q = G4UniformRand();
247 G4double deltaKinEnergy = minKinEnergy * maxKinEnergy
248 / (minKinEnergy * (1.0 - q) + maxKinEnergy * q);
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);
258 G4double sint = sqrt((1.0 - cost) * (1.0 + cost));
260 G4double phi = twopi * G4UniformRand() ;
262 G4ThreeVector deltaDirection(sint * cos(phi), sint * sin(phi), cost);
263 G4ThreeVector direction = dp->GetMomentumDirection();
264 deltaDirection.rotateUz(direction);
267 G4DynamicParticle* delta =
268 new G4DynamicParticle(
theElectron, deltaDirection, deltaKinEnergy);
270 vdp->push_back(delta);
273 kineticEnergy -= deltaKinEnergy;
274 G4ThreeVector finalP = direction * totMomentum - deltaDirection * deltaMomentum;
275 finalP = finalP.unit();
282 const G4MaterialCutsCouple* couple,
283 const G4DynamicParticle* dp,
288 G4double siga =
Dispersion(couple->GetMaterial(), dp, tmax, length);
289 G4double loss = meanLoss;
291 G4double twomeanLoss = meanLoss + meanLoss;
293 if (twomeanLoss < siga) {
296 loss = twomeanLoss * G4UniformRand();
297 x = (loss - meanLoss) / siga;
299 }
while (1.0 - 0.5 * x * x < G4UniformRand());
302 loss = G4RandGauss::shoot(meanLoss, siga);
304 }
while (0.0 > loss || loss > twomeanLoss);
311 const G4DynamicParticle* dp,
316 G4double tau = dp->GetKineticEnergy() /
mass;
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
331 G4double tau = kinEnergy /
mass;
332 return 2.0 * electron_mass_c2 * tau * (tau + 2.);