From a408f0ede8ce2b47ba8b534cc2d36b779ff9a27a Mon Sep 17 00:00:00 2001 From: Devin Date: Sat, 26 Sep 2026 00:13:38 +0200 Subject: [PATCH 01/12] Read a category from a per-model-point input directory The HH->bbWW single-lepton shapes hold, in the file for mass MX, a signal region for every mass (SR_SL_M300 ... SR_SL_M1000); only SR_SL_M${MX} belongs to the MX limit. Listing those directories as categories puts all ten into every card as separate bins, and nothing reports it, since every file has all of them. A configuration can now map a category to its input directory, with the model parameters substituted at read time: input_categories: SR_SL/res2b: SR_SL_M${MX}/res2b The category keeps one name for every mass -- bin names, per-category card directories, Process/Uncertainty scoping -- and only the read is redirected, in the maker's four histogram-path builders. A map entry for a category that is not configured is an error; configurations without a map read exactly as before. This covers cards that read the merged shapes directly. rebin_2d.py does not apply the map, so a card with a `preprocess:` step still needs its input under the category name. Checked on the real SL shapes (v2605a_TwoStageDNN_Shapes_v3, MX=500, Run3_Early, no preprocessing): the card has the six bins {e,mu} x SR_SL/{res2b,boosted,recovery}, and TT and signal in every bin equal the sum of the four sub-eras' SR_SL_M500 inputs. Co-Authored-By: Claude Opus 5.5 (1M context) --- common/tools.py | 52 ++++++++++++++++++++++++++++++++++++++++++++++++ dc_make/maker.py | 28 +++++++++++++++++++++----- 2 files changed, 75 insertions(+), 5 deletions(-) diff --git a/common/tools.py b/common/tools.py index 59aeee7..e0665e5 100644 --- a/common/tools.py +++ b/common/tools.py @@ -101,6 +101,58 @@ def base(self, category): return self.split(category)[0] +class InputCategories: + """Where a category's shapes live in the input file, when that is not its own name. + + A datacard category is one name for every model point: it names the bin, the + per-category datacard directory, and what Process/Uncertainty scoping matches + against. Some inputs do not store it under one name. The HH->bbWW single-lepton + shapes hold, in the file for mass MX, a region for every mass (SR_SL_M300 ... + SR_SL_M1000), and only SR_SL_M${MX} is the one the MX hypothesis selected. Listing + those as categories would put all ten into every datacard as separate bins -- every + file has all of them, so nothing would ever report one missing. + + So the configuration maps a category to its input directory, with the model's + parameters substituted at read time: + + input_categories: + SR_SL/res2b: SR_SL_M${MX}/res2b + + Only exact names are mapped. A sliced name (SR_SL/res2b_hmebox0) is read from shapes + a preprocessing step has already written under the category's own name, so it + passes through, as does every category of a configuration that declares no map. + """ + + def __init__(self, mapping=None): + mapping = mapping or {} + if not isinstance(mapping, dict): + raise RuntimeError( + f"input_categories must map category -> input directory, got {mapping!r}" + ) + self.mapping = {str(k): str(v) for k, v in mapping.items()} + + @classmethod + def fromConfig(cls, cfg): + return cls(cfg.get("input_categories")) + + def unknown(self, categories): + """Mapped names that are not among the configuration's categories -- a typo + there would otherwise leave the category read from its own name, silently.""" + return sorted(set(self.mapping) - set(categories)) + + def path(self, category, params): + """The directory holding `category` for model point `params` ({"MX": 500}).""" + pattern = self.mapping.get(category, category) + for name, value in (params or {}).items(): + pattern = pattern.replace("${" + name + "}", str(value)) + if "${" in pattern: + raise RuntimeError( + f"input category '{pattern}' for '{category}' still has an unresolved " + f"parameter after substituting {dict(params or {})}" + ) + return pattern + + class PackageWrapper: def __init__(self, import_fn): self._package = None diff --git a/dc_make/maker.py b/dc_make/maker.py index abfd482..d25c87e 100644 --- a/dc_make/maker.py +++ b/dc_make/maker.py @@ -12,6 +12,7 @@ resolveNegativeBins, getRelevantBins, CategoryNaming, + InputCategories, ) from .process import Process from .uncertainty import ( @@ -68,6 +69,15 @@ def __init__( # group the slices of one base category. self.naming = CategoryNaming.fromConfig(cfg) self.categories = list(cfg["categories"]) + # Where each category's shapes sit in the input, when that depends on the model + # point; see InputCategories. + self.input_categories = InputCategories.fromConfig(cfg) + unknown = self.input_categories.unknown(self.categories) + if unknown: + raise RuntimeError( + f"input_categories maps {unknown}, which are not configured categories " + f"({self.categories})" + ) self.signalFractionForRelevantBins = cfg["signalFractionForRelevantBins"] self.era_groups = cfg.get("era_groups", {}) @@ -289,6 +299,11 @@ def PPECC(self): self.processes.keys(), param_bins, self.eras, self.channels, self.categories ) + def histDir(self, channel, category, model_params): + """The input directory holding one bin's histograms. The category is the bin's + own name unless `input_categories` places it per model point.""" + return f"{channel}/{self.input_categories.path(category, model_params)}" + def getInputFile(self, era, model_params): file_name = self.model.getInputFileName(era, model_params) if file_name not in self.input_files: @@ -349,7 +364,7 @@ def getMVLnNValue(self, unc, process, era, channel, category, model_params): channel, category, model_params, - f"{channel}/{category}/{subp}", + f"{self.histDir(channel, category, model_params)}/{subp}", ) sub_yield = hist.Integral(1, hist.GetNbinsX() + 1) if sub_yield <= 0: @@ -444,7 +459,7 @@ def _getSubEraLnNVariedShapes( channel, category, model_params, - f"{channel}/{category}/{hist_name_suffix}", + f"{self.histDir(channel, category, model_params)}/{hist_name_suffix}", ) unc_value = self._getLnNValue( unc, process, proc_name_for_unc, sub_era, channel, category @@ -728,11 +743,12 @@ def _rawShape( "this bin." ) else: - hist_name = f"{channel}/{category}/{process.hist_name}" + hist_dir = self.histDir(channel, category, model_params) + hist_name = f"{hist_dir}/{process.hist_name}" hists = [] if process.subprocesses: for subp in process.subprocesses: - hist_name = f"{channel}/{category}/{subp}" + hist_name = f"{hist_dir}/{subp}" if unc_name and unc_scale: hist_name += f"_{unc_name}_{unc_scale}" subhist = self.readHist(file, hist_name) @@ -972,7 +988,9 @@ def hasNominalShape(self, process, era, channel, category): MC -- e.g. a standalone single-era limit for a sparse category/channel that only has signal statistics once combined with other eras.""" sub_eras = self.getSubEras(era) if self.isMetaEra(era) else [era] - hist_name = f"{channel}/{category}/{process.hist_name}" + hist_name = ( + f"{self.histDir(channel, category, process.params)}/{process.hist_name}" + ) for sub_era in sub_eras: _, file = self.getInputFile(sub_era, process.params) obj = self.readHist(file, hist_name) From c36b236b4c4a6f5a3be561ca03d42ba81625e38d Mon Sep 17 00:00:00 2001 From: Devin Date: Tue, 8 Sep 2026 22:35:16 +0200 Subject: [PATCH 02/12] Gate every background's effective entries per bin, and never exempt positivity Two changes to what a candidate binning has to satisfy, both about backgrounds that are present but not measured. min_bin_bkg_each_neff requires every non-negligible background to have at least this many effective MC entries in every bin. min_bin_bkg_neff already existed but constrains only the *summed* background, so it is satisfied by any one well-measured process -- TT alone was covering for DY at N_eff 0.1. The sliced axis has carried the per-process arm for a while (min_slice_bkg_each_neff); the mass axis, where essentially all fit bins live, had no per-process test at all. Failing bins are merged into their lower-signal neighbour, so the resolution given up is where the signal is thinnest. Defaults to 0.0, which is the previous behaviour exactly. Positivity is no longer exemptable, in either _bin_passes or _slice_passes. `exempt` exists so that a negligible process cannot veto every candidate split -- it is about the *magnitude* floor, not about sign. Letting an exempt process go negative is not merely an ugly bin: resolveNegativeBins in dc_make/maker.py rejects a negative bin holding >= signalFractionForRelevantBins of the signal unconditionally, ignoring allow_negative_bins_within_error, so the datacard fails to build at all. Observed at m1000 muMu/SR/recovery_dnn1, where ST was 4.24% of the slice -- just under min_bkg_frac 0.05, hence exempt -- and came out at -0.567 +- 0.789, killing the production. Measured on HH->bbWW Run3_Early: with n_slices 4 and bkg_per_bin 1.0, turning the new gate on at 2.0 gives a geometric-mean expected limit 0.988x the production binning -- better on the limit *and* on every validity metric at once. Per-process bins below N_eff 1 fall 1.5% -> 0.2%, below N_eff 2 fall 3.2% -> 1.1%, the median DY N_eff nearly doubles 16.3 -> 31.8, and DY bins below N_eff 4 fall 29.2% -> 18.2%, for 3.7% of the fit bins (6220 -> 5989). Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01We5E9ZkM9tanReeRE2WSSR --- bin_opt_2d/rebin_2d.py | 71 +++++++++++++++++++++++++++++++++++------- 1 file changed, 60 insertions(+), 11 deletions(-) diff --git a/bin_opt_2d/rebin_2d.py b/bin_opt_2d/rebin_2d.py index 9132a16..41ad31b 100644 --- a/bin_opt_2d/rebin_2d.py +++ b/bin_opt_2d/rebin_2d.py @@ -46,6 +46,7 @@ "min_slice_bkg_each_neff": 0.0, "min_bin_bkg_each": 0.01, "min_bin_bkg_neff": 4.0, + "min_bin_bkg_each_neff": 0.0, "bkg_per_bin": 5.0, "min_bkg_frac": 0.05, "min_signal": 0.5, @@ -336,15 +337,17 @@ def _slice_passes( if min_neff > 0 and total_error is not None: if effective_entries(total, total_error) < min_neff: return False - if min_each > 0 or min_proc_neff > 0: - for name, value in yields.items(): - if name in exempt: - continue - if value <= min_each: + for name, value in yields.items(): + # As in _bin_passes: exemption covers the magnitude floor, never the sign. + if value < 0: + return False + if name in exempt: + continue + if min_each > 0 and value <= min_each: + return False + if min_proc_neff > 0 and errors is not None: + if effective_entries(value, errors.get(name, 0.0)) < min_proc_neff: return False - if min_proc_neff > 0 and errors is not None: - if effective_entries(value, errors.get(name, 0.0)) < min_proc_neff: - return False return True @@ -368,7 +371,15 @@ def minor_backgrounds(bkg_hists_by_name, lo, hi, min_frac): return {name for name, value in yields.items() if value < min_frac * total} -def _bin_passes(yields, min_each, exempt=(), total_error=None, min_neff=0.0): +def _bin_passes( + yields, + min_each, + exempt=(), + total_error=None, + min_neff=0.0, + errors=None, + min_proc_neff=0.0, +): """Mass-bin validity: every *relevant* background must exceed min_each, and (when min_neff > 0) the summed background must be known to at least min_neff effective MC entries. @@ -392,10 +403,29 @@ def _bin_passes(yields, min_each, exempt=(), total_error=None, min_neff=0.0): if effective_entries(total, total_error) < min_neff: return False for name, value in yields.items(): + # Positivity is never exemptable. `exempt` exists so a negligible process + # cannot veto every candidate split -- it is about the *magnitude* floor, not + # about sign. Letting an exempt process go negative does not merely leave an + # ugly bin: resolveNegativeBins in dc_make/maker.py rejects a negative bin + # that holds >= signalFractionForRelevantBins of the signal *unconditionally*, + # ignoring allow_negative_bins_within_error, so the datacard fails to build at + # all. Observed at m1000 muMu/SR/recovery_dnn1, where ST was 4.24% of the + # slice -- just under min_bkg_frac 0.05, hence exempt -- and came out at + # -0.567 +- 0.789, killing the production. + if value < 0: + return False if name in exempt: continue if value <= min_each: return False + # Every background must be measured, not merely present: a yield known only + # to a few hundred percent is not a background estimate. _slice_passes has + # carried this arm for the sliced axis; the mass axis, where essentially all + # fit bins live, had no per-process test at all -- only the summed one, which + # any single well-measured process satisfies on its own. + if min_proc_neff > 0 and errors is not None: + if effective_entries(value, errors.get(name, 0.0)) < min_proc_neff: + return False return True @@ -551,7 +581,13 @@ def signal_quantile_ranges(sig_hist, n_bins, first_bin, last_bin): def merge_until_valid( - ranges, sig_hist, bkg_hists_by_name, min_each, exempt=(), min_neff=0.0 + ranges, + sig_hist, + bkg_hists_by_name, + min_each, + exempt=(), + min_neff=0.0, + min_proc_neff=0.0, ): """Merge adjacent ranges until every one satisfies the background gates. @@ -572,6 +608,8 @@ def merge_until_valid( exempt, _total_bkg_error(bkg_hists_by_name, lo, hi), min_neff, + _bkg_errors(bkg_hists_by_name, lo, hi) if min_proc_neff > 0 else None, + min_proc_neff, ): bad = i break @@ -599,6 +637,7 @@ def find_bins( min_each, min_frac=0.0, min_neff=0.0, + min_proc_neff=0.0, ): """Bins inside one slice: signal quantiles for the edges, background gates for the count. @@ -620,7 +659,7 @@ def find_bins( exempt = minor_backgrounds(bkg_hists_by_name, first_bin, last_bin, min_frac) ranges = signal_quantile_ranges(sig_hist, max_bins, first_bin, last_bin) return merge_until_valid( - ranges, sig_hist, bkg_hists_by_name, min_each, exempt, min_neff + ranges, sig_hist, bkg_hists_by_name, min_each, exempt, min_neff, min_proc_neff ) @@ -671,6 +710,7 @@ def discover_binning( sig_mode="sb", min_slice_bkg_each=0.0, min_slice_bkg_each_neff=0.0, + min_bin_bkg_each_neff=0.0, ): """bkg2d_by_name: {background_name: [hist per discovery era, ...]}. The list is usually a single era's own histogram (standalone limit) or all of a @@ -728,6 +768,7 @@ def discover_binning( min_bin_each, min_bkg_frac, min_bin_bkg_neff, + min_bin_bkg_each_neff, ) bin_ranges = extend_outer_edges(bin_ranges, 0, ny + 1) result.append({"x_range": (xlo, xhi), "y_ranges": bin_ranges}) @@ -1003,6 +1044,7 @@ def load2d(f, key): knobs["significance_mode"], knobs["min_slice_bkg_each"], knobs["min_slice_bkg_each_neff"], + knobs["min_bin_bkg_each_neff"], ) if slices is None: print( @@ -1285,6 +1327,13 @@ def run( float, "minimum effective MC entries of the summed background in a bin", ), + "min_bin_bkg_each_neff": ( + float, + "minimum effective MC entries of every non-negligible background in a bin. " + "The summed test above is satisfied by any one well-measured process, so " + "this is what makes each background individually measured; bins are merged " + "into their lower-signal neighbour until it holds", + ), "bkg_per_bin": ( float, "summed background to aim for per bin; overrides max_bins_per_slice when it " From 983bdbc58a8b28ccd953354cddcb71f39ef17e62 Mon Sep 17 00:00:00 2001 From: Devin Date: Sat, 26 Sep 2026 00:26:03 +0200 Subject: [PATCH 03/12] Move the dimension-agnostic binning code into common/binning_core.py Everything that decides *where* a boundary goes -- slice finding and the gates it applies, which processes are minor enough to exempt, the figures of merit, and the histogram plumbing -- moves out of bin_opt_2d/rebin_2d.py into common/binning_core.py, unchanged. The cell tables the HME-box search is built on need the same answers to the same questions, and importing them from one place is what keeps a second copy from drifting. What stays behind in rebin_2d.py is the part that knows what a "bin" is, plus BINNING_DEFAULTS and load_binning_config(), whose "unknown key raises" contract only works if the file that reads the knobs states its own valid keys. No behaviour change: the moved code is byte-identical to the lines it came from, and the only other change to rebin_2d.py is its import header. Carved out of e205797 on Binning_2D_DP, without the 1D binner that commit also added. Co-Authored-By: Claude Opus 5.5 (1M context) --- bin_opt_2d/rebin_2d.py | 477 +++------------------------------------- common/binning_core.py | 480 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 504 insertions(+), 453 deletions(-) create mode 100644 common/binning_core.py diff --git a/bin_opt_2d/rebin_2d.py b/bin_opt_2d/rebin_2d.py index 41ad31b..a5dca54 100644 --- a/bin_opt_2d/rebin_2d.py +++ b/bin_opt_2d/rebin_2d.py @@ -1,6 +1,5 @@ import array import json -import math import os import sys import yaml @@ -15,26 +14,40 @@ __package__ = pkg_dir_name from StatInference.common.tools import importROOT, CategoryNaming -from StatInference.common.param_parse import extractParameters, applyParameters -from StatInference.dc_make.model import Model +from StatInference.common.param_parse import applyParameters + +# The dimension-agnostic half of the binning: slice finding, the gates it applies, +# the figures of merit, and the histogram plumbing. See that module's docstring for +# what stays here and why. +from StatInference.common.binning_core import ( + BINNING_JSON, + SIGNIFICANCE_MODES, + _bkg_errors, + _bkg_yields, + _detach, + _total_bkg_error, + bin_budget, + bin_edges, + effective_entries, + extend_outer_edges, + find_slices, + get_hist, + load_config, + lookup_frozen, + minor_backgrounds, + open_input_file, + sum_hists, +) ROOT = importROOT() -# Written inside each era's own output directory, beside the shapes it describes, and -# read back by --binning to replay them. -BINNING_JSON = "binning.json" - # The knobs that decide the binning, with the values a configuration gets when it does # not say otherwise. They belong to an analysis rather than to this file -- what a slice # needs to be worth keeping depends on the sample sizes and the selection, and a test # configuration wants them at zero -- so every production should state its own in a # binning yaml -- see config/Datacards/binning_2d.yaml in HH_bbWW for an annotated # set. This file ships no configuration of its own, only these fallbacks. -# The figures of merit significance() implements. Named here so the yaml path can be -# checked against the same list the command line's choices= uses. -SIGNIFICANCE_MODES = ("sb", "asimov") - BINNING_DEFAULTS = { "n_slices": 4, "category_pattern": None, # None -> CategoryNaming's own neutral default @@ -86,291 +99,6 @@ def load_binning_config(path, overrides=None): return knobs -def lookup_frozen(frozen_binning, era, mass, channel, category): - """The recorded binning for one channel/category, or None if it was skipped then. - - A category absent from the record was skipped by the run that wrote it (too little - signal, or a background missing from a discovery era), so it is skipped again rather - than quietly re-optimised -- a replay that rebinned some categories and froze others - would be neither the old binning nor a new one. - """ - return ( - frozen_binning.get("binning", {}) - .get(era, {}) - .get(str(mass), {}) - .get(channel, {}) - .get(category) - ) - - -def load_config(config_path): - with open(config_path, "r") as f: - cfg = yaml.safe_load(f) - model = Model.fromConfig(cfg["model"]) - channels = cfg["channels"] - # Taken verbatim here; run() reduces them to the base categories the 2D input is - # actually keyed by, once it has the pattern to do it with. - categories = list(cfg["categories"]) - - # Several processes may carry is_signal (e.g. the bbWW(2l) and bbtautau decay - # modes of the same resonance, both scaled by the same signal strength). They - # are summed to form the discovery signal that steers the slice - # boundaries, so the binning is optimised for the total signal in the fit. - signal_hist_names = [] - mass_values = None - background_entries = [] - for entry in cfg["processes"]: - if type(entry) == str: - background_entries.append((entry, entry, [])) - continue - if entry.get("is_data", False): - continue - base_name = entry["process"] - hist_name = entry.get("hist_name", base_name) - if entry.get("is_signal", False): - signal_hist_names.append(hist_name) - if mass_values is None: - mass_values = entry["param_values"] - elif list(entry["param_values"]) != list(mass_values): - raise RuntimeError( - f"Signal {hist_name} has param_values {entry['param_values']}, " - f"which differ from {mass_values}; every signal must be defined " - "at the same mass points" - ) - else: - background_entries.append((base_name, hist_name, entry.get("channels", []))) - - if not signal_hist_names: - raise RuntimeError("No signal process found in config") - - return { - "model": model, - "channels": channels, - "categories": categories, - "era_groups": cfg.get("era_groups", {}), - # Declared by a configuration that reads sliced shapes, so that run() can check - # it against the pattern the shapes are actually written with. None when the - # configuration names no pattern, and then there is nothing to check. - "category_pattern": cfg.get("category_pattern"), - "signal_hist_name_patterns": signal_hist_names, - "signal_param_name": extractParameters(signal_hist_names[0])[0], - "mass_values": mass_values, - "background_entries": background_entries, - } - - -def open_input_file(input_dir, model, era, mass, param_name): - file_name = model.getInputFileName(era, {param_name: mass}) - full_path = os.path.join(input_dir, file_name) - f = ROOT.TFile.Open(full_path, "READ") - if f is None or f.IsZombie(): - raise RuntimeError(f"Cannot open file {full_path}") - return f - - -def _detach(h): - """Detach a histogram from its TFile and hand ownership to Python. - - SetDirectory(0) alone makes the histogram survive the file's close, but - leaves it owned by nobody -- so every 2D histogram read here (one per - systematic variation, per category, per mass) leaked for the lifetime of the - process. Across the ten-mass loop that grew without bound and got the job - SIGKILLed partway through the final mass, leaving a truncated ROOT file. - SetOwnership makes the object die with its last Python reference. - """ - h.SetDirectory(0) - ROOT.SetOwnership(h, True) - return h - - -def get_hist(f, path): - h = f.Get(path) - # TFile.Get() on a fully-missing nested path can return a PyROOT wrapper - # around a null C++ pointer, which is not `is None` but is falsy. - if not h: - return None - return _detach(h) - - -def sum_hists(hists): - total = None - for h in hists: - if h is None: - continue - if total is None: - total = _detach(h.Clone()) - else: - total.Add(h) - return total - - -def _integral(hist, lo, hi): - return ( - hist.Integral(lo, hi, 0, -1) - if hist.GetDimension() == 2 - else hist.Integral(lo, hi) - ) - - -def _integral_and_error(hist, lo, hi): - """(yield, MC statistical error) over [lo, hi].""" - err = array.array("d", [0.0]) - if hist.GetDimension() == 2: - value = hist.IntegralAndError(lo, hi, 0, -1, err) - else: - value = hist.IntegralAndError(lo, hi, err) - return value, err[0] - - -def _bkg_yields(bkg_hists_by_name, lo, hi): - """For each background, the summed yield across all discovery eras (combined - statistics). Individual eras are not required to individually clear the - threshold -- use allow_negative_bins_within_error (maker.py, per process/ - category) for backgrounds/categories where a specific era can land negative - in a bin that's fine once combined.""" - return { - name: sum(_integral(h, lo, hi) for h in hists) - for name, hists in bkg_hists_by_name.items() - } - - -def _bkg_errors(bkg_hists_by_name, lo, hi): - """Per background, the MC statistical error on the summed yield (eras added in - quadrature).""" - return { - name: math.sqrt(sum(_integral_and_error(h, lo, hi)[1] ** 2 for h in hists)) - for name, hists in bkg_hists_by_name.items() - } - - -def _total_bkg_error(bkg_hists_by_name, lo, hi): - errors = _bkg_errors(bkg_hists_by_name, lo, hi) - return math.sqrt(sum(e**2 for e in errors.values())) - - -def effective_entries(value, error): - """(yield / error)^2 -- the unweighted MC event count a weighted yield is worth. - - A DY slice holding a couple of very-high-weight aMC@NLO events can carry a - sizeable yield with an effective count far below 1, i.e. a background estimate - that is statistically compatible with almost anything. - """ - if error <= 0: - return float("inf") if value > 0 else 0.0 - return (value / error) ** 2 - - -def asimov_significance(s, b, b_err=0.0): - """Median discovery significance for a counting experiment (Cowan et al., - arXiv:1007.1727), with the background uncertainty folded in when given. - - Reduces to sqrt(2*((s+b)*ln(1+s/b) - s)) for b_err = 0. Unlike S/sqrt(B) this - stays valid when b is O(1), which is exactly the regime the high-mass boosted - slices live in -- there S/sqrt(B) reports significances of ~30 on ~1.5 - background events, and drives the slice boundary on that basis. - """ - if b <= 0 or s <= 0: - return 0.0 - if b_err <= 0: - return math.sqrt(max(2.0 * ((s + b) * math.log1p(s / b) - s), 0.0)) - var = b_err**2 - term1 = (s + b) * math.log(((s + b) * (b + var)) / (b * b + (s + b) * var)) - term2 = (b * b / var) * math.log1p(var * s / (b * (b + var))) - return math.sqrt(max(2.0 * (term1 - term2), 0.0)) - - -def significance(s, b, b_err=0.0, mode="sb"): - """Figure of merit steering the slice boundaries. - - mode="sb": S / sqrt(B + sigma_B^2). Folding sigma_B into plain S/sqrt(B) - stops a downward fluctuation of a statistics-starved background - (DY in the b-tagged muMu slices, effective MC count <1) from - inflating the apparent significance. Still assumes the Gaussian - regime, which breaks down for B of order a few. - mode="asimov": the Poisson-correct Asimov significance, valid at low B. - """ - if b is None or b <= 0: - return 0.0 - if mode == "asimov": - return asimov_significance(s, b, b_err) - denominator = b + b_err**2 - if denominator <= 0: - return 0.0 - return s / (denominator**0.5) - - -def _slice_passes( - yields, - min_sum, - total_error=None, - min_neff=0.0, - errors=None, - min_each=0.0, - min_proc_neff=0.0, - exempt=(), -): - """Slice validity: the summed background must clear min_sum, and be known - to better than min_neff effective MC entries -- a window whose background is - statistically undetermined must not be selectable at all. - - The per-process arms (min_each, min_proc_neff) are what _bin_passes already - does on the mass axis, applied here as well. Testing only the sum hides an - individual background that has fluctuated negative behind a large, well - measured neighbour: muMu/SR/res2b at m500 selected a slice holding - DY = -5.87 +- 8.92 (N_eff 0.43) because the summed background there was - +39.9 +- 9.0 (N_eff 22), clearing both summed tests comfortably. That slice - cannot be turned into a datacard at all -- a negative DY integral is rejected - outright by resolveNegativeBins -- so the binning has to not select it in the - first place. Worse, the significance being maximised is S/sqrt(B + sigma_B^2), - which a downward fluctuation *increases* by lowering B, so such a window is - mildly preferred rather than merely tolerated. - - `exempt` comes from minor_backgrounds() judged over the whole sliced-axis range of the - category, never over the candidate window: a background that is negligible in - this category must not be able to veto every boundary, and judging it inside - the window under test is circular (a background that is exactly zero there is - trivially below any fraction of the total). - """ - total = sum(yields.values()) - if total <= min_sum: - return False - if min_neff > 0 and total_error is not None: - if effective_entries(total, total_error) < min_neff: - return False - for name, value in yields.items(): - # As in _bin_passes: exemption covers the magnitude floor, never the sign. - if value < 0: - return False - if name in exempt: - continue - if min_each > 0 and value <= min_each: - return False - if min_proc_neff > 0 and errors is not None: - if effective_entries(value, errors.get(name, 0.0)) < min_proc_neff: - return False - return True - - -def minor_backgrounds(bkg_hists_by_name, lo, hi, min_frac): - """Backgrounds negligible over the *whole* [lo, hi] range, which may therefore - be exempted from the per-bin min_each floor. - - This must be judged once over the full slice, never inside the candidate bin - under test. Evaluating the fraction per bin is circular: a background that is - exactly zero in that bin is trivially below any fraction of the bin total, so - it is always exempted -- which is precisely the case the floor exists to - catch. That let a bin through with TT = 0 in a slice where TT is 94% of the - background (m300, eMu, res2b_dnn2). - """ - if min_frac <= 0: - return set() - yields = _bkg_yields(bkg_hists_by_name, lo, hi) - total = sum(yields.values()) - if total <= 0: - return set() - return {name for name, value in yields.items() if value < min_frac * total} - - def _bin_passes( yields, min_each, @@ -429,114 +157,6 @@ def _bin_passes( return True -def grow_slice( - sig_hist, - bkg_hists_by_name, - right, - first_bin, - min_sum, - min_neff=0.0, - sig_mode="sb", - min_each=0.0, - min_proc_neff=0.0, - exempt=(), -): - """Among every candidate [left, right] with summed background > min_sum, - pick the one maximizing S/sqrt(B + sigma_B^2) -- not just the first one that - clears it. This is what actually drives where the cut lands; the content - floor is only a validity gate. Falls back to first_bin (best effort) if no - candidate clears the floor anywhere.""" - best_left = None - best_sig = -1.0 - for left in range(right, first_bin - 1, -1): - bkg_y = _bkg_yields(bkg_hists_by_name, left, right) - b_err = _total_bkg_error(bkg_hists_by_name, left, right) - bkg_e = ( - _bkg_errors(bkg_hists_by_name, left, right) if min_proc_neff > 0 else None - ) - if not _slice_passes( - bkg_y, - min_sum, - b_err, - min_neff, - bkg_e, - min_each, - min_proc_neff, - exempt, - ): - continue - s = _integral(sig_hist, left, right) - b = sum(bkg_y.values()) - sig = significance(s, b, b_err, sig_mode) - if sig > best_sig: - best_sig = sig - best_left = left - return best_left if best_left is not None else first_bin - - -def find_slices( - sig_hist, - bkg_hists_by_name, - n_slices, - first_bin, - last_bin, - min_sum, - min_neff=0.0, - sig_mode="sb", - min_each=0.0, - min_proc_neff=0.0, - min_frac=0.0, -): - """Split [first_bin, last_bin] into exactly `n_slices` ranges (fixed count, - required so every mass point shares the same category list), scanning right - to left. Each slice (except the final, leftover one) is placed to maximize - signal significance among boundaries with summed background > min_sum. - - Which backgrounds are minor enough to be exempt from the per-process floors is - decided once here, over the whole [first_bin, last_bin] range, and held fixed - for every candidate window -- see _slice_passes on why it cannot be re-judged - per window. - """ - exempt = ( - minor_backgrounds(bkg_hists_by_name, first_bin, last_bin, min_frac) - if (min_each > 0 or min_proc_neff > 0) - else set() - ) - slices = [] - right = last_bin - for slice_idx in range(n_slices): - if slice_idx == n_slices - 1 or right <= first_bin: - if right < first_bin: - # The axis ran out before n_slices boundaries were placed. Reported as - # None rather than invented: the previous code emitted (0,0), (-1,-1), - # (-2,-2) here to keep the count, and ROOT reads any range whose upper - # bound is below its lower one as the *whole axis including overflow* -- - # so those became copies of the full plane, and the same events entered - # the fit two or three times as independent categories. discover_binning - # rejects the category instead. - slices.append(None) - continue - slices.append((first_bin, right)) - right = first_bin - 1 - continue - left = grow_slice( - sig_hist, - bkg_hists_by_name, - right, - first_bin, - min_sum, - min_neff, - sig_mode, - min_each, - min_proc_neff, - exempt, - ) - slices.append((left, right)) - right = left - 1 - slices.reverse() - return slices - - def signal_quantile_ranges(sig_hist, n_bins, first_bin, last_bin): """Split [first_bin, last_bin] into `n_bins` ranges each holding an equal share of the signal. @@ -663,39 +283,6 @@ def find_bins( ) -def extend_outer_edges(ranges, full_lo, full_hi): - """Widen the first/last range to swallow under/overflow (bin 0 / nbins+1), - so no events are silently dropped at the extremes of the axis. - - Only real ranges are widened. Widening a None (an exhausted slice) would invent one, - and widening past a None would put the underflow on a range that is not the outermost - -- so a list containing any None is left alone and rejected by the caller. - """ - ranges = list(ranges) - if any(r is None for r in ranges): - return ranges - ranges[0] = (full_lo, ranges[0][1]) - ranges[-1] = (ranges[-1][0], full_hi) - return ranges - - -def bin_budget(bkg_hists_by_name, lo, hi, max_bins_per_slice, bkg_per_bin): - """How many bins this slice can actually afford. - - A fixed max_bins_per_slice is applied blind to slice content: the high-mass boosted - slices hold ~1.5 total background events and were still being split into 10 - bins, i.e. ~0.15 events per bin. Capping at B_slice / bkg_per_bin ties the - binning to the statistics that are really there. Returns max_bins_per_slice - unchanged when bkg_per_bin <= 0 (feature off). - """ - if bkg_per_bin <= 0: - return max_bins_per_slice - total = sum(_bkg_yields(bkg_hists_by_name, lo, hi).values()) - if total <= 0: - return 1 - return max(1, min(max_bins_per_slice, int(total / bkg_per_bin))) - - def discover_binning( sig2d, bkg2d_by_name, @@ -775,22 +362,6 @@ def discover_binning( return result -def bin_edges(y_axis, y_ranges): - """Physical edges of the discovered y ranges, for booking the output TH1. - - The rebinned shapes used to be booked as n_bins over [0, n_bins], which threw - the axis scale away and left every plot labelled by bin index. The ranges are - contiguous and ordered, so the edges are each range's low edge plus the last - range's upper edge. extend_outer_edges() pushes the outer ranges into the - underflow/overflow bins, which have no finite edge of their own -- those are - clamped to the axis limits. - """ - n = y_axis.GetNbins() - edges = [y_axis.GetBinLowEdge(max(lo, 1)) for lo, _ in y_ranges] - edges.append(y_axis.GetBinUpEdge(min(y_ranges[-1][1], n))) - return edges - - def mkdir_titled(directory, path, title): """mkdir a nested path, putting `title` on the leaf directory only. diff --git a/common/binning_core.py b/common/binning_core.py new file mode 100644 index 0000000..5750181 --- /dev/null +++ b/common/binning_core.py @@ -0,0 +1,480 @@ +"""The parts of the shape rebinning that do not care how many axes a histogram has. + +Split out of bin_opt_2d/rebin_2d.py, which is where all of this was written and where +the comments' point of view comes from. It moved here unchanged so that anything else +that needs the *same* answers to the same questions -- what a slice has to contain to be +worth keeping, which processes are minor enough to exempt, which figure of merit decides +a boundary -- reads them from here, rather than from a second copy that would be a second +place for them to drift. + +`_integral()` branches on GetDimension(), so everything downstream of it reads a TH1 and +a TH2 alike. What stays behind in the binner is the part that knows what a "bin" is: the +2D driver treats a slice as a directory of mass bins. + +Nothing here is configuration. BINNING_DEFAULTS and load_binning_config() stay with the +binner, because their "unknown key raises" contract only works if the file that reads the +knobs states its own valid keys. +""" + +import array +import math +import os + +import yaml + +from StatInference.common.tools import importROOT +from StatInference.common.param_parse import extractParameters +from StatInference.dc_make.model import Model + +ROOT = importROOT() + + +# Written inside each era's own output directory, beside the shapes it describes, and +# read back by --binning to replay them. +BINNING_JSON = "binning.json" + +# The figures of merit significance() implements. Named here so the yaml path can be +# checked against the same list the command line's choices= uses. +SIGNIFICANCE_MODES = ("sb", "asimov") + + +def lookup_frozen(frozen_binning, era, mass, channel, category): + """The recorded binning for one channel/category, or None if it was skipped then. + + A category absent from the record was skipped by the run that wrote it (too little + signal, or a background missing from a discovery era), so it is skipped again rather + than quietly re-optimised -- a replay that rebinned some categories and froze others + would be neither the old binning nor a new one. + """ + return ( + frozen_binning.get("binning", {}) + .get(era, {}) + .get(str(mass), {}) + .get(channel, {}) + .get(category) + ) + + +def load_config(config_path): + with open(config_path, "r") as f: + cfg = yaml.safe_load(f) + model = Model.fromConfig(cfg["model"]) + channels = cfg["channels"] + # Taken verbatim here; run() reduces them to the base categories the 2D input is + # actually keyed by, once it has the pattern to do it with. + categories = list(cfg["categories"]) + + # Several processes may carry is_signal (e.g. the bbWW(2l) and bbtautau decay + # modes of the same resonance, both scaled by the same signal strength). They + # are summed to form the discovery signal that steers the slice + # boundaries, so the binning is optimised for the total signal in the fit. + signal_hist_names = [] + mass_values = None + background_entries = [] + for entry in cfg["processes"]: + if type(entry) == str: + background_entries.append((entry, entry, [])) + continue + if entry.get("is_data", False): + continue + base_name = entry["process"] + hist_name = entry.get("hist_name", base_name) + if entry.get("is_signal", False): + signal_hist_names.append(hist_name) + if mass_values is None: + mass_values = entry["param_values"] + elif list(entry["param_values"]) != list(mass_values): + raise RuntimeError( + f"Signal {hist_name} has param_values {entry['param_values']}, " + f"which differ from {mass_values}; every signal must be defined " + "at the same mass points" + ) + else: + background_entries.append((base_name, hist_name, entry.get("channels", []))) + + if not signal_hist_names: + raise RuntimeError("No signal process found in config") + + return { + "model": model, + "channels": channels, + "categories": categories, + "era_groups": cfg.get("era_groups", {}), + # Declared by a configuration that reads sliced shapes, so that run() can check + # it against the pattern the shapes are actually written with. None when the + # configuration names no pattern, and then there is nothing to check. + "category_pattern": cfg.get("category_pattern"), + "signal_hist_name_patterns": signal_hist_names, + "signal_param_name": extractParameters(signal_hist_names[0])[0], + "mass_values": mass_values, + "background_entries": background_entries, + } + + +def open_input_file(input_dir, model, era, mass, param_name): + file_name = model.getInputFileName(era, {param_name: mass}) + full_path = os.path.join(input_dir, file_name) + f = ROOT.TFile.Open(full_path, "READ") + if f is None or f.IsZombie(): + raise RuntimeError(f"Cannot open file {full_path}") + return f + + +def _detach(h): + """Detach a histogram from its TFile and hand ownership to Python. + + SetDirectory(0) alone makes the histogram survive the file's close, but + leaves it owned by nobody -- so every 2D histogram read here (one per + systematic variation, per category, per mass) leaked for the lifetime of the + process. Across the ten-mass loop that grew without bound and got the job + SIGKILLed partway through the final mass, leaving a truncated ROOT file. + SetOwnership makes the object die with its last Python reference. + """ + h.SetDirectory(0) + ROOT.SetOwnership(h, True) + return h + + +def get_hist(f, path): + h = f.Get(path) + # TFile.Get() on a fully-missing nested path can return a PyROOT wrapper + # around a null C++ pointer, which is not `is None` but is falsy. + if not h: + return None + return _detach(h) + + +def sum_hists(hists): + total = None + for h in hists: + if h is None: + continue + if total is None: + total = _detach(h.Clone()) + else: + total.Add(h) + return total + + +def _integral(hist, lo, hi): + return ( + hist.Integral(lo, hi, 0, -1) + if hist.GetDimension() == 2 + else hist.Integral(lo, hi) + ) + + +def _integral_and_error(hist, lo, hi): + """(yield, MC statistical error) over [lo, hi].""" + err = array.array("d", [0.0]) + if hist.GetDimension() == 2: + value = hist.IntegralAndError(lo, hi, 0, -1, err) + else: + value = hist.IntegralAndError(lo, hi, err) + return value, err[0] + + +def _bkg_yields(bkg_hists_by_name, lo, hi): + """For each background, the summed yield across all discovery eras (combined + statistics). Individual eras are not required to individually clear the + threshold -- use allow_negative_bins_within_error (maker.py, per process/ + category) for backgrounds/categories where a specific era can land negative + in a bin that's fine once combined.""" + return { + name: sum(_integral(h, lo, hi) for h in hists) + for name, hists in bkg_hists_by_name.items() + } + + +def _bkg_errors(bkg_hists_by_name, lo, hi): + """Per background, the MC statistical error on the summed yield (eras added in + quadrature).""" + return { + name: math.sqrt(sum(_integral_and_error(h, lo, hi)[1] ** 2 for h in hists)) + for name, hists in bkg_hists_by_name.items() + } + + +def _total_bkg_error(bkg_hists_by_name, lo, hi): + errors = _bkg_errors(bkg_hists_by_name, lo, hi) + return math.sqrt(sum(e**2 for e in errors.values())) + + +def effective_entries(value, error): + """(yield / error)^2 -- the unweighted MC event count a weighted yield is worth. + + A DY slice holding a couple of very-high-weight aMC@NLO events can carry a + sizeable yield with an effective count far below 1, i.e. a background estimate + that is statistically compatible with almost anything. + """ + if error <= 0: + return float("inf") if value > 0 else 0.0 + return (value / error) ** 2 + + +def asimov_significance(s, b, b_err=0.0): + """Median discovery significance for a counting experiment (Cowan et al., + arXiv:1007.1727), with the background uncertainty folded in when given. + + Reduces to sqrt(2*((s+b)*ln(1+s/b) - s)) for b_err = 0. Unlike S/sqrt(B) this + stays valid when b is O(1), which is exactly the regime the high-mass boosted + slices live in -- there S/sqrt(B) reports significances of ~30 on ~1.5 + background events, and drives the slice boundary on that basis. + """ + if b <= 0 or s <= 0: + return 0.0 + if b_err <= 0: + return math.sqrt(max(2.0 * ((s + b) * math.log1p(s / b) - s), 0.0)) + var = b_err**2 + term1 = (s + b) * math.log(((s + b) * (b + var)) / (b * b + (s + b) * var)) + term2 = (b * b / var) * math.log1p(var * s / (b * (b + var))) + return math.sqrt(max(2.0 * (term1 - term2), 0.0)) + + +def significance(s, b, b_err=0.0, mode="sb"): + """Figure of merit steering the slice boundaries. + + mode="sb": S / sqrt(B + sigma_B^2). Folding sigma_B into plain S/sqrt(B) + stops a downward fluctuation of a statistics-starved background + (DY in the b-tagged muMu slices, effective MC count <1) from + inflating the apparent significance. Still assumes the Gaussian + regime, which breaks down for B of order a few. + mode="asimov": the Poisson-correct Asimov significance, valid at low B. + """ + if b is None or b <= 0: + return 0.0 + if mode == "asimov": + return asimov_significance(s, b, b_err) + denominator = b + b_err**2 + if denominator <= 0: + return 0.0 + return s / (denominator**0.5) + + +def _slice_passes( + yields, + min_sum, + total_error=None, + min_neff=0.0, + errors=None, + min_each=0.0, + min_proc_neff=0.0, + exempt=(), +): + """Slice validity: the summed background must clear min_sum, and be known + to better than min_neff effective MC entries -- a window whose background is + statistically undetermined must not be selectable at all. + + The per-process arms (min_each, min_proc_neff) are what _bin_passes already + does on the mass axis, applied here as well. Testing only the sum hides an + individual background that has fluctuated negative behind a large, well + measured neighbour: muMu/SR/res2b at m500 selected a slice holding + DY = -5.87 +- 8.92 (N_eff 0.43) because the summed background there was + +39.9 +- 9.0 (N_eff 22), clearing both summed tests comfortably. That slice + cannot be turned into a datacard at all -- a negative DY integral is rejected + outright by resolveNegativeBins -- so the binning has to not select it in the + first place. Worse, the significance being maximised is S/sqrt(B + sigma_B^2), + which a downward fluctuation *increases* by lowering B, so such a window is + mildly preferred rather than merely tolerated. + + `exempt` comes from minor_backgrounds() judged over the whole sliced-axis range of the + category, never over the candidate window: a background that is negligible in + this category must not be able to veto every boundary, and judging it inside + the window under test is circular (a background that is exactly zero there is + trivially below any fraction of the total). + """ + total = sum(yields.values()) + if total <= min_sum: + return False + if min_neff > 0 and total_error is not None: + if effective_entries(total, total_error) < min_neff: + return False + for name, value in yields.items(): + # As in _bin_passes: exemption covers the magnitude floor, never the sign. + if value < 0: + return False + if name in exempt: + continue + if min_each > 0 and value <= min_each: + return False + if min_proc_neff > 0 and errors is not None: + if effective_entries(value, errors.get(name, 0.0)) < min_proc_neff: + return False + return True + + +def minor_backgrounds(bkg_hists_by_name, lo, hi, min_frac): + """Backgrounds negligible over the *whole* [lo, hi] range, which may therefore + be exempted from the per-bin min_each floor. + + This must be judged once over the full slice, never inside the candidate bin + under test. Evaluating the fraction per bin is circular: a background that is + exactly zero in that bin is trivially below any fraction of the bin total, so + it is always exempted -- which is precisely the case the floor exists to + catch. That let a bin through with TT = 0 in a slice where TT is 94% of the + background (m300, eMu, res2b_dnn2). + """ + if min_frac <= 0: + return set() + yields = _bkg_yields(bkg_hists_by_name, lo, hi) + total = sum(yields.values()) + if total <= 0: + return set() + return {name for name, value in yields.items() if value < min_frac * total} + + +def grow_slice( + sig_hist, + bkg_hists_by_name, + right, + first_bin, + min_sum, + min_neff=0.0, + sig_mode="sb", + min_each=0.0, + min_proc_neff=0.0, + exempt=(), +): + """Among every candidate [left, right] with summed background > min_sum, + pick the one maximizing S/sqrt(B + sigma_B^2) -- not just the first one that + clears it. This is what actually drives where the cut lands; the content + floor is only a validity gate. Falls back to first_bin (best effort) if no + candidate clears the floor anywhere.""" + best_left = None + best_sig = -1.0 + for left in range(right, first_bin - 1, -1): + bkg_y = _bkg_yields(bkg_hists_by_name, left, right) + b_err = _total_bkg_error(bkg_hists_by_name, left, right) + bkg_e = ( + _bkg_errors(bkg_hists_by_name, left, right) if min_proc_neff > 0 else None + ) + if not _slice_passes( + bkg_y, + min_sum, + b_err, + min_neff, + bkg_e, + min_each, + min_proc_neff, + exempt, + ): + continue + s = _integral(sig_hist, left, right) + b = sum(bkg_y.values()) + sig = significance(s, b, b_err, sig_mode) + if sig > best_sig: + best_sig = sig + best_left = left + return best_left if best_left is not None else first_bin + + +def find_slices( + sig_hist, + bkg_hists_by_name, + n_slices, + first_bin, + last_bin, + min_sum, + min_neff=0.0, + sig_mode="sb", + min_each=0.0, + min_proc_neff=0.0, + min_frac=0.0, +): + """Split [first_bin, last_bin] into exactly `n_slices` ranges (fixed count, + required so every mass point shares the same category list), scanning right + to left. Each slice (except the final, leftover one) is placed to maximize + signal significance among boundaries with summed background > min_sum. + + Which backgrounds are minor enough to be exempt from the per-process floors is + decided once here, over the whole [first_bin, last_bin] range, and held fixed + for every candidate window -- see _slice_passes on why it cannot be re-judged + per window. + """ + exempt = ( + minor_backgrounds(bkg_hists_by_name, first_bin, last_bin, min_frac) + if (min_each > 0 or min_proc_neff > 0) + else set() + ) + slices = [] + right = last_bin + for slice_idx in range(n_slices): + if slice_idx == n_slices - 1 or right <= first_bin: + if right < first_bin: + # The axis ran out before n_slices boundaries were placed. Reported as + # None rather than invented: the previous code emitted (0,0), (-1,-1), + # (-2,-2) here to keep the count, and ROOT reads any range whose upper + # bound is below its lower one as the *whole axis including overflow* -- + # so those became copies of the full plane, and the same events entered + # the fit two or three times as independent categories. discover_binning + # rejects the category instead. + slices.append(None) + continue + slices.append((first_bin, right)) + right = first_bin - 1 + continue + left = grow_slice( + sig_hist, + bkg_hists_by_name, + right, + first_bin, + min_sum, + min_neff, + sig_mode, + min_each, + min_proc_neff, + exempt, + ) + slices.append((left, right)) + right = left - 1 + slices.reverse() + return slices + + +def extend_outer_edges(ranges, full_lo, full_hi): + """Widen the first/last range to swallow under/overflow (bin 0 / nbins+1), + so no events are silently dropped at the extremes of the axis. + + Only real ranges are widened. Widening a None (an exhausted slice) would invent one, + and widening past a None would put the underflow on a range that is not the outermost + -- so a list containing any None is left alone and rejected by the caller. + """ + ranges = list(ranges) + if any(r is None for r in ranges): + return ranges + ranges[0] = (full_lo, ranges[0][1]) + ranges[-1] = (ranges[-1][0], full_hi) + return ranges + + +def bin_budget(bkg_hists_by_name, lo, hi, max_bins_per_slice, bkg_per_bin): + """How many bins this slice can actually afford. + + A fixed max_bins_per_slice is applied blind to slice content: the high-mass boosted + slices hold ~1.5 total background events and were still being split into 10 + bins, i.e. ~0.15 events per bin. Capping at B_slice / bkg_per_bin ties the + binning to the statistics that are really there. Returns max_bins_per_slice + unchanged when bkg_per_bin <= 0 (feature off). + """ + if bkg_per_bin <= 0: + return max_bins_per_slice + total = sum(_bkg_yields(bkg_hists_by_name, lo, hi).values()) + if total <= 0: + return 1 + return max(1, min(max_bins_per_slice, int(total / bkg_per_bin))) + + +def bin_edges(y_axis, y_ranges): + """Physical edges of the discovered y ranges, for booking the output TH1. + + The rebinned shapes used to be booked as n_bins over [0, n_bins], which threw + the axis scale away and left every plot labelled by bin index. The ranges are + contiguous and ordered, so the edges are each range's low edge plus the last + range's upper edge. extend_outer_edges() pushes the outer ranges into the + underflow/overflow bins, which have no finite edge of their own -- those are + clamped to the axis limits. + """ + n = y_axis.GetNbins() + edges = [y_axis.GetBinLowEdge(max(lo, 1)) for lo, _ in y_ranges] + edges.append(y_axis.GetBinUpEdge(min(y_ranges[-1][1], n))) + return edges From 30301577c7c83d8f9db3143a2dc61af3a5aa467b Mon Sep 17 00:00:00 2001 From: Devin Date: Sat, 26 Sep 2026 00:49:46 +0200 Subject: [PATCH 04/12] Check input_categories against the configured categories, not the resolved ones Once a category list is expanded against a rebinning record (the base "SR/res2b" becomes "SR/res2b_hmebox0"), the map's keys -- the names the configuration writes -- no longer appear in self.categories, and every key read as unknown. Main does not expand categories yet, so nothing is affected here; the HME-box binning PR adds the expansion. Co-Authored-By: Claude Opus 5.5 (1M context) --- dc_make/maker.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/dc_make/maker.py b/dc_make/maker.py index d25c87e..002c00b 100644 --- a/dc_make/maker.py +++ b/dc_make/maker.py @@ -70,13 +70,16 @@ def __init__( self.naming = CategoryNaming.fromConfig(cfg) self.categories = list(cfg["categories"]) # Where each category's shapes sit in the input, when that depends on the model - # point; see InputCategories. + # point; see InputCategories. Checked against the categories as the configuration + # writes them, not self.categories: once a category list is expanded against a + # rebinning record (base "SR/res2b" -> "SR/res2b_hmebox0"), the map's keys -- the + # configured names -- no longer appear in it, and every key would read as unknown. self.input_categories = InputCategories.fromConfig(cfg) - unknown = self.input_categories.unknown(self.categories) + unknown = self.input_categories.unknown(cfg["categories"]) if unknown: raise RuntimeError( f"input_categories maps {unknown}, which are not configured categories " - f"({self.categories})" + f"({list(cfg['categories'])})" ) self.signalFractionForRelevantBins = cfg["signalFractionForRelevantBins"] From dc8224fc4165333003ef3d9eeca5465eeb55a194 Mon Sep 17 00:00:00 2001 From: Devin Date: Sat, 26 Sep 2026 01:03:06 +0200 Subject: [PATCH 05/12] Add the hme_box binning strategy strategy: hme_box cuts one contiguous HME window per category and bins the DNN inside it, instead of cutting the DNN into slices and binning HME inside each. Everything outside the box is dropped, so the category becomes a single datacard bin (SR/res2b -> SR/res2b_hmebox0). Box and DNN binning are chosen together: a window is scored by the best DNN binning it admits (_box_value), not by its own integrated significance, which would only widen it to swallow signal and ignore the shape the fit uses. The DNN binning inside a box is an exact partition search (partition_dp) against the same Asimov Z^2 the greedy slice boundaries use, with the bin gates as a hard mask; _best_count() takes the best-scoring count up to max_bins_per_slice, and dp_min_bin_gain gives back any bin whose split raised Z^2 by less than that fraction of the category's total. The windows are scanned on a box_stride grid inside the box_signal_quantile range of the signal, then refined within +-box_refine. min_bkg_frac is judged once over the whole box, and the gates are those of _bin_passes, including min_bin_bkg_each_neff and the never-exempt positivity. common/binning_dp.py holds the prefix-sum cells that make this affordable (every window and candidate bin is an O(1) lookup instead of a ROOT integral), the vectorised box tables, and check_box_mask(), which re-derives a sample of the tables with the canonical _bin_passes on every run so the two spellings of the gate cannot drift. binning.json records the box layout as "y_range"/"x_ranges" beside the DNN-slice layout's "x_range"/"y_ranges", and replays either; both strategies now also record the objective (Z, Z^2 per slice, bin count) recomputed from the written ranges. CreateDatacardsTask labels the shape plots by the surviving axis read from that record, so box panels are labelled by the DNN rather than by HME. Default stays "greedy", which is unchanged by this commit. Carved out of Binning_2D_DP (92e0010 and the parts of 39f9511, 88f612d, 14287b5 it uses), without the "dp" strategy, the budget allocator, axis alternation, binning_exclude_processes, fill_negative_bins, fitted tails, per-process N_eff floors, exempt_scope, the dc_make base-category expansion and the plot book. Verified identical to Binning_2D_DP on v2605a_DNNOutputs_v3 Run3_Early, MX 500 and 1000: binning.json equal apart from the dropped knobs (all at their defaults there), and every one of 188928 hme_box and 755712 greedy rebinned histograms equal in edges, contents and errors. Co-Authored-By: Claude Opus 5.5 (1M context) --- bin_opt_2d/rebin_2d.py | 370 +++++++++++++++++---- common/binning_dp.py | 661 +++++++++++++++++++++++++++++++++++++ law/CreateDatacardsTask.py | 43 ++- 3 files changed, 1008 insertions(+), 66 deletions(-) create mode 100644 common/binning_dp.py diff --git a/bin_opt_2d/rebin_2d.py b/bin_opt_2d/rebin_2d.py index a5dca54..12df51e 100644 --- a/bin_opt_2d/rebin_2d.py +++ b/bin_opt_2d/rebin_2d.py @@ -1,7 +1,10 @@ import array import json +import math import os import sys + +import numpy as np import yaml if __name__ == "__main__": @@ -38,6 +41,14 @@ open_input_file, sum_hists, ) +from StatInference.common.binning_dp import ( + binning_objective, + build_cells, + box_tables, + check_box_mask, + partition_dp, + trim_by_marginal_gain, +) ROOT = importROOT() @@ -64,8 +75,34 @@ "min_bkg_frac": 0.05, "min_signal": 0.5, "significance_mode": "asimov", + # How the edges are found. "greedy" is the original: each slice boundary is the best + # one taken on its own, and the mass bins inside a slice are equal-signal quantiles + # merged until the gates pass. "hme_box" does the reverse: one HME window per + # category, with the DNN binned inside it by an exact partition search -- see + # _discover_hme_box(). Default stays "greedy" so adding this moves nothing until a + # configuration asks for it. + "strategy": "greedy", + # hme_box only: a bin must earn its place. The split that creates it has to raise + # this category's Z^2 by at least this fraction of the category's own achievable + # total, or it is given back. 0.0 leaves the count to max_bins_per_slice alone. + "dp_min_bin_gain": 0.0, + # hme_box only: the search grid for the box edges. A full scan is ny(ny+1)/2 windows + # and every one of them carries its own dynamic program over the binned axis, so the + # edges are scanned on a stride and then refined within +-box_refine of the winner. + "box_stride": 6, + "box_refine": 5, + # hme_box only: the box is searched only where the signal is. A window that starts + # beyond the resonance adds background and no signal; scanning it is wasted time, and + # including it in the grid costs resolution where the answer actually lies. + "box_signal_quantile": 0.001, } +# Named here so the yaml path can be checked against the same list the command line's +# choices= uses -- significance_mode already learned this lesson: it is read by a +# function that silently treats anything unrecognised as a default, so a typo in the +# yaml quietly changed the answer instead of failing. +BINNING_STRATEGIES = ("greedy", "hme_box") + def load_binning_config(path, overrides=None): """Merge the binning yaml over the defaults, then command-line overrides over that. @@ -96,6 +133,11 @@ def load_binning_config(path, overrides=None): f"{path or 'binning configuration'}: significance_mode " f"'{knobs['significance_mode']}' is not one of {sorted(SIGNIFICANCE_MODES)}." ) + if knobs["strategy"] not in BINNING_STRATEGIES: + raise RuntimeError( + f"{path or 'binning configuration'}: strategy '{knobs['strategy']}' is not " + f"one of {sorted(BINNING_STRATEGIES)}." + ) return knobs @@ -283,22 +325,7 @@ def find_bins( ) -def discover_binning( - sig2d, - bkg2d_by_name, - n_slices, - max_bins_per_slice, - min_slice_sum, - min_bin_each, - min_slice_bkg_neff=0.0, - min_bkg_frac=0.0, - min_bin_bkg_neff=0.0, - bkg_per_bin=0.0, - sig_mode="sb", - min_slice_bkg_each=0.0, - min_slice_bkg_each_neff=0.0, - min_bin_bkg_each_neff=0.0, -): +def _discover_greedy(sig2d, bkg2d_by_name, knobs): """bkg2d_by_name: {background_name: [hist per discovery era, ...]}. The list is usually a single era's own histogram (standalone limit) or all of a meta-era's sub-eras (combined limit) -- see --discovery-eras. @@ -306,6 +333,19 @@ def discover_binning( the same discovery reference's (already-summed) signal histogram: its x projection picks the significance-maximizing slice boundaries, and its y projection within each slice places the mass bin edges by signal quantile.""" + n_slices = knobs["n_slices"] + max_bins_per_slice = knobs["max_bins_per_slice"] + min_slice_sum = knobs["min_slice_bkg_sum"] + min_bin_each = knobs["min_bin_bkg_each"] + min_slice_bkg_neff = knobs["min_slice_bkg_neff"] + min_bkg_frac = knobs["min_bkg_frac"] + min_bin_bkg_neff = knobs["min_bin_bkg_neff"] + bkg_per_bin = knobs["bkg_per_bin"] + sig_mode = knobs["significance_mode"] + min_slice_bkg_each = knobs["min_slice_bkg_each"] + min_slice_bkg_each_neff = knobs["min_slice_bkg_each_neff"] + min_bin_bkg_each_neff = knobs["min_bin_bkg_each_neff"] + any_hist = next(iter(bkg2d_by_name.values()))[0] nx = any_hist.GetNbinsX() ny = any_hist.GetNbinsY() @@ -362,6 +402,146 @@ def discover_binning( return result +def _best_count(values, cap): + """The bin count at or below `cap` that actually scores best, or None if none is + feasible. + + Not the largest feasible count, which is the obvious choice and the wrong one. With + the background uncertainty folded into the figure of merit -- which is what + significance() does, and the reason it does not reward a background that fluctuated + low -- splitting a bin is no longer guaranteed to raise Z^2: each half carries a + larger relative sigma_B, and for a slice that is already statistics-limited the + split can cost more than the extra shape information gains. Taking the largest + feasible count would then hand back a binning worse than a coarser one that passes + the same gates. + + Choosing the argmax instead means every coarser feasible count is one of the + candidates, so the winner is at least as good as any of them. It also tends to spend + fewer bins than the budget allows, which is the right direction -- an extra bin that + buys nothing is pure MC-statistical exposure. + """ + best, best_n = -np.inf, None + for n in range(1, min(cap, len(values) - 1) + 1): + if values[n] > best: + best, best_n = values[n], n + return best_n + + +def _box_value(cells, y0, y1, knobs, check=False): + """What one HME box is worth: the best 1D binning of the DNN inside it. + + A box is scored by the binning it admits, not by its own integrated significance. + Those are different questions and they pick different windows -- scored on its own, a + box is rewarded for swallowing as much signal as it can and the DNN shape inside it + plays no part, which is the opposite of what the fit will do with it. + """ + mode = knobs["significance_mode"] + # which backgrounds count as negligible is decided once, over the whole box -- see + # minor_backgrounds() on why it cannot be re-judged inside each candidate bin + exempt = cells.exempt(1, cells.nx, y0, y1, knobs["min_bkg_frac"]) + score, valid = box_tables(cells, y0, y1, exempt, knobs, mode) + if check: + + def bin_passes(a, b, ylo, yhi, exempt_set): + return _bin_passes( + cells.yields(a, b, ylo, yhi), + knobs["min_bin_bkg_each"], + exempt_set, + cells.total_bkg_error(a, b, ylo, yhi), + knobs["min_bin_bkg_neff"], + cells.errors(a, b, ylo, yhi), + knobs["min_bin_bkg_each_neff"], + ) + + check_box_mask( + cells, + y0, + y1, + valid, + score, + exempt, + bin_passes, + mode, + ) + values, partitions = partition_dp(score, valid, knobs["max_bins_per_slice"]) + n = _best_count(values, knobs["max_bins_per_slice"]) + if n is None: + return float("-inf"), None + if knobs["dp_min_bin_gain"] > 0 and values[n] > 0: + n = trim_by_marginal_gain([values], [n], knobs["dp_min_bin_gain"] * values[n])[ + 0 + ] + return values[n], [(a + 1, b + 1) for a, b in partitions[n]] + + +def _discover_hme_box(sig2d, bkg2d_by_name, knobs): + """One HME window per category, with the DNN binned inside it. + + The greedy strategy cuts the DNN into slices that each become a datacard category and + gives every slice a shape along HME. This one does the reverse and keeps only one + category: a single contiguous HME box around the resonance, and a 1D DNN distribution + inside it. Everything outside the box is dropped -- that is what makes it a cut rather + than a slice, and it is why the box has to earn its acceptance loss in purity. + + Box and binning are chosen together; see _box_value. + """ + cells = build_cells(sig2d, bkg2d_by_name) + ny = cells.ny + # where the signal actually is -- see box_signal_quantile + cumulative = np.array( + [cells.signal(0, cells.nx + 1, 1, j) for j in range(1, ny + 1)] + ) + total = cumulative[-1] if len(cumulative) else 0.0 + if total <= 0: + return None + q = knobs["box_signal_quantile"] + lo_edge = max(1, int(np.searchsorted(cumulative, q * total))) + hi_edge = min(ny, int(np.searchsorted(cumulative, (1.0 - q) * total)) + 1) + + best_value, best_box = float("-inf"), None + stride = max(1, int(knobs["box_stride"])) + for a in range(lo_edge, hi_edge + 1, stride): + for b in range(a, hi_edge + 1, stride): + value, _ = _box_value(cells, a, b, knobs) + if value > best_value: + best_value, best_box = value, (a, b) + if best_box is None: + return None + refine = max(0, int(knobs["box_refine"])) + a0, b0 = best_box + for a in range(max(lo_edge, a0 - refine), min(hi_edge, a0 + refine) + 1): + for b in range(max(a, b0 - refine), min(hi_edge, b0 + refine) + 1): + value, _ = _box_value(cells, a, b, knobs) + if value > best_value: + best_value, best_box = value, (a, b) + + y0, y1 = best_box + value, bins = _box_value(cells, y0, y1, knobs, check=True) + if bins is None: + return None + # A box touching an end of the axis takes that end's overflow with it: the events are + # inside the selection the box stands for, and dropping them would make the written + # yield disagree with the cut the category claims to be. + y0 = 0 if y0 <= 1 else y0 + y1 = ny + 1 if y1 >= ny else y1 + return [{"y_range": (y0, y1), "x_ranges": bins}] + + +def discover_binning(sig2d, bkg2d_by_name, knobs): + """Find this category's binning, by whichever strategy the configuration names. + + bkg2d_by_name: {background_name: [hist per discovery era, ...]}. The list is usually + a single era's own histogram (standalone limit) or all of a meta-era's sub-eras + (combined limit). Yields are summed across whatever is in the list; see _bkg_yields(). + sig2d is the same discovery reference's already-summed signal histogram. + + Returns None when the category cannot be binned, which the caller reports as a skip. + """ + if knobs["strategy"] == "hme_box": + return _discover_hme_box(sig2d, bkg2d_by_name, knobs) + return _discover_greedy(sig2d, bkg2d_by_name, knobs) + + def mkdir_titled(directory, path, title): """mkdir a nested path, putting `title` on the leaf directory only. @@ -397,7 +577,7 @@ def format_var_range(lo, hi, var): return f"{lo:.2f} < {var} < {hi:.2f}" -def slice_ranges(x_axis, slices): +def slice_ranges(x_axis, y_axis, slices): """Physical edges of the discovered slices on the sliced axis, as [[lo, hi], ...]. The slice x_ranges are bin indices, and extend_outer_edges() has already pushed the @@ -406,39 +586,74 @@ def slice_ranges(x_axis, slices): the plots can say which selection each slice actually is; nothing downstream of the datacards needs it. """ - n = x_axis.GetNbins() ranges = [] for sl in slices: - lo, hi = sl["x_range"] + (lo, hi), _, slice_axis = slice_parts(sl) + axis = x_axis if slice_axis == "x" else y_axis + n = axis.GetNbins() ranges.append( [ - None if lo < 1 else x_axis.GetBinLowEdge(lo), - None if hi > n else x_axis.GetBinUpEdge(hi), + None if lo < 1 else axis.GetBinLowEdge(lo), + None if hi > n else axis.GetBinUpEdge(hi), ] ) return ranges -def slices_to_record(slices, x_axis, y_axis): +def _slice_record(sl, sel_edges, x_axis, y_axis): + """One category as plain data, in whichever layout produced it. + + The key names say which axis is the selection and which carries the bins, so a reader + -- human or replay -- never has to infer it from the strategy that happened to run. + """ + (lo, hi), bins, slice_axis = slice_parts(sl) + binned_axis = y_axis if slice_axis == "x" else x_axis + sel_key, bin_key = ( + ("x_range", "y_ranges") + if slice_axis == "x" + else ( + "y_range", + "x_ranges", + ) + ) + return { + sel_key: [lo, hi], + f"{slice_axis}_edges": sel_edges, + bin_key: [list(r) for r in bins], + ("y_edges" if slice_axis == "x" else "x_edges"): bin_edges(binned_axis, bins), + } + + +def slices_to_record(slices, x_axis, y_axis, objective=None): """The discovered structure as plain data, for binning.json. Both forms are written. The bin index ranges are what the code actually cuts on and are what a replay restores; the physical edges alongside them are what a human reads, and what makes the record meaningful next to a plot. + + `objective` is what the binning scored, recomputed from these ranges rather than + carried out of the optimiser -- so it is filled in for the greedy strategy too, which + never computes it, and two runs can be compared from their records alone. It is + written as a sibling of the slices and never read back: record_to_slices() takes only + the ranges, so a record written before this existed still replays. """ - return { + record = { "n_x_bins": x_axis.GetNbins(), "n_y_bins": y_axis.GetNbins(), "slices": [ - { - "x_range": list(sl["x_range"]), - "x_edges": x_edges, - "y_ranges": [list(r) for r in sl["y_ranges"]], - "y_edges": bin_edges(y_axis, sl["y_ranges"]), - } - for sl, x_edges in zip(slices, slice_ranges(x_axis, slices)) + _slice_record(sl, sel_edges, x_axis, y_axis) + for sl, sel_edges in zip(slices, slice_ranges(x_axis, y_axis, slices)) ], } + if objective is not None: + total, per_slice = objective + record["objective"] = { + "z": math.sqrt(max(total, 0.0)), + "z2": total, + "z2_per_slice": per_slice, + "n_bins": sum(len(slice_parts(sl)[1]) for sl in slices), + } + return record def record_to_slices(record, x_axis, y_axis, where): @@ -460,30 +675,52 @@ def record_to_slices(record, x_axis, y_axis, where): "to this input." ) return [ - { - "x_range": tuple(sl["x_range"]), - "y_ranges": [tuple(r) for r in sl["y_ranges"]], - } + ( + { + "y_range": tuple(sl["y_range"]), + "x_ranges": [tuple(r) for r in sl["x_ranges"]], + } + if "y_range" in sl + else { + "x_range": tuple(sl["x_range"]), + "y_ranges": [tuple(r) for r in sl["y_ranges"]], + } + ) for sl in record["slices"] ] +def slice_parts(sl): + """(selection range, bin ranges, which axis the selection is on) for one category. + + Two layouts reach this point. The DNN-sliced strategy cuts x into categories and gives + each a shape along y: {"x_range", "y_ranges"}. The HME-box strategy cuts y and gives + the one category a shape along x: {"y_range", "x_ranges"}. Everything downstream -- + the writer, the record, the objective -- needs the same three things from either, and + asking here is what keeps that code from having to know which strategy ran. + """ + if "y_range" in sl: + return sl["y_range"], sl["x_ranges"], "y" + return sl["x_range"], sl["y_ranges"], "x" + + def rebin_hist_2d(hist2d, slices, name, naming): """Given the discovered slice structure, produce one final TH1 per slice for this specific histogram (nominal or a systematic variation).""" outputs = [] for slice_idx, sl in enumerate(slices): - xlo, xhi = sl["x_range"] + (sel_lo, sel_hi), bin_ranges, slice_axis = slice_parts(sl) # ROOT reinterprets an inverted or negative range as the full axis including # under/overflow, silently and with a plausible-looking positive yield, so an # invalid range must never reach IntegralAndError below. Checked here rather # than trusted because this is the last point where it is still cheap to say so. - if not 0 <= xlo <= xhi: + if not 0 <= sel_lo <= sel_hi: raise RuntimeError( - f"invalid x bin range ({xlo}, {xhi}) for slice {slice_idx} of {name}; " - "ROOT would read this as the whole plane." + f"invalid {slice_axis} bin range ({sel_lo}, {sel_hi}) for slice " + f"{slice_idx} of {name}; ROOT would read this as the whole plane." ) - edges = array.array("d", bin_edges(hist2d.GetYaxis(), sl["y_ranges"])) + binned_axis = hist2d.GetYaxis() if slice_axis == "x" else hist2d.GetXaxis() + edges = array.array("d", bin_edges(binned_axis, bin_ranges)) # Detached for the same reason as the projections above: Write() targets # gDirectory regardless, so nothing needs these to stay attached. h = _detach( @@ -494,9 +731,14 @@ def rebin_hist_2d(hist2d, slices, name, naming): edges, ) ) - for bin_idx, (ylo, yhi) in enumerate(sl["y_ranges"], start=1): + for bin_idx, (blo, bhi) in enumerate(bin_ranges, start=1): err = array.array("d", [0.0]) - content = hist2d.IntegralAndError(xlo, xhi, ylo, yhi, err) + # the selection is on one axis and the bin on the other; IntegralAndError + # always wants (x_lo, x_hi, y_lo, y_hi) + if slice_axis == "x": + content = hist2d.IntegralAndError(sel_lo, sel_hi, blo, bhi, err) + else: + content = hist2d.IntegralAndError(blo, bhi, sel_lo, sel_hi, err) h.SetBinContent(bin_idx, content) h.SetBinError(bin_idx, err[0]) outputs.append(h) @@ -601,22 +843,7 @@ def load2d(f, key): frozen, disc_sig.GetXaxis(), disc_sig.GetYaxis(), where ) else: - slices = discover_binning( - disc_sig, - disc_bkg_by_name, - knobs["n_slices"], - knobs["max_bins_per_slice"], - knobs["min_slice_bkg_sum"], - knobs["min_bin_bkg_each"], - knobs["min_slice_bkg_neff"], - knobs["min_bkg_frac"], - knobs["min_bin_bkg_neff"], - knobs["bkg_per_bin"], - knobs["significance_mode"], - knobs["min_slice_bkg_each"], - knobs["min_slice_bkg_each_neff"], - knobs["min_bin_bkg_each_neff"], - ) + slices = discover_binning(disc_sig, disc_bkg_by_name, knobs) if slices is None: print( f" [skip] {channel}/{category} {param_name}={mass}: the " @@ -632,7 +859,7 @@ def load2d(f, key): # side-car path to hand a reader correctly and no key name for the two ends to agree # on -- anything that can open the shapes can already read it. naming = cfg["naming"] - ranges = slice_ranges(disc_sig.GetXaxis(), slices) + ranges = slice_ranges(disc_sig.GetXaxis(), disc_sig.GetYaxis(), slices) # One shared set of edges, applied to every source era in its own file. for source_era, in_file, out_file in sources: @@ -654,7 +881,21 @@ def load2d(f, key): out_file.cd(f"{channel}/{naming.name(category, slice_idx)}") h.Write(key) - return slices_to_record(slices, disc_sig.GetXaxis(), disc_sig.GetYaxis()) + # Scored from the ranges that were just written, for both strategies and for a replay + # alike, so the number in the record always describes the shapes beside it. Cheap + # next to the search, and the greedy path has no other way to report what it achieved. + # Only when there is something to score. A replay is allowed to proceed with a + # background missing from the input -- it applies recorded edges and asks no + # questions -- and a figure of merit computed against no background would be a + # number, but not a true one. + objective = ( + binning_objective( + build_cells(disc_sig, disc_bkg_by_name), slices, knobs["significance_mode"] + ) + if disc_bkg_by_name + else None + ) + return slices_to_record(slices, disc_sig.GetXaxis(), disc_sig.GetYaxis(), objective) def run( @@ -925,6 +1166,16 @@ def run( default=None, help=help_text, ) + parser.add_argument( + "--strategy", + required=False, + type=str, + default=None, + choices=list(BINNING_STRATEGIES), + help="how the edges are found: 'greedy' = each slice boundary taken on its own " + "and equal-signal quantiles inside it, 'hme_box' = one HME window per category " + "with the DNN binned inside it by an exact partition search", + ) parser.add_argument( "--significance-mode", required=False, @@ -938,6 +1189,7 @@ def run( overrides = {name: getattr(args, name) for name in knob_args} overrides["significance_mode"] = args.significance_mode + overrides["strategy"] = args.strategy knobs = load_binning_config(args.binning_config, overrides) frozen_binning = None diff --git a/common/binning_dp.py b/common/binning_dp.py new file mode 100644 index 0000000..ea0043a --- /dev/null +++ b/common/binning_dp.py @@ -0,0 +1,661 @@ +"""Prefix-sum cells over a 2D shape, so any candidate binning can be scored in O(1). + +The greedy binner asks its questions one candidate at a time and re-integrates the +histograms for each, which is affordable because it only ever looks at ~nx windows. The +HME-box search looks at ~ny^2 windows and runs an exact partition search over ~nx^2 +candidate DNN bins inside every one, so the integrals have to become table lookups or the +search is not worth running. +That is all this module does: it turns the TH2s into cumulative arrays and hands back the +same yields, errors and figure of merit the greedy path would have computed. + +It deliberately knows nothing about what a valid bin is. The gates live in +binning_core (and _bin_passes in the 2D binner), and are passed in as callables, so +there is exactly one definition of what a bin has to contain and this module cannot +drift from it. What it provides is the *arguments* those gates take -- the yields and +errors dicts -- built from the tables instead of from ROOT. + +Under- and overflow are carried, and that is load-bearing rather than tidy. _integral() +in binning_core reads a TH2 with Integral(lo, hi, 0, -1), and ROOT reads biny2 < biny1 as +"the whole y axis including under- and overflow", so every slice-level quantity in the +greedy path already includes the y under/overflow rows. A table that stopped at bin ny +would disagree with the greedy path by a fraction of an event, which is easily enough to +move a boundary and produce a difference nobody can explain. So the arrays are indexed +0..n+1 throughout, and a slice-level query spans y = 0..ny+1. + +Differencing prefix sums is not exact, and one of the ways it is inexact has teeth. +A window whose content is ~1e8 times smaller than the whole plane's -- routine for DY in +a narrow high-HME window -- loses about eight digits to cancellation, so a quantity that +is exactly 0.0 in ROOT comes back as -2e-16 here. That is numerically nothing and +physically nothing, but _bin_passes() tests `value < 0` to reject a background that has +gone negative, and it would reject 137 of the windows sampled below on the strength of a +rounding sign. So every yield is snapped to zero below a per-process epsilon scaled by +that process's own L1 norm, which is the size of the cancellation that can occur. The +epsilon lands around 1e-8 events for the largest process here, against a magnitude floor +(min_bin_bkg_each) of 0.01 -- six orders of margin, so the snap cannot mask a real +negative. + +The check that all of this is true is parity_check(), not inspection. It is the only +reason to trust anything built on top of this file. +""" + +import math + +import numpy as np + +from StatInference.common.binning_core import ( + _bkg_errors, + _bkg_yields, + _integral, + _total_bkg_error, + significance, +) + + +def _hist_arrays(hist): + """(values, variances) as (nx+2, ny+2) arrays indexed [binx, biny], bins 0..n+1. + + Read through GetBinContent/GetBinError rather than the internal buffer: the buffer + is faster but its layout is a ROOT implementation detail, and getting it wrong is + silent. This runs once per histogram per category, which is nothing next to the + search it feeds. + """ + nx = hist.GetNbinsX() + ny = hist.GetNbinsY() + values = np.empty((nx + 2, ny + 2), dtype=np.float64) + variances = np.empty((nx + 2, ny + 2), dtype=np.float64) + for bx in range(nx + 2): + for by in range(ny + 2): + values[bx, by] = hist.GetBinContent(bx, by) + err = hist.GetBinError(bx, by) + variances[bx, by] = err * err + return values, variances + + +def _prefix(a): + """Inclusive 2D prefix sums, shape (nx+3, ny+3), with P[i, j] = sum(a[:i, :j]).""" + p = np.zeros((a.shape[0] + 1, a.shape[1] + 1), dtype=np.float64) + p[1:, 1:] = a.cumsum(axis=0).cumsum(axis=1) + return p + + +def _cancellation_eps(a, rel=1e-12): + """How large a value this array's prefix sums can invent out of rounding. + + Scaled by the L1 norm rather than the sum, because the cancellation that matters + comes from negative-weight events: a process whose bins nearly cancel has a small + sum and a large L1, and it is the L1 that sets how much precision the difference of + two prefix sums can lose. rel is a few hundred times the double-precision epsilon, + which covers accumulating over the ~1e4 cells of one of these planes. + """ + return rel * float(np.abs(a).sum()) + + +def _snap(value, eps): + """Zero out a value that is indistinguishable from zero at this array's precision.""" + return 0.0 if abs(value) < eps else value + + +def _rect(p, x0, x1, y0, y1): + """Sum over the inclusive bin-index rectangle [x0..x1] x [y0..y1].""" + return float(p[x1 + 1, y1 + 1] - p[x0, y1 + 1] - p[x1 + 1, y0] + p[x0, y0]) + + +class Cells: + """Every quantity the binning gates and the figure of merit need, as O(1) lookups. + + Background yields are summed across discovery eras and their MC errors added in + quadrature, which is what _bkg_yields()/_bkg_errors() do -- the eras are summed + because the binning is derived from the combination it will be applied to. + """ + + def __init__(self, nx, ny, sig, bkg, var): + self.nx = nx + self.ny = ny + self.names = list(bkg) + self._sig = _prefix(sig) + self._bkg = {name: _prefix(bkg[name]) for name in self.names} + self._var = {name: _prefix(var[name]) for name in self.names} + # An empty background dict is reachable -- a replay whose input is missing a + # process it was derived with -- and sum() of nothing is the integer 0, which is + # not an array. Zeros keep every query defined and answering zero, which is the + # truthful answer when there is no background to integrate. + zero = np.zeros_like(sig) + bkg_tot = sum(bkg.values()) if bkg else zero + var_tot = sum(var.values()) if var else zero + self._bkg_tot = _prefix(bkg_tot) + self._var_tot = _prefix(var_tot) + self._eps_sig = _cancellation_eps(sig) + self._eps_bkg = {name: _cancellation_eps(bkg[name]) for name in self.names} + self._eps_tot = _cancellation_eps(bkg_tot) + self._eps_var = {name: _cancellation_eps(var[name]) for name in self.names} + self._eps_var_tot = _cancellation_eps(var_tot) + + # -- the whole y axis, under- and overflow included: what a slice-level query means + def full_y(self): + return 0, self.ny + 1 + + def signal(self, x0, x1, y0, y1): + return _snap(_rect(self._sig, x0, x1, y0, y1), self._eps_sig) + + def total_bkg(self, x0, x1, y0, y1): + return _snap(_rect(self._bkg_tot, x0, x1, y0, y1), self._eps_tot) + + def total_bkg_error(self, x0, x1, y0, y1): + return math.sqrt(max(_rect(self._var_tot, x0, x1, y0, y1), 0.0)) + + def variances(self, x0, x1, y0, y1): + """{background: squared MC statistical error}. Exposed because this, not the + error, is the quantity the tables actually add -- so it is the one a parity + check can bound the rounding of.""" + return {n: _rect(self._var[n], x0, x1, y0, y1) for n in self.names} + + def total_bkg_variance(self, x0, x1, y0, y1): + return _rect(self._var_tot, x0, x1, y0, y1) + + def yields(self, x0, x1, y0, y1): + """{background: yield}, the dict the gates take.""" + return { + n: _snap(_rect(self._bkg[n], x0, x1, y0, y1), self._eps_bkg[n]) + for n in self.names + } + + def errors(self, x0, x1, y0, y1): + """{background: MC statistical error}, the dict the gates take.""" + return { + n: math.sqrt(max(_rect(self._var[n], x0, x1, y0, y1), 0.0)) + for n in self.names + } + + def score(self, x0, x1, y0, y1, mode): + """The figure of merit for one cell, squared. + + Squared because Asimov Z^2 is what adds across independent counting bins, so a + partition's value is the sum of its cells' scores and the search can be a + dynamic program. significance() is binning_core's, so the cells are ranked by + exactly the quantity the greedy path maximises. + """ + s = self.signal(x0, x1, y0, y1) + b = self.total_bkg(x0, x1, y0, y1) + b_err = self.total_bkg_error(x0, x1, y0, y1) + return significance(s, b, b_err, mode) ** 2 + + def exempt(self, x0, x1, y0, y1, min_frac): + """minor_backgrounds() over a rectangle, for the range the caller is about to + subdivide -- never for the candidate sub-range itself, which is circular. + """ + if min_frac <= 0: + return set() + y = self.yields(x0, x1, y0, y1) + total = sum(y.values()) + if total <= 0: + return set() + return {n for n, v in y.items() if v < min_frac * total} + + +def build_cells(sig2d, bkg2d_by_name): + """Cells for one channel/category, from the same histograms discover_binning() reads. + + sig2d is the already-summed discovery signal; bkg2d_by_name is + {background: [hist per discovery era]}, summed here the way _bkg_yields() sums it. + """ + nx = sig2d.GetNbinsX() + ny = sig2d.GetNbinsY() + sig, _ = _hist_arrays(sig2d) + bkg = {} + var = {} + for name, hists in bkg2d_by_name.items(): + values = np.zeros((nx + 2, ny + 2), dtype=np.float64) + variances = np.zeros((nx + 2, ny + 2), dtype=np.float64) + for h in hists: + v, e2 = _hist_arrays(h) + values += v + variances += e2 + bkg[name] = values + var[name] = variances + return Cells(nx, ny, sig, bkg, var) + + +def parity_check(cells, sig2d, bkg2d_by_name, n_samples=200, seed=0): + """Assert the tables agree with the ROOT integrals they replace. + + Everything downstream is only as trustworthy as this. Over random x windows it + compares the quantities the gates and the figure of merit are built from -- + per-process yields and variances, the summed variance, and the signal -- against + binning_core's own _bkg_yields/_bkg_errors/_total_bkg_error/_integral, which is what + the greedy path calls. Slice level (full y, under- and overflow included), because + that is the query whose ROOT spelling is easy to get wrong. + + **The tolerance is absolute and scaled by the cancellation, not relative.** A + relative tolerance is the wrong test for a difference of prefix sums: a window + holding 1e-8 of the plane's content has already spent eight of its sixteen digits + before the comparison starts, so demanding a fixed relative agreement on it is + demanding that floating point not be floating point. What *can* be bounded is the + absolute rounding, and it is bounded by the same per-process epsilon Cells uses to + snap -- the L1 norm times a few hundred machine epsilons. Errors are compared as + variances for the same reason: the variance is what the tables add, and the square + root only inflates the relative error of an already-cancelled quantity. + + It also asserts the two agree on the *sign* of every yield, which is a separate + property from agreeing on its value and the one the binning actually turns on: + _bin_passes rejects a bin outright when a background is negative, so a yield that is + 0.0 in ROOT and -2e-16 here is not a rounding difference, it is a different binning. + + Raises on the first disagreement rather than returning a verdict: a table that is + subtly wrong must not be usable. + """ + rng = np.random.default_rng(seed) + y0, y1 = cells.full_y() + checked = 0 + for _ in range(n_samples): + lo = int(rng.integers(1, cells.nx + 1)) + hi = int(rng.integers(lo, cells.nx + 1)) + + got_y = cells.yields(lo, hi, y0, y1) + want_y = _bkg_yields(bkg2d_by_name, lo, hi) + got_v = cells.variances(lo, hi, y0, y1) + want_e = _bkg_errors(bkg2d_by_name, lo, hi) + for name in want_y: + _assert_close( + got_y[name], + want_y[name], + cells._eps_bkg[name], + f"yield[{name}] {lo}..{hi}", + ) + _assert_close( + got_v[name], + want_e[name] ** 2, + cells._eps_var[name], + f"variance[{name}] {lo}..{hi}", + ) + _assert_same_sign(got_y[name], want_y[name], f"yield[{name}] {lo}..{hi}") + + _assert_close( + cells.total_bkg_variance(lo, hi, y0, y1), + _total_bkg_error(bkg2d_by_name, lo, hi) ** 2, + cells._eps_var_tot, + f"total bkg variance {lo}..{hi}", + ) + _assert_close( + cells.signal(lo, hi, y0, y1), + _integral(sig2d, lo, hi), + cells._eps_sig, + f"signal {lo}..{hi}", + ) + checked += 1 + return checked + + +def _assert_same_sign(got, want, what): + """The gates branch on sign, so agreeing to 1e-9 is not the same as agreeing.""" + if (got < 0) != (want < 0): + raise AssertionError( + f"binning_dp parity failure on the sign of {what}: table says {got!r}, " + f"ROOT says {want!r}. _bin_passes rejects a bin whose background is " + "negative, so these two would not produce the same binning even though the " + "magnitudes agree. The zero-snapping in Cells is what should have prevented " + "this." + ) + + +def _assert_close(got, want, atol, what): + """Agreement to within the rounding the prefix sums can actually produce.""" + if abs(got - want) > atol: + raise AssertionError( + f"binning_dp parity failure on {what}: table says {got!r}, ROOT says " + f"{want!r}; they differ by {abs(got - want):.3e}, which exceeds the " + f"{atol:.3e} this quantity's cancellation can account for. The cell tables " + "do not reproduce the integrals the gates are defined on, so any binning " + "derived from them would be cut in the wrong places." + ) + + +NEG = -np.inf + + +def partition_dp(score, valid, max_parts): + """Best partition of [0..n-1] into exactly k contiguous valid ranges, for every k. + + Returns (values, partitions): values[k] is the total score of the best k-part + partition, or -inf if there is none, and partitions[k] is that partition as a list of + (lo, hi) index pairs. Index 0 is unused in both so that k reads as the part count. + + Every k is returned rather than just the requested one, because the caller needs the + whole curve and it costs nothing to keep: the bin count is chosen by comparing what + the partition scores at each count, and trim_by_marginal_gain() gives a bin back by + comparing values[k] against values[k-1]. Computing them one at a time would re-run the + same DP. + + Exact, not greedy. The score is additive over the parts -- Asimov Z^2 is what adds + across independent counting bins -- so the optimal k-part partition of a prefix is + built from an optimal (k-1)-part partition of a shorter prefix, and the usual + interval DP applies. O(max_parts * n^2). + + `valid` is a mask, not a penalty: an invalid range is unreachable rather than merely + expensive, so a partition containing one cannot be returned at any score. That is what + keeps the background gates hard constraints instead of preferences. + """ + n = score.shape[0] + max_parts = max(1, min(max_parts, n)) + masked = np.where(valid, score, NEG) + + dp = np.full((max_parts + 1, n + 1), NEG) + arg = np.full((max_parts + 1, n + 1), -1, dtype=int) + dp[0, 0] = 0.0 + for k in range(1, max_parts + 1): + for e in range(k, n + 1): + # candidate: the last part is [s .. e-1], the rest is a (k-1)-part prefix + cand = dp[k - 1, :e] + masked[:e, e - 1] + s = int(np.argmax(cand)) + if cand[s] > NEG: + dp[k, e] = cand[s] + arg[k, e] = s + + values = [NEG] * (max_parts + 1) + partitions = [None] * (max_parts + 1) + for k in range(1, max_parts + 1): + if dp[k, n] == NEG: + continue + values[k] = float(dp[k, n]) + parts = [] + e = n + for kk in range(k, 0, -1): + s = int(arg[kk, e]) + parts.append((s, e - 1)) + e = s + partitions[k] = list(reversed(parts)) + return values, partitions + + +def _self_test(trials=400, seed=0): + """Check partition_dp against exhaustive enumeration. + + It is a short dynamic program whose failure mode is silence: a subtly wrong + recurrence still returns a partition, and the binning it produces still looks like a + binning. On problems small enough to enumerate completely there is no reason to + settle for a plausible answer, so this compares against every partition there is. + + Pure numpy, no ROOT and no input files, so it runs anywhere: + python3 -m StatInference.common.binning_dp + The half of this module that does need ROOT is checked by parity_check() instead. + """ + import itertools + + def brute_partitions(n, k): + for cuts in itertools.combinations(range(1, n), k - 1): + bounds = (0,) + cuts + (n,) + parts = [(bounds[i], bounds[i + 1] - 1) for i in range(k)] + yield parts + + rng = np.random.default_rng(seed) + checked = 0 + for _ in range(trials): + n = int(rng.integers(2, 9)) + score = rng.normal(size=(n, n)) + valid = rng.random((n, n)) > rng.uniform(0.0, 0.6) + for a in range(n): + for b in range(n): + if b < a: + valid[a, b] = False + score[a, b] = 0.0 + values, partitions = partition_dp(score, valid, n) + for k in range(1, n + 1): + best = NEG + for part in brute_partitions(n, k): + if all(valid[a, b] for a, b in part): + total = sum(score[a, b] for a, b in part) + best = max(best, total) + got = values[k] + if best == NEG: + assert got == NEG, f"found a {k}-part partition where none exists" + else: + assert ( + abs(got - best) == 0 or abs(got - best) < 1e-9 + ), f"partition_dp n={n} k={k}: {got} != optimum {best}" + part = partitions[k] + assert len(part) == k + assert ( + part[0][0] == 0 and part[-1][1] == n - 1 + ), "does not cover the axis" + assert all( + part[i][1] + 1 == part[i + 1][0] for i in range(k - 1) + ), "parts are not contiguous" + assert all(valid[a, b] for a, b in part), "returned an invalid range" + checked += 1 + + return checked + + +def binning_objective(cells, slices, mode): + """(total Z^2, per-slice Z^2) for a finished binning. + + The quantity the search maximises, recomputed from the binning that was actually + written rather than carried out of the optimiser. That makes it meaningful for the + greedy strategy too -- which never computes it -- so the two can be compared from + their binning.json alone, and it means a discrepancy between what the optimiser + thought it achieved and what the shapes contain shows up as a discrepancy rather + than going unnoticed. + + Note this is evaluated on the ranges as recorded, i.e. after extend_outer_edges has + pushed the outermost ones into under/overflow, so it describes the shapes on disk. + """ + per_slice = [] + for sl in slices: + if "y_range" in sl: # HME box: selection on y, bins along x + ylo, yhi = sl["y_range"] + per_slice.append( + sum(cells.score(a, b, ylo, yhi, mode) for a, b in sl["x_ranges"]) + ) + else: # DNN slice: selection on x, bins along y + xlo, xhi = sl["x_range"] + per_slice.append( + sum(cells.score(xlo, xhi, a, b, mode) for a, b in sl["y_ranges"]) + ) + return sum(per_slice), per_slice + + +def _ranges_from_column(column, n, include_outer=False): + """M[i, j] = the rectangle sum for bins (i+1)..(j+1), from one column of a prefix sum. + + For a fixed x window, P[x1+1, :] - P[x0, :] is the running sum along y, and every + y range is a difference of two of its entries -- so the whole (n, n) table of + candidate ranges is one outer subtraction rather than n^2 lookups. + + With include_outer the first and last bins swallow the under- and overflow, which is + what extend_outer_edges() does to the outermost ranges before they are written. A + search that scores the ranges without them is optimising something slightly different + from what ends up in the datacard -- 0.11% of the background on the DNN axis of these + shapes, small but not nothing, and concentrated entirely in the two outermost slices. + """ + lo = column[1 : n + 1].copy() + hi = column[2 : n + 2].copy() + if include_outer: + lo[0] = column[0] + hi[-1] = column[n + 2] + return hi[None, :] - lo[:, None] + + +def _neff_matrix(value, error): + """effective_entries() over arrays, with its two degenerate cases kept. + + error <= 0 means the yield carries no MC uncertainty at all: infinite effective + entries if there is something there, none if there is not. Taking (v/e)^2 blindly + would make the first case a division by zero and the second a nan, and a nan + compares false against every threshold -- which would silently reject exactly the + empty bins the scalar gate accepts. + """ + with np.errstate(divide="ignore", invalid="ignore"): + ratio = np.where(error > 0, value / np.where(error > 0, error, 1.0), 0.0) ** 2 + degenerate = np.where(value > 0, np.inf, 0.0) + return np.where(error > 0, ratio, degenerate) + + +def _asimov_z2_matrix(s, b, b_err): + """asimov_significance(...)**2 over arrays, matching the scalar branch for branch.""" + out = np.zeros_like(s) + live = (b > 0) & (s > 0) + if not live.any(): + return out + var = b_err**2 + no_err = live & (b_err <= 0) + if no_err.any(): + ss, bb = s[no_err], b[no_err] + out[no_err] = np.maximum(2.0 * ((ss + bb) * np.log1p(ss / bb) - ss), 0.0) + with_err = live & (b_err > 0) + if with_err.any(): + ss, bb, vv = s[with_err], b[with_err], var[with_err] + term1 = (ss + bb) * np.log(((ss + bb) * (bb + vv)) / (bb * bb + (ss + bb) * vv)) + term2 = (bb * bb / vv) * np.log1p(vv * ss / (bb * (bb + vv))) + out[with_err] = np.maximum(2.0 * (term1 - term2), 0.0) + return out + + +def _sb_z2_matrix(s, b, b_err): + """significance(..., mode="sb")**2 over arrays: S^2 / (B + sigma_B^2).""" + denom = b + b_err**2 + ok = (b > 0) & (denom > 0) + with np.errstate(divide="ignore", invalid="ignore"): + return np.where(ok, s**2 / np.where(ok, denom, 1.0), 0.0) + + +def trim_by_marginal_gain(values_per_slice, counts, threshold): + """Give back every bin whose last split bought less than `threshold`. + + With Z^2 very nearly monotone in bin count, a search told only to maximise it will + always spend whatever it is given. So the binding constraint has to be the value of a + bin, not the size of a pot. A bin is kept only if the split that created it raised + this category's Z^2 by at least `threshold` -- expressed as a fraction of the + category's own achievable total, so it means "worth something on the scale that + matters here" rather than a yield in events. A bin that buys nothing is + MC-statistical exposure bought for nothing; this is what declines it. + + Stepping down one bin at a time and stopping at the first marginal gain that clears + the threshold is exact for a concave curve, which Z^2 in bin count is: each extra + split has less left to separate. It is not assumed -- a curve that is not concave + simply stops at its first qualifying step, which is still a count whose last bin paid + for itself. + """ + if threshold <= 0: + return list(counts) + trimmed = [] + for values, n in zip(values_per_slice, counts): + while n > 1 and values[n] - values[n - 1] < threshold: + n -= 1 + trimmed.append(n) + return trimmed + + +def _axis_ranges(prefix, a, b, n, axis, include_outer=False): + """All ranges along `axis`, at a fixed window [a, b] on the other one.""" + if axis == "y": + column = prefix[b + 1] - prefix[a] + else: + column = prefix[:, b + 1] - prefix[:, a] + return _ranges_from_column(column, n, include_outer) + + +def box_tables(cells, y0, y1, exempt, knobs, mode): + """(score, valid) over every candidate bin of the *sliced-on-y* layout. + + The window is a cut on y -- an HME box around the resonance -- and the bins run along + x, the DNN. The gates are _bin_passes() and the figure of merit significance(), as + whole-array operations: the scalar path -- building a yields dict and an errors dict + per candidate and calling _bin_passes -- costs about 80 microseconds a candidate, and + there are nx(nx+1)/2 of them per box, for every box the search tries. + + This is the one place that restates the gate rather than calling it, which is a real + risk: two spellings of the same rule are two rules as soon as one of them is edited. + It is contained by check_box_mask(), which re-evaluates the canonical _bin_passes on a + random sample of candidates and raises on the first disagreement -- so the fast path + is checked against the slow one on real inputs, on every run, rather than at review + time. + + The bins tile the whole x axis including its under/overflow, because everything inside + the box is kept and binned. The box itself is a selection: what falls outside it is + discarded, not swept into an outer bin, which is the difference between a cut and a + slice. + + `exempt` is the set of processes excused from the per-bin floors, judged once over + the whole box -- see minor_backgrounds(). Positivity is never excused -- see + _bin_passes. + """ + nx = cells.nx + ranges = lambda prefix: _axis_ranges(prefix, y0, y1, nx, "x", True) + + sig = ranges(cells._sig) + b_tot = ranges(cells._bkg_tot) + b_tot = np.where(np.abs(b_tot) < cells._eps_tot, 0.0, b_tot) + v_tot = ranges(cells._var_tot) + err = np.sqrt(np.maximum(v_tot, 0.0)) + + valid = np.triu(np.ones((nx, nx), dtype=bool)) + if knobs["min_bin_bkg_neff"] > 0: + valid &= _neff_matrix(b_tot, err) >= knobs["min_bin_bkg_neff"] + for name in cells.names: + b_p = ranges(cells._bkg[name]) + b_p = np.where(np.abs(b_p) < cells._eps_bkg[name], 0.0, b_p) + valid &= b_p >= 0 + if name in exempt: + continue + magnitude = b_p > knobs["min_bin_bkg_each"] + required = knobs["min_bin_bkg_each_neff"] + if required > 0: + v_p = ranges(cells._var[name]) + magnitude &= _neff_matrix(b_p, np.sqrt(np.maximum(v_p, 0.0))) >= required + valid &= magnitude + score = ( + _asimov_z2_matrix(sig, b_tot, err) + if mode == "asimov" + else _sb_z2_matrix(sig, b_tot, err) + ) + return np.where(valid, score, 0.0), valid + + +def check_box_mask( + cells, + y0, + y1, + valid, + score, + exempt, + bin_passes, + mode, + n_samples=100, + seed=0, +): + """The cross-check of box_tables() against the canonical scalar gate. + + Cheap -- a hundred calls against the thousands the tables replace -- and it is what + makes the fast path safe to trust: the tables are re-derived on a random sample with + _bin_passes and significance(), so if the two spellings of the rule ever disagree + about a candidate bin the run stops here instead of quietly producing a different + binning. + """ + nx = cells.nx + rng = np.random.default_rng(seed) + for _ in range(n_samples): + a = int(rng.integers(1, nx + 1)) + b = int(rng.integers(a, nx + 1)) + # the outermost bins carry the x under/overflow, because box_tables builds them + # with include_outer -- the reference has to ask the same question + lo = 0 if a == 1 else a + hi = nx + 1 if b == nx else b + want = bool(bin_passes(lo, hi, y0, y1, exempt)) + got = bool(valid[a - 1, b - 1]) + if got != want: + raise AssertionError( + f"binning_dp box gate disagrees with _bin_passes on DNN bins {a}..{b} of " + f"the HME box {y0}..{y1}: fast path says {'valid' if got else 'invalid'}, " + f"_bin_passes says {'valid' if want else 'invalid'}." + ) + if want: + expected = cells.score(lo, hi, y0, y1, mode) + if abs(expected - score[a - 1, b - 1]) > 1e-9 * max(abs(expected), 1.0): + raise AssertionError( + f"binning_dp box score disagrees with significance() on DNN bins " + f"{a}..{b} of box {y0}..{y1}: {score[a - 1, b - 1]!r} vs {expected!r}." + ) + + +if __name__ == "__main__": + n_part = _self_test() + print(f"partition_dp: {n_part} (n, k) cases vs exhaustive enumeration -- OK") diff --git a/law/CreateDatacardsTask.py b/law/CreateDatacardsTask.py index 1b09c9a..386a319 100644 --- a/law/CreateDatacardsTask.py +++ b/law/CreateDatacardsTask.py @@ -143,13 +143,42 @@ def run(self): base_dir_local.abspath, config, local_output.abspath ) - def plot_variable(self, variable): + @staticmethod + def _surviving_axis(base_dir): + """Which axis of the 2D input the rebinned shapes are binned along: 0 = x, 1 = y. + + Read off the binning record the preprocess step wrote beside the shapes. A DNN + slice selects on x and bins y; an HME box (a slice with a "y_range") selects on y + and bins x. With no record the shapes are the DNN-sliced kind, whose surviving + axis is y -- which is also what this assumed before boxes existed. + """ + import json + + path = os.path.join(base_dir, "binning.json") + try: + with open(path) as f: + record = json.load(f) + except (OSError, ValueError): + return 1 + node = record.get("binning", {}) + # era -> mass -> channel -> category -> {"slices": [...]} + while isinstance(node, dict) and "slices" not in node: + if not node: + return 1 + node = next(iter(node.values())) + slices = node.get("slices") if isinstance(node, dict) else None + return 0 if slices and "y_range" in slices[0] else 1 + + def plot_variable(self, variable, axis=1): """The variable whose axis the shapes are binned along, for histograms.yaml. - For input a 2D->1D rebinning produced, that is the y variable of the 2D entry -- - histograms.yaml records both as ``var_list: [x, y]``, so the axis metadata the - plotter needs is already described and does not have to be restated here. For - input that was always 1D, the variable is its own answer. + For input a 2D->1D rebinning produced, that is one of the two variables of the 2D + entry -- histograms.yaml records both as ``var_list: [x, y]``, so the axis + metadata the plotter needs is already described and does not have to be restated + here. `axis` picks which, from _surviving_axis(): the DNN-sliced shapes keep y + (HME), an HME box keeps x (the DNN). Taking y unconditionally labelled every box + panel "Deep HME Mass" under DNN-score bins. For input that was always 1D, the + variable is its own answer. """ try: import FLAF.Common.Setup as Setup @@ -157,7 +186,7 @@ def plot_variable(self, variable): hists = Setup.Setup(self.ana_path(), self.period, self.version).hists var_list = hists[variable].get("var_list") if var_list and len(var_list) > 1: - return var_list[1] + return var_list[axis] except Exception as e: print(f"Warning: no var_list for {variable} ({e}); plotting it as itself") return variable @@ -332,7 +361,7 @@ def plot_rebinned_shapes(self, base_dir, config, output_dir): # binning -- so no --rebin, which would coarsen it back to the histograms.yaml # grid and undo the whole point of the rebinning. "--var", - self.plot_variable(variable), + self.plot_variable(variable, self._surviving_axis(base_dir)), # The plotted shapes are the sum over every sub-era, so the label has to # name the combination: HistPlotter reads config/plot/.yaml for the # luminosity, and a sub-era's file states that sub-era's luminosity alone. From 4be3d0be4243abba3a6082890306a57d2a79f56d Mon Sep 17 00:00:00 2001 From: Devin Date: Sat, 26 Sep 2026 02:06:14 +0200 Subject: [PATCH 06/12] Keep the helpers in rebin_2d.py and name the strategy y_window The helpers moved to common/binning_core.py were moved to share them with a 1D binner that is not part of this change, so they go back where they were on main, unchanged except _slice_passes, which keeps the positivity rule. binning_dp.py no longer imports from the binner: the one figure of merit it needs is passed to build_cells, as the gate already is passed to check_window_mask. The development-time checks (parity_check, the partition self-test) and the Cells methods only they used are dropped; check_window_mask still re-derives a sample of every window's tables against _bin_passes and significance() on each run. The method is not specific to HH->bbWW's HME: it cuts one window on y and bins x inside it, whatever the two axes hold. So the strategy is y_window (was hme_box), its knobs are window_stride, window_refine and window_signal_quantile, and the prose describes x and y with HH->bbWW as the example. greedy stays the default. Checked against Binning_2D_DP on Run3_Early, v2605a_DNNOutputs_v3, MX 500 and 1000: the rebinned shapes are identical for y_window (188928 histograms) and greedy (755712), and binning.json is identical outside the knobs block, the reference's old names mapped to the new. Co-Authored-By: Claude Opus 5.5 (1M context) --- bin_opt_2d/rebin_2d.py | 583 ++++++++++++++++++++++++++++++++----- common/binning_core.py | 480 ------------------------------ common/binning_dp.py | 260 +++-------------- law/CreateDatacardsTask.py | 16 +- 4 files changed, 557 insertions(+), 782 deletions(-) delete mode 100644 common/binning_core.py diff --git a/bin_opt_2d/rebin_2d.py b/bin_opt_2d/rebin_2d.py index 12df51e..a71abd0 100644 --- a/bin_opt_2d/rebin_2d.py +++ b/bin_opt_2d/rebin_2d.py @@ -3,7 +3,6 @@ import math import os import sys - import numpy as np import yaml @@ -17,35 +16,13 @@ __package__ = pkg_dir_name from StatInference.common.tools import importROOT, CategoryNaming -from StatInference.common.param_parse import applyParameters - -# The dimension-agnostic half of the binning: slice finding, the gates it applies, -# the figures of merit, and the histogram plumbing. See that module's docstring for -# what stays here and why. -from StatInference.common.binning_core import ( - BINNING_JSON, - SIGNIFICANCE_MODES, - _bkg_errors, - _bkg_yields, - _detach, - _total_bkg_error, - bin_budget, - bin_edges, - effective_entries, - extend_outer_edges, - find_slices, - get_hist, - load_config, - lookup_frozen, - minor_backgrounds, - open_input_file, - sum_hists, -) +from StatInference.common.param_parse import extractParameters, applyParameters +from StatInference.dc_make.model import Model from StatInference.common.binning_dp import ( binning_objective, build_cells, - box_tables, - check_box_mask, + window_tables, + check_window_mask, partition_dp, trim_by_marginal_gain, ) @@ -53,12 +30,20 @@ ROOT = importROOT() +# Written inside each era's own output directory, beside the shapes it describes, and +# read back by --binning to replay them. +BINNING_JSON = "binning.json" + # The knobs that decide the binning, with the values a configuration gets when it does # not say otherwise. They belong to an analysis rather than to this file -- what a slice # needs to be worth keeping depends on the sample sizes and the selection, and a test # configuration wants them at zero -- so every production should state its own in a # binning yaml -- see config/Datacards/binning_2d.yaml in HH_bbWW for an annotated # set. This file ships no configuration of its own, only these fallbacks. +# The figures of merit significance() implements. Named here so the yaml path can be +# checked against the same list the command line's choices= uses. +SIGNIFICANCE_MODES = ("sb", "asimov") + BINNING_DEFAULTS = { "n_slices": 4, "category_pattern": None, # None -> CategoryNaming's own neutral default @@ -75,33 +60,33 @@ "min_bkg_frac": 0.05, "min_signal": 0.5, "significance_mode": "asimov", - # How the edges are found. "greedy" is the original: each slice boundary is the best - # one taken on its own, and the mass bins inside a slice are equal-signal quantiles - # merged until the gates pass. "hme_box" does the reverse: one HME window per - # category, with the DNN binned inside it by an exact partition search -- see - # _discover_hme_box(). Default stays "greedy" so adding this moves nothing until a - # configuration asks for it. + # How the edges are found. "greedy" is the original: x is cut into slices, each + # boundary the best one taken on its own, and the y bins inside a slice are + # equal-signal quantiles merged until the gates pass. "y_window" does the reverse: one + # window on y per category, with x binned inside it by an exact partition search -- + # see _discover_y_window(). Default stays "greedy" so adding this moves nothing until + # a configuration asks for it. "strategy": "greedy", - # hme_box only: a bin must earn its place. The split that creates it has to raise + # y_window only: a bin must earn its place. The split that creates it has to raise # this category's Z^2 by at least this fraction of the category's own achievable # total, or it is given back. 0.0 leaves the count to max_bins_per_slice alone. "dp_min_bin_gain": 0.0, - # hme_box only: the search grid for the box edges. A full scan is ny(ny+1)/2 windows + # y_window only: the search grid for the window edges. A full scan is ny(ny+1)/2 windows # and every one of them carries its own dynamic program over the binned axis, so the - # edges are scanned on a stride and then refined within +-box_refine of the winner. - "box_stride": 6, - "box_refine": 5, - # hme_box only: the box is searched only where the signal is. A window that starts - # beyond the resonance adds background and no signal; scanning it is wasted time, and - # including it in the grid costs resolution where the answer actually lies. - "box_signal_quantile": 0.001, + # edges are scanned on a stride and then refined within +-window_refine of the winner. + "window_stride": 6, + "window_refine": 5, + # y_window only: the window is searched only where the signal is. A window that + # starts beyond the signal adds background and no signal; scanning it is wasted time, + # and including it in the grid costs resolution where the answer actually lies. + "window_signal_quantile": 0.001, } # Named here so the yaml path can be checked against the same list the command line's # choices= uses -- significance_mode already learned this lesson: it is read by a # function that silently treats anything unrecognised as a default, so a typo in the # yaml quietly changed the answer instead of failing. -BINNING_STRATEGIES = ("greedy", "hme_box") +BINNING_STRATEGIES = ("greedy", "y_window") def load_binning_config(path, overrides=None): @@ -141,6 +126,291 @@ def load_binning_config(path, overrides=None): return knobs +def lookup_frozen(frozen_binning, era, mass, channel, category): + """The recorded binning for one channel/category, or None if it was skipped then. + + A category absent from the record was skipped by the run that wrote it (too little + signal, or a background missing from a discovery era), so it is skipped again rather + than quietly re-optimised -- a replay that rebinned some categories and froze others + would be neither the old binning nor a new one. + """ + return ( + frozen_binning.get("binning", {}) + .get(era, {}) + .get(str(mass), {}) + .get(channel, {}) + .get(category) + ) + + +def load_config(config_path): + with open(config_path, "r") as f: + cfg = yaml.safe_load(f) + model = Model.fromConfig(cfg["model"]) + channels = cfg["channels"] + # Taken verbatim here; run() reduces them to the base categories the 2D input is + # actually keyed by, once it has the pattern to do it with. + categories = list(cfg["categories"]) + + # Several processes may carry is_signal (e.g. the bbWW(2l) and bbtautau decay + # modes of the same resonance, both scaled by the same signal strength). They + # are summed to form the discovery signal that steers the slice + # boundaries, so the binning is optimised for the total signal in the fit. + signal_hist_names = [] + mass_values = None + background_entries = [] + for entry in cfg["processes"]: + if type(entry) == str: + background_entries.append((entry, entry, [])) + continue + if entry.get("is_data", False): + continue + base_name = entry["process"] + hist_name = entry.get("hist_name", base_name) + if entry.get("is_signal", False): + signal_hist_names.append(hist_name) + if mass_values is None: + mass_values = entry["param_values"] + elif list(entry["param_values"]) != list(mass_values): + raise RuntimeError( + f"Signal {hist_name} has param_values {entry['param_values']}, " + f"which differ from {mass_values}; every signal must be defined " + "at the same mass points" + ) + else: + background_entries.append((base_name, hist_name, entry.get("channels", []))) + + if not signal_hist_names: + raise RuntimeError("No signal process found in config") + + return { + "model": model, + "channels": channels, + "categories": categories, + "era_groups": cfg.get("era_groups", {}), + # Declared by a configuration that reads sliced shapes, so that run() can check + # it against the pattern the shapes are actually written with. None when the + # configuration names no pattern, and then there is nothing to check. + "category_pattern": cfg.get("category_pattern"), + "signal_hist_name_patterns": signal_hist_names, + "signal_param_name": extractParameters(signal_hist_names[0])[0], + "mass_values": mass_values, + "background_entries": background_entries, + } + + +def open_input_file(input_dir, model, era, mass, param_name): + file_name = model.getInputFileName(era, {param_name: mass}) + full_path = os.path.join(input_dir, file_name) + f = ROOT.TFile.Open(full_path, "READ") + if f is None or f.IsZombie(): + raise RuntimeError(f"Cannot open file {full_path}") + return f + + +def _detach(h): + """Detach a histogram from its TFile and hand ownership to Python. + + SetDirectory(0) alone makes the histogram survive the file's close, but + leaves it owned by nobody -- so every 2D histogram read here (one per + systematic variation, per category, per mass) leaked for the lifetime of the + process. Across the ten-mass loop that grew without bound and got the job + SIGKILLed partway through the final mass, leaving a truncated ROOT file. + SetOwnership makes the object die with its last Python reference. + """ + h.SetDirectory(0) + ROOT.SetOwnership(h, True) + return h + + +def get_hist(f, path): + h = f.Get(path) + # TFile.Get() on a fully-missing nested path can return a PyROOT wrapper + # around a null C++ pointer, which is not `is None` but is falsy. + if not h: + return None + return _detach(h) + + +def sum_hists(hists): + total = None + for h in hists: + if h is None: + continue + if total is None: + total = _detach(h.Clone()) + else: + total.Add(h) + return total + + +def _integral(hist, lo, hi): + return ( + hist.Integral(lo, hi, 0, -1) + if hist.GetDimension() == 2 + else hist.Integral(lo, hi) + ) + + +def _integral_and_error(hist, lo, hi): + """(yield, MC statistical error) over [lo, hi].""" + err = array.array("d", [0.0]) + if hist.GetDimension() == 2: + value = hist.IntegralAndError(lo, hi, 0, -1, err) + else: + value = hist.IntegralAndError(lo, hi, err) + return value, err[0] + + +def _bkg_yields(bkg_hists_by_name, lo, hi): + """For each background, the summed yield across all discovery eras (combined + statistics). Individual eras are not required to individually clear the + threshold -- use allow_negative_bins_within_error (maker.py, per process/ + category) for backgrounds/categories where a specific era can land negative + in a bin that's fine once combined.""" + return { + name: sum(_integral(h, lo, hi) for h in hists) + for name, hists in bkg_hists_by_name.items() + } + + +def _bkg_errors(bkg_hists_by_name, lo, hi): + """Per background, the MC statistical error on the summed yield (eras added in + quadrature).""" + return { + name: math.sqrt(sum(_integral_and_error(h, lo, hi)[1] ** 2 for h in hists)) + for name, hists in bkg_hists_by_name.items() + } + + +def _total_bkg_error(bkg_hists_by_name, lo, hi): + errors = _bkg_errors(bkg_hists_by_name, lo, hi) + return math.sqrt(sum(e**2 for e in errors.values())) + + +def effective_entries(value, error): + """(yield / error)^2 -- the unweighted MC event count a weighted yield is worth. + + A DY slice holding a couple of very-high-weight aMC@NLO events can carry a + sizeable yield with an effective count far below 1, i.e. a background estimate + that is statistically compatible with almost anything. + """ + if error <= 0: + return float("inf") if value > 0 else 0.0 + return (value / error) ** 2 + + +def asimov_significance(s, b, b_err=0.0): + """Median discovery significance for a counting experiment (Cowan et al., + arXiv:1007.1727), with the background uncertainty folded in when given. + + Reduces to sqrt(2*((s+b)*ln(1+s/b) - s)) for b_err = 0. Unlike S/sqrt(B) this + stays valid when b is O(1), which is exactly the regime the high-mass boosted + slices live in -- there S/sqrt(B) reports significances of ~30 on ~1.5 + background events, and drives the slice boundary on that basis. + """ + if b <= 0 or s <= 0: + return 0.0 + if b_err <= 0: + return math.sqrt(max(2.0 * ((s + b) * math.log1p(s / b) - s), 0.0)) + var = b_err**2 + term1 = (s + b) * math.log(((s + b) * (b + var)) / (b * b + (s + b) * var)) + term2 = (b * b / var) * math.log1p(var * s / (b * (b + var))) + return math.sqrt(max(2.0 * (term1 - term2), 0.0)) + + +def significance(s, b, b_err=0.0, mode="sb"): + """Figure of merit steering the slice boundaries. + + mode="sb": S / sqrt(B + sigma_B^2). Folding sigma_B into plain S/sqrt(B) + stops a downward fluctuation of a statistics-starved background + (DY in the b-tagged muMu slices, effective MC count <1) from + inflating the apparent significance. Still assumes the Gaussian + regime, which breaks down for B of order a few. + mode="asimov": the Poisson-correct Asimov significance, valid at low B. + """ + if b is None or b <= 0: + return 0.0 + if mode == "asimov": + return asimov_significance(s, b, b_err) + denominator = b + b_err**2 + if denominator <= 0: + return 0.0 + return s / (denominator**0.5) + + +def _slice_passes( + yields, + min_sum, + total_error=None, + min_neff=0.0, + errors=None, + min_each=0.0, + min_proc_neff=0.0, + exempt=(), +): + """Slice validity: the summed background must clear min_sum, and be known + to better than min_neff effective MC entries -- a window whose background is + statistically undetermined must not be selectable at all. + + The per-process arms (min_each, min_proc_neff) are what _bin_passes already + does on the mass axis, applied here as well. Testing only the sum hides an + individual background that has fluctuated negative behind a large, well + measured neighbour: muMu/SR/res2b at m500 selected a slice holding + DY = -5.87 +- 8.92 (N_eff 0.43) because the summed background there was + +39.9 +- 9.0 (N_eff 22), clearing both summed tests comfortably. That slice + cannot be turned into a datacard at all -- a negative DY integral is rejected + outright by resolveNegativeBins -- so the binning has to not select it in the + first place. Worse, the significance being maximised is S/sqrt(B + sigma_B^2), + which a downward fluctuation *increases* by lowering B, so such a window is + mildly preferred rather than merely tolerated. + + `exempt` comes from minor_backgrounds() judged over the whole sliced-axis range of the + category, never over the candidate window: a background that is negligible in + this category must not be able to veto every boundary, and judging it inside + the window under test is circular (a background that is exactly zero there is + trivially below any fraction of the total). + """ + total = sum(yields.values()) + if total <= min_sum: + return False + if min_neff > 0 and total_error is not None: + if effective_entries(total, total_error) < min_neff: + return False + for name, value in yields.items(): + # As in _bin_passes: exemption covers the magnitude floor, never the sign. + if value < 0: + return False + if name in exempt: + continue + if min_each > 0 and value <= min_each: + return False + if min_proc_neff > 0 and errors is not None: + if effective_entries(value, errors.get(name, 0.0)) < min_proc_neff: + return False + return True + + +def minor_backgrounds(bkg_hists_by_name, lo, hi, min_frac): + """Backgrounds negligible over the *whole* [lo, hi] range, which may therefore + be exempted from the per-bin min_each floor. + + This must be judged once over the full slice, never inside the candidate bin + under test. Evaluating the fraction per bin is circular: a background that is + exactly zero in that bin is trivially below any fraction of the bin total, so + it is always exempted -- which is precisely the case the floor exists to + catch. That let a bin through with TT = 0 in a slice where TT is 94% of the + background (m300, eMu, res2b_dnn2). + """ + if min_frac <= 0: + return set() + yields = _bkg_yields(bkg_hists_by_name, lo, hi) + total = sum(yields.values()) + if total <= 0: + return set() + return {name for name, value in yields.items() if value < min_frac * total} + + def _bin_passes( yields, min_each, @@ -199,6 +469,114 @@ def _bin_passes( return True +def grow_slice( + sig_hist, + bkg_hists_by_name, + right, + first_bin, + min_sum, + min_neff=0.0, + sig_mode="sb", + min_each=0.0, + min_proc_neff=0.0, + exempt=(), +): + """Among every candidate [left, right] with summed background > min_sum, + pick the one maximizing S/sqrt(B + sigma_B^2) -- not just the first one that + clears it. This is what actually drives where the cut lands; the content + floor is only a validity gate. Falls back to first_bin (best effort) if no + candidate clears the floor anywhere.""" + best_left = None + best_sig = -1.0 + for left in range(right, first_bin - 1, -1): + bkg_y = _bkg_yields(bkg_hists_by_name, left, right) + b_err = _total_bkg_error(bkg_hists_by_name, left, right) + bkg_e = ( + _bkg_errors(bkg_hists_by_name, left, right) if min_proc_neff > 0 else None + ) + if not _slice_passes( + bkg_y, + min_sum, + b_err, + min_neff, + bkg_e, + min_each, + min_proc_neff, + exempt, + ): + continue + s = _integral(sig_hist, left, right) + b = sum(bkg_y.values()) + sig = significance(s, b, b_err, sig_mode) + if sig > best_sig: + best_sig = sig + best_left = left + return best_left if best_left is not None else first_bin + + +def find_slices( + sig_hist, + bkg_hists_by_name, + n_slices, + first_bin, + last_bin, + min_sum, + min_neff=0.0, + sig_mode="sb", + min_each=0.0, + min_proc_neff=0.0, + min_frac=0.0, +): + """Split [first_bin, last_bin] into exactly `n_slices` ranges (fixed count, + required so every mass point shares the same category list), scanning right + to left. Each slice (except the final, leftover one) is placed to maximize + signal significance among boundaries with summed background > min_sum. + + Which backgrounds are minor enough to be exempt from the per-process floors is + decided once here, over the whole [first_bin, last_bin] range, and held fixed + for every candidate window -- see _slice_passes on why it cannot be re-judged + per window. + """ + exempt = ( + minor_backgrounds(bkg_hists_by_name, first_bin, last_bin, min_frac) + if (min_each > 0 or min_proc_neff > 0) + else set() + ) + slices = [] + right = last_bin + for slice_idx in range(n_slices): + if slice_idx == n_slices - 1 or right <= first_bin: + if right < first_bin: + # The axis ran out before n_slices boundaries were placed. Reported as + # None rather than invented: the previous code emitted (0,0), (-1,-1), + # (-2,-2) here to keep the count, and ROOT reads any range whose upper + # bound is below its lower one as the *whole axis including overflow* -- + # so those became copies of the full plane, and the same events entered + # the fit two or three times as independent categories. discover_binning + # rejects the category instead. + slices.append(None) + continue + slices.append((first_bin, right)) + right = first_bin - 1 + continue + left = grow_slice( + sig_hist, + bkg_hists_by_name, + right, + first_bin, + min_sum, + min_neff, + sig_mode, + min_each, + min_proc_neff, + exempt, + ) + slices.append((left, right)) + right = left - 1 + slices.reverse() + return slices + + def signal_quantile_ranges(sig_hist, n_bins, first_bin, last_bin): """Split [first_bin, last_bin] into `n_bins` ranges each holding an equal share of the signal. @@ -325,6 +703,39 @@ def find_bins( ) +def extend_outer_edges(ranges, full_lo, full_hi): + """Widen the first/last range to swallow under/overflow (bin 0 / nbins+1), + so no events are silently dropped at the extremes of the axis. + + Only real ranges are widened. Widening a None (an exhausted slice) would invent one, + and widening past a None would put the underflow on a range that is not the outermost + -- so a list containing any None is left alone and rejected by the caller. + """ + ranges = list(ranges) + if any(r is None for r in ranges): + return ranges + ranges[0] = (full_lo, ranges[0][1]) + ranges[-1] = (ranges[-1][0], full_hi) + return ranges + + +def bin_budget(bkg_hists_by_name, lo, hi, max_bins_per_slice, bkg_per_bin): + """How many bins this slice can actually afford. + + A fixed max_bins_per_slice is applied blind to slice content: the high-mass boosted + slices hold ~1.5 total background events and were still being split into 10 + bins, i.e. ~0.15 events per bin. Capping at B_slice / bkg_per_bin ties the + binning to the statistics that are really there. Returns max_bins_per_slice + unchanged when bkg_per_bin <= 0 (feature off). + """ + if bkg_per_bin <= 0: + return max_bins_per_slice + total = sum(_bkg_yields(bkg_hists_by_name, lo, hi).values()) + if total <= 0: + return 1 + return max(1, min(max_bins_per_slice, int(total / bkg_per_bin))) + + def _discover_greedy(sig2d, bkg2d_by_name, knobs): """bkg2d_by_name: {background_name: [hist per discovery era, ...]}. The list is usually a single era's own histogram (standalone limit) or all of a @@ -427,19 +838,19 @@ def _best_count(values, cap): return best_n -def _box_value(cells, y0, y1, knobs, check=False): - """What one HME box is worth: the best 1D binning of the DNN inside it. +def _window_value(cells, y0, y1, knobs, check=False): + """What one window on y is worth: the best 1D binning of x inside it. - A box is scored by the binning it admits, not by its own integrated significance. + A window is scored by the binning it admits, not by its own integrated significance. Those are different questions and they pick different windows -- scored on its own, a - box is rewarded for swallowing as much signal as it can and the DNN shape inside it - plays no part, which is the opposite of what the fit will do with it. + window is rewarded for swallowing as much signal as it can and the shape along x + inside it plays no part, which is the opposite of what the fit will do with it. """ mode = knobs["significance_mode"] - # which backgrounds count as negligible is decided once, over the whole box -- see + # which backgrounds count as negligible is decided once, over the whole window -- see # minor_backgrounds() on why it cannot be re-judged inside each candidate bin exempt = cells.exempt(1, cells.nx, y0, y1, knobs["min_bkg_frac"]) - score, valid = box_tables(cells, y0, y1, exempt, knobs, mode) + score, valid = window_tables(cells, y0, y1, exempt, knobs, mode) if check: def bin_passes(a, b, ylo, yhi, exempt_set): @@ -453,7 +864,7 @@ def bin_passes(a, b, ylo, yhi, exempt_set): knobs["min_bin_bkg_each_neff"], ) - check_box_mask( + check_window_mask( cells, y0, y1, @@ -474,53 +885,55 @@ def bin_passes(a, b, ylo, yhi, exempt_set): return values[n], [(a + 1, b + 1) for a, b in partitions[n]] -def _discover_hme_box(sig2d, bkg2d_by_name, knobs): - """One HME window per category, with the DNN binned inside it. +def _discover_y_window(sig2d, bkg2d_by_name, knobs): + """One window on y per category, with x binned inside it. - The greedy strategy cuts the DNN into slices that each become a datacard category and - gives every slice a shape along HME. This one does the reverse and keeps only one - category: a single contiguous HME box around the resonance, and a 1D DNN distribution - inside it. Everything outside the box is dropped -- that is what makes it a cut rather - than a slice, and it is why the box has to earn its acceptance loss in purity. + The greedy strategy cuts x into slices that each become a datacard category and gives + every slice a shape along y. This one does the reverse and keeps only one category: a + single contiguous window on y around the signal, and a 1D distribution of x inside it. + Everything outside the window is dropped -- that is what makes it a cut rather than a + slice, and it is why the window has to earn its acceptance loss in purity. It suits a + y the signal peaks in (in HH->bbWW, x is the DNN score and y the HME mass), and it + keeps only one region of y, so a signal with two separated peaks would lose one. - Box and binning are chosen together; see _box_value. + Window and binning are chosen together; see _window_value. """ - cells = build_cells(sig2d, bkg2d_by_name) + cells = build_cells(sig2d, bkg2d_by_name, significance) ny = cells.ny - # where the signal actually is -- see box_signal_quantile + # where the signal actually is -- see window_signal_quantile cumulative = np.array( [cells.signal(0, cells.nx + 1, 1, j) for j in range(1, ny + 1)] ) total = cumulative[-1] if len(cumulative) else 0.0 if total <= 0: return None - q = knobs["box_signal_quantile"] + q = knobs["window_signal_quantile"] lo_edge = max(1, int(np.searchsorted(cumulative, q * total))) hi_edge = min(ny, int(np.searchsorted(cumulative, (1.0 - q) * total)) + 1) best_value, best_box = float("-inf"), None - stride = max(1, int(knobs["box_stride"])) + stride = max(1, int(knobs["window_stride"])) for a in range(lo_edge, hi_edge + 1, stride): for b in range(a, hi_edge + 1, stride): - value, _ = _box_value(cells, a, b, knobs) + value, _ = _window_value(cells, a, b, knobs) if value > best_value: best_value, best_box = value, (a, b) if best_box is None: return None - refine = max(0, int(knobs["box_refine"])) + refine = max(0, int(knobs["window_refine"])) a0, b0 = best_box for a in range(max(lo_edge, a0 - refine), min(hi_edge, a0 + refine) + 1): for b in range(max(a, b0 - refine), min(hi_edge, b0 + refine) + 1): - value, _ = _box_value(cells, a, b, knobs) + value, _ = _window_value(cells, a, b, knobs) if value > best_value: best_value, best_box = value, (a, b) y0, y1 = best_box - value, bins = _box_value(cells, y0, y1, knobs, check=True) + value, bins = _window_value(cells, y0, y1, knobs, check=True) if bins is None: return None - # A box touching an end of the axis takes that end's overflow with it: the events are - # inside the selection the box stands for, and dropping them would make the written + # A window touching an end of the axis takes that end's overflow with it: the events + # are inside the selection the window stands for, and dropping them would make the written # yield disagree with the cut the category claims to be. y0 = 0 if y0 <= 1 else y0 y1 = ny + 1 if y1 >= ny else y1 @@ -537,11 +950,27 @@ def discover_binning(sig2d, bkg2d_by_name, knobs): Returns None when the category cannot be binned, which the caller reports as a skip. """ - if knobs["strategy"] == "hme_box": - return _discover_hme_box(sig2d, bkg2d_by_name, knobs) + if knobs["strategy"] == "y_window": + return _discover_y_window(sig2d, bkg2d_by_name, knobs) return _discover_greedy(sig2d, bkg2d_by_name, knobs) +def bin_edges(y_axis, y_ranges): + """Physical edges of the discovered y ranges, for booking the output TH1. + + The rebinned shapes used to be booked as n_bins over [0, n_bins], which threw + the axis scale away and left every plot labelled by bin index. The ranges are + contiguous and ordered, so the edges are each range's low edge plus the last + range's upper edge. extend_outer_edges() pushes the outer ranges into the + underflow/overflow bins, which have no finite edge of their own -- those are + clamped to the axis limits. + """ + n = y_axis.GetNbins() + edges = [y_axis.GetBinLowEdge(max(lo, 1)) for lo, _ in y_ranges] + edges.append(y_axis.GetBinUpEdge(min(y_ranges[-1][1], n))) + return edges + + def mkdir_titled(directory, path, title): """mkdir a nested path, putting `title` on the leaf directory only. @@ -693,8 +1122,8 @@ def record_to_slices(record, x_axis, y_axis, where): def slice_parts(sl): """(selection range, bin ranges, which axis the selection is on) for one category. - Two layouts reach this point. The DNN-sliced strategy cuts x into categories and gives - each a shape along y: {"x_range", "y_ranges"}. The HME-box strategy cuts y and gives + Two layouts reach this point. The greedy strategy cuts x into categories and gives + each a shape along y: {"x_range", "y_ranges"}. The y_window strategy cuts y and gives the one category a shape along x: {"y_range", "x_ranges"}. Everything downstream -- the writer, the record, the objective -- needs the same three things from either, and asking here is what keeps that code from having to know which strategy ran. @@ -890,7 +1319,9 @@ def load2d(f, key): # number, but not a true one. objective = ( binning_objective( - build_cells(disc_sig, disc_bkg_by_name), slices, knobs["significance_mode"] + build_cells(disc_sig, disc_bkg_by_name, significance), + slices, + knobs["significance_mode"], ) if disc_bkg_by_name else None @@ -1173,8 +1604,8 @@ def run( default=None, choices=list(BINNING_STRATEGIES), help="how the edges are found: 'greedy' = each slice boundary taken on its own " - "and equal-signal quantiles inside it, 'hme_box' = one HME window per category " - "with the DNN binned inside it by an exact partition search", + "and equal-signal quantiles inside it, 'y_window' = one window on y per category " + "with x binned inside it by an exact partition search", ) parser.add_argument( "--significance-mode", diff --git a/common/binning_core.py b/common/binning_core.py deleted file mode 100644 index 5750181..0000000 --- a/common/binning_core.py +++ /dev/null @@ -1,480 +0,0 @@ -"""The parts of the shape rebinning that do not care how many axes a histogram has. - -Split out of bin_opt_2d/rebin_2d.py, which is where all of this was written and where -the comments' point of view comes from. It moved here unchanged so that anything else -that needs the *same* answers to the same questions -- what a slice has to contain to be -worth keeping, which processes are minor enough to exempt, which figure of merit decides -a boundary -- reads them from here, rather than from a second copy that would be a second -place for them to drift. - -`_integral()` branches on GetDimension(), so everything downstream of it reads a TH1 and -a TH2 alike. What stays behind in the binner is the part that knows what a "bin" is: the -2D driver treats a slice as a directory of mass bins. - -Nothing here is configuration. BINNING_DEFAULTS and load_binning_config() stay with the -binner, because their "unknown key raises" contract only works if the file that reads the -knobs states its own valid keys. -""" - -import array -import math -import os - -import yaml - -from StatInference.common.tools import importROOT -from StatInference.common.param_parse import extractParameters -from StatInference.dc_make.model import Model - -ROOT = importROOT() - - -# Written inside each era's own output directory, beside the shapes it describes, and -# read back by --binning to replay them. -BINNING_JSON = "binning.json" - -# The figures of merit significance() implements. Named here so the yaml path can be -# checked against the same list the command line's choices= uses. -SIGNIFICANCE_MODES = ("sb", "asimov") - - -def lookup_frozen(frozen_binning, era, mass, channel, category): - """The recorded binning for one channel/category, or None if it was skipped then. - - A category absent from the record was skipped by the run that wrote it (too little - signal, or a background missing from a discovery era), so it is skipped again rather - than quietly re-optimised -- a replay that rebinned some categories and froze others - would be neither the old binning nor a new one. - """ - return ( - frozen_binning.get("binning", {}) - .get(era, {}) - .get(str(mass), {}) - .get(channel, {}) - .get(category) - ) - - -def load_config(config_path): - with open(config_path, "r") as f: - cfg = yaml.safe_load(f) - model = Model.fromConfig(cfg["model"]) - channels = cfg["channels"] - # Taken verbatim here; run() reduces them to the base categories the 2D input is - # actually keyed by, once it has the pattern to do it with. - categories = list(cfg["categories"]) - - # Several processes may carry is_signal (e.g. the bbWW(2l) and bbtautau decay - # modes of the same resonance, both scaled by the same signal strength). They - # are summed to form the discovery signal that steers the slice - # boundaries, so the binning is optimised for the total signal in the fit. - signal_hist_names = [] - mass_values = None - background_entries = [] - for entry in cfg["processes"]: - if type(entry) == str: - background_entries.append((entry, entry, [])) - continue - if entry.get("is_data", False): - continue - base_name = entry["process"] - hist_name = entry.get("hist_name", base_name) - if entry.get("is_signal", False): - signal_hist_names.append(hist_name) - if mass_values is None: - mass_values = entry["param_values"] - elif list(entry["param_values"]) != list(mass_values): - raise RuntimeError( - f"Signal {hist_name} has param_values {entry['param_values']}, " - f"which differ from {mass_values}; every signal must be defined " - "at the same mass points" - ) - else: - background_entries.append((base_name, hist_name, entry.get("channels", []))) - - if not signal_hist_names: - raise RuntimeError("No signal process found in config") - - return { - "model": model, - "channels": channels, - "categories": categories, - "era_groups": cfg.get("era_groups", {}), - # Declared by a configuration that reads sliced shapes, so that run() can check - # it against the pattern the shapes are actually written with. None when the - # configuration names no pattern, and then there is nothing to check. - "category_pattern": cfg.get("category_pattern"), - "signal_hist_name_patterns": signal_hist_names, - "signal_param_name": extractParameters(signal_hist_names[0])[0], - "mass_values": mass_values, - "background_entries": background_entries, - } - - -def open_input_file(input_dir, model, era, mass, param_name): - file_name = model.getInputFileName(era, {param_name: mass}) - full_path = os.path.join(input_dir, file_name) - f = ROOT.TFile.Open(full_path, "READ") - if f is None or f.IsZombie(): - raise RuntimeError(f"Cannot open file {full_path}") - return f - - -def _detach(h): - """Detach a histogram from its TFile and hand ownership to Python. - - SetDirectory(0) alone makes the histogram survive the file's close, but - leaves it owned by nobody -- so every 2D histogram read here (one per - systematic variation, per category, per mass) leaked for the lifetime of the - process. Across the ten-mass loop that grew without bound and got the job - SIGKILLed partway through the final mass, leaving a truncated ROOT file. - SetOwnership makes the object die with its last Python reference. - """ - h.SetDirectory(0) - ROOT.SetOwnership(h, True) - return h - - -def get_hist(f, path): - h = f.Get(path) - # TFile.Get() on a fully-missing nested path can return a PyROOT wrapper - # around a null C++ pointer, which is not `is None` but is falsy. - if not h: - return None - return _detach(h) - - -def sum_hists(hists): - total = None - for h in hists: - if h is None: - continue - if total is None: - total = _detach(h.Clone()) - else: - total.Add(h) - return total - - -def _integral(hist, lo, hi): - return ( - hist.Integral(lo, hi, 0, -1) - if hist.GetDimension() == 2 - else hist.Integral(lo, hi) - ) - - -def _integral_and_error(hist, lo, hi): - """(yield, MC statistical error) over [lo, hi].""" - err = array.array("d", [0.0]) - if hist.GetDimension() == 2: - value = hist.IntegralAndError(lo, hi, 0, -1, err) - else: - value = hist.IntegralAndError(lo, hi, err) - return value, err[0] - - -def _bkg_yields(bkg_hists_by_name, lo, hi): - """For each background, the summed yield across all discovery eras (combined - statistics). Individual eras are not required to individually clear the - threshold -- use allow_negative_bins_within_error (maker.py, per process/ - category) for backgrounds/categories where a specific era can land negative - in a bin that's fine once combined.""" - return { - name: sum(_integral(h, lo, hi) for h in hists) - for name, hists in bkg_hists_by_name.items() - } - - -def _bkg_errors(bkg_hists_by_name, lo, hi): - """Per background, the MC statistical error on the summed yield (eras added in - quadrature).""" - return { - name: math.sqrt(sum(_integral_and_error(h, lo, hi)[1] ** 2 for h in hists)) - for name, hists in bkg_hists_by_name.items() - } - - -def _total_bkg_error(bkg_hists_by_name, lo, hi): - errors = _bkg_errors(bkg_hists_by_name, lo, hi) - return math.sqrt(sum(e**2 for e in errors.values())) - - -def effective_entries(value, error): - """(yield / error)^2 -- the unweighted MC event count a weighted yield is worth. - - A DY slice holding a couple of very-high-weight aMC@NLO events can carry a - sizeable yield with an effective count far below 1, i.e. a background estimate - that is statistically compatible with almost anything. - """ - if error <= 0: - return float("inf") if value > 0 else 0.0 - return (value / error) ** 2 - - -def asimov_significance(s, b, b_err=0.0): - """Median discovery significance for a counting experiment (Cowan et al., - arXiv:1007.1727), with the background uncertainty folded in when given. - - Reduces to sqrt(2*((s+b)*ln(1+s/b) - s)) for b_err = 0. Unlike S/sqrt(B) this - stays valid when b is O(1), which is exactly the regime the high-mass boosted - slices live in -- there S/sqrt(B) reports significances of ~30 on ~1.5 - background events, and drives the slice boundary on that basis. - """ - if b <= 0 or s <= 0: - return 0.0 - if b_err <= 0: - return math.sqrt(max(2.0 * ((s + b) * math.log1p(s / b) - s), 0.0)) - var = b_err**2 - term1 = (s + b) * math.log(((s + b) * (b + var)) / (b * b + (s + b) * var)) - term2 = (b * b / var) * math.log1p(var * s / (b * (b + var))) - return math.sqrt(max(2.0 * (term1 - term2), 0.0)) - - -def significance(s, b, b_err=0.0, mode="sb"): - """Figure of merit steering the slice boundaries. - - mode="sb": S / sqrt(B + sigma_B^2). Folding sigma_B into plain S/sqrt(B) - stops a downward fluctuation of a statistics-starved background - (DY in the b-tagged muMu slices, effective MC count <1) from - inflating the apparent significance. Still assumes the Gaussian - regime, which breaks down for B of order a few. - mode="asimov": the Poisson-correct Asimov significance, valid at low B. - """ - if b is None or b <= 0: - return 0.0 - if mode == "asimov": - return asimov_significance(s, b, b_err) - denominator = b + b_err**2 - if denominator <= 0: - return 0.0 - return s / (denominator**0.5) - - -def _slice_passes( - yields, - min_sum, - total_error=None, - min_neff=0.0, - errors=None, - min_each=0.0, - min_proc_neff=0.0, - exempt=(), -): - """Slice validity: the summed background must clear min_sum, and be known - to better than min_neff effective MC entries -- a window whose background is - statistically undetermined must not be selectable at all. - - The per-process arms (min_each, min_proc_neff) are what _bin_passes already - does on the mass axis, applied here as well. Testing only the sum hides an - individual background that has fluctuated negative behind a large, well - measured neighbour: muMu/SR/res2b at m500 selected a slice holding - DY = -5.87 +- 8.92 (N_eff 0.43) because the summed background there was - +39.9 +- 9.0 (N_eff 22), clearing both summed tests comfortably. That slice - cannot be turned into a datacard at all -- a negative DY integral is rejected - outright by resolveNegativeBins -- so the binning has to not select it in the - first place. Worse, the significance being maximised is S/sqrt(B + sigma_B^2), - which a downward fluctuation *increases* by lowering B, so such a window is - mildly preferred rather than merely tolerated. - - `exempt` comes from minor_backgrounds() judged over the whole sliced-axis range of the - category, never over the candidate window: a background that is negligible in - this category must not be able to veto every boundary, and judging it inside - the window under test is circular (a background that is exactly zero there is - trivially below any fraction of the total). - """ - total = sum(yields.values()) - if total <= min_sum: - return False - if min_neff > 0 and total_error is not None: - if effective_entries(total, total_error) < min_neff: - return False - for name, value in yields.items(): - # As in _bin_passes: exemption covers the magnitude floor, never the sign. - if value < 0: - return False - if name in exempt: - continue - if min_each > 0 and value <= min_each: - return False - if min_proc_neff > 0 and errors is not None: - if effective_entries(value, errors.get(name, 0.0)) < min_proc_neff: - return False - return True - - -def minor_backgrounds(bkg_hists_by_name, lo, hi, min_frac): - """Backgrounds negligible over the *whole* [lo, hi] range, which may therefore - be exempted from the per-bin min_each floor. - - This must be judged once over the full slice, never inside the candidate bin - under test. Evaluating the fraction per bin is circular: a background that is - exactly zero in that bin is trivially below any fraction of the bin total, so - it is always exempted -- which is precisely the case the floor exists to - catch. That let a bin through with TT = 0 in a slice where TT is 94% of the - background (m300, eMu, res2b_dnn2). - """ - if min_frac <= 0: - return set() - yields = _bkg_yields(bkg_hists_by_name, lo, hi) - total = sum(yields.values()) - if total <= 0: - return set() - return {name for name, value in yields.items() if value < min_frac * total} - - -def grow_slice( - sig_hist, - bkg_hists_by_name, - right, - first_bin, - min_sum, - min_neff=0.0, - sig_mode="sb", - min_each=0.0, - min_proc_neff=0.0, - exempt=(), -): - """Among every candidate [left, right] with summed background > min_sum, - pick the one maximizing S/sqrt(B + sigma_B^2) -- not just the first one that - clears it. This is what actually drives where the cut lands; the content - floor is only a validity gate. Falls back to first_bin (best effort) if no - candidate clears the floor anywhere.""" - best_left = None - best_sig = -1.0 - for left in range(right, first_bin - 1, -1): - bkg_y = _bkg_yields(bkg_hists_by_name, left, right) - b_err = _total_bkg_error(bkg_hists_by_name, left, right) - bkg_e = ( - _bkg_errors(bkg_hists_by_name, left, right) if min_proc_neff > 0 else None - ) - if not _slice_passes( - bkg_y, - min_sum, - b_err, - min_neff, - bkg_e, - min_each, - min_proc_neff, - exempt, - ): - continue - s = _integral(sig_hist, left, right) - b = sum(bkg_y.values()) - sig = significance(s, b, b_err, sig_mode) - if sig > best_sig: - best_sig = sig - best_left = left - return best_left if best_left is not None else first_bin - - -def find_slices( - sig_hist, - bkg_hists_by_name, - n_slices, - first_bin, - last_bin, - min_sum, - min_neff=0.0, - sig_mode="sb", - min_each=0.0, - min_proc_neff=0.0, - min_frac=0.0, -): - """Split [first_bin, last_bin] into exactly `n_slices` ranges (fixed count, - required so every mass point shares the same category list), scanning right - to left. Each slice (except the final, leftover one) is placed to maximize - signal significance among boundaries with summed background > min_sum. - - Which backgrounds are minor enough to be exempt from the per-process floors is - decided once here, over the whole [first_bin, last_bin] range, and held fixed - for every candidate window -- see _slice_passes on why it cannot be re-judged - per window. - """ - exempt = ( - minor_backgrounds(bkg_hists_by_name, first_bin, last_bin, min_frac) - if (min_each > 0 or min_proc_neff > 0) - else set() - ) - slices = [] - right = last_bin - for slice_idx in range(n_slices): - if slice_idx == n_slices - 1 or right <= first_bin: - if right < first_bin: - # The axis ran out before n_slices boundaries were placed. Reported as - # None rather than invented: the previous code emitted (0,0), (-1,-1), - # (-2,-2) here to keep the count, and ROOT reads any range whose upper - # bound is below its lower one as the *whole axis including overflow* -- - # so those became copies of the full plane, and the same events entered - # the fit two or three times as independent categories. discover_binning - # rejects the category instead. - slices.append(None) - continue - slices.append((first_bin, right)) - right = first_bin - 1 - continue - left = grow_slice( - sig_hist, - bkg_hists_by_name, - right, - first_bin, - min_sum, - min_neff, - sig_mode, - min_each, - min_proc_neff, - exempt, - ) - slices.append((left, right)) - right = left - 1 - slices.reverse() - return slices - - -def extend_outer_edges(ranges, full_lo, full_hi): - """Widen the first/last range to swallow under/overflow (bin 0 / nbins+1), - so no events are silently dropped at the extremes of the axis. - - Only real ranges are widened. Widening a None (an exhausted slice) would invent one, - and widening past a None would put the underflow on a range that is not the outermost - -- so a list containing any None is left alone and rejected by the caller. - """ - ranges = list(ranges) - if any(r is None for r in ranges): - return ranges - ranges[0] = (full_lo, ranges[0][1]) - ranges[-1] = (ranges[-1][0], full_hi) - return ranges - - -def bin_budget(bkg_hists_by_name, lo, hi, max_bins_per_slice, bkg_per_bin): - """How many bins this slice can actually afford. - - A fixed max_bins_per_slice is applied blind to slice content: the high-mass boosted - slices hold ~1.5 total background events and were still being split into 10 - bins, i.e. ~0.15 events per bin. Capping at B_slice / bkg_per_bin ties the - binning to the statistics that are really there. Returns max_bins_per_slice - unchanged when bkg_per_bin <= 0 (feature off). - """ - if bkg_per_bin <= 0: - return max_bins_per_slice - total = sum(_bkg_yields(bkg_hists_by_name, lo, hi).values()) - if total <= 0: - return 1 - return max(1, min(max_bins_per_slice, int(total / bkg_per_bin))) - - -def bin_edges(y_axis, y_ranges): - """Physical edges of the discovered y ranges, for booking the output TH1. - - The rebinned shapes used to be booked as n_bins over [0, n_bins], which threw - the axis scale away and left every plot labelled by bin index. The ranges are - contiguous and ordered, so the edges are each range's low edge plus the last - range's upper edge. extend_outer_edges() pushes the outer ranges into the - underflow/overflow bins, which have no finite edge of their own -- those are - clamped to the axis limits. - """ - n = y_axis.GetNbins() - edges = [y_axis.GetBinLowEdge(max(lo, 1)) for lo, _ in y_ranges] - edges.append(y_axis.GetBinUpEdge(min(y_ranges[-1][1], n))) - return edges diff --git a/common/binning_dp.py b/common/binning_dp.py index ea0043a..ffa78d7 100644 --- a/common/binning_dp.py +++ b/common/binning_dp.py @@ -2,20 +2,20 @@ The greedy binner asks its questions one candidate at a time and re-integrates the histograms for each, which is affordable because it only ever looks at ~nx windows. The -HME-box search looks at ~ny^2 windows and runs an exact partition search over ~nx^2 -candidate DNN bins inside every one, so the integrals have to become table lookups or the -search is not worth running. +y_window search looks at ~ny^2 windows and runs an exact partition search over ~nx^2 +candidate bins of x inside every one, so the integrals have to become table lookups or +the search is not worth running. That is all this module does: it turns the TH2s into cumulative arrays and hands back the same yields, errors and figure of merit the greedy path would have computed. -It deliberately knows nothing about what a valid bin is. The gates live in -binning_core (and _bin_passes in the 2D binner), and are passed in as callables, so -there is exactly one definition of what a bin has to contain and this module cannot -drift from it. What it provides is the *arguments* those gates take -- the yields and -errors dicts -- built from the tables instead of from ROOT. +It deliberately knows nothing about what a valid bin is or how a bin is scored. The +gate (_bin_passes) and the figure of merit (significance) live in the 2D binner and are +passed in as callables, so there is exactly one definition of each and this module +cannot drift from it. What it provides is the *arguments* those gates take -- the +yields and errors dicts -- built from the tables instead of from ROOT. Under- and overflow are carried, and that is load-bearing rather than tidy. _integral() -in binning_core reads a TH2 with Integral(lo, hi, 0, -1), and ROOT reads biny2 < biny1 as +in the 2D binner reads a TH2 with Integral(lo, hi, 0, -1), and ROOT reads biny2 < biny1 as "the whole y axis including under- and overflow", so every slice-level quantity in the greedy path already includes the y under/overflow rows. A table that stopped at bin ny would disagree with the greedy path by a fraction of an event, which is easily enough to @@ -24,32 +24,25 @@ Differencing prefix sums is not exact, and one of the ways it is inexact has teeth. A window whose content is ~1e8 times smaller than the whole plane's -- routine for DY in -a narrow high-HME window -- loses about eight digits to cancellation, so a quantity that +a narrow window at high y -- loses about eight digits to cancellation, so a quantity that is exactly 0.0 in ROOT comes back as -2e-16 here. That is numerically nothing and physically nothing, but _bin_passes() tests `value < 0` to reject a background that has -gone negative, and it would reject 137 of the windows sampled below on the strength of a -rounding sign. So every yield is snapped to zero below a per-process epsilon scaled by -that process's own L1 norm, which is the size of the cancellation that can occur. The +gone negative, and on the real inputs it rejected 137 sampled windows on the strength +of a rounding sign. So every yield is snapped to zero below a per-process epsilon scaled +by that process's own L1 norm, which is the size of the cancellation that can occur. The epsilon lands around 1e-8 events for the largest process here, against a magnitude floor (min_bin_bkg_each) of 0.01 -- six orders of margin, so the snap cannot mask a real negative. -The check that all of this is true is parity_check(), not inspection. It is the only -reason to trust anything built on top of this file. +The check that all of this is true is check_window_mask(), not inspection: every +window the search settles on has a sample of its table entries re-derived with _bin_passes and +significance on the histograms themselves, and the run stops on any disagreement. """ import math import numpy as np -from StatInference.common.binning_core import ( - _bkg_errors, - _bkg_yields, - _integral, - _total_bkg_error, - significance, -) - def _hist_arrays(hist): """(values, variances) as (nx+2, ny+2) arrays indexed [binx, biny], bins 0..n+1. @@ -108,8 +101,9 @@ class Cells: because the binning is derived from the combination it will be applied to. """ - def __init__(self, nx, ny, sig, bkg, var): + def __init__(self, nx, ny, sig, bkg, var, significance): self.nx = nx + self._significance = significance self.ny = ny self.names = list(bkg) self._sig = _prefix(sig) @@ -130,10 +124,6 @@ def __init__(self, nx, ny, sig, bkg, var): self._eps_var = {name: _cancellation_eps(var[name]) for name in self.names} self._eps_var_tot = _cancellation_eps(var_tot) - # -- the whole y axis, under- and overflow included: what a slice-level query means - def full_y(self): - return 0, self.ny + 1 - def signal(self, x0, x1, y0, y1): return _snap(_rect(self._sig, x0, x1, y0, y1), self._eps_sig) @@ -143,15 +133,6 @@ def total_bkg(self, x0, x1, y0, y1): def total_bkg_error(self, x0, x1, y0, y1): return math.sqrt(max(_rect(self._var_tot, x0, x1, y0, y1), 0.0)) - def variances(self, x0, x1, y0, y1): - """{background: squared MC statistical error}. Exposed because this, not the - error, is the quantity the tables actually add -- so it is the one a parity - check can bound the rounding of.""" - return {n: _rect(self._var[n], x0, x1, y0, y1) for n in self.names} - - def total_bkg_variance(self, x0, x1, y0, y1): - return _rect(self._var_tot, x0, x1, y0, y1) - def yields(self, x0, x1, y0, y1): """{background: yield}, the dict the gates take.""" return { @@ -171,13 +152,13 @@ def score(self, x0, x1, y0, y1, mode): Squared because Asimov Z^2 is what adds across independent counting bins, so a partition's value is the sum of its cells' scores and the search can be a - dynamic program. significance() is binning_core's, so the cells are ranked by - exactly the quantity the greedy path maximises. + dynamic program. `significance` is the 2D binner's own, passed in, so the cells + are ranked by exactly the quantity the greedy path maximises. """ s = self.signal(x0, x1, y0, y1) b = self.total_bkg(x0, x1, y0, y1) b_err = self.total_bkg_error(x0, x1, y0, y1) - return significance(s, b, b_err, mode) ** 2 + return self._significance(s, b, b_err, mode) ** 2 def exempt(self, x0, x1, y0, y1, min_frac): """minor_backgrounds() over a rectangle, for the range the caller is about to @@ -192,7 +173,7 @@ def exempt(self, x0, x1, y0, y1, min_frac): return {n for n, v in y.items() if v < min_frac * total} -def build_cells(sig2d, bkg2d_by_name): +def build_cells(sig2d, bkg2d_by_name, significance): """Cells for one channel/category, from the same histograms discover_binning() reads. sig2d is the already-summed discovery signal; bkg2d_by_name is @@ -212,101 +193,7 @@ def build_cells(sig2d, bkg2d_by_name): variances += e2 bkg[name] = values var[name] = variances - return Cells(nx, ny, sig, bkg, var) - - -def parity_check(cells, sig2d, bkg2d_by_name, n_samples=200, seed=0): - """Assert the tables agree with the ROOT integrals they replace. - - Everything downstream is only as trustworthy as this. Over random x windows it - compares the quantities the gates and the figure of merit are built from -- - per-process yields and variances, the summed variance, and the signal -- against - binning_core's own _bkg_yields/_bkg_errors/_total_bkg_error/_integral, which is what - the greedy path calls. Slice level (full y, under- and overflow included), because - that is the query whose ROOT spelling is easy to get wrong. - - **The tolerance is absolute and scaled by the cancellation, not relative.** A - relative tolerance is the wrong test for a difference of prefix sums: a window - holding 1e-8 of the plane's content has already spent eight of its sixteen digits - before the comparison starts, so demanding a fixed relative agreement on it is - demanding that floating point not be floating point. What *can* be bounded is the - absolute rounding, and it is bounded by the same per-process epsilon Cells uses to - snap -- the L1 norm times a few hundred machine epsilons. Errors are compared as - variances for the same reason: the variance is what the tables add, and the square - root only inflates the relative error of an already-cancelled quantity. - - It also asserts the two agree on the *sign* of every yield, which is a separate - property from agreeing on its value and the one the binning actually turns on: - _bin_passes rejects a bin outright when a background is negative, so a yield that is - 0.0 in ROOT and -2e-16 here is not a rounding difference, it is a different binning. - - Raises on the first disagreement rather than returning a verdict: a table that is - subtly wrong must not be usable. - """ - rng = np.random.default_rng(seed) - y0, y1 = cells.full_y() - checked = 0 - for _ in range(n_samples): - lo = int(rng.integers(1, cells.nx + 1)) - hi = int(rng.integers(lo, cells.nx + 1)) - - got_y = cells.yields(lo, hi, y0, y1) - want_y = _bkg_yields(bkg2d_by_name, lo, hi) - got_v = cells.variances(lo, hi, y0, y1) - want_e = _bkg_errors(bkg2d_by_name, lo, hi) - for name in want_y: - _assert_close( - got_y[name], - want_y[name], - cells._eps_bkg[name], - f"yield[{name}] {lo}..{hi}", - ) - _assert_close( - got_v[name], - want_e[name] ** 2, - cells._eps_var[name], - f"variance[{name}] {lo}..{hi}", - ) - _assert_same_sign(got_y[name], want_y[name], f"yield[{name}] {lo}..{hi}") - - _assert_close( - cells.total_bkg_variance(lo, hi, y0, y1), - _total_bkg_error(bkg2d_by_name, lo, hi) ** 2, - cells._eps_var_tot, - f"total bkg variance {lo}..{hi}", - ) - _assert_close( - cells.signal(lo, hi, y0, y1), - _integral(sig2d, lo, hi), - cells._eps_sig, - f"signal {lo}..{hi}", - ) - checked += 1 - return checked - - -def _assert_same_sign(got, want, what): - """The gates branch on sign, so agreeing to 1e-9 is not the same as agreeing.""" - if (got < 0) != (want < 0): - raise AssertionError( - f"binning_dp parity failure on the sign of {what}: table says {got!r}, " - f"ROOT says {want!r}. _bin_passes rejects a bin whose background is " - "negative, so these two would not produce the same binning even though the " - "magnitudes agree. The zero-snapping in Cells is what should have prevented " - "this." - ) - - -def _assert_close(got, want, atol, what): - """Agreement to within the rounding the prefix sums can actually produce.""" - if abs(got - want) > atol: - raise AssertionError( - f"binning_dp parity failure on {what}: table says {got!r}, ROOT says " - f"{want!r}; they differ by {abs(got - want):.3e}, which exceeds the " - f"{atol:.3e} this quantity's cancellation can account for. The cell tables " - "do not reproduce the integrals the gates are defined on, so any binning " - "derived from them would be cut in the wrong places." - ) + return Cells(nx, ny, sig, bkg, var, significance) NEG = -np.inf @@ -366,65 +253,6 @@ def partition_dp(score, valid, max_parts): return values, partitions -def _self_test(trials=400, seed=0): - """Check partition_dp against exhaustive enumeration. - - It is a short dynamic program whose failure mode is silence: a subtly wrong - recurrence still returns a partition, and the binning it produces still looks like a - binning. On problems small enough to enumerate completely there is no reason to - settle for a plausible answer, so this compares against every partition there is. - - Pure numpy, no ROOT and no input files, so it runs anywhere: - python3 -m StatInference.common.binning_dp - The half of this module that does need ROOT is checked by parity_check() instead. - """ - import itertools - - def brute_partitions(n, k): - for cuts in itertools.combinations(range(1, n), k - 1): - bounds = (0,) + cuts + (n,) - parts = [(bounds[i], bounds[i + 1] - 1) for i in range(k)] - yield parts - - rng = np.random.default_rng(seed) - checked = 0 - for _ in range(trials): - n = int(rng.integers(2, 9)) - score = rng.normal(size=(n, n)) - valid = rng.random((n, n)) > rng.uniform(0.0, 0.6) - for a in range(n): - for b in range(n): - if b < a: - valid[a, b] = False - score[a, b] = 0.0 - values, partitions = partition_dp(score, valid, n) - for k in range(1, n + 1): - best = NEG - for part in brute_partitions(n, k): - if all(valid[a, b] for a, b in part): - total = sum(score[a, b] for a, b in part) - best = max(best, total) - got = values[k] - if best == NEG: - assert got == NEG, f"found a {k}-part partition where none exists" - else: - assert ( - abs(got - best) == 0 or abs(got - best) < 1e-9 - ), f"partition_dp n={n} k={k}: {got} != optimum {best}" - part = partitions[k] - assert len(part) == k - assert ( - part[0][0] == 0 and part[-1][1] == n - 1 - ), "does not cover the axis" - assert all( - part[i][1] + 1 == part[i + 1][0] for i in range(k - 1) - ), "parts are not contiguous" - assert all(valid[a, b] for a, b in part), "returned an invalid range" - checked += 1 - - return checked - - def binning_objective(cells, slices, mode): """(total Z^2, per-slice Z^2) for a finished binning. @@ -440,12 +268,12 @@ def binning_objective(cells, slices, mode): """ per_slice = [] for sl in slices: - if "y_range" in sl: # HME box: selection on y, bins along x + if "y_range" in sl: # y_window: selection on y, bins along x ylo, yhi = sl["y_range"] per_slice.append( sum(cells.score(a, b, ylo, yhi, mode) for a, b in sl["x_ranges"]) ) - else: # DNN slice: selection on x, bins along y + else: # greedy slice: selection on x, bins along y xlo, xhi = sl["x_range"] per_slice.append( sum(cells.score(xlo, xhi, a, b, mode) for a, b in sl["y_ranges"]) @@ -463,8 +291,9 @@ def _ranges_from_column(column, n, include_outer=False): With include_outer the first and last bins swallow the under- and overflow, which is what extend_outer_edges() does to the outermost ranges before they are written. A search that scores the ranges without them is optimising something slightly different - from what ends up in the datacard -- 0.11% of the background on the DNN axis of these - shapes, small but not nothing, and concentrated entirely in the two outermost slices. + from what ends up in the datacard -- 0.11% of the background on the DNN axis of the + HH->bbWW shapes, small but not nothing, and concentrated entirely in the two + outermost slices. """ lo = column[1 : n + 1].copy() hi = column[2 : n + 2].copy() @@ -553,29 +382,29 @@ def _axis_ranges(prefix, a, b, n, axis, include_outer=False): return _ranges_from_column(column, n, include_outer) -def box_tables(cells, y0, y1, exempt, knobs, mode): +def window_tables(cells, y0, y1, exempt, knobs, mode): """(score, valid) over every candidate bin of the *sliced-on-y* layout. - The window is a cut on y -- an HME box around the resonance -- and the bins run along - x, the DNN. The gates are _bin_passes() and the figure of merit significance(), as + The window is a cut on y around the signal, and the bins run along x. The gates are + _bin_passes() and the figure of merit significance(), as whole-array operations: the scalar path -- building a yields dict and an errors dict per candidate and calling _bin_passes -- costs about 80 microseconds a candidate, and - there are nx(nx+1)/2 of them per box, for every box the search tries. + there are nx(nx+1)/2 of them per window, for every window the search tries. This is the one place that restates the gate rather than calling it, which is a real risk: two spellings of the same rule are two rules as soon as one of them is edited. - It is contained by check_box_mask(), which re-evaluates the canonical _bin_passes on a + It is contained by check_window_mask(), which re-evaluates the canonical _bin_passes on a random sample of candidates and raises on the first disagreement -- so the fast path is checked against the slow one on real inputs, on every run, rather than at review time. The bins tile the whole x axis including its under/overflow, because everything inside - the box is kept and binned. The box itself is a selection: what falls outside it is + the window is kept and binned. The window itself is a selection: what falls outside it is discarded, not swept into an outer bin, which is the difference between a cut and a slice. `exempt` is the set of processes excused from the per-bin floors, judged once over - the whole box -- see minor_backgrounds(). Positivity is never excused -- see + the whole window -- see minor_backgrounds(). Positivity is never excused -- see _bin_passes. """ nx = cells.nx @@ -610,7 +439,7 @@ def box_tables(cells, y0, y1, exempt, knobs, mode): return np.where(valid, score, 0.0), valid -def check_box_mask( +def check_window_mask( cells, y0, y1, @@ -622,7 +451,7 @@ def check_box_mask( n_samples=100, seed=0, ): - """The cross-check of box_tables() against the canonical scalar gate. + """The cross-check of window_tables() against the canonical scalar gate. Cheap -- a hundred calls against the thousands the tables replace -- and it is what makes the fast path safe to trust: the tables are re-derived on a random sample with @@ -635,7 +464,7 @@ def check_box_mask( for _ in range(n_samples): a = int(rng.integers(1, nx + 1)) b = int(rng.integers(a, nx + 1)) - # the outermost bins carry the x under/overflow, because box_tables builds them + # the outermost bins carry the x under/overflow, because window_tables builds them # with include_outer -- the reference has to ask the same question lo = 0 if a == 1 else a hi = nx + 1 if b == nx else b @@ -643,19 +472,14 @@ def check_box_mask( got = bool(valid[a - 1, b - 1]) if got != want: raise AssertionError( - f"binning_dp box gate disagrees with _bin_passes on DNN bins {a}..{b} of " - f"the HME box {y0}..{y1}: fast path says {'valid' if got else 'invalid'}, " + f"binning_dp window gate disagrees with _bin_passes on x bins {a}..{b} " + f"of the y window {y0}..{y1}: fast path says {'valid' if got else 'invalid'}, " f"_bin_passes says {'valid' if want else 'invalid'}." ) if want: expected = cells.score(lo, hi, y0, y1, mode) if abs(expected - score[a - 1, b - 1]) > 1e-9 * max(abs(expected), 1.0): raise AssertionError( - f"binning_dp box score disagrees with significance() on DNN bins " - f"{a}..{b} of box {y0}..{y1}: {score[a - 1, b - 1]!r} vs {expected!r}." + f"binning_dp window score disagrees with significance() on x bins " + f"{a}..{b} of the y window {y0}..{y1}: {score[a - 1, b - 1]!r} vs {expected!r}." ) - - -if __name__ == "__main__": - n_part = _self_test() - print(f"partition_dp: {n_part} (n, k) cases vs exhaustive enumeration -- OK") diff --git a/law/CreateDatacardsTask.py b/law/CreateDatacardsTask.py index 386a319..b380639 100644 --- a/law/CreateDatacardsTask.py +++ b/law/CreateDatacardsTask.py @@ -147,10 +147,10 @@ def run(self): def _surviving_axis(base_dir): """Which axis of the 2D input the rebinned shapes are binned along: 0 = x, 1 = y. - Read off the binning record the preprocess step wrote beside the shapes. A DNN - slice selects on x and bins y; an HME box (a slice with a "y_range") selects on y - and bins x. With no record the shapes are the DNN-sliced kind, whose surviving - axis is y -- which is also what this assumed before boxes existed. + Read off the binning record the preprocess step wrote beside the shapes. A greedy + slice selects on x and bins y; a y_window category (a slice with a "y_range") + selects on y and bins x. With no record the shapes are the x-sliced kind, whose + surviving axis is y -- which is also what this assumed before windows existed. """ import json @@ -175,10 +175,10 @@ def plot_variable(self, variable, axis=1): For input a 2D->1D rebinning produced, that is one of the two variables of the 2D entry -- histograms.yaml records both as ``var_list: [x, y]``, so the axis metadata the plotter needs is already described and does not have to be restated - here. `axis` picks which, from _surviving_axis(): the DNN-sliced shapes keep y - (HME), an HME box keeps x (the DNN). Taking y unconditionally labelled every box - panel "Deep HME Mass" under DNN-score bins. For input that was always 1D, the - variable is its own answer. + here. `axis` picks which, from _surviving_axis(): x-sliced shapes keep y, a + y_window keeps x. Taking y unconditionally labelled every window panel with the y + variable under bins of x (in HH->bbWW, "Deep HME Mass" under DNN-score bins). For + input that was always 1D, the variable is its own answer. """ try: import FLAF.Common.Setup as Setup From 99c4383fccb3ab8b41066442625745d8417c3b8b Mon Sep 17 00:00:00 2001 From: Devin Date: Sat, 26 Sep 2026 03:39:53 +0200 Subject: [PATCH 07/12] Name the strategy window and let the configuration choose its axis A window on y with x binned inside it is one orientation of a general method, so the strategy is now "window", and window_axis (y, the default, or x) says which axis the window is cut on; the other is the one binned inside it. --window-axis sets it from the command line, and any other value is refused. The search is written with the window on y. For window_axis x, build_cells builds the cells transposed, the same search runs, and the result is returned in the input's axes: a selection on x with bins along y, the layout the greedy slices already use, so the writer, binning.json, replay and plot labels need no change. Checked on Run3_Early, v2605a_DNNOutputs_v3, against Binning_2D_DP (the reference the method was measured on): window with window_axis y reproduces its rebinned shapes exactly at MX 500 and 1000 (188928 histograms), and greedy does too (755712). For window_axis x, every 2D input at MX 500 was transposed and rebinned with window_axis x: the binning is the same as window_axis y on the original input, and the 94464 shapes agree, 66464 bit-identical and the rest within 1e-9 relative (the transposed sums round differently). Co-Authored-By: Claude Opus 5.5 (1M context) --- bin_opt_2d/rebin_2d.py | 86 ++++++++++++++++++++++++++------------ common/binning_dp.py | 30 ++++++++----- law/CreateDatacardsTask.py | 6 +-- 3 files changed, 82 insertions(+), 40 deletions(-) diff --git a/bin_opt_2d/rebin_2d.py b/bin_opt_2d/rebin_2d.py index a71abd0..fad3dda 100644 --- a/bin_opt_2d/rebin_2d.py +++ b/bin_opt_2d/rebin_2d.py @@ -62,21 +62,24 @@ "significance_mode": "asimov", # How the edges are found. "greedy" is the original: x is cut into slices, each # boundary the best one taken on its own, and the y bins inside a slice are - # equal-signal quantiles merged until the gates pass. "y_window" does the reverse: one - # window on y per category, with x binned inside it by an exact partition search -- - # see _discover_y_window(). Default stays "greedy" so adding this moves nothing until - # a configuration asks for it. + # equal-signal quantiles merged until the gates pass. "window" keeps one window per + # category on window_axis, with the other axis binned inside it by an exact partition + # search -- see _discover_window(). Default stays "greedy" so adding this moves nothing + # until a configuration asks for it. "strategy": "greedy", - # y_window only: a bin must earn its place. The split that creates it has to raise + # window only: the axis the window is cut on, "y" or "x". The other axis is the one + # binned inside it. + "window_axis": "y", + # window only: a bin must earn its place. The split that creates it has to raise # this category's Z^2 by at least this fraction of the category's own achievable # total, or it is given back. 0.0 leaves the count to max_bins_per_slice alone. "dp_min_bin_gain": 0.0, - # y_window only: the search grid for the window edges. A full scan is ny(ny+1)/2 windows + # window only: the search grid for the window edges. A full scan is ny(ny+1)/2 windows # and every one of them carries its own dynamic program over the binned axis, so the # edges are scanned on a stride and then refined within +-window_refine of the winner. "window_stride": 6, "window_refine": 5, - # y_window only: the window is searched only where the signal is. A window that + # window only: the window is searched only where the signal is. A window that # starts beyond the signal adds background and no signal; scanning it is wasted time, # and including it in the grid costs resolution where the answer actually lies. "window_signal_quantile": 0.001, @@ -86,7 +89,8 @@ # choices= uses -- significance_mode already learned this lesson: it is read by a # function that silently treats anything unrecognised as a default, so a typo in the # yaml quietly changed the answer instead of failing. -BINNING_STRATEGIES = ("greedy", "y_window") +BINNING_STRATEGIES = ("greedy", "window") +WINDOW_AXES = ("y", "x") def load_binning_config(path, overrides=None): @@ -123,6 +127,11 @@ def load_binning_config(path, overrides=None): f"{path or 'binning configuration'}: strategy '{knobs['strategy']}' is not " f"one of {sorted(BINNING_STRATEGIES)}." ) + if knobs["window_axis"] not in WINDOW_AXES: + raise RuntimeError( + f"{path or 'binning configuration'}: window_axis '{knobs['window_axis']}' is " + f"not one of {sorted(WINDOW_AXES)}." + ) return knobs @@ -839,12 +848,15 @@ def _best_count(values, cap): def _window_value(cells, y0, y1, knobs, check=False): - """What one window on y is worth: the best 1D binning of x inside it. + """What one window is worth: the best 1D binning of the other axis inside it. + + In the cells' own frame the window is always on y and the bins along x; + _discover_window builds the cells transposed when window_axis is x. A window is scored by the binning it admits, not by its own integrated significance. Those are different questions and they pick different windows -- scored on its own, a - window is rewarded for swallowing as much signal as it can and the shape along x - inside it plays no part, which is the opposite of what the fit will do with it. + window is rewarded for swallowing as much signal as it can and the shape inside it + plays no part, which is the opposite of what the fit will do with it. """ mode = knobs["significance_mode"] # which backgrounds count as negligible is decided once, over the whole window -- see @@ -885,20 +897,26 @@ def bin_passes(a, b, ylo, yhi, exempt_set): return values[n], [(a + 1, b + 1) for a, b in partitions[n]] -def _discover_y_window(sig2d, bkg2d_by_name, knobs): - """One window on y per category, with x binned inside it. +def _discover_window(sig2d, bkg2d_by_name, knobs): + """One window per category on window_axis, with the other axis binned inside it. The greedy strategy cuts x into slices that each become a datacard category and gives - every slice a shape along y. This one does the reverse and keeps only one category: a - single contiguous window on y around the signal, and a 1D distribution of x inside it. - Everything outside the window is dropped -- that is what makes it a cut rather than a - slice, and it is why the window has to earn its acceptance loss in purity. It suits a - y the signal peaks in (in HH->bbWW, x is the DNN score and y the HME mass), and it - keeps only one region of y, so a signal with two separated peaks would lose one. + every slice a shape along y. This one keeps only one category: a single contiguous + window around the signal on one axis, and a 1D distribution of the other axis inside + it. Everything outside the window is dropped -- that is what makes it a cut rather than + a slice, and it is why the window has to earn its acceptance loss in purity. It suits + a window axis the signal peaks in (in HH->bbWW, window_axis y is the HME mass and x the + DNN score), and it keeps only one region of that axis, so a signal with two separated + peaks would lose one. + + The search below is written with the window on y. For window_axis x the cells are + built transposed, so the same search runs, and the result is turned back into the + original axes before it is returned. Window and binning are chosen together; see _window_value. """ - cells = build_cells(sig2d, bkg2d_by_name, significance) + transpose = knobs["window_axis"] == "x" + cells = build_cells(sig2d, bkg2d_by_name, significance, transpose=transpose) ny = cells.ny # where the signal actually is -- see window_signal_quantile cumulative = np.array( @@ -937,6 +955,10 @@ def _discover_y_window(sig2d, bkg2d_by_name, knobs): # yield disagree with the cut the category claims to be. y0 = 0 if y0 <= 1 else y0 y1 = ny + 1 if y1 >= ny else y1 + if transpose: + # back to the input's axes: the window is on x and the bins run along y, which is + # the layout the greedy slices use, so everything downstream already reads it + return [{"x_range": (y0, y1), "y_ranges": bins}] return [{"y_range": (y0, y1), "x_ranges": bins}] @@ -950,8 +972,8 @@ def discover_binning(sig2d, bkg2d_by_name, knobs): Returns None when the category cannot be binned, which the caller reports as a skip. """ - if knobs["strategy"] == "y_window": - return _discover_y_window(sig2d, bkg2d_by_name, knobs) + if knobs["strategy"] == "window": + return _discover_window(sig2d, bkg2d_by_name, knobs) return _discover_greedy(sig2d, bkg2d_by_name, knobs) @@ -1122,9 +1144,10 @@ def record_to_slices(record, x_axis, y_axis, where): def slice_parts(sl): """(selection range, bin ranges, which axis the selection is on) for one category. - Two layouts reach this point. The greedy strategy cuts x into categories and gives - each a shape along y: {"x_range", "y_ranges"}. The y_window strategy cuts y and gives - the one category a shape along x: {"y_range", "x_ranges"}. Everything downstream -- + Two layouts reach this point: a selection on x with a shape along y, + {"x_range", "y_ranges"} -- the greedy slices, and a window with window_axis x -- and a + selection on y with a shape along x, {"y_range", "x_ranges"} -- a window with + window_axis y. Everything downstream -- the writer, the record, the objective -- needs the same three things from either, and asking here is what keeps that code from having to know which strategy ran. """ @@ -1604,8 +1627,16 @@ def run( default=None, choices=list(BINNING_STRATEGIES), help="how the edges are found: 'greedy' = each slice boundary taken on its own " - "and equal-signal quantiles inside it, 'y_window' = one window on y per category " - "with x binned inside it by an exact partition search", + "and equal-signal quantiles inside it, 'window' = one window per category on " + "--window-axis with the other axis binned inside it by an exact partition search", + ) + parser.add_argument( + "--window-axis", + required=False, + type=str, + default=None, + choices=list(WINDOW_AXES), + help="for strategy 'window': the axis the window is cut on (default y)", ) parser.add_argument( "--significance-mode", @@ -1621,6 +1652,7 @@ def run( overrides = {name: getattr(args, name) for name in knob_args} overrides["significance_mode"] = args.significance_mode overrides["strategy"] = args.strategy + overrides["window_axis"] = args.window_axis knobs = load_binning_config(args.binning_config, overrides) frozen_binning = None diff --git a/common/binning_dp.py b/common/binning_dp.py index ffa78d7..f2dd4a7 100644 --- a/common/binning_dp.py +++ b/common/binning_dp.py @@ -2,8 +2,9 @@ The greedy binner asks its questions one candidate at a time and re-integrates the histograms for each, which is affordable because it only ever looks at ~nx windows. The -y_window search looks at ~ny^2 windows and runs an exact partition search over ~nx^2 -candidate bins of x inside every one, so the integrals have to become table lookups or +window search looks at ~ny^2 windows and runs an exact partition search over ~nx^2 +candidate bins of x inside every one (in the cells' own frame, where the window is on y; +build_cells transposes the input when the window is on x), so the integrals have to become table lookups or the search is not worth running. That is all this module does: it turns the TH2s into cumulative arrays and hands back the same yields, errors and figure of merit the greedy path would have computed. @@ -173,11 +174,14 @@ def exempt(self, x0, x1, y0, y1, min_frac): return {n for n, v in y.items() if v < min_frac * total} -def build_cells(sig2d, bkg2d_by_name, significance): +def build_cells(sig2d, bkg2d_by_name, significance, transpose=False): """Cells for one channel/category, from the same histograms discover_binning() reads. sig2d is the already-summed discovery signal; bkg2d_by_name is {background: [hist per discovery era]}, summed here the way _bkg_yields() sums it. + + transpose swaps the two axes, so the cells' x is the input's y and vice versa. The + window search is written with the window on y; this is how it runs on x instead. """ nx = sig2d.GetNbinsX() ny = sig2d.GetNbinsY() @@ -193,6 +197,11 @@ def build_cells(sig2d, bkg2d_by_name, significance): variances += e2 bkg[name] = values var[name] = variances + if transpose: + nx, ny = ny, nx + sig = sig.T.copy() + bkg = {name: a.T.copy() for name, a in bkg.items()} + var = {name: a.T.copy() for name, a in var.items()} return Cells(nx, ny, sig, bkg, var, significance) @@ -268,12 +277,12 @@ def binning_objective(cells, slices, mode): """ per_slice = [] for sl in slices: - if "y_range" in sl: # y_window: selection on y, bins along x + if "y_range" in sl: # selection on y, bins along x ylo, yhi = sl["y_range"] per_slice.append( sum(cells.score(a, b, ylo, yhi, mode) for a, b in sl["x_ranges"]) ) - else: # greedy slice: selection on x, bins along y + else: # selection on x, bins along y xlo, xhi = sl["x_range"] per_slice.append( sum(cells.score(xlo, xhi, a, b, mode) for a, b in sl["y_ranges"]) @@ -385,7 +394,8 @@ def _axis_ranges(prefix, a, b, n, axis, include_outer=False): def window_tables(cells, y0, y1, exempt, knobs, mode): """(score, valid) over every candidate bin of the *sliced-on-y* layout. - The window is a cut on y around the signal, and the bins run along x. The gates are + The window is a cut on the cells' y around the signal, and the bins run along their + x (build_cells transposes the input when the window is on x). The gates are _bin_passes() and the figure of merit significance(), as whole-array operations: the scalar path -- building a yields dict and an errors dict per candidate and calling _bin_passes -- costs about 80 microseconds a candidate, and @@ -472,14 +482,14 @@ def check_window_mask( got = bool(valid[a - 1, b - 1]) if got != want: raise AssertionError( - f"binning_dp window gate disagrees with _bin_passes on x bins {a}..{b} " - f"of the y window {y0}..{y1}: fast path says {'valid' if got else 'invalid'}, " + f"binning_dp window gate disagrees with _bin_passes on bins {a}..{b} " + f"of the window {y0}..{y1}: fast path says {'valid' if got else 'invalid'}, " f"_bin_passes says {'valid' if want else 'invalid'}." ) if want: expected = cells.score(lo, hi, y0, y1, mode) if abs(expected - score[a - 1, b - 1]) > 1e-9 * max(abs(expected), 1.0): raise AssertionError( - f"binning_dp window score disagrees with significance() on x bins " - f"{a}..{b} of the y window {y0}..{y1}: {score[a - 1, b - 1]!r} vs {expected!r}." + f"binning_dp window score disagrees with significance() on bins " + f"{a}..{b} of the window {y0}..{y1}: {score[a - 1, b - 1]!r} vs {expected!r}." ) diff --git a/law/CreateDatacardsTask.py b/law/CreateDatacardsTask.py index b380639..c000bc5 100644 --- a/law/CreateDatacardsTask.py +++ b/law/CreateDatacardsTask.py @@ -148,8 +148,8 @@ def _surviving_axis(base_dir): """Which axis of the 2D input the rebinned shapes are binned along: 0 = x, 1 = y. Read off the binning record the preprocess step wrote beside the shapes. A greedy - slice selects on x and bins y; a y_window category (a slice with a "y_range") - selects on y and bins x. With no record the shapes are the x-sliced kind, whose + slice, or a window on x, selects on x and bins y; a window on y (a slice with a + "y_range") selects on y and bins x. With no record the shapes are the x-sliced kind, whose surviving axis is y -- which is also what this assumed before windows existed. """ import json @@ -176,7 +176,7 @@ def plot_variable(self, variable, axis=1): entry -- histograms.yaml records both as ``var_list: [x, y]``, so the axis metadata the plotter needs is already described and does not have to be restated here. `axis` picks which, from _surviving_axis(): x-sliced shapes keep y, a - y_window keeps x. Taking y unconditionally labelled every window panel with the y + window on y keeps x. Taking y unconditionally labelled every window panel with the y variable under bins of x (in HH->bbWW, "Deep HME Mass" under DNN-score bins). For input that was always 1D, the variable is its own answer. """ From de355041daef83a934ee9e8898ce94221a35b23e Mon Sep 17 00:00:00 2001 From: Devin Date: Sat, 26 Sep 2026 04:07:09 +0200 Subject: [PATCH 08/12] Shorten the comments and docstrings of the window strategy Comments and docstrings only: the code is unchanged (the syntax tree without docstrings is identical to the previous commit's). Co-Authored-By: Claude Opus 5.5 (1M context) --- bin_opt_2d/rebin_2d.py | 139 ++++++---------------- common/binning_dp.py | 230 ++++++++++--------------------------- law/CreateDatacardsTask.py | 15 +-- 3 files changed, 98 insertions(+), 286 deletions(-) diff --git a/bin_opt_2d/rebin_2d.py b/bin_opt_2d/rebin_2d.py index fad3dda..abee8fe 100644 --- a/bin_opt_2d/rebin_2d.py +++ b/bin_opt_2d/rebin_2d.py @@ -60,35 +60,23 @@ "min_bkg_frac": 0.05, "min_signal": 0.5, "significance_mode": "asimov", - # How the edges are found. "greedy" is the original: x is cut into slices, each - # boundary the best one taken on its own, and the y bins inside a slice are - # equal-signal quantiles merged until the gates pass. "window" keeps one window per - # category on window_axis, with the other axis binned inside it by an exact partition - # search -- see _discover_window(). Default stays "greedy" so adding this moves nothing - # until a configuration asks for it. + # "greedy": slice x, equal-signal bins along y. "window": one window per category on + # window_axis, the other axis binned inside it by an exact search (_discover_window). "strategy": "greedy", - # window only: the axis the window is cut on, "y" or "x". The other axis is the one - # binned inside it. + # window only: the axis the window is cut on, "y" or "x" "window_axis": "y", - # window only: a bin must earn its place. The split that creates it has to raise - # this category's Z^2 by at least this fraction of the category's own achievable - # total, or it is given back. 0.0 leaves the count to max_bins_per_slice alone. + # window only: a bin is kept only if it raises the category's Z^2 by at least this + # fraction of its total (0 = keep up to max_bins_per_slice) "dp_min_bin_gain": 0.0, - # window only: the search grid for the window edges. A full scan is ny(ny+1)/2 windows - # and every one of them carries its own dynamic program over the binned axis, so the - # edges are scanned on a stride and then refined within +-window_refine of the winner. + # window only: window edges are scanned on this stride, then refined within + # +-window_refine of the best "window_stride": 6, "window_refine": 5, - # window only: the window is searched only where the signal is. A window that - # starts beyond the signal adds background and no signal; scanning it is wasted time, - # and including it in the grid costs resolution where the answer actually lies. + # window only: edges are searched between these signal quantiles "window_signal_quantile": 0.001, } -# Named here so the yaml path can be checked against the same list the command line's -# choices= uses -- significance_mode already learned this lesson: it is read by a -# function that silently treats anything unrecognised as a default, so a typo in the -# yaml quietly changed the answer instead of failing. +# checked for the yaml as well as by the command line's choices= BINNING_STRATEGIES = ("greedy", "window") WINDOW_AXES = ("y", "x") @@ -387,7 +375,7 @@ def _slice_passes( if effective_entries(total, total_error) < min_neff: return False for name, value in yields.items(): - # As in _bin_passes: exemption covers the magnitude floor, never the sign. + # exemption covers the magnitude floor, never the sign (see _bin_passes) if value < 0: return False if name in exempt: @@ -452,26 +440,15 @@ def _bin_passes( if effective_entries(total, total_error) < min_neff: return False for name, value in yields.items(): - # Positivity is never exemptable. `exempt` exists so a negligible process - # cannot veto every candidate split -- it is about the *magnitude* floor, not - # about sign. Letting an exempt process go negative does not merely leave an - # ugly bin: resolveNegativeBins in dc_make/maker.py rejects a negative bin - # that holds >= signalFractionForRelevantBins of the signal *unconditionally*, - # ignoring allow_negative_bins_within_error, so the datacard fails to build at - # all. Observed at m1000 muMu/SR/recovery_dnn1, where ST was 4.24% of the - # slice -- just under min_bkg_frac 0.05, hence exempt -- and came out at - # -0.567 +- 0.789, killing the production. + # Positivity is never exempt: an exempt process that goes negative in a + # signal-relevant bin makes resolveNegativeBins reject the datacard. if value < 0: return False if name in exempt: continue if value <= min_each: return False - # Every background must be measured, not merely present: a yield known only - # to a few hundred percent is not a background estimate. _slice_passes has - # carried this arm for the sliced axis; the mass axis, where essentially all - # fit bins live, had no per-process test at all -- only the summed one, which - # any single well-measured process satisfies on its own. + # each background must be measured, not only the sum if min_proc_neff > 0 and errors is not None: if effective_entries(value, errors.get(name, 0.0)) < min_proc_neff: return False @@ -823,22 +800,10 @@ def _discover_greedy(sig2d, bkg2d_by_name, knobs): def _best_count(values, cap): - """The bin count at or below `cap` that actually scores best, or None if none is - feasible. - - Not the largest feasible count, which is the obvious choice and the wrong one. With - the background uncertainty folded into the figure of merit -- which is what - significance() does, and the reason it does not reward a background that fluctuated - low -- splitting a bin is no longer guaranteed to raise Z^2: each half carries a - larger relative sigma_B, and for a slice that is already statistics-limited the - split can cost more than the extra shape information gains. Taking the largest - feasible count would then hand back a binning worse than a coarser one that passes - the same gates. - - Choosing the argmax instead means every coarser feasible count is one of the - candidates, so the winner is at least as good as any of them. It also tends to spend - fewer bins than the budget allows, which is the right direction -- an extra bin that - buys nothing is pure MC-statistical exposure. + """The bin count at or below `cap` that scores best, or None if none is feasible. + + Not simply the largest: with the background error in the figure of merit, an extra + split can lower Z^2. """ best, best_n = -np.inf, None for n in range(1, min(cap, len(values) - 1) + 1): @@ -848,19 +813,13 @@ def _best_count(values, cap): def _window_value(cells, y0, y1, knobs, check=False): - """What one window is worth: the best 1D binning of the other axis inside it. - - In the cells' own frame the window is always on y and the bins along x; - _discover_window builds the cells transposed when window_axis is x. + """(Z^2, bins) of the best binning inside window y0..y1 (in the cells' frame). - A window is scored by the binning it admits, not by its own integrated significance. - Those are different questions and they pick different windows -- scored on its own, a - window is rewarded for swallowing as much signal as it can and the shape inside it - plays no part, which is the opposite of what the fit will do with it. + A window is scored by the binning it admits, not by its integrated significance, + which would favour windows that swallow signal and ignore the shape. """ mode = knobs["significance_mode"] - # which backgrounds count as negligible is decided once, over the whole window -- see - # minor_backgrounds() on why it cannot be re-judged inside each candidate bin + # negligible backgrounds are judged once, over the whole window exempt = cells.exempt(1, cells.nx, y0, y1, knobs["min_bkg_frac"]) score, valid = window_tables(cells, y0, y1, exempt, knobs, mode) if check: @@ -900,25 +859,16 @@ def bin_passes(a, b, ylo, yhi, exempt_set): def _discover_window(sig2d, bkg2d_by_name, knobs): """One window per category on window_axis, with the other axis binned inside it. - The greedy strategy cuts x into slices that each become a datacard category and gives - every slice a shape along y. This one keeps only one category: a single contiguous - window around the signal on one axis, and a 1D distribution of the other axis inside - it. Everything outside the window is dropped -- that is what makes it a cut rather than - a slice, and it is why the window has to earn its acceptance loss in purity. It suits - a window axis the signal peaks in (in HH->bbWW, window_axis y is the HME mass and x the - DNN score), and it keeps only one region of that axis, so a signal with two separated - peaks would lose one. - - The search below is written with the window on y. For window_axis x the cells are - built transposed, so the same search runs, and the result is turned back into the - original axes before it is returned. + Events outside the window are dropped. Suited to a window axis the signal peaks in + (HH->bbWW: HME on y, DNN on x); only one region of it is kept. - Window and binning are chosen together; see _window_value. + The search is written with the window on y; for window_axis x the cells are built + transposed and the result is mapped back to the input's axes. """ transpose = knobs["window_axis"] == "x" cells = build_cells(sig2d, bkg2d_by_name, significance, transpose=transpose) ny = cells.ny - # where the signal actually is -- see window_signal_quantile + # only search where the signal is cumulative = np.array( [cells.signal(0, cells.nx + 1, 1, j) for j in range(1, ny + 1)] ) @@ -950,14 +900,11 @@ def _discover_window(sig2d, bkg2d_by_name, knobs): value, bins = _window_value(cells, y0, y1, knobs, check=True) if bins is None: return None - # A window touching an end of the axis takes that end's overflow with it: the events - # are inside the selection the window stands for, and dropping them would make the written - # yield disagree with the cut the category claims to be. + # a window at an end of the axis includes that end's under/overflow y0 = 0 if y0 <= 1 else y0 y1 = ny + 1 if y1 >= ny else y1 if transpose: - # back to the input's axes: the window is on x and the bins run along y, which is - # the layout the greedy slices use, so everything downstream already reads it + # window on x, bins along y: the same layout as the greedy slices return [{"x_range": (y0, y1), "y_ranges": bins}] return [{"y_range": (y0, y1), "x_ranges": bins}] @@ -1052,11 +999,7 @@ def slice_ranges(x_axis, y_axis, slices): def _slice_record(sl, sel_edges, x_axis, y_axis): - """One category as plain data, in whichever layout produced it. - - The key names say which axis is the selection and which carries the bins, so a reader - -- human or replay -- never has to infer it from the strategy that happened to run. - """ + """One category as plain data; the keys say which axis is selected and which binned.""" (lo, hi), bins, slice_axis = slice_parts(sl) binned_axis = y_axis if slice_axis == "x" else x_axis sel_key, bin_key = ( @@ -1082,11 +1025,8 @@ def slices_to_record(slices, x_axis, y_axis, objective=None): are what a replay restores; the physical edges alongside them are what a human reads, and what makes the record meaningful next to a plot. - `objective` is what the binning scored, recomputed from these ranges rather than - carried out of the optimiser -- so it is filled in for the greedy strategy too, which - never computes it, and two runs can be compared from their records alone. It is - written as a sibling of the slices and never read back: record_to_slices() takes only - the ranges, so a record written before this existed still replays. + `objective` (the binning's Z^2) is written for information only; replay reads only + the ranges. """ record = { "n_x_bins": x_axis.GetNbins(), @@ -1144,12 +1084,8 @@ def record_to_slices(record, x_axis, y_axis, where): def slice_parts(sl): """(selection range, bin ranges, which axis the selection is on) for one category. - Two layouts reach this point: a selection on x with a shape along y, - {"x_range", "y_ranges"} -- the greedy slices, and a window with window_axis x -- and a - selection on y with a shape along x, {"y_range", "x_ranges"} -- a window with - window_axis y. Everything downstream -- - the writer, the record, the objective -- needs the same three things from either, and - asking here is what keeps that code from having to know which strategy ran. + Layouts: {"x_range", "y_ranges"} (greedy slices, window on x) and + {"y_range", "x_ranges"} (window on y). """ if "y_range" in sl: return sl["y_range"], sl["x_ranges"], "y" @@ -1333,13 +1269,8 @@ def load2d(f, key): out_file.cd(f"{channel}/{naming.name(category, slice_idx)}") h.Write(key) - # Scored from the ranges that were just written, for both strategies and for a replay - # alike, so the number in the record always describes the shapes beside it. Cheap - # next to the search, and the greedy path has no other way to report what it achieved. - # Only when there is something to score. A replay is allowed to proceed with a - # background missing from the input -- it applies recorded edges and asks no - # questions -- and a figure of merit computed against no background would be a - # number, but not a true one. + # the Z^2 of the written binning, for the record; skipped without backgrounds (a + # replay can run on input missing one) objective = ( binning_objective( build_cells(disc_sig, disc_bkg_by_name, significance), diff --git a/common/binning_dp.py b/common/binning_dp.py index f2dd4a7..a48d4c1 100644 --- a/common/binning_dp.py +++ b/common/binning_dp.py @@ -1,43 +1,21 @@ -"""Prefix-sum cells over a 2D shape, so any candidate binning can be scored in O(1). - -The greedy binner asks its questions one candidate at a time and re-integrates the -histograms for each, which is affordable because it only ever looks at ~nx windows. The -window search looks at ~ny^2 windows and runs an exact partition search over ~nx^2 -candidate bins of x inside every one (in the cells' own frame, where the window is on y; -build_cells transposes the input when the window is on x), so the integrals have to become table lookups or -the search is not worth running. -That is all this module does: it turns the TH2s into cumulative arrays and hands back the -same yields, errors and figure of merit the greedy path would have computed. - -It deliberately knows nothing about what a valid bin is or how a bin is scored. The -gate (_bin_passes) and the figure of merit (significance) live in the 2D binner and are -passed in as callables, so there is exactly one definition of each and this module -cannot drift from it. What it provides is the *arguments* those gates take -- the -yields and errors dicts -- built from the tables instead of from ROOT. - -Under- and overflow are carried, and that is load-bearing rather than tidy. _integral() -in the 2D binner reads a TH2 with Integral(lo, hi, 0, -1), and ROOT reads biny2 < biny1 as -"the whole y axis including under- and overflow", so every slice-level quantity in the -greedy path already includes the y under/overflow rows. A table that stopped at bin ny -would disagree with the greedy path by a fraction of an event, which is easily enough to -move a boundary and produce a difference nobody can explain. So the arrays are indexed -0..n+1 throughout, and a slice-level query spans y = 0..ny+1. - -Differencing prefix sums is not exact, and one of the ways it is inexact has teeth. -A window whose content is ~1e8 times smaller than the whole plane's -- routine for DY in -a narrow window at high y -- loses about eight digits to cancellation, so a quantity that -is exactly 0.0 in ROOT comes back as -2e-16 here. That is numerically nothing and -physically nothing, but _bin_passes() tests `value < 0` to reject a background that has -gone negative, and on the real inputs it rejected 137 sampled windows on the strength -of a rounding sign. So every yield is snapped to zero below a per-process epsilon scaled -by that process's own L1 norm, which is the size of the cancellation that can occur. The -epsilon lands around 1e-8 events for the largest process here, against a magnitude floor -(min_bin_bkg_each) of 0.01 -- six orders of margin, so the snap cannot mask a real -negative. - -The check that all of this is true is check_window_mask(), not inspection: every -window the search settles on has a sample of its table entries re-derived with _bin_passes and -significance on the histograms themselves, and the run stops on any disagreement. +"""Prefix-sum tables over a 2D shape, so the window search can score any bin in O(1). + +The window search tries ~ny^2 windows, each with an exact partition search over ~nx^2 +candidate bins, which is only affordable with table lookups instead of TH2 integrals. +The search is written with the window on y; build_cells transposes the input for a +window on x. + +The gate (_bin_passes) and the figure of merit (significance) stay in rebin_2d.py and are +passed in, so each has one definition. + +Notes: +- Under/overflow are included (arrays are indexed 0..n+1), to match what the TH2 + integrals in rebin_2d.py return. +- Differencing prefix sums can turn an exact 0 into -2e-16, which the `value < 0` + positivity gate would reject. Yields are snapped to 0 below an epsilon scaled by each + process's L1 norm (~1e-8 events, far below the 0.01 floor). +- check_window_mask() re-checks a sample of every chosen window against _bin_passes and + significance on each run. """ import math @@ -46,13 +24,7 @@ def _hist_arrays(hist): - """(values, variances) as (nx+2, ny+2) arrays indexed [binx, biny], bins 0..n+1. - - Read through GetBinContent/GetBinError rather than the internal buffer: the buffer - is faster but its layout is a ROOT implementation detail, and getting it wrong is - silent. This runs once per histogram per category, which is nothing next to the - search it feeds. - """ + """(values, variances) as (nx+2, ny+2) arrays indexed [binx, biny], bins 0..n+1.""" nx = hist.GetNbinsX() ny = hist.GetNbinsY() values = np.empty((nx + 2, ny + 2), dtype=np.float64) @@ -75,11 +47,8 @@ def _prefix(a): def _cancellation_eps(a, rel=1e-12): """How large a value this array's prefix sums can invent out of rounding. - Scaled by the L1 norm rather than the sum, because the cancellation that matters - comes from negative-weight events: a process whose bins nearly cancel has a small - sum and a large L1, and it is the L1 that sets how much precision the difference of - two prefix sums can lose. rel is a few hundred times the double-precision epsilon, - which covers accumulating over the ~1e4 cells of one of these planes. + Scaled by the L1 norm, not the sum: with negative weights a process can nearly cancel + to a small sum while its prefix sums still lose precision on the L1 scale. """ return rel * float(np.abs(a).sum()) @@ -97,9 +66,8 @@ def _rect(p, x0, x1, y0, y1): class Cells: """Every quantity the binning gates and the figure of merit need, as O(1) lookups. - Background yields are summed across discovery eras and their MC errors added in - quadrature, which is what _bkg_yields()/_bkg_errors() do -- the eras are summed - because the binning is derived from the combination it will be applied to. + Backgrounds are summed across discovery eras (errors in quadrature), as + _bkg_yields()/_bkg_errors() do. """ def __init__(self, nx, ny, sig, bkg, var, significance): @@ -110,10 +78,7 @@ def __init__(self, nx, ny, sig, bkg, var, significance): self._sig = _prefix(sig) self._bkg = {name: _prefix(bkg[name]) for name in self.names} self._var = {name: _prefix(var[name]) for name in self.names} - # An empty background dict is reachable -- a replay whose input is missing a - # process it was derived with -- and sum() of nothing is the integer 0, which is - # not an array. Zeros keep every query defined and answering zero, which is the - # truthful answer when there is no background to integrate. + # no backgrounds at all is possible (a replay whose input lacks them); keep arrays zero = np.zeros_like(sig) bkg_tot = sum(bkg.values()) if bkg else zero var_tot = sum(var.values()) if var else zero @@ -149,22 +114,15 @@ def errors(self, x0, x1, y0, y1): } def score(self, x0, x1, y0, y1, mode): - """The figure of merit for one cell, squared. - - Squared because Asimov Z^2 is what adds across independent counting bins, so a - partition's value is the sum of its cells' scores and the search can be a - dynamic program. `significance` is the 2D binner's own, passed in, so the cells - are ranked by exactly the quantity the greedy path maximises. - """ + """significance()^2 for one rectangle; Z^2 adds across bins, so a partition's + value is the sum of its bins'.""" s = self.signal(x0, x1, y0, y1) b = self.total_bkg(x0, x1, y0, y1) b_err = self.total_bkg_error(x0, x1, y0, y1) return self._significance(s, b, b_err, mode) ** 2 def exempt(self, x0, x1, y0, y1, min_frac): - """minor_backgrounds() over a rectangle, for the range the caller is about to - subdivide -- never for the candidate sub-range itself, which is circular. - """ + """minor_backgrounds() over the rectangle about to be subdivided.""" if min_frac <= 0: return set() y = self.yields(x0, x1, y0, y1) @@ -175,13 +133,10 @@ def exempt(self, x0, x1, y0, y1, min_frac): def build_cells(sig2d, bkg2d_by_name, significance, transpose=False): - """Cells for one channel/category, from the same histograms discover_binning() reads. - - sig2d is the already-summed discovery signal; bkg2d_by_name is - {background: [hist per discovery era]}, summed here the way _bkg_yields() sums it. + """Cells for one channel/category. - transpose swaps the two axes, so the cells' x is the input's y and vice versa. The - window search is written with the window on y; this is how it runs on x instead. + sig2d is the summed discovery signal; bkg2d_by_name is {background: [hist per era]}. + transpose swaps x and y, which is how the window search runs with the window on x. """ nx = sig2d.GetNbinsX() ny = sig2d.GetNbinsY() @@ -211,24 +166,12 @@ def build_cells(sig2d, bkg2d_by_name, significance, transpose=False): def partition_dp(score, valid, max_parts): """Best partition of [0..n-1] into exactly k contiguous valid ranges, for every k. - Returns (values, partitions): values[k] is the total score of the best k-part - partition, or -inf if there is none, and partitions[k] is that partition as a list of - (lo, hi) index pairs. Index 0 is unused in both so that k reads as the part count. + Returns (values, partitions): values[k] is the best k-part total score (-inf if none) + and partitions[k] its (lo, hi) index pairs; index 0 is unused. Every k is returned + because the bin count is chosen from the whole curve. - Every k is returned rather than just the requested one, because the caller needs the - whole curve and it costs nothing to keep: the bin count is chosen by comparing what - the partition scores at each count, and trim_by_marginal_gain() gives a bin back by - comparing values[k] against values[k-1]. Computing them one at a time would re-run the - same DP. - - Exact, not greedy. The score is additive over the parts -- Asimov Z^2 is what adds - across independent counting bins -- so the optimal k-part partition of a prefix is - built from an optimal (k-1)-part partition of a shorter prefix, and the usual - interval DP applies. O(max_parts * n^2). - - `valid` is a mask, not a penalty: an invalid range is unreachable rather than merely - expensive, so a partition containing one cannot be returned at any score. That is what - keeps the background gates hard constraints instead of preferences. + Exact interval DP, O(max_parts * n^2), valid because the score is additive over parts. + Invalid ranges are excluded outright, so the background gates are hard constraints. """ n = score.shape[0] max_parts = max(1, min(max_parts, n)) @@ -263,17 +206,9 @@ def partition_dp(score, valid, max_parts): def binning_objective(cells, slices, mode): - """(total Z^2, per-slice Z^2) for a finished binning. + """(total Z^2, per-slice Z^2) of a finished binning, from the recorded ranges. - The quantity the search maximises, recomputed from the binning that was actually - written rather than carried out of the optimiser. That makes it meaningful for the - greedy strategy too -- which never computes it -- so the two can be compared from - their binning.json alone, and it means a discrepancy between what the optimiser - thought it achieved and what the shapes contain shows up as a discrepancy rather - than going unnoticed. - - Note this is evaluated on the ranges as recorded, i.e. after extend_outer_edges has - pushed the outermost ones into under/overflow, so it describes the shapes on disk. + Computed for every strategy, so binnings can be compared from binning.json alone. """ per_slice = [] for sl in slices: @@ -291,18 +226,10 @@ def binning_objective(cells, slices, mode): def _ranges_from_column(column, n, include_outer=False): - """M[i, j] = the rectangle sum for bins (i+1)..(j+1), from one column of a prefix sum. - - For a fixed x window, P[x1+1, :] - P[x0, :] is the running sum along y, and every - y range is a difference of two of its entries -- so the whole (n, n) table of - candidate ranges is one outer subtraction rather than n^2 lookups. - - With include_outer the first and last bins swallow the under- and overflow, which is - what extend_outer_edges() does to the outermost ranges before they are written. A - search that scores the ranges without them is optimising something slightly different - from what ends up in the datacard -- 0.11% of the background on the DNN axis of the - HH->bbWW shapes, small but not nothing, and concentrated entirely in the two - outermost slices. + """M[i, j] = the sum over bins (i+1)..(j+1), from one column of a prefix sum. + + The whole (n, n) table of candidate ranges is one outer subtraction. With + include_outer the first and last bins also take the under/overflow. """ lo = column[1 : n + 1].copy() hi = column[2 : n + 2].copy() @@ -313,14 +240,8 @@ def _ranges_from_column(column, n, include_outer=False): def _neff_matrix(value, error): - """effective_entries() over arrays, with its two degenerate cases kept. - - error <= 0 means the yield carries no MC uncertainty at all: infinite effective - entries if there is something there, none if there is not. Taking (v/e)^2 blindly - would make the first case a division by zero and the second a nan, and a nan - compares false against every threshold -- which would silently reject exactly the - empty bins the scalar gate accepts. - """ + """effective_entries() over arrays. error <= 0 gives inf if value > 0, else 0, as + the scalar does (a nan would silently fail every threshold).""" with np.errstate(divide="ignore", invalid="ignore"): ratio = np.where(error > 0, value / np.where(error > 0, error, 1.0), 0.0) ** 2 degenerate = np.where(value > 0, np.inf, 0.0) @@ -356,21 +277,11 @@ def _sb_z2_matrix(s, b, b_err): def trim_by_marginal_gain(values_per_slice, counts, threshold): - """Give back every bin whose last split bought less than `threshold`. - - With Z^2 very nearly monotone in bin count, a search told only to maximise it will - always spend whatever it is given. So the binding constraint has to be the value of a - bin, not the size of a pot. A bin is kept only if the split that created it raised - this category's Z^2 by at least `threshold` -- expressed as a fraction of the - category's own achievable total, so it means "worth something on the scale that - matters here" rather than a yield in events. A bin that buys nothing is - MC-statistical exposure bought for nothing; this is what declines it. - - Stepping down one bin at a time and stopping at the first marginal gain that clears - the threshold is exact for a concave curve, which Z^2 in bin count is: each extra - split has less left to separate. It is not assumed -- a curve that is not concave - simply stops at its first qualifying step, which is still a count whose last bin paid - for itself. + """Give back every bin whose split raised Z^2 by less than `threshold`. + + Z^2 almost always grows with the bin count, so without this the search would use + every bin it is allowed. The caller sets `threshold` as a fraction of the category's + total, and the count is stepped down until the last bin pays for itself. """ if threshold <= 0: return list(counts) @@ -392,30 +303,14 @@ def _axis_ranges(prefix, a, b, n, axis, include_outer=False): def window_tables(cells, y0, y1, exempt, knobs, mode): - """(score, valid) over every candidate bin of the *sliced-on-y* layout. - - The window is a cut on the cells' y around the signal, and the bins run along their - x (build_cells transposes the input when the window is on x). The gates are - _bin_passes() and the figure of merit significance(), as - whole-array operations: the scalar path -- building a yields dict and an errors dict - per candidate and calling _bin_passes -- costs about 80 microseconds a candidate, and - there are nx(nx+1)/2 of them per window, for every window the search tries. - - This is the one place that restates the gate rather than calling it, which is a real - risk: two spellings of the same rule are two rules as soon as one of them is edited. - It is contained by check_window_mask(), which re-evaluates the canonical _bin_passes on a - random sample of candidates and raises on the first disagreement -- so the fast path - is checked against the slow one on real inputs, on every run, rather than at review - time. - - The bins tile the whole x axis including its under/overflow, because everything inside - the window is kept and binned. The window itself is a selection: what falls outside it is - discarded, not swept into an outer bin, which is the difference between a cut and a - slice. - - `exempt` is the set of processes excused from the per-bin floors, judged once over - the whole window -- see minor_backgrounds(). Positivity is never excused -- see - _bin_passes. + """(score, valid) for every candidate x bin inside the window y0..y1. + + _bin_passes() and significance() restated as array operations: calling them per + candidate is far too slow for the search. check_window_mask() guards against the two + versions drifting apart. + + `exempt` (processes excused from the per-bin floors) is judged once over the whole + window. Positivity is never excused. """ nx = cells.nx ranges = lambda prefix: _axis_ranges(prefix, y0, y1, nx, "x", True) @@ -461,21 +356,14 @@ def check_window_mask( n_samples=100, seed=0, ): - """The cross-check of window_tables() against the canonical scalar gate. - - Cheap -- a hundred calls against the thousands the tables replace -- and it is what - makes the fast path safe to trust: the tables are re-derived on a random sample with - _bin_passes and significance(), so if the two spellings of the rule ever disagree - about a candidate bin the run stops here instead of quietly producing a different - binning. - """ + """Re-derive a random sample of window_tables() with _bin_passes and significance(), + and stop the run on any disagreement.""" nx = cells.nx rng = np.random.default_rng(seed) for _ in range(n_samples): a = int(rng.integers(1, nx + 1)) b = int(rng.integers(a, nx + 1)) - # the outermost bins carry the x under/overflow, because window_tables builds them - # with include_outer -- the reference has to ask the same question + # the outermost bins include the x under/overflow, as in window_tables lo = 0 if a == 1 else a hi = nx + 1 if b == nx else b want = bool(bin_passes(lo, hi, y0, y1, exempt)) diff --git a/law/CreateDatacardsTask.py b/law/CreateDatacardsTask.py index c000bc5..88794e1 100644 --- a/law/CreateDatacardsTask.py +++ b/law/CreateDatacardsTask.py @@ -147,10 +147,8 @@ def run(self): def _surviving_axis(base_dir): """Which axis of the 2D input the rebinned shapes are binned along: 0 = x, 1 = y. - Read off the binning record the preprocess step wrote beside the shapes. A greedy - slice, or a window on x, selects on x and bins y; a window on y (a slice with a - "y_range") selects on y and bins x. With no record the shapes are the x-sliced kind, whose - surviving axis is y -- which is also what this assumed before windows existed. + Read from binning.json: a slice with a "y_range" (a window on y) is binned along x, + anything else along y. Without a record, y. """ import json @@ -172,13 +170,8 @@ def _surviving_axis(base_dir): def plot_variable(self, variable, axis=1): """The variable whose axis the shapes are binned along, for histograms.yaml. - For input a 2D->1D rebinning produced, that is one of the two variables of the 2D - entry -- histograms.yaml records both as ``var_list: [x, y]``, so the axis - metadata the plotter needs is already described and does not have to be restated - here. `axis` picks which, from _surviving_axis(): x-sliced shapes keep y, a - window on y keeps x. Taking y unconditionally labelled every window panel with the y - variable under bins of x (in HH->bbWW, "Deep HME Mass" under DNN-score bins). For - input that was always 1D, the variable is its own answer. + For 2D->1D rebinned input, the entry of the 2D variable's ``var_list: [x, y]`` + picked by `axis` (from _surviving_axis()); for 1D input, the variable itself. """ try: import FLAF.Common.Setup as Setup From ec8d016aa957243b62d7801d3d7e7f8c8227711b Mon Sep 17 00:00:00 2001 From: Devin Date: Sat, 26 Sep 2026 04:46:51 +0200 Subject: [PATCH 09/12] Read rebinning inputs from per-model-point directories rebin_2d.py reads each base category from where input_categories places it for the file's mass (e.g. SR/res2b from SR_SL_M${MX}/res2b) and writes it under the category's own name, so binning.json and the rebinned shapes are keyed the same way for every mass. Co-Authored-By: Claude Opus 5.5 (1M context) --- bin_opt_2d/rebin_2d.py | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/bin_opt_2d/rebin_2d.py b/bin_opt_2d/rebin_2d.py index abee8fe..adcf6bd 100644 --- a/bin_opt_2d/rebin_2d.py +++ b/bin_opt_2d/rebin_2d.py @@ -15,7 +15,7 @@ sys.path.append(base_dir) __package__ = pkg_dir_name -from StatInference.common.tools import importROOT, CategoryNaming +from StatInference.common.tools import importROOT, CategoryNaming, InputCategories from StatInference.common.param_parse import extractParameters, applyParameters from StatInference.dc_make.model import Model from StatInference.common.binning_dp import ( @@ -189,6 +189,8 @@ def load_config(config_path): # it against the pattern the shapes are actually written with. None when the # configuration names no pattern, and then there is nothing to check. "category_pattern": cfg.get("category_pattern"), + # where each base category is read from, if the input keys it per model point + "input_categories": InputCategories.fromConfig(cfg), "signal_hist_name_patterns": signal_hist_names, "signal_param_name": extractParameters(signal_hist_names[0])[0], "mass_values": mass_values, @@ -1164,12 +1166,15 @@ def process_category( # The resonance parameter is named by the configuration, not by this script -- it is # "MX" for bbWW but the slicing knows nothing about which parameter it is scanning. param_name = cfg["signal_param_name"] - prefix = f"{channel}/{category}/" + # read from where the input keeps this category for this mass (see InputCategories); + # written under the category's own name + in_dir = f"{channel}/{cfg['input_categories'].path(category, {param_name: mass})}" + prefix = f"{in_dir}/" # The key list comes from the first source era; a systematic that only some eras carry # is filled in from their nominal when the slices are written below. - cat_dir = sources[0][1].Get(f"{channel}/{category}") + cat_dir = sources[0][1].Get(in_dir) if not cat_dir: - print(f" [skip] {channel}/{category}: not found in {sources[0][1].GetName()}") + print(f" [skip] {in_dir}: not found in {sources[0][1].GetName()}") return signal_keys = [ @@ -1257,9 +1262,9 @@ def load2d(f, key): f"{channel}/{naming.name(category, slice_idx)}", format_var_range(*ranges[slice_idx], var=cfg["slice_var"]), ) - cat_dir = in_file.Get(f"{channel}/{category}") + cat_dir = in_file.Get(in_dir) if not cat_dir: - print(f" [skip] {source_era} {channel}/{category}: not in the input") + print(f" [skip] {source_era} {in_dir}: not in the input") continue for key in [k.GetName() for k in cat_dir.GetListOfKeys()]: hist2d = get_hist(in_file, prefix + key) From e9ff36fa901cffd3769b1fa84334c7aa1b7a981f Mon Sep 17 00:00:00 2001 From: Devin Date: Sat, 26 Sep 2026 04:47:23 +0200 Subject: [PATCH 10/12] Accept an input_categories key that is the base of a configured category A card that reads rebinned shapes lists the written names (SR/res2b_hmebox0) while its input_categories map is keyed by the base category (SR/res2b) that rebin_2d.py reads. The maker now accepts a key that is a configured category or the base of one, and still refuses any other, which catches typos. Sliced names are read under their own names. Co-Authored-By: Claude Opus 5.5 (1M context) --- dc_make/maker.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/dc_make/maker.py b/dc_make/maker.py index 002c00b..929ce38 100644 --- a/dc_make/maker.py +++ b/dc_make/maker.py @@ -74,12 +74,17 @@ def __init__( # writes them, not self.categories: once a category list is expanded against a # rebinning record (base "SR/res2b" -> "SR/res2b_hmebox0"), the map's keys -- the # configured names -- no longer appear in it, and every key would read as unknown. + # A key may also be the base of a configured sliced name: that map is for + # rebin_2d.py, and the sliced shapes it writes are read under their own names. self.input_categories = InputCategories.fromConfig(cfg) - unknown = self.input_categories.unknown(cfg["categories"]) + configured = list(cfg["categories"]) + unknown = self.input_categories.unknown( + configured + [self.naming.base(c) for c in configured] + ) if unknown: raise RuntimeError( - f"input_categories maps {unknown}, which are not configured categories " - f"({list(cfg['categories'])})" + f"input_categories maps {unknown}, which are neither configured " + f"categories nor their bases ({configured})" ) self.signalFractionForRelevantBins = cfg["signalFractionForRelevantBins"] From 06ad255302de803592bb1057122da0258a093905 Mon Sep 17 00:00:00 2001 From: Devin Date: Sat, 26 Sep 2026 04:40:37 +0200 Subject: [PATCH 11/12] Write only the file's own mass point's signal when rebinning The merged input for one mass carries every signal sample at every mass, but a rebinned file is only read for its own mass, since input_file_pattern names it. Other points' signal histograms, nominal and systematic, are now skipped, matched from the configuration's signal name patterns; backgrounds and data are untouched. Co-Authored-By: Claude Opus 5.5 (1M context) --- bin_opt_2d/rebin_2d.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/bin_opt_2d/rebin_2d.py b/bin_opt_2d/rebin_2d.py index adcf6bd..8be247e 100644 --- a/bin_opt_2d/rebin_2d.py +++ b/bin_opt_2d/rebin_2d.py @@ -2,6 +2,7 @@ import json import math import os +import re import sys import numpy as np import yaml @@ -1135,6 +1136,21 @@ def rebin_hist_2d(hist2d, slices, name, naming): return outputs +def other_point_signal(key, patterns, param_name, value): + """Whether `key` is a signal histogram (nominal or variation) of a model point other + than `value`. A rebinned file is only read for its own point, so those are dropped. + """ + placeholder = "${" + param_name + "}" + for pattern in patterns: + head, found, tail = pattern.partition(placeholder) + if not found: + continue + match = re.match(re.escape(head) + r"(\d+)" + re.escape(tail) + r"(?:_|$)", key) + if match and match.group(1) != str(value): + return True + return False + + def process_category( sources, channel, @@ -1267,6 +1283,10 @@ def load2d(f, key): print(f" [skip] {source_era} {in_dir}: not in the input") continue for key in [k.GetName() for k in cat_dir.GetListOfKeys()]: + if other_point_signal( + key, cfg["signal_hist_name_patterns"], param_name, mass + ): + continue hist2d = get_hist(in_file, prefix + key) if hist2d is None or hist2d.GetDimension() != 2: continue From 6824052fadc00b4eb63bb0d2e2a289bbeaecf3f5 Mon Sep 17 00:00:00 2001 From: Devin Date: Sat, 26 Sep 2026 04:41:32 +0200 Subject: [PATCH 12/12] Combine datacard configurations: per-config namespaces and a combined limit Configurations built from different inputs (the single- and double-lepton channels) cannot be summed histogram by histogram like a group era's sub-eras; they are combined card by card, in one law invocation: - datacard_config / datacard_tag on StatInferenceTask: a task can run a named configuration, with its products under //. An empty tag keeps every path as before. - A combination configuration (`members: {name: {config, hists_version, user_custom}}`). member_req() runs each member's chain under its own tag and histograms, and CombinedResonantLimitsTask combines their headline cards per mass, labelled by member, after checking they share the mass points and no bin. - A limit_plots curve may be a list of globs; in a combination each names its member, "SL:${ERA}/e/*.txt", and each member's own limit_plots are drawn as well. ResonantLimitsTask's per-mass combineCards and publish steps are factored out to be shared. Impacts on a combination are refused for now. Co-Authored-By: Claude Opus 5.5 (1M context) --- law/CombinedResonantLimitsTask.py | 84 +++++++++++++++++++++++++++++ law/CreateDatacardsTask.py | 4 +- law/MergedHists.py | 6 +++ law/PlotPullsAndImpactsTask.py | 8 ++- law/PlotResonantLimitsTask.py | 46 +++++++++++----- law/PreprocessShapesTask.py | 2 +- law/ResonantLimitsTask.py | 89 ++++++++++++++++++++----------- law/StatInferenceTask.py | 71 ++++++++++++++++++++++-- law/tasks.py | 2 + 9 files changed, 261 insertions(+), 51 deletions(-) create mode 100644 law/CombinedResonantLimitsTask.py diff --git a/law/CombinedResonantLimitsTask.py b/law/CombinedResonantLimitsTask.py new file mode 100644 index 0000000..400cff3 --- /dev/null +++ b/law/CombinedResonantLimitsTask.py @@ -0,0 +1,84 @@ +import law +import luigi +import shutil +import tempfile + +from dhi.tasks.resonant import MergeResonantLimits + +from .StatInferenceTask import StatInferenceTask +from .ResonantLimitsTask import ( + ResonantLimitsTask, + cards_by_mass, + combine_cards_per_mass, + publish_dir, +) + + +class CombinedResonantLimitsTask(StatInferenceTask): + """The limit of a combination configuration (`members:`): each member's headline + cards, combined per mass and labelled by member. Same-named nuisances correlate. + """ + + # Not a workflow itself; see ResonantLimitsTask. + workflow = luigi.Parameter(default=law.parameter.NO_STR) + + def requires(self): + if not self.is_combination(): + raise RuntimeError( + f"{self.datacard_config_path()} declares no 'members', so there is " + "nothing to combine; run ResonantLimitsTask for a single configuration." + ) + return { + name: self.member_req(name, ResonantLimitsTask) for name in self.members() + } + + def store_parts(self): + return (*self.version_parts(), self.__class__.__name__) + + def output(self): + return { + "limits": self.local_target("limits.npz"), + "datacards": law.LocalDirectoryTarget(self.datacards_dir("combined")), + } + + @staticmethod + def check_members(member_cards): + """The members must share their mass points and no bin.""" + masses = { + name: sorted(cards_by_mass(cards), key=int) + for name, cards in member_cards.items() + } + reference = next(iter(masses.values())) + if any(m != reference for m in masses.values()): + raise RuntimeError(f"the members disagree on the mass points: {masses}") + + owner = {} + for name, cards in member_cards.items(): + with open(cards[0]) as f: + # the first "bin" line lists each bin once (the second, per process) + bins = next(line.split()[1:] for line in f if line.startswith("bin ")) + for b in bins: + if owner.setdefault(b, name) != name: + raise RuntimeError( + f"bin '{b}' is in both '{owner[b]}' and '{name}'; the members' " + "channels and categories must not overlap" + ) + + def run(self): + member_cards = {} + for name, req in self.requires().items(): + cards = req.headline_cards() + if not cards: + raise RuntimeError(f"member '{name}' has no datacards to combine") + member_cards[name] = cards + self.check_members(member_cards) + + with tempfile.TemporaryDirectory() as staging: + combine_cards_per_mass(list(member_cards.items()), staging, self.cmssw_env) + combined = publish_dir(staging, self.output()["datacards"]) + + merged = yield MergeResonantLimits( + version=self.version, datacards=tuple(sorted(combined)) + ) + self.output()["limits"].parent.touch() + shutil.copy2(merged.path, self.output()["limits"].path) diff --git a/law/CreateDatacardsTask.py b/law/CreateDatacardsTask.py index 88794e1..a63fbdb 100644 --- a/law/CreateDatacardsTask.py +++ b/law/CreateDatacardsTask.py @@ -81,7 +81,9 @@ def output(self): # fs_default. Note that combine cannot read these directly: # ResonantLimitsTask mirrors them back to datacards_dir() before handing them to # dhi -- see ResonantLimitsTask.stage_datacards. - return self.output_dir_target(self.version, "Datacards", self.datacard_era) + return self.output_dir_target( + *self.version_parts(), "Datacards", self.datacard_era + ) def run(self): statInf_entry = self.global_params["StatInference"] diff --git a/law/MergedHists.py b/law/MergedHists.py index 41a5709..39f523a 100644 --- a/law/MergedHists.py +++ b/law/MergedHists.py @@ -19,6 +19,12 @@ class MergedHists(StatInferenceTask, law.ExternalTask): variable = luigi.Parameter(description="Hists_merged variable (directory) name") + # The input file is the same whichever configuration reads it. + exclude_params_req = StatInferenceTask.exclude_params_req | { + "datacard_config", + "datacard_tag", + } + def output(self): return self.remote_target( self.input_hists_version, diff --git a/law/PlotPullsAndImpactsTask.py b/law/PlotPullsAndImpactsTask.py index 771fac6..15ff254 100644 --- a/law/PlotPullsAndImpactsTask.py +++ b/law/PlotPullsAndImpactsTask.py @@ -42,10 +42,16 @@ class PlotPullsAndImpactsTask(DhiPlotMixin, StatInferenceTask): def requires(self): # ResonantLimitsTask is what writes the combined cards this reads. + if self.is_combination(): + # card_for()/card_family() only know a single configuration's layout + raise RuntimeError( + f"{self.datacard_config_path()} is a combination; impacts are not yet " + "supported for it. Run them on a member with --datacard-config." + ) return [ResonantLimitsTask.req(self)] def output(self): - return self.output_dir_target(self.version, "ImpactPlots") + return self.output_dir_target(*self.version_parts(), "ImpactPlots") def entry_eras(self, entry): """The eras an entry is drawn for. diff --git a/law/PlotResonantLimitsTask.py b/law/PlotResonantLimitsTask.py index fd41998..67fdc86 100644 --- a/law/PlotResonantLimitsTask.py +++ b/law/PlotResonantLimitsTask.py @@ -18,6 +18,7 @@ from .DhiPlotMixin import DhiPlotMixin from .StatInferenceTask import StatInferenceTask from .ResonantLimitsTask import ResonantLimitsTask +from .CombinedResonantLimitsTask import CombinedResonantLimitsTask class PlotResonantLimitsTask(DhiPlotMixin, StatInferenceTask): @@ -46,12 +47,17 @@ class PlotResonantLimitsTask(DhiPlotMixin, StatInferenceTask): def requires(self): # ResonantLimitsTask is what mirrors the cards to datacards_dir(), which is what # the globs below resolve against -- the plots cannot be built before it has run. - return [ResonantLimitsTask.req(self)] + if not self.is_combination(): + return [ResonantLimitsTask.req(self)] + # a combination also draws each member's own limit_plots + return [CombinedResonantLimitsTask.req(self)] + [ + self.member_req(name, PlotResonantLimitsTask) for name in self.members() + ] def output(self): # fs_default, like CreateDatacardsTask. Holds one # / sub-directory of plots plus plots.json naming them. - return self.output_dir_target(self.version, "LimitPlots") + return self.output_dir_target(*self.version_parts(), "LimitPlots") def _resolve_config_path(self, path): return path if os.path.isabs(path) else os.path.join(self.ana_path(), path) @@ -160,18 +166,21 @@ def datacard_sequences(self, entry, era): same cards the leading curve is drawn from. """ name = entry.get("name") or "limits" - base_dir = self.datacards_dir(era) sequences, names = [], [] for card in entry["datacards"]: - pattern = os.path.join( - base_dir, Template(card["glob"]).safe_substitute(ERA=era) - ) - if not glob.glob(pattern): - raise RuntimeError( - f"limit_plots entry '{name}': no datacard matches '{pattern}'. " - f"Check the glob against the cards under {base_dir}." - ) + # several globs: dhi combines the cards they match per mass + globs = card["glob"] if isinstance(card["glob"], list) else [card["glob"]] + patterns = [] + for g in globs: + base_dir, rel = self.glob_base(g, era, name) + pattern = os.path.join(base_dir, Template(rel).safe_substitute(ERA=era)) + if not glob.glob(pattern): + raise RuntimeError( + f"limit_plots entry '{name}': no datacard matches '{pattern}'. " + f"Check the glob against the cards under {base_dir}." + ) + patterns.append(pattern) label = card["name"] if "{" in label or "}" in label: # law brace-expands datacard_names, so "fb^{-1}" would arrive as "fb^-1". @@ -179,10 +188,23 @@ def datacard_sequences(self, entry, era): f"limit_plots entry '{name}': datacard name '{label}' contains braces, " "which law brace-expands. Put the information in the campaign label instead." ) - sequences.append((pattern,)) + sequences.append(tuple(patterns)) names.append(label) return sequences, names + def glob_base(self, card_glob, era, entry_name): + """(directory, relative glob) a limit_plots glob resolves against. In a + combination the glob names its member, "SL:${ERA}/e/*.txt".""" + if not self.is_combination(): + return self.datacards_dir(era), card_glob + member, sep, rel = card_glob.partition(":") + if not sep or member not in self.members(): + raise RuntimeError( + f"limit_plots entry '{entry_name}': glob '{card_glob}' must name a member " + f"as ':'; members are {sorted(self.members())}" + ) + return self.member_req(member, ResonantLimitsTask).datacards_dir(era), rel + def limits_npz(self, sequence): """The MergeResonantLimits .npz for a datacard sequence, produced if absent. diff --git a/law/PreprocessShapesTask.py b/law/PreprocessShapesTask.py index 01235a2..56ddbfa 100644 --- a/law/PreprocessShapesTask.py +++ b/law/PreprocessShapesTask.py @@ -63,7 +63,7 @@ def create_branch_map(self): def output(self): return self.output_dir_target( - self.version, "Hists_preprocessed", self.datacard_era + *self.version_parts(), "Hists_preprocessed", self.datacard_era ) def run(self): diff --git a/law/ResonantLimitsTask.py b/law/ResonantLimitsTask.py index 646601f..7dfb578 100644 --- a/law/ResonantLimitsTask.py +++ b/law/ResonantLimitsTask.py @@ -17,6 +17,45 @@ COMBINED_KEY = "__combined__" +def cards_by_mass(cards): + """{mass: card} from datacard paths named "..._.txt".""" + out = {} + for c in cards: + m = re.search(r"_(\d+)\.txt$", c) + if m: + out[m.group(1)] = c + return out + + +def combine_cards_per_mass(labelled_cards, staging, env): + """Write staging/combined_.txt per mass from [(label, [card, ...])]; the + unique label prefixes each card's bins.""" + by_label = [(label, cards_by_mass(cards)) for label, cards in labelled_cards] + masses = sorted({m for _, cards in by_label for m in cards}, key=int) + written = [] + for mass in masses: + args = [f"{label}={cards[mass]}" for label, cards in by_label if mass in cards] + out_file = os.path.join(staging, f"combined_{mass}.txt") + with open(out_file, "w") as f: + subprocess.run(["combineCards.py"] + args, env=env, stdout=f, check=True) + written.append(out_file) + return written + + +def publish_dir(staging, target): + """Replace directory `target` with the files in `staging`, returning their paths. + Only called once everything is staged, so a failure never looks complete to law.""" + if target.exists(): + target.remove() + target.touch() + published = [] + for name in sorted(os.listdir(staging)): + dest = os.path.join(target.path, name) + shutil.copy2(os.path.join(staging, name), dest) + published.append(dest) + return published + + class ResonantLimitsTask(StatInferenceTask): # Not a workflow itself. The parameter exists so that --workflow reaches the tasks # below that are -- law's req() only forwards parameters both tasks declare -- and @@ -24,7 +63,7 @@ class ResonantLimitsTask(StatInferenceTask): workflow = luigi.Parameter(default=law.parameter.NO_STR) def store_parts(self): - return (self.version, self.__class__.__name__, "combined") + return (*self.version_parts(), self.__class__.__name__, "combined") def get_eras(self): """Every era the configuration lists: each gets its own cards and its own limit.""" @@ -53,6 +92,18 @@ def output(self): "datacards": law.LocalDirectoryTarget(self.datacards_dir("combined")), } + def headline_cards(self): + """The per-mass cards the limit is set on: the cross-era combined ones, or the + single top-level era's.""" + top_level = self.get_top_level_eras() + if len(top_level) > 1: + return sorted( + glob.glob(os.path.join(self.datacards_dir("combined"), "*.txt")) + ) + return sorted( + glob.glob(os.path.join(self.datacards_dir(top_level[0]), "*.txt")) + ) + def stage_datacards(self, era, remote_target): """Mirror an era's datacards from fs_default to a stable local path, returned. @@ -85,13 +136,6 @@ def run(self): output_dir = self.stage_datacards(e, create_dc_br0.output()) era_cards[e] = glob.glob(os.path.join(output_dir, "*.txt")) - masses = set() - for e, cards in era_cards.items(): - for c in cards: - m = re.search(r"_(\d+)\.txt$", c) - if m: - masses.add(m.group(1)) - # Built into a staging directory and published only once every card is there. # Creating the outputs first and filling them as we go left a failed # combineCards halfway through looking like a complete task -- law's default @@ -102,29 +146,12 @@ def run(self): # The cross-era combination is the one place a group era and its members # cannot both appear: the group already *is* their combination, so a card # built from both would count those events twice. - for mass in sorted(masses): - combine_args = [] - for e in self.get_top_level_eras(): - for c in era_cards.get(e, []): - if c.endswith(f"_{mass}.txt"): - combine_args.append(f"{e}={c}") - break - - if combine_args: - cmd = ["combineCards.py"] + combine_args - out_file = os.path.join(staging, f"combined_{mass}.txt") - with open(out_file, "w") as f: - subprocess.run(cmd, env=self.cmssw_env, stdout=f, check=True) - - out_dc_dir = self.output()["datacards"] - if out_dc_dir.exists(): - out_dc_dir.remove() - out_dc_dir.touch() - combined_cards = [] - for name in sorted(os.listdir(staging)): - dest = os.path.join(out_dc_dir.path, name) - shutil.copy2(os.path.join(staging, name), dest) - combined_cards.append(dest) + combine_cards_per_mass( + [(e, era_cards.get(e, [])) for e in self.get_top_level_eras()], + staging, + self.cmssw_env, + ) + combined_cards = publish_dir(staging, self.output()["datacards"]) # One merge per era, never one merge over every era's cards together. dhi groups # the datacards it is given by the mass in their file name alone and runs diff --git a/law/StatInferenceTask.py b/law/StatInferenceTask.py index 0dab693..b7a1cff 100644 --- a/law/StatInferenceTask.py +++ b/law/StatInferenceTask.py @@ -2,6 +2,7 @@ import law import luigi import os +import re import shutil import yaml @@ -29,10 +30,26 @@ class StatInferenceTask(Task): description="version of the Hists_merged tree to read; defaults to --version", ) + # set per member by a combination configuration + datacard_config = luigi.Parameter( + default="", + description="datacard configuration to use instead of StatInference.config", + ) + # empty keeps a single configuration's paths unchanged + datacard_tag = luigi.Parameter( + default="", + description="sub-directory of --version holding this configuration's products", + ) + @property def input_hists_version(self): return self.hists_version or self.version + def version_parts(self): + """--version, then the tag ("combined" for a combination's own products).""" + tag = self.datacard_tag or ("combined" if self.is_combination() else "") + return (self.version, tag) if tag else (self.version,) + def output_dir_target(self, *path): """remote_dir_target() that also works when fs_default is a local directory. @@ -50,9 +67,8 @@ def output_dir_target(self, *path): return self.remote_dir_target(*path) def datacard_config_path(self): - return os.path.join( - self.ana_path(), self.global_params["StatInference"]["config"] - ) + path = self.datacard_config or self.global_params["StatInference"]["config"] + return path if os.path.isabs(path) else os.path.join(self.ana_path(), path) def get_config_data(self): # Cached per instance: requires()/workflow_requires() are re-entered many times @@ -62,6 +78,49 @@ def get_config_data(self): self._config_data = yaml.safe_load(f) return self._config_data + def is_combination(self): + """Whether the configuration combines others (`members:`) rather than building + datacards itself.""" + return bool(self.get_config_data().get("members")) + + def members(self): + """{name: {config, hists_version, user_custom, customisations}}, validated.""" + members = self.get_config_data().get("members") or {} + known = {"config", "hists_version", "user_custom", "customisations"} + for name, spec in members.items(): + if not isinstance(spec, dict) or "config" not in spec: + raise RuntimeError( + f"{self.datacard_config_path()}: member '{name}' must be a mapping " + "with at least 'config'" + ) + unknown = sorted(set(spec) - known) + if unknown: + raise RuntimeError( + f"{self.datacard_config_path()}: member '{name}' has unknown keys " + f"{unknown}; known are {sorted(known)}" + ) + if not re.fullmatch(r"[A-Za-z0-9_]+", str(name)) or name == "combined": + raise RuntimeError( + f"{self.datacard_config_path()}: member name '{name}' must be " + "alphanumeric and not 'combined' -- it names a directory and a " + "combineCards label" + ) + return members + + def member_req(self, name, cls, **kwargs): + """`cls` for member `name`: its own configuration, tag and histograms (FLAF + builds a Setup per user_custom). Everything else is inherited.""" + spec = self.members()[name] + params = dict( + datacard_config=spec["config"], + datacard_tag=name, + hists_version=spec.get("hists_version", self.hists_version), + user_custom=spec.get("user_custom", self.user_custom), + customisations=spec.get("customisations", self.customisations), + ) + params.update(kwargs) + return cls.req(self, **params) + @property def datacard_era(self): """The era this task's products belong to: its meta-era if it has one, else @@ -85,7 +144,7 @@ def store_parts(self): Identical to the old value for every task without a meta_era. """ - return (self.version, self.__class__.__name__, self.datacard_era) + return (*self.version_parts(), self.__class__.__name__, self.datacard_era) def datacards_dir(self, era): """Local directory holding an era's datacards. @@ -94,7 +153,9 @@ def datacards_dir(self, era): combine sees them; ResonantLimitsTask mirrors them here and everything downstream (dhi's --multi-datacards globbing, the overlay plots) resolves against this path. """ - return os.path.join(self.ana_data_path(), self.version, "Datacards", era) + return os.path.join( + self.ana_data_path(), *self.version_parts(), "Datacards", era + ) def preprocess_config(self): """The datacard configuration's `preprocess:` block, or None. diff --git a/law/tasks.py b/law/tasks.py index 8e07ebf..185e606 100644 --- a/law/tasks.py +++ b/law/tasks.py @@ -11,6 +11,7 @@ from .PreprocessShapesTask import PreprocessShapesTask from .CreateDatacardsTask import CreateDatacardsTask from .ResonantLimitsTask import ResonantLimitsTask +from .CombinedResonantLimitsTask import CombinedResonantLimitsTask from .PlotResonantLimitsTask import PlotResonantLimitsTask from .PlotPullsAndImpactsTask import PlotPullsAndImpactsTask from .ResonantLimitsAndHistPlotTask import ResonantLimitsAndHistPlotTask @@ -21,6 +22,7 @@ "PreprocessShapesTask", "CreateDatacardsTask", "ResonantLimitsTask", + "CombinedResonantLimitsTask", "PlotResonantLimitsTask", "PlotPullsAndImpactsTask", "ResonantLimitsAndHistPlotTask",