ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbb/python/stack_from_dc.py
Revision: 1.6
Committed: Tue Nov 6 16:06:03 2012 UTC (12 years, 6 months ago) by nmohr
Content type: text/x-python
Branch: MAIN
Changes since 1.5: +8 -4 lines
Log Message:
Add norm from best fit

File Contents

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