#include <BremsstrahlungSimulator.h>

Inheritance diagram for BremsstrahlungSimulator:

Public Member Functions
	BremsstrahlungSimulator (double photonEnergyCut, double photonFractECut, const RandomEngine *engine)
	Constructor. More...

	~BremsstrahlungSimulator ()
	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
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	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 26 of file BremsstrahlungSimulator.h.

Constructor & Destructor Documentation

BremsstrahlungSimulator::BremsstrahlungSimulator	(	double	photonEnergyCut,
		double	photonFractECut,
		const RandomEngine *	engine
	)

Constructor.

Definition at line 7 of file BremsstrahlungSimulator.cc.

References photonEnergy, and photonFractE.

                                                      : 
   MaterialEffectsSimulator(engine) 
 {
 
   // Set the minimal photon energy for a Brem from e+/-
   photonEnergy = photonEnergyCut;
   photonFractE = photonFractECut;
 
 }

BremsstrahlungSimulator::~BremsstrahlungSimulator ( )

inline

Default destructor.

Definition at line 36 of file BremsstrahlungSimulator.h.

36 {}

Member Function Documentation

XYZTLorentzVector BremsstrahlungSimulator::brem ( ParticlePropagator & p ) const

private

Compute Brem photon energy and angles, if any.

Definition at line 72 of file BremsstrahlungSimulator.cc.

References funct::cos(), create_public_lumi_plots::exp, RandomEngine::flatShoot(), gbteth(), create_public_lumi_plots::log, M_PI, phi, MaterialEffectsSimulator::random, funct::sin(), theta(), histoStyle::weight, and xmin.

Referenced by compute().

                                                           {
 
   // This is a simple version (a la PDG) of a Brem generator.
   // It replaces the buggy GEANT3 -> C++ former version.
   // Author : Patrick Janot - 25-Dec-2003
   double emass = 0.0005109990615;
   double xp=0;
   double weight = 0.;
   
   do {
     xp = xmin * std::exp ( -std::log(xmin) * random->flatShoot() );
     weight = 1. - xp + 3./4.*xp*xp;
   } while ( weight < random->flatShoot() );
   
   
   // 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(),emass,xp)*emass/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 BremsstrahlungSimulator::compute ( ParticlePropagator & Particle )

privatevirtual

Generate Bremsstrahlung photons.

Implements MaterialEffectsSimulator.

Definition at line 21 of file BremsstrahlungSimulator.cc.

References MaterialEffectsSimulator::_theUpdatedState, brem(), i, create_public_lumi_plots::log, max(), RawParticle::momentum(), photonEnergy, photonFractE, poisson(), MaterialEffectsSimulator::radLengths, RawParticle::rotate(), and xmin.

 {
 
   // 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...
   // 08/02/06 - pv: increase protection from 1 to 4 X0 for eta>4.8 region
   // if ( radLengths > 1. ) Particle.SetXYZT(0.,0.,0.,0.);
   if ( radLengths > 4. ) Particle.SetXYZT(0.,0.,0.,0.);
 
   // Hard brem probability with a photon Energy above photonEnergy.
   if (Particle.e()<photonEnergy) return;
   xmin = std::max(photonEnergy/Particle.e(),photonFractE);
   if ( xmin >=1. || xmin <=0. ) return;
 
   double bremProba = radLengths * ( 4./3. * std::log(1./xmin)
                   - 4./3. * (1.-xmin)
                   + 1./2. * (1.-xmin*xmin) );
 
   
   // 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;
 
     // Add a photon
     RawParticle thePhoton(22,brem(Particle));    
     thePhoton.rotate(rotY);
     thePhoton.rotate(rotZ);
     _theUpdatedState.push_back(thePhoton);
     
     // Update the original e+/-
     Particle -= thePhoton.momentum();
 
   } 
 }

double BremsstrahlungSimulator::gbteth	(	const double	ener,
		const double	partm,
		const double	efrac
	)		const

private

A universal angular distribution - still from GEANT.

Definition at line 106 of file BremsstrahlungSimulator.cc.

References beta, 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 BremsstrahlungSimulator::poisson ( double ymu )

private

Generate numbers according to a Poisson distribution of mean ymu.

Definition at line 126 of file BremsstrahlungSimulator.cc.

References create_public_lumi_plots::exp, RandomEngine::flatShoot(), n, mix_2012_Summer_inTimeOnly_cff::prob, 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;                                                        
   
 }