#include <BremsstrahlungSimulator.h>

Inheritance diagram for BremsstrahlungSimulator:

Public Member Functions
	BremsstrahlungSimulator (double photonEnergyCut, double photonFractECut, const RandomEngine *engine)
	Constructor.
	~BremsstrahlungSimulator ()
	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
double	photonEnergy
	The minimum photon energy to be radiated, in GeV.
double	photonFractE
	The minimum photon fractional energy (wrt that of the electron)
double	xmin
	The fractional photon energy cut (determined from the above two)

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.

{}

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(), funct::exp(), RandomEngine::flatShoot(), gbteth(), funct::log(), M_PI, phi, MaterialEffectsSimulator::random, funct::sin(), theta(), CommonMethods::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 ) [private, virtual]

Generate Bremsstrahlung photons.

Implements MaterialEffectsSimulator.

Definition at line 21 of file BremsstrahlungSimulator.cc.

References MaterialEffectsSimulator::_theUpdatedState, brem(), i, funct::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(), 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 BremsstrahlungSimulator::poisson ( double ymu ) [private]

Generate numbers according to a Poisson distribution of mean ymu.

Definition at line 126 of file BremsstrahlungSimulator.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;                                                        
  
}