#ifndef __JPHYSICS__JDIS__
#define __JPHYSICS__JDIS__

#include <cmath>

#include <TRandom3.h>

#include "JPhysics/JConstants.hh"

/**
 * \file
 * Deep-inelastic muon-nucleon scattering.
 *
 * \author mdejong
 */

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

namespace JPHYSICS {

  /**
   * Deep-inelastic muon-nucleon scattering.
   *
   * For cross section, see reference:\n
   * D.A. Timashkov and A.A. Petrukhin, 29th International Cosmic Ray Conference Pune (2005) 9, 89-92.
   *
   * For kinematics, see reference:\n
   * A. van Ginneken, Nucl.\ Instr.\ and Meth.\ A251 (1986) 21.
   */
  class JDIS {
  public:
    /**
     * Default constructor.
     */
    JDIS()
    {}


    /**
     * Get cross section.
     *
     * \param  E         muon energy   [GeV]
     * \return           cross section [cm^2]
     */
    inline double getCrossSection(const double E) const
    {
      const double x = log10(E*1.0e+1);

      if (E > JDIS_t::E0)
	return 0.35e-30 * pow(1.0 - JDIS_t::E0/E, 2.0) * exp(-2.10/x) * pow(10.0, 0.125*x);
      else
	return 0.0;
    }


    /**
     * Get probability of given energy fraction.
     *
     * \param  E         muon energy [GeV]
     * \param  v         energy fraction
     * \return           probability
     */
    double getP(const double E, const double v) const
    {
      const JDIS_t dis(E);

      return dis.getP(v);
    }
  

    /**
     * Get shower energy.
     *
     * \param  E         muon energy [GeV]
     * \return           shower energy [GeV]
     */
    double getE(const double E) const
    {
      const JDIS_t dis(E);

      return dis.getE();
    }


    /**
     * Get RMS of scattering angle.
     *
     * \param  E         muon energy [GeV]
     * \param  v         energy loss fraction
     * \return           angle [rad]
     */
    double getThetaRMS(const double E, const double v) const
    {
      const double u = 1.0 - v;

      if (E*v > 0.135)
	return (0.39 / (E*u)) * pow(sqrt(E)*u*v, 0.17) * (1.0 - 0.135/(E*v));
      else
	return  0.0;
    }


    /**
     * Get breakpoint.
     *
     * \param  E         muon energy [GeV]
     * \return           energy fraction
     */
    double getV(const double E) const
    {
      return JDIS_t::E1/E;
    }


    /**
     * Auxiliary helper class for kinematics of deep-inelastic muon-nucleon scattering at fixed muon energy.
     */
    class JDIS_t {
    public:
      /**
       * Constructor.
       *
       * \param  E         muon energy [GeV]
       */
      JDIS_t(const double E) :
	E (E),
	E2(E - MASS_MUON)
      {
	k2 = pow(E1, p0 - p1) / pow(1.0 - E1/E, 2);

	const double A = getA(E1/E, false) - getA(E0/E, false);
	const double B = getB(E2/E, false) - getB(E1/E, false);

	k1 = 1.0 / (A + k2 * B);
      }


      /**
       * Get probability of given energy fraction.
       *
       * \param  v         energy fraction
       * \return           probability
       */
      double getP(const double v) const
      {
	const double x = E*v;
	const double u = 1.0 - v; 

	if      (x > E2)
	  return 0.0;
	else if (x > E1)
	  return k1 * k2 * pow(x, p1) * u*u;
	else if (x > E0)
	  return k1 *      pow(x, p0);
	else
	  return 0.0;
      }
  

      /**
       * Get shower energy.
       *
       * \return           shower energy [GeV]
       */
      double getE() const
      {
	const double A = getA(E1/E) - getA(E0/E);
	const double B = getB(E2/E) - getB(E1/E);

	for ( ; ; ) {

	  double y = gRandom->Rndm() * (A + B);

	  if (y <= A) {

	    y += getA(E0/E);
	    y /= k1 * pow(E,p0) / (p0 + 1);

	    const double v = pow(y, 1.0 / (p0 + 1));

	    return E * v;

	  } else {

	    y -= A;
	    y *= (getb(E2/E) - getb(E1/E)) / B;
	    y +=  getb(E1/E);
	    y /= k1 * k2 * pow(E,p1) / (p1 + 1);
	  
	    const double v = pow(y, 1.0 / (p1 + 1));
	    const double u = 1.0 - v; 

	    if (gRandom->Rndm() < u*u) {
	      return E * v;
	    }
	  }
	}
      }

      const double E;                        //!< actual  energy [GeV]
      const double E2;                       //!< maximal energy [GeV]

      static constexpr double E0 = 0.144;    //!< minimal energy [GeV]
      static constexpr double E1 = 0.35;     //!< breakpoint     [GeV]

      static constexpr double p0 =  2.0;     //!< spectral index upto breakpoint
      static constexpr double p1 = -1.11;    //!< spectral index from breakpoint

    private:
      /**
       * Integral upto breakpoint.
       *
       * \param  v         energy fraction
       * \param  option    option (if false, suppress normalisation constant)
       * \return           integral
       */    
      double getA(const double v, const bool option = true) const
      {
	return (option ? k1 : 1.0)      * pow(E, p0) * pow(v, p0 + 1) * (1.0 / (p0 + 1));
      }


      /**
       * Integral from breakpoint.
       *
       * \param  v         energy fraction
       * \param  option    option (if false, suppress normalisation constant)
       * \return           integral
       */    
      double getB(const double v, const bool option = true) const
      {
	return (option ? k1 * k2 : 1.0) * pow(E, p1) * pow(v, p1 + 1) * (1.0 / (p1 + 1)  -
									 2*v / (p1 + 2)  +
									 v*v / (p1 + 3));
      }


      /**
       * Integral from breakpoint without suppression factor.
       *
       * \param  v         energy fraction
       * \return           integral
       */    
      double getb(const double v) const
      {
	return k1 * k2 * pow(E, p1) * pow(v, p1 + 1) / (p1 + 1);
      }

      double k1;                             //!< normalisation constant upto breakpoint
      double k2;                             //!< normalisation constant from breakpoint
    };
  };
}

#endif