1 print "\nrun : --- CAF analysis" 3 from optparse
import OptionParser
4 parser = OptionParser()
5 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)
6 parser.add_option(
"-i",
"--input_file", help=
"input file names", action=
"store", type=
"string", dest=
"input_files", default=
False)
7 parser.add_option(
"-o",
"--output_file", help=
"output file name", action=
"store", type=
"string", dest=
"output_file", default=
"inspect_dataset.root")
8 parser.add_option(
"-s",
"--skip_spills", help=
"don't compute PoT", action=
"store_true", dest=
"skip_spills", default=
False)
9 parser.add_option(
"-v",
"--verbose", help=
"turn on verbose mode", action=
"store_true", dest=
"verbose", default=
False)
10 (options, args) = parser.parse_args()
11 print "run : --- Options" 12 print "run : number of events: ",options.n_events
13 print "run : input files: ",options.input_files
14 print "run : output file name: ",options.output_file
15 print "run : skip spills: ",options.skip_spills
16 print "run : verbose mode: ",options.verbose
17 assert options.input_files,
"No input files specified, provide with -i" 18 options.mc_mode =
False 19 if (
"genie" in options.input_files)
or (
"cry" in options.input_files): options.mc_mode =
True 20 print "run : Monte Carlo mode: ",options.mc_mode
22 print "run : Importing ROOT" 24 ROOT.gROOT.SetBatch(1)
26 import caf_tools
as tools
27 from collections
import OrderedDict
30 print "run : --- Booking histograms" 33 figures[
"reco_neutrino_energy"] = ROOT.TH1F(
"reco_neutrino_energy",
"reco_neutrino_energy", 100, 0., 5.)
34 figures[
"reco_muon_length"] = ROOT.TH1F(
"reco_muon_length",
"reco_muon_length", 100, 0., 1800)
35 figures[
"reco_hadronic_energy"] = ROOT.TH1F(
"reco_hadronic_energy",
"reco_hadronic_energy", 100, 0., 2.5)
36 figures[
"reco_fit_hadronic_energy"] = ROOT.TH1F(
"reco_fit_hadronic_energy",
"reco_fit_hadronic_energy", 100, 0., 2.5)
38 figures[
"true_neutrino_energy"] = ROOT.TH1F(
"true_neutrino_energy",
"true_neutrino_energy", 100, 0., 5.)
39 figures[
"true_muon_energy"] = ROOT.TH1F(
"true_muon_energy",
"true_muon_energy", 100, 0., 5.)
40 figures[
"true_hadronic_energy"] = ROOT.TH1F(
"true_hadronic_energy",
"true_hadronic_energy", 100, 0., 3.)
41 figures[
"true_x"] = ROOT.TH1F(
"true_x",
"true_x", 100, 0., 1.5)
42 figures[
"true_y"] = ROOT.TH1F(
"true_y",
"true_y", 100, 0., 1.)
43 figures[
"true_w2"] = ROOT.TH1F(
"true_w2",
"true_w2", 100, 0., 4.)
44 figures[
"true_q2"] = ROOT.TH1F(
"true_q2",
"true_q2", 100, 0., 2.)
57 figure_keys = figures.keys()
61 "QE",
"Res",
"DIS",
"Coh",
"MEC",
67 figures[name] = figures[f].
Clone(name)
69 counters = OrderedDict()
70 counters[
"spills"] = 0
71 counters[
"events"] = 0
72 counters[
"pass spill"] = 0
73 counters[
"pass preselection"] = 0
74 counters[
"pass quality"] = 0
75 counters[
"pass energy"] = 0
76 counters[
"pass containment"] = 0
77 counters[
"pass loose cosmic rejection"] = 0
78 counters[
"pass cosmic rejection"] = 0
79 counters[
"pass remid"] = 0
80 counters[
"pass numu cuts"] = 0
81 counters[
"pass true cuts"] = 0
82 figures[
"counters"] = ROOT.TH1F(
"meta-counters",
"counters", len(counters.keys()), 0.5, len(counters.keys())+0.5)
84 figures[
"TotalPOT"] = ROOT.TH1F(
"meta-TotalPOT",
"TotalPOT", 1, 0., 1)
86 print "run : --- Loading CAF libraries" 87 ROOT.gSystem.Load(
"libCintex.so")
89 ROOT.gSystem.Load(
"$SRT_PUBLIC_CONTEXT/lib/Linux2.6-GCC-maxopt/libStandardRecord_dict.so")
91 print "run : --- Input file" 92 if options.input_files:
93 if "*" in options.input_files:
95 print "run : globing input files" 96 files = glob.glob(options.input_files)
98 files = options.input_files.split(
",")
99 print "run : using command line list of %i input files"%(len(files))
101 assert len(files),
"run : no input files found" 103 r_chain_name =
"recTree" 104 r_chain = ROOT.TChain(r_chain_name)
106 s_chain_name =
"spillTree" 107 s_chain = ROOT.TChain(s_chain_name)
109 for this_file
in files:
110 print "run : - %s"%this_file
111 r_chain.Add(this_file)
112 s_chain.Add(this_file)
115 assert False,
"run : no input file given" 116 print "run : All files open" 117 nevents = r_chain.GetEntries()
118 s_nevents = s_chain.GetEntries()
119 counters[
"spills"] = s_nevents
120 print "run : Dataset: %i file(s)"%len(files)
121 print "run : rec tree has {:,} events".
format(nevents)
122 print "run : spill tree has {:,} events".
format(s_nevents)
127 nevents = options.n_events
128 print "run : running on {:,} event(s)".
format(nevents)
137 "energy.numusimp.trkccE",
138 "energy.numusimp.hadcalE",
139 "energy.numusimp.hadtrkE",
149 for a
in tools.fa_nu_mu_vars: active_branches.append(a)
150 r_chain.SetBranchStatus(
"*",0)
151 for b
in active_branches: r_chain.SetBranchStatus(b,1)
153 print "run : --- Beginning Analysis" 154 print "run : Beginning event loop" 155 start_time = time.time()
156 if options.verbose: print_every_n_events = 1
157 elif nevents < 1000000: print_every_n_events = 1000
158 else: print_every_n_events = 10000
159 for event
in xrange(nevents):
161 counters[
"events"] += 1
163 if r_chain.GetEntry(event) <=0:
print "run : couldn't get rec entry: ",event
164 if event%print_every_n_events == 0:
165 dt = time.time() - start_time
166 print 'run : --- Event %10s (%4.1f%%), %6.1f s'%\
169 print "run : Run: %i, sub-run: %i, event number: %i, detector: %i"%\
170 (r_chain.rec.hdr.run,r_chain.rec.hdr.subrun,r_chain.rec.hdr.evt, r_chain.rec.hdr.det)
171 if (options.mc_mode): tools.dumpTruth(r_chain)
174 pass_numu_cuts, counters = tools.numuFACuts(r_chain, counters=counters)
175 if (
not pass_numu_cuts):
continue 178 counters[
"pass numu cuts"] += 1
184 values[
"reco_neutrino_energy"] = r_chain.rec.energy.numusimp.trkccE
185 values[
"reco_muon_length"] = r_chain.rec.sel.remid.len
186 values[
"reco_hadronic_energy"] = r_chain.rec.energy.numusimp.hadcalE+r_chain.rec.energy.numusimp.hadtrkE
187 values[
"reco_fit_hadronic_energy"] = r_chain.rec.energy.numusimp.trkccE-r_chain.rec.energy.mutrkE[r_chain.rec.sel.remid.bestidx].E
192 true_event_types = tools.DetermineTruthEventType(r_chain, verbose=options.verbose)
194 if (len(r_chain.rec.mc.nu) == 1):
195 true_neutrino = r_chain.rec.mc.nu[0]
198 for true_particle
in true_neutrino.prim:
200 if true_particle.pdg
in [13,-13]:
201 if (
not true_muon)
or (true_particle.p.E() > true_muon.p.E()): true_muon = true_particle
206 counters[
"pass true cuts"] += 1
209 values[
"true_neutrino_energy"] = true_neutrino.E
210 values[
"true_muon_energy"] = true_muon.p.E()
211 values[
"true_hadronic_energy"] = true_neutrino.E - true_muon.p.E()
212 values[
"true_x"] = true_neutrino.x
213 values[
"true_y"] = true_neutrino.y
214 values[
"true_w2"] = true_neutrino.W2
215 values[
"true_q2"] = true_neutrino.q2
217 true_event_types = [
""]
219 for figure
in figure_keys:
221 x_variable = figure.split(
"_vs")[0]
222 y_variable = figure.split(
"vs_")[1]
223 if (x_variable
not in values.keys()) :
continue 224 if (y_variable
not in values.keys()) :
continue 225 for event_type
in true_event_types:
226 figures[
"%s%s"%(figure,event_type)].
Fill(values[x_variable],values[y_variable],weight)
228 if (figure
not in values.keys()) :
continue 229 for event_type
in true_event_types:
230 figures[
"%s%s"%(figure,event_type)].
Fill(values[figure],weight)
232 print "run : --- Event loop finished after {:,} events".
format(nevents)
233 dt = time.time() - start_time
234 print "run : time: %.1f s, %.1f Hz"%(dt,nevents/dt)
236 bad_runs_file =
"/nova/app/users/tamsett/art/python/bad_timing_runs.txt" 237 bad_runs_list = set([
int(r.strip())
for r
in open(bad_runs_file,
"r").readlines()]) 238 print "run : found %i runs with bad timing"%len(bad_runs_list)
240 print "run : --- Looping over spill tree" 242 start_time = time.time()
243 if options.verbose: print_every_n_events = 1
244 elif s_nevents < 1000000: print_every_n_events = 1000
245 else: print_every_n_events = 10000
246 for event
in xrange(s_nevents):
247 if options.skip_spills:
break 248 if s_chain.GetEntry(event) <=0:
print "run : couldn't get spill entry: ",event
249 if event%print_every_n_events == 0:
250 dt = time.time() - start_time
251 print 'run : --- Event %10s (%4.1f%%), %6.1f s'%\
253 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)
254 trigger_int =
int(repr(s_chain.spill.trigger)[3:-1])
255 if (trigger_int != 0):
256 print "run : Non-beam trigger found: %s"%repr(s_chain.spill.trigger)
257 print "run : Run: %i, sub-run: %i, event number: %i"%(s_chain.spill.run,s_chain.spill.subrun,s_chain.spill.evt)
261 timing_dq = (s_chain.spill.run
not in bad_runs_list)
262 if not (timing_dq):
continue 264 spill = tools.standardSpillCut(s_chain.spill)
265 if not (spill):
continue 267 total_pot += s_chain.spill.spillpot
269 print "run : --- Event loop finished after {:,} events".
format(nevents)
270 print "run : total PoT: ",total_pot
271 dt = time.time() - start_time
272 print "run : time: %.1f s, %.1f Hz"%(dt,nevents/dt)
274 print "run : --- counters" 275 for i,counter
in enumerate(counters.keys()):
277 figures[
"counters"].
GetXaxis().SetBinLabel(i+1,counter)
278 print "run : %-25s: %s"%(counter,
"{:,}".
format(counters[counter]))
280 print "run : Writing to: %s"%(options.output_file)
281 saved_root_file = ROOT.TFile(options.output_file,
"RECREATE" )
283 for i,f
in enumerate(figures.keys()):
284 if ((
"true" in f)
and (
"reco" in f)):
292 figures[f].SetName(prefix+figures[f].
GetName())
294 print "run : Wrote %i histograms"%(len(figures.keys()))
295 saved_root_file.Close()
297 print "run : --- CAF Analysis finished\n"
std::string GetName(int i)
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
correl_xv GetXaxis() -> SetDecimals()
std::string format(const int32_t &value, const int &ndigits=8)
procfile open("FD_BRL_v0.txt")
cout<< "--"<< endl;for(Int_t iP=1;iP<=hyz->GetNbinsX();iP++){for(Int_t iC=1;iC<=hyz->GetNbinsY();iC++){if(hyv->GetBinContent(iP, iC)>-999){goal_hyv-> SetBinContent(iP, iC,-(dy[iP-1][iC-1]))