Feat: variation sum

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 `<name>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 `<name>Up@t`),
+  and `VariationSum(name, parts)`: `<name>Up/Down` = nominal + the sum of (part - nominal) for
+  every histogram with a part, the part `p` being the column variation `<name>Up%p` /
+  `<name>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

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
 
 

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]:

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

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)

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
 

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]

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]
+

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]
+

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")

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)

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,)))
+