#include <MuonBremsstrahlungSimulator.h>

Inheritance diagram for MuonBremsstrahlungSimulator:

Public Member Functions
const XYZTLorentzVector &	deltaP_BremMuon () const
const XYZTLorentzVector &	deltaP_BremPhoton () const
	MuonBremsstrahlungSimulator (const RandomEngine *engine, double A, double Z, double density, double radLen, double photonEnergyCut, double photonFractECut)
	Constructor.
	~MuonBremsstrahlungSimulator ()
	Default destructor.
Private Member Functions
XYZTLorentzVector	brem (ParticlePropagator &p) const
	Compute Brem photon energy and angles, if any.
void	compute (ParticlePropagator &Particle)
	Generate Bremsstrahlung photons.
double	gbteth (const double ener, const double partm, const double efrac) const
	A universal angular distribution - still from GEANT.
unsigned int	poisson (double ymu)
	Generate numbers according to a Poisson distribution of mean ymu.
Private Attributes
XYZTLorentzVector	brem_photon
double	bremProba
double	d
XYZTLorentzVector	deltaPMuon
TF1 *	f1
int	npar
double	photonEnergy
	The minimum photon energy to be radiated, in GeV.
double	photonFractE
	The minimum photon fractional energy (wrt that of the electron)
double	rand
double	xmax
double	xmin
	The fractional photon energy cut (determined from the above two)

Detailed Description

Definition at line 31 of file MuonBremsstrahlungSimulator.h.

Constructor & Destructor Documentation

MuonBremsstrahlungSimulator::MuonBremsstrahlungSimulator	(	const RandomEngine *	engine,
		double	A,
		double	Z,
		double	density,
		double	radLen,
		double	photonEnergyCut,
		double	photonFractECut
	)

Constructor.

Definition at line 22 of file MuonBremsstrahlungSimulator.cc.

References d, LogDebug, photonEnergy, and photonFractE.

                                                                                                        : 
  MaterialEffectsSimulator(engine,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;
}

MuonBremsstrahlungSimulator::~MuonBremsstrahlungSimulator ( ) [inline]

Default destructor.

Definition at line 41 of file MuonBremsstrahlungSimulator.h.

{}

Member Function Documentation

XYZTLorentzVector MuonBremsstrahlungSimulator::brem ( ParticlePropagator & p ) const [private]

Compute Brem photon energy and angles, if any.

Definition at line 131 of file MuonBremsstrahlungSimulator.cc.

References funct::cos(), f1, RandomEngine::flatShoot(), gbteth(), LogDebug, M_PI, phi, MaterialEffectsSimulator::random, funct::sin(), and theta().

Referenced by compute().

                                                              {

  // This is a simple version of a Muon Brem using Petrukhin model .
  //Ref: https://pdg.lbl.gov/2008/AtomicNuclearProperties/adndt.pdf 
  double mumass = 0.105658367;//mu mass  (GeV/c^2)
  double xp =0.;  
 //+++++++++++++++++++++++++++++++++++++++++++++++++++ 
  // Ref: Framework    
  //https://hypernews.cern.ch/HyperNews/CMS/get/swDevelopment/1969/1.html
  // https://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/CMSSW/DQM/DTMonitorClient/src/DTResolutionTest.cc?view=markup
  try {
 
    xp = f1->GetRandom();
  } catch (...) {
    
    edm::LogInfo ("MuonBremsstrahlungSimulator") << "Exception when fitting..."; 
  }
  LogDebug("MuonBremsstrahlungSimulator")<<  "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.e(),mumass,xp)*mumass/pp.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.e() * XYZTLorentzVector(stheta*cphi,stheta*sphi,ctheta,1.);
 
}

void MuonBremsstrahlungSimulator::compute ( ParticlePropagator & Particle ) [private, virtual]

Generate Bremsstrahlung photons.

Implements MaterialEffectsSimulator.

Definition at line 38 of file MuonBremsstrahlungSimulator.cc.

References MaterialEffectsSimulator::_theUpdatedState, MaterialEffectsSimulator::A, brem(), brem_photon, bremProba, d, deltaPMuon, MaterialEffectsSimulator::density, f1, i, LogDebug, max(), RawParticle::momentum(), npar, PetrukhinFunc(), photonEnergy, photonFractE, poisson(), MaterialEffectsSimulator::radLen, MaterialEffectsSimulator::radLengths, RawParticle::rotate(), xmax, xmin, and MaterialEffectsSimulator::Z.

{

  double NA = 6.022e+23;  //Avogadro's number

  if ( radLengths > 4. )  {
    Particle.SetXYZT(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.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.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 = poisson(bremProba);
  _theUpdatedState.reserve(nPhotons);
 
 
  if ( !nPhotons ) return;
 
  //Rotate to the lab frame
  double chi = Particle.theta();
  double psi = 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.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(22,brem(Particle));    
    if (thePhoton.E()>0.){

    thePhoton.rotate(rotY);
    thePhoton.rotate(rotZ);
    
    _theUpdatedState.push_back(thePhoton);
        
    // Update the original mu +/-
    deltaPMuon = Particle -= 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;
    
    } 
    
  }

}

const XYZTLorentzVector& MuonBremsstrahlungSimulator::deltaP_BremMuon ( ) const [inline]

Definition at line 45 of file MuonBremsstrahlungSimulator.h.

References deltaPMuon.

{ return deltaPMuon ; }

const XYZTLorentzVector& MuonBremsstrahlungSimulator::deltaP_BremPhoton ( ) const [inline]

Definition at line 48 of file MuonBremsstrahlungSimulator.h.

References brem_photon.

{ return brem_photon ; }

double MuonBremsstrahlungSimulator::gbteth	(	const double	ener,
		const double	partm,
		const double	efrac
	)		const `[private]`

A universal angular distribution - still from GEANT.

Definition at line 171 of file MuonBremsstrahlungSimulator.cc.

References beta, d, RandomEngine::flatShoot(), funct::log(), M_PI, MaterialEffectsSimulator::random, and MaterialEffectsSimulator::theZ().

Referenced by brem().

                                                          {
  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;
}

unsigned int MuonBremsstrahlungSimulator::poisson ( double ymu ) [private]

Generate numbers according to a Poisson distribution of mean ymu.

Definition at line 191 of file MuonBremsstrahlungSimulator.cc.

References funct::exp(), RandomEngine::flatShoot(), n, MaterialEffectsSimulator::random, and x.

Referenced by compute().

                                               {

  unsigned int n = 0;
  double prob = std::exp(-ymu);
  double proba = prob;
  double x = random->flatShoot();
  
  while ( proba <= x ) {
    prob *= ymu / double(++n);
    proba += prob;
  }
  
  return n;                                                        
  
}