ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbb/python/scaleFactorPlot.py
Revision: 1.2
Committed: Sun Mar 17 13:29:34 2013 UTC (12 years, 2 months ago) by bortigno
Content type: text/x-python
Branch: MAIN
Changes since 1.1: +65 -23 lines
Log Message:
@UPDATE

File Contents

# User Rev Content
1 bortigno 1.1 #!/usr/bin/env python
2     import re
3     from sys import argv, stdout, stderr, exit
4     from optparse import OptionParser
5    
6     # import ROOT with a fix to get batch mode (http://root.cern.ch/phpBB3/viewtopic.php?t=3198)
7     hasHelp = False
8     for X in ("-h", "-?", "--help"):
9     if X in argv:
10     hasHelp = True
11     argv.remove(X)
12     argv.append( '-b-' )
13     import ROOT
14     ROOT.gROOT.SetBatch(True)
15     ROOT.gSystem.Load("libHiggsAnalysisCombinedLimit.so")
16     argv.remove( '-b-' )
17     if hasHelp: argv.append("-h")
18    
19     parser = OptionParser(usage="usage: %prog [options] in.root \nrun with --help to get list of options")
20     parser.add_option("--vtol", "--val-tolerance", dest="vtol", default=0.30, type="float", help="Report nuisances whose value changes by more than this amount of sigmas")
21     parser.add_option("--stol", "--sig-tolerance", dest="stol", default=0.10, type="float", help="Report nuisances whose sigma changes by more than this amount")
22     parser.add_option("--vtol2", "--val-tolerance2", dest="vtol2", default=2.0, type="float", help="Report severely nuisances whose value changes by more than this amount of sigmas")
23     parser.add_option("--stol2", "--sig-tolerance2", dest="stol2", default=0.50, type="float", help="Report severely nuisances whose sigma changes by more than this amount")
24     parser.add_option("-a", "--all", dest="all", default=False, action="store_true", help="Print all nuisances, even the ones which are unchanged w.r.t. pre-fit values.")
25     parser.add_option("-A", "--absolute", dest="abs", default=False, action="store_true", help="Report also absolute values of nuisance values and errors, not only the ones normalized to the input sigma")
26     parser.add_option("-p", "--poi", dest="poi", default="r", type="string", help="Name of signal strength parameter (default is 'r' as per text2workspace.py)")
27     parser.add_option("-f", "--format", dest="format", default="text", type="string", help="Output format ('text', 'latex', 'twiki'")
28     parser.add_option("-g", "--histogram", dest="plotfile", default=None, type="string", help="If true, plot the pulls of the nuisances to the given file.")
29     parser.add_option("--sf", "--scale-factor", dest="plotsf", default=None, type="string", help="If true, plot the scale factors in a histograms.")
30    
31     (options, args) = parser.parse_args()
32     if len(args) == 0:
33     parser.print_usage()
34     exit(1)
35    
36     file = ROOT.TFile(args[0])
37     if file == None: raise RuntimeError, "Cannot open file %s" % args[0]
38     fit_s = file.Get("fit_s")
39     fit_b = file.Get("fit_b")
40     prefit = file.Get("nuisances_prefit")
41     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]
42     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]
43     if prefit == None or prefit.ClassName() != "RooArgSet": raise RuntimeError, "File %s does not contain the prefit nuisances 'nuisances_prefit'" % args[0]
44    
45     isFlagged = {}
46     table = {}
47     fpf_b = fit_b.floatParsFinal()
48     fpf_s = fit_s.floatParsFinal()
49     pulls = []
50     for i in range(fpf_s.getSize()):
51     nuis_s = fpf_s.at(i)
52     name = nuis_s.GetName();
53     nuis_b = fpf_b.find(name)
54     nuis_p = prefit.find(name)
55     row = []
56     flag = False;
57     mean_p, sigma_p = 0,0
58     if nuis_p == None:
59     if not options.abs: continue
60     row += [ "[%.2f, %.2f]" % (nuis_s.getMin(), nuis_s.getMax()) ]
61     else:
62     mean_p, sigma_p = (nuis_p.getVal(), nuis_p.getError())
63     if options.abs: row += [ "%.2f +/- %.2f" % (nuis_p.getVal(), nuis_p.getError()) ]
64     for fit_name, nuis_x in [('b', nuis_b), ('s',nuis_s)]:
65     if nuis_x == None:
66     row += [ " n/a " ]
67     else:
68     row += [ "%+.2f +/- %.2f" % (nuis_x.getVal(), nuis_x.getError()) ]
69     if nuis_p != None:
70     valShift = (nuis_x.getVal() - mean_p)/sigma_p
71     if fit_name == 'b':
72     pulls.append(valShift)
73     sigShift = nuis_x.getError()/sigma_p
74     if options.abs:
75     row[-1] += " (%+4.2fsig, %4.2f)" % (valShift, sigShift)
76     else:
77     row[-1] = " %+4.2f, %4.2f" % (valShift, sigShift)
78     if (abs(valShift) > options.vtol2 or abs(sigShift-1) > options.stol2):
79     isFlagged[(name,fit_name)] = 2
80     flag = True
81     elif (abs(valShift) > options.vtol or abs(sigShift-1) > options.stol):
82     if options.all: isFlagged[(name,fit_name)] = 1
83     flag = True
84     elif options.all:
85     flag = True
86     row += [ "%+4.2f" % fit_s.correlation(name, options.poi) ]
87     if flag or options.all: table[name] = row
88    
89     fmtstring = "%-40s %15s %15s %10s"
90     highlight = "*%s*"
91     morelight = "!%s!"
92     pmsub, sigsub = None, None
93     if options.format == 'text':
94     if options.abs:
95     fmtstring = "%-40s %15s %30s %30s %10s"
96     print fmtstring % ('name', 'pre fit', 'b-only fit', 's+b fit', 'rho')
97     else:
98     print fmtstring % ('name', 'b-only fit', 's+b fit', 'rho')
99     elif options.format == 'latex':
100     pmsub = (r"(\S+) \+/- (\S+)", r"$\1 \\pm \2$")
101     sigsub = ("sig", r"$\\sigma$")
102     highlight = "\\textbf{%s}"
103     morelight = "{{\\color{red}\\textbf{%s}}}"
104     if options.abs:
105     fmtstring = "%-40s & %15s & %30s & %30s & %6s \\\\"
106     print "\\begin{tabular}{|l|r|r|r|r|} \\hline ";
107     print (fmtstring % ('name', 'pre fit', '$b$-only fit', '$s+b$ fit', r'$\rho(\theta, \mu)$')), " \\hline"
108     else:
109     fmtstring = "%-40s & %15s & %15s & %6s \\\\"
110     print "\\begin{tabular}{|l|r|r|r|} \\hline ";
111     #what = r"$(x_\text{out} - x_\text{in})/\sigma_{\text{in}}$, $\sigma_{\text{out}}/\sigma_{\text{in}}$"
112     what = r"\Delta x/\sigma_{\text{in}}$, $\sigma_{\text{out}}/\sigma_{\text{in}}$"
113     print fmtstring % ('', '$b$-only fit', '$s+b$ fit', '')
114     print (fmtstring % ('name', what, what, r'$\rho(\theta, \mu)$')), " \\hline"
115     elif options.format == 'twiki':
116     pmsub = (r"(\S+) \+/- (\S+)", r"\1 ± \2")
117     sigsub = ("sig", r"σ")
118     highlight = "<b>%s</b>"
119     morelight = "<b style='color:red;'>%s</b>"
120     if options.abs:
121     fmtstring = "| <verbatim>%-40s</verbatim> | %-15s | %-30s | %-30s | %-15s |"
122     print "| *name* | *pre fit* | *b-only fit* | *s+b fit* | "
123     else:
124     fmtstring = "| <verbatim>%-40s</verbatim> | %-15s | %-15s | %-15s |"
125     print "| *name* | *b-only fit* | *s+b fit* | *corr.* |"
126     elif options.format == 'html':
127     pmsub = (r"(\S+) \+/- (\S+)", r"\1 &plusmn; \2")
128     sigsub = ("sig", r"&sigma;")
129     highlight = "<b>%s</b>"
130     morelight = "<strong>%s</strong>"
131     print """
132     <html><head><title>Comparison of nuisances</title>
133     <style type="text/css">
134     td, th { border-bottom: 1px solid black; padding: 1px 1em; }
135     td { font-family: 'Consolas', 'Courier New', courier, monospace; }
136     strong { color: red; font-weight: bolder; }
137     </style>
138     </head><body style="font-family: 'Verdana', sans-serif; font-size: 10pt;"><h1>Comparison of nuisances</h1>
139     <table>
140     """
141     if options.abs:
142     print "<tr><th>nuisance</th><th>pre fit</th><th>background fit </th><th>signal fit</th><th>correlation</th></tr>"
143     fmtstring = "<tr><td><tt>%-40s</tt> </td><td> %-15s </td><td> %-30s </td><td> %-30s </td><td> %-15s </td></tr>"
144     else:
145     what = "&Delta;x/&sigma;<sub>in</sub>, &sigma;<sub>out</sub>/&sigma;<sub>in</sub>";
146     print "<tr><th>nuisance</th><th>background fit<br/>%s </th><th>signal fit<br/>%s</th><th>&rho;(&mu;, &theta;)</tr>" % (what,what)
147     fmtstring = "<tr><td><tt>%-40s</tt> </td><td> %-15s </td><td> %-15s </td><td> %-15s </td></tr>"
148    
149     names = table.keys()
150     names.sort()
151     highlighters = { 1:highlight, 2:morelight };
152     for n in names:
153     v = table[n]
154     if options.format == "latex": n = n.replace(r"_", r"\_")
155     if pmsub != None: v = [ re.sub(pmsub[0], pmsub[1], i) for i in v ]
156     if sigsub != None: v = [ re.sub(sigsub[0], sigsub[1], i) for i in v ]
157     if (n,'b') in isFlagged: v[-3] = highlighters[isFlagged[(n,'b')]] % v[-3]
158     if (n,'s') in isFlagged: v[-2] = highlighters[isFlagged[(n,'s')]] % v[-2]
159     if options.abs:
160     print fmtstring % (n, v[0], v[1], v[2], v[3])
161     else:
162     print fmtstring % (n, v[0], v[1], v[2])
163    
164     if options.format == "latex":
165     print " \\hline\n\end{tabular}"
166     elif options.format == "html":
167     print "</table></body></html>"
168    
169     if options.plotfile:
170     import ROOT
171     ROOT.gROOT.SetStyle("Plain")
172     ROOT.gStyle.SetOptFit(1)
173     histogram = ROOT.TH1F("pulls", "Pulls", 60, -3, 3)
174     for pull in pulls:
175     histogram.Fill(pull)
176     canvas = ROOT.TCanvas("asdf", "asdf", 800, 800)
177     histogram.GetXaxis().SetTitle("pull")
178     histogram.SetTitle("Post-fit nuisance pull distribution")
179     histogram.SetMarkerStyle(20)
180     histogram.SetMarkerSize(2)
181     #histogram.Fit("gaus")
182     histogram.Draw("pe")
183     canvas.SaveAs(options.plotfile)
184    
185    
186 bortigno 1.2 ########### SCALE FACTOR PLOT ###############
187     ## If you fit scale factors this part of ##
188     ## will create a plt reading from the ##
189     ## nuisances. All the nuisance cointaning ##
190     ## "SF" in their name will be plotted ##
191     #############################################
192    
193 bortigno 1.1 import numpy
194     if options.plotsf:
195     n=0
196     labels=[]
197 bortigno 1.2 v_s=[] # dX/sigma_in values
198     v_b=[] # sigma_out/sigma_in values
199     rho=[] # rho
200     sys=[] # systematics. Not properly filled yet
201     x_position=[] # x position in the plot, basically the SF value
202     y_position=[] # y position in the plot. @TO IMPLEMENT: According to the number of entries per background the spacing is different.
203     # loop on all the nuisances
204 bortigno 1.1 for name in names:
205 bortigno 1.2 # take only the nuisances which have SF in their name
206 bortigno 1.1 if ('SF' in name ):
207     print name
208     labels.append(name)
209     n+=1
210 bortigno 1.2 # Take the values
211     # The usual order is: dX/sigma_in for background only - sigma_out/sigma_in for background only - dX/sigma_in for background+signal - sigma_out/sigma_in for background+signal - rho
212     v = table[name]
213     # forget about the flag
214 bortigno 1.1 v = [ re.sub('!','',i) for i in v ]
215     if pmsub != None: v = [ re.sub(pmsub[0], pmsub[1], i) for i in v ]
216     if sigsub != None: v = [ re.sub(sigsub[0], sigsub[1], i) for i in v ]
217     v_b.append(v[0].split(','))
218     v_s.append(v[1].split(','))
219     rho.append(v[2].split(','))
220     x_position.append(0.)
221 bortigno 1.2 y_position.append(n-1.)
222     sys.append(0.1) # need to decide how to quote the systematics and where to take them from
223    
224     # count the number of different channels
225     channels = [0.,0.,0.]
226     print labels
227     for label in labels:
228     if label.find('Zll') > 0. : channels[0] = 1.
229     if label.find('Wln') > 0. : channels[1] = 1.
230     if label.find('Znn') > 0. : channels[2] = 1.
231     # calculate the shift on the y position
232     shift=0.
233     for channel in channels: shift+=channel;
234     print shift
235     shift=1./(shift+1)
236     print shift
237     #shift the elements in the array
238     for i in range(0,len(y_position)): y_position[i]+=shift
239     print y_position
240     # clean the label
241     labels = [re.sub('CMS_vhbb_','',label) for label in labels ]
242     try:
243     labels = [re.sub('_Zll_SF_8TeV','',label) for label in labels ]
244     labels = [re.sub('_Wln_SF_8TeV','',label) for label in labels ]
245     labels = [re.sub('_Znn_SF_8TeV','',label) for label in labels ]
246     except:
247     print '@WARNING: No usual naming for datacard scale factors nuisances'
248     print labels
249    
250     #CMS_vhbb_TT_Zll_SF_8TeV ! +0.56, 0.13! ! +0.54, 0.13! -0.00
251     #CMS_vhbb_ZjHF_Zll_SF_8TeV +1.79, 1.14 +1.80, 1.13 +0.00
252     #CMS_vhbb_ZjLF_Zll_SF_8TeV ! +3.51, 0.94! ! +3.42, 0.95! -0.00
253    
254 bortigno 1.1 print n
255     print labels
256     print v_s
257     print v_b
258     print rho
259     sf=[]
260     sf_e=[]
261 bortigno 1.2 initial_uncertainty=0.2 # initial uncertainty. @TO FIX: this can go in a config or as input arg
262 bortigno 1.1 for i in v_b:
263     try:
264 bortigno 1.2 sf.append(1+initial_uncertainty*eval(i[0])) # calculate the actual value for the scale factors
265     sf_e.append(initial_uncertainty*eval(i[1])) # calculate the actula value for the uncertainties
266 bortigno 1.1 except:
267     print 'Problem evaluating the SF or the SF errors'
268     print sf
269     print sf_e
270     c = ROOT.TCanvas("c","c",600,600)
271    
272 bortigno 1.2
273     label_dictionary = {"TT":"t#bar{t}","ZjHF":"Z+b#bar{b}","ZjLF":"Z+udscg"}
274 bortigno 1.1 #labels = ["W+udscg","W+b#bar{b}","Z+udscg","Z+b#bar{b}","t#bar{t}"]
275     #d = numpy.array([0.94, 1.72, 1.10,1.08,1.01])
276     d = numpy.array(sf)
277     #e = numpy.array([0.02,0.16,0.02,0.04,0.02])
278     e = numpy.array(sf_e)
279     #for the moment just random systematics
280     sys_e = numpy.array(sys)
281     #TO FIX: Y position of the point.
282     p = numpy.array(y_position)
283     zero = numpy.array(x_position)
284    
285 bortigno 1.2 h2 = ROOT.TH2F("h2","",1,0.7,2.2,n,0,n)
286 bortigno 1.1 h2.GetXaxis().SetTitle("Scale factor")
287    
288     for i in range(n):
289 bortigno 1.2 h2.GetYaxis().SetBinLabel(i+1,label_dictionary[labels[i]])
290 bortigno 1.1
291     g = ROOT.TGraphErrors(n,d,p,e,zero)
292     g.SetFillColor(0)
293     g.SetLineColor(2)
294     g.SetLineWidth(3)
295     g.SetMarkerStyle(21)
296    
297     g2 = ROOT.TGraphErrors(n,d,p,sys_e,zero)
298     g2.SetFillColor(0)
299     #g2.SetLineColor(2);
300     g2.SetLineWidth(3);
301     g2.SetMarkerStyle(21);
302    
303     h2.Draw(); ROOT.gStyle.SetOptStat(0);
304     h2.GetXaxis().SetTitleSize(0.04);
305     h2.GetXaxis().SetLabelSize(0.04);
306     h2.GetYaxis().SetLabelSize(0.06);
307     #h2.SetFillStyle(4000)
308     c.SetFillStyle(4000)
309    
310     globalFitBand = ROOT.TBox(1.0, 0., 1.5, n);
311     globalFitBand.SetFillStyle(3013);
312     globalFitBand.SetFillColor(65);
313     globalFitBand.SetLineStyle(0);
314     #globalFitBand.Draw("same");
315 bortigno 1.2 globalFitLine = ROOT.TLine(1.0, 0., 1.0, n);
316 bortigno 1.1 globalFitLine.SetLineWidth(2);
317     globalFitLine.SetLineColor(214);#214
318     globalFitLine.Draw("same");
319    
320     #l2 = ROOT.TLegend(0.5, 0.82,0.92,0.95)
321 bortigno 1.2 l2 = ROOT.TLegend(0.55, 0.85,0.85,0.87)
322 bortigno 1.1 l2.SetLineWidth(2)
323     l2.SetBorderSize(0)
324     l2.SetFillColor(0)
325     l2.SetFillStyle(4000)
326     l2.SetTextFont(62)
327 bortigno 1.2 l2.AddEntry(g,"SF","p")
328     l2.AddEntry(g,"Stat.","l")
329     l2.AddEntry(g2,"Syst.","l")
330 bortigno 1.1 l2.SetTextSize(0.035)
331 bortigno 1.2 l2.SetNColumns(3)
332 bortigno 1.1 l2.Draw("same")
333    
334     g.Draw("P same")
335     g2.Draw("[] same")
336     ROOT.gPad.SetLeftMargin(0.2)
337     ROOT.gPad.Update()
338     c.Print("histo.pdf")
339