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=[]


        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 & EventForTraining == 0'%(cutcut,job.subcuts[subsample])        
            else:
                treeCut='%s & EventForTraining == 0'%(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])

            if job.type != 'DATA':
                if CuttedTree.GetEntries():
                    output.cd() 
                    CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax), weightF, "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':
                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.4
Committed:	Thu Oct 4 13:02:23 2012 UTC (12 years, 7 months ago) by nmohr
Content type:	text/x-python
Branch:	MAIN
Changes since 1.3:	+10 -0 lines
Log Message:	Plotting style
#	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.2
48			plot_path = self.config.get('Directories','plotpath')
49	nmohr	1.4	addOverFlow=eval(self.config.get('Plot_general','addOverFlow'))
50	peller	1.2
51			# define treeCut
52	peller	1.1	if job.type != 'DATA':
53	peller	1.2	if type(self.region)==str:
54			cutcut=self.config.get('Cuts',self.region)
55			elif type(self.region)==list:
56			#replace vars with other vars in the cutstring (used in DC writer)
57			cutcut=self.config.get('Cuts',self.region[0])
58			cutcut=cutcut.replace(self.region[1],self.region[2])
59			#print cutcut
60	peller	1.1	if subsample>-1:
61			treeCut='%s & %s & EventForTraining == 0'%(cutcut,job.subcuts[subsample])
62			else:
63			treeCut='%s & EventForTraining == 0'%(cutcut)
64			elif job.type == 'DATA':
65	peller	1.2	cutcut=self.config.get('Cuts',self.region)
66	peller	1.1	treeCut='%s'%(cutcut)
67	peller	1.2
68			# get and skim the Trees
69			output=TFile.Open(plot_path+'/tmp_plotCache_%s_%s.root'%(self.region,job.identifier),'recreate')
70	peller	1.1	input = TFile.Open(self.path+'/'+job.getpath(),'read')
71			Tree = input.Get(job.tree)
72	peller	1.2	output.cd()
73			CuttedTree=Tree.CopyTree(treeCut)
74
75			# get all Histos at once
76	peller	1.1	weightF=self.config.get('Weights',self.which_weightF)
77			for options in self.optionsList:
78			if subsample>-1:
79			name=job.subnames[subsample]
80			group=job.group[subsample]
81			else:
82			name=job.name
83			group=job.group
84			treeVar=options[0]
85			name=options[1]
86			nBins=int(options[3])
87			xMin=float(options[4])
88			xMax=float(options[5])
89
90			if job.type != 'DATA':
91			if CuttedTree.GetEntries():
92			output.cd()
93			CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax), weightF, "goff,e")
94			full=True
95			else:
96			full=False
97			elif job.type == 'DATA':
98			if options[11] == 'blind':
99			output.cd()
100			CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),treeVar+'<0', "goff,e")
101			else:
102			output.cd()
103			CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),'1', "goff,e")
104			full = True
105			if full:
106			hTree = ROOT.gDirectory.Get(name)
107			else:
108			output.cd()
109			hTree = ROOT.TH1F('%s'%name,'%s'%name,nBins,xMin,xMax)
110			hTree.Sumw2()
111			if job.type != 'DATA':
112			ScaleFactor = self.getScale(job,subsample)
113			if ScaleFactor != 0:
114			hTree.Scale(ScaleFactor)
115			#print '\t-->import %s\t Integral: %s'%(job.name,hTree.Integral())
116	nmohr	1.4	if addOverFlow:
117			uFlow = hTree.GetBinContent(0)+hTree.GetBinContent(1)
118			oFlow = hTree.GetBinContent(hTree.GetNbinsX()+1)+hTree.GetBinContent(hTree.GetNbinsX())
119			uFlowErr = ROOT.TMath.Sqrt(ROOT.TMath.Power(hTree.GetBinError(0),2)+ROOT.TMath.Power(hTree.GetBinError(1),2))
120			oFlowErr = ROOT.TMath.Sqrt(ROOT.TMath.Power(hTree.GetBinError(hTree.GetNbinsX()),2)+ROOT.TMath.Power(hTree.GetBinError(hTree.GetNbinsX()+1),2))
121			hTree.SetBinContent(1,uFlow)
122			hTree.SetBinContent(hTree.GetNbinsX(),oFlow)
123			hTree.SetBinError(1,uFlowErr)
124			hTree.SetBinError(hTree.GetNbinsX(),oFlowErr)
125	peller	1.1	hTree.SetDirectory(0)
126			input.Close()
127			hTreeList.append(hTree)
128			groupList.append(group)
129
130			return hTreeList, groupList
131
132
133			######################
134			def orderandadd(histos,typs,setup):
135			#ORDER AND ADD TOGETHER
136			ordnung=[]
137			ordnungtyp=[]
138			num=[0]*len(setup)
139			for i in range(0,len(setup)):
140			for j in range(0,len(histos)):
141			if typs[j] in setup[i]:
142			num[i]+=1
143			ordnung.append(histos[j])
144			ordnungtyp.append(typs[j])
145			del histos
146			del typs
147			histos=ordnung
148			typs=ordnungtyp
149			print typs
150			for k in range(0,len(num)):
151			for m in range(0,num[k]):
152			if m > 0:
153			#add
154			histos[k].Add(histos[k+1],1)
155			printc('magenta','','\t--> added %s to %s'%(typs[k],typs[k+1]))
156			del histos[k+1]
157			del typs[k+1]
158			del histos[len(setup):]
159			del typs[len(setup):]
160			return histos, typs