#pragma once

#include "constants.hh" 

#include "Hit.hh"
#include "Evt.hh"
#include "Det.hh"
#include "Spline.hh"
#include "Mat.hh"
#include "io_util.hh"

#include "JPhysics/JGeanz.hh"
//#include "JGeano.hh"

#include "JAAnet/JAAnetToolkit.hh"
#include "JGeometry3D/JRotation3D.hh"
#include "JGeometry3D/JPosition3D.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JGeometry3D/JAxis3D.hh"

#include "TASImage.h"
// def filter_hit_residuals ( hit_vector, trk, pred = lambda x : -100<x<900 ) :
    
//     r = ROOT.vector(ROOT.Hit)()
//     nk =0
//     for h in hit_vector :
//         res = h.t - trk.t - (trk.pos-h.pos).len() / v_light;
//         ok = pred(res)
//         if ok : r.push_back( h )

//     return r

using namespace JPP;

inline Hit hit_from_pmt( Pmt& pmt ){
  Hit h;
  pmt.dress( h );
  return h;
}

template<typename T>
void load_jpp_shower_pdf( T& pdf , string fn)
{
  pdf.load( fn.c_str() );
  //pdf.compile();

  cout << "loaded " << fn << endl;
}

template<typename T>
void add_jpp_shower_pdfs( T& pdf , std::vector<string> fn )
{
  for(std::vector<string>::const_iterator ifn = fn.begin(); ifn != fn.end(); ifn ++){

    T temp_pdf;
    load_jpp_shower_pdf<T::function_type>(temp_pdf, *ifn);
    
    try{ 
     pdf.add(temp_pdf);
    } catch(JLANG::JException e ){
      cout << "Failed adding: " << e << endl;
      throw e;
    }
    cout << "Added " << *ifn << endl;
  }
}

template<typename T>
void add_jpp_shower_pdfs( T& pdf0, T& pdf1 )
{
  
  try{
    pdf0.add(pdf1);
  } catch(JLANG::JException e ){
    cout << "Failed adding, " << e << endl;
    throw e;
  }
  cout << "Added " << endl;
  //Set to something else?
  //pdf0.setExceptionHandler( new JLANG::JBool<false> );
}

template<typename T>
void load_jpp_track_pdf( T& pdf , string fn)
{
  pdf.load( fn.c_str() );
  cout << "loaded " << fn << endl;
}


template<typename T>
void add_jpp_track_pdfs( T& pdf0, T& pdf1 )
{
  
  try{
    pdf0.add(pdf1);
  } catch(JLANG::JException e ){
    cout << "Failed adding, " << e << endl;
    throw e;
  }
  cout << "Added " << endl;
  //Set to something else?
  //pdf0.setExceptionHandler( new JLANG::JBool<false> );
}


/*!
* This function returns a parametrisation of the shower elongation
* for a number of steps. 
* geanz.getLength returns shower length for a given integrated probability.
*/
inline vector<double> getShowerDepths(const double E, const uint elong_steps) {
  vector<double> Z;
  
  for (uint i = 0; i != elong_steps; ++i) {
    Z.push_back( geanz.getLength( E, (i + 0.5) / (double) elong_steps ));  
  }
  return Z;

}

/*!
* This function returns the depth of the shower maximum
*/
inline double getShowerMaximum(const double E) {
  const double depth = geanz.getMaximum( E );
  return depth;
}

/*!
* This function returns the depth of the shower
*/
inline double getShowerLength(const double E, const double rn) {
  const double depth = geanz.getLength( E, rn );  
  return depth;
}

inline void trk_to_smax_pars( Trk& t, double* r )
{
  double smax_d = JPP::geanz.getLength( t.E, 0.5 );
  Vec pos_smax = t.pos + t.dir * smax_d;
  r[0] = pos_smax.x; r[1] = pos_smax.y; r[2] = pos_smax.z;
  r[3] = acos( t.dir.z );
  double phi = atan2( t.dir.y, t.dir.x );
  while( phi < 0 )    phi += 2*pi;
  while( phi > 2*pi ) phi -= 2*pi;
  r[4] = phi;
  r[5] = t.t + smax_d / c_light;
  r[6] = log10(t.E);
}
//https://www.osti.gov/pages/servlets/purl/1330371
//Particle showers evolve at ~c
inline void smax_pars_to_trk( const double* r, Trk& t )
{
  t.dir.set_angles( r[3],r[4] );
  try{
    t.E = (double)pow(10., r[6]); //Note that the power should be a double
  } catch( string except ){
    std::cout << "No energy: " << except << std::endl;
    t.E = -1;
  }
  double smax_d = JPP::geanz.getLength( t.E, 0.5 );
  t.pos = Vec(r[0], r[1], r[2]) - t.dir * smax_d;
  t.t   = r[5] - smax_d / c_light;
}

/*! Struct to compute parameters that describe relation between shower/track and a hit */
struct HitParams
{
  
  double t1;
  double dt;
  double D;
  double R;
  double cd;
  double theta;
  double phi;
  Vec    _trk_pos;
  Vec    _hit_pos;
  double t;
  
  HitParams(){};
  
  //Position constructor
  HitParams( const Vec& trk_pos ):
    _trk_pos( trk_pos ){}

  //Position + unhit constructor    
  HitParams( const Vec&    trk_pos,
	     const Vec&    hit_pos ):
    HitParams( trk_pos )
  {
    set_hit_pos( hit_pos );
  } 

  //Full unhit constructor
  HitParams( const Trk& track,
	     const Vec& hit_pos,
	     const Vec& hit_dir ):
    HitParams( track.pos,
	       hit_pos )
  {
    set_hit_dir( JPosition3D( _trk_pos.x, _trk_pos.y, _trk_pos.z ), getDirection( track ), hit_dir );
  }
  

  
  //Position + time constructor
  HitParams( const Vec&    trk_pos,
	     const double& trk_t):
    _trk_pos( trk_pos ),
    t  ( trk_t ){}
  
  //Position + hit constructor    
  HitParams( const Vec&    trk_pos,
	     const double& trk_t,
	     const Vec&    hit_pos,
	     const double& hit_t):
    HitParams( trk_pos, trk_t)
  {
    set_hit_pos( hit_pos );
    set_hit_t( hit_t );
  } 
  
  //Full hit constructor
  HitParams( const Trk&    track,
	     const Vec&    hit_pos,
	     const Vec&    hit_dir,
	     const double& hit_t):
    HitParams( track.pos,
	       track.t,
	       hit_pos,
	       hit_t   )
  {
    set_hit_dir( getPosition( track ), getDirection( track ), hit_dir );
  }

  void set_hit_pos( const Vec& hit_pos ){
    _hit_pos = hit_pos;
    D = ( _trk_pos - _hit_pos ).len();
  };
  
  void set_hit_t( const double& hit_t ){
    t1    = t + D / aa::v_light; //* getInverseSpeedOfLight() * getIndexOfRefraction();
    dt    = hit_t - t1;
  };
  
  void set_hit_dir( JPosition3D Jpos, const JDirection3D& Jdir, const Vec& hit_dir ){
    
    const JRotation3D Rot( Jdir );
    
    Jpos.rotate(Rot);
    
    JAxis3D axis( JPosition3D(  _hit_pos.x, _hit_pos.y, _hit_pos.z ),
		  JDirection3D( hit_dir.x, hit_dir.y, hit_dir.z ) );
    axis.transform( Rot, Jpos);
    
    R     = axis.getPosition().getX();  
    cd    = axis.getZ()/D;
    theta = axis.getTheta();
    phi   = fabs(axis.getPhi());
  }
  
  HitParams( const Vec&    trk_pos,
	     const double& trk_t,
	     const Hit&    hit ):
    HitParams( trk_pos,
	       trk_t,
	       Vec(hit.pos.x, hit.pos.y, hit.pos.z),
	       hit.t){}
  
  HitParams( const Trk& track,
	     const Hit& hit ):
    HitParams( track,
	       Vec(hit.pos.x, hit.pos.y, hit.pos.z),
	       Vec(hit.dir.x, hit.dir.y, hit.dir.z),
	       hit.t){}
  
  void print( std::ostream& out = std::cout ) const 
  {
    out << "D=" << D << " cd=" << cd << " theta="<< theta << " phi=" << phi << " dt=" << dt;
  }

  void dump(ostream& out = cout )
  {
    out << " Hit params object " << endl;
    out << " t1    " << t1 << endl;
    out << " dt    " << dt << endl;
    out << " R     " << R << endl;
    out << " D     " << D << endl;
    out << " cd    " << cd << endl;
    out << " theta " << theta << endl;
    out << " phi   " << phi << endl;
  }
  
  double getShowerDt()
  {
    return dt;
  }
  
  double getTrackDt()
  {
    //Change this to something correct
    //t1    = track.t + D * getInverseSpeedOfLight() * getIndexOfRefraction();
    return -9999; //hit_time - t1;
  }
  
};

ADD_PRINTING( HitParams );

inline const double get_shifted_D( double D, double cd, double shift ){
  return sqrt( D*D - 2.0*D*cd*shift + shift*shift );
};
inline const double sD_get_shifted_cd( double D, double shifted_D, double cd, double shift ){
  return ( D * cd - shift ) / shifted_D;
};
inline const double get_shifted_cd( double D, double cd, double shift ){
  return sD_get_shifted_cd( D, get_shifted_D( D, cd, shift ), cd, shift );
};
inline const double sD_get_shifted_dt( double D, double shifted_D, double dt, double shift ){
  return dt + ( D - shifted_D )  / aa::v_light - shift / aa::c_light;
};
inline const double get_shifted_dt( double D, double cd, double dt, double shift ){
  return sD_get_shifted_dt( D, get_shifted_D( D, cd, shift ), dt, shift );
};


inline const double get_unshifted_D( double shifted_D, double shifted_cd, double shift ){
  return sqrt( shifted_D*shifted_D + 2.0*shifted_D*shifted_cd*shift + shift*shift );
};
inline const double get_unshifted_cd( double shifted_D, double shifted_cd, double shift ){
  return ( shifted_D*shifted_cd + shift ) / get_unshifted_D( shifted_D, shifted_cd, shift );
};
inline const double get_unshifted_dt( double shifted_D, double shifted_cd, double shifted_dt, double shift ){
  return shifted_dt + ( shifted_D - get_unshifted_D( shifted_D, shifted_cd, shift ) )  / aa::v_light + shift / aa::c_light;
};


inline const double get_shower_max_D(double D, double cd, double E){
  return get_shifted_D( D, cd, geanz.getMaximum( E ) );
};
inline const double get_shower_max_cd(double D, double cd, double E){
  return get_shifted_cd( D, cd, geanz.getMaximum( E ) ); 
};

inline const double get_vertex_D(double shifted_D, double shifted_cd, double E){
  return         get_unshifted_D(       shifted_D,        shifted_cd, geanz.getMaximum( E ) );
};
inline const double get_vertex_cd(double shifted_D, double shifted_cd, double E){
  return         get_unshifted_cd(       shifted_D,        shifted_cd, geanz.getMaximum( E ) ); 
};
inline const double get_vertex_dt(double shifted_D, double shifted_cd, double shifted_dt, double E){
  return         get_unshifted_dt(       shifted_D,        shifted_cd,        shifted_dt, geanz.getMaximum( E ) ); 
};

inline const long long int emerge( double num,
				   uint decimal_place = 10){
  return std::floor( std::pow( 10, decimal_place ) * num );
};

inline const bool are_emerged_equal( double num1,
				     double num2,
				     uint decimal_place = 10,
				     bool ignore_small = true){
  if( ignore_small && log10( abs( num1 ) ) < -decimal_place && log10( abs( num2 ) ) < -decimal_place ){
    return true;
  }
  const int scale = std::pow( 10, -std::floor( log10( abs( num1 ) ) ) );
  return abs( emerge( scale * num1, decimal_place ) - emerge( scale * num2, decimal_place ) ) <= 1;
};

//Translate aanet angles to jpp angles for given irad
inline const vector< double > get_pmt_angles ( const double cd,
					       const double pmt_a,
					       const double irad ){
  
  //Vec dir = Vec( sqrt(( 1 - cd )*( 1 + cd )), 0, cd ).normalize();
  
  Vec dir = Vec( sin( acos( cd ) ), 0, cd );
  //Must be normalized because we would occasionally get rot_dir = { x, 0, y } where x is very small, and y is very slightly under -1.0. 
  const Vec rot_dir = ( Mat( dir ) * Vec(0, 0, 1).rotate_y( pmt_a ).rotate_z( irad ) ).normalize();
  
  return { rot_dir.theta(), rot_dir.phi() };
};



/*! Structure to store information 
    needed to compute the likelihood faster. */

struct ExtDomSum
{
  Vec pos ;                 // dom position
  vector< double > hit_ls;   
  HitParams pre_hp;
};


//Functions to test above

/*Distance from track vertex to hit*/
inline double calc_D(const Vec trk_pos,
		     const Vec hit_pos)
{
  return (trk_pos - hit_pos).len();
}

/*Length of projection of vertex hit vector on track axis*/
inline double calc_L(const Vec trk_pos,
		     const Vec hit_pos,
		     const Vec trk_dir)
{
  return trk_dir.dot(hit_pos - trk_pos);
}

/*Perpendicular distance from track axis to hit*/
inline double calc_R(const Vec trk_pos,
		     const Vec hit_pos,
		           Vec trk_dir)
{
  const Vec track_hit_perpendicular = hit_pos - trk_pos - trk_dir.normalize()*calc_L(trk_pos, hit_pos, trk_dir);
  return track_hit_perpendicular.len();
}

/*Angle of hit off direction of track*/
inline double calc_cd(const Vec trk_pos,
		      const Vec hit_pos,
		      const Vec trk_dir)
{
  return trk_dir.dot(hit_pos - trk_pos)/calc_D(trk_pos, hit_pos); 
}

/*Difference of expected light arrival time and hit time*/
inline double calc_dt(const Vec trk_pos,
		      const Vec hit_pos,
		      const double trk_t,
		      const double hit_t)
{
  return hit_t - trk_t - calc_D(trk_pos, hit_pos) / aa::v_light;
}

/*Difference of expected light arrival time and hit time*/
inline double calc_dt(const Trk& trk,
		      const Hit& hit)
{
  return calc_dt(trk.pos, hit.pos, trk.t, hit.t);
}

struct PyFuncs{

  PyFuncs(){};

  /*!
   * This function returns the depth of the shower
   */
  inline double getShowerLength( const double E, const double rn ) {
    const double depth = geanz.getLength( E, rn, 1e-6);  
    return depth;
  };

  inline double getShowerMax( const double E ) {
    return geanz.getMaximum( E );
  };
  /* // under construction...
     vector<char> getTASImageBuffer( TImage image ){
     
     char** buffer;
     int*   size;
     image.GetImageBuffer( buffer, size, TImage::kPng );
     
     vector<char> out(*size);
     for( uint ii = 0; ii != *size; ii++ ){
      out[ii] = buffer[ii];
      }
      return out;
      };
  */
};


class testclass
{
  testclass(){};
};

  
// inline ostream& operator<<(ostream& out, const Cherenkov_pars& p ) {}

inline bool sort_t( const Hit& h1, const Hit& h2 ){
  return h1.t < h2.t;
}