SimplifiedGeometryPropagator/plugins/Bremsstrahlung.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: 25-Dec-2003
0021 //
0022 // Revision: Class structure modified to match SimplifiedGeometryPropagator
0023 //           Fixed a bug in which particles could radiate more energy than they have
0024 //           S. Kurz, 29 May 2017
0025 //////////////////////////////////////////////////////////
0026
0027 namespace fastsim {
0028   //! Implementation of Bremsstrahlung from e+/e- in the tracker layers.
0029   /*!
0030         Computes the number, energy and angles of Bremsstrahlung photons emitted by electrons and positrons
0031         and modifies e+/e- particle accordingly.
0032     */
0033   class Bremsstrahlung : public InteractionModel {
0034   public:
0035     //! Constructor.
0036     Bremsstrahlung(const std::string& name, const edm::ParameterSet& cfg);
0037
0038     //! Default destructor.
0039     ~Bremsstrahlung() override { ; };
0040
0041     //! Perform the interaction.
0042     /*!
0043             \param particle The particle that interacts with the matter.
0044             \param layer The detector layer that interacts with the particle.
0045             \param secondaries Particles that are produced in the interaction (if any).
0046             \param random The Random Engine.
0047         */
0048     void interact(Particle& particle,
0049                   const SimplifiedGeometry& layer,
0050                   std::vector<std::unique_ptr<Particle> >& secondaries,
0051                   const RandomEngineAndDistribution& random) override;
0052
0053   private:
0054     //! Compute Brem photon energy and angles, if any.
0055     /*!
0056             \param particle The particle that interacts with the matter.
0057             \param xmin Minimum fraction of the particle's energy that has to be converted to a photon.
0058             \param random The Random Engine.
0059             \return Momentum 4-vector of a bremsstrahlung photon.
0060         */
0061     math::XYZTLorentzVector brem(Particle& particle, double xmin, const RandomEngineAndDistribution& random) const;
0062
0063     //! A universal angular distribution.
0064     /*!
0065             \param ener
0066             \param partm
0067             \param efrac
0068             \param random The Random Engine.
0069             \return Theta from universal distribution
0070         */
0071     double gbteth(const double ener,
0072                   const double partm,
0073                   const double efrac,
0074                   const RandomEngineAndDistribution& random) const;
0075
0076     double minPhotonEnergy_;          //!< Cut on minimum energy of bremsstrahlung photons
0077     double minPhotonEnergyFraction_;  //!< Cut on minimum fraction of particle's energy which has to be carried by photon
0078     double Z_;                        //!< Atomic number of material (usually silicon Z=14)
0079   };
0080 }  // namespace fastsim
0081
0082 fastsim::Bremsstrahlung::Bremsstrahlung(const std::string& name, const edm::ParameterSet& cfg)
0083     : fastsim::InteractionModel(name) {
0084   // Set the minimal photon energy for a Brem from e+/-
0085   minPhotonEnergy_ = cfg.getParameter<double>("minPhotonEnergy");
0086   minPhotonEnergyFraction_ = cfg.getParameter<double>("minPhotonEnergyFraction");
0087   // Material properties
0088   Z_ = cfg.getParameter<double>("Z");
0089 }
0090
0091 void fastsim::Bremsstrahlung::interact(fastsim::Particle& particle,
0092                                        const SimplifiedGeometry& layer,
0093                                        std::vector<std::unique_ptr<fastsim::Particle> >& secondaries,
0094                                        const RandomEngineAndDistribution& random) {
0095   // only consider electrons and positrons
0096   if (std::abs(particle.pdgId()) != 11) {
0097     return;
0098   }
0099
0100   double radLengths = layer.getThickness(particle.position(), particle.momentum());
0101   //
0102   // no material
0103   //
0104   if (radLengths < 1E-10) {
0105     return;
0106   }
0107
0108   // Protection : Just stop the electron if more than 1 radiation lengths.
0109   // This case corresponds to an electron entering the layer parallel to
0110   // the layer axis - no reliable simulation can be done in that case...
0111   if (radLengths > 4.) {
0112     particle.momentum().SetXYZT(0., 0., 0., 0.);
0113     return;
0114   }
0115
0116   // electron must have more energy than minimum photon energy
0117   if (particle.momentum().E() - particle.momentum().mass() < minPhotonEnergy_) {
0118     return;
0119   }
0120
0121   // Hard brem probability with a photon Energy above threshold.
0122   double xmin = std::max(minPhotonEnergy_ / particle.momentum().E(), minPhotonEnergyFraction_);
0123   if (xmin >= 1. || xmin <= 0.) {
0124     return;
0125   }
0126
0127   // probability to radiate a photon
0128   double bremProba =
0129       radLengths * (4. / 3. * std::log(1. / xmin) - 4. / 3. * (1. - xmin) + 1. / 2. * (1. - xmin * xmin));
0130
0131   // Number of photons to be radiated.
0132   unsigned int nPhotons = random.poissonShoot(bremProba);
0133   if (nPhotons == 0) {
0134     return;
0135   }
0136
0137   // Needed to rotate photons to the lab frame
0138   double theta = particle.momentum().Theta();
0139   double phi = particle.momentum().Phi();
0140
0141   // Calculate energy of these photons and add them to the event
0142   for (unsigned int i = 0; i < nPhotons; ++i) {
0143     // Throw momentum of the photon
0144     math::XYZTLorentzVector photonMom = brem(particle, xmin, random);
0145
0146     // Check that there is enough energy left.
0147     if (particle.momentum().E() - particle.momentum().mass() < photonMom.E())
0148       break;
0149
0150     // Rotate to the lab frame
0151     photonMom = ROOT::Math::RotationZ(phi) * (ROOT::Math::RotationY(theta) * photonMom);
0152
0153     // Add a photon
0154     secondaries.emplace_back(new fastsim::Particle(22, particle.position(), photonMom));
0155
0156     // Update the original e+/-
0157     particle.momentum() -= photonMom;
0158   }
0159 }
0160
0161 math::XYZTLorentzVector fastsim::Bremsstrahlung::brem(fastsim::Particle& particle,
0162                                                       double xmin,
0163                                                       const RandomEngineAndDistribution& random) const {
0164   double xp = 0;
0165   double weight = 0.;
0166
0167   do {
0168     xp = xmin * std::exp(-std::log(xmin) * random.flatShoot());
0169     weight = 1. - xp + 3. / 4. * xp * xp;
0170   } while (weight < random.flatShoot());
0171
0172   // Have photon energy. Now generate angles with respect to the z axis
0173   // defined by the incoming particle's momentum.
0174
0175   // Isotropic in phi
0176   const double phi = random.flatShoot() * 2. * M_PI;
0177   // theta from universal distribution
0178   const double theta = gbteth(particle.momentum().E(), fastsim::Constants::eMass, xp, random) *
0179                        fastsim::Constants::eMass / particle.momentum().E();
0180
0181   // Make momentum components
0182   double stheta = std::sin(theta);
0183   double ctheta = std::cos(theta);
0184   double sphi = std::sin(phi);
0185   double cphi = std::cos(phi);
0186
0187   return xp * particle.momentum().E() * math::XYZTLorentzVector(stheta * cphi, stheta * sphi, ctheta, 1.);
0188 }
0189
0190 double fastsim::Bremsstrahlung::gbteth(const double ener,
0191                                        const double partm,
0192                                        const double efrac,
0193                                        const RandomEngineAndDistribution& random) const {
0194   // Details on implementation here
0195   // http://www.dnp.fmph.uniba.sk/cernlib/asdoc/geant_html3/node299.html#GBTETH
0196   // http://svn.cern.ch/guest/AliRoot/tags/v3-07-03/GEANT321/gphys/gbteth.F
0197
0198   const double alfa = 0.625;
0199
0200   const double d = 0.13 * (0.8 + 1.3 / Z_) * (100.0 + (1.0 / ener)) * (1.0 + efrac);
0201   const double w1 = 9.0 / (9.0 + d);
0202   const double umax = ener * M_PI / partm;
0203   double u;
0204
0205   do {
0206     double beta = (random.flatShoot() <= w1) ? alfa : 3.0 * alfa;
0207     u = -std::log(random.flatShoot() * random.flatShoot()) / beta;
0208   } while (u >= umax);
0209
0210   return u;
0211 }
0212
0213 DEFINE_EDM_PLUGIN(fastsim::InteractionModelFactory, fastsim::Bremsstrahlung, "fastsim::Bremsstrahlung");