plotSysts.py
Go to the documentation of this file.
1 import sys, os
2 
3 from ROOT import *
4 from collections import OrderedDict
5 import argparse
6 
7 parser = argparse.ArgumentParser(description='Make some systs plots.')
8 parser.add_argument('inFile', help='The input file to run over',
9  type=str)
10 
11 parser.add_argument('-b', help='Batch mode... no graphics.',
12  action="store_true")
13 
14 args = parser.parse_args()
15 
16 class Decomp:
17  def __init__(self, dir, pot=None):
18  self.hists = OrderedDict()
19  for key in dir.GetListOfKeys():
20  hist = dir.Get(str(key.GetName()) + "/hist")
21  if not hist: # The thing didn't have a key/hist
22  continue
23  hist.SetName(key.GetName())
24  if pot:
25  potHist = dir.Get(str(key.GetName()) + "/pot")
26 
27  if not potHist:
28  raise Exception("Trying to scale to POT, but original POT not found.")
29  potOrig = potHist.GetBinContent(1)
30  if potOrig == 0.0:
31  print "Warning: spectrum", key.GetName(), "has 0 POT, cannot be scaled."
32  continue
33  hist.Scale(pot/potOrig)
34  self.hists[str(key.GetName())] = hist
35 
36  self.sig = self.getSig()
37  self.bkg = self.getBkg()
38  self.components = OrderedDict()
39  self.components["sig"] = self.sig
40  self.components["bkg"] = self.bkg
41 
42  self.errors = OrderedDict()
43  self.errors["sig"] = [0 for i in range(0, self.sig.GetNbinsX() + 2)]
44  self.errors["bkg"] = [0 for i in range(0, self.bkg.GetNbinsX() + 2)]
45 
46  def getSig(self):
47  self.hists["numu_comp"].Add(self.hists["antinumu_comp"])
48  return self.hists["numu_comp"]
49 
50  def getBkg(self):
51  exclude = ["numu_comp", "antinumu_comp", "nc_tot_comp", ]
52  components = [hist for key, hist in self.hists.items() if not key in exclude]
53  before = self.hists["nc_tot_comp"].Integral()
54  for component in components:
55  self.hists["nc_tot_comp"].Add(component)
56 
57  after = self.hists["nc_tot_comp"].Integral()
58  return self.hists["nc_tot_comp"]
59 
60  def addShiftInQuadrature(self, shiftUp, shiftDown):
61  for component, hist in self.components.items():
62  for iBin in range(0, hist.GetNbinsX() + 2):
63  nominal = hist.GetBinContent(iBin)
64  up = abs(nominal - shiftUp.components[component].GetBinContent(iBin))
65  down = abs(nominal - shiftDown.components[component].GetBinContent(iBin))
66  error = max(up, down)
67  self.errors[component][iBin] += error**2
68 
69  def setErrors(self):
70  for component, hist in self.components.items():
71  for iBin in range(0, hist.GetNbinsX() + 2):
72  hist.SetBinError(iBin, sqrt(self.errors[component][iBin]))
73 
74 
75 
76 
77 #fileName = "/local/nova/users/rocco/Decomposition/development/results/2014-10-24_all_plus_qe_separate/twoSample_mc.root"
78 #fileName = "/local/nova/users/rocco/Decomposition/Systs_dev/CAFAna/numu/files_FA14-11-25_mc_exclude_10385.root"
79 fileName = args.inFile
80 
81 f = TFile(fileName, "READ")
82 
83 
84 dirKeys = [k.GetName() for k in f.GetListOfKeys()]
85 
86 unshifted = OrderedDict()
87 
88 shifted = OrderedDict()
89 
90 
91 colors = OrderedDict()
92 colors["-2"] = kRed
93 colors["2"] = kRed
94 colors["-1"] = kBlue -4
95 colors["1"] = kBlue -4
96 
97 redundant = ["ReweightMaCCQEshape", \
98  "ReweightMaCCRESshape", \
99  "ReweightMvCCRESshape", \
100  "ReweightMaNCRESshape", \
101  "ReweightMvNCRESshape", \
102  "ReweightAhtBYshape", \
103  "ReweightBhtBYshape", \
104  "ReweightCV1uBYshape", \
105  "ReweightCV2uBYshape", \
106  "ReweightNormCCQE", \
107  "ReweightNormCCQEenu", \
108  "ReweightNormCCRES", \
109  "ReweightNormNCRES", \
110  "ReweightNC"]
111 
112 broken = ["ReweightMaCOHpi", \
113  "ReweightR0COHpi"]
114 
115 
116 large = ["ReweightMaCCQE", \
117  "ReweightMaCCRES", \
118  "ReweightMvCCRES", \
119  "ReweightMaNCRES", \
120  "ReweightMaNCEL"]
121 
122 # Effectively turn off small for now... not using it.
123 small = []
124 
125 """
126 small = [\
127 "ReweightAGKY_xF1pi", \
128 "ReweightAGKY_pT1pi", \
129 "ReweightBR1gamma", \
130 "ReweightCCQEMomDistroFGtoSF", \
131 "ReweightDISNuclMod", \
132 "ReweightFormZone", \
133 "ReweightEtaNCEL", \
134 "ReweightFrPiProd_pi", \
135 "ReweightNormDISCC", \
136 "ReweightRnubarnuCC", \
137 "ReweightRvbarnCC1pi", \
138 "ReweightRvbarnCC2pi", \
139 "ReweightRvbarnNC1pi", \
140 "ReweightRvbarnNC2pi", \
141 "ReweightRvbarpCC1pi", \
142 "ReweightRvbarpCC2pi", \
143 "ReweightRvbarpNC1pi", \
144 "ReweightRvbarpNC2pi", \
145 "ReweightTheta_Delta2Npi", \
146 "ReweightFrCEx_pi", \
147 "ReweightFrCEx_N", \
148 "ReweightFrAbs_N", \
149 "ReweightFrAbs_pi", \
150 "ReweightFrElas_pi", \
151 "ReweightFrPiProd_pi", \
152 "ReweightFrPiProd_N", \
153 "ReweightFrElas_N", \
154 "ReweightFrInel_N" \
155 ]
156 """
157 
158 
159 
160 def exclude(name, list):
161  exclude = False
162  for item in list:
163  if item == name:
164  exclude = True
165  break
166  return exclude
167 
168 
169 
170 for key in dirKeys:
171  if not "_" in key:
172  unshifted[key] = Decomp(f.Get(key))
173  else:
174  bits = key.split("_")
175  sample = bits[0]
176  syst = "_".join(bits[1:len(bits)-1])
177  shift = bits[len(bits)-1]
178  if not syst in shifted.keys():
179  shifted[syst] = OrderedDict()
180  if not sample in shifted[syst].keys():
181  shifted[syst][sample] = OrderedDict()
182 
183  shifted[syst][sample][shift] = Decomp(f.Get(key))
184 
185 
186 
187 
188 histLists = []
189 ratioLists = []
190 cansToSave = []
191 nSystsSoFar = 0
192 
193 cans = []
194 canNum = 0
195 canPositions = [[10, 10, 600, 450],[615, 10, 600, 450], [10, 480, 600, 450], [615, 480, 600, 450]]
196 
197 def addCan(can):
198  global canNum
199  global cans
200  cans.append(can)
201  canNum +=1
202 
203 def newCan(title):
204  global canNum
205  global cans
206  global canPositions
207 
208  cans.append(TCanvas("can" + str(canNum), title))
209  cans[canNum].SetWindowPosition(canPositions[canNum % 4][0], canPositions[canNum % 4][1])
210  cans[canNum].SetWindowSize(canPositions[canNum % 4][2], canPositions[canNum % 4][3])
211  canNum += 1
212  if canNum % 4 == 0:
213  posUp()
214  return cans[canNum-1]
215 
216 
217 def canUp():
218  for i in range(0, len(cans)):
219  cans[i].RedrawAxis()
220  cans[i].Update()
221 
222 
223 def saveCans(prefix = "", suffix=".pdf"):
224  for i in range(0, len(cans)):
225  cans[i].SaveAs(prefix + cans[i].GetTitle() + suffix)
226 
227 def posUp():
228  global canPositions
229  for i in range(0,len(canPositions)):
230  canPositions[i][0] += 30
231  canPositions[i][1] += 30
232 
233 
234 texts = []
235 for syst, sampleDict in shifted.items():
236  nSystsSoFar += 1
237  print syst
238 # print "N so far:", nSystsSoFar
239  for sample, shiftDict in sampleDict.items():
240 
241  hists = OrderedDict()
242  histLists.append(hists)
243  ratios = OrderedDict( )
244 
245  if not exclude(syst, redundant + broken + large):
246  unshifted[sample].addShiftInQuadrature(shiftDict["1"], shiftDict["-1"])
247 
248  for component in ["sig", "bkg"]:
249  c = newCan("_".join([sample, component,syst]))
250  hists[component] = OrderedDict()
251  ratios[component] = OrderedDict()
252  hists[component]["Nominal"] = unshifted[sample].components[component]
253  hists[component]["Nominal"].GetYaxis().CenterTitle()
254  legEntries = ["Nominal"]
255  ratioHists = []
256  for shift, spectrumDict in shiftDict.items():
257  spectrumDict.components[component].SetLineColor(colors[shift])
258  spectrumDict.components[component].SetLineStyle(2)
259  legEntries += [shift + " #sigma"]
260 
261 
262  hists[component][shift] = spectrumDict.components[component]
263  hists[component][shift].GetYaxis().CenterTitle()
264  ratio = TH1D(hists[component][shift])
265  ratio.Divide(hists[component]["Nominal"])
266  ratio.SetMarkerColor(ratio.GetLineColor())
267  ratio.GetYaxis().SetTitle("#frac{Shifted}{Nominal}")
268  ratio.GetYaxis().CenterTitle()
269  ratio.SetMarkerStyle(0)
270  ratio.SetLineStyle(0)
271  #ratio.SetLineColor(0)
272  for iBin in range(0, ratio.GetNbinsX()+2):
273  ratio.SetBinError(iBin,0.00000001)
274 
275  ratios[component][shift]= ratio
276 
277 
278  drawTwoPad(drawOverlay(hists[component].values(), opt="HIST", preserveStyle=True, legEntries=legEntries, axisMin=0.0001, legPos=[0.65, 0.5],legTextSize=0.055, legHeight=0.35, padFactor=1.7),
279  drawOverlay(ratios[component].values() , opt="E", preserveStyle=True, noLeg=True, axisMin=0.5, axisMax=1.5),
280  can=c, titleOffY=0.05, bottomFraction=0.35 )
281  if exclude(syst, broken):
282  c.cd()
283  t = TText(0.2, 0.8, "Excluded -- Broken")
284  t.SetTextColor(kRed)
285  texts.append(t)
286  t.Draw()
287  if exclude(syst, redundant):
288  c.cd()
289  t = TText(0.2, 0.8, "Excluded -- Redundant")
290  t.SetTextColor(kRed)
291  texts.append(t)
292  t.Draw()
293 
294  if exclude(syst, large):
295  c.cd()
296  t = TText(0.2, 0.8, "Large")
297  t.SetTextColor(kBlue)
298  texts.append(t)
299  t.Draw()
300  if exclude(syst, small):
301  c.cd()
302  t = TText(0.2, 0.8, "Excluded -- Small")
303  t.SetTextColor(kGreen + 2)
304  texts.append(t)
305  t.Draw()
306 
307 
308  cansToSave.append(c)
309 
310 
311 outDir = os.path.join("SystPlots", args.inFile.replace(".root", ""))
312 
313 
314 summedSystHistFile = TFile(
315  os.path.join(outDir,"summed_numu_genie_systs.root"), "RECREATE")
316 
317 
318 
319 errorHists = OrderedDict()
320 
321 
322 
323 for sample, decomp in unshifted.items():
324  decomp.setErrors()
325  errorHists[sample] = OrderedDict()
326  for component, hist in decomp.components.items():
327  c = canMan.newCan("_".join([sample, component, "errors"]))
328  summedSystHistFile.cd()
329  errorHists[sample][component] = TH1D("_".join([sample, component]),
330  ";Error;Energy", hist.GetNbinsX(),
331  hist.GetXaxis().GetXmin(), hist.GetXaxis().GetXmax())
332  for i in range(1, hist.GetNbinsX()+1):
333  weight = 1.0 + (hist.GetBinError(i) / float(hist.GetBinContent(i))) \
334  if hist.GetBinContent(i) > 0 else 1.0
335  errorHists[sample][component].SetBinContent(i, weight)
336 
337  errorHists[sample][component].Write()
338  hist.SetFillColor(kRed)
339  hist.SetFillStyle(3001)
340  #hist.Draw("")
341  hist.Draw("")
342  cansToSave.append(c)
343 
344 
345 summedSystHistFile.Close()
346 
347 
348 
349 canMan.canUp()
350 
351 
352 if not os.path.isdir(outDir): os.mkdir(outDir)
353 
354 for can in cansToSave:
355  can.SaveAs(os.path.join(outDir ,can.GetTitle() + ".pdf"))
356 
357 
358 
359 
def posUp()
Definition: plotSysts.py:227
keys
Reco plots.
Definition: caf_analysis.py:46
c1 Update()
correl_xv GetYaxis() -> SetDecimals()
def setErrors(self)
Definition: plotSysts.py:69
hmean SetLineStyle(2)
T sqrt(T number)
Definition: d0nt_math.hpp:156
def addShiftInQuadrature(self, shiftUp, shiftDown)
Definition: plotSysts.py:60
def newCan(title)
Definition: plotSysts.py:203
double Integral(const Spectrum &s, double pot)
gargamelle SetTitle("Gargamelle #nu_{e} CC data")
float abs(float number)
Definition: d0nt_math.hpp:39
def addCan(can)
Definition: plotSysts.py:197
def canUp()
Definition: canMan.py:39
def getSig(self)
Definition: plotSysts.py:46
hmean SetLineColor(4)
def newCan(title)
Definition: canMan.py:25
TGraph * join(TGraph *a, TGraph *b, int col)
Definition: contours.C:359
def saveCans(prefix="", suffix=".pdf")
Definition: plotSysts.py:223
def drawOverlay(listOfHists, name="", opt="", varOpts=[""], colors=[kRed, kBlue, kGreen, kMagenta, kOrange, kBlack, kViolet, lineStyles=[1], legEntries=[], logY=False, axisMin=0, axisMax=None, legPos="right", legTextSize=0, inputCan=None, legHeight=0.15, legWidth=0.2, supressZero=False, padFactor=1.1, lineWidths=[gStyle.GetHistLineWidth()], fillStyles=[0], preserveStyle=False, noLeg=False)
Definition: canMan.py:69
def canUp()
Definition: plotSysts.py:217
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]))
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
def __init__(self, dir, pot=None)
Definition: plotSysts.py:17
def exclude(name, list)
Definition: plotSysts.py:160
cosmicTree SaveAs("cosmicTree.root")
TCanvas * drawTwoPad(TCanvas *top, TCanvas *bottom, TCanvas *can=NULL, float titleOffY=0.035, float leftMargin=0.13, float topMargin=0.1, float bottomFraction=0.4)
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
def getBkg(self)
Definition: plotSysts.py:50
gm Write()