File indexing completed on 2024-04-06 12:11:21
0001 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0002 #include "FastSimulation/SimplifiedGeometryPropagator/interface/Particle.h"
0003 #include "FastSimulation/SimplifiedGeometryPropagator/interface/SimplifiedGeometry.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006
0007 #include <cmath>
0008 #include <memory>
0009
0010 #include <Math/RotationY.h>
0011 #include <Math/RotationZ.h>
0012
0013 #include "FastSimulation/SimplifiedGeometryPropagator/interface/InteractionModelFactory.h"
0014 #include "FastSimulation/SimplifiedGeometryPropagator/interface/InteractionModel.h"
0015 #include "FastSimulation/SimplifiedGeometryPropagator/interface/Constants.h"
0016 #include "DataFormats/Math/interface/LorentzVector.h"
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026 namespace fastsim {
0027
0028
0029
0030
0031 class PairProduction : public InteractionModel {
0032 public:
0033
0034 PairProduction(const std::string& name, const edm::ParameterSet& cfg);
0035
0036
0037 ~PairProduction() override { ; };
0038
0039
0040
0041
0042
0043
0044
0045
0046 void interact(fastsim::Particle& particle,
0047 const SimplifiedGeometry& layer,
0048 std::vector<std::unique_ptr<fastsim::Particle> >& secondaries,
0049 const RandomEngineAndDistribution& random) override;
0050
0051 private:
0052
0053
0054
0055
0056
0057
0058
0059
0060 double gbteth(double ener, double partm, double efrac, const RandomEngineAndDistribution& random) const;
0061
0062 double minPhotonEnergy_;
0063 double Z_;
0064 };
0065 }
0066
0067 fastsim::PairProduction::PairProduction(const std::string& name, const edm::ParameterSet& cfg)
0068 : fastsim::InteractionModel(name) {
0069
0070 minPhotonEnergy_ = cfg.getParameter<double>("photonEnergyCut");
0071
0072 Z_ = cfg.getParameter<double>("Z");
0073 }
0074
0075 void fastsim::PairProduction::interact(fastsim::Particle& particle,
0076 const SimplifiedGeometry& layer,
0077 std::vector<std::unique_ptr<fastsim::Particle> >& secondaries,
0078 const RandomEngineAndDistribution& random) {
0079 double eGamma = particle.momentum().e();
0080
0081
0082
0083 if (particle.pdgId() != 22) {
0084 return;
0085 }
0086
0087 double radLengths = layer.getThickness(particle.position(), particle.momentum());
0088
0089
0090
0091 if (radLengths < 1E-10) {
0092 return;
0093 }
0094
0095
0096
0097
0098 if (eGamma < minPhotonEnergy_) {
0099 return;
0100 }
0101
0102
0103
0104
0105 if (-std::log(random.flatShoot()) > (7. / 9.) * radLengths) {
0106 return;
0107 }
0108
0109 double xe = 0;
0110 double eMass = fastsim::Constants::eMass;
0111 double xm = eMass / eGamma;
0112 double weight = 0.;
0113
0114
0115 do {
0116 xe = random.flatShoot() * (1. - 2. * xm) + xm;
0117 weight = 1. - 4. / 3. * xe * (1. - xe);
0118 } while (weight < random.flatShoot());
0119
0120
0121 double eElectron = xe * eGamma;
0122 double tElectron = eElectron - eMass;
0123 double pElectron = std::sqrt(std::max((eElectron + eMass) * tElectron, 0.));
0124
0125
0126 double ePositron = eGamma - eElectron;
0127 double tPositron = ePositron - eMass;
0128 double pPositron = std::sqrt((ePositron + eMass) * tPositron);
0129
0130
0131 double phi = random.flatShoot() * 2. * M_PI;
0132 double sphi = std::sin(phi);
0133 double cphi = std::cos(phi);
0134
0135 double stheta1, stheta2, ctheta1, ctheta2;
0136
0137 if (eElectron > ePositron) {
0138 double theta1 = gbteth(eElectron, eMass, xe, random) * eMass / eElectron;
0139 stheta1 = std::sin(theta1);
0140 ctheta1 = std::cos(theta1);
0141 stheta2 = stheta1 * pElectron / pPositron;
0142 ctheta2 = std::sqrt(std::max(0., 1.0 - (stheta2 * stheta2)));
0143 } else {
0144 double theta2 = gbteth(ePositron, eMass, xe, random) * eMass / ePositron;
0145 stheta2 = std::sin(theta2);
0146 ctheta2 = std::cos(theta2);
0147 stheta1 = stheta2 * pPositron / pElectron;
0148 ctheta1 = std::sqrt(std::max(0., 1.0 - (stheta1 * stheta1)));
0149 }
0150
0151
0152 double thetaLab = particle.momentum().Theta();
0153 double phiLab = particle.momentum().Phi();
0154
0155
0156 secondaries.emplace_back(new fastsim::Particle(
0157 11,
0158 particle.position(),
0159 math::XYZTLorentzVector(pElectron * stheta1 * cphi, pElectron * stheta1 * sphi, pElectron * ctheta1, eElectron)));
0160 secondaries.back()->momentum() =
0161 ROOT::Math::RotationZ(phiLab) * (ROOT::Math::RotationY(thetaLab) * secondaries.back()->momentum());
0162
0163
0164 secondaries.emplace_back(new fastsim::Particle(
0165 -11,
0166 particle.position(),
0167 math::XYZTLorentzVector(
0168 -pPositron * stheta2 * cphi, -pPositron * stheta2 * sphi, pPositron * ctheta2, ePositron)));
0169 secondaries.back()->momentum() =
0170 ROOT::Math::RotationZ(phiLab) * (ROOT::Math::RotationY(thetaLab) * secondaries.back()->momentum());
0171
0172
0173 particle.momentum().SetXYZT(0., 0., 0., 0.);
0174 }
0175
0176 double fastsim::PairProduction::gbteth(const double ener,
0177 const double partm,
0178 const double efrac,
0179 const RandomEngineAndDistribution& random) const {
0180
0181
0182
0183
0184 const double alfa = 0.625;
0185
0186 const double d = 0.13 * (0.8 + 1.3 / Z_) * (100.0 + (1.0 / ener)) * (1.0 + efrac);
0187 const double w1 = 9.0 / (9.0 + d);
0188 const double umax = ener * M_PI / partm;
0189 double u;
0190
0191 do {
0192 double beta = (random.flatShoot() <= w1) ? alfa : 3.0 * alfa;
0193 u = -std::log(random.flatShoot() * random.flatShoot()) / beta;
0194 } while (u >= umax);
0195
0196 return u;
0197 }
0198
0199 DEFINE_EDM_PLUGIN(fastsim::InteractionModelFactory, fastsim::PairProduction, "fastsim::PairProduction");