caf_analysis.py
Go to the documentation of this file.
1 print "\nrun : --- CAF analysis"
2 ########################### Option parser
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
21 ################################## Setup ROOT
22 print "run : Importing ROOT"
23 import ROOT
24 ROOT.gROOT.SetBatch(1)
25 ################################## Imports
26 import caf_tools as tools
27 from collections import OrderedDict
28 import time
29 ################################## Book histograms
30 print "run : --- Booking histograms"
31 figures = {}
32 #### Reco plots
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)
37 #### True plots
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.)
45 #### Multi-dimensional plots
46 keys = figures.keys()
47 # dependent_variables = ["true_neutrino_energy"]
48 # for d_v in dependent_variables:
49 # for f in keys:
50 # if f == d_v: continue
51 # x_bins = tools.GetBins(figures[d_v].GetXaxis())
52 # y_bins = tools.GetBins(figures[f].GetXaxis())
53 # figures["%s_vs_%s"%(d_v,f)] = ROOT.TH2F("%s_vs_%s"%(d_v,f), "%s_vs_%s"%(d_v,f),
54 # len(x_bins)-1, x_bins[0], x_bins[-1],
55 # len(y_bins)-1, y_bins[0], y_bins[-1])
56 #### Store plot keys
57 figure_keys = figures.keys()
58 #### Sub-samples
59 if options.mc_mode:
60  int_types = [
61  "QE","Res","DIS", "Coh","MEC",
62  ]
63  keys = figures.keys()
64  for f in keys:
65  for t in int_types:
66  name = "%s-%s"%(f,t)
67  figures[name] = figures[f].Clone(name)
68 ################################## Initialise a counter dictionary
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)
83 #### Meta histograms
84 figures["TotalPOT"] = ROOT.TH1F("meta-TotalPOT", "TotalPOT", 1, 0., 1)
85 ################################## NOvA setup
86 print "run : --- Loading CAF libraries"
87 ROOT.gSystem.Load("libCintex.so")
88 ROOT.Cintex.Enable()
89 ROOT.gSystem.Load("$SRT_PUBLIC_CONTEXT/lib/Linux2.6-GCC-maxopt/libStandardRecord_dict.so")
90 ################################## load input file
91 print "run : --- Input file"
92 if options.input_files:
93  if "*" in options.input_files:
94  import glob
95  print "run : globing input files"
96  files = glob.glob(options.input_files)
97  else:
98  files = options.input_files.split(",")
99  print "run : using command line list of %i input files"%(len(files))
100  files.sort()
101  assert len(files),"run : no input files found"
102 
103  r_chain_name = "recTree"
104  r_chain = ROOT.TChain(r_chain_name)
105 
106  s_chain_name = "spillTree"
107  s_chain = ROOT.TChain(s_chain_name)
108 
109  for this_file in files:
110  print "run : - %s"%this_file
111  r_chain.Add(this_file)
112  s_chain.Add(this_file)
113  # assert r_chain.GetEntries(), "run : No entries in input file: %s, an error has probably occured"%this_file
114 else:
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)
123 # assert nevents,"run : No entries in any input file(s), an error has probably occurred"
124 for f in files:
125  print "run : - %s"%f
126 if options.n_events:
127  nevents = options.n_events
128  print "run : running on {:,} event(s)".format(nevents)
129 ################################## Optimisations
130 active_branches = [
131  "hdr.run",
132  "hdr.subrun",
133  "hdr.evt",
134  "hdr.det",
135  "sel.remid.len",
136  "sel.remid.bestidx",
137  "energy.numusimp.trkccE",
138  "energy.numusimp.hadcalE",
139  "energy.numusimp.hadtrkE",
140  "energy.mutrkE.E",
141  "mc.nu",
142  "mc.nnu",
143  "mc.nu.iscc",
144  "mc.nu.pdg",
145  "mc.nu.pdgorig",
146  "mc.nu.E",
147  "mc.nu.*",
148 ]
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)
152 ################################## Run analysis
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):
160  ################################## Begin loop
161  counters["events"] += 1
162  ################################## Get entries
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'%\
167  ("{:,}".format(event),( (float(event)/float(nevents))*100.), dt)
168  if options.verbose:
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)
172  ################################## Nu mu analysis cuts
173  # These are the standard nu mu analysis cuts used in first analysis
174  pass_numu_cuts, counters = tools.numuFACuts(r_chain, counters=counters)
175  if (not pass_numu_cuts): continue
176 
177  # if we get here then we've passed all the nu mu cuts
178  counters["pass numu cuts"] += 1
179 
180  # end of nu mu cuts
181  ################################## Fill reconstruction values
182  weight = 1.
183  values = {}
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
188  ################################## Truth studies
189  # if the '-m' option is set then we assume this is MC
190  # if so we run some simple truth analysis
191  if options.mc_mode:
192  true_event_types = tools.DetermineTruthEventType(r_chain, verbose=options.verbose)
193  # if there is more than one truth neutrino this event is too complicated
194  if (len(r_chain.rec.mc.nu) == 1):
195  true_neutrino = r_chain.rec.mc.nu[0]
196  true_muon = None
197 
198  for true_particle in true_neutrino.prim:
199  # select the most energetic muon
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
202 
203  # if we didn't find a true muon then this event is no good
204  if true_muon:
205  # if we get here then we've passed the true cuts
206  counters["pass true cuts"] += 1
207 
208  ################################## Fill truth values
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
216  else:
217  true_event_types = [""]
218  ################################## Fill histograms
219  for figure in figure_keys:
220  if "vs" in figure:
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)
227  else:
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)
231 ################################## Event loop finished
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)
235 ################################## Bad timing list
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)
239 ################################## Spill tree
240 print "run : --- Looping over spill tree"
241 total_pot = 0.
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'%\
252  ("{:,}".format(event),( (float(event)/float(s_nevents))*100.), dt)
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)
258  ################################## Analysis cuts
259  # the full analysis cuts should be applied here
260  # timing dq
261  timing_dq = (s_chain.spill.run not in bad_runs_list)
262  if not (timing_dq): continue
263 
264  spill = tools.standardSpillCut(s_chain.spill)
265  if not (spill): continue
266 
267  total_pot += s_chain.spill.spillpot
268 figures["TotalPOT"].SetBinContent(1,total_pot)
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)
273 ################################## print counters
274 print "run : --- counters"
275 for i,counter in enumerate(counters.keys()):
276  figures["counters"].SetBinContent(i+1,counters[counter])
277  figures["counters"].GetXaxis().SetBinLabel(i+1,counter)
278  print "run : %-25s: %s"%(counter,"{:,}".format(counters[counter]))
279 ################################## Save histograms
280 print "run : Writing to: %s"%(options.output_file)
281 saved_root_file = ROOT.TFile(options.output_file, "RECREATE" )
282 saved_root_file.cd()
283 for i,f in enumerate(figures.keys()):
284  if (("true" in f) and ("reco" in f)):
285  prefix = "mixed-"
286  elif ("true" in f):
287  prefix = "true-"
288  elif ("reco" in f):
289  prefix = "reco-"
290  else:
291  prefix = ""
292  figures[f].SetName(prefix+figures[f].GetName())
293  figures[f].Write()
294 print "run : Wrote %i histograms"%(len(figures.keys()))
295 saved_root_file.Close()
296 ################################## Done
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)
Definition: HexUtils.cpp:14
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]))
gm Write()