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()
                    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.5
Committed:	Fri Oct 5 11:58:18 2012 UTC (12 years, 7 months ago) by peller
Content type:	text/x-python
Branch:	MAIN
CVS Tags:	hcpPreAppFreeze
Changes since 1.4:	+19 -4 lines
Log Message:	plot and other fixes
#	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			CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),treeVar+'<0', "goff,e")
115			else:
116			output.cd()
117			CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),'1', "goff,e")
118			full = True
119			if full:
120			hTree = ROOT.gDirectory.Get(name)
121			else:
122			output.cd()
123			hTree = ROOT.TH1F('%s'%name,'%s'%name,nBins,xMin,xMax)
124			hTree.Sumw2()
125			if job.type != 'DATA':
126	peller	1.5	if 'BDT' in treeVar: ScaleFactor = self.getScale(job,subsample,MC_rescale_factor)
127			else: ScaleFactor = self.getScale(job,subsample)
128	peller	1.1	if ScaleFactor != 0:
129			hTree.Scale(ScaleFactor)
130			#print '\t-->import %s\t Integral: %s'%(job.name,hTree.Integral())
131	nmohr	1.4	if addOverFlow:
132			uFlow = hTree.GetBinContent(0)+hTree.GetBinContent(1)
133			oFlow = hTree.GetBinContent(hTree.GetNbinsX()+1)+hTree.GetBinContent(hTree.GetNbinsX())
134			uFlowErr = ROOT.TMath.Sqrt(ROOT.TMath.Power(hTree.GetBinError(0),2)+ROOT.TMath.Power(hTree.GetBinError(1),2))
135			oFlowErr = ROOT.TMath.Sqrt(ROOT.TMath.Power(hTree.GetBinError(hTree.GetNbinsX()),2)+ROOT.TMath.Power(hTree.GetBinError(hTree.GetNbinsX()+1),2))
136			hTree.SetBinContent(1,uFlow)
137			hTree.SetBinContent(hTree.GetNbinsX(),oFlow)
138			hTree.SetBinError(1,uFlowErr)
139			hTree.SetBinError(hTree.GetNbinsX(),oFlowErr)
140	peller	1.1	hTree.SetDirectory(0)
141			input.Close()
142			hTreeList.append(hTree)
143			groupList.append(group)
144
145			return hTreeList, groupList
146
147
148			######################
149			def orderandadd(histos,typs,setup):
150			#ORDER AND ADD TOGETHER
151			ordnung=[]
152			ordnungtyp=[]
153			num=[0]*len(setup)
154			for i in range(0,len(setup)):
155			for j in range(0,len(histos)):
156			if typs[j] in setup[i]:
157			num[i]+=1
158			ordnung.append(histos[j])
159			ordnungtyp.append(typs[j])
160			del histos
161			del typs
162			histos=ordnung
163			typs=ordnungtyp
164			print typs
165			for k in range(0,len(num)):
166			for m in range(0,num[k]):
167			if m > 0:
168			#add
169			histos[k].Add(histos[k+1],1)
170			printc('magenta','','\t--> added %s to %s'%(typs[k],typs[k+1]))
171			del histos[k+1]
172			del typs[k+1]
173			del histos[len(setup):]
174			del typs[len(setup):]
175			return histos, typs