EnergyLossSimulator.cc

Go to the documentation of this file.

#include "FastSimulation/MaterialEffects/interface/EnergyLossSimulator.h"
//#include "FastSimulation/Utilities/interface/RandomEngine.h"
#include "FastSimulation/Utilities/interface/LandauFluctuationGenerator.h"

#include <cmath>

EnergyLossSimulator::EnergyLossSimulator(double A, double Z, double density, double radLen)
    : MaterialEffectsSimulator(A, Z, density, radLen) {
  theGenerator = new LandauFluctuationGenerator();
}

EnergyLossSimulator::~EnergyLossSimulator() { delete theGenerator; }

void EnergyLossSimulator::compute(ParticlePropagator& Particle, RandomEngineAndDistribution const* random) {
  //  FamosHistos* myHistos = FamosHistos::instance();

  // double gamma_e = 0.577215664901532861;  // Euler constant

  // The thickness in cm
  double thick = radLengths * radLenIncm();

  // This is a simple version (a la PDG) of a dE/dx generator.
  // It replaces the buggy GEANT3 -> C++ former version.
  // Author : Patrick Janot - 8-Jan-2004

  double p2 = Particle.particle().Vect().Mag2();
  double verySmallP2 = 0.0001;
  if (p2 <= verySmallP2) {
    deltaP.SetXYZT(0., 0., 0., 0.);
    return;
  }
  double m2 = Particle.particle().mass() * Particle.particle().mass();
  double e2 = p2 + m2;

  double beta2 = p2 / e2;
  double gama2 = e2 / m2;

  double charge2 = Particle.particle().charge() * Particle.particle().charge();

  // Energy loss spread in GeV
  double eSpread = 0.1536E-3 * charge2 * (theZ() / theA()) * rho() * thick / beta2;

  // Most probable energy loss (from the integrated Bethe-Bloch equation)
  mostProbableLoss = eSpread * (log(2. * eMass() * beta2 * gama2 * eSpread / (excitE() * excitE())) - beta2 + 0.200);

  // This one can be needed on output (but is not used internally)
  // meanEnergyLoss = 2.*eSpread * ( log ( 2.*eMass()*beta2*gama2 /excitE() ) - beta2 );

  // Generate the energy loss with Landau fluctuations
  double dedx = mostProbableLoss + eSpread * theGenerator->landau(random);

  // Compute the new energy and momentum
  double aBitAboveMass = Particle.particle().mass() * 1.0001;
  double newE = std::max(aBitAboveMass, Particle.particle().e() - dedx);
  //  double newE = std::max(Particle.particle().mass(),Particle.particle().e()-dedx);
  double fac = std::sqrt((newE * newE - m2) / p2);

  // Update the momentum
  deltaP.SetXYZT(Particle.particle().Px() * (1. - fac),
                 Particle.particle().Py() * (1. - fac),
                 Particle.particle().Pz() * (1. - fac),
                 Particle.particle().E() - newE);
  Particle.particle().setMomentum(
      Particle.particle().Px() * fac, Particle.particle().Py() * fac, Particle.particle().Pz() * fac, newE);
}