diff --git a/README.md b/README.md index 690d918..03b7d4f 100644 --- a/README.md +++ b/README.md @@ -1,18 +1,18 @@ # ShapeSmith -From CROWN ntuples to histograms, background estimates, combine datacards, fits, plots and -ML training folds. Columnar (uproot, pandas, numexpr, numpy), no ROOT needed except for -`shapesmith fit`, which calls combine in a CMSSW subshell. +From CROWN ntuples to histograms, background estimates, combine datacards, fits, plots, +ML training folds and measurements (correctionlib payloads). Columnar (uproot, pandas, +numexpr, numpy), no ROOT needed except for `shapesmith fit`, which calls combine in a CMSSW +subshell. -An analysis is a Python package that builds a `shapesmith.model.Analysis` (samples, -selections, processes, regions, categories, variables, systematics, style) from a small run +An analysis is a Python package that builds a `shapesmith.model.Analysis` from a small run YAML; ShapeSmith runs the steps. Example analysis: [BBTauTauAnalysis-ShapeSmith](https://github.com/KIT-CMS/BBTauTauAnalysis-ShapeSmith). ## Install ```bash -source /cvmfs/sft.cern.ch/lcg/views/LCG_108/x86_64-el9-gcc15-opt/setup.sh # Python 3.12, uproot, pandas, mplhep, XRootD +source /cvmfs/sft.cern.ch/lcg/views/LCG_108/x86_64-el9-gcc15-opt/setup.sh # Python 3.12, uproot, pandas, mplhep, correctionlib, XRootD python3 -m venv --system-site-packages ~/.venvs/shapesmith && source ~/.venvs/shapesmith/bin/activate pip install -e ".[test]" # editable: changes in this checkout are live pytest # tests on synthetic ntuples, no network @@ -21,48 +21,125 @@ pytest # tests on synthetic ntuples, no network ## Steps ```bash -shapesmith example-config run.yaml # annotated run configuration -shapesmith validate -c run.yaml # build + validate the analysis, list the columns each channel needs -shapesmith skim -c run.yaml [--channels mt] [--samples TT,SingleMuon] [--workers 8] # ntuples (+ friends) -> Parquet -shapesmith hist -c run.yaml [--control] [--regions all] [--skip-systematics] [--processes ztt,data] # histograms (ROOT + JSON index) -shapesmith estimate -c run.yaml [--control] # data-driven processes (fake factors / ABCD), embedding variations +shapesmith validate -c run.yaml # build + validate the analysis, list processes and needed columns per channel +shapesmith skim -c run.yaml [--channels mt] [--samples TT,SingleMuon] [--force] [--workers 8] # ntuples (+ friends) -> Parquet +shapesmith hist -c run.yaml [--control] [--regions all] [--skip-systematics] [--processes ZTT,data] # histograms (ROOT + JSON index) +shapesmith estimate -c run.yaml [--control] # estimated processes and variations (Channel.estimators) +shapesmith measure -c run.yaml [--suggest-binning] [--merge] # Analysis.measurement -> /// shapesmith plot -c run.yaml [--control] [--region nominal] [--blind] [--log] shapesmith sync -c run.yaml # combine-style shape files per channel shapesmith datacards -c run.yaml [--min-background 1.0] [--no-systematics] -shapesmith fit -c run.yaml [--final-states mt,all] [--skip-combine] +shapesmith fit -c run.yaml [--final-states mt,all] [--skip-combine] # limits, significance, best fit shapesmith ml-export -c run.yaml # Feather training folds -shapesmith inspect output/shapes.root # what a histogram file contains -shapesmith inventory datasets.json nicks.txt inventory.txt # production inventory (nick + DBS path) +shapesmith inspect output/shapes.root [--unchanged] # what a histogram file contains; variations equal to their nominal ``` -Two stages: `skim` reads the ntuples once (loose selection, all needed columns, normalisation -weight `norm_weight = xsec / (nevents * generator_weight) * sign(genWeight)`) into -`///*.parquet` plus a `manifest.json`; everything else works on the -skims. `skim`, `hist` and `ml-export` write a `versions.json` (versions, git hashes of core, -analysis and sample database, full run configuration) next to their outputs. - -For named controls, `hist --regions same_sign,anti_iso` fills those regions, including -signal contamination; `--regions all` expands each channel's defined regions plus nominal. -Without this option the existing nominal/estimator booking applies. Use `--force` to refill -an existing histogram file; partial refills preserve other regions, variables, channels and -processes. `plot --region same_sign` reads that region and writes beneath -`plots//same_sign/`, with the region in the plot label. A nominal background -estimate is never substituted into another region. Unblinded plots require actual data -histograms; only explicit `--blind` uses Asimov data. - -A sample's own `cut` (e.g. one generator-level part of an inclusive sample) is applied at skim -time on top of the loose selection; its `norm_weight` stays that of the whole sample. Only that -sample reads the cut's columns and only its manifest records the cut, so the skims of the other -samples stay reusable. - -Skim manifests record cuts, required columns and normalization. Reuse after changing this -contract fails with `skim --force` guidance; legacy manifests without recorded cuts also -require regeneration. The manifest is written before the first file of a sample and after -every finished file (`completed`, with per-file event counts), so an interrupted skim resumes -with only the unfinished files; `--force` starts the record afresh. Parquet files and manifests -are replaced atomically through temporary files. A region weight that replaces a named -baseline weight keeps its MC-only applicability. Newly added weights, such as a fake factor, -also apply to data. +`skim` reads the ntuples once into `///*.parquet` plus a +`manifest.json`; everything else works on the skims. `skim`, `hist`, `measure` and `ml-export` +record what produced their outputs in `/versions/.json` (versions, +git hashes of core, analysis and sample database, the full run configuration and its hash). + +## The analysis model + +Everything that differs by channel lives in its `Channel`: samples, cuts, processes, regions, +categories, control variables, variations and estimators. + +- `Sample(nick, group, kind, xsec, nevents, generator_weight, cut)`: kind `data`, `mc` or + `embedding`; `norm_weight = xsec / (nevents * generator_weight) * sign(genWeight)` for MC. + A sample's own `cut` (e.g. one generator-level part of an inclusive sample) is applied at skim + time; its normalisation stays that of the whole sample. +- `Process(name, group, role, plot_group, selection)`: role `data`, `signal`, `background` or + `auxiliary` (booked in the nominal region only, for an estimator; never in datacards, plots or + the ML export). A process owns its whole weight set. +- `Region(name, replace_cuts, add_weights, replace_weights)`: cuts replaced by name, weights + replaced where a process carries them, weights added to every process. +- `WeightVariation(name, replace_weights, applies_to, regions)` and + `ColumnVariation(name, suffix | derived, applies_to, regions, groups)`: a column variation reads + column `c` as `c + suffix` where that branch exists (a CROWN shift), or replaces it by + `derived[c]`, an expression of nominal columns; `groups` restricts it to the samples of some + 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`). + +**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 +the process weights with the region's replacements, then the region's added weights; a weight +variation replaces weights in place (and does not apply to a process lacking one of them), a +column variation rewrites every expression. Event weight: `(norm_weight * lumi) * product`, +`lumi` for MC only. + +**Booking**: the signal in the nominal region, other processes in the nominal and the estimator +regions (or `--regions`). A variation is filled where its `applies_to` and `regions` allow; a +column variation only where it changes an expression of the histogram (elsewhere it equals the +nominal). `DataMinus` builds its output for every column variation of data or a subtracted +process, taking the nominal of inputs that lack it; weight variations are not propagated, ABCD +is nominal only. `--skip-systematics` skips weight and column variations. `hist` always refills +the requested scopes and keeps the other histograms of its file. + +## Skims + +Per sample, the skim keeps the columns of every expression its processes can use and the +shifted branches of its column variations. An event is kept if the skim cuts pass nominally or +under any column variation of its kind (hist re-applies the varied skim cuts). Where the channel +declares CROWN shifts for a sample (its kind and group), every file must carry a shifted branch for each of them, +and every shifted branch `c__X` of a needed column must be declared; otherwise the skim fails. + +The manifest records the skim contract. A stored skim is reused when its skim cuts, +normalisation and sample cut are unchanged, its recorded column variations and friends contain +the current ones, and its Parquet schema holds the needed columns; so one skim made with all +variations and friends serves every subset. Incompatible skims are listed together, with the +`skim --force --samples ...` fix. A new friend version gets a new friend base (a friend rewritten +in place is not detected). The manifest is written before the first file of a sample and after +every finished file, so an interrupted skim resumes with only the unfinished files. Parquet +files and manifests are replaced atomically. + +## Samples + +`shapesmith.samples.read_sample_list` reads a KingMaker sample list copied unchanged (one nick per +line); `normalisation(database, nicks)` looks every nick up in the KingMaker sample database +(`datasets.json`) and reports all missing nicks at once. The database renames nicks now and +then, so an older production needs the database checkout it was produced with. + +## Measurements + +`Analysis.measurement` is an object with a `name` and `run(context)`; `shapesmith measure` calls +it with a `MeasureContext` (configuration, analysis, channels, output directory, +`events(query)`, `provenance()`, and the flags `suggest_binning` and `merge`, the latter to combine +the results of earlier runs into one payload). Measurements: + +- `shapesmith.measurements.tau_id_es`: tau-ID scale factor and energy scale of embedded taus, the + chain of smhtt_ul `tauID_SFs_dev`. One run per working-point combination fills the mt channel + and the mm control region, writes the shapes in the input format of MorphingTauID2017 and runs + per category in CMSSW the datacards, T2W, a 2D likelihood scan, close-up 1D scans of both POIs + and the MultiDimFit singles fit, the result (`results.json`, plots). As in the predecessor, the + close-up scans and the singles fit start at the lowest grid point of the step before, within its + 2 sigma intervals widened by a margin. The energy scale is a grid of template column variations + in units of 0.1 % (`grid.py`); `--merge` writes the correctionlib payload of all combinations. + A 1 sigma interval at the fit range or a failed crossing is an error. +- `shapesmith.measurements.fake_factors`: the fake-factor measurement (the SM method of + TauFakeFactors). An analysis sets `Analysis.measurement = FakeFactorMeasurement(legs)`, built + from its region and process names and per-channel tables (`Binned`, `Fit`, `Split`, + `ProcessFF`, `DrSr`, `DataScale`, `Fractions`, `Leg`). `shapesmith measure` writes + `fake_factors_.json.gz` and `FF_corrections_.json.gz` (the conventions of the CROWN + fake-factor friend), the intermediate DR->SR payload, `measurement_.json` (every ratio, fit + input and result) and plots; `--suggest-binning` prints equipopulated edges instead. The + histograms (`hist.py`: centre of mass, MC-suppressed errors), fits (`fit.py`: bands, + SystMCShift, the compatibility-with-1 reset, sparsify) and the kernel are checked against ROOT + and TauFakeFactors (`tests/reference/make_fake_factor_reference.py`); the normative algorithm is + the "Algorithm reference" of the design spec. + +Building blocks: + +- `shapesmith.payloads`: correctionlib schema-v2 builders; the provenance is JSON in + `CorrectionSet.description`; payloads are checked (schema, evaluator, no `$schema` key, so + correctionlib 2.6 reads them) and written gzip compressed. +- `shapesmith.measurements.smoothing`: TauFakeFactors' kernel smoothing without ROOT, a port of + ROOT's `TGraphSmooth::SmoothKern` checked against ROOT and FF_Updated + (`tests/reference/make_smoothing_reference.py`). +- `shapesmith.cmssw`: run commands in a CMSSW environment (`combine.cmssw_dir`), started from a + clean environment so that the Python environment of ShapeSmith does not shadow CMSSW's. ## Run configuration @@ -94,5 +171,5 @@ combine: {cmssw_dir: /work/USER/CMSSW_14_1_9, scram_arch: el9_amd64_gcc12} - `workers > 1` uses a spawn-based process pool (fork deadlocks with uproot/pyarrow threads). Scripts that call ShapeSmith functions with more than one worker therefore need the usual `if __name__ == "__main__":` guard; `--workers 1` runs everything inline for debugging. -- Systematics are weight variations from ntuple columns plus lnN normalisation uncertainties; - NN scores and fake factors come from CROWN friend trees. +- `tests/golden/` holds the mini-dataset outputs of core `aa5daa5`; the rebuilt core reproduces + them bitwise (`tests/test_golden.py`). diff --git a/environment.yml b/environment.yml index a29d712..dcbd6bb 100644 --- a/environment.yml +++ b/environment.yml @@ -6,7 +6,8 @@ dependencies: - numexpr>=2.8 - pandas>=2.1 - pyarrow>=15 - - numpy>=1.26 + - numpy>=2.0 + - correctionlib>=2.6 - mplhep>=0.3.30 - matplotlib>=3.8 - typer>=0.12 diff --git a/pyproject.toml b/pyproject.toml index b5c1535..897e48a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "shapesmith" -description = "From CROWN ntuples to histograms, background estimates, datacards, fits and ML training folds - without ROOT." +description = "From CROWN ntuples to histograms, background estimates, datacards, fits, ML training folds and measurements - without ROOT." readme = "README.md" requires-python = ">=3.11" license = {text = "MIT"} @@ -14,7 +14,8 @@ dependencies = [ "numexpr>=2.8", "pandas>=2.1", "pyarrow>=15", - "numpy>=1.26", + "numpy>=2.0", + "correctionlib>=2.6", "mplhep>=0.3.30", "matplotlib>=3.8", "typer>=0.12", diff --git a/src/shapesmith/cli.py b/src/shapesmith/cli.py index 455bb38..9301b89 100644 --- a/src/shapesmith/cli.py +++ b/src/shapesmith/cli.py @@ -1,4 +1,4 @@ -"""Command-line interface: one sub-command per analysis step (Spec §13).""" +"""Command-line interface: one sub-command per analysis step.""" from __future__ import annotations import logging @@ -11,13 +11,13 @@ from rich.table import Table from shapesmith import __version__ -from shapesmith.config import RunConfig, load_analysis, load_config, write_example_config, write_versions +from shapesmith.config import RunConfig, load_analysis, load_config +from shapesmith.provenance import write_versions -app = typer.Typer(help=f"ShapeSmith {__version__}: CROWN ntuples -> skims -> histograms -> datacards, fits, plots, ML folds.", no_args_is_help=True) +app = typer.Typer(help=f"ShapeSmith {__version__}: CROWN ntuples -> skims -> histograms -> datacards, fits, plots, ML folds, measurements.", no_args_is_help=True) ConfigOption = typer.Option(..., "--config", "-c", help="run configuration YAML") ChannelsOption = typer.Option(None, "--channels", help="comma separated subset of the configured channels") -ForceOption = typer.Option(False, "--force", help="recreate existing outputs") WorkersOption = typer.Option(None, "--workers", help="override `workers` of the configuration (1 = inline, best for debugging)") SetOption = typer.Option(None, "--set", "-s", help="override a configuration entry, e.g. -s switches.jet_fakes=ff -s workers=4 (repeatable, values parsed as YAML)") @@ -46,100 +46,87 @@ def _final_states(channels: list[str]) -> dict[str, list[str]]: return states -@app.command("example-config") -def example_config(path: Path = typer.Argument(..., help="where to write the example YAML")): - """Write an example run configuration.""" - write_example_config(path) - typer.echo(f"wrote {path}") - - -@app.command() -def inventory(database: Path = typer.Argument(..., help="KingMaker datasets.json that still knows the nicks by name"), nicks: Path = typer.Argument(..., help="text file, one produced nick per line"), output: Path = typer.Argument(..., help="inventory to write: nick and DBS path per line")): - """Write a production inventory: for every produced nick its DBS dataset path, the key that survives database renames.""" - from shapesmith.samples import write_inventory - - names = [line.strip() for line in nicks.read_text().splitlines() if line.strip()] - entries = write_inventory(output, database, names) - typer.echo(f"wrote {len(entries)} entries to {output}") +def _split(value: Optional[str]) -> Optional[list[str]]: + return value.split(",") if value else None @app.command() def validate(config: Path = ConfigOption, channels: Optional[str] = ChannelsOption, overrides: Optional[list[str]] = SetOption): - """Load and validate the analysis; report the columns every channel needs.""" - from shapesmith.skim import required_columns + """Load and validate the analysis; report processes, selections and the columns every channel needs.""" + from shapesmith.skim import needed_columns cfg, analysis, selected = _setup(config, channels, overrides=overrides) console = Console() console.print(f"[bold]{analysis.name}[/bold] era {analysis.era}, {analysis.lumi_pb:g} pb^-1, switches {dict(cfg.switches)}") - console.print(f"channels: {', '.join(selected)} samples: {len(analysis.samples)} ({', '.join(f'{kind} {sum(s.kind == kind for s in analysis.samples)}' for kind in sorted({s.kind for s in analysis.samples}))})") - table = Table(title="processes", show_edge=False) - for column in ("key", "name", "kind", "group", "plot group", *[f"samples {c}" for c in selected]): - table.add_column(column) - for process in analysis.processes: - table.add_row(process.key, process.name, process.kind, process.group, process.plot_group, *[str(len(analysis.samples_for(process.group, c))) for c in selected]) - console.print(table) - for channel in selected: - ch = analysis.channel(channel) - console.print(f"{channel}: baseline cuts {list(ch.baseline.cuts)}, weights {list(ch.baseline.weights)}, regions {[r.name for r in ch.regions]}") - console.print(f"categories: {[c.name for c in analysis.categories] or 'none'} control variables: {len(analysis.control_variables)}") - console.print(f"systematics: {len(analysis.weight_variations)} weight variations, {len(analysis.lnn)} lnN estimator: {analysis.estimator.name + ' -> ' + analysis.estimator.output if analysis.estimator else 'none'}") - for channel in selected: - typer.echo(f"{channel}: {len(required_columns(analysis, channel))} columns required") + for name in selected: + channel = analysis.channel(name) + kinds = sorted({s.kind for s in channel.samples}) + console.print(f"[bold]{name}[/bold]: samples {len(channel.samples)} ({', '.join(f'{k} {sum(s.kind == k for s in channel.samples)}' for k in kinds)})") + table = Table(show_edge=False) + for column in ("process", "role", "group", "plot group", "samples"): + table.add_column(column) + for process in channel.processes: + table.add_row(process.name, process.role, process.group, process.plot_group, str(len(channel.samples_of(process.group)))) + console.print(table) + console.print(f" cuts {list(channel.cuts)}, regions {[r.name for r in channel.regions]}") + console.print(f" categories {[c.name for c in channel.categories] or 'none'}, control variables {len(channel.variables)}, variations {len(channel.variations)}") + console.print(f" estimators {[type(e).__name__ for e in channel.estimators] or 'none'}, backgrounds {list(channel.backgrounds())}") + columns = set().union(*(needed_columns(channel, s) for s in channel.samples)) + typer.echo(f"{name}: {len(columns)} columns required") typer.echo(f"analysis {analysis.name} is valid") @app.command() -def skim(config: Path = ConfigOption, channels: Optional[str] = ChannelsOption, overrides: Optional[list[str]] = SetOption, force: bool = ForceOption, samples: Optional[str] = typer.Option(None, "--samples", help="comma separated sample groups or nick prefixes (default: all)"), workers: Optional[int] = WorkersOption): +def skim(config: Path = ConfigOption, channels: Optional[str] = ChannelsOption, overrides: Optional[list[str]] = SetOption, force: bool = typer.Option(False, "--force", help="skim again, ignoring stored skims"), samples: Optional[str] = typer.Option(None, "--samples", help="comma separated sample groups or nick prefixes (default: all)"), workers: Optional[int] = WorkersOption): """Stage 1: ntuples (+ friends) -> Parquet skims.""" from shapesmith.skim import run_skim cfg, analysis, selected = _setup(config, channels, workers=workers, overrides=overrides) - write_versions(cfg, cfg.skim_dir) - results = run_skim(cfg, analysis, selected, force=force, samples=samples.split(",") if samples else None) + write_versions(cfg, analysis.name, cfg.skim_dir) + results = run_skim(cfg, analysis, selected, force=force, samples=_split(samples)) typer.echo(f"{sum(1 for r in results if not r.skipped)} files skimmed, {sum(1 for r in results if r.skipped)} reused, {sum(r.n_out for r in results)} events kept") @app.command() -def hist(config: Path = ConfigOption, channels: Optional[str] = ChannelsOption, overrides: Optional[list[str]] = SetOption, control: bool = typer.Option(False, "--control", help="control variables instead of NN categories"), variables: Optional[str] = typer.Option(None, help="comma separated control variables"), regions: Optional[str] = typer.Option(None, "--regions", help="comma separated regions, or 'all' (default: nominal plus estimator regions)"), skip_systematics: bool = typer.Option(False, "--skip-systematics"), processes: Optional[str] = typer.Option(None, help="comma separated process keys"), force: bool = ForceOption, workers: Optional[int] = WorkersOption): - """Stage 2: skims -> histograms (control_shapes.root or shapes.root).""" - from shapesmith.histograms import run_hist +def hist(config: Path = ConfigOption, channels: Optional[str] = ChannelsOption, overrides: Optional[list[str]] = SetOption, control: bool = typer.Option(False, "--control", help="control variables instead of NN categories"), variables: Optional[str] = typer.Option(None, help="comma separated control variables"), regions: Optional[str] = typer.Option(None, "--regions", help="comma separated regions, or 'all' (default: nominal plus estimator regions)"), skip_systematics: bool = typer.Option(False, "--skip-systematics", help="no weight or column variations"), processes: Optional[str] = typer.Option(None, help="comma separated process names"), workers: Optional[int] = WorkersOption): + """Stage 2: skims -> histograms (control_shapes.root or shapes.root); other histograms in the file are kept.""" + from shapesmith.fill import run_hist cfg, analysis, selected = _setup(config, channels, workers=workers, overrides=overrides) - write_versions(cfg, cfg.output_dir) - hset = run_hist( - cfg, - analysis, - selected, - control, - variables.split(",") if variables else None, - not skip_systematics, - processes.split(",") if processes else None, - _shapes_path(cfg, control), - force, - regions.split(",") if regions else None, - ) - typer.echo(f"{len(hset)} histograms in {_shapes_path(cfg, control)}") + write_versions(cfg, analysis.name, cfg.output_dir) + output = _shapes_path(cfg, control) + hset = run_hist(cfg, analysis, selected, control, _split(variables), not skip_systematics, _split(processes), output, _split(regions)) + typer.echo(f"{len(hset)} histograms in {output}") @app.command() def estimate(config: Path = ConfigOption, channels: Optional[str] = ChannelsOption, overrides: Optional[list[str]] = SetOption, control: bool = typer.Option(False, "--control")): - """Add the estimated jet-fake process (and embedding variations) to the histogram file.""" - from shapesmith.estimators import run_estimate - from shapesmith.histograms import HistogramSet + """Stage 3: add the estimated processes and variations of every channel to the histogram file.""" + from shapesmith.estimates import run_estimates + from shapesmith.histogram import HistogramSet cfg, analysis, selected = _setup(config, channels, overrides=overrides) path = _shapes_path(cfg, control) hset = HistogramSet.load(path) - added = run_estimate(hset, analysis, selected) + added = run_estimates(hset, analysis, selected) hset.save(path) typer.echo(f"{len(added)} histograms added to {path}") +@app.command() +def measure(config: Path = ConfigOption, channels: Optional[str] = ChannelsOption, overrides: Optional[list[str]] = SetOption, suggest_binning: bool = typer.Option(False, "--suggest-binning", help="print proposed bin edges instead of measuring"), merge: bool = typer.Option(False, "--merge", help="combine the results of the earlier runs (e.g. one per working point) into one payload"), workers: Optional[int] = WorkersOption): + """Run the analysis's measurement; results in ///.""" + from shapesmith.measurements import run_measure + + cfg, analysis, selected = _setup(config, channels, workers=workers, overrides=overrides) + typer.echo(f"measurement output: {run_measure(cfg, analysis, selected, suggest_binning, merge)}") + + @app.command() def sync(config: Path = ConfigOption, channels: Optional[str] = ChannelsOption, overrides: Optional[list[str]] = SetOption): """Write combine-style shape files per channel.""" - from shapesmith.histograms import HistogramSet - from shapesmith.sync import run_sync + from shapesmith.histogram import HistogramSet + from shapesmith.shapes import run_sync cfg, analysis, selected = _setup(config, channels, overrides=overrides) for path in run_sync(HistogramSet.load(_shapes_path(cfg, False)), analysis, selected, cfg.output_dir): @@ -150,7 +137,7 @@ def sync(config: Path = ConfigOption, channels: Optional[str] = ChannelsOption, def datacards(config: Path = ConfigOption, channels: Optional[str] = ChannelsOption, overrides: Optional[list[str]] = SetOption, systematics: bool = typer.Option(True, "--systematics/--no-systematics"), min_background: Optional[float] = typer.Option(1.0, help="merge bins below this background yield; negative disables")): """Write text datacards and their shapes file per final state (each channel and the combination).""" from shapesmith.datacards import run_datacards - from shapesmith.histograms import HistogramSet + from shapesmith.histogram import HistogramSet cfg, analysis, selected = _setup(config, channels, overrides=overrides) threshold = None if min_background is None or min_background < 0 else min_background @@ -161,11 +148,11 @@ def datacards(config: Path = ConfigOption, channels: Optional[str] = ChannelsOpt @app.command() def fit(config: Path = ConfigOption, channels: Optional[str] = ChannelsOption, overrides: Optional[list[str]] = SetOption, final_states: Optional[str] = typer.Option(None, help="comma separated; default: each channel and 'all'"), skip_combine: bool = typer.Option(False, "--skip-combine", help="only collect existing combine outputs")): """Run combine (expected limits, significance, best fit) and summarise the results.""" - from shapesmith.fit import run_fit + from shapesmith.limits import run_limits cfg, _, selected = _setup(config, channels, overrides=overrides) - states = final_states.split(",") if final_states else list(_final_states(selected)) - results = run_fit(cfg, cfg.output_dir / "datacards", states, skip_combine) + states = _split(final_states) or list(_final_states(selected)) + results = run_limits(cfg, cfg.output_dir / "datacards", states, skip_combine) typer.echo((cfg.output_dir / "datacards" / "limits.md").read_text()) typer.echo(f"{len(results)} final states summarised") @@ -173,11 +160,11 @@ def fit(config: Path = ConfigOption, channels: Optional[str] = ChannelsOption, o @app.command() def plot(config: Path = ConfigOption, channels: Optional[str] = ChannelsOption, overrides: Optional[list[str]] = SetOption, control: bool = typer.Option(False, "--control"), category: Optional[str] = typer.Option(None), variables: Optional[str] = typer.Option(None), region: str = typer.Option("nominal", "--region", help="histogram region to plot"), blind: bool = typer.Option(False, "--blind"), log: bool = typer.Option(False, "--log"), signal_scale: Optional[float] = typer.Option(None), normalize_by_bin_width: bool = typer.Option(False, "--normalize-by-bin-width")): """Prefit stack plots of control variables (--control) or NN categories.""" - from shapesmith.histograms import HistogramSet + from shapesmith.histogram import HistogramSet from shapesmith.plotting.stack import run_plot cfg, analysis, selected = _setup(config, channels, overrides=overrides) - files = run_plot(HistogramSet.load(_shapes_path(cfg, control)), analysis, selected, control, category, variables.split(",") if variables else None, cfg.output_dir / "plots", region=region, blind=blind, log=log, signal_scale=signal_scale, normalize_by_bin_width=normalize_by_bin_width) + files = run_plot(HistogramSet.load(_shapes_path(cfg, control)), analysis, selected, control, category, _split(variables), cfg.output_dir / "plots", region=region, blind=blind, log=log, signal_scale=signal_scale, normalize_by_bin_width=normalize_by_bin_width) typer.echo(f"{len(files)} files written to {cfg.output_dir / 'plots'}") @@ -188,19 +175,23 @@ def ml_export(config: Path = ConfigOption, channels: Optional[str] = ChannelsOpt cfg, analysis, selected = _setup(config, channels, workers=workers, overrides=overrides) if cfg.ml_dir is not None: - write_versions(cfg, cfg.ml_dir) + write_versions(cfg, analysis.name, cfg.ml_dir) for path in run_ml_export(cfg, analysis, selected): typer.echo(f"wrote {path}") @app.command() -def inspect(shapes: Path = typer.Argument(..., help="a shapes ROOT file written by `hist`")): +def inspect(shapes: Path = typer.Argument(..., help="a shapes ROOT file written by `hist`"), unchanged: bool = typer.Option(False, "--unchanged", help="list the variations that are bitwise equal to their nominal")): """List the processes, regions, variations and variables in a histogram file.""" - from shapesmith.histograms import HistogramSet + from shapesmith.histogram import HistogramSet, unchanged_variations hset = HistogramSet.load(shapes) + if unchanged: + for key in unchanged_variations(hset): + typer.echo(key.path) + return for field in ("channel", "category", "process", "region", "variation", "variable"): - values = sorted({getattr(key, field) for key in hset.keys()}) + values = sorted({getattr(key, field) for key in hset}) typer.echo(f"{field:10s} ({len(values)}): {', '.join(values[:40])}{' ...' if len(values) > 40 else ''}") typer.echo(f"{len(hset)} histograms") diff --git a/src/shapesmith/cmssw.py b/src/shapesmith/cmssw.py new file mode 100644 index 0000000..ca0d0d6 --- /dev/null +++ b/src/shapesmith/cmssw.py @@ -0,0 +1,55 @@ +"""Run commands inside a CMSSW environment (combine, CombineHarvester) in a bash subshell. + +The subshell starts from a clean environment: the Python environment ShapeSmith runs in (e.g. an LCG view) would +otherwise shadow CMSSW's libraries and Python packages. With a timeout, the whole process group of the subshell is +killed when it expires (a fit can iterate without end). +""" +from __future__ import annotations + +import logging +import os +import signal +import subprocess +from pathlib import Path +from typing import Sequence + +from shapesmith.config import CombineConfig + +logger = logging.getLogger(__name__) +KEPT_VARIABLES = ("HOME", "USER", "LOGNAME", "TMPDIR", "X509_USER_PROXY") # the environment a CMSSW subshell starts with, besides PATH + + +def script(combine: CombineConfig, cwd: Path, commands: Sequence[str]) -> str: + """The bash script: CMSSW environment of `combine.cmssw_dir`, then the commands in `cwd`, stopping at the first failure.""" + return "\n".join( + [ + "set -e", + f"export SCRAM_ARCH={combine.scram_arch}", + "source /cvmfs/cms.cern.ch/cmsset_default.sh", + f"cd {combine.cmssw_dir}/src", + "eval $(scramv1 runtime -sh)", + f"cd {Path(cwd).resolve()}", + *commands, + ] + ) + + +def clean_environment() -> dict[str, str]: + return {"PATH": "/usr/bin:/bin", **{name: os.environ[name] for name in KEPT_VARIABLES if name in os.environ}} + + +def run(commands: Sequence[str], combine: CombineConfig | None, cwd: Path, timeout: float | None = None) -> None: + """Raise CalledProcessError on a failing command, TimeoutExpired after `timeout` seconds.""" + if combine is None: + raise ValueError("the run configuration needs `combine` (cmssw_dir) to run CMSSW commands") + logger.info(f"running {len(commands)} CMSSW commands in {cwd}") + arguments = ["bash", "-c", script(combine, cwd, commands)] + process = subprocess.Popen(arguments, env=clean_environment(), start_new_session=True) + try: + code = process.wait(timeout=timeout) + except subprocess.TimeoutExpired: + os.killpg(process.pid, signal.SIGKILL) + process.wait() + raise + if code: + raise subprocess.CalledProcessError(code, arguments) diff --git a/src/shapesmith/config.py b/src/shapesmith/config.py index a727fa3..7182d53 100644 --- a/src/shapesmith/config.py +++ b/src/shapesmith/config.py @@ -1,22 +1,16 @@ -"""Run configuration (YAML -> pydantic) and loading of the analysis object (Spec §13).""" +"""Run configuration (YAML -> pydantic) and loading of the analysis object.""" from __future__ import annotations import importlib -import json import os -import subprocess -from datetime import datetime, timezone from pathlib import Path -from typing import Any, Iterable, Literal +from typing import Any, Iterable import yaml - -from shapesmith import __version__ from pydantic import BaseModel, ConfigDict, Field -from shapesmith.model import Analysis - -SampleKind = Literal["data", "mc", "embedding"] +from shapesmith.model import Analysis, SampleKind +from shapesmith.validate import validate class FriendConfig(BaseModel): @@ -88,58 +82,7 @@ def load_analysis(config: RunConfig) -> Analysis: module_name, _, function_name = config.analysis.partition(":") if not function_name: raise ValueError(f"analysis must be 'module:function', got {config.analysis!r}") - module = importlib.import_module(module_name) - build = getattr(module, function_name) + build = getattr(importlib.import_module(module_name), function_name) analysis = build(config) - analysis.validate() + validate(analysis) return analysis - - -def git_hash(path: Path) -> str | None: - """HEAD commit of the git checkout containing `path`, or None if there is none (or git is missing).""" - path = Path(path) - directory = path.parent if path.is_file() else path - try: - return subprocess.run(["git", "-C", str(directory), "rev-parse", "HEAD"], capture_output=True, text=True, check=True).stdout.strip() - except (OSError, subprocess.CalledProcessError): - return None - - -def write_versions(config: RunConfig, directory: Path) -> Path: - """Record what produced the contents of `directory`: versions, git hashes and the full run configuration.""" - module = importlib.import_module(config.analysis.partition(":")[0]) - record = { - "written": datetime.now(timezone.utc).isoformat(timespec="seconds"), - "shapesmith": {"version": __version__, "git": git_hash(Path(__file__))}, - "analysis": {"module": config.analysis, "git": git_hash(Path(module.__file__))}, - "sample_database": {"path": str(config.sample_database), "git": git_hash(config.sample_database)} if config.sample_database else None, - "config": json.loads(config.model_dump_json()), - } - directory = Path(directory) - directory.mkdir(parents=True, exist_ok=True) - path = directory / "versions.json" - path.write_text(json.dumps(record, indent=2)) - return path - - -EXAMPLE = """\ -analysis: my_analysis.analysis:build # module:function returning an Analysis -era: "2018" -channels: [et, mt, tt] -switches: {jet_fakes: mc, embedding: false} -ntuples: - server: root://cmsdcache-kit-disk.gridka.de - base: /store/user/USER/CROWN/ntuples/TAG/CROWNRun - friends: - - {base: /store/user/USER/CROWN/ntuples/TAG/CROWNFriends/nn_v1, applies_to: [mc, data, embedding]} -skim_dir: /ceph/USER/shapesmith/TAG -output_dir: output/TAG_v1 -ml_dir: /ceph/USER/shapesmith_ml/v1 -sample_database: ../../KingMaker_sample_database/nanoAOD_v15/datasets.json # relative paths: relative to this file -workers: 16 -combine: {cmssw_dir: /work/USER/CMSSW_14_1_9, scram_arch: el9_amd64_gcc12} -""" - - -def write_example_config(path: str | Path) -> None: - Path(path).write_text(EXAMPLE) diff --git a/src/shapesmith/datacards.py b/src/shapesmith/datacards.py index 9fdbfa0..80596cd 100644 --- a/src/shapesmith/datacards.py +++ b/src/shapesmith/datacards.py @@ -1,31 +1,32 @@ -"""Native combine datacards: text card + shapes file, optional rebinning (Spec §10.2). +"""Native combine datacards: text card + shapes file, optional rebinning. -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 -weight variation present, `* autoMCStats 0`. +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`. """ from __future__ import annotations +from dataclasses import dataclass from pathlib import Path -import uproot +from shapesmith.histogram import NOMINAL_VARIATION, HistKey, Histogram, HistogramSet +from shapesmith.model import NOMINAL, Analysis, Category, LnN +from shapesmith.shapes import shape_name, write_shapes -from shapesmith.histogram import Histogram -from shapesmith.histograms import NOMINAL_REGION, NOMINAL_VARIATION, HistKey, HistogramSet -from shapesmith.model import Analysis, LnN + +@dataclass +class Bin: + channel: str + name: str + data: Histogram + processes: dict[str, Histogram] # the processes with a positive rate, rebinned, in card order + variations: dict[str, dict[str, Histogram]] # process -> variation -> rebinned histogram def bin_name(analysis: Analysis, channel: str, category_index: int) -> str: return f"htt_{channel}_{category_index + 1}_{analysis.era}" -def shape_processes(analysis: Analysis) -> tuple[str, ...]: - """Processes that carry shape systematics: simulated backgrounds and the signal.""" - excluded = {analysis.estimator.output} if analysis.estimator else set() - excluded |= {p.name for p in analysis.processes_of_kind("embedding")} - return (*(p for p in analysis.backgrounds() if p not in excluded), analysis.signal) - - def rebin_edges(total_background: Histogram, min_background: float | None) -> list[float]: """Merge bins from the right until every bin holds at least `min_background`.""" edges = list(total_background.edges) @@ -44,6 +45,33 @@ def rebin_edges(total_background: Histogram, min_background: float | None) -> li return edges +def collect_bin(hset: HistogramSet, analysis: Analysis, channel_name: str, index: int, category: Category, min_background: float | None) -> Bin | None: + """The rebinned data, processes and variations of one bin; None without data or processes.""" + channel = analysis.channel(channel_name) + variable = category.variable.name + + def key_of(process: str, variation: str = NOMINAL_VARIATION) -> HistKey: + return HistKey(channel_name, category.name, process, NOMINAL, variation, variable) + + nominal = {p: hset[key_of(p)] for p in (analysis.signal, *channel.backgrounds()) if key_of(p) in hset} + if not nominal or key_of(channel.data()) not in hset: + return None + total = None + for process, h in nominal.items(): + if process != analysis.signal: + total = h.copy() if total is None else total.add(h) + edges = rebin_edges(total, min_background) if total is not None else list(hset[key_of(channel.data())].edges) + result = Bin(channel_name, bin_name(analysis, channel_name, index), hset[key_of(channel.data())].rebin(edges), {}, {}) + for process, h in nominal.items(): + rebinned = h.rebin(edges) + if rebinned.sum() <= 0.0: + 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} + return result + + def _value(lnn: LnN) -> str: if isinstance(lnn.value, tuple): down, up = lnn.value @@ -51,91 +79,72 @@ def _value(lnn: LnN) -> str: return f"{lnn.value:g}" -def lnn_lines(analysis: Analysis, bins: list[tuple[str, str]], processes_per_bin: dict[str, list[str]]) -> list[str]: +def lnn_lines(analysis: Analysis, bins: list[Bin]) -> list[str]: """One ` lnN ...` line per concrete nuisance name ($CHANNEL/$ERA expanded), columns = (bin, process) pairs.""" lines = [] for lnn in analysis.lnn: - concrete_names = sorted({lnn.name.replace("$CHANNEL", channel).replace("$ERA", analysis.era) for channel, _ in bins}) + concrete_names = sorted({lnn.name.replace("$CHANNEL", b.channel).replace("$ERA", analysis.era) for b in bins}) for name in concrete_names: columns = [] - for channel, bin_ in bins: - this_name = lnn.name.replace("$CHANNEL", channel).replace("$ERA", analysis.era) - for process in processes_per_bin[bin_]: - applies = this_name == name and (lnn.channels is None or channel in lnn.channels) and ("*" in lnn.processes or process in lnn.processes) + for b in bins: + this_name = lnn.name.replace("$CHANNEL", b.channel).replace("$ERA", analysis.era) + for process in b.processes: + applies = this_name == name and (lnn.channels is None or b.channel in lnn.channels) and ("*" in lnn.processes or process in lnn.processes) columns.append(_value(lnn) if applies else "-") if any(column != "-" for column in columns): lines.append(f"{name} lnN " + " ".join(columns)) return lines -def write_datacard(hset: HistogramSet, analysis: Analysis, channels: list[str], output_dir: Path, systematics: bool = True, min_background: float | None = 1.0) -> Path: - output_dir = Path(output_dir) - (output_dir / "common").mkdir(parents=True, exist_ok=True) - shapes_file = output_dir / "common" / f"htt_input_{analysis.era}.root" - processes = (analysis.signal, *analysis.backgrounds()) - bins: list[tuple[str, str]] = [] - processes_per_bin: dict[str, list[str]] = {} - observations: dict[str, float] = {} - variations_per_bin: dict[str, dict[str, set[str]]] = {} - with uproot.recreate(shapes_file) as f: - for channel in channels: - for index, category in enumerate(analysis.categories): - bin_ = bin_name(analysis, channel, index) - variable = category.variable.name - - def key_of(process, variation=NOMINAL_VARIATION): - return HistKey(channel, category.name, process, NOMINAL_REGION, variation, variable) - - nominal = {p: hset.get(key_of(p)) for p in processes if hset.has(key_of(p))} - if not nominal or not hset.has(key_of("data")): - continue - total = None - for process, h in nominal.items(): - if process == analysis.signal: - continue - total = h.copy() if total is None else total.add(h) - edges = rebin_edges(total, min_background) if total is not None else list(hset.get(key_of("data")).edges) - data = hset.get(key_of("data")).rebin(edges) - f[f"{bin_}/data_obs"] = data.to_root("data_obs") - observations[bin_] = data.sum() - kept = [] - variations_per_bin[bin_] = {} - for process, h in nominal.items(): - rebinned = h.rebin(edges) - if rebinned.sum() <= 0.0: - continue - kept.append(process) - f[f"{bin_}/{process}"] = rebinned.to_root(process) - variations_per_bin[bin_][process] = set() - for key in hset.keys(channel=channel, category=category.name, process=process, region=NOMINAL_REGION, variable=variable): - if key.variation != NOMINAL_VARIATION: - f[f"{bin_}/{process}_{key.variation}"] = hset.get(key).rebin(edges).to_root(f"{process}_{key.variation}") - variations_per_bin[bin_][process].add(key.variation) - bins.append((channel, bin_)) - processes_per_bin[bin_] = kept +def shape_lines(bins: list[Bin]) -> list[str]: + columns = [(b, process) for b in bins for process in b.processes] + names = sorted({v[: -len("Up")] for b in bins for variations in b.variations.values() for v in variations if v.endswith("Up")}) + lines = [] + for name in names: + entries = ["1" if {f"{name}Up", f"{name}Down"} <= set(b.variations[process]) else "-" for b, process in columns] + if "1" in entries: + lines.append(f"{name} shape " + " ".join(entries)) + return lines + + +def card_text(analysis: Analysis, bins: list[Bin], systematics: bool) -> str: + backgrounds = list(dict.fromkeys(p for b in bins for p in analysis.channel(b.channel).backgrounds())) + index = {name: i for i, name in enumerate((analysis.signal, *backgrounds))} + columns = [(b.name, process) for b in bins for process in b.processes] lines = ["imax * number of bins", "jmax * number of processes minus 1", "kmax * number of nuisance parameters", "-" * 60] - for _, bin_ in bins: - lines.append(f"shapes * {bin_} common/htt_input_{analysis.era}.root {bin_}/$PROCESS {bin_}/$PROCESS_$SYSTEMATIC") + lines += [f"shapes * {b.name} common/htt_input_{analysis.era}.root {b.name}/$PROCESS {b.name}/$PROCESS_$SYSTEMATIC" for b in bins] lines.append("-" * 60) - lines.append("bin " + " ".join(bin_ for _, bin_ in bins)) - lines.append("observation " + " ".join(f"{observations[bin_]:.1f}" for _, bin_ in bins)) + lines.append("bin " + " ".join(b.name for b in bins)) + lines.append("observation " + " ".join(f"{b.data.sum():.1f}" for b in bins)) lines.append("-" * 60) - columns = [(bin_, process) for _, bin_ in bins for process in processes_per_bin[bin_]] - lines.append("bin " + " ".join(bin_ for bin_, _ in columns)) + lines.append("bin " + " ".join(name for name, _ in columns)) lines.append("process " + " ".join(process for _, process in columns)) - lines.append("process " + " ".join("0" if process == analysis.signal else str(processes.index(process)) for _, process in columns)) + lines.append("process " + " ".join(str(index[process]) for _, process in columns)) lines.append("rate " + " ".join("-1" for _ in columns)) lines.append("-" * 60) if systematics: - lines += lnn_lines(analysis, bins, processes_per_bin) - systematic_names = sorted({v[: -len("Up")] for per_bin in variations_per_bin.values() for vs in per_bin.values() for v in vs if v.endswith("Up")}) - for name in systematic_names: - entries = ["1" if {f"{name}Up", f"{name}Down"} <= variations_per_bin[bin_].get(process, set()) else "-" for bin_, process in columns] - if any(entry == "1" for entry in entries): - lines.append(f"{name} shape " + " ".join(entries)) + lines += lnn_lines(analysis, bins) + shape_lines(bins) lines.append("* autoMCStats 0") + return "\n".join(lines) + "\n" + + +def write_datacard(hset: HistogramSet, analysis: Analysis, channels: list[str], output_dir: Path, systematics: bool = True, min_background: float | None = 1.0) -> Path: + output_dir = Path(output_dir) + bins = [ + b + for channel in channels + for index, category in enumerate(analysis.channel(channel).categories) + if (b := collect_bin(hset, analysis, channel, index, category, min_background)) is not None + ] + entries = [] + for b in bins: + entries.append((b.name, "data_obs", b.data)) + for process, h in b.processes.items(): + entries.append((b.name, process, h)) + entries += [(b.name, shape_name(process, variation), varied) for variation, varied in b.variations[process].items()] + write_shapes(output_dir / "common" / f"htt_input_{analysis.era}.root", entries) card = output_dir / "combined.txt" - card.write_text("\n".join(lines) + "\n") + card.write_text(card_text(analysis, bins, systematics)) return card diff --git a/src/shapesmith/estimates.py b/src/shapesmith/estimates.py new file mode 100644 index 0000000..6dc4040 --- /dev/null +++ b/src/shapesmith/estimates.py @@ -0,0 +1,121 @@ +"""Stage 3: processes estimated from histograms, run per channel in the order of `Channel.estimators`. + +DataMinus: data minus the subtracted processes in its region, for the nominal and every column variation filled on +data or a subtracted process there; an input without that variation contributes its nominal (under the fill rule it +is unchanged). Weight variations are not propagated. Optionally negative bins are clipped. +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`. +""" +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 + +logger = logging.getLogger(__name__) + + +def data_minus(hset: HistogramSet, channel: Channel, category: str, variable: str, region: str, subtract: tuple[str, ...], variation: str = NOMINAL_VARIATION) -> Histogram | None: + """data[region] - sum of the given processes[region] under `variation` (nominal where an input lacks it); + None if the data histogram is missing.""" + + def get(process: str) -> Histogram | None: + key = HistKey(channel.name, category, process, region, variation, variable) + if key not in hset: + key = HistKey(channel.name, category, process, region, NOMINAL_VARIATION, variable) + return hset.get(key) + + data = get(channel.data()) + if data is None: + return None + result = data.copy() + for process in subtract: + h = get(process) + if h is None: + logger.warning(f"{channel.name}_{category}/{process}#{region}#{variable} missing, not subtracted") + continue + result.add(h, -1.0) + return result + + +def categories_and_variables(hset: HistogramSet, channel: Channel, region: str) -> list[tuple[str, str]]: + return sorted({(k.category, k.variable) for k in hset.select(channel=channel.name, process=channel.data(), region=region, variation=NOMINAL_VARIATION)}) + + +def estimate_data_minus(hset: HistogramSet, channel: Channel, estimator: DataMinus) -> list[HistKey]: + column_variations = {v.name for v in channel.variations if isinstance(v, ColumnVariation)} + inputs = (channel.data(), *estimator.subtract) + added = [] + for category, variable in categories_and_variables(hset, channel, estimator.region): + present = {k.variation for process in inputs for k in hset.select(channel=channel.name, category=category, process=process, region=estimator.region, variable=variable)} + for variation in (NOMINAL_VARIATION, *sorted(present & column_variations)): + h = data_minus(hset, channel, category, variable, estimator.region, estimator.subtract, variation).scale(estimator.scale) + key = HistKey(channel.name, category, estimator.output, NOMINAL, variation, variable) + hset[key] = clip_negative_bins(h) if estimator.clip_negative else h + added.append(key) + return added + + +def clip_negative_bins(h: Histogram) -> Histogram: + """Set negative bins to zero and rescale to the original integral (zero everything if the integral is <= 0).""" + total = h.sum() + positive = float(h.values[h.values > 0].sum()) + h.values[h.values < 0] = 0.0 + if total <= 0.0: + logger.warning(f"integral {total:.3f} <= 0 after subtraction, histogram set to zero") + h.values[:] = 0.0 + h.variances[:] = 0.0 + elif positive > 0.0 and positive != total: + h.scale(total / positive) + return h + + +def estimate_abcd(hset: HistogramSet, channel: Channel, estimator: ABCD) -> list[HistKey]: + added = [] + for category, variable in categories_and_variables(hset, channel, estimator.b): + b, c, d = (data_minus(hset, channel, category, variable, region, estimator.subtract) for region in (estimator.b, estimator.c, estimator.d)) + if b is None or c is None or d is None: + logger.warning(f"{channel.name}/{category}/{variable}: ABCD control region missing, skipped") + continue + yield_c, yield_d = c.sum(), d.sum() + if yield_c <= 0.0 or yield_d <= 0.0: + logger.warning(f"{channel.name}/{category}/{variable}: C={yield_c:.2f}, D={yield_d:.2f}, {estimator.output} set to zero") + factor = 0.0 + else: + factor = yield_c / yield_d + key = HistKey(channel.name, category, estimator.output, NOMINAL, NOMINAL_VARIATION, variable) + hset[key] = clip_negative_bins(b.scale(factor)) + added.append(key) + return added + + +def estimate_template_shift(hset: HistogramSet, channel: Channel, estimator: TemplateShift) -> list[HistKey]: + added = [] + for key in hset.select(channel=channel.name, process=estimator.process, region=NOMINAL): + if key.variation != NOMINAL_VARIATION and not is_template(key.variation): + continue + template = HistKey(channel.name, key.category, estimator.template, NOMINAL, NOMINAL_VARIATION, key.variable) + if template not in hset: + logger.warning(f"{template.path} missing, no {estimator.name} variation") + continue + for direction, sign in (("Up", 1.0), ("Down", -1.0)): + name = f"{estimator.name}{direction}" + variation = name if key.variation == NOMINAL_VARIATION else on_template(name, key.variation) + varied = HistKey(channel.name, key.category, estimator.process, NOMINAL, variation, key.variable) + hset[varied] = hset[key].copy().add(hset[template], sign * estimator.fraction) + added.append(varied) + return added + + +ESTIMATES = {DataMinus: estimate_data_minus, ABCD: estimate_abcd, TemplateShift: estimate_template_shift} + + +def run_estimates(hset: HistogramSet, analysis: Analysis, channels: list[str]) -> list[HistKey]: + added: list[HistKey] = [] + for name in channels: + channel = analysis.channel(name) + for estimator in channel.estimators: + added += ESTIMATES[type(estimator)](hset, channel, estimator) + return added diff --git a/src/shapesmith/estimators/__init__.py b/src/shapesmith/estimators/__init__.py deleted file mode 100644 index 3bc6420..0000000 --- a/src/shapesmith/estimators/__init__.py +++ /dev/null @@ -1,42 +0,0 @@ -"""Derived processes computed from histograms of the estimation regions (Spec §9).""" -from __future__ import annotations - -import logging - -from shapesmith.histogram import Histogram -from shapesmith.histograms import HistKey, HistogramSet -from shapesmith.model import Analysis - -logger = logging.getLogger(__name__) - - -def data_minus(hset: HistogramSet, channel: str, category: str, variable: str, region: str, subtract: tuple[str, ...]) -> Histogram | None: - """data[region] - sum of the given processes[region]; None if the data histogram is missing.""" - data_key = HistKey(channel, category, "data", region, "Nominal", variable) - if not hset.has(data_key): - return None - result = hset.get(data_key).copy() - for process in subtract: - key = HistKey(channel, category, process, region, "Nominal", variable) - if not hset.has(key): - logger.warning(f"{key.path} missing, not subtracted") - continue - result.add(hset.get(key), -1.0) - return result - - -def categories_and_variables(hset: HistogramSet, channel: str, region: str) -> list[tuple[str, str]]: - return sorted({(k.category, k.variable) for k in hset.keys(channel=channel, process="data", region=region)}) - - -def run_estimate(hset: HistogramSet, analysis: Analysis, channels: list[str]) -> list[HistKey]: - from shapesmith.estimators import abcd, embedding, fake_factors # local import: the modules import this package - - added: list[HistKey] = [] - for channel in channels: - if analysis.estimator is not None: - estimate = {"fake_factors": fake_factors.estimate, "abcd": abcd.estimate}[analysis.estimator.name] - added += estimate(hset, analysis, channel) - if analysis.processes_of_kind("embedding"): - added += embedding.ttbar_contamination(hset, analysis, channel) - return added diff --git a/src/shapesmith/estimators/abcd.py b/src/shapesmith/estimators/abcd.py deleted file mode 100644 index 744c4b0..0000000 --- a/src/shapesmith/estimators/abcd.py +++ /dev/null @@ -1,51 +0,0 @@ -"""ABCD method: QCD(A) = (data - MC)(B) * (data - MC)(C).sum() / (data - MC)(D).sum(). - -A = signal region (OS, isolated), B = OS anti-isolated, C = SS isolated, D = SS anti-isolated. -""" -from __future__ import annotations - -import logging - -from shapesmith.estimators import categories_and_variables, data_minus -from shapesmith.histogram import Histogram -from shapesmith.histograms import NOMINAL_REGION, NOMINAL_VARIATION, HistKey, HistogramSet -from shapesmith.model import Analysis - -logger = logging.getLogger(__name__) - - -def clip_negative_bins(h: Histogram) -> Histogram: - """Set negative bins to zero and rescale to the original integral (zero everything if the integral is <= 0).""" - total = h.sum() - positive = float(h.values[h.values > 0].sum()) - h.values[h.values < 0] = 0.0 - if total <= 0.0: - logger.warning(f"integral {total:.3f} <= 0 after subtraction, histogram set to zero") - h.values[:] = 0.0 - h.variances[:] = 0.0 - elif positive > 0.0 and positive != total: - h.scale(total / positive) - return h - - -def estimate(hset: HistogramSet, analysis: Analysis, channel: str) -> list[HistKey]: - estimator = analysis.estimator - regions = estimator.regions - added = [] - for category, variable in categories_and_variables(hset, channel, regions["B"]): - b = data_minus(hset, channel, category, variable, regions["B"], estimator.subtract) - c = data_minus(hset, channel, category, variable, regions["C"], estimator.subtract) - d = data_minus(hset, channel, category, variable, regions["D"], estimator.subtract) - if b is None or c is None or d is None: - logger.warning(f"{channel}/{category}/{variable}: ABCD control region missing, skipped") - continue - yield_c, yield_d = c.sum(), d.sum() - if yield_c <= 0.0 or yield_d <= 0.0: - logger.warning(f"{channel}/{category}/{variable}: C={yield_c:.2f}, D={yield_d:.2f}, QCD set to zero") - factor = 0.0 - else: - factor = yield_c / yield_d - key = HistKey(channel, category, estimator.output, NOMINAL_REGION, NOMINAL_VARIATION, variable) - hset.add(key, clip_negative_bins(b.scale(factor))) - added.append(key) - return added diff --git a/src/shapesmith/estimators/embedding.py b/src/shapesmith/estimators/embedding.py deleted file mode 100644 index 199f353..0000000 --- a/src/shapesmith/estimators/embedding.py +++ /dev/null @@ -1,27 +0,0 @@ -"""ttbar contamination variation of the embedded sample: EMB -/+ 10 % of the genuine-tau ttbar histogram.""" -from __future__ import annotations - -import logging - -from shapesmith.histograms import NOMINAL_REGION, NOMINAL_VARIATION, HistKey, HistogramSet -from shapesmith.model import Analysis - -logger = logging.getLogger(__name__) -TTBAR_GENUINE = "TTT" -FRACTION = 0.1 - - -def ttbar_contamination(hset: HistogramSet, analysis: Analysis, channel: str) -> list[HistKey]: - added = [] - for emb in analysis.processes_of_kind("embedding"): - for key in hset.keys(channel=channel, process=emb.name, region=NOMINAL_REGION, variation=NOMINAL_VARIATION): - ttbar_key = HistKey(channel, key.category, TTBAR_GENUINE, NOMINAL_REGION, NOMINAL_VARIATION, key.variable) - if not hset.has(ttbar_key): - logger.warning(f"{ttbar_key.path} missing, no embedding ttbar variation") - continue - for direction, sign in (("Down", -1.0), ("Up", +1.0)): - varied = hset.get(key).copy().add(hset.get(ttbar_key), sign * FRACTION) - new_key = HistKey(channel, key.category, emb.name, NOMINAL_REGION, f"CMS_htt_emb_ttbar_{analysis.era}{direction}", key.variable) - hset.add(new_key, varied) - added.append(new_key) - return added diff --git a/src/shapesmith/estimators/fake_factors.py b/src/shapesmith/estimators/fake_factors.py deleted file mode 100644 index c4ae1ef..0000000 --- a/src/shapesmith/estimators/fake_factors.py +++ /dev/null @@ -1,20 +0,0 @@ -"""Fake-factor method: jetFakes = FF-weighted anti-isolated data minus genuine-tau and lepton-fake MC.""" -from __future__ import annotations - -from shapesmith.estimators import categories_and_variables, data_minus -from shapesmith.histograms import NOMINAL_REGION, NOMINAL_VARIATION, HistKey, HistogramSet -from shapesmith.model import Analysis - - -def estimate(hset: HistogramSet, analysis: Analysis, channel: str) -> list[HistKey]: - estimator = analysis.estimator - region = estimator.regions["anti_iso"] - added = [] - for category, variable in categories_and_variables(hset, channel, region): - h = data_minus(hset, channel, category, variable, region, estimator.subtract) - if h is None: - continue - key = HistKey(channel, category, estimator.output, NOMINAL_REGION, NOMINAL_VARIATION, variable) - hset.add(key, h) - added.append(key) - return added diff --git a/src/shapesmith/events.py b/src/shapesmith/events.py new file mode 100644 index 0000000..da89cb3 --- /dev/null +++ b/src/shapesmith/events.py @@ -0,0 +1,74 @@ +"""The one rule that turns a process in a region under a variation into cuts and weights, and loads its events. + +Cuts: the channel cuts with the region's replacements, the process cuts and the skim cuts. The skim cuts are a no-op +on nominal rows; under a variation they keep out the events a wider skim stored for other variations. +Weights, in order: the process weights, with the region's replacements where the process carries the name, then the +region's added weights. A weight variation replaces weights in place and does not apply to a process lacking one of +them; a column variation rewrites every cut and weight (expressions.shift). +Event weight: (norm_weight * lumi) * product of the weights, with lumi only for simulated (mc) processes. +""" +from __future__ import annotations + +from dataclasses import dataclass +from typing import Iterable + +import numpy as np +import pandas as pd + +from shapesmith.config import RunConfig +from shapesmith.expressions import columns_of, mask, product, shift +from shapesmith.model import NOMINAL, Analysis, Channel, ColumnVariation, Process, Region, Variation, WeightVariation +from shapesmith.store import SKIM_COLUMNS, read_skims, schema + + +@dataclass(frozen=True) +class Query: + channel: str + process: str + region: str = NOMINAL + variation: Variation | None = None + + +@dataclass +class Events: + frame: pd.DataFrame + weights: np.ndarray + + +def select(channel: Channel, process: Process, region: Region, variation: Variation | None = None, available: set[str] = frozenset()) -> tuple[list[str], list[str]] | None: + """Cut and weight expressions of `process` in `region` under `variation`; None if the variation does not apply.""" + weights = {name: region.replace_weights.get(name, expr) for name, expr in process.selection.weights.items()} + weights.update(region.add_weights) + if isinstance(variation, WeightVariation): + if any(name not in weights for name in variation.replace_weights): + return None + weights.update(variation.replace_weights) + cuts = [*{**channel.cuts, **region.replace_cuts}.values(), *process.selection.cuts.values(), *channel.skim.values()] + weights = list(weights.values()) + if isinstance(variation, ColumnVariation): + cuts = [shift(expr, variation, available) for expr in cuts] + weights = [shift(expr, variation, available) for expr in weights] + return cuts, weights + + +def lumi(analysis: Analysis, channel: Channel, process: Process) -> float: + return analysis.lumi_pb if channel.kind_of(process) == "mc" else 1.0 + + +def event_weights(frame: pd.DataFrame, weights: Iterable[str], lumi_pb: float) -> np.ndarray: + return frame["norm_weight"].to_numpy() * lumi_pb * product(frame, weights) + + +def load(config: RunConfig, analysis: Analysis, query: Query, columns: Iterable[str] = ()) -> Events: + """The selected events of a query (one frame, in sample order) with their weights; `columns` are read in addition.""" + channel = analysis.channel(query.channel) + process = channel.process(query.process) + nicks = [s.nick for s in channel.samples_of(process.group)] + available = schema(config.skim_dir, channel.name, nicks) if isinstance(query.variation, ColumnVariation) else set() + selection = select(channel, process, channel.region(query.region), query.variation, available) + if selection is None: + raise ValueError(f"{query.channel}/{query.process}: variation {query.variation.name} does not apply") + cuts, weights = selection + frame = read_skims(config.skim_dir, channel.name, nicks, set(SKIM_COLUMNS) | columns_of(cuts + weights) | set(columns)) + selected = frame[mask(frame, cuts)].reset_index(drop=True) + return Events(selected, event_weights(selected, weights, lumi(analysis, channel, process))) diff --git a/src/shapesmith/expressions.py b/src/shapesmith/expressions.py index cd5c79b..9128b3c 100644 --- a/src/shapesmith/expressions.py +++ b/src/shapesmith/expressions.py @@ -1,15 +1,14 @@ -"""Evaluation of cut and weight expressions (pandas.eval syntax, numexpr engine), Spec §6.""" +"""Evaluation of cut and weight expressions (pandas.eval syntax, numexpr engine) and their column variations.""" from __future__ import annotations -import re -from typing import Mapping - import logging +import re +from typing import Iterable import numpy as np import pandas as pd -from shapesmith.model import Region, Selection, WeightVariation +from shapesmith.model import ColumnVariation FUNCTIONS = frozenset( "abs sqrt exp log log10 sin cos tan arctan2 arcsin arccos arctan sinh cosh tanh expm1 log1p".split() @@ -17,7 +16,6 @@ KEYWORDS = frozenset({"and", "or", "not", "True", "False", "inf", "nan"}) _IDENTIFIER = re.compile(r"(? set[str]: return {token for token in _IDENTIFIER.findall(expr) if token not in FUNCTIONS and token not in KEYWORDS} +def columns_of(exprs: Iterable[str]) -> set[str]: + columns: set[str] = set() + for expr in set(exprs): # the selections of many regions and variations repeat the same expressions + columns |= columns_in(expr) + return columns + + def evaluate(frame: pd.DataFrame, expr: str) -> np.ndarray: """Evaluate `expr` on `frame`; scalars are broadcast to one value per row.""" missing = sorted(columns_in(expr) - set(frame.columns)) @@ -50,57 +55,31 @@ def evaluate(frame: pd.DataFrame, expr: str) -> np.ndarray: return np.asarray(result) -def mask(frame: pd.DataFrame, cuts: Mapping[str, str]) -> np.ndarray: - """Logical AND of all cuts (all True for an empty mapping).""" +def mask(frame: pd.DataFrame, cuts: Iterable[str]) -> np.ndarray: + """Logical AND of all cuts (all True for none).""" result = np.ones(len(frame), dtype=bool) - for expr in cuts.values(): + for expr in cuts: result &= evaluate(frame, expr).astype(bool) return result -def weight(frame: pd.DataFrame, weights: Mapping[str, str]) -> np.ndarray: - """Product of all weights as float64 (all 1.0 for an empty mapping).""" +def product(frame: pd.DataFrame, weights: Iterable[str]) -> np.ndarray: + """Product of all weights as float64, in the given order (all 1.0 for none).""" result = np.ones(len(frame), dtype=np.float64) - for expr in weights.values(): + for expr in weights: result *= evaluate(frame, expr).astype(np.float64) return result -def apply_region(selection: Selection, region: Region) -> Selection: - """Replace the named cuts of `selection` and append the region weights.""" - cuts = dict(selection.cuts) - for name, expr in region.replace_cuts.items(): - if name not in cuts: - raise KeyError(f"region {region.name}: cut {name!r} not in selection {sorted(cuts)}") - cuts[name] = expr - weights = dict(selection.weights) - weights.update(region.add_weights) - return Selection(cuts=cuts, weights=weights) - - -def apply_weight_variation(selection: Selection, variation: WeightVariation) -> Selection: - """Replace the named weights of `selection` (weights not present are an error).""" - weights = dict(selection.weights) - for name, expr in variation.replace_weights.items(): - if name not in weights: - raise KeyError(f"variation {variation.name}: weight {name!r} not in selection {sorted(weights)}") - weights[name] = expr - return Selection(cuts=dict(selection.cuts), weights=weights) - - -def apply_column_variation(expr: str, suffix: str, available: set[str]) -> str: - """Rename every column `c` in `expr` to `c + suffix` if that shifted column exists.""" +def shift(expr: str, variation: ColumnVariation, available: set[str]) -> str: + """`expr` under a column variation: a derived column becomes its expression, a column `c` whose shifted branch + `c + suffix` is in `available` is read shifted, and every other column stays nominal.""" def rename(match: re.Match) -> str: name = match.group(0) - shifted = name + suffix - return shifted if name in available and shifted in available else name + if name in variation.derived: + return f"({variation.derived[name]})" + shifted = name + variation.suffix + return shifted if variation.suffix and shifted in available else name return _IDENTIFIER.sub(rename, expr) - - -def selection_columns(selection: Selection) -> set[str]: - columns: set[str] = set() - for expr in list(selection.cuts.values()) + list(selection.weights.values()): - columns |= columns_in(expr) - return columns diff --git a/src/shapesmith/fill.py b/src/shapesmith/fill.py new file mode 100644 index 0000000..0568538 --- /dev/null +++ b/src/shapesmith/fill.py @@ -0,0 +1,184 @@ +"""Stage 2: skims -> histograms per process, region, category, variable and variation. + +Bookings: the signal in the nominal region, auxiliary processes in the nominal region without variations, every +other process in the nominal region and the estimator regions (or the explicitly requested regions). A variation is +filled for the processes whose sample kind it applies to, in its regions; a column variation only where it changes +an expression of the histogram (elsewhere it equals the nominal, and consumers fall back to that). +""" +from __future__ import annotations + +import logging +from dataclasses import dataclass +from pathlib import Path +from typing import NamedTuple + +import numpy as np + +from shapesmith.config import RunConfig +from shapesmith.events import event_weights, lumi, select +from shapesmith.expressions import columns_of, evaluate, mask, shift +from shapesmith.histogram import INCLUSIVE, NOMINAL_VARIATION, HistKey, Histogram, HistogramSet +from shapesmith.model import ABCD, NOMINAL, Analysis, Channel, ColumnVariation, DataMinus, Process, Variable, Variation, applies +from shapesmith.parallel import run_jobs +from shapesmith.store import SKIM_COLUMNS, read_skims, schema + +logger = logging.getLogger(__name__) + + +@dataclass(frozen=True) +class Booking: + process: Process + regions: tuple[str, ...] + variations: tuple[Variation, ...] + + +@dataclass(frozen=True) +class Target: + category: str + cut: str | None + variable: Variable + + +def estimator_regions(channel: Channel) -> tuple[str, ...]: + regions = [] + for estimator in channel.estimators: + if isinstance(estimator, DataMinus): + regions.append(estimator.region) + elif isinstance(estimator, ABCD): + regions += [estimator.b, estimator.c, estimator.d] + return tuple(dict.fromkeys(regions)) + + +def resolve_regions(channel: Channel, regions: list[str] | None) -> tuple[str, ...] | None: + """An explicit region request; `all` expands to the nominal and every region of the channel.""" + if regions is None: + return None + names = (NOMINAL, *(region.name for region in channel.regions)) if "all" in regions else tuple(regions) + resolved = tuple(dict.fromkeys(names)) + for name in resolved: + channel.region(name) # validates with a channel-specific error + return resolved + + +def bookings(channel: Channel, systematics: bool, regions: list[str] | None = None) -> list[Booking]: + requested = resolve_regions(channel, regions) + result = [] + for process in channel.processes: + if process.role == "auxiliary": + if requested is None or NOMINAL in requested: + result.append(Booking(process, (NOMINAL,), ())) + continue + if requested is not None: + process_regions = requested + elif process.role == "signal": + process_regions = (NOMINAL,) + else: + process_regions = tuple(dict.fromkeys((NOMINAL, *estimator_regions(channel)))) + kind = channel.kind_of(process) + variations = tuple(v for v in channel.variations if applies(v, kind, process.group)) if systematics else () + result.append(Booking(process, process_regions, variations)) + return result + + +def targets(channel: Channel, control: bool, variables: list[str] | None) -> list[Target]: + if control: + names = variables or list(channel.variables) + return [Target(INCLUSIVE, None, channel.variables[name]) for name in names] # KeyError for unknown variables + return [Target(c.name, c.cut, c.variable) for c in channel.categories] + + +class Fill(NamedTuple): + """The histograms of one region and variation of a booking: its cuts and weights, and per target the category + cut and the variable expression.""" + + region: str + variation: Variation | None + cuts: list[str] + weights: list[str] + targets: list[tuple[Target, str | None, str]] + + +def _target_exprs(target: Target, variation: Variation | None, available: set[str]) -> tuple[str | None, str]: + if not isinstance(variation, ColumnVariation): + return target.cut, target.variable.expr + return (shift(target.cut, variation, available) if target.cut else None), shift(target.variable.expr, variation, available) + + +def plan_booking(channel: Channel, booking: Booking, targets_: list[Target], available: set[str]) -> list[Fill]: + """What to fill: every region with its applicable variations; a column variation only for the targets whose + cuts, weights, category or variable it changes.""" + fills = [] + for region_name in booking.regions: + region = channel.region(region_name) + nominal = select(channel, booking.process, region) + for variation in (None, *(v for v in booking.variations if v.regions is None or region_name in v.regions)): + selection = select(channel, booking.process, region, variation, available) + if selection is None: + logger.debug(f"{booking.process.name}: variation {variation.name} not applicable, skipped") + continue + exprs = [(target, *_target_exprs(target, variation, available)) for target in targets_] + if isinstance(variation, ColumnVariation) and selection == nominal: + exprs = [(target, cut, expr) for target, cut, expr in exprs if (cut, expr) != (target.cut, target.variable.expr)] + if exprs: + fills.append(Fill(region_name, variation, *selection, exprs)) + return fills + + +def fill_booking(config: RunConfig, analysis: Analysis, channel_name: str, booking: Booking, targets_: list[Target]) -> dict[HistKey, Histogram]: + """All histograms of one booking, from one frame of the process's skims.""" + channel = analysis.channel(channel_name) + process = booking.process + nicks = [s.nick for s in channel.samples_of(process.group)] + fills = plan_booking(channel, booking, targets_, schema(config.skim_dir, channel_name, nicks)) + exprs = [] + for f in fills: + exprs += f.cuts + f.weights + [e for _, cut, var in f.targets for e in (cut, var) if e] + frame = read_skims(config.skim_dir, channel_name, nicks, set(SKIM_COLUMNS) | columns_of(exprs)) + lumi_pb = lumi(analysis, channel, process) + result: dict[HistKey, Histogram] = {} + for f in fills: + selected = frame[mask(frame, f.cuts)] + w = event_weights(selected, f.weights, lumi_pb) + for target, cut, expr in f.targets: + in_category = mask(selected, [cut]) if cut else np.ones(len(selected), dtype=bool) + values = evaluate(selected, expr).astype(np.float64) + finite = np.isfinite(values) & np.isfinite(w) + if (in_category & ~finite).any(): + logger.warning(f"{process.name}/{f.region}/{target.variable.name}: {(in_category & ~finite).sum()} events with NaN/inf skipped") + keep = in_category & finite + key = HistKey(channel_name, target.category, process.name, f.region, f.variation.name if f.variation else NOMINAL_VARIATION, target.variable.name) + result[key] = Histogram.fill(target.variable.edges, values[keep], w[keep]) + return result + + +def _fill_job(args) -> dict[HistKey, Histogram]: + return fill_booking(*args) + + +def run_hist(config: RunConfig, analysis: Analysis, channels: list[str] | None, control: bool, variables: list[str] | None, systematics: bool, processes: list[str] | None, output: Path, regions: list[str] | None = None) -> HistogramSet: + """Fill every booked histogram of the requested channels into `output` (+ .json index). + + The requested scopes (channel, process, region, category, variable) are always filled; histograms of other + scopes already in `output` are kept.""" + output = Path(output) + jobs, scopes = [], set() + for channel_name in channels or config.channels: + channel = analysis.channel(channel_name) + targets_ = targets(channel, control, variables) + for booking in bookings(channel, systematics, regions): + if processes and booking.process.name not in processes: + continue + jobs.append((config, analysis, channel_name, booking, targets_)) + scopes |= {(channel_name, booking.process.name, region, t.category, t.variable.name) for region in booking.regions for t in targets_} + hset = HistogramSet() + if output.exists(): + previous = HistogramSet.load(output) + hset.update((key, h) for key, h in previous.items() if (key.channel, key.process, key.region, key.category, key.variable) not in scopes) + logger.info(f"{output}: filling {len(scopes)} histogram scopes, keeping {len(hset)} other histograms") + logger.info(f"filling {len(jobs)} process bookings with {config.workers} workers") + for job, result, error in run_jobs(_fill_job, jobs, config.workers): + if error is not None: + raise RuntimeError(f"{job[2]}/{job[3].process.name}: {error}") from error + hset.update(result) + hset.save(output) + return hset diff --git a/src/shapesmith/histogram.py b/src/shapesmith/histogram.py index 991b1bc..a4a9a12 100644 --- a/src/shapesmith/histogram.py +++ b/src/shapesmith/histogram.py @@ -1,4 +1,5 @@ -"""A minimal weighted 1D histogram on numpy arrays with ROOT TH1D conversion through uproot. +"""A minimal weighted 1D histogram on numpy arrays with ROOT TH1D conversion through uproot, and the keyed set +of histograms the stages exchange (one ROOT file plus a JSON index). boost-histogram/hist are deliberately not used: in LCG_108 (boost_histogram 1.3.2, numpy 2.1) the axis `edges` property returns constant values, and uproot writes exactly those corrupted edges. numpy plus @@ -6,11 +7,29 @@ """ from __future__ import annotations -from dataclasses import dataclass +import json +from dataclasses import asdict, dataclass +from pathlib import Path import numpy as np +import uproot from uproot.writing.identify import to_TAxis, to_TH1x +NOMINAL_VARIATION = "Nominal" +INCLUSIVE = "inclusive" # the category of control-variable histograms +TEMPLATE_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 + + +def on_template(variation: str, template: str) -> str: + """The name of `variation` applied on top of the template variation `template`. It ends in the template's name, + so datacards never read it as a shape.""" + return f"{variation}{TEMPLATE_SEPARATOR}{template}" + @dataclass class Histogram: @@ -76,3 +95,72 @@ def to_root(self, name: str): @classmethod def from_root(cls, th) -> "Histogram": return cls(th.axis().edges(), th.values(), th.variances()) + + +@dataclass(frozen=True, order=True) +class HistKey: + channel: str + category: str + process: str + region: str + variation: str + variable: str + + @property + def directory(self) -> str: + return f"{self.channel}_{self.category}" + + @property + def object_name(self) -> str: + return f"{self.process}#{self.region}#{self.variation}#{self.variable}" + + @property + def path(self) -> str: + return f"{self.directory}/{self.object_name}" + + @classmethod + def parse(cls, path: str) -> "HistKey": + directory, name = path.split("/", 1) + channel, category = directory.split("_", 1) + process, region, variation, variable = name.split("#") + return cls(channel, category, process, region, variation, variable) + + +class HistogramSet(dict): + """HistKey -> Histogram, saved as one ROOT file (`_/###`) + plus a JSON index next to it.""" + + def select(self, **fields) -> list[HistKey]: + """The keys whose fields equal the given values, sorted.""" + return sorted(key for key in self if all(getattr(key, name) == value for name, value in fields.items())) + + def save(self, path: Path) -> None: + path = Path(path) + path.parent.mkdir(parents=True, exist_ok=True) + items = sorted(self.items()) + with uproot.recreate(path) as f: + for key, h in items: + f[key.path] = h.to_root(key.object_name) + index = [{**asdict(key), "sum": h.sum(), "bins": len(h.values)} for key, h in items] + path.with_suffix(".json").write_text(json.dumps(index, indent=1)) + + @classmethod + def load(cls, path: Path) -> "HistogramSet": + hset = cls() + with uproot.open(path) as f: + for name, classname in f.classnames(recursive=True, cycle=False).items(): + if classname.startswith("TH1"): + hset[HistKey.parse(name)] = Histogram.from_root(f[name]) + return hset + + +def unchanged_variations(hset: HistogramSet) -> list[HistKey]: + """The variation histograms that are bitwise equal to their nominal (a declared variation that has no effect).""" + result = [] + for key in hset.select(): + if key.variation == NOMINAL_VARIATION: + continue + nominal = hset.get(HistKey(key.channel, key.category, key.process, key.region, NOMINAL_VARIATION, key.variable)) + if nominal is not None and np.array_equal(hset[key].values, nominal.values) and np.array_equal(hset[key].variances, nominal.variances): + result.append(key) + return result diff --git a/src/shapesmith/histograms.py b/src/shapesmith/histograms.py deleted file mode 100644 index 29f5d02..0000000 --- a/src/shapesmith/histograms.py +++ /dev/null @@ -1,273 +0,0 @@ -"""Stage 2: skims -> histograms per process, region, category, variable and variation (Spec §8). - -Weight rules: - * MC: norm_weight * lumi_pb * baseline weights * process weights * region weights - * embedding: norm_weight (=1) * process weights * region weights - * data: region weights only (e.g. the fake factor in the anti-isolated region) -A weight variation naming a weight the process does not carry is skipped for that process. -""" -from __future__ import annotations - -import json -import logging -from dataclasses import asdict, dataclass -from pathlib import Path -from typing import Iterator - -import numpy as np -import uproot - -from shapesmith.config import RunConfig -from shapesmith.expressions import apply_region, columns_in, evaluate, mask, weight -from shapesmith.histogram import Histogram -from shapesmith.io.skims import read_skims -from shapesmith.parallel import run_jobs -from shapesmith.model import Analysis, Process, Variable, WeightVariation -from shapesmith.skim import SKIM_COLUMNS - -logger = logging.getLogger(__name__) - -NOMINAL_REGION = "nominal" -NOMINAL_VARIATION = "Nominal" -INCLUSIVE = "inclusive" - - -@dataclass(frozen=True, order=True) -class HistKey: - channel: str - category: str - process: str - region: str - variation: str - variable: str - - @property - def directory(self) -> str: - return f"{self.channel}_{self.category}" - - @property - def object_name(self) -> str: - return f"{self.process}#{self.region}#{self.variation}#{self.variable}" - - @property - def path(self) -> str: - return f"{self.directory}/{self.object_name}" - - @classmethod - def parse(cls, path: str) -> "HistKey": - directory, name = path.split("/", 1) - channel, category = directory.split("_", 1) - process, region, variation, variable = name.split("#") - return cls(channel, category, process, region, variation, variable) - - -class HistogramSet: - """HistKey -> Histogram with ROOT persistence (uproot) and a JSON index.""" - - def __init__(self): - self._hists: dict[HistKey, Histogram] = {} - - def add(self, key: HistKey, h: Histogram) -> None: - self._hists[key] = h - - def get(self, key: HistKey) -> Histogram: - return self._hists[key] - - def has(self, key: HistKey) -> bool: - return key in self._hists - - def keys(self, **filters) -> list[HistKey]: - return sorted(k for k in self._hists if all(getattr(k, name) == value for name, value in filters.items())) - - def items(self) -> Iterator[tuple[HistKey, Histogram]]: - return iter(sorted(self._hists.items(), key=lambda kv: kv[0])) - - def update(self, other: "HistogramSet") -> None: - self._hists.update(other._hists) - - def __len__(self) -> int: - return len(self._hists) - - def save(self, path: Path) -> None: - path = Path(path) - path.parent.mkdir(parents=True, exist_ok=True) - with uproot.recreate(path) as f: - for key, h in self.items(): - f[key.path] = h.to_root(key.object_name) - index = [{**asdict(key), "sum": h.sum(), "bins": len(h.values)} for key, h in self.items()] - path.with_suffix(".json").write_text(json.dumps(index, indent=1)) - - @classmethod - def load(cls, path: Path) -> "HistogramSet": - hset = cls() - with uproot.open(path) as f: - for name, classname in f.classnames(recursive=True, cycle=False).items(): - if classname.startswith("TH1"): - hset.add(HistKey.parse(name), Histogram.from_root(f[name])) - return hset - - -@dataclass(frozen=True) -class Booking: - process: Process - regions: tuple[str, ...] - variations: tuple[WeightVariation, ...] - - -@dataclass(frozen=True) -class Target: - category: str - cut: str | None - variable: Variable - - -def resolve_regions(analysis: Analysis, channel_name: str, regions: list[str] | None) -> tuple[str, ...] | None: - """Resolve an explicit region request for one channel. - - ``None`` retains the historical process-dependent booking. ``all`` is - expanded per channel because channels need not define the same regions. - """ - if regions is None: - return None - channel = analysis.channel(channel_name) - names = (NOMINAL_REGION, *(region.name for region in channel.regions)) if "all" in regions else tuple(regions) - resolved = tuple(dict.fromkeys(names)) - for name in resolved: - channel.region(name) # validates with a channel-specific error - return resolved - - -def bookings(analysis: Analysis, channel_name: str, systematics: bool, regions: list[str] | None = None) -> list[Booking]: - requested_regions = resolve_regions(analysis, channel_name, regions) - estimator_regions = tuple(analysis.estimator.regions.values()) if analysis.estimator else () - result = [] - for process in analysis.processes: - if not analysis.samples_for(process.group, channel_name): - continue - process_regions = requested_regions - if process_regions is None: - process_regions = (NOMINAL_REGION,) if process.kind == "signal" else tuple(dict.fromkeys((NOMINAL_REGION, *estimator_regions))) - variations = analysis.weight_variations if systematics and process.kind not in ("data", "embedding") else () - result.append(Booking(process, process_regions, tuple(variations))) - return result - - -def targets(analysis: Analysis, control: bool, variables: list[str] | None) -> list[Target]: - if control: - names = variables or list(analysis.control_variables) - return [Target(INCLUSIVE, None, analysis.control_variables[name]) for name in names] # KeyError for unknown variables - return [Target(c.name, c.cut, c.variable) for c in analysis.categories] - - -def _needed_columns(analysis: Analysis, channel_name: str, booking: Booking, targets_: list[Target], kinds: set[str]) -> set[str]: - """Columns to read for this process: cuts always, weights only for the sample kinds that use them.""" - channel = analysis.channel(channel_name) - process_selection = booking.process.selection_for(channel_name) - columns = set(SKIM_COLUMNS) - for expr in list(channel.baseline.cuts.values()) + list(process_selection.cuts.values()): - columns |= columns_in(expr) - if kinds == {"mc"}: - for expr in channel.baseline.weights.values(): - columns |= columns_in(expr) - for variation in booking.variations: - for expr in variation.replace_weights.values(): - columns |= columns_in(expr) - if kinds <= {"mc", "embedding"}: - for expr in process_selection.weights.values(): - columns |= columns_in(expr) - for region_name in booking.regions: - region = channel.region(region_name) - region_weights = [expr for name, expr in region.add_weights.items() if kinds == {"mc"} or name not in channel.baseline.weights] - for expr in list(region.replace_cuts.values()) + region_weights: - columns |= columns_in(expr) - for target in targets_: - columns |= columns_in(target.variable.expr) | (columns_in(target.cut) if target.cut else set()) - return columns - - -def _varied(weights: dict[str, str], variation: WeightVariation | None) -> dict[str, str] | None: - """Weights with the variation applied, or None if the variation names a weight not present.""" - if variation is None: - return weights - if any(name not in weights for name in variation.replace_weights): - return None - return {**weights, **variation.replace_weights} - - -def fill_booking(config: RunConfig, analysis: Analysis, channel_name: str, booking: Booking, targets_: list[Target]) -> dict[HistKey, Histogram]: - channel = analysis.channel(channel_name) - process = booking.process - samples = analysis.samples_for(process.group, channel_name) - kinds = {s.kind for s in samples} - frame = read_skims(config.skim_dir, channel_name, [s.nick for s in samples], _needed_columns(analysis, channel_name, booking, targets_, kinds)) - lumi = analysis.lumi_pb if kinds == {"mc"} else 1.0 - process_selection = process.selection_for(channel_name) - process_mask = mask(frame, process_selection.cuts) - result: dict[HistKey, Histogram] = {} - for region_name in booking.regions: - region = channel.region(region_name) - cuts = apply_region(channel.baseline, region).cuts - selected = frame[process_mask & mask(frame, cuts)] - weights: dict[str, str] = {} - if kinds == {"mc"}: - weights.update(channel.baseline.weights) - if kinds <= {"mc", "embedding"}: - weights.update(process_selection.weights) - weights.update({name: expr for name, expr in region.add_weights.items() if kinds == {"mc"} or name not in channel.baseline.weights}) - for variation in (None, *booking.variations): - varied = _varied(weights, variation) - if varied is None: - logger.debug(f"{process.name}: variation {variation.name} not applicable, skipped") - continue - w = selected["norm_weight"].to_numpy() * lumi * weight(selected, varied) - for target in targets_: - in_category = mask(selected, {"category": target.cut}) if target.cut else np.ones(len(selected), dtype=bool) - values = evaluate(selected, target.variable.expr).astype(np.float64) - finite = np.isfinite(values) & np.isfinite(w) - if (in_category & ~finite).any(): - logger.warning(f"{process.name}/{region_name}/{target.variable.name}: {(in_category & ~finite).sum()} events with NaN/inf skipped") - keep = in_category & finite - key = HistKey(channel_name, target.category, process.name, region_name, variation.name if variation else NOMINAL_VARIATION, target.variable.name) - result[key] = Histogram.fill(target.variable.edges, values[keep], w[keep]) - return result - - -def _fill_job(args) -> dict[HistKey, Histogram]: - return fill_booking(*args) - - -def run_hist(config: RunConfig, analysis: Analysis, channels: list[str] | None, control: bool, variables: list[str] | None, systematics: bool, processes: list[str] | None, output: Path, force: bool = False, regions: list[str] | None = None) -> HistogramSet: - """Fill every booked histogram of the requested channels and save them to `output` (+ .json index).""" - output = Path(output) - selected_channels = list(channels or config.channels) - targets_ = targets(analysis, control, variables) - jobs = [] - refill_scope: set[tuple[str, str, str, str, str]] = set() - for channel_name in selected_channels: - for booking in bookings(analysis, channel_name, systematics, regions): - if processes and booking.process.key not in processes: - continue - jobs.append((config, analysis, channel_name, booking, targets_)) - refill_scope.update( - (channel_name, booking.process.name, region, target.category, target.variable.name) - for region in booking.regions - for target in targets_ - ) - hset = HistogramSet() - if output.exists(): - if not force: - logger.info(f"{output} exists, loading (use --force to refill)") - return HistogramSet.load(output) - for key, h in HistogramSet.load(output).items(): # keep what this call does not refill - scope = (key.channel, key.process, key.region, key.category, key.variable) - if scope not in refill_scope: - hset.add(key, h) - logger.info(f"{output}: refilling {len(refill_scope)} histogram scopes, keeping {len(hset)} unrelated histograms") - logger.info(f"filling {len(jobs)} process bookings with {config.workers} workers") - for job, result, error in run_jobs(_fill_job, jobs, config.workers): - if error is not None: - raise RuntimeError(f"{job[2]}/{job[3].process.name}: {error}") from error - for key, h in result.items(): - hset.add(key, h) - hset.save(output) - return hset diff --git a/src/shapesmith/io/__init__.py b/src/shapesmith/io/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/src/shapesmith/io/discovery.py b/src/shapesmith/io/discovery.py deleted file mode 100644 index 8a956a5..0000000 --- a/src/shapesmith/io/discovery.py +++ /dev/null @@ -1,61 +0,0 @@ -"""Locate CROWN ntuples and their friend trees (Spec §7 step 1).""" -from __future__ import annotations - -import os -from dataclasses import dataclass - -from shapesmith.config import NtupleConfig - - -@dataclass(frozen=True) -class NtupleFile: - path: str # full path or xrootd URL of the main ntuple - friends: tuple[str, ...] # full paths/URLs of the friend files, same order as the config - nick: str - channel: str - basename: str - - -def join_url(server: str, path: str) -> str: - """'root://server' + '/store/...' -> 'root://server//store/...'; local paths pass through.""" - if not server: - return path - return f"{server.rstrip('/')}/{path}" - - -def list_root_files(server: str, directory: str) -> list[str]: - """Sorted basenames of the .root files in `directory` (xrootd when `server` is set).""" - if server: - from XRootD import client # optional dependency, only needed for remote access - from XRootD.client.flags import DirListFlags - - status, listing = client.FileSystem(server).dirlist(directory, DirListFlags.STAT) - if not status.ok or listing is None: - raise FileNotFoundError(f"cannot list {join_url(server, directory)}: {status.message}") - names = [entry.name for entry in listing] - else: - if not os.path.isdir(directory): - raise FileNotFoundError(f"cannot list {directory}") - names = os.listdir(directory) - return sorted(name for name in names if name.endswith(".root")) - - -def discover(ntuples: NtupleConfig, era: str, nick: str, channel: str, kind: str) -> list[NtupleFile]: - """Main files of one sample and channel with their friend files; missing friends are an error.""" - main_dir = f"{ntuples.base.rstrip('/')}/{era}/{nick}/{channel}" - basenames = list_root_files(ntuples.server, main_dir) - friend_dirs = [f"{f.base.rstrip('/')}/{era}/{nick}/{channel}" for f in ntuples.friends if kind in f.applies_to] - friend_listings = [set(list_root_files(ntuples.server, directory)) for directory in friend_dirs] - missing = [f"{directory}/{name}" for directory, names in zip(friend_dirs, friend_listings) for name in basenames if name not in names] - if missing: - raise FileNotFoundError("missing friend files:\n" + "\n".join(missing)) - return [ - NtupleFile( - path=join_url(ntuples.server, f"{main_dir}/{name}"), - friends=tuple(join_url(ntuples.server, f"{directory}/{name}") for directory in friend_dirs), - nick=nick, - channel=channel, - basename=name, - ) - for name in basenames - ] diff --git a/src/shapesmith/io/ntuples.py b/src/shapesmith/io/ntuples.py deleted file mode 100644 index aaff31a..0000000 --- a/src/shapesmith/io/ntuples.py +++ /dev/null @@ -1,80 +0,0 @@ -"""Read columns from a CROWN ntuple and its friend trees with uproot (Spec §7 steps 2-3).""" -from __future__ import annotations - -import json - -import pandas as pd -import uproot - -from shapesmith.io.discovery import NtupleFile - -TREE = "ntuple" - - -class MissingColumnsError(ValueError): - def __init__(self, missing: list[str], path: str): - self.missing = sorted(missing) - self.path = path - super().__init__(f"{path}: missing columns: {', '.join(self.missing)}") - - def __reduce__(self): # picklable with its two arguments, so that worker processes can hand it back - return (MissingColumnsError, (self.missing, self.path)) - - -def available_columns(ntuple: NtupleFile) -> dict[str, str]: - """column name -> file that provides it (the main file wins over friends).""" - result: dict[str, str] = {} - for path in (*ntuple.friends, ntuple.path): # main file last so that it overrides - with uproot.open(path) as handle: - for name in handle[TREE].keys(): - result[name] = path - return result - - -def read_ntuple(ntuple: NtupleFile, columns: set[str], optional: set[str] = frozenset()) -> tuple[pd.DataFrame, dict]: - """The requested columns of one ntuple (+ friends) and the CROWN metadata of the main file, opening every file once. - - Remote opens cost seconds, so columns are listed and read from the same handle. Columns in `optional` may be - absent (they are left out); any other missing column raises MissingColumnsError. - """ - paths = (*ntuple.friends, ntuple.path) # main file last so that it overrides friends - handles = [uproot.open(path) for path in paths] - try: - sources: dict[str, int] = {} - for index, handle in enumerate(handles): - for name in handle[TREE].keys(): - sources[name] = index - missing = set(columns) - set(sources) - if missing - set(optional): - raise MissingColumnsError(sorted(missing - set(optional)), ntuple.path) - by_handle: dict[int, list[str]] = {} - for column in sorted(set(columns) & set(sources)): - by_handle.setdefault(sources[column], []).append(column) - frames = [pd.DataFrame(handles[index][TREE].arrays(names, library="np")) for index, names in by_handle.items()] - main = handles[-1] - metadata = json.loads(str(main["metadata"])) if "metadata" in main else {} - finally: - for handle in handles: - handle.close() - lengths = {len(frame) for frame in frames} - if len(lengths) > 1: - raise ValueError(f"{ntuple.path}: friend trees have different numbers of entries {sorted(lengths)}") - return (pd.concat(frames, axis=1) if frames else pd.DataFrame()), metadata - - -def read_columns(ntuple: NtupleFile, columns: set[str]) -> pd.DataFrame: - """The requested columns of one ntuple (+ friends) as a DataFrame; dtypes as stored.""" - return read_ntuple(ntuple, columns)[0] - - -def read_metadata(path: str) -> dict: - """The CROWN `metadata` JSON block (empty dict when absent).""" - with uproot.open(path) as handle: - if "metadata" not in handle: - return {} - return json.loads(str(handle["metadata"])) - - -def num_entries(path: str) -> int: - with uproot.open(path) as handle: - return handle[TREE].num_entries diff --git a/src/shapesmith/fit.py b/src/shapesmith/limits.py similarity index 71% rename from src/shapesmith/fit.py rename to src/shapesmith/limits.py index 5dd4b4a..8f9d93b 100644 --- a/src/shapesmith/fit.py +++ b/src/shapesmith/limits.py @@ -1,9 +1,7 @@ -"""Run combine inside the CMSSW environment and collect expected limits, significance and best-fit r (Spec §10.3).""" +"""The limit chain: combine on the datacards (expected limits, significance, best-fit r) and its summary.""" from __future__ import annotations import json -import logging -import subprocess from pathlib import Path import matplotlib @@ -12,37 +10,21 @@ import matplotlib.pyplot as plt # noqa: E402 import uproot # noqa: E402 +from shapesmith import cmssw # noqa: E402 from shapesmith.config import RunConfig # noqa: E402 -logger = logging.getLogger(__name__) - QUANTILES = {0.025: "exp_m2", 0.16: "exp_m1", 0.5: "exp_median", 0.84: "exp_p1", 0.975: "exp_p2", -1.0: "observed"} COMMON = "-m {mass} --setParameterRanges r=-40,40 -t -1" -def combine_script(cmssw_dir: str, scram_arch: str, card_dir: Path, mass: str = "125") -> str: - return "\n".join( - [ - "set -e", - f"export SCRAM_ARCH={scram_arch}", - "source /cvmfs/cms.cern.ch/cmsset_default.sh", - f"cd {cmssw_dir}/src", - "eval $(scramv1 runtime -sh)", - f"cd {Path(card_dir).resolve()}", - f"text2workspace.py combined.txt -o workspace.root -m {mass}", - f"combine -M AsymptoticLimits -d workspace.root {COMMON.format(mass=mass)} --expectSignal 0 -n .Limit", - f"combine -M Significance -d workspace.root {COMMON.format(mass=mass)} --expectSignal 1 -n .Significance", - f"combine -M MultiDimFit --algo singles -d workspace.root {COMMON.format(mass=mass)} --expectSignal 1 -n .Fit", - ] - ) - - -def run_combine(config: RunConfig, card_dir: Path) -> None: - if config.combine is None: - raise ValueError("RunConfig.combine (cmssw_dir) is required to run combine") - script = combine_script(config.combine.cmssw_dir, config.combine.scram_arch, card_dir) - logger.info(f"running combine in {card_dir}") - subprocess.run(["bash", "-c", script], check=True) +def limit_commands(mass: str = "125") -> list[str]: + common = COMMON.format(mass=mass) + return [ + f"text2workspace.py combined.txt -o workspace.root -m {mass}", + f"combine -M AsymptoticLimits -d workspace.root {common} --expectSignal 0 -n .Limit", + f"combine -M Significance -d workspace.root {common} --expectSignal 1 -n .Significance", + f"combine -M MultiDimFit --algo singles -d workspace.root {common} --expectSignal 1 -n .Fit", + ] def _read(path: Path) -> list[tuple[float, float, float | None]]: @@ -107,13 +89,13 @@ def write_summary(results: dict[str, dict], output_dir: Path) -> tuple[Path, Pat return json_path, md_path, pdf_path -def run_fit(config: RunConfig, datacard_dir: Path, final_states: list[str], skip_combine: bool = False) -> dict: +def run_limits(config: RunConfig, datacard_dir: Path, final_states: list[str], skip_combine: bool = False) -> dict: datacard_dir = Path(datacard_dir) results = {} for name in final_states: card_dir = datacard_dir / name if not skip_combine: - run_combine(config, card_dir) + cmssw.run(limit_commands(), config.combine, card_dir) results[name] = collect(card_dir) write_summary(results, datacard_dir) return results diff --git a/src/shapesmith/measurements/__init__.py b/src/shapesmith/measurements/__init__.py new file mode 100644 index 0000000..07851b5 --- /dev/null +++ b/src/shapesmith/measurements/__init__.py @@ -0,0 +1,51 @@ +"""Measurements: plugins on the same event query and histogram tools as the analysis stages. + +An analysis sets `Analysis.measurement` to an object with a `name` and a `run(context)` method; +`shapesmith measure` runs it with its output directory `///`. The measurement reads the skims +itself (`context.events`), possibly in several passes, and writes its results there. +""" +from __future__ import annotations + +from dataclasses import dataclass +from pathlib import Path +from typing import Iterable, Protocol + +from shapesmith.config import RunConfig +from shapesmith.events import Events, Query, load +from shapesmith.model import Analysis +from shapesmith.provenance import record, write_versions + + +@dataclass(frozen=True) +class MeasureContext: + config: RunConfig + analysis: Analysis + channels: list[str] + output: Path + suggest_binning: bool = False # print proposed bin edges instead of measuring + merge: bool = False # combine the results of earlier runs (e.g. one per switch combination) into one payload + + def events(self, query: Query, columns: Iterable[str] = ()) -> Events: + return load(self.config, self.analysis, query, columns) + + def provenance(self, **extra) -> dict: + """What a payload records about its origin: the versions record with the analysis name, the era and `extra`.""" + result = record(self.config) + result["analysis"]["name"] = self.analysis.name + return {**result, "era": self.analysis.era, **extra} + + +class Measurement(Protocol): + name: str + + def run(self, context: MeasureContext) -> None: ... + + +def run_measure(config: RunConfig, analysis: Analysis, channels: list[str], suggest_binning: bool = False, merge: bool = False) -> Path: + if analysis.measurement is None: + raise ValueError(f"analysis {analysis.name} defines no measurement") + output = Path(config.output_dir) / analysis.measurement.name / analysis.era + if not suggest_binning: + write_versions(config, analysis.name, output) + analysis.measurement.run(MeasureContext(config, analysis, channels, output, suggest_binning, merge)) + return output diff --git a/src/shapesmith/measurements/fake_factors/__init__.py b/src/shapesmith/measurements/fake_factors/__init__.py new file mode 100644 index 0000000..e9025a3 --- /dev/null +++ b/src/shapesmith/measurements/fake_factors/__init__.py @@ -0,0 +1,13 @@ +"""The fake-factor measurement (the SM method of TauFakeFactors): model in model.py, the chain in measure.py.""" +from shapesmith.measurements.fake_factors.model import ( # noqa: F401 + Binned, + DataScale, + DrSr, + Equipopulated, + FakeFactorMeasurement, + Fit, + Fractions, + Leg, + ProcessFF, + Split, +) diff --git a/src/shapesmith/measurements/fake_factors/binning.py b/src/shapesmith/measurements/fake_factors/binning.py new file mode 100644 index 0000000..e3be7c2 --- /dev/null +++ b/src/shapesmith/measurements/fake_factors/binning.py @@ -0,0 +1,60 @@ +"""`measure --suggest-binning`: equipopulated bin edges from data, as TauFakeFactors' adjust_binning.py +(FF_Updated _equipopulated_binned_variable and get_binning). + +Per category: the |weight|-weighted quantiles (inverted CDF) of the variable in [low, high], the outer edges set to +low and high, rounded, the `add_left` edges put in front. The data are those of the quantity's region: the SR-like +region of the process for its fake factors and non-closures (the data/MC scale region for an MC fake factor), the +orthogonal SR-like region for the DR->SR correction and its non-closures, the fraction region for the fractions. +""" +from __future__ import annotations + +from typing import TYPE_CHECKING, Iterator + +import numpy as np + +if TYPE_CHECKING: + from shapesmith.measurements.fake_factors.measure import EventSource + from shapesmith.measurements.fake_factors.model import Binned, Equipopulated, FakeFactorMeasurement, Split + + +def equipopulated(x: np.ndarray, weights: np.ndarray, options: Equipopulated, n_bins: int) -> list[float]: + inside = (x >= options.low) & (x <= options.high) + edges = np.quantile(x[inside], np.linspace(0, 1, n_bins + 1), weights=np.abs(weights[inside]), method="inverted_cdf") if inside.any() else np.array([options.low, options.high]) + edges[0], edges[-1] = options.low, options.high + return [*options.add_left, *(round(edge, options.rounding) for edge in np.round(edges, options.rounding).tolist())] + + +def quantities(measurement: FakeFactorMeasurement, channel: str) -> Iterator[tuple[str, Binned, Split, str]]: + """(name, quantity, split, region) of every quantity with equipopulated options.""" + for leg in measurement.legs[channel]: + for process in leg.processes: + name = process.name + leg.suffix + region = process.scale.sr_like if process.scale is not None else process.sr_like + yield f"{name}_fake_factors", process.fake_factors, process.split, region + yield from ((f"{name}_non_closure_{b.variable}", b, process.split, region) for b in process.non_closures) + if process.dr_sr is not None: + dr_sr = process.dr_sr + yield f"{name}_DR_SR", dr_sr.correction, process.split, dr_sr.sr_like + yield from ((f"{name}_for_DRtoSR_non_closure_{b.variable}", b, process.split, dr_sr.sr_like) for b in dr_sr.non_closures) + yield f"process_fractions{leg.suffix}", leg.fractions.binned, leg.fractions.split, leg.fractions.region + + +def suggestions(measurement: FakeFactorMeasurement, source: EventSource) -> list[tuple[str, int, list[float], list[float]]]: + """(quantity, category, suggested edges, current edges).""" + result = [] + for name, binned, split, region in quantities(measurement, source.channel): + if binned.equipopulated is None: + continue + events = source.events(source.data, region) + split_values = events.frame[split.variable].to_numpy() + for category in split.categories: + selected = (split_values >= split.edges[category]) & (split_values < split.edges[category + 1]) + x = events.frame[binned.variable].to_numpy()[selected] + suggested = equipopulated(x, events.weights[selected], binned.equipopulated, binned.equipopulated.n_bins[category]) + result.append((name, category, suggested, list(binned.edges[category]))) + return result + + +def print_suggestions(measurement: FakeFactorMeasurement, source: EventSource) -> None: + for name, category, suggested, current in suggestions(measurement, source): + print(f"{source.channel} {name} category {category}: suggested {suggested}, current {current}") diff --git a/src/shapesmith/measurements/fake_factors/fit.py b/src/shapesmith/measurements/fake_factors/fit.py new file mode 100644 index 0000000..c01ddd9 --- /dev/null +++ b/src/shapesmith/measurements/fake_factors/fit.py @@ -0,0 +1,139 @@ +"""From a measured ratio to the stored step functions with their uncertainty variations (TauFakeFactors' +smooth_function, FF_Updated helper/ff_functions.py sha256 31f89add...). + +The variations are StatShift (the statistical band), SystMCShift (the ratio with every subtracted process shifted by ++-1 sigma, fitted the same way: both sides smoothed at their own centres of mass, user decision U13; TauFakeFactors +smooths the Down side of the corrections at the bin centres, because a deepcopy drops the centres of mass), +SystBandHigh/Low (the curve with bandwidth x1.5/x0.5, Down its mirror around the +nominal) and SystBandAsym (Up the x1.5 and Down the x0.5 curve). A smoothed ratio that is compatible with 1 +(statistical_check) is reset to 1; smoothed curves are then sparsified. The spec's "Algorithm reference" states +every rule. +""" +from __future__ import annotations + +import math +from dataclasses import dataclass, replace + +import numpy as np + +from shapesmith.measurements import smoothing +from shapesmith.measurements.fake_factors.hist import Hist +from shapesmith.measurements.fake_factors.model import Fit + +STAT, MC, BAND_HIGH, BAND_LOW, BAND_ASYM = "StatShift", "SystMCShift", "SystBandHigh", "SystBandLow", "SystBandAsym" +BANDWIDTH_FACTORS = {BAND_HIGH: 1.5, BAND_LOW: 0.5} +RESET_P_VALUE = 0.05 # a smoothed ratio with a larger p-value of the compatibility with 1 is set to 1 +SPARSIFY_THRESHOLD = 0.01 # a stored bin edge per 1 % of the y range of the variation + + +@dataclass(frozen=True) +class Step: + """values[i] on [edges[i], edges[i + 1]).""" + + edges: np.ndarray + values: np.ndarray + + +@dataclass(frozen=True) +class Fitted: + """The fitted ratio: the curve at full resolution (plots) and the stored steps (payload), with variations.""" + + curve: Step + variations: dict[str, np.ndarray] # on curve.edges + stored: Step + stored_variations: dict[str, Step] + p_value: float | None = None # of the compatibility with 1 (smoothed fits only) + reset: bool = False # compatible with 1, set to 1 + + +def fit(ratio: Hist, fit: Fit, mc_shifted: tuple[Hist, Hist] | None = None, stat_sigma: float = 1.0) -> Fitted: + """The ratio as binwise or smoothed steps with all variations; `mc_shifted` are the ratios with the subtracted + processes shifted up and down (None: nothing is subtracted, SystMCShift is the nominal).""" + measured = curve(ratio, fit, stat_sigma) + nominal = measured.nominal + variations = {f"{STAT}Up": measured.up, f"{STAT}Down": measured.down} + shifted = [curve(h, fit, stat_sigma).nominal for h in mc_shifted] if mc_shifted else [nominal, nominal] + variations.update({f"{MC}Up": shifted[0], f"{MC}Down": shifted[1]}) + if fit.kind == "binwise": + variations.update({f"{band}{direction}": nominal for band in (BAND_HIGH, BAND_LOW, BAND_ASYM) for direction in ("Up", "Down")}) + steps = {name: Step(measured.edges, values) for name, values in variations.items()} + return Fitted(Step(measured.edges, nominal), variations, Step(measured.edges, nominal), steps) + bands = {band: curve(ratio, replace(fit, bandwidth=fit.bandwidth * factor), stat_sigma).nominal for band, factor in BANDWIDTH_FACTORS.items()} + for band, values in bands.items(): + variations[f"{band}Up"] = (values - nominal) + nominal + variations[f"{band}Down"] = (nominal - values) + nominal + variations.update({f"{BAND_ASYM}Up": bands[BAND_HIGH], f"{BAND_ASYM}Down": bands[BAND_LOW]}) + p_value = compatibility_with_one(ratio) + reset = p_value > RESET_P_VALUE + if reset: + nominal, variations = _reset(ratio, nominal, variations, mc_shifted, stat_sigma) + stored = Step(*sparsify(measured.edges, nominal)) + stored_variations = {name: Step(*sparsify(measured.edges, values)) for name, values in variations.items()} + return Fitted(Step(measured.edges, nominal), variations, stored, stored_variations, p_value, reset) + + +def curve(ratio: Hist, fit: Fit, stat_sigma: float) -> smoothing.Curve: + """The full-resolution curve of a ratio with its statistical band, at the centres of mass.""" + if fit.kind == "binwise": + return smoothing.binwise(ratio.values, ratio.errors, ratio.edges, stat_sigma) + return smoothing.smooth(ratio.means, ratio.values, ratio.errors, ratio.edges, fit.bandwidth, fit.binwise_left, fit.binwise_right, stat_sigma) + + +def compatibility_with_one(ratio: Hist) -> float: + """TauFakeFactors' statistical_check: the pulls (y - 1) / sigma with the data-only (suppressed) errors; the + smallest two-sided p-value of the pull sum over all windows of adjacent bins, corrected for the number of bins.""" + errors = ratio.suppressed_errors.clip(min=1e-6) + pulls = (ratio.values - 1.0) / errors + n = len(pulls) + p_min = 1.0 + for width in range(1, n + 1): + for start in range(n - width + 1): + z = abs(np.sum(pulls[start : start + width]) / np.sqrt(width)) + p_min = min(p_min, math.erfc(z / math.sqrt(2.0))) # 2 * (1 - Phi(z)) + return 1 - (1 - p_min) ** n + + +def _reset(ratio: Hist, nominal: np.ndarray, variations: dict[str, np.ndarray], mc_shifted: tuple[Hist, Hist] | None, stat_sigma: float) -> tuple[np.ndarray, dict[str, np.ndarray]]: + """A ratio compatible with 1: nominal and band variations 1, the statistical band that of the inverse-variance + mean, SystMCShift the constant fitted to each shifted ratio (1 without shifted ratios).""" + errors = ratio.suppressed_errors.clip(min=1e-6) + inclusive = 1.0 / np.sqrt(np.sum(1.0 / errors**2)) + ones = np.ones_like(nominal) + result = {} + for name, values in variations.items(): + if name.startswith((BAND_ASYM, BAND_HIGH, BAND_LOW)): + result[name] = ones + elif name.startswith(STAT): + result[name] = np.full_like(values, 1.0 + (stat_sigma if name.endswith("Up") else -stat_sigma) * inclusive) + else: + result[name] = values + constants = [constant_fit(h) for h in mc_shifted] if mc_shifted else [1.0, 1.0] + result[f"{MC}Up"], result[f"{MC}Down"] = (np.full_like(nominal, constant) for constant in constants) + return ones, result + + +def constant_fit(hist: Hist) -> float: + """ROOT's pol0 chi2 fit of the bin values: the inverse-variance weighted mean of the bins with an error; 1 if + there is none (the fit fails).""" + errors = hist.errors + used = errors > 0 + if not used.any(): + return 1.0 + weights = 1.0 / errors[used] ** 2 + return float(np.sum(weights * hist.values[used]) / np.sum(weights)) + + +def sparsify(edges: np.ndarray, values: np.ndarray, threshold: float = SPARSIFY_THRESHOLD) -> tuple[np.ndarray, np.ndarray]: + """TauFakeFactors' sparsify: keep a bin edge wherever the cumulative |change| of the values passes another + `threshold` times their range, and the first and the last bin; constant values become one bin.""" + cumulative = np.concatenate(([0], np.cumsum(np.abs(np.diff(values))))) + step = (values.max() - values.min()) * threshold + if step == 0: + return np.array([edges[0], edges[-1]]), np.array([values[0]]) + kept = np.unique(np.searchsorted(cumulative, np.arange(0, cumulative[-1], step))) + kept = kept[kept < len(values)] + if 0 not in kept: + kept = np.insert(kept, 0, 0) + if len(values) - 1 not in kept: + kept = np.append(kept, len(values) - 1) + return np.append(edges[kept], edges[-1]), values[kept] diff --git a/src/shapesmith/measurements/fake_factors/hist.py b/src/shapesmith/measurements/fake_factors/hist.py new file mode 100644 index 0000000..7e42e7a --- /dev/null +++ b/src/shapesmith/measurements/fake_factors/hist.py @@ -0,0 +1,130 @@ +"""The histogram of the fake-factor measurement: TauFakeFactors' patched ROOT TH1D, without ROOT. + +Besides ROOT's bin contents and Sumw2 (`values`, `variances`) a histogram carries per bin the weighted event count +and mean of x (`counts`, `means`: the centre of mass, the x position of the bin in the fits) and TauFakeFactors' +"base" values with their errors, of which `suppressed_errors` keep only the errors of the minuend of every +subtraction (the data errors, the input of the compatibility check). Every operation follows ROOT and FF_Updated +helper/hooks_and_patches.py (sha256 79acd651..., the state of the SM 2018 payload): `add` and `divide` are its +patched TH1.Add/Divide, `shifted` is its AddError, `scaled` is ROOT's unpatched TH1.Scale (contents and Sumw2 +only). tests/reference/make_fake_factor_reference.py checks the port against ROOT with the patches. +""" +from __future__ import annotations + +from dataclasses import dataclass, replace + +import numpy as np + + +@dataclass(frozen=True) +class Hist: + edges: np.ndarray + values: np.ndarray + variances: np.ndarray + counts: np.ndarray + means: np.ndarray + base_values: np.ndarray + base_errors: np.ndarray + suppressed_errors: np.ndarray + + @classmethod + def fill(cls, edges, x: np.ndarray, weights: np.ndarray) -> Hist: + """ROOT's TH1D::Fill per event (bins [low, high), in event order) and TauFakeFactors' centre of mass + (np.histogram, whose last bin includes its upper edge).""" + edges = np.asarray(edges, dtype=np.float64) + x, weights = np.asarray(x), np.asarray(weights, dtype=np.float64) + n = len(edges) - 1 + index = np.searchsorted(edges, x, side="right") - 1 + inside = (index >= 0) & (index < n) + values = np.bincount(index[inside], weights=weights[inside], minlength=n) + variances = np.bincount(index[inside], weights=weights[inside] ** 2, minlength=n) + counts, _ = np.histogram(x, bins=edges, weights=weights) + sums, _ = np.histogram(x, bins=edges, weights=x * weights) + centres = 0.5 * (edges[:-1] + edges[1:]) + means = np.array([s / c if c != 0 else centre for s, c, centre in zip(sums, counts, centres)]) + errors = np.sqrt(variances) + return cls(edges, values, variances, counts, means, values.copy(), errors, errors.copy()) + + @property + def errors(self) -> np.ndarray: + return np.sqrt(self.variances) + + @property + def centres(self) -> np.ndarray: + return 0.5 * (self.edges[:-1] + self.edges[1:]) + + def add(self, other: Hist, factor: float = 1.0) -> Hist: + """self + factor * other; a subtraction (factor < 0) keeps the suppressed errors of self.""" + counts, means = self._centre_of_mass(other, factor, self.counts + factor * other.counts) + base_errors = np.sqrt(self.base_errors**2 + (factor * other.base_errors) ** 2) + suppressed = self.suppressed_errors.copy() if factor < 0 else np.sqrt(self.suppressed_errors**2 + (factor * other.suppressed_errors) ** 2) + return Hist( + self.edges, self.values + factor * other.values, self.variances + factor * factor * other.variances, + counts, means, self.base_values + factor * other.base_values, base_errors, suppressed, + ) + + def divide(self, other: Hist) -> Hist: + """self / other bin by bin, 0 where other is 0 (ROOT's TH1::Divide); the centre of mass of both combined.""" + c0, c1 = self.values, other.values + nonzero = c1 != 0 + values = np.divide(c0, c1, out=np.zeros_like(c0), where=nonzero) + c1sq = c1 * c1 + variances = np.divide(self.variances * c1sq + other.variances * c0 * c0, c1sq * c1sq, out=np.zeros_like(c0), where=nonzero) + ratio_counts = np.divide(self.counts, other.counts, out=np.zeros_like(self.counts), where=other.counts != 0) + counts, means = self._centre_of_mass(other, 1.0, ratio_counts) + base_values = np.divide(self.base_values, other.base_values, out=np.zeros_like(self.base_values), where=other.base_values != 0) + return Hist( + self.edges, values, variances, counts, means, base_values, + _relative_sum(base_values, (self.base_errors, self.base_values), (other.base_errors, other.base_values)), + _relative_sum(base_values, (self.suppressed_errors, self.base_values), (other.suppressed_errors, other.base_values)), + ) + + def scaled(self, factor: float) -> Hist: + """ROOT's TH1::Scale: contents and Sumw2 only (TauFakeFactors does not patch it).""" + return replace(self, values=self.values * factor, variances=self.variances * (factor * factor)) + + def shifted(self, sigmas: float) -> Hist: + """Every bin moved by `sigmas` times its error and clipped at 0 (TauFakeFactors' AddError); the error of a + clipped bin is 0, and so is its count, whose mean moves to the bin centre. The base values are not clipped.""" + errors = self.errors + values = np.maximum(0.0, self.values + sigmas * errors) + variances = np.where(values > 0, errors, 0.0) ** 2 + moved = sigmas * (abs(sigmas) * errors) + counts = self.counts + moved + means = np.divide(self.counts * self.means + moved * self.means, counts, out=self.centres.copy(), where=counts != 0) + means = np.where(counts < 0, self.centres, _clip_means(means, self.edges)) + return Hist( + self.edges, values, variances, np.maximum(0.0, counts), means, + self.base_values + sigmas * self.base_errors, self.base_errors.copy(), self.suppressed_errors.copy(), + ) + + def with_negative_bins_zeroed(self) -> Hist: + """Negative contents set to 0 (errors, centre of mass and base values unchanged).""" + return replace(self, values=np.maximum(self.values, 0.0)) + + def with_empty_bins_set_to_one(self) -> Hist: + """Bins <= 0 become 1 +- 1 with the bin centre as their mean and count 1 (TauFakeFactors' + calculate_non_closure_correction); the base values are unchanged.""" + empty = self.values <= 0.0 + return replace( + self, values=np.where(empty, 1.0, self.values), variances=np.where(empty, 1.0, self.variances), + counts=np.where(empty, 1.0, self.counts), means=np.where(empty, self.centres, self.means), + ) + + def _centre_of_mass(self, other: Hist, factor: float, counts: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + """TauFakeFactors' calc_center_of_mass: the means weighted with the counts of both, inside the bin.""" + total = self.counts + factor * other.counts + weighted = self.counts * self.means + factor * other.counts * other.means + means = np.divide(weighted, total, out=self.centres.copy(), where=total != 0) + return counts, np.where(total != 0, _clip_means(means, self.edges), means) + + +def _clip_means(means: np.ndarray, edges: np.ndarray) -> np.ndarray: + """Means kept strictly inside their bins.""" + low, high = edges[:-1], edges[1:] + return np.maximum(np.minimum(means, np.nextafter(high, low)), np.nextafter(low, high)) + + +def _relative_sum(values: np.ndarray, *parts: tuple[np.ndarray, np.ndarray]) -> np.ndarray: + """|values| times the quadratic sum of the relative errors error / value of the parts (0 where value is 0).""" + squares = sum(np.divide(error, value, out=np.zeros_like(error), where=value != 0) ** 2 for error, value in parts) + return np.abs(values) * np.sqrt(squares) diff --git a/src/shapesmith/measurements/fake_factors/measure.py b/src/shapesmith/measurements/fake_factors/measure.py new file mode 100644 index 0000000..9cea16e --- /dev/null +++ b/src/shapesmith/measurements/fake_factors/measure.py @@ -0,0 +1,277 @@ +"""The fake-factor measurement of one channel, in TauFakeFactors' order (FF_Updated ff_calculation.py and +ff_corrections.py): + +1. per leg the QCD and ttbar fake factors and the process fractions -> fake_factors_.json.gz; +2. per QCD leg the orthogonal fake factors and their non-closures in the DR->SR regions, then the DR->SR correction + (the intermediate payloads -> FF_for_DRtoSR_.json.gz); +3. per process the non-closures in table order, each measured with the fake factors, the DR->SR correction (QCD) and + the earlier non-closures applied per event -> FF_corrections_.json.gz. + +Fake factors and corrections are applied per event with correctionlib, from the payload models of this run, as the +CROWN friend applies them. A "yield" is the histogram of the target process minus the subtracted processes, nominal +and with every subtracted histogram shifted by +-1 sigma (SystMCShift). +""" +from __future__ import annotations + +import json +import logging +from dataclasses import dataclass, replace +from typing import TYPE_CHECKING, Callable + +import correctionlib.schemav2 as cs +import numpy as np + +from shapesmith import payloads +from shapesmith.events import Events, Query +from shapesmith.measurements.fake_factors import payload, plots +from shapesmith.measurements.fake_factors.binning import print_suggestions +from shapesmith.measurements.fake_factors.fit import Fitted, fit +from shapesmith.measurements.fake_factors.hist import Hist + +if TYPE_CHECKING: + from shapesmith.measurements import MeasureContext + from shapesmith.measurements.fake_factors.model import Binned, FakeFactorMeasurement, Fractions, Leg, ProcessFF, Split + +logger = logging.getLogger(__name__) + +Yields = tuple[Hist, Hist, Hist] # nominal, subtracted shifted up, subtracted shifted down + + +@dataclass(frozen=True) +class Ratio: + """A measured ratio of the nominal yields, the ratios with shifted subtraction (None without subtraction) and the + nominal yields themselves (the closure plots).""" + + nominal: Hist + mc_shifted: tuple[Hist, Hist] | None + numerator: Hist + denominator: Hist + + +class EventSource: + """The selected events of the processes of one channel per region, each read once.""" + + def __init__(self, context: MeasureContext, channel: str, columns: set[str]): + self.channel = channel + self.data = context.analysis.channel(channel).data() + self._context, self._columns, self._events = context, columns, {} + + def events(self, process: str, region: str) -> Events: + if (process, region) not in self._events: + self._events[process, region] = self._context.events(Query(self.channel, process, region), self._columns) + return self._events[process, region] + + def hist(self, process: str, region: str, variable: str, edges, split: Split, category: int, weights: np.ndarray | None = None) -> Hist: + """`variable` in one category of the split, with the event weights or the given per-event weights.""" + events = self.events(process, region) + values = events.frame[split.variable].to_numpy() + selected = (values >= split.edges[category]) & (values < split.edges[category + 1]) + weights = events.weights if weights is None else weights + return Hist.fill(edges, events.frame[variable].to_numpy()[selected], weights[selected]) + + def weighted(self, process: str, region: str, corrections: list[Callable]) -> np.ndarray: + """The event weights times the corrections evaluated per event, multiplied in order.""" + events = self.events(process, region) + weights = events.weights + for correction in corrections: + weights = weights * correction(events.frame) + return weights + + +def evaluator(correction: cs.Correction) -> Callable: + """The nominal of a payload correction per event, with the inputs as float (as TauFakeFactors and the friend).""" + evaluate = payloads.correction_set([correction], {}).to_evaluator()[correction.name].evaluate + columns = [variable.name for variable in correction.inputs[:-1]] # the last input is syst + return lambda frame: evaluate(*(np.asarray(frame[c], dtype=np.float32).astype(np.float64) for c in columns), "nominal") + + +class ChannelMeasurement: + """The measurement of one channel: the payload corrections and a record of every measured ratio and fit.""" + + def __init__(self, measurement: FakeFactorMeasurement, source: EventSource): + self.measurement, self.source = measurement, source + self.fake_factors: list[cs.Correction] = [] + self.intermediate: list[cs.Correction] = [] # the orthogonal fake factors and their non-closures + self.corrections: list[cs.Correction] = [] + self.compounds: list[cs.CompoundCorrection] = [] + self.record: dict[str, list[dict]] = {} + + def measure(self, legs: tuple[Leg, ...]) -> None: + for leg in legs: + self.fake_factors += [self.process_fake_factors(p.name + leg.suffix, p, p.sr_like, p.ar_like, p.subtract) for p in leg.processes] + self.fake_factors.append(self.fractions(f"process_fractions{leg.suffix}", leg.fractions)) + by_name = {c.name: c for c in self.fake_factors} + for leg in legs: + for process in leg.processes: + name = process.name + leg.suffix + applied = [by_name[f"{name}_fake_factors"]] + if process.dr_sr is not None: + applied.append(self.dr_sr(name, process)) + closures = self.non_closures(name, process, process.non_closures, process.sr_like, process.ar_like, process.subtract, applied) + self.corrections += closures + if closures: + self.compounds.append(payload.compound(name, process.non_closures, process.split)) + + def dr_sr(self, name: str, process: ProcessFF) -> cs.Correction: + """The DR->SR correction, measured with the orthogonal fake factors and their non-closures.""" + dr_sr = process.dr_sr + orthogonal = self.process_fake_factors(name, process, dr_sr.sr_like, dr_sr.ar_like, dr_sr.subtract, stage="for_DRtoSR/") + closures = self.non_closures(name, process, dr_sr.non_closures, dr_sr.sr_like, dr_sr.ar_like, dr_sr.subtract, [orthogonal], stage="for_DRtoSR/") + self.intermediate += [orthogonal, *closures] + correction = self.closure(name, "DR_SR", process, dr_sr.correction, dr_sr.sr, dr_sr.ar, dr_sr.subtract, [orthogonal, *closures]) + self.corrections.append(correction) + return correction + + def process_fake_factors(self, name: str, process: ProcessFF, sr_like: str, ar_like: str, subtract: tuple[str, ...], stage: str = "") -> cs.Correction: + """{name}_fake_factors, the target yield in `sr_like` over that in `ar_like`, times the data/MC factor of an + MC fake factor.""" + scale = None if process.scale is None else self.data_scale(process) + if scale is not None: + self.record[f"{stage}{name}_data_scale"] = [{"factor": scale}] + + def measure(category: int) -> Ratio: + ratio = self.ratio(process, process.fake_factors, category, (sr_like, ar_like), subtract, from_data=scale is None) + return ratio if scale is None else replace(ratio, nominal=ratio.nominal.scaled(scale)) + + fitted = self.fitted(f"{stage}{name}_fake_factors", process.fake_factors, process.split, measure) + return payload.fake_factors(name, process.fake_factors, process.split, fitted, self.measurement.version) + + def closure(self, name: str, correction_name: str, process: ProcessFF, binned: Binned, sr_like: str, ar_like: str, subtract: tuple[str, ...], applied: list[cs.Correction], stage: str = "") -> cs.Correction: + """{name}_{correction_name}_correction: the target yield in `sr_like` over that in `ar_like` weighted per event + with the `applied` fake factors and corrections. From data, bins <= 0 are set to 1 +- 1.""" + from_data = process.scale is None + evaluators = [evaluator(correction) for correction in applied] + weights = {p: self.source.weighted(p, ar_like, evaluators) for p in (process.target, *subtract)} + + def measure(category: int) -> Ratio: + ratio = self.ratio(process, binned, category, (sr_like, ar_like), subtract, from_data, weights) + if not from_data: + return ratio + return replace(ratio, nominal=ratio.nominal.with_empty_bins_set_to_one(), mc_shifted=tuple(h.with_empty_bins_set_to_one() for h in ratio.mc_shifted)) + + fitted = self.fitted(f"{stage}{name}_{correction_name}", binned, process.split, measure) + return payload.correction(name, correction_name, binned, process.split, fitted, self.measurement.version) + + def non_closures(self, name: str, process: ProcessFF, non_closures: tuple[Binned, ...], sr_like: str, ar_like: str, subtract: tuple[str, ...], applied: list[cs.Correction], stage: str = "") -> list[cs.Correction]: + """The non-closures in table order, each on top of `applied` and the earlier ones.""" + result = [] + for binned in non_closures: + result.append(self.closure(name, f"non_closure_{binned.variable}", process, binned, sr_like, ar_like, subtract, applied + result, stage)) + return result + + def ratio(self, process: ProcessFF, binned: Binned, category: int, regions: tuple[str, str], subtract: tuple[str, ...], from_data: bool, ar_weights: dict[str, np.ndarray] | None = None) -> Ratio: + """The yield in the SR-like over the (weighted) yield in the AR-like region; from data also the ratios of the + shifted yields.""" + sr_like, ar_like = regions + numerator = self.yields(process.target, subtract, sr_like, binned, process.split, category) + denominator = self.yields(process.target, subtract, ar_like, binned, process.split, category, ar_weights) + nominal, up, down = (n.divide(d) for n, d in zip(numerator, denominator)) + return Ratio(nominal, (up, down) if from_data else None, numerator[0], denominator[0]) + + def yields(self, target: str, subtract: tuple[str, ...], region: str, binned: Binned, split: Split, category: int, weights: dict[str, np.ndarray] | None = None) -> Yields: + edges = binned.edges[category] + + def hist(process: str) -> Hist: + return self.source.hist(process, region, binned.variable, edges, split, category, weights[process] if weights else None) + + nominal = up = down = hist(target) + for process in subtract: + h = hist(process) + nominal, up, down = nominal.add(h, -1.0), up.add(h.shifted(1.0), -1.0), down.add(h.shifted(-1.0), -1.0) + return nominal, up, down + + def data_scale(self, process: ProcessFF) -> float: + """The global data/MC factor of an MC fake factor, from the event totals of its four scale regions (summed and + subtracted in order, as TauFakeFactors' one-bin histograms).""" + scale, data = process.scale, self.source.data + + def total(name: str, region: str) -> float: + weights = self.source.events(name, region).weights + return float(np.cumsum(weights)[-1]) if len(weights) else 0.0 + + def minus(value: float, names, region: str) -> float: + for name in names: + value -= total(name, region) + return value + + def data_minus_mc(region: str, same_sign: str) -> float: + qcd = max(0.0, minus(total(data, same_sign), scale.subtract, same_sign)) + return minus(total(data, region), (p for p in scale.subtract if p != process.target), region) - qcd + + ratio = data_minus_mc(scale.sr_like, scale.sr_like_same_sign) / data_minus_mc(scale.ar_like, scale.ar_like_same_sign) + return ratio / (total(process.target, scale.sr_like) / total(process.target, scale.ar_like)) + + def fractions(self, name: str, fractions: Fractions) -> cs.Correction: + """process_fractions{s}: QCD and ttbar over their sum, nominal and with each of them shifted by +-1 sigma.""" + values: dict[str, dict[str, list[np.ndarray]]] = {} + for category in fractions.split.categories: + hists = self.fraction_hists(fractions, category) + variations = {"nominal": hists} + for varied in hists: + variations[f"frac_{varied}_up"] = {**hists, varied: hists[varied].shifted(1.0)} + variations[f"frac_{varied}_down"] = {**hists, varied: hists[varied].shifted(-1.0)} + for variation, shifted in variations.items(): + for process, h in shifted.items(): + total = h + for other in (p for p in shifted if p != process): + total = total.add(shifted[other]) + values.setdefault(variation, {}).setdefault(process, []).append(h.divide(total).values) + self.record.setdefault(name, []).append({ + "category": category, "edges": list(fractions.binned.edges[category]), + "yields": {p: h.values.tolist() for p, h in hists.items()}, "fractions": {p: values["nominal"][p][-1].tolist() for p in hists}, + }) + return payload.fractions(name, fractions.binned, fractions.split, values, self.measurement.version) + + def fraction_hists(self, fractions: Fractions, category: int) -> dict[str, Hist]: + """QCD = data - all MC with negative bins 0, and ttbar, in one category.""" + def hist(process: str) -> Hist: + return self.source.hist(process, fractions.region, fractions.binned.variable, fractions.binned.edges[category], fractions.split, category) + + qcd = hist(self.source.data) + for process in fractions.subtract: + qcd = qcd.add(hist(process), -1.0) + return {"QCD": qcd.with_negative_bins_zeroed(), "ttbar_J": hist(fractions.ttbar)} + + def fitted(self, name: str, binned: Binned, split: Split, measure: Callable[[int], Ratio]) -> list[Fitted]: + """The fit per category of the ratio `measure(category)`, recorded under `name`.""" + result = [] + for category in split.categories: + ratio = measure(category) + fitted = fit(ratio.nominal, binned.fits[category], ratio.mc_shifted, self.measurement.stat_sigma) + self.record.setdefault(name, []).append(_record(binned.variable, split, category, ratio, fitted)) + result.append(fitted) + return result + + +def run(measurement: FakeFactorMeasurement, context: MeasureContext) -> None: + for channel in context.channels: + source = EventSource(context, channel, measurement.columns(channel)) + if context.suggest_binning: + print_suggestions(measurement, source) + continue + logger.info(f"{channel}: measuring the fake factors") + result = ChannelMeasurement(measurement, source) + result.measure(measurement.legs[channel]) + provenance, output = context.provenance(channel=channel), context.output + payloads.write(payloads.correction_set(result.fake_factors, provenance), output / f"fake_factors_{channel}.json.gz") + payloads.write(payloads.correction_set(result.corrections, provenance, result.compounds), output / f"FF_corrections_{channel}.json.gz") + payloads.write(payloads.correction_set(result.intermediate, provenance), output / f"FF_for_DRtoSR_{channel}.json.gz") + (output / f"measurement_{channel}.json").write_text(json.dumps(result.record, indent=1)) + plots.plot_record(result.record, output / "plots" / channel) + logger.info(f"{channel}: payloads written to {output}") + + +def _points(h: Hist) -> dict: + return {"x": h.means.tolist(), "y": h.values.tolist(), "errors": h.errors.tolist()} + + +def _record(variable: str, split: Split, category: int, ratio: Ratio, fitted: Fitted) -> dict: + """The fit inputs (the ratio and the shifted ratios as points at their centres of mass, the yields) and results.""" + return { + "variable": variable, "category": category, "split": [split.edges[category], split.edges[category + 1]], + "edges": ratio.nominal.edges.tolist(), **_points(ratio.nominal), + "mc_shifted": [_points(h) for h in ratio.mc_shifted] if ratio.mc_shifted else None, + "numerator": ratio.numerator.values.tolist(), "denominator": ratio.denominator.values.tolist(), + "p_value": fitted.p_value, "reset": fitted.reset, + "curve": {"edges": fitted.curve.edges.tolist(), "nominal": fitted.curve.values.tolist(), **{k: v.tolist() for k, v in fitted.variations.items()}}, + } diff --git a/src/shapesmith/measurements/fake_factors/model.py b/src/shapesmith/measurements/fake_factors/model.py new file mode 100644 index 0000000..7fb2964 --- /dev/null +++ b/src/shapesmith/measurements/fake_factors/model.py @@ -0,0 +1,157 @@ +"""The fake-factor measurement as data: the SM method of TauFakeFactors (QCD and ttbar fake factors, their fractions, +the QCD DR->SR correction and the non-closure corrections), built by an analysis from its regions and tables. + +Regions and processes are names in the analysis channel; every quantity is measured per category of one split +variable (n_jets), in bins of its own variable. +""" +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import Literal, Mapping + + +@dataclass(frozen=True) +class Fit: + """How a measured ratio becomes the stored function: `binwise` keeps the bins; `smoothed` is the kernel + smoothing with `bandwidth`, where the first `binwise_left` and the last `binwise_right` bins keep their values.""" + + kind: Literal["binwise", "smoothed"] + bandwidth: float = 0.0 + binwise_left: int = 0 + binwise_right: int = 0 + + +@dataclass(frozen=True) +class Equipopulated: + """The --suggest-binning options of a quantity (TauFakeFactors' equipopulated_binning_options): `n_bins` per + category over [low, high], edges rounded to `rounding` digits, `add_left` edges put in front.""" + + n_bins: tuple[int, ...] + low: float + high: float + rounding: int = 2 + add_left: tuple[float, ...] = () + + +@dataclass(frozen=True) +class Split: + """Category i holds the events with edges[i] <= variable < edges[i + 1]; the payload bins on the same edges.""" + + variable: str + edges: tuple[float, ...] + + @property + def categories(self) -> range: + return range(len(self.edges) - 1) + + +@dataclass(frozen=True) +class Binned: + """A quantity in bins of `variable`: edges[i] and fits[i] in category i of the split (fractions have no fits).""" + + variable: str + edges: tuple[tuple[float, ...], ...] + fits: tuple[Fit, ...] = () + equipopulated: Equipopulated | None = None + + +@dataclass(frozen=True) +class DrSr: + """The QCD DR->SR correction. Orthogonal fake factors (the QCD variable, bins and fits) from data in + `sr_like`/`ar_like`, with their own `non_closures` there; the correction is the SR-like yield over the + orthogonal-FF weighted AR-like yield in `sr`/`ar`. Every step subtracts `subtract`.""" + + sr_like: str + ar_like: str + sr: str + ar: str + subtract: tuple[str, ...] + correction: Binned + non_closures: tuple[Binned, ...] = () + + +@dataclass(frozen=True) +class DataScale: + """The global data/MC factor of an MC fake factor: (data - the MC but the target - QCD)(sr_like) / (...)(ar_like) + over the same ratio of the target, from event totals; QCD is data minus all MC in the same-sign version of each + region. `subtract` lists all MC processes, the target among them, in subtraction order.""" + + sr_like: str + ar_like: str + sr_like_same_sign: str + ar_like_same_sign: str + subtract: tuple[str, ...] + + +@dataclass(frozen=True) +class ProcessFF: + """The fake factors of one process, `target`(sr_like) / `target`(ar_like) with target = the `target` process + minus the `subtract` processes: from data (QCD), or from MC scaled by `scale` (ttbar). The non-closures are + measured in the same regions in table order, each on top of the fake factors, the DR->SR correction and the + earlier non-closures.""" + + name: str # the payload process: QCD or ttbar + target: str + subtract: tuple[str, ...] + sr_like: str + ar_like: str + split: Split + fake_factors: Binned + non_closures: tuple[Binned, ...] = () + dr_sr: DrSr | None = None + scale: DataScale | None = None + + +@dataclass(frozen=True) +class Fractions: + """The QCD and ttbar fractions in `region`: QCD = data - `subtract` (negative bins 0), ttbar = `ttbar`.""" + + region: str + subtract: tuple[str, ...] + ttbar: str + split: Split + binned: Binned + + +@dataclass(frozen=True) +class Leg: + """The fake factors of one hadronic tau: payload names carry `suffix` ("" or "_subleading").""" + + suffix: str + qcd: ProcessFF + ttbar: ProcessFF + fractions: Fractions + + @property + def processes(self) -> tuple[ProcessFF, ProcessFF]: + return (self.qcd, self.ttbar) + + +@dataclass(frozen=True) +class FakeFactorMeasurement: + """`shapesmith measure` for the fake factors of every channel: the payloads fake_factors_.json.gz and + FF_corrections_.json.gz in /fake_factors//.""" + + legs: Mapping[str, tuple[Leg, ...]] # channel -> legs + version: int = 1 # Correction.version of the payloads + stat_sigma: float = 1.0 + name: str = field(default="fake_factors", init=False) + + def columns(self, channel: str) -> set[str]: + """The columns every quantity of the channel is measured and evaluated in.""" + columns = set() + for leg in self.legs[channel]: + columns |= {leg.fractions.split.variable, leg.fractions.binned.variable} + for process in leg.processes: + binned = [process.fake_factors, *process.non_closures] + if process.dr_sr is not None: + binned += [process.dr_sr.correction, *process.dr_sr.non_closures] + columns |= {process.split.variable, *(b.variable for b in binned)} + return columns + + def run(self, context) -> None: + if context.merge: + raise ValueError("fake_factors: --merge is not supported, one run measures every channel") + from shapesmith.measurements.fake_factors.measure import run # plots and correctionlib load only for a run + + run(self, context) diff --git a/src/shapesmith/measurements/fake_factors/payload.py b/src/shapesmith/measurements/fake_factors/payload.py new file mode 100644 index 0000000..22c5f71 --- /dev/null +++ b/src/shapesmith/measurements/fake_factors/payload.py @@ -0,0 +1,77 @@ +"""The correctionlib models of the fake-factor payloads, in the conventions the CROWN fake-factor friend reads. + +| correction | inputs | syst keys (default: the nominal) | +|---|---|---| +| `{p}{s}_fake_factors` | variable, split, syst | `{p}{s}{StatShift,SystMCShift,SystBandHigh,SystBandLow,SystBandAsym}{Up,Down}` | +| `process_fractions{s}` | process, variable, split, syst | `process_fractions{s}{nominal,frac_QCD_up/down,frac_ttbar_J_up/down}`; processes QCD, ttbar | +| `QCD{s}_DR_SR_correction` | variable, split, syst | `QCD{s}_DR_SR_Corr{StatShift,SystMCShift,SystBandAsym}{Up,Down}`, `nominal` | +| `{p}{s}_non_closure__correction` | variable, split, syst | `{p}{s}_non_closure__Corr...` and the global `{p}{s}_non_closure_Corr...`, `nominal` | +| `{p}{s}_compound_correction` | the non-closure variables in order, split, syst | the product of the non-closures of {p}{s} | + +{p} is QCD or ttbar, {s} the leg suffix. Every correction is binned in the split categories first, then in its +variable (flow clamp), below a `syst` category. +""" +from __future__ import annotations + +from typing import Mapping, Sequence + +import correctionlib.schemav2 as cs + +from shapesmith import payloads +from shapesmith.measurements.fake_factors.fit import BAND_ASYM, MC, STAT, Fitted, Step +from shapesmith.measurements.fake_factors.model import Binned, Split + +CORRECTION_VARIATIONS = tuple(f"{name}{direction}" for name in (STAT, MC, BAND_ASYM) for direction in ("Up", "Down")) +FRACTION_PROCESSES = {"QCD": "QCD", "ttbar_J": "ttbar"} # TauFakeFactors process -> payload category + + +def _by_category(split: Split, variable: str, steps: Sequence[Step]) -> cs.Binning: + return payloads.binning(split.variable, split.edges, [payloads.binning(variable, step.edges, step.values) for step in steps]) + + +def _inputs(*names: str) -> list[cs.Variable]: + return [payloads.variable(name) for name in names] + [payloads.variable("syst", "string")] + + +def fake_factors(name: str, binned: Binned, split: Split, fitted: Sequence[Fitted], version: int) -> cs.Correction: + """`name` = {p}{s}; the stored fit per category.""" + keys = fitted[0].stored_variations + content = {name + key: _by_category(split, binned.variable, [f.stored_variations[key] for f in fitted]) for key in keys} + data = payloads.category("syst", content, default=_by_category(split, binned.variable, [f.stored for f in fitted])) + return payloads.correction(f"{name}_fake_factors", _inputs(binned.variable, split.variable), data, version, output=payloads.variable(f"{name}_ff")) + + +def fractions(name: str, binned: Binned, split: Split, values: Mapping[str, Mapping[str, Sequence]], version: int) -> cs.Correction: + """`name` = process_fractions{s}; values[variation][process] holds the fraction values per category, variation + nominal or frac__up/down.""" + + def by_process(variation: str) -> cs.Category: + return payloads.category("process", { + label: _by_category(split, binned.variable, [Step(edges, category) for edges, category in zip(binned.edges, values[variation][process])]) + for process, label in FRACTION_PROCESSES.items() + }) + + data = payloads.category("syst", {name + variation: by_process(variation) for variation in values}, default=by_process("nominal")) + return payloads.correction(name, [payloads.variable("process", "string"), *_inputs(binned.variable, split.variable)], data, version, output=payloads.variable("fraction")) + + +def correction(process: str, correction_name: str, binned: Binned, split: Split, fitted: Sequence[Fitted], version: int) -> cs.Correction: + """`process` = {p}{s}, `correction_name` DR_SR or non_closure_; a non-closure repeats its keys under the + global non_closure name, which the compound correction shifts as one.""" + names = [f"_{correction_name}"] + (["_non_closure"] if correction_name.startswith("non_closure") else []) + nominal = _by_category(split, binned.variable, [f.stored for f in fitted]) + content = { + f"{process}{name}_Corr{variation}": _by_category(split, binned.variable, [f.stored_variations[variation] for f in fitted]) + for variation in CORRECTION_VARIATIONS for name in names + } + content["nominal"] = nominal + full_name = f"{process}_{correction_name}_correction" + data = payloads.category("syst", content, default=nominal) + return payloads.correction(full_name, _inputs(binned.variable, split.variable), data, version, output=payloads.variable(full_name)) + + +def compound(process: str, non_closures: Sequence[Binned], split: Split) -> cs.CompoundCorrection: + """The product of the non-closure corrections of {p}{s} = `process`, in table order.""" + name = f"{process}_compound_correction" + stack = [f"{process}_non_closure_{binned.variable}_correction" for binned in non_closures] + return payloads.compound(name, _inputs(*(b.variable for b in non_closures), split.variable), stack, output=payloads.variable(name)) diff --git a/src/shapesmith/measurements/fake_factors/plots.py b/src/shapesmith/measurements/fake_factors/plots.py new file mode 100644 index 0000000..c9b4720 --- /dev/null +++ b/src/shapesmith/measurements/fake_factors/plots.py @@ -0,0 +1,65 @@ +"""Plots of the fake-factor measurement from its record: every fitted ratio with its curve and bands per category, +the SR-like yield against the (weighted) AR-like yield below it (the closure), and the fractions.""" +from __future__ import annotations + +from pathlib import Path + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt # noqa: E402 +import mplhep as hep # noqa: E402 +import numpy as np # noqa: E402 + +hep.style.use("CMS") + + +def plot_record(record: dict[str, list[dict]], directory: Path) -> list[Path]: + directory.mkdir(parents=True, exist_ok=True) + paths = [] + for name, entries in record.items(): + for entry in entries: + path = directory / f"{name.replace('/', '_')}_category{entry.get('category', 0)}.png" + if "curve" in entry: + paths.append(_plot_fit(name, entry, path)) + elif "fractions" in entry: + paths.append(_plot_fractions(name, entry, path)) + return paths + + +def _plot_fit(name: str, entry: dict, path: Path) -> Path: + figure, (top, bottom) = plt.subplots(2, 1, figsize=(10, 10), sharex=True, gridspec_kw={"height_ratios": (2, 1)}) + curve = entry["curve"] + edges = np.asarray(curve["edges"]) + for band, colour in (("SystBandAsym", "#f89c20"), ("StatShift", "#5790fc")): + top.stairs(np.asarray(curve[f"{band}Up"]), edges, baseline=np.asarray(curve[f"{band}Down"]), fill=True, alpha=0.4, color=colour, label=band) + top.stairs(np.asarray(curve["nominal"]), edges, baseline=None, color="black", label="fit") + top.errorbar(entry["x"], entry["y"], yerr=entry["errors"], fmt="o", color="black", label="measured") + top.set_ylabel("ratio") + top.set_title(f"{name} {entry['split']}" + (" (reset to 1)" if entry["reset"] else ""), fontsize=14) + top.legend(fontsize=12) + bin_edges = entry["edges"] + closure = not name.endswith("_fake_factors") # a correction's denominator is the prediction from the AR-like region + bottom.stairs(entry["numerator"], bin_edges, color="black", label="SR-like") + bottom.stairs(entry["denominator"], bin_edges, color="#e42536", label="AR-like x fake factors" if closure else "AR-like") + bottom.set_ylabel("yield") + bottom.set_xlabel(entry["variable"]) + bottom.legend(fontsize=12) + figure.savefig(path) + plt.close(figure) + return path + + +def _plot_fractions(name: str, entry: dict, path: Path) -> Path: + figure, ax = plt.subplots(figsize=(10, 8)) + bottom = np.zeros(len(entry["edges"]) - 1) + for (process, values), colour in zip(entry["fractions"].items(), ("#5790fc", "#e42536")): + ax.stairs(bottom + np.asarray(values), entry["edges"], baseline=bottom, fill=True, color=colour, label=process) + bottom = bottom + np.asarray(values) + ax.set_ylim(0, 1.3) + ax.set_ylabel("fraction") + ax.set_title(f"{name} category {entry['category']}", fontsize=14) + ax.legend(fontsize=12) + figure.savefig(path) + plt.close(figure) + return path diff --git a/src/shapesmith/measurements/smoothing.py b/src/shapesmith/measurements/smoothing.py new file mode 100644 index 0000000..44a1d8b --- /dev/null +++ b/src/shapesmith/measurements/smoothing.py @@ -0,0 +1,106 @@ +"""Kernel smoothing of binned ratios as TauFakeFactors does it, without ROOT. + +The reference is FF_Updated `helper/ff_functions.py` (`_smooth_function`, `_get_index_and_slices`; sha256 31f89add..., +the state that produced the SM 2018 payload): 100 samples per bin from Normal(y, stat_sigma * sigma) with seed 19, +each smoothed with ROOT's TGraphSmooth::SmoothKern (normal kernel) on the grid linspace(x[start], x[end], 100 n + 1) +without its last point; the nominal is the mean and the statistical band the standard deviation over the samples. +A hybrid keeps the first or last bins binwise (TauFakeFactors' `binwise#[0]+smoothed`: the first bin). +tests/reference/make_smoothing_reference.py checks the port against ROOT and FF_Updated. +""" +from __future__ import annotations + +import math +from dataclasses import dataclass + +import numpy as np + +KERNEL_SCALE = 0.3706506 # ROOT's normal kernel: quartiles at +-0.25 * bandwidth +SAMPLES = 100 +SEED = 19 + + +@dataclass(frozen=True) +class Curve: + """A piecewise constant function: `nominal[i]` on [edges[i], edges[i + 1]), with its statistical band.""" + + edges: np.ndarray + nominal: np.ndarray + up: np.ndarray + down: np.ndarray + + +def smooth_kern(x: np.ndarray, y: np.ndarray, x_out: np.ndarray, bandwidth: float) -> np.ndarray: + """ROOT's TGraphSmooth::SmoothKern with the normal kernel (BDRksmooth), at the sorted points `x_out`. + + A transcription of ROOT's loop, not a textbook kernel smoother: the window of an output point starts at a running + index that only moves forward and may still include the last point below x0 - cutoff (ROOT sums it), and ends at + the last point <= x0 + cutoff; where no point contributes, the result is 0. The sums run in ROOT's order with + libm's exp, so that the result is bitwise ROOT's. `y` may have a second axis (samples), smoothed in one go. + """ + x_out = np.asarray(x_out, dtype=np.float64) + if np.any(np.diff(x_out) < 0): + raise ValueError("x_out must be sorted") + order = np.argsort(x, kind="stable") + x, y = np.asarray(x, dtype=np.float64)[order], np.asarray(y, dtype=np.float64)[order] + width = KERNEL_SCALE * bandwidth + cutoff = 4 * width + result = np.zeros((len(x_out), *y.shape[1:])) + start = int(np.searchsorted(x, x_out[0] - cutoff, side="left")) + points = x.tolist() + for j, x0 in enumerate(x_out.tolist()): + end = int(np.searchsorted(x, x0 + cutoff, side="right")) + numerator, denominator = 0.0, 0.0 + for i in range(start, end): + distance = abs(points[i] - x0) / width + weight = math.exp(-0.5 * distance * distance) + numerator = numerator + weight * y[i] + denominator += weight + if denominator > 0: + result[j] = numerator / denominator + start = max(start, int(np.searchsorted(x, x0 - cutoff, side="left")) - 1) + return result + + +def binwise(y: np.ndarray, errors: np.ndarray, bin_edges: np.ndarray, stat_sigma: float = 1.0) -> Curve: + """The bins themselves as the curve, with their statistical errors; negative values are set to 0.""" + y, errors, bin_edges = (np.asarray(a, dtype=np.float64) for a in (y, errors, bin_edges)) + return _clipped(bin_edges, y, y + stat_sigma * errors, y - stat_sigma * errors, bin_edges) + + +def smooth(x: np.ndarray, y: np.ndarray, errors: np.ndarray, bin_edges: np.ndarray, bandwidth: float, binwise_left: int = 0, binwise_right: int = 0, stat_sigma: float = 1.0) -> Curve: + """The smoothed curve of a histogram with values `y`, symmetric `errors`, bin `x` positions (the centres of mass) + and `bin_edges`; the first `binwise_left` and the last `binwise_right` bins keep their bin values (a hybrid). + Negative values are set to 0, the curve covers the full bin range.""" + x, y, errors, bin_edges = (np.asarray(a, dtype=np.float64) for a in (x, y, errors, bin_edges)) + start, end = binwise_left, len(y) - 1 - binwise_right # the first and the last smoothed bin + rng = np.random.RandomState(SEED) # the stream of TauFakeFactors' np.random.seed(19), without global state + samples = np.array([rng.normal(mu, stat_sigma * sigma, SAMPLES) for mu, sigma in zip(y, errors)]) + samples[samples < 0.0] = 0.0 + n = 100 * len(y) + grid = np.linspace(x[start], x[end], n + 1) + smoothed = smooth_kern(x, samples, grid[:n], bandwidth) + mean, std = smoothed.mean(axis=1), smoothed.std(axis=1) + if not binwise_left and not binwise_right: + return _clipped(grid, mean, mean + std, mean - std, bin_edges) + left = slice(None, start) + right = slice(-binwise_right, None) if binwise_right else slice(0, 0) # counted from the end, as in TauFakeFactors + overlap = slice(None, -1) if binwise_right else slice(None, None) # the grid ends where the right bins begin + parts = [(bin_edges[left], y[left], stat_sigma * errors[left])] + if binwise_left: # the smoothed part begins at the lower edge of bin `start`, with its first value + parts.append((bin_edges[start : start + 1], mean[:1], std[:1])) + parts.append((grid[overlap], mean[overlap], std[overlap])) + if binwise_right: # and ends at the lower edge of bin `end + 1`, with its last value + parts.append((bin_edges[end + 1 : end + 2], mean[-1:], std[-1:])) + parts.append((bin_edges[right], y[right], stat_sigma * errors[right])) + edges, nominal, band = (np.concatenate(column) for column in zip(*parts)) + return _clipped(edges, nominal, nominal + band, nominal - band, bin_edges) + + +def _clipped(edges: np.ndarray, nominal: np.ndarray, up: np.ndarray, down: np.ndarray, bin_edges: np.ndarray) -> Curve: + """Negative values set to 0; the first and last value extended to the outer bin edges where the curve stops short.""" + nominal, up, down = (np.where(values < 0, 0.0, values) for values in (nominal, up, down)) + if edges[0] != bin_edges[0]: + edges, nominal, up, down = np.insert(edges, 0, bin_edges[0]), *(np.insert(v, 0, v[0]) for v in (nominal, up, down)) + if edges[-1] != bin_edges[-1]: + edges, nominal, up, down = np.append(edges, bin_edges[-1]), *(np.append(v, v[-1]) for v in (nominal, up, down)) + return Curve(edges, nominal, up, down) diff --git a/src/shapesmith/measurements/tau_id_es/__init__.py b/src/shapesmith/measurements/tau_id_es/__init__.py new file mode 100644 index 0000000..4401c80 --- /dev/null +++ b/src/shapesmith/measurements/tau_id_es/__init__.py @@ -0,0 +1,155 @@ +"""Tau-ID scale factor and energy scale of embedded taus, per decay-mode (and pT) category. + +The chain of the predecessor (smhtt_ul tauID_SFs_dev) on ShapeSmith histograms. One run per working-point +combination fills the mt channel and its mm control region (the analysis builds them for `vsjet_wp` and `vsele_wp`), +adds the estimates, writes the shapes in the MorphingTauID2017 format and runs per category in CMSSW (combine.py); +the result is the singles fit. `--merge` combines the runs into the correctionlib payload. + +Output, below /tau_id_es//: +- _/: shapes.root (histograms and estimates), synced/ (MorphingTauID2017 inputs), fits/ (per + category the cards, workspace, scans and singles fit), results.json, plots/; +- DeepTau2018v2p5_id_es_embeddingUL.json.gz, written by --merge. +A result with a problem (an interval at the fit range, a failed crossing, a CMSSW step that failed or ran longer than +STEP_TIMEOUT) is an error: the run raises after writing results.json, and --merge refuses it. +""" +from __future__ import annotations + +import json +import logging +import subprocess +from concurrent.futures import ThreadPoolExecutor +from dataclasses import asdict, dataclass +from pathlib import Path + +from shapesmith import cmssw, payloads +from shapesmith.config import CombineConfig +from shapesmith.estimates import run_estimates +from shapesmith.fill import run_hist +from shapesmith.measurements import MeasureContext +from shapesmith.measurements.tau_id_es import combine, plots +from shapesmith.measurements.tau_id_es.combine import Interval +from shapesmith.measurements.tau_id_es.grid import CATEGORIES, CHANNEL, CONTROL_CATEGORY, CONTROL_CHANNEL, SIGNAL, grid_name +from shapesmith.measurements.tau_id_es.payload import correction_set, payload_name +from shapesmith.measurements.tau_id_es.synced import write_synced +from shapesmith.model import Analysis + +logger = logging.getLogger(__name__) +STEP_TIMEOUT = 3600.0 # seconds per CMSSW step of a category; a singles fit can iterate without end + + +@dataclass(frozen=True) +class TauIdEsMeasurement: + vsjet_wp: str + vsele_wp: str + grid: tuple[int, ...] # the ES shifts in 0.1 % (grid.grid); each is a template variation of EMB in mt + name: str = "tau_id_es" + + def run(self, context: MeasureContext) -> None: + if context.suggest_binning: + raise ValueError("tau_id_es: --suggest-binning is not supported, the m_vis bins are fixed in the analysis") + if context.merge: + merge(context) + else: + measure(self, context) + + +def check(analysis: Analysis, measurement: TauIdEsMeasurement) -> None: + """The analysis has the channels, categories, signal and grid variations the chain expects.""" + problems = [f"channel {name} is missing" for name in (CHANNEL, CONTROL_CHANNEL) if name not in analysis.channels] + if not problems: + channel = analysis.channel(CHANNEL) + problems += [f"unknown mt category {c.name} (known: {', '.join(CATEGORIES)})" for c in channel.categories if c.name not in CATEGORIES] + if [c.name for c in analysis.channel(CONTROL_CHANNEL).categories] != [CONTROL_CATEGORY]: + problems.append(f"the {CONTROL_CHANNEL} channel needs exactly the category {CONTROL_CATEGORY}") + if SIGNAL not in channel.backgrounds(): + problems.append(f"process {SIGNAL} is missing in {CHANNEL}") + declared = {v.name for v in channel.variations} + problems += [f"grid variation {grid_name(shift)} is missing in {CHANNEL}" for shift in measurement.grid if grid_name(shift) not in declared] + if problems: + raise ValueError("tau_id_es: " + "; ".join(problems)) + + +def fit_category(combine_config: CombineConfig | None, category: str, synced: Path, directory: Path, era: str, grid: tuple[int, ...]) -> dict: + """Datacards, workspace, 2D scan, close-up scans, singles fit and postfit shapes of one category below + `directory`/fits; its result and problems.""" + work_dir = Path(directory) / "fits" + work_dir.mkdir(parents=True, exist_ok=True) + cmssw.run([f"{combine.morphing_command(category, synced, era, grid)} > morphing_{category}.log 2>&1"], combine_config, work_dir, STEP_TIMEOUT) + cards = combine.card_dir(work_dir, category) + combine.add_rate_parameters(cards, category) + cmssw.run([f"{combine.workspace_command(category)} > workspace.log 2>&1", f"{combine.scan_command(category, grid)} > scan.log 2>&1"], combine_config, cards, STEP_TIMEOUT) + scan = combine.read_scan(combine.output_path(cards, f"scan_2D_{category}"), category) + closeup = combine.closeup_inputs(scan, grid) + cmssw.run([f"{combine.closeup_command(category, poi, *closeup)} > closeup_{poi}.log 2>&1" for poi in combine.pois(category)], combine_config, cards, STEP_TIMEOUT) + profiles = tuple(combine.read_profile(combine.output_path(cards, f"closeup_{poi}"), poi) for poi in combine.pois(category)) + start, r_range, es_range = combine.singles_inputs(profiles, grid) + cmssw.run([f"{combine.singles_command(category, start, r_range, es_range)} > singles.log 2>&1", f"{combine.postfit_command(category)} > postfit.log 2>&1"], combine_config, cards, STEP_TIMEOUT) + sf, es = combine.read_singles(combine.output_path(cards, f"singles_{category}"), category) + problems = combine.interval_problems("r", sf, r_range) + combine.interval_problems("ES", es, es_range) + return {"sf": asdict(sf), "es": asdict(es), **_inputs(start, r_range, es_range), "closeup": _inputs(*closeup), "problems": problems} + + +def _inputs(start: tuple[float, float], r_range: tuple[float, float], es_range: tuple[float, float]) -> dict: + return {"start": start, "ranges": {"r": r_range, "ES": es_range}} + + +def fit_or_failure(combine_config: CombineConfig | None, category: str, synced: Path, directory: Path, era: str, grid: tuple[int, ...]) -> dict: + """fit_category, or a result without values whose problem is the failed CMSSW step (see its log).""" + try: + return fit_category(combine_config, category, synced, directory, era, grid) + except (subprocess.CalledProcessError, subprocess.TimeoutExpired) as error: + return {"sf": None, "es": None, "problems": [f"CMSSW step failed: {error}"]} + + +def plot_category(directory: Path, category: str, fit: dict) -> None: + """The 2D scan with the singles result and the postfit shapes of one fitted category.""" + cards = combine.card_dir(Path(directory) / "fits", category) + scan = combine.read_scan(combine.output_path(cards, f"scan_2D_{category}"), category) + plots.plot_scan(scan, Interval(**fit["sf"]), Interval(**fit["es"]), category, Path(directory) / "plots" / "scans" / category) + plots.plot_postfit(cards / "postfit_shapes.root", Path(directory) / "plots" / "postfit" / category) + + +def measure(measurement: TauIdEsMeasurement, context: MeasureContext) -> Path: + analysis, config = context.analysis, context.config + check(analysis, measurement) + directory = context.output / f"{measurement.vsjet_wp}_{measurement.vsele_wp}" + shapes = directory / "shapes.root" + hset = run_hist(config, analysis, [CHANNEL, CONTROL_CHANNEL], False, None, True, None, shapes) + run_estimates(hset, analysis, [CHANNEL, CONTROL_CHANNEL]) + hset.save(shapes) + synced = write_synced(hset, analysis, directory / "synced") + categories = [c.name for c in analysis.channel(CHANNEL).categories] + for category in categories: + plots.plot_control(hset, analysis, category, directory / "plots" / "control" / category) + with ThreadPoolExecutor(max_workers=max(1, config.workers)) as pool: # each category runs in its own CMSSW process + fits = list(pool.map(lambda category: fit_or_failure(config.combine, category, synced, directory, analysis.era, measurement.grid), categories)) + for category, fit in zip(categories, fits): + if fit["sf"] is not None: + plot_category(directory, category, fit) + record = { + "vsjet_wp": measurement.vsjet_wp, + "vsele_wp": measurement.vsele_wp, + "categories": dict(zip(categories, fits)), + "provenance": context.provenance(vsjet_wp=measurement.vsjet_wp, vsele_wp=measurement.vsele_wp), + } + (directory / "results.json").write_text(json.dumps(record, indent=2)) + problems = [f"{category}: {problem}" for category, fit in record["categories"].items() for problem in fit["problems"]] + if problems: + raise ValueError(f"tau_id_es {measurement.vsjet_wp}/{measurement.vsele_wp}, see {directory / 'results.json'}:\n" + "\n".join(problems)) + return directory + + +def merge(context: MeasureContext) -> Path: + """The payload of every working-point combination measured below the output directory.""" + records = [json.loads(path.read_text()) for path in sorted(Path(context.output).glob("*/results.json"))] + if not records: + raise ValueError(f"tau_id_es: no results.json below {context.output}; run the measurement per working point first") + problems = [f"{r['vsjet_wp']}_{r['vsele_wp']}/{category}: {problem}" for r in records for category, fit in r["categories"].items() for problem in fit["problems"]] + if problems: + raise ValueError("tau_id_es: results with problems, not merged:\n" + "\n".join(problems)) + results = {(r["vsjet_wp"], r["vsele_wp"]): {category: (Interval(**fit["sf"]), Interval(**fit["es"])) for category, fit in r["categories"].items()} for r in records} + sources = {f"{r['vsjet_wp']}_{r['vsele_wp']}": r["provenance"] for r in records} + cset = correction_set(results, context.provenance(working_points=sorted(sources), measurements=sources)) + path = payloads.write(cset, Path(context.output) / payload_name(context.analysis.era)) + logger.info(f"payload of {len(results)} working-point combinations: {path}") + return path diff --git a/src/shapesmith/measurements/tau_id_es/combine.py b/src/shapesmith/measurements/tau_id_es/combine.py new file mode 100644 index 0000000..509133f --- /dev/null +++ b/src/shapesmith/measurements/tau_id_es/combine.py @@ -0,0 +1,239 @@ +"""The CMSSW part of the measurement, per category, with the options of the predecessor (smhtt_ul tauID_SFs_dev, +tau_id_es_measurement/emb_tau_id_sfs_ul.sh): MorphingTauID2017 datacards plus the r_DY_incl rate parameters, the +multiSignalModel workspace, the 2D likelihood scan over (r_EMB, ES), the close-up 1D scans of each POI with the other +profiled and the MultiDimFit singles fit, which is the result. As in the predecessor (plot_2D_scan.py), each step +starts at the lowest grid point of the previous one and its ranges are the 2 sigma intervals of that step, widened +by a margin: the profiles of the 2D scan give the close-up ranges (the predecessor took them from 1D scans over the +full range, some with margins set by hand), the close-up scans give the singles ranges. + +The POIs: r_EMB_ (the scale factor) and ES_ (the energy scale in percent). +""" +from __future__ import annotations + +from dataclasses import dataclass +from pathlib import Path + +import numpy as np +import uproot + +from shapesmith.measurements.tau_id_es.grid import CHANNEL, CONTROL_CHANNEL, SIGNAL +from shapesmith.measurements.tau_id_es.synced import POSTFIX + +MASS = "125" +R_RANGE = (0.1, 2.9) # of the workspace +SCAN_R_RANGE = (0.11, 2.89) +SCAN_POINTS = 289 +CLOSEUP_POINTS = 17 +LEVEL = 6.18 / 2 # delta NLL of the predecessor's 2 sigma intervals (2 delta NLL = 6.18, the 2 sigma level of two parameters) +CLOSEUP_MARGIN = (0.2, 2.0) # (r, ES): widening of the 2D scan's intervals to the close-up ranges +SINGLES_MARGIN = (0.05, 0.5) # (r, ES): widening of the close-up scans' intervals to the singles ranges +ROBUST_FIT = "--robustFit=1 --setRobustFitAlgo=Minuit2 --X-rtd FITTER_NEW_CROSSING_ALGO --X-rtd FITTER_NEVER_GIVE_UP" +SCAN_FIT = f"{ROBUST_FIT} --cminFallbackAlgo Minuit2,Migrad,0:0.001 --cminFallbackAlgo Minuit2,Migrad,1:0.01 --cminPreScan" +Pair = tuple[float, float] + + +@dataclass(frozen=True) +class Interval: + best: float + low: float + high: float + + +@dataclass(frozen=True) +class Scan: + """The grid points of the 2D scan (r, ES, delta NLL).""" + + r: np.ndarray + es: np.ndarray + dnll: np.ndarray + + def profile(self, axis: int) -> Profile: + """The profile of r (axis 0) or ES (axis 1) on the grid: the lowest delta NLL over the other POI.""" + values = (self.r, self.es)[axis] + dnll = np.where(np.isfinite(self.dnll), self.dnll, np.inf) + grid = np.unique(values) + return Profile.of(grid, np.array([dnll[values == value].min() for value in grid])) + + +@dataclass(frozen=True) +class Profile: + """A 1D likelihood scan: sorted grid values and their delta NLL, relative to the fit or, if a grid point lies + lower (a fit that missed the minimum), to that point.""" + + values: np.ndarray + dnll: np.ndarray + + @classmethod + def of(cls, values: np.ndarray, dnll: np.ndarray) -> Profile: + finite = np.isfinite(dnll) + order = np.argsort(values[finite], kind="stable") + values, dnll = values[finite][order], dnll[finite][order] + return cls(values, dnll - min(0.0, dnll.min())) + + @property + def minimum(self) -> float: + """The lowest grid point, the start value of the next step (the predecessor's workaround for failed fits).""" + return float(self.values[np.argmin(self.dnll)]) + + def crossings(self, level: float) -> tuple[float | None, float | None]: + """The outermost crossings of `level` on each side of the minimum, linearly interpolated; None if there is + none (plot_2D_scan.extract_confidence_interval).""" + x, y, best = self.values, self.dnll, int(np.argmin(self.dnll)) + + def crossing(i: int) -> float: # on the segment from grid point i to i + 1 + return float(x[i] + (level - y[i]) * (x[i + 1] - x[i]) / (y[i + 1] - y[i])) + + low = [i for i in range(best) if y[i] > level > y[i + 1]] + high = [i for i in range(best, len(x) - 1) if y[i] < level < y[i + 1]] + return (crossing(low[0]) if low else None), (crossing(high[-1]) if high else None) + + def window(self, margin: float, bounds: Pair) -> Pair: + """The 2 sigma interval (an end of the scan where it does not cross), widened by `margin`, within `bounds`.""" + low, high = self.crossings(LEVEL) + low, high = self.values[0] if low is None else low, self.values[-1] if high is None else high + return max(float(low) - margin, bounds[0]), min(float(high) + margin, bounds[1]) + + +def pois(category: str) -> tuple[str, str]: + return f"r_{SIGNAL}_{category}", f"ES_{category}" + + +def es_range(grid: tuple[int, ...]) -> tuple[float, float]: + """The ES range of the datacards (percent): the ends of the grid.""" + return min(grid) / 10, max(grid) / 10 + + +def scan_es_range(grid: tuple[int, ...]) -> tuple[float, float]: + low, high = es_range(grid) + return low + 0.1, high - 0.1 + + +def morphing_command(category: str, synced_dir: Path, era: str, grid: tuple[int, ...]) -> str: + """MorphingTauID2017 as in the predecessor; it writes output//htt_mt_/ below the working + directory. Its ES loop runs in float from es_min in steps of the precision and skips 0; the accumulated rounding + drops the last grid point (with the +-20 %, 0.2 % grid of the predecessor: +20.0).""" + steps = set(np.diff(sorted(grid + (0,)))) + if len(steps) != 1: + raise ValueError(f"MorphingTauID2017 needs a uniform grid, got steps {sorted(steps)}") + low, high = es_range(grid) + options = { + "base_path": Path(synced_dir).resolve(), "input_folder_mt": CHANNEL, "input_folder_mm": CONTROL_CHANNEL, + "real_data": "true", "classic_bbb": "false", "binomial_bbb": "false", "jetfakes": 0, "embedding": 1, + "verbose": "false", "postfix": POSTFIX, "use_control_region": "true", "auto_rebin": "false", + "rebin_categories": "false", "manual_rebin_for_yields": "false", "categories": category, "era": era, + "tes_precision": steps.pop() / 10, "es_min": low, "es_max": high, "output": category, + } + return "MorphingTauID2017 " + " ".join(f"--{name}={value}" for name, value in options.items()) + + +def card_dir(work_dir: Path, category: str) -> Path: + return Path(work_dir) / "output" / category / f"htt_{CHANNEL}_{category}" + + +def add_rate_parameters(cards: Path, category: str) -> None: + """The common normalisation r_DY_incl_ of the Z->ll processes in both regions: EMB, ZL and ZJ in mt, + the embedded muons in mm.""" + line = f"r_DY_incl_{category} rateParam * {{}} 1.0 [0.5,1.5]\n" + for card in sorted(Path(cards).glob(f"htt_{CHANNEL}_*.txt")): + processes = [f"{SIGNAL}_{category}", "ZL", "ZJ"] if "ZL" in card.read_text() else [f"{SIGNAL}_{category}", "ZJ"] + with card.open("a") as handle: + handle.writelines(line.format(process) for process in processes) + for card in sorted(Path(cards).glob(f"htt_{CONTROL_CHANNEL}_*.txt")): + with card.open("a") as handle: + handle.write(line.format("MUEMB")) + + +def workspace_command(category: str) -> str: + r = pois(category)[0] + return (f"combineTool.py -M T2W -i . -o workspace.root -m {MASS} -P HiggsAnalysis.CombinedLimit.PhysicsModel:multiSignalModel" + f" --PO \"map=^.*/{SIGNAL}_{category}:{r}[1,{R_RANGE[0]},{R_RANGE[1]}]\"") + + +def _ranges(category: str, r_range: Pair, es: Pair) -> str: + r, es_poi = pois(category) + return f"--setParameterRanges {r}={r_range[0]!r},{r_range[1]!r}:{es_poi}={es[0]!r},{es[1]!r}" + + +def scan_command(category: str, grid: tuple[int, ...]) -> str: + r, es = pois(category) + return (f"combineTool.py -M MultiDimFit -n .scan_2D_{category} -d workspace.root --setParameters {r}=1.0,{es}=0.0" + f" {_ranges(category, SCAN_R_RANGE, scan_es_range(grid))} {SCAN_FIT}" + f" --redefineSignalPOIs {r},{es} --floatOtherPOIs=1 --points={SCAN_POINTS} --algo grid -m {MASS}" + " --alignEdges=1 --cminDefaultMinimizerStrategy 0") + + +def closeup_command(category: str, poi: str, start: Pair, r_range: Pair, es: Pair) -> str: + """The close-up scan of `poi`, the other POI profiled.""" + r, es_poi = pois(category) + return (f"combineTool.py -M MultiDimFit -n .closeup_{poi} -d workspace.root --setParameters {es_poi}={start[1]!r},{r}={start[0]!r}" + f" {_ranges(category, r_range, es)} {SCAN_FIT} --redefineSignalPOIs {poi} --algo grid -m {MASS}" + f" --cminDefaultMinimizerStrategy 0 --floatOtherPOIs=1 --points={CLOSEUP_POINTS} --alignEdges=1") + + +def singles_command(category: str, start: Pair, r_range: Pair, es: Pair) -> str: + r, es_poi = pois(category) + return (f"combineTool.py -M MultiDimFit -n .singles_{category} -d workspace.root" + f" --setParameters {es_poi}={start[1]!r},{r}={start[0]!r} {_ranges(category, r_range, es)} {ROBUST_FIT}" + " --cminFallbackAlgo Minuit2,Migrad,0:0.001,Minuit2,Migrad,0:0.01 --cminPreScan" + f" --redefineSignalPOIs {r},{es_poi} --floatOtherPOIs=1 --algo singles -m {MASS} --saveFitResult" + " --cminDefaultMinimizerStrategy=1") + + +def postfit_command(category: str) -> str: + return (f"PostFitShapesFromWorkspace -m {MASS} -w workspace.root --output postfit_shapes.root" + f" -f multidimfit.singles_{category}.root:fit_mdf --postfit") + + +def output_path(cards: Path, name: str) -> Path: + return Path(cards) / f"higgsCombine.{name}.MultiDimFit.mH{MASS}.root" + + +def read_scan(path: Path, category: str) -> Scan: + """The grid points of the 2D scan (its first entry is the fit).""" + r, es = pois(category) + with uproot.open(path) as f: + rows = f["limit"].arrays([r, es, "deltaNLL"], library="np") + return Scan(rows[r][1:].astype(float), rows[es][1:].astype(float), rows["deltaNLL"][1:].astype(float)) + + +def read_profile(path: Path, poi: str) -> Profile: + """A close-up scan (its first entry is the fit).""" + with uproot.open(path) as f: + rows = f["limit"].arrays([poi, "deltaNLL"], library="np") + return Profile.of(rows[poi][1:].astype(float), rows["deltaNLL"][1:].astype(float)) + + +def closeup_inputs(scan: Scan, grid: tuple[int, ...]) -> tuple[Pair, Pair, Pair]: + """Start values, r range and ES range of the close-up scans, from the profiles of the 2D scan.""" + return _next_inputs((scan.profile(0), scan.profile(1)), CLOSEUP_MARGIN, grid) + + +def singles_inputs(closeups: tuple[Profile, Profile], grid: tuple[int, ...]) -> tuple[Pair, Pair, Pair]: + """Start values, r range and ES range of the singles fit, from the close-up scans of r and ES.""" + return _next_inputs(closeups, SINGLES_MARGIN, grid) + + +def _next_inputs(profiles: tuple[Profile, Profile], margins: Pair, grid: tuple[int, ...]) -> tuple[Pair, Pair, Pair]: + start = (profiles[0].minimum, profiles[1].minimum) + r_range, es = (p.window(margin, bounds) for p, margin, bounds in zip(profiles, margins, (SCAN_R_RANGE, scan_es_range(grid)))) + return start, r_range, es + + +def read_singles(path: Path, category: str) -> tuple[Interval, Interval]: + """The singles fit: best fit and 1 sigma interval of r and ES (the limit tree holds the best fit, then the low + and high crossing of each POI with the others at their best fit).""" + names = pois(category) + with uproot.open(path) as f: + rows = f["limit"].arrays(list(names), library="np") + return tuple(Interval(float(rows[name][0]), float(rows[name][1:].min()), float(rows[name][1:].max())) for name in names) + + +def interval_problems(name: str, interval: Interval, fit_range: tuple[float, float]) -> list[str]: + """An interval that failed (an end equals the best fit) or ends at the fit range is a problem.""" + tolerance = 1e-5 * (fit_range[1] - fit_range[0]) + problems = [] + if not interval.low < interval.best < interval.high: + problems.append(f"{name}: no interval around the best fit {interval.best:g} ({interval.low:g}, {interval.high:g})") + if interval.low - fit_range[0] < tolerance or fit_range[1] - interval.high < tolerance: + problems.append(f"{name}: interval [{interval.low:g}, {interval.high:g}] at the fit range [{fit_range[0]:g}, {fit_range[1]:g}]") + return problems diff --git a/src/shapesmith/measurements/tau_id_es/grid.py b/src/shapesmith/measurements/tau_id_es/grid.py new file mode 100644 index 0000000..0c2738d --- /dev/null +++ b/src/shapesmith/measurements/tau_id_es/grid.py @@ -0,0 +1,54 @@ +"""The tau energy-scale grid and the naming conventions of MorphingTauID2017. + +A grid point is an energy-scale shift in integer units of 0.1 %: shift -198 scales the tau four-momentum by 0.802. +Its column variation is named `es-198`, and MorphingTauID2017 reads its shape under the mass string "-19.8" (the +shift in percent with one decimal). Integers keep the three names in exact correspondence. +""" +from __future__ import annotations + +import re + +SIGNAL = "EMB" # the embedded signal of the mt channel, one shape per grid point +CHANNEL = "mt" +CONTROL_CHANNEL = "mm" # the Z->mumu control region, one category +CONTROL_CATEGORY = "control_region" + +# MorphingTauID2017's categories: decay modes and tau pT bin in GeV (None: every tau above 20 GeV) +CATEGORIES = { + "DM0": ((0,), None), + "DM1": ((1,), None), + "DM1011": ((10, 11), None), + "DM0_PT20_40": ((0,), (20.0, 40.0)), + "DM1_PT20_40": ((1,), (20.0, 40.0)), + "DM1011_PT20_40": ((10, 11), (20.0, 40.0)), + "DM0_PT40_200": ((0,), (40.0, 200.0)), + "DM1_PT40_200": ((1,), (40.0, 200.0)), + "DM1011_PT40_200": ((10, 11), (40.0, 200.0)), +} + +_NAME = re.compile(r"es([+-]\d+)") + + +def grid(low: int, high: int, step: int) -> tuple[int, ...]: + """The shifts from `low` to `high` (inclusive) in steps of `step`, without the nominal 0.""" + return tuple(shift for shift in range(low, high + 1, step) if shift != 0) + + +def grid_name(shift: int) -> str: + return f"es{shift:+d}" + + +def shift_of(name: str) -> int | None: + """The shift of a grid variation name, None for any other name.""" + match = _NAME.fullmatch(name) + return int(match.group(1)) if match else None + + +def factor(shift: int) -> float: + """The scale factor of the tau four-momentum.""" + return (1000 + shift) / 1000 + + +def mass(shift: int) -> str: + """MorphingTauID2017's mass string of a shift: the percentage with one decimal.""" + return f"{shift / 10:.1f}" diff --git a/src/shapesmith/measurements/tau_id_es/payload.py b/src/shapesmith/measurements/tau_id_es/payload.py new file mode 100644 index 0000000..8bd11eb --- /dev/null +++ b/src/shapesmith/measurements/tau_id_es/payload.py @@ -0,0 +1,113 @@ +"""The correctionlib payload of the measurement, with the node structure of the predecessor's +friends/create_xpog_json_v15.py (the file the embedding configuration of CROWN reads), plus provenance. + +- DeepTau2018v2p5VSjet(pt, dm, genmatch, wp, wp_VSe, syst, flag): the scale factor r; flag "pt" in the pT bins of + the categories, flag "dm" per decay mode for every pT from 20 GeV. +- tau_energy_scale and tau_energy_scale_dm_binned(pt, eta, dm, genmatch, id, wp, wp_VSe, syst): the energy scale + factor (100 + ES) / 100 of the fitted ES in percent, per decay mode and pT bin or per decay mode, for |eta| < 2.5. +Every other genmatch gives 1.0. DM10 and DM11 both take the DM1011 result; syst up/down are the ends of the 1 sigma +interval. +""" +from __future__ import annotations + +from typing import Callable, Mapping + +import correctionlib.schemav2 as cs + +from shapesmith import payloads +from shapesmith.measurements.tau_id_es.combine import Interval +from shapesmith.measurements.tau_id_es.grid import CATEGORIES + +DECAY_MODES = (0, 1, 10, 11) +OTHER_GENMATCH = (0, 1, 2, 3, 4, 6) +WORKING_POINTS = ("VVVLoose", "VVLoose", "VLoose", "Loose", "Medium", "Tight", "VTight", "VVTight") +WP_ORDER = {wp: index for index, wp in enumerate(WORKING_POINTS)} +MAX_PT = 100000.0 # the upper edge of the per-decay-mode bins +TAU_ID = "DeepTau2018v2p5" + +# per working-point combination (vsjet, vsele) and category: the intervals of the scale factor and of ES (percent) +Results = Mapping[tuple[str, str], Mapping[str, tuple[Interval, Interval]]] + + +def es_factor(percent: float) -> float: + return (100 + percent) / 100 + + +def _syst(interval: Interval, convert: Callable[[float], float]) -> cs.Category: + return payloads.category("syst", {"nom": convert(interval.best), "up": convert(interval.high), "down": convert(interval.low)}) + + +def _pt_bins(dm: int, results: Mapping[str, tuple[Interval, Interval]], pt_binned: bool) -> list[tuple[tuple[float, float], str]]: + """The (pT bin, category) of a decay mode: its pT categories, or its inclusive category from 20 GeV.""" + bins = [(pt, name) for name, (modes, pt) in CATEGORIES.items() if dm in modes and (pt is not None) == pt_binned and name in results] + if not pt_binned: + return [((20.0, MAX_PT), name) for _, name in bins] + return sorted(bins) + + +def _by_decay_mode(results: Mapping[str, tuple[Interval, Interval]], index: int, pt_binned: bool, convert: Callable[[float], float], wrap: Callable = lambda node: node) -> cs.Category: + content = {} + for dm in DECAY_MODES: + bins = _pt_bins(dm, results, pt_binned) + if bins: + edges = [bins[0][0][0]] + [pt[1] for pt, _ in bins] + content[dm] = wrap(payloads.binning("pt", edges, [_syst(results[name][index], convert) for _, name in bins])) + return payloads.category("dm", content) + + +def _genuine_taus(results: Results, node: Callable[[Mapping[str, tuple[Interval, Interval]]], cs.Category]) -> cs.Category: + """1.0 unless genmatch is 5; then per working point, vsEle working point and decay mode.""" + by_wp: dict[str, dict] = {} + for vsjet, vsele in sorted(results, key=lambda pair: (WP_ORDER[pair[0]], WP_ORDER[pair[1]])): + by_wp.setdefault(vsjet, {})[vsele] = node(results[vsjet, vsele]) + wps = payloads.category("wp", {vsjet: payloads.category("wp_VSe", by_vsele) for vsjet, by_vsele in by_wp.items()}) + return payloads.category("genmatch", {**{genmatch: 1.0 for genmatch in OTHER_GENMATCH}, 5: wps}) + + +def _inputs(results: Results, *names: str) -> list[cs.Variable]: + vsjet = ",".join(sorted({pair[0] for pair in results}, key=WP_ORDER.get)) + vsele = ",".join(sorted({pair[1] for pair in results}, key=WP_ORDER.get)) + known = { + "pt": ("real", "Reconstructed tau pT"), + "eta": ("real", "Reconstructed tau eta"), + "dm": ("int", "Reconstructed tau decay mode: 0, 1, 10, 11"), + "genmatch": ("int", "genmatch: 0 or 6 = unmatched or jet, 1 or 3 = electron, 2 or 4 = muon, 5 = real tau"), + "id": ("string", f"{TAU_ID}:{TAU_ID}"), + "wp": ("string", f"{TAU_ID}VSjet working point: {vsjet}"), + "wp_VSe": ("string", f"{TAU_ID}VSe working point: {vsele}"), + "syst": ("string", "Systematic variations: 'nom', 'up', 'down'"), + "flag": ("string", "Flag: 'pt' = DM-pT-dependent SFs(20,40,200), 'dm' = DM-dependent SFs with pt>=20"), + } + return [payloads.variable(name, *known[name]) for name in names] + + +def _abs_eta(node: cs.Binning) -> cs.Transform: + return payloads.transform("eta", "abs(x)", payloads.binning("eta", [0.0, 2.5], [node])) + + +def scale_factors(results: Results) -> cs.Correction: + data = payloads.category("flag", { + "pt": _genuine_taus(results, lambda r: _by_decay_mode(r, 0, True, float)), + "dm": _genuine_taus(results, lambda r: _by_decay_mode(r, 0, False, float)), + }) + output = payloads.variable("sf", "real", "DM-pT-dependent or DM-dependent scale factor") + return payloads.correction(f"{TAU_ID}VSjet", _inputs(results, "pt", "dm", "genmatch", "wp", "wp_VSe", "syst", "flag"), data, version=0, + description="dm-pt-dependent tau ID scale factor for tau embedded samples", output=output) + + +def energy_scale(results: Results, pt_binned: bool) -> cs.Correction: + data = payloads.category("id", {TAU_ID: _genuine_taus(results, lambda r: _by_decay_mode(r, 1, pt_binned, es_factor, _abs_eta))}) + inputs = _inputs(results, "pt", "eta", "dm", "genmatch", "id", "wp", "wp_VSe", "syst") + if pt_binned: + output = payloads.variable("tes", "real", "DM-pT-dependent tau energy scale") + return payloads.correction("tau_energy_scale", inputs, data, version=0, description="dm-pt-dependent tau ES scale factor for tau embedded samples", output=output) + output = payloads.variable("tes", "real", "DM-dependent tau energy scale") + return payloads.correction("tau_energy_scale_dm_binned", inputs, data, version=0, description="dm-dependent tau ES scale factor for tau embedded samples", output=output) + + +def correction_set(results: Results, provenance: Mapping) -> cs.CorrectionSet: + return payloads.correction_set([scale_factors(results), energy_scale(results, True), energy_scale(results, False)], provenance) + + +def payload_name(era: str) -> str: + return f"{TAU_ID}_id_es_embedding{era}UL.json.gz" diff --git a/src/shapesmith/measurements/tau_id_es/plots.py b/src/shapesmith/measurements/tau_id_es/plots.py new file mode 100644 index 0000000..8f6d578 --- /dev/null +++ b/src/shapesmith/measurements/tau_id_es/plots.py @@ -0,0 +1,107 @@ +"""Plots of the measurement: the 2D likelihood scan per category, the postfit shapes of the singles fit, and prefit +control plots of m_vis with a family of energy-scale grid points.""" +from __future__ import annotations + +from pathlib import Path + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt # noqa: E402 +import mplhep as hep # noqa: E402 +import numpy as np # noqa: E402 +import uproot # noqa: E402 + +from shapesmith.histogram import NOMINAL_VARIATION, HistKey, HistogramSet # noqa: E402 +from shapesmith.measurements.tau_id_es.combine import Interval, Scan # noqa: E402 +from shapesmith.measurements.tau_id_es.grid import CHANNEL, SIGNAL, grid_name # noqa: E402 +from shapesmith.model import NOMINAL, Analysis # noqa: E402 + +hep.style.use("CMS") +ES_FAMILY = (-200, -100, 100, 200) # grid points drawn in the control plots, if present +LEVELS = {2.30: "68 % CL", 5.99: "95 % CL"} # 2 delta NLL of two parameters + + +def _colors(n: int) -> list: + """Distinct fill colours for up to 20 stacked processes (the default cycle repeats after 10).""" + return [plt.cm.tab20(index % 20) for index in range(n)] + + +def _save(fig, path: Path) -> list[Path]: + path.parent.mkdir(parents=True, exist_ok=True) + paths = [path.with_suffix(ext) for ext in (".pdf", ".png")] + for p in paths: + fig.savefig(p, bbox_inches="tight") + plt.close(fig) + return paths + + +def plot_scan(scan: Scan, sf: Interval, es: Interval, category: str, path: Path) -> list[Path]: + """2 delta NLL of the scan grid with the 68 % and 95 % contours and the singles result.""" + finite = np.isfinite(scan.dnll) + r, e, z = scan.r[finite], scan.es[finite], 2 * (scan.dnll[finite] - scan.dnll[finite].min()) + fig, ax = plt.subplots(figsize=(10, 9)) + filled = ax.tricontourf(r, e, np.minimum(z, 30.0), levels=30, cmap="viridis") + lines = ax.tricontour(r, e, z, levels=sorted(LEVELS), colors=["white", "orange"]) + ax.clabel(lines, fmt={level: label for level, label in LEVELS.items()}, fontsize=14) + ax.errorbar([sf.best], [es.best], xerr=[[sf.best - sf.low], [sf.high - sf.best]], yerr=[[es.best - es.low], [es.high - es.best]], fmt="o", color="red", label="singles fit") + ax.set_xlabel(f"$r_{{EMB}}$ ({category})") + ax.set_ylabel("tau ES shift [%]") + ax.legend(loc="upper right") + fig.colorbar(filled, ax=ax, label=r"$2\,\Delta$NLL") + return _save(fig, path) + + +def plot_postfit(shapes: Path, output_dir: Path) -> list[Path]: + """Data and the stacked postfit processes of every bin of the singles fit, with the data/prediction ratio.""" + written = [] + with uproot.open(shapes) as f: + for directory in (name for name in f.keys(cycle=False, recursive=False) if name.endswith("_postfit")): + processes = [name for name in f[directory].keys(cycle=False) if not name.startswith(("Total", "data_obs"))] + edges = f[f"{directory}/data_obs"].axis().edges() + stack = [f[f"{directory}/{name}"].values() for name in processes] + data = f[f"{directory}/data_obs"] + total = f[f"{directory}/TotalProcs"] + fig, (ax, rax) = plt.subplots(2, 1, figsize=(10, 10), gridspec_kw={"height_ratios": [3, 1], "hspace": 0.06}, sharex=True) + hep.histplot(stack, bins=edges, stack=True, histtype="fill", label=processes, color=_colors(len(stack)), ax=ax) + hep.histplot(data.values(), bins=edges, yerr=np.sqrt(data.variances()), histtype="errorbar", color="black", label="Data", ax=ax) + ax.set_ylabel("Events") + ax.set_ylim(0, 1.6 * data.values().max()) + ax.legend(ncol=3, fontsize=14, frameon=False, loc="upper right") + ax.text(0.04, 0.95, directory, transform=ax.transAxes, va="top", fontsize=16) + with np.errstate(divide="ignore", invalid="ignore"): + rax.errorbar(0.5 * (edges[1:] + edges[:-1]), data.values() / total.values(), yerr=np.sqrt(data.variances()) / total.values(), fmt="o", color="black") + rax.axhline(1.0, color="gray") + rax.set_ylim(0.8, 1.2) + rax.set_ylabel("Data / fit") + rax.set_xlabel(r"$m_{vis}$ [GeV]") + written += _save(fig, Path(output_dir) / directory) + return written + + +def plot_control(hset: HistogramSet, analysis: Analysis, category: str, path: Path) -> list[Path]: + """Prefit m_vis of one category: data, the stacked backgrounds with the nominal embedded signal, and the total + prediction with the embedded signal at a few grid points.""" + channel = analysis.channel(CHANNEL) + variable = next(c.variable for c in channel.categories if c.name == category) + + def get(process: str, variation: str = NOMINAL_VARIATION): + return hset.get(HistKey(CHANNEL, category, process, NOMINAL, variation, variable.name)) + + backgrounds = [(name, h) for name in channel.backgrounds() if (h := get(name)) is not None] + edges, widths = np.asarray(variable.edges), np.diff(variable.edges) + total = np.sum([h.values for _, h in backgrounds], axis=0) + signal = get(SIGNAL) + fig, ax = plt.subplots(figsize=(10, 8)) + hep.histplot([h.values / widths for _, h in backgrounds], bins=edges, stack=True, histtype="fill", label=[name for name, _ in backgrounds], color=_colors(len(backgrounds)), ax=ax) + for shift in ES_FAMILY: + shifted = get(SIGNAL, grid_name(shift)) + if shifted is not None: + hep.histplot((total - signal.values + shifted.values) / widths, bins=edges, histtype="step", linewidth=2, label=f"ES {shift / 10:+.1f} %", ax=ax) + data = get(channel.data()) + hep.histplot(data.values / widths, bins=edges, yerr=np.sqrt(data.variances) / widths, histtype="errorbar", color="black", label="Data", ax=ax) + ax.set_ylabel("dN / dm$_{vis}$ [1/GeV]") + ax.set_xlabel(r"$m_{vis}$ [GeV]") + ax.legend(ncol=2, fontsize=14, frameon=False) + ax.text(0.04, 0.95, f"{CHANNEL}, {category}", transform=ax.transAxes, va="top", fontsize=16) + return _save(fig, path) diff --git a/src/shapesmith/measurements/tau_id_es/synced.py b/src/shapesmith/measurements/tau_id_es/synced.py new file mode 100644 index 0000000..acb96b3 --- /dev/null +++ b/src/shapesmith/measurements/tau_id_es/synced.py @@ -0,0 +1,41 @@ +"""The shapes in the input format of MorphingTauID2017 (that of the predecessor's convert_to_synced_shapes.py). + +One file per channel, `/htt_.inputs-sm-Run-TauID_ES.root`, with one directory per category, +`_`, holding data_obs and every background, nominal as `` and varied as +`_`. The embedded signal of the mt channel is `EMB__` per grid point, the +nominal with mass "0.0", and its variations `EMB___`. +""" +from __future__ import annotations + +from pathlib import Path + +from shapesmith.histogram import TEMPLATE_SEPARATOR, HistogramSet +from shapesmith.measurements.tau_id_es.grid import SIGNAL, mass, shift_of +from shapesmith.model import Analysis +from shapesmith.shapes import shape_name, synced_entries, write_shapes + +POSTFIX = "-TauID_ES" + + +def shapes_path(directory: Path, channel: str, era: str) -> Path: + return Path(directory) / channel / f"htt_{channel}.inputs-sm-Run{era}{POSTFIX}.root" + + +def signal_name(category: str, variation: str) -> str: + """MorphingTauID2017's name of the embedded signal under a variation (a grid point, one on top of it, or none).""" + systematic, _, template = variation.partition(TEMPLATE_SEPARATOR) + shift = shift_of(template or systematic) + if shift is None: + return shape_name(f"{SIGNAL}_{category}_{mass(0)}", variation) + return f"{SIGNAL}_{category}_{mass(shift)}" + (f"_{systematic}" if template else "") + + +def _name(category: str, process: str, variation: str) -> str: + return signal_name(category, variation) if process == SIGNAL else shape_name(process, variation) + + +def write_synced(hset: HistogramSet, analysis: Analysis, directory: Path) -> Path: + """The shapes file of every channel of the analysis under `directory`; returns `directory`.""" + for channel in analysis.channels: + write_shapes(shapes_path(directory, channel, analysis.era), synced_entries(hset, analysis, channel, _name)) + return Path(directory) diff --git a/src/shapesmith/ml_export.py b/src/shapesmith/ml_export.py index 2a2953f..bbf43ea 100644 --- a/src/shapesmith/ml_export.py +++ b/src/shapesmith/ml_export.py @@ -1,4 +1,4 @@ -"""Training folds from the skims, in the Feather layout of the previous smhtt_ul export (Spec §12). +"""Training folds from the skims, in the Feather layout of the previous smhtt_ul export. Columns are 5-level tuples: ("Event","id"), ("Event","event"), ("Labels",