VHbb/python/gethistofromtree.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


def getScale(job,path,config,rescale,subsample=-1):
    anaTag=config.get('Analysis','tag')
    input = TFile.Open(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(job.lumi)*xsec*sf/(0.46502*CountWithPU.GetBinContent(1)+0.53498*CountWithPU2011B.GetBinContent(1))*rescale/float(job.split)
    elif anaTag == '8TeV':
        theScale = float(job.lumi)*xsec*sf/(CountWithPU.GetBinContent(1))*rescale/float(job.split)
    return theScale 

def getHistoFromTree(job,path,config,options,rescale=1,subsample=-1,which_weightF='weightF'):

    #print job.getpath()
    #print options
    treeVar=options[0]
    if subsample>-1:
        name=job.subnames[subsample]
        group=job.group[subsample]
    else:
        name=job.name
        group=job.group

    #title=job.plotname()
    nBins=int(options[3])
    xMin=float(options[4])
    xMax=float(options[5])
    #addOverFlow=eval(config.get('Plot_general','addOverFlow'))
    addOverFlow = False

    if job.type != 'DATA':
    
        if type(options[7])==str:
            cutcut=config.get('Cuts',options[7])
        elif type(options[7])==list:
            cutcut=config.get('Cuts',options[7][0])
            cutcut=cutcut.replace(options[7][1],options[7][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=config.get('Cuts',options[8])
        treeCut='%s'%(cutcut)


    input = TFile.Open(path+'/'+job.getpath(),'read')

    Tree = input.Get(job.tree)
    #Tree=tmpTree.CloneTree()
    #Tree.SetDirectory(0)
    
    #Tree=tmpTree.Clone()
    weightF=config.get('Weights',which_weightF)
    #hTree = ROOT.TH1F('%s'%name,'%s'%title,nBins,xMin,xMax)
    #hTree.SetDirectory(0)
    #hTree.Sumw2()
    #print 'drawing...'
    if job.type != 'DATA':
        #print treeCut
        #print job.name
        if Tree.GetEntries():
            Tree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),'(%s)*(%s)' %(treeCut,weightF), "goff,e")
            full=True
        else:
            full=False
    elif job.type == 'DATA':
    
        if len(options)>10:
            if options[11] == 'blind':
                treeCut = treeCut + '&'+treeVar+'<0'
    
    
        Tree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),treeCut, "goff,e")
        full = True
    if full:
        hTree = ROOT.gDirectory.Get(name)
    else:
        hTree = ROOT.TH1F('%s'%name,'%s'%name,nBins,xMin,xMax)
        hTree.Sumw2()
    #print job.name + ' Sumw2', hTree.GetEntries()

    if job.type != 'DATA':
        ScaleFactor = getScale(job,path,config,rescale,subsample)
        if ScaleFactor != 0:
            hTree.Scale(ScaleFactor)
    
    if addOverFlow:
            print 'Adding overflow'
            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)
              
            
    print '\t-->import %s\t Integral: %s'%(job.name,hTree.Integral())
            
    hTree.SetDirectory(0)
    input.Close()  
    
    return hTree, group
    

######################


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.16
Committed:	Tue Oct 9 21:17:34 2012 UTC (12 years, 7 months ago) by peller
Content type:	text/x-python
Branch:	MAIN
CVS Tags:	hcpPreAppFreeze
Changes since 1.15:	+2 -1 lines
Log Message:	only DC writer
#	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	nmohr	1.3	from BetterConfigParser import BetterConfigParser
9	peller	1.1	import sys
10
11
12	nmohr	1.10	def getScale(job,path,config,rescale,subsample=-1):
13			anaTag=config.get('Analysis','tag')
14	nmohr	1.9	input = TFile.Open(path+'/'+job.getpath())
15	peller	1.1	CountWithPU = input.Get("CountWithPU")
16			CountWithPU2011B = input.Get("CountWithPU2011B")
17			#print lumi*xsecs[i]/hist.GetBinContent(1)
18	peller	1.7
19			if subsample>-1:
20			xsec=float(job.xsec[subsample])
21			sf=float(job.sf[subsample])
22			else:
23			xsec=float(job.xsec)
24			sf=float(job.sf)
25
26
27	nmohr	1.6	theScale = 1.
28			if anaTag == '7TeV':
29	peller	1.7	theScale = float(job.lumi)xsecsf/(0.46502CountWithPU.GetBinContent(1)+0.53498CountWithPU2011B.GetBinContent(1))*rescale/float(job.split)
30	nmohr	1.6	elif anaTag == '8TeV':
31	peller	1.7	theScale = float(job.lumi)xsecsf/(CountWithPU.GetBinContent(1))*rescale/float(job.split)
32	nmohr	1.6	return theScale
33	peller	1.1
34	peller	1.11	def getHistoFromTree(job,path,config,options,rescale=1,subsample=-1,which_weightF='weightF'):
35	peller	1.8
36			#print job.getpath()
37			#print options
38	peller	1.1	treeVar=options[0]
39	peller	1.7	if subsample>-1:
40			name=job.subnames[subsample]
41			group=job.group[subsample]
42			else:
43			name=job.name
44			group=job.group
45
46	peller	1.1	#title=job.plotname()
47			nBins=int(options[3])
48			xMin=float(options[4])
49			xMax=float(options[5])
50	peller	1.16	#addOverFlow=eval(config.get('Plot_general','addOverFlow'))
51			addOverFlow = False
52	peller	1.1
53			if job.type != 'DATA':
54	peller	1.12
55			if type(options[7])==str:
56			cutcut=config.get('Cuts',options[7])
57			elif type(options[7])==list:
58			cutcut=config.get('Cuts',options[7][0])
59			cutcut=cutcut.replace(options[7][1],options[7][2])
60	peller	1.13	#print cutcut
61	peller	1.8	if subsample>-1:
62	peller	1.7	treeCut='%s & %s & EventForTraining == 0'%(cutcut,job.subcuts[subsample])
63			else:
64			treeCut='%s & EventForTraining == 0'%(cutcut)
65	peller	1.1
66			elif job.type == 'DATA':
67	peller	1.5	cutcut=config.get('Cuts',options[8])
68	peller	1.8	treeCut='%s'%(cutcut)
69	peller	1.5
70	peller	1.1
71	nmohr	1.9	input = TFile.Open(path+'/'+job.getpath(),'read')
72	peller	1.1
73			Tree = input.Get(job.tree)
74			#Tree=tmpTree.CloneTree()
75			#Tree.SetDirectory(0)
76
77			#Tree=tmpTree.Clone()
78	peller	1.11	weightF=config.get('Weights',which_weightF)
79	peller	1.1	#hTree = ROOT.TH1F('%s'%name,'%s'%title,nBins,xMin,xMax)
80			#hTree.SetDirectory(0)
81			#hTree.Sumw2()
82			#print 'drawing...'
83			if job.type != 'DATA':
84			#print treeCut
85			#print job.name
86			if Tree.GetEntries():
87			Tree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),'(%s)*(%s)' %(treeCut,weightF), "goff,e")
88			full=True
89			else:
90			full=False
91			elif job.type == 'DATA':
92
93			if len(options)>10:
94			if options[11] == 'blind':
95			treeCut = treeCut + '&'+treeVar+'<0'
96
97
98			Tree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),treeCut, "goff,e")
99			full = True
100			if full:
101			hTree = ROOT.gDirectory.Get(name)
102			else:
103			hTree = ROOT.TH1F('%s'%name,'%s'%name,nBins,xMin,xMax)
104			hTree.Sumw2()
105			#print job.name + ' Sumw2', hTree.GetEntries()
106
107			if job.type != 'DATA':
108	nmohr	1.10	ScaleFactor = getScale(job,path,config,rescale,subsample)
109	peller	1.1	if ScaleFactor != 0:
110			hTree.Scale(ScaleFactor)
111	nmohr	1.15
112			if addOverFlow:
113			print 'Adding overflow'
114			uFlow = hTree.GetBinContent(0)+hTree.GetBinContent(1)
115			oFlow = hTree.GetBinContent(hTree.GetNbinsX()+1)+hTree.GetBinContent(hTree.GetNbinsX())
116			uFlowErr = ROOT.TMath.Sqrt(ROOT.TMath.Power(hTree.GetBinError(0),2)+ROOT.TMath.Power(hTree.GetBinError(1),2))
117			oFlowErr = ROOT.TMath.Sqrt(ROOT.TMath.Power(hTree.GetBinError(hTree.GetNbinsX()),2)+ROOT.TMath.Power(hTree.GetBinError(hTree.GetNbinsX()+1),2))
118			hTree.SetBinContent(1,uFlow)
119			hTree.SetBinContent(hTree.GetNbinsX(),oFlow)
120			hTree.SetBinError(1,uFlowErr)
121			hTree.SetBinError(hTree.GetNbinsX(),oFlowErr)
122
123	peller	1.1
124			print '\t-->import %s\t Integral: %s'%(job.name,hTree.Integral())
125
126			hTree.SetDirectory(0)
127	peller	1.7	input.Close()
128
129			return hTree, group
130	peller	1.1
131
132			######################
133
134
135
136			def orderandadd(histos,typs,setup):
137			#ORDER AND ADD TOGETHER
138
139			ordnung=[]
140			ordnungtyp=[]
141			num=[0]*len(setup)
142			for i in range(0,len(setup)):
143			for j in range(0,len(histos)):
144	peller	1.2	if typs[j] in setup[i]:
145	peller	1.1	num[i]+=1
146			ordnung.append(histos[j])
147			ordnungtyp.append(typs[j])
148
149			del histos
150			del typs
151
152			histos=ordnung
153			typs=ordnungtyp
154
155	peller	1.14	print typs
156
157	peller	1.1	for k in range(0,len(num)):
158			for m in range(0,num[k]):
159			if m > 0:
160
161			#add
162			histos[k].Add(histos[k+1],1)
163	peller	1.2	printc('magenta','','\t--> added %s to %s'%(typs[k],typs[k+1]))
164	peller	1.1	del histos[k+1]
165			del typs[k+1]
166
167			del histos[len(setup):]
168			del typs[len(setup):]
169
170			return histos, typs
171
172