Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:11:16

0001 ///////////////////////////////////////////////
0002 //MuonBremsstrahlungSimulator
0003 // Description: Implementation of Muon bremsstrahlung using the Petrukhin model
0004 //Authors :Sandro Fonseca de Souza and Andre Sznajder (UERJ/Brazil)
0005 // Date: 23-Nov-2010
0006 //////////////////////////////////////////////////////////
0007 #include "FastSimulation/MaterialEffects/interface/MuonBremsstrahlungSimulator.h"
0008 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "FastSimulation/MaterialEffects/interface/PetrukhinModel.h"
0011 #include "FastSimulation/Particle/interface/makeParticle.h"
0012 
0013 #include <cmath>
0014 #include <string>
0015 #include <iostream>
0016 
0017 #include "TF1.h"
0018 
0019 //////////////////////////////////////////////////////////////////////////////////////
0020 MuonBremsstrahlungSimulator::MuonBremsstrahlungSimulator(
0021     double A, double Z, double density, double radLen, double photonEnergyCut, double photonFractECut)
0022     : MaterialEffectsSimulator(A, Z, density, radLen) {
0023   // Set the minimal photon energy for a Brem from mu+/-
0024   photonEnergy = photonEnergyCut;
0025   photonFractE = photonFractECut;
0026   d = 0.;  //distance
0027   LogDebug("MuonBremsstrahlungSimulator") << "Starting the MuonBremsstrahlungSimulator" << std::endl;
0028 }
0029 
0030 ///////////////////////////////////////////////////////////////////
0031 void MuonBremsstrahlungSimulator::compute(ParticlePropagator& Particle, RandomEngineAndDistribution const* random) {
0032   double NA = 6.022e+23;  //Avogadro's number
0033 
0034   if (radLengths > 4.) {
0035     Particle.particle().setMomentum(0., 0., 0., 0.);
0036     deltaPMuon.SetXYZT(0., 0., 0., 0.);
0037     brem_photon.SetXYZT(0., 0., 0., 0.);
0038   }
0039 
0040   // Hard brem probability with a photon Energy above photonEnergy.
0041 
0042   double EMuon = Particle.particle().e();  //Muon Energy
0043   if (EMuon < photonEnergy)
0044     return;
0045   xmin = std::max(photonEnergy / EMuon, photonFractE);  //fraction of muon's energy transferred to the photon
0046 
0047   //  xmax = photonEnergy/Particle.particle().e();
0048   if (xmin >= 1. || xmin <= 0.)
0049     return;
0050 
0051   xmax = 1.;
0052   npar = 3;  //Number of parameters
0053 
0054   // create TF1 using a free C function
0055   f1 = new TF1("f1", PetrukhinFunc, xmin, xmax, npar);
0056   //Setting parameters
0057   f1->SetParameters(EMuon, A, Z);
0058   /////////////////////////////////////////////////////////////////////////
0059   //d = distance for several materials
0060   //X0 = radLen
0061   //d = radLengths * X0(for tracker,yoke,ECAL and HCAL)
0062   d = radLengths * radLen;
0063   //Integration
0064   bremProba = density * d * (NA / A) * (f1->Integral(0., 1.));
0065 
0066   // Number of photons to be radiated.
0067   unsigned int nPhotons = random->poissonShoot(bremProba);
0068   _theUpdatedState.reserve(nPhotons);
0069 
0070   if (!nPhotons)
0071     return;
0072 
0073   //Rotate to the lab frame
0074   double chi = Particle.particle().theta();
0075   double psi = Particle.particle().phi();
0076   RawParticle::RotationZ rotZ(psi);
0077   RawParticle::RotationY rotY(chi);
0078 
0079   // Energy of these photons
0080   for (unsigned int i = 0; i < nPhotons; ++i) {
0081     // Check that there is enough energy left.
0082     if (Particle.particle().e() < photonEnergy)
0083       break;
0084     LogDebug("MuonBremsstrahlungSimulator") << "MuonBremsstrahlungSimulator parameters:" << std::endl;
0085     LogDebug("MuonBremsstrahlungSimulator") << "xmin-> " << xmin << std::endl;
0086     LogDebug("MuonBremsstrahlungSimulator") << "Atomic Weight-> " << A << std::endl;
0087     LogDebug("MuonBremsstrahlungSimulator") << "Density-> " << density << std::endl;
0088     LogDebug("MuonBremsstrahlungSimulator") << "Distance-> " << d << std::endl;
0089     LogDebug("MuonBremsstrahlungSimulator") << "bremProba->" << bremProba << std::endl;
0090     LogDebug("MuonBremsstrahlungSimulator") << "nPhotons->" << nPhotons << std::endl;
0091     LogDebug("MuonBremsstrahlungSimulator") << " Muon_Energy-> " << EMuon << std::endl;
0092     LogDebug("MuonBremsstrahlungSimulator") << "X0-> " << radLen << std::endl;
0093     LogDebug("MuonBremsstrahlungSimulator") << " radLengths-> " << radLengths << std::endl;
0094 
0095     // Add a photon
0096     RawParticle thePhoton = makeParticle(Particle.particleDataTable(), 22, brem(Particle, random));
0097     if (thePhoton.E() > 0.) {
0098       thePhoton.rotate(rotY);
0099       thePhoton.rotate(rotZ);
0100 
0101       _theUpdatedState.push_back(thePhoton);
0102 
0103       // Update the original mu +/-
0104       deltaPMuon = Particle.particle().momentum() -= thePhoton.momentum();
0105       // Information of brem photon
0106       brem_photon.SetXYZT(thePhoton.Px(), thePhoton.Py(), thePhoton.Pz(), thePhoton.E());
0107 
0108       LogDebug("MuonBremsstrahlungSimulator") << " Muon Bremsstrahlung: photon_energy-> " << thePhoton.E() << std::endl;
0109       LogDebug("MuonBremsstrahlungSimulator") << "photon_px->" << thePhoton.Px() << std::endl;
0110       LogDebug("MuonBremsstrahlungSimulator") << "photon_py->" << thePhoton.Py() << std::endl;
0111       LogDebug("MuonBremsstrahlungSimulator") << "photon_pz->" << thePhoton.Pz() << std::endl;
0112     }
0113   }
0114 }
0115 ///////////////////////////////////////////////////////////////////////////////////////////
0116 XYZTLorentzVector MuonBremsstrahlungSimulator::brem(ParticlePropagator& pp,
0117                                                     RandomEngineAndDistribution const* random) const {
0118   // This is a simple version of a Muon Brem using Petrukhin model .
0119   //Ref: http://pdg.lbl.gov/2008/AtomicNuclearProperties/adndt.pdf
0120   double mumass = 0.105658367;  //mu mass  (GeV/c^2)
0121 
0122   // RANDOM_NUMBER_ERROR
0123   // Random number should be generated by the engines from the
0124   // RandomNumberGeneratorService. This appears to use the global
0125   // engine in ROOT. This is not thread safe unless the module using
0126   // it is a one module and declares a shared resource and all
0127   // other modules using it also declare the same shared resource.
0128   // This also breaks replay.
0129   // It might be that this is never used because of the std::cout
0130   // statement 2 lines below. I've never seen or heard complaints
0131   // about that vebosity ....
0132   double xp = f1->GetRandom();
0133   LogDebug("MuonBremsstrahlungSimulator") << "MuonBremsstrahlungSimulator: xp->" << xp << std::endl;
0134   std::cout << "MuonBremsstrahlungSimulator: xp->" << xp << std::endl;
0135 
0136   // Have photon energy. Now generate angles with respect to the z axis
0137   // defined by the incoming particle's momentum.
0138 
0139   // Isotropic in phi
0140   const double phi = random->flatShoot() * 2 * M_PI;
0141   // theta from universal distribution
0142   const double theta = gbteth(pp.particle().e(), mumass, xp, random) * mumass / pp.particle().e();
0143 
0144   // Make momentum components
0145   double stheta = std::sin(theta);
0146   double ctheta = std::cos(theta);
0147   double sphi = std::sin(phi);
0148   double cphi = std::cos(phi);
0149 
0150   return xp * pp.particle().e() * XYZTLorentzVector(stheta * cphi, stheta * sphi, ctheta, 1.);
0151 }
0152 //////////////////////////////////////////////////////////////////////////////////////////////
0153 double MuonBremsstrahlungSimulator::gbteth(const double ener,
0154                                            const double partm,
0155                                            const double efrac,
0156                                            RandomEngineAndDistribution const* random) const {
0157   const double alfa = 0.625;
0158 
0159   const double d = 0.13 * (0.8 + 1.3 / theZ()) * (100.0 + (1.0 / ener)) * (1.0 + efrac);
0160   const double w1 = 9.0 / (9.0 + d);
0161   const double umax = ener * M_PI / partm;
0162   double u;
0163 
0164   do {
0165     double beta = (random->flatShoot() <= w1) ? alfa : 3.0 * alfa;
0166     u = -std::log(random->flatShoot() * random->flatShoot()) / beta;
0167   } while (u >= umax);
0168 
0169   return u;
0170 }