From 0cfb68130ec3d2017a38a1e39c55e9e84df5251e Mon Sep 17 00:00:00 2001 From: tolange 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 aae3ddc..78a04ef 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 86db70a..38f13cf 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 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 78a04ef..748f77b 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 "" 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 + "" 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 "" 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 0c83393..fcba284 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 26c2622..63fe295 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
---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"
-- GitLab From 6acacf1f6b8ccfc61d0e6ccc760a748062829b58 Mon Sep 17 00:00:00 2001 From: Marcel R 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 748f77b..d2e797e 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