Running the flowchart

The "flowchart" is the statistical pipeline of the Z′ → bb̄ TLA search: it starts from per-region m(jj) histograms and ends with the 95% CL limit on the coupling g_q versus m_Z′ (the Brazilian band). Everything lives in src/flowcharts/ and is driven from the command line by scripts/run_flowchart.jl, which mirrors notebooks/5-flow_chart.ipynb stage by stage.

Inputs

The pipeline never touches ntuples. It reads ROOT files of per-region m(jj) histograms produced by ensemble_mjj_task (see the private schema page for the ntuples themselves):

file	content
data/flowchart_histograms/data23_<4,20>percent.nv9.root	data, keys SR_nominal, C1_nominal, C0_nominal, D1_nominal, D0_nominal
data/flowchart_histograms/Zprime_bb_mRp<MASS>_mjj.nv9.root	signal per mass, SR_nominal plus SR_<systematic>__1up / __1down
data/sr_lib/*.jl	symbolic-regression background models: a ladder of functions f_c2, f_c4, … of increasing complexity
data/100percent_histograms/	the same for the full 2023 data set (--lumi-pct 100)

Regions are defined by |y*| and the number of GN1 b-tags:

region	|y*|	b-tags	role
SR	< 0.8	2	signal region
C1	< 0.8	1	background model: fit by symbolic regression, downscaled to the SR yield
C0	< 0.8	0	shape cross-check
D1	≥ 0.8	1	shape cross-check
D0	≥ 0.8	0	shape cross-check

The background is data-driven from C1, which is why every systematic in the likelihood is a signal systematic.

Quick start

# full pipeline on the 4% slice with the validated background function (complexity index 2)
julia --project scripts/run_flowchart.jl --lumi-pct 4 --sr-idx 2 --plot-output-dir plots/my_run

# only the limits (dependencies are pulled in automatically, cached stages are reused)
julia --project scripts/run_flowchart.jl --lumi-pct 4 --stages limits --plot-output-dir plots/my_run

# recompute one stage from scratch
julia --project scripts/run_flowchart.jl --lumi-pct 4 --stages spurious --redo --plot-output-dir plots/my_run

# two masses only, quick look
julia --project scripts/run_flowchart.jl --lumi-pct 4 --masses 200,300 --plot-output-dir plots/scratch_200_300

julia --project scripts/run_flowchart.jl --help

Julia 1.11 (julia) and 1.13 (julia +1.13, which picks up Manifest-v1.13.toml) are both supported. Run jobs sequentially: the full 4% pipeline takes of the order of an hour and the toy stages are memory hungry.

Configuration

Each luminosity slice has validated defaults, collected in FlowchartConfig:

--lumi-pct	data file	background model	sr_idx
4	data23_4percent.nv9.root	test_4pct_06_29_80to620GeV.jl	2 (complexity 4)
20	data23_20percent.nv9.root	test_20pct_06_29_80to620GeV.jl	9 (complexity 18)
100	data23_100percent.nv9.root	test_100pct_08_17_80to620GeV.jl	9 (not yet validated)

The two model files have different ladders

sr_idx is an index into the file's sr_tuple_list, and the 4% and 20% files contain different complexities at the same index. Always quote the complexity, not the index, when comparing slices.

Two binnings are in play and must not be confused:

1. SR-search binning, uniform 1 GeV over 80–620 GeV. Only used when fitting the symbolic-regression background (--refit), so the result behaves like a smooth PDF.

2. Likelihood binning (cfg.fit_BINEDGES), the mass-resolution bins [90, 102, 114, …, 611] from calculate_res_bins. Every fit and plot uses it. Override with --fit-bin-edges 98:12:502 or an explicit list.

litehf_fit evaluates the background PDF and the signal templates on a uniform 1 GeV grid spanning the likelihood range and sums them into the resolution bins ("fine + rebinby"), so the background yield is a proper integral without per-bin width factors.

Stages

All stages run by default (--stages all), in this canonical order; --stages takes a subset and adds the dependencies automatically. Aliases in parentheses.

stage	what it does	output
residuals	pulls of several ladder rungs fit to the downscaled C1 (run_residual_plots)	residual_C1.svg
chi2pval (chi2)	χ² and p-value vs complexity (run_chi2_pval_plot)	complexity_pval_C1.svg
chi2pval_atlas (atlasref)	same with the ATLAS 4p/5p dijet functions as reference stars	complexity_pval_ATLAS_C1.svg
crdivergence (crdiv, shape)	KL/JS/Hellinger shape divergence of every CR to the SR	cr_shape_divergence.svg
bonly	background-only fit to the SR data (run_bonly_fit)	bonly_fit_SRidx<i>.svg
spurious (ss)	background-only toy fits per mass, fitted-μ bias (run_spurious_test)	SS_<m>GeV_SRidx<i>.svg, ss_summary_SRidx<i>.svg, spurious_SRidx<i>.csv
injection (si)	signal-injection linearity (run_injection_test)	SI_scan_SRidx<i>.svg
limits	S+B fit to data with all systematics, CLs scan, Brazilian bands (run_limit_scan)	limit_plot*_SRidx<i>.svg
llrscan (cls, llr)	profile-likelihood-ratio scan −2lnλ(μ) per mass	llrscan_<m>GeV_SRidx<i>.svg
pulls (pull)	nuisance pulls and post-fit plots from the limit fits (make_pull_plots)	pulls_<m>GeV_*.svg, postfit_<m>GeV_*.svg, pulls_SRidx<i>_minos.csv
impacts (ranking)	CMS-style impact ranking per mass (run_impacts)	impacts_<m>GeV_SRidx<i>.svg
syst (templates)	nominal vs ±1σ overlay of every signal systematic (dump_syst_templates)	syst_templates/

Opt-in only (never part of all): chi2pval_atlas_poisson (atlastoys) repeats the χ² ladder on --atlas-toys N Poisson toys of C1 downsampled to the SR yield, so the points get SR-statistics error bars (run_chi2_pval_atlas_poisson_plot).

Dependencies: limits ← spurious; pulls, llrscan ← limits, spurious; impacts ← spurious. The spurious stage supplies the per-mass spurious-signal means that become the spurious-signal nuisance in every downstream likelihood.

The likelihood model

Per mass point, litehf_fit maximises a binned Poisson likelihood with

• background: the chosen symbolic-regression function of C1, with a free yield NB and free shape parameters, entering as a normalised PDF;

• signal: the MC template scaled to the slice luminosity, times the signal strength μ;

• nuisances (48, all on the signal): 45 jet / b-tag shape systematics read as SR_<name>__1up/down, two flat yield systematics (luminosity ±2%, photon-ID ±4%), and the spurious-signal term, see SystematicErrors and _load_signal_systematics.

Limits use the asymptotic CLs method with the q̃μ test statistic (get_T_qmutilde, CLs_scan): the expected bands come from the μ′=0 Asimov data set, the observed limit from the real-data test statistic against the same Asimov-calibrated distributions. μ limits are converted to g_q with gq_from_mu (signal generated at g_q = 0.1, cross-section ∝ g_q², statistics-only extrapolation to the full 2023 luminosity for the partial slices).

Caching

Every stage serialises its result to <plot-output-dir>/.flowchart_cache/<stage>_SRidx<idx>_<tag>.jls, where the tag digests the model file name, sr_idx, the likelihood binning and (for mass-dependent stages) the mass subset. On a rerun a cached stage is loaded and skipped; --redo recomputes only the stages you name, while their dependencies still come from the cache.

What the cache does not know

• The tag does not include the systematics content. After changing the nuisance list use --redo or a fresh output directory, or the run silently reuses pre-change results.

• A --masses <subset> run rewrites the shared summary plots and CSVs (ss_summary, spurious_SRidx<i>.csv, …) even though its checkpoint is separately keyed; a later full run reloads from cache and never restores them. Use a separate output directory for subsets.

• Checkpoints written before September 2026 use a Julia-version-dependent tag; they are adopted automatically when the same Julia version that wrote them runs the directory once.

Reading the results back

• spurious_SRidx<i>.csv: per-mass mean/std of the fitted μ; the test passes when |mean/std| < 0.5.

• pulls_SRidx<i>_minos.csv and pulls_summary_SRidx<i>_minos.csv: per-nuisance pulls, MINOS errors and over-constraint flags (average error < 0.6).

• complexity_pval_ATLAS_poisson_C1.csv: toy-ensemble χ²/p-value per function.

• Limits have no CSV. Deserialize the limits checkpoint:

using BjetTLA, Serialization, Glob
ck  = only(glob("limits_SRidx2_*.jls", "plots/my_run/.flowchart_cache"))
lim = deserialize(ck)                      # (; exp_lim, exp_lim_statonly, obs_lim, CLs_scans, llr_scans, masses, fitstats)
cfg = BjetTLA.FlowchartConfig(; lumi_pct=4, plot_out_dir="plots/my_run")
gq_median = BjetTLA.gq_from_mu(cfg, lim.exp_lim[3])   # bands: [−2σ, −1σ, median, +1σ, +2σ]; NaN = no crossing

Re-fitting the background model

--refit --refit-iters 400 runs the symbolic-regression search on the downscaled C1 (uniform 1 GeV binning) and overwrites the configured model file (export_sr_model); the standalone scripts/fit_and_export_sr.jl --region C1 does the same to a file of your choice. The search is stochastic: judge a new draw by the residuals, chi2pval and spurious stages before adopting it.

Using the stages from Julia

using BjetTLA
cfg      = BjetTLA.FlowchartConfig(; lumi_pct=4, plot_out_dir="plots/my_run")
inp      = BjetTLA.load_flowchart_inputs(cfg)
sr_model = BjetTLA.load_sr_model(cfg)          # includes generated code: advances the world age

Base.invokelatest() do                          # so every call below sees the new functions
    BjetTLA.run_residual_plots(cfg, inp, sr_model; idxs=1:2:9)
    SS_means, SS_stds = BjetTLA.run_spurious_test(cfg, inp, sr_model; N=400)
    res = BjetTLA.run_limit_scan(cfg, inp, sr_model, SS_means; compute_fitstats=true)
    BjetTLA.make_pull_plots(cfg, res.fitstats; masses=res.masses)
end

load_sr_model includes the model file, so call it at top level and enter the consuming code through Base.invokelatest, exactly as run_flowchart.jl does.

API

Configuration and inputs

BjetTLA.FlowchartConfig Type

FlowchartConfig(; lumi_pct, plot_out_dir=nothing, base_path=nothing,
                  sr_lib_dir=<data/sr_lib>, sr_jl_file=nothing, sr_idx=nothing,
                  fit_bin_edges=nothing)

Build the run configuration for a given luminosity slice. lumi_pct must be 4, 20 or 100; the per-slice defaults (data file, histogram dir, SR model, sr_idx) mirror the notebook. Any of plot_out_dir, base_path, sr_jl_file, sr_idx, fit_bin_edges can be overridden. base_path=nothing picks the slice's own histogram directory (data/flowchart_histograms for 4%/20%, data/100percent_histograms for 100%).

fit_bin_edges is the likelihood (resolution) binning used for every fit and plot. The default is the mass-resolution binning _DEFAULT_FIT_BINEDGES = [90, …, 611] (~12–30 GeV bins, widening with mⱼⱼ), comfortably inside the 80→620 SR-fit range. It may be non-uniform.

The litehf S+B fits use the "fine + rebinby" scheme: data and signal/systematic templates are histogrammed on a uniform 1 GeV grid spanning this range (_fine_grid(cfg)), the SR background PDF is evaluated there, and litehf_fit's rebinby sums the fine predictions into fit_BINEDGES. So the background yield is a proper fine-grained integral with no per-bin width factor (the 1 GeV grid is uniform), and the MC-noisy 1 GeV signal templates are summed to the resolution bins by the same rebin, keeping the pulls well-conditioned. The residual/χ²/syst stages instead fit/display directly on fit_BINEDGES (stable, well-populated), which is why refitted_smooth_hist_sr / BinnedDistributionFit are width-aware. The symbolic-regression background fit is done separately, on sr_binedges = uniform 1 GeV over 80→620 (so the SR model is a smooth PDF, consistent with the 1 GeV _fine_grid evaluation above); the shipped *_80to620GeV.jl files were produced this way in notebooks/5-flow_chart.ipynb. See the –refit path in run_flowchart.jl.

source

BjetTLA.load_flowchart_inputs Function

load_flowchart_inputs(cfg; region="C1") -> NamedTuple

Load the data SR + control-region histograms, their rebinned / downscaled variants, the per-mass signal nominal templates (all_sig_Nss), and the Ns bundle used by the spurious / injection tests. Mirrors notebook cells "Load Histograms", "SR Search for candidates", and the "SS test" setup.

source

BjetTLA.load_sr_model Function

load_sr_model(cfg) -> sr_tuple_list

Include the precomputed SR model file into Main (it defines module SymRegModel, whose functions resolve as Main.SymRegModel.f_c*) and return its sr_tuple_list.

source

BjetTLA.export_sr_model Function

export_sr_model(target_downscaled, out_jl; niterations=400)

Fit a symbolic-regression background model to target_downscaled and write the functional forms to out_jl as a SymRegModel module (the LOAD_SR=false branch of the notebook). Returns out_jl.

source

BjetTLA.calculate_res_bins Function

calculate_res_bins(; max_mjj=500, min_mjj=50)

Mass-resolution bin edges over [min_mjj, max_mjj]. The binning is always seeded at min_mjj (the lowest edge sits exactly at min_mjj), then grows upward in steps of the local dijet mass resolution mjj_res(m) = m·sqrt(a/m + (b/m)^2 + c) until the last edge passes max_mjj. Returns integer edges. (The seed used to be a separate argument; it is now fixed to min_mjj by design.)

source

BjetTLA._load_signal_systematics Function

Load a signal's nominal + per-systematic up/down templates on the fit binning.

source

BjetTLA._limit_fit_bounds Function

Bounds for the limit S+B fit: μ≥0, |θ|<6 for each (jet + spurious) nuisance, NB≥0, shape ps free. Shared by run_limit_scan and run_llr_scan so the two fit the identical model.

source

Stages

BjetTLA.run_residual_plots Function

run_residual_plots(cfg, inp, sr_model; idxs=1:5) -> Figure

Stage residuals: refit each SR-ladder rung in idxs to the control region downscaled to the SR yield (on cfg.fit_BINEDGES) and overlay the per-bin pulls with χ²/ndof and p-value per rung. Writes residual_<region>.svg.

source

BjetTLA.run_chi2_pval_plot Function

run_chi2_pval_plot(cfg, inp, sr_model; idx_range=1:5) -> Figure

Stage chi2pval: χ² and p-value of every SR-ladder rung in idx_range refit to the downscaled control region, against its complexity (chi2s_on_sr). Writes complexity_pval_<region>.svg.

source

BjetTLA.run_chi2_pval_atlas_plot Function

run_chi2_pval_atlas_plot(cfg, inp, sr_model; idx_range=1:9)

Like run_chi2_pval_plot, but overlays the standard ATLAS dijet background functions ATLAS_f_4p_xnorm and ATLAS_f_5p_xnorm as ⭐ χ² reference points (on the χ² axis, at the fitted χ²). Each is fit to the same downscaled-CR SR template (identical _refit_smooth procedure as the SR ladder) and placed at its SR-template complexity — nodes(B)+nodes(E) over the {+,−,×} basis: 4p → 8, 5p → 14 — so you can read off where the ATLAS forms' χ² lands versus the symbolic-regression candidates. Writes complexity_pval_ATLAS_<region>.svg.

source

BjetTLA.run_chi2_pval_atlas_poisson_plot Function

run_chi2_pval_atlas_poisson_plot(cfg, inp, sr_model; idx_range=1:9, ntoys=100)

Ensemble version of run_chi2_pval_atlas_plot. Instead of the single deterministic downscale_to(C1, SR) template, draw ntoys Poisson toys of the control region downsampled to the SR event count (resample_poisson(C1, integral(SR)) — per-bin Poisson around the SR-normalized C1 shape, so the total yield fluctuates like a real counting experiment) and fit every SR-ladder rung plus both ATLAS dijet forms to each toy. The plotted χ²/p-value is the ensemble mean, the error bar its standard deviation.

Why: the deterministic template has C1's (≈14× smaller at 4%) statistical error but is judged with SR-sized Neyman errors √obs, so its χ² is one particular realisation with no notion of spread — a 0.5σ shuffle of the points can reorder the ladder or move ATLAS 4p above/below a neighbouring rung. The toy ensemble gives each point the sampling error it would have at SR statistics, which is what you need to say whether two functions genuinely differ. Note this measures the SR-statistics spread only: all toys share the same underlying C1 shape, so a systematic C1→SR shape difference is common to the whole ensemble and does not enter the error bars.

Writes complexity_pval_ATLAS_poisson_<region>.svg + complexity_pval_ATLAS_poisson_<region>.csv via draw_chi2_pval_atlas_poisson, and returns the per-toy χ² matrices, so the same figure and CSV can be reproduced later from the checkpoint without re-fitting anything.

source

BjetTLA.run_cr_divergence_plot Function

run_cr_divergence_plot(cfg)

Shape-divergence diagnostic: compare the area-normalized mⱼⱼ shape of each control region (C1, C0, D1, D0) to the SR, on the fit binning, annotating each with KL/JS/Hellinger to the SR (the metric once used to pick which CR models the SR background). One mass-independent figure, cr_shape_divergence.svg. Loads the regions straight from cfg.data_file — no SR fit needed.

source

BjetTLA.run_bonly_fit Function

run_bonly_fit(cfg, inp, sr_model)

Background-only fit to the SR data: the SR-function background (free NB + shape ps) is fit to the signal-region data with the signal strength fixed to μ=0 and no systematics — i.e. the post-fit plot without any signal, so it is a single mass-independent figure bonly_fit_SRidx<idx>.svg. Uses the same "fine (1 GeV) + rebinby" evaluation as the S+B limit fits, so it is directly comparable.

source

BjetTLA.run_spurious_test Function

run_spurious_test(cfg, inp, sr_model; N=250, debug_mask=1:10, poisson=true) -> (SS_means, SS_stds)

Per signal mass, resample background-only toys, fit S+B, and summarize the fitted-μ bias. Saves a per-mass plot_resample_SS and the plot_SS_scan summary. Returns the mean/std of fitted μ per mass; SS_means feeds the spurious nuisance in the limit fit.

source

BjetTLA.run_injection_test Function

run_injection_test(cfg, inp, sr_model; factors=[1,4,9,16,25], N=200, ntasks=5) -> SI_results

Inject signal at the given strength factors and recover μ, testing response linearity. Saves the plot_SI_scan figure and returns the raw scan results.

source

BjetTLA.run_limit_scan Function

run_limit_scan(cfg, inp, sr_model, SS_means; adaptive_mu=true, nσ_max=5, n_mu=20,
               mu_range=15:0.5:20, debug_mask=1:10, sys_names=jet_syst_names,
               maxiters=60_000, cls_maxiters=150_000) -> NamedTuple

For each mass: build the per-mass systematics (jet variations + spurious nuisance from SS_means), fit S+B to the real SR data, and run the CLs scan. With adaptive_mu=true (default) the scan range is seeded per mass from that fit's σ_μ as 0 → nσ_max·σ_μinn_musteps, so it brackets the ±2σ expected band at uniform resolution for every mass; setadaptive_mu=falseto use the fixedmu_rangefor all masses.maxiterscaps the main S+B fit;cls_maxiterscaps the CLs-scan profile/Asimov fits (threaded down tolitehf_freeze_fit, where the convergence warnings originate).

With compute_statonly=true (default) each mass ALSO gets a second, systematics-free fit + CLs scan (no jet nuisances, no spurious-signal term), yielding _statonly companions of both limit plots. With compute_llr=true each mass ALSO gets a profile-likelihood-ratio scan q(μ)=−2lnλ(μ) on the real-data S+B fit over the same adaptive μ grid (feeds the :llrscan per-mass plots); these are nothing when compute_llr=false. Saves limit_plot[_mu][_statonly]_SRidx<idx>.svg and returns (; sols_sys, CLs_scans, llr_scans, exp_lim, exp_lim_statonly, masses, fig, fitstats).

source

BjetTLA.run_llr_scan Function

run_llr_scan(cfg, inp, sr_model, SS_means; debug_mask=1:10, maxiters=60_000,
             cls_maxiters=150_000, adaptive_mu=true, nσ_max=5, n_mu=20, progress_offset=0)

Standalone per-mass profile-likelihood-ratio scan: build the per-mass systematics, fit S+B to the real SR data (identical model + bounds as run_limit_scan), and scan q(μ)=−2lnλ(μ) over the adaptive μ grid — WITHOUT the CLs scan or MINOS. This is the :llrscan fallback used when the cached limits checkpoint predates compute_llr; it computes only what the llrscan plots need and never rewrites the limits checkpoint. Returns the vector of (μ grid, −2lnλ, μ̂) tuples (one per mass, nothing on a non-convergent fit).

source

BjetTLA.make_llrscan_plots Function

make_llrscan_plots(cfg, llr_scans; masses, progress_offset=0)

Draw one profile-likelihood-ratio scan per mass (llrscan_<mass>GeV_SRidx<idx>.svg) from the (μ grid, −2lnλ, μ̂) tuples produced by run_limit_scan(...; compute_llr=true). A nothing entry (scan not computed or failed to converge) is skipped with a warning.

source

BjetTLA.make_pull_plots Function

make_pull_plots(cfg, fitstats; masses=mass_points, tag="_minos")

Per-mass nuisance pull plots and the matching post-fit data/prediction plots. fitstats is a vector of FitStats (produced by run_limit_scan(...; compute_fitstats=true), which can be cached and reloaded). Pass a vector of sols instead to build the FitStats on the fly (use_minos controls how). Writes pulls_<mass>GeV_*.svg and postfit_<mass>GeV_*.svg.

source

BjetTLA.run_impacts Function

run_impacts(cfg, inp, sr_model, SS_means; debug_mask=1:10, maxiters=60_000,
            freeze_maxiters=150_000, use_minos=true, progress_offset=0)

CMS-style impact/ranking plot per mass. For each mass: fit S+B to the real SR data (identical model

• bounds as run_limit_scan, incl. jet systematics + spurious nuisance), then for every nuisance θ_i

compute the POST-fit impact on the POI Δµ_i^± = µ̂(θ_i = θ̂_i ± σ_i) − µ̂ by freezing θ_i at its post-fit ±1σ (MINOS) and re-profiling (litehf_freeze_theta_fit). Writesimpacts_<mass>GeV_SRidx<idx>.svg (pull + constraint on one x-axis, Δµ bars ranked by |Δµ| on the other). Heavy: 2×N_nuisance refits per mass on top of the FitStats/MINOS build.

source

BjetTLA.dump_syst_templates Function

dump_syst_templates(cfg; masses=mass_points, sys_names=jet_syst_names) -> Dict

For every signal mass, plot nominal vs ±1σ for every systematic and save under <plot_out_path>/syst_templates/. Returns the loaded histograms keyed by mass.

source

Likelihood and test statistics

BjetTLA.SystematicErrors Type

SystematicErrors(sig_var_names, sig_var_hists_hi, sig_var_hists_lo, mu_spur=0.0)
SystematicErrors(mu_spur)

The nuisance-parameter content of a litehf_fit. Every entry of sig_var_names becomes a Gaussian-constrained shape systematic on the signal template, interpolated piecewise-linearly between the _hi/_lo histograms (which must share the signal's binning). mu_spur is the spurious-signal mean from run_spurious_test; when non-zero it adds one more constrained nuisance ss_theta whose ±1σ template is ±mu_spur × signal. The background is data-driven, so it carries no systematics here. SystematicErrors(0) means "no systematics at all".

source

BjetTLA.litehf_fit Function

litehf_fit(hist_data, func_bkg, p0, hist_sig; systematics=nothing, rebinby=ResBin_4p,
           lb=nothing, ub=nothing, maxiters=10^4, b_only=false) -> NamedTuple

Binned signal-plus-background maximum-likelihood fit (Migrad, via LiteHF), the likelihood at the heart of the spurious-signal, injection, limit and pull stages.

Model, evaluated on the fine binning of hist_data (1 GeV in the flowchart) and then summed into the likelihood bins rebinby (an integer factor or a vector of edges) before the Poisson term:

expected = μ · S(θ) + ss_theta · mu_spur · S + NB · f(x; ps) / ∫f

• func_bkg(x, ps) is the symbolic-regression background shape (an SR-ladder f_c<N>), p0 its starting parameters; NB is a free yield, so the background enters as a normalised PDF.

• hist_sig is the nominal signal template; systematics::SystematicErrors adds the Gaussian-constrained shape nuisances sys_thetas and the spurious-signal nuisance ss_theta.

• lb/ub bound the flat parameter vector [μ, sys_thetas…, ss_theta…, NB, ps…]; by default only NB ≥ 0. The limit fits use _limit_fit_bounds (μ ≥ 0, |θ| < 6).

• b_only=true drops the signal block entirely (see litehf_bonly_fit).

Returns (; objective, original, u, hist_data, hist_pdf1, hist_pdf2, hist_pdf_bkg, hist_pdf_spur, NLL, fit_args, EB): the NLL at the minimum, the Minuit object, the best-fit ComponentArray u (fields mu, sys_thetas, ss_theta, NB, ps), the post-fit background (hist_pdf1, including the spurious component), signal (hist_pdf2) and pure-PDF (hist_pdf_bkg) predictions on the fine binning, plus the NLL closure and the inputs (fit_args) needed by litehf_freeze_fit, generate_asimov and FitStats.

source

BjetTLA.litehf_bonly_fit Function

litehf_bonly_fit(hist_data, func_bkg, p0; kwargs...)

Background-only litehf_fit: only NB and the shape parameters ps float (no μ, no nuisances). kwargs are forwarded (rebinby, lb, ub, maxiters).

source

BjetTLA.litehf_freeze_fit Function

litehf_freeze_fit(nominal_fit, mu; constrain_p1=false, maxiters=150_000, strategy=2)

Conditional (profile) fit: re-minimise the NLL of a litehf_fit result with the signal strength frozen at mu, all nuisances, NB and the shape parameters floating from the nominal best fit and inside the nominal bounds. Returns the Optimization solution (.objective is the conditional NLL, .u the profiled parameters). This is the building block of get_lnLRtilde, generate_asimov and every CLs / LLR scan point.

source

BjetTLA.litehf_freeze_theta_fit Function

litehf_freeze_theta_fit(nominal_fit, theta_idx, theta_val; maxiters=150_000, strategy=2) -> μ̂

Refit with ONE nuisance frozen at theta_val (µ and all other parameters float) and return the resulting best-fit µ. The mirror of litehf_freeze_fit, which freezes µ. theta_idx indexes the combined nuisance vector [sys_thetas; ss_theta] (1-based). Building block of the CMS-style impacts in run_impacts: Δµ_i^± = µ̂(θ_i = θ̂_i ± σ_i) − µ̂.

source

BjetTLA.get_lnLRtilde Function

get_lnLRtilde(nominal_fit; maxiters=150_000) -> lnLRtilde(μ)

Return the log profile-likelihood ratio ln λ̃(μ) = ln L(μ, θ̂̂(μ)) − ln L(μ̂, θ̂) of a litehf_fit result (Cowan et al., arXiv:1007.1727 eq. 10): the denominator is the unconditional maximum, replaced by the μ=0 conditional maximum when μ̂ < 0. Each call costs one litehf_freeze_fit.

source

BjetTLA.get_T_qmutilde Function

get_T_qmutilde(nominal_fit; maxiters=150_000) -> q̃(μ)

Return the upper-limit test statistic q̃_μ = −2 ln λ̃(μ) for μ ≥ μ̂ and 0 otherwise (arXiv:1007.1727 eq. 16), built on get_lnLRtilde. Evaluate it on the real-data fit for the observed limit and on the μ'=0 Asimov fit for the expected bands, see CLs_expected_and_observed.

source

BjetTLA.generate_asimov Function

generate_asimov(mu, nominal_fit; maxiters=150_000) -> Hist1D

Asimov data set for signal strength mu: profile the nuisances at mu with litehf_freeze_fit and return the expected counts of the fitted model on the fine binning of the original data (no fluctuations). generate_asimov(0, fit) is the background-only Asimov used to calibrate the expected CLs bands in CLs_scan.

source

BjetTLA.CLs_expected_and_observed Function

CLs_expected_and_observed(qA_f, qO_f, mu) -> (bands::Vector, obs::Float64)

The five EXPECTED CLs bands and the OBSERVED CLs at mu, from one shared Asimov test-statistic evaluation. qA_f / qO_f are pre-built get_T_qmutilde closures for the Asimov and the REAL data respectively (hoisted by the caller so the closure setup isn't repeated per µ).

Expected bands come from expected_pvalue on the µ'=0 Asimov, exactly as before.

The observed leg uses the real-data test statistic q_obs(µ) against the SAME Asimov-calibrated distributions. IMPORTANT: LiteHF's pvalue(d, v) = Φ(-(v - d.shift)) takes v in the SHIFTED variable that expected_pvaluefeeds it (where the median corresponds to v = 0), NOT √q itself. For the q̃_µ test statistic the observed value enters as (Cowan et al. arXiv:1007.1727 eqs. 65-66, identical to pyhf'sAsymptoticCalculator.teststatisticforqtilde)

v = √q_obs − √q_A                     if √q_obs ≤ √q_A   (q̃_µ ≤ µ²/σ²)
v = (q_obs − q_A) / (2 √q_A)          otherwise          (q̃_µ >  µ²/σ², data below background)

The first branch is fixed by two limiting requirements, both asserted in the test suite:

• q_obs == q_A (data == background-only Asimov) -> CLs_obs == the MEDIAN expected band.

• q_obs == 0 (µ <= µ̂, where q̃ is defined as 0) -> CLs_obs ≈ 0.5, i.e. NO exclusion.

Passing √q_obs directly instead silently collapses case 2 onto the median band, making the observed limit exactly equal the expected median wherever µ̂ exceeds the crossing. The second branch matters whenever the data fluctuate BELOW the background (observed limit stronger than the median): there q̃_µ follows the non-central-χ² tail of eq. 65, and using the first branch for it over-estimates the observed CLs (weaker exclusion) by an amount growing with (q_obs − q_A).

Pass qO_f === nothing to skip the observed leg (returns NaN for it).

Non-finite test statistics are mapped to CLs = 0: a µ whose freeze-fit fails at high µ is physically maximally excluded, and a stray NaN would otherwise abort the crossing search.

source

BjetTLA.CLs_scan Function

CLs_scan(sol, range; search_crossing=true, maxiters=150_000, max_extend=400, observed=true)
    -> (exp_lims, mus, cls_bands, obs_lim, cls_obs)

Asymptotic CLs scan for one mass point from a litehf_fit result sol on the real data. Builds the μ'=0 Asimov data set (generate_asimov), refits it with the same model, and at every μ of range evaluates the five expected CLs bands [−2σ, −1σ, median, +1σ, +2σ] and the observed CLs (CLs_expected_and_observed). With search_crossing=true the grid is extended upward (by the step of range, at most max_extend points) until both the +2σ band and the observed curve have dropped below 0.05, then each curve's 0.05 crossing is interpolated (find_zero_crossing).

Returns the five expected 95% CL limits on μ (nothing where a band never crosses), the μ grid, the per-μ band vectors, the observed limit (nothing when observed=false) and the per-μ observed CLs. With search_crossing=false only cls_bands is returned. Convert μ limits to g_q with gq_from_mu.

source

BjetTLA.find_zero_crossing Function

find_zero_crossing(x, y) -> Union{Float64, Nothing}

Linearly interpolated position of the first sign change of y along x; nothing if there is none. Pairs containing a NaN never count as a crossing. Used to read the 95% CL limit off a CLs scan (find_zero_crossing(mus, cls .- 0.05)).

source

BjetTLA.gq_stat_scale Function

gq_stat_scale(cfg) -> Float64

Statistics-only luminosity extrapolation factor √(L_2023 / L_slice) by which a μ limit measured on the cfg.lumi_pct slice is divided to estimate the full-2023 limit (limit ∝ 1/√L). Exactly 1 for the 100% slice.

source

BjetTLA.gq_from_mu Function

gq_from_mu(cfg, mu) -> g_q

Convert a 95% CL limit on the signal strength μ into a limit on the Z′ coupling g_q: the signal is generated at g_q = 0.1 and σ ∝ g_q², so g_q = 0.1 √(μ / gq_stat_scale(cfg))after the stat-only luminosity extrapolation of gq_stat_scale. Broadcasts over vectors; NaN in means NaN out (no CLs crossing). Mirrors the notebook'sto_gq.

source

Toy ensembles

BjetTLA.resample_poisson Function

resample_poisson(h, target_integral) -> Hist1D

Background-only toy: renormalise the template h (a control-region shape) to target_integral events and fluctuate every bin independently with a Poisson distribution, so the total yield fluctuates too (the physically correct model for a counting experiment; resample instead draws a fixed number of events and under-estimates the toy variance).

source

BjetTLA.resample Function

resample(h, n) -> Hist1D

Draw exactly n events from the shape of h (multinomial, total pinned to n) and histogram them on the same binning. Legacy toy model; prefer resample_poisson.

source

BjetTLA.resample_SS Function

resample_SS(sig_hist, data_hists, sr_model; factor=0.0, region=:N_C1, sr_idx=2, N=200,
            is_machine=false, maxiters=10_000, poisson=true, rebinby=1)
    -> (fitinfo::Vector, oneshot_fitinfo)

Spurious-signal / injection toy ensemble. Draw N background-only toys of the control region getproperty(data_hists, region) normalised to the signal-region yield integral(data_hists.N_A) (Poisson per bin by default, see resample_poisson), inject factor × sig_hist into each and fit them with lf_signal_injection. Also returns the single "oneshot" fit to the control region deterministically downscaled to the SR yield. The spread of fi.u.mu over fitinfo is the spurious-signal (factor = 0) or recovered-strength (factor > 0) distribution.

source

BjetTLA.lf_signal_injection Function

lf_signal_injection(sig_hist, bkg_hist, sr_model; sr_idx, factor=0.0, rebinby=1,
                    is_machine=false, systematics=nothing, inj_hist=nothing, maxiters=10_000)

Fit bkg_hist + factor × (inj_hist or sig_hist) with litehf_fit, using rung sr_idx of sr_model (an exported SR tuple list, or an MLJ machine when is_machine=true) as the background shape. factor=0 is a pure background-only toy fit (spurious-signal test); factor>0 injects signal (linearity test). Returns the litehf_fit result.

source

BjetTLA.lf_scan Function

lf_scan(sig, target, sr_model; sr_idx, factors, N=100, region=:N_C1, is_machine=false, rebinby=1)
    -> (factors, mean_mu, std_mu)

Signal-injection linearity scan: for every injected strength in factors, run N toy fits via resample_SS and return the mean and standard deviation of the fitted μ (non-finite toy fits are dropped with a warning). A linear response has mean_mu ≈ factors.

source

Background-model goodness of fit

BjetTLA._refit_smooth Function

_refit_smooth(hist, sr_f, sr_p0) -> (; hist, p)

Extended-Poisson fit of sr_f(x, p) (shape starting at sr_p0) to hist, with a DECOUPLED linear normalization N (the total yield, started at the data integral) carrying the scale and sr_f supplying the shape: pred_i = N · ŝ_i, where ŝ is sr_f(cᵢ,p)·wᵢ normalized to unit sum. The shape normalization is max-stabilized (divide by the peak before summing) so that an overall scale factor which underflows Float64 — e.g. the |p2|^p1 prefactor of a constant-base SR form — cancels numerically instead of giving 0/0 = NaN. That underflow (with the old ExtendPdf+norms path) is what sent constant-base rungs (c2/c10 on 4%) to Inf χ².

Keeping N decoupled and well-started is what makes the fit CONVERGE: the SR forms bury their scale nonlinearly in the shape params (no clean linear knob), so a norm-free "direct" fit cannot recover the ~14× downscale factor and every rung mis-fits (c4 χ² 1.2 → 1100+). So the normalization d.o.f. is needed here too, not only in the litehf S+B hybrid fit — just made robust. For a non-degenerate form the MLE is identical to the old ExtendPdf result (same N·ŝ model), so the residual and χ²-p-value plots are unchanged except that the previously-broken rungs now render finite. Shared by refitted_smooth_hist_sr (residual/χ² ladder) and chi2_pval_of_func (external refs).

source

BjetTLA.refitted_smooth_hist_sr Function

refitted_smooth_hist_sr(hist, sr_model; idx=5, is_machine=true) -> (; sol=(; complexity), hist)

Refit rung idx of the SR ladder to hist with _refit_smooth and return the fitted prediction as a histogram (plus the rung's complexity). sr_model is either an exported sr_tuple_list (is_machine=false, the flowchart case) or the MLJ machine of fit_sr_model.

source

BjetTLA.chi2s_on_sr Function

chi2s_on_sr(sr_model, hist, idx_range; rebinby, is_machine=true) -> Vector{Float64}

Neyman χ² (σ = √observed) of every SR-ladder rung in idx_range refit to hist rebinned by rebinby — the goodness-of-fit ladder behind the χ²/p-value-vs-complexity plot (run_chi2_pval_plot). The number of free parameters per rung is 1 + length(p0).

source

BjetTLA.chi2_pval_of_func Function

chi2_pval_of_func(hist, func, p0; rebinby=1, nparas=length(p0)) -> (; chi2, pval, ndof)

Fit an arbitrary background func(x, p) to hist with _refit_smooth and return its goodness of fit using the SAME Neyman χ² (σ = √observed) and ndof convention as chi2s_on_sr / the χ²-p-value plot. nparas is the free-parameter count (default length(p0) — for the ATLAS *_xnorm forms the leading parameter is the overall norm, so this equals 1 + #shape params, matching the SR ladder's 1 + length(p0)). Used to overlay external reference functions on the complexity/p-value plot.

source

BjetTLA.ATLAS_f_4p_xnorm Function

ATLAS_f_4p_xnorm(x, p)

Standard ATLAS 4-parameter dijet function in GeV, x normalised internally by √s = 13600: p₁ (1−x)^{p₂} x^{p₃ + p₄ ln x}. ATLAS_f_3p_xnorm, ATLAS_f_5p_xnorm and ATLAS_f_6p_xnorm drop / add the ln² x, ln³ x terms. Used as the ⭐ reference points in run_chi2_pval_atlas_plot. The un-suffixed ATLAS_f_4p etc. expect x already divided by √s.

source

Fit summaries and plots

BjetTLA.FitStats Type

FitStats

Serializable summary of one litehf_fit, on the likelihood (resolution) binning: the data, post-fit background (incl. spurious component) and signal histograms; mu_fit with its symmetric / MINOS _lo/_hi errors; the nuisance pulls thetas with errors and sys_names (jet systematics, luminosity, photon-ID, then "Spurious signal"); the bin edges bes; the SR complexity; and ndf, χ2ndf, p_val of the post-fit χ² (Neyman, empty bins floored at 1). Constrained nuisances do not count as free parameters in ndf. Built with FitStats(fitinfo; use_minos=true, minos_cl=1, b_only=false); b_only=true blanks the signal and nuisance fields (they become nothing). Consumed by draw_postfit, draw_pulls, draw_bonly, draw_impacts and the pulls CSVs of make_pull_plots.

source

BjetTLA.draw_postfit Function

draw_postfit(fs::FitStats; lumi=930, mass="", signal_label="Signal", save_path=nothing,
             ymin_r=0.5, ymax_r=1.5) -> Figure

Post-fit plot of an S+B fit: data (per GeV) with the fitted background and background+signal overlaid, plus a Data/Prediction panel whose error bars carry the data's Poisson error only. The ATLAS label quotes μ_fit with its asymmetric errors and χ²/ndf. lumi in pb⁻¹.

source

BjetTLA.draw_bonly Function

draw_bonly(fs::FitStats; lumi=930, save_path=nothing, ymin_r=0.5, ymax_r=1.5)

Background-only fit plot: SR data + fitted background PDF + a data/prediction ratio panel. Like draw_postfit but with NO signal component (the fit is μ-fixed-at-0), so it is a single mass-independent figure. Expects a FitStats built from a background-only litehf_fit (FitStats(fitinfo; b_only=true)), whose hist_sig is ~0.

source

BjetTLA.draw_pulls Function

draw_pulls(fs::FitStats; save_path=nothing, mass="") -> Figure

Nuisance-parameter pull plot: one row per entry of fs.sys_names showing (θ̂ − θ₀)/Δθ with its post-fit (MINOS) error over the ±1σ / ±2σ prior bands. An error bar visibly shorter than 1 means the fit constrains that nuisance beyond its prior ("over-constrained").

source

BjetTLA.draw_impacts Function

draw_impacts(fs::FitStats; impacts_hi, impacts_lo, mass="", save_path=nothing)

CMS-style impact/ranking plot: nuisances ranked by their impact on µ (largest at top). Each row shows the PULL (θ̂−θ₀)/Δθ with its post-fit error (black marker, over the ±1σ/±2σ prior bands, read on the bottom axis) and the IMPACT Δµ^± bars (blue = +1σ, red = −1σ, read on the top axis). impacts_hi/impacts_lo are the Δµ from freezing each nuisance at θ̂±σ (see run_impacts); they are drawn rescaled onto the pull axis so bands→bars→marker layer correctly, with a twin top axis carrying the true Δµ scale.

source

BjetTLA.plot_gq_limit Function

plot_gq_limit(cfg, masses, exp_lim; stat_only=false) -> Figure

Brazilian-band expected exclusion in g_q vs m_Z′ from the five per-band μ limits exp_lim ([−2σ, −1σ, median, +1σ, +2σ], each a vector over masses), converted with gq_from_mu. stat_only only changes the annotation.

source

BjetTLA.plot_mu_limit Function

Brazilian exclusion plot in raw signal strength μ (no luminosity rescaling, no μ -> g_q). stat_only only changes the annotation.

source