File indexing completed on 2024-04-06 12:11:20
0001 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0002 #include "FastSimulation/SimplifiedGeometryPropagator/interface/Particle.h"
0003 #include "FastSimulation/SimplifiedGeometryPropagator/interface/SimplifiedGeometry.h"
0004 #include "FastSimulation/SimplifiedGeometryPropagator/interface/Constants.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007
0008 #include <cmath>
0009 #include <memory>
0010
0011 #include <Math/AxisAngle.h>
0012
0013 #include "FastSimulation/SimplifiedGeometryPropagator/interface/InteractionModelFactory.h"
0014 #include "FastSimulation/SimplifiedGeometryPropagator/interface/InteractionModel.h"
0015 #include "DataFormats/Math/interface/LorentzVector.h"
0016 #include "DataFormats/Math/interface/Vector3D.h"
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027 typedef math::XYZVector XYZVector;
0028
0029 namespace fastsim {
0030
0031
0032
0033
0034 class MultipleScattering : public InteractionModel {
0035 public:
0036
0037 MultipleScattering(const std::string& name, const edm::ParameterSet& cfg);
0038
0039
0040 ~MultipleScattering() override { ; };
0041
0042
0043
0044
0045
0046
0047
0048
0049 void interact(Particle& particle,
0050 const SimplifiedGeometry& layer,
0051 std::vector<std::unique_ptr<Particle> >& secondaries,
0052 const RandomEngineAndDistribution& random) override;
0053
0054 private:
0055
0056 XYZVector orthogonal(const XYZVector& aVector) const;
0057
0058 double minPt_;
0059 double radLenInCm_;
0060 };
0061 }
0062
0063 fastsim::MultipleScattering::MultipleScattering(const std::string& name, const edm::ParameterSet& cfg)
0064 : fastsim::InteractionModel(name) {
0065
0066 minPt_ = cfg.getParameter<double>("minPt");
0067
0068 radLenInCm_ = cfg.getParameter<double>("radLen");
0069 }
0070
0071 void fastsim::MultipleScattering::interact(fastsim::Particle& particle,
0072 const SimplifiedGeometry& layer,
0073 std::vector<std::unique_ptr<fastsim::Particle> >& secondaries,
0074 const RandomEngineAndDistribution& random) {
0075
0076
0077
0078 if (particle.charge() == 0) {
0079 return;
0080 }
0081
0082 double radLengths = layer.getThickness(particle.position(), particle.momentum());
0083
0084
0085
0086 if (radLengths < 1E-10) {
0087 return;
0088 }
0089
0090
0091 if (particle.momentum().Pt() < minPt_) {
0092 return;
0093 }
0094
0095 double p2 = particle.momentum().Vect().Mag2();
0096 double m2 = particle.momentum().mass2();
0097 double e = std::sqrt(p2 + m2);
0098
0099
0100 double pbeta = p2 / e;
0101
0102
0103
0104 double theta0 = 0.0136 / pbeta * particle.charge() * std::sqrt(2. * radLengths) * (1. + 0.038 * std::log(radLengths));
0105
0106
0107 double theta = random.gaussShoot(0., theta0);
0108
0109 double phi = 2. * M_PI * random.flatShoot();
0110
0111
0112 ROOT::Math::AxisAngle rotation1(orthogonal(particle.momentum().Vect()), theta);
0113 ROOT::Math::AxisAngle rotation2(particle.momentum().Vect(), phi);
0114
0115 XYZVector rotated = rotation2((rotation1(particle.momentum().Vect())));
0116 particle.momentum().SetXYZT(rotated.X(), rotated.Y(), rotated.Z(), particle.momentum().E());
0117
0118
0119
0120
0121 double xp = (cos(phi) * theta / 2. + random.gaussShoot(0., theta0) / sqrt(12.)) * radLengths * radLenInCm_;
0122 double yp = (sin(phi) * theta / 2. + random.gaussShoot(0., theta0) / sqrt(12.)) * radLengths * radLenInCm_;
0123
0124
0125 XYZVector normal;
0126 if (layer.isForward()) {
0127 normal = XYZVector(0., 0., particle.momentum().Pz() > 0. ? 1. : -1.);
0128 } else {
0129 double norm = particle.position().Rho();
0130 normal = XYZVector(particle.position().X() / norm, particle.position().Y() / norm, 0.);
0131 }
0132 XYZVector tangent = orthogonal(normal);
0133
0134 XYZVector delta = xp * tangent + yp * normal.Cross(tangent);
0135
0136 particle.position().SetXYZT(particle.position().X() + delta.X(),
0137 particle.position().Y() + delta.Y(),
0138 particle.position().Z() + delta.Z(),
0139 particle.position().T());
0140
0141
0142
0143 if (!layer.isForward() && std::abs(layer.getGeomProperty() - particle.position().Rho()) > 1e-2) {
0144
0145 double scalePos = layer.getGeomProperty() / particle.position().Rho();
0146 particle.position().SetXYZT(particle.position().X() * scalePos,
0147 particle.position().Y() * scalePos,
0148 particle.position().Z(),
0149 particle.position().T());
0150
0151
0152 if (std::abs(layer.getGeomProperty() - particle.position().Rho()) > 1e-2) {
0153 throw cms::Exception("fastsim::MultipleScattering") << "particle no longer on layer's surface";
0154 }
0155 }
0156 }
0157
0158 XYZVector fastsim::MultipleScattering::orthogonal(const XYZVector& aVector) const {
0159 double x = fabs(aVector.X());
0160 double y = fabs(aVector.Y());
0161 double z = fabs(aVector.Z());
0162
0163 if (x < y)
0164 return x < z ? XYZVector(0., aVector.Z(), -aVector.Y()) : XYZVector(aVector.Y(), -aVector.X(), 0.);
0165 else
0166 return y < z ? XYZVector(-aVector.Z(), 0., aVector.X()) : XYZVector(aVector.Y(), -aVector.X(), 0.);
0167 }
0168
0169 DEFINE_EDM_PLUGIN(fastsim::InteractionModelFactory, fastsim::MultipleScattering, "fastsim::MultipleScattering");