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, os
from mvainfos import mvainfo
#from gethistofromtree import getHistoFromTree, orderandadd
from Ratio import getRatio
from optparse import OptionParser
from HistoMaker import HistoMaker, orderandadd

#CONFIGURE
argv = sys.argv
parser = OptionParser()
parser.add_option("-P", "--path", dest="path", default="",
                      help="path to samples")
parser.add_option("-R", "--reg", dest="region", default="",
                      help="region to plot")
parser.add_option("-C", "--config", dest="config", default=[], action="append",
                      help="configuration file")
(opts, args) = parser.parse_args(argv)
if opts.config =="":
        opts.config = "config"
print opts.config
config = BetterConfigParser()
config.read(opts.config)
anaTag = config.get("Analysis","tag")

path = opts.path
region = opts.region

plotConfig = BetterConfigParser()
plotConfig.read('vhbbPlotDef.ini')

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

section='Plot:%s'%region

Normalize = eval(config.get(section,'Normalize'))
log = eval(config.get(section,'log'))
blind = eval(config.get(section,'blind'))

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

#options = plot.split(',')

mass = config.get(section,'Signal')

vars = (config.get(section, 'vars')).split(',')

names = [plotConfig.get('plotDef:%s'%x,'relPath') for x in vars]
nBins = [eval(plotConfig.get('plotDef:%s'%x,'nBins')) for x in vars]
xMin = [eval(plotConfig.get('plotDef:%s'%x,'min')) for x in vars]
xMax = [eval(plotConfig.get('plotDef:%s'%x,'max')) for x in vars]
xAxis = [plotConfig.get('plotDef:%s'%x,'xAxis') for x in vars]

for p in range(0,len(names)):
    if '<mass>' in names[p]:
        newp= names[p].replace('<mass>',mass)
        names[p]=newp
    print names[p]

data = config.get(section,'Datas')
if config.has_option(section, 'Datacut'):
    datacut=config.get(section, 'Datacut')
else:
    datacut = region

options=[]

if blind: blindopt='blind'
else: blindopt = 'noblind'

for i in range(0,len(vars)):
    options.append([names[i],'',xAxis[i],nBins[i],xMin[i],xMax[i],'%s_%s.pdf'%(region,vars[i]),region,datacut,mass,data,blindopt])


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

samples=config.get('Plot_general','samples')
samples=samples.split(',')

colorDict=eval(config.get('Plot_general','colorDict'))
#color=color.split(',')


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


#GETALL AT ONCE

Plotter=HistoMaker(path,config,region,options)
 
#print '\nProducing Plot of %s\n'%vars[v]
Lhistos = [[] for _ in range(0,len(vars))]
Ltyps = [[] for _ in range(0,len(vars))]
Ldatas = [[] for _ in range(0,len(vars))]
Ldatatyps = [[] for _ in range(0,len(vars))]
Ldatanames = [[] for _ in range(0,len(vars))]


#Find out Lumi:
for job in info:
    if job.name in data: lumi_data=float(job.lumi)

Plotter.lumi=lumi_data

for job in info:
    if eval(job.active):
        if job.subsamples:
            for subsample in range(0,len(job.subnames)):
                
                if job.subnames[subsample] in samples:
                    hTempList, typList = Plotter.getHistoFromTree(job,subsample)
                    for v in range(0,len(vars)):
                        Lhistos[v].append(hTempList[v])
                        Ltyps[v].append(Group[job.subnames[subsample]])

        else:
            if job.name in samples:
                #print job.getpath()
                hTempList, typList = Plotter.getHistoFromTree(job)
                for v in range(0,len(vars)):
                    Lhistos[v].append(hTempList[v])
                    Ltyps[v].append(Group[job.name])

            elif job.name in data:
                #print 'DATA'
                hTemp, typ = Plotter.getHistoFromTree(job)
                for v in range(0,len(vars)):
                    Ldatas[v].append(hTemp[v])
                    Ldatatyps[v].append(typ[v])
                    Ldatanames[v].append(job.name)


for v in range(0,len(vars)):

    histos = Lhistos[v]
    typs = Ltyps[v]
    datas = Ldatas[v]
    datatyps = Ldatatyps[v]
    datanames= Ldatanames[v]

    ROOT.gROOT.SetStyle("Plain")
    
    c = ROOT.TCanvas(vars[v],'', 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()

    oben.cd()

    allStack = ROOT.THStack(vars[v],'')     
    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

    #ORDER AND ADD TOGETHER
    #print typs
    #print setup
    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(colorDict[setup[i]]))
        histos[i].SetLineColor(1)
        allStack.Add(histos[i])
        l.AddEntry(histos[j],typs[j],'F')

    d1 = ROOT.TH1F('noData','noData',nBins[v],xMin[v],xMax[v])
    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('')
    allStack.GetYaxis().SetTitle('Counts')
    allStack.GetXaxis().SetRangeUser(xMin[v],xMax[v])
    allStack.GetYaxis().SetRangeUser(0,20000)
    Ymax = max(allStack.GetMaximum(),d1.GetMaximum())*1.3
    allStack.SetMaximum(Ymax)
    allStack.SetMinimum(0.1)
    c.Update()
    if log:
        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[v],xMin[v],xMax[v])
    allMC.Sumw2()
    for bin in range(0,nBins[v]):
        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")
    t.SetTextSize(0.03)
    t.DrawLatex(0.13,0.79,"#sqrt{s} =  %s, L = %s fb^{-1}"%(anaTag,(float(lumi_data)/1000.)))

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

    ratio, error, ksScore, chiScore = getRatio(d1,allMC,xMin[v],xMax[v])
    ksScore = allMC.KolmogorovTest( d1 )
    chiScore = allMC.Chi2Test( d1 , "UWCHI2/NDF")
    print ksScore
    print chiScore
    ratio.SetStats(0)
    ratio.GetYaxis().SetRangeUser(0,2)
    ratio.GetYaxis().SetNdivisions(502,0)
    ratio.GetXaxis().SetTitle(xAxis[v])
    ratio.Draw("E1")
    ratio.SetTitle("")
    m_one_line = ROOT.TLine(xMin[v],1,xMax[v],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.07)
    t.DrawLatex(0.12,0.9,"K_{s}: %.2f"%(ksScore))
    t.DrawLatex(0.12,0.4,"#chi_{#nu}^{2}: %.2f"%(chiScore))

    name = '%s/%s' %(config.get('Directories','plotpath'),options[v][6])
    c.Print(name)

    os.system('rm %s/tmp_plotCache_%s*'%(config.get('Directories','plotpath'),region))
    print 'i am done!\n'

sys.exit(0)
Revision:	1.13
Committed:	Wed Oct 3 09:51:40 2012 UTC (12 years, 7 months ago) by nmohr
Content type:	text/x-python
Branch:	MAIN
Changes since 1.12:	+2 -2 lines
Log Message:	Float
#	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.11	import sys, os
11	peller	1.1	from mvainfos import mvainfo
12	peller	1.10	#from gethistofromtree import getHistoFromTree, orderandadd
13	peller	1.1	from Ratio import getRatio
14	peller	1.8	from optparse import OptionParser
15	peller	1.10	from HistoMaker import HistoMaker, orderandadd
16	peller	1.1
17	peller	1.7	#CONFIGURE
18			argv = sys.argv
19			parser = OptionParser()
20			parser.add_option("-P", "--path", dest="path", default="",
21			help="path to samples")
22	peller	1.10	parser.add_option("-R", "--reg", dest="region", default="",
23			help="region to plot")
24	peller	1.7	parser.add_option("-C", "--config", dest="config", default=[], action="append",
25			help="configuration file")
26			(opts, args) = parser.parse_args(argv)
27			if opts.config =="":
28			opts.config = "config"
29			print opts.config
30			config = BetterConfigParser()
31			config.read(opts.config)
32			anaTag = config.get("Analysis","tag")
33
34			path = opts.path
35	peller	1.10	region = opts.region
36	peller	1.7
37	peller	1.10	plotConfig = BetterConfigParser()
38			plotConfig.read('vhbbPlotDef.ini')
39	peller	1.7
40	peller	1.10	#get locations:
41			Wdir=config.get('Directories','Wdir')
42	peller	1.1
43	peller	1.10	section='Plot:%s'%region
44	peller	1.1
45	peller	1.10	Normalize = eval(config.get(section,'Normalize'))
46			log = eval(config.get(section,'log'))
47			blind = eval(config.get(section,'blind'))
48
49			infofile = open(path+'/samples.info','r')
50			info = pickle.load(infofile)
51			infofile.close()
52	peller	1.1
53	peller	1.10	#options = plot.split(',')
54	peller	1.1
55	peller	1.10	mass = config.get(section,'Signal')
56	peller	1.1
57	peller	1.10	vars = (config.get(section, 'vars')).split(',')
58	peller	1.1
59	peller	1.10	names = [plotConfig.get('plotDef:%s'%x,'relPath') for x in vars]
60			nBins = [eval(plotConfig.get('plotDef:%s'%x,'nBins')) for x in vars]
61			xMin = [eval(plotConfig.get('plotDef:%s'%x,'min')) for x in vars]
62			xMax = [eval(plotConfig.get('plotDef:%s'%x,'max')) for x in vars]
63			xAxis = [plotConfig.get('plotDef:%s'%x,'xAxis') for x in vars]
64	peller	1.5
65	peller	1.10	for p in range(0,len(names)):
66			if '<mass>' in names[p]:
67			newp= names[p].replace('<mass>',mass)
68			names[p]=newp
69			print names[p]
70	peller	1.5
71	peller	1.10	data = config.get(section,'Datas')
72			if config.has_option(section, 'Datacut'):
73			datacut=config.get(section, 'Datacut')
74			else:
75			datacut = region
76	peller	1.1
77	peller	1.10	options=[]
78	peller	1.1
79	peller	1.10	if blind: blindopt='blind'
80			else: blindopt = 'noblind'
81	peller	1.1
82	peller	1.10	for i in range(0,len(vars)):
83			options.append([names[i],'',xAxis[i],nBins[i],xMin[i],xMax[i],'%s_%s.pdf'%(region,vars[i]),region,datacut,mass,data,blindopt])
84	peller	1.1
85
86	peller	1.11	setup=config.get('Plot_general','setup')
87	peller	1.1	setup=setup.split(',')
88
89	peller	1.11	samples=config.get('Plot_general','samples')
90	peller	1.8	samples=samples.split(',')
91
92	peller	1.11	colorDict=eval(config.get('Plot_general','colorDict'))
93	peller	1.8	#color=color.split(',')
94	peller	1.1
95
96			weightF=config.get('Weights','weightF')
97	peller	1.6	Group = eval(config.get('LimitGeneral','Group'))
98	peller	1.1
99
100	peller	1.10	#GETALL AT ONCE
101	peller	1.1
102	peller	1.11	Plotter=HistoMaker(path,config,region,options)
103	peller	1.10
104			#print '\nProducing Plot of %s\n'%vars[v]
105			Lhistos = [[] for _ in range(0,len(vars))]
106			Ltyps = [[] for _ in range(0,len(vars))]
107			Ldatas = [[] for _ in range(0,len(vars))]
108			Ldatatyps = [[] for _ in range(0,len(vars))]
109			Ldatanames = [[] for _ in range(0,len(vars))]
110	peller	1.1
111	peller	1.12
112			#Find out Lumi:
113			for job in info:
114	nmohr	1.13	if job.name in data: lumi_data=float(job.lumi)
115	peller	1.12
116			Plotter.lumi=lumi_data
117
118	peller	1.6	for job in info:
119			if eval(job.active):
120			if job.subsamples:
121			for subsample in range(0,len(job.subnames)):
122
123	peller	1.8	if job.subnames[subsample] in samples:
124	peller	1.10	hTempList, typList = Plotter.getHistoFromTree(job,subsample)
125			for v in range(0,len(vars)):
126			Lhistos[v].append(hTempList[v])
127			Ltyps[v].append(Group[job.subnames[subsample]])
128	peller	1.6
129			else:
130	peller	1.8	if job.name in samples:
131	peller	1.6	#print job.getpath()
132	peller	1.10	hTempList, typList = Plotter.getHistoFromTree(job)
133			for v in range(0,len(vars)):
134			Lhistos[v].append(hTempList[v])
135			Ltyps[v].append(Group[job.name])
136	peller	1.6
137			elif job.name in data:
138			#print 'DATA'
139	peller	1.10	hTemp, typ = Plotter.getHistoFromTree(job)
140			for v in range(0,len(vars)):
141			Ldatas[v].append(hTemp[v])
142			Ldatatyps[v].append(typ[v])
143			Ldatanames[v].append(job.name)
144
145
146			for v in range(0,len(vars)):
147
148			histos = Lhistos[v]
149			typs = Ltyps[v]
150			datas = Ldatas[v]
151			datatyps = Ldatatyps[v]
152			datanames= Ldatanames[v]
153	peller	1.6
154	peller	1.10	ROOT.gROOT.SetStyle("Plain")
155	peller	1.1
156	peller	1.10	c = ROOT.TCanvas(vars[v],'', 700, 600)
157			c.SetFillStyle(4000)
158			c.SetFrameFillStyle(1000)
159			c.SetFrameFillColor(0)
160
161			oben = ROOT.TPad('oben','oben',0,0.3 ,1.0,1.0)
162			oben.SetBottomMargin(0)
163			oben.SetFillStyle(4000)
164			oben.SetFrameFillStyle(1000)
165			oben.SetFrameFillColor(0)
166			unten = ROOT.TPad('unten','unten',0,0.0,1.0,0.3)
167			unten.SetTopMargin(0.)
168			unten.SetBottomMargin(0.35)
169			unten.SetFillStyle(4000)
170			unten.SetFrameFillStyle(1000)
171			unten.SetFrameFillColor(0)
172
173			oben.Draw()
174			unten.Draw()
175
176			oben.cd()
177
178			allStack = ROOT.THStack(vars[v],'')
179			l = ROOT.TLegend(0.75, 0.63, 0.88, 0.88)
180			MC_integral=0
181			MC_entries=0
182
183			for histo in histos:
184			MC_integral+=histo.Integral()
185			#MC_entries+=histo.GetEntries()
186			print "\033[1;32m\n\tMC integral = %s\033[1;m"%MC_integral
187
188			#ORDER AND ADD TOGETHER
189	peller	1.11	#print typs
190			#print setup
191	peller	1.10	histos, typs = orderandadd(histos,typs,setup)
192
193
194			k=len(histos)
195			for j in range(0,k):
196			#print histos[j].GetBinContent(1)
197			i=k-j-1
198			histos[i].SetFillColor(int(colorDict[setup[i]]))
199			histos[i].SetLineColor(1)
200			allStack.Add(histos[i])
201			l.AddEntry(histos[j],typs[j],'F')
202
203			d1 = ROOT.TH1F('noData','noData',nBins[v],xMin[v],xMax[v])
204			datatitle=''
205			for i in range(0,len(datas)):
206			d1.Add(datas[i],1)
207			if i ==0:
208			datatitle=datanames[i]
209			else:
210			datatitle=datatitle+ ' + '+datanames[i]
211			print "\033[1;32m\n\tDATA integral = %s\033[1;m"%d1.Integral()
212			flow = d1.GetEntries()-d1.Integral()
213			if flow > 0:
214			print "\033[1;31m\tU/O flow: %s\033[1;m"%flow
215			l.AddEntry(d1,datatitle,'PL')
216
217			if Normalize:
218			stackscale=d1.Integral()/MC_integral
219			stackhists=allStack.GetHists()
220			for blabla in stackhists:
221			blabla.Scale(stackscale)
222
223			allStack.SetTitle()
224			allStack.Draw("hist")
225			allStack.GetXaxis().SetTitle('')
226			allStack.GetYaxis().SetTitle('Counts')
227			allStack.GetXaxis().SetRangeUser(xMin[v],xMax[v])
228			allStack.GetYaxis().SetRangeUser(0,20000)
229			Ymax = max(allStack.GetMaximum(),d1.GetMaximum())*1.3
230			allStack.SetMaximum(Ymax)
231			allStack.SetMinimum(0.1)
232			c.Update()
233			if log:
234			ROOT.gPad.SetLogy()
235			ROOT.gPad.SetTicks(1,1)
236			allStack.Draw("hist")
237			d1.SetMarkerStyle(21)
238			d1.Draw("P,E1,X0,same")
239			l.SetFillColor(0)
240			l.SetBorderSize(0)
241			l.Draw()
242
243
244			allMC=ROOT.TH1F('allMC','allMC',nBins[v],xMin[v],xMax[v])
245			allMC.Sumw2()
246			for bin in range(0,nBins[v]):
247			allMC.SetBinContent(bin,allStack.GetStack().Last().GetBinContent(bin))
248
249			t = ROOT.TLatex()
250			t.SetNDC()
251			t.SetTextAlign(12)
252			t.SetTextSize(0.04)
253			t.DrawLatex(0.13,0.85,"CMS Preliminary")
254			t.SetTextSize(0.03)
255	nmohr	1.13	t.DrawLatex(0.13,0.79,"#sqrt{s} = %s, L = %s fb^{-1}"%(anaTag,(float(lumi_data)/1000.)))
256	peller	1.10
257			unten.cd()
258			ROOT.gPad.SetTicks(1,1)
259
260			ratio, error, ksScore, chiScore = getRatio(d1,allMC,xMin[v],xMax[v])
261			ksScore = allMC.KolmogorovTest( d1 )
262			chiScore = allMC.Chi2Test( d1 , "UWCHI2/NDF")
263			print ksScore
264			print chiScore
265			ratio.SetStats(0)
266			ratio.GetYaxis().SetRangeUser(0,2)
267			ratio.GetYaxis().SetNdivisions(502,0)
268			ratio.GetXaxis().SetTitle(xAxis[v])
269			ratio.Draw("E1")
270			ratio.SetTitle("")
271			m_one_line = ROOT.TLine(xMin[v],1,xMax[v],1)
272			m_one_line.SetLineStyle(7)
273			m_one_line.SetLineColor(4)
274			m_one_line.Draw("Same")
275
276			t = ROOT.TLatex()
277			t.SetNDC()
278			t.SetTextAlign(12)
279			t.SetTextSize(0.07)
280			t.DrawLatex(0.12,0.9,"K_{s}: %.2f"%(ksScore))
281			t.DrawLatex(0.12,0.4,"#chi_{#nu}^{2}: %.2f"%(chiScore))
282
283			name = '%s/%s' %(config.get('Directories','plotpath'),options[v][6])
284			c.Print(name)
285	peller	1.11
286			os.system('rm %s/tmp_plotCache_%s*'%(config.get('Directories','plotpath'),region))
287	peller	1.10	print 'i am done!\n'
288	peller	1.11
289	peller	1.5	sys.exit(0)