caf_analysis.py

Go to the documentation of this file.

print "\nrun  : --- CAF analysis"
########################### Option parser
from optparse import OptionParser
parser = OptionParser()
parser.add_option("-n", "--number_of_events", help="only run over a certain number of events", action="store",      type="int",         dest="n_events",        default=False)
parser.add_option("-i", "--input_file",       help="input file names",                         action="store",      type="string",      dest="input_files",     default=False)
parser.add_option("-o", "--output_file",      help="output file name",     action="store", type="string", dest="output_file", default="inspect_dataset.root")
parser.add_option("-s", "--skip_spills",      help="don't compute PoT",    action="store_true",           dest="skip_spills", default=False)
parser.add_option("-v", "--verbose",          help="turn on verbose mode", action="store_true",           dest="verbose",     default=False)
(options, args) = parser.parse_args()
print "run  : --- Options"
print "run  :     number of events: ",options.n_events
print "run  :     input files:      ",options.input_files
print "run  :     output file name: ",options.output_file
print "run  :     skip spills:      ",options.skip_spills
print "run  :     verbose mode:     ",options.verbose
assert options.input_files, "No input files specified, provide with -i"
options.mc_mode = False
if ("genie" in options.input_files) or ("cry" in options.input_files): options.mc_mode = True
print "run  :     Monte Carlo mode:  ",options.mc_mode
################################## Setup ROOT
print "run  :     Importing ROOT"
import ROOT
ROOT.gROOT.SetBatch(1)
################################## Imports
import caf_tools as tools
from collections import OrderedDict
import time
################################## Book histograms
print "run  : --- Booking histograms"
figures = {}
#### Reco plots
figures["reco_neutrino_energy"]      = ROOT.TH1F("reco_neutrino_energy",      "reco_neutrino_energy",      100,   0.,   5.)
figures["reco_muon_length"]          = ROOT.TH1F("reco_muon_length",          "reco_muon_length",          100,   0.,   1800)
figures["reco_hadronic_energy"]      = ROOT.TH1F("reco_hadronic_energy",      "reco_hadronic_energy",      100,   0.,   2.5)
figures["reco_fit_hadronic_energy"]  = ROOT.TH1F("reco_fit_hadronic_energy",  "reco_fit_hadronic_energy",  100,   0.,   2.5)
#### True plots
figures["true_neutrino_energy"] = ROOT.TH1F("true_neutrino_energy", "true_neutrino_energy", 100,   0.,   5.)
figures["true_muon_energy"]     = ROOT.TH1F("true_muon_energy",     "true_muon_energy",     100,   0.,   5.)
figures["true_hadronic_energy"] = ROOT.TH1F("true_hadronic_energy", "true_hadronic_energy", 100,   0.,   3.)
figures["true_x"]               = ROOT.TH1F("true_x",               "true_x",               100,   0.,   1.5)
figures["true_y"]               = ROOT.TH1F("true_y",               "true_y",               100,   0.,   1.)
figures["true_w2"]              = ROOT.TH1F("true_w2",              "true_w2",              100,   0.,   4.)
figures["true_q2"]              = ROOT.TH1F("true_q2",              "true_q2",              100,   0.,   2.)
#### Multi-dimensional plots
keys = figures.keys()
# dependent_variables = ["true_neutrino_energy"]
# for d_v in dependent_variables:
#     for f in keys:
#         if f == d_v: continue
#         x_bins = tools.GetBins(figures[d_v].GetXaxis())
#         y_bins = tools.GetBins(figures[f].GetXaxis())
#         figures["%s_vs_%s"%(d_v,f)] = ROOT.TH2F("%s_vs_%s"%(d_v,f), "%s_vs_%s"%(d_v,f), 
#                                                 len(x_bins)-1, x_bins[0], x_bins[-1], 
#                                                 len(y_bins)-1, y_bins[0], y_bins[-1])
#### Store plot keys
figure_keys = figures.keys()
#### Sub-samples
if options.mc_mode:
    int_types = [
                 "QE","Res","DIS", "Coh","MEC",
                 ]
    keys = figures.keys()
    for f in keys:
        for t in int_types:
            name = "%s-%s"%(f,t)
            figures[name] = figures[f].Clone(name)
################################## Initialise a counter dictionary
counters = OrderedDict()
counters["spills"]                      = 0
counters["events"]                      = 0
counters["pass spill"]                  = 0
counters["pass preselection"]           = 0
counters["pass quality"]                = 0
counters["pass energy"]                 = 0
counters["pass containment"]            = 0
counters["pass loose cosmic rejection"] = 0
counters["pass cosmic rejection"]       = 0
counters["pass remid"]                  = 0
counters["pass numu cuts"]              = 0
counters["pass true cuts"]              = 0
figures["counters"] = ROOT.TH1F("meta-counters",  "counters",  len(counters.keys()),  0.5,  len(counters.keys())+0.5)
#### Meta histograms
figures["TotalPOT"]  = ROOT.TH1F("meta-TotalPOT", "TotalPOT",              1,   0.,   1)
################################## NOvA setup
print "run  : --- Loading CAF libraries"
ROOT.gSystem.Load("libCintex.so")
ROOT.Cintex.Enable()
ROOT.gSystem.Load("$SRT_PUBLIC_CONTEXT/lib/Linux2.6-GCC-maxopt/libStandardRecord_dict.so")
################################## load input file
print "run  : --- Input file"
if options.input_files:
    if "*" in options.input_files:
        import glob
        print "run  :     globing input files"
        files = glob.glob(options.input_files)
    else:
        files = options.input_files.split(",")
        print "run  :     using command line list of %i input files"%(len(files))
    files.sort()
    assert len(files),"run  :     no input files found"
    
    r_chain_name = "recTree"
    r_chain = ROOT.TChain(r_chain_name)

    s_chain_name = "spillTree"
    s_chain = ROOT.TChain(s_chain_name)
    
    for this_file in files:
        print "run  :         - %s"%this_file
        r_chain.Add(this_file)
        s_chain.Add(this_file)
        # assert r_chain.GetEntries(), "run  :     No entries in input file: %s, an error has probably occured"%this_file
else:
    assert False,"run  :     no input file given"
print "run  :     All files open"
nevents   = r_chain.GetEntries()
s_nevents = s_chain.GetEntries()
counters["spills"] = s_nevents
print "run  :     Dataset: %i file(s)"%len(files)
print "run  :         rec tree has {:,} events".format(nevents)
print "run  :       spill tree has {:,} events".format(s_nevents)
# assert nevents,"run  : No entries in any input file(s), an error has probably occurred"
for f in files:
    print "run  :         - %s"%f
if options.n_events: 
    nevents = options.n_events
    print "run  :     running on {:,} event(s)".format(nevents)
################################## Optimisations
active_branches = [
                   "hdr.run",
                   "hdr.subrun",
                   "hdr.evt",
                   "hdr.det",
                   "sel.remid.len",
                   "sel.remid.bestidx",
                   "energy.numusimp.trkccE",
                   "energy.numusimp.hadcalE",
                   "energy.numusimp.hadtrkE",
                   "energy.mutrkE.E",
                   "mc.nu",
                   "mc.nnu",
                   "mc.nu.iscc",
                   "mc.nu.pdg",
                   "mc.nu.pdgorig",
                   "mc.nu.E",
                   "mc.nu.*",
]
for a in tools.fa_nu_mu_vars: active_branches.append(a)
r_chain.SetBranchStatus("*",0)
for b in active_branches: r_chain.SetBranchStatus(b,1)
################################## Run analysis
print "run  : --- Beginning Analysis"
print "run  :     Beginning event loop"
start_time = time.time()
if options.verbose:     print_every_n_events = 1
elif nevents < 1000000: print_every_n_events = 1000
else:                   print_every_n_events = 10000
for event in xrange(nevents):
    ################################## Begin loop
    counters["events"] += 1
    ################################## Get entries
    if r_chain.GetEntry(event) <=0: print "run  :    couldn't get rec entry: ",event
    if event%print_every_n_events == 0:
        dt = time.time() - start_time
        print 'run  : --- Event %10s (%4.1f%%), %6.1f s'%\
         ("{:,}".format(event),( (float(event)/float(nevents))*100.), dt)
        if options.verbose:
            print "run  :     Run: %i, sub-run: %i, event number: %i, detector: %i"%\
                (r_chain.rec.hdr.run,r_chain.rec.hdr.subrun,r_chain.rec.hdr.evt, r_chain.rec.hdr.det)
            if (options.mc_mode): tools.dumpTruth(r_chain)
    ################################## Nu mu analysis cuts
    # These are the standard nu mu analysis cuts used in first analysis
    pass_numu_cuts, counters = tools.numuFACuts(r_chain, counters=counters)
    if (not pass_numu_cuts): continue
    
    # if we get here then we've passed all the nu mu cuts
    counters["pass numu cuts"] += 1

    # end of nu mu cuts
    ################################## Fill reconstruction values
    weight = 1.
    values = {}
    values["reco_neutrino_energy"]     = r_chain.rec.energy.numusimp.trkccE
    values["reco_muon_length"]         = r_chain.rec.sel.remid.len
    values["reco_hadronic_energy"]     = r_chain.rec.energy.numusimp.hadcalE+r_chain.rec.energy.numusimp.hadtrkE
    values["reco_fit_hadronic_energy"] = r_chain.rec.energy.numusimp.trkccE-r_chain.rec.energy.mutrkE[r_chain.rec.sel.remid.bestidx].E
    ################################## Truth studies
    # if the '-m' option is set then we assume this is MC
    # if so we run some simple truth analysis
    if options.mc_mode:
        true_event_types = tools.DetermineTruthEventType(r_chain, verbose=options.verbose)
        # if there is more than one truth neutrino this event is too complicated
        if (len(r_chain.rec.mc.nu) == 1):
            true_neutrino = r_chain.rec.mc.nu[0]
            true_muon     = None

            for true_particle in true_neutrino.prim:
                # select the most energetic muon
                if true_particle.pdg in [13,-13]:
                    if (not true_muon) or (true_particle.p.E() > true_muon.p.E()): true_muon = true_particle

            # if we didn't find a true muon then this event is no good
            if true_muon:
                # if we get here then we've passed the true cuts
                counters["pass true cuts"] += 1

                ################################## Fill truth values
                values["true_neutrino_energy"] = true_neutrino.E
                values["true_muon_energy"]     = true_muon.p.E()
                values["true_hadronic_energy"] = true_neutrino.E - true_muon.p.E()
                values["true_x"]               = true_neutrino.x
                values["true_y"]               = true_neutrino.y
                values["true_w2"]              = true_neutrino.W2
                values["true_q2"]              = true_neutrino.q2
    else:
        true_event_types = [""]
    ################################## Fill histograms
    for figure in figure_keys:
        if "vs" in figure:
            x_variable = figure.split("_vs")[0]
            y_variable = figure.split("vs_")[1]
            if (x_variable not in values.keys()) : continue
            if (y_variable not in values.keys()) : continue
            for event_type in true_event_types:
                figures["%s%s"%(figure,event_type)].Fill(values[x_variable],values[y_variable],weight)
        else:
            if (figure not in values.keys()) : continue
            for event_type in true_event_types:
                figures["%s%s"%(figure,event_type)].Fill(values[figure],weight)
################################## Event loop finished
print "run  : --- Event loop finished after {:,} events".format(nevents)
dt = time.time() - start_time
print "run  :     time: %.1f s, %.1f Hz"%(dt,nevents/dt)
################################## Bad timing list
bad_runs_file  = "/nova/app/users/tamsett/art/python/bad_timing_runs.txt"
bad_runs_list  = set([int(r.strip()) for r in open(bad_runs_file,"r").readlines()])
print "run  : found %i runs with bad timing"%len(bad_runs_list)
################################## Spill tree
print "run  : --- Looping over spill tree"
total_pot = 0.
start_time = time.time()
if options.verbose:       print_every_n_events = 1
elif s_nevents < 1000000: print_every_n_events = 1000
else:                     print_every_n_events = 10000
for event in xrange(s_nevents):
    if options.skip_spills: break
    if s_chain.GetEntry(event) <=0: print "run  :    couldn't get spill entry: ",event        
    if event%print_every_n_events == 0:
        dt = time.time() - start_time
        print 'run  : --- Event %10s (%4.1f%%), %6.1f s'%\
         ("{:,}".format(event),( (float(event)/float(s_nevents))*100.), dt)
        if options.verbose: print "run  :     Run: %i, sub-run: %i, event number: %i"%(s_chain.spill.run,s_chain.spill.subrun,s_chain.spill.evt)
    trigger_int = int(repr(s_chain.spill.trigger)[3:-1])
    if (trigger_int != 0):
        print "run  : Non-beam trigger found: %s"%repr(s_chain.spill.trigger)
        print "run  :     Run: %i, sub-run: %i, event number: %i"%(s_chain.spill.run,s_chain.spill.subrun,s_chain.spill.evt)
    ################################## Analysis cuts
    # the full analysis cuts should be applied here
    # timing dq
    timing_dq = (s_chain.spill.run not in bad_runs_list)
    if not (timing_dq): continue
    
    spill = tools.standardSpillCut(s_chain.spill)
    if not (spill): continue
    
    total_pot += s_chain.spill.spillpot     
figures["TotalPOT"].SetBinContent(1,total_pot)
print "run  : --- Event loop finished after {:,} events".format(nevents)
print "run  :     total PoT: ",total_pot
dt = time.time() - start_time
print "run  :     time: %.1f s, %.1f Hz"%(dt,nevents/dt)
################################## print counters
print "run  : --- counters"
for i,counter in enumerate(counters.keys()):
    figures["counters"].SetBinContent(i+1,counters[counter])
    figures["counters"].GetXaxis().SetBinLabel(i+1,counter)
    print "run  :     %-25s: %s"%(counter,"{:,}".format(counters[counter]))
################################## Save histograms
print "run  :  Writing to: %s"%(options.output_file)
saved_root_file = ROOT.TFile(options.output_file, "RECREATE" )
saved_root_file.cd()
for i,f in enumerate(figures.keys()):
    if (("true" in f) and ("reco" in f)):
        prefix = "mixed-"
    elif ("true" in f):
        prefix = "true-"
    elif ("reco" in f):
        prefix = "reco-"
    else:
        prefix = ""
    figures[f].SetName(prefix+figures[f].GetName())
    figures[f].Write()
print "run  :  Wrote %i histograms"%(len(figures.keys()))
saved_root_file.Close()
################################## Done
print "run  : --- CAF Analysis finished\n"