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,os

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,MC_rescale_factor=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))*MC_rescale_factor/float(job.split)
        elif anaTag == '8TeV':
            theScale = float(self.lumi)*xsec*sf/(CountWithPU.GetBinContent(1))*MC_rescale_factor/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'))
        BDT_add_cut='EventForTraining == 0'


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

        scratchDir = os.environ["TMPDIR"]
        #scratchDir = '/shome/peller/'
        # 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(scratchDir+'/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)
        input.Close()
        del input
        # 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 'RTight' in treeVar or 'RMed' 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+'<90. || '+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 'RTight' in treeVar or 'RMed' in treeVar:
                    if TrainFlag:
                        MC_rescale_factor=2.
                        print 'I RESCALE BY 2.0'
                    else: MC_rescale_factor = 1.
                    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)
            hTreeList.append(hTree)
            groupList.append(group)
         
        output.Close()
        del output
        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.10
Committed:	Tue Oct 23 16:39:55 2012 UTC (12 years, 6 months ago) by peller
Content type:	text/x-python
Branch:	MAIN
CVS Tags:	hcpApproval, HCP_unblinding
Changes since 1.9:	+5 -3 lines
Log Message:	memory leak fix
#	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	peller	1.7	import sys,os
10	peller	1.1
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	peller	1.8	def getScale(self,job,subsample=-1,MC_rescale_factor=1):
22	peller	1.1	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.8	theScale = float(self.lumi)xsecsf/(0.46502CountWithPU.GetBinContent(1)+0.53498CountWithPU2011B.GetBinContent(1))*MC_rescale_factor/float(job.split)
36	peller	1.1	elif anaTag == '8TeV':
37	peller	1.8	theScale = float(self.lumi)xsecsf/(CountWithPU.GetBinContent(1))*MC_rescale_factor/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			BDT_add_cut='EventForTraining == 0'
50
51	peller	1.2
52			plot_path = self.config.get('Directories','plotpath')
53	nmohr	1.4	addOverFlow=eval(self.config.get('Plot_general','addOverFlow'))
54	peller	1.2
55	peller	1.7	scratchDir = os.environ["TMPDIR"]
56			#scratchDir = '/shome/peller/'
57	peller	1.2	# define treeCut
58	peller	1.1	if job.type != 'DATA':
59	peller	1.2	if type(self.region)==str:
60			cutcut=self.config.get('Cuts',self.region)
61			elif type(self.region)==list:
62			#replace vars with other vars in the cutstring (used in DC writer)
63			cutcut=self.config.get('Cuts',self.region[0])
64			cutcut=cutcut.replace(self.region[1],self.region[2])
65			#print cutcut
66	peller	1.1	if subsample>-1:
67	peller	1.5	treeCut='%s & %s'%(cutcut,job.subcuts[subsample])
68	peller	1.1	else:
69	peller	1.5	treeCut='%s'%(cutcut)
70	peller	1.1	elif job.type == 'DATA':
71	peller	1.2	cutcut=self.config.get('Cuts',self.region)
72	peller	1.1	treeCut='%s'%(cutcut)
73	peller	1.2
74			# get and skim the Trees
75	peller	1.7	output=TFile.Open(scratchDir+'/tmp_plotCache_%s_%s.root'%(self.region,job.identifier),'recreate')
76	peller	1.1	input = TFile.Open(self.path+'/'+job.getpath(),'read')
77			Tree = input.Get(job.tree)
78	peller	1.2	output.cd()
79			CuttedTree=Tree.CopyTree(treeCut)
80	peller	1.10	input.Close()
81			del input
82	peller	1.2	# get all Histos at once
83	peller	1.1	weightF=self.config.get('Weights',self.which_weightF)
84			for options in self.optionsList:
85			if subsample>-1:
86			name=job.subnames[subsample]
87			group=job.group[subsample]
88			else:
89			name=job.name
90			group=job.group
91			treeVar=options[0]
92			name=options[1]
93			nBins=int(options[3])
94			xMin=float(options[4])
95			xMax=float(options[5])
96
97	peller	1.5	#options
98
99	peller	1.1	if job.type != 'DATA':
100			if CuttedTree.GetEntries():
101	peller	1.5
102	peller	1.8	if 'RTight' in treeVar or 'RMed' in treeVar: drawoption = '(%s)*(%s)'%(weightF,BDT_add_cut)
103	peller	1.5	else: drawoption = '%s'%(weightF)
104	peller	1.1	output.cd()
105	peller	1.5	CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax), drawoption, "goff,e")
106	peller	1.1	full=True
107			else:
108			full=False
109			elif job.type == 'DATA':
110			if options[11] == 'blind':
111			output.cd()
112	peller	1.6	if treeVar == 'H.mass':
113	nmohr	1.9	CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),treeVar+'<90. \|\| '+treeVar + '>150.' , "goff,e")
114	peller	1.6	else:
115			CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),treeVar+'<0', "goff,e")
116
117	peller	1.1	else:
118			output.cd()
119			CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),'1', "goff,e")
120			full = True
121			if full:
122			hTree = ROOT.gDirectory.Get(name)
123			else:
124			output.cd()
125			hTree = ROOT.TH1F('%s'%name,'%s'%name,nBins,xMin,xMax)
126			hTree.Sumw2()
127			if job.type != 'DATA':
128	peller	1.8	if 'RTight' in treeVar or 'RMed' in treeVar:
129			if TrainFlag:
130			MC_rescale_factor=2.
131			print 'I RESCALE BY 2.0'
132			else: MC_rescale_factor = 1.
133			ScaleFactor = self.getScale(job,subsample,MC_rescale_factor)
134	peller	1.5	else: ScaleFactor = self.getScale(job,subsample)
135	peller	1.1	if ScaleFactor != 0:
136			hTree.Scale(ScaleFactor)
137			#print '\t-->import %s\t Integral: %s'%(job.name,hTree.Integral())
138	nmohr	1.4	if addOverFlow:
139			uFlow = hTree.GetBinContent(0)+hTree.GetBinContent(1)
140			oFlow = hTree.GetBinContent(hTree.GetNbinsX()+1)+hTree.GetBinContent(hTree.GetNbinsX())
141			uFlowErr = ROOT.TMath.Sqrt(ROOT.TMath.Power(hTree.GetBinError(0),2)+ROOT.TMath.Power(hTree.GetBinError(1),2))
142			oFlowErr = ROOT.TMath.Sqrt(ROOT.TMath.Power(hTree.GetBinError(hTree.GetNbinsX()),2)+ROOT.TMath.Power(hTree.GetBinError(hTree.GetNbinsX()+1),2))
143			hTree.SetBinContent(1,uFlow)
144			hTree.SetBinContent(hTree.GetNbinsX(),oFlow)
145			hTree.SetBinError(1,uFlowErr)
146			hTree.SetBinError(hTree.GetNbinsX(),oFlowErr)
147	peller	1.1	hTree.SetDirectory(0)
148			hTreeList.append(hTree)
149			groupList.append(group)
150	peller	1.10
151			output.Close()
152			del output
153	peller	1.1	return hTreeList, groupList
154
155
156			######################
157			def orderandadd(histos,typs,setup):
158			#ORDER AND ADD TOGETHER
159			ordnung=[]
160			ordnungtyp=[]
161			num=[0]*len(setup)
162			for i in range(0,len(setup)):
163			for j in range(0,len(histos)):
164			if typs[j] in setup[i]:
165			num[i]+=1
166			ordnung.append(histos[j])
167			ordnungtyp.append(typs[j])
168			del histos
169			del typs
170			histos=ordnung
171			typs=ordnungtyp
172			print typs
173			for k in range(0,len(num)):
174			for m in range(0,num[k]):
175			if m > 0:
176			#add
177			histos[k].Add(histos[k+1],1)
178			printc('magenta','','\t--> added %s to %s'%(typs[k],typs[k+1]))
179			del histos[k+1]
180			del typs[k+1]
181			del histos[len(setup):]
182			del typs[len(setup):]
183			return histos, typs