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

# This script shows how to use the event weights to compute event rates
# for a given flux, and the neutrino effective area.
# 
# This example runs on a jgandalf output file, which already contains the
# MC information in aanet format. So there is no need to convert it, and
# we can immediately use the aanet trie


import aa
from ROOT import *
gStyle.SetOptStat(0)
from math import *

# Open the event files and get the total number of generated events
# from the header.

infile = "../data/mc5.1.numuCC.nohits.aa.root"

f = EventFile(infile)

print f.header
ngen = float ( f.header.get_field("genvol","numberOfEvents") ) 


# Make a ROOT Canvas
c = TCanvas("c","c",300,900)
c.Divide(1,3)


# Define the histograms, and their axis labels..

hevts  = TH1D("hevts","events for flux of 1e-4 (E/GeV)^{-2} GeV^{-1} m^{-2} s^{-1};log(E_{#nu});events yr^{-1}", 24, 1,7 )
haeff  = TH1D("haeff","effective area;log(E_{#nu});m^2", 24, 1,7 )
hcheck = TH1D("hcheck","event rate computed with Aeff;log(E_{#nu});events yr^{-1}", 24, 1,7 )


#-----------------------------------------------------------------------------------------
# plot the number of events/yr for a flux of 1e-8 GeV-1 cm-2 s-1 into the hevts histogram
#-----------------------------------------------------------------------------------------

c.cd(1)
# note that you have access to E, simply because root has opened the file.
E.Draw("log10(mc_trks[0].E)>>+hevts","w[1]*1e-4*mc_trks[0].E**-2*"+str(1.0/ngen) )

#-----------------------------------------------------------------------------------------
# The effective is related to the response to a flat flux (i.e. 1/E as we have log bins)
# See http://cherenkov.nl:8888/notebooks/aanet_tutorial/05%20-%20Effective%20Area.ipynb
# for details.
#-----------------------------------------------------------------------------------------

c.cd(2)
E.Draw("log10(mc_trks[0].E)>>+haeff","w[1]/mc_trks[0].E","goff")
haeff.Scale ( 1 / ( ngen * log(10) * haeff.GetBinWidth(1) * 3600.0 * 24 * 365 * 4*pi ))
gPad.SetLogy()
haeff.Draw()

#----------------------------------------------------------------------------------------
# If you work out the flux in each (logarithmic bin) and multiply with the effective
# Area, you should get the event rate again...
#----------------------------------------------------------------------------------------

c.cd(3)
for b in range (1, hcheck.GetNbinsX()+1) :

    E    = 10**hcheck.GetBinCenter(b)
    flux = 1e-4 * E **-2
    Aeff = haeff.GetBinContent(b)

    # how many neutrinos in this bin in a year?
    Nnus = flux * ( 10**hcheck.GetXaxis().GetBinUpEdge(b) - 10**hcheck.GetXaxis().GetBinLowEdge(b) )
    Nnus *= 4*pi * 3600 * 24 * 365
    hcheck.SetBinContent(b, Nnus * Aeff )

hcheck.Draw()
print hevts.Integral(), hcheck.Integral()