#ifndef __JPHYSICS__JGEANX__
#define __JPHYSICS__JGEANX__

#include <cmath>

#include "JMath/JMathSupportkit.hh"
#include "JPhysics/JConstants.hh"

/**
 * \file
 * Photon emission profile EM-shower.
 * \author mdejong
 */

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

namespace JPHYSICS {

  /**
   * Probability density function of photon emission from EM-shower as a function of cosine of the emission angle.
   *
   *          \f[ P(\cos(\theta))  =  c \times e^{b \times |\cos(\theta) - 1/n|^a} \f]
   * 
   * where \f$ c \f$ is a normalisation constant such that the integral of \f$ P \f$ over the full solid angle is one.
   *
   * The parametrisation is taken from reference:
   * R. Mirani, "Parametrisation of EM-showers in the ANTARES detector volume.",
   * Doctoral thesis in computational physics, University of Amsterdam.
   */
  class JGeanx 
  {
  public:
    /**
     * Constructor.
     * 
     * \param  __a     power
     * \param  __b     exponential slope
     * \param  __n     index of refraction
     */
    JGeanx(const double __a,
	   const double __b,
	   const double __n = getIndexOfRefractionPhase()) :
      a(__a),
      b(__b),
      n(__n)
    {
      const double y = evaluate(-1.0, +1.0);

      c  = 1.0 / (y * 2*PI);
    }


    /**
     * Number of photons from EM-shower as a function of emission angle.
     * The integral over full solid angle is normalised to one.
     *
     * \param  ct      cosine angle of emmision
     * \return         d^2P/dcos()dphi
     */
    double operator()(const double ct) const
    {
      return c * evaluate(ct);
    }


    /**
     * Integral number of photons from EM-shower between two emission angles.
     * The integral over full solid angle is normalised to one.
     *
     * \param  xmin    minimal cosine angle of emmision
     * \param  xmax    maximal cosine angle of emmision
     * \return         dnpe/dphi 
     */
    double operator()(const double xmin, 
		      const double xmax) const
    {
      return c * evaluate(xmin, xmax);
    }


    /**
     * Functional dependence.
     * In this, the normalisation constant c = 1.
     *
     * \param  ct      cosine angle of emmision
     * \return         f()
     */
    double evaluate(const double ct) const
    {
      const double x = fabs(ct - 1.0/n);

      return exp(b * pow(x,a));
    }


    /**
     * Integral.
     * In this, the normalisation constant <tt>c = 1</tt>.
     *
     * \param  xmin    minimal cosine angle of emmision
     * \param  xmax    maximal cosine angle of emmision
     * \return         integral f() x dcos()
     */
    double evaluate(const double xmin, 
		    const double xmax) const
    {
      using namespace std;
      using namespace JPP;

      const double x  = 1.0 / n;
      
      const double ai = 1.0 / a;
      const double bi = 1.0 / b;
      
      const double xl = pow(fabs(x - xmin), a) * -b;
      const double xr = pow(fabs(x - xmax), a) * -b;
      
      const double gp = tgamma(ai);
      
      double yl = Gamma(ai, xl) * gp;
      double yr = Gamma(ai, xr) * gp;
      
      if (xmin > x) 
	yl = -yl;
      if (xmax < x)
	yr = -yr;
      
      return ai * pow(-bi,ai) * (yl + yr);
    }


    const double a;     //!< power
    const double b;     //!< slope
    const double n;     //!< index of refraction

  private:
    double c;           //!< normalisation constant
  };


  /**
   * Function object for the number of photons from EM-shower as a function of emission angle.
   */
  static const JGeanx geanx(0.35, -5.40);
}

#endif