#ifndef __GRADTOOLS__
#define __GRADTOOLS__

#include "Params.hh"
#include "Trk.hh"
#include "CascadePdf.hh"

std::pair< Trk, Trk > delta_trks( const Trk &trk, uint ipar, double delta );
std::pair< Trk, Trk > delta_trks_e_frac( const Trk &trk, const double fraction, uint ipar, double delta );
std::pair< Trk, Trk > delta_trks_e( const Trk &trk, const double len, uint ipar, double delta );

std::pair< Paramst, Paramst > delta_pst( const Paramst &pst, uint ipar, double delta );
std::pair< Params, Params > delta_ps( const Params &ps, uint ipar, double delta );

vector< vector< double > > num_dX_dTs( const Trk &trk, const Pmt &unhit, double delta);
vector< vector< double > > num_dX_dTs( const Trk &trk, const Hit &hit, double delta);
vector< vector< double > > num_dTe_dTs( const Trk &trk, const double len, double delta);

int get_trk_hit_hash( Trk &trk, Hit &hit );
int get_trk_unhit_hash( Trk &trk, Pmt &unhit );

//template< typename pdf_type>
//vector<double> num_dnpe_dtdXs( pdf_type *pdf, double D, double cd, double theta, double phi, double t, double E, double delta );

class numFuncs
{

public:

  numFuncs(){};

  template< typename pdf_type>
  vector<double> num_dnpe_dtdXs( pdf_type *pdf, double E, double D, double cd, double theta, double phi, double t, double delta )
  {
    vector<double> ds( 6, 0 );

    if( E < delta ) E = delta;
    if( D < delta ) D = delta;
    if( cd >  1.0 - delta ) cd =  1.0 - delta;
    if( cd < -1.0 + delta ) cd = -1.0 + delta;

    const vector< double > pars = { E, D, cd, theta, phi, t };

    // E, x, y, z, dphi, dtheta, t
    for( uint ii = 0; ii != 6; ii++ ){

      vector< double > pars_l = pars, pars_h = pars;
      pars_l[ii] -= delta/2.;
      pars_h[ii] += delta/2.;
      
      //cout << "ii : " << ii << endl;
      //cout << std::setprecision(10) << "pars_l: 0 " << pars_l[0] << ", 1 " << pars_l[1] << ", 2 " << pars_l[2] << ", 3 " << pars_l[3] << ", 4 " << pars_l[4] << ", 5 " << pars_l[5] << endl;
      //cout << std::setprecision(10) << "pars_h: 0  " << pars_h[0] << ", 1 " << pars_h[1] << ", 2 " << pars_h[2] << ", 3 " << pars_h[3] << ", 4 " << pars_h[4] << ", 5 " << pars_h[5] << endl;

      ds[ii] = ( pdf->dnpe_dt( Paramst( pars_h[0],
					pars_h[1],
					pars_h[2],
					pars_h[3],
					pars_h[4], 
					pars_h[5] ) )
		 - pdf->dnpe_dt( Paramst( pars_l[0],
					  pars_l[1], 
					  pars_l[2], 
					  pars_l[3], 
					  pars_l[4],
					  pars_l[5] ) ) ) / delta;
      //cout << std::setprecision(10) << "ds[ii]: " << ds[ii] << endl;
      //cin.ignore();
    }

    return ds;

  };

  template< typename pdf_type>
  vector<double> num_dnpe_dtdXs_ps( pdf_type *pdf, Paramst ps, double delta ){
    return num_dnpe_dtdXs( pdf, ps.E, ps.d, ps.z, ps.theta, ps.phi, ps.t, delta );
  };

  //template< typename pdf_type>
  //vector<double> num_dnpe_dXs( pdf_type *pdf, double E, double D, double cd, double theta, double phi, double t, double delta );

  //template< typename pdf_type>
  //vector<double> num_dnpe_dXs_ps( pdf_type *pdf, Paramst ps, double delta );

  //template< typename pdf_type>
  //vector<double> num_dnpe_dXs( pdf_type *pdf,double E,  double D, double cd, double theta, double phi, double delta );

  //template< typename pdf_type>
  //vector<double> num_dnpe_dXs_ps( pdf_type *pdf, Params ps, double delta );

  template< typename pdf_type>
  vector<double> num_dnpe_dXs( pdf_type *pdf, double E,  double D, double cd, double theta, double phi, double delta )
  {
    vector<double> ds( 6, 0 );

    if( E < delta ) E = delta;
    if( D < delta ) D = delta;
    if( cd >  1.0 - delta ) cd =  1.0 - delta;
    if( cd < -1.0 + delta ) cd = -1.0 + delta;

    const vector< double > pars = { E, D, cd, theta, phi };

    // x, y, z, dphi, dtheta, E
    for( uint ii = 0; ii != 5; ii++ ){	
      vector< double > pars_l = pars, pars_h = pars;
      pars_l[ii] -= delta/2.;
      pars_h[ii] += delta/2.;

      ds[ii] = ( pdf->npe( Params( pars_h[0],
				   pars_h[1],
				   pars_h[2], 
				   pars_h[3], 
				   pars_h[4] ) )
		 - pdf->npe( Params( pars_l[0],
				     pars_l[1], 
				     pars_l[2], 
				     pars_l[3], 
				     pars_l[4] ) ) ) / delta;
    }

    return ds;

  };

  //template< typename pdf_type>
  //vector<double> num_dnpe_dXs( pdf_type *pdf, double E, double D, double cd, double theta, double phi, double t, double delta );


  //extern template<> vector<double> num_dnpe_dXs<CascadePdf>( CascadePdf *pdf, double E, double D, double cd, double theta, double phi, double t, double delta );

  template< typename pdf_type>
  vector<double> num_dnpe_dXs( pdf_type *pdf, double E, double D, double cd, double theta, double phi, double t, double delta )
  {
    vector<double> ds( 6, 0 );

    if( E < delta ) E = delta;
    if( D < delta ) D = delta;
    if( cd >  1.0 - delta ) cd =  1.0 - delta;
    if( cd < -1.0 + delta ) cd = -1.0 + delta;

    const vector< double > pars = { E, D, cd, theta, phi, t };

    // E, D, cd, theta, phi, t
    for( uint ii = 0; ii != 6; ii++ ){	

      vector< double > pars_l = pars, pars_h = pars;
      pars_l[ii] -= delta/2.;
      pars_h[ii] += delta/2.;
      ds[ii] = ( pdf->npe( Paramst( pars_h[0],
				    pars_h[1], 
				    pars_h[2], 
				    pars_h[3],
				    pars_h[4],
				    pars_h[5] ) )
		 - pdf->npe( Paramst( pars_l[0],
				      pars_l[1],
				      pars_l[2],
				      pars_l[3], 
				      pars_l[4],
				      pars_l[5] ) ) ) / delta;
    }

    return ds;

  };

  template< typename pdf_type>
  vector<double> num_dnpe_dtdXs( pdf_type *pdf, Paramst ps, double delta ){
    return num_dnpe_dtdXs( pdf, ps.E, ps.d, ps.z, ps.theta, ps.phi, ps.t, delta );
  }

  template< typename pdf_type>
  vector<double> num_dnpe_dXs( pdf_type *pdf, Paramst ps, double delta ){
    return num_dnpe_dXs( pdf, ps.E, ps.d, ps.z, ps.theta, ps.phi, ps.t, delta );
  }

  template< typename pdf_type>
  vector<double> num_dnpe_dXs( pdf_type *pdf, Params ps, double delta ){
    return num_dnpe_dXs( pdf, ps.E, ps.d, ps.z, ps.theta, ps.phi, delta );
  }
  
};

#endif