#ifndef PARAMSHH_INCLUDED
#define PARAMSHH_INCLUDED

#include "constants.hh" 

#include "../utilJpp.hh"
#include "Hit.hh"
#include "Vec.hh"
#include "io_util.hh"

//! Paramaters describing shower-pmt relation
struct Params
{
  double E, d, z, theta, phi;

  Params(): 
    E(0.0), d(0.0), z(0.0), theta(0.0), phi(0.0) {}

  Params(double E, double d, double z, double theta, double phi) :
    E(E), d(d), z(z), theta(theta), phi(phi) {}

  void init( const Trk& trk, const Vec& pos, const Vec& dir )
  {
    // trk.dir.normalize();
    // dir.normalize();

    const Vec    v = pos - trk.pos;
    const double y = v.dot( trk.dir );              
    double u = dir.dot( trk.dir );

    if ( u > 1.0 ) u = 1.0;
    if ( u < -1.0 ) u = -1.0;

    if ( trk.E > 1e14 )   E = 1e14;
    else if ( trk.E < 1 ) E = 1;
    else                  E = trk.E;

    d = v.len();
    if ( d == 0.0 ) z = 0.0; //Pick arbitrary direction
    else            z = y/d;
    if ( z >  1.0 ) z =  1.0;
    if ( z < -1.0 ) z = -1.0;
    theta = acos(u);
    phi   = angle_between ( v - y * trk.dir, dir - u * trk.dir );
    if ( phi != phi ) phi = 0.0;
    assert( E == E );
    assert( d == d );
    assert( z == z );
    assert( phi == phi );
    if( theta != theta ){
      cout << "u: " << u << endl;
      cout << trk << endl;
      cout << pos << endl;
      cout << dir << endl;
    }
    assert( theta == theta );

  }
  Params( const Trk& trk, const Vec& pos, const Vec& dir ) { init( trk, pos, dir ); }
  Params( const Trk& trk, const Hit& hit ) { init( trk, hit.pos, hit.dir ); }
  Params( const Trk& trk, const Pmt& pmt ) { init( trk, pmt.pos, pmt.dir ); }

  // copy
  Params( const Params& p2 ) { E = p2.E; d = p2.d; z = p2.z; theta = p2.theta; phi = p2.phi; }

  //! compute parameters corresponding to moving 'length' along the shower direction
  Params propagate( double length ) const
  {
    Params r = *this;
    const double Z0 = z * d;
    const double X2 = d*d - Z0*Z0; // X remains unchanged
    const double Z1 = Z0 - length;
    r.d = sqrt(Z1*Z1 + X2);
    r.z = Z1 / r.d;

    if ( r.z > 1.0 ) r.z = 1.0;
    if ( r.z < -1.0 ) r.z = -1.0;

    return r;
  }


  void print( std::ostream& out = std::cout ) const 
  {
    out << " E=" << E << " d=" << d << " z=" << z << " theta="<< theta << " phi=" << phi << endl;
  }

};



//template Params::Params(const Trk&, const Pmt&); // template instantiation

ADD_PRINTING( Params );

//! Parameters, including time (residual) 

struct Paramst : public Params
{
  double t;

  Paramst( const Params& p, double t ) : Params(p), t(t) {}

  Paramst() :
    Params(), t(0.0) {}

  Paramst(double E, double d, double z, double theta, double phi , double t) :
    Params( E,d,z,theta,phi), t(t) {}

  Paramst( const Trk& trk, const Hit& hit ) : 
    Params( trk, hit ), t( hit.t - ( trk.t + d / aa::v_light ) ) {}

  void init_t( const Trk &trk, const Hit &hit ){
    init( trk, hit.pos, hit.dir );
    t = hit.t - ( trk.t + d / aa::v_light );
  };

  // copy
  Paramst( const Paramst& p2 ) : Paramst( p2.E, p2.d, p2.z, p2.theta, p2.phi, p2.t ) {}

  Paramst propagate( double length ) const 
  {
    Params r = Params::propagate( length ); // slicing
    return Paramst( r, t - length / aa::c_light - ( r.d - d ) / aa::v_light ) ;
  }

  void print( std::ostream& out = std::cout ) const 
  {
    Params::print(); out << " t=" << t ;
  }

};

ADD_PRINTING( Paramst );

Trk propagate( Trk& trk, double length );

Trk get_elong_trk( const Trk &trk, const double fraction, const double precision = 1e-6);

#endif