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 #include <TMath.h>
0019 #include <TF1.h>
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033 namespace fastsim {
0034
0035
0036
0037
0038
0039 class MuonBremsstrahlung : public InteractionModel {
0040 public:
0041
0042 MuonBremsstrahlung(const std::string &name, const edm::ParameterSet &cfg);
0043
0044
0045 ~MuonBremsstrahlung() override { ; };
0046
0047
0048
0049
0050
0051
0052
0053
0054 void interact(Particle &particle,
0055 const SimplifiedGeometry &layer,
0056 std::vector<std::unique_ptr<Particle> > &secondaries,
0057 const RandomEngineAndDistribution &random) override;
0058
0059 private:
0060
0061
0062
0063
0064
0065
0066
0067 math::XYZTLorentzVector brem(Particle &particle, double xmin, const RandomEngineAndDistribution &random) const;
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077 double gbteth(const double ener,
0078 const double partm,
0079 const double efrac,
0080 const RandomEngineAndDistribution &random) const;
0081
0082
0083 static double PetrukhinFunc(double *x, double *p);
0084
0085 TF1 *Petrfunc;
0086 double minPhotonEnergy_;
0087 double minPhotonEnergyFraction_;
0088 double density_;
0089 double radLenInCm_;
0090 double A_;
0091 double Z_;
0092 };
0093 }
0094
0095 fastsim::MuonBremsstrahlung::MuonBremsstrahlung(const std::string &name, const edm::ParameterSet &cfg)
0096 : fastsim::InteractionModel(name) {
0097
0098 minPhotonEnergy_ = cfg.getParameter<double>("minPhotonEnergy");
0099 minPhotonEnergyFraction_ = cfg.getParameter<double>("minPhotonEnergyFraction");
0100
0101 A_ = cfg.getParameter<double>("A");
0102 Z_ = cfg.getParameter<double>("Z");
0103 density_ = cfg.getParameter<double>("density");
0104 radLenInCm_ = cfg.getParameter<double>("radLen");
0105 }
0106
0107 void fastsim::MuonBremsstrahlung::interact(fastsim::Particle &particle,
0108 const SimplifiedGeometry &layer,
0109 std::vector<std::unique_ptr<fastsim::Particle> > &secondaries,
0110 const RandomEngineAndDistribution &random) {
0111
0112 if (std::abs(particle.pdgId()) != 13) {
0113 return;
0114 }
0115
0116 double radLengths = layer.getThickness(particle.position(), particle.momentum());
0117
0118
0119
0120 if (radLengths < 1E-10) {
0121 return;
0122 }
0123
0124
0125
0126
0127 if (radLengths > 4.) {
0128 particle.momentum().SetXYZT(0., 0., 0., 0.);
0129 return;
0130 }
0131
0132
0133 if (particle.momentum().E() - particle.momentum().mass() < minPhotonEnergy_) {
0134 return;
0135 }
0136
0137
0138 double xmin = std::max(minPhotonEnergy_ / particle.momentum().E(), minPhotonEnergyFraction_);
0139
0140 if (xmin >= 1. || xmin <= 0.) {
0141 return;
0142 }
0143
0144
0145 double xmax = 1.;
0146
0147
0148 Petrfunc = new TF1("Petrfunc", PetrukhinFunc, xmin, xmax, 3);
0149
0150 Petrfunc->SetParameters(particle.momentum().E(), A_, Z_);
0151
0152
0153
0154 double distance = radLengths * radLenInCm_;
0155
0156
0157
0158 double bremProba = density_ * distance * (fastsim::Constants::NA / A_) * (Petrfunc->Integral(xmin, xmax));
0159 if (bremProba < 1E-10) {
0160 return;
0161 }
0162
0163
0164 unsigned int nPhotons = random.poissonShoot(bremProba);
0165 if (nPhotons == 0) {
0166 return;
0167 }
0168
0169
0170 double theta = particle.momentum().Theta();
0171 double phi = particle.momentum().Phi();
0172
0173
0174 for (unsigned int i = 0; i < nPhotons; ++i) {
0175
0176 math::XYZTLorentzVector photonMom = brem(particle, xmin, random);
0177
0178
0179 if (particle.momentum().E() - particle.momentum().mass() < photonMom.E())
0180 break;
0181
0182
0183 photonMom = ROOT::Math::RotationZ(phi) * (ROOT::Math::RotationY(theta) * photonMom);
0184
0185
0186 secondaries.emplace_back(new fastsim::Particle(22, particle.position(), photonMom));
0187
0188
0189 particle.momentum() -= photonMom;
0190 }
0191 }
0192
0193 math::XYZTLorentzVector fastsim::MuonBremsstrahlung::brem(fastsim::Particle &particle,
0194 double xmin,
0195 const RandomEngineAndDistribution &random) const {
0196
0197
0198 double xp = Petrfunc->GetRandom();
0199
0200
0201
0202
0203
0204 const double phi = random.flatShoot() * 2. * M_PI;
0205
0206 const double theta = gbteth(particle.momentum().E(), fastsim::Constants::muMass, xp, random) *
0207 fastsim::Constants::muMass / particle.momentum().E();
0208
0209
0210 double stheta = std::sin(theta);
0211 double ctheta = std::cos(theta);
0212 double sphi = std::sin(phi);
0213 double cphi = std::cos(phi);
0214
0215 return xp * particle.momentum().E() * math::XYZTLorentzVector(stheta * cphi, stheta * sphi, ctheta, 1.);
0216 }
0217
0218 double fastsim::MuonBremsstrahlung::gbteth(const double ener,
0219 const double partm,
0220 const double efrac,
0221 const RandomEngineAndDistribution &random) const {
0222
0223
0224
0225
0226 const double alfa = 0.625;
0227
0228 const double d = 0.13 * (0.8 + 1.3 / Z_) * (100.0 + (1.0 / ener)) * (1.0 + efrac);
0229 const double w1 = 9.0 / (9.0 + d);
0230 const double umax = ener * M_PI / partm;
0231 double u;
0232
0233 do {
0234 double beta = (random.flatShoot() <= w1) ? alfa : 3.0 * alfa;
0235 u = -std::log(random.flatShoot() * random.flatShoot()) / beta;
0236 } while (u >= umax);
0237
0238 return u;
0239 }
0240
0241 double fastsim::MuonBremsstrahlung::PetrukhinFunc(double *x, double *p) {
0242
0243 double nu = x[0];
0244
0245
0246 double E = p[0];
0247 double A = p[1];
0248 double Z = p[2];
0249
0250
0251
0252
0253
0254
0255
0256
0257
0258
0259
0260
0261
0262
0263
0264
0265
0266
0267
0268
0269
0270
0271
0272
0273
0274
0275
0276
0277
0278
0279
0280
0281
0282
0283
0284
0285 double B = 183.;
0286 double Bl = 1429.;
0287 double ee = 1.64872;
0288 double Z13 = pow(Z, -1. / 3.);
0289 double Z23 = pow(Z, -2. / 3.);
0290
0291
0292 double emass = 0.0005109990615;
0293 double mumass = 0.105658367;
0294
0295 double alpha = 0.00729735;
0296 double constant = 1.85736e-30;
0297
0298
0299 if (nu * E >= E - mumass)
0300 return 0;
0301
0302 double Dn = 1.54 * (pow(A, 0.27));
0303 double Dnl = pow(Dn, (1. - 1. / Z));
0304
0305 double delta = (mumass * mumass * nu) / (2. * E * (1. - nu));
0306
0307 double Phi_n = TMath::Log(B * Z13 * (mumass + delta * (Dnl * ee - 2)) / (Dnl * (emass + delta * ee * B * Z13)));
0308 if (Phi_n < 0)
0309 Phi_n = 0;
0310
0311 double Phi_e =
0312 TMath::Log((Bl * Z23 * mumass) / (1. + delta * mumass / (emass * emass * ee)) / (emass + delta * ee * Bl * Z23));
0313 if (Phi_e < 0 || nu * E >= E / (1. + (mumass * mumass / (2. * emass * E))))
0314 Phi_e = 0;
0315
0316
0317 double f = 16. / 3. * alpha * constant * Z * (Z * Phi_n + Phi_e) * (1. / nu) * (1. - nu + 0.75 * nu * nu);
0318
0319 return f;
0320 }
0321
0322 DEFINE_EDM_PLUGIN(fastsim::InteractionModelFactory, fastsim::MuonBremsstrahlung, "fastsim::MuonBremsstrahlung");