ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbb/python/tree_stack.py
Revision: 1.32
Committed: Fri Oct 19 15:44:07 2012 UTC (12 years, 6 months ago) by peller
Content type: text/x-python
Branch: MAIN
Changes since 1.31: +212 -100 lines
Log Message:
BDT plots from TH files

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 array import array
7 from BetterConfigParser import BetterConfigParser
8 import sys, os
9 from mvainfos import mvainfo
10 #from gethistofromtree import getHistoFromTree, orderandadd
11 from optparse import OptionParser
12 from HistoMaker import HistoMaker, orderandadd
13 from copy import copy,deepcopy
14 from StackMaker import StackMaker
15 from math import sqrt
16
17 #CONFIGURE
18 argv = sys.argv
19 parser = OptionParser()
20 parser.add_option("-P", "--path", dest="path", default="",
21 help="path to samples")
22 parser.add_option("-R", "--reg", dest="region", default="",
23 help="region to plot")
24 parser.add_option("-C", "--config", dest="config", default=[], action="append",
25 help="configuration file")
26 (opts, args) = parser.parse_args(argv)
27 if opts.config =="":
28 opts.config = "config"
29 print opts.config
30 config = BetterConfigParser()
31 config.read(opts.config)
32
33 path = opts.path
34 region = opts.region
35
36
37 #get locations:
38 Wdir=config.get('Directories','Wdir')
39 samplesinfo=config.get('Directories','samplesinfo')
40 #limitpath=config.get('Directories','limits')
41 limitpath=path
42 section='Plot:%s'%region
43
44 infofile = open(samplesinfo,'r')
45 info = pickle.load(infofile)
46 infofile.close()
47
48 if 'vhbb_TH_BDT' in region:
49 #-----------Histo from TH File------------------------------------
50 if 'Zee' in region: d='Zee'
51 elif 'Zmm' in region: d='Zmm'
52 if 'LowPt' in region:
53 var='BDT8_RMed'
54 newregion='LowPt_%s'%d
55 elif 'HighPt' in region:
56 var='BDT8_RTight'
57 newregion='HighPt_%s'%d
58 elif 'HighPtLooseBTag' in region:
59 var='BDT8_RTightLooseBTag'
60 newregion='HighPtLooseBTag_%s'%d
61 Stack=StackMaker(config,var,newregion,True)
62
63 log = eval(config.get('Plot:%s'%newregion,'log'))
64
65 if log:
66 setup = config.get('Plot_general','setupLog').split(',')
67 else:
68 setup = config.get('Plot_general','setup').split(',')
69 Dict = eval(config.get('LimitGeneral','Dict'))
70
71 setup.remove('DYc')
72
73 sys_BDT= eval(config.get('LimitGeneral','sys_BDT'))
74 systematicsnaming8TeV = eval(config.get('LimitGeneral','systematicsnaming8TeV'))
75 systs=[systematicsnaming8TeV[s] for s in sys_BDT]
76 if eval(config.get('LimitGeneral','weightF_sys')): systs.append('UEPS')
77 input = ROOT.TFile.Open(limitpath+'/'+region+'.root','read')
78
79 lumi=0
80 for job in info:
81 if job.name == d:
82 lumi=job.lumi
83 break
84 else: pass
85
86 histos = []
87 typs = []
88
89 #systs=[]
90
91 shapesUp = [[] for _ in range(0,len(setup))]
92 shapesDown = [[] for _ in range(0,len(setup))]
93 for s in setup:
94 if 'ZH' in s:
95 Overlay=copy(input.Get(Dict[s]))
96 histos.append(input.Get(Dict[s]))
97 typs.append(s)
98 for syst in systs:
99 shapesUp[setup.index(s)].append(input.Get(Dict[s]+syst+'Up'))
100 shapesDown[setup.index(s)].append(input.Get(Dict[s]+syst+'Down'))
101
102 #print shapesUp
103
104 #calculate the Errors
105 counter = 0
106 total=[]
107 errUp=[]
108 errDown=[]
109 print 'total bins %s'%histos[0].GetNbinsX()
110 for h in range(0,len(histos)):
111 if counter == 0:
112 Error = ROOT.TGraphAsymmErrors(histos[h])
113 for bin in range(1,histos[h].GetNbinsX()+1):
114 if counter == 0 and h == 0:
115 total.append(0)
116 errUp.append(0)
117 errDown.append(0)
118 point=histos[h].GetXaxis().GetBinCenter(bin)
119 total[bin-1]+=histos[h].GetBinContent(bin)
120 for i in range(0,len(shapesUp[h])):
121 errUp[bin-1]+=(shapesUp[h][i].GetBinContent(bin)-histos[h].GetBinContent(bin))**2
122 errDown[bin-1]+=(shapesDown[h][i].GetBinContent(bin)-histos[h].GetBinContent(bin))**2
123 errUp[bin-1]+=(histos[h].GetBinError(bin))**2
124 errDown[bin-1]+=(histos[h].GetBinError(bin))**2
125
126 Error.SetPoint(bin-1,point,1)
127 counter += 1
128
129 for bin in range(0,len(total)):
130 if not total[bin] == 0:
131 Error.SetPointEYlow(bin,sqrt(errDown[bin])/total[bin])
132 Error.SetPointEYhigh(bin,sqrt(errUp[bin])/total[bin])
133
134 #if counter == 0: Errors=Error
135 #else: Errors.Add(Error)
136 #set the errors
137
138 #get the Ratios
139
140
141
142
143 datas=[input.Get('data_obs')]
144 datatyps = [None]
145 datanames=[d]
146
147 #histos.append(copy(Overlay))
148 #typs.append('ZH')
149
150 Stack.histos = histos
151 Stack.typs = typs
152 Stack.datas = datas
153 Stack.datatyps = datatyps
154 Stack.datanames= datanames
155 Stack.overlay = Overlay
156 Stack.AddErrors=Error
157 Stack.lumi = lumi
158 Stack.doPlot()
159
160 print 'i am done!\n'
161 #-------------------------------------------------
162
163
164 else:
165 #----------Histo from trees------------
166 vars = (config.get(section, 'vars')).split(',')
167
168 if 'ZLight' in region or 'TTbar' in region or 'Zbb' in region: SignalRegion = False
169 else:
170 SignalRegion = True
171 print 'You are in the Signal Region!'
172
173 data = config.get(section,'Datas')
174
175 samples=config.get('Plot_general','samples')
176 samples=samples.split(',')
177
178 weightF=config.get('Weights','weightF')
179 Group = eval(config.get('Plot_general','Group'))
180
181 #GETALL AT ONCE
182 options = []
183 Stacks = []
184 for i in range(len(vars)):
185 Stacks.append(StackMaker(config,vars[i],region,SignalRegion))
186 options.append(Stacks[i].options)
187
188 Plotter=HistoMaker(path,config,region,options)
189
190 #print '\nProducing Plot of %s\n'%vars[v]
191 Lhistos = [[] for _ in range(0,len(vars))]
192 Ltyps = [[] for _ in range(0,len(vars))]
193 Ldatas = [[] for _ in range(0,len(vars))]
194 Ldatatyps = [[] for _ in range(0,len(vars))]
195 Ldatanames = [[] for _ in range(0,len(vars))]
196
197 #Find out Lumi:
198 lumicounter=0.
199 lumi=0.
200 for job in info:
201 if job.name in data:
202 lumi+=float(job.lumi)
203 lumicounter+=1.
204
205 if lumicounter > 0:
206 lumi=lumi/lumicounter
207
208 Plotter.lumi=lumi
209 mass = Stacks[0].mass
210
211 for job in info:
212 if eval(job.active):
213 if job.subsamples:
214 for subsample in range(0,len(job.subnames)):
215
216 if job.subnames[subsample] in samples:
217 hTempList, typList = Plotter.getHistoFromTree(job,subsample)
218 for v in range(0,len(vars)):
219 Lhistos[v].append(hTempList[v])
220 Ltyps[v].append(Group[job.subnames[subsample]])
221 print job.subnames[subsample]
222
223 else:
224 if job.name in samples:
225 if job.name == mass:
226 print job.name
227 hTempList, typList = Plotter.getHistoFromTree(job)
228 for v in range(0,len(vars)):
229 if SignalRegion:
230 Lhistos[v].append(hTempList[v])
231 Ltyps[v].append(Group[job.name])
232 Overlaylist= deepcopy(hTempList)
233
234 else:
235 print job.name
236 hTempList, typList = Plotter.getHistoFromTree(job)
237 for v in range(0,len(vars)):
238 Lhistos[v].append(hTempList[v])
239 Ltyps[v].append(Group[job.name])
240
241 elif job.name in data:
242 #print 'DATA'
243 hTemp, typ = Plotter.getHistoFromTree(job)
244 for v in range(0,len(vars)):
245 Ldatas[v].append(hTemp[v])
246 Ldatatyps[v].append(typ[v])
247 Ldatanames[v].append(job.name)
248
249 for v in range(0,len(vars)):
250
251 histos = Lhistos[v]
252 typs = Ltyps[v]
253 Stacks[v].histos = Lhistos[v]
254 Stacks[v].typs = Ltyps[v]
255 Stacks[v].datas = Ldatas[v]
256 Stacks[v].datatyps = Ldatatyps[v]
257 Stacks[v].datanames= Ldatanames[v]
258 Stacks[v].overlay = Overlaylist[v]
259 Stacks[v].lumi = lumi
260 Stacks[v].doPlot()
261 print 'i am done!\n'
262 #----------------------------------------------------
263 sys.exit(0)