#!/usr/local/python/python-2.7/bin/python -i


import aa
from ROOT import *
from math import log10


f = EventFile("/afs/in2p3.fr/home/h/heijboer/vici/gSeaGen/test_outputs/dus.1.aa.root")

h_nu_energy = TH1D("h_nu_energy","log of energy of neutrinos, weighted weith E^{-2} spectrum",
                   28,1,8 )

h_nmuons = TH1D("h_nmuons","number of muons (unweighted)", 10,0,10)

h_yy = TH2D("h_yy","computed vs generated y variable",
            30,0,1, 30,0,1 )


# print all the header information
print f.header

# get the number of generated neutrino interactions from the header
ngen = f.header.ngen()


for evt in f :

 
    nu = evt.primary_neutrino()


    flux = 1e-4 * nu.E**-2 # E-2 flux with units : GeV-1 cm-2, sr-1 s-1
    
    h_nu_energy.Fill( log10 ( nu.E ) ,
                      evt.w[1] / ngen * flux )  # weight the event 
    

    # get all muons

    all_muons = [ trk for trk in evt.mc_trks if abs( trk.type ) == 13 ]

    h_nmuons.Fill( len( all_muons ) )
        

    # compute ('bjorken') y variable

    y = 1 - evt.leading_lepton().E / evt.primary_neutrino().E


    # compare with the one from the w2list

    h_yy.Fill( y, evt.w2list[8] )



c = TCanvas("c","c",900,900)
c.Divide(2,2)

c.cd(1);
h_nu_energy.Draw()

c.cd(2);
h_nmuons.Draw()

c.cd(3);
h_yy.Draw("colz")