File indexing completed on 2023-03-17 11:00:47
0001 #include "FastSimulation/MaterialEffects/interface/PairProductionSimulator.h"
0002 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0003 #include "FastSimulation/Particle/interface/makeParticle.h"
0004
0005 #include <cmath>
0006
0007 PairProductionSimulator::PairProductionSimulator(double photonEnergyCut) {
0008
0009 photonEnergy = std::max(0.100, photonEnergyCut);
0010 }
0011
0012 void PairProductionSimulator::compute(ParticlePropagator& Particle, RandomEngineAndDistribution const* random) {
0013 double eGamma = Particle.particle().e();
0014
0015
0016 if (eGamma >= photonEnergy) {
0017
0018
0019
0020
0021
0022 if (-std::log(random->flatShoot()) <= (7. / 9.) * radLengths) {
0023 double xe = 0;
0024 double xm = eMass() / eGamma;
0025 double weight = 0.;
0026
0027
0028 do {
0029 xe = random->flatShoot() * (1. - 2. * xm) + xm;
0030 weight = 1. - 4. / 3. * xe * (1. - xe);
0031 } while (weight < random->flatShoot());
0032
0033 double eElectron = xe * eGamma;
0034 double tElectron = eElectron - eMass();
0035 double pElectron = std::sqrt(std::max((eElectron + eMass()) * tElectron, 0.));
0036
0037 double ePositron = eGamma - eElectron;
0038 double tPositron = ePositron - eMass();
0039 double pPositron = std::sqrt((ePositron + eMass()) * tPositron);
0040
0041
0042 double phi = random->flatShoot() * 2. * M_PI;
0043 double sphi = std::sin(phi);
0044 double cphi = std::cos(phi);
0045
0046 double stheta1, stheta2, ctheta1, ctheta2;
0047
0048 if (eElectron > ePositron) {
0049 double theta1 = gbteth(eElectron, eMass(), xe, random) * eMass() / eElectron;
0050 stheta1 = std::sin(theta1);
0051 ctheta1 = std::cos(theta1);
0052 stheta2 = stheta1 * pElectron / pPositron;
0053 ctheta2 = std::sqrt(std::max(0., 1.0 - (stheta2 * stheta2)));
0054 } else {
0055 double theta2 = gbteth(ePositron, eMass(), xe, random) * eMass() / ePositron;
0056 stheta2 = std::sin(theta2);
0057 ctheta2 = std::cos(theta2);
0058 stheta1 = stheta2 * pPositron / pElectron;
0059 ctheta1 = std::sqrt(std::max(0., 1.0 - (stheta1 * stheta1)));
0060 }
0061
0062 double chi = Particle.particle().theta();
0063 double psi = Particle.particle().phi();
0064 RawParticle::RotationZ rotZ(psi);
0065 RawParticle::RotationY rotY(chi);
0066
0067 _theUpdatedState.reserve(2);
0068 _theUpdatedState.clear();
0069
0070
0071 _theUpdatedState.emplace_back(makeParticle(
0072 Particle.particleDataTable(),
0073 +11,
0074 XYZTLorentzVector(pElectron * stheta1 * cphi, pElectron * stheta1 * sphi, pElectron * ctheta1, eElectron)));
0075 _theUpdatedState[0].rotate(rotY);
0076 _theUpdatedState[0].rotate(rotZ);
0077
0078
0079 _theUpdatedState.emplace_back(makeParticle(
0080 Particle.particleDataTable(),
0081 -11,
0082 XYZTLorentzVector(-pPositron * stheta2 * cphi, -pPositron * stheta2 * sphi, pPositron * ctheta2, ePositron)));
0083 _theUpdatedState[1].rotate(rotY);
0084 _theUpdatedState[1].rotate(rotZ);
0085 }
0086 }
0087 }
0088
0089 double PairProductionSimulator::gbteth(double ener,
0090 double partm,
0091 double efrac,
0092 RandomEngineAndDistribution const* random) {
0093 const double alfa = 0.625;
0094
0095 double d = 0.13 * (0.8 + 1.3 / theZ()) * (100.0 + (1.0 / ener)) * (1.0 + efrac);
0096 double w1 = 9.0 / (9.0 + d);
0097 double umax = ener * M_PI / partm;
0098 double u;
0099
0100 do {
0101 double beta;
0102 if (random->flatShoot() <= w1)
0103 beta = alfa;
0104 else
0105 beta = 3.0 * alfa;
0106 u = -(std::log(random->flatShoot() * random->flatShoot())) / beta;
0107 } while (u >= umax);
0108
0109 return u;
0110 }