From 8a10da3b161c3c07fb6d1879a8f72133bd59db41 Mon Sep 17 00:00:00 2001 From: Devin Date: Thu, 13 Aug 2026 23:21:31 +0200 Subject: [PATCH 01/27] Add HistRebinTask, limit plotting and meta-era support Adds an in-chain rebinning and limit-plotting path to the datacard tasks: - HistRebinTask (law/tasks.py) derives 2D DNN x HME bin edges from the shapes themselves via dc_make/hist_rebin_2d.py, driven by the datacard config's `binning:` block, and writes the rebinned histograms CreateDatacardsTask builds from. With no `binning:` block the task drops out of the graph and the merged histograms are read directly, which is what the 1D DNN-score configurations need. - Sliced category names are produced and parsed by common/tools.py:CategoryNaming from a pattern the datacard configuration supplies (`binning: category_pattern`), so datacard bins and written histograms cannot fall out of step. - ResonantLimitsTask / PlotResonantLimitsTask run combine per mass point and draw the overlays declared in the config's `limit_plots`, including the banded single-curve plots. - Meta-eras are handled through the normal `eras:` list plus `era_groups:` rather than a separate mechanism (dc_make/maker.py, dc_make/uncertainty.py). - Diagnostic plots of the rebinned shapes: dc_make/plot_rebinned.py. The luminosity projection is built inline in PlotResonantLimitsTask, since dhi applies the scaling factor itself and its own --lumi-scale exists only on the single-curve plotter, where it overwrites the measured curve. common/tools.py:resolveNegativeBins and dc_make/process.py are left exactly as they are on main: negative-bin acceptance policy is deliberately out of scope here. The sliced categories this task produces are sparse enough that a few summed DY shapes go negative within their statistical error and are rejected by the existing rule, which is a binning problem to be fixed in the binning, not by loosening acceptance. With the Run3_Early input, four nominal shapes are affected (all DY, all in SR): eE/SR/boosted_dnn1 4-5 bins integral -0.99 +- 3.54 eE/SR/boosted_dnn2 1 bin integral -0.02 +- 0.18 muMu/SR/boosted_dnn2 1 bin integral -0.21 +- 1.34 muMu/SR/res2b_dnn2 4 bins integral -5.87 +- 8.92 No shape systematic variation is affected. Two of the four are single-bin categories, where no rebinning of the HME axis can help and the DNN slicing itself has to change. A process may now be scoped to a subset of categories (Process.categories, matched as a prefix so "SR/boosted" covers every DNN slice of it), and a process whose `subprocesses` name other *datacard* processes replaces them wherever it applies (DatacardMaker.mergedAwayIn / processInBin). That is what lets the bbWW DL configuration collapse TT/DY/ST/VV into one TotalBkg template in the boosted slices, where DY on its own has a median effective MC count of 5.8 and goes negative in several slices, while leaving every other category untouched. The constituents keep their own configuration and simply do not enter the merged bins. Because the merged shape is assembled from subprocess histograms by path, it also picks up DY in eMu, which DY's own `channels` list had excluded from that channel's background entirely. Inert for existing configurations: every `subprocesses` list in x_hh_bbtautau_run2.yaml names sample-level histograms (WW, WZ, ZZ, ...), none of which is a datacard process, so nothing is ever absorbed there, and an empty `categories` (the default) means every category as before. Verified by building the Run3_Early SR/boosted datacards: 59 files, no negative bins anywhere, and all 96 TotalBkg rates (32 per channel) equal the raw TT+DY+ST+VV histogram sums to 1e-4 -- nothing was clamped or donor balanced. The DNN slice boundaries gain per-process floors (min_dnn_bkg_each, min_dnn_bkg_neff), the analogues of what _mass_passes already applied on the mass axis. _dnn_passes tested only the summed background, which never binds in practice -- over the 232 eE/muMu SR slices of Run3_Early the summed N_eff never once fell below the configured 4, while DY individually sat below it in 42% of them. A background that has fluctuated negative is invisible inside a healthy total, and the significance being maximised, S/sqrt(B + sigma_B^2), actively prefers it because a lower B raises the score: muMu/SR/res2b at m500 chose left edge 68 (Z 0.834, DY -5.87 +- 8.92) over 67 (Z 0.775, DY +7.31). The resulting slice cannot be made into a datacard at all. min_dnn_bkg_each: 0.01 removes all 7 negative-background slices in Run3_Early for -0.16% of the total combined Asimov Z (39.103 -> 39.042); worst single category muMu/res2b at m500, -2.3%. min_dnn_bkg_neff is left at 0: requiring those backgrounds to be measured rather than merely positive costs 35% of the Z there, because the whole signal-like region of that category holds only N_eff 5.2 of DY and cannot be split three ways above a threshold of 4. Backgrounds below min_bkg_frac of the category total are exempt from both, as on the mass axis, so a negligible process cannot veto every boundary. The DNN selection each slice stands for is the slice directory's own title, not a side-car dnn_slices.json. It is otherwise nowhere in the output -- the sliced categories are named by index and the surviving axis is the HME one -- so the plotter cannot re-derive it: the edges exist only on the pre-rebin 2D x axis, and plot_rebinned.py runs in the default env for FLAF's PlotKit, without ROOT. Carrying the label inside the artifact removes the requirement that the reader be handed exactly the right directory, which was satisfied only because CreateDatacardsTask (not HistRebinTask) invokes the plotter with self.input()[0].parent; a wrong path silently dropped every label. format_var_range moves to hist_rebin_2d.py accordingly -- the producer formats once, the consumer displays a string -- and using the directory's own title leaves no key name for the two ends to agree on at all. It takes the axis name from the configuration's `slice_var` rather than defaulting to "DNN". Directories are now created a level at a time (mkdir_titled). TDirectory.mkdir() given a slashed path puts the title on the *first* level and hands the rest of the path down as the sub-levels' titles, which is why the output had "muMu" titled "muMu/SR/res2b_dnn0"; the parents now carry their own names. common/tools.py keeps only what the rebinned output's layout genuinely shares, as one class: CategoryNaming, which both writes a sliced category name and reads it back. It replaces a format string that was authoritative for writing while four hardcoded regexes parsed it back (maker.py:getCategoryGroups, plot_rebinned.py x3) -- changing the format used to break grouping silently, the regex stopping to match and every slice becoming its own base category rather than raising. Both directions are now derived from one pattern, and that pattern is configuration rather than code: `binning: category_pattern`, defaulting to "{base_category}_slice{slice_idx}". Nothing in this repository assumes the sliced axis is a DNN score; HH->bbWW's own configuration pins "{base_category}_dnn{slice_idx}" and `slice_var: DNN`, which is what keeps the names in this analysis's output unchanged. "_slice" rather than "_cat" for the neutral default because hand-written category names do use "_cat" (x_hh_bbtautau_run2.yaml has "res1b_cat3_masswindow"), and a default that could collide with an unsliced name would have split() claiming it as a slice. The pattern is a BINNING_PARAM like the rest, so HistRebinTask resolves it once and hands the same value to create_datacards.py and plot_rebinned.py, on the same rule --n-dnn-slices already followed: a CLI override cannot leave a reader deriving names the rebinned files do not contain. Verified by re-running the rebin for all four sub-eras and rebuilding: 40 files, no dnn_slices.json, slice labels matching the old JSON edge-for-edge, datacard .txt byte-identical and all 40992 shape histograms unchanged against the previously validated build. Rebased onto the black formatting of #21, so this diff is code only. Co-Authored-By: Claude Opus 5 --- README.md | 41 ++ common/tools.py | 105 ++++ dc_make/create_datacards.py | 22 + dc_make/hist_rebin_2d.py | 1087 +++++++++++++++++++++++++++++++++++ dc_make/maker.py | 719 +++++++++++++++++++---- dc_make/plot_rebinned.py | 566 ++++++++++++++++++ dc_make/process.py | 17 + dc_make/uncertainty.py | 13 +- law/tasks.py | 1086 ++++++++++++++++++++++++++++++++-- 9 files changed, 3481 insertions(+), 175 deletions(-) create mode 100644 dc_make/hist_rebin_2d.py create mode 100644 dc_make/plot_rebinned.py diff --git a/README.md b/README.md index 9d89ce8..e3cb58b 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,46 @@ # StatInference +## Two ways bins get decided + +Two independent things in this repository choose the binning of the shapes that go into +the datacards. They do not share code, and an analysis uses one or the other. + +**1. In-chain rebinning — `dc_make/hist_rebin_2d.py`.** A production step, run by +`HistRebinTask` in `law/tasks.py`: it derives the bin edges from the shapes themselves +(no fits, no limits) and writes the rebinned histograms the datacards are built from, so +`CreateDatacardsTask` cannot run without it. Its knobs come from the `binning:` block of +the datacard configuration — see the annotated block in the HH→bbWW Run 3 configuration, +`config/Datacards/x_hh_bbww_DL_run3.yaml` in the analysis area. It cuts each base +category (`SR/res2b`) into per-slice categories whose names come from that block's +`category_pattern`, e.g. `{base_category}_dnn{slice_idx}` -- the pattern is the +analysis's choice, and nothing in this repository assumes the sliced axis is a DNN +score. `common/tools.py:CategoryNaming` both writes those names and parses them back +from the one pattern, and every consumer of the configuration's `categories` list +expands them through it, so the datacard bins and the written histograms cannot fall out +of step. + +**2. Offline binning optimisation — `bin_opt/`.** A search harness, documented below: it +builds candidate binnings, runs limits with combine for each, and ranks them. Its product +is a `hist_bins` JSON, applied at datacard time by `dc_make/binner.py`. It is not part of +the datacard chain and exposes no importable API — every module is a script driven by +`bin_opt/bin_optimization.yaml`. + +Which one an analysis is on is visible in its configuration: the in-chain path has a +`binning:` block and leaves `hist_bins` unset (as HH→bbWW's `config/global.yaml` does), +while the `bin_opt` path sets `hist_bins` and has no `binning:` block (as +`config/x_hh_bbtautau_run2.yaml` does). + +That difference is also what decides how `categories:` is read. With a `binning:` block +(or an explicit `--n-dnn-slices`), `dc_make` expands the listed base categories into the +sliced names above, because those are the directories `HistRebinTask` wrote. Without one, +the categories are used exactly as listed — which is what an input file that is already +1D and binned contains. Pass `--n-dnn-slices 0` to force the latter for a configuration +that does carry a `binning:` block. + +In the law chain the same block decides the graph: with no `binning:` block +`HistRebinTask` drops out entirely and `CreateDatacardsTask` reads the merged histograms +straight from the `Hists_merged` tree, since there is no step in between. + ## How to run binning optimisation on lxplus Open two separate LXPLUS terminals: one for **server side** and the other for **worker side**. Set up the environment and proxy in analysis area as usual on both terminals. diff --git a/common/tools.py b/common/tools.py index 930b6df..40664ab 100644 --- a/common/tools.py +++ b/common/tools.py @@ -1,7 +1,112 @@ import math +import re + +from string import Formatter + import numpy as np +class CategoryNaming: + """Names the sub-categories a 2D -> categorized-1D transformation cuts a base + category into, and reads those names back. + + The pattern is configuration rather than code: nothing here assumes the sliced axis + is a DNN score, or that the analysis producing the 2D shapes has one at all. Both + directions are built from that one pattern, so the code that writes a category name + and the code that parses it back cannot drift apart -- which is the failure this + replaces, where a format string wrote the names and hardcoded regexes read them, and + changing the format silently stopped the regexes matching, leaving every slice as its + own base category instead of raising. + """ + + # Deliberately neutral: an analysis that wants its discriminant in the name says so in + # its configuration (HH->bbWW pins "{base_category}_dnn{slice_idx}"), rather than every + # analysis inheriting one analysis's choice from here. "_slice" rather than "_cat" + # because hand-written category names do use "_cat" -- x_hh_bbtautau_run2.yaml has + # "res1b_cat3_masswindow" -- and a default that can collide with an unsliced name would + # have split() claiming it as a slice of something. + default_pattern = "{base_category}_slice{slice_idx}" + + # Sub-expression per placeholder. base_category is greedy so a base name that itself + # ends in something the pattern could match still resolves to the last slice index, + # which is the one the pattern wrote. + placeholders = { + "base_category": r"(?P.+)", + "slice_idx": r"(?P\d+)", + } + + def __init__(self, pattern=None): + self.pattern = pattern or self.default_pattern + regex = "" + seen = set() + for literal, field, _, conversion in Formatter().parse(self.pattern): + regex += re.escape(literal) + if field is None: + continue + if field not in self.placeholders: + raise RuntimeError( + f"category pattern '{self.pattern}': unknown placeholder " + f"'{{{field}}}'; known are " + + ", ".join("'{%s}'" % p for p in sorted(self.placeholders)) + ) + if field in seen: + raise RuntimeError( + f"category pattern '{self.pattern}': '{{{field}}}' appears more than " + "once, so a name built from it cannot be read back unambiguously" + ) + if conversion is not None: + raise RuntimeError( + f"category pattern '{self.pattern}': conversion '!{conversion}' on " + f"'{{{field}}}' would not survive the round trip back to a category" + ) + seen.add(field) + regex += self.placeholders[field] + missing = sorted(set(self.placeholders) - seen) + if missing: + raise RuntimeError( + f"category pattern '{self.pattern}' must use every placeholder; missing " + + ", ".join("'{%s}'" % m for m in missing) + ) + self._regex = re.compile("^" + regex + "$") + + @classmethod + def fromConfig(cls, cfg): + """From a datacard configuration's ``binning:`` block, or the default pattern + when it declares none (or when there is no ``binning:`` block at all, i.e. the + input is already 1D and nothing was ever sliced).""" + return cls((cfg.get("binning") or {}).get("category_pattern")) + + def name(self, base_category, slice_idx): + return self.pattern.format(base_category=base_category, slice_idx=slice_idx) + + def expand(self, base_categories, n_slices): + """Base category names -> the per-slice names that exist in the rebinned files. + + Every consumer of the datacard configuration's ``categories`` list needs the same + expansion (maker.py's datacard bins, plot_rebinned.py's panels), and hist_rebin_2d.py + writes with the same pattern, so it is done in one place. + """ + return [ + self.name(base, idx) for base in base_categories for idx in range(n_slices) + ] + + def split(self, category): + """Inverse of name(): "SR/res2b_dnn2" -> ("SR/res2b", 2). + + An unsliced name returns (category, None) rather than raising -- a configuration + with no ``binning:`` block has categories that were never sliced, and they are + legitimately their own base. + """ + match = self._regex.match(category) + if not match: + return category, None + return match.group("base_category"), int(match.group("slice_idx")) + + def base(self, category): + """The base category a slice belongs to, or the name itself if unsliced.""" + return self.split(category)[0] + + class PackageWrapper: def __init__(self, import_fn): self._package = None diff --git a/dc_make/create_datacards.py b/dc_make/create_datacards.py index 192e984..b1c26e7 100644 --- a/dc_make/create_datacards.py +++ b/dc_make/create_datacards.py @@ -34,6 +34,26 @@ default=None, help="parameter values to run only certain masses", ) + parser.add_argument( + "--n-dnn-slices", + required=False, + type=int, + default=None, + help="DNN slices each base category was cut into by HistRebinTask; " + "defaults to the config's binning block. Pass the value that " + "task actually used, so the datacard bins match the input files. " + "0 -- or an unset value with no binning block in the config -- means the " + "input is already binned and the config's category names are used verbatim", + ) + parser.add_argument( + "--category-pattern", + required=False, + type=str, + default=None, + help="pattern HistRebinTask named those slices with; defaults to the config's " + "binning block. Pass the value that task actually used, so the datacard bins " + "carry the names the input files do", + ) for param in DatacardMaker.customizeble_parameters: parser.add_argument( @@ -68,6 +88,8 @@ args.input, hist_bins=hist_bins, param_values=param_values, + n_dnn_slices=args.n_dnn_slices, + category_pattern=args.category_pattern, **kwargs, ) maker.createDatacards(args.output) diff --git a/dc_make/hist_rebin_2d.py b/dc_make/hist_rebin_2d.py new file mode 100644 index 0000000..1fbeac5 --- /dev/null +++ b/dc_make/hist_rebin_2d.py @@ -0,0 +1,1087 @@ +import array +import math +import os +import sys +import yaml + +if __name__ == "__main__": + file_dir = os.path.dirname(os.path.abspath(__file__)) + pkg_dir = os.path.dirname(file_dir) + base_dir = os.path.dirname(pkg_dir) + pkg_dir_name = os.path.split(pkg_dir)[1] + if base_dir not in sys.path: + sys.path.append(base_dir) + __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 + +ROOT = importROOT() + + +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"] + # cfg["categories"] holds base category names ("SR/res2b"), which is exactly what + # the raw 2D histograms are keyed by. The per-slice names this script writes are + # built by the CategoryNaming run() installs below. + categories = list(dict.fromkeys(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 DNN 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, + "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 DNN 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 _dnn_passes( + yields, + min_sum, + total_error=None, + min_neff=0.0, + errors=None, + min_each=0.0, + min_proc_neff=0.0, + exempt=(), +): + """DNN-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 _mass_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 DNN 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 + if min_each > 0 or min_proc_neff > 0: + for name, value in yields.items(): + if name in exempt: + continue + if 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 _mass_passes(yields, min_each, exempt=(), total_error=None, min_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. + + `exempt` is the set of backgrounds judged negligible across the whole slice by + minor_backgrounds(); they are excused from min_each so that a process worth + 0.4% of the yield -- and known only to a few hundred percent -- cannot veto + every candidate split, which would make find_mass_bins() back off all the way + to a single bin and discard the mass shape in exactly the high-significance + slices that matter most. + + The min_neff arm is the same gate _dnn_passes() applies to slice boundaries. + Applying it only there left the mass bins *inside* each slice ungated, and + since that is where essentially every fit bin lives, the median per-bin + effective background count came out at ~2.8 against a slice threshold of 4. + Note it constrains only the *summed* background, so it is satisfied by any one + well-measured process; min_each/exempt is what protects the individual ones. + """ + total = sum(yields.values()) + 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(): + if name in exempt: + continue + if value <= min_each: + return False + return True + + +def grow_dnn_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 _dnn_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_dnn_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 _dnn_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: + lo = first_bin if right >= first_bin else right + slices.append((lo, right if right >= lo else lo)) + right = lo - 1 + continue + left = grow_dnn_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. + + This is what puts the bins where the resonance is. The HME axis is 0-1500 GeV + for every mass point, so the signal occupies a narrow window whose position + moves with MX while the axis does not; binning must follow the signal rather + than the axis. Equal-signal quantiles do that automatically -- bin edges + cluster wherever dS/dm is large (the peak) and a single wide bin absorbs the + long empty stretches on either side. + + The CDF clamps negative bin contents to zero so it stays monotonic; a + statistical undershoot in a signal MC bin must not move an edge backwards. + """ + n_avail = last_bin - first_bin + 1 + n_bins = max(1, min(n_bins, n_avail)) + if n_bins == 1: + return [(first_bin, last_bin)] + + cumulative = [] + running = 0.0 + for b in range(first_bin, last_bin + 1): + running += max(sig_hist.GetBinContent(b), 0.0) + cumulative.append(running) + total = running + if total <= 0: + return [(first_bin, last_bin)] + + ranges = [] + lo = first_bin + for k in range(1, n_bins): + target = k * total / n_bins + b = lo + while b < last_bin and cumulative[b - first_bin] < target: + b += 1 + # every one of the n_bins-k ranges still to come needs at least one bin + b = max(lo, min(b, last_bin - (n_bins - k))) + ranges.append((lo, b)) + lo = b + 1 + ranges.append((lo, last_bin)) + return ranges + + +def merge_until_valid( + ranges, sig_hist, bkg_hists_by_name, min_each, exempt=(), min_neff=0.0 +): + """Merge adjacent ranges until every one satisfies the background gates. + + The signal quantiles decide where the edges want to be; this decides how many + of them the background statistics can actually support. A failing range is + merged into whichever neighbour holds *less* signal, so the dense bins around + the peak -- the ones carrying the discrimination -- are the last to be given + up. Terminates because every step removes one range, ending at the single + full-range bin, which has nothing left to fail against. + """ + ranges = list(ranges) + while len(ranges) > 1: + bad = None + for i, (lo, hi) in enumerate(ranges): + if not _mass_passes( + _bkg_yields(bkg_hists_by_name, lo, hi), + min_each, + exempt, + _total_bkg_error(bkg_hists_by_name, lo, hi), + min_neff, + ): + bad = i + break + if bad is None: + break + if bad == 0: + other = 1 + elif bad == len(ranges) - 1: + other = bad - 1 + else: + left_sig = sig_hist.Integral(*ranges[bad - 1]) + right_sig = sig_hist.Integral(*ranges[bad + 1]) + other = bad - 1 if left_sig <= right_sig else bad + 1 + first, second = min(bad, other), max(bad, other) + ranges[first : second + 1] = [(ranges[first][0], ranges[second][1])] + return ranges + + +def find_mass_bins( + sig_hist, + bkg_hists_by_name, + max_bins, + first_bin, + last_bin, + min_each, + min_frac=0.0, + min_neff=0.0, +): + """Mass bins inside one DNN slice: signal quantiles for the edges, background + gates for the count. + + The previous version scanned right-to-left growing each bin leftward until the + backgrounds cleared their floors, mirroring find_dnn_slices(). That is right + for the DNN axis (signal-like = high score, so resolution belongs at the top) + and wrong for HME, where the signal sits at HME ~ MX and both tails are empty. + Starting from the top of the axis spent the bin budget on the empty + high-HME tail and left the entire resonance peak in the single leftover bin: + at MX=600 in muMu/res2b_dnn3, one bin covered 0-710 GeV holding 97.5% of the + signal while five bins shared the 710-1500 GeV region holding 2.5%. The fit + then had no shape to work with in the only region where signal and background + differ. + """ + # which backgrounds count as negligible is decided once, over the whole slice + 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 + ) + + +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.""" + ranges = list(ranges) + ranges[0] = (full_lo, ranges[0][1]) + ranges[-1] = (ranges[-1][0], full_hi) + return ranges + + +def mass_bin_budget(bkg_hists_by_name, lo, hi, max_mass_bins, bkg_per_mass_bin): + """How many mass bins this DNN slice can actually afford. + + A fixed max_mass_bins 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_mass_bin ties the + binning to the statistics that are really there. Returns max_mass_bins + unchanged when bkg_per_mass_bin <= 0 (feature off). + """ + if bkg_per_mass_bin <= 0: + return max_mass_bins + total = sum(_bkg_yields(bkg_hists_by_name, lo, hi).values()) + if total <= 0: + return 1 + return max(1, min(max_mass_bins, int(total / bkg_per_mass_bin))) + + +def discover_binning( + sig2d, + bkg2d_by_name, + n_dnn_slices, + max_mass_bins, + min_dnn_sum, + min_mass_each, + min_bkg_neff=0.0, + min_bkg_frac=0.0, + min_mass_bkg_neff=0.0, + bkg_per_mass_bin=0.0, + sig_mode="sb", + min_dnn_bkg_each=0.0, + min_dnn_bkg_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 + meta-era's sub-eras (combined limit) -- see HistRebinTask.get_discovery_eras(). + Yields are summed across whatever's in the list; see _bkg_yields(). sig2d is + the same discovery reference's (already-summed) signal histogram: its x + projection picks the significance-maximizing DNN slice boundaries, and its y + projection within each slice places the mass bin edges by signal quantile.""" + any_hist = next(iter(bkg2d_by_name.values()))[0] + nx = any_hist.GetNbinsX() + ny = any_hist.GetNbinsY() + dnn_slices = find_dnn_slices( + sig2d, + bkg2d_by_name, + n_dnn_slices, + 1, + nx, + min_dnn_sum, + min_bkg_neff, + sig_mode, + min_dnn_bkg_each, + min_dnn_bkg_neff, + min_bkg_frac, + ) + dnn_slices = extend_outer_edges(dnn_slices, 0, nx + 1) + + result = [] + for xlo, xhi in dnn_slices: + # ProjectionY attaches its result to gDirectory (here, the open output + # file); detach so these die with the loop iteration instead of piling up + # in the output file's in-memory object list. + bkg_y_by_name = { + name: [ + _detach(h.ProjectionY(f"_disc_{name}_{xlo}_{xhi}_{i}_y", xlo, xhi, "e")) + for i, h in enumerate(hists) + ] + for name, hists in bkg2d_by_name.items() + } + sig_y = _detach(sig2d.ProjectionY(f"_disc_sig_{xlo}_{xhi}_y", xlo, xhi, "e")) + n_bins = mass_bin_budget(bkg_y_by_name, 1, ny, max_mass_bins, bkg_per_mass_bin) + mass_ranges = find_mass_bins( + sig_y, + bkg_y_by_name, + n_bins, + 1, + ny, + min_mass_each, + min_bkg_frac, + min_mass_bkg_neff, + ) + mass_ranges = extend_outer_edges(mass_ranges, 0, ny + 1) + result.append({"x_range": (xlo, xhi), "y_ranges": mass_ranges}) + return result + + +def mass_bin_edges(y_axis, y_ranges): + """Physical HME 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 HME 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. + + TDirectory.mkdir() given a slashed path applies the title to the *first* level and + hands the rest of the path down as the sub-levels' titles, so every parent ends up + labelled with a stale fragment of whichever slice was created first -- the output + currently has "muMu" titled "muMu/SR/res2b_dnn0". Walking the components keeps the + parents clean and puts the label where it belongs. + """ + parts = path.split("/") + for part in parts[:-1]: + directory = directory.GetDirectory(part) or directory.mkdir(part) + return directory.mkdir(parts[-1], title) + + +def format_var_range(lo, hi, var): + """One slice's edges on the sliced axis as a selection label, e.g. "1.20 < DNN < 4.50". + + `var` names that axis and comes from the configuration -- this script does not assume + the analysis slices on a DNN score. + + Formatted here, by the code that discovered the edges, so the label travels with the + histograms and nothing downstream has to re-derive it. Both edges open means the slice + covers the whole axis, i.e. there is no selection to state -- that returns an empty + string rather than a vacuous label. + """ + if lo is None and hi is None: + return "" + if lo is None: + return f"{var} < {hi:.2f}" + if hi is None: + return f"{var} > {lo:.2f}" + return f"{lo:.2f} < {var} < {hi:.2f}" + + +def slice_ranges(x_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 + outermost ones into underflow/overflow -- those have no finite edge, so they are + recorded as null and read back as an open-ended selection. Written out by run() so + the plots can say which DNN 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"] + ranges.append( + [ + None if lo < 1 else x_axis.GetBinLowEdge(lo), + None if hi > n else x_axis.GetBinUpEdge(hi), + ] + ) + return ranges + + +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"] + edges = array.array("d", mass_bin_edges(hist2d.GetYaxis(), sl["y_ranges"])) + # Detached for the same reason as the projections above: Write() targets + # gDirectory regardless, so nothing needs these to stay attached. + h = _detach( + ROOT.TH1D( + naming.name(name, slice_idx), + name, + len(edges) - 1, + edges, + ) + ) + for bin_idx, (ylo, yhi) in enumerate(sl["y_ranges"], start=1): + err = array.array("d", [0.0]) + content = hist2d.IntegralAndError(xlo, xhi, ylo, yhi, err) + h.SetBinContent(bin_idx, content) + h.SetBinError(bin_idx, err[0]) + outputs.append(h) + return outputs + + +def process_category( + in_file, + out_file, + channel, + category, + cfg, + mass, + era, + discovery_files, + n_dnn_slices, + max_mass_bins, + min_dnn_sum, + min_mass_each, + min_signal, + min_bkg_neff=0.0, + min_bkg_frac=0.0, + min_mass_bkg_neff=0.0, + bkg_per_mass_bin=0.0, + sig_mode="sb", + min_dnn_bkg_each=0.0, + min_dnn_bkg_neff=0.0, +): + prefix = f"{channel}/{category}/" + cat_dir = in_file.Get(f"{channel}/{category}") + if not cat_dir: + print(f" [skip] {channel}/{category}: not found in {in_file.GetName()}") + return + keys = [k.GetName() for k in cat_dir.GetListOfKeys()] + + signal_keys = [ + applyParameters(pattern, {cfg["signal_param_name"]: mass}) + for pattern in cfg["signal_hist_name_patterns"] + ] + background_names = [ + hist_name + for (base_name, hist_name, allowed_channels) in cfg["background_entries"] + if not allowed_channels or channel in allowed_channels + ] + + def load2d(f, key): + return get_hist(f, prefix + key) + + disc_sig = sum_hists( + [load2d(f, key) for f in discovery_files for key in signal_keys] + ) + disc_bkg_by_name = {} + for bkg_key in background_names: + per_era = [load2d(f, bkg_key) for f in discovery_files] + per_era = [h for h in per_era if h is not None] + if len(per_era) == len(discovery_files): + disc_bkg_by_name[bkg_key] = per_era + + sig_integral = disc_sig.Integral() if disc_sig is not None else 0 + if disc_sig is None or sig_integral < min_signal or not disc_bkg_by_name: + if not disc_bkg_by_name: + reason = "no background histograms found in all discovery eras" + else: + reason = f"signal too small for discovery ({sig_integral} < {min_signal})" + print(f" [skip] {channel}/{category} MX={mass}: {reason}, skipping") + return + + slices = discover_binning( + disc_sig, + disc_bkg_by_name, + n_dnn_slices, + max_mass_bins, + min_dnn_sum, + min_mass_each, + min_bkg_neff, + min_bkg_frac, + min_mass_bkg_neff, + bkg_per_mass_bin, + sig_mode, + min_dnn_bkg_each, + min_dnn_bkg_neff, + ) + # The DNN selection each slice stands for is otherwise nowhere in the output: the + # sliced categories are named by index and the surviving axis is the HME one. It is + # the slice directory's own title, so it travels with the histograms and there is no + # 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) + for slice_idx in range(len(slices)): + mkdir_titled( + out_file, + f"{channel}/{naming.name(category, slice_idx)}", + format_var_range(*ranges[slice_idx], var=cfg["slice_var"]), + ) + + for key in keys: + hist2d = load2d(in_file, key) + if hist2d is None or hist2d.GetDimension() != 2: + continue + rebinned = rebin_hist_2d(hist2d, slices, key, naming) + for slice_idx, h in enumerate(rebinned): + out_cat = naming.name(category, slice_idx) + out_file.cd(f"{channel}/{out_cat}") + h.Write(key) + + +def run( + input_dir, + output_dir, + config_path, + era, + discovery_eras, + n_dnn_slices, + max_mass_bins, + min_dnn_sum, + min_mass_each, + min_signal, + min_bkg_neff=0.0, + min_bkg_frac=0.0, + min_mass_bkg_neff=0.0, + bkg_per_mass_bin=0.0, + sig_mode="sb", + min_dnn_bkg_each=0.0, + min_dnn_bkg_neff=0.0, + category_pattern=None, + slice_var="x", +): + cfg = load_config(config_path) + # Taken from the command line rather than from the configuration just read, so that an + # override reaches the names actually written. HistRebinTask resolves the value once + # and hands the same one to every reader, so the written names and the datacard bins + # cannot come from different patterns. + cfg["naming"] = CategoryNaming(category_pattern) + cfg["slice_var"] = slice_var + model = cfg["model"] + + for mass in cfg["mass_values"]: + file_name = model.getInputFileName(era, {cfg["signal_param_name"]: mass}) + in_file = open_input_file(input_dir, model, era, mass, cfg["signal_param_name"]) + discovery_files = [ + open_input_file(input_dir, model, disc_era, mass, cfg["signal_param_name"]) + for disc_era in discovery_eras + ] + + # output_dir is already period-scoped (HistRebinTask.output() = HistRebin/), + # so strip the leading "/" that getInputFileName() adds -- otherwise the era + # ends up doubled in the output path. + era_prefix = era + "/" + out_rel = ( + file_name[len(era_prefix) :] + if file_name.startswith(era_prefix) + else file_name + ) + out_path = os.path.join(output_dir, out_rel) + os.makedirs(os.path.dirname(out_path), exist_ok=True) + out_file = ROOT.TFile.Open(out_path, "RECREATE") + + print(f"Rebinning era={era} MX={mass} -> {out_path}") + for channel in cfg["channels"]: + for category in cfg["categories"]: + process_category( + in_file, + out_file, + channel, + category, + cfg, + mass, + era, + discovery_files, + n_dnn_slices, + max_mass_bins, + min_dnn_sum, + min_mass_each, + min_signal, + min_bkg_neff, + min_bkg_frac, + min_mass_bkg_neff, + bkg_per_mass_bin, + sig_mode, + min_dnn_bkg_each, + min_dnn_bkg_neff, + ) + + out_file.Close() + in_file.Close() + for f in discovery_files: + f.Close() + + +if __name__ == "__main__": + import argparse + + parser = argparse.ArgumentParser( + description="Rebin 2D (DNN x HME) histograms into significance-sliced 1D shapes." + ) + parser.add_argument( + "--input", + required=True, + type=str, + help="base directory containing //.root", + ) + parser.add_argument( + "--output", + required=True, + type=str, + help="output base directory, mirrors --input layout", + ) + parser.add_argument( + "--config", required=True, type=str, help="datacard configuration yaml" + ) + parser.add_argument( + "--era", required=True, type=str, help="era to rebin and write out" + ) + parser.add_argument( + "--discovery-eras", + required=False, + type=str, + default=None, + help="comma-separated eras summed to discover bin edges (defaults to --era alone)", + ) + parser.add_argument( + "--category-pattern", + required=False, + type=str, + default=None, + help="pattern naming the categories a base category is sliced into, e.g. " + f"'{CategoryNaming.default_pattern}' (the default). Must use both " + "{base_category} and {slice_idx}; every reader of these shapes is handed the " + "same pattern, and it is what they parse the names back with", + ) + parser.add_argument( + "--slice-var", + required=False, + type=str, + default="x", + help="name of the sliced axis, used only to label each slice directory with the " + "selection it stands for (e.g. 'DNN' -> '0.80 < DNN < 1.00')", + ) + parser.add_argument( + "--n-dnn-slices", + required=False, + type=int, + default=4, + help="fixed number of DNN slices per category (must be the same for every mass point)", + ) + parser.add_argument( + "--max-mass-bins", + required=False, + type=int, + default=10, + help="mass bins to aim for inside each DNN slice; backed off until achievable", + ) + parser.add_argument( + "--min-dnn-bkg-sum", + required=False, + type=float, + default=1.0, + help="minimum summed-background yield required in a DNN slice", + ) + parser.add_argument( + "--min-mass-bkg-each", + required=False, + type=float, + default=0.01, + help="minimum yield required of every individual background in a mass bin", + ) + parser.add_argument( + "--min-bkg-neff", + required=False, + type=float, + default=0.0, + help="minimum effective MC entries ((sum/err)^2) required of the summed " + "background for a DNN slice boundary to be selectable", + ) + parser.add_argument( + "--min-bkg-frac", + required=False, + type=float, + default=0.0, + help="backgrounds contributing less than this fraction of the total are " + "exempt from --min-mass-bkg-each, so a statistically starved minor " + "process cannot collapse a slice to a single mass bin", + ) + parser.add_argument( + "--min-signal", + required=False, + type=float, + default=0.5, + help="minimum discovery signal yield required to slice a category at all " + "(below this, e.g. boosted at low MX, the category is skipped entirely)", + ) + parser.add_argument( + "--min-mass-bkg-neff", + required=False, + type=float, + default=0.0, + help="minimum effective MC entries required of the summed background in " + "every mass bin. The --min-bkg-neff analogue for the bins inside a " + "slice, which is where essentially every fit bin lives", + ) + parser.add_argument( + "--bkg-per-mass-bin", + required=False, + type=float, + default=0.0, + help="target summed-background yield per mass bin; caps the bin count at " + "B_slice/this instead of always using --max-mass-bins. 0 disables, " + "restoring the fixed --max-mass-bins for every slice", + ) + parser.add_argument( + "--min-dnn-bkg-each", + required=False, + type=float, + default=0.0, + help="minimum yield required of every non-negligible background in a DNN " + "slice. The --min-mass-bkg-each analogue for the slice boundaries " + "themselves, which were previously gated on the summed background only: " + "a background that had fluctuated negative was hidden inside a healthy " + "total, and the resulting slice could not be made into a datacard. " + "Backgrounds below --min-bkg-frac of the category total are exempt. " + "0 disables", + ) + parser.add_argument( + "--min-dnn-bkg-neff", + required=False, + type=float, + default=0.0, + help="minimum effective MC entries required of every non-negligible " + "background in a DNN slice (same exemption as --min-dnn-bkg-each). Much " + "stronger than --min-dnn-bkg-each and correspondingly expensive: on " + "Run3_Early a threshold of 4 cost 35% of the combined Asimov Z in " + "muMu/SR/res2b at m500, against 2.3% for the yield floor alone. 0 disables", + ) + parser.add_argument( + "--significance-mode", + required=False, + type=str, + default="sb", + choices=["sb", "asimov"], + help="figure of merit for DNN slice boundaries: 'sb' = S/sqrt(B+sigmaB^2), " + "'asimov' = Poisson-correct Asimov significance (valid at low B)", + ) + args = parser.parse_args() + + discovery_eras = ( + args.discovery_eras.split(",") if args.discovery_eras else [args.era] + ) + + run( + args.input, + args.output, + args.config, + args.era, + discovery_eras, + args.n_dnn_slices, + args.max_mass_bins, + args.min_dnn_bkg_sum, + args.min_mass_bkg_each, + args.min_signal, + args.min_bkg_neff, + args.min_bkg_frac, + args.min_mass_bkg_neff, + args.bkg_per_mass_bin, + args.significance_mode, + args.min_dnn_bkg_each, + args.min_dnn_bkg_neff, + category_pattern=args.category_pattern, + slice_var=args.slice_var, + ) diff --git a/dc_make/maker.py b/dc_make/maker.py index a5b6902..274c7e5 100644 --- a/dc_make/maker.py +++ b/dc_make/maker.py @@ -11,6 +11,7 @@ importROOT, resolveNegativeBins, getRelevantBins, + CategoryNaming, ) from .process import Process from .uncertainty import ( @@ -18,6 +19,8 @@ UncertaintyType, UncertaintyScale, MultiValueLnNUncertainty, + LnNUncertainty, + ShapeUncertainty, ) from .model import Model from .binner import Binner @@ -29,7 +32,14 @@ class DatacardMaker: customizeble_parameters = ["eras", "channels", "categories"] def __init__( - self, cfg_file, input_path, hist_bins=None, param_values=None, **kwargs + self, + cfg_file, + input_path, + hist_bins=None, + param_values=None, + n_dnn_slices=None, + category_pattern=None, + **kwargs, ): self.cb = CombineHarvester() @@ -53,9 +63,35 @@ def __init__( self.analysis = cfg["analysis"] self.eras = cfg["eras"] self.channels = cfg["channels"] - self.categories = cfg["categories"] + # For input produced by HistRebinTask the configuration lists base categories + # ("SR/res2b") while the datacard bins are the per-DNN-slice names it wrote. + # n_dnn_slices comes from the caller when given -- HistRebinTask's resolved value + # is authoritative, since a command-line override would otherwise leave this + # deriving a different slice count than the rebinned files contain. + # + # No slice count anywhere (no `binning:` block and no caller value) means the + # input did not come from HistRebinTask -- it is already binned, and its + # categories are exactly the ones the configuration lists. An explicit 0 says the + # same for a configuration that does carry a `binning:` block. + n_slices = n_dnn_slices + if n_slices is None: + n_slices = (cfg.get("binning") or {}).get("n_dnn_slices") + # Same rule for the pattern that names those slices: the caller's resolved value + # wins, otherwise the configuration's own. + self.naming = ( + CategoryNaming(category_pattern) + if category_pattern + else CategoryNaming.fromConfig(cfg) + ) + self.categories = ( + self.naming.expand(cfg["categories"], int(n_slices)) + if n_slices + else list(cfg["categories"]) + ) self.signalFractionForRelevantBins = cfg["signalFractionForRelevantBins"] + self.era_groups = cfg.get("era_groups", {}) + self.bins = [] for era, channel, cat in self.ECC(): bin = self.getBin(era, channel, cat, return_index=False) @@ -127,6 +163,7 @@ def __init__( # print(f"Using hist_bins: {self.hist_binner.hist_bins}") self.input_files = {} + self._merged_away = {} self.shapes = {} self.signal_hists_by_key = {} @@ -141,6 +178,54 @@ def getBin(self, era, channel, category, return_name=True, return_index=True): return index return (index, name) + def mergedAwayIn(self, channel, category): + """Process names absorbed by an active merged process in this bin. + + A process declared with `subprocesses` that name other *datacard* processes + replaces them wherever it applies -- see MinorBkg in the bbWW DL config, + which merges DY/ST/VV in the boosted slices where DY alone has no usable MC + statistics. Suppressing the constituents here is what stops the merge from + double counting, and it means they keep their own configuration unchanged + instead of needing a mirror-image category list to carve the merged bins + back out. + + Inert for every existing configuration: the `subprocesses` lists in + x_hh_bbtautau_run2.yaml name sample-level histograms (WW, WZ, ZZ ...), none + of which is a datacard process, so nothing is ever absorbed there. + """ + key = (channel, category) + if key not in self._merged_away: + absorbed = set() + for p in self.processes.values(): + if not p.subprocesses: + continue + if p.name not in self.channel_processes[channel]: + continue + if not p.appliesToCategory(category): + continue + absorbed |= {s for s in p.subprocesses if s in self.processes} + self._merged_away[key] = absorbed + return self._merged_away[key] + + def processInBin(self, name, channel, category): + """Whether process `name` enters the datacard for this (channel, category).""" + if name not in self.channel_processes[channel]: + return False + if not self.processes[name].appliesToCategory(category): + return False + return name not in self.mergedAwayIn(channel, category) + + def getSubEras(self, era): + """Get sub-eras for a given era. If era is a meta-era, return its sub-eras. + Otherwise return [era].""" + if self.isMetaEra(era): + return self.era_groups[era] + return [era] + + def isMetaEra(self, era): + """Check if era is a meta-era.""" + return era in self.era_groups + def cbCopy(self, param_str, process, era, channel, category): bin_idx, bin_name = self.getBin(era, channel, category) return self.cb.cp().mass([param_str]).process([process]).bin([bin_name]) @@ -148,6 +233,18 @@ def cbCopy(self, param_str, process, era, channel, category): def ECC(self): return itertools.product(self.eras, self.channels, self.categories) + def getCategoryGroups(self): + """{"SR/res2b": ["SR/res2b_dnn0", ...]} -- the slices of one base category. + + Slices of the same base category are the natural unit for a per-category + breakdown: they are one physical selection cut into pieces, not independent + categories. + """ + groups = {} + for cat in self.categories: + groups.setdefault(self.naming.base(cat), []).append(cat) + return groups + def PPECC(self): param_bins = list(self.param_bins.keys()) if not self.model.param_dependent_bkg: @@ -166,50 +263,251 @@ def getInputFile(self, era, model_params): self.input_files[file_name] = file return file_name, self.input_files[file_name] - def getMultiValueLnUnc( - self, unc, unc_name, process, era, channel, category, model_params - ): # , unc_name=None, unc_scale=None) - file_name, file = self.getInputFile(era, model_params) - hist_name = f"{channel}/{category}/{process.hist_name}" - if unc.getUncertaintyForProcess(process.name) != None: - return unc.getUncertaintyForProcess(process.name) - elif process.subprocesses: - unc_value_tot_down = 0.0 - unc_value_tot_up = 0.0 - yield_value_tot = 0.0 - for subp in process.subprocesses: - hist_name = f"{channel}/{category}/{subp}" - subhist = file.Get(hist_name) - # newhist = self.hist_binner.applyBinning(era, channel, category, model_params, subhist) - if subhist == None: - raise RuntimeError( - f"Cannot find histogram {hist_name} in {file.GetName()}" + def _getLnNValue(self, unc, process, proc_name_for_unc, sub_era, channel, category): + if isinstance(unc, MultiValueLnNUncertainty): + return unc.getUncertaintyForProcess( + proc_name_for_unc, sub_era, channel, category + ) + if unc.appliesTo(process, sub_era, channel, category): + return unc.value + return None + + def _applyLnNToHist(self, hist, unc_value, direction): + scaled = hist.Clone() + if isinstance(unc_value, dict): + factor = 1 + unc_value[direction] + elif direction == UncertaintyScale.Up: + factor = 1 + unc_value + else: + factor = 1 - unc_value + scaled.Scale(factor) + scaled.SetDirectory(0) + return scaled + + def _loadBinnedHist(self, file, era, channel, category, model_params, hist_name): + hist = file.Get(hist_name) + if hist is None: + raise RuntimeError(f"Cannot find histogram {hist_name} in {file.GetName()}") + binned = self.hist_binner.applyBinning( + era, channel, category, model_params, hist + ) + binned.SetDirectory(0) + return binned + + def _getSubEraLnNVariedShapes( + self, unc, process, sub_era, channel, category, model_params + ): + file_name, file = self.getInputFile(sub_era, model_params) + hist_names = ( + [(subp, subp) for subp in process.subprocesses] + if process.subprocesses + else [(process.hist_name, process.name)] + ) + up_hist = None + down_hist = None + applies = False + + for hist_name_suffix, proc_name_for_unc in hist_names: + hist = self._loadBinnedHist( + file, + sub_era, + channel, + category, + model_params, + f"{channel}/{category}/{hist_name_suffix}", + ) + unc_value = self._getLnNValue( + unc, process, proc_name_for_unc, sub_era, channel, category + ) + if unc_value is not None: + applies = True + sub_up = self._applyLnNToHist(hist, unc_value, UncertaintyScale.Up) + sub_down = self._applyLnNToHist(hist, unc_value, UncertaintyScale.Down) + else: + sub_up = hist.Clone() + sub_down = hist.Clone() + sub_up.SetDirectory(0) + sub_down.SetDirectory(0) + + if up_hist is None: + up_hist = sub_up + down_hist = sub_down + else: + up_hist.Add(sub_up) + down_hist.Add(sub_down) + + if process.scale != 1: + up_hist.Scale(process.scale) + down_hist.Scale(process.scale) + return up_hist, down_hist, applies + + def getMetaEraLnNShapeUnc(self, unc, process, era, channel, category, model_params): + if not self.isMetaEra(era): + return None + + nominal_shape = self.getShape(process, era, channel, category, model_params) + combined_up = None + combined_down = None + any_applies = False + + for sub_era in self.getSubEras(era): + up, down, applies = self._getSubEraLnNVariedShapes( + unc, process, sub_era, channel, category, model_params + ) + if applies: + any_applies = True + if combined_up is None: + combined_up = up.Clone() + combined_down = down.Clone() + else: + combined_up.Add(up) + combined_down.Add(down) + + if not any_applies: + return None + return nominal_shape, { + UncertaintyScale.Up: combined_up, + UncertaintyScale.Down: combined_down, + } + + def _canIgnoreLnNShape(self, nominal_shape, shapes): + nom_int = nominal_shape.Integral() + if nom_int == 0: + return True + up_frac = (shapes[UncertaintyScale.Up].Integral() - nom_int) / nom_int + down_frac = (shapes[UncertaintyScale.Down].Integral() - nom_int) / nom_int + return abs(up_frac) < self.ignorelnNThr and abs(down_frac) < self.ignorelnNThr + + def _addMetaEraLnNAsShapeUnc( + self, unc_name, proc, param_str, process, era, channel, category, model_params + ): + unc = self.uncertainties[unc_name] + shape_result = self.getMetaEraLnNShapeUnc( + unc, process, era, channel, category, model_params + ) + if shape_result is None: + return False + nominal_shape, shapes = shape_result + if self._canIgnoreLnNShape(nominal_shape, shapes): + print( + f"Ignoring uncertainty {unc_name} for {proc} in {era} {channel} {category}" + ) + return False + + cb_copy = self.cbCopy(param_str, proc, era, channel, category) + cb_copy.AddSyst( + self.cb, + unc_name, + UncertaintyType.shape.name, + ShapeUncertainty(unc_name).valueToMap(), + ) + shape_set = False + + def setShape(syst): + nonlocal shape_set + print(f"Setting unc shape for {syst}") + if shape_set: + raise RuntimeError("Shape already set") + syst.set_shapes( + shapes[UncertaintyScale.Up], + shapes[UncertaintyScale.Down], + nominal_shape, + ) + shape_set = True + + self.cbCopy(param_str, proc, era, channel, category).syst_name( + [unc_name] + ).ForEachSyst(setShape) + return True + + def getCombinedShape( + self, + process, + era, + channel, + category, + model_params, + unc_name=None, + unc_scale=None, + ): + """Combine histograms from multiple sub-eras for a meta-era. + For meta-eras, this sums histograms from constituent sub-eras. + For regular eras, delegates to getShape.""" + if not self.isMetaEra(era): + # Regular era - just get the shape normally + return self.getShape( + process, era, channel, category, model_params, unc_name, unc_scale + ) + + sub_eras = self.getSubEras(era) + + if process.is_asimov_data: + # Build the combined asimov sum from each background's own combined + # (already negative-bin-resolved, with that background's own + # tolerance) shape -- not from raw per-sub-era background shapes + # summed then checked under data_obs's own (untolerant) settings. + # This mirrors how a real era builds asimov data: by summing + # already-resolved per-process shapes, never raw ones. + combined_hist = None + for bkg_proc in self.processes.values(): + if bkg_proc.is_background: + if not self.processInBin(bkg_proc.name, channel, category): + continue + bkg_hist = self.getCombinedShape( + bkg_proc, era, channel, category, model_params ) - axis = subhist.GetXaxis() - yield_subproc = subhist.Integral(1, axis.GetNbins() + 1) - unc_value = unc.getUncertaintyForProcess(subp) - if unc_value != None: - if yield_subproc == 0: + if bkg_hist is None: continue - # print(unc_value) - if isinstance(unc_value, dict): - unc_value_tot_up += ( - unc_value[UncertaintyScale.Up] * yield_subproc - ) - unc_value_tot_down += ( - unc_value[UncertaintyScale.Down] * yield_subproc - ) + if combined_hist is None: + combined_hist = bkg_hist.Clone() else: - unc_value_tot_up += unc_value * yield_subproc - unc_value_tot_down -= unc_value * yield_subproc - yield_value_tot += yield_subproc - if unc_value_tot_up != 0.0 and unc_value_tot_down != 0: - return { - UncertaintyScale.Down: unc_value_tot_down / yield_value_tot, - UncertaintyScale.Up: unc_value_tot_up / yield_value_tot, - } - return None - return None + combined_hist.Add(bkg_hist) + if combined_hist is None: + raise RuntimeError("Cannot create asimov data histogram") + return combined_hist + + # Meta-era: combine histograms from all sub-eras. Negative-bin + # validation is deferred until after summing (below) rather than + # applied per sub-era here -- a sub-era can dip negative on its own + # statistics while the combined shape is fine, and only the combined + # shape is what actually goes into the meta-era datacard. + combined_hist = None + + for sub_era in sub_eras: + sub_hist = self.getShape( + process, + sub_era, + channel, + category, + model_params, + unc_name, + unc_scale, + skip_negative_bin_check=True, + ) + if sub_hist is None: + continue + if combined_hist is None: + combined_hist = sub_hist.Clone() + else: + combined_hist.Add(sub_hist) + + needs_check = combined_hist is not None and not ( + process.is_signal and not (unc_name and unc_scale) + ) + if needs_check: + self.resolveOrRaiseNegativeBins( + combined_hist, + process, + era, + channel, + category, + model_params, + unc_name, + unc_scale, + discovery_eras=sub_eras, + ) + + return combined_hist def getShape( self, @@ -220,9 +518,26 @@ def getShape( model_params, unc_name=None, unc_scale=None, + skip_negative_bin_check=False, ): + # Handle meta-eras by combining sub-era shapes + if self.isMetaEra(era): + return self.getCombinedShape( + process, era, channel, category, model_params, unc_name, unc_scale + ) + file_name, file = self.getInputFile(era, model_params) - key = (file_name, process.name, era, channel, category, unc_name, unc_scale) + key = ( + file_name, + process.name, + era, + channel, + category, + unc_name, + unc_scale, + skip_negative_bin_check, + ) + if key not in self.shapes: if process.is_data and (unc_name is not None or unc_scale is not None): raise RuntimeError("Cannot apply uncertainty to the data process") @@ -230,10 +545,15 @@ def getShape( hist = None for bkg_proc in self.processes.values(): if bkg_proc.is_background: - if bkg_proc.name not in self.channel_processes[channel]: + if not self.processInBin(bkg_proc.name, channel, category): continue bkg_hist = self.getShape( - bkg_proc, era, channel, category, model_params + bkg_proc, + era, + channel, + category, + model_params, + skip_negative_bin_check=skip_negative_bin_check, ) if hist is None: hist = bkg_hist.Clone() @@ -298,65 +618,137 @@ def getShape( ), ) self.signal_hists_by_key.setdefault(key_sig, []).append(hist) - else: - param_str = ( - self.model.paramStr(model_params) if model_params else "*" - ) - key_sig = ( - era, - channel, - category, - ( - param_str - if not self.keep_all_signal_hypothesis_into_single_datacard - else "*" - ), - ) - signals = self.signal_hists_by_key.get(key_sig, []) - relevant_bins = getRelevantBins( + elif not skip_negative_bin_check: + self.resolveOrRaiseNegativeBins( + hist, + process, era, channel, category, - signals, - self.signalFractionForRelevantBins, - ) - solution = resolveNegativeBins( - hist, - relevant_bins=relevant_bins, - allow_zero_integral=process.allow_zero_integral, - allow_negative_bins_within_error=process.allow_negative_bins_within_error, - max_n_sigma_for_negative_bins=process.max_n_sigma_for_negative_bins, - allow_negative_integral=process.allow_negative_integral, + model_params, + unc_name, + unc_scale, ) - - if not solution.accepted: - axis = hist.GetXaxis() - bins_edges = [ - str(axis.GetBinLowEdge(n)) - for n in range(1, axis.GetNbins() + 2) - ] - bin_values = [ - str(hist.GetBinContent(n)) - for n in range(1, axis.GetNbins() + 1) - ] - bin_errors = [ - str(hist.GetBinError(n)) - for n in range(1, axis.GetNbins() + 1) - ] - print(f'bins_edges: [ {", ".join(bins_edges)} ]') - print(f'bin_values: [ {", ".join(bin_values)} ]') - print(f'bin_errors: [ {", ".join(bin_errors)} ]') - raise RuntimeError( - f"Negative bins found in histogram for {channel}/{category}/{process.hist_name}" - + ( - f" (syst {unc_name}{unc_scale})" - if unc_name and unc_scale - else "" - ) - ) self.shapes[key] = hist return self.shapes[key] + def resolveOrRaiseNegativeBins( + self, + hist, + process, + era, + channel, + category, + model_params, + unc_name=None, + unc_scale=None, + discovery_eras=None, + ): + """Validate/rebalance negative bins in-place on `hist`, raising if the + result isn't accepted. `discovery_eras`, when given (meta-era combined + shapes), unions relevant-signal-bin lookups across those real sub-eras + instead of the single `era` -- signal shapes are cached per real + sub-era, never under the meta-era name itself.""" + param_str = self.model.paramStr(model_params) if model_params else "*" + key_param = ( + param_str + if not self.keep_all_signal_hypothesis_into_single_datacard + else "*" + ) + lookup_eras = discovery_eras if discovery_eras else [era] + signals = [] + for lookup_era in lookup_eras: + signals.extend( + self.signal_hists_by_key.get( + (lookup_era, channel, category, key_param), [] + ) + ) + relevant_bins = getRelevantBins( + era, + channel, + category, + signals, + self.signalFractionForRelevantBins, + ) + solution = resolveNegativeBins( + hist, + relevant_bins=relevant_bins, + allow_zero_integral=process.allow_zero_integral, + allow_negative_bins_within_error=process.allow_negative_bins_within_error, + max_n_sigma_for_negative_bins=process.max_n_sigma_for_negative_bins, + allow_negative_integral=process.allow_negative_integral, + ) + + final_integral = sum( + hist.GetBinContent(n) for n in range(1, hist.GetNbinsX() + 1) + ) + is_degenerate = not solution.accepted or final_integral <= 0 + + if is_degenerate and unc_name and unc_scale: + # A shape systematic variation that can't be resolved into a + # valid (positive-integral) histogram -- whether flagged directly + # by resolveNegativeBins, or only zero/negative after its donor + # balancing happened to cancel out the whole shape -- is + # inherently unusable for combine's shape + # morphing (it requires a nonzero norm for every variation). This + # is a low-statistics artifact of the up/down reweighting, not a + # real central-value problem, so fall back to the nominal shape: + # i.e. treat the systematic as having no effect in this bin. + nominal = self.getShape(process, era, channel, category, model_params) + for n in range(1, hist.GetNbinsX() + 1): + hist.SetBinContent(n, nominal.GetBinContent(n)) + hist.SetBinError(n, nominal.GetBinError(n)) + return + + if not solution.accepted: + axis = hist.GetXaxis() + bins_edges = [ + str(axis.GetBinLowEdge(n)) for n in range(1, axis.GetNbins() + 2) + ] + bin_values = [ + str(hist.GetBinContent(n)) for n in range(1, axis.GetNbins() + 1) + ] + bin_errors = [ + str(hist.GetBinError(n)) for n in range(1, axis.GetNbins() + 1) + ] + print(f'bins_edges: [ {", ".join(bins_edges)} ]') + print(f'bin_values: [ {", ".join(bin_values)} ]') + print(f'bin_errors: [ {", ".join(bin_errors)} ]') + raise RuntimeError( + f"Negative bins found in histogram for {channel}/{category}/{process.hist_name}" + + (f" (syst {unc_name}{unc_scale})" if unc_name and unc_scale else "") + ) + + def getSignalProcessForParams(self, model_params): + """Signal Process matching model_params, or None. Used to gate a + param-dependent background on whether the signal hypothesis it's + being evaluated for actually has a shape in a given era/channel/ + category -- backgrounds are looked up per-MX (param_dependent_bkg), + so a category HistRebinTask skipped for that MX has no background + histograms either, not just no signal.""" + for p in self.processes.values(): + if p.is_signal and p.params == model_params: + return p + return None + + def hasNominalShape(self, process, era, channel, category): + """Whether process's nominal shape exists for (era, channel, category), + without raising. Used to skip a signal (and its per-mass background + counterpart) where a specific era+category+mass genuinely has no signal + 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}" + for sub_era in sub_eras: + _, file = self.getInputFile(sub_era, process.params) + obj = file.Get(hist_name) + # 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 obj:` (not `is not None`) is the correct null check. + if obj and obj.InheritsFrom("TH1"): + return True + return False + def addProcess(self, proc, era, channel, category): bin_idx, bin_name = self.getBin(era, channel, category) process = self.processes[proc] @@ -386,7 +778,6 @@ def add(model_params, param_str, process_name): def setShape(p): nonlocal shape_set - print(f"Setting shape for {p}") if shape_set: raise RuntimeError("Shape already set") p.set_shape(shape, True) @@ -399,6 +790,11 @@ def setShape(p): cb_copy.ForEachProc(setShape) if process.is_signal: + if not self.hasNominalShape(process, era, channel, category): + print( + f"Skipping {process.name} in {era}/{channel}/{category}: no signal shape found" + ) + return model_params = process.params param_str = self.model.paramStr(model_params) if self.keep_all_signal_hypothesis_into_single_datacard: @@ -413,15 +809,25 @@ def setShape(p): self.base_of[actual_proc_name] = process.name elif self.model.param_dependent_bkg: + # One copy of this process per distinct signal *parameter point*, not per + # signal process: several signal processes (e.g. the bbWW and bbtautau + # decay modes) share the same mass grid, and adding the copy once per + # process would set the same shape twice ("Shape already set"). + seen_params = set() for signal_proc in self.processes.values(): if not signal_proc.is_signal: continue + if not self.hasNominalShape(signal_proc, era, channel, category): + continue model_params = signal_proc.params param_str = ( self.model.paramStr(model_params) if not self.keep_all_signal_hypothesis_into_single_datacard else "*" ) + if param_str in seen_params: + continue + seen_params.add(param_str) add(model_params, param_str, proc) self.param_of[(param_str, proc)] = model_params self.base_of[proc] = proc @@ -433,16 +839,45 @@ def addUncertainty(self, unc_name): isMVLnUnc = isinstance(unc, MultiValueLnNUncertainty) for proc, param_str, era, channel, category in self.PPECC(): - if proc not in self.channel_processes[channel]: + if not self.processInBin(proc, channel, category): continue process = self.processes[proc] if process.is_data: continue model_params = self.param_bins.get(param_str, None) + if not process.hasCompatibleModelParams( + model_params, self.model.param_dependent_bkg + ): + continue + if process.is_signal: + if not self.hasNominalShape(process, era, channel, category): + continue + elif self.model.param_dependent_bkg and model_params is not None: + signal_proc = self.getSignalProcessForParams(model_params) + if signal_proc is not None and not self.hasNominalShape( + signal_proc, era, channel, category + ): + continue + + if self.isMetaEra(era) and isinstance( + unc, (LnNUncertainty, MultiValueLnNUncertainty) + ): + self._addMetaEraLnNAsShapeUnc( + unc_name, + proc, + param_str, + process, + era, + channel, + category, + model_params, + ) + continue + if isMVLnUnc: - unc_value = self.getMultiValueLnUnc( - unc, unc_name, process, era, channel, category, model_params - ) # , unc_name=None, unc_scale=None + unc_value = unc.getUncertaintyForProcess( + process.name, era, channel, category + ) uncApplies = ( unc_value != None @@ -451,10 +886,6 @@ def addUncertainty(self, unc_name): ) if not uncApplies: continue - if not process.hasCompatibleModelParams( - model_params, self.model.param_dependent_bkg - ): - continue nominal_shape = None shapes = {} @@ -521,10 +952,27 @@ def setShape(syst): process = self.processes[base_name] if process.is_data: continue + if not process.hasCompatibleModelParams( + params, self.model.param_dependent_bkg + ): + continue + + if self.isMetaEra(era) and unc.type == UncertaintyType.lnN: + self._addMetaEraLnNAsShapeUnc( + unc_name, + proc_name, + param_str, + process, + era, + channel, + category, + params, + ) + continue if isMVLnUnc: - unc_value = self.getMultiValueLnUnc( - unc, unc_name, process, era, channel, category, params + unc_value = unc.getUncertaintyForProcess( + process.name, era, channel, category ) uncApplies = ( (unc_value is not None) @@ -533,10 +981,6 @@ def setShape(syst): ) if not uncApplies: continue - if not process.hasCompatibleModelParams( - params, self.model.param_dependent_bkg - ): - continue nominal_shape = None shapes = {} @@ -606,13 +1050,26 @@ def writeDatacards(self, output): return background_names = [n for n, p in self.processes.items() if p.is_background] + + # Group the signal processes by parameter point. Several signal processes can + # share a mass (e.g. the bbWW and bbtautau decay modes of the same resonance); + # they must go into the *same* datacard so the fit scales them with a common + # signal strength, rather than yielding a separate limit per decay mode. + signals_by_param = {} for proc_name, process in self.processes.items(): if not process.is_signal: continue - processes = [proc_name] + background_names - param_list = [self.model.paramStr(process.params)] + key = self.model.paramStr(process.params) + signals_by_param.setdefault(key, []).append(proc_name) + + for param_str, signal_names in signals_by_param.items(): + processes = list(signal_names) + background_names + param_list = [param_str] if not self.model.param_dependent_bkg: param_list.append("*") + # Named after the primary (first configured) signal, so a single-signal + # config keeps exactly the file names it produced before. + proc_name = signal_names[0] dc_file = os.path.join(output, f"datacard_{proc_name}.txt") shape_file = os.path.join(output, f"{proc_name}.root") @@ -626,6 +1083,34 @@ def writeDatacards(self, output): param_list ).process(processes).WriteDatacard(tmp_dc_file, tmp_shape_file) + # Same breakdown by base category (all its DNN slices, all channels), + # for per-category limits alongside the per-channel ones. + for base_cat, slice_cats in self.getCategoryGroups().items(): + bin_names = [ + self.getBin(subera, subchannel, cat, return_index=False) + for subchannel in self.channels + for cat in slice_cats + ] + selected = ( + self.cb.cp() + .era([subera]) + .bin(bin_names) + .mass(param_list) + .process(processes) + ) + # A base category can be absent for a given mass hypothesis (e.g. + # boosted at low MX, where HistRebinTask found too little signal + # to slice it) -- there is no card to write then. + if len(selected.bin_set()) == 0: + continue + cat_dir = os.path.join( + output, subera, "categories", base_cat.replace("/", "_") + ) + os.makedirs(cat_dir, exist_ok=True) + selected.WriteDatacard( + os.path.join(cat_dir, f"datacard_{proc_name}.txt"), shape_file + ) + self.cb.cp().mass(param_list).process(processes).WriteDatacard( dc_file, shape_file ) @@ -634,13 +1119,13 @@ def createDatacards(self, output, verbose=1): try: for era, channel, category in self.ECC(): for name, p in self.processes.items(): - if name not in self.channel_processes[channel]: + if not self.processInBin(name, channel, category): continue if p.is_signal: self.addProcess(name, era, channel, category) for era, channel, category in self.ECC(): for name, p in self.processes.items(): - if name not in self.channel_processes[channel]: + if not self.processInBin(name, channel, category): continue if not p.is_signal: self.addProcess(name, era, channel, category) diff --git a/dc_make/plot_rebinned.py b/dc_make/plot_rebinned.py new file mode 100644 index 0000000..56719d1 --- /dev/null +++ b/dc_make/plot_rebinned.py @@ -0,0 +1,566 @@ +"""Stacked plots of the rebinned 1D shapes that actually enter the datacards. + +HistPlotTask plots the pre-rebinning histograms straight out of HistMergerTask, one +era at a time. That is the wrong object to look at once HistRebinTask has sliced the +2D (DNN x HME) plane into per-DNN-slice 1D mass shapes: what the fit sees is the +*rebinned* shape, summed over the meta-era's sub-eras exactly as maker.py's +getCombinedShape does at datacard-build time. + +This script reproduces HistPlotTask's look (same PlotKit Plotter, same +config/plot/*.yaml style files, same signal-overlay convention) but takes the +rebinned files as input and sums the sub-eras first, so every panel corresponds +one-to-one with a datacard bin block. + +One file per base category and mass, laid out as a grid: DNN slices across, +lepton channels down. The fit sees these bins together, so they are easier to judge +together -- a slice that looks reasonable in eMu and pathological in eE is obvious +side by side and invisible in separate files. +""" + +import os +import sys +import yaml + +if __name__ == "__main__": + file_dir = os.path.dirname(os.path.abspath(__file__)) + pkg_dir = os.path.dirname(file_dir) + base_dir = os.path.dirname(pkg_dir) + pkg_dir_name = os.path.split(pkg_dir)[1] + if base_dir not in sys.path: + sys.path.append(base_dir) + __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 + +ROOT = importROOT() + +# Fallback colours, used only when a process is missing from processes.yaml so that +# a config gap degrades to an ugly plot rather than a crash. +FALLBACK_COLORS = { + "TT": "kAzure+1", + "DY": "kOrange-3", + "ST": "kGreen+1", + "VV": "kViolet-4", + "W": "kRed-7", + "SingleHiggs": "kGray+1", +} +FALLBACK_SIGNAL_COLOR = "kRed" + + +def load_config(config_path, n_dnn_slices=None, category_pattern=None): + """Like hist_rebin_2d.load_config, but expands the base categories into the + *sliced* names -- those are the directories that exist in the rebinned files and + the ones that map onto datacard bins.""" + with open(config_path, "r") as f: + cfg = yaml.safe_load(f) + + # HistRebinTask's resolved value is authoritative when the caller passes it; a + # command-line override would otherwise leave this deriving a different slice count + # than the rebinned files contain. No slice count anywhere (or an explicit 0) means + # the input is already binned and its categories are the ones listed -- same rule as + # DatacardMaker, so the panels and the datacard bins cannot disagree. + n_slices = n_dnn_slices + if n_slices is None: + n_slices = (cfg.get("binning") or {}).get("n_dnn_slices") + # Same rule for the pattern naming those slices: the caller's resolved value wins, + # otherwise the configuration's own. + naming = ( + CategoryNaming(category_pattern) + if category_pattern + else CategoryNaming.fromConfig(cfg) + ) + + # multiple signal processes are allowed; each is drawn as its own stack entry + signal_hist_names = [] + mass_values = None + background_entries = [] + for entry in cfg["processes"]: + if type(entry) == str: + background_entries.append((entry, [])) + continue + if entry.get("is_data", False): + continue + hist_name = entry.get("hist_name", entry["process"]) + if entry.get("is_signal", False): + signal_hist_names.append(hist_name) + if mass_values is None: + mass_values = entry["param_values"] + else: + background_entries.append((hist_name, entry.get("channels", []))) + + if not signal_hist_names: + raise RuntimeError("No signal process found in config") + + return { + "model": Model.fromConfig(cfg["model"]), + "channels": cfg["channels"], + "naming": naming, + "categories": ( + naming.expand(cfg["categories"], int(n_slices)) + if n_slices + else list(cfg["categories"]) + ), + "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 get_process_styles(ana_path, period, version, process_names): + """Plot label + colour per process, taken from the analysis processes.yaml so + the rebinned plots use the same colours as HistPlotTask.""" + styles = {} + try: + import FLAF.Common.Setup as Setup + + setup = Setup.Setup(ana_path, period, version) + parent = setup.parent_processes + for name in process_names: + if name in parent: + proc = parent[name] + styles[name] = (proc.get("name", name), proc["color"]) + except Exception as e: + print( + f"Warning: could not load process styles from Setup ({e}); using fallbacks" + ) + + for name in process_names: + if name not in styles or styles[name][1] == "kBlack": + # kBlack is what Setup assigns to meta-process members it was not asked + # to plot (every signal mass but the three in `to_plot`), which would be + # invisible against the axis -- give the overlaid signal its own colour. + is_signal = name not in FALLBACK_COLORS + color = FALLBACK_SIGNAL_COLOR if is_signal else FALLBACK_COLORS[name] + label = styles[name][0] if name in styles else name + styles[name] = (label, color) + return styles + + +def build_hist_cfg(x_title, y_title, log_y): + """PlotKit reads binning/axis metadata out of histograms.yaml keyed by variable + name. The rebinned axis is a bin index with no config entry (and needs none -- + the binning is already final), so hand it a synthetic one.""" + return { + "rebinned": { + "x_title": x_title, + "y_title": y_title, + "use_log_y": log_y, + # Enough headroom for the legend without leaving the stack squashed into + # the bottom decade (HistPlotter's 2000x is tuned for a full-page canvas). + "max_y_sf": 100.0 if log_y else 1.5, + "divide_by_bin_width": False, + } + } + + +class PanelBackend: + """Draws a PlotKit StackSpec into an axes we already own. + + PlotKit's MplhepBackend makes its own figure in _new_figure() and saves it in + render_stacked(). Swapping just those two things out lets every DNN slice share + one canvas while the drawing code -- and therefore the styling -- stays exactly + the one HistPlotTask uses. + """ + + def __init__(self, ax, draw_cms_label, draw_legend): + from FLAF.PlotKit.backends import MplhepBackend + + self._ax = ax + self._draw_cms_label = draw_cms_label + self._draw_legend = draw_legend + + outer = self + + class _Panel(MplhepBackend): + def _new_figure(self, spec): + import mplhep as hep + + class _NoSave: + def savefig(self, *args, **kwargs): + pass + + class _NoClose: + @staticmethod + def close(*args, **kwargs): + pass + + return _NoClose(), hep, _NoSave(), outer._ax, None + + def _cms_label(self, hep, ax, spec): + # Only the leftmost panel carries "CMS Simulation" + the lumi. + if outer._draw_cms_label: + super()._cms_label(hep, ax, spec) + + self._impl = _Panel() + + def render(self, spec): + self._impl.render_stacked(spec, None) + if not self._draw_legend: + legend = self._ax.get_legend() + if legend is not None: + legend.remove() + + +def plot_channel_slice_grid( + rows, + plotter_cfg, + out_path, + want_data, + signal_scale, + title, +): + """One figure per base category: DNN slices across, lepton channels down. + + `rows` is [(channel, [panel or None, ...]), ...], each list holding one entry per + DNN slice in slice order. A missing slice keeps its column slot as None instead of + shifting the row left -- dnn0 has to sit above dnn0 across channels, since making + the columns comparable is the whole reason for sharing a canvas. + + Each panel keeps its own legend and channel label: the signal entry carries that + panel's yield, so it is per-panel information rather than boilerplate. + """ + import matplotlib + + if matplotlib.get_backend().lower() not in ("agg", "pdf", "svg", "ps"): + matplotlib.use("Agg") + import matplotlib.pyplot as plt + import mplhep as hep + + from FLAF.PlotKit.plotters.stacked import StackedPlotter + from FLAF.PlotKit.rootcompat import tlatex_to_mpl as _to_mathtext + + plt.style.use(hep.style.CMS) + n_rows = len(rows) + n_cols = max(len(panels) for _, panels in rows) + fig, axes = plt.subplots( + n_rows, n_cols, figsize=(7.0 * n_cols, 7.5 * n_rows), squeeze=False + ) + + lumi_text = plotter_cfg.text_box("lumi_text").get("text", "") + drawn_cms_label = False + for r, (channel, panels) in enumerate(rows): + for c in range(n_cols): + ax = axes[r][c] + panel = panels[c] if c < len(panels) else None + if panel is None: + # HistRebinTask produced no shapes here (e.g. a channel whose only + # backgrounds are excluded). Blank the cell rather than closing the gap. + ax.set_axis_off() + continue + slice_label, hists, custom = panel + backend = PanelBackend( + ax, draw_cms_label=not drawn_cms_label, draw_legend=True + ) + drawn_cms_label = True + sp = StackedPlotter(plotter_cfg, backend._impl) + spec = sp.build_spec("rebinned", hists, want_data, custom, signal_scale) + # The lumi belongs to the whole canvas, and inside a single narrow panel it + # collides with "CMS Simulation" -- draw it once at figure level instead. + spec.lumi_text = "" + backend.render(spec) + # Right-aligned: the CMS label the first panel draws is left-aligned above + # the axes, and the lumi that would otherwise sit on the right was just + # blanked, so this is the only thing on that side. + ax.set_title(slice_label, loc="right", fontsize=16) + if c > 0: + ax.set_ylabel("") + + fig.suptitle(title, y=1.0, fontsize=20) + if lumi_text: + fig.text( + 1.0, 1.0, _to_mathtext(lumi_text), ha="right", va="bottom", fontsize=18 + ) + fig.tight_layout() + os.makedirs(os.path.dirname(out_path), exist_ok=True) + fig.savefig(out_path, bbox_inches="tight") + plt.close(fig) + + +def sum_over_eras(files, hist_path): + total = None + for f in files: + h = f.Get(hist_path) + if not h: + continue + h.SetDirectory(0) + if total is None: + total = h.Clone() + total.SetDirectory(0) + else: + total.Add(h) + return total + + +def build_category_panel(files, channel, category, cfg, mass, styles, era_label_cfg): + """(hists, custom) for one DNN slice, or None if the slice has no shapes. + + A slice is missing whenever HistRebinTask skipped the whole base category for + lack of discovery signal (e.g. boosted at MX=300).""" + signal_keys = [ + applyParameters(pattern, {cfg["signal_param_name"]: mass}) + for pattern in cfg["signal_hist_name_patterns"] + ] + prefix = f"{channel}/{category}/" + + hists = {} + for bkg_key, allowed_channels in cfg["background_entries"]: + if allowed_channels and channel not in allowed_channels: + continue + h = sum_over_eras(files, prefix + bkg_key) + if h is None: + continue + label, color = styles.get(bkg_key, (bkg_key, "kGray")) + hists[bkg_key] = (h, label, color, "backgrounds") + + if not hists: + return None + + for signal_key in signal_keys: + sig = sum_over_eras(files, prefix + signal_key) + if sig is not None and sig.Integral() > 0: + label, color = styles.get(signal_key, (signal_key, FALLBACK_SIGNAL_COLOR)) + hists[signal_key] = (sig, label, color, "signals") + + # category is "SR/res2b_dnn0"; the region prefix is already in the label boxes + custom = { + "cat_text": category.split("/")[-1], + "ch_text": era_label_cfg.get("channel_text", {}).get(channel, channel), + "customreg_text": era_label_cfg.get("customregion_text", {}).get( + category.split("/")[0], "" + ), + # PlotKit already draws "Simulation" itself when want_data is False, so the + # scope must not repeat it -- HistPlotter passes "CMS simulation" here and + # ends up rendering "CMS Simulation simulation". Leaving the scope empty + # defers to scope_text in the era's plot config. + "datasim_text": "CMS ", + "scope_text": "", + } + return hists, custom + + +def slice_title(files, channel, category): + """Panel title: the sliced category name plus the DNN selection it stands for. + + The selection is the slice directory's own title, set by hist_rebin_2d.py. Every era + holds the same edges -- they are discovered once from the summed discovery eras -- so + the first file carrying the slice wins. Input that was never sliced (a configuration + with no `binning:` block) has no title of its own and simply carries no label. + """ + name = category.split("/")[-1] + for f in files: + slice_dir = f.Get(f"{channel}/{category}") + if not slice_dir: + continue + # ROOT defaults a directory's title to its own name, so only a title that + # differs from it is a selection hist_rebin_2d.py actually wrote. + label = slice_dir.GetTitle() + return f"{name} {label}" if label and label != name else name + return name + + +def group_categories_by_base(categories, naming): + """{"SR/res2b": ["SR/res2b_dnn0", ...]} in slice order -- the slices of one base + category are what share a canvas.""" + groups = {} + for cat in categories: + groups.setdefault(naming.base(cat), []).append(cat) + return groups + + +def run( + input_dir, + output_dir, + config_path, + eras, + era_label, + ana_path, + version, + signal_scale, + log_y, + masses=None, + n_dnn_slices=None, + category_pattern=None, +): + from FLAF.PlotKit.config import PlotConfig + + cfg = load_config(config_path, n_dnn_slices, category_pattern) + model = cfg["model"] + mass_values = masses if masses else cfg["mass_values"] + + page_cfg = os.path.join(ana_path, "config", "plot", "cms_stacked.yaml") + page_cfg_custom = os.path.join(ana_path, "config", "plot", f"{era_label}.yaml") + if not os.path.exists(page_cfg_custom): + print(f"Warning: {page_cfg_custom} not found, falling back to {eras[0]}") + page_cfg_custom = os.path.join(ana_path, "config", "plot", f"{eras[0]}.yaml") + with open(page_cfg_custom, "r") as f: + era_label_cfg = yaml.safe_load(f) + + process_names = [name for name, _ in cfg["background_entries"]] + process_names += [ + applyParameters(pattern, {cfg["signal_param_name"]: m}) + for pattern in cfg["signal_hist_name_patterns"] + for m in mass_values + ] + styles = get_process_styles(ana_path, eras[0], version, process_names) + + hist_cfg = build_hist_cfg( + x_title="m_{X}^{HME} (GeV)", y_title="Events", log_y=log_y + ) + plotter_cfg = PlotConfig(page_cfg, page_cfg_custom, hist_cfg) + groups = group_categories_by_base(cfg["categories"], cfg["naming"]) + + n_made, n_skipped = 0, 0 + for mass in mass_values: + rel = model.getInputFileName(eras[0], {cfg["signal_param_name"]: mass}) + # getInputFileName() prefixes the era; input_dir holds one such tree per era + var_rel = rel[len(eras[0]) + 1 :] if rel.startswith(eras[0] + "/") else rel + + files = [] + for era in eras: + path = os.path.join(input_dir, era, var_rel) + f = ROOT.TFile.Open(path, "READ") + if f is None or f.IsZombie(): + raise RuntimeError(f"Cannot open rebinned file {path}") + files.append(f) + + for base_cat, slice_cats in groups.items(): + # All lepton channels of a base category share one canvas, so the DNN + # slices can be compared across channels at a glance instead of by + # flipping between files. + rows = [] + for channel in cfg["channels"]: + panels = [] + for category in slice_cats: + built = build_category_panel( + files, channel, category, cfg, mass, styles, era_label_cfg + ) + if built is None: + panels.append(None) + continue + hists, custom = built + panels.append( + ( + slice_title(files, channel, category), + hists, + custom, + ) + ) + rows.append((channel, panels)) + + if not any(panel for _, panels in rows for panel in panels): + # HistRebinTask skipped this base category outright (e.g. boosted + # at low MX): no slices exist in any channel, so there is nothing to draw. + n_skipped += 1 + continue + + out_path = os.path.join(output_dir, f"m{mass}", f"{base_cat}.pdf") + plot_channel_slice_grid( + rows, + plotter_cfg, + out_path, + want_data=False, + signal_scale=signal_scale, + title=f"{base_cat} $M_X$ = {mass} GeV", + ) + n_made += 1 + + for f in files: + f.Close() + + print( + f"plot_rebinned: {n_made} canvases written to {output_dir} ({n_skipped} empty categories skipped)" + ) + + +if __name__ == "__main__": + import argparse + + parser = argparse.ArgumentParser( + description="Stacked plots of the merged, rebinned shapes entering the datacards." + ) + parser.add_argument( + "--input", + required=True, + type=str, + help="base directory holding //.root for each era in --eras", + ) + parser.add_argument( + "--output", required=True, type=str, help="output directory for the pdfs" + ) + parser.add_argument( + "--config", required=True, type=str, help="datacard configuration yaml" + ) + parser.add_argument( + "--eras", + required=True, + type=str, + help="comma-separated eras to sum before plotting", + ) + parser.add_argument( + "--era-label", + required=True, + type=str, + help="name used for config/plot/