ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbb/python/myutils/HistoMaker.py
Revision: 1.1
Committed: Wed Jan 16 16:35:40 2013 UTC (12 years, 3 months ago) by peller
Content type: text/x-python
Branch: MAIN
CVS Tags: workingVersionAfterHCP
Log Message:
reorganized the whole repository. Macros im myutils, config files in subdirectories. Config file split in parts. Path config file restructured. Moved all path options to the path config. Changed the code accordingly.

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