VHbb/python/tree_stack.py

#!/usr/bin/env python
from samplesclass import sample
from printcolor import printc
import pickle
import ROOT 
from ROOT import TFile, TTree
import ROOT
from array import array
from BetterConfigParser import BetterConfigParser
import sys
from mvainfos import mvainfo
from gethistofromtree import getHistoFromTree, orderandadd
from Ratio import getRatio

#warnings.filterwarnings( action='ignore', category=RuntimeWarning, message='creating converter.*' )


#load config
config = BetterConfigParser()
config.read('./config7TeV_ZZ')

#get locations:
Wdir=config.get('Directories','Wdir')


Normalize=False

path = sys.argv[1]
var = sys.argv[2]

if 'bb' in var or 'Light' in var or 'Top' in var:
    Normalize=True

Normalize=False


plot=config.get('Plot',var)

infofile = open(path+'/samples.info','r')
info = pickle.load(infofile)
infofile.close()

options = plot.split(',')
name=options[1]
title = options[2]
nBins=int(options[3])
xMin=float(options[4])
xMax=float(options[5])

mass=options[9]
data=options[10]

blind=options[11]

setup=config.get('Plot','setup')
setup=setup.split(',')

color=config.get('Plot','color')
color=color.split(',')


weightF=config.get('Weights','weightF')
Group = eval(config.get('LimitGeneral','Group'))


print '\nProducing Plot of %s\n'%title


histos = []
typs = []
datas = []
datatyps =[]
datanames=[]
'''
for job in info:
    if job.type == 'BKG':
        #print 'MC'
        hTemp, typ = getHistoFromTree(job,options,1)
        histos.append(hTemp)
        typs.append(typ)
    elif job.type == 'SIG' and job.name == mass:
        hTemp, typ = getHistoFromTree(job,options,1)
        histos.append(hTemp)
        typs.append(typ)    
    elif job.name in data:
        #print 'DATA'
        hTemp, typ = getHistoFromTree(job,options)
        datas.append(hTemp)
        datatyps.append(typ)
        datanames.append(job.name)
'''
for job in info:
    if eval(job.active):
        if job.subsamples:
            for subsample in range(0,len(job.subnames)):
                
                if job.subnames[subsample] in setup:
                    hTemp, typ = getHistoFromTree(job,options,1,subsample)
                    histos.append(hTemp)
                    typs.append(Group[job.subnames[subsample]])


        else:
            if job.name in setup:
                #print job.getpath()
                hTemp, typ = getHistoFromTree(job,options,1)
                histos.append(hTemp)
                typs.append(Group[job.name])

            elif job.name in data:
                #print 'DATA'
                hTemp, typ = getHistoFromTree(job,options)
                datas.append(hTemp)
                datatyps.append(typ)
                datanames.append(job.name)


ROOT.gROOT.SetStyle("Plain")
#import TdrStyles
#TdrStyles.tdrStyle()
c = ROOT.TCanvas(name,title, 700, 600)
c.SetFillStyle(4000)
c.SetFrameFillStyle(1000)
c.SetFrameFillColor(0)

oben = ROOT.TPad('oben','oben',0,0.3 ,1.0,1.0)
oben.SetBottomMargin(0)
oben.SetFillStyle(4000)
oben.SetFrameFillStyle(1000)
oben.SetFrameFillColor(0)
unten = ROOT.TPad('unten','unten',0,0.0,1.0,0.3)
unten.SetTopMargin(0.)
unten.SetBottomMargin(0.35)
unten.SetFillStyle(4000)
unten.SetFrameFillStyle(1000)
unten.SetFrameFillColor(0)

oben.Draw()
unten.Draw()
#c.Divide(1,2)

#c.cd(1)
oben.cd()

allStack = ROOT.THStack(name,title)     
l = ROOT.TLegend(0.75, 0.63, 0.88, 0.88)
MC_integral=0
MC_entries=0

for histo in histos:
    MC_integral+=histo.Integral()
    #MC_entries+=histo.GetEntries()
print "\033[1;32m\n\tMC integral = %s\033[1;m"%MC_integral
#flow = MC_entries-MC_integral
#if flow > 0:
#    print "\033[1;31m\tU/O flow: %s\033[1;m"%flow    

#ORDER AND ADD TOGETHER

#histos, typs = orderandadd(histos,typs,setup)


k=len(histos)
for j in range(0,k):
    #print histos[j].GetBinContent(1)
    i=k-j-1
    histos[i].SetFillColor(int(color[i]))
    histos[i].SetLineColor(1)
    allStack.Add(histos[i])
    l.AddEntry(histos[j],typs[j],'F')
    

d1 = ROOT.TH1F('noData','noData',nBins,xMin,xMax)
datatitle=''
for i in range(0,len(datas)):
    d1.Add(datas[i],1)
    if i ==0:
        datatitle=datanames[i]
    else:
        datatitle=datatitle+ ' + '+datanames[i]
print "\033[1;32m\n\tDATA integral = %s\033[1;m"%d1.Integral()
flow = d1.GetEntries()-d1.Integral()
if flow > 0:
    print "\033[1;31m\tU/O flow: %s\033[1;m"%flow
l.AddEntry(d1,datatitle,'PL')

if Normalize:
    stackscale=d1.Integral()/MC_integral
    stackhists=allStack.GetHists()
    for blabla in stackhists:
        blabla.Scale(stackscale)


allStack.SetTitle()
allStack.Draw("hist")
allStack.GetXaxis().SetTitle(title)
allStack.GetYaxis().SetTitle('Counts')
allStack.GetXaxis().SetRangeUser(xMin,xMax)
allStack.GetYaxis().SetRangeUser(0,20000)
Ymax = max(allStack.GetMaximum(),d1.GetMaximum())*1.3
allStack.SetMaximum(Ymax)
allStack.SetMinimum(0.1)
c.Update()
if config.get('Plot','logy') == '1':
    ROOT.gPad.SetLogy()
ROOT.gPad.SetTicks(1,1)
allStack.Draw("hist")
d1.SetMarkerStyle(21)
d1.Draw("P,E1,X0,same")
l.SetFillColor(0)
l.SetBorderSize(0)
l.Draw()


allMC=ROOT.TH1F('allMC','allMC',nBins,xMin,xMax)
allMC.Sumw2()
for bin in range(0,nBins):
    allMC.SetBinContent(bin,allStack.GetStack().Last().GetBinContent(bin))


t = ROOT.TLatex()
t.SetNDC()
t.SetTextAlign(12)
t.SetTextSize(0.04)
t.DrawLatex(0.13,0.85,"CMS Preliminary")#, BDT Shape")
t.SetTextSize(0.03)
t.DrawLatex(0.13,0.79,"#sqrt{s} = 8 TeV, L = 5.0 fb^{-1}")
#t.DrawLatex(0.13,0.74,"Z(ll)H(b#bar{b})")


unten.cd()
ROOT.gPad.SetTicks(1,1)

ratio, error, ksScore, chiScore = getRatio(d1,allMC,xMin,xMax)
ksScore = allMC.KolmogorovTest( d1 )
chiScore = allMC.Chi2Test( d1 , "UWCHI2/NDF")
print ksScore
print chiScore
ratio.SetStats(0)
ratio.GetYaxis().SetRangeUser(0.5,1.5)
ratio.GetYaxis().SetNdivisions(502,0)
ratio.GetYaxis().SetLabelSize(0.2)
ratio.GetYaxis().SetTitleSize(0.2)
ratio.GetYaxis().SetTitleOffset(0.2)
ratio.GetXaxis().SetLabelColor(10)
ratio.Draw("E1")
ratio.SetTitle("")
m_one_line = ROOT.TLine(xMin,1,xMax,1)
m_one_line.SetLineStyle(7)
m_one_line.SetLineColor(4)
m_one_line.Draw("Same")

t = ROOT.TLatex()
t.SetNDC()
t.SetTextAlign(12)
t.SetTextSize(0.15)
t.DrawLatex(0.12,0.8,"K_{s}: %.2f"%(ksScore))
t.DrawLatex(0.12,0.25,"#chi_{#nu}^{2}: %.2f"%(chiScore))

name = '%s/%s' %(config.get('Directories','plotpath'),options[6])
c.Print(name)
print 'i am done!\n'
sys.exit(0)
Revision:	1.6
Committed:	Fri Aug 10 09:36:16 2012 UTC (12 years, 9 months ago) by peller
Content type:	text/x-python
Branch:	MAIN
Changes since 1.5:	+35 -5 lines
Log Message:	update
#	User	Rev	Content
1	peller	1.1	#!/usr/bin/env python
2			from samplesclass import sample
3			from printcolor import printc
4			import pickle
5			import ROOT
6			from ROOT import TFile, TTree
7			import ROOT
8			from array import array
9	nmohr	1.3	from BetterConfigParser import BetterConfigParser
10	peller	1.1	import sys
11			from mvainfos import mvainfo
12			from gethistofromtree import getHistoFromTree, orderandadd
13			from Ratio import getRatio
14
15			#warnings.filterwarnings( action='ignore', category=RuntimeWarning, message='creating converter.*' )
16
17
18			#load config
19	nmohr	1.3	config = BetterConfigParser()
20	peller	1.6	config.read('./config7TeV_ZZ')
21	peller	1.1
22			#get locations:
23			Wdir=config.get('Directories','Wdir')
24
25
26	peller	1.5	Normalize=False
27	peller	1.1
28			path = sys.argv[1]
29			var = sys.argv[2]
30
31	peller	1.2	if 'bb' in var or 'Light' in var or 'Top' in var:
32			Normalize=True
33	peller	1.1
34	peller	1.5	Normalize=False
35
36
37	peller	1.1	plot=config.get('Plot',var)
38
39			infofile = open(path+'/samples.info','r')
40			info = pickle.load(infofile)
41			infofile.close()
42
43			options = plot.split(',')
44			name=options[1]
45			title = options[2]
46			nBins=int(options[3])
47			xMin=float(options[4])
48			xMax=float(options[5])
49
50			mass=options[9]
51			data=options[10]
52
53			blind=options[11]
54
55			setup=config.get('Plot','setup')
56			setup=setup.split(',')
57
58			color=config.get('Plot','color')
59			color=color.split(',')
60
61
62			weightF=config.get('Weights','weightF')
63	peller	1.6	Group = eval(config.get('LimitGeneral','Group'))
64	peller	1.1
65
66			print '\nProducing Plot of %s\n'%title
67
68
69			histos = []
70			typs = []
71			datas = []
72			datatyps =[]
73			datanames=[]
74	peller	1.6	'''
75	peller	1.1	for job in info:
76			if job.type == 'BKG':
77			#print 'MC'
78	peller	1.5	hTemp, typ = getHistoFromTree(job,options,1)
79	peller	1.1	histos.append(hTemp)
80			typs.append(typ)
81			elif job.type == 'SIG' and job.name == mass:
82	peller	1.5	hTemp, typ = getHistoFromTree(job,options,1)
83	peller	1.1	histos.append(hTemp)
84			typs.append(typ)
85			elif job.name in data:
86			#print 'DATA'
87			hTemp, typ = getHistoFromTree(job,options)
88			datas.append(hTemp)
89			datatyps.append(typ)
90			datanames.append(job.name)
91	peller	1.6	'''
92			for job in info:
93			if eval(job.active):
94			if job.subsamples:
95			for subsample in range(0,len(job.subnames)):
96
97			if job.subnames[subsample] in setup:
98			hTemp, typ = getHistoFromTree(job,options,1,subsample)
99			histos.append(hTemp)
100			typs.append(Group[job.subnames[subsample]])
101
102
103
104			else:
105			if job.name in setup:
106			#print job.getpath()
107			hTemp, typ = getHistoFromTree(job,options,1)
108			histos.append(hTemp)
109			typs.append(Group[job.name])
110
111			elif job.name in data:
112			#print 'DATA'
113			hTemp, typ = getHistoFromTree(job,options)
114			datas.append(hTemp)
115			datatyps.append(typ)
116			datanames.append(job.name)
117
118
119
120	peller	1.1
121
122
123			ROOT.gROOT.SetStyle("Plain")
124	peller	1.6	#import TdrStyles
125			#TdrStyles.tdrStyle()
126	nmohr	1.4	c = ROOT.TCanvas(name,title, 700, 600)
127			c.SetFillStyle(4000)
128			c.SetFrameFillStyle(1000)
129			c.SetFrameFillColor(0)
130
131			oben = ROOT.TPad('oben','oben',0,0.3 ,1.0,1.0)
132			oben.SetBottomMargin(0)
133			oben.SetFillStyle(4000)
134			oben.SetFrameFillStyle(1000)
135			oben.SetFrameFillColor(0)
136			unten = ROOT.TPad('unten','unten',0,0.0,1.0,0.3)
137			unten.SetTopMargin(0.)
138			unten.SetBottomMargin(0.35)
139			unten.SetFillStyle(4000)
140			unten.SetFrameFillStyle(1000)
141			unten.SetFrameFillColor(0)
142	peller	1.1
143			oben.Draw()
144			unten.Draw()
145			#c.Divide(1,2)
146
147			#c.cd(1)
148			oben.cd()
149
150			allStack = ROOT.THStack(name,title)
151			l = ROOT.TLegend(0.75, 0.63, 0.88, 0.88)
152			MC_integral=0
153			MC_entries=0
154
155			for histo in histos:
156			MC_integral+=histo.Integral()
157			#MC_entries+=histo.GetEntries()
158			print "\033[1;32m\n\tMC integral = %s\033[1;m"%MC_integral
159			#flow = MC_entries-MC_integral
160			#if flow > 0:
161			# print "\033[1;31m\tU/O flow: %s\033[1;m"%flow
162
163			#ORDER AND ADD TOGETHER
164
165	peller	1.6	#histos, typs = orderandadd(histos,typs,setup)
166	peller	1.1
167
168			k=len(histos)
169			for j in range(0,k):
170			#print histos[j].GetBinContent(1)
171			i=k-j-1
172			histos[i].SetFillColor(int(color[i]))
173	peller	1.5	histos[i].SetLineColor(1)
174	peller	1.1	allStack.Add(histos[i])
175			l.AddEntry(histos[j],typs[j],'F')
176
177
178			d1 = ROOT.TH1F('noData','noData',nBins,xMin,xMax)
179			datatitle=''
180			for i in range(0,len(datas)):
181			d1.Add(datas[i],1)
182			if i ==0:
183			datatitle=datanames[i]
184			else:
185			datatitle=datatitle+ ' + '+datanames[i]
186			print "\033[1;32m\n\tDATA integral = %s\033[1;m"%d1.Integral()
187			flow = d1.GetEntries()-d1.Integral()
188			if flow > 0:
189			print "\033[1;31m\tU/O flow: %s\033[1;m"%flow
190			l.AddEntry(d1,datatitle,'PL')
191	peller	1.2
192			if Normalize:
193			stackscale=d1.Integral()/MC_integral
194			stackhists=allStack.GetHists()
195			for blabla in stackhists:
196			blabla.Scale(stackscale)
197
198
199
200	peller	1.1	allStack.SetTitle()
201			allStack.Draw("hist")
202			allStack.GetXaxis().SetTitle(title)
203			allStack.GetYaxis().SetTitle('Counts')
204			allStack.GetXaxis().SetRangeUser(xMin,xMax)
205			allStack.GetYaxis().SetRangeUser(0,20000)
206			Ymax = max(allStack.GetMaximum(),d1.GetMaximum())*1.3
207			allStack.SetMaximum(Ymax)
208			allStack.SetMinimum(0.1)
209			c.Update()
210			if config.get('Plot','logy') == '1':
211			ROOT.gPad.SetLogy()
212			ROOT.gPad.SetTicks(1,1)
213			allStack.Draw("hist")
214			d1.SetMarkerStyle(21)
215			d1.Draw("P,E1,X0,same")
216			l.SetFillColor(0)
217			l.SetBorderSize(0)
218			l.Draw()
219
220
221			allMC=ROOT.TH1F('allMC','allMC',nBins,xMin,xMax)
222			allMC.Sumw2()
223			for bin in range(0,nBins):
224			allMC.SetBinContent(bin,allStack.GetStack().Last().GetBinContent(bin))
225
226
227
228
229			t = ROOT.TLatex()
230			t.SetNDC()
231			t.SetTextAlign(12)
232			t.SetTextSize(0.04)
233			t.DrawLatex(0.13,0.85,"CMS Preliminary")#, BDT Shape")
234			t.SetTextSize(0.03)
235	peller	1.5	t.DrawLatex(0.13,0.79,"#sqrt{s} = 8 TeV, L = 5.0 fb^{-1}")
236	peller	1.1	#t.DrawLatex(0.13,0.74,"Z(ll)H(b#bar{b})")
237
238
239			unten.cd()
240			ROOT.gPad.SetTicks(1,1)
241
242			ratio, error, ksScore, chiScore = getRatio(d1,allMC,xMin,xMax)
243	nmohr	1.4	ksScore = allMC.KolmogorovTest( d1 )
244			chiScore = allMC.Chi2Test( d1 , "UWCHI2/NDF")
245	peller	1.1	print ksScore
246			print chiScore
247			ratio.SetStats(0)
248			ratio.GetYaxis().SetRangeUser(0.5,1.5)
249			ratio.GetYaxis().SetNdivisions(502,0)
250			ratio.GetYaxis().SetLabelSize(0.2)
251			ratio.GetYaxis().SetTitleSize(0.2)
252			ratio.GetYaxis().SetTitleOffset(0.2)
253			ratio.GetXaxis().SetLabelColor(10)
254			ratio.Draw("E1")
255			ratio.SetTitle("")
256			m_one_line = ROOT.TLine(xMin,1,xMax,1)
257			m_one_line.SetLineStyle(7)
258			m_one_line.SetLineColor(4)
259			m_one_line.Draw("Same")
260
261	peller	1.2	t = ROOT.TLatex()
262			t.SetNDC()
263			t.SetTextAlign(12)
264			t.SetTextSize(0.15)
265	nmohr	1.4	t.DrawLatex(0.12,0.8,"K_{s}: %.2f"%(ksScore))
266			t.DrawLatex(0.12,0.25,"#chi_{#nu}^{2}: %.2f"%(chiScore))
267	peller	1.1
268			name = '%s/%s' %(config.get('Directories','plotpath'),options[6])
269	nmohr	1.3	c.Print(name)
270	peller	1.5	print 'i am done!\n'
271			sys.exit(0)