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):
        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'))

        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 '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.7
Committed:	Thu Oct 11 13:40:18 2012 UTC (12 years, 7 months ago) by peller
Content type:	text/x-python
Branch:	MAIN
Changes since 1.6:	+4 -2 lines
Log Message:	plotting cache
#	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			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	peller	1.7	scratchDir = os.environ["TMPDIR"]
61			#scratchDir = '/shome/peller/'
62	peller	1.2	# define treeCut
63	peller	1.1	if job.type != 'DATA':
64	peller	1.2	if type(self.region)==str:
65			cutcut=self.config.get('Cuts',self.region)
66			elif type(self.region)==list:
67			#replace vars with other vars in the cutstring (used in DC writer)
68			cutcut=self.config.get('Cuts',self.region[0])
69			cutcut=cutcut.replace(self.region[1],self.region[2])
70			#print cutcut
71	peller	1.1	if subsample>-1:
72	peller	1.5	treeCut='%s & %s'%(cutcut,job.subcuts[subsample])
73	peller	1.1	else:
74	peller	1.5	treeCut='%s'%(cutcut)
75	peller	1.1	elif job.type == 'DATA':
76	peller	1.2	cutcut=self.config.get('Cuts',self.region)
77	peller	1.1	treeCut='%s'%(cutcut)
78	peller	1.2
79			# get and skim the Trees
80	peller	1.7	output=TFile.Open(scratchDir+'/tmp_plotCache_%s_%s.root'%(self.region,job.identifier),'recreate')
81	peller	1.1	input = TFile.Open(self.path+'/'+job.getpath(),'read')
82			Tree = input.Get(job.tree)
83	peller	1.2	output.cd()
84			CuttedTree=Tree.CopyTree(treeCut)
85
86			# get all Histos at once
87	peller	1.1	weightF=self.config.get('Weights',self.which_weightF)
88			for options in self.optionsList:
89			if subsample>-1:
90			name=job.subnames[subsample]
91			group=job.group[subsample]
92			else:
93			name=job.name
94			group=job.group
95			treeVar=options[0]
96			name=options[1]
97			nBins=int(options[3])
98			xMin=float(options[4])
99			xMax=float(options[5])
100
101	peller	1.5	#options
102
103	peller	1.1	if job.type != 'DATA':
104			if CuttedTree.GetEntries():
105	peller	1.5
106			if 'BDT' in treeVar: drawoption = '(%s)*(%s)'%(weightF,BDT_add_cut)
107			else: drawoption = '%s'%(weightF)
108	peller	1.1	output.cd()
109	peller	1.5	CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax), drawoption, "goff,e")
110	peller	1.1	full=True
111			else:
112			full=False
113			elif job.type == 'DATA':
114			if options[11] == 'blind':
115			output.cd()
116	peller	1.6	if treeVar == 'H.mass':
117			CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),treeVar+'<80. \|\| '+treeVar + '>150.' , "goff,e")
118			else:
119			CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),treeVar+'<0', "goff,e")
120
121	peller	1.1	else:
122			output.cd()
123			CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),'1', "goff,e")
124			full = True
125			if full:
126			hTree = ROOT.gDirectory.Get(name)
127			else:
128			output.cd()
129			hTree = ROOT.TH1F('%s'%name,'%s'%name,nBins,xMin,xMax)
130			hTree.Sumw2()
131			if job.type != 'DATA':
132	peller	1.5	if 'BDT' in treeVar: ScaleFactor = self.getScale(job,subsample,MC_rescale_factor)
133			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