SimplifiedGeometryPropagator/plugins/EnergyLoss.cc

0001 #include "FastSimulation/Utilities/interface/LandauFluctuationGenerator.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 "FastSimulation/SimplifiedGeometryPropagator/interface/InteractionModelFactory.h"
0011 #include "FastSimulation/SimplifiedGeometryPropagator/interface/InteractionModel.h"
0012 #include "FastSimulation/SimplifiedGeometryPropagator/interface/Constants.h"
0013 #include "DataFormats/Math/interface/LorentzVector.h"
0014
0015 ///////////////////////////////////////////////
0016 // Author: Patrick Janot
0017 // Date: 8-Jan-2004
0018 //
0019 // Revision: Class structure modified to match SimplifiedGeometryPropagator
0020 //           Fixed a bug in which particles could deposit more energy than they have
0021 //           S. Kurz, 29 May 2017
0022 //////////////////////////////////////////////////////////
0023
0024 namespace fastsim {
0025   //! Implementation of most probable energy loss by ionization in the tracker layers.
0026   /*!
0027         Computes the most probable energy loss by ionization from a charged particle in the tracker layer,
0028         smears it with Landau fluctuations and returns the particle with modified energy.
0029         The deposited energy is assigned with a produced SimHit (if active material hit).
0030         \sa TrackerSimHitProducer
0031     */
0032   class EnergyLoss : public InteractionModel {
0033   public:
0034     //! Constructor.
0035     EnergyLoss(const std::string& name, const edm::ParameterSet& cfg);
0036
0037     //! Default destructor.
0038     ~EnergyLoss() override { ; };
0039
0040     //! Perform the interaction.
0041     /*!
0042             \param particle The particle that interacts with the matter.
0043             \param layer The detector layer that interacts with the particle.
0044             \param secondaries Particles that are produced in the interaction (if any).
0045             \param random The Random Engine.
0046         */
0047     void interact(fastsim::Particle& particle,
0048                   const SimplifiedGeometry& layer,
0049                   std::vector<std::unique_ptr<fastsim::Particle> >& secondaries,
0050                   const RandomEngineAndDistribution& random) override;
0051
0052   private:
0053     LandauFluctuationGenerator theGenerator;  //!< Generator to do Landau fluctuation
0054     double minMomentum_;                      //!< Minimum momentum of incoming (charged) particle
0055     double density_;                          //!< Density of material (usually silicon rho=2.329)
0056     double radLenInCm_;                       //!< Radiation length of material (usually silicon X0=9.360)
0057     double A_;                                //!< Atomic weight of material (usually silicon A=28.0855)
0058     double Z_;                                //!< Atomic number of material (usually silicon Z=14)
0059   };
0060 }  // namespace fastsim
0061
0062 fastsim::EnergyLoss::EnergyLoss(const std::string& name, const edm::ParameterSet& cfg)
0063     : fastsim::InteractionModel(name), theGenerator(LandauFluctuationGenerator()) {
0064   // Set the minimal momentum
0065   minMomentum_ = cfg.getParameter<double>("minMomentumCut");
0066   // Material properties
0067   A_ = cfg.getParameter<double>("A");
0068   Z_ = cfg.getParameter<double>("Z");
0069   density_ = cfg.getParameter<double>("density");
0070   radLenInCm_ = cfg.getParameter<double>("radLen");
0071 }
0072
0073 void fastsim::EnergyLoss::interact(fastsim::Particle& particle,
0074                                    const SimplifiedGeometry& layer,
0075                                    std::vector<std::unique_ptr<fastsim::Particle> >& secondaries,
0076                                    const RandomEngineAndDistribution& random) {
0077   // Reset the energy deposit in the layer
0078   particle.setEnergyDeposit(0);
0079
0080   //
0081   // no material
0082   //
0083   double radLengths = layer.getThickness(particle.position(), particle.momentum());
0084   if (radLengths < 1E-10) {
0085     return;
0086   }
0087
0088   //
0089   // only charged particles
0090   //
0091   if (particle.charge() == 0) {
0092     return;
0093   }
0094
0095   //
0096   // minimum momentum
0097   //
0098   double p2 = particle.momentum().Vect().Mag2();
0099   if (p2 < minMomentum_ * minMomentum_) {
0100     return;
0101   }
0102
0103   // Mean excitation energy (in GeV)
0104   double excitE = 12.5E-9 * Z_;
0105
0106   // The thickness in cm
0107   double thick = radLengths * radLenInCm_;
0108
0109   // This is a simple version (a la PDG) of a dE/dx generator.
0110   // It replaces the buggy GEANT3 -> C++ former version.
0111   // Author : Patrick Janot - 8-Jan-2004
0112
0113   double m2 = particle.momentum().mass() * particle.momentum().mass();
0114   double e2 = p2 + m2;
0115
0116   double beta2 = p2 / e2;
0117   double gama2 = e2 / m2;
0118
0119   double charge2 = particle.charge() * particle.charge();
0120
0121   // Energy loss spread in GeV
0122   double eSpread = 0.1536E-3 * charge2 * (Z_ / A_) * density_ * thick / beta2;
0123
0124   // Most probable energy loss (from the integrated Bethe-Bloch equation)
0125   double mostProbableLoss =
0126       eSpread * (log(2. * fastsim::Constants::eMass * beta2 * gama2 * eSpread / (excitE * excitE)) - beta2 + 0.200);
0127
0128   // Generate the energy loss with Landau fluctuations
0129   double dedx = mostProbableLoss + eSpread * theGenerator.landau(&random);
0130
0131   // Compute the new energy and momentum
0132   double newE = particle.momentum().e() - dedx;
0133
0134   // Particle is stopped
0135   double eDiff2 = newE * newE - m2;
0136   if (eDiff2 < 0) {
0137     particle.momentum().SetXYZT(0., 0., 0., 0.);
0138     // The energy is deposited in the detector
0139     // Assigned with SimHit (if active layer) -> see TrackerSimHitProducer
0140     particle.setEnergyDeposit(particle.momentum().e() - particle.momentum().mass());
0141     return;
0142   }
0143
0144   // Relative change in momentum
0145   double fac = std::sqrt(eDiff2 / p2);
0146
0147   // The energy is deposited in the detector
0148   // Assigned with SimHit (if active layer) -> see TrackerSimHitProducer
0149   particle.setEnergyDeposit(dedx);
0150
0151   // Update the momentum
0152   particle.momentum().SetXYZT(
0153       particle.momentum().Px() * fac, particle.momentum().Py() * fac, particle.momentum().Pz() * fac, newE);
0154 }
0155
0156 DEFINE_EDM_PLUGIN(fastsim::InteractionModelFactory, fastsim::EnergyLoss, "fastsim::EnergyLoss");