| 1 |
import sys,os
|
| 2 |
import pickle
|
| 3 |
import ROOT
|
| 4 |
from array import array
|
| 5 |
from printcolor import printc
|
| 6 |
from BetterConfigParser import BetterConfigParser
|
| 7 |
from TreeCache import TreeCache
|
| 8 |
from math import sqrt
|
| 9 |
from copy import copy
|
| 10 |
|
| 11 |
class HistoMaker:
|
| 12 |
def __init__(self, samples, path, config, optionsList,GroupDict=None):
|
| 13 |
self.path = path
|
| 14 |
self.config = config
|
| 15 |
self.optionsList = optionsList
|
| 16 |
self.nBins = optionsList[0]['nBins']
|
| 17 |
self.lumi=0.
|
| 18 |
self.cuts = []
|
| 19 |
for options in optionsList:
|
| 20 |
self.cuts.append(options['cut'])
|
| 21 |
#print self.cuts
|
| 22 |
#self.tc = TreeCache(self.cuts,samples,path)
|
| 23 |
self.tc = TreeCache(self.cuts,samples,path,config)
|
| 24 |
self._rebin = False
|
| 25 |
self.mybinning = None
|
| 26 |
self.GroupDict=GroupDict
|
| 27 |
self.calc_rebin_flag = False
|
| 28 |
VHbbNameSpace=config.get('VHbbNameSpace','library')
|
| 29 |
ROOT.gSystem.Load(VHbbNameSpace)
|
| 30 |
|
| 31 |
def get_histos_from_tree(self,job,cutOverWrite=None):
|
| 32 |
if self.lumi == 0:
|
| 33 |
raise Exception("You're trying to plot with no lumi")
|
| 34 |
|
| 35 |
hTreeList=[]
|
| 36 |
|
| 37 |
#get the conversion rate in case of BDT plots
|
| 38 |
TrainFlag = eval(self.config.get('Analysis','TrainFlag'))
|
| 39 |
BDT_add_cut='EventForTraining == 0'
|
| 40 |
|
| 41 |
|
| 42 |
plot_path = self.config.get('Directories','plotpath')
|
| 43 |
addOverFlow=eval(self.config.get('Plot_general','addOverFlow'))
|
| 44 |
|
| 45 |
# get all Histos at once
|
| 46 |
CuttedTree = self.tc.get_tree(job,'1')
|
| 47 |
for options in self.optionsList:
|
| 48 |
name=job.name
|
| 49 |
if self.GroupDict is None:
|
| 50 |
group=job.group
|
| 51 |
else:
|
| 52 |
group=self.GroupDict[job.name]
|
| 53 |
treeVar=options['var']
|
| 54 |
name=options['name']
|
| 55 |
if self._rebin or self.calc_rebin_flag:
|
| 56 |
nBins = self.nBins
|
| 57 |
else:
|
| 58 |
nBins = int(options['nBins'])
|
| 59 |
xMin=float(options['xMin'])
|
| 60 |
xMax=float(options['xMax'])
|
| 61 |
weightF=options['weight']
|
| 62 |
if cutOverWrite:
|
| 63 |
treeCut=cutOverWrite
|
| 64 |
else:
|
| 65 |
treeCut='%s'%(options['cut'])
|
| 66 |
|
| 67 |
#options
|
| 68 |
|
| 69 |
if job.type != 'DATA':
|
| 70 |
if CuttedTree.GetEntries():
|
| 71 |
if 'RTight' in treeVar or 'RMed' in treeVar:
|
| 72 |
drawoption = '(%s)*(%s & %s)'%(weightF,treeCut,BDT_add_cut)
|
| 73 |
#print drawoption
|
| 74 |
else:
|
| 75 |
drawoption = '(%s)*(%s)'%(weightF,treeCut)
|
| 76 |
CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax), drawoption, "goff,e")
|
| 77 |
full=True
|
| 78 |
else:
|
| 79 |
full=False
|
| 80 |
elif job.type == 'DATA':
|
| 81 |
if options['blind']:
|
| 82 |
if treeVar == 'H.mass':
|
| 83 |
CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),' (%(var)s <90. || %(var)s > 150.) & %(cut)s' %options, "goff,e")
|
| 84 |
else:
|
| 85 |
CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),'%(var)s < 0. & %(cut)s'%options, "goff,e")
|
| 86 |
|
| 87 |
else:
|
| 88 |
CuttedTree.Draw('%s>>%s(%s,%s,%s)' %(treeVar,name,nBins,xMin,xMax),'%s' %treeCut, "goff,e")
|
| 89 |
full = True
|
| 90 |
if full:
|
| 91 |
hTree = ROOT.gDirectory.Get(name)
|
| 92 |
else:
|
| 93 |
hTree = ROOT.TH1F('%s'%name,'%s'%name,nBins,xMin,xMax)
|
| 94 |
hTree.Sumw2()
|
| 95 |
if job.type != 'DATA':
|
| 96 |
if 'RTight' in treeVar or 'RMed' in treeVar:
|
| 97 |
if TrainFlag:
|
| 98 |
MC_rescale_factor=2.
|
| 99 |
#print 'I RESCALE BY 2.0'
|
| 100 |
else:
|
| 101 |
MC_rescale_factor = 1.
|
| 102 |
ScaleFactor = self.tc.get_scale(job,self.config,self.lumi)*MC_rescale_factor
|
| 103 |
else:
|
| 104 |
ScaleFactor = self.tc.get_scale(job,self.config,self.lumi)
|
| 105 |
if ScaleFactor != 0:
|
| 106 |
hTree.Scale(ScaleFactor)
|
| 107 |
#print '\t-->import %s\t Integral: %s'%(job.name,hTree.Integral())
|
| 108 |
if addOverFlow:
|
| 109 |
uFlow = hTree.GetBinContent(0)+hTree.GetBinContent(1)
|
| 110 |
oFlow = hTree.GetBinContent(hTree.GetNbinsX()+1)+hTree.GetBinContent(hTree.GetNbinsX())
|
| 111 |
uFlowErr = ROOT.TMath.Sqrt(ROOT.TMath.Power(hTree.GetBinError(0),2)+ROOT.TMath.Power(hTree.GetBinError(1),2))
|
| 112 |
oFlowErr = ROOT.TMath.Sqrt(ROOT.TMath.Power(hTree.GetBinError(hTree.GetNbinsX()),2)+ROOT.TMath.Power(hTree.GetBinError(hTree.GetNbinsX()+1),2))
|
| 113 |
hTree.SetBinContent(1,uFlow)
|
| 114 |
hTree.SetBinContent(hTree.GetNbinsX(),oFlow)
|
| 115 |
hTree.SetBinError(1,uFlowErr)
|
| 116 |
hTree.SetBinError(hTree.GetNbinsX(),oFlowErr)
|
| 117 |
hTree.SetDirectory(0)
|
| 118 |
gDict = {}
|
| 119 |
if self._rebin:
|
| 120 |
gDict[group] = self.mybinning.rebin(hTree)
|
| 121 |
del hTree
|
| 122 |
else:
|
| 123 |
#print 'not rebinning %s'%job.name
|
| 124 |
gDict[group] = hTree
|
| 125 |
hTreeList.append(gDict)
|
| 126 |
CuttedTree.IsA().Destructor(CuttedTree)
|
| 127 |
del CuttedTree
|
| 128 |
return hTreeList
|
| 129 |
|
| 130 |
@property
|
| 131 |
def rebin(self):
|
| 132 |
return self._rebin
|
| 133 |
|
| 134 |
@property
|
| 135 |
def rebin(self, value):
|
| 136 |
if self._rebin and value:
|
| 137 |
return True
|
| 138 |
elif self._rebin and not value:
|
| 139 |
self.nBins = self.norebin_nBins
|
| 140 |
self._rebin = False
|
| 141 |
elif not self._rebin and value:
|
| 142 |
if self.mybinning is None:
|
| 143 |
raise Exception('define rebinning first')
|
| 144 |
else:
|
| 145 |
self.nBins = self.rebin_nBins
|
| 146 |
self._rebin = True
|
| 147 |
return True
|
| 148 |
elif not self._rebin and not self.value:
|
| 149 |
return False
|
| 150 |
|
| 151 |
def calc_rebin(self, bg_list, nBins_start=1000, tolerance=0.25):
|
| 152 |
self.calc_rebin_flag = True
|
| 153 |
self.norebin_nBins = copy(self.nBins)
|
| 154 |
self.rebin_nBins = nBins_start
|
| 155 |
self.nBins = nBins_start
|
| 156 |
i=0
|
| 157 |
#add all together:
|
| 158 |
print '\n\t...calculating rebinning...'
|
| 159 |
for job in bg_list:
|
| 160 |
htree = self.get_histos_from_tree(job)[0].values()[0]
|
| 161 |
if not i:
|
| 162 |
totalBG = copy(htree)
|
| 163 |
else:
|
| 164 |
totalBG.Add(htree,1)
|
| 165 |
del htree
|
| 166 |
i+=1
|
| 167 |
ErrorR=0
|
| 168 |
ErrorL=0
|
| 169 |
TotR=0
|
| 170 |
TotL=0
|
| 171 |
binR=self.rebin_nBins
|
| 172 |
binL=1
|
| 173 |
rel=1.0
|
| 174 |
#---- from right
|
| 175 |
while rel > tolerance:
|
| 176 |
TotR+=totalBG.GetBinContent(binR)
|
| 177 |
ErrorR=sqrt(ErrorR**2+totalBG.GetBinError(binR)**2)
|
| 178 |
binR-=1
|
| 179 |
if not TotR == 0 and not ErrorR == 0:
|
| 180 |
rel=ErrorR/TotR
|
| 181 |
#print rel
|
| 182 |
#print 'upper bin is %s'%binR
|
| 183 |
|
| 184 |
#---- from left
|
| 185 |
rel=1.0
|
| 186 |
while rel > tolerance:
|
| 187 |
TotL+=totalBG.GetBinContent(binL)
|
| 188 |
ErrorL=sqrt(ErrorL**2+totalBG.GetBinError(binL)**2)
|
| 189 |
binL+=1
|
| 190 |
if not TotL == 0 and not ErrorL == 0:
|
| 191 |
rel=ErrorL/TotL
|
| 192 |
#print rel
|
| 193 |
#it's the lower edge
|
| 194 |
binL+=1
|
| 195 |
#print 'lower bin is %s'%binL
|
| 196 |
|
| 197 |
inbetween=binR-binL
|
| 198 |
stepsize=int(inbetween)/(int(self.norebin_nBins)-2)
|
| 199 |
modulo = int(inbetween)%(int(self.norebin_nBins)-2)
|
| 200 |
|
| 201 |
#print 'stepsize %s'% stepsize
|
| 202 |
#print 'modulo %s'%modulo
|
| 203 |
binlist=[binL]
|
| 204 |
for i in range(0,int(self.norebin_nBins)-3):
|
| 205 |
binlist.append(binlist[-1]+stepsize)
|
| 206 |
binlist[-1]+=modulo
|
| 207 |
binlist.append(binR)
|
| 208 |
binlist.append(self.rebin_nBins+1)
|
| 209 |
#print 'binning set to %s'%binlist
|
| 210 |
self.mybinning = Rebinner(int(self.norebin_nBins),array('d',[-1.0]+[totalBG.GetBinLowEdge(i) for i in binlist]),True)
|
| 211 |
self._rebin = True
|
| 212 |
print '\t > rebinning is set <\n'
|
| 213 |
|
| 214 |
@staticmethod
|
| 215 |
def orderandadd(histo_dicts,setup):
|
| 216 |
ordered_histo_dict = {}
|
| 217 |
for sample in setup:
|
| 218 |
nSample = 0
|
| 219 |
for histo_dict in histo_dicts:
|
| 220 |
if histo_dict.has_key(sample):
|
| 221 |
if nSample == 0:
|
| 222 |
ordered_histo_dict[sample] = histo_dict[sample].Clone()
|
| 223 |
else:
|
| 224 |
printc('magenta','','\t--> added %s to %s'%(sample,sample))
|
| 225 |
ordered_histo_dict[sample].Add(histo_dict[sample])
|
| 226 |
nSample += 1
|
| 227 |
del histo_dicts
|
| 228 |
return ordered_histo_dict
|
| 229 |
|
| 230 |
class Rebinner:
|
| 231 |
def __init__(self,nBins,lowedgearray,active=True):
|
| 232 |
self.lowedgearray=lowedgearray
|
| 233 |
self.nBins=nBins
|
| 234 |
self.active=active
|
| 235 |
def rebin(self, histo):
|
| 236 |
if not self.active: return histo
|
| 237 |
#print histo.Integral()
|
| 238 |
ROOT.gDirectory.Delete('hnew')
|
| 239 |
histo.Rebin(self.nBins,'hnew',self.lowedgearray)
|
| 240 |
binhisto=ROOT.gDirectory.Get('hnew')
|
| 241 |
#print binhisto.Integral()
|
| 242 |
newhisto=ROOT.TH1F('new','new',self.nBins,self.lowedgearray[0],self.lowedgearray[-1])
|
| 243 |
newhisto.Sumw2()
|
| 244 |
for bin in range(1,self.nBins+1):
|
| 245 |
newhisto.SetBinContent(bin,binhisto.GetBinContent(bin))
|
| 246 |
newhisto.SetBinError(bin,binhisto.GetBinError(bin))
|
| 247 |
newhisto.SetName(binhisto.GetName())
|
| 248 |
newhisto.SetTitle(binhisto.GetTitle())
|
| 249 |
#print newhisto.Integral()
|
| 250 |
del histo
|
| 251 |
del binhisto
|
| 252 |
return copy(newhisto)
|
