#ifndef AHRECUTILINCD
#define AHRECUTILINCD

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

#include <unordered_map>

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

struct DomSum
{
  uint id ;
  Vec pos ;               /*! dom position */

  vector< std::pair<uint, Vec> >  hit_pmts_id;  
  vector< std::pair<uint, Vec> >  unhit_pmts_id;

  vector< Vec >  hit_pmts;   /*! directions of pmts with at least a hit */
  vector< Vec >  unhit_pmts; /*! directions of the unhit pmts */

  DomSum() { //cout << "domsum ctor" << endl; 
  }
  ~DomSum()  { //cout  << "domsum dtor" << endl; 
  }

};

/*! Structure to store precomputed information for very
    fast evaluation of the likelihood vs energy and direction
    for a given position */

struct Fdom
{
  Fdom( const TH1D& h ) : hist0(h) {}
  Fdom() { //{cout << "fdom ctor" << endl;
  }

  Vec dir;                 // normalized vector from shower to dom
  TH1D hist0;              // sum of z-slices (mu vs z) for unhit pmts 
  Spline spline0;

  vector<TH1D> hit_hists;  // z-slices for the hit pmts
  vector<Spline> hit_splines;


  ~Fdom()  { //cout  << "fdom dtor" << endl; 
  }

};

inline bool first_sort (pair< double, double > i, pair< double, double > j) { return (i.first<j.first); }

inline double poisson(const double k, const double lambda, const double softening = 1e-12)
{
  return exp(k * log(lambda + softening) - lgamma(k + 1.0) - lambda);
};

inline vector<Hit> filter_hits( const vector<Hit>& input,
				const Trk& trk,
				double min_residual  = -100,
				double max_residual  = 900,
				double max_distance  = 1e100,
				double min_hit_angle = 0 )
{
  vector<Hit> result;

  foreach( h, input )
    {
      const double l = (trk.pos-h.pos).len();
      if ( l > max_distance ) continue;

      const double r = h.t - trk.t - l / v_light;
      if ( r < min_residual || r > max_residual ) continue; 

      const double a = angle_between( h.dir, h.pos - trk.pos );
      if ( a < min_hit_angle ) continue; 

      result.push_back( h );
    }

  return result;
}

/*! This function select all hits within a residuals window and a certain distance. It also orders the hits with
 *  by PMT id so it can also return the first hit on every PMT. */
inline vector<Hit> select_hits( const vector<Hit>& hits,
				const Trk& trk,
				double min_residual  = -100,
				double max_residual  = 900,
				double max_distance  = 1e100,
				bool first_only = false )
{

  vector<Hit> result;

  std::map<int, vector<const Hit*> > M;

  for( auto& h : hits ) {

    const double l = (trk.pos - h.pos).len();
    if ( l > max_distance ) continue;

    const double r = h.t - trk.t - l / v_light;
    if ( r < min_residual || r > max_residual ) continue; 

    M[ h.pmt_id ].push_back( &h );
  }


  foreach_map( pmt_id, hitlist , M )
  {
    (void) pmt_id; // stop warnings about unused map key
    std::sort( hitlist.begin(), hitlist.end() , [] (const Hit* a, const Hit* b) { return a->t < b->t; } );

    foreach( h, hitlist ) {
      result.push_back( *h );
      if ( first_only ) break;
    }

  }
  return result;

}

inline double log_without_error( double p ) {

  double logp = 0;
  if ( p < 1e-14 ) { // check for extremely small probability; this gives infinite likelihoods 
    logp = log( 1e-14 ); 
  }
  else {
    logp = log( p ); 
  }

  return logp;

}

inline bool earlier( const Hit* a, const Hit* b )
{
  return a->t < b->t;
}

/*! merge a list of hits. Keep adding hits together as long the 
    distance to the next hits is < gate */

inline vector<Hit*> greedy_merge( vector<Hit*>& v , double gate = 300 )
{
   if ( v.size() < 2 ) return v;
   vector<Hit*> r;
   
   for(unsigned i=0; i< v.size(); )
   {
     double T = v[i]->t + gate;
  
     for ( unsigned j= i+1 ; j < v.size()+1 ; ++j )
        {
           if ( j==v.size() || v[j]->t > T )
             {
	       r.push_back( v[i] );
	       i=j;
	       break;
             }
           else // add j to i
             {
               v[i]->a += v[j]->a;
               T = v[j]->t + gate;
             }  
        }
   }
   return r;
}

/*! merge hits on the same PMT (or DOM) within gate (ns) */

inline vector<Hit> merge_hits(vector<Hit>& hits, double gate = 300 , bool merge_on_doms = false ) 
{
  vector<Hit> r;

  map<int, vector<Hit*> > m;

  if (merge_on_doms) foreach( h, hits ) m[h.dom_id].push_back(&h);
  else               foreach( h, hits ) m[h.pmt_id].push_back(&h);

  foreach_map( k, v, m )   {
    (void) k;
    std::sort( v.begin(), v.end(), earlier ); 
    vector<Hit*> y = greedy_merge( v , gate );
    foreach( h, y ) r.push_back( *h ) ;
  }  

  return r;
}

inline vector<Hit> get_coincidences(vector<Hit> &hits,
                                    int level = 2,
                                    double dt = 20)
{
  vector<Hit> r;

  std::unordered_map<int, vector<Hit *>> m;

  foreach (h, hits)
    m[h.dom_id].push_back(&h);

  foreach_map(domid, hits_on_dom, m)
  {
    (void)domid;

    if (hits_on_dom.size() < 2)
      continue;

    std::sort(hits_on_dom.begin(), hits_on_dom.end(), earlier);

    unsigned int i = 0;

    while (i < hits_on_dom.size())
    {
      unsigned int j = i + 1;

      while (j < hits_on_dom.size() &&
             hits_on_dom[j]->t - hits_on_dom[i]->t < dt) ++j;
      // j points to the first hit not in coincidence with i

      if (j > i + level - 1)
      {
        Hit h(*hits_on_dom[i]);

        for (unsigned int k = i + 1; k < j; ++k)
        {
          h.a   += hits_on_dom[k]->a;
          h.tot += hits_on_dom[k]->tot;
          if (hits_on_dom[k]->type == -1)
            h.type = -1;
        }

        r.push_back(h);
        i = j;
      }
      else
      {
        i++;
      }
    }
  }

  return r;
}

inline double hit_to_dom_ratio( Det& det, 
				vector<Hit>& hits,
				Vec pos,
				double radius = 200 )
{
  int ndoms =0;
  int nhits =0;
  
  foreach_map( domid, dom, det.doms )
    {
      (void) domid;
      if ( (dom.pos - pos).len() < radius ) ndoms++;
    }

  foreach( hit, hits )
    {
      if ( (hit.pos - pos).len() < radius ) nhits++;
    }
  
  if ( ndoms == 0 ) return 0;
  return double(nhits)/ndoms;
}

/*! structure to compute the parameters describing the relation between 
    a track (shower) and PMT */

struct Cherenkov_pars
{
  double t;              ///< time of arrival
  double d_track;        ///< distance the muon travels along the track
  double d_photon;       ///< distance the photon must travel
  double d_closest;      ///< distance of closest approach between track and pmt
  double cos_theta;      ///< theta = angle between PMT and track direction.
  double cos_phi;        ///< phi = pi(0) means pmt looks at(away) track in perpendicular plain
  double cos_angle;      ///< angle between photon and pmt direction
  Vec k;
  double l;

  string __str__()
  {
    return form( "%4.3f  %4.3f  %4.3f  %4.3f  %4.3f  %4.3f",
		 t, d_track, d_closest, cos_theta, cos_phi, cos_angle );
  }

  Cherenkov_pars( const Trk& trk , const Hit& hit ) : Cherenkov_pars( trk, hit.pos, hit.dir ) {}

  Cherenkov_pars( const Trk& trk , const Vec& p, const Vec& w )
  {
    init_track( trk, p, w );
  }


  void init_track(  const Trk& trk , const Vec& p, const Vec& w )
  {
    const Vec v = p - trk.pos;
    l = v.dot(trk.dir);
    k = v - trk.dir * l;
    
    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 );
    t         = 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 );
    cos_theta = w.dot( trk.dir );

    if ( cos_theta > 1 ) cos_theta = 1;
    if ( cos_theta <-1 ) cos_theta =-1;
    
    const double len_w_plane = sqrt( 1 - cos_theta * cos_theta );
    
    if ( d_closest * len_w_plane < 1e-8 ) cos_phi = 1;
    else                                  cos_phi   = w.dot( k ) / ( d_closest * len_w_plane ) ;
    
    if ( cos_phi < -1 ) cos_phi = -1;
    if ( cos_phi > 1  ) cos_phi = 1;
  }
  

  ClassDefNV( Cherenkov_pars ,1 );
};

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