diff --git a/README.md b/README.md index 9d89ce8..1af06a9 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,61 @@ # 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. **Neither is +part of the datacard chain**: both are offline pre-steps that produce input the chain then +consumes, so `CreateDatacardsTask` reads whatever shapes it is pointed at and does no +binning of its own. + +**1. Shape-driven 2D→1D binning — `bin_opt_2d/rebin_2d.py`.** Derives bin edges from the +shapes themselves (no fits, no limits): it cuts the x axis of a 2D input into slices that +become datacard categories, and rebins y into the bins of each slice. It writes the +rebinned shapes as `///.root`, plus the +`binning.json` recording the edges it chose, beside them. + +It runs as the datacard configuration's **preprocessing step**, not as part of this +repository's logic. `PreprocessShapesTask` runs whatever `preprocess:` names, supplying +`--input`, `--output`, `--era` and `--config` (the datacard configuration); nothing here +knows what the step does, and a +configuration that declares no `preprocess:` block skips the task entirely and reads the +merged histograms as they are: + +```yaml +preprocess: + script: StatInference/bin_opt_2d/rebin_2d.py + args: [--binning-config, config/Datacards/binning_2d.yaml] +``` + +The knobs live with the analysis (`config/Datacards/binning_2d.yaml`) because they are +analysis configuration; the derived `binning.json` lives with the shapes it produced, +because it is a product of the run rather than an input to it. + +Which era is being produced decides how it is derived. A plain era is binned on its own +statistics, for a standalone limit; an era that is a key of `era_groups:` is binned on its +members' summed statistics, and those edges are then applied to each member separately -- +kept in their own files, because the datacard step sums them and a per-era lnN can only be +built by scaling a sub-era's own shape. + +Each base category (`SR/res2b`) becomes per-slice categories named by the `category_pattern` +knob, e.g. `{base_category}_dnn{slice_idx}`. The pattern is the analysis's choice — nothing +here assumes the sliced axis is a DNN score. `common/tools.py:CategoryNaming` both writes +those names and parses them back from that one pattern. + +**2. Limit-driven 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`. Every module is a +script driven by `bin_opt/bin_optimization.yaml`, and it exposes no importable API. + +Which one an analysis is on is visible in its configuration: the `bin_opt` path sets +`hist_bins` (as `config/x_hh_bbtautau_run2.yaml` does), while an analysis reading rebinned +shapes leaves `hist_bins` unset and instead lists the sliced `categories:` and the +`category_pattern` that names them (as HH→bbWW's `config/Datacards/x_hh_bbww_DL_run3.yaml` +does). An analysis using neither simply lists the categories its input already has. + +`categories:` is always taken verbatim — it is the set of directories that exist in the +input shapes, and nothing in `dc_make` derives or expands it. + ## 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/bin_opt_2d/rebin_2d.py b/bin_opt_2d/rebin_2d.py new file mode 100644 index 0000000..7bb8f2f --- /dev/null +++ b/bin_opt_2d/rebin_2d.py @@ -0,0 +1,1218 @@ +import array +import json +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() + + +# Written inside each era's own output directory, beside the shapes it describes, and +# read back by --binning to replay them. +BINNING_JSON = "binning.json" + +# The knobs that decide the binning, with the values a configuration gets when it does +# not say otherwise. They belong to an analysis rather than to this file -- what a slice +# needs to be worth keeping depends on the sample sizes and the selection, and a test +# configuration wants them at zero -- so every production should state its own in a +# binning yaml. See bin_opt_2d/binning.yaml for the annotated HH->bbWW set. +BINNING_DEFAULTS = { + "n_slices": 4, + "category_pattern": None, # None -> CategoryNaming's own neutral default + "slice_var": "x", + "max_bins_per_slice": 10, + "min_slice_bkg_sum": 1.0, + "min_slice_bkg_neff": 4.0, + "min_slice_bkg_each": 0.01, + "min_slice_bkg_each_neff": 0.0, + "min_bin_bkg_each": 0.01, + "min_bin_bkg_neff": 4.0, + "bkg_per_bin": 5.0, + "min_bkg_frac": 0.05, + "min_signal": 0.5, + "significance_mode": "asimov", +} + + +def load_binning_config(path, overrides=None): + """Merge the binning yaml over the defaults, then command-line overrides over that. + + Unknown keys raise rather than being ignored: a misspelled floor that silently does + nothing is the failure mode worth spending an exception on. + """ + declared = {} + if path: + with open(path, "r") as f: + declared = yaml.safe_load(f) or {} + unknown = set(declared) - set(BINNING_DEFAULTS) + if unknown: + raise RuntimeError( + f"{path}: unknown binning key(s) {sorted(unknown)}. " + f"Known keys: {sorted(BINNING_DEFAULTS)}." + ) + knobs = dict(BINNING_DEFAULTS) + knobs.update(declared) + for key, value in (overrides or {}).items(): + if value is not None: + knobs[key] = value + return knobs + + +def lookup_frozen(frozen_binning, era, mass, channel, category): + """The recorded binning for one channel/category, or None if it was skipped then. + + A category absent from the record was skipped by the run that wrote it (too little + signal, or a background missing from a discovery era), so it is skipped again rather + than quietly re-optimised -- a replay that rebinned some categories and froze others + would be neither the old binning nor a new one. + """ + return ( + frozen_binning.get("binning", {}) + .get(era, {}) + .get(str(mass), {}) + .get(channel, {}) + .get(category) + ) + + +def load_config(config_path): + with open(config_path, "r") as f: + cfg = yaml.safe_load(f) + model = Model.fromConfig(cfg["model"]) + channels = cfg["channels"] + # Taken verbatim here; run() reduces them to the base categories the 2D input is + # actually keyed by, once it has the pattern to do it with. + categories = list(cfg["categories"]) + + # Several processes may carry is_signal (e.g. the bbWW(2l) and bbtautau decay + # modes of the same resonance, both scaled by the same signal strength). They + # are summed to form the discovery signal that steers the slice + # boundaries, so the binning is optimised for the total signal in the fit. + signal_hist_names = [] + mass_values = None + background_entries = [] + for entry in cfg["processes"]: + if type(entry) == str: + background_entries.append((entry, entry, [])) + continue + if entry.get("is_data", False): + continue + base_name = entry["process"] + hist_name = entry.get("hist_name", base_name) + if entry.get("is_signal", False): + signal_hist_names.append(hist_name) + if mass_values is None: + mass_values = entry["param_values"] + elif list(entry["param_values"]) != list(mass_values): + raise RuntimeError( + f"Signal {hist_name} has param_values {entry['param_values']}, " + f"which differ from {mass_values}; every signal must be defined " + "at the same mass points" + ) + else: + background_entries.append((base_name, hist_name, entry.get("channels", []))) + + if not signal_hist_names: + raise RuntimeError("No signal process found in config") + + return { + "model": model, + "channels": channels, + "categories": categories, + "era_groups": cfg.get("era_groups", {}), + "signal_hist_name_patterns": signal_hist_names, + "signal_param_name": extractParameters(signal_hist_names[0])[0], + "mass_values": mass_values, + "background_entries": background_entries, + } + + +def open_input_file(input_dir, model, era, mass, param_name): + file_name = model.getInputFileName(era, {param_name: mass}) + full_path = os.path.join(input_dir, file_name) + f = ROOT.TFile.Open(full_path, "READ") + if f is None or f.IsZombie(): + raise RuntimeError(f"Cannot open file {full_path}") + return f + + +def _detach(h): + """Detach a histogram from its TFile and hand ownership to Python. + + SetDirectory(0) alone makes the histogram survive the file's close, but + leaves it owned by nobody -- so every 2D histogram read here (one per + systematic variation, per category, per mass) leaked for the lifetime of the + process. Across the ten-mass loop that grew without bound and got the job + SIGKILLed partway through the final mass, leaving a truncated ROOT file. + SetOwnership makes the object die with its last Python reference. + """ + h.SetDirectory(0) + ROOT.SetOwnership(h, True) + return h + + +def get_hist(f, path): + h = f.Get(path) + # TFile.Get() on a fully-missing nested path can return a PyROOT wrapper + # around a null C++ pointer, which is not `is None` but is falsy. + if not h: + return None + return _detach(h) + + +def sum_hists(hists): + total = None + for h in hists: + if h is None: + continue + if total is None: + total = _detach(h.Clone()) + else: + total.Add(h) + return total + + +def _integral(hist, lo, hi): + return ( + hist.Integral(lo, hi, 0, -1) + if hist.GetDimension() == 2 + else hist.Integral(lo, hi) + ) + + +def _integral_and_error(hist, lo, hi): + """(yield, MC statistical error) over [lo, hi].""" + err = array.array("d", [0.0]) + if hist.GetDimension() == 2: + value = hist.IntegralAndError(lo, hi, 0, -1, err) + else: + value = hist.IntegralAndError(lo, hi, err) + return value, err[0] + + +def _bkg_yields(bkg_hists_by_name, lo, hi): + """For each background, the summed yield across all discovery eras (combined + statistics). Individual eras are not required to individually clear the + threshold -- use allow_negative_bins_within_error (maker.py, per process/ + category) for backgrounds/categories where a specific era can land negative + in a bin that's fine once combined.""" + return { + name: sum(_integral(h, lo, hi) for h in hists) + for name, hists in bkg_hists_by_name.items() + } + + +def _bkg_errors(bkg_hists_by_name, lo, hi): + """Per background, the MC statistical error on the summed yield (eras added in + quadrature).""" + return { + name: math.sqrt(sum(_integral_and_error(h, lo, hi)[1] ** 2 for h in hists)) + for name, hists in bkg_hists_by_name.items() + } + + +def _total_bkg_error(bkg_hists_by_name, lo, hi): + errors = _bkg_errors(bkg_hists_by_name, lo, hi) + return math.sqrt(sum(e**2 for e in errors.values())) + + +def effective_entries(value, error): + """(yield / error)^2 -- the unweighted MC event count a weighted yield is worth. + + A DY slice holding a couple of very-high-weight aMC@NLO events can carry a + sizeable yield with an effective count far below 1, i.e. a background estimate + that is statistically compatible with almost anything. + """ + if error <= 0: + return float("inf") if value > 0 else 0.0 + return (value / error) ** 2 + + +def asimov_significance(s, b, b_err=0.0): + """Median discovery significance for a counting experiment (Cowan et al., + arXiv:1007.1727), with the background uncertainty folded in when given. + + Reduces to sqrt(2*((s+b)*ln(1+s/b) - s)) for b_err = 0. Unlike S/sqrt(B) this + stays valid when b is O(1), which is exactly the regime the high-mass boosted + slices live in -- there S/sqrt(B) reports significances of ~30 on ~1.5 + background events, and drives the slice boundary on that basis. + """ + if b <= 0 or s <= 0: + return 0.0 + if b_err <= 0: + return math.sqrt(max(2.0 * ((s + b) * math.log1p(s / b) - s), 0.0)) + var = b_err**2 + term1 = (s + b) * math.log(((s + b) * (b + var)) / (b * b + (s + b) * var)) + term2 = (b * b / var) * math.log1p(var * s / (b * (b + var))) + return math.sqrt(max(2.0 * (term1 - term2), 0.0)) + + +def significance(s, b, b_err=0.0, mode="sb"): + """Figure of merit steering the slice boundaries. + + mode="sb": S / sqrt(B + sigma_B^2). Folding sigma_B into plain S/sqrt(B) + stops a downward fluctuation of a statistics-starved background + (DY in the b-tagged muMu slices, effective MC count <1) from + inflating the apparent significance. Still assumes the Gaussian + regime, which breaks down for B of order a few. + mode="asimov": the Poisson-correct Asimov significance, valid at low B. + """ + if b is None or b <= 0: + return 0.0 + if mode == "asimov": + return asimov_significance(s, b, b_err) + denominator = b + b_err**2 + if denominator <= 0: + return 0.0 + return s / (denominator**0.5) + + +def _slice_passes( + yields, + min_sum, + total_error=None, + min_neff=0.0, + errors=None, + min_each=0.0, + min_proc_neff=0.0, + exempt=(), +): + """Slice validity: the summed background must clear min_sum, and be known + to better than min_neff effective MC entries -- a window whose background is + statistically undetermined must not be selectable at all. + + The per-process arms (min_each, min_proc_neff) are what _bin_passes already + does on the mass axis, applied here as well. Testing only the sum hides an + individual background that has fluctuated negative behind a large, well + measured neighbour: muMu/SR/res2b at m500 selected a slice holding + DY = -5.87 +- 8.92 (N_eff 0.43) because the summed background there was + +39.9 +- 9.0 (N_eff 22), clearing both summed tests comfortably. That slice + cannot be turned into a datacard at all -- a negative DY integral is rejected + outright by resolveNegativeBins -- so the binning has to not select it in the + first place. Worse, the significance being maximised is S/sqrt(B + sigma_B^2), + which a downward fluctuation *increases* by lowering B, so such a window is + mildly preferred rather than merely tolerated. + + `exempt` comes from minor_backgrounds() judged over the whole sliced-axis range of the + category, never over the candidate window: a background that is negligible in + this category must not be able to veto every boundary, and judging it inside + the window under test is circular (a background that is exactly zero there is + trivially below any fraction of the total). + """ + total = sum(yields.values()) + if total <= min_sum: + return False + if min_neff > 0 and total_error is not None: + if effective_entries(total, total_error) < min_neff: + return False + 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 _bin_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_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 _slice_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_slice( + sig_hist, + bkg_hists_by_name, + right, + first_bin, + min_sum, + min_neff=0.0, + sig_mode="sb", + min_each=0.0, + min_proc_neff=0.0, + exempt=(), +): + """Among every candidate [left, right] with summed background > min_sum, + pick the one maximizing S/sqrt(B + sigma_B^2) -- not just the first one that + clears it. This is what actually drives where the cut lands; the content + floor is only a validity gate. Falls back to first_bin (best effort) if no + candidate clears the floor anywhere.""" + best_left = None + best_sig = -1.0 + for left in range(right, first_bin - 1, -1): + bkg_y = _bkg_yields(bkg_hists_by_name, left, right) + b_err = _total_bkg_error(bkg_hists_by_name, left, right) + bkg_e = ( + _bkg_errors(bkg_hists_by_name, left, right) if min_proc_neff > 0 else None + ) + if not _slice_passes( + bkg_y, + min_sum, + b_err, + min_neff, + bkg_e, + min_each, + min_proc_neff, + exempt, + ): + continue + s = _integral(sig_hist, left, right) + b = sum(bkg_y.values()) + sig = significance(s, b, b_err, sig_mode) + if sig > best_sig: + best_sig = sig + best_left = left + return best_left if best_left is not None else first_bin + + +def find_slices( + sig_hist, + bkg_hists_by_name, + n_slices, + first_bin, + last_bin, + min_sum, + min_neff=0.0, + sig_mode="sb", + min_each=0.0, + min_proc_neff=0.0, + min_frac=0.0, +): + """Split [first_bin, last_bin] into exactly `n_slices` ranges (fixed count, + required so every mass point shares the same category list), scanning right + to left. Each slice (except the final, leftover one) is placed to maximize + signal significance among boundaries with summed background > min_sum. + + Which backgrounds are minor enough to be exempt from the per-process floors is + decided once here, over the whole [first_bin, last_bin] range, and held fixed + for every candidate window -- see _slice_passes on why it cannot be re-judged + per window. + """ + exempt = ( + minor_backgrounds(bkg_hists_by_name, first_bin, last_bin, min_frac) + if (min_each > 0 or min_proc_neff > 0) + else set() + ) + slices = [] + right = last_bin + for slice_idx in range(n_slices): + if slice_idx == n_slices - 1 or right <= first_bin: + lo = first_bin if right >= first_bin else right + slices.append((lo, right if right >= lo else lo)) + right = lo - 1 + continue + left = grow_slice( + sig_hist, + bkg_hists_by_name, + right, + first_bin, + min_sum, + min_neff, + sig_mode, + min_each, + min_proc_neff, + exempt, + ) + slices.append((left, right)) + right = left - 1 + slices.reverse() + return slices + + +def signal_quantile_ranges(sig_hist, n_bins, first_bin, last_bin): + """Split [first_bin, last_bin] into `n_bins` ranges each holding an equal share + of the signal. + + This is what puts the bins where the resonance is. The axis is the same range at + every mass point (bbWW's HME runs 0-1500 GeV), 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 _bin_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_bins( + sig_hist, + bkg_hists_by_name, + max_bins, + first_bin, + last_bin, + min_each, + min_frac=0.0, + min_neff=0.0, +): + """Bins inside one slice: signal quantiles for the edges, background gates for + the count. + + The two axes need opposite rules, which is why this is not find_slices(). The + previous version scanned right-to-left growing each bin leftward until the + backgrounds cleared their floors, mirroring find_slices(). That is right for the + sliced axis, where the signal piles up at one end so resolution belongs there, and + wrong for an axis whose signal sits in the middle with both tails empty -- as + bbWW's HME does, peaking at HME ~ MX. + Starting from the top of the axis spent the bin budget on the empty + upper 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 bin_budget(bkg_hists_by_name, lo, hi, max_bins_per_slice, bkg_per_bin): + """How many bins this slice can actually afford. + + A fixed max_bins_per_slice is applied blind to slice content: the high-mass boosted + slices hold ~1.5 total background events and were still being split into 10 + bins, i.e. ~0.15 events per bin. Capping at B_slice / bkg_per_bin ties the + binning to the statistics that are really there. Returns max_bins_per_slice + unchanged when bkg_per_bin <= 0 (feature off). + """ + if bkg_per_bin <= 0: + return max_bins_per_slice + total = sum(_bkg_yields(bkg_hists_by_name, lo, hi).values()) + if total <= 0: + return 1 + return max(1, min(max_bins_per_slice, int(total / bkg_per_bin))) + + +def discover_binning( + sig2d, + bkg2d_by_name, + n_slices, + max_bins_per_slice, + min_slice_sum, + min_bin_each, + min_slice_bkg_neff=0.0, + min_bkg_frac=0.0, + min_bin_bkg_neff=0.0, + bkg_per_bin=0.0, + sig_mode="sb", + min_slice_bkg_each=0.0, + min_slice_bkg_each_neff=0.0, +): + """bkg2d_by_name: {background_name: [hist per discovery era, ...]}. The list + is usually a single era's own histogram (standalone limit) or all of a + meta-era's sub-eras (combined limit) -- see --discovery-eras. + 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 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() + slices = find_slices( + sig2d, + bkg2d_by_name, + n_slices, + 1, + nx, + min_slice_sum, + min_slice_bkg_neff, + sig_mode, + min_slice_bkg_each, + min_slice_bkg_each_neff, + min_bkg_frac, + ) + slices = extend_outer_edges(slices, 0, nx + 1) + + result = [] + for xlo, xhi in 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 = bin_budget(bkg_y_by_name, 1, ny, max_bins_per_slice, bkg_per_bin) + bin_ranges = find_bins( + sig_y, + bkg_y_by_name, + n_bins, + 1, + ny, + min_bin_each, + min_bkg_frac, + min_bin_bkg_neff, + ) + bin_ranges = extend_outer_edges(bin_ranges, 0, ny + 1) + result.append({"x_range": (xlo, xhi), "y_ranges": bin_ranges}) + return result + + +def bin_edges(y_axis, y_ranges): + """Physical edges of the discovered y ranges, for booking the output TH1. + + The rebinned shapes used to be booked as n_bins over [0, n_bins], which threw + the axis scale away and left every plot labelled by bin index. The ranges are + contiguous and ordered, so the edges are each range's low edge plus the last + range's upper edge. extend_outer_edges() pushes the outer ranges into the + underflow/overflow bins, which have no finite edge of their own -- those are + clamped to the axis limits. + """ + n = y_axis.GetNbins() + edges = [y_axis.GetBinLowEdge(max(lo, 1)) for lo, _ in y_ranges] + edges.append(y_axis.GetBinUpEdge(min(y_ranges[-1][1], n))) + return edges + + +def mkdir_titled(directory, path, title): + """mkdir a nested path, putting `title` on the leaf directory only. + + 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 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 slices_to_record(slices, x_axis, y_axis): + """The discovered structure as plain data, for binning.json. + + Both forms are written. The bin index ranges are what the code actually cuts on and + are what a replay restores; the physical edges alongside them are what a human reads, + and what makes the record meaningful next to a plot. + """ + return { + "n_x_bins": x_axis.GetNbins(), + "n_y_bins": y_axis.GetNbins(), + "slices": [ + { + "x_range": list(sl["x_range"]), + "x_edges": x_edges, + "y_ranges": [list(r) for r in sl["y_ranges"]], + "y_edges": bin_edges(y_axis, sl["y_ranges"]), + } + for sl, x_edges in zip(slices, slice_ranges(x_axis, slices)) + ], + } + + +def record_to_slices(record, x_axis, y_axis, where): + """Inverse of slices_to_record(): what discover_binning() would have returned. + + The axis sizes are checked rather than trusted. Bin indices only mean anything against + the axes they were found on, so a binning.json replayed over input with a different + binning would otherwise cut the shapes in silently wrong places -- which is exactly + what freezing a binning is supposed to rule out. + """ + for name, axis, recorded in ( + ("x", x_axis, record["n_x_bins"]), + ("y", y_axis, record["n_y_bins"]), + ): + if axis.GetNbins() != recorded: + raise RuntimeError( + f"{where}: recorded binning was found on a {recorded}-bin {name} axis, " + f"but the input has {axis.GetNbins()}. The binning.json does not belong " + "to this input." + ) + return [ + { + "x_range": tuple(sl["x_range"]), + "y_ranges": [tuple(r) for r in sl["y_ranges"]], + } + for sl in record["slices"] + ] + + +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", 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( + sources, + channel, + category, + cfg, + mass, + era, + discovery_files, + knobs, + frozen=None, +): + """Rebin one channel/category, returning the binning it used (None if skipped). + + The edges are discovered once from `discovery_files` -- every source era summed -- and + then applied to each source era separately, so a group era's members are cut on the + combination's statistics but stay in their own files. The datacard step still sums + them, which is what keeps a per-era lnN expressible: scaling one sub-era's own shape + needs that sub-era to still exist as a shape. + + `sources` is [(source_era, in_file, out_file)]; for a plain era there is one. + + With `frozen` given the edges are replayed from a previous run's binning.json and no + optimisation happens at all; the discovery files are then only read for the axes. + """ + min_signal = knobs["min_signal"] + prefix = f"{channel}/{category}/" + # The key list comes from the first source era; a systematic that only some eras carry + # is filled in from their nominal by sum_over_sources(). + cat_dir = sources[0][1].Get(f"{channel}/{category}") + if not cat_dir: + print(f" [skip] {channel}/{category}: not found in {sources[0][1].GetName()}") + return + + 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 + + where = f"{era}/MX={mass}/{channel}/{category}" + if frozen is not None: + slices = record_to_slices( + frozen, disc_sig.GetXaxis(), disc_sig.GetYaxis(), where + ) + else: + slices = discover_binning( + disc_sig, + disc_bkg_by_name, + knobs["n_slices"], + knobs["max_bins_per_slice"], + knobs["min_slice_bkg_sum"], + knobs["min_bin_bkg_each"], + knobs["min_slice_bkg_neff"], + knobs["min_bkg_frac"], + knobs["min_bin_bkg_neff"], + knobs["bkg_per_bin"], + knobs["significance_mode"], + knobs["min_slice_bkg_each"], + knobs["min_slice_bkg_each_neff"], + ) + # The selection each slice stands for is otherwise nowhere in the output: the sliced + # categories are named by index and the surviving axis is the rebinned 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) + + # One shared set of edges, applied to every source era in its own file. + for source_era, in_file, out_file in sources: + 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"]), + ) + cat_dir = in_file.Get(f"{channel}/{category}") + if not cat_dir: + print(f" [skip] {source_era} {channel}/{category}: not in the input") + continue + for key in [k.GetName() for k in cat_dir.GetListOfKeys()]: + hist2d = get_hist(in_file, prefix + key) + if hist2d is None or hist2d.GetDimension() != 2: + continue + for slice_idx, h in enumerate(rebin_hist_2d(hist2d, slices, key, naming)): + out_file.cd(f"{channel}/{naming.name(category, slice_idx)}") + h.Write(key) + + return slices_to_record(slices, disc_sig.GetXaxis(), disc_sig.GetYaxis()) + + +def run( + input_dir, + output_dir, + config_path, + era, + knobs, + frozen_binning=None, +): + """Produce one era of the datacard configuration's `eras:` list. + + Which era it is decides everything. A plain era is binned on its own statistics, for + a standalone limit. An era that is a key of `era_groups:` is binned on its members' + summed statistics -- a combination supports finer bins than any single era can -- and + those edges are then applied to each member separately. + + Output layout is "///.root", plus the + binning.json that produced it. One era per --output directory, since the same source + era carries different edges under different target eras and they must not overwrite + each other. The members are kept in their own files rather than summed here: the + datacard step sums them, and a per-era lnN can only be built by scaling a sub-era's + own shape, which a pre-summed shape no longer has. + """ + cfg = load_config(config_path) + # The pattern that names the sliced categories comes from the binning configuration + # alongside the edges themselves: it is the writer's choice, and every reader of these + # shapes recovers the base category from the same pattern in the datacard config. + cfg["naming"] = CategoryNaming(knobs["category_pattern"]) + cfg["slice_var"] = knobs["slice_var"] + # The datacard configuration lists the sliced names this script *writes* + # ("SR/res2b_dnn0"); the 2D input it *reads* is keyed by the base categories those + # slices are cut from ("SR/res2b"). Recover them with the same pattern that names + # them, so the two ends cannot disagree about which is which. A configuration that + # lists base names already is unchanged by this -- base() returns an unsliced name + # as-is. + cfg["categories"] = list( + dict.fromkeys(cfg["naming"].base(c) for c in cfg["categories"]) + ) + model = cfg["model"] + + # A group era is built from its members; a plain era is built from itself. + source_eras = cfg["era_groups"].get(era, [era]) + if era in cfg["era_groups"]: + print(f"{era} is a group of {source_eras}: binning on their summed statistics") + else: + print(f"{era} is a standalone era: binning on its own statistics") + + record = { + "slice_var": knobs["slice_var"], + "category_pattern": cfg["naming"].pattern, + "knobs": {k: v for k, v in sorted(knobs.items())}, + "binning": {}, + } + by_mass = record["binning"].setdefault(era, {}) + + for mass in cfg["mass_values"]: + # One handle per source era, used both to discover the edges and to write the + # shapes -- so the binning is derived from exactly the statistics it is applied to. + sources = [] + for src in source_eras: + in_file = open_input_file( + input_dir, model, src, mass, cfg["signal_param_name"] + ) + # getInputFileName() yields "//.root", which is + # the layout input_file_pattern resolves against. --output is one era's own + # directory, so the era being produced is not repeated inside it. + out_path = os.path.join( + output_dir, + model.getInputFileName(src, {cfg["signal_param_name"]: mass}), + ) + os.makedirs(os.path.dirname(out_path), exist_ok=True) + sources.append((src, in_file, ROOT.TFile.Open(out_path, "RECREATE"))) + discovery_files = [f for _, f, _ in sources] + + print( + f"Rebinning {era} MX={mass} from {len(sources)} source era(s) -> " + f"{os.path.join(output_dir, era)}" + ) + by_channel = by_mass.setdefault(str(mass), {}) + for channel in cfg["channels"]: + for category in cfg["categories"]: + frozen = None + if frozen_binning is not None: + frozen = lookup_frozen(frozen_binning, era, mass, channel, category) + used = process_category( + sources, + channel, + category, + cfg, + mass, + era, + discovery_files, + knobs, + frozen=frozen, + ) + if used is not None: + by_channel.setdefault(channel, {})[category] = used + + for _, in_file, out_file in sources: + out_file.Close() + in_file.Close() + + # Written beside the shapes it describes, in the same output as them: the binning is + # a product of this run, not configuration, so it travels with what it produced and a + # reader never has to work out which binning a set of shapes came from. + os.makedirs(output_dir, exist_ok=True) + json_path = os.path.join(output_dir, BINNING_JSON) + with open(json_path, "w") as f: + json.dump(record, f, indent=2, sort_keys=True) + print(f"Wrote {json_path}") + + +if __name__ == "__main__": + import argparse + + parser = argparse.ArgumentParser( + description="Rebin 2D histograms into significance-sliced 1D shapes.\n\n" + "A standalone pre-step, not part of the datacard chain: it writes a " + "'//.root' tree in the same layout HistMergerTask " + "produces, so a production of these is consumed by pointing the chain's " + "--hists-version at it. It also writes " + BINNING_JSON + ", which --binning " + "replays to reproduce a binning instead of re-deriving one.", + formatter_class=argparse.RawDescriptionHelpFormatter, + ) + 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( + "--binning-config", + required=False, + type=str, + default=None, + help="binning yaml holding the knobs below; anything it does not set takes the " + "built-in default, and an explicit flag overrides both", + ) + parser.add_argument( + "--binning", + required=False, + type=str, + default=None, + help=f"a previous run's {BINNING_JSON}. Given one, the recorded edges are applied " + "as-is and nothing is optimised -- this is how a binning is frozen and a " + "production reproduced. The knobs are then unused", + ) + parser.add_argument( + "--era", + required=True, + type=str, + help="the era of the datacard configuration's `eras:` list to produce. A plain " + "era is binned on its own statistics; an era that is a key of `era_groups:` is " + "binned on its members' summed statistics and its shapes are that sum. Each is " + "written under its own era name, since a combination supports finer bins than " + "any single era and the two must not overwrite each other", + ) + + # Knob overrides. All default to None so that "not given" is distinguishable from + # "given the same value as the default", which is what lets --binning-config win. + knob_args = { + "n_slices": (int, "fixed number of slices per category (same for every mass)"), + "category_pattern": ( + str, + "pattern naming the sliced categories, e.g. " + "'{base_category}_dnn{slice_idx}'; must use both {base_category} and " + "{slice_idx}, and the datacard configuration reading these shapes must " + "declare the same one", + ), + "slice_var": ( + str, + "name of the sliced axis, used to label each slice directory with the " + "selection it stands for (e.g. 'DNN' -> '0.80 < DNN < 1.00')", + ), + "max_bins_per_slice": (int, "bins to aim for inside each slice"), + "min_slice_bkg_sum": ( + float, + "minimum summed-background yield required in a slice", + ), + "min_slice_bkg_neff": ( + float, + "minimum effective MC entries of the summed background for a slice boundary " + "to be selectable", + ), + "min_slice_bkg_each": ( + float, + "minimum yield required of every non-negligible background in a slice", + ), + "min_slice_bkg_each_neff": ( + float, + "minimum effective MC entries of every non-negligible background in a slice; " + "much stronger than the yield floor and correspondingly expensive", + ), + "min_bin_bkg_each": ( + float, + "minimum yield required of every non-negligible background in a bin", + ), + "min_bin_bkg_neff": ( + float, + "minimum effective MC entries of the summed background in a bin", + ), + "bkg_per_bin": ( + float, + "summed background to aim for per bin; overrides max_bins_per_slice when it " + "binds", + ), + "min_bkg_frac": ( + float, + "backgrounds below this fraction of the category total are exempt from the " + "per-background floors, so a negligible process cannot veto every boundary", + ), + "min_signal": (float, "minimum signal integral for a category to be rebinned"), + } + for name, (typ, help_text) in knob_args.items(): + parser.add_argument( + f"--{name.replace('_', '-')}", + required=False, + type=typ, + default=None, + help=help_text, + ) + parser.add_argument( + "--significance-mode", + required=False, + type=str, + default=None, + choices=["sb", "asimov"], + help="figure of merit for slice boundaries: 'sb' = S/sqrt(B+sigmaB^2), " + "'asimov' = Poisson-correct Asimov significance (valid at low B)", + ) + args = parser.parse_args() + + overrides = {name: getattr(args, name) for name in knob_args} + overrides["significance_mode"] = args.significance_mode + knobs = load_binning_config(args.binning_config, overrides) + + frozen_binning = None + if args.binning: + with open(args.binning, "r") as f: + frozen_binning = json.load(f) + print(f"Replaying the binning recorded in {args.binning}; not optimising") + + run( + args.input, + args.output, + args.config, + args.era, + knobs, + frozen_binning=frozen_binning, + ) diff --git a/common/tools.py b/common/tools.py index 930b6df..87362a8 100644 --- a/common/tools.py +++ b/common/tools.py @@ -1,7 +1,117 @@ 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 top-level ``category_pattern``, or the default + when it declares none -- a configuration whose categories were never sliced has + nothing to parse back, and any name is then its own base. + + A configuration reading shapes from bin_opt_2d/rebin_2d.py must repeat the pattern + that run used, since that is the only way its sliced category names can be taken + apart again. + """ + return cls(cfg.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, and bin_opt_2d/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..2f5bfac 100644 --- a/dc_make/create_datacards.py +++ b/dc_make/create_datacards.py @@ -34,7 +34,6 @@ default=None, help="parameter values to run only certain masses", ) - for param in DatacardMaker.customizeble_parameters: parser.add_argument( f"--{param}", diff --git a/dc_make/maker.py b/dc_make/maker.py index a5b6902..ac19660 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,12 @@ 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, + **kwargs, ): self.cb = CombineHarvester() @@ -53,9 +61,17 @@ def __init__( self.analysis = cfg["analysis"] self.eras = cfg["eras"] self.channels = cfg["channels"] - self.categories = cfg["categories"] + # The categories are exactly the ones the configuration lists -- the datacard bins + # are whatever directories the input shapes contain, and nothing here derives them. + # For input that a 2D->1D rebinning produced that means the sliced names + # ("SR/res2b_dnn0"), and `category_pattern` is how they are taken apart again to + # group the slices of one base category. + self.naming = CategoryNaming.fromConfig(cfg) + self.categories = 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 +143,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 +158,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 +213,58 @@ def cbCopy(self, param_str, process, era, channel, category): def ECC(self): return itertools.product(self.eras, self.channels, self.categories) + @staticmethod + def lnNIsEraDependent(unc): + """Whether an lnN's value depends on which era it is applied to. + + Only an era-dependent lnN needs the meta-era shape treatment. One that applies + uniformly -- no `eras:` anywhere in its entry -- scales the whole summed yield by + the same factor, which is what an lnN already says, so turning it into a template + buys nothing and states the same uncertainty a different way: combine morphs a + shape and takes an lnN as an exact log-normal. + + For MultiValueLnNUncertainty the eras are read off the value keys, not off + `unc.eras`. Uncertainty.fromConfig reassigns its `args` dict inside the sub-entry + loop and passes the result to the constructor, so a multi-value entry carries + whatever the *last* sub-entry happened to scope by -- top_mass ends up with + processes ('ST',) and lumi_1_13p6TeV with the 2023 eras. Nothing reads those today + (addUncertainty decides multi-value applicability from getUncertaintyForProcess), + but they are not to be trusted. + """ + if isinstance(unc, MultiValueLnNUncertainty): + # key is (processes, eras, channels, categories), built in fromConfig. + return any(key[1] for key in unc.values) + return bool(unc.eras) + + def uncAppliesTo(self, unc, process, era, channel, category): + """Whether an uncertainty applies, with a meta-era inheriting its sub-eras'. + + Uncertainty.appliesTo matches the era it is given against the uncertainty's own + `eras:` list, and a meta-era is not in it -- CMS_scale_j_2022 is scoped to + Run3_2022, and Run3_Early is not that string. Answering False is right at that + level; the meta-era is a fact this class knows and uncertainty.py does not. + + A nuisance scoped to one sub-era does apply to the combination: it varies that + era's part of the summed shape. Without this every per-era shape systematic was + dropped from a meta-era card silently -- 64 of 76 for HH->bbWW, including the jet + energy scale and resolution. Per-era lnN escaped only because they are intercepted + for meta-eras earlier, in _addMetaEraLnNAsShapeUnc. + """ + eras = self.getSubEras(era) if self.isMetaEra(era) else [era] + return any(unc.appliesTo(process, e, channel, category) for e in eras) + + 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,52 +283,281 @@ 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 + + @staticmethod + def readHist(file, hist_name): + """A histogram from a file, owned by Python rather than by ROOT. + + TFile.Get() hands back an object ROOT keeps alive for the lifetime of the file, so + dropping the Python reference frees nothing. That is affordable when the objects + are the small 1D shapes a datacard bin is made of; it is not when they are the 2D + inputs a binning is cut from, where one era's build walked through ~15 GB. Taking + ownership lets refcounting free each one after it has been sliced, while anything + still referenced -- a cached shape -- stays alive as usual. + + Returns None if the path is missing: TFile.Get() on a fully-missing nested path + returns a PyROOT wrapper around a null pointer, which is not `is None` but is + falsy, so `if obj:` is the correct check. + """ + obj = file.Get(hist_name) + if not obj: + return None + ROOT.SetOwnership(obj, True) + return obj + + def _loadBinnedHist(self, file, era, channel, category, model_params, hist_name): + hist = self.readHist(file, 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 - def getShape( + # 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: + # A per-era nuisance varies its own era and leaves the rest alone, so a sub-era + # it is not scoped to contributes its *nominal*. Two things depend on that: the + # variation only exists in its own era's file, so asking the others for it + # raises; and the summed Up shape has to carry the full yield, or its ratio to + # the nominal -- which is what the fit reads -- would be meaningless. + applies = unc_name is None or self.uncAppliesTo( + self.uncertainties[unc_name], process, sub_era, channel, category + ) + sub_hist = self._rawShape( + process, + sub_era, + channel, + category, + model_params, + unc_name if applies else None, + unc_scale if applies else None, + ) + 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 _rawShape( self, process, era, @@ -221,8 +567,26 @@ def getShape( unc_name=None, unc_scale=None, ): + """One process's shape in one bin, as read, with no negative-bin validation. + + Separate from getShape() for getCombinedShape()'s sake: a meta-era's sub-eras are + summed before anything is checked, because a sub-era can dip negative on its own + statistics while the combined shape -- the one that actually enters the datacard -- + is fine. Validating per sub-era would reject shapes that are good, and would also + change them: resolveNegativeBins rebalances in place, so clamping each sub-era and + then summing does not give the same shape as summing and then clamping. + """ 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, + ) + 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,8 +594,13 @@ 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 + # Validated, not raw: asimov data is the sum of the + # per-process shapes as they enter the fit, each resolved + # under its own tolerance, never raw shapes summed and then + # checked under data_obs's own (untolerant) settings. The + # meta-era branch of getCombinedShape mirrors this. bkg_hist = self.getShape( bkg_proc, era, channel, category, model_params ) @@ -249,7 +618,7 @@ def getShape( hist_name = f"{channel}/{category}/{subp}" if unc_name and unc_scale: hist_name += f"_{unc_name}_{unc_scale}" - subhist = file.Get(hist_name) + subhist = self.readHist(file, hist_name) if subhist == None: raise RuntimeError( f"Cannot find histogram {hist_name} in {file.GetName()}" @@ -262,7 +631,7 @@ def getShape( else: if unc_name and unc_scale: hist_name += f"_{unc_name}_{unc_scale}" - hist = file.Get(hist_name) + hist = self.readHist(file, hist_name) if hist == None: raise RuntimeError( f"Cannot find histogram {hist_name} in {file.GetName()}" @@ -298,65 +667,188 @@ 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( - 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, - ) + self.shapes[key] = hist + return self.shapes[key] - 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 _isUnvalidatedSignal(self, process, unc_name, unc_scale): + """Whether this is a nominal signal shape, which is never negative-bin checked. + + Those are collected into signal_hists_by_key to define the relevant (signal + carrying) bins that the check itself consults, so they are the input to the rule + rather than subject to it. + """ + return process.is_signal and not (unc_name and unc_scale) + + def getShape( + self, + process, + era, + channel, + category, + model_params, + unc_name=None, + unc_scale=None, + ): + """One process's shape in one bin, negative-bin validated -- what a datacard bin + is built from. + + A meta-era is summed from its sub-eras first (getCombinedShape), which validates + the sum rather than the parts; see _rawShape() for why. + """ + if self.isMetaEra(era): + return self.getCombinedShape( + process, era, channel, category, model_params, unc_name, unc_scale + ) + + raw = self._rawShape( + process, era, channel, category, model_params, unc_name, unc_scale + ) + if raw is None or self._isUnvalidatedSignal(process, unc_name, unc_scale): + return raw + + # Cached separately from the raw form, and validated on a copy: the check + # rebalances in place, and the raw shape is still needed unmodified by + # getCombinedShape, which sums the sub-eras before checking anything. + key = ( + self.model.getInputFileName(era, model_params), + process.name, + era, + channel, + category, + unc_name, + unc_scale, + "validated", + ) + if key not in self.shapes: + hist = raw.Clone() + hist.SetDirectory(0) + self.resolveOrRaiseNegativeBins( + hist, + process, + era, + channel, + category, + model_params, + unc_name, + unc_scale, + ) 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 the rebinning 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 = self.readHist(file, hist_name) + if obj is not None 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 +878,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 +890,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,48 +909,111 @@ 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 else: add(None, "*", proc) + def warnUnusedUncertainties(self): + """Name the configured uncertainties that reached no bin of this card. + + addUncertainty walks the process/era/channel/category product and adds nothing + where the uncertainty does not apply, so one that matches nothing at all is not an + error anywhere -- it simply never appears, and the card is quietly missing a + nuisance. That is worth saying out loud: the usual cause is a `processes:` or + `categories:` pattern that matches no name, and the usual reason for that is a + pattern written with a placeholder the maker never substitutes (a signal's ${MX} + is already resolved to the mass being built by the time it is matched, so it has + to be written as a ^-anchored regex). + + Some entries legitimately reach nothing -- an era-scoped nuisance in a card for a + different era -- so this warns rather than raises. + """ + used = set(self.cb.syst_name_set()) + unused = [name for name in self.uncertainties if name not in used] + if unused: + print( + f"WARNING: {len(unused)} configured uncertainty(ies) were added to no bin " + f"of this card and are absent from it: {', '.join(sorted(unused))}. " + "Check their processes/eras/channels/categories patterns against the " + "names actually in the card." + ) + return unused + def addUncertainty(self, unc_name): unc = self.uncertainties[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)) + and self.lnNIsEraDependent(unc) + ): + 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 if isMVLnUnc - else unc.appliesTo(process, era, channel, category) + else self.uncAppliesTo(unc, process, era, channel, category) ) if not uncApplies: continue - if not process.hasCompatibleModelParams( - model_params, self.model.param_dependent_bkg - ): - continue nominal_shape = None shapes = {} @@ -495,7 +1054,7 @@ def addUncertainty(self, unc_name): ) cb_copy = self.cbCopy(param_str, proc, era, channel, category) cb_copy.AddSyst(self.cb, unc_name, unc_to_apply.type.name, systMap) - if unc_to_apply.type == UncertaintyType.shape: + if isinstance(unc_to_apply, ShapeUncertainty): shape_set = False def setShape(syst): @@ -521,22 +1080,37 @@ 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 isinstance( + unc, (LnNUncertainty, MultiValueLnNUncertainty) + ): + 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) if isMVLnUnc - else unc.appliesTo(process, era, channel, category) + else self.uncAppliesTo(unc, process, era, channel, category) ) if not uncApplies: continue - if not process.hasCompatibleModelParams( - params, self.model.param_dependent_bkg - ): - continue nominal_shape = None shapes = {} @@ -572,7 +1146,7 @@ def setShape(syst): cb_copy = self.cbCopy(param_str, proc_name, era, channel, category) cb_copy.AddSyst(self.cb, unc_name, unc_to_apply.type.name, systMap) - if unc_to_apply.type == UncertaintyType.shape: + if isinstance(unc_to_apply, ShapeUncertainty): def setShape(syst): syst.set_shapes( @@ -606,13 +1180,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 +1213,34 @@ def writeDatacards(self, output): param_list ).process(processes).WriteDatacard(tmp_dc_file, tmp_shape_file) + # Same breakdown by base category (all its 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 the rebinning 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 +1249,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) @@ -648,6 +1263,7 @@ def createDatacards(self, output, verbose=1): for unc_name in self.uncertainties.keys(): print(f"adding uncertainty: {unc_name}") self.addUncertainty(unc_name) + self.warnUnusedUncertainties() if self.autoMCStats["apply"]: self.cb.SetAutoMCStats( self.cb, diff --git a/dc_make/process.py b/dc_make/process.py index 70613c5..a6dd4a9 100644 --- a/dc_make/process.py +++ b/dc_make/process.py @@ -17,6 +17,7 @@ def __init__( max_n_sigma_for_negative_bins=1, allow_negative_integral=False, channels=[], + categories=[], ): self.name = name self.hist_name = hist_name @@ -32,6 +33,7 @@ def __init__( self.max_n_sigma_for_negative_bins = max_n_sigma_for_negative_bins self.allow_negative_integral = allow_negative_integral self.channels = channels + self.categories = categories if is_data and is_signal: raise RuntimeError("Data and signal flags cannot be set simultaneously") if is_asimov_data and not is_data: @@ -44,6 +46,18 @@ def __init__( else: self.type = "background" + def appliesToCategory(self, category): + """Whether this process contributes to `category`. Empty list = all of them. + + Entries match as prefixes of the datacard category name, so "SR/boosted" + covers SR/boosted_dnn0..3 without listing the slices: how many slices + a category is cut into is a binning parameter, and the process list should + not have to change when it does. + """ + if not self.categories: + return True + return any(category.startswith(prefix) for prefix in self.categories) + def __str__(self): str_rep = ( f"Process({self.name}, type={self.type}, subprocesses={self.subprocesses}" @@ -85,6 +99,7 @@ def fromConfig(entry, model): ) max_n_sigma_for_negative_bins = entry.get("max_n_sigma_for_negative_bins", 1) channels = entry.get("channels", []) + categories = entry.get("categories", []) if type(scale) == str: scale = eval(scale) if "param_values" not in entry: @@ -105,6 +120,7 @@ def fromConfig(entry, model): max_n_sigma_for_negative_bins=max_n_sigma_for_negative_bins, allow_negative_integral=allow_negative_integral, channels=channels, + categories=categories, ) ] @@ -133,6 +149,7 @@ def fromConfig(entry, model): allow_negative_integral=allow_negative_integral, params=param_dict, channels=channels, + categories=categories, ) ) return processes diff --git a/dc_make/uncertainty.py b/dc_make/uncertainty.py index 1142237..4981f7c 100644 --- a/dc_make/uncertainty.py +++ b/dc_make/uncertainty.py @@ -293,11 +293,18 @@ def checkValue(self): # if not isinstance(processes, tuple) or not all(isinstance(p, str) for p in processes): # raise ValueError(f"Invalid processes list: {processes}. Must be a list of strings.") - def getUncertaintyForProcess(self, process): + def getUncertaintyForProcess(self, process, era=None, channel=None, category=None): for key in self.values.keys(): processes, eras, channels, categories = key - if process in processes: - return self.values[key] + if processes and process not in processes: + continue + if era is not None and eras and era not in eras: + continue + if channel is not None and channels and channel not in channels: + continue + if category is not None and categories and category not in categories: + continue + return self.values[key] return None def valueToMap(self, unc_value, digits=3): diff --git a/law/CreateDatacardsTask.py b/law/CreateDatacardsTask.py new file mode 100644 index 0000000..e8f969c --- /dev/null +++ b/law/CreateDatacardsTask.py @@ -0,0 +1,352 @@ +import contextlib +import law +import luigi +import os + +from string import Template + +from FLAF.RunKit.run_tools import ps_call +from FLAF.run_tools.law_customizations import HTCondorWorkflow, copy_param + +from StatInference.common.tools import CategoryNaming, importROOT + +from .PreprocessShapesTask import PreprocessShapesTask +from .StatInferenceTask import StatInferenceTask + + +class CreateDatacardsTask(StatInferenceTask, HTCondorWorkflow, law.LocalWorkflow): + max_runtime = copy_param(HTCondorWorkflow.max_runtime, 2.0) + n_cpus = copy_param(HTCondorWorkflow.n_cpus, 1) + + # If set, build datacards for this meta-era instead of self.period. self.period must + # then be one of its real sub-eras -- it is only used to construct a valid FLAF Setup + # (Setup.getGlobal requires a real, known period; a meta-era name isn't one). + meta_era = luigi.Parameter(default="") + + # Stacked plots of the shapes going into the cards. This task is where they belong: + # it is the only one holding every sub-era's histograms at once, which is what the + # merged shape (and hence each datacard bin) is built from. + make_plots = luigi.BoolParameter(default=True) + + @property + def datacard_era(self): + """The era datacards are actually built for: self.meta_era if set, else self.period.""" + return self.meta_era or self.period + + def get_sub_periods(self): + """Real periods whose histograms feed into datacard_era: its constituent + sub-eras for a meta-era, otherwise just [self.period].""" + return self.get_era_groups().get(self.datacard_era, [self.period]) + + def input_hist_reqs(self): + """{key: task} for the shapes the cards are built from. + + The configuration's `preprocess:` step if it declares one, otherwise the merged + histograms directly. Either way what arrives is "//.root", + so nothing below here knows which it got. + """ + if self.preprocess_config(): + return { + "PreprocessShapes": PreprocessShapesTask.req( + self, meta_era=self.meta_era, branches=() + ) + } + return { + f"MergedHists_{era}_{variable}": req + for (era, variable), req in self.merged_hist_reqs( + self.get_sub_periods() + ).items() + } + + def workflow_requires(self): + return self.input_hist_reqs() + + def requires(self): + return list(self.input_hist_reqs().values()) + + def create_branch_map(self): + return {0: None} + + def output(self): + # fs_default. Note that combine cannot read these directly: + # ResonantLimitsTask mirrors them back to datacards_dir() before handing them to + # dhi -- see ResonantLimitsTask.stage_datacards. + return self.output_dir_target(self.version, "Datacards", self.datacard_era) + + def run(self): + statInf_entry = self.global_params["StatInference"] + config = self.datacard_config_path() + # ${ERA} in hist_bins names the era the cards are for, the same way the model's + # input_file_pattern does. Each era has its own binning -- derived from its own + # statistics, or its members' summed -- so one configuration serves them all. + hist_bins_rel = statInf_entry.get("hist_bins") + hist_bins = ( + os.path.join( + self.ana_path(), + Template(hist_bins_rel).safe_substitute(ERA=self.datacard_era), + ) + if hist_bins_rel + else None + ) + param_values = statInf_entry.get("param_values", []) + create_datacards_py = os.path.join( + self.ana_path(), "StatInference", "dc_make", "create_datacards.py" + ) + with contextlib.ExitStack() as stack: + if self.preprocess_config(): + # PreprocessShapesTask already wrote the "//.root" + # layout, so its output directory is the base input_file_pattern resolves + # against -- nothing to stage. + reqs = self.input_hist_reqs()["PreprocessShapes"] + base_dir_local = stack.enter_context(reqs.output().localize("r")) + else: + targets = { + key: req.output() + for key, req in self.merged_hist_reqs( + self.get_sub_periods() + ).items() + } + self.check_inputs(targets) + base_dir_local = stack.enter_context(self.stage_inputs(targets)) + local_output = stack.enter_context(self.output().localize("w")) + cmd = [ + "python3", + create_datacards_py, + "--input", + base_dir_local.abspath, + "--output", + local_output.abspath, + "--config", + config, + "--eras", + self.datacard_era, + ] + if hist_bins: + cmd += ["--hist-bins", hist_bins] + if len(param_values) > 0: + param_values_str = ",".join(str(v) for v in param_values) + cmd += ["--param_values", param_values_str] + ps_call(cmd, env=self.cmssw_env, verbose=1) + + if self.make_plots: + self.plot_rebinned_shapes( + base_dir_local.abspath, config, local_output.abspath + ) + + def plot_variable(self, variable): + """The variable whose axis the shapes are binned along, for histograms.yaml. + + For input a 2D->1D rebinning produced, that is the y variable of the 2D entry -- + histograms.yaml records both as ``var_list: [x, y]``, so the axis metadata the + plotter needs is already described and does not have to be restated here. For + input that was always 1D, the variable is its own answer. + """ + try: + import FLAF.Common.Setup as Setup + + hists = Setup.Setup(self.ana_path(), self.period, self.version).hists + var_list = hists[variable].get("var_list") + if var_list and len(var_list) > 1: + return var_list[1] + except Exception as e: + print(f"Warning: no var_list for {variable} ({e}); plotting it as itself") + return variable + + @staticmethod + def _stitch_grid(panels, out_path): + """Lay rendered panels out on one page: a row per channel, a column per slice. + + The slices of a base category are one physical selection cut into pieces, and the + fit sees them together -- a slice that looks reasonable in eMu and pathological in + eE is obvious side by side and invisible in separate files. Only page geometry + happens here; every panel is drawn by HistPlotter.py, so there is still one + plotting implementation. + + `panels` is [[path or None, ...] per column] per row; a None leaves its cell blank, + which is what keeps column N under column N when a slice was skipped. + """ + from pypdf import PageObject, PdfReader, PdfWriter, Transformation + + pages = { + p: PdfReader(p).pages[0] for row in panels for p in row if p is not None + } + if not pages: + return False + w = max(float(p.mediabox.width) for p in pages.values()) + h = max(float(p.mediabox.height) for p in pages.values()) + n_rows, n_cols = len(panels), max(len(r) for r in panels) + + sheet = PageObject.create_blank_page(width=w * n_cols, height=h * n_rows) + for r, row in enumerate(panels): + for c, path in enumerate(row): + if path is None: + continue + # PDF origin is bottom-left, so the first row goes at the top. + sheet.merge_transformed_page( + pages[path], + Transformation().translate(c * w, (n_rows - 1 - r) * h), + ) + writer = PdfWriter() + writer.add_page(sheet) + with open(out_path, "wb") as f: + writer.write(f) + return True + + @staticmethod + def _present_keys(shape_file, cfg): + """The (channel, region, category) triples the summed shape file actually holds.""" + ROOT = importROOT() + present = set() + f = ROOT.TFile.Open(shape_file) + if not f or f.IsZombie(): + return present + try: + for category in cfg["categories"]: + region, _, cat = category.rpartition("/") + for channel in cfg["channels"]: + path = "/".join(x for x in (channel, region, cat) if x) + if f.Get(path): + present.add((channel, region, cat)) + finally: + f.Close() + return present + + def _stitch_grids(self, cfg, plots_dir, variable, panel_of_key): + """One grid per base category: its slices across, the channels down.""" + naming = CategoryNaming.fromConfig(cfg) + bases = {} + for category in cfg["categories"]: + region, _, cat = category.rpartition("/") + base, slice_idx = naming.split(cat) + bases.setdefault((region, base), {})[slice_idx] = cat + for (region, base), slices in bases.items(): + # None for an unsliced category, which is a single-column grid of channels. + order = sorted(slices, key=lambda i: (i is None, i)) + panels = [ + [ + (lambda p: p if p and os.path.exists(p) else None)( + panel_of_key.get(f"{channel}:{slices[i]}:{region}") + ) + for i in order + ] + for channel in cfg["channels"] + ] + out = os.path.join( + plots_dir, f"{variable}_{region.replace('/', '_')}_{base}_grid.pdf" + ) + if self._stitch_grid(panels, out): + print(f"Wrote grid {out}") + + def plot_rebinned_shapes(self, base_dir, config, output_dir): + """Stacked plots of the shapes the cards are built from, one per datacard bin. + + Runs FLAF's own HistPlotter.py -- the script HistPlotTask uses -- rather than a + plotter of our own, so these come out in the same style as every other plot in the + analysis and there is one implementation to maintain. It is agnostic about which + categories exist: it plots whatever `channel:category:region` keys it is handed, + so the sliced names need nothing on the FLAF side. + + Runs in the default env rather than cmssw_env: the plotter needs PlotKit + (matplotlib/mplhep), same as HistPlotTask. + """ + import glob + + cfg = self.get_config_data() + plots_dir = os.path.join(output_dir, "plots") + os.makedirs(plots_dir, exist_ok=True) + plotter = os.path.join(self.ana_path(), "FLAF", "Analysis", "HistPlotter.py") + + # The datacard bin is the sum over the sub-eras (getCombinedShape does the same at + # card-build time), so the sub-era files are hadd'ed into the one file the plotter + # reads. Plotting them separately would show something the fit never sees. + for variable in self.get_required_variables(): + inputs = [ + p + for era in self.get_sub_periods() + for p in glob.glob( + os.path.join(base_dir, era, variable, f"{variable}.root") + ) + ] + if not inputs: + continue + summed = os.path.join(plots_dir, f"_summed_{variable}.root") + try: + ps_call(["hadd", "-f", summed] + inputs, verbose=1) + except Exception as e: + print( + f"WARNING: could not sum the shapes for {variable}, " + f"datacards are unaffected: {e}" + ) + continue + + # Only the categories this mass point actually has. The configuration lists + # every datacard bin, but the preprocessing gates drop a category where the + # signal is too small -- boosted at a low mass, say -- and it is then absent + # from the shapes. HistPlotter exits non-zero on the first key it cannot find, + # so passing the full list costs every panel and grid that would have come + # after the gap, not just the missing one. + present = self._present_keys(summed, cfg) + keys, outputs = [], [] + for category in cfg["categories"]: + # "SR/res2b_dnn0" -> region "SR", category "res2b_dnn0": HistPlotter + # navigates channel -> region -> category. + region, _, cat = category.rpartition("/") + for channel in cfg["channels"]: + if (channel, region, cat) not in present: + continue + keys.append(f"{channel}:{cat}:{region}") + outputs.append( + os.path.join( + plots_dir, f"{variable}_{channel}_{region}_{cat}.pdf" + ) + ) + if not keys: + print(f"WARNING: no shapes to plot for {variable}") + continue + cmd = [ + "python3", + plotter, + "--inFile", + summed, + "--all_outFiles", + ",".join(outputs), + "--all_keys", + ",".join(keys), + "--globalConfig", + os.path.join( + self.ana_path(), + self.global_params["analysis_config_area"], + "global.yaml", + ), + # The surviving axis of the rebinned shapes, and already at its final + # binning -- so no --rebin, which would coarsen it back to the histograms.yaml + # grid and undo the whole point of the rebinning. + "--var", + self.plot_variable(variable), + "--year", + self.period, + "--ana_path", + self.ana_path(), + "--period", + self.period, + "--LAWrunVersion", + self.version, + "--wantSignals", + # The slices span ~5 decades in yield (the low-significance slice holds + # most of the background), so a linear axis hides everything but slice 0. + "--wantLogScale", + "y", + ] + try: + ps_call(cmd, verbose=1) + self._stitch_grids(cfg, plots_dir, variable, dict(zip(keys, outputs))) + except Exception as e: + # The datacards are the product that matters; a plotting failure should be + # loud but must not leave the task looking broken with valid cards on disk. + print( + f"WARNING: shape plotting failed for {variable}, " + f"datacards are unaffected: {e}" + ) + finally: + if os.path.exists(summed): + os.remove(summed) diff --git a/law/DhiPlotMixin.py b/law/DhiPlotMixin.py new file mode 100644 index 0000000..b0d6388 --- /dev/null +++ b/law/DhiPlotMixin.py @@ -0,0 +1,122 @@ +import law +import luigi +import os +import shutil + +from FLAF.RunKit.run_tools import ps_call + + +class DhiPlotMixin: + """Running a dhi plot task from one of ours, and keeping what it drew. + + dhi writes its plots under inference/data/store, which is what makes an + already-drawn plot cheap to re-request but leaves the products somewhere other than + the rest of the chain's. Every task here follows the same three steps: work out where + dhi will write, run it as a subprocess, copy the result into our own output. + + The dhi tasks are invoked as subprocesses rather than yielded as dynamic dependencies. + luigi round-trips a dynamic dependency through to_str_params() / from_str_params(), + and dhi declares parameters whose unset defaults do not survive that -- lumi_scale is + a FloatParameter(default=None) that serialises to "None" and then raises on + float("None"). Going through the command line keeps dhi unmodified. + """ + + # Forwarded as "--remove-output 0,a,y": depth 0 drops the plot itself and nothing + # below it, so no fit is recomputed. Needed because much of the plot styling is + # significant=False and so does not move the output path -- changing it would + # otherwise leave the old plot in place and look like the change had no effect. + redraw = luigi.BoolParameter( + default=False, + significant=False, + description="discard existing plots and draw them again", + ) + + @staticmethod + def known_plot_params(task_cls): + """The parameter names a dhi plot task accepts. + + get_params() rather than get_param_names(), which drops the insignificant ones -- + that is most of the plot styling (x_min, campaign, parameters_per_page, ...). + """ + return {param_name for param_name, _ in task_cls.get_params()} + + def validate_plot_params(self, task_cls, params, context): + """Reject a plot_params key the dhi task would silently ignore.""" + known = self.known_plot_params(task_cls) + for key in params: + if key not in known: + raise RuntimeError( + f"{context}: '{key}' is not a {task_cls.__name__} parameter." + ) + + def plot_targets(self, task_cls, spec): + """Where the dhi task will write, without scheduling it. + + Constructing the task is safe and is the only honest way to learn its output + paths -- they are built from a hash of the datacards plus several parameters, + which is not something to reimplement here. flatten() because a dhi task's + output() may be a single target, a list, or a dict of both (PlotPullsAndImpacts + returns {"plots": [...], "plot_data": ...}). + """ + return law.util.flatten(task_cls(**spec).output()) + + def plot_command(self, task_cls, spec, extra_args=()): + """`law run ...` for a spec. + + ``extra_args`` is appended verbatim. It exists for the arguments that belong to a + task upstream of the one being run rather than to it -- law namespaces those as + ---, so they are not parameters of task_cls and cannot go in the + spec. Sending the per-parameter fits of PullsAndImpacts to HTCondor is the case + that needs it: --PullsAndImpacts-workflow htcondor. + """ + cmd = ["law", "run", task_cls.__name__] + for key, value in spec.items(): + flag = "--" + key.replace("_", "-") + if isinstance(value, bool): + # luigi bool parameters are set by presence, not by value + if value: + cmd.append(flag) + elif key == "multi_datacards": + # colon between datacard sequences, comma within one + cmd += [flag, ":".join(",".join(seq) for seq in value)] + elif isinstance(value, (tuple, list)) and all( + isinstance(v, (tuple, list)) and len(v) == 2 for v in value + ): + # dhi's MultiCSVParameter pairs, e.g. parameter_values ("r", 6.8) -> r=6.8 + cmd += [flag, ",".join(f"{a}={b}" for a, b in value)] + elif isinstance(value, (tuple, list)): + cmd += [flag, ",".join(str(v) for v in value)] + else: + cmd += [flag, str(value)] + cmd += [str(a) for a in extra_args] + if self.redraw: + cmd += ["--remove-output", "0,a,y"] + return cmd + + def draw_plot(self, task_cls, spec, name, era, dest_dir, extra_args=()): + """Run a dhi plot task, check it actually drew, copy the result into ``dest_dir``. + + Returns the basenames copied, for the manifest. + """ + # cwd is pinned to the analysis root: dhi's resolve_datacards() takes a + # different branch when the process happens to sit inside a configured + # datacards_run2 directory. + ps_call( + self.plot_command(task_cls, spec, extra_args), + cwd=self.ana_path(), + verbose=1, + ) + + basenames = [] + for target in self.plot_targets(task_cls, spec): + # luigi's retcode defaults return 0 even when a task fails, and no + # [retcode] section overrides them here -- so a clean exit is not + # evidence that anything was drawn. The file is. + if not os.path.exists(target.path): + raise RuntimeError( + f"'{name}' ({era}): law exited cleanly but {target.path} was not " + "written. See the output above for the task that actually failed." + ) + shutil.copy2(target.path, dest_dir) + basenames.append(os.path.basename(target.path)) + return basenames diff --git a/law/MergedHists.py b/law/MergedHists.py new file mode 100644 index 0000000..41a5709 --- /dev/null +++ b/law/MergedHists.py @@ -0,0 +1,30 @@ +import law +import luigi + +from .StatInferenceTask import StatInferenceTask + + +class MergedHists(StatInferenceTask, law.ExternalTask): + """One merged histogram file, addressed by path under --hists-version. + + External on purpose. The limit-setting chain consumes histograms it does not make, so + a missing input is reported as a missing dependency rather than pulling the whole + AnaTuple production graph (AnaTupleFileListTask -> HistFromNtupleProducerTask -> + HistMergerTask) into the run and rebuilding its branch maps just to discover the file + is already there. + + The path mirrors HistMergerTask.output() in FLAF/Analysis/tasks.py -- that is the + contract between the two halves, and the only thing this chain needs from FLAF. + """ + + variable = luigi.Parameter(description="Hists_merged variable (directory) name") + + def output(self): + return self.remote_target( + self.input_hists_version, + "Hists_merged", + self.period, + self.variable, + f"{self.variable}.root", + fs=self.fs_HistTuple, + ) diff --git a/law/PlotPullsAndImpactsTask.py b/law/PlotPullsAndImpactsTask.py new file mode 100644 index 0000000..4d95484 --- /dev/null +++ b/law/PlotPullsAndImpactsTask.py @@ -0,0 +1,314 @@ +import glob +import json +import law +import luigi +import os +import re + +from string import Template + +from FLAF.RunKit.run_tools import ps_call +from dhi.tasks.pulls_impacts import MergePullsAndImpacts, PlotPullsAndImpacts +from dhi.tasks.resonant import MergeResonantLimits + +from .DhiPlotMixin import DhiPlotMixin +from .StatInferenceTask import StatInferenceTask +from .ResonantLimitsTask import ResonantLimitsTask + + +class PlotPullsAndImpactsTask(DhiPlotMixin, StatInferenceTask): + """The nuisance pull and impact plots, declared in the datacard configuration. + + Each entry of the config's ``impact_plots`` block becomes one dhi PlotPullsAndImpacts + task per mass point it lists. Like the limit plots, the results are copied to + fs_default beside the rest of the chain's products; dhi keeps its own copy under + inference/data/store, which is what makes a redraw cheap. + + The input is the combined card ResonantLimitsTask already writes -- one per mass, + every era in it, at ``//Datacards/combined/combined_.txt``. Their + shape paths are absolute, which is what makes them usable here at all: PullsAndImpacts + runs combine in a temporary directory. + + This is a task of its own rather than something the limit chain drags in because it is + far more expensive than everything before it. PullsAndImpacts is a per-parameter + workflow: roughly two combine fits per nuisance per mass point, so a card with 76 + nuisances costs ~150 fits for a single mass. + """ + + workflow = luigi.Parameter(default=law.parameter.NO_STR) + + def requires(self): + # ResonantLimitsTask is what writes the combined cards this reads. + return [ResonantLimitsTask.req(self)] + + def output(self): + return self.output_dir_target(self.version, "ImpactPlots") + + def entry_eras(self, entry): + """The eras an entry is drawn for. + + An entry with a `glob:` addresses cards inside one era's directory, so it is + drawn once per era the configuration lists. The default combined card already + contains every era, so it is drawn once and labelled "combined". + """ + return self.get_all_eras() if entry.get("glob") else ["combined"] + + def card_for(self, entry, era, mass): + """The single datacard an entry's plot is built from. + + PullsAndImpacts fits one card, not a glob -- every parameter gets its own + combine job against that one workspace -- so an entry's `glob:` and mass have to + resolve to exactly one file. + """ + name = entry.get("name") or "impacts" + pattern = entry.get("glob") + if not pattern: + # The combined card ResonantLimitsTask writes: one per mass, every era in it. + base = self.datacards_dir("combined") + candidates = glob.glob(os.path.join(base, "combined_*.txt")) + else: + base = self.datacards_dir(era) + candidates = glob.glob( + os.path.join(base, Template(pattern).safe_substitute(ERA=era)) + ) + + # Same convention dhi's datacard_pattern uses: the mass is the trailing number. + matched = [ + p + for p in candidates + if (m := re.search(r"_(\d+)\.txt$", os.path.basename(p))) + and int(m.group(1)) == int(mass) + ] + if len(matched) != 1: + found = sorted(os.path.basename(p) for p in candidates) + raise RuntimeError( + f"impact_plots entry '{name}' ({era}, mass {mass}): expected exactly one " + f"datacard, found {len(matched)}. Looked in {base} for " + f"{pattern or 'combined_*.txt'}; it holds: {', '.join(found) or 'nothing'}." + ) + return matched[0] + + def card_family(self, entry, era): + """The datacard set an entry's card belongs to, as a MergeResonantLimits glob. + + The limit is merged over a family of per-mass cards, so the family is what has to + be named to look one up -- the same glob the entry selects its own card from, + without the mass filter. + """ + pattern = entry.get("glob") + if pattern: + return ( + os.path.join( + self.datacards_dir(era), Template(pattern).safe_substitute(ERA=era) + ), + ) + + # The combined card's own family, combined_*.txt, is a datacard set nobody has + # limits for -- asking for them starts a fresh workspace and fit per mass. The + # limits that do exist are the ones the limit plots merged, whose "Combined" curve + # globs an era's cards; with a single top-level era that is the same measurement + # the combined card makes. + eras = self.get_top_level_eras() + if len(eras) != 1: + raise RuntimeError( + f"impact_plots: poi_value 'limit' needs one top-level era to take the " + f"limit from, but the configuration has {eras}. Give the entry a number " + "instead." + ) + return (os.path.join(self.datacards_dir(eras[0]), "*.txt"),) + + def expected_limit(self, entry, era, mass): + """The expected limit on r at a mass, from the limits already computed for it. + + The path comes from dhi's own task rather than being assembled here: the store + directory is a hash of the resolved datacard set. MergeResonantLimits is run if it + has not been already -- ResonantLimitsTask has produced the per-mass limits this + merges, so this is a merge and not a fit. + """ + import numpy as np + + sequence = self.card_family(entry, era) + target = MergeResonantLimits( + version=self.version, datacards=tuple(sequence) + ).output() + if not os.path.exists(target.path): + ps_call( + [ + "law", + "run", + "MergeResonantLimits", + "--version", + self.version, + "--datacards", + ",".join(sequence), + ], + cwd=self.ana_path(), + verbose=1, + ) + if not os.path.exists(target.path): + raise RuntimeError( + f"impact_plots: no merged limits at {target.path} for {list(sequence)}, " + "so poi_value: limit cannot be resolved. Set a number instead." + ) + data = np.load(target.path, allow_pickle=True)["data"] + for row in data: + if int(row["mhh"]) == int(mass): + return float(row["limit"]) + raise RuntimeError( + f"impact_plots: the merged limits at {target.path} carry no mass {mass}; " + f"they hold {sorted({int(r['mhh']) for r in data})}." + ) + + # The parameter the ranking is of. dhi's own default for a resonant search. + poi_name = "r" + + def poi_value(self, entry, era, mass): + """The value of r the Asimov dataset is built at, or None to leave dhi's default. + + This matters more than it looks. dhi builds the Asimov at r=1 unless told + otherwise, and r=1 is not a meaningful reference at either end of this scan: where + the limit is ~7 it is an almost invisible signal, so the ranking is really the + background-only one, and where the limit is ~0.16 it is several times more signal + than could be excluded, which pins r and collapses every impact towards zero. + Ranking at the limit puts the fit where the measurement actually is. + + Note this cannot go through dhi's `parameter_values`: for a resonant search + hh_model is NO_STR, and POITask then hard-codes both the joined parameter values + ('""') and the output postfix (r=1.0), so the value would be silently dropped. It + has to reach combine as --expectSignal through PullsAndImpacts' custom_args, which + is significant and so does change the fit's output path. + """ + value = entry.get("poi_value") + if value is None: + return None + if value == "limit": + return self.expected_limit(entry, era, mass) + return float(value) + + def build_plot_spec(self, entry, era, mass): + """PlotPullsAndImpacts parameters for an entry and mass point. + + Kept as a plain dict so the same values can both construct the task (to learn + where it will write) and be rendered into a command line (to make it write). + """ + name = entry.get("name") or "impacts" + plot_params = dict(entry.get("plot_params") or {}) + self.validate_plot_params( + PlotPullsAndImpacts, plot_params, f"impact_plots entry '{name}'" + ) + + # The combined cards carry a few hundred autoMCStats bins, so mc_stats puts ~460 + # parameters on the plot. parameters_per_page defaults to -1, meaning a single + # page, and the result is an unreadable hairline strip rather than an error. + if ( + plot_params.get("mc_stats") + and plot_params.get("parameters_per_page", -1) < 1 + ): + raise RuntimeError( + f"impact_plots entry '{name}': mc_stats puts every autoMCStats bin on the " + "plot, and parameters_per_page defaults to -1 (one page), which is " + "unreadable. Set parameters_per_page (25 works) or drop mc_stats." + ) + + return dict( + version=self.version, + datacards=(self.card_for(entry, era, mass),), + mass=float(mass), + # PullsAndImpacts is a plain POITask, so hh_model defaults to the + # non-resonant model_default -- which would fit r together with kl, kt, CV and + # C2V, and build its own workspace under a hh_model__model_default/ store path + # rather than reusing the one the resonant chain already made. dhi's own + # resonant tasks pin hh_model to NO_STR for the same reason + # (dhi/tasks/resonant.py:35); allow_empty_hh_model is already True on the + # PullsAndImpacts base, so this is simply saying the analysis is resonant. + # An entry may still override it through plot_params. + **{"hh_model": law.NO_STR, **plot_params}, + ) + + def report_dropped_parameters(self, entry, era, mass, spec): + """Warn about nuisances robustHesse removed from the fit. + + robustHesse drops a parameter it cannot invert, logs "Dropping from the + hessian" and carries on successfully; the nuisance is then simply absent from the + plot and the merged JSON, with nothing on the plot marking its absence. So the + card's own nuisance list is the reference: anything in it that the merged JSON + does not carry was dropped. + """ + merged = MergePullsAndImpacts( + **{ + k: v + for k, v in spec.items() + if k in self.known_plot_params(MergePullsAndImpacts) + } + ).output() + if not os.path.exists(merged.path): + return + with open(merged.path) as f: + fitted = {p["name"] for p in json.load(f).get("params", [])} + + card_params = set() + with open(self.card_for(entry, era, mass)) as f: + for line in f: + parts = line.split() + if len(parts) >= 2 and parts[1] in ("shape", "lnN"): + card_params.add(parts[0]) + + dropped = sorted(card_params - fitted) + if dropped: + name = entry.get("name") or "impacts" + print( + f"WARNING: impact_plots entry '{name}' ({era}, mass {mass}): " + f"{len(dropped)} " + f"nuisance(s) are in the card but not in the fit, so they are missing " + f"from the plot without being marked: {', '.join(dropped)}. " + "With --method robust this is robustHesse dropping what it could not " + "invert; the ranking is not complete." + ) + + def run(self): + entries = self.get_config_data().get("impact_plots", []) + if not entries: + raise RuntimeError( + f"{self.datacard_config_path()} declares no 'impact_plots' block, so " + "there is nothing to plot." + ) + + with self.output().localize("w") as local_output: + produced = [] + for entry in entries: + name = entry.get("name") or "impacts" + masses = entry.get("masses") + if not masses: + raise RuntimeError( + f"impact_plots entry '{name}': no 'masses'. Each one costs about " + "two combine fits per nuisance, so they are listed explicitly " + "rather than defaulting to every mass in the model." + ) + + for era in self.entry_eras(entry): + for mass in masses: + rel = os.path.join(name, era, str(mass)) + dest_dir = os.path.join(local_output.abspath, rel) + os.makedirs(dest_dir, exist_ok=True) + + spec = self.build_plot_spec(entry, era, mass) + extra = list(entry.get("dhi_args") or ()) + poi = self.poi_value(entry, era, mass) + if poi is not None: + extra.append( + "--PullsAndImpacts-custom-args=" + f"--expectSignal={poi:.4g}" + ) + basenames = self.draw_plot( + PlotPullsAndImpacts, + spec, + name, + f"{era} mass {mass}", + dest_dir, + extra_args=extra, + ) + self.report_dropped_parameters(entry, era, mass, spec) + produced.extend(os.path.join(rel, b) for b in basenames) + + with open(os.path.join(local_output.abspath, "impacts.json"), "w") as f: + json.dump(sorted(produced), f, indent=2) diff --git a/law/PlotResonantLimitsTask.py b/law/PlotResonantLimitsTask.py new file mode 100644 index 0000000..c706db7 --- /dev/null +++ b/law/PlotResonantLimitsTask.py @@ -0,0 +1,320 @@ +import glob +import json +import law +import luigi +import math +import os +import re + +from string import Template +from FLAF.RunKit.run_tools import ps_call +from dhi.config import campaign_labels, campaign_lumis +from dhi.tasks.resonant import ( + MergeResonantLimits, + PlotMultipleResonantLimits, + PlotResonantLimits, +) + +from .DhiPlotMixin import DhiPlotMixin +from .StatInferenceTask import StatInferenceTask +from .ResonantLimitsTask import ResonantLimitsTask + + +class PlotResonantLimitsTask(DhiPlotMixin, StatInferenceTask): + """The limit plots, declared in the datacard configuration. + + Each entry of the config's ``limit_plots`` block becomes one dhi + PlotMultipleResonantLimits task: a list of datacard globs with the labels they should + carry, optional external limit curves, and the plot styling. Globs are relative to the + era's datacards directory and may use ``${ERA}``; the per-channel and per-category + sub-directories they address come from DatacardMaker.writeDatacards. + + An entry with ``bands: true`` additionally gets one dhi PlotResonantLimits plot per + curve -- the standard single-curve plot with the +-1/+-2 sigma bands, which the overlay + cannot show. External limits and luminosity projections are separate curves rather than + datacards, so they appear only on the overlay. + + The finished plots are copied to fs_default alongside the rest of the chain's + products; see DhiPlotMixin for how a dhi plot task is run and its output collected. + """ + + workflow = luigi.Parameter(default=law.parameter.NO_STR) + + def requires(self): + # ResonantLimitsTask is what mirrors the cards to datacards_dir(), which is what + # the globs below resolve against -- the plots cannot be built before it has run. + return [ResonantLimitsTask.req(self)] + + def output(self): + # fs_default, like CreateDatacardsTask. Holds one + # / sub-directory of plots plus plots.json naming them. + return self.output_dir_target(self.version, "LimitPlots") + + def _resolve_config_path(self, path): + return path if os.path.isabs(path) else os.path.join(self.ana_path(), path) + + def common_plot_params(self, entry, era): + """The dhi parameters an entry's plots share, whatever kind they are. + + Validated against PlotMultipleResonantLimits, which is the subclass and so carries + every parameter the single-curve PlotResonantLimits has plus its own -- one check + covers both plot kinds, and the band builder drops what does not apply. + """ + name = entry.get("name") or "limits" + + params = { + # Distinguishes the plot files of entries that differ only in labelling; + # not significant, so it does not move the output directory. + "plot_postfix": name, + # Per era, not per entry: the same entry is built for every top-level era and + # each one carries a different luminosity. plot_params may still override it. + "campaign": self.get_campaign(era), + } + plot_params = entry.get("plot_params") or {} + self.validate_plot_params( + PlotMultipleResonantLimits, plot_params, f"limit_plots entry '{name}'" + ) + params.update(plot_params) + + # dhi turns an unset campaign into None and simply draws no luminosity label, so + # a missing entry would silently produce an unlabelled plot. Refuse instead. + if not params.get("campaign"): + raise RuntimeError( + f"limit_plots entry '{name}': no campaign for era '{era}'. Add it to the " + f"'campaigns' map in {self.datacard_config_path()} (a key of " + "dhi.config.campaign_labels), or set it in this entry's plot_params." + ) + if params["campaign"] not in campaign_labels: + # dhi falls back to drawing the key itself, which is a usable escape hatch for + # a one-off label but is also exactly what a typo looks like. + print( + f"WARNING: limit_plots entry '{name}': campaign '{params['campaign']}' is " + "not in dhi.config.campaign_labels; it will be drawn verbatim as the " + "luminosity label." + ) + + return params + + def build_plot_spec(self, entry, era, extra_external=()): + """PlotMultipleResonantLimits parameters for a ``limit_plots`` entry and era. + + Kept as a plain dict so the same values can both construct the task (to learn + where it will write) and be rendered into a command line (to make it write). + + ``extra_external`` holds external-limit files generated for this run (the + luminosity projections). They come first, so a projection of our own curve is + drawn before any fixed reference such as the Run 2 result. + """ + sequences, names = self.datacard_sequences(entry, era) + external = tuple(extra_external) + tuple( + self._resolve_config_path(p) for p in entry.get("external_limits") or [] + ) + + return dict( + version=self.version, + multi_datacards=tuple(sequences), + datacard_names=tuple(names), + external_limits=external, + **self.common_plot_params(entry, era), + ) + + def build_band_specs(self, entry, era): + """PlotResonantLimits parameters, one per curve, for an entry with ``bands: true``. + + The band plot is drawn from a single datacard set, so the multi-curve parameters + (multi_datacards, datacard_names, external_limits, colors, markers) have no meaning + here and are dropped by filtering against the parameters PlotResonantLimits + actually declares. Everything else -- axes, xsec, campaign -- is shared with the + overlay by construction, so the two plots of an entry cannot drift apart. + """ + if not entry.get("bands"): + return [] + + sequences, names = self.datacard_sequences(entry, era) + shared = self.common_plot_params(entry, era) + known = {param_name for param_name, _ in PlotResonantLimits.get_params()} + + specs = [] + for sequence, label in zip(sequences, names): + spec = {key: value for key, value in shared.items() if key in known} + spec.update( + version=self.version, + datacards=tuple(sequence), + # Entry name and curve label both. All three entries of the HH->bbWW + # config draw "*.txt", so dhi would hash them into one output directory; + # the postfix is what keeps their band plots from overwriting each other. + # dhi's join_postfix strips whatever is not [a-zA-Z0-9._+-], so a label + # like "e#mu" needs no sanitising here. + plot_postfix=f"{shared['plot_postfix']}_{label}", + ) + specs.append(spec) + return specs + + def datacard_sequences(self, entry, era): + """(datacard glob per curve, legend label per curve) for an entry, validated. + + Shared by the plot itself and by the luminosity projection, which has to scale the + same cards the leading curve is drawn from. + """ + name = entry.get("name") or "limits" + base_dir = self.datacards_dir(era) + + sequences, names = [], [] + for card in entry["datacards"]: + pattern = os.path.join( + base_dir, Template(card["glob"]).safe_substitute(ERA=era) + ) + if not glob.glob(pattern): + raise RuntimeError( + f"limit_plots entry '{name}': no datacard matches '{pattern}'. " + f"Check the glob against the cards under {base_dir}." + ) + label = card["name"] + if "{" in label or "}" in label: + # law brace-expands datacard_names, so "fb^{-1}" would arrive as "fb^-1". + raise RuntimeError( + f"limit_plots entry '{name}': datacard name '{label}' contains braces, " + "which law brace-expands. Put the information in the campaign label instead." + ) + sequences.append((pattern,)) + names.append(label) + return sequences, names + + def limits_npz(self, sequence): + """The MergeResonantLimits .npz for a datacard sequence, produced if absent. + + The path is taken from dhi's own task rather than assembled here: the store + directory is a hash of the resolved datacard set, which is not something to + reimplement. Passing the glob is equivalent to passing the resolved paths -- + dhi resolves it in modify_param_values before hashing. + """ + target = MergeResonantLimits( + version=self.version, datacards=tuple(sequence) + ).output() + if not os.path.exists(target.path): + ps_call( + [ + "law", + "run", + "MergeResonantLimits", + "--version", + self.version, + "--datacards", + ",".join(sequence), + ], + cwd=self.ana_path(), + verbose=1, + ) + if not os.path.exists(target.path): + raise RuntimeError( + f"MergeResonantLimits produced no limits for {list(sequence)}; expected " + f"{target.path}." + ) + return target.path + + def make_lumi_projection(self, entry, era, projection, sequence, out_dir): + """Write a dhi external-limits JSON scaling this entry's own limit curve to a + different integrated luminosity, and return its path. + + Regenerated from the current limits on every run rather than read from a + checked-in file: the projection is a function of the measured curve, so a stored + copy silently goes stale the moment the binning, the datacards or the fits change. + """ + name = entry.get("name") or "limits" + label = projection["name"] + campaign = self.get_campaign(era) + # Defaults to the luminosity dhi prints on the plot for this campaign, so the + # curve cannot be scaled from a different number than the one drawn beside it. + lumi_now = projection.get("lumi_now", campaign_lumis.get(campaign)) + if not lumi_now: + raise RuntimeError( + f"limit_plots entry '{name}': lumi_projection '{label}' has no lumi_now " + f"and campaign '{campaign}' is not in dhi.config.campaign_lumis." + ) + + out_path = os.path.join( + out_dir, re.sub(r"[^\w.-]+", "_", f"{name}_{label}").strip("_") + ".json" + ) + + import numpy as np + + # The same purely statistical scaling dhi's own --lumi-scale applies + # (plots/limits.py): limit(L_target) = limit(L_now) * sqrt(L_now / L_target). + # It assumes the systematics scale away with the data, which they do not, so + # this curve is optimistic -- an extrapolation of the statistical reach, not a + # projected result. --lumi-scale itself is not usable here: it exists only on + # plot_limit_scan, and it overwrites the measured curve rather than adding a + # second one, whereas the point is to draw both side by side. + # + # Only "factor" is written, not pre-scaled limits: dhi's read_external_limits + # (tasks/resonant.py) multiplies by factor * scale itself, so applying it here + # too would double-count it. + data = np.load(self.limits_npz(sequence), allow_pickle=True)["data"] + entry_json = [ + { + "name": label, + "scan_parameter": "mhh", + "factor": math.sqrt(lumi_now / projection["lumi_target"]), + "limits": {repr(float(r["mhh"])): float(r["limit"]) for r in data}, + } + ] + with open(out_path, "w") as f: + json.dump(entry_json, f, indent=2) + f.write("\n") + return out_path + + def run(self): + entries = self.get_config_data().get("limit_plots", []) + if not entries: + raise RuntimeError( + f"{self.datacard_config_path()} declares no 'limit_plots' block, so there " + "is nothing to plot." + ) + + with self.output().localize("w") as local_output: + produced = [] + for era in self.get_eras_to_plot(): + dest_dir = os.path.join(local_output.abspath, era) + os.makedirs(dest_dir, exist_ok=True) + + for entry in entries: + name = entry.get("name") or "limits" + + # Projections are written into the output directory before the plot + # runs: they are --external-limits inputs, and they are worth keeping + # next to the plot as the record of what was actually drawn. + sequence = self.datacard_sequences(entry, era)[0][0] + projections = [ + self.make_lumi_projection(entry, era, proj, sequence, dest_dir) + for proj in entry.get("lumi_projections") or [] + ] + + specs = [ + ( + PlotMultipleResonantLimits, + self.build_plot_spec( + entry, era, extra_external=projections + ), + ) + ] + specs += [ + (PlotResonantLimits, s) + for s in self.build_band_specs(entry, era) + ] + + for task_cls, spec in specs: + produced.extend( + os.path.join(era, b) + for b in self.draw_plot(task_cls, spec, name, era, dest_dir) + ) + produced.extend( + os.path.join(era, os.path.basename(p)) for p in projections + ) + + with open(os.path.join(local_output.abspath, "plots.json"), "w") as f: + json.dump(sorted(produced), f, indent=2) + + def get_eras_to_plot(self): + """Every era the configuration lists -- each has its own limit, so each gets its + own plots.""" + return self.get_all_eras() diff --git a/law/PreprocessShapesTask.py b/law/PreprocessShapesTask.py new file mode 100644 index 0000000..441137b --- /dev/null +++ b/law/PreprocessShapesTask.py @@ -0,0 +1,106 @@ +import contextlib +import law +import luigi +import os + +from string import Template + +from FLAF.RunKit.run_tools import ps_call +from FLAF.run_tools.law_customizations import HTCondorWorkflow, copy_param + +from .StatInferenceTask import StatInferenceTask + + +class PreprocessShapesTask(StatInferenceTask, HTCondorWorkflow, law.LocalWorkflow): + """Run the datacard configuration's `preprocess:` step over the merged histograms. + + A hook, not a rebinning. The configuration names a script and any arguments it wants; + this task supplies --input, --output, --era and --config and knows nothing else about + what the step does. --config is the datacard configuration, which any step working on + these shapes needs anyway: it is where the processes, the model and the categories + are. HH->bbWW plugs in bin_opt_2d/rebin_2d.py to cut its 2D DNN-vs-HME shapes + into per-slice 1D ones, but an analysis that needs some other transformation writes its + own script, and one that needs none declares no `preprocess:` block at all -- then this + task is never scheduled and CreateDatacardsTask reads the merged histograms directly. + + Whatever the step writes must be laid out as "//.root", which + is what the model's input_file_pattern resolves against and what the merged histograms + already look like. That is the whole contract: the datacard step cannot tell whether it + is reading preprocessed shapes or raw ones. + """ + + max_runtime = copy_param(HTCondorWorkflow.max_runtime, 4.0) + n_cpus = copy_param(HTCondorWorkflow.n_cpus, 1) + + # As on CreateDatacardsTask: preprocessing a meta-era needs its members, and + # self.period must stay a real era because FLAF's Setup resolves it against config/. + meta_era = luigi.Parameter(default="") + + @property + def datacard_era(self): + return self.meta_era or self.period + + def get_sub_periods(self): + return self.get_era_groups().get(self.datacard_era, [self.period]) + + def input_hist_reqs(self): + return { + f"MergedHists_{era}_{variable}": req + for (era, variable), req in self.merged_hist_reqs( + self.get_sub_periods() + ).items() + } + + def workflow_requires(self): + return self.input_hist_reqs() + + def requires(self): + return list(self.input_hist_reqs().values()) + + def create_branch_map(self): + return {0: None} + + def output(self): + return self.output_dir_target( + self.version, "Hists_preprocessed", self.datacard_era + ) + + def run(self): + cfg = self.preprocess_config() + if not cfg: + raise RuntimeError( + f"{self.datacard_config_path()} declares no 'preprocess' block, so this " + "task has nothing to run. CreateDatacardsTask reads the merged histograms " + "directly in that case and does not require it." + ) + script = os.path.join(self.ana_path(), cfg["script"]) + with contextlib.ExitStack() as stack: + targets = { + key: req.output() + for key, req in self.merged_hist_reqs(self.get_sub_periods()).items() + } + self.check_inputs(targets) + base_dir_local = stack.enter_context(self.stage_inputs(targets)) + local_output = stack.enter_context(self.output().localize("w")) + cmd = [ + "python3", + "-u", + script, + "--input", + base_dir_local.abspath, + "--output", + local_output.abspath, + "--era", + self.datacard_era, + "--config", + self.datacard_config_path(), + ] + # ${ERA} in an argument names the era being produced, the same way the model's + # input_file_pattern does. Relative paths are resolved against the analysis + # area so a configuration can point at its own files without absolute paths. + for arg in cfg.get("args", []): + arg = Template(str(arg)).safe_substitute(ERA=self.datacard_era) + if arg.startswith("config/"): + arg = os.path.join(self.ana_path(), arg) + cmd.append(arg) + ps_call(cmd, env=self.cmssw_env, verbose=1) diff --git a/law/ResonantLimitsAndHistPlotTask.py b/law/ResonantLimitsAndHistPlotTask.py new file mode 100644 index 0000000..4903513 --- /dev/null +++ b/law/ResonantLimitsAndHistPlotTask.py @@ -0,0 +1,45 @@ +import law +import luigi + +from FLAF.Analysis.tasks import HistPlotTask + +from .StatInferenceTask import StatInferenceTask +from .ResonantLimitsTask import ResonantLimitsTask + + +class ResonantLimitsAndHistPlotTask(StatInferenceTask): + """The limits, and the plots of the histograms that went into them, from one + `law run`. + + Depends on ResonantLimitsTask, not PlotResonantLimitsTask: this task is about getting + the limits together with their input histograms, and the limit plots the configuration + declares in `limit_plots` are asked for separately, by running PlotResonantLimitsTask. + """ + + workflow = luigi.Parameter(default=law.parameter.NO_STR) + + def get_plottable_eras(self): + """Real production periods, for FLAF's HistPlotTask. + + Deliberately not the same set as ResonantLimitsTask.get_eras(), which returns the + eras that carry their own datacards -- for a configuration with an `era_groups:` + entry that is the meta-era, and this is its members. The two are complements: + limits are set on the combination, while the input histograms can only be plotted + for periods that really exist, since HistPlotTask resolves `period` against + config//. + """ + era_groups = self.get_era_groups() + eras = self.get_config_data().get("eras", [self.period]) + return [e for e in eras if e not in era_groups] + + def requires(self): + reqs = [ResonantLimitsTask.req(self)] + for e in self.get_plottable_eras(): + reqs.append(HistPlotTask.req(self, period=e)) + return reqs + + def output(self): + return self.local_target("dummy.txt") + + def run(self): + self.output().touch() diff --git a/law/ResonantLimitsTask.py b/law/ResonantLimitsTask.py new file mode 100644 index 0000000..3389a22 --- /dev/null +++ b/law/ResonantLimitsTask.py @@ -0,0 +1,110 @@ +import glob +import law +import luigi +import os +import re +import shutil +import subprocess + +from dhi.tasks.resonant import MergeResonantLimits + +from .StatInferenceTask import StatInferenceTask +from .CreateDatacardsTask import CreateDatacardsTask + + +class ResonantLimitsTask(StatInferenceTask): + workflow = luigi.Parameter(default=law.parameter.NO_STR) + + def store_parts(self): + return (self.version, self.__class__.__name__, "combined") + + def get_eras(self): + """Every era the configuration lists: each gets its own cards and its own limit.""" + return self.get_all_eras() + + def _create_datacards_req(self, era, **kwargs): + era_groups = self.get_era_groups() + if era in era_groups: + return CreateDatacardsTask.req( + self, period=era_groups[era][0], meta_era=era, **kwargs + ) + return CreateDatacardsTask.req(self, period=era, **kwargs) + + def requires(self): + return [self._create_datacards_req(e, branches=()) for e in self.get_eras()] + + def output(self): + return { + "limits": self.local_target("limits.npz"), + "datacards": law.LocalDirectoryTarget(self.datacards_dir("combined")), + } + + def stage_datacards(self, era, remote_target): + """Mirror an era's datacards from fs_default to a stable local path, returned. + + CreateDatacardsTask writes to fs_default, but the cards must be real local files + by the time combine sees them: MergeResonantLimits shells out to combine, and + PlotMultipleResonantLimits resolves --multi-datacards by globbing the filesystem + outside law entirely. A localize() scratch dir will not do -- run() yields, which + suspends and re-enters it, so the scratch dir is gone before dhi reads anything. + + The mirror path is the one CreateDatacardsTask used to write to directly, so + existing --multi-datacards invocations keep working unchanged. + """ + local_dir = self.datacards_dir(era) + local_target = law.LocalDirectoryTarget(local_dir) + # Re-staged on every run, including the re-entry after the yield below. The + # copy is small (~30 MB) and unconditional refresh is what keeps the mirror + # from going stale when CreateDatacardsTask reruns within the same version. + if local_target.exists(): + local_target.remove() + local_target.touch() + remote_target.copy_to_local(local_target) + return local_dir + + def run(self): + datacards = [] + eras = self.get_eras() + era_cards = {} + + for e in eras: + create_dc_br0 = self._create_datacards_req(e, branch=0, branches=()) + output_dir = self.stage_datacards(e, create_dc_br0.output()) + cards = glob.glob(os.path.join(output_dir, "*.txt")) + era_cards[e] = cards + datacards.extend(cards) + + limits = yield MergeResonantLimits( + version=self.version, datacards=tuple(datacards) + ) + print(f"Merged limits: {limits}") + + self.output()["limits"].parent.touch() + shutil.copy2(limits.path, self.output()["limits"].path) + + out_dc_dir = self.output()["datacards"] + out_dc_dir.touch() + + masses = set() + for e, cards in era_cards.items(): + for c in cards: + m = re.search(r"_(\d+)\.txt$", c) + if m: + masses.add(m.group(1)) + + # The cross-era combination is the one place a group era and its members cannot + # both appear: the group already *is* their combination, so a card built from both + # would count those events twice. Every era still gets its own limit above. + for mass in masses: + combine_args = [] + for e in self.get_top_level_eras(): + for c in era_cards.get(e, []): + if c.endswith(f"_{mass}.txt"): + combine_args.append(f"{e}={c}") + break + + if combine_args: + cmd = ["combineCards.py"] + combine_args + out_file = os.path.join(out_dc_dir.path, f"combined_{mass}.txt") + with open(out_file, "w") as f: + subprocess.run(cmd, env=self.cmssw_env, stdout=f, check=True) diff --git a/law/StatInferenceTask.py b/law/StatInferenceTask.py new file mode 100644 index 0000000..f4d0946 --- /dev/null +++ b/law/StatInferenceTask.py @@ -0,0 +1,216 @@ +import contextlib +import law +import luigi +import os +import shutil +import yaml + +from string import Template +from FLAF.run_tools.law_customizations import Task + + +class StatInferenceTask(Task): + """Shared datacard-configuration access for the limit-setting chain. + + Everything downstream of the merged histograms -- rebinning, datacards, limits, limit + plots -- is driven by the datacard configuration named in global.yaml's + ``StatInference.config``, not by global.yaml's own variable lists. This base class is + the single place that file is read, so the four tasks below cannot disagree about + which eras, masses or variables the analysis consists of. + """ + + # The Hists_merged tree to read. Separate from `version`, which names what this chain + # *writes*: re-binning or re-fitting under a new version must not require the input + # histograms to be reproduced (or copied) under that name. Not significant -- it + # identifies an input, not a product, and law's store paths are about products. + hists_version = luigi.Parameter( + default="", + significant=False, + description="version of the Hists_merged tree to read; defaults to --version", + ) + + @property + def input_hists_version(self): + return self.hists_version or self.version + + def output_dir_target(self, *path): + """remote_dir_target() that also works when fs_default is a local directory. + + FLAF's remote_dir_target() handles a str fs and a remote fs but not a + law.LocalFileSystem: it falls through to WLCGDirectoryTarget(), which raises + "fs must be a RemoteFileSystem instance". Its sibling remote_target() (the file + variant) already carries the branch mirrored here. The CI configures fs_default + as a plain local path, so every directory output in this chain has to cope -- + without this, output() raises and luigi reports the task as having "error in + complete() method". + """ + fs = self.fs_default + if isinstance(fs, law.LocalFileSystem): + return law.LocalDirectoryTarget(os.path.join(*path), fs=fs) + return self.remote_dir_target(*path) + + def datacard_config_path(self): + return os.path.join( + self.ana_path(), self.global_params["StatInference"]["config"] + ) + + def get_config_data(self): + # Cached per instance: requires()/workflow_requires() are re-entered many times + # during graph construction and each call would otherwise re-parse the yaml. + if getattr(self, "_config_data", None) is None: + with open(self.datacard_config_path(), "r") as f: + self._config_data = yaml.safe_load(f) + return self._config_data + + def store_parts(self): + return (self.version, self.__class__.__name__, self.period) + + def datacards_dir(self, era): + """Local directory holding an era's datacards. + + The cards are produced on fs_default but must be real local files by the time + combine sees them; ResonantLimitsTask mirrors them here and everything downstream + (dhi's --multi-datacards globbing, the overlay plots) resolves against this path. + """ + return os.path.join(self.ana_data_path(), self.version, "Datacards", era) + + def preprocess_config(self): + """The datacard configuration's `preprocess:` block, or None. + + None means the chain reads the merged histograms as they are -- an analysis that + needs no transformation declares nothing and PreprocessShapesTask never runs. + """ + return self.get_config_data().get("preprocess") or None + + def get_era_groups(self): + return self.get_config_data().get("era_groups", {}) + + def get_all_eras(self): + """Every era the configuration lists, each of which gets its own datacards, its + own limit and its own plots. + + A group era and its members all appear here: they are separate measurements of + different datasets, not a double count. Only combining them into one card would + be -- see get_top_level_eras(). + """ + return self.get_config_data().get("eras", [self.period]) + + def get_top_level_eras(self): + """The eras that may be *combined* with each other: group eras plus any standalone + real era. + + Members of a group are excluded because the group already is their combination, so + a card built from both would count those events twice. This is a strictly smaller + set than get_all_eras(), and only the cross-era combination uses it. + """ + grouped = {e for sub in self.get_era_groups().values() for e in sub} + return [e for e in self.get_all_eras() if e not in grouped] + + def get_campaign(self, era): + """dhi campaign key for an era, from the config's ``campaigns`` map. + + It cannot be derived from the era name. 'Run3_Early' happens to lowercase into + dhi's 'run3_early', but 'Run3_2022' would give 'run3_2022', which is not a key of + dhi.config.campaign_labels -- and dhi draws an unknown key verbatim rather than + failing, so a guessed value ends up printed on the plot as the luminosity label. + """ + return self.get_config_data().get("campaigns", {}).get(era) + + def get_masses(self): + """Every parameter value the configuration declares, across all processes. + + Not restricted to signals: with ``param_dependent_bkg`` a background can be + parameterised too, and its histograms live in the same per-mass input file. + """ + masses = set() + for proc in self.get_config_data().get("processes", []): + if not isinstance(proc, dict): + # A bare string entry names a background directly and carries no + # settings -- the shape dc_make's load_config functions also accept. + continue + for value in proc.get("param_values") or []: + masses.add(value) + if not masses: + raise RuntimeError( + f"No process in {self.datacard_config_path()} declares param_values, so " + "there is no mass point to build inputs for." + ) + return sorted(masses) + + def get_required_variables(self): + """The Hists_merged variables this chain reads: one per mass point (the 2D + shape the datacards are built from), derived from the datacard config's + input_file_pattern. + + This is the *complete* input list -- there is no intersection with global.yaml's + histTuple_flavor variable list anywhere downstream, so a mass point present here + is read whether or not the active flavour happens to mention it. + """ + data = self.get_config_data() + model = data["model"] + pattern = model["input_file_pattern"] + param_name = model["parameters"][0] + + variables = set() + for mass in self.get_masses(): + rel = Template(pattern).safe_substitute({"ERA": "ERA", param_name: mass}) + # "//.root": the variable is the directory holding + # the file, not the file itself. + variables.add(os.path.basename(os.path.dirname(rel))) + return sorted(variables) + + def merged_hist_reqs(self, eras): + """{(era, variable): MergedHists} -- every merged input file for those eras. + + Read by CreateDatacardsTask over the sub-periods of the era it builds cards for. + """ + # Deferred: MergedHists subclasses this class, so importing it at module level + # would be circular. + from .MergedHists import MergedHists + + return { + (era, variable): MergedHists.req(self, period=era, variable=variable) + for era in eras + for variable in self.get_required_variables() + } + + @contextlib.contextmanager + def stage_inputs(self, targets): + """Assemble the required merged histograms into the "// + .root" layout that Model.getInputFileName resolves input_file_pattern + against, and yield the directory holding it. + + Each input is localized individually so only the files this task actually reads + cross the network -- the merged tree holds every variable of the active + histTuple_flavor (~106 for 'default') across every era. + """ + with contextlib.ExitStack() as stack: + staging = law.LocalDirectoryTarget(is_tmp=True) + staging.touch() + stack.callback(lambda: staging.remove(silent=True)) + + for (era, variable), target in targets.items(): + dest_dir = os.path.join(staging.abspath, era, variable) + os.makedirs(dest_dir, exist_ok=True) + local_inp = stack.enter_context(target.localize("r")) + shutil.copy2( + local_inp.abspath, os.path.join(dest_dir, f"{variable}.root") + ) + + yield staging + + def check_inputs(self, targets): + """Report every missing input at once, naming the version they were looked for + under. law already refuses to run with an incomplete MergedHists dependency; this + makes the condor log self-explanatory when the task is forced anyway, instead of + failing on whichever file ROOT happened to open first.""" + missing = sorted( + f"{era}/{var}" for (era, var), t in targets.items() if not t.exists() + ) + if missing: + raise RuntimeError( + f"{len(missing)} of {len(targets)} merged histograms are missing under " + f"hists_version='{self.input_hists_version}' " + f"(/{self.input_hists_version}/Hists_merged///.root): " + + ", ".join(missing) + ) diff --git a/law/tasks.py b/law/tasks.py index 385a965..8e07ebf 100644 --- a/law/tasks.py +++ b/law/tasks.py @@ -1,175 +1,27 @@ -import law -import luigi -import os - -from FLAF.RunKit.run_tools import ps_call -from FLAF.run_tools.law_customizations import ( - Task, - HTCondorWorkflow, - copy_param, -) -from FLAF.Analysis.tasks import HistMergerTask, HistPlotTask -from dhi.tasks.resonant import MergeResonantLimits - - -class CreateDatacardsTask(Task, HTCondorWorkflow, law.LocalWorkflow): - max_runtime = copy_param(HTCondorWorkflow.max_runtime, 2.0) - n_cpus = copy_param(HTCondorWorkflow.n_cpus, 1) - - def workflow_requires(self): - return {"HistMerger": HistMergerTask.req(self, branches=())} - - def requires(self): - merge_map = HistMergerTask.req(self, branch=-1, branches=()).create_branch_map() - - return [ - HistMergerTask.req(self, branch=br, branches=(br,)) - for br in merge_map.keys() - ] - - def create_branch_map(self): - return {0: None} - - def output(self): - path = os.path.join( - self.ana_data_path(), self.version, "Datacards", self.period - ) - return law.LocalDirectoryTarget(path) - - def run(self): - statInf_entry = self.global_params["StatInference"] - config = os.path.join(self.ana_path(), statInf_entry["config"]) - hist_bins = os.path.join(self.ana_path(), statInf_entry["hist_bins"]) - param_values = statInf_entry.get("param_values", []) - create_datacards_py = os.path.join( - self.ana_path(), "StatInference", "dc_make", "create_datacards.py" - ) - base_input_dir_remote = self.input()[0].parent.parent.parent - with base_input_dir_remote.localize("r") as base_dir_local: - cmd = [ - "python3", - create_datacards_py, - "--input", - base_dir_local.abspath, - "--output", - self.output().abspath, - "--config", - config, - "--hist-bins", - hist_bins, - "--eras", - self.period, - ] - if len(param_values) > 0: - param_values_str = ",".join(str(v) for v in param_values) - cmd += ["--param_values", param_values_str] - ps_call(cmd, env=self.cmssw_env, verbose=1) - - -class ResonantLimitsTask(Task): - workflow = luigi.Parameter(default=law.parameter.NO_STR) - - def store_parts(self): - return (self.version, self.__class__.__name__, "combined") - - def get_eras(self): - statInf_entry = self.global_params["StatInference"] - config = os.path.join(self.ana_path(), statInf_entry["config"]) - import yaml - - with open(config, "r") as f: - data = yaml.safe_load(f) - return data.get("eras", [self.period]) - - def requires(self): - return [ - CreateDatacardsTask.req(self, period=e, branches=()) - for e in self.get_eras() - ] - - def output(self): - return { - "limits": self.local_target("limits.npz"), - "datacards": law.LocalDirectoryTarget( - os.path.join( - self.ana_data_path(), self.version, "Datacards", "combined" - ) - ), - } - - def run(self): - datacards = [] - eras = self.get_eras() - era_cards = {} - import glob - import re - - for e in eras: - create_dc_br0 = CreateDatacardsTask.req( - self, period=e, branch=0, branches=() - ) - output_dir = create_dc_br0.output().abspath - cards = glob.glob(os.path.join(output_dir, "*.txt")) - era_cards[e] = cards - datacards.extend(cards) - - limits = yield MergeResonantLimits( - version=self.version, datacards=tuple(datacards) - ) - print(f"Merged limits: {limits}") - - import shutil - - self.output()["limits"].parent.touch() - shutil.copy2(limits.path, self.output()["limits"].path) - - out_dc_dir = self.output()["datacards"] - out_dc_dir.touch() - - masses = set() - for e, cards in era_cards.items(): - for c in cards: - m = re.search(r"_(\d+)\.txt$", c) - if m: - masses.add(m.group(1)) - - for mass in masses: - combine_args = [] - for e in eras: - for c in era_cards[e]: - if c.endswith(f"_{mass}.txt"): - combine_args.append(f"{e}={c}") - break - - if combine_args: - import subprocess - - cmd = ["combineCards.py"] + combine_args - out_file = os.path.join(out_dc_dir.path, f"combined_{mass}.txt") - with open(out_file, "w") as f: - subprocess.run(cmd, env=self.cmssw_env, stdout=f, check=True) - - -class ResonantLimitsAndHistPlotTask(Task): - workflow = luigi.Parameter(default=law.parameter.NO_STR) - - def get_eras(self): - statInf_entry = self.global_params["StatInference"] - config = os.path.join(self.ana_path(), statInf_entry["config"]) - import yaml - - with open(config, "r") as f: - data = yaml.safe_load(f) - return data.get("eras", [self.period]) - - def requires(self): - reqs = [ResonantLimitsTask.req(self)] - for e in self.get_eras(): - reqs.append(HistPlotTask.req(self, period=e)) - return reqs - - def output(self): - return self.local_target("dummy.txt") - - def run(self): - self.output().touch() +"""Every task in one namespace. + +Kept as a module of its own so an analysis can go on naming +``StatInference.law.tasks`` in its law.cfg ``[modules]`` list; each task lives in the +file beside this one that carries its name. law discovers tasks by walking +``Task.__subclasses__()``, so importing them here is what puts them in the index. +""" + +from .StatInferenceTask import StatInferenceTask +from .MergedHists import MergedHists +from .PreprocessShapesTask import PreprocessShapesTask +from .CreateDatacardsTask import CreateDatacardsTask +from .ResonantLimitsTask import ResonantLimitsTask +from .PlotResonantLimitsTask import PlotResonantLimitsTask +from .PlotPullsAndImpactsTask import PlotPullsAndImpactsTask +from .ResonantLimitsAndHistPlotTask import ResonantLimitsAndHistPlotTask + +__all__ = [ + "StatInferenceTask", + "MergedHists", + "PreprocessShapesTask", + "CreateDatacardsTask", + "ResonantLimitsTask", + "PlotResonantLimitsTask", + "PlotPullsAndImpactsTask", + "ResonantLimitsAndHistPlotTask", +]