Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:26:48

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");