File indexing completed on 2023-03-17 11:00:46
0001
0002 #include "FastSimulation/MaterialEffects/interface/BremsstrahlungSimulator.h"
0003 #include "FastSimulation/Particle/interface/makeParticle.h"
0004 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0005
0006 #include <cmath>
0007
0008 BremsstrahlungSimulator::BremsstrahlungSimulator(double photonEnergyCut, double photonFractECut) {
0009
0010 photonEnergy = photonEnergyCut;
0011 photonFractE = photonFractECut;
0012 }
0013
0014 void BremsstrahlungSimulator::compute(ParticlePropagator& Particle, RandomEngineAndDistribution const* random) {
0015
0016
0017
0018
0019
0020 if (radLengths > 4.)
0021 Particle.particle().setMomentum(0., 0., 0., 0.);
0022
0023
0024 if (Particle.particle().e() < photonEnergy)
0025 return;
0026 xmin = std::max(photonEnergy / Particle.particle().e(), photonFractE);
0027 if (xmin >= 1. || xmin <= 0.)
0028 return;
0029
0030 double bremProba =
0031 radLengths * (4. / 3. * std::log(1. / xmin) - 4. / 3. * (1. - xmin) + 1. / 2. * (1. - xmin * xmin));
0032
0033
0034 unsigned int nPhotons = poisson(bremProba, random);
0035 _theUpdatedState.reserve(nPhotons);
0036
0037 if (!nPhotons)
0038 return;
0039
0040
0041 double chi = Particle.particle().theta();
0042 double psi = Particle.particle().phi();
0043 RawParticle::RotationZ rotZ(psi);
0044 RawParticle::RotationY rotY(chi);
0045
0046
0047 for (unsigned int i = 0; i < nPhotons; ++i) {
0048
0049 if (Particle.particle().e() < photonEnergy)
0050 break;
0051
0052
0053 RawParticle thePhoton = makeParticle(Particle.particleDataTable(), 22, brem(Particle, random));
0054 thePhoton.rotate(rotY);
0055 thePhoton.rotate(rotZ);
0056 _theUpdatedState.push_back(thePhoton);
0057
0058
0059 Particle.particle().momentum() -= thePhoton.momentum();
0060 }
0061 }
0062
0063 XYZTLorentzVector BremsstrahlungSimulator::brem(ParticlePropagator& pp,
0064 RandomEngineAndDistribution const* random) const {
0065
0066
0067
0068 double emass = 0.0005109990615;
0069 double xp = 0;
0070 double weight = 0.;
0071
0072 do {
0073 xp = xmin * std::exp(-std::log(xmin) * random->flatShoot());
0074 weight = 1. - xp + 3. / 4. * xp * xp;
0075 } while (weight < random->flatShoot());
0076
0077
0078
0079
0080
0081 const double phi = random->flatShoot() * 2 * M_PI;
0082
0083 const double theta = gbteth(pp.particle().e(), emass, xp, random) * emass / pp.particle().e();
0084
0085
0086 double stheta = std::sin(theta);
0087 double ctheta = std::cos(theta);
0088 double sphi = std::sin(phi);
0089 double cphi = std::cos(phi);
0090
0091 return xp * pp.particle().e() * XYZTLorentzVector(stheta * cphi, stheta * sphi, ctheta, 1.);
0092 }
0093
0094 double BremsstrahlungSimulator::gbteth(const double ener,
0095 const double partm,
0096 const double efrac,
0097 RandomEngineAndDistribution const* random) const {
0098 const double alfa = 0.625;
0099
0100 const double d = 0.13 * (0.8 + 1.3 / theZ()) * (100.0 + (1.0 / ener)) * (1.0 + efrac);
0101 const double w1 = 9.0 / (9.0 + d);
0102 const double umax = ener * M_PI / partm;
0103 double u;
0104
0105 do {
0106 double beta = (random->flatShoot() <= w1) ? alfa : 3.0 * alfa;
0107 u = -std::log(random->flatShoot() * random->flatShoot()) / beta;
0108 } while (u >= umax);
0109
0110 return u;
0111 }
0112
0113 unsigned int BremsstrahlungSimulator::poisson(double ymu, RandomEngineAndDistribution const* random) {
0114 unsigned int n = 0;
0115 double prob = std::exp(-ymu);
0116 double proba = prob;
0117 double x = random->flatShoot();
0118
0119 while (proba <= x) {
0120 prob *= ymu / double(++n);
0121 proba += prob;
0122 }
0123
0124 return n;
0125 }