SimplifiedGeometryPropagator/plugins/PairProduction.cc

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 // Author: Patrick Janot
0020 // Date: 24-Dec-2003
0021 //
0022 // Revision: Class structure modified to match SimplifiedGeometryPropagator
0023 //           S. Kurz, 29 May 2017
0024 //////////////////////////////////////////////////////////
0025
0026 namespace fastsim {
0027   //! Computes the probability for photons to convert into an e+e- pair in the tracker layer.
0028   /*!
0029         In case, it returns a list of two Secondaries (e+ and e-).
0030     */
0031   class PairProduction : public InteractionModel {
0032   public:
0033     //! Constructor.
0034     PairProduction(const std::string& name, const edm::ParameterSet& cfg);
0035
0036     //! Default destructor.
0037     ~PairProduction() override { ; };
0038
0039     //! Perform the interaction.
0040     /*!
0041             \param particle The particle that interacts with the matter.
0042             \param layer The detector layer that interacts with the particle.
0043             \param secondaries Particles that are produced in the interaction (if any).
0044             \param random The Random Engine.
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     //! A universal angular distribution.
0053     /*!
0054             \param ener
0055             \param partm
0056             \param efrac
0057             \param random The Random Engine.
0058             \return Theta from universal distribution
0059         */
0060     double gbteth(double ener, double partm, double efrac, const RandomEngineAndDistribution& random) const;
0061
0062     double minPhotonEnergy_;  //!< Cut on minimum energy of photons
0063     double Z_;                //!< Atomic number of material (usually silicon Z=14)
0064   };
0065 }  // namespace fastsim
0066
0067 fastsim::PairProduction::PairProduction(const std::string& name, const edm::ParameterSet& cfg)
0068     : fastsim::InteractionModel(name) {
0069   // Set the minimal photon energy for possible conversion
0070   minPhotonEnergy_ = cfg.getParameter<double>("photonEnergyCut");
0071   // Material properties
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   // only consider photons
0082   //
0083   if (particle.pdgId() != 22) {
0084     return;
0085   }
0086
0087   double radLengths = layer.getThickness(particle.position(), particle.momentum());
0088   //
0089   // no material
0090   //
0091   if (radLengths < 1E-10) {
0092     return;
0093   }
0094
0095   //
0096   // The photon has enough energy to create a pair
0097   //
0098   if (eGamma < minPhotonEnergy_) {
0099     return;
0100   }
0101
0102   //
0103   // Probability to convert is 7/9*(dx/X0)
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   // Generate electron energy between emass and eGamma-emass
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   // the electron
0121   double eElectron = xe * eGamma;
0122   double tElectron = eElectron - eMass;
0123   double pElectron = std::sqrt(std::max((eElectron + eMass) * tElectron, 0.));
0124
0125   // the positron
0126   double ePositron = eGamma - eElectron;
0127   double tPositron = ePositron - eMass;
0128   double pPositron = std::sqrt((ePositron + eMass) * tPositron);
0129
0130   // Generate angles
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   //Rotate to the lab frame
0152   double thetaLab = particle.momentum().Theta();
0153   double phiLab = particle.momentum().Phi();
0154
0155   // Add a electron
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   // Add a positron
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   // The photon converted
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   // Details on implementation here
0181   // http://www.dnp.fmph.uniba.sk/cernlib/asdoc/geant_html3/node299.html#GBTETH
0182   // http://svn.cern.ch/guest/AliRoot/tags/v3-07-03/GEANT321/gphys/gbteth.F
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");