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

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

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


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