ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbb/python/stack_from_dc.py
Revision: 1.17
Committed: Tue Apr 23 19:16:45 2013 UTC (12 years ago) by nmohr
Content type: text/x-python
Branch: MAIN
CVS Tags: lhcp_UnblindFix
Changes since 1.16: +22 -10 lines
Log Message:
Changes for latest cards

File Contents

# Content
1 #!/usr/bin/env python
2 import pickle
3 import sys, os
4 from optparse import OptionParser
5 from copy import copy,deepcopy
6 from math import sqrt
7 import math
8 from HiggsAnalysis.CombinedLimit.DatacardParser import *
9 from HiggsAnalysis.CombinedLimit.ShapeTools import *
10
11 # import ROOT with a fix to get batch mode (http://root.cern.ch/phpBB3/viewtopic.php?t=3198)
12 hasHelp = False
13 argv = sys.argv
14 for X in ("-h", "-?", "--help"):
15 if X in argv:
16 hasHelp = True
17 argv.remove(X)
18 argv.append( '-b-' )
19 import ROOT
20 from myutils import StackMaker, BetterConfigParser
21
22 ROOT.gROOT.SetBatch(True)
23 ROOT.gSystem.Load("libHiggsAnalysisCombinedLimit.so")
24 argv.remove( '-b-' )
25 if hasHelp: argv.append("-h")
26
27 #CONFIGURE
28 parser = OptionParser()
29 parser.add_option("-D", "--datacard", dest="dc", default="",
30 help="Datacard to be plotted")
31 parser.add_option("-B", "--bin", dest="bin", default="",
32 help="DC bin to plot")
33 parser.add_option("-M", "--mlfit", dest="mlfit", default="",
34 help="mlfit file for nuisances")
35 parser.add_option("-F", "--fitresult", dest="fit", default="s",
36 help="Fit result to be used, 's' (signal+background) or 'b' (background only), default is 's'")
37 parser.add_option("-C", "--config", dest="config", default=[], action="append",
38 help="configuration file")
39 parser.add_option("-V", "--variable", dest="var", default="",
40 help="variable to be fitted")
41
42 (opts, args) = parser.parse_args(argv)
43
44 def readBestFit(theFile):
45 file = ROOT.TFile(theFile)
46 if file == None: raise RuntimeError, "Cannot open file %s" % theFile
47 fit_s = file.Get("fit_s")
48 fit_b = file.Get("fit_b")
49 prefit = file.Get("nuisances_prefit")
50 if fit_s == None or fit_s.ClassName() != "RooFitResult": raise RuntimeError, "File %s does not contain the output of the signal fit 'fit_s'" % args[0]
51 if fit_b == None or fit_b.ClassName() != "RooFitResult": raise RuntimeError, "File %s does not contain the output of the background fit 'fit_b'" % args[0]
52 if prefit == None or prefit.ClassName() != "RooArgSet": raise RuntimeError, "File %s does not contain the prefit nuisances 'nuisances_prefit'" % args[0]
53
54 isFlagged = {}
55 table = {}
56 fpf_b = fit_b.floatParsFinal()
57 fpf_s = fit_s.floatParsFinal()
58 nuiVariation = {}
59 for i in range(fpf_s.getSize()):
60 nuis_s = fpf_s.at(i)
61 name = nuis_s.GetName();
62 nuis_b = fpf_b.find(name)
63 nuis_p = prefit.find(name)
64 if nuis_p != None:
65 mean_p, sigma_p = (nuis_p.getVal(), nuis_p.getError())
66 for fit_name, nuis_x in [('b', nuis_b), ('s',nuis_s)]:
67 if nuis_p != None:
68 valShift = (nuis_x.getVal() - mean_p)/sigma_p
69 #sigShift = nuis_x.getError()/sigma_p
70 print fit_name, name
71 print valShift
72 nuiVariation['%s_%s'%(fit_name,name)] = valShift
73 #print valShift
74 return nuiVariation
75
76 def getBestFitShapes(procs,theShapes,shapeNui,theBestFit,DC,setup,opts,Dict):
77 b = opts.bin
78 for p in procs:
79 counter = 0
80 nom = theShapes[p].Clone()
81 for (lsyst,nofloat,pdf,pdfargs,errline) in DC.systs:
82 if errline[b][p] == 0: continue
83 if ("shape" in pdf):
84 if shapeNui[p+lsyst] > 0.:
85 theVari = 'Up'
86 else:
87 theVari = 'Down'
88 bestNuiVar = theShapes[p+lsyst+theVari].Clone()
89 bestNuiVar.Add(nom,-1.)
90 #print p,lsyst,abs(shapeNui[p+lsyst]),bestNuiVar.Integral()
91 bestNuiVar.Scale(abs(shapeNui[p+lsyst]))
92 if counter == 0:
93 bestNui = bestNuiVar.Clone()
94 else:
95 bestNui.Add(bestNuiVar)
96 counter +=1
97 nom.Add(bestNui)
98 #nom.Scale(theBestFit[p])
99 nom.Scale(theShapes[p].Integral()/nom.Integral()*theBestFit[p])
100 nBins = nom.GetNbinsX()
101 for bin in range(1,nBins+1):
102 nom.SetBinError(bin,theShapes[p].GetBinError(bin))
103 theShapes['%s_%s'%(opts.fit,p)] = nom.Clone()
104 histos = []
105 typs = []
106 sigCount = 0
107 signalList = ['ZH','WH']
108 #signalList = ['VVb']
109 for s in setup:
110 if s in signalList:
111 if sigCount ==0:
112 Overlay=copy(theShapes[Dict[s]])
113 else:
114 Overlay.Add(theShapes[Dict[s]])
115 sigCount += 1
116 else:
117 histos.append(theShapes['%s_%s'%(opts.fit,Dict[s])])
118 typs.append(s)
119 return histos,Overlay,typs
120
121
122 def drawFromDC():
123 config = BetterConfigParser()
124 config.read(opts.config)
125 print opts.config
126 dataname = ''
127 if 'Zmm' in opts.bin: dataname = 'Zmm'
128 elif 'Zee' in opts.bin: dataname = 'Zee'
129 elif 'Wmunu' in opts.bin: dataname = 'Wmn'
130 elif 'Wenu' in opts.bin: dataname = 'Wen'
131 elif 'Znunu' in opts.bin: dataname = 'Znn'
132 elif 'Wtn' in opts.bin: dataname = 'Wtn'
133
134 print 'Variable printing'
135 print opts.var
136 if(opts.var == ''):
137 var = 'BDT'
138 if dataname == 'Zmm' or dataname == 'Zee': var = 'BDT_Zll'
139 elif dataname == 'Wmn' or dataname == 'Wen': var = 'BDT_Wln'
140 elif dataname == 'Znn':
141 if 'HighPt' in opts.bin: var = 'BDT_ZnnHighPt'
142 if 'LowPt' in opts.bin: var = 'BDT_ZnnLowPt'
143 if 'LowCSV' in opts.bin: var = 'BDT_ZnnLowCSV'
144 if dataname == '' or var == 'BDT': raise RuntimeError, "Did not recognise mode or var from %s" % opts.bin
145 else:
146 var = opts.var
147
148 region = 'BDT'
149 ws_var = config.get('plotDef:%s'%var,'relPath')
150 ws_var = ROOT.RooRealVar(ws_var,ws_var,-1.,1.)
151 blind = eval(config.get('Plot:%s'%region,'blind'))
152 Stack=StackMaker(config,var,region,True)
153
154 preFit = False
155 addName = 'PostFit_%s' %(opts.fit)
156 if not opts.mlfit:
157 addName = 'PreFit'
158 preFit = True
159
160 Stack.options['pdfName'] = '%s_%s_%s.pdf' %(var,opts.bin,addName)
161
162 log = eval(config.get('Plot:%s'%region,'log'))
163
164 setup = config.get('Plot_general','setup').split(',')
165 if dataname == 'Zmm' or dataname == 'Zee':
166 try:
167 setup.remove('W1b')
168 setup.remove('W2b')
169 setup.remove('Wlight')
170 setup.remove('WH')
171 except:
172 print '@INFO: Wb / Wligh / WH not present in the datacard'
173 if not dataname == 'Znn':
174 setup.remove('QCD')
175 Stack.setup = setup
176
177 Dict = eval(config.get('LimitGeneral','Dict'))
178 lumi = eval(config.get('Plot_general','lumi'))
179
180 options = copy(opts)
181 options.dataname = "data_obs"
182 options.mass = 0
183 options.format = "%8.3f +/- %6.3f"
184 options.channel = opts.bin
185 options.excludeSyst = []
186 options.norm = False
187 options.stat = False
188 options.bin = True # fake that is a binary output, so that we parse shape lines
189 options.out = "tmp.root"
190 options.fileName = args[0]
191 options.cexpr = False
192 options.fixpars = False
193 options.libs = []
194 options.verbose = 0
195 options.poisson = 0
196 options.nuisancesToExclude = []
197 options.noJMax = None
198 theBinning = ROOT.RooFit.Binning(Stack.nBins,Stack.xMin,Stack.xMax)
199
200 file = open(opts.dc, "r")
201 os.chdir(os.path.dirname(opts.dc))
202 DC = parseCard(file, options)
203 if not DC.hasShapes: DC.hasShapes = True
204 MB = ShapeBuilder(DC, options)
205 theShapes = {}
206 theSyst = {}
207 nuiVar = {}
208 if opts.mlfit:
209 nuiVar = readBestFit(opts.mlfit)
210 if not opts.bin in DC.bins: raise RuntimeError, "Cannot open find %s in bins %s of %s" % (opts.bin,DC.bins,opts.dc)
211 for b in DC.bins:
212 if options.channel != None and (options.channel != b): continue
213 exps = {}
214 expNui = {}
215 shapeNui = {}
216 for (p,e) in DC.exp[b].items(): # so that we get only self.DC.processes contributing to this bin
217 exps[p] = [ e, [] ]
218 expNui[p] = [ e, [] ]
219 for (lsyst,nofloat,pdf,pdfargs,errline) in DC.systs:
220 if pdf in ('param', 'flatParam'): continue
221 # begin skip systematics
222 skipme = False
223 for xs in options.excludeSyst:
224 if re.search(xs, lsyst):
225 skipme = True
226 if skipme: continue
227 # end skip systematics
228 counter = 0
229 for p in DC.exp[b].keys(): # so that we get only self.DC.processes contributing to this bin
230 if errline[b][p] == 0: continue
231 if p == 'QCD' and not dataname == 'Znn': continue
232 if pdf == 'gmN':
233 exps[p][1].append(1/sqrt(pdfargs[0]+1));
234 elif pdf == 'gmM':
235 exps[p][1].append(errline[b][p]);
236 elif type(errline[b][p]) == list:
237 kmax = max(errline[b][p][0], errline[b][p][1], 1.0/errline[b][p][0], 1.0/errline[b][p][1]);
238 exps[p][1].append(kmax-1.);
239 elif pdf == 'lnN':
240 exps[p][1].append(max(errline[b][p], 1.0/errline[b][p])-1.);
241 if not nuiVar.has_key('%s_%s'%(opts.fit,lsyst)):
242 nui = 0.
243 else:
244 nui= nuiVar['%s_%s'%(opts.fit,lsyst)]
245 expNui[p][1].append(abs(1-errline[b][p])*nui);
246 elif ("shape" in pdf):
247 #print 'shape %s %s: %s'%(pdf,p,lsyst)
248 s0 = MB.getShape(b,p)
249 sUp = MB.getShape(b,p,lsyst+"Up")
250 sDown = MB.getShape(b,p,lsyst+"Down")
251 if (s0.InheritsFrom("RooDataHist")):
252 s0 = ROOT.RooAbsData.createHistogram(s0,p,ws_var,theBinning)
253 s0.SetName(p)
254 sUp = ROOT.RooAbsData.createHistogram(sUp,p+lsyst+'Up',ws_var,theBinning)
255 sUp.SetName(p+lsyst+'Up')
256 sDown = ROOT.RooAbsData.createHistogram(sDown,p+lsyst+'Down',ws_var,theBinning)
257 sDown.SetName(p+lsyst+'Down')
258 theShapes[p] = s0.Clone()
259 theShapes[p+lsyst+'Up'] = sUp.Clone()
260 theShapes[p+lsyst+'Down'] = sDown.Clone()
261 if not nuiVar.has_key('%s_%s'%(opts.fit,lsyst)):
262 nui = 0.
263 else:
264 nui= nuiVar['%s_%s'%(opts.fit,lsyst)]
265 shapeNui[p+lsyst] = nui
266 if not 'CMS_vhbb_stat' in lsyst:
267 if counter == 0:
268 theSyst[lsyst] = s0.Clone()
269 theSyst[lsyst+'Up'] = sUp.Clone()
270 theSyst[lsyst+'Down'] = sDown.Clone()
271 else:
272 theSyst[lsyst].Add(s0)
273 theSyst[lsyst+'Up'].Add(sUp.Clone())
274 theSyst[lsyst+'Down'].Add(sDown.Clone())
275 counter += 1
276
277 procs = DC.exp[b].keys(); procs.sort()
278 if not dataname == 'Znn' and 'QCD' in procs:
279 procs.remove('QCD')
280 fmt = ("%%-%ds " % max([len(p) for p in procs]))+" "+options.format;
281 #Compute norm uncertainty and best fit
282 theNormUncert = {}
283 theBestFit = {}
284 for p in procs:
285 relunc = sqrt(sum([x*x for x in exps[p][1]]))
286 print fmt % (p, exps[p][0], exps[p][0]*relunc)
287 theNormUncert[p] = relunc
288 absBestFit = sum([x for x in expNui[p][1]])
289 theBestFit[p] = 1.+absBestFit
290
291 histos = []
292 typs = []
293
294 setup2=copy(setup)
295
296 shapesUp = [[] for _ in range(0,len(setup2))]
297 shapesDown = [[] for _ in range(0,len(setup2))]
298
299 sigCount = 0
300 signalList = ['ZH','WH']
301 #signalList = ['VVb']
302 for p in procs:
303 b = opts.bin
304 for s in setup:
305 if not Dict[s] == p: continue
306 if s in signalList:
307 if sigCount ==0:
308 Overlay=copy(theShapes[Dict[s]])
309 else:
310 Overlay.Add(theShapes[Dict[s]])
311 sigCount += 1
312 else:
313 histos.append(theShapes[Dict[s]])
314 typs.append(s)
315 for (lsyst,nofloat,pdf,pdfargs,errline) in DC.systs:
316 if errline[b][p] == 0: continue
317 if ("shape" in pdf) and not 'CMS_vhbb_stat' in lsyst:
318 print 'syst %s'%lsyst
319 shapesUp[setup2.index(s)].append(theShapes[Dict[s]+lsyst+'Up'])
320 shapesDown[setup2.index(s)].append(theShapes[Dict[s]+lsyst+'Down'])
321
322 #-------------
323 #Compute absolute uncertainty from shapes
324 counter = 0
325 for (lsyst,nofloat,pdf,pdfargs,errline) in DC.systs:
326 sumErr = 0
327 for p in procs:
328 sumErr += errline[b][p]
329 if ("shape" in pdf) and not 'CMS_vhbb_stat' in lsyst and not sumErr == 0:
330 theSystUp = theSyst[lsyst+'Up'].Clone()
331 theSystUp.Add(theSyst[lsyst].Clone(),-1.)
332 theSystUp.Multiply(theSystUp)
333 theSystDown = theSyst[lsyst+'Down'].Clone()
334 theSystDown.Add(theSyst[lsyst].Clone(),-1.)
335 theSystDown.Multiply(theSystDown)
336 if counter == 0:
337 theAbsSystUp = theSystUp.Clone()
338 theAbsSystDown = theSystDown.Clone()
339 else:
340 theAbsSystUp.Add(theSystUp.Clone())
341 theAbsSystDown.Add(theSystDown.Clone())
342 counter +=1
343
344 #-------------
345 #Best fit for shapes
346 if not preFit:
347 histos, Overlay, typs = getBestFitShapes(procs,theShapes,shapeNui,theBestFit,DC,setup,opts,Dict)
348
349 counter = 0
350 errUp=[]
351 total=[]
352 errDown=[]
353 nBins = histos[0].GetNbinsX()
354 print 'total bins %s'%nBins
355 Error = ROOT.TGraphAsymmErrors(histos[0])
356 theTotalMC = histos[0].Clone()
357 for h in range(1,len(histos)):
358 theTotalMC.Add(histos[h])
359
360 total = [[]]*nBins
361 errUp = [[]]*nBins
362 errDown = [[]]*nBins
363 for bin in range(1,nBins+1):
364 binError = theTotalMC.GetBinError(bin)
365 if math.isnan(binError):
366 binError = 0.
367 total[bin-1]=theTotalMC.GetBinContent(bin)
368 #Stat uncertainty of the MC outline
369 errUp[bin-1] = [binError]
370 errDown[bin-1] = [binError]
371 #Relative norm uncertainty of the individual MC
372 for h in range(0,len(histos)):
373 errUp[bin-1].append(histos[h].GetBinContent(bin)*theNormUncert[histos[h].GetName()])
374 errDown[bin-1].append(histos[h].GetBinContent(bin)*theNormUncert[histos[h].GetName()])
375 #Shape uncertainty of the MC
376 for bin in range(1,nBins+1):
377 #print sqrt(theSystUp.GetBinContent(bin))
378 errUp[bin-1].append(sqrt(theAbsSystUp.GetBinContent(bin)))
379 errDown[bin-1].append(sqrt(theAbsSystDown.GetBinContent(bin)))
380
381
382 #Add all in quadrature
383 totErrUp=[sqrt(sum([x**2 for x in bin])) for bin in errUp]
384 totErrDown=[sqrt(sum([x**2 for x in bin])) for bin in errDown]
385
386 #Make TGraph with errors
387 for bin in range(1,nBins+1):
388 if not total[bin-1] == 0:
389 point=histos[0].GetXaxis().GetBinCenter(bin)
390 Error.SetPoint(bin-1,point,1)
391 Error.SetPointEYlow(bin-1,totErrDown[bin-1]/total[bin-1])
392 print 'down %s'%(totErrDown[bin-1]/total[bin-1])
393 Error.SetPointEYhigh(bin-1,totErrUp[bin-1]/total[bin-1])
394 print 'up %s'%(totErrUp[bin-1]/total[bin-1])
395
396 #-----------------------
397 #Read data
398 data0 = MB.getShape(opts.bin,'data_obs')
399 if (data0.InheritsFrom("RooDataHist")):
400 data0 = ROOT.RooAbsData.createHistogram(data0,'data_obs',ws_var,theBinning)
401 data0.SetName('data_obs')
402 datas=[data0]
403 datatyps = [None]
404 datanames=[dataname]
405
406
407 if blind and 'BDT' in var:
408 for bin in range(10,datas[0].GetNbinsX()+1):
409 datas[0].SetBinContent(bin,0)
410
411 histos.append(copy(Overlay))
412 if 'ZH' in signalList and 'WH' in signalList:
413 typs.append('VH')
414 Stack.setup.remove('ZH')
415 if 'WH' in Stack.setup: Stack.setup.remove('WH')
416 Stack.setup.insert(0,'VH')
417 elif 'ZH' in signalList:
418 typs.append('ZH')
419 elif 'WH' in signalList:
420 typs.append('WH')
421 elif 'VVb' in signalList:
422 typs.append('VVb')
423 print Stack.setup
424
425 Stack.histos = histos
426 Stack.typs = typs
427 Stack.datas = datas
428 Stack.datatyps = datatyps
429 Stack.datanames= datanames
430 Stack.overlay = Overlay
431 Stack.AddErrors=Error
432 Stack.lumi = lumi
433 Stack.doPlot()
434
435 print 'i am done!\n'
436 #-------------------------------------------------
437
438
439 if __name__ == "__main__":
440 drawFromDC()
441 sys.exit(0)