SimplifiedGeometryPropagator/plugins/MuonBremsstrahlung.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 #include <TMath.h>
0019 #include <TF1.h>
0020
0021 ///////////////////////////////////////////////
0022 // Authors: Sandro Fonseca de Souza and Andre Sznajder (UERJ/Brazil)
0023 // Date: 23-Nov-2010
0024 //
0025 // Revision: Class structure modified to match SimplifiedGeometryPropagator
0026 //           S. Sekmen, 18 May 2017
0027 //
0028 // Revision: Code very buggy, PetrukhinFunc return negative values, double bremProba wasn't properly defined etc.
0029 //           Should be all fixed by now
0030 //           S. Kurz, 23 May 2017
0031 //////////////////////////////////////////////////////////
0032
0033 namespace fastsim {
0034   //! Implementation of Bremsstrahlung from mu+/mu- in the tracker layers based on a Petrukhin Model (nuclear screening correction).
0035   /*!
0036         Computes the number, energy and angles of Bremsstrahlung photons emitted by muons
0037         and modifies mu+/mu- particle accordingly.
0038     */
0039   class MuonBremsstrahlung : public InteractionModel {
0040   public:
0041     //! Constructor.
0042     MuonBremsstrahlung(const std::string &name, const edm::ParameterSet &cfg);
0043
0044     //! Default destructor.
0045     ~MuonBremsstrahlung() override { ; };
0046
0047     //! Perform the interaction.
0048     /*!
0049             \param particle The particle that interacts with the matter.
0050             \param layer The detector layer that interacts with the particle.
0051             \param secondaries Particles that are produced in the interaction (if any).
0052             \param random The Random Engine.
0053         */
0054     void interact(Particle &particle,
0055                   const SimplifiedGeometry &layer,
0056                   std::vector<std::unique_ptr<Particle> > &secondaries,
0057                   const RandomEngineAndDistribution &random) override;
0058
0059   private:
0060     //! Compute Brem photon energy and angles, if any.
0061     /*!
0062             \param particle The particle that interacts with the matter.
0063             \param xmin Minimum fraction of the particle's energy that has to be converted to a photon.
0064             \param random The Random Engine.
0065             \return Momentum 4-vector of a bremsstrahlung photon.
0066         */
0067     math::XYZTLorentzVector brem(Particle &particle, double xmin, const RandomEngineAndDistribution &random) const;
0068
0069     //! A universal angular distribution.
0070     /*!
0071             \param ener
0072             \param partm
0073             \param efrac
0074             \param random The Random Engine.
0075             \return Theta from universal distribution
0076         */
0077     double gbteth(const double ener,
0078                   const double partm,
0079                   const double efrac,
0080                   const RandomEngineAndDistribution &random) const;
0081
0082     //! Petrukhin Function: Returns cross section using nuclear-electron screening correction from G4 style
0083     static double PetrukhinFunc(double *x, double *p);
0084
0085     TF1 *Petrfunc;                    //!< The Petrukhin Function
0086     double minPhotonEnergy_;          //!< Cut on minimum energy of bremsstrahlung photons
0087     double minPhotonEnergyFraction_;  //!< Cut on minimum fraction of particle's energy which has to be carried by photon
0088     double density_;                  //!< Density of material (usually silicon rho=2.329)
0089     double radLenInCm_;               //!< Radiation length of material (usually silicon X0=9.360)
0090     double A_;                        //!< Atomic weight of material (usually silicon A=28.0855)
0091     double Z_;                        //!< Atomic number of material (usually silicon Z=14)
0092   };
0093 }  // namespace fastsim
0094
0095 fastsim::MuonBremsstrahlung::MuonBremsstrahlung(const std::string &name, const edm::ParameterSet &cfg)
0096     : fastsim::InteractionModel(name) {
0097   // Set the minimal photon energy for a Brem from mu+/-
0098   minPhotonEnergy_ = cfg.getParameter<double>("minPhotonEnergy");
0099   minPhotonEnergyFraction_ = cfg.getParameter<double>("minPhotonEnergyFraction");
0100   // Material properties
0101   A_ = cfg.getParameter<double>("A");
0102   Z_ = cfg.getParameter<double>("Z");
0103   density_ = cfg.getParameter<double>("density");
0104   radLenInCm_ = cfg.getParameter<double>("radLen");
0105 }
0106
0107 void fastsim::MuonBremsstrahlung::interact(fastsim::Particle &particle,
0108                                            const SimplifiedGeometry &layer,
0109                                            std::vector<std::unique_ptr<fastsim::Particle> > &secondaries,
0110                                            const RandomEngineAndDistribution &random) {
0111   // only consider muons
0112   if (std::abs(particle.pdgId()) != 13) {
0113     return;
0114   }
0115
0116   double radLengths = layer.getThickness(particle.position(), particle.momentum());
0117   //
0118   // no material
0119   //
0120   if (radLengths < 1E-10) {
0121     return;
0122   }
0123
0124   // Protection : Just stop the electron if more than 1 radiation lengths.
0125   // This case corresponds to an electron entering the layer parallel to
0126   // the layer axis - no reliable simulation can be done in that case...
0127   if (radLengths > 4.) {
0128     particle.momentum().SetXYZT(0., 0., 0., 0.);
0129     return;
0130   }
0131
0132   // muon must have more energy than minimum photon energy
0133   if (particle.momentum().E() - particle.momentum().mass() < minPhotonEnergy_) {
0134     return;
0135   }
0136
0137   // Min fraction of muon's energy transferred to the photon
0138   double xmin = std::max(minPhotonEnergy_ / particle.momentum().E(), minPhotonEnergyFraction_);
0139   // Hard brem probability with a photon Energy above threshold.
0140   if (xmin >= 1. || xmin <= 0.) {
0141     return;
0142   }
0143
0144   // Max fraction of muon's energy transferred to the photon
0145   double xmax = 1.;
0146
0147   // create TF1 using a free C function
0148   Petrfunc = new TF1("Petrfunc", PetrukhinFunc, xmin, xmax, 3);
0149   // Setting parameters
0150   Petrfunc->SetParameters(particle.momentum().E(), A_, Z_);
0151   // d = distance for several materials
0152   // X0 = radLen
0153   // d = radLengths * X0 (for tracker, yoke, ECAL and HCAL)
0154   double distance = radLengths * radLenInCm_;
0155
0156   // Integration
0157   // Fixed previous version which used Petrfunc->Integral(0.,1.) -> does not make sense
0158   double bremProba = density_ * distance * (fastsim::Constants::NA / A_) * (Petrfunc->Integral(xmin, xmax));
0159   if (bremProba < 1E-10) {
0160     return;
0161   }
0162
0163   // Number of photons to be radiated.
0164   unsigned int nPhotons = random.poissonShoot(bremProba);
0165   if (nPhotons == 0) {
0166     return;
0167   }
0168
0169   //Rotate to the lab frame
0170   double theta = particle.momentum().Theta();
0171   double phi = particle.momentum().Phi();
0172
0173   // Energy of these photons
0174   for (unsigned int i = 0; i < nPhotons; ++i) {
0175     // Throw momentum of the photon
0176     math::XYZTLorentzVector photonMom = brem(particle, xmin, random);
0177
0178     // Check that there is enough energy left.
0179     if (particle.momentum().E() - particle.momentum().mass() < photonMom.E())
0180       break;
0181
0182     // Rotate to the lab frame
0183     photonMom = ROOT::Math::RotationZ(phi) * (ROOT::Math::RotationY(theta) * photonMom);
0184
0185     // Add a photon
0186     secondaries.emplace_back(new fastsim::Particle(22, particle.position(), photonMom));
0187
0188     // Update the original e+/-
0189     particle.momentum() -= photonMom;
0190   }
0191 }
0192
0193 math::XYZTLorentzVector fastsim::MuonBremsstrahlung::brem(fastsim::Particle &particle,
0194                                                           double xmin,
0195                                                           const RandomEngineAndDistribution &random) const {
0196   // This is a simple version of a Muon Brem using Petrukhin model.
0197   // Ref: http://pdg.lbl.gov/2008/AtomicNuclearProperties/adndt.pdf
0198   double xp = Petrfunc->GetRandom();
0199
0200   // Have photon energy. Now generate angles with respect to the z axis
0201   // defined by the incoming particle's momentum.
0202
0203   // Isotropic in phi
0204   const double phi = random.flatShoot() * 2. * M_PI;
0205   // theta from universal distribution
0206   const double theta = gbteth(particle.momentum().E(), fastsim::Constants::muMass, xp, random) *
0207                        fastsim::Constants::muMass / particle.momentum().E();
0208
0209   // Make momentum components
0210   double stheta = std::sin(theta);
0211   double ctheta = std::cos(theta);
0212   double sphi = std::sin(phi);
0213   double cphi = std::cos(phi);
0214
0215   return xp * particle.momentum().E() * math::XYZTLorentzVector(stheta * cphi, stheta * sphi, ctheta, 1.);
0216 }
0217
0218 double fastsim::MuonBremsstrahlung::gbteth(const double ener,
0219                                            const double partm,
0220                                            const double efrac,
0221                                            const RandomEngineAndDistribution &random) const {
0222   // Details on implementation here
0223   // http://www.dnp.fmph.uniba.sk/cernlib/asdoc/geant_html3/node299.html#GBTETH
0224   // http://svn.cern.ch/guest/AliRoot/tags/v3-07-03/GEANT321/gphys/gbteth.F
0225
0226   const double alfa = 0.625;
0227
0228   const double d = 0.13 * (0.8 + 1.3 / Z_) * (100.0 + (1.0 / ener)) * (1.0 + efrac);
0229   const double w1 = 9.0 / (9.0 + d);
0230   const double umax = ener * M_PI / partm;
0231   double u;
0232
0233   do {
0234     double beta = (random.flatShoot() <= w1) ? alfa : 3.0 * alfa;
0235     u = -std::log(random.flatShoot() * random.flatShoot()) / beta;
0236   } while (u >= umax);
0237
0238   return u;
0239 }
0240
0241 double fastsim::MuonBremsstrahlung::PetrukhinFunc(double *x, double *p) {
0242   // Function independent variable
0243   double nu = x[0];  // fraction of muon's energy transferred to the photon
0244
0245   // Parameters
0246   double E = p[0];  // Muon Energy (in GeV)
0247   double A = p[1];  // Atomic weight
0248   double Z = p[2];  // Atomic number
0249
0250   /*
0251         //////////////////////////////////////////////////
0252         // Function of Muom Brem using  nuclear screening correction
0253         // Ref: http://pdg.lbl.gov/2008/AtomicNuclearProperties/adndt.pdf
0254         // http://geant4.cern.ch/G4UsersDocuments/UsersGuides/PhysicsReferenceManual/html/node48.html
0255
0256         // Physical constants
0257         double B = 182.7;
0258         double ee = sqrt(2.7181) ; // sqrt(e)
0259         double ZZ=  pow( Z,-1./3.); // Z^-1/3
0260         ///////////////////////////////////////////////////////////////////
0261         double emass = 0.0005109990615;  // electron mass (GeV/c^2)
0262         double mumass = 0.105658367;//mu mass  (GeV/c^2)
0263
0264         double re = 2.817940285e-13;// Classical electron radius (Units: cm)
0265         double alpha = 1./137.03599976; // fine structure constant
0266         double Dn = 1.54* (pow(A,0.27));
0267         double constant =  pow((2.0 * Z * emass/mumass * re ),2.0);
0268         //////////////////////////////////////////////
0269
0270         double delta = (mumass * mumass * nu) /(2.* E * (1.- nu));
0271
0272         double Delta_n = TMath::Log(Dn / (1.+ delta *( Dn * ee -2.)/ mumass)); //nuclear screening correction
0273
0274         double Phi = TMath::Log((B * mumass * ZZ / emass)/ (1.+ delta * ee * B * ZZ  / emass)) - Delta_n;//phi(delta)
0275
0276         // Diff. Cross Section for Muon Brem from a screened nuclear (Equation 16: REF: LBNL-44742)
0277         double f = alpha * constant *(4./3.-4./3.*nu + nu*nu)*Phi/nu;
0278     */
0279
0280   //////////////////////////////////////////////////
0281   // Function for Muon Brem Xsec from G4
0282   //////////////////////////////////////////////////
0283
0284   // Physical constants
0285   double B = 183.;
0286   double Bl = 1429.;
0287   double ee = 1.64872;            // sqrt(e)
0288   double Z13 = pow(Z, -1. / 3.);  // Z^-1/3
0289   double Z23 = pow(Z, -2. / 3.);  // Z^-2/3
0290
0291   // Original values of paper
0292   double emass = 0.0005109990615;  // electron mass (GeV/c^2)
0293   double mumass = 0.105658367;     // muon mass  (GeV/c^2)
0294   // double re = 2.817940285e-13;     // Classical electron radius (Units: cm)
0295   double alpha = 0.00729735;      // 1./137.03599976; // fine structure constant
0296   double constant = 1.85736e-30;  // pow( ( emass / mumass * re ) , 2.0);
0297
0298   // Use nomenclature from reference -> Follow those formula step by step
0299   if (nu * E >= E - mumass)
0300     return 0;
0301
0302   double Dn = 1.54 * (pow(A, 0.27));
0303   double Dnl = pow(Dn, (1. - 1. / Z));
0304
0305   double delta = (mumass * mumass * nu) / (2. * E * (1. - nu));
0306
0307   double Phi_n = TMath::Log(B * Z13 * (mumass + delta * (Dnl * ee - 2)) / (Dnl * (emass + delta * ee * B * Z13)));
0308   if (Phi_n < 0)
0309     Phi_n = 0;
0310
0311   double Phi_e =
0312       TMath::Log((Bl * Z23 * mumass) / (1. + delta * mumass / (emass * emass * ee)) / (emass + delta * ee * Bl * Z23));
0313   if (Phi_e < 0 || nu * E >= E / (1. + (mumass * mumass / (2. * emass * E))))
0314     Phi_e = 0;
0315
0316   // Diff. Cross Section for Muon Brem from G4 (without NA/A factor)
0317   double f = 16. / 3. * alpha * constant * Z * (Z * Phi_n + Phi_e) * (1. / nu) * (1. - nu + 0.75 * nu * nu);
0318
0319   return f;
0320 }
0321
0322 DEFINE_EDM_PLUGIN(fastsim::InteractionModelFactory, fastsim::MuonBremsstrahlung, "fastsim::MuonBremsstrahlung");