-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy pathLatestBinning.py
224 lines (184 loc) · 9.27 KB
/
LatestBinning.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
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
import sys
from differentials.observable import Observable
from collections import namedtuple
from copy import deepcopy
from math import sqrt
import numpy
class AttrDict(dict):
def __init__(self, *args, **kwargs):
super(AttrDict, self).__init__(*args, **kwargs)
self.__dict__ = self
def normalize(l, normalization=1.0):
s = float(sum(l))
return [ e/s*normalization for e in l ]
########################################
# Physical quantities
########################################
YR4_ggF_n3lo = 4.852E+01
YR4_VBF = 3.779E+00
YR4_WH = 1.369E+00
YR4_ZH = 8.824E-01
YR4_ttH = 5.065E-01
YR4_bbH = 4.863E-01
YR4_tH_t_ch = 7.426E-02
YR4_tH_s_ch = 2.875E-03
YR4_tH_W_associated = 0.000E+00
YR4_totalXS = YR4_ggF_n3lo + YR4_VBF + YR4_WH + YR4_ZH + YR4_ttH + YR4_bbH + YR4_tH_t_ch + YR4_tH_s_ch + YR4_tH_W_associated
YR4_xH = YR4_totalXS - YR4_ggF_n3lo
YR4_totalXS_uncertainty = 2.5
# YR4_totalXS = 55.70628722 # pb
# YR4_ggHXS = 48.58 # pb
SM_BR_hgg = 2.270E-03
SM_BR_hzz = 0.02641
SM_ratio_of_BRs = SM_BR_hgg / SM_BR_hzz
SM_ratio_of_BRs_unc_up = +0.00223
SM_ratio_of_BRs_unc_down = -0.00220
Observable.YR4_totalXS = YR4_totalXS
########################################
# Binning and SM differentials
########################################
# ----------------
# Observables other than pth
binning_ptjet = [ -1000.0, 30.0, 55.0, 95.0, 120.0, 200.0, 13000.0 ]
binning_yh = [ 0.0, 0.15, 0.3, 0.6, 0.9, 1.2, 2.5 ]
binning_njets = [ 0., 1.0, 2.0, 3.0, 4.0, 100.0 ]
shape_ptjet = normalize([ 78.53302819, 20.52882251, 13.40260032, 3.79521629, 4.44346963, 1.62691503])
shape_yh = normalize([ 8.88641092, 8.86086463, 17.34621208, 16.38134718, 15.14497828, 51.87741436 ])
shape_njets = normalize([ 78.53302819, 30.8512632, 9.1331588, 2.41446376, 1.39813801 ])
obs_njets = Observable( name = 'njets', title = 'N_{jets}', shape = shape_njets, binning = binning_njets )
obs_yh = Observable( name = 'yh', title = '|y_{H}|', shape = shape_yh, binning = binning_yh, lastBinIsOverflow=False )
obs_ptjet = Observable( name = 'ptjet', title = 'p_{T}^{jet}', shape = shape_ptjet, binning = binning_ptjet )
hzz_binMerging_ptjet = [ 0, 1, 2, [3,4,5] ]
hzz_binMerging_njets = [ 0, 1, 2, [3,4] ]
obs_ptjet_hzzBinning = obs_ptjet.mergeBins(hzz_binMerging_ptjet)
obs_njets_hzzBinning = obs_njets.mergeBins(hzz_binMerging_njets)
# ----------------
# pth
binning_pth = [ 0.0, 15.0, 30.0, 45.0, 80.0, 120.0, 200.0, 350.0, 600., 10000.0 ]
shape_pth_smH = [
2.74660743e+01, 2.94496162e+01, 1.96934684e+01, 2.44876793e+01,
1.17284547e+01, 7.16436043e+00, 2.05821120e+00, 2.63714145e-01,
1.84732619e-02]
shape_pth_ggH = [
2.70051430e+01, 2.81966864e+01, 1.79667435e+01, 2.03254009e+01,
8.39956789e+00, 4.53454528e+00, 1.19397255e+00, 1.39795879e-01,
7.51079245e-03
]
shape_pth_xH = [ s_tot-s_ggH for s_tot, s_ggH in zip(shape_pth_smH, shape_pth_ggH) ]
shape_pth_smH = normalize(shape_pth_smH)
shape_pth_ggH = normalize(shape_pth_ggH)
shape_pth_xH = normalize(shape_pth_xH)
obs_pth_smH = Observable( name = 'pth_smH', title = 'p_{T}^{H}', shape = shape_pth_smH, binning = binning_pth )
obs_pth_ggH = Observable( name = 'pth_ggH', title = 'p_{T}^{H} (ggH)', shape = shape_pth_ggH, binning = binning_pth )
obs_pth_ggH.YR4_totalXS = YR4_ggF_n3lo
obs_pth_xH = Observable( name = 'pth_xH', title = 'p_{T}^{H} (xH)', shape = shape_pth_xH, binning = binning_pth )
obs_pth_xH.YR4_totalXS = YR4_xH
# xH individual shapes
ttH_npz = numpy.load('suppliedInput/fromVittorio/FullPhaseSpace_spectrumNNLOPS_ttH_Pt.npz')
shape_pth_ttH = normalize(ttH_npz['spectrum'])
obs_pth_ttH = Observable( name = 'pth_ttH', title = 'p_{T}^{H} (ttH)', shape = shape_pth_ttH, binning = binning_pth )
obs_pth_ttH.YR4_totalXS = YR4_ttH
# Different binning schemes for HZZ and Hbb
hzz_binMerging_pth = [ 0, 1, [2,3], [4,5], [6,7,8] ]
obs_pth_smH_hzzBinning = obs_pth_smH.mergeBins(hzz_binMerging_pth)
obs_pth_ggH_hzzBinning = obs_pth_ggH.mergeBins(hzz_binMerging_pth)
obs_pth_smH_hbbBinning = deepcopy(obs_pth_smH)
obs_pth_ggH_hbbBinning = deepcopy(obs_pth_ggH)
# !!!: Dropping up to 350 now!! 200-350 no longer included; careful when multiplying this with a scan!
for i in xrange(7):
obs_pth_smH_hbbBinning.drop_first_bin()
obs_pth_ggH_hbbBinning.drop_first_bin()
#____________________________________________________________________
# Convenience tuples
obstuple_pth_smH = AttrDict(
hgg = obs_pth_smH,
hzz = obs_pth_smH_hzzBinning,
combination = obs_pth_smH,
combWithHbb = obs_pth_smH,
hbb = obs_pth_smH_hbbBinning
)
obstuple_pth_ggH = AttrDict(
hgg = obs_pth_ggH,
hzz = obs_pth_ggH_hzzBinning,
combination = obs_pth_ggH,
combWithHbb = obs_pth_ggH,
hbb = obs_pth_ggH_hbbBinning
)
obstuple_njets = AttrDict(
hgg = obs_njets,
hzz = obs_njets_hzzBinning,
combination = obs_njets,
combWithHbb = obs_njets,
)
obstuple_ptjet = AttrDict(
hgg = obs_ptjet,
hzz = obs_ptjet_hzzBinning,
combination = obs_ptjet,
combWithHbb = obs_ptjet,
)
obstuple_rapidity = AttrDict(
hgg = obs_yh,
hzz = obs_yh,
combination = obs_yh,
combWithHbb = obs_yh,
)
#____________________________________________________________________
# xH uncertainty, inclusive
# Uncertainties per mode, all in percentages of total XS
# first one is scale, second PDF, third alpha_s
uncs_VBF = [ 0.35, 2.1, 0.05 ]
uncs_WH = [ 0.6, 1.7, 0.9 ]
uncs_ZH = [ 3.4, 1.3, 0.9 ]
uncs_ttH = [ 7.5, 3.0, 2.0 ]
uncs_bbH = [ 22.0 ]
uncs_tH_t_ch = [ 10.6, 3.5, 1.2 ]
uncs_tH_s_ch = [ 2.1, 2.2, 0.2 ]
uncs_tH_W_associated = [ 5.8, 6.1, 1.5 ]
unc_squared_per_mode = lambda uncs, xs: sum([ (0.01*unc * xs)**2 for unc in uncs ])
tot_unc_squared = 0.0
tot_unc_squared += unc_squared_per_mode(uncs_VBF, YR4_VBF)
tot_unc_squared += unc_squared_per_mode(uncs_WH, YR4_WH)
tot_unc_squared += unc_squared_per_mode(uncs_ZH, YR4_ZH)
tot_unc_squared += unc_squared_per_mode(uncs_ttH, YR4_ttH)
tot_unc_squared += unc_squared_per_mode(uncs_bbH, YR4_bbH)
tot_unc_squared += unc_squared_per_mode(uncs_tH_t_ch, YR4_tH_t_ch)
tot_unc_squared += unc_squared_per_mode(uncs_tH_s_ch, YR4_tH_s_ch)
tot_unc_squared += unc_squared_per_mode(uncs_tH_W_associated, YR4_tH_W_associated)
xH_unc_inclusive = sqrt(tot_unc_squared)
xH_unc_inclusive_fraction = xH_unc_inclusive / YR4_xH
uncs_ggF = [ 5.65, 3.2 ]
tot_unc_squared += unc_squared_per_mode(uncs_ggF, YR4_ggF_n3lo)
smH_unc_inclusive = sqrt(tot_unc_squared)
smH_unc_inclusive_fraction = smH_unc_inclusive / YR4_totalXS
#____________________________________________________________________
# smH uncertainties
add_quad = lambda *args: sqrt(sum([ x**2 for x in args ]))
def add_unc(obs, unc_file):
xs = obs.crosssection()
shape_unc_perc = list(numpy.load(unc_file)['uncetainty'])
# Add in quadrature with incl xs uncertainty from YR4
# Do this all in 'percent-space'
unc_perc = [ add_quad(err, smH_unc_inclusive_fraction) for err in shape_unc_perc ]
obs.unc_fraction = unc_perc
add_unc(obs_yh, 'suppliedInput/fromVittorio/sm_uncertainties/uncertaintyFullPhaseSpaceCombination_AbsRapidity.npz')
add_unc(obs_ptjet, 'suppliedInput/fromVittorio/sm_uncertainties/uncertaintyFullPhaseSpaceCombination_Jet2p5Pt0.npz')
add_unc(obs_njets, 'suppliedInput/fromVittorio/sm_uncertainties/uncertaintyFullPhaseSpaceCombination_Njets2p5.npz')
add_unc(obs_pth_smH, 'suppliedInput/fromVittorio/sm_uncertainties/uncertaintyFullPhaseSpaceCombination_Pt.npz')
add_unc(obs_pth_ggH, 'suppliedInput/fromVittorio/sm_uncertainties/uncertaintyFullPhaseSpaceCombination_ggH_Pt.npz')
if __name__ == "__main__":
print 'Unc from VBF: {0:.2e} , xs = {1:.2e}'.format(sqrt(unc_squared_per_mode(uncs_VBF, YR4_VBF)), YR4_VBF)
print 'Unc from WH: {0:.2e} , xs = {1:.2e}'.format(sqrt(unc_squared_per_mode(uncs_WH, YR4_WH)), YR4_WH)
print 'Unc from ZH: {0:.2e} , xs = {1:.2e}'.format(sqrt(unc_squared_per_mode(uncs_ZH, YR4_ZH)), YR4_ZH)
print 'Unc from ttH: {0:.2e} , xs = {1:.2e}'.format(sqrt(unc_squared_per_mode(uncs_ttH, YR4_ttH)), YR4_ttH)
print 'Unc from bbH: {0:.2e} , xs = {1:.2e}'.format(sqrt(unc_squared_per_mode(uncs_bbH, YR4_bbH)), YR4_bbH)
print 'Unc from tH_t_ch: {0:.2e} , xs = {1:.2e}'.format(sqrt(unc_squared_per_mode(uncs_tH_t_ch, YR4_tH_t_ch)), YR4_tH_t_ch)
print 'Unc from tH_s_ch: {0:.2e} , xs = {1:.2e}'.format(sqrt(unc_squared_per_mode(uncs_tH_s_ch, YR4_tH_s_ch)), YR4_tH_s_ch)
print 'Unc from tH_W_associated: {0:.2e} , xs = {1:.2e}'.format(sqrt(unc_squared_per_mode(uncs_tH_W_associated, YR4_tH_W_associated)), YR4_tH_W_associated)
print
print 'YR4_xH: ', YR4_xH
print 'xH_unc_inclusive: ', xH_unc_inclusive
print 'xH_unc_inclusive_fraction:', xH_unc_inclusive_fraction
print 'smH_unc_inclusive: ', smH_unc_inclusive
print 'smH_unc_inclusive_fraction:', smH_unc_inclusive_fraction
# print list(numpy.load('suppliedInput/fromVittorio/sm_uncertainties/uncertaintyFullPhaseSpaceCombination_AbsRapidity.npz')['uncetainty'])
# print obs_yh.unc_fraction