3 Plots for energy standard candles studies. 5 Requires ntuple(s) created by the EnergyStandardCandles_module.cc in this package. 7 Original author: J. Wolcott <jwolcott@fnal.gov> 21 ROOT.gROOT.SetBatch(
True)
26 ROOT.gSystem.Load(
"libCalibration_dict")
28 sys.path.append(
"/nova/app/users/jwolcott/analysis")
29 from plot_tools.NtupleUtils
import TreeProcessor
30 from plot_tools
import HistoTools
34 TREE_NAME =
"stdcandles/MichelE" 36 OUT_LOCATION =
"/nova/ana/users/jwolcott/stdcandles" 37 OUTPUT_TYPES = [
"eps",
"png",
"root"]
39 TRUE_ORIGIN_PDGS = [11, 13, 111, 211, 2112, 2212, 0]
41 TRUE_ORIGIN_COLORS = {
53 111:
"#gamma from #pi^{0}",
55 2112:
"#gamma from neutron capture",
62 1000170350:
"^{35}Cl",
63 1000220480:
"^{48}Ti",
68 1000160320: ROOT.kMagenta+4,
69 1000170350: ROOT.kBlue,
70 1000170370: ROOT.kRed,
71 1000220480: ROOT.kOrange+2,
72 1000220490: ROOT.kGreen+2,
76 PARENTFRONTZ_RANGE =
None 80 MICHELID_RANGE = (-9,
float(
"inf"))
85 ENERGY_BINS = [0.5*i
for i
in range(120)]
86 ZCOORD_BINS = [0,] +
range(10, 1160, 50)
90 if "pc_" in NTUPLE_FILELIST:
92 elif "MC" in NTUPLE_FILELIST:
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)
115 self._tree_methods.append( self.
AddPOT )
119 return "MC" in self.name
or "PC" in self.name
123 for name, plotobj
in self.michel_histos.iteritems():
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)) )
131 if NEUTRON_REWEIGHT != 1:
132 print >> sys.stderr,
"\nWARNING: neutron captures being reweighted by factor:", NEUTRON_REWEIGHT
134 if PARENTFRONTZ_RANGE
is not None:
135 print >> sys.stderr,
"\nWARNING: requiring parent tracks to have z coordinates in range", PARENTFRONTZ_RANGE
137 if MICHELID_RANGE
is not None:
138 print "\nNote: Michel PID range:", MICHELID_RANGE
148 tfile = self.fChain.GetCurrentFile()
149 POT_param = tfile.Get(
"stdcandles/POT")
151 print >> sys.stderr,
"WARNING: no POT information for file:", tfile.GetName()
153 self.
POT += POT_param.GetVal()
158 print "Total POT:", self.
POT 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)
170 elif hasattr(plotobj,
"GetHists"):
171 for obj
in obj.GetHists():
172 if hasattr(obj,
"GetListOfFunctions"):
173 objs_to_tag.append(obj)
176 if not hasattr(plotobj,
"GetListOfFunctions"):
177 print >> sys.stderr,
"WARNING: can't store POT in plot '%s'" % plotname
180 objs_to_tag.append(plotobj)
183 for obj
in objs_to_tag:
184 param = ROOT.TParameter(
"double")(
"POT", self.
POT)
185 ROOT.SetOwnership(param,
False)
186 obj.GetListOfFunctions().
Add(param)
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)
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))
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)
212 for plotname, plot
in self.michel_histos.iteritems():
213 if "TrueMichel" in plotname
or isinstance(plot, ROOT.TH2D):
216 for pdg
in TRUE_ORIGIN_PDGS:
218 h.SetName(h.GetName() +
"_TruePID_%d" % pdg)
219 h.SetLineColor(TRUE_ORIGIN_COLORS[pdg])
228 michel_candidate = tree.michel_candidate
230 if PARENTFRONTZ_RANGE
is not None and not (PARENTFRONTZ_RANGE[0] <= michel_candidate.parentFrontZ <= PARENTFRONTZ_RANGE[1]):
233 if MICHELID_RANGE
is not None and not (MICHELID_RANGE[0] <= michel_candidate.MID < MICHELID_RANGE[1]):
236 energy =
sum(michel_candidate.hitEnergies)
237 nHits = len(michel_candidate.hitEnergies)
242 "MichelEnergy": energy,
243 "MichelPID": michel_candidate.MID,
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,
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"):
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)
269 pdg =
abs(contrib.pdg)
270 if pdg
not in TRUE_ORIGIN_PDGS:
273 origin_breakdown[pdg] = origin_breakdown.setdefault(pdg, 0) + contrib.contribEFrac
275 if hasattr(contrib,
"nCaptureNucleus")
and contrib.nCaptureNucleus != 0:
276 self.
neutron_captures[contrib.nCaptureNucleus] = self.neutron_captures.setdefault(contrib.nCaptureNucleus, 0) + 1
278 total_frac =
sum(origin_breakdown.values())
280 assert abs(total_frac - 1) < 1e-5,
"total fraction not 1! sum: %f" % total_frac
282 self.michel_origins.append(origin_breakdown)
284 for pdg, frac
in origin_breakdown.iteritems():
287 wgt = NEUTRON_REWEIGHT
288 for plotname, val
in vals.iteritems():
302 for plotname, val
in vals.iteritems():
316 hists = self.michel_histos.values()
319 for coll
in self.michel_histos_by_pdg.itervalues():
320 hists += coll.values()
324 axis_variable = [getattr(p,
"Get%saxis" % axis)().IsVariableBinSize()
for axis
in (
"X",
"Y")]
325 if any(axis_variable):
343 nEvtsInBin_int =
int(n_evts_in_bin)
344 nEvtsInBin_frac = n_evts_in_bin - nEvtsInBin_int
348 if nCurrEvts >= nEvtsInBin_int +
int(random.random() < nEvtsInBin_frac):
349 for pdg
in TRUE_ORIGIN_PDGS:
352 frac = fracs[pdg] /
float(nCurrEvts)
360 for pdg, frac
in breakdown.iteritems():
362 fracs[pdg] = fracs.setdefault(pdg, 0) + frac
370 ROOT.SetOwnership(twoD_plot,
False)
371 canvas = ROOT.TCanvas()
373 twoD_plot.Draw(
"colz")
374 prof = twoD_plot.ProfileX()
378 proj = twoD_plot.ProjectionY(
"TrueELostToBrem")
383 key_order = NUCLEUS_LABELS.keys()
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])
389 pie = ROOT.TPie(
"NeutronCaptureNucleus",
"", n_nuclei, vals, colors)
390 for i, pdg
in enumerate(key_order):
391 pie.SetEntryLabel(i, NUCLEUS_LABELS[pdg])
393 pie.SetLabelFormat(
"#splitline{%txt}{%perc}");
397 for plotname, plot
in self.michel_histos.iteritems():
398 if "True" in plotname:
401 ROOT.SetOwnership(plot,
False)
403 legend = ROOT.TLegend(0.55, 0.6, 0.83, 0.9)
404 ROOT.SetOwnership(legend,
False)
408 ROOT.SetOwnership(total_exclusive,
False)
411 total_exclusive.SetLineColor(ROOT.kBlack)
412 total_exclusive.SetLineStyle(3)
414 h_by_pdg = ROOT.THStack()
415 ROOT.SetOwnership(h_by_pdg,
False)
416 hists_to_consider = {
"Total": plot,
"Total-exclusive": total_exclusive}
418 for pdg
in [
"Total",
"Total-exclusive"] + TRUE_ORIGIN_PDGS:
419 h = hists_to_consider[pdg]
421 h_by_pdg.Add(h,
"hist")
423 legend_entry =
"Total" 424 elif pdg ==
"Total-exclusive":
425 legend_entry =
"Total of exclusive categories" 427 legend_entry = LEGEND_ENTRIES[pdg]
428 legend.AddEntry(h, legend_entry,
"l")
430 h_by_pdg.Draw(
"nostack")
431 HistoTools.SanitizeStackLabels(h_by_pdg)
435 if "MichelPID" in plotname:
436 print "fixing Michel PID plot legend" 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)
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
457 legend.AddEntry(h, LEGEND_ENTRIES[pdg],
"f")
459 self.michel_origins.sort(
463 sum(100**(len(TRUE_ORIGIN_PDGS)-i)*(bkdown[pdg]
if pdg
in bkdown
else 0)
for i, pdg
in enumerate(TRUE_ORIGIN_PDGS))
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")
477 HistoTools.SanitizeStackLabels(st)
486 def __init__(self, name, path, mc_norm_factor=AREA_NORM):
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 ),
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()
518 elif hasattr(input,
"GetHists"):
520 assert stack,
"Could not determine the stack to pull the sum out of!" 525 for plot_config
in COMBINATION_PLOTS:
526 for source
in (
"data",
"mc"):
527 collection = getattr(self,
"%s_plots_in" % source)
530 source_infix = source
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)
537 f = ROOT.TFile(plotfile_name)
538 keys = f.GetListOfKeys()
540 raise Exception(
"Don't know which of the %d objects in file you want: %s" % (len(objs),
list(keys)))
542 raise Exception(
"File has nothing in it: '%s'" % plotname_file)
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)
548 collection[plot_config.path] = obj
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)()
556 raise Exception(
"Don't know how to make '%s' plot (no method '%s' defined)" % (plot_config.name, method_name))
560 assert isinstance(mc_canvas, ROOT.TCanvas)
561 mc_stack = [o
for o
in mc_canvas.GetListOfPrimitives()
if isinstance(o, ROOT.THStack)][0]
563 legend = ROOT.TLegend(0.65, 0.7, 0.85, 0.9)
564 ROOT.SetOwnership(legend,
False)
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")
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)
579 h.SetLineColor(ROOT.kRed)
580 legend.AddEntry(h,
"MC",
"l")
582 out_stack.Add(data_plot,
"pe")
586 ROOT.SetOwnership(c,
False)
587 out_stack.Draw(
"nostack")
588 HistoTools.SanitizeStackLabels(out_stack)
595 assert data_plot
and hasattr(mc_plot,
"Scale")
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()
611 assert scale_factor
is not None,
"Could not understand your MC normalization factor: %s" % mc_norm_factor
613 mc_plot.Scale(scale_factor)
624 stack = DataMCPlots.GetMCDataStack(canvas)
625 mc_plot = stack.GetHists()[0]
626 data_plot = stack.GetHists()[1]
628 DataMCPlots.NormalizeMCPlot(data_plot, mc_plot, mc_norm_factor)
632 setattr(self,
"Make_%s_Plot" % plot_name, types.MethodType(f, self))
638 michel_config = [config
for config
in COMBINATION_PLOTS
if config.name ==
"MichelPID"][0]
640 leg = [o
for o
in canvas.GetListOfPrimitives()
if isinstance(o, ROOT.TLegend)][0]
644 stack = DataMCPlots.GetMCDataStack(canvas)
645 mc_plot = stack.GetHists()[0]
646 data_plot = stack.GetHists()[1]
648 DataMCPlots.NormalizeMCPlot(data_plot, mc_plot, michel_config.mc_norm_factor)
650 self.
plots_out[michel_config.name] = canvas
655 for name, plotobj
in self.plots_out.iteritems():
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)) )
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",
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()
674 if sys.argv[1] ==
"comp":
679 NTUPLE_FILELIST = NTUPLE_FILELISTS[sys.argv[1]]
def RegisterBranches(self)
def Make_Default_Plot_Method(self, plot_name, plot_path, mc_norm_factor)
def MakeOverlayMCPlots(self)
def FillMichelEPlots(self)
def __init__(self, testing=False)
def NormalizeMCPlot(data_plot, mc_plot, mc_norm_factor=PlotConfig.AREA_NORM)
def __init__(self, name, path, mc_norm_factor=AREA_NORM)
def Make_MichelPID_Plot(self)
def Make_DataMC_Overlay(self, in_varname)
def PrepareMichelEPlots(self, chain)
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)
def FinalizeMichelEPlots(self)
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
def MakeNeutronCapturePlot(self)
def MakeTrueOriginBreakdownPlot(self)
static void Add(TH3D *h, const int bx, const int by, const int bz, const double w)
def GetMCDataStack(input)
def FillMichelBreakdownPlots(self, histo_stack)
def RegisterMethods(self)