SimplifiedGeometryPropagator/plugins/MultipleScattering.cc

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 // Author: Patrick Janot
0020 // Date: 8-Jan-2004
0021 //
0022 // Revision: Class structure modified to match SimplifiedGeometryPropagator
0023 //           Fixed a bug in which particles could no longer be on the layer after scattering
0024 //           S. Kurz, 29 May 2017
0025 //////////////////////////////////////////////////////////
0026
0027 typedef math::XYZVector XYZVector;
0028
0029 namespace fastsim {
0030   //! Implementation of multiple scattering in the tracker layers.
0031   /*!
0032         Computes the direction change by multiple scattering of a charged particle (assumes constant properties of material).
0033     */
0034   class MultipleScattering : public InteractionModel {
0035   public:
0036     //! Constructor.
0037     MultipleScattering(const std::string& name, const edm::ParameterSet& cfg);
0038
0039     //! Default destructor.
0040     ~MultipleScattering() override { ; };
0041
0042     //! Perform the interaction.
0043     /*!
0044             \param particle The particle that interacts with the matter.
0045             \param layer The detector layer that interacts with the particle.
0046             \param secondaries Particles that are produced in the interaction (if any).
0047             \param random The Random Engine.
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     //! Return an orthogonal vector.
0056     XYZVector orthogonal(const XYZVector& aVector) const;
0057
0058     double minPt_;       //!< Cut on minimum pT of particle
0059     double radLenInCm_;  //!< Radiation length of material (usually silicon X0=9.360)
0060   };
0061 }  // namespace fastsim
0062
0063 fastsim::MultipleScattering::MultipleScattering(const std::string& name, const edm::ParameterSet& cfg)
0064     : fastsim::InteractionModel(name) {
0065   // Set the minimal pT for interaction
0066   minPt_ = cfg.getParameter<double>("minPt");
0067   // Material properties
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   // only charged particles
0077   //
0078   if (particle.charge() == 0) {
0079     return;
0080   }
0081
0082   double radLengths = layer.getThickness(particle.position(), particle.momentum());
0083   //
0084   // no material
0085   //
0086   if (radLengths < 1E-10) {
0087     return;
0088   }
0089
0090   // particle must have minimum pT
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   // This is p*beta
0100   double pbeta = p2 / e;
0101
0102   // Average multiple scattering angle from Moliere radius
0103   // The sqrt(2) factor is because of the *space* angle
0104   double theta0 = 0.0136 / pbeta * particle.charge() * std::sqrt(2. * radLengths) * (1. + 0.038 * std::log(radLengths));
0105
0106   // Generate multiple scattering space angle perpendicular to the particle motion
0107   double theta = random.gaussShoot(0., theta0);
0108   // Plus a random rotation angle around the particle motion
0109   double phi = 2. * M_PI * random.flatShoot();
0110
0111   // The two rotations
0112   ROOT::Math::AxisAngle rotation1(orthogonal(particle.momentum().Vect()), theta);
0113   ROOT::Math::AxisAngle rotation2(particle.momentum().Vect(), phi);
0114   // Rotate!
0115   XYZVector rotated = rotation2((rotation1(particle.momentum().Vect())));
0116   particle.momentum().SetXYZT(rotated.X(), rotated.Y(), rotated.Z(), particle.momentum().E());
0117
0118   // Generate mutiple scattering displacements in cm (assuming the detectors
0119   // are silicon only to determine the thickness) in the directions orthogonal
0120   // to the vector normal to the surface
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   // Determine a unitary vector tangent to the surface
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   // The total displacement
0134   XYZVector delta = xp * tangent + yp * normal.Cross(tangent);
0135   // Translate!
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   // Make sure particle is still physically on the layer (particle moved on a tangent but a barrel layer is bent)
0142   // Happens only in very very few cases
0143   if (!layer.isForward() && std::abs(layer.getGeomProperty() - particle.position().Rho()) > 1e-2) {
0144     // Scale positions so the radii agree
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     // Add a protection in case something goes wrong
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");