MuonBremsstrahlungSimulator.cc

Go to the documentation of this file.

//MuonBremsstrahlungSimulator
// Description: Implementation of Muon bremsstrahlung using the Petrukhin model
//Authors :Sandro Fonseca de Souza and Andre Sznajder (UERJ/Brazil)
// Date: 23-Nov-2010
#include "FastSimulation/MaterialEffects/interface/MuonBremsstrahlungSimulator.h"
#include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FastSimulation/MaterialEffects/interface/PetrukhinModel.h"
#include "FastSimulation/Particle/interface/makeParticle.h"

#include <cmath>
#include <string>
#include <iostream>

#include "TF1.h"

MuonBremsstrahlungSimulator::MuonBremsstrahlungSimulator(
    double A, double Z, double density, double radLen, double photonEnergyCut, double photonFractECut)
    : MaterialEffectsSimulator(A, Z, density, radLen) {
  // Set the minimal photon energy for a Brem from mu+/-
  photonEnergy = photonEnergyCut;
  photonFractE = photonFractECut;
  d = 0.;  //distance
  LogDebug("MuonBremsstrahlungSimulator") << "Starting the MuonBremsstrahlungSimulator" << std::endl;
}

void MuonBremsstrahlungSimulator::compute(ParticlePropagator& Particle, RandomEngineAndDistribution const* random) {
  double NA = 6.022e+23;  //Avogadro's number

  if (radLengths > 4.) {
    Particle.particle().setMomentum(0., 0., 0., 0.);
    deltaPMuon.SetXYZT(0., 0., 0., 0.);
    brem_photon.SetXYZT(0., 0., 0., 0.);
  }

  // Hard brem probability with a photon Energy above photonEnergy.

  double EMuon = Particle.particle().e();  //Muon Energy
  if (EMuon < photonEnergy)
    return;
  xmin = std::max(photonEnergy / EMuon, photonFractE);  //fraction of muon's energy transferred to the photon

  //  xmax = photonEnergy/Particle.particle().e();
  if (xmin >= 1. || xmin <= 0.)
    return;

  xmax = 1.;
  npar = 3;  //Number of parameters

  // create TF1 using a free C function
  f1 = new TF1("f1", PetrukhinFunc, xmin, xmax, npar);
  //Setting parameters
  f1->SetParameters(EMuon, A, Z);
  //d = distance for several materials
  //X0 = radLen
  //d = radLengths * X0(for tracker,yoke,ECAL and HCAL)
  d = radLengths * radLen;
  //Integration
  bremProba = density * d * (NA / A) * (f1->Integral(0., 1.));

  // Number of photons to be radiated.
  unsigned int nPhotons = random->poissonShoot(bremProba);
  _theUpdatedState.reserve(nPhotons);

  if (!nPhotons)
    return;

  //Rotate to the lab frame
  double chi = Particle.particle().theta();
  double psi = Particle.particle().phi();
  RawParticle::RotationZ rotZ(psi);
  RawParticle::RotationY rotY(chi);

  // Energy of these photons
  for (unsigned int i = 0; i < nPhotons; ++i) {
    // Check that there is enough energy left.
    if (Particle.particle().e() < photonEnergy)
      break;
    LogDebug("MuonBremsstrahlungSimulator") << "MuonBremsstrahlungSimulator parameters:" << std::endl;
    LogDebug("MuonBremsstrahlungSimulator") << "xmin-> " << xmin << std::endl;
    LogDebug("MuonBremsstrahlungSimulator") << "Atomic Weight-> " << A << std::endl;
    LogDebug("MuonBremsstrahlungSimulator") << "Density-> " << density << std::endl;
    LogDebug("MuonBremsstrahlungSimulator") << "Distance-> " << d << std::endl;
    LogDebug("MuonBremsstrahlungSimulator") << "bremProba->" << bremProba << std::endl;
    LogDebug("MuonBremsstrahlungSimulator") << "nPhotons->" << nPhotons << std::endl;
    LogDebug("MuonBremsstrahlungSimulator") << " Muon_Energy-> " << EMuon << std::endl;
    LogDebug("MuonBremsstrahlungSimulator") << "X0-> " << radLen << std::endl;
    LogDebug("MuonBremsstrahlungSimulator") << " radLengths-> " << radLengths << std::endl;

    // Add a photon
    RawParticle thePhoton = makeParticle(Particle.particleDataTable(), 22, brem(Particle, random));
    if (thePhoton.E() > 0.) {
      thePhoton.rotate(rotY);
      thePhoton.rotate(rotZ);

      _theUpdatedState.push_back(thePhoton);

      // Update the original mu +/-
      deltaPMuon = Particle.particle().momentum() -= thePhoton.momentum();
      // Information of brem photon
      brem_photon.SetXYZT(thePhoton.Px(), thePhoton.Py(), thePhoton.Pz(), thePhoton.E());

      LogDebug("MuonBremsstrahlungSimulator") << " Muon Bremsstrahlung: photon_energy-> " << thePhoton.E() << std::endl;
      LogDebug("MuonBremsstrahlungSimulator") << "photon_px->" << thePhoton.Px() << std::endl;
      LogDebug("MuonBremsstrahlungSimulator") << "photon_py->" << thePhoton.Py() << std::endl;
      LogDebug("MuonBremsstrahlungSimulator") << "photon_pz->" << thePhoton.Pz() << std::endl;
    }
  }
}
XYZTLorentzVector MuonBremsstrahlungSimulator::brem(ParticlePropagator& pp,
                                                    RandomEngineAndDistribution const* random) const {
  // This is a simple version of a Muon Brem using Petrukhin model .
  //Ref: http://pdg.lbl.gov/2008/AtomicNuclearProperties/adndt.pdf
  double mumass = 0.105658367;  //mu mass  (GeV/c^2)

  // RANDOM_NUMBER_ERROR
  // Random number should be generated by the engines from the
  // RandomNumberGeneratorService. This appears to use the global
  // engine in ROOT. This is not thread safe unless the module using
  // it is a one module and declares a shared resource and all
  // other modules using it also declare the same shared resource.
  // This also breaks replay.
  // It might be that this is never used because of the std::cout
  // statement 2 lines below. I've never seen or heard complaints
  // about that vebosity ....
  double xp = f1->GetRandom();
  LogDebug("MuonBremsstrahlungSimulator") << "MuonBremsstrahlungSimulator: xp->" << xp << std::endl;
  std::cout << "MuonBremsstrahlungSimulator: xp->" << xp << std::endl;

  // 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(pp.particle().e(), mumass, xp, random) * mumass / pp.particle().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 * pp.particle().e() * XYZTLorentzVector(stheta * cphi, stheta * sphi, ctheta, 1.);
}
double MuonBremsstrahlungSimulator::gbteth(const double ener,
                                           const double partm,
                                           const double efrac,
                                           RandomEngineAndDistribution const* random) const {
  const double alfa = 0.625;

  const double d = 0.13 * (0.8 + 1.3 / theZ()) * (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;
}