fastsim::PairProduction Class Reference

Computes the probability for photons to convert into an e+e- pair in the tracker layer. More...

Inheritance diagram for fastsim::PairProduction:

Public Member Functions
void	interact (fastsim::Particle &particle, const SimplifiedGeometry &layer, std::vector< std::unique_ptr< fastsim::Particle > > &secondaries, const RandomEngineAndDistribution &random) override
	Perform the interaction. More...
	PairProduction (const std::string &name, const edm::ParameterSet &cfg)
	Constructor. More...
	~PairProduction () 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
double	gbteth (double ener, double partm, double efrac, const RandomEngineAndDistribution &random) const
	A universal angular distribution. More...

Private Attributes
double	minPhotonEnergy_
	Cut on minimum energy of photons. More...
double	Z_
	Atomic number of material (usually silicon Z=14) More...

Detailed Description

Computes the probability for photons to convert into an e+e- pair in the tracker layer.

In case, it returns a list of two Secondaries (e+ and e-).

Definition at line 34 of file PairProduction.cc.

Constructor & Destructor Documentation

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

Constructor.

Definition at line 68 of file PairProduction.cc.

References edm::ParameterSet::getParameter(), minPhotonEnergy_, and Z_.

    : fastsim::InteractionModel(name)
{
    // Set the minimal photon energy for possible conversion 
    minPhotonEnergy_ = cfg.getParameter<double>("photonEnergyCut");
    // Material properties
    Z_ = cfg.getParameter<double>("Z");
}

fastsim::PairProduction::~PairProduction ( )

inlineoverride

Default destructor.

Definition at line 41 of file PairProduction.cc.

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

{;};

Member Function Documentation

double fastsim::PairProduction::gbteth ( double  ener, double  partm, double  efrac, const RandomEngineAndDistribution &  random  ) const

private

A universal angular distribution.

Parameters

ener
partm
efrac
random	The Random Engine.

Returns

Theta from universal distribution

Definition at line 182 of file PairProduction.cc.

References pfBoostedDoubleSVAK8TagInfos_cfi::beta, edmIntegrityCheck::d, DEFINE_EDM_PLUGIN, RandomEngineAndDistribution::flatShoot(), cmsBatch::log, M_PI, and Z_.

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

{
    // Details on implementation here
    // http://www.dnp.fmph.uniba.sk/cernlib/asdoc/geant_html3/node299.html#GBTETH
    // http://svn.cern.ch/guest/AliRoot/tags/v3-07-03/GEANT321/gphys/gbteth.F

    const double alfa = 0.625;
    
    const double d = 0.13*(0.8+1.3/Z_)*(100.0+(1.0/ener))*(1.0+efrac);
    const double w1 = 9.0/(9.0+d);
    const double umax = ener*M_PI/partm;
    double u;
    
    do 
    {
        double beta = (random.flatShoot()<=w1) ? alfa : 3.0*alfa;
        u = -std::log(random.flatShoot()*random.flatShoot())/beta;
    }while (u>=umax);

    return u;
}

void fastsim::PairProduction::interact ( fastsim::Particle &  particle, const SimplifiedGeometry &  layer, std::vector< std::unique_ptr< fastsim::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 78 of file PairProduction.cc.

References funct::cos(), fastsim::Constants::eMass, RandomEngineAndDistribution::flatShoot(), gbteth(), fastsim::SimplifiedGeometry::getThickness(), cmsBatch::log, M_PI, SiStripPI::max, minPhotonEnergy_, fastsim::Particle::momentum(), fastsim::Particle::pdgId(), fastsim::Particle::position(), funct::sin(), and mathSSE::sqrt().

Referenced by ~PairProduction().

{

    double eGamma = particle.momentum().e();
    //
    // only consider photons
    //
    if(particle.pdgId()!=22)
    {
    return;
    }
    
    double radLengths = layer.getThickness(particle.position(),particle.momentum());
    //
    // no material
    //
    if(radLengths < 1E-10)
    {
    return;
    }

    //
    // The photon has enough energy to create a pair
    //
    if (eGamma < minPhotonEnergy_)
    {
    return;
    }

    //
    // Probability to convert is 7/9*(dx/X0)
    //
    if(-std::log(random.flatShoot()) > (7./9.)*radLengths)
    {
    return;
    }

    double xe=0;
    double eMass = fastsim::Constants::eMass;
    double xm=eMass/eGamma;
    double weight = 0.;

    // Generate electron energy between emass and eGamma-emass
    do{
        xe = random.flatShoot()*(1.-2.*xm) + xm;
        weight = 1. - 4./3.*xe*(1.-xe);
    }while(weight < random.flatShoot());

    // the electron
    double eElectron = xe * eGamma;
    double tElectron = eElectron-eMass;
    double pElectron = std::sqrt(std::max((eElectron+eMass)*tElectron,0.));

    // the positron
    double ePositron = eGamma-eElectron;
    double tPositron = ePositron-eMass;
    double pPositron = std::sqrt((ePositron+eMass)*tPositron);

    // Generate angles
    double phi  = random.flatShoot()*2.*M_PI;
    double sphi = std::sin(phi);
    double cphi = std::cos(phi);

    double stheta1, stheta2, ctheta1, ctheta2;

    if(eElectron > ePositron){
        double theta1 = gbteth(eElectron,eMass,xe,random)*eMass/eElectron;
        stheta1 = std::sin(theta1);
        ctheta1 = std::cos(theta1);
        stheta2 = stheta1*pElectron/pPositron;
        ctheta2 = std::sqrt(std::max(0.,1.0-(stheta2*stheta2)));
    }else{
        double theta2 = gbteth(ePositron,eMass,xe,random)*eMass/ePositron;
        stheta2 = std::sin(theta2);
        ctheta2 = std::cos(theta2);
        stheta1 = stheta2*pPositron/pElectron;
        ctheta1 = std::sqrt(std::max(0.,1.0-(stheta1*stheta1)));
    }

    //Rotate to the lab frame
    double thetaLab = particle.momentum().Theta();
    double phiLab = particle.momentum().Phi();
    
    // Add a electron
    secondaries.emplace_back(new fastsim::Particle(11,particle.position(),
        math::XYZTLorentzVector(pElectron*stheta1*cphi,
        pElectron*stheta1*sphi,
        pElectron*ctheta1,
        eElectron)));
    secondaries.back()->momentum() = ROOT::Math::RotationZ(phiLab)*(ROOT::Math::RotationY(thetaLab)*secondaries.back()->momentum());

    // Add a positron
    secondaries.emplace_back(new fastsim::Particle(-11,particle.position(),
        math::XYZTLorentzVector(-pPositron*stheta2*cphi,
        -pPositron*stheta2*sphi,
        pPositron*ctheta2,
        ePositron)));
    secondaries.back()->momentum() = ROOT::Math::RotationZ(phiLab)*(ROOT::Math::RotationY(thetaLab)*secondaries.back()->momentum());
    
    // The photon converted
    particle.momentum().SetXYZT(0.,0.,0.,0.);
}

Member Data Documentation

double fastsim::PairProduction::minPhotonEnergy_

private

Cut on minimum energy of photons.

Definition at line 63 of file PairProduction.cc.

Referenced by interact(), and PairProduction().

double fastsim::PairProduction::Z_

private

Atomic number of material (usually silicon Z=14)

Definition at line 64 of file PairProduction.cc.

Referenced by gbteth(), and PairProduction().