Fast implementation of class JRadiation. More...

#include <JRadiationSource.hh>

Inheritance diagram for JPHYSICS::JRadiationFunction:

Public Member Functions
	JRadiationFunction (const JRadiation &radiation, const unsigned int number_of_bins, const double Emin, const double Emax)
	Constructor. More...

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

virtual double	CalculateACoeff (const double E) const override
	A Coefficient. More...

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

double	SigmaEErad (const double E, const double eps) 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...

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

double	getThetaRMSfromBrems (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	getThetaRMSfromEErad (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	getThetaRMSfromGNrad (const double E, const double v) const
	Get RMS of scattering angle for photo-nuclear shower. More...

Protected Member Functions
virtual double	IntegralofG (const double E, const double eps) const

double	GofZEvrho (const double E, const double v, const double r) 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
JFunction1D_t	sigmaEE

JFunction1D_t	sigmaGN

JFunction1D_t	Acoeff

JFunction2D_t	integral

const double	Z

const double	A

const double	Dn

const double	DnP

const int	steps

const double	EminBrems

const double	EminEErad

const double	EminGNrad

Private Types
typedef JTOOLS::JGridPolint1Function1D_t	JFunction1D_t

typedef JTOOLS::JMultiFunction < JFunction1D_t, JTOOLS::JMapList < JTOOLS::JPolint1FunctionalGridMap > >	JFunction2D_t

Detailed Description

Fast implementation of class JRadiation.

In this, the methods

JRadiation::TotalCrossSectionEErad(...), JRadiation::TotalCrossSectionGNrad(...) and JRadiation::IntegralofG(...)

are replaced by lookup tables.

Definition at line 35 of file JRadiationSource.hh.

Member Typedef Documentation

typedef JTOOLS::JGridPolint1Function1D_t JPHYSICS::JRadiationFunction::JFunction1D_t

private

Definition at line 39 of file JRadiationSource.hh.

typedef JTOOLS::JMultiFunction<JFunction1D_t, JTOOLS::JMapList<JTOOLS::JPolint1FunctionalGridMap> > JPHYSICS::JRadiationFunction::JFunction2D_t

private

Definition at line 41 of file JRadiationSource.hh.

Constructor & Destructor Documentation

JPHYSICS::JRadiationFunction::JRadiationFunction	(	const JRadiation &	radiation,
		const unsigned int	number_of_bins,
		const double	Emin,
		const double	Emax
	)

inlineexplicit

Constructor.

Parameters

radiation	JRadiation object
number_of_bins	number of bins
Emin	minimal muon energy [GeV]
Emax	maximal muon energy [GeV]

Acoeff

Definition at line 52 of file JRadiationSource.hh.

                                                          :
       JRadiation(radiation)
     {
       using namespace JTOOLS;
 
       const double xmin = log(Emin);
       const double xmax = log(Emax);
 
     
       integral.configure(make_grid(number_of_bins, xmin, xmax));
 
       for (JFunction2D_t::iterator i = integral.begin(); i != integral.end(); ++i) {
 
         const double x    = i->getX();
         const double E    = exp(x);
 
         const double ymin = log(EminEErad/E);
         const double ymax = 0.0;
 
         if (ymax > ymin) {
 
           i->getY().configure(make_grid(number_of_bins, ymin, ymax));
           
           for (JFunction1D_t::iterator j = i->getY().begin(); j != i->getY().end(); ++j) {
             
             const double y   = j->getX();
             const double eps = exp(y) * E;
             const double z   = JRadiation::IntegralofG(E, eps);
 
             j->getY() = z;
           }
         }
       }
     
       integral.compile();
       integral.setExceptionHandler(new JFunction1D_t::JDefaultResult(0.0));
 
 
       sigmaEE.configure(make_grid(number_of_bins, xmin, xmax));
 
       for (JFunction1D_t::iterator i = sigmaEE.begin(); i != sigmaEE.end(); ++i) {
 
         const double E = exp(i->getX());
         const double y = JRadiation::TotalCrossSectionEErad(E);
 
         i->getY() = y;
       }
                     
       sigmaEE.compile();
 
 
       sigmaGN.configure(make_grid(number_of_bins, xmin, xmax));
 
       for (JFunction1D_t::iterator i = sigmaGN.begin(); i != sigmaGN.end(); ++i) {
 
         const double E = exp(i->getX());
         const double y = JRadiation::TotalCrossSectionGNrad(E);
 
         i->getY() = y;
       }
                     
       sigmaGN.compile();
         
     /**
        Acoeff*/
      Acoeff.configure(make_grid(number_of_bins, xmin, xmax));
      for (JFunction1D_t::iterator i = Acoeff.begin(); i != Acoeff.end(); ++i) {
      
         const double E = exp(i->getX());
         const double y = JRadiation::CalculateACoeff(E);
 
         i->getY() = y;                                                                                                                                        
     }
     
   Acoeff.compile();                                                                                                                                    
 
   }

Member Function Documentation

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

inlineoverridevirtual

Pair production cross section.

Parameters

E	muon energy [GeV]

Returns: cross section [m^2/g]

Reimplemented from JPHYSICS::JRadiation.

Definition at line 140 of file JRadiationSource.hh.

     {
       const double x = log(E);
 
       if (x >= sigmaEE. begin()->getX() &&
           x <= sigmaEE.rbegin()->getX()) {
 
         try {
           return sigmaEE(x);
         }
         catch(std::exception& error) {
           std::cerr << "JRadiation::TotalCrossSectionEErad() " << E << ' ' << error.what() << std::endl;
         }
       }
 
       return JRadiation::TotalCrossSectionEErad(E);
     }

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

inlineoverridevirtual

A Coefficient.

Parameters

E	muon energy [GeV]

Returns: Ionization A parameter [GeV/m]

Reimplemented from JPHYSICS::JRadiation.

Definition at line 163 of file JRadiationSource.hh.

   {
     const double x = log(E);
 
     if (x >= Acoeff. begin()->getX() &&
         x <= Acoeff.rbegin()->getX()) {
 
       try {
         return Acoeff(x);
       }
       catch(std::exception& error) {
         std::cerr << "JRadiation::ACoeff() " << E << ' ' << error.what() << std::endl;
       }
     }
     return JRadiation::CalculateACoeff(E);
   }

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

inlineoverridevirtual

Photo-nuclear cross section.

Parameters

E	muon energy [GeV]

Returns: cross section [m^2/g]

Reimplemented from JPHYSICS::JRadiation.

Definition at line 187 of file JRadiationSource.hh.

     {
       const double x = log(E);
 
       if (x >= sigmaGN. begin()->getX() &&
           x <= sigmaGN.rbegin()->getX()) {
 
         try {
           return sigmaGN(x);
         }
         catch(std::exception& error) {
           std::cerr << "JRadiation::TotalCrossSectionGNrad() " << E << ' ' << error.what() << std::endl;
         }
       }
 
       return JRadiation::TotalCrossSectionGNrad(E);
     }

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

inlineprotectedvirtual

Reimplemented from JPHYSICS::JRadiation.

Definition at line 206 of file JRadiationSource.hh.

     {
       try {
         
         const double x = log(E); 
         const double y = log(eps/E); 
         const double z = integral(x, y);
 
         return z;
       }
       catch(std::exception& error) {
         std::cerr << "JRadiation::IntegralofG() " << E << ' ' << eps << ' ' << error.what() << std::endl;
       }
 
       return JRadiation::IntegralofG(E, eps);
     }

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

inlineinherited

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

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

inlineinherited

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

inlinevirtualinherited

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

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

inlineinherited

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::getThetaRMSfromBrems	(	const double	E,
		const double	v
	)		const

inlineinherited

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

inlineinherited

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::getThetaRMSfromEErad	(	const double	E,
		const double	v
	)		const

inlineinherited

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

inlineinherited

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::getThetaRMSfromGNrad	(	const double	E,
		const double	v
	)		const

inlineinherited

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

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

inlineprotectedinherited

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

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