MaterialEffects/src/BremsstrahlungSimulator.cc

0001 //FAMOS Headers
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   // Set the minimal photon energy for a Brem from e+/-
0010   photonEnergy = photonEnergyCut;
0011   photonFractE = photonFractECut;
0012 }
0013
0014 void BremsstrahlungSimulator::compute(ParticlePropagator& Particle, RandomEngineAndDistribution const* random) {
0015   // Protection : Just stop the electron if more than 1 radiation lengths.
0016   // This case corresponds to an electron entering the layer parallel to
0017   // the layer axis - no reliable simulation can be done in that case...
0018   // 08/02/06 - pv: increase protection from 1 to 4 X0 for eta>4.8 region
0019   // if ( radLengths > 1. ) Particle.particle().setMomentum(0.,0.,0.,0.);
0020   if (radLengths > 4.)
0021     Particle.particle().setMomentum(0., 0., 0., 0.);
0022
0023   // Hard brem probability with a photon Energy above photonEnergy.
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   // Number of photons to be radiated.
0034   unsigned int nPhotons = poisson(bremProba, random);
0035   _theUpdatedState.reserve(nPhotons);
0036
0037   if (!nPhotons)
0038     return;
0039
0040   //Rotate to the lab frame
0041   double chi = Particle.particle().theta();
0042   double psi = Particle.particle().phi();
0043   RawParticle::RotationZ rotZ(psi);
0044   RawParticle::RotationY rotY(chi);
0045
0046   // Energy of these photons
0047   for (unsigned int i = 0; i < nPhotons; ++i) {
0048     // Check that there is enough energy left.
0049     if (Particle.particle().e() < photonEnergy)
0050       break;
0051
0052     // Add a photon
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     // Update the original e+/-
0059     Particle.particle().momentum() -= thePhoton.momentum();
0060   }
0061 }
0062
0063 XYZTLorentzVector BremsstrahlungSimulator::brem(ParticlePropagator& pp,
0064                                                 RandomEngineAndDistribution const* random) const {
0065   // This is a simple version (a la PDG) of a Brem generator.
0066   // It replaces the buggy GEANT3 -> C++ former version.
0067   // Author : Patrick Janot - 25-Dec-2003
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   // Have photon energy. Now generate angles with respect to the z axis
0078   // defined by the incoming particle's momentum.
0079
0080   // Isotropic in phi
0081   const double phi = random->flatShoot() * 2 * M_PI;
0082   // theta from universal distribution
0083   const double theta = gbteth(pp.particle().e(), emass, xp, random) * emass / pp.particle().e();
0084
0085   // Make momentum components
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 }