File indexing completed on 2024-04-06 12:11:16
0001
0002
0003
0004
0005
0006
0007 #include "FastSimulation/MaterialEffects/interface/MuonBremsstrahlungSimulator.h"
0008 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "FastSimulation/MaterialEffects/interface/PetrukhinModel.h"
0011 #include "FastSimulation/Particle/interface/makeParticle.h"
0012
0013 #include <cmath>
0014 #include <string>
0015 #include <iostream>
0016
0017 #include "TF1.h"
0018
0019
0020 MuonBremsstrahlungSimulator::MuonBremsstrahlungSimulator(
0021 double A, double Z, double density, double radLen, double photonEnergyCut, double photonFractECut)
0022 : MaterialEffectsSimulator(A, Z, density, radLen) {
0023
0024 photonEnergy = photonEnergyCut;
0025 photonFractE = photonFractECut;
0026 d = 0.;
0027 LogDebug("MuonBremsstrahlungSimulator") << "Starting the MuonBremsstrahlungSimulator" << std::endl;
0028 }
0029
0030
0031 void MuonBremsstrahlungSimulator::compute(ParticlePropagator& Particle, RandomEngineAndDistribution const* random) {
0032 double NA = 6.022e+23;
0033
0034 if (radLengths > 4.) {
0035 Particle.particle().setMomentum(0., 0., 0., 0.);
0036 deltaPMuon.SetXYZT(0., 0., 0., 0.);
0037 brem_photon.SetXYZT(0., 0., 0., 0.);
0038 }
0039
0040
0041
0042 double EMuon = Particle.particle().e();
0043 if (EMuon < photonEnergy)
0044 return;
0045 xmin = std::max(photonEnergy / EMuon, photonFractE);
0046
0047
0048 if (xmin >= 1. || xmin <= 0.)
0049 return;
0050
0051 xmax = 1.;
0052 npar = 3;
0053
0054
0055 f1 = new TF1("f1", PetrukhinFunc, xmin, xmax, npar);
0056
0057 f1->SetParameters(EMuon, A, Z);
0058
0059
0060
0061
0062 d = radLengths * radLen;
0063
0064 bremProba = density * d * (NA / A) * (f1->Integral(0., 1.));
0065
0066
0067 unsigned int nPhotons = random->poissonShoot(bremProba);
0068 _theUpdatedState.reserve(nPhotons);
0069
0070 if (!nPhotons)
0071 return;
0072
0073
0074 double chi = Particle.particle().theta();
0075 double psi = Particle.particle().phi();
0076 RawParticle::RotationZ rotZ(psi);
0077 RawParticle::RotationY rotY(chi);
0078
0079
0080 for (unsigned int i = 0; i < nPhotons; ++i) {
0081
0082 if (Particle.particle().e() < photonEnergy)
0083 break;
0084 LogDebug("MuonBremsstrahlungSimulator") << "MuonBremsstrahlungSimulator parameters:" << std::endl;
0085 LogDebug("MuonBremsstrahlungSimulator") << "xmin-> " << xmin << std::endl;
0086 LogDebug("MuonBremsstrahlungSimulator") << "Atomic Weight-> " << A << std::endl;
0087 LogDebug("MuonBremsstrahlungSimulator") << "Density-> " << density << std::endl;
0088 LogDebug("MuonBremsstrahlungSimulator") << "Distance-> " << d << std::endl;
0089 LogDebug("MuonBremsstrahlungSimulator") << "bremProba->" << bremProba << std::endl;
0090 LogDebug("MuonBremsstrahlungSimulator") << "nPhotons->" << nPhotons << std::endl;
0091 LogDebug("MuonBremsstrahlungSimulator") << " Muon_Energy-> " << EMuon << std::endl;
0092 LogDebug("MuonBremsstrahlungSimulator") << "X0-> " << radLen << std::endl;
0093 LogDebug("MuonBremsstrahlungSimulator") << " radLengths-> " << radLengths << std::endl;
0094
0095
0096 RawParticle thePhoton = makeParticle(Particle.particleDataTable(), 22, brem(Particle, random));
0097 if (thePhoton.E() > 0.) {
0098 thePhoton.rotate(rotY);
0099 thePhoton.rotate(rotZ);
0100
0101 _theUpdatedState.push_back(thePhoton);
0102
0103
0104 deltaPMuon = Particle.particle().momentum() -= thePhoton.momentum();
0105
0106 brem_photon.SetXYZT(thePhoton.Px(), thePhoton.Py(), thePhoton.Pz(), thePhoton.E());
0107
0108 LogDebug("MuonBremsstrahlungSimulator") << " Muon Bremsstrahlung: photon_energy-> " << thePhoton.E() << std::endl;
0109 LogDebug("MuonBremsstrahlungSimulator") << "photon_px->" << thePhoton.Px() << std::endl;
0110 LogDebug("MuonBremsstrahlungSimulator") << "photon_py->" << thePhoton.Py() << std::endl;
0111 LogDebug("MuonBremsstrahlungSimulator") << "photon_pz->" << thePhoton.Pz() << std::endl;
0112 }
0113 }
0114 }
0115
0116 XYZTLorentzVector MuonBremsstrahlungSimulator::brem(ParticlePropagator& pp,
0117 RandomEngineAndDistribution const* random) const {
0118
0119
0120 double mumass = 0.105658367;
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131
0132 double xp = f1->GetRandom();
0133 LogDebug("MuonBremsstrahlungSimulator") << "MuonBremsstrahlungSimulator: xp->" << xp << std::endl;
0134 std::cout << "MuonBremsstrahlungSimulator: xp->" << xp << std::endl;
0135
0136
0137
0138
0139
0140 const double phi = random->flatShoot() * 2 * M_PI;
0141
0142 const double theta = gbteth(pp.particle().e(), mumass, xp, random) * mumass / pp.particle().e();
0143
0144
0145 double stheta = std::sin(theta);
0146 double ctheta = std::cos(theta);
0147 double sphi = std::sin(phi);
0148 double cphi = std::cos(phi);
0149
0150 return xp * pp.particle().e() * XYZTLorentzVector(stheta * cphi, stheta * sphi, ctheta, 1.);
0151 }
0152
0153 double MuonBremsstrahlungSimulator::gbteth(const double ener,
0154 const double partm,
0155 const double efrac,
0156 RandomEngineAndDistribution const* random) const {
0157 const double alfa = 0.625;
0158
0159 const double d = 0.13 * (0.8 + 1.3 / theZ()) * (100.0 + (1.0 / ener)) * (1.0 + efrac);
0160 const double w1 = 9.0 / (9.0 + d);
0161 const double umax = ener * M_PI / partm;
0162 double u;
0163
0164 do {
0165 double beta = (random->flatShoot() <= w1) ? alfa : 3.0 * alfa;
0166 u = -std::log(random->flatShoot() * random->flatShoot()) / beta;
0167 } while (u >= umax);
0168
0169 return u;
0170 }