fastsim::MultipleScattering Class Reference

Implementation of multiple scattering in the tracker layers. More...

Inheritance diagram for fastsim::MultipleScattering:

Public Member Functions
void	interact (Particle &particle, const SimplifiedGeometry &layer, std::vector< std::unique_ptr< Particle > > &secondaries, const RandomEngineAndDistribution &random) override
	Perform the interaction. More...
	MultipleScattering (const std::string &name, const edm::ParameterSet &cfg)
	Constructor. More...
	~MultipleScattering () override
	Default destructor. More...
Public Member Functions inherited from fastsim::InteractionModel
const std::string	getName ()
	Return (unique) name of this interaction. More...
	InteractionModel (std::string name)
	Constructor. More...
virtual void	registerProducts (edm::ProducesCollector) const
	In case interaction produces and stores content in the event (e.g. TrackerSimHits). More...
virtual void	storeProducts (edm::Event &iEvent)
	In case interaction produces and stores content in the event (e.g. TrackerSimHits). More...
virtual	~InteractionModel ()
	Default destructor. More...

Private Member Functions
XYZVector	orthogonal (const XYZVector &aVector) const
	Return an orthogonal vector. More...

Private Attributes
double	minPt_
	Cut on minimum pT of particle. More...
double	radLenInCm_
	Radiation length of material (usually silicon X0=9.360) More...

Detailed Description

Implementation of multiple scattering in the tracker layers.

Computes the direction change by multiple scattering of a charged particle (assumes constant properties of material).

Definition at line 34 of file MultipleScattering.cc.

Constructor & Destructor Documentation

MultipleScattering()

fastsim::MultipleScattering::MultipleScattering	(	const std::string &	name,
		const edm::ParameterSet &	cfg
	)

Constructor.

Definition at line 63 of file MultipleScattering.cc.

    : fastsim::InteractionModel(name) {
  // Set the minimal pT for interaction
  minPt_ = cfg.getParameter<double>("minPt");
  // Material properties
  radLenInCm_ = cfg.getParameter<double>("radLen");
}

References looper::cfg, minPt_, and radLenInCm_.

~MultipleScattering()

fastsim::MultipleScattering::~MultipleScattering ( )

inlineoverride

Default destructor.

Definition at line 40 of file MultipleScattering.cc.

{ ; };

Member Function Documentation

interact()

void fastsim::MultipleScattering::interact ( fastsim::Particle &  particle, const SimplifiedGeometry &  layer, std::vector< std::unique_ptr< Particle > > &  secondaries, const RandomEngineAndDistribution &  random  )

overridevirtual

Perform the interaction.

Parameters

particle	The particle that interacts with the matter.
layer	The detector layer that interacts with the particle.
secondaries	Particles that are produced in the interaction (if any).
random	The Random Engine.

Implements fastsim::InteractionModel.

Definition at line 71 of file MultipleScattering.cc.

                                                                                      {
  //
  // only charged particles
  //
  if (particle.charge() == 0) {
    return;
  }
 
  double radLengths = layer.getThickness(particle.position(), particle.momentum());
  //
  // no material
  //
  if (radLengths < 1E-10) {
    return;
  }
 
  // particle must have minimum pT
  if (particle.momentum().Pt() < minPt_) {
    return;
  }
 
  double p2 = particle.momentum().Vect().Mag2();
  double m2 = particle.momentum().mass2();
  double e = std::sqrt(p2 + m2);
 
  // This is p*beta
  double pbeta = p2 / e;
 
  // Average multiple scattering angle from Moliere radius
  // The sqrt(2) factor is because of the *space* angle
  double theta0 = 0.0136 / pbeta * particle.charge() * std::sqrt(2. * radLengths) * (1. + 0.038 * std::log(radLengths));
 
  // Generate multiple scattering space angle perpendicular to the particle motion
  double theta = random.gaussShoot(0., theta0);
  // Plus a random rotation angle around the particle motion
  double phi = 2. * M_PI * random.flatShoot();
 
  // The two rotations
  ROOT::Math::AxisAngle rotation1(orthogonal(particle.momentum().Vect()), theta);
  ROOT::Math::AxisAngle rotation2(particle.momentum().Vect(), phi);
  // Rotate!
  XYZVector rotated = rotation2((rotation1(particle.momentum().Vect())));
  particle.momentum().SetXYZT(rotated.X(), rotated.Y(), rotated.Z(), particle.momentum().E());
 
  // Generate mutiple scattering displacements in cm (assuming the detectors
  // are silicon only to determine the thickness) in the directions orthogonal
  // to the vector normal to the surface
  double xp = (cos(phi) * theta / 2. + random.gaussShoot(0., theta0) / sqrt(12.)) * radLengths * radLenInCm_;
  double yp = (sin(phi) * theta / 2. + random.gaussShoot(0., theta0) / sqrt(12.)) * radLengths * radLenInCm_;
 
  // Determine a unitary vector tangent to the surface
  XYZVector normal;
  if (layer.isForward()) {
    normal = XYZVector(0., 0., particle.momentum().Pz() > 0. ? 1. : -1.);
  } else {
    double norm = particle.position().Rho();
    normal = XYZVector(particle.position().X() / norm, particle.position().Y() / norm, 0.);
  }
  XYZVector tangent = orthogonal(normal);
  // The total displacement
  XYZVector delta = xp * tangent + yp * normal.Cross(tangent);
  // Translate!
  particle.position().SetXYZT(particle.position().X() + delta.X(),
                              particle.position().Y() + delta.Y(),
                              particle.position().Z() + delta.Z(),
                              particle.position().T());
 
  // Make sure particle is still physically on the layer (particle moved on a tangent but a barrel layer is bent)
  // Happens only in very very few cases
  if (!layer.isForward() && std::abs(layer.getGeomProperty() - particle.position().Rho()) > 1e-2) {
    // Scale positions so the radii agree
    double scalePos = layer.getGeomProperty() / particle.position().Rho();
    particle.position().SetXYZT(particle.position().X() * scalePos,
                                particle.position().Y() * scalePos,
                                particle.position().Z(),
                                particle.position().T());
 
    // Add a protection in case something goes wrong
    if (std::abs(layer.getGeomProperty() - particle.position().Rho()) > 1e-2) {
      throw cms::Exception("fastsim::MultipleScattering") << "particle no longer on layer's surface";
    }
  }
}

References funct::abs(), fastsim::Particle::charge(), funct::cos(), dumpMFGeometry_cfg::delta, MillePedeFileConverter_cfg::e, Exception, RandomEngineAndDistribution::flatShoot(), RandomEngineAndDistribution::gaussShoot(), fastsim::SimplifiedGeometry::getGeomProperty(), fastsim::SimplifiedGeometry::getThickness(), fastsim::SimplifiedGeometry::isForward(), dqm-mbProfile::log, M_PI, fastsim::Particle::momentum(), p2, fastsim::Particle::position(), funct::sin(), mathSSE::sqrt(), and theta().

orthogonal()

XYZVector fastsim::MultipleScattering::orthogonal ( const XYZVector &  aVector ) const

private

Return an orthogonal vector.

Definition at line 158 of file MultipleScattering.cc.

                                                                              {
  double x = fabs(aVector.X());
  double y = fabs(aVector.Y());
  double z = fabs(aVector.Z());
 
  if (x < y)
    return x < z ? XYZVector(0., aVector.Z(), -aVector.Y()) : XYZVector(aVector.Y(), -aVector.X(), 0.);
  else
    return y < z ? XYZVector(-aVector.Z(), 0., aVector.X()) : XYZVector(aVector.Y(), -aVector.X(), 0.);
}

Member Data Documentation

minPt_

double fastsim::MultipleScattering::minPt_

private

Cut on minimum pT of particle.

Definition at line 58 of file MultipleScattering.cc.

Referenced by MultipleScattering().

radLenInCm_

double fastsim::MultipleScattering::radLenInCm_

private

Radiation length of material (usually silicon X0=9.360)

Definition at line 59 of file MultipleScattering.cc.

Referenced by MultipleScattering().