MaterialEffects/src/EnergyLossSimulator.cc

0001 #include "FastSimulation/MaterialEffects/interface/EnergyLossSimulator.h"
0002 //#include "FastSimulation/Utilities/interface/RandomEngine.h"
0003 #include "FastSimulation/Utilities/interface/LandauFluctuationGenerator.h"
0004
0005 #include <cmath>
0006
0007 EnergyLossSimulator::EnergyLossSimulator(double A, double Z, double density, double radLen)
0008     : MaterialEffectsSimulator(A, Z, density, radLen) {
0009   theGenerator = new LandauFluctuationGenerator();
0010 }
0011
0012 EnergyLossSimulator::~EnergyLossSimulator() { delete theGenerator; }
0013
0014 void EnergyLossSimulator::compute(ParticlePropagator& Particle, RandomEngineAndDistribution const* random) {
0015   //  FamosHistos* myHistos = FamosHistos::instance();
0016
0017   // double gamma_e = 0.577215664901532861;  // Euler constant
0018
0019   // The thickness in cm
0020   double thick = radLengths * radLenIncm();
0021
0022   // This is a simple version (a la PDG) of a dE/dx generator.
0023   // It replaces the buggy GEANT3 -> C++ former version.
0024   // Author : Patrick Janot - 8-Jan-2004
0025
0026   double p2 = Particle.particle().Vect().Mag2();
0027   double verySmallP2 = 0.0001;
0028   if (p2 <= verySmallP2) {
0029     deltaP.SetXYZT(0., 0., 0., 0.);
0030     return;
0031   }
0032   double m2 = Particle.particle().mass() * Particle.particle().mass();
0033   double e2 = p2 + m2;
0034
0035   double beta2 = p2 / e2;
0036   double gama2 = e2 / m2;
0037
0038   double charge2 = Particle.particle().charge() * Particle.particle().charge();
0039
0040   // Energy loss spread in GeV
0041   double eSpread = 0.1536E-3 * charge2 * (theZ() / theA()) * rho() * thick / beta2;
0042
0043   // Most probable energy loss (from the integrated Bethe-Bloch equation)
0044   mostProbableLoss = eSpread * (log(2. * eMass() * beta2 * gama2 * eSpread / (excitE() * excitE())) - beta2 + 0.200);
0045
0046   // This one can be needed on output (but is not used internally)
0047   // meanEnergyLoss = 2.*eSpread * ( log ( 2.*eMass()*beta2*gama2 /excitE() ) - beta2 );
0048
0049   // Generate the energy loss with Landau fluctuations
0050   double dedx = mostProbableLoss + eSpread * theGenerator->landau(random);
0051
0052   // Compute the new energy and momentum
0053   double aBitAboveMass = Particle.particle().mass() * 1.0001;
0054   double newE = std::max(aBitAboveMass, Particle.particle().e() - dedx);
0055   //  double newE = std::max(Particle.particle().mass(),Particle.particle().e()-dedx);
0056   double fac = std::sqrt((newE * newE - m2) / p2);
0057
0058   // Update the momentum
0059   deltaP.SetXYZT(Particle.particle().Px() * (1. - fac),
0060                  Particle.particle().Py() * (1. - fac),
0061                  Particle.particle().Pz() * (1. - fac),
0062                  Particle.particle().E() - newE);
0063   Particle.particle().setMomentum(
0064       Particle.particle().Px() * fac, Particle.particle().Py() * fac, Particle.particle().Pz() * fac, newE);
0065 }