Implementation of Bremsstrahlung from mu+/mu- in the tracker layers based on a Petrukhin Model (nuclear screening correction). More...

Inheritance diagram for fastsim::MuonBremsstrahlung:

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

	MuonBremsstrahlung (const std::string &name, const edm::ParameterSet &cfg)
	Constructor. More...

	~MuonBremsstrahlung () 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::ProducesCollector) 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...

Static Private Member Functions
static double	PetrukhinFunc (double x, double p)
	Petrukhin Function: Returns cross section using nuclear-electron screening correction from G4 style. More...

Private Attributes
double	A_
	Atomic weight of material (usually silicon A=28.0855) More...

double	density_
	Density of material (usually silicon rho=2.329) More...

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

TF1 *	Petrfunc
	The Petrukhin Function. More...

double	radLenInCm_
	Radiation length of material (usually silicon X0=9.360) More...

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

Detailed Description

Implementation of Bremsstrahlung from mu+/mu- in the tracker layers based on a Petrukhin Model (nuclear screening correction).

Computes the number, energy and angles of Bremsstrahlung photons emitted by muons and modifies mu+/mu- particle accordingly.

Definition at line 39 of file MuonBremsstrahlung.cc.

Constructor & Destructor Documentation

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

Constructor.

Definition at line 95 of file MuonBremsstrahlung.cc.

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

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

fastsim::MuonBremsstrahlung::~MuonBremsstrahlung ( )

inlineoverride

Default destructor.

Definition at line 45 of file MuonBremsstrahlung.cc.

45 { ; };

Member Function Documentation

math::XYZTLorentzVector fastsim::MuonBremsstrahlung::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 193 of file MuonBremsstrahlung.cc.

References funct::cos(), RandomEngineAndDistribution::flatShoot(), M_PI, fastsim::Particle::momentum(), fastsim::Constants::muMass, funct::sin(), and theta().

                                                                                                            {
   // This is a simple version of a Muon Brem using Petrukhin model.
   // Ref: http://pdg.lbl.gov/2008/AtomicNuclearProperties/adndt.pdf
   double xp = Petrfunc->GetRandom();
 
   // 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::muMass, xp, random) *
                        fastsim::Constants::muMass / 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::MuonBremsstrahlung::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 218 of file MuonBremsstrahlung.cc.

References HLT_FULL_cff::beta, ztail::d, RandomEngineAndDistribution::flatShoot(), log, and M_PI.

                                                                                             {
   // 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::MuonBremsstrahlung::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 107 of file MuonBremsstrahlung.cc.

References funct::abs(), HLT_FULL_cff::distance, fastsim::SimplifiedGeometry::getThickness(), mps_fire::i, SiStripPI::max, fastsim::Particle::momentum(), fastsim::Constants::NA, fastsim::Particle::pdgId(), PetrukhinFunc(), RandomEngineAndDistribution::poissonShoot(), fastsim::Particle::position(), theta(), hlt_dqm_clientPB-live_cfg::xmax, and hlt_dqm_clientPB-live_cfg::xmin.

                                                                                       {
   // only consider muons
   if (std::abs(particle.pdgId()) != 13) {
     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;
   }
 
   // muon must have more energy than minimum photon energy
   if (particle.momentum().E() - particle.momentum().mass() < minPhotonEnergy_) {
     return;
   }
 
   // Min fraction of muon's energy transferred to the photon
   double xmin = std::max(minPhotonEnergy_ / particle.momentum().E(), minPhotonEnergyFraction_);
   // Hard brem probability with a photon Energy above threshold.
   if (xmin >= 1. || xmin <= 0.) {
     return;
   }
 
   // Max fraction of muon's energy transferred to the photon
   double xmax = 1.;
 
   // create TF1 using a free C function
   Petrfunc = new TF1("Petrfunc", PetrukhinFunc, xmin, xmax, 3);
   // Setting parameters
   Petrfunc->SetParameters(particle.momentum().E(), A_, Z_);
   // d = distance for several materials
   // X0 = radLen
   // d = radLengths * X0 (for tracker, yoke, ECAL and HCAL)
   double distance = radLengths * radLenInCm_;
 
   // Integration
   // Fixed previous version which used Petrfunc->Integral(0.,1.) -> does not make sense
   double bremProba = density_ * distance * (fastsim::Constants::NA / A_) * (Petrfunc->Integral(xmin, xmax));
   if (bremProba < 1E-10) {
     return;
   }
 
   // Number of photons to be radiated.
   unsigned int nPhotons = random.poissonShoot(bremProba);
   if (nPhotons == 0) {
     return;
   }
 
   //Rotate to the lab frame
   double theta = particle.momentum().Theta();
   double phi = particle.momentum().Phi();
 
   // Energy of these photons
   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;
   }
 }

double fastsim::MuonBremsstrahlung::PetrukhinFunc	(	double *	x,
		double *	p
	)

staticprivate

Petrukhin Function: Returns cross section using nuclear-electron screening correction from G4 style.

Definition at line 241 of file MuonBremsstrahlung.cc.

References funct::A, alpha, TtFullHadDaughter::B, CommonMethods::delta(), validate-o2o-wbm::f, funct::pow(), and BeamSpotPI::Z.

                                                                     {
   // Function independent variable
   double nu = x[0];  // fraction of muon's energy transferred to the photon
 
   // Parameters
   double E = p[0];  // Muon Energy (in GeV)
   double A = p[1];  // Atomic weight
   double Z = p[2];  // Atomic number
 
   /*
         // Function of Muom Brem using  nuclear screening correction
         // Ref: http://pdg.lbl.gov/2008/AtomicNuclearProperties/adndt.pdf
         // http://geant4.cern.ch/G4UsersDocuments/UsersGuides/PhysicsReferenceManual/html/node48.html
 
         // Physical constants
         double B = 182.7;
         double ee = sqrt(2.7181) ; // sqrt(e)
         double ZZ=  pow( Z,-1./3.); // Z^-1/3
         double emass = 0.0005109990615;  // electron mass (GeV/c^2)
         double mumass = 0.105658367;//mu mass  (GeV/c^2)
 
         double re = 2.817940285e-13;// Classical electron radius (Units: cm)
         double alpha = 1./137.03599976; // fine structure constant
         double Dn = 1.54* (pow(A,0.27));
         double constant =  pow((2.0 * Z * emass/mumass * re ),2.0);
 
         double delta = (mumass * mumass * nu) /(2.* E * (1.- nu)); 
 
         double Delta_n = TMath::Log(Dn / (1.+ delta *( Dn * ee -2.)/ mumass)); //nuclear screening correction 
 
         double Phi = TMath::Log((B * mumass * ZZ / emass)/ (1.+ delta * ee * B * ZZ  / emass)) - Delta_n;//phi(delta)
 
         // Diff. Cross Section for Muon Brem from a screened nuclear (Equation 16: REF: LBNL-44742)
         double f = alpha * constant *(4./3.-4./3.*nu + nu*nu)*Phi/nu;
     */
 
   // Function for Muon Brem Xsec from G4
 
   // Physical constants
   double B = 183.;
   double Bl = 1429.;
   double ee = 1.64872;            // sqrt(e)
   double Z13 = pow(Z, -1. / 3.);  // Z^-1/3
   double Z23 = pow(Z, -2. / 3.);  // Z^-2/3
 
   // Original values of paper
   double emass = 0.0005109990615;  // electron mass (GeV/c^2)
   double mumass = 0.105658367;     // muon mass  (GeV/c^2)
   // double re = 2.817940285e-13;     // Classical electron radius (Units: cm)
   double alpha = 0.00729735;      // 1./137.03599976; // fine structure constant
   double constant = 1.85736e-30;  // pow( ( emass / mumass * re ) , 2.0);
 
   // Use nomenclature from reference -> Follow those formula step by step
   if (nu * E >= E - mumass)
     return 0;
 
   double Dn = 1.54 * (pow(A, 0.27));
   double Dnl = pow(Dn, (1. - 1. / Z));
 
   double delta = (mumass * mumass * nu) / (2. * E * (1. - nu));
 
   double Phi_n = TMath::Log(B * Z13 * (mumass + delta * (Dnl * ee - 2)) / (Dnl * (emass + delta * ee * B * Z13)));
   if (Phi_n < 0)
     Phi_n = 0;
 
   double Phi_e =
       TMath::Log((Bl * Z23 * mumass) / (1. + delta * mumass / (emass * emass * ee)) / (emass + delta * ee * Bl * Z23));
   if (Phi_e < 0 || nu * E >= E / (1. + (mumass * mumass / (2. * emass * E))))
     Phi_e = 0;
 
   // Diff. Cross Section for Muon Brem from G4 (without NA/A factor)
   double f = 16. / 3. * alpha * constant * Z * (Z * Phi_n + Phi_e) * (1. / nu) * (1. - nu + 0.75 * nu * nu);
 
   return f;
 }