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
#	User	Rev	Content
1	algomez	1.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 ' )