MuonBremsstrahlung.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"

#include <TMath.h>
#include <TF1.h>

// Authors: Sandro Fonseca de Souza and Andre Sznajder (UERJ/Brazil)
// Date: 23-Nov-2010
//
// Revision: Class structure modified to match SimplifiedGeometryPropagator
//           S. Sekmen, 18 May 2017
//
// Revision: Code very buggy, PetrukhinFunc return negative values, double bremProba wasn't properly defined etc.
//           Should be all fixed by now
//           S. Kurz, 23 May 2017

namespace fastsim {

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

    ~MuonBremsstrahlung() 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;

    static double PetrukhinFunc(double *x, double *p);

    TF1 *Petrfunc;                    
    double minPhotonEnergy_;          
    double minPhotonEnergyFraction_;  
    double density_;                  
    double radLenInCm_;               
    double A_;                        
    double Z_;                        
  };
}  // namespace fastsim

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

void fastsim::MuonBremsstrahlung::interact(fastsim::Particle &particle,
                                           const SimplifiedGeometry &layer,
                                           std::vector<std::unique_ptr<fastsim::Particle> > &secondaries,
                                           const RandomEngineAndDistribution &random) {
  // only consider muons
  if (std::abs(particle.pdgId()) != 13) {
    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;
  }

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

  // Min fraction of muon's energy transferred to the photon
  double xmin = std::max(minPhotonEnergy_ / particle.momentum().E(), minPhotonEnergyFraction_);
  // Hard brem probability with a photon Energy above threshold.
  if (xmin >= 1. || xmin <= 0.) {
    return;
  }

  // Max fraction of muon's energy transferred to the photon
  double xmax = 1.;

  // create TF1 using a free C function
  Petrfunc = new TF1("Petrfunc", PetrukhinFunc, xmin, xmax, 3);
  // Setting parameters
  Petrfunc->SetParameters(particle.momentum().E(), A_, Z_);
  // d = distance for several materials
  // X0 = radLen
  // d = radLengths * X0 (for tracker, yoke, ECAL and HCAL)
  double distance = radLengths * radLenInCm_;

  // Integration
  // Fixed previous version which used Petrfunc->Integral(0.,1.) -> does not make sense
  double bremProba = density_ * distance * (fastsim::Constants::NA / A_) * (Petrfunc->Integral(xmin, xmax));
  if (bremProba < 1E-10) {
    return;
  }

  // Number of photons to be radiated.
  unsigned int nPhotons = random.poissonShoot(bremProba);
  if (nPhotons == 0) {
    return;
  }

  //Rotate to the lab frame
  double theta = particle.momentum().Theta();
  double phi = particle.momentum().Phi();

  // Energy of these photons
  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;
  }
}

math::XYZTLorentzVector fastsim::MuonBremsstrahlung::brem(fastsim::Particle &particle,
                                                          double xmin,
                                                          const RandomEngineAndDistribution &random) const {
  // This is a simple version of a Muon Brem using Petrukhin model.
  // Ref: http://pdg.lbl.gov/2008/AtomicNuclearProperties/adndt.pdf
  double xp = Petrfunc->GetRandom();

  // 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::muMass, xp, random) *
                       fastsim::Constants::muMass / 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::MuonBremsstrahlung::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;
}

double fastsim::MuonBremsstrahlung::PetrukhinFunc(double *x, double *p) {
  // Function independent variable
  double nu = x[0];  // fraction of muon's energy transferred to the photon

  // Parameters
  double E = p[0];  // Muon Energy (in GeV)
  double A = p[1];  // Atomic weight
  double Z = p[2];  // Atomic number

  /*
        // Function of Muom Brem using  nuclear screening correction
        // Ref: http://pdg.lbl.gov/2008/AtomicNuclearProperties/adndt.pdf
        // http://geant4.cern.ch/G4UsersDocuments/UsersGuides/PhysicsReferenceManual/html/node48.html

        // Physical constants
        double B = 182.7;
        double ee = sqrt(2.7181) ; // sqrt(e)
        double ZZ=  pow( Z,-1./3.); // Z^-1/3
        double emass = 0.0005109990615;  // electron mass (GeV/c^2)
        double mumass = 0.105658367;//mu mass  (GeV/c^2)

        double re = 2.817940285e-13;// Classical electron radius (Units: cm)
        double alpha = 1./137.03599976; // fine structure constant
        double Dn = 1.54* (pow(A,0.27));
        double constant =  pow((2.0 * Z * emass/mumass * re ),2.0);

        double delta = (mumass * mumass * nu) /(2.* E * (1.- nu)); 

        double Delta_n = TMath::Log(Dn / (1.+ delta *( Dn * ee -2.)/ mumass)); //nuclear screening correction 

        double Phi = TMath::Log((B * mumass * ZZ / emass)/ (1.+ delta * ee * B * ZZ  / emass)) - Delta_n;//phi(delta)

        // Diff. Cross Section for Muon Brem from a screened nuclear (Equation 16: REF: LBNL-44742)
        double f = alpha * constant *(4./3.-4./3.*nu + nu*nu)*Phi/nu;
    */

  // Function for Muon Brem Xsec from G4

  // Physical constants
  double B = 183.;
  double Bl = 1429.;
  double ee = 1.64872;            // sqrt(e)
  double Z13 = pow(Z, -1. / 3.);  // Z^-1/3
  double Z23 = pow(Z, -2. / 3.);  // Z^-2/3

  // Original values of paper
  double emass = 0.0005109990615;  // electron mass (GeV/c^2)
  double mumass = 0.105658367;     // muon mass  (GeV/c^2)
  // double re = 2.817940285e-13;     // Classical electron radius (Units: cm)
  double alpha = 0.00729735;      // 1./137.03599976; // fine structure constant
  double constant = 1.85736e-30;  // pow( ( emass / mumass * re ) , 2.0);

  // Use nomenclature from reference -> Follow those formula step by step
  if (nu * E >= E - mumass)
    return 0;

  double Dn = 1.54 * (pow(A, 0.27));
  double Dnl = pow(Dn, (1. - 1. / Z));

  double delta = (mumass * mumass * nu) / (2. * E * (1. - nu));

  double Phi_n = TMath::Log(B * Z13 * (mumass + delta * (Dnl * ee - 2)) / (Dnl * (emass + delta * ee * B * Z13)));
  if (Phi_n < 0)
    Phi_n = 0;

  double Phi_e =
      TMath::Log((Bl * Z23 * mumass) / (1. + delta * mumass / (emass * emass * ee)) / (emass + delta * ee * Bl * Z23));
  if (Phi_e < 0 || nu * E >= E / (1. + (mumass * mumass / (2. * emass * E))))
    Phi_e = 0;

  // Diff. Cross Section for Muon Brem from G4 (without NA/A factor)
  double f = 16. / 3. * alpha * constant * Z * (Z * Phi_n + Phi_e) * (1. / nu) * (1. - nu + 0.75 * nu * nu);

  return f;
}

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