ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbb/python/tree_stack.py
Revision: 1.27
Committed: Thu Oct 11 13:40:18 2012 UTC (12 years, 7 months ago) by peller
Content type: text/x-python
Branch: MAIN
Changes since 1.26: +0 -1 lines
Log Message:
plotting cache

File Contents

# Content
1 #!/usr/bin/env python
2 from samplesclass import sample
3 from printcolor import printc
4 import pickle
5 import ROOT
6 from ROOT import TFile, TTree
7 import ROOT
8 from array import array
9 from BetterConfigParser import BetterConfigParser
10 import sys, os
11 from mvainfos import mvainfo
12 #from gethistofromtree import getHistoFromTree, orderandadd
13 from Ratio import getRatio
14 from optparse import OptionParser
15 from HistoMaker import HistoMaker, orderandadd
16 import TdrStyles
17 from copy import deepcopy
18
19 #CONFIGURE
20 argv = sys.argv
21 parser = OptionParser()
22 parser.add_option("-P", "--path", dest="path", default="",
23 help="path to samples")
24 parser.add_option("-R", "--reg", dest="region", default="",
25 help="region to plot")
26 parser.add_option("-C", "--config", dest="config", default=[], action="append",
27 help="configuration file")
28 (opts, args) = parser.parse_args(argv)
29 if opts.config =="":
30 opts.config = "config"
31 print opts.config
32 config = BetterConfigParser()
33 config.read(opts.config)
34 anaTag = config.get("Analysis","tag")
35
36
37 path = opts.path
38 region = opts.region
39
40 plotConfig = BetterConfigParser()
41 plotConfig.read('vhbbPlotDef.ini')
42
43 #get locations:
44 Wdir=config.get('Directories','Wdir')
45
46 section='Plot:%s'%region
47
48 Normalize = eval(config.get(section,'Normalize'))
49 log = eval(config.get(section,'log'))
50 blind = eval(config.get(section,'blind'))
51
52 infofile = open(path+'/samples.info','r')
53 info = pickle.load(infofile)
54 infofile.close()
55
56 #options = plot.split(',')
57
58 mass = config.get(section,'Signal')
59
60 vars = (config.get(section, 'vars')).split(',')
61
62 names = [plotConfig.get('plotDef:%s'%x,'relPath') for x in vars]
63 nBins = [eval(plotConfig.get('plotDef:%s'%x,'nBins')) for x in vars]
64 xMin = [eval(plotConfig.get('plotDef:%s'%x,'min')) for x in vars]
65 xMax = [eval(plotConfig.get('plotDef:%s'%x,'max')) for x in vars]
66 xAxis = [plotConfig.get('plotDef:%s'%x,'xAxis') for x in vars]
67
68 for p in range(0,len(names)):
69 if '<mass>' in names[p]:
70 newp= names[p].replace('<mass>',mass)
71 names[p]=newp
72 print names[p]
73
74 if 'ZLight' in region or 'TTbar' in region or 'Zbb' in region: SignalRegion = False
75 else:
76 SignalRegion = True
77 print 'You are in the Signal Region!'
78
79 data = config.get(section,'Datas')
80 if config.has_option(section, 'Datacut'):
81 datacut=config.get(section, 'Datacut')
82 else:
83 datacut = region
84
85 options=[]
86
87 if blind: blindopt='blind'
88 else: blindopt = 'noblind'
89
90 for i in range(0,len(vars)):
91 options.append([names[i],'',xAxis[i],nBins[i],xMin[i],xMax[i],'%s_%s.pdf'%(region,vars[i]),region,datacut,mass,data,blindopt])
92
93
94 setup=config.get('Plot_general','setup')
95 if log:
96 setup=config.get('Plot_general','setupLog')
97 setup=setup.split(',')
98
99 samples=config.get('Plot_general','samples')
100 samples=samples.split(',')
101
102 colorDict=eval(config.get('Plot_general','colorDict'))
103 typLegendDict=eval(config.get('Plot_general','typLegendDict'))
104 #color=color.split(',')
105
106
107 weightF=config.get('Weights','weightF')
108 Group = eval(config.get('Plot_general','Group'))
109
110
111 #GETALL AT ONCE
112
113 Plotter=HistoMaker(path,config,region,options)
114
115 #print '\nProducing Plot of %s\n'%vars[v]
116 Lhistos = [[] for _ in range(0,len(vars))]
117 Ltyps = [[] for _ in range(0,len(vars))]
118 Ldatas = [[] for _ in range(0,len(vars))]
119 Ldatatyps = [[] for _ in range(0,len(vars))]
120 Ldatanames = [[] for _ in range(0,len(vars))]
121
122 def myText(txt="CMS Preliminary",ndcX=0,ndcY=0,size=0.8):
123 ROOT.gPad.Update()
124 text = ROOT.TLatex()
125 text.SetNDC()
126 text.SetTextColor(ROOT.kBlack)
127 text.SetTextSize(text.GetTextSize()*size)
128 text.DrawLatex(ndcX,ndcY,txt)
129 return text
130
131
132 #Find out Lumi:
133 lumicounter=0.
134 lumi=0.
135 for job in info:
136 if job.name in data:
137 lumi+=float(job.lumi)
138 lumicounter+=1.
139
140 if lumicounter > 0:
141 lumi=lumi/lumicounter
142
143 Plotter.lumi=lumi
144
145 for job in info:
146 if eval(job.active):
147 if job.subsamples:
148 for subsample in range(0,len(job.subnames)):
149
150 if job.subnames[subsample] in samples:
151 hTempList, typList = Plotter.getHistoFromTree(job,subsample)
152 for v in range(0,len(vars)):
153 Lhistos[v].append(hTempList[v])
154 Ltyps[v].append(Group[job.subnames[subsample]])
155 print job.subnames[subsample]
156
157 else:
158 if job.name in samples:
159 if job.name == mass:
160 print job.name
161 hTempList, typList = Plotter.getHistoFromTree(job)
162 for v in range(0,len(vars)):
163 if SignalRegion:
164 Lhistos[v].append(hTempList[v])
165 Ltyps[v].append(Group[job.name])
166 Overlaylist= deepcopy(hTempList)
167
168 else:
169 print job.name
170 hTempList, typList = Plotter.getHistoFromTree(job)
171 for v in range(0,len(vars)):
172 Lhistos[v].append(hTempList[v])
173 Ltyps[v].append(Group[job.name])
174
175 elif job.name in data:
176 #print 'DATA'
177 hTemp, typ = Plotter.getHistoFromTree(job)
178 for v in range(0,len(vars)):
179 Ldatas[v].append(hTemp[v])
180 Ldatatyps[v].append(typ[v])
181 Ldatanames[v].append(job.name)
182
183 if not SignalRegion: setup.remove('ZH')
184
185 print setup
186
187 for v in range(0,len(vars)):
188
189 histos = Lhistos[v]
190 typs = Ltyps[v]
191 datas = Ldatas[v]
192 datatyps = Ldatatyps[v]
193 datanames= Ldatanames[v]
194
195 Overlay = Overlaylist[v]
196
197 TdrStyles.tdrStyle()
198
199 c = ROOT.TCanvas(vars[v],'', 600, 600)
200 c.SetFillStyle(4000)
201 c.SetFrameFillStyle(1000)
202 c.SetFrameFillColor(0)
203
204 oben = ROOT.TPad('oben','oben',0,0.3 ,1.0,1.0)
205 oben.SetBottomMargin(0)
206 oben.SetFillStyle(4000)
207 oben.SetFrameFillStyle(1000)
208 oben.SetFrameFillColor(0)
209 unten = ROOT.TPad('unten','unten',0,0.0,1.0,0.3)
210 unten.SetTopMargin(0.)
211 unten.SetBottomMargin(0.35)
212 unten.SetFillStyle(4000)
213 unten.SetFrameFillStyle(1000)
214 unten.SetFrameFillColor(0)
215
216 oben.Draw()
217 unten.Draw()
218
219 oben.cd()
220
221
222
223 allStack = ROOT.THStack(vars[v],'')
224 l = ROOT.TLegend(0.63, 0.60,0.92,0.92)
225 l.SetLineWidth(2)
226 l.SetBorderSize(0)
227 l.SetFillColor(0)
228 l.SetFillStyle(4000)
229 l.SetTextFont(62)
230 l.SetTextSize(0.035)
231 MC_integral=0
232 MC_entries=0
233
234 for histo in histos:
235 MC_integral+=histo.Integral()
236 #MC_entries+=histo.GetEntries()
237 print "\033[1;32m\n\tMC integral = %s\033[1;m"%MC_integral
238
239 #ORDER AND ADD TOGETHER
240 #print typs
241 #print setup
242 histos, typs = orderandadd(histos,typs,setup)
243
244
245 k=len(histos)
246
247 for j in range(0,k):
248 #print histos[j].GetBinContent(1)
249 i=k-j-1
250 histos[i].SetFillColor(int(colorDict[setup[i]]))
251 histos[i].SetLineColor(1)
252 allStack.Add(histos[i])
253
254 d1 = ROOT.TH1F('noData','noData',nBins[v],xMin[v],xMax[v])
255 datatitle='Data'
256 addFlag = ''
257 if 'Zee' in datanames and 'Zmm' in datanames:
258 addFlag = 'Z(l^{-}l^{+})H(b#bar{b})'
259 elif 'Zee' in datanames:
260 addFlag = 'Z(e^{-}e^{+})H(b#bar{b})'
261 elif 'Zmm' in datanames:
262 addFlag = 'Z(#mu^{-}#mu^{+})H(b#bar{b})'
263 for i in range(0,len(datas)):
264 d1.Add(datas[i],1)
265 print "\033[1;32m\n\tDATA integral = %s\033[1;m"%d1.Integral()
266 flow = d1.GetEntries()-d1.Integral()
267 if flow > 0:
268 print "\033[1;31m\tU/O flow: %s\033[1;m"%flow
269
270 Overlay.SetLineColor(2)
271 Overlay.SetLineWidth(2)
272 Overlay.SetFillColor(0)
273 Overlay.SetFillStyle(4000)
274 Overlay.SetNameTitle('Overlay','Overlay')
275
276 l.AddEntry(d1,datatitle,'P')
277 for j in range(0,k):
278 l.AddEntry(histos[j],typLegendDict[typs[j]],'F')
279 l.AddEntry(Overlay,typLegendDict['Overlay'],'L')
280
281 if Normalize:
282 if MC_integral != 0: stackscale=d1.Integral()/MC_integral
283 Overlay.Scale(stackscale)
284 stackhists=allStack.GetHists()
285 for blabla in stackhists:
286 if MC_integral != 0: blabla.Scale(stackscale)
287
288 allMC=allStack.GetStack().Last().Clone()
289
290 allStack.SetTitle()
291 allStack.Draw("hist")
292 allStack.GetXaxis().SetTitle('')
293 yTitle = 'Entries'
294 if not '/' in yTitle:
295 yAppend = '%s' %(allStack.GetXaxis().GetBinWidth(1))
296 yTitle = '%s / %s' %(yTitle, yAppend)
297 allStack.GetYaxis().SetTitle(yTitle)
298 allStack.GetXaxis().SetRangeUser(xMin[v],xMax[v])
299 allStack.GetYaxis().SetRangeUser(0,20000)
300 theErrorGraph = ROOT.TGraphErrors(allMC)
301 theErrorGraph.SetFillColor(ROOT.kGray+3)
302 theErrorGraph.SetFillStyle(3013)
303 theErrorGraph.Draw('SAME2')
304 l.AddEntry(theErrorGraph,"MC uncert. (stat.)","fl")
305 Ymax = max(allStack.GetMaximum(),d1.GetMaximum())*1.7
306 if log:
307 allStack.SetMinimum(0.05)
308 Ymax = Ymax*ROOT.TMath.Power(10,1.6*(ROOT.TMath.Log(1.6*(Ymax/0.1))/ROOT.TMath.Log(10)))*(0.6*0.1)
309 ROOT.gPad.SetLogy()
310 allStack.SetMaximum(Ymax)
311 c.Update()
312 ROOT.gPad.SetTicks(1,1)
313 #allStack.Draw("hist")
314 d1.Draw("E,same")
315 l.SetFillColor(0)
316 l.SetBorderSize(0)
317 l.Draw()
318
319 Overlay.Draw('hist,same')
320
321 tPrel = myText("CMS Preliminary",0.17,0.88,1.04)
322 tLumi = myText("#sqrt{s} = %s, L = %s fb^{-1}"%(anaTag,(float(lumi)/1000.)),0.17,0.83)
323 tAddFlag = myText(addFlag,0.17,0.78)
324
325 unten.cd()
326 ROOT.gPad.SetTicks(1,1)
327
328 ratio, error = getRatio(d1,allMC,xMin[v],xMax[v])
329 ksScore = d1.KolmogorovTest( allMC )
330 chiScore = d1.Chi2Test( allMC , "UWCHI2/NDF")
331 print ksScore
332 print chiScore
333 ratio.SetStats(0)
334 ratio.GetXaxis().SetTitle(xAxis[v])
335 ratioError = ROOT.TGraphErrors(error)
336 ratioError.SetFillColor(ROOT.kGray+3)
337 ratioError.SetFillStyle(3013)
338 ratio.Draw("E1")
339 ratioError.Draw('SAME2')
340 ratio.Draw("E1SAME")
341 ratio.SetTitle("")
342 m_one_line = ROOT.TLine(xMin[v],1,xMax[v],1)
343 m_one_line.SetLineStyle(ROOT.kDashed)
344 m_one_line.Draw("Same")
345
346 if not blind:
347 tKsChi = myText("#chi_{#nu}^{2} = %.3f K_{s} = %.3f"%(chiScore,ksScore),0.17,0.9,1.5)
348 t0 = ROOT.TText()
349 t0.SetTextSize(ROOT.gStyle.GetLabelSize()*2.4)
350 t0.SetTextFont(ROOT.gStyle.GetLabelFont())
351 if not log:
352 t0.DrawTextNDC(0.1059,0.96, "0")
353
354 name = '%s/%s' %(config.get('Directories','plotpath'),options[v][6])
355 c.Print(name)
356
357 print 'i am done!\n'
358
359 sys.exit(0)