Auxiliary class for the calculation of the muon radiative cross sections. More...

#include <JRadiation.hh>

Inheritance diagram for JPHYSICS::JRadiation:

Public Member Functions
virtual	~JRadiation ()
	Virtual desctructor. More...

	JRadiation (const double z, const double a, const int integrsteps, const double eminBrems, const double eminEErad, const double eminGNrad)
	Constructor. More...

double	SigmaEErad (const double E, const double eps) const
	Pair production cross section. More...

virtual double	TotalCrossSectionEErad (const double E) const
	Pair production cross section. More...

double	TotalCrossSectionBrems (const double E) const
	Bremsstrahlung cross section. More...

virtual double	SigmaGNrad (const double E, const double eps) const
	Photo-nuclear cross section. More...

virtual double	TotalCrossSectionGNrad (const double E) const
	Photo-nuclear cross section. More...

double	EfromBrems (const double E) const
	Bremsstrahlung shower energy. More...

double	ThetaRMSfromBrems (const double E, const double v) const
	Get RMS of scattering angle for Bremsstrahlung. More...

double	EfromEErad (const double E) const
	Pair production shower energy. More...

double	ThetaRMSfromEErad (const double E, const double v) const
	Get RMS of scattering angle for pair production. More...

double	EfromGNrad (const double E) const
	Photo-nuclear shower energy. More...

double	ThetaRMSfromGNrad (const double E, const double v) const
	Get RMS of scattering angle for photo-nuclear shower. More...

virtual double	CalculateACoeff (double E) const
	Ionization a parameter. More...

Protected Member Functions
double	GofZEvrho (const double E, const double v, const double r) const

virtual double	IntegralofG (const double E, const double eps) const

Static Protected Member Functions
static double	sigmaGammaPparam (const double eps)

static double	PhiofEofepsparam (const double E, const double eps)

static double	le ()

static double	r0 ()

static double	Astar ()

static double	B ()

static double	BP ()

Protected Attributes
const double	Z

const double	A

const double	Dn

const double	DnP

const int	steps

const double	EminBrems

const double	EminEErad

const double	EminGNrad

Detailed Description

Auxiliary class for the calculation of the muon radiative cross sections.

Author: P. Kooijman

Muon scattering: A. van Ginneken, Nucl. Instr. and Meth. A251 (1986) 21.

Definition at line 34 of file JRadiation.hh.

Constructor & Destructor Documentation

virtual JPHYSICS::JRadiation::~JRadiation ( )

inlinevirtual

Virtual desctructor.

Definition at line 40 of file JRadiation.hh.

41 {}

JPHYSICS::JRadiation::JRadiation	(	const double	z,
		const double	a,
		const int	integrsteps,
		const double	eminBrems,
		const double	eminEErad,
		const double	eminGNrad
	)

inline

Constructor.

Parameters

z	number of protons
a	number of nucleons
integrsteps	number of integration steps
eminBrems	energy threshold Bremsstrahlung [GeV]
eminEErad	energy threshold pair production [GeV]
eminGNrad	energy threshold photo-production [GeV]

Definition at line 54 of file JRadiation.hh.

                                        :
       Z(z),
       A(a),
       Dn  (1.54*pow(A,0.27)),
       DnP (pow(Dn,(1-1/Z))),
       steps(integrsteps),
       EminBrems(eminBrems),
       EminEErad(eminEErad),
       EminGNrad(eminGNrad)
     {}

Member Function Documentation

double JPHYSICS::JRadiation::SigmaEErad	(	const double	E,
		const double	eps
	)		const

inline

Pair production cross section.

Parameters

E	muon energy [GeV]
eps	shower energy [GeV]

Returns: cross section [m^2/g]

Definition at line 78 of file JRadiation.hh.

     {
       //routine to calculate dsigma/de in cm^2/GeV for radiating a photon of energy eps
   
       if(eps<4.*MASS_ELECTRON)return 0.;
       if(eps>E-0.75*exp(1.0)*pow(Z,1./3.)) return 0.;
       double zeta;
       if(E<35.*MASS_MUON)
         {
           zeta = 0.;
         }
       else
         {
           zeta = (0.073*log((E/MASS_MUON)/(1.+1.95e-5*pow(Z,2./3.)*E/MASS_MUON))-0.26);
           if(zeta<0.)
             {
               zeta=0.;
             }
           else
             {
               zeta = zeta/(0.058*log((E/MASS_MUON)/(1.+5.3e-5*pow(Z,1./3.)*E/MASS_MUON))-0.14);
             }
         }
       double integ = IntegralofG(E,eps);
       //cout<< eps<<" integ "<< integ<<endl;
       return integ*(4/(3.*M_PI))*(Z*(Z+zeta)/A)*AVOGADRO*pow((ALPHA_ELECTRO_MAGNETIC*r0()),2.)*(1-eps/E)/eps;
     }

virtual double JPHYSICS::JRadiation::TotalCrossSectionEErad ( const double E ) const

inlinevirtual

Pair production cross section.

Parameters

E	muon energy [GeV]

Returns: cross section [m^2/g]

Reimplemented in JPHYSICS::JRadiationFunction.

Definition at line 114 of file JRadiation.hh.

     {
       //cout<<" entering TotcrsecEErad"<<endl;
       // radiation goes approx as 1/E so sample in log E
       // calculate the total cross section as S(dsigma/de)*de
       // start at E and go backwards for precision
       double sum(0.);
       double dle = log(E/EminEErad)/steps;
       //cout << "dle "<< dle<<endl;
       for(int i=0;i<steps+1;i++)
         {
           double factor = 1.0;
           if(i==0 || i==steps)factor=0.5;
           //cout<<"i"<<i<<endl;
           double eps = EminEErad*exp(i*dle);sum += factor*eps*SigmaEErad( E, eps);
           //cout<<"eps "<<eps<<" sum "<<sum<<endl;
         }
       return sum*dle;
     }

double JPHYSICS::JRadiation::TotalCrossSectionBrems ( const double E ) const

inline

Bremsstrahlung cross section.

Parameters

E	muon energy [GeV]

Returns: cross section [m^2/g]

Definition at line 142 of file JRadiation.hh.

     {
       //  double delta = (MASS_MUON*MASS_MUON*eps)/(2*E*(E-eps));
       double delta = (MASS_MUON*MASS_MUON)/(2*E);
       //double v = eps/E;
       double Phin = log((B()*pow(Z,-1./3.)*(MASS_MUON+delta*(DnP*sqrt(exp(1.0))-2.)))/(DnP*(MASS_ELECTRON+delta*sqrt(exp(1.0))*B()*pow(Z,-1./3.))));
       if(Phin<0.)Phin=0.;
       double Phie = log((BP()*pow(Z,-2./3.)*MASS_MUON)/((1+delta*MASS_MUON/(MASS_ELECTRON*MASS_ELECTRON*sqrt(exp(1.0))))*(MASS_ELECTRON+delta*sqrt(exp(1.0))*BP()*pow(Z,-2./3.))));
       if(Phie<0.)Phie=0.;
       //if(eps>E/(1+MASS_MUON*MASS_MUON/(2.*MASS_ELECTRON*E)))Phie=0.;
       double sig = (16./3.)*ALPHA_ELECTRO_MAGNETIC*AVOGADRO*pow((MASS_ELECTRON/MASS_MUON*r0()),2.0)*Z*(Z*Phin+Phie)/(A);
       double epsint = log((E-MASS_MUON)/EminBrems)-(E-MASS_MUON-EminBrems)/E+0.375*(pow(E-MASS_MUON,2.)-pow(EminBrems,2.0))/pow(E,2.0);
       return epsint*sig;//"cross section" in m^2/g multiplied by density and inverted gives mean free path 
     }

virtual double JPHYSICS::JRadiation::SigmaGNrad	(	const double	E,
		const double	eps
	)		const

inlinevirtual

Photo-nuclear cross section.

Parameters

E	muon energy [GeV]
eps	shower energy [GeV]

Returns: cross section [m^2/g]

Definition at line 165 of file JRadiation.hh.

     {
      
       //minimum radiated energy is 0.2 GeV, not very critical formulae become invalid for very much lower
       //result is given in m^2/g GeV
       if (eps<0.2) return 0.;
       if (eps>E-MASS_PROTON) return 0.;
       double Aeff;
       if (A==1)
         {Aeff = 1.;}
       else
         {Aeff = (0.22*A + 0.78 * pow(A,0.89));}
       double sigmaGammaPofeps = (49.2 + 11.1 * log(eps) + 151.8/sqrt(eps))*pow(10,-34.);
       double Psiofeps = ALPHA_ELECTRO_MAGNETIC/PI*Aeff*AVOGADRO/A*sigmaGammaPofeps/eps;
       double Denom = 1+eps/LAMBDA*(1+LAMBDA/(2*MASS_PROTON)+eps/LAMBDA);
       double epsoverE = eps/E;
       double Numerator = E * E * (1 - epsoverE) / (MASS_MUON * MASS_MUON) * (1 + MASS_MUON * MASS_MUON * epsoverE * epsoverE / (LAMBDA * LAMBDA * (1 - epsoverE)));
       double Factor = 1 - epsoverE + epsoverE * epsoverE / 2 * (1 + 2 * MASS_MUON * MASS_MUON / (LAMBDA * LAMBDA));
       double PhiofEofeps = epsoverE - 1 + Factor * log (Numerator / Denom);
       //cout<<"PhiofEofeps "<<PhiofEofeps<<" Psiofeps "<<Psiofeps<<endl;
       return Psiofeps*PhiofEofeps;
     }

virtual double JPHYSICS::JRadiation::TotalCrossSectionGNrad ( const double E ) const

inlinevirtual

Photo-nuclear cross section.

Parameters

E	muon energy [GeV]

Returns: cross section [m^2/g]

Reimplemented in JPHYSICS::JRadiationFunction.

Definition at line 195 of file JRadiation.hh.

     {
       double epsmin = 0.2;
       double epsmax = E-MASS_PROTON/2;
       double sum(0.);
       double dle = log(epsmax/epsmin)/steps;
       for(int i=0;i<steps+1;i++)
         {
           float factor = 1.0;
           if(i==0 || i==steps)factor=0.5;
           //cout<<"i"<<i<<endl;
           double eps = epsmin*exp(i*dle);sum += factor*eps*SigmaGNrad( E, eps);
           //cout<<"eps "<<eps<<" sum "<<sum<<endl;
         }
 
       return sum*dle;
     }

double JPHYSICS::JRadiation::EfromBrems ( const double E ) const

inline

Bremsstrahlung shower energy.

Parameters

E	muon energy [GeV]

Returns: shower energy [GeV]

Definition at line 220 of file JRadiation.hh.

     {
       //double EminBrems(0.01);
       //generate according to 1/k from minimum energy to E
       double Er = 0.0;
       for (int i = 1000; i != 0; --i) {
         Er = EminBrems*exp(gRandom->Rndm()*log(E/EminBrems));
         //check on the extra factor (1-v+3/4v^2)
         if(gRandom->Rndm()<(1.-Er/E+0.75*pow(Er/E,2))) break;
       }
       return Er;
     }

double JPHYSICS::JRadiation::ThetaRMSfromBrems	(	const double	E,
		const double	v
	)		const

inline

Get RMS of scattering angle for Bremsstrahlung.

Parameters

E	muon energy [GeV]
v	energy loss fraction

Returns: angle [rad]

Definition at line 241 of file JRadiation.hh.

     {
       using namespace std;
       
       const double precision = 1.0e-3;
 
       const double k1 = 0.092 * pow(E, -1.0/3.0);
       const double k2 = 0.052 * pow(E, -1.0) * pow(Z, -1.0/4.0);
       const double k3 = 0.220 * pow(E, -0.92);
       const double k4 = 0.260 * pow(E, -0.91);
 
       double rms = 0.0;
     
       if (v <= 0.5) {
 
         rms = max(min(k1*sqrt(v), k2), k3*v);
 
       } else {
 
         const double n  = 0.81 * sqrt(E) / (sqrt(E) + 1.8);
 
         double k5 = 0.0;
 
         for (double vmin = 0.5, vmax = 1.0; ; ) {
 
           const double v = 0.5 * (vmin + vmax);
 
           const double y = k4 * pow(v, 1.0 + n) * pow(1.0 - v, -n);
 
           if (abs(y - 0.2) <= precision) {
 
             k5 = y * pow(1.0 - v, 1.0/2.0);
 
             break;
           }
 
           if (y < 0.2)
             vmin = v;
           else
             vmax = v;
         }
 
         rms = k4 * pow(v, 1.0 + n) * pow(1.0 - v, -n);
 
         if (rms >= 0.2) {
           rms = k5 * pow(1.0 - v, -1.0/2.0);
         }
       }
 
       return rms;
     }

double JPHYSICS::JRadiation::EfromEErad ( const double E ) const

inline

Pair production shower energy.

Parameters

E	muon energy [GeV]

Returns: shower energy [GeV]

Definition at line 300 of file JRadiation.hh.

     {
       //generate according to 1/k from minimum energy to E
 
       const double eps =0.2;
       const double IntGmax = IntegralofG(E,eps)*2.0;
 
       double Er = 0.0;
       for (int i = 1000; i != 0; --i)
         {
           Er = EminEErad*exp(gRandom->Rndm()*log(E/EminEErad));
           //check on the extra factor, (1-v) and IntofG
           double factor = (1.-Er/E)*IntegralofG(E,Er)/IntGmax;
           if(gRandom->Rndm()<factor) break;
         }
       return Er;
     }

double JPHYSICS::JRadiation::ThetaRMSfromEErad	(	const double	E,
		const double	v
	)		const

inline

Get RMS of scattering angle for pair production.

Parameters

E	muon energy [GeV]
v	energy loss fraction

Returns: angle [rad]

Definition at line 326 of file JRadiation.hh.

     {
       using namespace std;
 
       const double a = 8.9e-4;
       const double b = 1.5e-5;
       const double c = 0.032;
       const double d = 1.0;
       const double e = 0.1;
 
       const double n = -1.0;
 
       const double u = v - 2.0*MASS_ELECTRON/E;
 
       if (u > 0.0) 
         return (2.3 + log(E)) * (1.0/E) * pow(1.0 - v, n) * (u*u) * (1.0/(v*v)) * 
           min(a * pow(v, 1.0/4.0) * (1.0 + b*E) + c*v/(v+d), e);
       else
         return 0.0;
     }

double JPHYSICS::JRadiation::EfromGNrad ( const double E ) const

inline

Photo-nuclear shower energy.

Parameters

E	muon energy [GeV]

Returns: shower energy [GeV]

Definition at line 354 of file JRadiation.hh.

     {
       //generate according to 1/k from minimum energy to E
 
       double cmax = sigmaGammaPparam(EminGNrad); 
       if (cmax < sigmaGammaPparam(E))
         cmax=sigmaGammaPparam(E);
 
       double Pmax = PhiofEofepsparam(E,EminGNrad);
 
       double Er = 0.0;
       for (int i = 1000; i != 0; --i) {
         Er = EminGNrad*exp(gRandom->Rndm()*log(E/EminGNrad));
         const double factor = PhiofEofepsparam(E,Er)*sigmaGammaPparam(Er)/(cmax*Pmax);
         if (gRandom->Rndm() < factor) return Er;
       }
 
       return Er;
     }

double JPHYSICS::JRadiation::ThetaRMSfromGNrad	(	const double	E,
		const double	v
	)		const

inline

Get RMS of scattering angle for photo-nuclear shower.

Parameters

E	muon energy [GeV]
v	energy loss fraction

Returns: angle [rad]

Definition at line 382 of file JRadiation.hh.

     {
       return 0.0;
     }

virtual double JPHYSICS::JRadiation::CalculateACoeff ( double E ) const

inlinevirtual

Ionization a parameter.

Parameters

E	muon energy [GeV]

Returns: ionization coefficient [GeV*m^2/g]

Reimplemented in JPHYSICS::JRadiationFunction.

Definition at line 394 of file JRadiation.hh.

     {
       const double E2      = E*E;
       const double beta    = sqrt(E2 - MASS_MUON*MASS_MUON) / E;
       const double beta2   = beta*beta;
       const double gamma   = E/MASS_MUON;
       const double gamma2  = gamma*gamma;
 
       const double p2      = E2 - MASS_MUON*MASS_MUON;
       const double EMaxT   = 2.*MASS_ELECTRON*p2 / (MASS_ELECTRON*MASS_ELECTRON + MASS_MUON*MASS_MUON + 2.*MASS_ELECTRON*E);
       const double EMaxT2  = EMaxT*EMaxT;
       const JSter& coeff   = getSterCoefficient(Z);
       const double I2      = coeff.I*coeff.I;           // in GeV^2
       const double X       = log10(beta*gamma);
 
       double delta = 0.0;
 
       if (coeff.X0 < X && X < coeff.X1) {
         delta = 4.6052*X + coeff.a*pow(coeff.X1-X,coeff.m) + coeff.C;
       }
 
       if (X > coeff.X1) {
         delta = 4.6052*X + coeff.C;
       }
 
       double acoeff = (2*PI*AVOGADRO*ALPHA_ELECTRO_MAGNETIC*ALPHA_ELECTRO_MAGNETIC*le()*le())*Z/A*(MASS_ELECTRON/beta2)*(log(2.*MASS_ELECTRON*beta2*gamma2*EMaxT/I2)-2.*beta2+0.25*(EMaxT2/E2)-delta);
 
       return acoeff;
     }

double JPHYSICS::JRadiation::GofZEvrho	(	const double	E,
		const double	v,
		const double	r
	)		const

inlineprotected

Definition at line 425 of file JRadiation.hh.

     {
       const double ksi = MASS_MUON*MASS_MUON*v*v*(1-r*r)/(4*MASS_ELECTRON*MASS_ELECTRON*(1-v));//
       const double b   = v*v/(2*(1-v));//
       const double Be  = ((2+r*r)*(1+b)+ksi*(3+r*r))*log(1+1/ksi)+(1-r*r-b)/(1+ksi)-(3+r*r);
       const double Bm  = ((1+r*r)*(1+3*b/2)-(1+2*b)*(1-r*r)/ksi)*log(1+ksi)+ksi*(1-r*r-b)/(1+ksi)+(1+2*b)*(1-r*r);
       const double Ye  = (5-r*r+4*b*(1+r*r))/(2*(1+3*b)*log(3+1/ksi)-r*r-2*b*(2-r*r));
       const double Ym  = (4+r*r+3*b*(1+r*r))/((1+r*r)*(1.5+2*b)*log(3+ksi)+1-1.5*r*r);
       const double Le  = log((Astar()*pow(Z,-1./3.)*sqrt((1+ksi)*(1+Ye)))/
                              (1.+(2.*MASS_ELECTRON*sqrt(exp(1.))*Astar()*pow(Z,-1./3.)*(1+ksi)*(1+Ye))/(E*v*(1-r*r))));
       const double Lm  = log(((MASS_MUON/MASS_ELECTRON)*Astar()*pow(Z,-1./3.)*sqrt((1.+1./ksi)*(1.+Ym)))/
                              (1.+(2.*MASS_ELECTRON*sqrt(exp(1.))*Astar()*pow(Z,-1./3.)*(1+ksi)*(1+Ym))/(E*v*(1-r*r))));
       double Phie = Be*Le; if(Phie<0.)Phie=0.;
       double Phim = Bm*Lm; if(Phim<0.)Phim=0.;
       return Phie+(MASS_ELECTRON*MASS_ELECTRON/(MASS_MUON*MASS_MUON))*Phim;
     }

virtual double JPHYSICS::JRadiation::IntegralofG	(	const double	E,
		const double	eps
	)		const

inlineprotectedvirtual

Reimplemented in JPHYSICS::JRadiationFunction.

Definition at line 444 of file JRadiation.hh.

     {
       const double EP   = E-eps;
       const double v    = eps/E;
       const double tmin = log((4*MASS_ELECTRON/eps+12*pow(MASS_MUON,2)/(E*EP)*(1-4*MASS_ELECTRON/eps))
                               /(1+(1-6*pow(MASS_MUON,2)/(E*EP))*sqrt(1-4*MASS_ELECTRON/eps)));
 
       if (tmin < 0.0) {
         const double dt   = -tmin/steps;
         double sum(0.);
         for (int i = 0;i<steps+1;i++)
           {
             double fac = 1.0;if(i==0 || i==steps)fac=0.5;
             double t = tmin+i*dt;
             double r = 1.-exp(t);
             sum+=fac*GofZEvrho(E,v,r)*exp(t);
           } 
         return sum*dt;
       } else
         return 0.0;
     }

static double JPHYSICS::JRadiation::sigmaGammaPparam ( const double eps )

inlinestaticprotected

Definition at line 467 of file JRadiation.hh.

     { 
       return (49.2 + 11.1 * log(eps) + 151.8/sqrt(eps));
     }

static double JPHYSICS::JRadiation::PhiofEofepsparam	(	const double	E,
		const double	eps
	)

inlinestaticprotected

Definition at line 472 of file JRadiation.hh.

     {
       const double Denom = 1+eps/LAMBDA*(1+LAMBDA/(2*MASS_PROTON)+eps/LAMBDA);
       const double epsoverE = eps/E;
       const double Numerator = E * E * (1 - epsoverE) / (MASS_MUON * MASS_MUON) * (1 + MASS_MUON * MASS_MUON * epsoverE * epsoverE /
                                                                                    (LAMBDA * LAMBDA * (1 - epsoverE)));
       const double Factor = 1 - epsoverE + epsoverE * epsoverE / 2 * (1 + 2 * MASS_MUON * MASS_MUON / (LAMBDA * LAMBDA));
       return (epsoverE - 1 + Factor * log (Numerator / Denom));
     }

static double JPHYSICS::JRadiation::le ( )

inlinestaticprotected

Definition at line 483 of file JRadiation.hh.

483 { return 3.8616E-13;}

static double JPHYSICS::JRadiation::r0 ( )

inlinestaticprotected

Definition at line 484 of file JRadiation.hh.

484 { return 2.817940e-15; }

static double JPHYSICS::JRadiation::Astar ( )

inlinestaticprotected

Definition at line 485 of file JRadiation.hh.

485 { return 183.0; }

static double JPHYSICS::JRadiation::B ( )

inlinestaticprotected

Definition at line 486 of file JRadiation.hh.

486 { return 183.0; }

static double JPHYSICS::JRadiation::BP ( )

inlinestaticprotected

Definition at line 487 of file JRadiation.hh.

487 { return 1429.0; }

Member Data Documentation

const double JPHYSICS::JRadiation::Z

protected

Definition at line 489 of file JRadiation.hh.

const double JPHYSICS::JRadiation::A

protected

Definition at line 490 of file JRadiation.hh.

const double JPHYSICS::JRadiation::Dn

protected

Definition at line 491 of file JRadiation.hh.

const double JPHYSICS::JRadiation::DnP

protected

Definition at line 492 of file JRadiation.hh.

const int JPHYSICS::JRadiation::steps

protected

Definition at line 493 of file JRadiation.hh.

const double JPHYSICS::JRadiation::EminBrems

protected

Definition at line 494 of file JRadiation.hh.

const double JPHYSICS::JRadiation::EminEErad

protected

Definition at line 495 of file JRadiation.hh.

const double JPHYSICS::JRadiation::EminGNrad

protected

Definition at line 496 of file JRadiation.hh.

The documentation for this class was generated from the following file:

JRadiation.hh

Public Member Functions

Protected Member Functions

Static Protected Member Functions

Protected Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation