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
import TdrStyles

#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")

TrainFlag = eval(config.get('Analysis','TrainFlag'))

if TrainFlag:
    MC_rescale_factor=2.
    print 'I RESCALE BY 2.0'
else: MC_rescale_factor = 1.

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('Plot_general','Group'))


#GETALL AT ONCE

Plotter=HistoMaker(path,config,region,options,MC_rescale_factor)

#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))]

def myText(txt="CMS Preliminary",ndcX=0,ndcY=0,size=0.8):
    ROOT.gPad.Update()
    text = ROOT.TLatex()
    text.SetNDC()
    text.SetTextColor(ROOT.kBlack)
    text.SetTextSize(text.GetTextSize()*size)
    text.DrawLatex(ndcX,ndcY,txt)
    return text


#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]

    TdrStyles.tdrStyle()
    
    c = ROOT.TCanvas(vars[v],'', 600, 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:
        if MC_integral != 0:    stackscale=d1.Integral()/MC_integral
        stackhists=allStack.GetHists()
        for blabla in stackhists:
                if MC_integral != 0: blabla.Scale(stackscale)
    
    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))
        allMC.SetBinError(bin,allStack.GetStack().Last().GetBinError(bin))

    allStack.SetTitle()
    allStack.Draw("hist")
    allStack.GetXaxis().SetTitle('')
    yTitle = 'Entries'
    if not '/' in yTitle:
            yAppend = '%s' %(allStack.GetXaxis().GetBinWidth(1)) 
            yTitle = '%s / %s' %(yTitle, yAppend)
    allStack.GetYaxis().SetTitle(yTitle)
    allStack.GetXaxis().SetRangeUser(xMin[v],xMax[v])
    allStack.GetYaxis().SetRangeUser(0,20000)
    theErrorGraph = ROOT.TGraphErrors(allMC)
    theErrorGraph.SetFillColor(ROOT.kGray+3)
    theErrorGraph.SetFillStyle(3013)
    theErrorGraph.Draw('SAME2')
    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.Draw("E0same")
    l.SetFillColor(0)
    l.SetBorderSize(0)
    l.Draw()


    tPrel = myText("CMS Preliminary",0.17,0.88,1.04)
    tLumi = myText("#sqrt{s} =  %s, L = %s fb^{-1}"%(anaTag,(float(lumi_data)/1000.)),0.17,0.83)

    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")

    tKsChi = myText("#chi_{#nu}^{2} = %.3f K_{s} = %.3f"%(chiScore,ksScore),0.17,0.9,1.5)

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