VHbb/python/tree_stack.py

#!/usr/bin/env python
from samplesclass import sample
from printcolor import printc
import pickle
import ROOT 
from array import array
from BetterConfigParser import BetterConfigParser
import sys, os
from mvainfos import mvainfo
#from gethistofromtree import getHistoFromTree, orderandadd
from optparse import OptionParser
from HistoMaker import HistoMaker, orderandadd
from copy import copy,deepcopy
from StackMaker import StackMaker
from math import sqrt

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

path = opts.path
region = opts.region


#get locations:
Wdir=config.get('Directories','Wdir')
samplesinfo=config.get('Directories','samplesinfo')
#limitpath=config.get('Directories','limits')
limitpath=path
section='Plot:%s'%region

infofile = open(samplesinfo,'r')
info = pickle.load(infofile)
infofile.close()

if 'vhbb_TH_BDT' in region:
#-----------Histo from TH File------------------------------------
    if 'Zee' in region: d='Zee'
    elif 'Zmm' in region: d='Zmm'
    if 'LowPt' in region:
        var='BDT8_RMed'
        newregion='LowPt_%s'%d
    elif 'HighPtLooseBTag' in region:
        var='BDT8_RTightLooseBTag'
        newregion='HighPtLooseBTag_%s'%d
    elif 'HighPt' in region:
        var='BDT8_RTight'
        newregion='HighPt_%s'%d
    Stack=StackMaker(config,var,newregion,True)

    log = eval(config.get('Plot:%s'%newregion,'log'))

    if log:
        setup = config.get('Plot_general','setupLog').split(',')
    else:
        setup = config.get('Plot_general','setup').split(',')
    Dict = eval(config.get('LimitGeneral','Dict'))

    setup.remove('DYc')

    sys_BDT= eval(config.get('LimitGeneral','sys_BDT'))
    systematicsnaming8TeV = eval(config.get('LimitGeneral','systematicsnaming8TeV'))
    systs=[systematicsnaming8TeV[s] for s in sys_BDT]
    if eval(config.get('LimitGeneral','weightF_sys')): systs.append('UEPS')
    input = ROOT.TFile.Open(limitpath+'/'+region+'.root','read')

    lumi=0
    for job in info:
        if job.name == d:
            lumi=job.lumi
            break
        else: pass

    histos = []
    typs = []

    #systs=[] 

    setup2=copy(setup)
    setup2.remove('ZH')

    shapesUp = [[] for _ in range(0,len(setup2))]
    shapesDown = [[] for _ in range(0,len(setup2))]
    for s in setup:
        if 'ZH' == s:
            Overlay=copy(input.Get(Dict[s]))
        else:
            histos.append(input.Get(Dict[s]))
            typs.append(s)
            print s
            for syst in systs:
                shapesUp[setup2.index(s)].append(input.Get(Dict[s]+syst+'Up'))
                shapesDown[setup2.index(s)].append(input.Get(Dict[s]+syst+'Down'))

    #print shapesUp

    #calculate the Errors
    counter = 0
    total=[]
    errUp=[]
    errDown=[]
    print 'total bins %s'%histos[0].GetNbinsX()
    for h in range(0,len(histos)):
        if counter == 0:
            Error = ROOT.TGraphAsymmErrors(histos[h])
        for bin in range(1,histos[h].GetNbinsX()+1):
            if counter == 0 and h == 0:
                total.append(0)
                errUp.append(0)
                errDown.append(0)
            point=histos[h].GetXaxis().GetBinCenter(bin)
            total[bin-1]+=histos[h].GetBinContent(bin)
            for i in range(0,len(shapesUp[h])):
                errUp[bin-1]+=(shapesUp[h][i].GetBinContent(bin)-histos[h].GetBinContent(bin))**2
                errDown[bin-1]+=(shapesDown[h][i].GetBinContent(bin)-histos[h].GetBinContent(bin))**2
            errUp[bin-1]+=(histos[h].GetBinError(bin))**2
            errDown[bin-1]+=(histos[h].GetBinError(bin))**2
            
            Error.SetPoint(bin-1,point,1)
        counter += 1
        
    for bin in range(0,len(total)):
        if not total[bin] == 0: 
            Error.SetPointEYlow(bin,sqrt(errDown[bin])/total[bin])
            Error.SetPointEYhigh(bin,sqrt(errUp[bin])/total[bin])

        #if counter == 0: Errors=Error
        #else: Errors.Add(Error)
    #set the errors

    #get the Ratios


    datas=[input.Get('data_obs')]
    datatyps = [None]
    datanames=[d] 


    histos.append(copy(Overlay))
    typs.append('ZH')
    #histos.append(copy(Overlay))
    #typs.append('ZH')

    Stack.histos = histos
    Stack.typs = typs
    Stack.datas = datas
    Stack.datatyps = datatyps
    Stack.datanames= datanames
    Stack.overlay = Overlay
    Stack.AddErrors=Error
    Stack.lumi = lumi
    Stack.doPlot()

    print 'i am done!\n'
#-------------------------------------------------


else:
#----------Histo from trees------------
    vars = (config.get(section, 'vars')).split(',')

    if 'ZLight' in region or 'TTbar' in region or 'Zbb' in region: SignalRegion = False
    else:
        SignalRegion = True
        print 'You are in the Signal Region!'

    data = config.get(section,'Datas')

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

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

    #GETALL AT ONCE
    options = []
    Stacks = []
    for i in range(len(vars)):
        Stacks.append(StackMaker(config,vars[i],region,SignalRegion))
        options.append(Stacks[i].options)

    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:
    lumicounter=0.
    lumi=0.
    for job in info:
        if job.name in data:
            lumi+=float(job.lumi)
            lumicounter+=1.

    if lumicounter > 0:
        lumi=lumi/lumicounter

    Plotter.lumi=lumi
    mass = Stacks[0].mass

    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]])
                            print job.subnames[subsample]

            else:
                if job.name in samples:
                    if job.name == mass:
                        print job.name
                        hTempList, typList = Plotter.getHistoFromTree(job)
                        for v in range(0,len(vars)):
                            if SignalRegion:
                                Lhistos[v].append(hTempList[v])
                                Ltyps[v].append(Group[job.name])
                            Overlaylist= deepcopy(hTempList)
                                                                                                                                 
                    else:
                        print job.name
                        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]
        Stacks[v].histos = Lhistos[v]
        Stacks[v].typs = Ltyps[v]
        Stacks[v].datas = Ldatas[v]
        Stacks[v].datatyps = Ldatatyps[v]
        Stacks[v].datanames= Ldatanames[v]
        Stacks[v].overlay = Overlaylist[v]
        Stacks[v].lumi = lumi
        Stacks[v].doPlot()
        print 'i am done!\n'
#----------------------------------------------------
sys.exit(0)
Revision:	1.33
Committed:	Mon Oct 22 14:33:42 2012 UTC (12 years, 6 months ago) by peller
Content type:	text/x-python
Branch:	MAIN
Changes since 1.32:	+19 -11 lines
Log Message:	bin trafo
#	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 array import array
7	nmohr	1.3	from BetterConfigParser import BetterConfigParser
8	peller	1.11	import sys, os
9	peller	1.1	from mvainfos import mvainfo
10	peller	1.10	#from gethistofromtree import getHistoFromTree, orderandadd
11	peller	1.8	from optparse import OptionParser
12	peller	1.10	from HistoMaker import HistoMaker, orderandadd
13	peller	1.32	from copy import copy,deepcopy
14	nmohr	1.31	from StackMaker import StackMaker
15	peller	1.32	from math import sqrt
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	peller	1.17
33	peller	1.7	path = opts.path
34	peller	1.10	region = opts.region
35	peller	1.7
36	peller	1.32
37	peller	1.10	#get locations:
38			Wdir=config.get('Directories','Wdir')
39	nmohr	1.28	samplesinfo=config.get('Directories','samplesinfo')
40	peller	1.32	#limitpath=config.get('Directories','limits')
41			limitpath=path
42	peller	1.10	section='Plot:%s'%region
43	peller	1.1
44	nmohr	1.28	infofile = open(samplesinfo,'r')
45	peller	1.10	info = pickle.load(infofile)
46			infofile.close()
47	peller	1.1
48	peller	1.32	if 'vhbb_TH_BDT' in region:
49			#-----------Histo from TH File------------------------------------
50			if 'Zee' in region: d='Zee'
51			elif 'Zmm' in region: d='Zmm'
52			if 'LowPt' in region:
53			var='BDT8_RMed'
54			newregion='LowPt_%s'%d
55	peller	1.33	elif 'HighPtLooseBTag' in region:
56			var='BDT8_RTightLooseBTag'
57			newregion='HighPtLooseBTag_%s'%d
58	peller	1.32	elif 'HighPt' in region:
59			var='BDT8_RTight'
60			newregion='HighPt_%s'%d
61			Stack=StackMaker(config,var,newregion,True)
62
63			log = eval(config.get('Plot:%s'%newregion,'log'))
64
65			if log:
66			setup = config.get('Plot_general','setupLog').split(',')
67			else:
68			setup = config.get('Plot_general','setup').split(',')
69			Dict = eval(config.get('LimitGeneral','Dict'))
70
71			setup.remove('DYc')
72
73			sys_BDT= eval(config.get('LimitGeneral','sys_BDT'))
74			systematicsnaming8TeV = eval(config.get('LimitGeneral','systematicsnaming8TeV'))
75			systs=[systematicsnaming8TeV[s] for s in sys_BDT]
76			if eval(config.get('LimitGeneral','weightF_sys')): systs.append('UEPS')
77			input = ROOT.TFile.Open(limitpath+'/'+region+'.root','read')
78
79			lumi=0
80			for job in info:
81			if job.name == d:
82			lumi=job.lumi
83			break
84			else: pass
85
86			histos = []
87			typs = []
88
89			#systs=[]
90
91	peller	1.33	setup2=copy(setup)
92			setup2.remove('ZH')
93
94			shapesUp = [[] for _ in range(0,len(setup2))]
95			shapesDown = [[] for _ in range(0,len(setup2))]
96	peller	1.32	for s in setup:
97	peller	1.33	if 'ZH' == s:
98	peller	1.32	Overlay=copy(input.Get(Dict[s]))
99	peller	1.33	else:
100			histos.append(input.Get(Dict[s]))
101			typs.append(s)
102			print s
103			for syst in systs:
104			shapesUp[setup2.index(s)].append(input.Get(Dict[s]+syst+'Up'))
105			shapesDown[setup2.index(s)].append(input.Get(Dict[s]+syst+'Down'))
106	peller	1.32
107			#print shapesUp
108
109			#calculate the Errors
110			counter = 0
111			total=[]
112			errUp=[]
113			errDown=[]
114			print 'total bins %s'%histos[0].GetNbinsX()
115			for h in range(0,len(histos)):
116			if counter == 0:
117			Error = ROOT.TGraphAsymmErrors(histos[h])
118			for bin in range(1,histos[h].GetNbinsX()+1):
119			if counter == 0 and h == 0:
120			total.append(0)
121			errUp.append(0)
122			errDown.append(0)
123			point=histos[h].GetXaxis().GetBinCenter(bin)
124			total[bin-1]+=histos[h].GetBinContent(bin)
125			for i in range(0,len(shapesUp[h])):
126			errUp[bin-1]+=(shapesUp[h][i].GetBinContent(bin)-histos[h].GetBinContent(bin))**2
127			errDown[bin-1]+=(shapesDown[h][i].GetBinContent(bin)-histos[h].GetBinContent(bin))**2
128			errUp[bin-1]+=(histos[h].GetBinError(bin))**2
129			errDown[bin-1]+=(histos[h].GetBinError(bin))**2
130
131			Error.SetPoint(bin-1,point,1)
132			counter += 1
133
134			for bin in range(0,len(total)):
135			if not total[bin] == 0:
136			Error.SetPointEYlow(bin,sqrt(errDown[bin])/total[bin])
137			Error.SetPointEYhigh(bin,sqrt(errUp[bin])/total[bin])
138
139			#if counter == 0: Errors=Error
140			#else: Errors.Add(Error)
141			#set the errors
142
143			#get the Ratios
144
145
146
147
148			datas=[input.Get('data_obs')]
149			datatyps = [None]
150			datanames=[d]
151
152	peller	1.33
153			histos.append(copy(Overlay))
154			typs.append('ZH')
155	peller	1.32	#histos.append(copy(Overlay))
156			#typs.append('ZH')
157
158			Stack.histos = histos
159			Stack.typs = typs
160			Stack.datas = datas
161			Stack.datatyps = datatyps
162			Stack.datanames= datanames
163			Stack.overlay = Overlay
164			Stack.AddErrors=Error
165			Stack.lumi = lumi
166			Stack.doPlot()
167
168			print 'i am done!\n'
169			#-------------------------------------------------
170	peller	1.1
171
172	peller	1.25	else:
173	peller	1.32	#----------Histo from trees------------
174			vars = (config.get(section, 'vars')).split(',')
175	peller	1.1
176	peller	1.32	if 'ZLight' in region or 'TTbar' in region or 'Zbb' in region: SignalRegion = False
177			else:
178			SignalRegion = True
179			print 'You are in the Signal Region!'
180
181			data = config.get(section,'Datas')
182
183			samples=config.get('Plot_general','samples')
184			samples=samples.split(',')
185
186			weightF=config.get('Weights','weightF')
187			Group = eval(config.get('Plot_general','Group'))
188
189			#GETALL AT ONCE
190			options = []
191			Stacks = []
192			for i in range(len(vars)):
193			Stacks.append(StackMaker(config,vars[i],region,SignalRegion))
194			options.append(Stacks[i].options)
195
196			Plotter=HistoMaker(path,config,region,options)
197
198			#print '\nProducing Plot of %s\n'%vars[v]
199			Lhistos = [[] for _ in range(0,len(vars))]
200			Ltyps = [[] for _ in range(0,len(vars))]
201			Ldatas = [[] for _ in range(0,len(vars))]
202			Ldatatyps = [[] for _ in range(0,len(vars))]
203			Ldatanames = [[] for _ in range(0,len(vars))]
204
205			#Find out Lumi:
206			lumicounter=0.
207			lumi=0.
208			for job in info:
209			if job.name in data:
210			lumi+=float(job.lumi)
211			lumicounter+=1.
212
213			if lumicounter > 0:
214			lumi=lumi/lumicounter
215
216			Plotter.lumi=lumi
217			mass = Stacks[0].mass
218
219			for job in info:
220			if eval(job.active):
221			if job.subsamples:
222			for subsample in range(0,len(job.subnames)):
223
224			if job.subnames[subsample] in samples:
225			hTempList, typList = Plotter.getHistoFromTree(job,subsample)
226			for v in range(0,len(vars)):
227			Lhistos[v].append(hTempList[v])
228			Ltyps[v].append(Group[job.subnames[subsample]])
229			print job.subnames[subsample]
230	peller	1.12
231	peller	1.32	else:
232			if job.name in samples:
233			if job.name == mass:
234			print job.name
235			hTempList, typList = Plotter.getHistoFromTree(job)
236			for v in range(0,len(vars)):
237			if SignalRegion:
238			Lhistos[v].append(hTempList[v])
239			Ltyps[v].append(Group[job.name])
240			Overlaylist= deepcopy(hTempList)
241
242			else:
243			print job.name
244			hTempList, typList = Plotter.getHistoFromTree(job)
245			for v in range(0,len(vars)):
246	peller	1.25	Lhistos[v].append(hTempList[v])
247			Ltyps[v].append(Group[job.name])
248	peller	1.32
249			elif job.name in data:
250			#print 'DATA'
251			hTemp, typ = Plotter.getHistoFromTree(job)
252	peller	1.25	for v in range(0,len(vars)):
253	peller	1.32	Ldatas[v].append(hTemp[v])
254			Ldatatyps[v].append(typ[v])
255			Ldatanames[v].append(job.name)
256
257			for v in range(0,len(vars)):
258
259			histos = Lhistos[v]
260			typs = Ltyps[v]
261			Stacks[v].histos = Lhistos[v]
262			Stacks[v].typs = Ltyps[v]
263			Stacks[v].datas = Ldatas[v]
264			Stacks[v].datatyps = Ldatatyps[v]
265			Stacks[v].datanames= Ldatanames[v]
266			Stacks[v].overlay = Overlaylist[v]
267			Stacks[v].lumi = lumi
268			Stacks[v].doPlot()
269			print 'i am done!\n'
270			#----------------------------------------------------
271	peller	1.5	sys.exit(0)