#ifndef __JPHYSICS__JPDFTOOLKIT__
#define __JPHYSICS__JPDFTOOLKIT__

#include <cmath>

#include "JTools/JRange.hh"

#include "JPhysics/JConstants.hh"
#include "JPhysics/JPhysicsToolkit.hh"


/**
 * \file
 * Auxiliary methods for PDF calculations.
 * \author mdejong, bjung
 */

namespace JPHYSICS {}
namespace JPP { using namespace JPHYSICS; }

namespace JPHYSICS {

  using JTOOLS::JRange;

  
  static const double DELTARAY_TMIN = 0.000915499; //!< Minimum allowed delta-ray kinetic energy [GeV]
  static const double DELTARAY_TMAX = 1.0e10;      //!< Maximum allowed delta-ray kinetic energy [GeV]
  

  /**
   * Get minimal wavelength for PDF evaluations.
   *
   * \return            wavelength of light [nm]
   */
  inline double getMinimalWavelength()
  {
    return 300.0;
  }


  /**
   * Get maximal wavelength for PDF evaluations.
   *
   * \return            wavelength of light [nm]
   */
  inline double getMaximalWavelength()
  {
    return 700.0;
  }
  

  /**
   * Number of Cherenkov photons per unit track length and per unit wavelength.
   *
   * \param  lambda     wavelength of light [nm]
   * \param  n          index of refraction
   * \return            number of photons per unit track length and per unit wavelength [m^-1 nm^-1]
   */
  inline double cherenkov(const double lambda,
			  const double n)
  {
    const double x = n*lambda;
    
    return 1.0e9 * 2 * PI * ALPHA_ELECTRO_MAGNETIC * (n*n - 1.0) / (x*x);
  }


  /**
   * Get minimum delta-ray kinetic energy.
   *
   * \return     minimum delta-ray kinetic energy [GeV]
   */
  inline double getDeltaRayTmin()
  {
    const double Emin = MASS_ELECTRON / getSinThetaC(); // delta-ray Cherenkov threshold [GeV]

    return getKineticEnergy(Emin, MASS_ELECTRON);
  }


  /**
   * Get maximum delta-ray kinetic energy for given lepton energy and mass.\n
   * This formula is taken from reference
   * https://pdg.lbl.gov/2020/reviews/rpp2020-rev-passage-particles-matter.pdf
   * \n\n
   * N.B.: This function should not be used for electrons.\n
   *       For electrons use `0.5 * getKineticEnergy(E, MASS_ELECTRON)` instead.
   *
   * \param  E             particle energy                  [GeV]
   * \param  M             particle mass                    [GeV]
   * \return               maximum delta-ray kinetic energy [GeV]
   */
  inline double getDeltaRayTmax(const double E,
				const double M)
  {
    const double ratio = MASS_ELECTRON / M;

    const double gamma = E / M;
    const double beta  = sqrt((1.0 + 1.0/gamma) * (1.0 - 1.0/gamma));    

    return ( (2.0*MASS_ELECTRON*beta*beta*gamma*gamma) /
	     (1.0 + 2.0*gamma*ratio + ratio*ratio) );
  }			     


  /**
   * Get equivalent EM-shower energy due to delta-rays per unit track length\n
   * for an ionising particle with given energy and given mass.\n\n
   *
   * N.B: For this function a form factor \f$F(T) = \frac{1}{2E^2}T^2 - \frac{\beta^2}{T_{\rm max}} + 1\f$ is assumed, where\n:
   *      - \f$T\f$ corresponds to the delta-ray kinetic energy,\n
   *      - \f$T_{\rm max}\f$ to the maximum delta-ray kinetic energy,\n
   *      - \f$E\f$ to the ionising particle's energy and
   *      - \f$\beta\f$ to the corresponding particle velocity, relative to the speed of light.
   *
   * \param  E      particle energy                  [GeV]
   * \param  M      particle mass                    [GeV]
   * \param  Tmin   minimum delta-ray kinetic energy [GeV]
   * \param  Tmax   maximum delta-ray kinetic energy [GeV]
   * \param  Z      atomic number                    [unit]
   * \param  A      atomic mass                      [g/mol]
   * \return        equivalent energy loss           [GeV/m]
   */
  inline double getDeltaRays(const double E,
			     const double M,
			     const double Tmin,
			     const double Tmax,
			     const double Z,
			     const double A)
  {      
    if (Tmax > Tmin) {

      const double K = 0.307075;                            // [MeV mol^-1 cm^2]

      const double gamma = E / M;
      const double beta  = sqrt((1.0 + 1.0/gamma) * (1 - 1.0/gamma));
      
      const double a = 0.25/(E*E);
      const double b = beta*beta/Tmax;
      const double c = 1.0;

      const double W = 0.5 * K * (Z/A) * (1.0/(beta*beta)); // [MeV g^-1 cm^2]

      const double sT = Tmax + Tmin;
      const double dT = Tmax - Tmin;
      const double rT = Tmax / Tmin;
	
      const double weight = W * (a*sT*dT - b*dT + c*log(rT));

      return weight * DENSITY_SEA_WATER * 1.0e-1;           // [GeV m^-1]
	
    } else {

      return 0.0;
    }
  }


  /**
   * Get equivalent EM-shower energy due to delta-rays per unit track length\n
   * for an ionising particle with given energy and given mass and for a given form factor.\n
   * The template parameter corresponds to a class which contains an `operator()(const double)`\n
   * to compute the form factor corresponding to a given delta-ray kinetic energy.
   *
   * \param  E      particle energy                                     [GeV]
   * \param  M      particle mass                                       [GeV]
   * \param  Tmin   minimum delta-ray kinetic energy                    [GeV]
   * \param  Tmax   maximum delta-ray kinetic energy                    [GeV]
   * \param  Z      atomic number                                       [unit]
   * \param  A      atomic mass                                         [g/mol]
   * \param  F      form factor functor
   * \param  N      number of points for numeric integration
   * \return        equivalent energy loss                              [GeV/m]
   */
  template<class JFormFactor_t>
  inline double getDeltaRays(const double         E,
			     const double         M,
			     const double         Tmin,
			     const double         Tmax,
			     const double         Z,
			     const double         A,
			     const JFormFactor_t& F,
			     const int            N = 1000000)
  {
    using namespace JPP;

    if (Tmax > Tmin) {

      const double K = 0.307075;                            // [MeV mol^-1 cm^2]

      const double gamma = E / M;
      const double beta  = sqrt((1.0 + 1.0/gamma) * (1 - 1.0/gamma));

      const double W = 0.5 * K * (Z/A) * (1.0/(beta*beta)); // [MeV g^-1 cm^2]

      const double xmin = log(Tmin);
      const double xmax = log(Tmax);
      const double dx   = (xmax - xmin) / ((double) N);
	
      double weight = 0.0;

      for (double x = xmin; x <= xmax; x += dx) {

	const double T = exp(x);
	const double y = W * F(T) * dx;                     // [MeV g^-1 cm^-2]
	  
	weight += y;
      }

      return weight * DENSITY_SEA_WATER * 1.0e-1;           // [GeV m^-1]
	
    } else {

      return 0.0;
    }
  }


  /**
   * Equivalent EM-shower energy due to delta-rays per unit electron track length.\n\n
   *
   * N.B.: This function assumes a medium with atomic number Z=10 and atomic mass A=18\n
   *       and a form factor \f$F(T) = \frac{1}{2E^2}T^2 - \frac{\beta^2}{T_{\rm max}} + 1\f$, where\n:
   *       - \f$T\f$ corresponds to the delta-ray kinetic energy,\n
   *       - \f$T_{\rm max}\f$ to the maximum delta-ray kinetic energy,\n
   *       - \f$E\f$ to the ionising particle's energy and
   *       - \f$\beta\f$ to the corresponding particle velocity, relative to the speed of light.
   *
   * \param  E       electron energy                            [GeV]
   * \param  T_GeV   allowed kinetic energy range of delta-rays [GeV]
   * \return         equivalent energy loss                     [GeV/m]
   */  
  inline double getDeltaRaysFromElectron(const double         E,
					 const JRange<double> T_GeV = JRange<double>(DELTARAY_TMIN, DELTARAY_TMAX))
  {
    static const double Z = 10.0; // atomic number [unit]
    static const double A = 18.0; // atomic mass   [g/mol]
    
    const double Tmin = T_GeV.constrain(getDeltaRayTmin());
    const double Tmax = T_GeV.constrain(0.5 * getKineticEnergy(E, MASS_ELECTRON));
    
    return getDeltaRays(E, MASS_ELECTRON, Tmin, Tmax, Z, A);
  }
  
  
  /**
   * Equivalent EM-shower energy due to delta-rays per unit muon track length.\n\n
   * 
   * N.B.: This function assumes a medium with atomic number Z=10 and atomic mass A=18\n
   *       and a form factor \f$F(T) = \frac{1}{2E^2}T^2 - \frac{\beta^2}{T_{\rm max}} + 1\f$, where\n:
   *       - \f$T\f$ corresponds to the delta-ray kinetic energy,\n
   *       - \f$T_{\rm max}\f$ to the maximum delta-ray kinetic energy,\n
   *       - \f$E\f$ to the ionising particle's energy and
   *       - \f$\beta\f$ to the corresponding particle velocity, relative to the speed of light.
   *
   * \param  E       muon energy                                [GeV]
   * \param  T_GeV   allowed kinetic energy range of delta-rays [GeV]
   * \return         equivalent energy loss                     [GeV/m]
   */  
  inline double getDeltaRaysFromMuon(const double         E,
				     const JRange<double> T_GeV = JRange<double>(DELTARAY_TMIN, DELTARAY_TMAX))
  {
    static const double Z = 10.0; // atomic number [unit]
    static const double A = 18.0; // atomic mass   [g/mol]
    
    const double Tmin = T_GeV.constrain(getDeltaRayTmin());
    const double Tmax = T_GeV.constrain(getDeltaRayTmax(E, MASS_MUON));
    
    return getDeltaRays(E, MASS_MUON, Tmin, Tmax, Z, A);
  }


  /**
   * Equivalent EM-shower energy due to delta-rays per unit muon track length.\n\n
   * 
   * N.B.: This function assumes a medium with atomic number Z=10 and atomic mass A=18\n
   *       and a form factor \f$F(T) = \frac{1}{2E^2}T^2 - \frac{\beta^2}{T_{\rm max}} + 1\f$, where\n:
   *       - \f$T\f$ corresponds to the delta-ray kinetic energy,\n
   *       - \f$T_{\rm max}\f$ to the maximum delta-ray kinetic energy,\n
   *       - \f$E\f$ to the ionising particle's energy and
   *       - \f$\beta\f$ to the corresponding particle velocity, relative to the speed of light.
   *
   * \param  E       muon energy                                [GeV]
   * \return         equivalent energy loss                     [GeV/m]
   */  
  inline double getDeltaRaysFromMuonFit(const double E)
  {
    static const double a =  3.195e-01;
    static const double b =  3.373e-01;
    static const double c = -2.731e-02;
    static const double d =  1.610e-03;
    static const double Emin = 0.13078; // [GeV]
     
    if (E > Emin) {
      
      const double x = log10(E);                      //
      const double y = a + x*(b + x*(c + x*(d)));     // [MeV g^-1 cm^2]
      
      if (y > 0.0) {
	return y * DENSITY_SEA_WATER * 1.0e-1;        // [GeV/m]
      }
    }
    
    return 0.0;
  }

  
  /**
   * Equivalent EM-shower energy due to delta-rays per unit tau track length.\n\n
   *
   * N.B.: This function assumes a medium with atomic number Z=10 and atomic mass A=18\n
   *       and a form factor \f$F(T) = \frac{1}{2E^2}T^2 - \frac{\beta^2}{T_{\rm max}} + 1\f$, where\n:
   *       - \f$T\f$ corresponds to the delta-ray kinetic energy,\n
   *       - \f$T_{\rm max}\f$ to the maximum delta-ray kinetic energy,\n
   *       - \f$E\f$ to the ionising particle's energy and
   *       - \f$\beta\f$ to the corresponding particle velocity, relative to the speed of light.
   *
   * \param  E       tau energy                                 [GeV]
   * \param  T_GeV   allowed kinetic energy range of delta-rays [GeV]
   * \return         equivalent energy loss                     [GeV/m]
   */  
  inline double getDeltaRaysFromTau(const double         E,
				    const JRange<double> T_GeV = JRange<double>(DELTARAY_TMIN, DELTARAY_TMAX))				    
  {
    static const double Z = 10.0; // atomic number [unit]
    static const double A = 18.0; // atomic mass   [g/mol]
    
    const double Tmin = T_GeV.constrain(getDeltaRayTmin());
    const double Tmax = T_GeV.constrain(getDeltaRayTmax(E, MASS_TAU));
    
    return getDeltaRays(E, MASS_TAU, Tmin, Tmax, Z, A);
  }


  /**
   * Equivalent EM-shower energy due to delta-rays per unit tau track length.\n\n
   *
   * N.B.: This function assumes a medium with atomic number Z=10 and atomic mass A=18\n
   *       and a form factor \f$F(T) = \frac{1}{2E^2}T^2 - \frac{\beta^2}{T_{\rm max}} + 1\f$, where\n:
   *       - \f$T\f$ corresponds to the delta-ray kinetic energy,\n
   *       - \f$T_{\rm max}\f$ to the maximum delta-ray kinetic energy,\n
   *       - \f$E\f$ to the ionising particle's energy and
   *       - \f$\beta\f$ to the corresponding particle velocity, relative to the speed of light.
   *
   * \param  E       tau energy                                 [GeV]
   * \return         equivalent energy loss                     [GeV/m]
   */  
  inline double getDeltaRaysFromTauFit(const double E)
  {
    static const double a = -2.178e-01;
    static const double b =  4.995e-01;
    static const double c = -3.906e-02;
    static const double d =  1.615e-03;
    static const double Emin = 2.19500; // [GeV]
       
    if (E > Emin) {
   
      const double x = log10(E);                      //
      const double y = a + x*(b + x*(c + x*(d)));     // [MeV g^-1 cm^2]
   
      if (y > 0) {
	return y * DENSITY_SEA_WATER * 1.0e-1;        // [GeV/m]
      }
    }
   
    return 0.0;
  }


  /**
   * Emission profile of photons from delta-rays.
   *
   * Profile is taken from reference ANTARES-SOFT-2002-015, J.\ Brunner (fig.\ 3).
   *
   * \param  x      cosine emission angle
   * \return        probability
   */
  inline double getDeltaRayProbability(const double x)
  {
    //return 1.0 / (4.0 * PI);
    return 0.188 * exp(-1.25 * pow(fabs(x - 0.90), 1.30));
  }


  /**
   * Rayleigh cross section.
   *
   * \param  n          index of refraction
   * \param  lambda     wavelength of light [nm]
   * \return            cross section       [cm^2]
   */
  inline const double getRayleighCrossSection(const double n,
					      const double lambda)
  {
    static const double d = 0.36;                    // size of H2O molecule [nm]
    static const double U = PI*PI*PI*PI*PI*2.0/3.0;
    static const double V = d*d*d*d*d*d;

    const double W     = (n*n - 1.0) / (n*n + 2.0);
    const double sigma = 1.0e-14 * U*V*W*W / (lambda*lambda*lambda*lambda);   // [cm^2]

    return sigma;
  }


  /**
   * Rayleigh scattering length.
   *
   * \param  n          index of refraction
   * \param  lambda     wavelength of light [nm]
   * \return            scattering length   [m]
   */
  inline const double getRayleighScatteringLength(const double n,
						  const double lambda)
  {
    static const double amu = 18.01528; // H20 mass in Atomic mass units

    const double sigma = getRayleighCrossSection(n, lambda);
    const double ls    = 1.0e-2 / (DENSITY_SEA_WATER * AVOGADRO * sigma / amu);     // [m]

    return ls;
  }
}

#endif