ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbb/python/stack_from_dc.py
Revision: 1.7
Committed: Tue Nov 6 17:17:56 2012 UTC (12 years, 6 months ago) by nmohr
Content type: text/x-python
Branch: MAIN
Changes since 1.6: +8 -4 lines
Log Message:
Compatibility with Wln

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 nmohr 1.7 if ("shape" in pdf) and not 'CMS_vhbb_stat' in lsyst:
73 nmohr 1.4 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 nmohr 1.7 if p == 'QCD': continue
183 nmohr 1.1 if pdf == 'gmN':
184     exps[p][1].append(1/sqrt(pdfargs[0]+1));
185     elif pdf == 'gmM':
186     exps[p][1].append(errline[b][p]);
187     elif type(errline[b][p]) == list:
188     kmax = max(errline[b][p][0], errline[b][p][1], 1.0/errline[b][p][0], 1.0/errline[b][p][1]);
189     exps[p][1].append(kmax-1.);
190     elif pdf == 'lnN':
191     exps[p][1].append(max(errline[b][p], 1.0/errline[b][p])-1.);
192     if not nuiVar.has_key('%s_%s'%(opts.fit,lsyst)):
193     nui = 0.
194     else:
195     nui= nuiVar['%s_%s'%(opts.fit,lsyst)]
196     expNui[p][1].append(abs(1-errline[b][p])*nui);
197 nmohr 1.7 elif ("shape" in pdf) and not 'CMS_vhbb_stat' in lsyst:
198     #print 'shape %s %s: %s'%(pdf,p,lsyst)
199 nmohr 1.1 s0 = MB.getShape(b,p)
200     sUp = MB.getShape(b,p,lsyst+"Up")
201     sDown = MB.getShape(b,p,lsyst+"Down")
202     if (s0.InheritsFrom("RooDataHist")):
203     s0 = ROOT.RooAbsData.createHistogram(s0,ws_var)
204     s0.SetName(p)
205     sUp = ROOT.RooAbsData.createHistogram(sUp,ws_var)
206     sUp.SetName(p+lsyst+'Up')
207     sDown = ROOT.RooAbsData.createHistogram(sDown,ws_var)
208     sDown.SetName(p+lsyst+'Down')
209     theShapes[p] = s0.Clone()
210     theShapes[p+lsyst+'Up'] = sUp.Clone()
211     theShapes[p+lsyst+'Down'] = sDown.Clone()
212     if not nuiVar.has_key('%s_%s'%(opts.fit,lsyst)):
213     nui = 0.
214     else:
215     nui= nuiVar['%s_%s'%(opts.fit,lsyst)]
216     shapeNui[p] = nui
217     if counter == 0:
218     theSyst[lsyst] = s0.Clone()
219     theSyst[lsyst+'Up'] = sUp.Clone()
220     theSyst[lsyst+'Down'] = sDown.Clone()
221     else:
222     theSyst[lsyst].Add(s0)
223     theSyst[lsyst+'Up'].Add(sUp.Clone())
224     theSyst[lsyst+'Down'].Add(sDown.Clone())
225     counter += 1
226     procs = DC.exp[b].keys(); procs.sort()
227 nmohr 1.7 if 'QCD' in procs:
228     procs.remove('QCD')
229 nmohr 1.1 fmt = ("%%-%ds " % max([len(p) for p in procs]))+" "+options.format;
230     #Compute norm uncertainty and best fit
231     theNormUncert = {}
232     theBestFit = {}
233     for p in procs:
234     relunc = sqrt(sum([x*x for x in exps[p][1]]))
235     print fmt % (p, exps[p][0], exps[p][0]*relunc)
236     theNormUncert[p] = relunc
237     absBestFit = sum([x for x in expNui[p][1]])
238     theBestFit[p] = 1.+absBestFit
239    
240     histos = []
241     typs = []
242    
243     setup2=copy(setup)
244    
245     shapesUp = [[] for _ in range(0,len(setup2))]
246     shapesDown = [[] for _ in range(0,len(setup2))]
247    
248     for p in procs:
249 nmohr 1.2 b = opts.bin
250 nmohr 1.1 for s in setup:
251     if not Dict[s] == p: continue
252     if 'ZH' == s:
253     Overlay=copy(theShapes[Dict[s]])
254     else:
255     histos.append(theShapes[Dict[s]])
256     typs.append(s)
257     print s
258     for (lsyst,nofloat,pdf,pdfargs,errline) in DC.systs:
259     if errline[b][p] == 0: continue
260 nmohr 1.7 if ("shape" in pdf) and not 'CMS_vhbb_stat' in lsyst:
261 nmohr 1.1 print 'syst %s'%lsyst
262     shapesUp[setup2.index(s)].append(theShapes[Dict[s]+lsyst+'Up'])
263     shapesDown[setup2.index(s)].append(theShapes[Dict[s]+lsyst+'Down'])
264    
265     #-------------
266     #Compute absolute uncertainty from shapes
267     counter = 0
268     for (lsyst,nofloat,pdf,pdfargs,errline) in DC.systs:
269 nmohr 1.7 if ("shape" in pdf) and not 'CMS_vhbb_stat' in lsyst:
270 nmohr 1.1 theSystUp = theSyst[lsyst+'Up'].Clone()
271     theSystUp.Add(theSyst[lsyst].Clone(),-1.)
272     theSystUp.Multiply(theSystUp)
273     theSystDown = theSyst[lsyst+'Down'].Clone()
274     theSystDown.Add(theSyst[lsyst].Clone(),-1.)
275     theSystDown.Multiply(theSystDown)
276     if counter == 0:
277     theAbsSystUp = theSystUp.Clone()
278     theAbsSystDown = theSystDown.Clone()
279     else:
280     theAbsSystUp.Add(theSystUp.Clone())
281     theAbsSystDown.Add(theSystDown.Clone())
282     counter +=1
283    
284     #-------------
285     #Best fit for shapes
286 nmohr 1.4 if not preFit:
287 nmohr 1.6 histos, typs = getBestFitShapes(procs,theShapes,shapeNui,theBestFit,DC,setup,opts,Dict)
288 nmohr 1.1
289     counter = 0
290     errUp=[]
291     total=[]
292     errDown=[]
293     nBins = histos[0].GetNbinsX()
294     print 'total bins %s'%nBins
295     Error = ROOT.TGraphAsymmErrors(histos[0])
296     theTotalMC = histos[0].Clone()
297     for h in range(1,len(histos)):
298     theTotalMC.Add(histos[h])
299    
300     total = [[]]*nBins
301     errUp = [[]]*nBins
302     errDown = [[]]*nBins
303     for bin in range(1,nBins+1):
304     binError = theTotalMC.GetBinError(bin)
305     if math.isnan(binError):
306     binError = 0.
307     total[bin-1]=theTotalMC.GetBinContent(bin)
308     #Stat uncertainty of the MC outline
309     errUp[bin-1] = [binError]
310     errDown[bin-1] = [binError]
311     #Relative norm uncertainty of the individual MC
312     for h in range(0,len(histos)):
313     errUp[bin-1].append(histos[h].GetBinContent(bin)*theNormUncert[histos[h].GetName()])
314     errDown[bin-1].append(histos[h].GetBinContent(bin)*theNormUncert[histos[h].GetName()])
315     #Shape uncertainty of the MC
316     for bin in range(1,nBins+1):
317     #print sqrt(theSystUp.GetBinContent(bin))
318     errUp[bin-1].append(sqrt(theAbsSystUp.GetBinContent(bin)))
319     errDown[bin-1].append(sqrt(theAbsSystDown.GetBinContent(bin)))
320    
321    
322     #Add all in quadrature
323     totErrUp=[sqrt(sum([x**2 for x in bin])) for bin in errUp]
324     totErrDown=[sqrt(sum([x**2 for x in bin])) for bin in errDown]
325    
326     #Make TGraph with errors
327     for bin in range(1,nBins+1):
328     if not total[bin-1] == 0:
329     point=histos[0].GetXaxis().GetBinCenter(bin)
330     Error.SetPoint(bin-1,point,1)
331     Error.SetPointEYlow(bin-1,totErrDown[bin-1]/total[bin-1])
332     print 'down %s'%(totErrDown[bin-1]/total[bin-1])
333     Error.SetPointEYhigh(bin-1,totErrUp[bin-1]/total[bin-1])
334     print 'up %s'%(totErrUp[bin-1]/total[bin-1])
335    
336     #-----------------------
337     #Read data
338 nmohr 1.2 data0 = MB.getShape(opts.bin,'data_obs')
339 nmohr 1.1 if (data0.InheritsFrom("RooDataHist")):
340     data0 = ROOT.RooAbsData.createHistogram(data0,ws_var)
341     data0.SetName('data_obs')
342     datas=[data0]
343     datatyps = [None]
344     datanames=[dataname]
345    
346    
347     if blind:
348     for bin in range(10,datas[0].GetNbinsX()+1):
349     datas[0].SetBinContent(bin,0)
350    
351     histos.append(copy(Overlay))
352     typs.append('ZH')
353    
354     Stack.histos = histos
355     Stack.typs = typs
356     Stack.datas = datas
357     Stack.datatyps = datatyps
358     Stack.datanames= datanames
359     Stack.overlay = Overlay
360     Stack.AddErrors=Error
361     Stack.lumi = lumi
362     Stack.doPlot()
363    
364     print 'i am done!\n'
365     #-------------------------------------------------
366    
367    
368     if __name__ == "__main__":
369     drawFromDC()
370     sys.exit(0)