text: ///////////////////////////////////////////////
  //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;
  }
  
  
  