// Root macro 
// run this with: root $AADIR/aa.C residuals.C
//
// see residuals.py for same example in python.


#include "../util/foreach.hh"


void fill_residuals( vector<Hit>& hits, 
		     Trk& nu,
		     TH1D& hist, double offset = 0 )
{
  for (int i=0; i<hits.size(); i++)
    {
      Hit& hit = hits[i];
      double d = (hit.pos - nu.pos).len();
      double r = hit.t - ( nu.t + d / v_light ) - offset;
      hist.Fill(r);
    }
}



TH1D hres1("hres1","hres1;residual(ns)",125,-50,200);
TH1D hres2("hres2","hres2",125,-50,200);


void residuals()
{
  EventFile f( "/afs/in2p3.fr/throng/km3net/dat/benchmark/JTE/v7/genhen_v7r6_nueCC_1.km3_v5r3.jte_v7.root" );
  
  string detfile = "/afs/in2p3.fr/home/throng/km3net/detectors/"+f.header.get_line("detector");
  Det det(detfile);
 
  // foreach does not work in CINT (will be better with cling/c++11!). 
  
  for (EventFile::iterator ievt=f.begin(); ievt!=f.end(); ++ievt)
    {
      Evt& evt = *ievt;

      det.apply(evt);
      
      if (f.index > 300) break;
      if (evt.mc_hits.size() < 1000 ) continue;

      double offset = evt.mc_t - evt.t.AsDouble()*1e9;

      // shift tracks to match detector
      
      for(int i = 0; i < evt.mc_trks.size() ;++i)
	{
	  evt.mc_trks[i].pos += f.header.coord_origin();
	}
					  
      Trk& nu = evt.mc_trks[0];

      cout << evt.hits.size() << endl;
      
      fill_residuals( evt.mc_hits, nu, hres1 );
      fill_residuals( evt.hits,    nu, hres2 , offset);
      
    }

  hres1.Scale( 1.0 / hres1.Integral() );
  hres2.Scale( 1.0 / hres2.Integral() );

  hres1.SetLineColor(2);
  hres2.SetLineColor(4);
  hres1.Draw();
  hres2.Draw("same");
  
}