#!/usr/bin/env python


""" 
    This example shows how to compute residuals of hits with resprect
    To monte carlo tracks. It also illustrates the use of C++ code
    fom python and the Timer utility.

    The hits that have passed though jtriggerefficieny have a time
    offset. This is because the trigger simulation must first put the
    event in a timeslice somewhere. We have to deal with a few different
    quantities.
    hits[].t is the time in nanoseconds of the hit with respect to
             the beginning of the timeslice. (between 0 and 1e8)
   
    evt.t    is a ROOT::TTimeStamp representing the start of the 
             timeslice in which the event occured. It is present in
             both data and MC. The exact time, Te, of a hit is therefore
             evt.t + hit.t
             
    evt.mc_t Is the MonteCarlo time. It corresponds to the Te a hit
             with t=0 will have after passing through triggerefficiency.
   
    Note all this is not applicable to the mc_hits which are not 
    effected at all by the trigger simulation.

"""



import aa, ROOT
from ROOT import EventFile, Det, cxx, TH1D, Timer, time_at
ROOT.gStyle.SetOptStat(0)


infile = aa.aadir+"/data/mc5.1.numuCC.ps10.aa.root"

f = EventFile( infile )


detfile = aa.aadir+'/data/KM3NeT_-00000001_20171212.detx'

det = Det( detfile )
 


# If we do the histogram filling in c++, it will be a lot faster.

use_cxx = True

if use_cxx :

    cxx("""
    void fill_residuals( vector<Hit>& hits, 
                         Trk& nu,
                         TH1D& hist, double offset = 0 )
    {
       for(auto& hit : hits )
       {
          double residual = hit.t - time_at( nu, hit.pos ) -offset ; // time_at from Trk_util.hh
          hist.Fill(residual);
       }
    }
    """)

    from ROOT import fill_residuals

else : # but it works exactly the same in python.

    def fill_residuals( hits, nu, hist, offset =0 ) :

        for hit in hits : 
            residual = hit.t - time_at( nu, hit.pos) - offset 
            hist.Fill( residual )



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


for evt in f :
    
    print(evt)    

    with Timer ("event loop") :

        # set the hit positions using the detector

        det.apply(evt)


        # -- we need two corrections 1) the offset of the hits
        #    and 2) the spatial offset of the monte-carlo tracks with 
        #    respect to the detector coordinates.
        
        offset = evt.mc_t - evt.t.AsDouble()*1e9

        for t in evt.mc_trks :
            t.pos += f.header.coord_origin();


        # -- ready to copute some residuals.

        nu = evt.neutrino()

        with Timer("fill mc_hits") :
            fill_residuals( evt.mc_hits, nu , hres1 )

        with Timer("fill hits") :
            fill_residuals( evt.hits   , nu , hres2, offset )
    

Timer.report()


for h in [hres1, hres2] : h.Scale( 1.0 / h.Integral() )
hres1.SetLineColor(2)
hres2.SetLineColor(4)
hres1.Draw()
hres2.Draw("same")
ROOT.gPad.Update()

aa.i() # interactive: start a prompt