MaterialEffects/src/PairProductionSimulator.cc

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   // Set the minimal photon energy for possible conversion
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   // The photon has enough energy to create a pair
0016   if (eGamma >= photonEnergy) {
0017     // This is a simple version (a la PDG) of a photon conversion generator.
0018     // It replaces the buggy GEANT3 -> C++ former version.
0019     // Author : Patrick Janot - 7-Jan-2004
0020
0021     // Probability to convert is 7/9*(dx/X0)
0022     if (-std::log(random->flatShoot()) <= (7. / 9.) * radLengths) {
0023       double xe = 0;
0024       double xm = eMass() / eGamma;
0025       double weight = 0.;
0026
0027       // Generate electron energy between emass and eGamma-emass
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       // Generate angles
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       // The electron
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       // The positron
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 }