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

# User Rev Content
1 nmohr 1.1 #!/usr/bin/env python
2     import pickle
3     import sys, os
4     from optparse import OptionParser
5 peller 1.14 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 bortigno 1.16
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 peller 1.14 from myutils import StackMaker, BetterConfigParser
21    
22     ROOT.gROOT.SetBatch(True)
23     ROOT.gSystem.Load("libHiggsAnalysisCombinedLimit.so")
24 bortigno 1.16 argv.remove( '-b-' )
25     if hasHelp: argv.append("-h")
26 nmohr 1.1
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 bortigno 1.16 parser.add_option("-V", "--variable", dest="var", default="",
40     help="variable to be fitted")
41    
42 nmohr 1.1 (opts, args) = parser.parse_args(argv)
43    
44     def readBestFit(theFile):
45     file = ROOT.TFile(theFile)
46 nmohr 1.4 if file == None: raise RuntimeError, "Cannot open file %s" % theFile
47 nmohr 1.1 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 nmohr 1.6 def getBestFitShapes(procs,theShapes,shapeNui,theBestFit,DC,setup,opts,Dict):
77 nmohr 1.4 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 nmohr 1.11 if ("shape" in pdf):
84     if shapeNui[p+lsyst] > 0.:
85 nmohr 1.4 theVari = 'Up'
86     else:
87     theVari = 'Down'
88     bestNuiVar = theShapes[p+lsyst+theVari].Clone()
89     bestNuiVar.Add(nom,-1.)
90 nmohr 1.11 #print p,lsyst,abs(shapeNui[p+lsyst]),bestNuiVar.Integral()
91     bestNuiVar.Scale(abs(shapeNui[p+lsyst]))
92 nmohr 1.4 if counter == 0:
93     bestNui = bestNuiVar.Clone()
94     else:
95     bestNui.Add(bestNuiVar)
96     counter +=1
97 nmohr 1.10 nom.Add(bestNui)
98 nmohr 1.6 #nom.Scale(theBestFit[p])
99     nom.Scale(theShapes[p].Integral()/nom.Integral()*theBestFit[p])
100 nmohr 1.4 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 nmohr 1.8 sigCount = 0
107 nmohr 1.17 signalList = ['ZH','WH']
108     #signalList = ['VVb']
109 nmohr 1.4 for s in setup:
110 nmohr 1.17 if s in signalList:
111 nmohr 1.8 if sigCount ==0:
112     Overlay=copy(theShapes[Dict[s]])
113     else:
114     Overlay.Add(theShapes[Dict[s]])
115     sigCount += 1
116 nmohr 1.4 else:
117     histos.append(theShapes['%s_%s'%(opts.fit,Dict[s])])
118     typs.append(s)
119 nmohr 1.8 return histos,Overlay,typs
120 nmohr 1.4
121 nmohr 1.1
122     def drawFromDC():
123     config = BetterConfigParser()
124     config.read(opts.config)
125 bortigno 1.16 print opts.config
126 nmohr 1.8 dataname = ''
127     if 'Zmm' in opts.bin: dataname = 'Zmm'
128     elif 'Zee' in opts.bin: dataname = 'Zee'
129 nmohr 1.9 elif 'Wmunu' in opts.bin: dataname = 'Wmn'
130     elif 'Wenu' in opts.bin: dataname = 'Wen'
131 nmohr 1.8 elif 'Znunu' in opts.bin: dataname = 'Znn'
132 nmohr 1.17 elif 'Wtn' in opts.bin: dataname = 'Wtn'
133 nmohr 1.8
134 bortigno 1.16 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 nmohr 1.1 region = 'BDT'
149     ws_var = config.get('plotDef:%s'%var,'relPath')
150 nmohr 1.12 ws_var = ROOT.RooRealVar(ws_var,ws_var,-1.,1.)
151 nmohr 1.1 blind = eval(config.get('Plot:%s'%region,'blind'))
152     Stack=StackMaker(config,var,region,True)
153    
154 nmohr 1.4 preFit = False
155 nmohr 1.6 addName = 'PostFit_%s' %(opts.fit)
156 nmohr 1.4 if not opts.mlfit:
157 nmohr 1.6 addName = 'PreFit'
158 nmohr 1.4 preFit = True
159    
160 nmohr 1.15 Stack.options['pdfName'] = '%s_%s_%s.pdf' %(var,opts.bin,addName)
161 nmohr 1.6
162 nmohr 1.1 log = eval(config.get('Plot:%s'%region,'log'))
163    
164     setup = config.get('Plot_general','setup').split(',')
165 nmohr 1.8 if dataname == 'Zmm' or dataname == 'Zee':
166 bortigno 1.16 try:
167 nmohr 1.17 setup.remove('W1b')
168     setup.remove('W2b')
169 bortigno 1.16 setup.remove('Wlight')
170     setup.remove('WH')
171     except:
172     print '@INFO: Wb / Wligh / WH not present in the datacard'
173 nmohr 1.17 if not dataname == 'Znn':
174     setup.remove('QCD')
175     Stack.setup = setup
176 bortigno 1.16
177 nmohr 1.1 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 nmohr 1.2 options.channel = opts.bin
185 nmohr 1.1 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 nmohr 1.12 theBinning = ROOT.RooFit.Binning(Stack.nBins,Stack.xMin,Stack.xMax)
199 nmohr 1.1
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 nmohr 1.4 nuiVar = {}
208 nmohr 1.1 if opts.mlfit:
209     nuiVar = readBestFit(opts.mlfit)
210 nmohr 1.8 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 nmohr 1.1 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 nmohr 1.17 if p == 'QCD' and not dataname == 'Znn': continue
232 nmohr 1.1 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 nmohr 1.11 elif ("shape" in pdf):
247 nmohr 1.7 #print 'shape %s %s: %s'%(pdf,p,lsyst)
248 nmohr 1.1 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 nmohr 1.12 s0 = ROOT.RooAbsData.createHistogram(s0,p,ws_var,theBinning)
253 nmohr 1.1 s0.SetName(p)
254 nmohr 1.12 sUp = ROOT.RooAbsData.createHistogram(sUp,p+lsyst+'Up',ws_var,theBinning)
255 nmohr 1.1 sUp.SetName(p+lsyst+'Up')
256 nmohr 1.12 sDown = ROOT.RooAbsData.createHistogram(sDown,p+lsyst+'Down',ws_var,theBinning)
257 nmohr 1.1 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 nmohr 1.11 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 nmohr 1.1 procs = DC.exp[b].keys(); procs.sort()
278 nmohr 1.17 if not dataname == 'Znn' and 'QCD' in procs:
279 nmohr 1.7 procs.remove('QCD')
280 nmohr 1.1 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 nmohr 1.8 sigCount = 0
300 nmohr 1.17 signalList = ['ZH','WH']
301     #signalList = ['VVb']
302 nmohr 1.1 for p in procs:
303 nmohr 1.2 b = opts.bin
304 nmohr 1.1 for s in setup:
305     if not Dict[s] == p: continue
306 nmohr 1.17 if s in signalList:
307 nmohr 1.8 if sigCount ==0:
308     Overlay=copy(theShapes[Dict[s]])
309     else:
310     Overlay.Add(theShapes[Dict[s]])
311     sigCount += 1
312 nmohr 1.1 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 nmohr 1.7 if ("shape" in pdf) and not 'CMS_vhbb_stat' in lsyst:
318 nmohr 1.1 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 nmohr 1.11 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 nmohr 1.1 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 nmohr 1.4 if not preFit:
347 nmohr 1.8 histos, Overlay, typs = getBestFitShapes(procs,theShapes,shapeNui,theBestFit,DC,setup,opts,Dict)
348 nmohr 1.1
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 nmohr 1.2 data0 = MB.getShape(opts.bin,'data_obs')
399 nmohr 1.1 if (data0.InheritsFrom("RooDataHist")):
400 nmohr 1.12 data0 = ROOT.RooAbsData.createHistogram(data0,'data_obs',ws_var,theBinning)
401 nmohr 1.1 data0.SetName('data_obs')
402     datas=[data0]
403     datatyps = [None]
404     datanames=[dataname]
405    
406    
407 nmohr 1.17 if blind and 'BDT' in var:
408 nmohr 1.1 for bin in range(10,datas[0].GetNbinsX()+1):
409     datas[0].SetBinContent(bin,0)
410    
411     histos.append(copy(Overlay))
412 nmohr 1.17 if 'ZH' in signalList and 'WH' in signalList:
413 nmohr 1.8 typs.append('VH')
414     Stack.setup.remove('ZH')
415 nmohr 1.17 if 'WH' in Stack.setup: Stack.setup.remove('WH')
416 nmohr 1.8 Stack.setup.insert(0,'VH')
417 nmohr 1.17 elif 'ZH' in signalList:
418 nmohr 1.8 typs.append('ZH')
419 nmohr 1.17 elif 'WH' in signalList:
420 nmohr 1.8 typs.append('WH')
421 nmohr 1.17 elif 'VVb' in signalList:
422     typs.append('VVb')
423     print Stack.setup
424 nmohr 1.1
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)