Skip to content

Commit b6fe97a

Browse files
andrzejnovakandrzejnovak
committed
missing scripts
1 parent 5da5ae2 commit b6fe97a

File tree

4 files changed

+175
-10
lines changed

4 files changed

+175
-10
lines changed

‎README.md‎

Lines changed: 31 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -242,7 +242,7 @@ combineTool.py -M MultiDimFit -d model_combined.root --algo singles --cminDefaul
242242
# Split year form full fit
243243
text2workspace.py -P HiggsAnalysis.CombinedLimit.PhysicsModel:multiSignalModel --PO verbose --PO 'map=.*/zcc:z[1,0,2]' --PO 'map=.*2016/*hcc*:r16[1,-200,300]' --PO 'map=.*2017/*hcc*:r17[1,-200,200]' --PO 'map=.*2018/*hcc*:r18[1,-200,200]' model_combined.txt
244244

245-
combineTool.py -M MultiDimFit -d model_combined.root --algo singles --cminDefaultMinimizerStrategy 0 --setParameters z=1,r16=1,r17=1,r18=1 --redefineSignalPOIs r16 -n .r16 --freezeParameters z
245+
combineTool.py -M MultiDimFit -d model_combined.root --algo singles --cminDefaultMinimizerStrategy 0 --setParameters z=1,r16=1,r17=1,r18=1 --redefineSignalPOIs r16 -n .r16 --freezeParameters z
246246
combineTool.py -M MultiDimFit -d model_combined.root --algo singles --cminDefaultMinimizerStrategy 0 --setParameters z=1,r16=1,r17=1,r18=1 --redefineSignalPOIs r17 -n .r17 --freezeParameters z
247247
combineTool.py -M MultiDimFit -d model_combined.root --algo singles --cminDefaultMinimizerStrategy 0 --setParameters z=1,r16=1,r17=1,r18=1 --redefineSignalPOIs r18 -n .r18 --freezeParameters z
248248

@@ -275,6 +275,33 @@ combineTool.py -M AsymptoticLimits -m 125 -d model_combined.root --freezeParamet
275275
combineTool.py -M AsymptoticLimits -m 125 -d model_combined.root --freezeParameters z,r16,r17 --setParameters z=1,r16=1,r17=1,r18=1 --redefineSignalPOIs r18 -n .r18 &
276276
```
277277

278+
### Unblind Z
279+
280+
```bash
281+
bash build.sh
282+
combine -M Significance model_combined.root --setParameters r=1,z=1 --redefineSignalPOIs z
283+
284+
text2workspace.py -P HiggsAnalysis.CombinedLimit.PhysicsModel:multiSignalModel --PO verbose --PO 'map=.*/*hcc*:r[1,-500,500]' --PO 'map=.*2016/zcc:z16[1,0,3]' --PO 'map=.*2017/zcc:z17[1,0,3]' --PO 'map=.*2018/zcc:z18[1,0,3]' model_combined.txt
285+
286+
combine -M Significance model_combined.root --setParameters r=1,z16=1,z17=1,z18=1 --freezeParameters z17,z18 --redefineSignalPOIs z16 -n .z16 &
287+
combine -M Significance model_combined.root --setParameters r=1,z16=1,z17=1,z18=1 --freezeParameters z16,z18 --redefineSignalPOIs z17 -n .z17 &
288+
combine -M Significance model_combined.root --setParameters r=1,z16=1,z17=1,z18=1 --freezeParameters z16,z17 --redefineSignalPOIs z18 -n .z18 &
289+
```
290+
291+
### ToysF Z
292+
293+
```bash
294+
bash build.sh
295+
combine -M Significance model_combined.root --setParameters r=1,z=1 --redefineSignalPOIs z -t -1 --toysFrequentist
296+
297+
text2workspace.py -P HiggsAnalysis.CombinedLimit.PhysicsModel:multiSignalModel --PO verbose --PO 'map=.*/*hcc*:r[1,-500,500]' --PO 'map=.*2016/zcc:z16[1,0,3]' --PO 'map=.*2017/zcc:z17[1,0,3]' --PO 'map=.*2018/zcc:z18[1,0,3]' model_combined.txt
298+
299+
combine -M Significance model_combined.root --setParameters r=1,z16=1,z17=1,z18=1 --freezeParameters z17,z18 --redefineSignalPOIs z16 -t -1 --toysFrequentist -n .z16 &
300+
combine -M Significance model_combined.root --setParameters r=1,z16=1,z17=1,z18=1 --freezeParameters z16,z18 --redefineSignalPOIs z17 -t -1 --toysFrequentist -n .z17 &
301+
combine -M Significance model_combined.root --setParameters r=1,z16=1,z17=1,z18=1 --freezeParameters z16,z17 --redefineSignalPOIs z18 -t -1 --toysFrequentist -n .z18 &
302+
```
303+
304+
278305

279306
### Uncertainties split stat+syst
280307
```
@@ -302,7 +329,7 @@ combine -M MultiDimFit higgsCombine.postfitZ.MultiDimFit.mH120.root -n .Ztotal -
302329
combine -M MultiDimFit higgsCombine.postfitZ.MultiDimFit.mH120.root -n .Znoexpsyst --algo grid --snapshotName MultiDimFit --setParameterRanges z=0,2 --setParameters r=1,z=1 --redefineSignalPOIs z --freezeParameters 'rgx{(?!.*_EW$)(?!.*_NLO$)(?!.*_th_scale.pt$).*}' &
303330
combine -M MultiDimFit higgsCombine.postfitZ.MultiDimFit.mH120.root -n .Zfreezeall --algo grid --snapshotName MultiDimFit --setParameterRanges z=0,2 --setParameters r=1,z=1 --redefineSignalPOIs z --freezeParameters allConstrainedNuisances &
304331
### Replace X year
305-
plot1DScan.py higgsCombine.Ztotal.MultiDimFit.mH120.root --POI z --main-label "Total Uncert." --others higgsCombine.Znoexpsyst.MultiDimFit.mH120.root:"Theory Uncert.":4 higgsCombine.Zfreezeall.MultiDimFit.mH120.root:"Statistical Only":2 --output plots/breakdown_Z_2017 --y-max 22 --y-cut 10 --breakdown "exp,thy,stat"
332+
plot1DScan.py higgsCombine.Ztotal.MultiDimFit.mH120.root --POI z --main-label "Total Uncert." --others higgsCombine.Znoexpsyst.MultiDimFit.mH120.root:"Theory Uncert.":4 higgsCombine.Zfreezeall.MultiDimFit.mH120.root:"Statistical Only":2 --y-max 22 --y-cut 10 --breakdown "syst,theo,stat" --logo-sub Preliminary --translate ../breakdown.json --output plots/breakdown_Z_201X
306333
```
307334
Higgs 2016 - shifted range
308335
```
@@ -312,7 +339,7 @@ combine -M MultiDimFit higgsCombinepostfit.MultiDimFit.mH120.root -n total --alg
312339
combine -M MultiDimFit higgsCombinepostfit.MultiDimFit.mH120.root -n noexpsyst --algo grid --snapshotName MultiDimFit --setParameterRanges r=-100,300 --setParameters r=1,z=1 --freezeParameters 'rgx{(?!.*_EW$)(?!.*_NLO$)(?!.*_th_scale.pt$)(?!z$).*}' &
313340
combine -M MultiDimFit higgsCombinepostfit.MultiDimFit.mH120.root -n freezeall --algo grid --snapshotName MultiDimFit --setParameterRanges r=-100,300 --setParameters r=1,z=1 --freezeParameters=allConstrainedNuisances,z &
314341
###
315-
plot1DScan.py higgsCombinetotal.MultiDimFit.mH120.root --POI r --main-label "Total Uncert." --others higgsCombinenoexpsyst.MultiDimFit.mH120.root:"Theory Uncert.":4 higgsCombinefreezeall.MultiDimFit.mH120.root:"Statistical Only":2 --output plots/breakdown_H_2016 --y-max 22 --y-cut 10 --breakdown "exp,thy,stat"
342+
plot1DScan.py higgsCombinetotal.MultiDimFit.mH120.root --POI r --main-label "Total Uncert." --others higgsCombinenoexpsyst.MultiDimFit.mH120.root:"Theory Uncert.":4 higgsCombinefreezeall.MultiDimFit.mH120.root:"Statistical Only":2 --y-max 22 --y-cut 10 --breakdown "syst,theo,stat" --logo-sub Preliminary --translate ../breakdown.json --output plots/breakdown_H_2016
316343
```
317344

318345
Higgs
@@ -323,7 +350,7 @@ combine -M MultiDimFit higgsCombine.postfit.MultiDimFit.mH120.root -n .total --a
323350
combine -M MultiDimFit higgsCombine.postfit.MultiDimFit.mH120.root -n .noexpsyst --algo grid --snapshotName MultiDimFit --setParameterRanges r=-100,100 --setParameters r=1,z=1 --freezeParameters 'rgx{(?!.*_EW$)(?!.*_NLO$)(?!.*_th_scale.pt$)(?!z$).*}' &
324351
combine -M MultiDimFit higgsCombine.postfit.MultiDimFit.mH120.root -n .freezeall --algo grid --snapshotName MultiDimFit --setParameterRanges r=-100,100 --setParameters r=1,z=1 --freezeParameters=allConstrainedNuisances,z &
325352
### Replace X year
326-
plot1DScan.py higgsCombine.total.MultiDimFit.mH120.root --POI r --main-label "Total Uncert." --others higgsCombine.noexpsyst.MultiDimFit.mH120.root:"Theory Uncert.":4 higgsCombine.freezeall.MultiDimFit.mH120.root:"Statistical Only":2 --output plots/breakdown_H_201X --y-max 22 --y-cut 10 --breakdown "exp,thy,stat"
353+
plot1DScan.py higgsCombine.total.MultiDimFit.mH120.root --POI r --main-label "Total Uncert." --others higgsCombine.noexpsyst.MultiDimFit.mH120.root:"Theory Uncert.":4 higgsCombine.freezeall.MultiDimFit.mH120.root:"Statistical Only":2 --y-max 22 --y-cut 10 --breakdown "syst,theo,stat" --logo-sub Preliminary --translate ../breakdown.json --output plots/breakdown_H_201X
327354
```
328355

329356
### Running GoFs

‎plot.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -52,7 +52,7 @@ def lite_plot(
5252
scaleH=None,
5353
stack_style=0,
5454
style_set=style_set_A,
55-
prelim=False,
55+
prelim=True,
5656
format='png',
5757
fitDiag=None, # Add signal strength labels
5858
bkgUnc=None, # External bkg uncertainty

‎plotTF.py‎

Lines changed: 8 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -45,8 +45,8 @@ def base_plot(TF, fptpts, fmsdpts, frhopts, ax=None, rho=False):
4545
coff = 0.3
4646
zmin, zmax = np.floor(10 * np.min(TF)) / 10, np.ceil(10 * np.max(TF)) / 10
4747
zmin, zmax = zmin + 0.001, zmax - 0.001
48-
clim = np.round(np.min([abs(zmin - 1), abs(zmax - 1)]), 1)
49-
clim = np.clip(clim, 0.2, 0.4)
48+
clim = np.round(np.mean([abs(zmin - 1), abs(zmax - 1)]), 1)
49+
clim = np.clip(clim, 0.2, 0.5)
5050
levels = np.linspace(1 - clim, 1 + clim, 500)
5151

5252
if np.min(TF) < 1 - clim and np.max(TF) > 1 + clim:
@@ -65,6 +65,9 @@ def base_plot(TF, fptpts, fmsdpts, frhopts, ax=None, rho=False):
6565
c_extend = 'max'
6666
else:
6767
c_extend = 'neither'
68+
69+
print "X", np.min(TF), np.max(TF)
70+
print "lims:", zmin, zmax, clim, _extend, c_extend
6871

6972
pmesh = plt.contourf(frhopts if rho else fmsdpts, fptpts, TF, cmap='RdBu_r', levels=levels, extend=_extend, corner_mask=False)
7073
cax = hep.make_square_add_cbar(ax, pad=0.2, size=0.5)
@@ -255,7 +258,7 @@ def plotTF_ratio(in_ratio, mask, region, args=None, zrange=None):
255258

256259
ax = singleTF(tf_res, rho=args.rho)
257260
ax.set_title("Residual (Data/MC) TF", x=0, ha='left', fontsize='small')
258-
hep.cms.label(ax=ax, loc=2, year=args.year)
261+
hep.cms.label(ax=ax, loc=2, year=args.year, data=not args.isMC)
259262
ax.figure.savefig('{}/TF_data_{}.png'.format(args.output_folder, args.year), dpi=300, bbox_inches="tight")
260263
ax.figure.savefig('{}/TF_data_{}.pdf'.format(args.output_folder, args.year), transparent=True, bbox_inches="tight")
261264

@@ -272,13 +275,13 @@ def plotTF_ratio(in_ratio, mask, region, args=None, zrange=None):
272275

273276
ax = singleTF(tf_MC, rho=args.rho)
274277
ax.set_title("Tagger Response TF", x=0, ha='left', fontsize='small')
275-
hep.cms.label(ax=ax, loc=2, year=args.year)
278+
hep.cms.label(ax=ax, loc=2, year=args.year, data=not args.isMC)
276279
ax.figure.savefig('{}/TF_MC_{}.png'.format(args.output_folder, args.year), dpi=300, bbox_inches="tight")
277280
ax.figure.savefig('{}/TF_MC_{}.pdf'.format(args.output_folder, args.year), transparent=True, bbox_inches="tight")
278281

279282
ax = combinedTF(tf_MC, tf_res, rho=args.rho)
280283
ax.set_title("Effective Transfer Factor", x=0, ha='left', fontsize='small')
281-
hep.cms.label(ax=ax, loc=2, year=args.year)
284+
hep.cms.label(ax=ax, loc=2, year=args.year, data=not args.isMC)
282285
ax.figure.savefig('{}/TF_eff_{}.png'.format(args.output_folder, args.year), dpi=300, bbox_inches="tight")
283286
ax.figure.savefig('{}/TF_eff_{}.pdf'.format(args.output_folder, args.year), transparent=True, bbox_inches="tight")
284287
# except:

‎resampler.py‎

Lines changed: 135 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,135 @@
1+
#!/usr/bin/env python
2+
import argparse
3+
import numpy as np
4+
import tqdm
5+
import pickle
6+
import gzip
7+
8+
import ROOT
9+
ROOT.PyConfig.IgnoreCommandLineOptions = True
10+
ROOT.gROOT.SetBatch(True)
11+
ROOT.gEnv.SetValue("RooFit.Banner=0")
12+
ROOT.RooMsgService.instance().setGlobalKillBelow(ROOT.RooFit.WARNING)
13+
14+
15+
def _RooAbsCollection__iter__(self):
16+
it = self.iterator()
17+
obj = it.Next()
18+
while obj != None: # noqa: E711
19+
yield obj
20+
obj = it.Next()
21+
22+
23+
def _RooAbsCollection_assign(self, other):
24+
if self == other:
25+
return
26+
for el in self:
27+
if not hasattr(el, 'setVal'):
28+
continue
29+
theirs = other.find(el)
30+
if not theirs:
31+
continue
32+
el.setVal(theirs.getVal())
33+
el.setError(theirs.getError())
34+
el.setAsymError(theirs.getErrorLo(), theirs.getErrorHi())
35+
el.setAttribute("Constant", theirs.isConstant())
36+
37+
38+
def _RooArgList_fromiter(cls, iterable, silent=False):
39+
items = cls()
40+
for item in iterable:
41+
items.add(item, silent)
42+
return items
43+
44+
45+
ROOT.RooAbsCollection.assign = _RooAbsCollection_assign
46+
ROOT.RooAbsCollection.__iter__ = _RooAbsCollection__iter__
47+
ROOT.RooArgList.fromiter = classmethod(_RooArgList_fromiter)
48+
49+
50+
def to_numpy(hinput):
51+
if "<class 'ROOT.TH1" in str(type(hinput)):
52+
sumw = np.zeros(hinput.GetNbinsX())
53+
binning = np.zeros(sumw.size + 1)
54+
name = hinput.GetName()
55+
for i in range(1, sumw.size + 1):
56+
sumw[i-1] = hinput.GetBinContent(i)
57+
binning[i-1] = hinput.GetXaxis().GetBinLowEdge(i)
58+
binning[i] = hinput.GetXaxis().GetBinUpEdge(i)
59+
return (sumw, binning, name)
60+
else:
61+
raise TypeError("Don't know how to convert %r to numpy" % type(hinput))
62+
63+
64+
def resample_fit_result(args):
65+
fws, wsname = args.workspace.split(':')
66+
fin = ROOT.TFile.Open(fws)
67+
w = fin.Get(wsname)
68+
69+
ffit, fitname = args.fit.split(':')
70+
fin2 = ROOT.TFile.Open(ffit)
71+
if fitname != "prefit":
72+
fit = fin2.Get(fitname)
73+
else:
74+
fit = fin2.Get("fit_s")
75+
76+
model = w.obj(args.model)
77+
obs = model.GetObservables()[args.observable]
78+
pdf = model.GetPdf()
79+
params = pdf.getParameters(model.GetObservables())
80+
# It would seem necessary to do this but sometimes the fit
81+
# constant parameters are somehow the wrong value w.r.t. the model
82+
# params.assign(fit.constPars())
83+
if fitname != "prefit":
84+
params.assign(fit.floatParsFinal())
85+
else:
86+
params.assign(fit.floatParsInit())
87+
88+
all_processes = []
89+
cat = pdf.indexCat()
90+
for iCat in range(cat.numBins('')):
91+
cat.setBin(iCat)
92+
catpdf = pdf.getPdf(cat.getLabel())
93+
addpdf = [p for p in catpdf.pdfList() if p.dependsOn(model.GetObservables())]
94+
if len(addpdf) != 1:
95+
raise RuntimeError("I failed to parse the model structure :(")
96+
addpdf = addpdf[0]
97+
if not addpdf.dependsOn(obs):
98+
continue
99+
for proc, norm in zip(addpdf.pdfList(), addpdf.coefList()):
100+
all_processes.append((proc, norm))
101+
102+
nominal = {}
103+
for proc, norm in all_processes:
104+
h = proc.createHistogram("", obs)
105+
h.SetDirectory(0)
106+
h.Scale(norm.getVal())
107+
nominal[proc.GetName()] = to_numpy(h)
108+
del h
109+
110+
varied = {name: np.zeros((args.samples, ) + h[0].shape) for name, h in nominal.items()}
111+
for i in tqdm.trange(args.samples):
112+
params.assignValueOnly(fit.randomizePars())
113+
for proc, norm in all_processes:
114+
h = proc.createHistogram("", obs)
115+
h.Scale(norm.getVal())
116+
h.SetDirectory(0)
117+
varied[proc.GetName()][i] = to_numpy(h)[0]
118+
del h
119+
120+
return nominal, varied, args.samples
121+
122+
123+
if __name__ == '__main__':
124+
parser = argparse.ArgumentParser(description="Resample fit result and save as a pickled dictionary of numpy arrays for subsequent processing.")
125+
parser.add_argument("-w", "--workspace", metavar="ROOTFILE:WORKSPACE", help="Workspace to load (e.g. card.root:w)", required=True)
126+
parser.add_argument("-f", "--fit", metavar="ROOTFILE:FIT_NAME", help="Fit result to load (e.g. fitDiagnostics.root:fit_s)", required=True)
127+
parser.add_argument("-m", "--model", help="Model to load (default: %(default)s)", default="ModelConfig")
128+
parser.add_argument("-o", "--observable", help="Observable (default: %(default)s)", default="x")
129+
parser.add_argument("-n", "--samples", help="Number of samples (default: %(default)s)", type=int, default=200)
130+
parser.add_argument("--out", help="Output filename. Suggested suffix: .pkl.gz", type=str, required=True)
131+
132+
args = parser.parse_args()
133+
arrays = resample_fit_result(args)
134+
with gzip.open(args.out, 'wb') as fout:
135+
pickle.dump(arrays, fout)

0 commit comments

Comments
 (0)