#include "../utilJpp.hh"
#include "CascadePdf.hh"


CascadePdf::CascadePdf( string shower_scattered_pdf,// = "/sps/km3net/users/jseneca/jpp/data/J14p.dat", //Scattered one first
			string shower_direct_pdf   ,// = "/sps/km3net/users/jseneca/jpp/data/J13p.dat", 
			double sigma_blur          ,// = 0, 
			uint elong_steps           ,
			double R_bg                ,// background rate in [Hz]
			double t_start             ,// start of hit time residual window 
			double t_end,
			double E_extra_sample ):
  nsamples( elong_steps ), // /= 1 ) : // 0 : do nothing, 1 = move to shower max, > 1 : sample
  light_type( total ),
  background_rate( R_bg ),
  window_start( t_start ),
  window_end( t_end ),
  energy_fraction_sample( E_extra_sample )
{
  
  cout << shower_scattered_pdf << endl;
  cout << shower_direct_pdf << endl;
  cout << sigma_blur << endl;
  cout << elong_steps << endl;
  cout << R_bg << endl;
  cout << t_start << endl;
  cout << t_end << endl;
  cout << E_extra_sample << endl;

  Timer tim("constructPDF", true);

  _pdf_funcs.resize(3);
    
  print ("loading scattered photon table.");
  _pdf_funcs[1].load( shower_scattered_pdf.c_str() ) ;
    
  print ("loading direct photon table.");
  _pdf_funcs[2].load( shower_direct_pdf.c_str() );
    
  print("adding direct and scattered tables.");
  _pdf_funcs[0] = _pdf_funcs[1];
    
  for( auto& p : _pdf_funcs ) 
    p.setExceptionHandler( new typename PDF_Function1D_t::JDefaultResult( JResultPDF<double>( 0.0, 0.0, 0.0, 0.0 ) ) );
    
  add_jpp_shower_pdfs( _pdf_funcs[0], _pdf_funcs[2]        );
    
  print("computing table of total npe.");
  _npe_funcs.resize(3);
  for( int i=0; i<3; i++) {
    _npe_funcs[i] = JCascadeNpe_t( _pdf_funcs[i] );
    _npe_funcs[i].setExceptionHandler( new JFunction<double,double>::JDefaultResult( JMATH::zero ) );
  }
    
  print("smearing combined table with sigma=",sigma_blur );
  if ( sigma_blur > 0 ){
    _pdf_funcs[0].blur( sigma_blur ); // only blur total : save some time
  }
    
 
  print("using", elong_steps, "elongation steps");

  print("background rate:", R_bg, "Hz");
  print("time window:", t_start, t_end );
  print("energy fraction sample", energy_fraction_sample);

  print("CascadePdf constructed in", tim.value(), "sec");

  set_light_type( light_type );
} 


JCascadePdf_t::result_type CascadePdf::eval1(const Paramst& p, JCascadePdf_t* f ) 
{ 
  JCascadePdf_t::result_type r;

  {
    Timer _("evalpdf");
    r =  (*f)( p.d, p.z, p.theta, p.phi, p.t ) * p.E ;
  }

  check(r.v, p );
  check(r.f, p );
  return r;
}

JCascadeNpe_t::result_type CascadePdf::eval1(const Params & p, JCascadeNpe_t* f ) 
{ 
  Timer _("evalnpe");
  return check( (*f)( p.d, p.z, p.theta, p.phi      ) * p.E , p );
}

JCascadePdf_t::result_type CascadePdf::eval1_opt(const Paramst& p, JCascadePdf_t* f ) 
{ 
  return (*f)( p.d, p.z, p.theta, p.phi, p.t ) * p.E ;
}

JCascadeNpe_t::result_type CascadePdf::eval1_opt(const Params & p, JCascadeNpe_t* f ) 
{ 
  return (*f)( p.d, p.z, p.theta, p.phi ) * p.E;
}

JCascadePdf_t::result_type CascadePdf::eval1(double d, double z, double theta, double phi, double t, double E ) 
{ 
  Timer _("evalpdf");
  return (*_pdf_func)( d, z, theta, phi, t ) * E ;
}


JCascadeNpe_t::result_type CascadePdf::eval1(double d, double z, double theta, double phi, double E ) 
{ 
  Timer _("evalnpe");
  return (*_npe_func)( d, z, theta, phi ) * E;
}

/*
  JCascadePdf_t::result_type CascadePdf::deval1(const Paramst& p, JppGrad* f ) 
  { 
  Timer _("evalgrad");
  return (*f).dnpe_dXs( p.d, p.z, p.theta, p.phi ) * p.E ;
  }

  JCascadePdf_t::result_type CascadePdf::ddeval1(const Paramst& p, JppGrad* f ) 
  { 
  Timer _("evalgrad2nd");
  return (*f).dnpe_dtdXs( p.d, p.z, p.theta, p.phi, p.t ) * p.E ;
  }
*/