#ifndef TRK_UTIL_HH_INCLUDED
#define TRK_UTIL_HH_INCLUDED

#include "Trk.hh"
#include "constants.hh"



 /*! \brief  Compute complete information for a Cherenkov photon 
    emited from the Trk and hitting a PMT. 

    \param trk         track.
    \param p           input: position of the PMT       				   
    \param w		       input direction of the PMT 					   
    \param time	       output: arrival time of the photon (i.e. hit-time)		   
    \param d_track	   output: distance along the track where the photon is emmited	   
    \param d_photon	   output: distance of the photon (from track to pmt)		   
    \param d_closest	 output: distance between the PMT and the track (closest approach)
    \param cos_angle	 output: cos of angle of impact of photon wrt PMT direction	   
    \param phot_dir 	 if the pointer is non-zero, the photon direction is returned.    
  */ 

inline void cherenkov_pars( const Trk& trk, 
  const Vec& p,              
  const Vec& w,              
  double& time,              
  double& d_track,           
  double& d_photon,          
  double& d_closest,         
  double& cos_angle,         
  Vec* phot_dir=0 )
{ 
  const Vec v = p - trk.pos;
  const double l = v.dot(trk.dir);
  double k2 = v.dot(v) - l*l;
  d_closest = k2 > 0 ? sqrt(k2) : 0 ;

  d_photon  = d_closest/ sin(cherenkov_angle );
  d_track   = l - d_closest / tan( cherenkov_angle );
  time      = trk.t + d_track / c_light + d_photon /v_light;

  Vec vphot =  v - trk.dir* d_track;
  vphot /= vphot.len();

  cos_angle = vphot.dot( w ) ;

  if (phot_dir)  *phot_dir = vphot;
}

 /*! Reset both the Track time and position so that 
    the position is as close as possible to p */

inline void move_trk_to( Trk& trk, const Vec& p) 
{
  const Vec v    = p - trk.pos;
  const double l = v.dot(trk.dir);
  trk.pos += trk.dir * l;
  trk.t   += l / c_light ;
}

  /*! Return the time the Cherenkov light reaches point p */

inline double time_at(const Trk& trk, const Vec& p)
{
  const Vec v = p - trk.pos;
  const double l = v.dot(trk.dir);
  double k2 = v.dot(v) - l*l;
  if (k2<0) k2=0;

  const double d_closest = sqrt(k2);
  const double d_photon  = d_closest / sin( cherenkov_angle );
  const double d_track   = l - d_closest / tan ( cherenkov_angle );
  return trk.t + d_track / c_light + d_photon /v_light;
}

#endif