MaterialEffects/src/MultipleScatteringSimulator.cc

0001 #include "FastSimulation/MaterialEffects/interface/MultipleScatteringSimulator.h"
0002 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0003
0004 #include <cmath>
0005
0006 MultipleScatteringSimulator::MultipleScatteringSimulator(double A, double Z, double density, double radLen)
0007     : MaterialEffectsSimulator(A, Z, density, radLen) {
0008   sqr12 = std::sqrt(12.);
0009 }
0010
0011 void MultipleScatteringSimulator::compute(ParticlePropagator& Particle, RandomEngineAndDistribution const* random) {
0012   double p2 = Particle.particle().Vect().Mag2();
0013   double m2 = Particle.particle().mass() * Particle.particle().mass();
0014   double e = std::sqrt(p2 + m2);
0015
0016   double pbeta = p2 / e;  // This is p*beta
0017
0018   // Average multiple scattering angle from Moliere radius
0019   // The sqrt(2) factor is because of the *space* angle
0020   double theta0 =
0021       0.0136 / pbeta * Particle.particle().charge() * std::sqrt(2. * radLengths) * (1. + 0.038 * std::log(radLengths));
0022
0023   // Generate multiple scattering space angle perpendicular to the particle motion
0024   double theta = random->gaussShoot(0., theta0);
0025   // Plus a random rotation angle around the particle motion
0026   double phi = 2. * 3.14159265358979323 * random->flatShoot();
0027   // The two rotations
0028   RawParticle::Rotation rotation1(orthogonal(Particle.particle().Vect()), theta);
0029   RawParticle::Rotation rotation2(Particle.particle().Vect(), phi);
0030   // Rotate!
0031   Particle.particle().rotate(rotation1);
0032   Particle.particle().rotate(rotation2);
0033
0034   // Generate mutiple scattering displacements in cm (assuming the detectors
0035   // are silicon only to determine the thickness) in the directions orthogonal
0036   // to the vector normal to the surface
0037   double xp = (cos(phi) * theta / 2. + random->gaussShoot(0., theta0) / sqr12) * radLengths * radLenIncm();
0038   double yp = (sin(phi) * theta / 2. + random->gaussShoot(0., theta0) / sqr12) * radLengths * radLenIncm();
0039
0040   // Determine a unitary vector tangent to the surface
0041   // This tangent vector is unitary because "normal" is
0042   // either (0,0,1) in the Endcap  or (x,y,0) in the Barrel !
0043   XYZVector normal(theNormalVector.x(), theNormalVector.y(), theNormalVector.z());
0044   XYZVector tangent = orthogonal(normal);
0045   // The total displacement
0046   XYZVector delta = xp * tangent + yp * normal.Cross(tangent);
0047   // Translate!
0048   Particle.particle().translate(delta);
0049 }