JPHYSICS::JRadiationFunction Class Reference

Fast implementation of class JRadiation. More...

#include <JRadiationFunction.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.
virtual double	TotalCrossSectionEErad (const double E) const override
	Pair production cross section.
virtual double	TotalCrossSectionGNrad (const double E) const override
	Photo-nuclear cross section.
virtual double	CalculateACoeff (const double E) const override
	Ionization a parameter.
double	SigmaEErad (const double E, const double eps) const
	Pair production cross section.
double	TotalCrossSectionBrems (const double E) const
	Bremsstrahlung cross section.
virtual double	SigmaGNrad (const double E, const double eps) const
	Photo-nuclear cross section.
double	EfromBrems (const double E) const
	Bremsstrahlung shower energy.
double	ThetaRMSfromBrems (const double E, const double v) const
	Get RMS of scattering angle for Bremsstrahlung.
double	EfromEErad (const double E) const
	Pair production shower energy.
double	ThetaRMSfromEErad (const double E, const double v) const
	Get RMS of scattering angle for pair production.
double	EfromGNrad (const double E) const
	Photo-nuclear shower energy.
double	ThetaRMSfromGNrad (const double E, const double v) const
	Get RMS of scattering angle for photo-nuclear shower.

Protected Member Functions
virtual double	IntegralofG (const double E, const double eps) const override
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,
• JRadiation::CalculateACoeff and
• JRadiation::IntegralofG

are reimplemented using lookup tables.

Definition at line 32 of file JRadiationFunction.hh.

Member Typedef Documentation

JFunction1D_t

JTOOLS::JGridPolint1Function1D_t JPHYSICS::JRadiationFunction::JFunction1D_t

private

Definition at line 35 of file JRadiationFunction.hh.

JFunction2D_t

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

private

Definition at line 37 of file JRadiationFunction.hh.

Constructor & Destructor Documentation

JRadiationFunction()

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

inline

Constructor.

Parameters

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

Definition at line 48 of file JRadiationFunction.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.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

TotalCrossSectionEErad()

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 135 of file JRadiationFunction.hh.

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

TotalCrossSectionGNrad()

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 158 of file JRadiationFunction.hh.

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

CalculateACoeff()

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

inlineoverridevirtual

Ionization a parameter.

Parameters

E	muon energy [GeV]

Returns

ionization coefficient [GeV*m^2/g]

Reimplemented from JPHYSICS::JRadiation.

Definition at line 181 of file JRadiationFunction.hh.

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

IntegralofG()

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

inlineoverrideprotectedvirtual

Reimplemented from JPHYSICS::JRadiation.

Definition at line 198 of file JRadiationFunction.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) {}
 
      return JRadiation::IntegralofG(E, eps);
    }

SigmaEErad()

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

TotalCrossSectionBrems()

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

SigmaGNrad()

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

EfromBrems()

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

inlineinherited

Bremsstrahlung shower energy.

Parameters

E	muon energy [GeV]

Returns

shower energy [GeV]

Definition at line 222 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;
    }

ThetaRMSfromBrems()

double JPHYSICS::JRadiation::ThetaRMSfromBrems ( 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 243 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;
    }

EfromEErad()

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

ThetaRMSfromEErad()

double JPHYSICS::JRadiation::ThetaRMSfromEErad ( 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 328 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;
    }

EfromGNrad()

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

ThetaRMSfromGNrad()

double JPHYSICS::JRadiation::ThetaRMSfromGNrad ( 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 384 of file JRadiation.hh.

    {
      return 0.0;
    }

GofZEvrho()

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

inlineprotectedinherited

Definition at line 427 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;
    }

sigmaGammaPparam()

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

inlinestaticprotectedinherited

Definition at line 469 of file JRadiation.hh.

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

PhiofEofepsparam()

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

inlinestaticprotectedinherited

Definition at line 474 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));
    }

le()

static double JPHYSICS::JRadiation::le ( )

inlinestaticprotectedinherited

Definition at line 485 of file JRadiation.hh.

{ return 3.8616E-13;}

r0()

static double JPHYSICS::JRadiation::r0 ( )

inlinestaticprotectedinherited

Definition at line 486 of file JRadiation.hh.

{ return 2.817940e-15; }

Astar()

static double JPHYSICS::JRadiation::Astar ( )

inlinestaticprotectedinherited

Definition at line 487 of file JRadiation.hh.

{ return 183.0; }

B()

static double JPHYSICS::JRadiation::B ( )

inlinestaticprotectedinherited

Definition at line 488 of file JRadiation.hh.

{ return 183.0; }

BP()

static double JPHYSICS::JRadiation::BP ( )

inlinestaticprotectedinherited

Definition at line 489 of file JRadiation.hh.

{ return 1429.0; }

Member Data Documentation

sigmaEE

JFunction1D_t JPHYSICS::JRadiationFunction::sigmaEE

protected

Definition at line 214 of file JRadiationFunction.hh.

sigmaGN

JFunction1D_t JPHYSICS::JRadiationFunction::sigmaGN

protected

Definition at line 215 of file JRadiationFunction.hh.

Acoeff

JFunction1D_t JPHYSICS::JRadiationFunction::Acoeff

protected

Definition at line 216 of file JRadiationFunction.hh.

integral

JFunction2D_t JPHYSICS::JRadiationFunction::integral

protected

Definition at line 217 of file JRadiationFunction.hh.

Z

const double JPHYSICS::JRadiation::Z

protectedinherited

Definition at line 491 of file JRadiation.hh.

A

const double JPHYSICS::JRadiation::A

protectedinherited

Definition at line 492 of file JRadiation.hh.

Dn

const double JPHYSICS::JRadiation::Dn

protectedinherited

Definition at line 493 of file JRadiation.hh.

DnP

const double JPHYSICS::JRadiation::DnP

protectedinherited

Definition at line 494 of file JRadiation.hh.

steps

const int JPHYSICS::JRadiation::steps

protectedinherited

Definition at line 495 of file JRadiation.hh.

EminBrems

const double JPHYSICS::JRadiation::EminBrems

protectedinherited

Definition at line 496 of file JRadiation.hh.

EminEErad

const double JPHYSICS::JRadiation::EminEErad

protectedinherited

Definition at line 497 of file JRadiation.hh.

EminGNrad

const double JPHYSICS::JRadiation::EminGNrad

protectedinherited

Definition at line 498 of file JRadiation.hh.

---

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

• JRadiationFunction.hh