#ifndef __JPHYSICS__JGEANT__
#define __JPHYSICS__JGEANT__

#include <cmath>
#include <algorithm>

#include "JPhysics/JConstants.hh"
#include "JIO/JSerialisable.hh"
#include "JPhysics/JGeant_t.hh"
#include "JPhysics/JGeanx.hh"

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

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

namespace JPHYSICS {

  using JIO::JReader;
  using JIO::JWriter;

  /**
   * Function object for the probability density function of photon emission from EM-shower 
   * as a function of the index of refraction and the cosine of the emission angle.
   */
  class JGeant : 
    public JGeant_t
  {
  public:
    /**
     * Default constructor.
     */
    JGeant() 
    {}


    /**
     * Constructor.
     * 
     * \param  geanx   PDF of EM-shower
     * \param  dx      step size for interpolation of function values
     */
    JGeant(const JGeanx& geanx,
	   const double  dx = 0.0001) 
    {
      using JTOOLS::PI;

      const double n0 = geanx.n - 0.02;
      const double n1 = geanx.n + 0.02;

      const JGeanx g0(geanx.a, geanx.b, n0);
      const JGeanx g1(geanx.a, geanx.b, n1);

      const double y0 = g0.evaluate(-1.0, +1.0);
      const double y1 = g1.evaluate(-1.0, +1.0);

      const double c  = 1.0 / (y0 * 2*PI);

      // following range corresponds to n = 1 to n = infinity

      for (double x = -2.0; x < 1.0; x += dx) {
	(*this)[x] = c * g0.evaluate(x + 1.0/n0);
      }

      this->compile();

      // linear approximation of dependence normalisation constant on index of refraction
      
      this->a1 = (y1 - y0) / y0 / (n1 - n0);
      this->a0 = 1.0 + a1*n0;
    }
  };


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

#endif