// #include "Hit.hh"
#include "util.hh" // recutil
#include "constants.hh"

#include "../definitions/reconstruction.hh"


struct cascade_summary
{
  // aashower
  int nhits_aashower;
  int ntot_aashower;
  double goldP_aashower;

  int nhits_aashower_early_10;
  int ntot_aashower_early_10;

  int nhits_aashower_early_20;
  int ntot_aashower_early_20;

  int nhits_aashower_coin_early_10;
  int ntot_aashower_coin_early_10;

  int nhits_aashower_coin_early_20;
  int ntot_aashower_coin_early_20;

  int nhitscrkv_aashower;
  int ntotcrkv_aashower;

  // qstrategy
  int nhits_qstrat;
  int ntot_qstrat;
  double goldP_qstrat;

  vector<double> inertia_tensor_eigenvalues;
  double inertia_ratio;

  // reco
  double angle_aa_jpp; // degrees
  double d_aa_jpp; // meters

};

const double crkv_time_hypothesis = 15;

// aashowerfit hypotheses
const double res_min_aashower = -100; // [ns]
const double res_max_aashower =  900; // [ns]
const double d_max_aashower   =  800; // [m]

// QStrategy hypotheses
const double res_min_qstrat  = -12;
const double res_max_qstrat  =  25;
const double d_max_qstrat    = 130;


//
// Gold parameter: https://drive.google.com/file/d/0BxzGGqhzVJmZYlk1clJVQldYbWs/view?resourcekey=0-iBI2eci7kQa9LB_-Bj0d6Q
// formula 10
//
double calculate_goldP( vector<Hit> hits, Trk trk )
{
  double result = 0;

  if (hits.size() == 0) return -1;

  for ( auto & h : hits )
  { 
    double d = (h.pos - trk.pos).len();
    double res = h.t - trk.t - d/v_light;
    // TODO add time resolution dom?
    double tsigma_dom = 1;
    result+= exp( - ( pow(res,2) / pow(tsigma_dom*1.5,2) ) ) / tsigma_dom;
  }

  return result/hits.size();
}

bool fill_qstrat( cascade_summary& s, vector<Hit> hits, Trk trk )
{

  TMatrixD tensor(3,3);

  for ( auto & h : hits )
  {
    s.nhits_qstrat+=1;
    s.ntot_qstrat+=h.tot;

    Vec d = (h.pos - trk.pos);

    double D = d.len();
    vector<double> r_k = {d.x, d.y, d.z};

    for (int i = 0; i < 3; i++)
    {
      for (int j = 0; j < 3; j++)
      {
        if ( i == j ) TMatrixDRow(tensor,i)(j) += pow(D,2) - r_k[i]*r_k[j];
        else TMatrixDRow(tensor,i)(j) += - r_k[i]*r_k[j];
      }
    }
  }

  TVectorD eigenval;
  TMatrixD vectors = tensor.EigenVectors(eigenval);

  double min_value = eigenval[0];
  double sum = 0;

  for (int i = 0; i < 3; i++)
  {
    s.inertia_tensor_eigenvalues.push_back( eigenval[i] );
    sum+=eigenval[i];
    if (eigenval[i] < min_value) min_value = eigenval[i];
  }

  s.inertia_ratio = min_value/sum;

  return true;

}

bool fill_aashower( cascade_summary& s, vector<Hit> aashower_hits, Trk trk )
{

  for ( auto & h : aashower_hits )
  {
    s.nhits_aashower+=1;
    s.ntot_aashower+=h.tot;
    if ( is_hit_from_casc(trk,h, 0) )
    {
      s.nhitscrkv_aashower+=1;
      s.ntotcrkv_aashower +=h.tot;
    }

    double d = (h.pos - trk.pos).len();
    double res = h.t - trk.t - d/v_light;

    if ( res < -10 )
    {
      s.nhits_aashower_early_10+=1;
      s.ntot_aashower_early_10+=h.tot;
    }
    if ( res < -20 )
    {
      s.nhits_aashower_early_20+=1;
      s.ntot_aashower_early_20+=h.tot;
    }

  }

  return true;

}

bool fill_coin( cascade_summary& s, vector<Hit> hits, Trk trk )
{

  for ( auto & h : hits )
  {

    double d = (h.pos - trk.pos).len();
    double res = h.t - trk.t - d/v_light;

    if ( res < -10 )
    {
      s.nhits_aashower_coin_early_10+=1;
      s.ntot_aashower_coin_early_10+=h.tot;
    }
    if ( res < -20 )
    {
      s.nhits_aashower_coin_early_20+=1;
      s.ntot_aashower_coin_early_20+=h.tot;
    }

  }

  return true;

}

inline bool summarize_cascade( Evt& evt,
                               cascade_summary& s ) 
{
  s.nhits_aashower =  0; s.ntot_aashower = 0; s.nhitscrkv_aashower = 0; s.ntotcrkv_aashower = 0; s.nhits_qstrat = 0; s.ntot_qstrat = 0;s.goldP_aashower = 0; s.goldP_qstrat = 0;
  s.nhits_aashower_early_10 =  0;s.nhits_aashower_early_20 =  0;s.ntot_aashower_early_10 =  0;s.ntot_aashower_early_20 =  0;
  s.nhits_aashower_coin_early_10 = 0;s.nhits_aashower_coin_early_20 = 0;s.ntot_aashower_coin_early_10 = 0;s.ntot_aashower_coin_early_20 = 0;
  s.inertia_tensor_eigenvalues.clear();
  s.inertia_ratio = 2;
  s.d_aa_jpp = -1; s.angle_aa_jpp = -1;

  Trk* best = best_reco_track( evt, AANET_RECONSTRUCTION_TYPE );
  if (!best) 
  {
    cout << "ERROR: event has no aashowerfit track";
    return true;
  }

  Trk& rectrk = *best;

  vector<Hit> aashower_hits = filter_hits( evt.hits, rectrk, res_min_aashower, res_max_aashower, d_max_aashower);
  vector<Hit> coin_hits     = get_coincidences( evt.hits, 2, 20 );
  vector<Hit> qstrat_hits   = filter_hits( coin_hits, rectrk, res_min_qstrat, res_max_qstrat, d_max_qstrat );  

  if ( qstrat_hits.empty() )
  {
    s.inertia_ratio = 3;
    return true;
  }

  fill_aashower( s, aashower_hits, rectrk );
  fill_qstrat  ( s, qstrat_hits, rectrk );
  fill_coin    ( s, coin_hits, rectrk );

  s.goldP_aashower = calculate_goldP(aashower_hits, rectrk); 
  s.goldP_qstrat = calculate_goldP(qstrat_hits, rectrk); 


  // reco extra options
  Trk* best_jpp = best_reco_track( evt, JPP_RECONSTRUCTION_TYPE );
  if (!best_jpp) 
  {
    cout << "ERROR: event has no Jpp track";
  }
  else
  {
    Trk& rectrk_jpp = *best_jpp;

    s.d_aa_jpp = (rectrk_jpp.pos - rectrk.pos).len();
    s.angle_aa_jpp = angle_between( rectrk_jpp.dir, rectrk.dir )*TMath::RadToDeg();

  }

  return true;
}