From 0cfb68130ec3d2017a38a1e39c55e9e84df5251e Mon Sep 17 00:00:00 2001
From: tolange <torben.lange@cern.ch>
Date: Thu, 2 Sep 2021 18:12:45 +0200
Subject: [PATCH 1/3] Added unblinding feature to PlotMultipleUpperLimits
---
dhi/plots/limits.py | 49 ++++++++++++++++++++++++++++++++++++---------
dhi/tasks/limits.py | 2 +-
2 files changed, 41 insertions(+), 10 deletions(-)
diff --git a/dhi/plots/limits.py b/dhi/plots/limits.py
index aae3ddcb..78a04ef6 100644
--- a/dhi/plots/limits.py
+++ b/dhi/plots/limits.py
@@ -23,7 +23,6 @@ from dhi.plots.util import (
get_graph_points, get_y_range, get_contours, fill_hist_from_points, infer_binning_from_grid,
)
-
colors = colors.root
@@ -266,6 +265,7 @@ def plot_limit_scans(
scan_parameter,
names,
expected_values,
+ unblinded = False,
theory_values=None,
y_log=False,
x_min=None,
@@ -299,6 +299,9 @@ def plot_limit_scans(
When *show_points* is *True*, the central scan points are drawn on top of the interpolated
curve. When *paper* is *True*, certain plot configurations are adjusted for use in publications.
+ If unblinded is *True*, *expected_values* should also contain the observed limits,
+ so that both the expected and observed limits can be drawn.
+
Example: https://cms-hh.web.cern.ch/tools/inference/tasks/limits.html#multiple-limits-on-poi-vs-scan-parameter
"""
import plotlib.root as r
@@ -311,13 +314,14 @@ def plot_limit_scans(
_ev = {key: _ev[key] for key in _ev.dtype.names}
_expected_values.append(_ev)
expected_values = _expected_values
-
# input checks
n_graphs = len(expected_values)
assert n_graphs >= 1
assert len(names) == n_graphs
assert all(scan_parameter in ev for ev in expected_values)
assert all("limit" in ev for ev in expected_values)
+ if unblinded:
+ assert all("observed" in ev for ev in expected_values)
scan_values = expected_values[0][scan_parameter]
has_thy = theory_values is not None
has_thy_err = False
@@ -352,10 +356,17 @@ def plot_limit_scans(
r.setup_hist(h_dummy, pad=pad, props={"LineWidth": 0})
draw_objs.append((h_dummy, "HIST"))
+ g_exp_dummy = None
+ g_obs_dummy = None
# central values
for i, (ev, col, ms) in enumerate(zip(expected_values[::-1], color_sequence[:n_graphs][::-1],
marker_sequence[:n_graphs][::-1])):
mask = ~np.isnan(ev["limit"])
+ limit_obs_values = None
+ if unblinded:
+ mask_obs = ~np.isnan(ev["observed"])
+ mask = np.logical_and(mask,mask_obs)
+ limit_obs_values = ev["observed"][mask]
limit_values = ev["limit"][mask]
scan_values = ev[scan_parameter][mask]
n_nans = (~mask).sum()
@@ -363,13 +374,27 @@ def plot_limit_scans(
print("WARNING: found {} NaN(s) in limit values at index {}".format(n_nans,
len(expected_values) - 1 - i))
- g_exp = create_tgraph(mask.sum(), scan_values, limit_values)
- r.setup_graph(g_exp, props={"LineWidth": 2, "MarkerStyle": ms, "MarkerSize": 1.2},
- color=colors[col])
- draw_objs.append((g_exp, "SAME,CP" if show_points else "SAME,C"))
name = names[n_graphs - i - 1]
- legend_entries.insert(0, (g_exp, to_root_latex(br_hh_names.get(name, name)),
- "LP" if show_points else "L"))
+ if unblinded:
+ g_exp = create_tgraph(mask.sum(), scan_values, limit_values)
+ g_obs = create_tgraph(mask.sum(), scan_values, limit_obs_values)
+ r.setup_graph(g_obs, props={"LineWidth": 2, "MarkerStyle": ms, "MarkerSize": 1.2},
+ color=colors[col])
+ r.setup_graph(g_exp, props={"LineWidth": 2, "MarkerStyle": ms, "MarkerSize": 1.2, "LineStyle": 2},
+ color=colors[col])
+ draw_objs.append((g_exp, "SAME,C"))
+ draw_objs.append((g_obs, "SAME,CP" if show_points else "SAME,C"))
+ legend_entries.insert(0, (g_obs, to_root_latex(br_hh_names.get(name, name)),
+ "LP" if show_points else "L"))
+ g_exp_dummy = g_exp.Clone()
+ g_obs_dummy = g_obs.Clone()
+ else:
+ g_exp = create_tgraph(mask.sum(), scan_values, limit_values)
+ r.setup_graph(g_exp, props={"LineWidth": 2, "MarkerStyle": ms, "MarkerSize": 1.2},
+ color=colors[col])
+ draw_objs.append((g_exp, "SAME,CP" if show_points else "SAME,C"))
+ legend_entries.insert(0, (g_exp, to_root_latex(br_hh_names.get(name, name)),
+ "LP" if show_points else "L"))
y_max_value = max(y_max_value, max(limit_values))
y_min_value = min(y_min_value, min(limit_values))
# print excluded ranges
@@ -388,7 +413,13 @@ def plot_limit_scans(
y_min, y_max, _ = get_y_range(y_min_value, y_max_value, y_min, y_max, log=y_log)
h_dummy.SetMinimum(y_min)
h_dummy.SetMaximum(y_max)
-
+ if unblinded:
+ r.apply_properties(g_exp_dummy, {"LineColor": colors.black})
+ r.apply_properties(g_obs_dummy, {"LineColor": colors.black})
+ legend_entries.insert(0, (g_obs_dummy, to_root_latex("observed"),
+ "L"))
+ legend_entries.insert(0, (g_exp_dummy, to_root_latex("expected"),
+ "L"))
# draw the theory prediction
if has_thy:
scan_values_thy = theory_values[scan_parameter]
diff --git a/dhi/tasks/limits.py b/dhi/tasks/limits.py
index 86db70a9..38f13cfd 100644
--- a/dhi/tasks/limits.py
+++ b/dhi/tasks/limits.py
@@ -328,7 +328,6 @@ class PlotUpperLimits(UpperLimitsBase, POIPlotTask):
class PlotMultipleUpperLimits(PlotUpperLimits, MultiDatacardTask):
- unblinded = None
@classmethod
def modify_param_values(cls, params):
@@ -431,6 +430,7 @@ class PlotMultipleUpperLimits(PlotUpperLimits, MultiDatacardTask):
scan_parameter=self.scan_parameter,
names=names,
expected_values=limit_values,
+ unblinded = self.unblinded,
theory_values=thy_values,
x_min=self.get_axis_limit("x_min"),
x_max=self.get_axis_limit("x_max"),
--
GitLab
From 2befe7eaf1fd138989046f2898f587bee171b196 Mon Sep 17 00:00:00 2001
From: Marcel R <github.riga@icloud.com>
Date: Tue, 7 Sep 2021 20:50:35 +0200
Subject: [PATCH 2/3] Align multi limit plot functions.
---
dhi/plots/limits.py | 127 ++++++++++--------
dhi/tasks/limits.py | 34 ++++-
.../plotmultipleupperlimits_param_tab.md | 2 +-
3 files changed, 101 insertions(+), 62 deletions(-)
diff --git a/dhi/plots/limits.py b/dhi/plots/limits.py
index 78a04ef6..748f77b7 100644
--- a/dhi/plots/limits.py
+++ b/dhi/plots/limits.py
@@ -23,6 +23,7 @@ from dhi.plots.util import (
get_graph_points, get_y_range, get_contours, fill_hist_from_points, infer_binning_from_grid,
)
+
colors = colors.root
@@ -265,7 +266,7 @@ def plot_limit_scans(
scan_parameter,
names,
expected_values,
- unblinded = False,
+ observed_values=None,
theory_values=None,
y_log=False,
x_min=None,
@@ -283,7 +284,9 @@ def plot_limit_scans(
Creates a plot showing multiple upper limit scans of a *poi* over a *scan_parameter* and saves
it at *paths*. *expected_values* should be a list of mappings to lists of values or a record
array with keys "<scan_parameter>" and "limit". Each mapping in *expected_values* will result in
- a different curve. When *theory_values* is set, it should have a similar format with keys
+ a different curve. When *observed_values* is set, it should have a similar format with keys
+ "<scan_parameter>" and "limit", where the i-th element is corresponding to the i-th element in
+ *expected_values*. When *theory_values* is set, it should have a similar format with keys
"<scan_parameter>" and "xsec", and optionally "xsec_p1" and "xsec_m1". *names* denote the names
of limit curves shown in the legend. When a name is found to be in dhi.config.br_hh_names, its
value is used as a label instead.
@@ -299,29 +302,32 @@ def plot_limit_scans(
When *show_points* is *True*, the central scan points are drawn on top of the interpolated
curve. When *paper* is *True*, certain plot configurations are adjusted for use in publications.
- If unblinded is *True*, *expected_values* should also contain the observed limits,
- so that both the expected and observed limits can be drawn.
-
Example: https://cms-hh.web.cern.ch/tools/inference/tasks/limits.html#multiple-limits-on-poi-vs-scan-parameter
"""
import plotlib.root as r
ROOT = import_ROOT()
# convert record arrays to dicts mapping to arrays
- _expected_values = []
- for _ev in expected_values:
- if isinstance(_ev, np.ndarray):
- _ev = {key: _ev[key] for key in _ev.dtype.names}
- _expected_values.append(_ev)
- expected_values = _expected_values
+ def check_values(values):
+ _values = []
+ for v in values:
+ if isinstance(v, np.ndarray):
+ v = {key: v[key] for key in v.dtype.names}
+ assert "limit" in v
+ assert scan_parameter in v
+ _values.append(v)
+ return _values
+
# input checks
+ expected_values = check_values(expected_values)
n_graphs = len(expected_values)
assert n_graphs >= 1
assert len(names) == n_graphs
- assert all(scan_parameter in ev for ev in expected_values)
- assert all("limit" in ev for ev in expected_values)
- if unblinded:
- assert all("observed" in ev for ev in expected_values)
+ has_obs = False
+ if observed_values:
+ assert len(observed_values) == n_graphs
+ observed_values = check_values(observed_values)
+ has_obs = True
scan_values = expected_values[0][scan_parameter]
has_thy = theory_values is not None
has_thy_err = False
@@ -361,65 +367,78 @@ def plot_limit_scans(
# central values
for i, (ev, col, ms) in enumerate(zip(expected_values[::-1], color_sequence[:n_graphs][::-1],
marker_sequence[:n_graphs][::-1])):
+ name = names[n_graphs - i - 1]
+
+ # expected graph
mask = ~np.isnan(ev["limit"])
- limit_obs_values = None
- if unblinded:
- mask_obs = ~np.isnan(ev["observed"])
- mask = np.logical_and(mask,mask_obs)
- limit_obs_values = ev["observed"][mask]
limit_values = ev["limit"][mask]
scan_values = ev[scan_parameter][mask]
n_nans = (~mask).sum()
if n_nans:
- print("WARNING: found {} NaN(s) in limit values at index {}".format(n_nans,
- len(expected_values) - 1 - i))
-
- name = names[n_graphs - i - 1]
- if unblinded:
- g_exp = create_tgraph(mask.sum(), scan_values, limit_values)
- g_obs = create_tgraph(mask.sum(), scan_values, limit_obs_values)
- r.setup_graph(g_obs, props={"LineWidth": 2, "MarkerStyle": ms, "MarkerSize": 1.2},
- color=colors[col])
- r.setup_graph(g_exp, props={"LineWidth": 2, "MarkerStyle": ms, "MarkerSize": 1.2, "LineStyle": 2},
- color=colors[col])
- draw_objs.append((g_exp, "SAME,C"))
- draw_objs.append((g_obs, "SAME,CP" if show_points else "SAME,C"))
- legend_entries.insert(0, (g_obs, to_root_latex(br_hh_names.get(name, name)),
- "LP" if show_points else "L"))
- g_exp_dummy = g_exp.Clone()
- g_obs_dummy = g_obs.Clone()
- else:
- g_exp = create_tgraph(mask.sum(), scan_values, limit_values)
- r.setup_graph(g_exp, props={"LineWidth": 2, "MarkerStyle": ms, "MarkerSize": 1.2},
- color=colors[col])
- draw_objs.append((g_exp, "SAME,CP" if show_points else "SAME,C"))
- legend_entries.insert(0, (g_exp, to_root_latex(br_hh_names.get(name, name)),
- "LP" if show_points else "L"))
+ print("WARNING: found {} NaN(s) in expected limit values at index {}".format(n_nans,
+ n_graphs - 1 - i))
+ g_exp = create_tgraph(mask.sum(), scan_values, limit_values)
+ r.setup_graph(g_exp, props={"LineWidth": 2, "MarkerStyle": ms, "MarkerSize": 1.2,
+ "LineStyle": 2 if has_obs else 1}, color=colors[col])
+ draw_objs.append((g_exp, "SAME,CP" if show_points and not has_obs else "SAME,C"))
+ legend_entries.insert(0, (g_exp, to_root_latex(br_hh_names.get(name, name)),
+ "LP" if show_points and not has_obs else "L"))
y_max_value = max(y_max_value, max(limit_values))
y_min_value = min(y_min_value, min(limit_values))
- # print excluded ranges
- print_excluded_ranges(scan_parameter, poi + " " + name,
+
+ # print expected excluded ranges
+ print_excluded_ranges(scan_parameter, poi + " " + name + " (expected)",
scan_values,
limit_values,
theory_values[scan_parameter] if has_thy else None,
theory_values["xsec"] if has_thy else None,
)
+ # observed graph
+ if has_obs:
+ ov = observed_values[n_graphs - i - 1]
+ obs_mask = ~np.isnan(ov["limit"])
+ obs_limit_values = ov["limit"][obs_mask]
+ obs_scan_values = ov[scan_parameter][obs_mask]
+ n_nans = (~obs_mask).sum()
+ if n_nans:
+ print("WARNING: found {} NaN(s) in observed limit values at index {}".format(n_nans,
+ n_graphs - 1 - i))
+ g_obs = create_tgraph(obs_mask.sum(), obs_scan_values, obs_limit_values)
+ r.setup_graph(g_obs, props={"LineWidth": 2, "MarkerStyle": ms, "MarkerSize": 1.2},
+ color=colors[col])
+ draw_objs.append((g_obs, "SAME,CP" if show_points else "SAME,C"))
+ legend_entries[0] = (g_obs, to_root_latex(br_hh_names.get(name, name)),
+ "LP" if show_points else "L")
+ y_max_value = max(y_max_value, max(obs_limit_values))
+ y_min_value = min(y_min_value, min(obs_limit_values))
+
+ # print observed excluded ranges
+ print_excluded_ranges(scan_parameter, poi + " " + name + " (observed)",
+ obs_scan_values,
+ obs_limit_values,
+ theory_values[scan_parameter] if has_thy else None,
+ theory_values["xsec"] if has_thy else None,
+ )
+
+ # add additional legend entries to distinguish expected and observed lines
+ if has_obs:
+ g_exp_dummy = g_exp.Clone()
+ g_obs_dummy = g_obs.Clone()
+ r.apply_properties(g_exp_dummy, {"LineColor": colors.black})
+ r.apply_properties(g_obs_dummy, {"LineColor": colors.black})
+ legend_entries.append((g_exp_dummy, "expected", "L"))
+ legend_entries.append((g_obs_dummy, "observed", "L"))
+
# get theory prediction limits
if has_thy:
y_min_value = min(y_min_value, min(theory_values["xsec_m1" if has_thy_err else "xsec"]))
- # set limits
+ # set axis limits
y_min, y_max, _ = get_y_range(y_min_value, y_max_value, y_min, y_max, log=y_log)
h_dummy.SetMinimum(y_min)
h_dummy.SetMaximum(y_max)
- if unblinded:
- r.apply_properties(g_exp_dummy, {"LineColor": colors.black})
- r.apply_properties(g_obs_dummy, {"LineColor": colors.black})
- legend_entries.insert(0, (g_obs_dummy, to_root_latex("observed"),
- "L"))
- legend_entries.insert(0, (g_exp_dummy, to_root_latex("expected"),
- "L"))
+
# draw the theory prediction
if has_thy:
scan_values_thy = theory_values[scan_parameter]
diff --git a/dhi/tasks/limits.py b/dhi/tasks/limits.py
index 0c833930..fcba2846 100644
--- a/dhi/tasks/limits.py
+++ b/dhi/tasks/limits.py
@@ -354,7 +354,6 @@ class PlotUpperLimits(UpperLimitsBase, POIPlotTask):
class PlotMultipleUpperLimits(PlotUpperLimits, MultiDatacardTask):
-
@classmethod
def modify_param_values(cls, params):
params = PlotUpperLimits.modify_param_values.__func__.__get__(cls)(params)
@@ -443,11 +442,22 @@ class PlotMultipleUpperLimits(PlotUpperLimits, MultiDatacardTask):
if self.datacard_names:
names = list(self.datacard_names)
- # reoder if requested
+ # reorder if requested
if self.datacard_order:
limit_values = [limit_values[i] for i in self.datacard_order]
names = [names[i] for i in self.datacard_order]
+ # prepare observed values
+ obs_values = None
+ if self.unblinded:
+ obs_values = [
+ {
+ self.scan_parameter: _limit_values[self.scan_parameter],
+ "limit": _limit_values["observed"],
+ }
+ for _limit_values in limit_values
+ ]
+
# call the plot function
self.call_plot_func(
"dhi.plots.limits.plot_limit_scans",
@@ -456,7 +466,7 @@ class PlotMultipleUpperLimits(PlotUpperLimits, MultiDatacardTask):
scan_parameter=self.scan_parameter,
names=names,
expected_values=limit_values,
- unblinded = self.unblinded,
+ observed_values=obs_values,
theory_values=thy_values,
x_min=self.get_axis_limit("x_min"),
x_max=self.get_axis_limit("x_max"),
@@ -474,8 +484,6 @@ class PlotMultipleUpperLimits(PlotUpperLimits, MultiDatacardTask):
class PlotMultipleUpperLimitsByModel(PlotUpperLimits, MultiHHModelTask):
- unblinded = None
-
allow_empty_hh_model = True
def requires(self):
@@ -568,11 +576,22 @@ class PlotMultipleUpperLimitsByModel(PlotUpperLimits, MultiHHModelTask):
if self.hh_model_names:
names = list(self.hh_model_names)
- # reoder if requested
+ # reorder if requested
if self.hh_model_order:
limit_values = [limit_values[i] for i in self.hh_model_order]
names = [names[i] for i in self.hh_model_order]
+ # prepare observed values
+ obs_values = None
+ if self.unblinded:
+ obs_values = [
+ {
+ self.scan_parameter: _limit_values[self.scan_parameter],
+ "limit": _limit_values["observed"],
+ }
+ for _limit_values in limit_values
+ ]
+
# call the plot function
self.call_plot_func(
"dhi.plots.limits.plot_limit_scans",
@@ -581,6 +600,7 @@ class PlotMultipleUpperLimitsByModel(PlotUpperLimits, MultiHHModelTask):
scan_parameter=self.scan_parameter,
names=names,
expected_values=limit_values,
+ observed_values=obs_values,
theory_values=thy_values,
x_min=self.get_axis_limit("x_min"),
x_max=self.get_axis_limit("x_max"),
@@ -829,7 +849,7 @@ class PlotUpperLimitsAtPoint(POIPlotTask, SnapshotUser, MultiDatacardTask, BoxPl
for d, label in zip(data, self.extra_labels):
d["label"] = label
- # reoder if requested
+ # reorder if requested
if self.datacard_order:
data = [data[i] for i in self.datacard_order]
diff --git a/docs/content/snippets/plotmultipleupperlimits_param_tab.md b/docs/content/snippets/plotmultipleupperlimits_param_tab.md
index 26c26221..63fe295d 100644
--- a/docs/content/snippets/plotmultipleupperlimits_param_tab.md
+++ b/docs/content/snippets/plotmultipleupperlimits_param_tab.md
@@ -2,6 +2,6 @@ The `PlotMultipleUpperLimits` task collects the fit results from multiple `Merge
<div class="dhi_parameter_table">
---8<-- "content/snippets/parameters.md@-2,21,19,16,54,14,66,12,17,18,74,10,11,3,4,22,23,5,6,7,8,9"
+--8<-- "content/snippets/parameters.md@-2,21,19,16,54,14,66,12,48,17,18,74,10,11,3,4,22,23,5,6,7,8,9"
</div>
--
GitLab
From 6acacf1f6b8ccfc61d0e6ccc760a748062829b58 Mon Sep 17 00:00:00 2001
From: Marcel R <github.riga@icloud.com>
Date: Tue, 7 Sep 2021 20:52:43 +0200
Subject: [PATCH 3/3] Align multi limit plot functions.
---
dhi/plots/limits.py | 6 ++----
1 file changed, 2 insertions(+), 4 deletions(-)
diff --git a/dhi/plots/limits.py b/dhi/plots/limits.py
index 748f77b7..d2e797eb 100644
--- a/dhi/plots/limits.py
+++ b/dhi/plots/limits.py
@@ -362,8 +362,6 @@ def plot_limit_scans(
r.setup_hist(h_dummy, pad=pad, props={"LineWidth": 0})
draw_objs.append((h_dummy, "HIST"))
- g_exp_dummy = None
- g_obs_dummy = None
# central values
for i, (ev, col, ms) in enumerate(zip(expected_values[::-1], color_sequence[:n_graphs][::-1],
marker_sequence[:n_graphs][::-1])):
@@ -387,7 +385,7 @@ def plot_limit_scans(
y_min_value = min(y_min_value, min(limit_values))
# print expected excluded ranges
- print_excluded_ranges(scan_parameter, poi + " " + name + " (expected)",
+ print_excluded_ranges(scan_parameter, "{}, {}, expected".format(poi, name),
scan_values,
limit_values,
theory_values[scan_parameter] if has_thy else None,
@@ -414,7 +412,7 @@ def plot_limit_scans(
y_min_value = min(y_min_value, min(obs_limit_values))
# print observed excluded ranges
- print_excluded_ranges(scan_parameter, poi + " " + name + " (observed)",
+ print_excluded_ranges(scan_parameter, "{}, {}, observed".format(poi, name),
obs_scan_values,
obs_limit_values,
theory_values[scan_parameter] if has_thy else None,
--
GitLab