#!/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")