Bremsstrahlung.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: 25-Dec-2003
//
// Revision: Class structure modified to match SimplifiedGeometryPropagator
//           Fixed a bug in which particles could radiate more energy than they have
//           S. Kurz, 29 May 2017


namespace fastsim
{

    class Bremsstrahlung : public InteractionModel
    {
        public:
        Bremsstrahlung(const std::string & name, const edm::ParameterSet & cfg);

        ~Bremsstrahlung() override{;};


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

        math::XYZTLorentzVector brem(Particle & particle, double xmin, const RandomEngineAndDistribution & random) const;
        

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

        double minPhotonEnergy_;  
        double minPhotonEnergyFraction_;  
        double Z_;  
    };
}

fastsim::Bremsstrahlung::Bremsstrahlung(const std::string & name,const edm::ParameterSet & cfg)
    : fastsim::InteractionModel(name)
{
    // Set the minimal photon energy for a Brem from e+/-
    minPhotonEnergy_ = cfg.getParameter<double>("minPhotonEnergy");
    minPhotonEnergyFraction_ = cfg.getParameter<double>("minPhotonEnergyFraction");
    // Material properties
    Z_ = cfg.getParameter<double>("Z");
}


void fastsim::Bremsstrahlung::interact(fastsim::Particle & particle, const SimplifiedGeometry & layer,std::vector<std::unique_ptr<fastsim::Particle> > & secondaries,const RandomEngineAndDistribution & random)
{
    // only consider electrons and positrons
    if(std::abs(particle.pdgId())!=11)
    {
    return;
    }
    
    double radLengths = layer.getThickness(particle.position(),particle.momentum());
    //
    // no material
    //
    if(radLengths < 1E-10)
    {
    return;
    }

    // Protection : Just stop the electron if more than 1 radiation lengths.
    // This case corresponds to an electron entering the layer parallel to 
    // the layer axis - no reliable simulation can be done in that case...
    if(radLengths > 4.) 
    {
    particle.momentum().SetXYZT(0.,0.,0.,0.);
    return;
    }

    // electron must have more energy than minimum photon energy
    if(particle.momentum().E() - particle.momentum().mass() < minPhotonEnergy_)
    {
    return;
    }

    // Hard brem probability with a photon Energy above threshold.
    double xmin = std::max(minPhotonEnergy_/particle.momentum().E(), minPhotonEnergyFraction_);
    if(xmin >=1. || xmin <=0.) 
    {
    return;
    }

    // probability to radiate a photon
    double bremProba = radLengths * (4./3. * std::log(1./xmin)
                      - 4./3. * (1.-xmin)
                      + 1./2. * (1.-xmin*xmin));
    
  
    // Number of photons to be radiated.
    unsigned int nPhotons = random.poissonShoot(bremProba);
    if(nPhotons == 0) 
    {
    return;
    }

    // Needed to rotate photons to the lab frame
    double theta = particle.momentum().Theta();
    double phi = particle.momentum().Phi();
    double m2dontchange = particle.momentum().mass()*particle.momentum().mass();
    
    // Calculate energy of these photons and add them to the event
    for(unsigned int i=0; i<nPhotons; ++i) 
    {
        // Throw momentum of the photon
        math::XYZTLorentzVector photonMom = brem(particle, xmin, random);

        // Check that there is enough energy left.
        if(particle.momentum().E() - particle.momentum().mass() < photonMom.E()) break;

        // Rotate to the lab frame
        photonMom = ROOT::Math::RotationZ(phi) * (ROOT::Math::RotationY(theta) * photonMom);

        // Add a photon
        secondaries.emplace_back(new fastsim::Particle(22, particle.position(), photonMom));        
        
        // Update the original e+/-
        particle.momentum() -= photonMom;
        
        // Reset mass to original, since the above codes for decay e->e+gamma, ignoring proton
        particle.momentum().SetXYZT(particle.momentum().px(),particle.momentum().py(), particle.momentum().pz(),sqrt(pow(particle.momentum().P(),2)+m2dontchange));
    }
}   


math::XYZTLorentzVector
fastsim::Bremsstrahlung::brem(fastsim::Particle & particle, double xmin, const RandomEngineAndDistribution & random) const 
{
    double xp=0;
    double weight = 0.;
  
    do{
        xp = xmin * std::exp ( -std::log(xmin) * random.flatShoot() );
        weight = 1. - xp + 3./4.*xp*xp;
    }while(weight < random.flatShoot());
  
  
    // Have photon energy. Now generate angles with respect to the z axis 
    // defined by the incoming particle's momentum.

    // Isotropic in phi
    const double phi = random.flatShoot() * 2. * M_PI;
    // theta from universal distribution
    const double theta = gbteth(particle.momentum().E(), fastsim::Constants::eMass, xp, random)
                            * fastsim::Constants::eMass / particle.momentum().E();
  
    // Make momentum components
    double stheta = std::sin(theta);
    double ctheta = std::cos(theta);
    double sphi   = std::sin(phi);
    double cphi   = std::cos(phi);
  
    return xp * particle.momentum().E() * math::XYZTLorentzVector(stheta*cphi, stheta*sphi, ctheta, 1.);  
}

double
fastsim::Bremsstrahlung::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::Bremsstrahlung,
    "fastsim::Bremsstrahlung"
    );