#include <MuonBremsstrahlungSimulator.h>

Inheritance diagram for MuonBremsstrahlungSimulator:

Public Member Functions
const XYZTLorentzVector &	deltaP_BremMuon () const

const XYZTLorentzVector &	deltaP_BremPhoton () const

	MuonBremsstrahlungSimulator (double A, double Z, double density, double radLen, double photonEnergyCut, double photonFractECut)
	Constructor. More...

	~MuonBremsstrahlungSimulator () override
	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 (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...

virtual void	save ()
	Used by NuclearInteractionSimulator to save last sampled event. 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, RandomEngineAndDistribution const *)
	Compute the material effect (calls the sub class) More...

virtual	~MaterialEffectsSimulator ()

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

void	compute (ParticlePropagator &Particle, RandomEngineAndDistribution const *) override
	Generate Bremsstrahlung photons. More...

double	gbteth (const double ener, const double partm, const double efrac, RandomEngineAndDistribution const *) 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

int	theClosestChargedDaughterId

GlobalVector	theNormalVector

double	Z

Detailed Description

Definition at line 29 of file MuonBremsstrahlungSimulator.h.

Constructor & Destructor Documentation

MuonBremsstrahlungSimulator::MuonBremsstrahlungSimulator	(	double	A,
		double	Z,
		double	density,
		double	radLen,
		double	photonEnergyCut,
		double	photonFractECut
	)

Constructor.

Definition at line 20 of file MuonBremsstrahlungSimulator.cc.

References d, LogDebug, photonEnergy, fastSimProducer_cff::photonEnergyCut, and photonFractE.

     : MaterialEffectsSimulator(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 ( )

inlineoverride

Default destructor.

Definition at line 36 of file MuonBremsstrahlungSimulator.h.

36 {}

Member Function Documentation

XYZTLorentzVector MuonBremsstrahlungSimulator::brem	(	ParticlePropagator &	p,
		RandomEngineAndDistribution const *	random
	)		const

private

Compute Brem photon energy and angles, if any.

Definition at line 116 of file MuonBremsstrahlungSimulator.cc.

References funct::cos(), gather_cfg::cout, RawParticle::e(), f1, RandomEngineAndDistribution::flatShoot(), gbteth(), LogDebug, M_PI, BaseParticlePropagator::particle(), phi, 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)
 
   // RANDOM_NUMBER_ERROR
   // Random number should be generated by the engines from the
   // RandomNumberGeneratorService. This appears to use the global
   // engine in ROOT. This is not thread safe unless the module using
   // it is a one module and declares a shared resource and all
   // other modules using it also declare the same shared resource.
   // This also breaks replay.
   // It might be that this is never used because of the std::cout
   // statement 2 lines below. I've never seen or heard complaints
   // about that vebosity ....
   double xp = f1->GetRandom();
   LogDebug("MuonBremsstrahlungSimulator") << "MuonBremsstrahlungSimulator: xp->" << xp << std::endl;
   std::cout << "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.particle().e(), mumass, xp, random) * mumass / pp.particle().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.particle().e() * XYZTLorentzVector(stheta * cphi, stheta * sphi, ctheta, 1.);
 }

void MuonBremsstrahlungSimulator::compute	(	ParticlePropagator &	Particle,
		RandomEngineAndDistribution const *	random
	)

overrideprivatevirtual

Generate Bremsstrahlung photons.

Implements MaterialEffectsSimulator.

Definition at line 31 of file MuonBremsstrahlungSimulator.cc.

References MaterialEffectsSimulator::_theUpdatedState, MaterialEffectsSimulator::A, brem(), brem_photon, bremProba, d, deltaPMuon, MaterialEffectsSimulator::density, RawParticle::e(), RawParticle::E(), f1, mps_fire::i, LogDebug, makeParticle(), SiStripPI::max, RawParticle::momentum(), fastsim::Constants::NA, npar, BaseParticlePropagator::particle(), ParticlePropagator::particleDataTable(), PetrukhinFunc(), RawParticle::phi(), photonEnergy, photonFractE, RandomEngineAndDistribution::poissonShoot(), RawParticle::Px(), RawParticle::Py(), RawParticle::Pz(), MaterialEffectsSimulator::radLen, MaterialEffectsSimulator::radLengths, RawParticle::rotate(), RawParticle::setMomentum(), RawParticle::theta(), xmax, xmin, and MaterialEffectsSimulator::Z.

                                                                                                                  {
   double NA = 6.022e+23;  //Avogadro's number
 
   if (radLengths > 4.) {
     Particle.particle().setMomentum(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.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.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 = random->poissonShoot(bremProba);
   _theUpdatedState.reserve(nPhotons);
 
   if (!nPhotons)
     return;
 
   //Rotate to the lab frame
   double chi = Particle.particle().theta();
   double psi = Particle.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.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 = makeParticle(Particle.particleDataTable(), 22, brem(Particle, random));
     if (thePhoton.E() > 0.) {
       thePhoton.rotate(rotY);
       thePhoton.rotate(rotZ);
 
       _theUpdatedState.push_back(thePhoton);
 
       // Update the original mu +/-
       deltaPMuon = Particle.particle().momentum() -= 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 39 of file MuonBremsstrahlungSimulator.h.

References deltaPMuon.

39 { return deltaPMuon; }

MuonBremsstrahlungSimulator::deltaPMuon

XYZTLorentzVector deltaPMuon

Definition: MuonBremsstrahlungSimulator.h:72

const XYZTLorentzVector& MuonBremsstrahlungSimulator::deltaP_BremPhoton ( ) const

inline

Definition at line 42 of file MuonBremsstrahlungSimulator.h.

References brem_photon.

42 { return brem_photon; }

MuonBremsstrahlungSimulator::brem_photon

XYZTLorentzVector brem_photon

Definition: MuonBremsstrahlungSimulator.h:74

double MuonBremsstrahlungSimulator::gbteth	(	const double	ener,
		const double	partm,
		const double	efrac,
		RandomEngineAndDistribution const *	random
	)		const

private

A universal angular distribution - still from GEANT.

Definition at line 153 of file MuonBremsstrahlungSimulator.cc.

References zMuMuMuonUserData::beta, d, RandomEngineAndDistribution::flatShoot(), dqm-mbProfile::log, M_PI, 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;
 }