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

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 42 of file MuonBremsstrahlung.cc.

Constructor & Destructor Documentation

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

Constructor.

Definition at line 98 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 49 of file MuonBremsstrahlung.cc.

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

49 {;};

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 209 of file MuonBremsstrahlung.cc.

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

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

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

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

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

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

References A_, funct::abs(), brem(), density_, SoftLeptonByDistance_cfi::distance, fastsim::SimplifiedGeometry::getThickness(), mps_fire::i, SiStripPI::max, minPhotonEnergy_, minPhotonEnergyFraction_, fastsim::Particle::momentum(), fastsim::Constants::NA, fastsim::Particle::pdgId(), Petrfunc, PetrukhinFunc(), RandomEngineAndDistribution::poissonShoot(), fastsim::Particle::position(), radLenInCm_, theta(), TrackerOfflineValidation_Dqm_cff::xmax, TrackerOfflineValidation_Dqm_cff::xmin, and Z_.

Referenced by ~MuonBremsstrahlung().

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

References patCaloMETCorrections_cff::A, alpha, TtFullHadDaughter::B, DEFINE_EDM_PLUGIN, delta, f, funct::pow(), and DOFs::Z.

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

                                                             {
  
     // 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;
 }