diff --git a/README.md b/README.md index 03b7d4f..4844118 100644 --- a/README.md +++ b/README.md @@ -60,8 +60,13 @@ categories, control variables, variations and estimators. groups (a CROWN shift produced for some samples only). Names ending in `Up`/`Down` need their partner and become datacard shapes; other names are templates. - Estimators, run in order by `estimate`: `DataMinus(output, region, subtract, scale, - clip_negative)`, `ABCD(output, b, c, d, subtract)` and `TemplateShift(name, process, template, - fraction)`, which also varies every template variation `t` of the process (as `Up@t`). + clip_negative)`, `ABCD(output, b, c, d, subtract)`, `TemplateShift(name, process, template, + fraction)`, which also varies every template variation `t` of the process (as `Up@t`), + and `VariationSum(name, parts)`: `Up/Down` = nominal + the sum of (part - nominal) for + every histogram with a part, the part `p` being the column variation `Up%p` / + `Down%p` (`histogram.part_of`), e.g. one nuisance from CROWN shifts produced per decay + mode (exact where no event is changed by two parts). Parts never become datacard shapes; place + the sum after the estimators whose outputs should carry it. **The event rule** (`shapesmith.events`, used by hist, the ML export and measurements): cuts are the channel cuts with the region's replacements, the process cuts and the skim cuts; weights are diff --git a/src/shapesmith/datacards.py b/src/shapesmith/datacards.py index 80596cd..571fcac 100644 --- a/src/shapesmith/datacards.py +++ b/src/shapesmith/datacards.py @@ -2,14 +2,15 @@ One card per final state lists every bin (channel x category); processes with a non-positive rate in a bin are left out of that bin. Systematics: lnN from Analysis.lnn, `shape` for every Up/Down pair of variations present (names -without a direction are templates and never become shape lines), `* autoMCStats 0`. +without a direction are templates and never become shape lines; parts of summed variations are left out of the shapes +file), `* autoMCStats 0`. """ from __future__ import annotations from dataclasses import dataclass from pathlib import Path -from shapesmith.histogram import NOMINAL_VARIATION, HistKey, Histogram, HistogramSet +from shapesmith.histogram import NOMINAL_VARIATION, HistKey, Histogram, HistogramSet, is_part from shapesmith.model import NOMINAL, Analysis, Category, LnN from shapesmith.shapes import shape_name, write_shapes @@ -68,7 +69,7 @@ def key_of(process: str, variation: str = NOMINAL_VARIATION) -> HistKey: continue result.processes[process] = rebinned keys = hset.select(channel=channel_name, category=category.name, process=process, region=NOMINAL, variable=variable) - result.variations[process] = {key.variation: hset[key].rebin(edges) for key in keys if key.variation != NOMINAL_VARIATION} + result.variations[process] = {key.variation: hset[key].rebin(edges) for key in keys if key.variation != NOMINAL_VARIATION and not is_part(key.variation)} return result diff --git a/src/shapesmith/estimates.py b/src/shapesmith/estimates.py index 6dc4040..c192cb5 100644 --- a/src/shapesmith/estimates.py +++ b/src/shapesmith/estimates.py @@ -6,13 +6,18 @@ ABCD: (data - MC)(B) * (data - MC)(C).sum() / (data - MC)(D).sum(), nominal only, negative bins clipped. TemplateShift: `process` +/- `fraction` * `template`, e.g. the ttbar contamination of the embedded sample; also on top of every template variation of `process`. +VariationSum: `name`Up/Down = nominal + the sum of (part - nominal) over its parts, for every histogram with a part; +place it after the estimators whose outputs should get the summed variation (DataMinus builds its output for the +parts as for any column variation). """ from __future__ import annotations import logging -from shapesmith.histogram import NOMINAL_VARIATION, HistKey, Histogram, HistogramSet, is_template, on_template -from shapesmith.model import ABCD, NOMINAL, Analysis, Channel, ColumnVariation, DataMinus, TemplateShift +from dataclasses import replace + +from shapesmith.histogram import NOMINAL_VARIATION, HistKey, Histogram, HistogramSet, is_template, on_template, part_of +from shapesmith.model import ABCD, NOMINAL, Analysis, Channel, ColumnVariation, DataMinus, TemplateShift, VariationSum logger = logging.getLogger(__name__) @@ -109,7 +114,27 @@ def estimate_template_shift(hset: HistogramSet, channel: Channel, estimator: Tem return added -ESTIMATES = {DataMinus: estimate_data_minus, ABCD: estimate_abcd, TemplateShift: estimate_template_shift} +def estimate_variation_sum(hset: HistogramSet, channel: Channel, estimator: VariationSum) -> list[HistKey]: + """The summed variation of every nominal histogram of the channel that has a part; the variances stay the + nominal ones.""" + added = [] + for key in hset.select(channel=channel.name, variation=NOMINAL_VARIATION): + for direction in ("Up", "Down"): + name = f"{estimator.name}{direction}" + parts = [replace(key, variation=part_of(name, part)) for part in estimator.parts] + parts = [part for part in parts if part in hset] + if not parts: + continue + nominal = hset[key] + varied = nominal.copy() + for part in parts: + varied.values += hset[part].values - nominal.values + hset[replace(key, variation=name)] = varied + added.append(replace(key, variation=name)) + return added + + +ESTIMATES = {DataMinus: estimate_data_minus, ABCD: estimate_abcd, TemplateShift: estimate_template_shift, VariationSum: estimate_variation_sum} def run_estimates(hset: HistogramSet, analysis: Analysis, channels: list[str]) -> list[HistKey]: diff --git a/src/shapesmith/histogram.py b/src/shapesmith/histogram.py index a4a9a12..8daff12 100644 --- a/src/shapesmith/histogram.py +++ b/src/shapesmith/histogram.py @@ -18,11 +18,18 @@ NOMINAL_VARIATION = "Nominal" INCLUSIVE = "inclusive" # the category of control-variable histograms TEMPLATE_SEPARATOR = "@" +PART_SEPARATOR = "%" def is_template(variation: str) -> bool: - """A variation without a direction (e.g. one point of an energy-scale grid), not one on top of a template.""" - return variation != NOMINAL_VARIATION and not variation.endswith(("Up", "Down")) and TEMPLATE_SEPARATOR not in variation + """A variation without a direction (e.g. one point of an energy-scale grid), not one on top of a template and + not a part of a summed variation.""" + return ( + variation != NOMINAL_VARIATION + and not variation.endswith(("Up", "Down")) + and TEMPLATE_SEPARATOR not in variation + and PART_SEPARATOR not in variation + ) def on_template(variation: str, template: str) -> str: @@ -31,6 +38,17 @@ def on_template(variation: str, template: str) -> str: return f"{variation}{TEMPLATE_SEPARATOR}{template}" +def part_of(variation: str, part: str) -> str: + """The name of the part `part` of the summed variation `variation` (a VariationSum). It ends in the part's name, + so datacards and template shifts never read it.""" + return f"{variation}{PART_SEPARATOR}{part}" + + +def is_part(variation: str) -> bool: + """A part of a summed variation: an input of the VariationSum, never written to a shapes file.""" + return PART_SEPARATOR in variation + + @dataclass class Histogram: edges: np.ndarray diff --git a/src/shapesmith/model.py b/src/shapesmith/model.py index d991058..2992370 100644 --- a/src/shapesmith/model.py +++ b/src/shapesmith/model.py @@ -175,7 +175,18 @@ class TemplateShift: fraction: float -Estimator = DataMinus | ABCD | TemplateShift +@dataclass(frozen=True) +class VariationSum: + """Variation `name`Up/Down of every process with a histogram of one of its parts: the nominal plus the sum over + `parts` of the part's difference to the nominal, e.g. one nuisance from CROWN shifts produced per decay mode. + The part `p` of a direction is the column variation `histogram.part_of(f"{name}{direction}", p)`. Exact where no + event is changed by two parts; a process without a part histogram gets no variation.""" + + name: str + parts: tuple[str, ...] + + +Estimator = DataMinus | ABCD | TemplateShift | VariationSum @dataclass(frozen=True) diff --git a/src/shapesmith/shapes.py b/src/shapesmith/shapes.py index 3e0a574..2f2b9c7 100644 --- a/src/shapesmith/shapes.py +++ b/src/shapesmith/shapes.py @@ -6,7 +6,7 @@ import uproot -from shapesmith.histogram import NOMINAL_VARIATION, HistKey, Histogram, HistogramSet +from shapesmith.histogram import NOMINAL_VARIATION, HistKey, Histogram, HistogramSet, is_part from shapesmith.model import NOMINAL, Analysis @@ -29,8 +29,8 @@ def _synced_name(category: str, process: str, variation: str) -> str: def synced_entries(hset: HistogramSet, analysis: Analysis, channel_name: str, name: Callable[[str, str, str], str] = _synced_name) -> list[tuple[str, str, Histogram]]: - """data_obs and every nominal-region histogram of the signal and the backgrounds, per category, each named - `name(category, process, variation)`.""" + """data_obs and every nominal-region histogram of the signal and the backgrounds but the parts of summed + variations, per category, each named `name(category, process, variation)`.""" channel = analysis.channel(channel_name) entries = [] for category in channel.categories: @@ -41,6 +41,8 @@ def synced_entries(hset: HistogramSet, analysis: Analysis, channel_name: str, na entries.append((folder, "data_obs", hset[data_key])) for process in (*channel.backgrounds(), analysis.signal): for key in hset.select(channel=channel_name, category=category.name, process=process, region=NOMINAL, variable=variable): + if is_part(key.variation): + continue entries.append((folder, name(category.name, process, key.variation), hset[key])) return entries diff --git a/src/shapesmith/validate.py b/src/shapesmith/validate.py index f24737d..9b0ff44 100644 --- a/src/shapesmith/validate.py +++ b/src/shapesmith/validate.py @@ -1,7 +1,8 @@ """Consistency checks of an Analysis; every problem is collected and reported at once.""" from __future__ import annotations -from shapesmith.model import NOMINAL, ROLES, SAMPLE_KINDS, Analysis, AnalysisError, Channel, ColumnVariation, DataMinus, TemplateShift, WeightVariation +from shapesmith.histogram import part_of +from shapesmith.model import NOMINAL, ROLES, SAMPLE_KINDS, Analysis, AnalysisError, Channel, ColumnVariation, DataMinus, TemplateShift, VariationSum, WeightVariation def _duplicates(what: str, names: list[str]) -> list[str]: @@ -68,10 +69,19 @@ def _variation_problems(channel: Channel, region_names: set[str]) -> list[str]: def _estimator_problems(channel: Channel, region_names: set[str]) -> list[str]: roles = {p.name: p.role for p in channel.processes} problems = [] + column_variations = {v.name for v in channel.variations if isinstance(v, ColumnVariation)} for estimator in channel.estimators: if isinstance(estimator, TemplateShift): problems += [f"template shift {estimator.name}: unknown process {name}" for name in (estimator.process, estimator.template) if name not in roles] continue + if isinstance(estimator, VariationSum): + if _partner(estimator.name) is not None: + problems.append(f"variation sum {estimator.name}: the name carries a direction") + if not estimator.parts: + problems.append(f"variation sum {estimator.name}: no parts") + parts = [part_of(f"{estimator.name}{d}", p) for p in estimator.parts for d in ("Up", "Down")] + problems += [f"variation sum {estimator.name}: part {part} is not a column variation" for part in parts if part not in column_variations] + continue regions = (estimator.region,) if isinstance(estimator, DataMinus) else (estimator.b, estimator.c, estimator.d) problems += [f"estimator {estimator.output}: region {r} does not exist" for r in regions if r not in region_names] problems += [f"estimator {estimator.output}: subtracts unknown process {name}" for name in estimator.subtract if name not in roles] diff --git a/tests/test_datacards.py b/tests/test_datacards.py index 0dcce7b..e99fdfb 100644 --- a/tests/test_datacards.py +++ b/tests/test_datacards.py @@ -1,7 +1,7 @@ import uproot from shapesmith.datacards import bin_name, rebin_edges, run_datacards, write_datacard -from shapesmith.histogram import HistKey, Histogram, HistogramSet +from shapesmith.histogram import HistKey, Histogram, HistogramSet, part_of from tests.mini_analysis import build, build_embedding @@ -78,3 +78,19 @@ def test_auxiliary_processes_and_templates_stay_out_of_the_card(tmp_path): processes = [line for line in text.splitlines() if line.startswith("process")][0].split()[1:] assert "ZTT" not in processes and "EMB" in processes # ZTT is the auxiliary template of the embedding assert "emb1p002" not in text + + +def test_parts_of_a_summed_variation_stay_out_of_the_card_and_the_shapes(tmp_path): + hset = _hset() + for category in ("sig", "bkg"): + hset[HistKey("mt", category, "EMB", "nominal", "Nominal", "score")] = _hist([3, 2, 1, 0.5]) + for d in ("Up", "Down"): + hset[HistKey("mt", category, "EMB", "nominal", part_of(f"CMS_tes{d}", "dm10"), "score")] = _hist([3, 2, 1, 0.6]) + hset[HistKey("mt", category, "EMB", "nominal", f"CMS_tes{d}", "score")] = _hist([3, 2, 1, 0.6]) + card = write_datacard(hset, build_embedding(), ["mt"], tmp_path / "mt") + text = card.read_text() + assert "CMS_tes shape" in " ".join(text.split()) and "%" not in text + with uproot.open(card.parent / "common" / "htt_input_2018.root") as f: + names = f.keys(recursive=True, cycle=False) + assert "htt_mt_1_2018/EMB_CMS_tesUp" in names and not [name for name in names if "%" in name] + diff --git a/tests/test_estimates.py b/tests/test_estimates.py index b6e6f36..5cda152 100644 --- a/tests/test_estimates.py +++ b/tests/test_estimates.py @@ -3,8 +3,8 @@ import pytest from shapesmith.estimates import clip_negative_bins, run_estimates -from shapesmith.histogram import INCLUSIVE, HistKey, Histogram, HistogramSet -from shapesmith.model import ABCD +from shapesmith.histogram import INCLUSIVE, HistKey, Histogram, HistogramSet, part_of +from shapesmith.model import ABCD, ColumnVariation, DataMinus, VariationSum from tests.mini_analysis import build, build_embedding @@ -107,3 +107,40 @@ def test_template_shift_also_shifts_every_template_variation(): assert hset[_key("EMB", "nominal", variation="CMS_htt_emb_ttbar_2018Up@es+2")].values.tolist() == pytest.approx([11.2, 9.4]) assert hset[_key("EMB", "nominal", variation="CMS_htt_emb_ttbar_2018Down@es+2")].values.tolist() == pytest.approx([10.8, 8.6]) assert not hset.select(process="EMB", variation="CMS_htt_emb_ttbar_2018Up@CMS_tesUp") + + +def _with_tes_parts(analysis): + """The tau ES of the embedding as one nuisance from two parts, e.g. CROWN shifts per decay mode.""" + channel = analysis.channel("mt") + parts = tuple( + ColumnVariation(part_of(f"CMS_tes{d}", dm), f"__tes_{dm}{d}", applies_to=("mc", "embedding")) for dm in ("dm10", "dm11") for d in ("Up", "Down") + ) + variations = tuple(v for v in channel.variations if not v.name.startswith("CMS_tes")) + parts + estimators = (DataMinus("jetFakes", "anti_iso", ("EMB", "ZL")), VariationSum("CMS_tes", ("dm10", "dm11"))) + return dataclasses.replace(analysis, channels={"mt": dataclasses.replace(channel, variations=variations, estimators=estimators)}) + + +def test_variation_sum_adds_the_differences_of_its_parts(): + hset = HistogramSet() + hset[_key("EMB", "nominal")] = _hist([10.0, 10.0], [1.0, 1.0]) + hset[_key("EMB", "nominal", variation=part_of("CMS_tesUp", "dm10"))] = _hist([11.0, 10.0], [2.0, 2.0]) + hset[_key("EMB", "nominal", variation=part_of("CMS_tesUp", "dm11"))] = _hist([10.0, 13.0], [2.0, 2.0]) + hset[_key("EMB", "nominal", variation=part_of("CMS_tesDown", "dm10"))] = _hist([9.0, 10.0]) + hset[_key("ZL", "nominal")] = _hist([5.0, 5.0]) # no part: no summed variation + run_estimates(hset, _with_tes_parts(build_embedding()), ["mt"]) + up = hset[_key("EMB", "nominal", variation="CMS_tesUp")] + assert up.values.tolist() == [11.0, 13.0] and up.variances.tolist() == [1.0, 1.0] # the nominal variances + assert hset[_key("EMB", "nominal", variation="CMS_tesDown")].values.tolist() == [9.0, 10.0] # the one part present + assert not hset.select(process="ZL", variation="CMS_tesUp") + + +def test_variation_sum_reaches_the_data_minus_output_through_its_parts(): + hset = HistogramSet() + hset[_key("data", "anti_iso")] = _hist([10.0, 20.0]) + hset[_key("EMB", "anti_iso")] = _hist([1.0, 2.0]) + hset[_key("EMB", "anti_iso", variation=part_of("CMS_tesUp", "dm10"))] = _hist([2.0, 2.0]) + hset[_key("EMB", "anti_iso", variation=part_of("CMS_tesUp", "dm11"))] = _hist([1.0, 4.0]) + run_estimates(hset, _with_tes_parts(build_embedding()), ["mt"]) + # jetFakes = data - EMB, nominal (9, 18), parts (8, 18) and (9, 16) + assert hset[_key("jetFakes", "nominal", variation="CMS_tesUp")].values.tolist() == [8.0, 16.0] + diff --git a/tests/test_histogram.py b/tests/test_histogram.py index 0b09a43..b270b65 100644 --- a/tests/test_histogram.py +++ b/tests/test_histogram.py @@ -2,7 +2,7 @@ import pytest import uproot -from shapesmith.histogram import HistKey, Histogram, HistogramSet, is_template, on_template, unchanged_variations +from shapesmith.histogram import HistKey, Histogram, HistogramSet, is_template, on_template, part_of, unchanged_variations def test_fill_add_scale_sum(): @@ -54,3 +54,4 @@ def test_template_names(): assert is_template("es+2") and is_template("es-198") assert not is_template("Nominal") and not is_template("CMS_tesUp") and not is_template("CMS_tesDown") assert on_template("CMS_ttbarUp", "es-2") == "CMS_ttbarUp@es-2" and not is_template("CMS_ttbarUp@es-2") + assert part_of("CMS_tesUp", "dm10") == "CMS_tesUp%dm10" and not is_template("CMS_tesUp%dm10") diff --git a/tests/test_shapes.py b/tests/test_shapes.py index 9bc3f60..bed2b4b 100644 --- a/tests/test_shapes.py +++ b/tests/test_shapes.py @@ -1,6 +1,6 @@ import uproot -from shapesmith.histogram import HistKey, Histogram, HistogramSet +from shapesmith.histogram import HistKey, Histogram, HistogramSet, part_of from shapesmith.shapes import run_sync, write_shapes from tests.mini_analysis import build @@ -23,6 +23,7 @@ def test_run_sync_layout(tmp_path): hset[HistKey("mt", category, "data", "anti_iso", "Nominal", "score")] = _hist(1.0) hset[HistKey("mt", category, "ZTT", "nominal", "Nominal", "score")] = _hist(2.0) hset[HistKey("mt", category, "ZTT", "nominal", "CMS_puUp", "score")] = _hist(2.1) + hset[HistKey("mt", category, "ZTT", "nominal", part_of("CMS_tesUp", "dm10"), "score")] = _hist(2.2) # a part, never synced hset[HistKey("mt", category, "ZTT", "anti_iso", "Nominal", "score")] = _hist(0.2) hset[HistKey("mt", category, "HH", "nominal", "Nominal", "score")] = _hist(0.1) hset[HistKey("mt", category, "jetFakes", "nominal", "Nominal", "score")] = _hist(0.8) diff --git a/tests/test_validate.py b/tests/test_validate.py index 9977af9..cee72fe 100644 --- a/tests/test_validate.py +++ b/tests/test_validate.py @@ -3,6 +3,7 @@ import pytest from shapesmith import model +from shapesmith.histogram import part_of from shapesmith.validate import validate from tests.mini_analysis import build, build_embedding @@ -105,3 +106,19 @@ def test_the_signal_is_a_process_with_the_signal_role(): def test_dataclasses_are_frozen(): with pytest.raises(dataclasses.FrozenInstanceError): build().signal = "other" + + +def test_validate_checks_variation_sums(): + analysis = build_embedding() + channel = analysis.channel("mt") + parts = tuple(model.ColumnVariation(part_of(f"CMS_tes{d}", "dm10"), f"__tes_dm10{d}") for d in ("Up", "Down")) + ok = model.VariationSum("CMS_tes", ("dm10",)) + validate(_with_channel(analysis, variations=channel.variations + parts, estimators=(ok,))) + for estimator, problem in ( + (model.VariationSum("CMS_tesUp", ("dm10",)), "the name carries a direction"), + (model.VariationSum("CMS_tes", ()), "no parts"), + (model.VariationSum("CMS_tes", ("dm11",)), r"part CMS_tesUp%dm11 is not a column variation"), + ): + with pytest.raises(model.AnalysisError, match=problem): + validate(_with_channel(analysis, variations=channel.variations + parts, estimators=(estimator,))) +