TreeAnalyzer/test/plotLimits.py

#!/usr/bin/env python
#____________________________________________________________
#
#
# Francisco Yumiceva
# yumiceva@fnal.gov
#
# Fermilab, 2008
#
#____________________________________________________________         

from ROOT import *
import sys
import os
import sys
from array import array

if os.path.isfile('tdrstyle.C'):
    gROOT.ProcessLine('.L tdrstyle.C')
    ROOT.setTDRStyle()
    print "Found tdrstyle.C file, using this style."
    HasCMSStyle = True
    if os.path.isfile('CMSTopStyle.cc'):
        gROOT.ProcessLine('.L CMSTopStyle.cc+')
        style = CMSTopStyle()
        style.setupICHEPv1()
        print "Found CMSTopStyle.cc file, use TOP style if requested in xml file."


# configuration
sample = sys.argv[1] # mu or e or combined
#sample = "e"
WType = "Gh"
calculatemethod = sys.argv[2] #"bayesian" #"cls" #"bayesian"
ExcludePoints = [] #2300, 2500]

exp_filename = {}
obs_filename = {}
colors = {}

nominal = 'nominal'

exp_filename['nominal'] = "bayesian_limits_expected_fourtop_test.txt"
#exp_filename['nominal'] = "bayesian_limits_expected_stjet_fourtop.txt"
obs_filename['nominal'] = "bayesian_limits_observed_fourtop_test.txt"
#obs_filename['nominal'] = "bayesian_limits_observed_stjet_fourtop.txt"
colors['nominal'] = 1

#exp_filename['only_stat'] = "res_bayesian_nosyst/bayesian_limits_expected_wprh.txt"
#obs_filename['only_stat'] = "res_bayesian_nosyst/bayesian_limits_observed_wprh.txt"
#colors['only_stat'] = 12

#exp_filename['only_rate'] = "res_bayesian_onlyrate/bayesian_limits_expected_wprh.txt"
#obs_filename['only_rate'] = "res_bayesian_onlyrate/bayesian_limits_observed_wprh.txt"
#colors['only_rate'] = ROOT.kRed

#exp_filename['no_btag'] = "res_bayesian_nobtag/bayesian_limits_expected_wprh.txt"
#obs_filename['no_btag'] = "res_bayesian_nobtag/bayesian_limits_observed_wprh.txt"
#colors['no_btag'] = ROOT.kBlue

#exp_filename['no_btag_ttbarmatch'] = "res_bayesian_nobtag_match/bayesian_limits_expected_wprh.txt"
#obs_filename['no_btag_ttbarmatch'] = "res_bayesian_nobtag_match/bayesian_limits_observed_wprh.txt"
#colors['no_btag_ttbarmatch'] = 6

#exp_filename['no_btag_ttbarmatch_wjets'] = "res_bayesian_nobtag_match_wjets/bayesian_limits_expected_wprh.txt"
#obs_filename['no_btag_ttbarmatch_wjets'] = "res_bayesian_nobtag_match_wjets/bayesian_limits_observed_wprh.txt"
#colors['no_btag_ttbarmatch_wjets'] = 7

#exp_filename['no_btag_ttbarmatch_wjets_nlo'] = "res_bayesian_nobtag_match_wjets_nlo/bayesian_limits_expected_wprh.txt"
#obs_filename['no_btag_ttbarmatch_wjets_nlo'] = "res_bayesian_nobtag_match_wjets_nlo/bayesian_limits_observed_wprh.txt"
#colors['no_btag_ttbarmatch_wjets_nlo'] = 9

#exp_filename['WbbWccWqq'] = "res_bayesian/bayesian_limits_expected_wprh.txt"
#obs_filename['WbbWccWqq'] = "res_bayesian/bayesian_limits_observed_wprh.txt"
#colors['WbbWccWqq'] = ROOT.kBlue


if calculatemethod == "cls":
    exp_filename = "cls_limits_expected_wprh.txt"
    obs_filename = "cls_limits_observed_wprh.txt"


# W' cross sections
kFactor = 1.0
xsec = {}
# right-handed
xsec[500] = 0.18182*kFactor;
xsec[700] = 0.012131*kFactor;
xsec[900] = 0.0011484*kFactor;
xsec[1000] = 0.00038029*kFactor;
xsec[1100] = 0.00012830*kFactor;

# read input text files

g_exp = {}
g_sigma1 = {}
g_sigma2 = {}
g_obs = {}


# theory
#xlist = [500, 700 ,900, 1000, 1100]
#ylist = [0.00025442, 0.000099366, 0.000022059, 0.0000098699, 0.0000043429]
xlist = xsec.keys()
xlist.sort()
ylist = []
for ii in xlist:
    ylist.append(xsec[ii])
    
x_theory = array('d')
y_theory = array('d')
x_theory.fromlist( xlist)
y_theory.fromlist( ylist)

g_theory = TGraph( len(x_theory), x_theory, y_theory)
g_theory.SetLineColor(2)
g_theory.SetLineWidth(2)


filekeys = exp_filename.keys()

for test in filekeys:
    
    # read expected
    exp_file = open(exp_filename[test])

    expected = {}
    expected_1sig_up = {}
    expected_1sig_low= {}
    expected_2sig_up = {}
    expected_2sig_low= {}

    for line in exp_file:
        if line.startswith("#"): continue
        list = line.split()
        point = int(list[0])
        skip = False
        for i in ExcludePoints:
            if point == i: skip = True
        if not skip:
            expected[point] = float(list[1])*xsec[point]
            expected_1sig_low[point] = float(list[4])*xsec[point]
            expected_1sig_up[point] = float(list[5])*xsec[point]
            expected_2sig_low[point] = float(list[2])*xsec[point]
            expected_2sig_up[point] = float(list[3])*xsec[point]

    exp_file.close()
#print expected_1sig_low
    print expected

    observed = {}

    obs_file = open(obs_filename[test])
    for line in obs_file:
        if line.startswith("#"): continue
        list = line.split()
        point = int(list[0])
        skip = False
        for i in ExcludePoints:
            if point == i: skip = True
        if not skip:
            observed[point] = float(list[1])*xsec[point]

    obs_file.close()
    print observed

    # create arrays
    x_exp = array('d')
    x_exp_onesigma = array('d')
    y_exp = array('d')
    y_exp_onesigma = array('d')
    y_exp_twosigma = array('d')

    thekeys = expected.keys()
    thekeys.sort()
    tmplist = []
    tmplist1 = []
    tmplist2 = []
    
    for i in thekeys:
        tmplist.append(expected[i])
        tmplist1.append(expected_1sig_low[i])
        tmplist2.append(expected_2sig_low[i])
    thekeys2 = expected.keys() 
    thekeys2.sort()
    thekeys2.reverse()
    for i in thekeys2:
        tmplist1.append(expected_1sig_up[i])
        tmplist2.append(expected_2sig_up[i])


    x_exp_onesigma.fromlist( thekeys + thekeys2 )
    x_exp.fromlist(thekeys)
    y_exp.fromlist( tmplist )
    y_exp_onesigma.fromlist(tmplist1)
    y_exp_twosigma.fromlist(tmplist2)

    print x_exp
    x_obs = array('d')
    y_obs =array('d')

    thekeys = observed.keys()
    thekeys.sort()
    tmplist = []
    for i in thekeys:
        tmplist.append(observed[i])

    x_obs.fromlist( thekeys )
    y_obs.fromlist( tmplist )


    # create graphs
    g_exp[test] = TGraph( len(x_exp), x_exp, y_exp)
    g_exp[test].SetLineStyle(kDashed)
    g_exp[test].SetLineColor( colors[test] )
    g_exp[test].SetLineWidth(2)
  
    g_sigma1[test] = TGraph( len(x_exp_onesigma), x_exp_onesigma, y_exp_onesigma)
    g_sigma1[test].SetFillColor(ROOT.kGreen)
  
    g_sigma2[test] = TGraph( len(x_exp_onesigma), x_exp_onesigma, y_exp_twosigma)
    g_sigma2[test].SetFillColor(ROOT.kYellow)

    g_obs[test] = TGraph( len(x_obs), x_obs, y_obs)
    g_obs[test].SetMarkerStyle(20)
    g_obs[test].SetMarkerColor( colors[test] )
    g_obs[test].SetLineColor( colors[test] )


# create frame
frame = TH2F("frame", "frame", 10, 490, x_obs[len(x_obs)-1]+10, 100, 0, 0.1)
frame.GetXaxis().SetTitle("Gh(t#bar{t}) mass [GeV]")
frame.GetYaxis().SetTitle("#sigma#upointBR(Gh#rightarrow t#bar{t}) [pb]")
frame.SetNdivisions(505)
      
# plot
cv = TCanvas("cv","cv",700,700)
frame.Draw()


g_sigma2[nominal].Draw("F")
g_sigma1[nominal].Draw("F")

# legend
aleg = TLegend(0.20,0.15,0.75,0.45)
aleg.SetMargin(0.12)
aleg.SetTextSize(0.025)
aleg.SetFillColor(10)
aleg.SetBorderSize(0)
aleg.SetLineStyle(0)
aleg.SetFillStyle(0)

    
for test in filekeys:
    
    g_exp[test].Draw("c")
    g_obs[test].Draw("Pl")
    
    aleg.AddEntry( g_obs[test], "95% Observed "+test,"P")
    aleg.AddEntry( g_exp[test], "95% Expected "+test,"L")

aleg.AddEntry( g_sigma1[nominal], "#pm1#sigma Expected","F")
aleg.AddEntry( g_sigma2[nominal], "#pm2#sigma Expected","F")
aleg.AddEntry( g_theory, "Gh","L")
aleg.Draw()

g_theory.Draw("c")
   
sbanner = "#splitline{CMS Preliminary}{4.7 fb^{-1} at #sqrt{s}=7TeV #mu+jets}"
if sample == "combined": sbanner = "#splitline{CMS Preliminary}{4.7 fb^{-1} at #sqrt{s}=7TeV #mu+jets,e+jets}"
if sample == "e": sbanner = "#splitline{CMS Preliminary}{4.7 fb^{-1} at #sqrt{s}=7TeV e+jets}"

banner = TLatex(0.55,0.85, sbanner);
banner.SetNDC();
banner.SetTextSize(0.035);
banner.Draw();

cv.SetLogy();

# Print
outname = "fourtop"+WType+"_" + calculatemethod +"_" + "Limits_"+sample+"_test.png"
#outname = "fourtopStjet"+WType+"_" + calculatemethod +"_" + "Limits_"+sample+".png"
#outname = "fourtop_" + calculatemethod +"_" + "Limits_"+sample+".png"

cv.Print(outname)

raw_input( 'Press ENTER to continue\n ' )
Revision:	1.1
Committed:	Thu Mar 1 17:52:55 2012 UTC (13 years, 2 months ago) by algomez
Content type:	text/x-python
Branch:	MAIN
CVS Tags:	HEAD
Log Message:	produce nice plots for limits
#	Content
1	#!/usr/bin/env python
2	#____________________________________________________________
3	#
4	#
5	# Francisco Yumiceva
6	# yumiceva@fnal.gov
7	#
8	# Fermilab, 2008
9	#
10	#____________________________________________________________
11
12	from ROOT import *
13	import sys
14	import os
15	import sys
16	from array import array
17
18	if os.path.isfile('tdrstyle.C'):
19	gROOT.ProcessLine('.L tdrstyle.C')
20	ROOT.setTDRStyle()
21	print "Found tdrstyle.C file, using this style."
22	HasCMSStyle = True
23	if os.path.isfile('CMSTopStyle.cc'):
24	gROOT.ProcessLine('.L CMSTopStyle.cc+')
25	style = CMSTopStyle()
26	style.setupICHEPv1()
27	print "Found CMSTopStyle.cc file, use TOP style if requested in xml file."
28
29
30	# configuration
31	sample = sys.argv[1] # mu or e or combined
32	#sample = "e"
33	WType = "Gh"
34	calculatemethod = sys.argv[2] #"bayesian" #"cls" #"bayesian"
35	ExcludePoints = [] #2300, 2500]
36
37	exp_filename = {}
38	obs_filename = {}
39	colors = {}
40
41	nominal = 'nominal'
42
43	exp_filename['nominal'] = "bayesian_limits_expected_fourtop_test.txt"
44	#exp_filename['nominal'] = "bayesian_limits_expected_stjet_fourtop.txt"
45	obs_filename['nominal'] = "bayesian_limits_observed_fourtop_test.txt"
46	#obs_filename['nominal'] = "bayesian_limits_observed_stjet_fourtop.txt"
47	colors['nominal'] = 1
48
49	#exp_filename['only_stat'] = "res_bayesian_nosyst/bayesian_limits_expected_wprh.txt"
50	#obs_filename['only_stat'] = "res_bayesian_nosyst/bayesian_limits_observed_wprh.txt"
51	#colors['only_stat'] = 12
52
53	#exp_filename['only_rate'] = "res_bayesian_onlyrate/bayesian_limits_expected_wprh.txt"
54	#obs_filename['only_rate'] = "res_bayesian_onlyrate/bayesian_limits_observed_wprh.txt"
55	#colors['only_rate'] = ROOT.kRed
56
57	#exp_filename['no_btag'] = "res_bayesian_nobtag/bayesian_limits_expected_wprh.txt"
58	#obs_filename['no_btag'] = "res_bayesian_nobtag/bayesian_limits_observed_wprh.txt"
59	#colors['no_btag'] = ROOT.kBlue
60
61	#exp_filename['no_btag_ttbarmatch'] = "res_bayesian_nobtag_match/bayesian_limits_expected_wprh.txt"
62	#obs_filename['no_btag_ttbarmatch'] = "res_bayesian_nobtag_match/bayesian_limits_observed_wprh.txt"
63	#colors['no_btag_ttbarmatch'] = 6
64
65	#exp_filename['no_btag_ttbarmatch_wjets'] = "res_bayesian_nobtag_match_wjets/bayesian_limits_expected_wprh.txt"
66	#obs_filename['no_btag_ttbarmatch_wjets'] = "res_bayesian_nobtag_match_wjets/bayesian_limits_observed_wprh.txt"
67	#colors['no_btag_ttbarmatch_wjets'] = 7
68
69	#exp_filename['no_btag_ttbarmatch_wjets_nlo'] = "res_bayesian_nobtag_match_wjets_nlo/bayesian_limits_expected_wprh.txt"
70	#obs_filename['no_btag_ttbarmatch_wjets_nlo'] = "res_bayesian_nobtag_match_wjets_nlo/bayesian_limits_observed_wprh.txt"
71	#colors['no_btag_ttbarmatch_wjets_nlo'] = 9
72
73	#exp_filename['WbbWccWqq'] = "res_bayesian/bayesian_limits_expected_wprh.txt"
74	#obs_filename['WbbWccWqq'] = "res_bayesian/bayesian_limits_observed_wprh.txt"
75	#colors['WbbWccWqq'] = ROOT.kBlue
76
77
78	if calculatemethod == "cls":
79	exp_filename = "cls_limits_expected_wprh.txt"
80	obs_filename = "cls_limits_observed_wprh.txt"
81
82
83	# W' cross sections
84	kFactor = 1.0
85	xsec = {}
86	# right-handed
87	xsec[500] = 0.18182*kFactor;
88	xsec[700] = 0.012131*kFactor;
89	xsec[900] = 0.0011484*kFactor;
90	xsec[1000] = 0.00038029*kFactor;
91	xsec[1100] = 0.00012830*kFactor;
92
93	# read input text files
94
95	g_exp = {}
96	g_sigma1 = {}
97	g_sigma2 = {}
98	g_obs = {}
99
100
101	# theory
102	#xlist = [500, 700 ,900, 1000, 1100]
103	#ylist = [0.00025442, 0.000099366, 0.000022059, 0.0000098699, 0.0000043429]
104	xlist = xsec.keys()
105	xlist.sort()
106	ylist = []
107	for ii in xlist:
108	ylist.append(xsec[ii])
109
110	x_theory = array('d')
111	y_theory = array('d')
112	x_theory.fromlist( xlist)
113	y_theory.fromlist( ylist)
114
115	g_theory = TGraph( len(x_theory), x_theory, y_theory)
116	g_theory.SetLineColor(2)
117	g_theory.SetLineWidth(2)
118
119
120	filekeys = exp_filename.keys()
121
122	for test in filekeys:
123
124	# read expected
125	exp_file = open(exp_filename[test])
126
127	expected = {}
128	expected_1sig_up = {}
129	expected_1sig_low= {}
130	expected_2sig_up = {}
131	expected_2sig_low= {}
132
133	for line in exp_file:
134	if line.startswith("#"): continue
135	list = line.split()
136	point = int(list[0])
137	skip = False
138	for i in ExcludePoints:
139	if point == i: skip = True
140	if not skip:
141	expected[point] = float(list[1])*xsec[point]
142	expected_1sig_low[point] = float(list[4])*xsec[point]
143	expected_1sig_up[point] = float(list[5])*xsec[point]
144	expected_2sig_low[point] = float(list[2])*xsec[point]
145	expected_2sig_up[point] = float(list[3])*xsec[point]
146
147	exp_file.close()
148	#print expected_1sig_low
149	print expected
150
151	observed = {}
152
153	obs_file = open(obs_filename[test])
154	for line in obs_file:
155	if line.startswith("#"): continue
156	list = line.split()
157	point = int(list[0])
158	skip = False
159	for i in ExcludePoints:
160	if point == i: skip = True
161	if not skip:
162	observed[point] = float(list[1])*xsec[point]
163
164	obs_file.close()
165	print observed
166
167	# create arrays
168	x_exp = array('d')
169	x_exp_onesigma = array('d')
170	y_exp = array('d')
171	y_exp_onesigma = array('d')
172	y_exp_twosigma = array('d')
173
174	thekeys = expected.keys()
175	thekeys.sort()
176	tmplist = []
177	tmplist1 = []
178	tmplist2 = []
179
180	for i in thekeys:
181	tmplist.append(expected[i])
182	tmplist1.append(expected_1sig_low[i])
183	tmplist2.append(expected_2sig_low[i])
184	thekeys2 = expected.keys()
185	thekeys2.sort()
186	thekeys2.reverse()
187	for i in thekeys2:
188	tmplist1.append(expected_1sig_up[i])
189	tmplist2.append(expected_2sig_up[i])
190
191
192	x_exp_onesigma.fromlist( thekeys + thekeys2 )
193	x_exp.fromlist(thekeys)
194	y_exp.fromlist( tmplist )
195	y_exp_onesigma.fromlist(tmplist1)
196	y_exp_twosigma.fromlist(tmplist2)
197
198	print x_exp
199	x_obs = array('d')
200	y_obs =array('d')
201
202	thekeys = observed.keys()
203	thekeys.sort()
204	tmplist = []
205	for i in thekeys:
206	tmplist.append(observed[i])
207
208	x_obs.fromlist( thekeys )
209	y_obs.fromlist( tmplist )
210
211
212	# create graphs
213	g_exp[test] = TGraph( len(x_exp), x_exp, y_exp)
214	g_exp[test].SetLineStyle(kDashed)
215	g_exp[test].SetLineColor( colors[test] )
216	g_exp[test].SetLineWidth(2)
217
218	g_sigma1[test] = TGraph( len(x_exp_onesigma), x_exp_onesigma, y_exp_onesigma)
219	g_sigma1[test].SetFillColor(ROOT.kGreen)
220
221	g_sigma2[test] = TGraph( len(x_exp_onesigma), x_exp_onesigma, y_exp_twosigma)
222	g_sigma2[test].SetFillColor(ROOT.kYellow)
223
224	g_obs[test] = TGraph( len(x_obs), x_obs, y_obs)
225	g_obs[test].SetMarkerStyle(20)
226	g_obs[test].SetMarkerColor( colors[test] )
227	g_obs[test].SetLineColor( colors[test] )
228
229
230	# create frame
231	frame = TH2F("frame", "frame", 10, 490, x_obs[len(x_obs)-1]+10, 100, 0, 0.1)
232	frame.GetXaxis().SetTitle("Gh(t#bar{t}) mass [GeV]")
233	frame.GetYaxis().SetTitle("#sigma#upointBR(Gh#rightarrow t#bar{t}) [pb]")
234	frame.SetNdivisions(505)
235
236	# plot
237	cv = TCanvas("cv","cv",700,700)
238	frame.Draw()
239
240
241	g_sigma2[nominal].Draw("F")
242	g_sigma1[nominal].Draw("F")
243
244	# legend
245	aleg = TLegend(0.20,0.15,0.75,0.45)
246	aleg.SetMargin(0.12)
247	aleg.SetTextSize(0.025)
248	aleg.SetFillColor(10)
249	aleg.SetBorderSize(0)
250	aleg.SetLineStyle(0)
251	aleg.SetFillStyle(0)
252
253
254	for test in filekeys:
255
256	g_exp[test].Draw("c")
257	g_obs[test].Draw("Pl")
258
259	aleg.AddEntry( g_obs[test], "95% Observed "+test,"P")
260	aleg.AddEntry( g_exp[test], "95% Expected "+test,"L")
261
262	aleg.AddEntry( g_sigma1[nominal], "#pm1#sigma Expected","F")
263	aleg.AddEntry( g_sigma2[nominal], "#pm2#sigma Expected","F")
264	aleg.AddEntry( g_theory, "Gh","L")
265	aleg.Draw()
266
267	g_theory.Draw("c")
268
269	sbanner = "#splitline{CMS Preliminary}{4.7 fb^{-1} at #sqrt{s}=7TeV #mu+jets}"
270	if sample == "combined": sbanner = "#splitline{CMS Preliminary}{4.7 fb^{-1} at #sqrt{s}=7TeV #mu+jets,e+jets}"
271	if sample == "e": sbanner = "#splitline{CMS Preliminary}{4.7 fb^{-1} at #sqrt{s}=7TeV e+jets}"
272
273	banner = TLatex(0.55,0.85, sbanner);
274	banner.SetNDC();
275	banner.SetTextSize(0.035);
276	banner.Draw();
277
278	cv.SetLogy();
279
280	# Print
281	outname = "fourtop"+WType+"_" + calculatemethod +"_" + "Limits_"+sample+"_test.png"
282	#outname = "fourtopStjet"+WType+"_" + calculatemethod +"_" + "Limits_"+sample+".png"
283	#outname = "fourtop_" + calculatemethod +"_" + "Limits_"+sample+".png"
284
285	cv.Print(outname)
286
287	raw_input( 'Press ENTER to continue\n ' )