diff --git a/dhi/plots/limits.py b/dhi/plots/limits.py index aae3ddcbc22bc3a21065cb79c63376a643d52adc..d2e797eb6ff9fb9fcc3a1efc369c2ee0f1988bb9 100644 --- a/dhi/plots/limits.py +++ b/dhi/plots/limits.py @@ -266,6 +266,7 @@ def plot_limit_scans( scan_parameter, names, expected_values, + 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. @@ -305,19 +308,26 @@ def plot_limit_scans( 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) + 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 @@ -355,36 +365,74 @@ 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_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)) - + 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}, - color=colors[col]) - draw_objs.append((g_exp, "SAME,CP" if show_points else "SAME,C")) - name = names[n_graphs - i - 1] + 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 else "L")) + "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, "{}, {}, expected".format(poi, name), 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, "{}, {}, observed".format(poi, name), + 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) diff --git a/dhi/tasks/limits.py b/dhi/tasks/limits.py index c0b680143024836e9f0361235c1ac6f3665b63cb..fcba284678a032f044e7b4f928746a1c3732c159 100644 --- a/dhi/tasks/limits.py +++ b/dhi/tasks/limits.py @@ -354,8 +354,6 @@ class PlotUpperLimits(UpperLimitsBase, POIPlotTask): class PlotMultipleUpperLimits(PlotUpperLimits, MultiDatacardTask): - unblinded = None - @classmethod def modify_param_values(cls, params): params = PlotUpperLimits.modify_param_values.__func__.__get__(cls)(params) @@ -444,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", @@ -457,6 +466,7 @@ class PlotMultipleUpperLimits(PlotUpperLimits, MultiDatacardTask): 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"), @@ -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 26c2622165c46edb50bdf98f861174671b2473f6..63fe295d94676e2685296003022842b8c1cd9ace 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"