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   double m2dontchange = particle.momentum().mass() * particle.momentum().mass();
0141
0142   // Calculate energy of these photons and add them to the event
0143   for (unsigned int i = 0; i < nPhotons; ++i) {
0144     // Throw momentum of the photon
0145     math::XYZTLorentzVector photonMom = brem(particle, xmin, random);
0146
0147     // Check that there is enough energy left.
0148     if (particle.momentum().E() - particle.momentum().mass() < photonMom.E())
0149       break;
0150
0151     // Rotate to the lab frame
0152     photonMom = ROOT::Math::RotationZ(phi) * (ROOT::Math::RotationY(theta) * photonMom);
0153
0154     // Add a photon
0155     secondaries.emplace_back(new fastsim::Particle(22, particle.position(), photonMom));
0156
0157     // Update the original e+/-
0158     particle.momentum() -= photonMom;
0159     particle.momentum().SetXYZT(particle.momentum().px(),
0160                                 particle.momentum().py(),
0161                                 particle.momentum().pz(),
0162                                 sqrt(pow(particle.momentum().P(), 2) + m2dontchange));
0163   }
0164 }
0165
0166 math::XYZTLorentzVector fastsim::Bremsstrahlung::brem(fastsim::Particle& particle,
0167                                                       double xmin,
0168                                                       const RandomEngineAndDistribution& random) const {
0169   double xp = 0;
0170   double weight = 0.;
0171
0172   do {
0173     xp = xmin * std::exp(-std::log(xmin) * random.flatShoot());
0174     weight = 1. - xp + 3. / 4. * xp * xp;
0175   } while (weight < random.flatShoot());
0176
0177   // Have photon energy. Now generate angles with respect to the z axis
0178   // defined by the incoming particle's momentum.
0179
0180   // Isotropic in phi
0181   const double phi = random.flatShoot() * 2. * M_PI;
0182   // theta from universal distribution
0183   const double theta = gbteth(particle.momentum().E(), fastsim::Constants::eMass, xp, random) *
0184                        fastsim::Constants::eMass / particle.momentum().E();
0185
0186   // Make momentum components
0187   double stheta = std::sin(theta);
0188   double ctheta = std::cos(theta);
0189   double sphi = std::sin(phi);
0190   double cphi = std::cos(phi);
0191
0192   return xp * particle.momentum().E() * math::XYZTLorentzVector(stheta * cphi, stheta * sphi, ctheta, 1.);
0193 }
0194
0195 double fastsim::Bremsstrahlung::gbteth(const double ener,
0196                                        const double partm,
0197                                        const double efrac,
0198                                        const RandomEngineAndDistribution& random) const {
0199   // Details on implementation here
0200   // http://www.dnp.fmph.uniba.sk/cernlib/asdoc/geant_html3/node299.html#GBTETH
0201   // http://svn.cern.ch/guest/AliRoot/tags/v3-07-03/GEANT321/gphys/gbteth.F
0202
0203   const double alfa = 0.625;
0204
0205   const double d = 0.13 * (0.8 + 1.3 / Z_) * (100.0 + (1.0 / ener)) * (1.0 + efrac);
0206   const double w1 = 9.0 / (9.0 + d);
0207   const double umax = ener * M_PI / partm;
0208   double u;
0209
0210   do {
0211     double beta = (random.flatShoot() <= w1) ? alfa : 3.0 * alfa;
0212     u = -std::log(random.flatShoot() * random.flatShoot()) / beta;
0213   } while (u >= umax);
0214
0215   return u;
0216 }
0217
0218 DEFINE_EDM_PLUGIN(fastsim::InteractionModelFactory, fastsim::Bremsstrahlung, "fastsim::Bremsstrahlung");