std_candles.py
Go to the documentation of this file.
1 """
2  std_candles.py:
3  Plots for energy standard candles studies.
4 
5  Requires ntuple(s) created by the EnergyStandardCandles_module.cc in this package.
6 
7  Original author: J. Wolcott <jwolcott@fnal.gov>
8  Date: December 2015
9 
10 """
11 
12 import array
13 import itertools
14 import numbers
15 import os, os.path
16 import random
17 import sys
18 import types
19 
20 import ROOT
21 ROOT.gROOT.SetBatch(True)
22 
23 # the ntuples contain custom objects. can't load them without this.
24 # (cafe does basically the same thing for CAFs.)
25 import PyCintex
26 ROOT.gSystem.Load("libCalibration_dict") # in the same package as this macro, so should be ok without modifying LD_LIBRARY_PATH further
27 
28 sys.path.append("/nova/app/users/jwolcott/analysis") # so that when I forget to add this to a package later somebody else can still find it
29 from plot_tools.NtupleUtils import TreeProcessor
30 from plot_tools import HistoTools
31 
32 TESTING = False
33 
34 TREE_NAME = "stdcandles/MichelE"
35 
36 OUT_LOCATION = "/nova/ana/users/jwolcott/stdcandles"
37 OUTPUT_TYPES = ["eps", "png", "root"]
38 
39 TRUE_ORIGIN_PDGS = [11, 13, 111, 211, 2112, 2212, 0]
40 #TRUE_ORIGIN_PDGS = [11, 2112, 0]
41 TRUE_ORIGIN_COLORS = {
42  11: ROOT.kRed+1,
43  13: ROOT.kBlue+1,
44  211: ROOT.kGreen+2,
45  111: ROOT.kViolet+2,
46  2112: ROOT.kOrange-8,
47  2212: ROOT.kBlack,
48  0: ROOT.kGray,
49 }
50 LEGEND_ENTRIES = {
51  11: "Michel e^{#pm}",
52  13: "#mu^{#pm}",
53  111: "#gamma from #pi^{0}",
54  211: "#pi^{#pm}",
55  2112: "#gamma from neutron capture",
56  2212: "p^{+}",
57  0: "Other",
58 }
59 
60 NUCLEUS_LABELS = {
61  0: "Other",
62  1000170350: "^{35}Cl",
63  1000220480: "^{48}Ti",
64 }
65 
66 NUCLEUS_COLORS = {
67  0: ROOT.kBlack,
68  1000160320: ROOT.kMagenta+4,
69  1000170350: ROOT.kBlue,
70  1000170370: ROOT.kRed,
71  1000220480: ROOT.kOrange+2,
72  1000220490: ROOT.kGreen+2,
73 }
74 
75 #PARENTFRONTZ_RANGE = (400, float("inf"))
76 PARENTFRONTZ_RANGE = None
77 
78 NEUTRON_REWEIGHT = 1 # 0.58
79 
80 MICHELID_RANGE = (-9, float("inf")) # keep everything except the junk
81 #MICHELID_RANGE = (0, float("inf")) # good quality Michels only
82 #MICHELID_RANGE = (-6, -2) # neutron capture rich sample
83 
84 
85 ENERGY_BINS = [0.5*i for i in range(120)] # + range(40, 100)
86 ZCOORD_BINS = [0,] + range(10, 1160, 50)
87 
88 class StdCandlesPlotter(TreeProcessor):
89  def __init__(self, testing=False):
90  if "pc_" in NTUPLE_FILELIST:
91  suffix = "PC"
92  elif "MC" in NTUPLE_FILELIST:
93  suffix = "MC"
94  else:
95  suffix = "data"
96  name = "StdCandles_" + suffix
97  max_nevts = 1000 if testing else None
98  super(StdCandlesPlotter, self).__init__(tree_name=TREE_NAME, file_list=NTUPLE_FILELIST, name=name, max_nevts=max_nevts)
99 
100  self.POT = 0
101 
102  def RegisterBranches(self):
104  "run",
105  "subRun",
106  "spillNum",
107  "michel_candidate",
108  ]
109 
110  def RegisterMethods(self):
111  self._initialize_methods += [ self.Warnings, self.PrepareMichelEPlots ]
112  self._process_methods.append( self.FillMichelEPlots )
113  self._finalize_methods += [ self.FinalizeMichelEPlots, self.StorePOT ]
114 
115  self._tree_methods.append( self.AddPOT )
116 
117  # --------------------
118  def IsMC(self):
119  return "MC" in self.name or "PC" in self.name
120 
121  def SavePlots(self):
122  c = ROOT.TCanvas()
123  for name, plotobj in self.michel_histos.iteritems():
124  c.Clear()
125  plotobj.Draw()
126  obj_to_print = plotobj if isinstance(plotobj, ROOT.TCanvas) else c
127  for out_type in OUTPUT_TYPES:
128  obj_to_print.Print( os.path.join(OUT_LOCATION, self.name.replace("StdCandles_", ""), "MichelEPlots", "%s.%s" % (name, out_type)) )
129 
130  def Warnings(self, tree):
131  if NEUTRON_REWEIGHT != 1:
132  print >> sys.stderr, "\nWARNING: neutron captures being reweighted by factor:", NEUTRON_REWEIGHT
133 
134  if PARENTFRONTZ_RANGE is not None:
135  print >> sys.stderr, "\nWARNING: requiring parent tracks to have z coordinates in range", PARENTFRONTZ_RANGE
136 
137  if MICHELID_RANGE is not None:
138  print "\nNote: Michel PID range:", MICHELID_RANGE
139 
140  print
141 
142  # --------------------
143 
144  # --------------------
145 
146  def AddPOT(self):
147 # print dir(self.fChain)
148  tfile = self.fChain.GetCurrentFile()
149  POT_param = tfile.Get("stdcandles/POT")
150  if not POT_param:
151  print >> sys.stderr, "WARNING: no POT information for file:", tfile.GetName()
152  return
153  self.POT += POT_param.GetVal()
154 
155  def StorePOT(self):
156  if not self.IsMC():
157  self.POT *= 1e12 # since the data ART files store it in units of 10^12
158  print "Total POT:", self.POT
159  objs_to_tag = []
160  for plotname, plotobj in self.michel_histos.iteritems():
161  if hasattr(plotobj, "GetListOfPrimitives"):
162  for obj in plotobj.GetListOfPrimitives():
163  if hasattr(obj, "GetListOfFunctions"):
164  objs_to_tag.append(obj)
165  elif hasattr(obj, "GetHists"):
166  for o in obj.GetHists():
167  if hasattr(o, "GetListOfFunctions"):
168  objs_to_tag.append(o)
169  continue
170  elif hasattr(plotobj, "GetHists"):
171  for obj in obj.GetHists():
172  if hasattr(obj, "GetListOfFunctions"):
173  objs_to_tag.append(obj)
174  continue
175 
176  if not hasattr(plotobj, "GetListOfFunctions"):
177  print >> sys.stderr, "WARNING: can't store POT in plot '%s'" % plotname
178  continue
179 
180  objs_to_tag.append(plotobj)
181 
182 # print "Saving POT in objs:", objs_to_tag
183  for obj in objs_to_tag:
184  param = ROOT.TParameter("double")("POT", self.POT)
185  ROOT.SetOwnership(param, False) # the plotobj's destructor will want to delete this, so don't let it disappear when it goes out of scope
186  obj.GetListOfFunctions().Add(param)
187 
188 
189  def PrepareMichelEPlots(self, chain):
190  self.michel_histos = {}
191 
192  self.michel_histos["MichelEnergy"] = ROOT.TH1D("michel_energy", ";Candidate Michel electron E_{reco} (MeV);Events/MeV;", len(ENERGY_BINS)-1, array.array('d', ENERGY_BINS))
193  self.michel_histos["MichelPID"] = ROOT.TH1D("michel_pid", ";Michel PID score;Events/MeV;", 120, -6, 6)
194 # self.michel_histos["MichelPIDVsDeltaT"] = ROOT.TH2D("michel_pid_vs_deltaT", ";#Delta T (ns);Michel PID score;", 20, 0, 18000, 120, -6, 6)
195 # self.michel_histos["MichelPIDVsDeltaT"].SetOption("colz")
196  self.michel_histos["DeltaT"] = ROOT.TH1D("time_diff", ";Candidate Michel #Delta_{t} (ns);Events;", 20, 0, 18000)
197  self.michel_histos["NHits"] = ROOT.TH1D("num_hits", ";N_{hit} in Michel candidates (MeV);Events;", 15, 0, 15)
198  self.michel_histos["HitE"] = ROOT.TH1D("energy_per_hit", ";All E_{hit} in Michel candidates (MeV);Events;", 100, 0, 20)
199  self.michel_histos["MeanHitE"] = ROOT.TH1D("mean_energy_per_hit", ";Candidate Michel electron E_{reco}/N_{hits} (MeV);Events;", 100, 0, 20)
200  self.michel_histos["HitTimeRMS"] = ROOT.TH1D("hit_time_RMS", ";RMS of hit times within candidate Michel (ns);Events;", 50, 0, 25)
201  self.michel_histos["ParentFrontZ"] = ROOT.TH1D("parent_front_Z", ";z-coord of parent track start (cm);Events/cm;", len(ZCOORD_BINS)-1, array.array('d', ZCOORD_BINS))
202 
203  if self.IsMC():
204  self.michel_histos["TrueMichelEnergy"] = ROOT.TH1D("true_michel_energy", ";Michel electron E_{true} (MeV);Events;", 100, 0, 100)
205  self.michel_histos["TrueMissingE_NonBrem"] = ROOT.TH1D("true_missing_E_nonbrem", ";Michel E_{true} - E_{'captured'} - Brem. E_{true} (MeV);Events;", 100, 0, 100)
206  self.michel_histos["TrueMissingE"] = ROOT.TH1D("true_missing_E", ";Michel E_{true} - E_{'captured'} - Brem. E_{true} - E_{inactive} (MeV);Events;", 100, -50, 50)
207  self.michel_histos["TrueDeadMaterialE"] = ROOT.TH1D("true_E_in_dead_materials", ";Estimated true E_{Michel} deposited in inactive materials (MeV);Events;", 100, 0, 100)
208  self.michel_histos["TrueMichelFracCaptured"] = ROOT.TH1D("michel_frac_captured", ";Fraction of true Michels' energy captured;Events;", 101, 0, 1.01)
209  self.michel_histos["TrueELostToBremVsRecoE"] = ROOT.TH2D("michel_E_brem", ";Michel E_{reco};True Michel energy lost to Brem. (MeV);", 100, 0, 100, 100, 0, 100)
210 
212  for plotname, plot in self.michel_histos.iteritems():
213  if "TrueMichel" in plotname or isinstance(plot, ROOT.TH2D):
214  continue
215  self.michel_histos_by_pdg[plotname] = {}
216  for pdg in TRUE_ORIGIN_PDGS:
217  h = ROOT.TH1D(plot)
218  h.SetName(h.GetName() + "_TruePID_%d" % pdg)
219  h.SetLineColor(TRUE_ORIGIN_COLORS[pdg])
220  self.michel_histos_by_pdg[plotname][pdg] = h
221 
222  self.michel_origins = []
224 
225  def FillMichelEPlots(self):
226  tree = self.fChain # save some typing
227 
228  michel_candidate = tree.michel_candidate # slightly less typing, and one less reference for Python to traverse
229 
230  if PARENTFRONTZ_RANGE is not None and not (PARENTFRONTZ_RANGE[0] <= michel_candidate.parentFrontZ <= PARENTFRONTZ_RANGE[1]):
231  return
232 
233  if MICHELID_RANGE is not None and not (MICHELID_RANGE[0] <= michel_candidate.MID < MICHELID_RANGE[1]):
234  return
235 
236  energy = sum(michel_candidate.hitEnergies)
237  nHits = len(michel_candidate.hitEnergies)
238  if energy <= 0:
239  return
240 
241  vals = {
242  "MichelEnergy": energy,
243  "MichelPID": michel_candidate.MID,
244  "NHits": nHits,
245  "HitE": michel_candidate.hitEnergies,
246  "DeltaT": michel_candidate.deltaT,
247  "MeanHitE": energy/nHits,
248  "HitTimeRMS": michel_candidate.t_RMS,
249  "ParentFrontZ": michel_candidate.parentFrontZ,
250  }
251 
252  wgt = 1
253 
254 # print tree.energy
255 # print len(michel_candidate.eContribs)
256  if self.IsMC() and len(michel_candidate.eContribs) > 0:
257  origin_breakdown = {}
258  for contrib in michel_candidate.eContribs:
259  if hasattr(contrib, "ElostToBrem"):
260  self.michel_histos["TrueMichelFracCaptured"].Fill(contrib.fracCaptured)
261  self.michel_histos["TrueMichelEnergy"].Fill(contrib.trueFourP.E())
262 
263  if contrib.contribEFrac > 0.9:
264  self.michel_histos["TrueMissingE_NonBrem"].Fill(contrib.trueFourP.E() - contrib.trueFourP.E()*contrib.fracCaptured - contrib.ElostToBrem)
265  self.michel_histos["TrueMissingE"].Fill(contrib.trueFourP.E() - contrib.trueFourP.E()*contrib.fracCaptured - contrib.ElostToBrem - contrib.ElostToDeadMaterials)
266  self.michel_histos["TrueDeadMaterialE"].Fill(contrib.ElostToDeadMaterials)
267  self.michel_histos["TrueELostToBremVsRecoE"].Fill(energy, contrib.ElostToBrem)
268 
269  pdg = abs(contrib.pdg)
270  if pdg not in TRUE_ORIGIN_PDGS:
271  pdg = 0
272 # print pdg, contrib.contribEFrac
273  origin_breakdown[pdg] = origin_breakdown.setdefault(pdg, 0) + contrib.contribEFrac
274 
275  if hasattr(contrib, "nCaptureNucleus") and contrib.nCaptureNucleus != 0:
276  self.neutron_captures[contrib.nCaptureNucleus] = self.neutron_captures.setdefault(contrib.nCaptureNucleus, 0) + 1
277 
278  total_frac = sum(origin_breakdown.values())
279 # print origin_breakdown, total_frac
280  assert abs(total_frac - 1) < 1e-5, "total fraction not 1! sum: %f" % total_frac
281 
282  self.michel_origins.append(origin_breakdown)
283 
284  for pdg, frac in origin_breakdown.iteritems():
285  if frac > 0.75:
286  if pdg == 2112:
287  wgt = NEUTRON_REWEIGHT
288  for plotname, val in vals.iteritems():
289  if plotname not in self.michel_histos:
290  continue
291  try:
292  vs = list(val)
293  except TypeError:
294  vs = [val,]
295  for v in vs:
296  self.michel_histos_by_pdg[plotname][pdg].Fill(v, wgt)
297 
298 # if set(origin_breakdown.keys()) != set([11,]):
299 # return
300 
301 # print tree.energy
302  for plotname, val in vals.iteritems():
303  if plotname not in self.michel_histos:
304  continue
305  try:
306  vs = list(val)
307  except TypeError:
308  vs = [val,]
309  for v in vs:
310  self.michel_histos[plotname].Fill(v, wgt)
311 
312 # self.michel_histos["MichelPIDVsDeltaT"].Fill(michel_candidate.deltaT, michel_candidate.MID)
313 
314 
316  hists = self.michel_histos.values()
317  if self.IsMC():
318 # print "adding hists:", self.michel_histos_by_pdg.values()
319  for coll in self.michel_histos_by_pdg.itervalues():
320  hists += coll.values()
321 # print coll, hists
322  for p in hists:
323 # print p
324  axis_variable = [getattr(p, "Get%saxis" % axis)().IsVariableBinSize() for axis in ("X", "Y")]
325  if any(axis_variable):
326  p.Scale(1, "width")
327 # print "Scaled plot:", p.GetName()
328 
329  if self.IsMC():
330  self.MakeBremPlots()
331  self.MakeOverlayMCPlots()
334 
335  # ---------------------
336 
337  def FillMichelBreakdownPlots(self, histo_stack):
338  fracs = {}
339  nCurrEvts = 0
340  nBin = 1
341  n_bins = min(200, len(self.michel_origins))
342  n_evts_in_bin = float(len(self.michel_origins)) / n_bins
343  nEvtsInBin_int = int(n_evts_in_bin)
344  nEvtsInBin_frac = n_evts_in_bin - nEvtsInBin_int
345  for breakdown in self.michel_origins:
346  # if this is all the events that should go into this bin,
347  # average and insert into the plot
348  if nCurrEvts >= nEvtsInBin_int + int(random.random() < nEvtsInBin_frac):
349  for pdg in TRUE_ORIGIN_PDGS:
350  frac = 0
351  if pdg in fracs:
352  frac = fracs[pdg] / float(nCurrEvts)
353  histo_stack[pdg].SetBinContent(nBin, frac)
354 # print " PDG %d: setting bin %d to have content %f" % (pdg, nBin, frac)
355  fracs = {}
356  nCurrEvts = 0
357  nBin += 1
358 
359 # print breakdown
360  for pdg, frac in breakdown.iteritems():
361 # LOG_DEBUG("StandardCandlePlots") << " PDG " << pdgFracPair.first << " contributed fraction " << pdgFracPair.second << " to Michel";
362  fracs[pdg] = fracs.setdefault(pdg, 0) + frac
363 # LOG_DEBUG("StandardCandlePlots") << " fracs is now: ";
364 # for (const auto & f : fracs)
365 # LOG_DEBUG("StandardCandlePlots") << " PDG " << f.first << ": " << f.second;
366  nCurrEvts += 1;
367 
368  def MakeBremPlots(self):
369  twoD_plot = self.michel_histos["TrueELostToBremVsRecoE"]
370  ROOT.SetOwnership(twoD_plot, False)
371  canvas = ROOT.TCanvas()
372  canvas.cd()
373  twoD_plot.Draw("colz")
374  prof = twoD_plot.ProfileX()
375  prof.Draw("same")
376  self.michel_histos["TrueELostToBremVsRecoE"] = canvas
377 
378  proj = twoD_plot.ProjectionY("TrueELostToBrem")
379  self.michel_histos["TrueELostToBrem"] = proj
380 
382 # print self.neutron_captures
383  key_order = NUCLEUS_LABELS.keys()
384  vals = dict( [(k, self.neutron_captures[k]) for k in key_order if k in self.neutron_captures])
385  vals[0] = sum([v for k, v in self.neutron_captures.iteritems() if k not in vals])
386  colors = array.array('i', [NUCLEUS_COLORS[pdg] for pdg in key_order])
387  vals = array.array('d', [vals[k] for k in key_order if k in vals])
388  n_nuclei = len(vals)
389  pie = ROOT.TPie("NeutronCaptureNucleus", "", n_nuclei, vals, colors)
390  for i, pdg in enumerate(key_order):
391  pie.SetEntryLabel(i, NUCLEUS_LABELS[pdg])
392  pie.SetRadius(0.2)
393  pie.SetLabelFormat("#splitline{%txt}{%perc}");
394  self.michel_histos["NeutronCaptureNucleus"] = pie
395 
397  for plotname, plot in self.michel_histos.iteritems():
398  if "True" in plotname:
399  continue
400 
401  ROOT.SetOwnership(plot, False)
402 
403  legend = ROOT.TLegend(0.55, 0.6, 0.83, 0.9)
404  ROOT.SetOwnership(legend, False)
405 
406  keys = self.michel_histos_by_pdg[plotname].keys()
407  total_exclusive = self.michel_histos_by_pdg[plotname][keys[0]].Clone("%s_sum_exclusive" % plotname)
408  ROOT.SetOwnership(total_exclusive, False)
409  for k in keys[1:]:
410  total_exclusive.Add(self.michel_histos_by_pdg[plotname][k])
411  total_exclusive.SetLineColor(ROOT.kBlack)
412  total_exclusive.SetLineStyle(3)
413 
414  h_by_pdg = ROOT.THStack()
415  ROOT.SetOwnership(h_by_pdg, False)
416  hists_to_consider = {"Total": plot, "Total-exclusive": total_exclusive}
417  hists_to_consider.update(self.michel_histos_by_pdg[plotname])
418  for pdg in ["Total", "Total-exclusive"] + TRUE_ORIGIN_PDGS:
419  h = hists_to_consider[pdg]
420 # print " adding histogram:", h
421  h_by_pdg.Add(h, "hist")
422  if pdg == "Total":
423  legend_entry = "Total"
424  elif pdg == "Total-exclusive":
425  legend_entry = "Total of exclusive categories"
426  else:
427  legend_entry = LEGEND_ENTRIES[pdg]
428  legend.AddEntry(h, legend_entry, "l")
429  c = ROOT.TCanvas()
430  h_by_pdg.Draw("nostack")
431  HistoTools.SanitizeStackLabels(h_by_pdg)
432 
433  # legend winds up overlapping with distributions in this one plot.
434  # not a wonderful place for special cases but will move if it gets unwieldy
435  if "MichelPID" in plotname:
436  print "fixing Michel PID plot legend"
437  legend.SetX1(0.2)
438  legend.SetX2(0.5)
439  legend.Draw()
440 
441  self.michel_histos[plotname] = c
442 
444  legend = ROOT.TLegend(0.55, 0.5, 0.83, 0.8)
445  legend.SetFillColor(ROOT.kWhite)
446  legend.SetBorderSize(1)
447  ROOT.SetOwnership(legend, False)
448 
449  n_bins = min(200, len(self.michel_origins))
450  origin_by_pdg = {}
451  for pdg in TRUE_ORIGIN_PDGS:
452  h = ROOT.TH1D("true_origin_%d" % pdg, ";Fraction of events;Fraction of true energy;", n_bins, 0, 1)
453  h.SetLineColor(TRUE_ORIGIN_COLORS[pdg])
454  h.SetFillColor(TRUE_ORIGIN_COLORS[pdg])
455  origin_by_pdg[pdg] = h
456 
457  legend.AddEntry(h, LEGEND_ENTRIES[pdg], "f")
458 
459  self.michel_origins.sort(
460  # yikes.
461  # basically: sort by amount of each PDG in priority order of TRUE_ORIGIN_PDGS
462  key=lambda bkdown:
463  sum(100**(len(TRUE_ORIGIN_PDGS)-i)*(bkdown[pdg] if pdg in bkdown else 0) for i, pdg in enumerate(TRUE_ORIGIN_PDGS))
464  )
465 
466  self.FillMichelBreakdownPlots(origin_by_pdg)
467 
468  st = ROOT.THStack()
469  ROOT.SetOwnership(st, False)
470  for pdg in TRUE_ORIGIN_PDGS:
471  ROOT.SetOwnership(origin_by_pdg[pdg], False)
472  st.Add(origin_by_pdg[pdg], "hist")
473 
474  c = ROOT.TCanvas()
475  st.Draw()
476  legend.Draw()
477  HistoTools.SanitizeStackLabels(st)
478  self.michel_histos["TrueMichelOriginBreakdown"] = c
479 
480 # -----------------
481 
483  AREA_NORM = object()
484  PEAK_NORM = object()
485  POT_NORM = object()
486  def __init__(self, name, path, mc_norm_factor=AREA_NORM):
487  self.name = name
488  self.path = path
489  self.mc_norm_factor = mc_norm_factor
490 
491 COMBINATION_PLOTS = [
492  PlotConfig( name="MichelEnergy", path="MichelEPlots/MichelEnergy", mc_norm_factor=PlotConfig.POT_NORM ),
493  PlotConfig( name="MichelPID", path="MichelEPlots/MichelPID", mc_norm_factor=PlotConfig.POT_NORM ),
494  PlotConfig( name="DeltaT", path="MichelEPlots/DeltaT", mc_norm_factor=PlotConfig.POT_NORM ),
495  PlotConfig( name="NHits", path="MichelEPlots/NHits", mc_norm_factor=PlotConfig.POT_NORM ),
496  PlotConfig( name="HitE", path="MichelEPlots/HitE", mc_norm_factor=PlotConfig.POT_NORM ),
497  PlotConfig( name="HitTimeRMS", path="MichelEPlots/HitTimeRMS", mc_norm_factor=PlotConfig.POT_NORM ),
498  PlotConfig (name="MeanHitE", path="MichelEPlots/MeanHitE", mc_norm_factor=PlotConfig.POT_NORM ),
499  PlotConfig (name="ParentFrontZ", path="MichelEPlots/ParentFrontZ", mc_norm_factor=PlotConfig.POT_NORM ),
500 ]
501 
503  def __init__(self):
504  self.data_plots_in = {}
505  self.mc_plots_in = {}
506 
507  self.plots_out = {}
508 
509  self.LoadPlots()
510 
511  @staticmethod
512  def GetMCDataStack(input):
513  stack = None
514  if hasattr(input, "GetListOfPrimitives"):
515  stacks = [o for o in input.GetListOfPrimitives() if hasattr(o, "GetHists")]
516  assert len(stacks) == 1, "Too many THStacks in collection '%s' for me to deduce which one you want!" % input.GetName()
517  stack = stacks[0]
518  elif hasattr(input, "GetHists"):
519  stack = input
520  assert stack, "Could not determine the stack to pull the sum out of!"
521 
522  return stack
523 
524  def LoadPlots(self):
525  for plot_config in COMBINATION_PLOTS:
526  for source in ("data", "mc"):
527  collection = getattr(self, "%s_plots_in" % source)
528 
529  # get the capitalization right. should've been more consistent
530  source_infix = source
531  if source == "mc":
532  source_infix = "MC"
533  plotfile_name = os.path.join(OUT_LOCATION, source_infix, "%s.root" % plot_config.path)
534  if not os.path.isfile(plotfile_name) :
535  raise Exception("Could not find necessary input file: '%s'" % plotfile_name)
536 
537  f = ROOT.TFile(plotfile_name)
538  keys = f.GetListOfKeys()
539  if len(keys) > 1:
540  raise Exception("Don't know which of the %d objects in file you want: %s" % (len(objs), list(keys)))
541  elif len(keys) == 0:
542  raise Exception("File has nothing in it: '%s'" % plotname_file)
543 
544  obj = keys.First().ReadObj()
545  assert obj, "Object read from file ('%s') is null!" % plotfile_name
546  if hasattr(obj, "SetDirectory"):
547  obj.SetDirectory(None) # otherwise they disappear when the file is closed
548  collection[plot_config.path] = obj
549 
550  def MakePlots(self):
551  for plot_config in COMBINATION_PLOTS:
552  method_name = "Make_%s_Plot" % plot_config.name
553  if hasattr(self, method_name) or self.Make_Default_Plot_Method(plot_config.name, plot_config.path, plot_config.mc_norm_factor):
554  getattr(self, method_name)()
555  else:
556  raise Exception("Don't know how to make '%s' plot (no method '%s' defined)" % (plot_config.name, method_name))
557 
558  def Make_DataMC_Overlay(self, in_varname):
559  mc_canvas = self.mc_plots_in[in_varname]
560  assert isinstance(mc_canvas, ROOT.TCanvas)
561  mc_stack = [o for o in mc_canvas.GetListOfPrimitives() if isinstance(o, ROOT.THStack)][0]
562 
563  legend = ROOT.TLegend(0.65, 0.7, 0.85, 0.9)
564  ROOT.SetOwnership(legend, False)
565 
566  data_canvas = self.data_plots_in[in_varname]
567  data_plot = [o for o in data_canvas.GetListOfPrimitives() if isinstance(o, ROOT.TH1)][0]
568  data_plot.SetMarkerSize(0)
569  legend.AddEntry(data_plot, "Data", "le")
570 
571  # we're just going to do data-MC total
572  # comparisons (to avoid over-complicating the figure).
573  out_stack = ROOT.THStack(mc_stack)
574  ROOT.SetOwnership(out_stack, False)
575  for h in list(out_stack.GetHists()):
576  if "TruePID" in h.GetName() or "sum_exclusive" in h.GetName():
577  out_stack.GetHists().Remove(h)
578  else:
579  h.SetLineColor(ROOT.kRed)
580  legend.AddEntry(h, "MC", "l")
581 
582  out_stack.Add(data_plot, "pe")
583 
584  c = ROOT.TCanvas()
585  c.cd()
586  ROOT.SetOwnership(c, False)
587  out_stack.Draw("nostack")
588  HistoTools.SanitizeStackLabels(out_stack)
589  legend.Draw()
590  return c
591 
592  @staticmethod
593  def NormalizeMCPlot(data_plot, mc_plot, mc_norm_factor=PlotConfig.AREA_NORM):
594 
595  assert data_plot and hasattr(mc_plot, "Scale")
596 
597  scale_factor = None
598  if isinstance(mc_norm_factor, numbers.Number) and mc_norm_factor >= 0:
599  scale_factor = mc_norm_factor
600  elif mc_norm_factor is PlotConfig.AREA_NORM:
601  scale_factor = data_plot.Integral() / mc_plot.Integral()
602  elif mc_norm_factor is PlotConfig.PEAK_NORM:
603  scale_factor = data_plot.GetMaximum() / mc_plot.GetMaximum()
604  elif mc_norm_factor is PlotConfig.POT_NORM:
605  data_pot = data_plot.GetListOfFunctions().FindObject("POT")
606  mc_pot = mc_plot.GetListOfFunctions().FindObject("POT")
607  assert hasattr(data_pot, "GetVal"), "Couldn't get POT from data plot: '%s'" % data_plot.GetName()
608  assert hasattr(mc_pot, "GetVal"), "Could't get POT values from MC plot: '%s'" % mc_plot.GetName()
609  scale_factor = data_pot.GetVal() / mc_pot.GetVal()
610 
611  assert scale_factor is not None, "Could not understand your MC normalization factor: %s" % mc_norm_factor
612 
613  mc_plot.Scale(scale_factor)
614 
615  # this will cover most cases.
616  # however, if you need especially fine-grained control,
617  # you can define a Make_<your plot name>_Plot(self)
618  # in which you can do whatever you want
619  def Make_Default_Plot_Method(self, plot_name, plot_path, mc_norm_factor):
620  def f(self):
621  canvas = self.Make_DataMC_Overlay(plot_path)
622 
623  # normalize...
624  stack = DataMCPlots.GetMCDataStack(canvas)
625  mc_plot = stack.GetHists()[0]
626  data_plot = stack.GetHists()[1]
627 
628  DataMCPlots.NormalizeMCPlot(data_plot, mc_plot, mc_norm_factor)
629 
630  self.plots_out[plot_name] = canvas
631 
632  setattr(self, "Make_%s_Plot" % plot_name, types.MethodType(f, self))
633 
634  return True
635 
637  # don't use the stock method because the legend winds up overlapping with the distribution
638  michel_config = [config for config in COMBINATION_PLOTS if config.name == "MichelPID"][0]
639  canvas = self.Make_DataMC_Overlay(michel_config.path)
640  leg = [o for o in canvas.GetListOfPrimitives() if isinstance(o, ROOT.TLegend)][0]
641  leg.SetX1(0.2)
642  leg.SetX2(0.5)
643 
644  stack = DataMCPlots.GetMCDataStack(canvas)
645  mc_plot = stack.GetHists()[0]
646  data_plot = stack.GetHists()[1]
647 
648  DataMCPlots.NormalizeMCPlot(data_plot, mc_plot, michel_config.mc_norm_factor)
649 
650  self.plots_out[michel_config.name] = canvas
651 
652 
653  def SavePlots(self):
654  c = ROOT.TCanvas()
655  for name, plotobj in self.plots_out.iteritems():
656  c.Clear()
657  plotobj.Draw()
658  obj_to_print = plotobj if isinstance(plotobj, ROOT.TCanvas) else c
659  for out_type in OUTPUT_TYPES:
660  obj_to_print.Print( os.path.join(OUT_LOCATION, "data+MC", "MichelEPlots", "%s.%s" % (name, out_type)) )
661 
662 ################################
663 NTUPLE_FILELISTS = {
664  "test": "/nova/ana/users/jwolcott/ntuple_lists/test_list_MC.txt",
665  "data": "/nova/ana/users/jwolcott/art/data/std_candles/playlist-jwolcott_data_ND_FA_small-20151208.txt",
666  "MC": "/nova/ana/users/jwolcott/art/mc/std_candles/playlist-jwolcott_MC_ND_FA_small-20151208.txt",
667  "PC": "/nova/ana/users/jwolcott/art/mc/pc_muplus_extraEM/playlist-jwolcott_pc_muplus_extraEM_small-20160108.txt",
668 }
669 
670 if __name__ == "__main__":
671  assert len(sys.argv) == 2 and (sys.argv[1] in NTUPLE_FILELISTS or sys.argv[1] == "comp"), \
672  "This script must be given exactly one argument (the name of the ntuple list to use).\nAvailable options: %s" % NTUPLE_FILELISTS.keys()
673 
674  if sys.argv[1] == "comp":
675  pl = DataMCPlots()
676  pl.MakePlots()
677  pl.SavePlots()
678  else:
679  NTUPLE_FILELIST = NTUPLE_FILELISTS[sys.argv[1]]
680  pl = StdCandlesPlotter(testing=TESTING)
681  pl.LoopTree()
682  pl.SavePlots()
keys
Reco plots.
Definition: caf_analysis.py:46
def Make_Default_Plot_Method(self, plot_name, plot_path, mc_norm_factor)
Definition: std_candles.py:619
def __init__(self, testing=False)
Definition: std_candles.py:89
void abs(TH1 *hist)
def NormalizeMCPlot(data_plot, mc_plot, mc_norm_factor=PlotConfig.AREA_NORM)
Definition: std_candles.py:593
def __init__(self, name, path, mc_norm_factor=AREA_NORM)
Definition: std_candles.py:486
Definition: novas.h:112
def Make_MichelPID_Plot(self)
Definition: std_candles.py:636
def Make_DataMC_Overlay(self, in_varname)
Definition: std_candles.py:558
def PrepareMichelEPlots(self, chain)
Definition: std_candles.py:189
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
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]))
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
Double_t sum
Definition: plot.C:31
def FillMichelBreakdownPlots(self, histo_stack)
Definition: std_candles.py:337