VHbb/python/HistoMaker.py

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

class HistoMaker:
    def __init__(self, path, config, region, optionsList,rescale=1,which_weightF='weightF'):
        self.path = path
        self.config = config
        self.optionsList = optionsList
        self.rescale = rescale
        self.which_weightF=which_weightF
        self.region = region
        self.lumi=0.

    def getScale(self,job,subsample=-1):
        anaTag=self.config.get('Analysis','tag')
        input = TFile.Open(self.path+'/'+job.getpath())
        CountWithPU = input.Get("CountWithPU")
        CountWithPU2011B = input.Get("CountWithPU2011B")
        #print lumi*xsecs[i]/hist.GetBinContent(1)
        if subsample>-1:
            xsec=float(job.xsec[subsample])
            sf=float(job.sf[subsample])
        else:
            xsec=float(job.xsec)
            sf=float(job.sf)
        theScale = 1.
        if anaTag == '7TeV':
            theScale = float(self.lumi)*xsec*sf/(0.46502*CountWithPU.GetBinContent(1)+0.53498*CountWithPU2011B.GetBinContent(1))*self.rescale/float(job.split)
        elif anaTag == '8TeV':
            theScale = float(self.lumi)*xsec*sf/(CountWithPU.GetBinContent(1))*self.rescale/float(job.split)
        return theScale 


    def getHistoFromTree(self,job,subsample=-1):
        if self.lumi == 0: raise Exception("You're trying to plot with no lumi")
         
        hTreeList=[]
        groupList=[]

        #get the conversion rate in case of BDT plots
        TrainFlag = eval(self.config.get('Analysis','TrainFlag'))
        if TrainFlag:
            MC_rescale_factor=2.
            print 'I RESCALE BY 2.0'
        else: MC_rescale_factor = 1.

        BDT_add_cut='EventForTraining == 0'


        plot_path = self.config.get('Directories','plotpath')
        addOverFlow=eval(self.config.get('Plot_general','addOverFlow'))

        # define treeCut
        if job.type != 'DATA':
            if type(self.region)==str:
                cutcut=self.config.get('Cuts',self.region)
            elif type(self.region)==list:
                #replace vars with other vars in the cutstring (used in DC writer)
                cutcut=self.config.get('Cuts',self.region[0])
                cutcut=cutcut.replace(self.region[1],self.region[2])
                #print cutcut
            if subsample>-1:
                treeCut='%s & %s'%(cutcut,job.subcuts[subsample])        
            else:
                treeCut='%s'%(cutcut)
        elif job.type == 'DATA':
            cutcut=self.config.get('Cuts',self.region)
            treeCut='%s'%(cutcut)

        # get and skim the Trees
        output=TFile.Open(plot_path+'/tmp_plotCache_%s_%s.root'%(self.region,job.identifier),'recreate')
        input = TFile.Open(self.path+'/'+job.getpath(),'read')
        Tree = input.Get(job.tree)
        output.cd()
        CuttedTree=Tree.CopyTree(treeCut)
    
        # get all Histos at once
        weightF=self.config.get('Weights',self.which_weightF)
        for options in self.optionsList:
            if subsample>-1:
                name=job.subnames[subsample]
                group=job.group[subsample]
            else:
                name=job.name
                group=job.group
            treeVar=options[0]
            name=options[1]
            nBins=int(options[3])
            xMin=float(options[4])
            xMax=float(options[5])

            #options

            if job.type != 'DATA':
                if CuttedTree.GetEntries():
                    
                    if 'BDT' in treeVar: drawoption = '(%s)*(%s)'%(weightF,BDT_add_cut)
                    else: drawoption = '%s'%(weightF)
                    output.cd() 
                    CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax), drawoption, "goff,e")
                    full=True
                else:
                    full=False
            elif job.type == 'DATA':
                if options[11] == 'blind':
                    output.cd()
                    if treeVar == 'H.mass':
                        CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),treeVar+'<80. || '+treeVar + '>150.' , "goff,e")
                    else:
                        CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),treeVar+'<0', "goff,e")

                else:
                    output.cd()
                    CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),'1', "goff,e")
                full = True
            if full:
                hTree = ROOT.gDirectory.Get(name)
            else:
                output.cd()
                hTree = ROOT.TH1F('%s'%name,'%s'%name,nBins,xMin,xMax)
                hTree.Sumw2()
            if job.type != 'DATA':
                if 'BDT' in treeVar: ScaleFactor = self.getScale(job,subsample,MC_rescale_factor)
                else: ScaleFactor = self.getScale(job,subsample)
                if ScaleFactor != 0:
                    hTree.Scale(ScaleFactor)
            #print '\t-->import %s\t Integral: %s'%(job.name,hTree.Integral())
            if addOverFlow:
                uFlow = hTree.GetBinContent(0)+hTree.GetBinContent(1)
                oFlow = hTree.GetBinContent(hTree.GetNbinsX()+1)+hTree.GetBinContent(hTree.GetNbinsX())
                uFlowErr = ROOT.TMath.Sqrt(ROOT.TMath.Power(hTree.GetBinError(0),2)+ROOT.TMath.Power(hTree.GetBinError(1),2))
                oFlowErr = ROOT.TMath.Sqrt(ROOT.TMath.Power(hTree.GetBinError(hTree.GetNbinsX()),2)+ROOT.TMath.Power(hTree.GetBinError(hTree.GetNbinsX()+1),2))
                hTree.SetBinContent(1,uFlow)
                hTree.SetBinContent(hTree.GetNbinsX(),oFlow)
                hTree.SetBinError(1,uFlowErr)
                hTree.SetBinError(hTree.GetNbinsX(),oFlowErr)
            hTree.SetDirectory(0)
            input.Close()
            hTreeList.append(hTree)
            groupList.append(group)
            
        return hTreeList, groupList
        

######################
def orderandadd(histos,typs,setup):
#ORDER AND ADD TOGETHER
    ordnung=[]
    ordnungtyp=[]
    num=[0]*len(setup)
    for i in range(0,len(setup)):
        for j in range(0,len(histos)):
            if typs[j] in setup[i]:
                num[i]+=1
                ordnung.append(histos[j])
                ordnungtyp.append(typs[j])
    del histos
    del typs
    histos=ordnung
    typs=ordnungtyp
    print typs
    for k in range(0,len(num)):
        for m in range(0,num[k]):
            if m > 0:
                #add
                histos[k].Add(histos[k+1],1)
                printc('magenta','','\t--> added %s to %s'%(typs[k],typs[k+1]))
                del histos[k+1]
                del typs[k+1]
    del histos[len(setup):]
    del typs[len(setup):]
    return histos, typs
Revision:	1.6
Committed:	Thu Oct 11 08:39:53 2012 UTC (12 years, 7 months ago) by peller
Content type:	text/x-python
Branch:	MAIN
Changes since 1.5:	+5 -1 lines
Log Message:	plotting
#	User	Rev	Content
1	peller	1.1	from samplesclass import sample
2			from printcolor import printc
3			import pickle
4			import ROOT
5			from ROOT import TFile, TTree
6			import ROOT
7			from array import array
8			from BetterConfigParser import BetterConfigParser
9			import sys
10
11			class HistoMaker:
12	peller	1.2	def __init__(self, path, config, region, optionsList,rescale=1,which_weightF='weightF'):
13	peller	1.1	self.path = path
14			self.config = config
15			self.optionsList = optionsList
16			self.rescale = rescale
17			self.which_weightF=which_weightF
18	peller	1.2	self.region = region
19	peller	1.3	self.lumi=0.
20	peller	1.1
21			def getScale(self,job,subsample=-1):
22			anaTag=self.config.get('Analysis','tag')
23			input = TFile.Open(self.path+'/'+job.getpath())
24			CountWithPU = input.Get("CountWithPU")
25			CountWithPU2011B = input.Get("CountWithPU2011B")
26			#print lumi*xsecs[i]/hist.GetBinContent(1)
27			if subsample>-1:
28			xsec=float(job.xsec[subsample])
29			sf=float(job.sf[subsample])
30			else:
31			xsec=float(job.xsec)
32			sf=float(job.sf)
33			theScale = 1.
34			if anaTag == '7TeV':
35	peller	1.3	theScale = float(self.lumi)xsecsf/(0.46502CountWithPU.GetBinContent(1)+0.53498CountWithPU2011B.GetBinContent(1))*self.rescale/float(job.split)
36	peller	1.1	elif anaTag == '8TeV':
37	peller	1.3	theScale = float(self.lumi)xsecsf/(CountWithPU.GetBinContent(1))*self.rescale/float(job.split)
38	peller	1.1	return theScale
39
40
41			def getHistoFromTree(self,job,subsample=-1):
42	peller	1.3	if self.lumi == 0: raise Exception("You're trying to plot with no lumi")
43
44	peller	1.1	hTreeList=[]
45			groupList=[]
46
47	peller	1.5	#get the conversion rate in case of BDT plots
48			TrainFlag = eval(self.config.get('Analysis','TrainFlag'))
49			if TrainFlag:
50			MC_rescale_factor=2.
51			print 'I RESCALE BY 2.0'
52			else: MC_rescale_factor = 1.
53
54			BDT_add_cut='EventForTraining == 0'
55
56	peller	1.2
57			plot_path = self.config.get('Directories','plotpath')
58	nmohr	1.4	addOverFlow=eval(self.config.get('Plot_general','addOverFlow'))
59	peller	1.2
60			# define treeCut
61	peller	1.1	if job.type != 'DATA':
62	peller	1.2	if type(self.region)==str:
63			cutcut=self.config.get('Cuts',self.region)
64			elif type(self.region)==list:
65			#replace vars with other vars in the cutstring (used in DC writer)
66			cutcut=self.config.get('Cuts',self.region[0])
67			cutcut=cutcut.replace(self.region[1],self.region[2])
68			#print cutcut
69	peller	1.1	if subsample>-1:
70	peller	1.5	treeCut='%s & %s'%(cutcut,job.subcuts[subsample])
71	peller	1.1	else:
72	peller	1.5	treeCut='%s'%(cutcut)
73	peller	1.1	elif job.type == 'DATA':
74	peller	1.2	cutcut=self.config.get('Cuts',self.region)
75	peller	1.1	treeCut='%s'%(cutcut)
76	peller	1.2
77			# get and skim the Trees
78			output=TFile.Open(plot_path+'/tmp_plotCache_%s_%s.root'%(self.region,job.identifier),'recreate')
79	peller	1.1	input = TFile.Open(self.path+'/'+job.getpath(),'read')
80			Tree = input.Get(job.tree)
81	peller	1.2	output.cd()
82			CuttedTree=Tree.CopyTree(treeCut)
83
84			# get all Histos at once
85	peller	1.1	weightF=self.config.get('Weights',self.which_weightF)
86			for options in self.optionsList:
87			if subsample>-1:
88			name=job.subnames[subsample]
89			group=job.group[subsample]
90			else:
91			name=job.name
92			group=job.group
93			treeVar=options[0]
94			name=options[1]
95			nBins=int(options[3])
96			xMin=float(options[4])
97			xMax=float(options[5])
98
99	peller	1.5	#options
100
101	peller	1.1	if job.type != 'DATA':
102			if CuttedTree.GetEntries():
103	peller	1.5
104			if 'BDT' in treeVar: drawoption = '(%s)*(%s)'%(weightF,BDT_add_cut)
105			else: drawoption = '%s'%(weightF)
106	peller	1.1	output.cd()
107	peller	1.5	CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax), drawoption, "goff,e")
108	peller	1.1	full=True
109			else:
110			full=False
111			elif job.type == 'DATA':
112			if options[11] == 'blind':
113			output.cd()
114	peller	1.6	if treeVar == 'H.mass':
115			CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),treeVar+'<80. \|\| '+treeVar + '>150.' , "goff,e")
116			else:
117			CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),treeVar+'<0', "goff,e")
118
119	peller	1.1	else:
120			output.cd()
121			CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),'1', "goff,e")
122			full = True
123			if full:
124			hTree = ROOT.gDirectory.Get(name)
125			else:
126			output.cd()
127			hTree = ROOT.TH1F('%s'%name,'%s'%name,nBins,xMin,xMax)
128			hTree.Sumw2()
129			if job.type != 'DATA':
130	peller	1.5	if 'BDT' in treeVar: ScaleFactor = self.getScale(job,subsample,MC_rescale_factor)
131			else: ScaleFactor = self.getScale(job,subsample)
132	peller	1.1	if ScaleFactor != 0:
133			hTree.Scale(ScaleFactor)
134			#print '\t-->import %s\t Integral: %s'%(job.name,hTree.Integral())
135	nmohr	1.4	if addOverFlow:
136			uFlow = hTree.GetBinContent(0)+hTree.GetBinContent(1)
137			oFlow = hTree.GetBinContent(hTree.GetNbinsX()+1)+hTree.GetBinContent(hTree.GetNbinsX())
138			uFlowErr = ROOT.TMath.Sqrt(ROOT.TMath.Power(hTree.GetBinError(0),2)+ROOT.TMath.Power(hTree.GetBinError(1),2))
139			oFlowErr = ROOT.TMath.Sqrt(ROOT.TMath.Power(hTree.GetBinError(hTree.GetNbinsX()),2)+ROOT.TMath.Power(hTree.GetBinError(hTree.GetNbinsX()+1),2))
140			hTree.SetBinContent(1,uFlow)
141			hTree.SetBinContent(hTree.GetNbinsX(),oFlow)
142			hTree.SetBinError(1,uFlowErr)
143			hTree.SetBinError(hTree.GetNbinsX(),oFlowErr)
144	peller	1.1	hTree.SetDirectory(0)
145			input.Close()
146			hTreeList.append(hTree)
147			groupList.append(group)
148
149			return hTreeList, groupList
150
151
152			######################
153			def orderandadd(histos,typs,setup):
154			#ORDER AND ADD TOGETHER
155			ordnung=[]
156			ordnungtyp=[]
157			num=[0]*len(setup)
158			for i in range(0,len(setup)):
159			for j in range(0,len(histos)):
160			if typs[j] in setup[i]:
161			num[i]+=1
162			ordnung.append(histos[j])
163			ordnungtyp.append(typs[j])
164			del histos
165			del typs
166			histos=ordnung
167			typs=ordnungtyp
168			print typs
169			for k in range(0,len(num)):
170			for m in range(0,num[k]):
171			if m > 0:
172			#add
173			histos[k].Add(histos[k+1],1)
174			printc('magenta','','\t--> added %s to %s'%(typs[k],typs[k+1]))
175			del histos[k+1]
176			del typs[k+1]
177			del histos[len(setup):]
178			del typs[len(setup):]
179			return histos, typs