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)
    
        # 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+'<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 '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)
            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.8
Committed:	Fri Oct 12 09:09:01 2012 UTC (12 years, 7 months ago) by peller
Content type:	text/x-python
Branch:	MAIN
CVS Tags:	hcpPreApp
Changes since 1.7:	+10 -10 lines
Log Message:	BDT fix in HistoMaker
#	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
81			# get all Histos at once
82	peller	1.1	weightF=self.config.get('Weights',self.which_weightF)
83			for options in self.optionsList:
84			if subsample>-1:
85			name=job.subnames[subsample]
86			group=job.group[subsample]
87			else:
88			name=job.name
89			group=job.group
90			treeVar=options[0]
91			name=options[1]
92			nBins=int(options[3])
93			xMin=float(options[4])
94			xMax=float(options[5])
95
96	peller	1.5	#options
97
98	peller	1.1	if job.type != 'DATA':
99			if CuttedTree.GetEntries():
100	peller	1.5
101	peller	1.8	if 'RTight' in treeVar or 'RMed' in treeVar: drawoption = '(%s)*(%s)'%(weightF,BDT_add_cut)
102	peller	1.5	else: drawoption = '%s'%(weightF)
103	peller	1.1	output.cd()
104	peller	1.5	CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax), drawoption, "goff,e")
105	peller	1.1	full=True
106			else:
107			full=False
108			elif job.type == 'DATA':
109			if options[11] == 'blind':
110			output.cd()
111	peller	1.6	if treeVar == 'H.mass':
112			CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),treeVar+'<80. \|\| '+treeVar + '>150.' , "goff,e")
113			else:
114			CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),treeVar+'<0', "goff,e")
115
116	peller	1.1	else:
117			output.cd()
118			CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),'1', "goff,e")
119			full = True
120			if full:
121			hTree = ROOT.gDirectory.Get(name)
122			else:
123			output.cd()
124			hTree = ROOT.TH1F('%s'%name,'%s'%name,nBins,xMin,xMax)
125			hTree.Sumw2()
126			if job.type != 'DATA':
127	peller	1.8	if 'RTight' in treeVar or 'RMed' in treeVar:
128			if TrainFlag:
129			MC_rescale_factor=2.
130			print 'I RESCALE BY 2.0'
131			else: MC_rescale_factor = 1.
132			ScaleFactor = self.getScale(job,subsample,MC_rescale_factor)
133	peller	1.5	else: ScaleFactor = self.getScale(job,subsample)
134	peller	1.1	if ScaleFactor != 0:
135			hTree.Scale(ScaleFactor)
136			#print '\t-->import %s\t Integral: %s'%(job.name,hTree.Integral())
137	nmohr	1.4	if addOverFlow:
138			uFlow = hTree.GetBinContent(0)+hTree.GetBinContent(1)
139			oFlow = hTree.GetBinContent(hTree.GetNbinsX()+1)+hTree.GetBinContent(hTree.GetNbinsX())
140			uFlowErr = ROOT.TMath.Sqrt(ROOT.TMath.Power(hTree.GetBinError(0),2)+ROOT.TMath.Power(hTree.GetBinError(1),2))
141			oFlowErr = ROOT.TMath.Sqrt(ROOT.TMath.Power(hTree.GetBinError(hTree.GetNbinsX()),2)+ROOT.TMath.Power(hTree.GetBinError(hTree.GetNbinsX()+1),2))
142			hTree.SetBinContent(1,uFlow)
143			hTree.SetBinContent(hTree.GetNbinsX(),oFlow)
144			hTree.SetBinError(1,uFlowErr)
145			hTree.SetBinError(hTree.GetNbinsX(),oFlowErr)
146	peller	1.1	hTree.SetDirectory(0)
147			input.Close()
148			hTreeList.append(hTree)
149			groupList.append(group)
150
151			return hTreeList, groupList
152
153
154			######################
155			def orderandadd(histos,typs,setup):
156			#ORDER AND ADD TOGETHER
157			ordnung=[]
158			ordnungtyp=[]
159			num=[0]*len(setup)
160			for i in range(0,len(setup)):
161			for j in range(0,len(histos)):
162			if typs[j] in setup[i]:
163			num[i]+=1
164			ordnung.append(histos[j])
165			ordnungtyp.append(typs[j])
166			del histos
167			del typs
168			histos=ordnung
169			typs=ordnungtyp
170			print typs
171			for k in range(0,len(num)):
172			for m in range(0,num[k]):
173			if m > 0:
174			#add
175			histos[k].Add(histos[k+1],1)
176			printc('magenta','','\t--> added %s to %s'%(typs[k],typs[k+1]))
177			del histos[k+1]
178			del typs[k+1]
179			del histos[len(setup):]
180			del typs[len(setup):]
181			return histos, typs