forked from CERN-PH-CMG/cmg-cmssw
-
Notifications
You must be signed in to change notification settings - Fork 5
/
calcSyst.py
executable file
·145 lines (97 loc) · 3.5 KB
/
calcSyst.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
#!/usr/bin/env python
#import re, sys, os, os.path
import glob, os, sys
from math import hypot
from ROOT import *
from readYields import getYield
def getHnames(fname,tdir):
tfile = TFile(fname,"READ")
tfile.cd(tdir)
hnames = []
for key in gDirectory.GetListOfKeys():
obj = key.ReadObj()
hnames.append(obj.GetName())
tfile.Close()
return hnames
def getSystHist(tfile, hname, syst = "Xsec"):
upName = hname + '_' + syst + '-Up'
dnName = hname + '_' + syst + '-Down'
#print tfile, hname, upName, dnName
hNorm = tfile.Get(hname)
hUp = tfile.Get(upName)
hDown = tfile.Get(dnName)
if not hUp or not hDown:
print 'No systematics found!'
return 1
hSyst = hNorm.Clone(hNorm.GetName() + '_' + syst + '_syst')
hPlus = hNorm.Clone(hNorm.GetName() + '_' + syst + '_plus')
hPlus.Add(hUp,-1)
hMinus = hNorm.Clone(hNorm.GetName() + '_' + syst + '_minus')
hMinus.Add(hDown,-1)
# find maximum deviations
for xbin in range(1,hSyst.GetNbinsX()+1):
for ybin in range(1,hSyst.GetNbinsY()+1):
# reset bins
hSyst.SetBinContent(xbin,ybin,0)
hSyst.SetBinError(xbin,ybin,0)
maxDev = 0
maxErr = 0
# fill maximum deviation
if abs(hPlus.GetBinContent(xbin,ybin)) > abs(hMinus.GetBinContent(xbin,ybin)):
maxDev = abs(hPlus.GetBinContent(xbin,ybin))
#maxErr = abs(hPlus.GetBinError(xbin,ybin))
else:
maxDev = abs(hMinus.GetBinContent(xbin,ybin))
#maxErr = abs(hMinus.GetBinError(xbin,ybin))
if hNorm.GetBinContent(xbin,ybin) > 0:
maxDev /= hNorm.GetBinContent(xbin,ybin)
#maxErr = hypot(maxErr,hNorm.GetBinError(xbin,ybin))
hSyst.SetBinContent(xbin,ybin,maxDev)
hSyst.SetBinError(xbin,ybin,maxErr)
return hSyst
def makeSystHists(fileList):
# filter
#fileList = [fname for fname in fileList if 'NB3' not in fname]
hnames = ["T1tttt_Scan"] # process name
#hnames = getHnames(fileList[0],'SR_MB') # get process names from file
#print 'Found these hists:', hnames
systNames = ["Xsec"]
bindirs = ['SR_MB','CR_MB','SR_SB','CR_SB']
for fname in fileList:
tfile = TFile(fname,"UPDATE")
for bindir in bindirs:
for hname in hnames:
for syst in systNames:
hSyst = getSystHist(tfile, bindir+'/'+ hname, syst)
tfile.cd(bindir)
hSyst.Write()
'''
# create Syst folder structure
if not tfile.GetDirectory(bindir+"/Syst"):
tfile.mkdir(bindir+"/Syst")
for hname in hnames:
for syst in systNames:
tfile.cd(bindir+"/Syst")
hSyst = getSystHist(tfile, bindir+'/'+ hname, syst)
hSyst.Write()
else:
print 'Already found syst'
'''
tfile.Close()
return 1
if __name__ == "__main__":
## remove '-b' option
_batchMode = False
if '-b' in sys.argv:
sys.argv.remove('-b')
_batchMode = True
if len(sys.argv) > 1:
pattern = sys.argv[1]
print '# pattern is', pattern
else:
print "No pattern given!"
exit(0)
# find files matching pattern
fileList = glob.glob(pattern+"*.root")
makeSystHists(fileList)
print 'Finished'