#ifndef AHRECUTILINCD
#define AHRECUTILINCD

#include "Hit.hh"
#include "Evt.hh"
#include "Det.hh"
#include "Spline.hh"
#include "Mat.hh"
#include "EventFile.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 );
};



constexpr int BADRATE =15; //! bit-number for use in pmt.status word


/*! set bit 15 of the pmt status flags if the count-rate is out of bounds */

inline int flag_bad_rates( Det& det, double min_rate = 2000, double max_rate = 1.95e6)
{
    const unsigned int mask = 1 << BADRATE;
    int n = 0;
    foreach_map( pmt_id, pmt, det.pmts ) 
    {
        bool ok = pmt->rate > min_rate && pmt->rate < max_rate;
        pmt->status &= ~mask;
        if (!ok) {
          pmt->status |= mask;
          //print("status is now", pmt->status );
          n++;
        }
    }
  return n;   
}

/*! Only keep the hits who's pmt status == 0 */

inline void filter_hits_by_status( Det& det, vector<Hit>& hits )
{
    remove_element_if( hits, [&](Hit& h){ return (det.doms[h.dom_id].pmts[h.channel_id].status != 0); } );
}

// #-----------------------------------------------
// # summary slice reading 
// #-----------------------------------------------


/*! Helper for reading timeslice info into the Pmt::rate fields of a Det */

#include "km3net-dataformat/definitions/root.hh" // from km3net-dataformat/tools

struct RateReader {

  TTree* sumtree;
  KM3NETDAQ::JDAQSummaryslice* sumslice;
  EventFile* eventfile;

  RateReader( EventFile& f ) : eventfile(&f)
  {
    prepare_summary_frame( eventfile->rootfile() );
    sumtree = (TTree*) eventfile->rootfile()->Get(TTREE_ONLINE_SUMMARYSLICE); // aka "KM3NET_SUMMARYSLICE"

    if (!sumtree) fatal("failed to get summary slice tree from file.")
    sumslice = new KM3NETDAQ::JDAQSummaryslice();

    sumtree->SetBranchAddress( TBRANCH_ONLINE_SUMMARYSLICE, &sumslice );
    int n = sumtree->BuildIndex("frame_index"); // this works, eventhough frame_index is procected; thank god.  

    if ( n < 1 ) fatal("building index for Summary Slice tree failed ", n );

    // this loads the full tree in memory -- hopefully for fast random access
    // it works -- loading time is from 0.1 s to 0.002 sec. 
    // But not sure it fully fits in memory for data-file.

    sumtree->LoadBaskets(); // default = 2 GB maximum
  }

  void operator() ( Det& det ) {
    int n = sumtree->GetEntryWithIndex( eventfile->evt.frame_index );

    if ( n < 1 ) fatal("failed to read entry with index ", eventfile->evt.frame_index, n );

    read ( det, *sumslice );
    // print("got frame for index", eventfile->evt.frame_index );
  }
};

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