File indexing completed on 2024-04-06 12:11:20
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
0017
0018
0019
0020
0021
0022
0023
0024 namespace fastsim {
0025
0026
0027
0028
0029
0030
0031
0032 class EnergyLoss : public InteractionModel {
0033 public:
0034
0035 EnergyLoss(const std::string& name, const edm::ParameterSet& cfg);
0036
0037
0038 ~EnergyLoss() override { ; };
0039
0040
0041
0042
0043
0044
0045
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;
0054 double minMomentum_;
0055 double density_;
0056 double radLenInCm_;
0057 double A_;
0058 double Z_;
0059 };
0060 }
0061
0062 fastsim::EnergyLoss::EnergyLoss(const std::string& name, const edm::ParameterSet& cfg)
0063 : fastsim::InteractionModel(name), theGenerator(LandauFluctuationGenerator()) {
0064
0065 minMomentum_ = cfg.getParameter<double>("minMomentumCut");
0066
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
0078 particle.setEnergyDeposit(0);
0079
0080
0081
0082
0083 double radLengths = layer.getThickness(particle.position(), particle.momentum());
0084 if (radLengths < 1E-10) {
0085 return;
0086 }
0087
0088
0089
0090
0091 if (particle.charge() == 0) {
0092 return;
0093 }
0094
0095
0096
0097
0098 double p2 = particle.momentum().Vect().Mag2();
0099 if (p2 < minMomentum_ * minMomentum_) {
0100 return;
0101 }
0102
0103
0104 double excitE = 12.5E-9 * Z_;
0105
0106
0107 double thick = radLengths * radLenInCm_;
0108
0109
0110
0111
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
0122 double eSpread = 0.1536E-3 * charge2 * (Z_ / A_) * density_ * thick / beta2;
0123
0124
0125 double mostProbableLoss =
0126 eSpread * (log(2. * fastsim::Constants::eMass * beta2 * gama2 * eSpread / (excitE * excitE)) - beta2 + 0.200);
0127
0128
0129 double dedx = mostProbableLoss + eSpread * theGenerator.landau(&random);
0130
0131
0132 double newE = particle.momentum().e() - dedx;
0133
0134
0135 double eDiff2 = newE * newE - m2;
0136 if (eDiff2 < 0) {
0137 particle.momentum().SetXYZT(0., 0., 0., 0.);
0138
0139
0140 particle.setEnergyDeposit(particle.momentum().e() - particle.momentum().mass());
0141 return;
0142 }
0143
0144
0145 double fac = std::sqrt(eDiff2 / p2);
0146
0147
0148
0149 particle.setEnergyDeposit(dedx);
0150
0151
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");