#ifndef __JCDFTABLE1D__
#define __JCDFTABLE1D__

#include <memory>

#include "JPhysics/JPDFTable.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"


namespace JPHYSICS {

  namespace {
    using JPHYSICS::JCDFTable;
    using JPHYSICS::JFunctionTransformer;
    
    using JLANG::JMapList;
    
    using JTOOLS::PI;
    using JTOOLS::JElement2D;
    using JTOOLS::JPolintFunction1D;
    using JTOOLS::JGrid;
    using JTOOLS::JGridBounds;
    using JTOOLS::JPolint1FunctionalGridMap;
    using JTOOLS::JMultiFunction;
    using JTOOLS::JFunctionTransformerInterface;
  }


  /**
   * Custom class for CDF table in 1 dimension (only number of photo-electrons).
   */
  template<class JArgument_t, class JResult_t>
  class JCDFTable1D :
    public JPolintFunction1D<1, JElement2D<JArgument_t,JResult_t>, JGrid>
  {
  public:

    typedef JArgument_t                                            argument_type;
    typedef JResult_t                                              result_type;
    typedef JElement2D<argument_type, result_type>                 element_type;

    typedef JPolintFunction1D<1, element_type, JGrid>              JFunction1D_t;

    typedef JFunctionTransformerInterface<1, argument_type>        transformer_type;
  

    /**
     * Default constructor.
     */
    JCDFTable1D() :
      JFunction1D_t()
    {}


    /**
     * Constructor.
     *
     * \param  cdf             CDF            (4D table)
     * \param  number_of_bins  number of bins (1D table)
     * \param  safety_factor   safety factor
     */
    template<class JFunction1D_s, class JMapList_s>
    JCDFTable1D(const JCDFTable<JFunction1D_s, JMapList_s>& cdf,
		const int    number_of_bins,
		const double safety_factor = 1.1) :
      JFunction1D_t()
    {

      // extract the weight function

      try {

	typedef JFunctionTransformer<3, argument_type> JFunctionTransformer_t;

	const JFunctionTransformer_t& object = dynamic_cast<const JFunctionTransformer_t&>(*(cdf.transformer));

	transformer.reset(new JFunctionTransformer<1, argument_type>(object.transformer));
      }
      catch(const std::exception& error) {

	transformer.reset(new transformer_type());
      }


      // build the 1D table

      configure(JGridBounds<argument_type>(number_of_bins, cdf.intensity.begin()->first, cdf.intensity.rbegin()->first));

      for (typename JFunction1D_t::iterator j = this->begin(); j != this->end(); ++j) {

	const double R  = j->getKey(); 
	const double W  = transformer->getWeight(R);

	j->setValue(0.0);

	for (double theta = 0.0; theta <= PI; theta += 0.01 * PI) {

	  try {
	    
	    const double y = safety_factor * cdf.getNPE(R, theta, 0.99999*PI) / W;
	    
	    if (y > j->getValue())
	      j->setValue(y);
	  }
	  catch(const std::exception& error) {
	    std::cerr << error.what() << std::endl;
	  }
	}
      }

      this->compile();
    }


    /**
     * Get number of photo-electrons.
     *
     * \param  R          distance between muon and PMT [m]
     * \return            number of photo-electrons
     */
    double getNPE(const double R) const
    {
      const double buffer[] = { R };

      const double y = this->evaluate(buffer);
      const double W = transformer->getWeight(R);

      return y * W;
    }


    std::auto_ptr<transformer_type> transformer;
  };
}

#endif