Implementation of Bremsstrahlung from e+/e- in the tracker layers. More...

Inheritance diagram for fastsim::Bremsstrahlung:

Public Member Functions
	Bremsstrahlung (const std::string &name, const edm::ParameterSet &cfg)
	Constructor. More...

void	interact (Particle &particle, const SimplifiedGeometry &layer, std::vector< std::unique_ptr< Particle > > &secondaries, const RandomEngineAndDistribution &random) override
	Perform the interaction. More...

	~Bremsstrahlung () override
	Default destructor. More...

Public Member Functions inherited from fastsim::InteractionModel
const std::string	getName ()
	Return (unique) name of this interaction. More...

	InteractionModel (std::string name)
	Constructor. More...

virtual void	registerProducts (edm::ProducerBase &producer) const
	In case interaction produces and stores content in the event (e.g. TrackerSimHits). More...

virtual void	storeProducts (edm::Event &iEvent)
	In case interaction produces and stores content in the event (e.g. TrackerSimHits). More...

virtual	~InteractionModel ()
	Default destructor. More...

Private Member Functions
math::XYZTLorentzVector	brem (Particle &particle, double xmin, const RandomEngineAndDistribution &random) const
	Compute Brem photon energy and angles, if any. More...

double	gbteth (const double ener, const double partm, const double efrac, const RandomEngineAndDistribution &random) const
	A universal angular distribution. More...

Private Attributes
double	minPhotonEnergy_
	Cut on minimum energy of bremsstrahlung photons. More...

double	minPhotonEnergyFraction_
	Cut on minimum fraction of particle's energy which has to be carried by photon. More...

double	Z_
	Atomic number of material (usually silicon Z=14) More...

Detailed Description

Implementation of Bremsstrahlung from e+/e- in the tracker layers.

Computes the number, energy and angles of Bremsstrahlung photons emitted by electrons and positrons and modifies e+/e- particle accordingly.

Definition at line 36 of file Bremsstrahlung.cc.

Constructor & Destructor Documentation

fastsim::Bremsstrahlung::Bremsstrahlung	(	const std::string &	name,
		const edm::ParameterSet &	cfg
	)

Constructor.

Definition at line 83 of file Bremsstrahlung.cc.

References edm::ParameterSet::getParameter(), minPhotonEnergy_, minPhotonEnergyFraction_, and Z_.

     : fastsim::InteractionModel(name)
 {
     // Set the minimal photon energy for a Brem from e+/-
     minPhotonEnergy_ = cfg.getParameter<double>("minPhotonEnergy");
     minPhotonEnergyFraction_ = cfg.getParameter<double>("minPhotonEnergyFraction");
     // Material properties
     Z_ = cfg.getParameter<double>("Z");
 }

fastsim::Bremsstrahlung::~Bremsstrahlung ( )

inlineoverride

Default destructor.

Definition at line 43 of file Bremsstrahlung.cc.

References brem(), gbteth(), interact(), random, and TrackerOfflineValidation_Dqm_cff::xmin.

43 {;};

Member Function Documentation

math::XYZTLorentzVector fastsim::Bremsstrahlung::brem	(	fastsim::Particle &	particle,
		double	xmin,
		const RandomEngineAndDistribution &	random
	)		const

private

Compute Brem photon energy and angles, if any.

Parameters

particle	The particle that interacts with the matter.
xmin	Minimum fraction of the particle's energy that has to be converted to a photon.
random	The Random Engine.

Returns: Momentum 4-vector of a bremsstrahlung photon.

Definition at line 172 of file Bremsstrahlung.cc.

References funct::cos(), fastsim::Constants::eMass, JetChargeProducer_cfi::exp, RandomEngineAndDistribution::flatShoot(), gbteth(), cmsBatch::log, M_PI, fastsim::Particle::momentum(), random, funct::sin(), and theta().

Referenced by interact(), and ~Bremsstrahlung().

 {
     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(particle.momentum().E(), fastsim::Constants::eMass, xp, random)
                             * fastsim::Constants::eMass / particle.momentum().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 * particle.momentum().E() * math::XYZTLorentzVector(stheta*cphi, stheta*sphi, ctheta, 1.);  
 }

double fastsim::Bremsstrahlung::gbteth	(	const double	ener,
		const double	partm,
		const double	efrac,
		const RandomEngineAndDistribution &	random
	)		const

private

A universal angular distribution.

Parameters

ener
partm
efrac
random	The Random Engine.

Returns: Theta from universal distribution

Definition at line 202 of file Bremsstrahlung.cc.

References beta, edmIntegrityCheck::d, DEFINE_EDM_PLUGIN, RandomEngineAndDistribution::flatShoot(), cmsBatch::log, M_PI, and Z_.

Referenced by brem(), and ~Bremsstrahlung().

 {
     // Details on implementation here
     // http://www.dnp.fmph.uniba.sk/cernlib/asdoc/geant_html3/node299.html#GBTETH
     // http://svn.cern.ch/guest/AliRoot/tags/v3-07-03/GEANT321/gphys/gbteth.F
 
     const double alfa = 0.625;
     
     const double d = 0.13*(0.8+1.3/Z_)*(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;
 }

void fastsim::Bremsstrahlung::interact	(	fastsim::Particle &	particle,
		const SimplifiedGeometry &	layer,
		std::vector< std::unique_ptr< Particle > > &	secondaries,
		const RandomEngineAndDistribution &	random
	)

overridevirtual

Perform the interaction.

Parameters

particle	The particle that interacts with the matter.
layer	The detector layer that interacts with the particle.
secondaries	Particles that are produced in the interaction (if any).
random	The Random Engine.

Implements fastsim::InteractionModel.

Definition at line 94 of file Bremsstrahlung.cc.

References funct::abs(), brem(), fastsim::SimplifiedGeometry::getThickness(), mps_fire::i, cmsBatch::log, SiStripPI::max, minPhotonEnergy_, minPhotonEnergyFraction_, fastsim::Particle::momentum(), fastsim::Particle::pdgId(), RandomEngineAndDistribution::poissonShoot(), fastsim::Particle::position(), theta(), and TrackerOfflineValidation_Dqm_cff::xmin.

Referenced by ~Bremsstrahlung().

 {
     // only consider electrons and positrons
     if(std::abs(particle.pdgId())!=11)
     {
     return;
     }
     
     double radLengths = layer.getThickness(particle.position(),particle.momentum());
     //
     // no material
     //
     if(radLengths < 1E-10)
     {
     return;
     }
 
     // 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...
     if(radLengths > 4.) 
     {
     particle.momentum().SetXYZT(0.,0.,0.,0.);
     return;
     }
 
     // electron must have more energy than minimum photon energy
     if(particle.momentum().E() - particle.momentum().mass() < minPhotonEnergy_)
     {
     return;
     }
 
     // Hard brem probability with a photon Energy above threshold.
     double xmin = std::max(minPhotonEnergy_/particle.momentum().E(), minPhotonEnergyFraction_);
     if(xmin >=1. || xmin <=0.) 
     {
     return;
     }
 
     // probability to radiate a photon
     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 = random.poissonShoot(bremProba);
     if(nPhotons == 0) 
     {
     return;
     }
 
     // Needed to rotate photons to the lab frame
     double theta = particle.momentum().Theta();
     double phi = particle.momentum().Phi();
     
     // Calculate energy of these photons and add them to the event
     for(unsigned int i=0; i<nPhotons; ++i) 
     {
         // Throw momentum of the photon
         math::XYZTLorentzVector photonMom = brem(particle, xmin, random);
 
         // Check that there is enough energy left.
         if(particle.momentum().E() - particle.momentum().mass() < photonMom.E()) break;
 
         // Rotate to the lab frame
         photonMom = ROOT::Math::RotationZ(phi) * (ROOT::Math::RotationY(theta) * photonMom);
 
         // Add a photon
         secondaries.emplace_back(new fastsim::Particle(22, particle.position(), photonMom));        
         
         // Update the original e+/-
         particle.momentum() -= photonMom;
     }
 }   