#ifndef __CASCADE_SHOWER_PDF__
#define __CASCADE_SHOWER_PDF__

#include <iostream>
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JNPETable.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JPhysics/JGeanz.hh"

#include "../utilJpp.hh"
#include "Hit.hh"
#include "Trk.hh"
#include "Det.hh"
#include "io_util.hh"
#include "rootutil.hh"

#include "TRandom.h"
#include "Params.hh"

using JPHYSICS::geanz; // longitudinal shower profile function

//-----------------------
// Jpp type incantations.
//--------------------

using namespace JPP;

typedef JPolintFunction1D<1, JPolintElement2S<double, double>, JCollection, JResultPDF<double> >    PDF_Function1D_t;
//typedef JPolintFunction1D<1, JElement2D<double, double>, JCollection, double >                      CDF_Function1D_t;

typedef JMAPLIST<JPolint1FunctionalMap,
                 JPolint1FunctionalMap,
                 JPolint1FunctionalGridMap,
                 JPolint1FunctionalGridMap>::maplist            JMapList_t;
typedef JPDFTable<PDF_Function1D_t, JMapList_t>                 JCascadePdf_t;
typedef JNPETable<double, double, JMapList_t>                   JCascadeNpe_t;

template< typename ...Ps> 
inline double check_ (double x, const char* fun, int line, Ps ...ps  ) {
  if ( std::isfinite( x )) return x;
  fatal (" non-finite number ", x, " in ", fun,", line:",line, "\n", ps...);
}
#define check(X, ...) check_( X, __FUNCTION__, __LINE__, __VA_ARGS__ )


struct CascadePdf
{
  enum light_type_t { total = 0, scattered = 1, direct = 2 };

  vector<double> _lengths;
  vector<JCascadePdf_t> _pdf_funcs;
  vector<JCascadeNpe_t> _npe_funcs;
  JCascadePdf_t* _pdf_func;
  JCascadeNpe_t* _npe_func;

  light_type_t light_type;
  int nsamples;  
  double background_rate;
  double window_start;
  double window_end;
  double energy_fraction_sample;

  CascadePdf( string shower_scattered_pdf = "/pbs/throng/km3net/src/Jpp/master/data/J14p.dat", //Scattered one first
              string shower_direct_pdf    = "/pbs/throng/km3net/src/Jpp/master/data/J13p.dat", 
              double sigma_blur           = 0, 
              uint elong_steps            = 0, //0 : do nothing, 1 = move to shower max, > 1 : sample
              double R_bg                 = 0,
              double t_start              = -100,
              double t_end                = 900,
              double E_extra_sample       = 0.0 );

  
  //! set the pdf/npe function to direct,scattered or total
  void set_light_type( light_type_t type )
  {
    _pdf_func = &_pdf_funcs[ type ];
    _npe_func = &_npe_funcs[ type ];
  }

  void set_nsamples( unsigned n ) 
  { 
    nsamples = n; 
  }

  void set_background( double rate )
  {
    background_rate = rate;
  }

  void set_energy_fraction( double fraction )
  {
    energy_fraction_sample = fraction;
  }


  double get_k40_npe( double time_window ) 
  {

    if ( time_window > 0 ) {
      return background_rate*1e-9*time_window;
    }
    else {
      return 0; // for residuals smaller than -100, we don't add background
    }

  }
  
  void set_lengths( const double E ){
    _lengths.resize(1);
    if      (nsamples==0) _lengths[0] = 0.0;
    else if (nsamples==1) _lengths[0] = geanz.getMaximum(E);
    else {
      _lengths.resize(nsamples);
      for( int ii = 0; ii < nsamples; ii++) _lengths[ii] = geanz.getLength( E, (ii + 0.5)/( (double) nsamples ) );
    };
  };

  JCascadePdf_t::result_type eval1(const Paramst& p, JCascadePdf_t* f ); 
  JCascadeNpe_t::result_type eval1(const Params & p, JCascadeNpe_t* f );
  JCascadePdf_t::result_type eval1_opt(const Paramst& p, JCascadePdf_t* f ); 
  JCascadeNpe_t::result_type eval1_opt(const Params & p, JCascadeNpe_t* f );
  JCascadePdf_t::result_type eval1(double D, double cd, double theta, double phi, double t, double E) ;
  JCascadeNpe_t::result_type eval1(double D, double cd, double theta, double phi, double E);

  //! Evaluate the 'elongated' function.
  template<typename F, typename Pars>
  typename F::result_type eval(F* func , const Pars& p )
  {
    if (nsamples==0) return eval1( p                                    , func );
    if (nsamples==1) return eval1( p.propagate( geanz.getMaximum(p.E) ) , func );

    typename F::result_type r;

    // cout << endl << "r " << typeid(r).name() << endl;

    for( int i = 0 ; i < nsamples ; i++) 
      {
	double l = geanz.getLength( p.E, (i + 0.5)/( (double) nsamples ) );
	auto x = eval1( p.propagate( l ) , func ) / (double) nsamples;

	if (i==0) r = x; // the horror
	else      r+= x;
      }

    if ( energy_fraction_sample > 0 ) 
      {
	Pars p_extra(p);
	p_extra.E = p_extra.E*energy_fraction_sample;
	r+= eval1( p_extra , func );
      }
    return r;
  }

  //! Evaluate the 'elongated' function.
  template<typename F, typename Pars>
  typename F::result_type eval_opt(F* func , const Pars& p )
  {  
    typename F::result_type r;
     
    uint ii = 0;
    for( vector<double>::iterator il = _lengths.begin(); il != _lengths.end(); il++, ii++) 
      {
	auto x = eval1_opt( p.propagate( *il ), func ) / (double) nsamples;
	 
	if (ii==0) r = x;
	else       r+= x;

      }

    return r;
  }

  // ---------------------------------------------------------------------
  // npe density function 
  // ---------------------------------------------------------------------

  double dN_dt_( const Paramst& p )
  {
    // cout << "dN_dt "; p.print(); 
    return eval( _pdf_func, p  ).f + get_k40_npe(1);
  }
  
  // anything you can pass to Paramst, you can also pass to dN_dt
  template<typename ...Ts> double dN_dt( Ts ...args ) { return dN_dt_( Paramst(args...)); }


  // ---------------------------------------------------------------------
  // Cummulative number of npe up to time t
  // ---------------------------------------------------------------------

  double N_( const Paramst& p )
  {
    return eval( _pdf_func, p  ).v + get_k40_npe( p.t - window_start );
  }

  template<typename ...Ts> double N( Ts ...args ) { return N_( Paramst(args...)); }


  // ---------------------------------------------------------------------
  // Pdf of the first hit time (normalized)
  // -------------------------------------------------------------------

  double dP1_dt_( const Paramst& p  )
  {
    auto z = eval( _pdf_func, p );
    double ntot = Ntot( p );

    /*
    cout << "z.f: " << z.f << endl;
    cout << "z.f (n): " << (z.f + get_k40_npe(1) ) << endl;
    cout << "z.v (N): " << exp( -z.v - get_k40_npe( p.t - window_start ) ) << endl;
    cout << "c: " <<  1 - exp( -ntot ) << endl;
    */

    return (z.f + get_k40_npe(1) ) * exp( -z.v - get_k40_npe( p.t - window_start ) ) / ( 1 - exp( -ntot ) );
  }

  // ---------------------------------------------------------------------
  // Pdf of the first hit time
  // -------------------------------------------------------------------

  double dP1_dt_a_( const Paramst& p  )
  {
    auto z = eval( _pdf_func, p );
    //cout << "p: " << p << endl;
    //cout << "z.f: " << z.f << endl;
    //cout << "z.v: " << z.v << endl;
    //cout << "p.t: " << p.t << endl;
    return ( z.f + get_k40_npe(1) ) * exp( -z.v - get_k40_npe( p.t - window_start ) );
  }

  template<typename ...Ts> double dP1_dt  ( Ts ...args ) { return dP1_dt_  ( Paramst(args...)); }
  template<typename ...Ts> double dP1_dt_a( Ts ...args ) { return dP1_dt_a_( Paramst(args...)); }
  //template<typename ...Ts> double dP1_dt2( Ts ...args ) { return dP1_dt_2( Paramst(args...)); }

  // ---------------------------------------------------------------------
  // Pdf of the first hit time (normalized) for double showers
  // -------------------------------------------------------------------
  
  double dP1_dt_double( const Paramst& p1, const Paramst& p2  )
  {
    auto z1 = eval( _pdf_func, p1 );
    auto z2 = eval( _pdf_func, p2 );
    double ntot = Ntot( p1 ) + Ntot( p2 );
    return (z1.f + z2.f + get_k40_npe(1) ) * exp( - z1.v - z2.v - get_k40_npe( p1.t - window_start ) ) / ( 1 - exp( -ntot ) );
  }
  
  
  // ---------------------------------------------------------------------
  // Cdf of the first hit time 
  // -------------------------------------------------------------------

  double P1_( const Paramst& p  )
  {
    // cout << "P1 "; p.print(); 
    return ( 1 - exp( -N( p ) ) ) / ( 1 - exp( -Ntot( p ) ) );
  }

  template<typename ...Ts> double P1( Ts ...args ) { return P1_( Paramst(args...)); }


  
  // ---------------------------------------------------------------------
  // total number of photons, integrated over all time
  // ---------------------------------------------------------------------

  double Ntot_( const Params& p  ) 
  {
    // cout << "Ntot "; p.print(); 
    return eval( _npe_func, p  ) + get_k40_npe( window_end - window_start );
  }

  template<typename ...Ts> double Ntot( Ts ...args ) { return Ntot_( Params(args...)); };

  double dnpe_dt( const Paramst &p );
  double npe    ( const Paramst &p );
  double npe    ( const Params  &p );
  //double npe_V  ( const Paramst  &p );

  double dnpe_dt ( double E, double D, double cd, double theta, double phi, double t );
  double npe     ( double E, double D, double cd, double theta, double phi, double t );
  double npe     ( double E, double D, double cd, double theta, double phi );

};


#endif