#ifndef __JPHYSICS__JRADIATIONFUNCTION__
#define __JPHYSICS__JRADIATIONFUNCTION__

#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JTools/JMultiFunction.hh"
#include "JTools/JGrid.hh"

#include "JPhysics/JRadiation.hh"


/**
 * \author mdejong
 */

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

namespace JPHYSICS {

  /**
   * Fast implementation of class JRadiation.
   *
   * In this, the methods
   *   - JRadiation::TotalCrossSectionEErad,
   *   - JRadiation::TotalCrossSectionGNrad,
   *   - JRadiation::CalculateACoeff and
   *   - JRadiation::IntegralofG
   *
   * are reimplemented using lookup tables.
   */
  class JRadiationFunction : 
    public JRadiation
  {
    typedef JTOOLS::JGridPolint1Function1D_t                                                       JFunction1D_t;
    typedef JTOOLS::JMultiFunction<JFunction1D_t, 
				   JTOOLS::JMapList<JTOOLS::JPolint1FunctionalGridMap> >           JFunction2D_t;

  public:
    /**
     * Constructor.
     *
     * \param  radiation       JRadiation object
     * \param  number_of_bins  number of bins
     * \param  Emin            minimal muon energy [GeV]
     * \param  Emax            maximal muon energy [GeV]
     */
    JRadiationFunction(const JRadiation&  radiation,
		       const unsigned int number_of_bins,
		       const double       Emin,
		       const double       Emax) :
      JRadiation(radiation)
    {
      using namespace JTOOLS;

      const double xmin = log(Emin);
      const double xmax = log(Emax);

    
      integral.configure(make_grid(number_of_bins, xmin, xmax));

      for (JFunction2D_t::iterator i = integral.begin(); i != integral.end(); ++i) {

	const double x    = i->getX();
	const double E    = exp(x);

	const double ymin = log(EminEErad/E);
	const double ymax = 0.0;

	if (ymax > ymin) {

	  i->getY().configure(make_grid(number_of_bins, ymin, ymax));
	  
	  for (JFunction1D_t::iterator j = i->getY().begin(); j != i->getY().end(); ++j) {
	    
	    const double y   = j->getX();
	    const double eps = exp(y) * E;
	    const double z   = JRadiation::IntegralofG(E, eps);

	    j->getY() = z;
	  }
	}
      }
    
      integral.compile();
      integral.setExceptionHandler(new JFunction1D_t::JDefaultResult(0.0));


      sigmaEE.configure(make_grid(number_of_bins, xmin, xmax));

      for (JFunction1D_t::iterator i = sigmaEE.begin(); i != sigmaEE.end(); ++i) {

	const double E = exp(i->getX());
	const double y = JRadiation::TotalCrossSectionEErad(E);

	i->getY() = y;
      }
    		    
      sigmaEE.compile();


      sigmaGN.configure(make_grid(number_of_bins, xmin, xmax));

      for (JFunction1D_t::iterator i = sigmaGN.begin(); i != sigmaGN.end(); ++i) {

	const double E = exp(i->getX());
	const double y = JRadiation::TotalCrossSectionGNrad(E);

	i->getY() = y;
      }
		    
      sigmaGN.compile();


      Acoeff.configure(make_grid(number_of_bins, xmin, xmax));

      for (JFunction1D_t::iterator i = Acoeff.begin(); i != Acoeff.end(); ++i) {
     
        const double E = exp(i->getX());
        const double y = JRadiation::CalculateACoeff(E);

        i->getY() = y;
      }
    
      Acoeff.compile();
    }
       

    /**
     * Pair production cross section.
     *
     * \param  E     muon    energy [GeV]
     * \return       cross section [m^2/g]
     */
    virtual double TotalCrossSectionEErad(const double E) const override 
    {
      const double x = log(E);

      if (x >= sigmaEE. begin()->getX() &&
	  x <= sigmaEE.rbegin()->getX()) {

	try {
	  return sigmaEE(x);
	}
	catch(std::exception& error) {}
      }

      return JRadiation::TotalCrossSectionEErad(E);
    }
  

    /**
     * Photo-nuclear cross section.
     *
     * \param  E     muon    energy [GeV]
     * \return       cross section [m^2/g]
     */
    virtual double TotalCrossSectionGNrad(const double E) const override 
    {
      const double x = log(E);

      if (x >= sigmaGN. begin()->getX() &&
	  x <= sigmaGN.rbegin()->getX()) {

	try {
	  return sigmaGN(x);
	}
	catch(std::exception& error) {}
      }

      return JRadiation::TotalCrossSectionGNrad(E);
    }


    /**
     * Ionization a parameter.
     *
     * \param  E       muon energy [GeV]
     * \return         ionization coefficient [GeV*m^2/g]
     */
    virtual double CalculateACoeff(const double E) const override 
    {
      const double x = log(E);

      if (x >= Acoeff. begin()->getX() &&
	  x <= Acoeff.rbegin()->getX()) {

	try {
	  return Acoeff(x);
	}
	catch(std::exception& error) {}
      }

      return JRadiation::CalculateACoeff(E);
    }

  protected:
    virtual double IntegralofG(const double E,
			       const double eps) const
    {
      try {
	
	const double x = log(E); 
	const double y = log(eps/E); 
	const double z = integral(x, y);

	return z;
      }
      catch(std::exception& error) {}

      return JRadiation::IntegralofG(E, eps);
    }

    JFunction1D_t sigmaEE;
    JFunction1D_t sigmaGN;
    JFunction1D_t Acoeff;
    JFunction2D_t integral;
  };
}

#endif