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::ProducerBase &producer) 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 37 of file MultipleScattering.cc.

Constructor & Destructor Documentation

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

Constructor.

Definition at line 64 of file MultipleScattering.cc.

References edm::ParameterSet::getParameter(), minPt_, and radLenInCm_.

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

fastsim::MultipleScattering::~MultipleScattering ( )

inlineoverride

Default destructor.

Definition at line 44 of file MultipleScattering.cc.

References interact(), orthogonal(), and random.

{;};

Member Function Documentation

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 74 of file MultipleScattering.cc.

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

Referenced by ~MultipleScattering().

{
    //
    // 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";
        }
    }
}

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

private

Return an orthogonal vector.

Definition at line 172 of file MultipleScattering.cc.

References DEFINE_EDM_PLUGIN.

Referenced by interact(), and ~MultipleScattering().

{ 
    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

double fastsim::MultipleScattering::minPt_

private

Cut on minimum pT of particle.

Definition at line 59 of file MultipleScattering.cc.

Referenced by interact(), and MultipleScattering().

double fastsim::MultipleScattering::radLenInCm_

private

Radiation length of material (usually silicon X0=9.360)

Definition at line 60 of file MultipleScattering.cc.

Referenced by interact(), and MultipleScattering().