PairProduction.cc

Go to the documentation of this file.

#include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
#include "FastSimulation/SimplifiedGeometryPropagator/interface/Particle.h"
#include "FastSimulation/SimplifiedGeometryPropagator/interface/SimplifiedGeometry.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include <cmath>
#include <memory>

#include <Math/RotationY.h>
#include <Math/RotationZ.h>

#include "FastSimulation/SimplifiedGeometryPropagator/interface/InteractionModelFactory.h"
#include "FastSimulation/SimplifiedGeometryPropagator/interface/InteractionModel.h"
#include "FastSimulation/SimplifiedGeometryPropagator/interface/Constants.h"
#include "DataFormats/Math/interface/LorentzVector.h"


// Author: Patrick Janot
// Date: 24-Dec-2003
//
// Revision: Class structure modified to match SimplifiedGeometryPropagator
//           S. Kurz, 29 May 2017


namespace fastsim
{

    class PairProduction : public InteractionModel
    {
        public:
        PairProduction(const std::string & name,const edm::ParameterSet & cfg);
        
        ~PairProduction() override{;};


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

        double gbteth(double ener,double partm,double efrac, const RandomEngineAndDistribution & random) const;
        
        double minPhotonEnergy_; 
        double Z_; 
    };
}

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


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

    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.);
}

double
fastsim::PairProduction::gbteth(const double ener,
                const double partm,
                const double efrac,
                const RandomEngineAndDistribution & random) const 
{
    // 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;
}


DEFINE_EDM_PLUGIN(
    fastsim::InteractionModelFactory,
    fastsim::PairProduction,
    "fastsim::PairProduction"
    );