#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
   *
   * Muon scattering:
   * A. van Ginneken, Nucl.\ Instr.\ and Meth.\ A251 (1986) 21. 
   */
  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;
    }
    
  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;
  };


  /**
   * Auxiliary data structure for handling member methods of class JRadiation.
   */
  struct JRadiationSource_t  {
    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 const JRadiationSource_t EErad_t = { &JRadiation::TotalCrossSectionEErad, &JRadiation::EfromEErad, &JRadiation::ThetaRMSfromEErad };
  static const JRadiationSource_t Brems_t = { &JRadiation::TotalCrossSectionBrems, &JRadiation::EfromBrems, &JRadiation::ThetaRMSfromBrems };
  static const JRadiationSource_t GNrad_t = { &JRadiation::TotalCrossSectionGNrad, &JRadiation::EfromGNrad, &JRadiation::ThetaRMSfromGNrad };
}

#endif