Back to home page

Project CMSSW displayed by LXR

 
 

    


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

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