From 47ff640649aee654a8ecb2de42c27d6525d5fe83 Mon Sep 17 00:00:00 2001 From: Luca Lavezzo Date: Sun, 1 Mar 2026 17:36:58 -0500 Subject: [PATCH 1/5] store group-nuisance mapping in fitresult; add support to plot group vars --- bin/rabbit_fit.py | 6 ++ bin/rabbit_plot_hists.py | 141 ++++++++++++++++++++++++++++++++++++++- 2 files changed, 144 insertions(+), 3 deletions(-) diff --git a/bin/rabbit_fit.py b/bin/rabbit_fit.py index 44893f1..7be75d4 100755 --- a/bin/rabbit_fit.py +++ b/bin/rabbit_fit.py @@ -487,6 +487,12 @@ def main(): "procs": ifitter.indata.procs, "pois": ifitter.poi_model.pois, "nois": ifitter.parms[ifitter.poi_model.npoi :][indata.noiidxs], + # Persist nuisance/group bookkeeping directly in fitresults + # so downstream plotting can reconstruct groupings without the input card. + "systs": ifitter.indata.systs, + "systgroups": ifitter.indata.systgroups, + "systgroupidxs": ifitter.indata.systgroupidxs, + "systsnoconstraint": ifitter.indata.systsnoconstraint, } with workspace.Workspace( diff --git a/bin/rabbit_plot_hists.py b/bin/rabbit_plot_hists.py index ce4e352..b35e58a 100755 --- a/bin/rabbit_plot_hists.py +++ b/bin/rabbit_plot_hists.py @@ -24,6 +24,38 @@ logger = None +def _decode_str(x): + return x.decode() if isinstance(x, (bytes, np.bytes_)) else str(x) + + +def _build_group_to_nuisance_map(meta): + # Prefer mapping persisted directly in fitresult meta, fallback to nested input meta. + systs = meta.get("systs") + groups = meta.get("systgroups") + groupidxs = meta.get("systgroupidxs") + + if systs is None or groups is None or groupidxs is None: + nested = meta.get("meta_info_input", {}) + systs = nested.get("systs") + groups = nested.get("systgroups") + groupidxs = nested.get("systgroupidxs") + + if systs is None or groups is None or groupidxs is None: + return {} + + systs = [_decode_str(x) for x in np.array(systs)] + groups = [_decode_str(x) for x in np.array(groups)] + + out = {} + for g, idxs in zip(groups, np.array(groupidxs, dtype=object)): + nuis = [] + for idx in np.array(idxs, dtype=int): + if 0 <= idx < len(systs): + nuis.append(systs[idx]) + out[g] = nuis + return out + + def parseArgs(): # choices for legend padding @@ -331,6 +363,30 @@ def parseArgs(): Additional varNames can be specified to add variations from the nominal input. """, ) + parser.add_argument( + "--varGroupNames", + type=str, + nargs="*", + default=None, + help=( + "Variation group names to build on-the-fly from nuisance variations " + "using fitresult meta mapping (systs/systgroups/systgroupidxs)." + ), + ) + parser.add_argument( + "--varGroupLabels", + type=str, + nargs="*", + default=None, + help="Label(s) for --varGroupNames.", + ) + parser.add_argument( + "--varGroupColors", + type=str, + nargs="*", + default=None, + help="Color(s) for --varGroupNames.", + ) parser.add_argument( "--varLabels", type=str, @@ -1174,6 +1230,9 @@ def make_plots( varFilesFitTypes=None, varMarkers=None, varNames=None, + varGroupNames=None, + varGroupLabels=None, + varGroupColors=None, varLabels=None, varColors=None, binwnorm=None, @@ -1214,13 +1273,16 @@ def make_plots( l if p not in args.suppressProcsLabel else None for l, p in zip(labels, procs) ] - if varNames is not None: + if varNames is not None or varGroupNames is not None: # take the first variations from the varFiles, empty if no varFiles are provided if len(varFilesFitTypes) == 1: varFilesFitTypes = varFilesFitTypes * len(varResults) hists_down = [] hists_up = [] + plot_var_names = [] if varNames is None else list(varNames) + plot_var_labels = [] if varLabels is None else list(varLabels) + plot_var_colors = [] if varColors is None else list(varColors) for r, t in zip(varResults, varFilesFitTypes): h = r[f"hist_{t}_inclusive"].get() @@ -1237,12 +1299,11 @@ def make_plots( hists_up.append(hist_up) # take the next variations from the nominal input file - if len(varNames) > len(varResults): + if varNames is not None and len(varNames) > len(varResults): # variations from the nominal input file hist_var = result[ f"hist_{fittype}_inclusive_variations{'_correlated' if args.correlatedVariations else ''}" ].get() - hists_down.extend( [ hist_var[{"downUpVar": 0, "vars": n}].project( @@ -1259,6 +1320,65 @@ def make_plots( for n in varNames[len(varResults) :] ] ) + + # grouped variations built on-the-fly from nuisance-level variations + if varGroupNames is not None and len(varGroupNames): + hist_var = result[ + f"hist_{fittype}_inclusive_variations{'_correlated' if args.correlatedVariations else ''}" + ].get() + axis_names = [a.name for a in axes] + var_axis_entries = {_decode_str(x) for x in hist_var.axes["vars"]} + group_to_nuis = _build_group_to_nuisance_map(kwopts.get("meta", {})) + h_nom = result[f"hist_{fittype}_inclusive"].get().project(*axis_names) + + for ig, gname in enumerate(varGroupNames): + members = group_to_nuis.get(gname, []) + members = [n for n in members if n in var_axis_entries] + if len(members) == 0: + logger.warning( + f"Group '{gname}' has no matching nuisances in variations axis; skipping." + ) + continue + + sigma2 = np.zeros_like(h_nom.values(), dtype=float) + for n in members: + hdn = hist_var[{"downUpVar": 0, "vars": n}].project(*axis_names) + hup = hist_var[{"downUpVar": 1, "vars": n}].project(*axis_names) + contrib = 0.5 * (hup.values() - hdn.values()) + sigma2 += contrib * contrib + + sigma = np.sqrt(sigma2) + h_up = h_nom.copy() + h_down = h_nom.copy() + h_up.values()[...] = h_nom.values() + sigma + h_down.values()[...] = h_nom.values() - sigma + hists_up.append(h_up) + hists_down.append(h_down) + + plot_var_names.append(gname) + if varGroupLabels is not None and ig < len(varGroupLabels): + plot_var_labels.append(varGroupLabels[ig]) + else: + plot_var_labels.append(gname) + if varGroupColors is not None and ig < len(varGroupColors): + plot_var_colors.append(varGroupColors[ig]) + + # Ensure metadata arrays align with the effective number of plotted variations. + n_total = len(hists_up) + if len(plot_var_names) < n_total: + plot_var_names.extend([f"var{i}" for i in range(len(plot_var_names), n_total)]) + if len(plot_var_labels) < n_total: + plot_var_labels.extend(plot_var_names[len(plot_var_labels) : n_total]) + if len(plot_var_colors) < n_total: + default_cols = [ + colormaps["tab10" if n_total < 10 else "tab20"](i) + for i in range(n_total) + ] + plot_var_colors.extend(default_cols[len(plot_var_colors) : n_total]) + + varNames = plot_var_names + varLabels = plot_var_labels + varColors = plot_var_colors else: hists_down = None hists_up = None @@ -1418,6 +1538,7 @@ def main(): varFiles = args.varFiles varNames = args.varNames + varGroupNames = args.varGroupNames varLabels = args.varLabels varColors = args.varColors if varNames is not None: @@ -1434,6 +1555,17 @@ def main(): colormaps["tab10" if len(varNames) < 10 else "tab20"](i) for i in range(len(varNames)) ] + varGroupLabels = args.varGroupLabels + varGroupColors = args.varGroupColors + if varGroupNames is not None: + if varGroupLabels is not None and len(varGroupLabels) != len(varGroupNames): + raise ValueError( + "Must specify the same number of args for --varGroupNames and --varGroupLabels" + ) + if varGroupColors is not None and len(varGroupColors) != len(varGroupNames): + raise ValueError( + "Must specify the same number of args for --varGroupNames and --varGroupColors" + ) fittype = "prefit" if args.prefit else "postfit" @@ -1470,6 +1602,9 @@ def main(): meta=meta, fittype=fittype, varNames=varNames, + varGroupNames=varGroupNames, + varGroupLabels=varGroupLabels, + varGroupColors=varGroupColors, varLabels=varLabels, varColors=varColors, varMarkers=args.varMarkers, From 388283ec7962fbb7e639c007cc50bb9009191748 Mon Sep 17 00:00:00 2001 From: Luca Lavezzo Date: Sun, 1 Mar 2026 17:43:16 -0500 Subject: [PATCH 2/5] lint --- bin/rabbit_plot_hists.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/bin/rabbit_plot_hists.py b/bin/rabbit_plot_hists.py index b35e58a..5481ec8 100755 --- a/bin/rabbit_plot_hists.py +++ b/bin/rabbit_plot_hists.py @@ -1366,7 +1366,9 @@ def make_plots( # Ensure metadata arrays align with the effective number of plotted variations. n_total = len(hists_up) if len(plot_var_names) < n_total: - plot_var_names.extend([f"var{i}" for i in range(len(plot_var_names), n_total)]) + plot_var_names.extend( + [f"var{i}" for i in range(len(plot_var_names), n_total)] + ) if len(plot_var_labels) < n_total: plot_var_labels.extend(plot_var_names[len(plot_var_labels) : n_total]) if len(plot_var_colors) < n_total: From afb2757eccecab53554b2628ab059b4835e1a33b Mon Sep 17 00:00:00 2001 From: Luca Lavezzo Date: Sun, 1 Mar 2026 18:11:27 -0500 Subject: [PATCH 3/5] plot using cov rather than the direct nuisances --- bin/rabbit_fit.py | 2 -- bin/rabbit_plot_hists.py | 68 +++++++++++----------------------------- 2 files changed, 18 insertions(+), 52 deletions(-) diff --git a/bin/rabbit_fit.py b/bin/rabbit_fit.py index 7be75d4..847e43c 100755 --- a/bin/rabbit_fit.py +++ b/bin/rabbit_fit.py @@ -487,8 +487,6 @@ def main(): "procs": ifitter.indata.procs, "pois": ifitter.poi_model.pois, "nois": ifitter.parms[ifitter.poi_model.npoi :][indata.noiidxs], - # Persist nuisance/group bookkeeping directly in fitresults - # so downstream plotting can reconstruct groupings without the input card. "systs": ifitter.indata.systs, "systgroups": ifitter.indata.systgroups, "systgroupidxs": ifitter.indata.systgroupidxs, diff --git a/bin/rabbit_plot_hists.py b/bin/rabbit_plot_hists.py index 5481ec8..ade432c 100755 --- a/bin/rabbit_plot_hists.py +++ b/bin/rabbit_plot_hists.py @@ -24,38 +24,6 @@ logger = None -def _decode_str(x): - return x.decode() if isinstance(x, (bytes, np.bytes_)) else str(x) - - -def _build_group_to_nuisance_map(meta): - # Prefer mapping persisted directly in fitresult meta, fallback to nested input meta. - systs = meta.get("systs") - groups = meta.get("systgroups") - groupidxs = meta.get("systgroupidxs") - - if systs is None or groups is None or groupidxs is None: - nested = meta.get("meta_info_input", {}) - systs = nested.get("systs") - groups = nested.get("systgroups") - groupidxs = nested.get("systgroupidxs") - - if systs is None or groups is None or groupidxs is None: - return {} - - systs = [_decode_str(x) for x in np.array(systs)] - groups = [_decode_str(x) for x in np.array(groups)] - - out = {} - for g, idxs in zip(groups, np.array(groupidxs, dtype=object)): - nuis = [] - for idx in np.array(idxs, dtype=int): - if 0 <= idx < len(systs): - nuis.append(systs[idx]) - out[g] = nuis - return out - - def parseArgs(): # choices for legend padding @@ -1321,33 +1289,33 @@ def make_plots( ] ) - # grouped variations built on-the-fly from nuisance-level variations + # grouped variations built on-the-fly from grouped global impacts + # (covariance-based decomposition). if varGroupNames is not None and len(varGroupNames): - hist_var = result[ - f"hist_{fittype}_inclusive_variations{'_correlated' if args.correlatedVariations else ''}" - ].get() axis_names = [a.name for a in axes] - var_axis_entries = {_decode_str(x) for x in hist_var.axes["vars"]} - group_to_nuis = _build_group_to_nuisance_map(kwopts.get("meta", {})) h_nom = result[f"hist_{fittype}_inclusive"].get().project(*axis_names) + name_impacts_grouped = f"hist_{fittype}_inclusive_global_impacts_grouped" + + if name_impacts_grouped not in result.keys(): + raise KeyError( + f"--varGroupNames requires '{name_impacts_grouped}' in the selected result/channel. " + "Run rabbit_fit with --computeHistImpacts to produce grouped impact histograms." + ) + h_impacts_grouped = result[name_impacts_grouped].get() + impact_entries = {_decode_str(x) for x in h_impacts_grouped.axes["impacts"]} for ig, gname in enumerate(varGroupNames): - members = group_to_nuis.get(gname, []) - members = [n for n in members if n in var_axis_entries] - if len(members) == 0: + if gname not in impact_entries: logger.warning( - f"Group '{gname}' has no matching nuisances in variations axis; skipping." + f"Group '{gname}' has no matching entry in grouped impacts axis; skipping." ) continue + sigma = ( + h_impacts_grouped[{"impacts": gname}] + .project(*axis_names) + .values() + ) - sigma2 = np.zeros_like(h_nom.values(), dtype=float) - for n in members: - hdn = hist_var[{"downUpVar": 0, "vars": n}].project(*axis_names) - hup = hist_var[{"downUpVar": 1, "vars": n}].project(*axis_names) - contrib = 0.5 * (hup.values() - hdn.values()) - sigma2 += contrib * contrib - - sigma = np.sqrt(sigma2) h_up = h_nom.copy() h_down = h_nom.copy() h_up.values()[...] = h_nom.values() + sigma From 45a785e8734563b4703b2d536f9ccf6753252def Mon Sep 17 00:00:00 2001 From: Luca Lavezzo Date: Sun, 1 Mar 2026 18:13:05 -0500 Subject: [PATCH 4/5] lint --- bin/rabbit_plot_hists.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/bin/rabbit_plot_hists.py b/bin/rabbit_plot_hists.py index ade432c..fb339ce 100755 --- a/bin/rabbit_plot_hists.py +++ b/bin/rabbit_plot_hists.py @@ -1311,9 +1311,7 @@ def make_plots( ) continue sigma = ( - h_impacts_grouped[{"impacts": gname}] - .project(*axis_names) - .values() + h_impacts_grouped[{"impacts": gname}].project(*axis_names).values() ) h_up = h_nom.copy() From 881e8babc32714818ead43bd0955e0796e45b35d Mon Sep 17 00:00:00 2001 From: Luca Lavezzo Date: Mon, 2 Mar 2026 15:13:47 -0500 Subject: [PATCH 5/5] CI be damned --- bin/rabbit_plot_hists.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bin/rabbit_plot_hists.py b/bin/rabbit_plot_hists.py index fb339ce..a3ae402 100755 --- a/bin/rabbit_plot_hists.py +++ b/bin/rabbit_plot_hists.py @@ -1302,7 +1302,7 @@ def make_plots( "Run rabbit_fit with --computeHistImpacts to produce grouped impact histograms." ) h_impacts_grouped = result[name_impacts_grouped].get() - impact_entries = {_decode_str(x) for x in h_impacts_grouped.axes["impacts"]} + impact_entries = {x for x in h_impacts_grouped.axes["impacts"]} for ig, gname in enumerate(varGroupNames): if gname not in impact_entries: