ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbb/python/HistoMaker.py
Revision: 1.11
Committed: Thu Nov 22 14:13:41 2012 UTC (12 years, 5 months ago) by nmohr
Content type: text/x-python
Branch: MAIN
Changes since 1.10: +1 -0 lines
Log Message:
Fix open input file

File Contents

# Content
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,os
10
11 class HistoMaker:
12 def __init__(self, path, config, region, optionsList,rescale=1,which_weightF='weightF'):
13 self.path = path
14 self.config = config
15 self.optionsList = optionsList
16 self.rescale = rescale
17 self.which_weightF=which_weightF
18 self.region = region
19 self.lumi=0.
20
21 def getScale(self,job,subsample=-1,MC_rescale_factor=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 theScale = float(self.lumi)*xsec*sf/(0.46502*CountWithPU.GetBinContent(1)+0.53498*CountWithPU2011B.GetBinContent(1))*MC_rescale_factor/float(job.split)
36 elif anaTag == '8TeV':
37 theScale = float(self.lumi)*xsec*sf/(CountWithPU.GetBinContent(1))*MC_rescale_factor/float(job.split)
38 input.Close()
39 return theScale
40
41
42 def getHistoFromTree(self,job,subsample=-1):
43 if self.lumi == 0: raise Exception("You're trying to plot with no lumi")
44
45 hTreeList=[]
46 groupList=[]
47
48 #get the conversion rate in case of BDT plots
49 TrainFlag = eval(self.config.get('Analysis','TrainFlag'))
50 BDT_add_cut='EventForTraining == 0'
51
52
53 plot_path = self.config.get('Directories','plotpath')
54 addOverFlow=eval(self.config.get('Plot_general','addOverFlow'))
55
56 scratchDir = os.environ["TMPDIR"]
57 #scratchDir = '/shome/peller/'
58 # define treeCut
59 if job.type != 'DATA':
60 if type(self.region)==str:
61 cutcut=self.config.get('Cuts',self.region)
62 elif type(self.region)==list:
63 #replace vars with other vars in the cutstring (used in DC writer)
64 cutcut=self.config.get('Cuts',self.region[0])
65 cutcut=cutcut.replace(self.region[1],self.region[2])
66 #print cutcut
67 if subsample>-1:
68 treeCut='%s & %s'%(cutcut,job.subcuts[subsample])
69 else:
70 treeCut='%s'%(cutcut)
71 elif job.type == 'DATA':
72 cutcut=self.config.get('Cuts',self.region)
73 treeCut='%s'%(cutcut)
74
75 # get and skim the Trees
76 output=TFile.Open(scratchDir+'/tmp_plotCache_%s_%s.root'%(self.region,job.identifier),'recreate')
77 input = TFile.Open(self.path+'/'+job.getpath(),'read')
78 Tree = input.Get(job.tree)
79 output.cd()
80 CuttedTree=Tree.CopyTree(treeCut)
81 input.Close()
82 del input
83 # get all Histos at once
84 weightF=self.config.get('Weights',self.which_weightF)
85 for options in self.optionsList:
86 if subsample>-1:
87 name=job.subnames[subsample]
88 group=job.group[subsample]
89 else:
90 name=job.name
91 group=job.group
92 treeVar=options[0]
93 name=options[1]
94 nBins=int(options[3])
95 xMin=float(options[4])
96 xMax=float(options[5])
97
98 #options
99
100 if job.type != 'DATA':
101 if CuttedTree.GetEntries():
102
103 if 'RTight' in treeVar or 'RMed' in treeVar: drawoption = '(%s)*(%s)'%(weightF,BDT_add_cut)
104 else: drawoption = '%s'%(weightF)
105 output.cd()
106 CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax), drawoption, "goff,e")
107 full=True
108 else:
109 full=False
110 elif job.type == 'DATA':
111 if options[11] == 'blind':
112 output.cd()
113 if treeVar == 'H.mass':
114 CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),treeVar+'<90. || '+treeVar + '>150.' , "goff,e")
115 else:
116 CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),treeVar+'<0', "goff,e")
117
118 else:
119 output.cd()
120 CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),'1', "goff,e")
121 full = True
122 if full:
123 hTree = ROOT.gDirectory.Get(name)
124 else:
125 output.cd()
126 hTree = ROOT.TH1F('%s'%name,'%s'%name,nBins,xMin,xMax)
127 hTree.Sumw2()
128 if job.type != 'DATA':
129 if 'RTight' in treeVar or 'RMed' in treeVar:
130 if TrainFlag:
131 MC_rescale_factor=2.
132 print 'I RESCALE BY 2.0'
133 else: MC_rescale_factor = 1.
134 ScaleFactor = self.getScale(job,subsample,MC_rescale_factor)
135 else: ScaleFactor = self.getScale(job,subsample)
136 if ScaleFactor != 0:
137 hTree.Scale(ScaleFactor)
138 #print '\t-->import %s\t Integral: %s'%(job.name,hTree.Integral())
139 if addOverFlow:
140 uFlow = hTree.GetBinContent(0)+hTree.GetBinContent(1)
141 oFlow = hTree.GetBinContent(hTree.GetNbinsX()+1)+hTree.GetBinContent(hTree.GetNbinsX())
142 uFlowErr = ROOT.TMath.Sqrt(ROOT.TMath.Power(hTree.GetBinError(0),2)+ROOT.TMath.Power(hTree.GetBinError(1),2))
143 oFlowErr = ROOT.TMath.Sqrt(ROOT.TMath.Power(hTree.GetBinError(hTree.GetNbinsX()),2)+ROOT.TMath.Power(hTree.GetBinError(hTree.GetNbinsX()+1),2))
144 hTree.SetBinContent(1,uFlow)
145 hTree.SetBinContent(hTree.GetNbinsX(),oFlow)
146 hTree.SetBinError(1,uFlowErr)
147 hTree.SetBinError(hTree.GetNbinsX(),oFlowErr)
148 hTree.SetDirectory(0)
149 hTreeList.append(hTree)
150 groupList.append(group)
151
152 output.Close()
153 del output
154 return hTreeList, groupList
155
156
157 ######################
158 def orderandadd(histos,typs,setup):
159 #ORDER AND ADD TOGETHER
160 ordnung=[]
161 ordnungtyp=[]
162 num=[0]*len(setup)
163 for i in range(0,len(setup)):
164 for j in range(0,len(histos)):
165 if typs[j] in setup[i]:
166 num[i]+=1
167 ordnung.append(histos[j])
168 ordnungtyp.append(typs[j])
169 del histos
170 del typs
171 histos=ordnung
172 typs=ordnungtyp
173 print typs
174 for k in range(0,len(num)):
175 for m in range(0,num[k]):
176 if m > 0:
177 #add
178 histos[k].Add(histos[k+1],1)
179 printc('magenta','','\t--> added %s to %s'%(typs[k],typs[k+1]))
180 del histos[k+1]
181 del typs[k+1]
182 del histos[len(setup):]
183 del typs[len(setup):]
184 return histos, typs