JRadiation.hh

Go to the documentation of this file.

#ifndef __JPHYSICS__JRADIATION__
#define __JPHYSICS__JRADIATION__
 
#include <cmath>
 
#include <TRandom3.h>
 
#include "JLang/JCC.hh"
#include "JPhysics/JConstants.hh"
#include "JPhysics/JIonization.hh"
 
/**
 * \file
 * Muon radiative cross sections.
 * \author P. Kooijman
 */
 
namespace JPHYSICS {}
namespace JPP { using namespace JPHYSICS; }
 
namespace JPHYSICS {
 
  namespace {
    static const double LAMBDA = 0.65;       // [GeV]
  }
 
  /**
   * Auxiliary class for the calculation of the muon radiative cross sections.
   * \author P. Kooijman
   *
   * Bremsstrahlung and pair production cross sections are taken from <a href="https://pdg.lbl.gov/index.html">Particle Data Group</a>. 
   *
   * Gamma-nucleon cross section:
   * Phys.\ Lett.\ B\ Vol.\ 697-3 (2011) 184-193.
   */
  class JRadiation {
  public:
    
    /**
     * Virtual desctructor.
     */
    virtual ~JRadiation()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  z            number of protons
     * \param  a            number of nucleons
     * \param  integrsteps  number of integration steps
     * \param  eminBrems    energy threshold Bremsstrahlung  [GeV]
     * \param  eminEErad    energy threshold pair  production [GeV]
     * \param  eminGNrad    energy threshold photo-production [GeV]
     */
    JRadiation(const double z, 
               const double a,
               const int    integrsteps,
               const double eminBrems,
               const double eminEErad,
               const double eminGNrad) :
      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)
    {}
 
 
    /**
     * Pair production cross section.
     *
     * \param  E     muon   energy [GeV]
     * \param  eps   shower energy [GeV]
     * \return       cross section [m^2/g]
     */
    double SigmaEErad(const double E,
                      const double eps) const
    {
      //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;
    }
 
 
    /**
     * Pair production cross section.
     *
     * \param  E     muon    energy [GeV]
     * \return       cross section [m^2/g]
     */
    virtual double TotalCrossSectionEErad(const double E) const
    {
      //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;
    }
 
 
    /**
     * Bremsstrahlung cross section.
     *
     * \param  E     muon    energy [GeV]
     * \return       cross section [m^2/g]
     */
    
    double TotalCrossSectionBrems(const double E) const
    {
      //  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 
    }
 
 
    /**
     * Photo-nuclear cross section.
     *
     * \param  E     muon   energy [GeV]
     * \param  eps   shower energy [GeV]
     * \return       cross section [m^2/g]
     */
    virtual double SigmaGNrad(const double E, const double eps) const
    {
     
      //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;
    }
 
 
    /**
     * Photo-nuclear cross section.
     *
     * \param  E     muon   energy [GeV]
     * \return       cross section [m^2/g]
     */
    virtual double TotalCrossSectionGNrad(const double E) const
    {
      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;
    }
 
 
    /**
     * Bremsstrahlung shower energy.
     *
     * \param  E     muon   energy [GeV]
     * \return       shower energy [GeV]
     */
    double EfromBrems(const double E) const
    {
      //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;
    }
 
 
    /**
     * Get RMS of scattering angle for Bremsstrahlung.
     *
     * \param  E       muon energy [GeV]
     * \param  v       energy loss fraction
     * \return         angle [rad]
     */
    double ThetaRMSfromBrems(const double E, const double v) const
    {
      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;
    }
 
 
    /**
     * Pair production shower energy.
     *
     * \param  E     muon   energy [GeV]
     * \return       shower energy [GeV]
     */
    double EfromEErad(const double E) const
    {
      //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;
    }
 
 
    /**
     * Get RMS of scattering angle for pair production.
     *
     * \param  E       muon energy [GeV]
     * \param  v       energy loss fraction
     * \return         angle [rad]
     */
    double ThetaRMSfromEErad(const double E, const double v) const
    {
      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;
    }
 
 
    /**
     * Photo-nuclear shower energy.
     *
     * \param  E     muon   energy [GeV]
     * \return       shower energy [GeV]
     */
    double EfromGNrad(const double E) const
    {
      //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;
    }
 
 
    /**
     * Get RMS of scattering angle for photo-nuclear shower.
     *
     * \param  E       muon energy [GeV]
     * \param  v       energy loss fraction
     * \return         angle [rad]
     */
    double ThetaRMSfromGNrad(const double E, const double v) const
    {
      return 0.0;
    }
 
 
    /**
     * Ionization a parameter.
     *
     * \param  E       muon energy [GeV]
     * \return         ionization coefficient [GeV*m^2/g]
     */    
    virtual double CalculateACoeff(double E) const
    {
      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;
    }
 
 
    /**
     * Auxiliary data structure for handling member methods of class JRadiation.
     */
    struct radiation_type  {
      double (JRadiation::*sigma)(const double) const;                 //!< total cross section method
      double (JRadiation::*eloss)(const double) const;                 //!< energy loss method
      double (JRadiation::*theta)(const double, const double) const;   //!< scattering angle method
    };
 
 
    static constexpr radiation_type EErad_t = { &JRadiation::TotalCrossSectionEErad, &JRadiation::EfromEErad, &JRadiation::ThetaRMSfromEErad };
    static constexpr radiation_type Brems_t = { &JRadiation::TotalCrossSectionBrems, &JRadiation::EfromBrems, &JRadiation::ThetaRMSfromBrems };
    static constexpr radiation_type GNrad_t = { &JRadiation::TotalCrossSectionGNrad, &JRadiation::EfromGNrad, &JRadiation::ThetaRMSfromGNrad };
    
  protected:
    double GofZEvrho(const double E, 
                     const double v, 
                     const double r) const
    {
      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 IntegralofG(const double E,
                               const double eps) const
    {
      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 sigmaGammaPparam(const double eps)
    { 
      return (49.2 + 11.1 * log(eps) + 151.8/sqrt(eps));
    }
 
    static double PhiofEofepsparam(const double E,const double eps)
    {
      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 le   () { return 3.8616E-13;}
    static double r0   () { return 2.817940e-15; }
    static double Astar() { return 183.0; }
    static double B    () { return 183.0; }
    static double BP   () { return 1429.0; }
 
    const double Z;
    const double A;
    const double Dn;
    const double DnP;
    const int    steps;
    const double EminBrems;
    const double EminEErad;
    const double EminGNrad;
  };
}
 
#endif