#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. More...

	~MuonBremsstrahlungSimulator ()
	Default destructor. More...

Public Member Functions inherited from MaterialEffectsSimulator
RHEP_const_iter	beginDaughters () const
	Returns const iterator to the beginning of the daughters list. More...

int	closestDaughterId ()
	The id of the closest charged daughter (filled for nuclear interactions only) More...

double	eMass () const
	Electron mass in GeV/c2. More...

RHEP_const_iter	endDaughters () const
	Returns const iterator to the end of the daughters list. More...

double	excitE () const
	Mean excitation energy (in GeV) More...

	MaterialEffectsSimulator (const RandomEngine *engine, double A=28.0855, double Z=14.0000, double density=2.329, double radLen=9.360)

unsigned	nDaughters () const
	Returns the number of daughters. More...

XYZVector	orthogonal (const XYZVector &) const
	A vector orthogonal to another one (because it's not in XYZTLorentzVector) More...

double	radLenIncm () const
	One radiation length in cm. More...

double	rho () const
	Density in g/cm3. More...

void	setNormalVector (const GlobalVector &normal)
	Sets the vector normal to the surface traversed. More...

double	theA () const
	A. More...

double	theZ () const
	Z. More...

void	updateState (ParticlePropagator &myTrack, double radlen)
	Compute the material effect (calls the sub class) More...

virtual	~MaterialEffectsSimulator ()

Private Member Functions
XYZTLorentzVector	brem (ParticlePropagator &p) const
	Compute Brem photon energy and angles, if any. More...

void	compute (ParticlePropagator &Particle)
	Generate Bremsstrahlung photons. More...

double	gbteth (const double ener, const double partm, const double efrac) const
	A universal angular distribution - still from GEANT. More...

unsigned int	poisson (double ymu)
	Generate numbers according to a Poisson distribution of mean ymu. More...

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. More...

double	photonFractE
	The minimum photon fractional energy (wrt that of the electron) More...

double	rand

double	xmax

double	xmin
	The fractional photon energy cut (determined from the above two) More...

Additional Inherited Members
Public Types inherited from MaterialEffectsSimulator
typedef std::vector < RawParticle > ::const_iterator	RHEP_const_iter

Protected Attributes inherited from MaterialEffectsSimulator
std::vector< RawParticle >	_theUpdatedState

double	A

double	density

double	radLen

double	radLengths

const RandomEngine *	random

int	theClosestChargedDaughterId

GlobalVector	theNormalVector

double	Z

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.

41 {}

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: http://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
   // http://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 )

privatevirtual

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.

45 { return deltaPMuon ; }

MuonBremsstrahlungSimulator::deltaPMuon

XYZTLorentzVector deltaPMuon

Definition: MuonBremsstrahlungSimulator.h:82

const XYZTLorentzVector& MuonBremsstrahlungSimulator::deltaP_BremPhoton ( ) const

inline

Definition at line 48 of file MuonBremsstrahlungSimulator.h.

References brem_photon.

48 { return brem_photon ; }

MuonBremsstrahlungSimulator::brem_photon

XYZTLorentzVector brem_photon

Definition: MuonBremsstrahlungSimulator.h:84

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(), create_public_lumi_plots::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 create_public_lumi_plots::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;                                                        
   
 }