#!/usr/bin/env python3
"""Reproducible public-data companion tests for END Theory 2.0.

This program tests the experimental analysis contracts in dossiers D02-D06.
It does not infer that END Theory is validated because the manuscript does not
freeze experiment-specific END coefficients.  Outputs are deterministic for
the declared random seed and are written as CSV, JSON, and PNG files.
"""

from __future__ import annotations

import argparse
import gzip
import hashlib
import json
import math
import pickle
import platform
import sys
from dataclasses import dataclass
from pathlib import Path

import h5py
import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy
import uproot
from scipy.linalg import expm
from scipy.optimize import brentq, minimize, minimize_scalar
from scipy.signal import butter, correlate, correlation_lags, sosfiltfilt, welch
from scipy.special import gammaln
from scipy.stats import chi2, norm


SEED = 20260808
INK = "#152238"
BLUE = "#2E74B5"
TEAL = "#1B7F79"
GOLD = "#B57A13"
RED = "#A33A3A"
GRAY = "#64748B"
LIGHT = "#E8EEF5"


def sha256(path: Path) -> str:
    h = hashlib.sha256()
    with path.open("rb") as stream:
        for block in iter(lambda: stream.read(1024 * 1024), b""):
            h.update(block)
    return h.hexdigest()


def write_csv(path: Path, rows: list[dict]) -> None:
    pd.DataFrame(rows).to_csv(path, index=False)


def json_ready(value):
    if isinstance(value, dict):
        return {str(k): json_ready(v) for k, v in value.items()}
    if isinstance(value, (list, tuple)):
        return [json_ready(v) for v in value]
    if isinstance(value, (np.integer,)):
        return int(value)
    if isinstance(value, (np.floating,)):
        return float(value)
    if isinstance(value, (np.bool_,)):
        return bool(value)
    if isinstance(value, np.ndarray):
        return value.tolist()
    return value


def configure_plotting() -> None:
    plt.rcParams.update(
        {
            "figure.dpi": 180,
            "savefig.dpi": 220,
            "font.family": "DejaVu Sans",
            "font.size": 9,
            "axes.labelcolor": INK,
            "axes.edgecolor": "#94A3B8",
            "axes.titlecolor": INK,
            "xtick.color": INK,
            "ytick.color": INK,
            "grid.color": "#CBD5E1",
            "grid.alpha": 0.45,
            "legend.frameon": False,
        }
    )


# ---------------------------------------------------------------------------
# GWOSC D02: inter-detector timing and frequency-coherence diagnostic


def whiten_strain(x: np.ndarray, fs: float, event_index: int) -> np.ndarray:
    mask = np.ones(len(x), dtype=bool)
    margin = int(0.5 * fs)
    mask[max(0, event_index - margin) : min(len(x), event_index + margin)] = False
    offsource = x[mask]
    f_psd, psd = welch(
        offsource,
        fs=fs,
        nperseg=int(4 * fs),
        noverlap=int(2 * fs),
        average="median",
    )
    f = np.fft.rfftfreq(len(x), 1.0 / fs)
    psd_i = np.interp(f, f_psd, psd, left=psd[0], right=psd[-1])
    # Absolute normalization cancels in normalized correlation.
    return np.fft.irfft(np.fft.rfft(x) / np.sqrt(psd_i / 2.0 * fs), n=len(x))


def estimate_lag_samples(h: np.ndarray, l: np.ndarray, fs: float) -> tuple[int, np.ndarray, np.ndarray]:
    cross = correlate(h, -l, mode="full", method="fft")
    lag_samples = correlation_lags(len(h), len(l), mode="full")
    keep = np.abs(lag_samples / fs) <= 0.015
    selected = int(lag_samples[keep][np.argmax(cross[keep])])
    denom = math.sqrt(float(np.dot(h, h) * np.dot(l, l)))
    return selected, lag_samples[keep] / fs, cross[keep] / denom


def run_gwosc(data_root: Path, out: Path) -> dict:
    fs = 4096.0
    gps_start = 1126259447.0
    event_gps = 1126259462.4
    event_index = int(round((event_gps - gps_start) * fs))
    h_path = data_root / "gwosc" / "H-H1_GWOSC_4KHZ_R1-1126259447-32.txt.gz"
    l_path = data_root / "gwosc" / "L-L1_GWOSC_4KHZ_R1-1126259447-32.txt.gz"
    with gzip.open(h_path, "rt") as stream:
        h_raw = np.loadtxt(stream)
    with gzip.open(l_path, "rt") as stream:
        l_raw = np.loadtxt(stream)
    if len(h_raw) != 131072 or len(l_raw) != 131072:
        raise ValueError("GWOSC sample count does not match 32 s at 4096 Hz")
    h_white = whiten_strain(h_raw, fs, event_index)
    l_white = whiten_strain(l_raw, fs, event_index)

    duration = 0.12
    nwin = int(round(duration * fs))
    lo = event_index - nwin // 2
    hi = lo + nwin
    bands = [(35, 70), (70, 120), (120, 200), (200, 350)]
    all_bands = [(35, 350)] + bands
    rng = np.random.default_rng(SEED)
    exclude = int(0.75 * fs)
    centers = np.r_[
        np.arange(nwin // 2 + 10, event_index - exclude - nwin // 2),
        np.arange(event_index + exclude + nwin // 2, len(h_white) - nwin // 2 - 10),
    ]

    rows = []
    broad_cross = None
    broad_lags = None
    for low, high in all_bands:
        sos = butter(4, [low, high], btype="bandpass", fs=fs, output="sos")
        h = sosfiltfilt(sos, h_white)
        l = sosfiltfilt(sos, l_white)
        h_event = h[lo:hi]
        l_event = l[lo:hi]
        lag_samples, lag_axis, cross_axis = estimate_lag_samples(h_event, l_event, fs)
        if (low, high) == (35, 350):
            broad_cross = cross_axis
            broad_lags = lag_axis

        # Injection/recovery: average the aligned event strains to form a
        # detector-independent signal proxy, inject it into independent
        # off-source noise windows, and re-estimate the lag.
        signal_proxy = (h_event + np.roll(-l_event, lag_samples)) / 2.0
        recovered = []
        for _ in range(500):
            ch = int(rng.choice(centers))
            cl = int(rng.choice(centers))
            nh = h[ch - nwin // 2 : ch - nwin // 2 + nwin]
            nl = l[cl - nwin // 2 : cl - nwin // 2 + nwin]
            h_inj = nh + signal_proxy
            l_inj = nl - np.roll(signal_proxy, -lag_samples)
            k, _, _ = estimate_lag_samples(h_inj, l_inj, fs)
            recovered.append(1000.0 * k / fs)
        q = np.percentile(recovered, [2.5, 16, 50, 84, 97.5])
        rows.append(
            {
                "f_low_hz": low,
                "f_high_hz": high,
                "f_geometric_hz": math.sqrt(low * high),
                "lag_ms": 1000.0 * lag_samples / fs,
                "bootstrap_p2_5_ms": q[0],
                "bootstrap_p16_ms": q[1],
                "bootstrap_median_ms": q[2],
                "bootstrap_p84_ms": q[3],
                "bootstrap_p97_5_ms": q[4],
            }
        )

    write_csv(out / "gwosc_subband_timing.csv", rows)
    sub = rows[1:]
    x = np.log(np.array([r["f_geometric_hz"] for r in sub]) / 100.0)
    y = np.array([r["lag_ms"] for r in sub])
    sigma = np.array(
        [max(r["lag_ms"] - r["bootstrap_p16_ms"], r["bootstrap_p84_ms"] - r["lag_ms"], 1.0 / fs * 1000) for r in sub]
    )
    design = np.column_stack([np.ones(len(x)), x])
    weight = np.diag(1.0 / sigma**2)
    cov = np.linalg.inv(design.T @ weight @ design)
    coeff = cov @ design.T @ weight @ y
    errors = np.sqrt(np.diag(cov))
    chi_sq = float(np.sum(((y - design @ coeff) / sigma) ** 2))

    official = {"central_ms": 6.9, "minus_ms": 0.4, "plus_ms": 0.5}
    broad = rows[0]
    within_official = (
        official["central_ms"] - official["minus_ms"]
        <= broad["lag_ms"]
        <= official["central_ms"] + official["plus_ms"]
    )

    fig, axes = plt.subplots(1, 2, figsize=(9.2, 3.7))
    ax = axes[0]
    ax.plot(1000 * broad_lags, broad_cross, color=BLUE, lw=1.8)
    ax.axvline(broad["lag_ms"], color=RED, lw=1.2, label=f"Recovered {broad['lag_ms']:.3f} ms")
    ax.axvspan(6.5, 7.4, color=GOLD, alpha=0.18, label="GWOSC published interval")
    ax.set(xlabel="H1 relative lag [ms]", ylabel="Normalized correlation", title="GW150914 H1/L1 timing")
    ax.grid(True)
    ax.legend(fontsize=8)

    ax = axes[1]
    fmid = np.array([r["f_geometric_hz"] for r in sub])
    ax.errorbar(fmid, y, yerr=sigma, fmt="o", color=TEAL, capsize=3, label="Sub-band recovery")
    fplot = np.geomspace(fmid.min() * 0.9, fmid.max() * 1.1, 100)
    ax.plot(fplot, coeff[0] + coeff[1] * np.log(fplot / 100), color=INK, lw=1.3, label="Weighted log-frequency fit")
    ax.set_xscale("log")
    ax.set(
        xlabel="Frequency [Hz]",
        ylabel="H1 relative lag [ms]",
        title="Frequency-coherence diagnostic",
    )
    ax.grid(True, which="both")
    ax.legend(fontsize=8)
    fig.tight_layout()
    fig.savefig(out / "figure_gwosc_timing.png", bbox_inches="tight")
    plt.close(fig)

    return {
        "dataset": "GW150914-v3 4 kHz 32 s H1/L1 strain",
        "samples_per_detector": len(h_raw),
        "broad_band_hz": [35, 350],
        "window_s": duration,
        "recovered_h1_relative_lag_ms": broad["lag_ms"],
        "bootstrap_68_interval_ms": [broad["bootstrap_p16_ms"], broad["bootstrap_p84_ms"]],
        "published_h1_relative_lag_ms": official,
        "within_published_interval": within_official,
        "log_frequency_slope_ms_per_ln_f_over_100hz": coeff[1],
        "log_frequency_slope_standard_error_ms": errors[1],
        "slope_z": coeff[1] / errors[1],
        "fit_chi2": chi_sq,
        "fit_dof": len(y) - 2,
        "inference_scope": "detector-baseline timing diagnostic; not a cosmological A_alpha bound",
        "files": {h_path.name: sha256(h_path), l_path.name: sha256(l_path)},
    }


# ---------------------------------------------------------------------------
# CERN D03 / Eq. 21.4: four-lepton scalar-template likelihood


@dataclass
class CERNFit:
    signal_yield: float
    background_yield: float
    slope_per_gev: float
    nll: float


def cern_probabilities(edges: np.ndarray, slope: float, mass: float = 125.0, sigma: float = 2.5):
    signal = np.diff(norm.cdf(edges, loc=mass, scale=sigma))
    signal = signal / signal.sum()
    a = edges[:-1] - edges[0]
    b = edges[1:] - edges[0]
    if abs(slope) < 1e-9:
        background = b - a
    else:
        background = (np.exp(slope * b) - np.exp(slope * a)) / slope
    background = background / background.sum()
    return signal, background


def fit_cern_counts(counts: np.ndarray, edges: np.ndarray, fixed_signal: float | None = None) -> CERNFit:
    def nll_free(par):
        if fixed_signal is None:
            s, b, slope = par
        else:
            s = fixed_signal
            b, slope = par
        sp, bp = cern_probabilities(edges, slope)
        mu = np.clip(s * sp + b * bp, 1e-12, None)
        return float(np.sum(mu - counts * np.log(mu) + gammaln(counts + 1)))

    if fixed_signal is None:
        initial = [max(1.0, counts.sum() / 3), max(1.0, counts.sum() * 2 / 3), 0.0]
        bounds = [(0, 1000), (1e-6, 1000), (-0.2, 0.2)]
    else:
        initial = [max(1.0, counts.sum() - fixed_signal), 0.0]
        bounds = [(1e-6, 1000), (-0.2, 0.2)]
    result = minimize(nll_free, initial, method="L-BFGS-B", bounds=bounds, options={"ftol": 1e-12, "maxiter": 2000})
    if not result.success:
        raise RuntimeError(f"CERN likelihood fit failed: {result.message}")
    if fixed_signal is None:
        s, b, slope = result.x
    else:
        s = fixed_signal
        b, slope = result.x
    return CERNFit(float(s), float(b), float(slope), float(result.fun))


def extract_four_lepton_masses(files: list[Path]) -> tuple[np.ndarray, dict]:
    stages = {"input_preselected": 0, "trigger": 0, "pt": 0, "tight_id": 0, "flavour": 0, "charge": 0, "isolation": 0}
    masses = []
    branches = [
        "lep_pt",
        "lep_eta",
        "lep_phi",
        "lep_E",
        "lep_charge",
        "lep_type",
        "lep_isTightID",
        "lep_ptcone30",
        "lep_etcone20",
        "lep_trigMatched",
        "trigE",
        "trigM",
    ]
    for path in files:
        tree = uproot.open(path)["mini"]
        data = tree.arrays(branches, library="ak")
        for i in range(len(data)):
            stages["input_preselected"] += 1
            if not (bool(data.trigE[i]) or bool(data.trigM[i])):
                continue
            stages["trigger"] += 1
            if len(data.lep_pt[i]) < 4:
                continue
            pt = np.asarray(data.lep_pt[i][:4], dtype=float) / 1000.0
            if not (pt[0] > 20 and pt[1] > 15 and pt[2] > 10):
                continue
            stages["pt"] += 1
            if not np.all(np.asarray(data.lep_isTightID[i][:4], dtype=bool)):
                continue
            stages["tight_id"] += 1
            lepton_type = np.abs(np.asarray(data.lep_type[i][:4], dtype=int))
            if int(lepton_type.sum()) not in (44, 48, 52):
                continue
            stages["flavour"] += 1
            if int(np.asarray(data.lep_charge[i][:4]).sum()) != 0:
                continue
            stages["charge"] += 1
            iso = (
                np.asarray(data.lep_ptcone30[i][:4]) / np.asarray(data.lep_pt[i][:4]) < 0.30
            ) & (
                np.asarray(data.lep_etcone20[i][:4]) / np.asarray(data.lep_pt[i][:4]) < 0.30
            )
            if not np.all(iso):
                continue
            stages["isolation"] += 1
            eta = np.asarray(data.lep_eta[i][:4], dtype=float)
            phi = np.asarray(data.lep_phi[i][:4], dtype=float)
            energy = np.asarray(data.lep_E[i][:4], dtype=float) / 1000.0
            px = pt * np.cos(phi)
            py = pt * np.sin(phi)
            pz = pt * np.sinh(eta)
            mass2 = energy.sum() ** 2 - px.sum() ** 2 - py.sum() ** 2 - pz.sum() ** 2
            masses.append(math.sqrt(max(0.0, float(mass2))))
    return np.asarray(masses), stages


def run_cern(data_root: Path, out: Path) -> dict:
    files = sorted((data_root / "cern").glob("data_?.4lep.root"))
    if len(files) != 4:
        raise FileNotFoundError("Expected four ATLAS data_?.4lep.root files")
    masses, cutflow = extract_four_lepton_masses(files)
    edges = np.arange(105.0, 160.0 + 2.5, 2.5)
    counts, _ = np.histogram(masses, bins=edges)
    best = fit_cern_counts(counts, edges)
    null = fit_cern_counts(counts, edges, fixed_signal=0.0)
    q0 = max(0.0, 2.0 * (null.nll - best.nll))
    local_z = math.sqrt(q0)

    def profile_delta(signal_yield: float, threshold: float = 1.0) -> float:
        return 2.0 * (fit_cern_counts(counts, edges, fixed_signal=signal_yield).nll - best.nll) - threshold

    if best.signal_yield > 0 and profile_delta(0.0) > 0:
        low68 = brentq(profile_delta, 0.0, best.signal_yield)
    else:
        low68 = 0.0
    hi = max(best.signal_yield + 1.0, 2.0)
    while profile_delta(hi) < 0:
        hi *= 1.8
    high68 = brentq(profile_delta, best.signal_yield, hi)

    # Frozen injection/recovery campaign under the fitted background model.
    rng = np.random.default_rng(SEED)
    signal_prob, background_prob = cern_probabilities(edges, best.slope_per_gev)
    injected_yield = 10.0
    expectation = injected_yield * signal_prob + best.background_yield * background_prob
    recovered = []
    for _ in range(500):
        toy = rng.poisson(expectation)
        recovered.append(fit_cern_counts(toy, edges).signal_yield)
    rec = np.percentile(recovered, [2.5, 16, 50, 84, 97.5])

    centers = 0.5 * (edges[:-1] + edges[1:])
    signal_prob, background_prob = cern_probabilities(edges, best.slope_per_gev)
    background = best.background_yield * background_prob
    signal = best.signal_yield * signal_prob
    hist_rows = [
        {
            "m_low_gev": edges[i],
            "m_high_gev": edges[i + 1],
            "observed": counts[i],
            "background_fit": background[i],
            "signal_fit": signal[i],
            "total_fit": background[i] + signal[i],
        }
        for i in range(len(counts))
    ]
    write_csv(out / "cern_four_lepton_histogram.csv", hist_rows)
    write_csv(out / "cern_selected_four_lepton_masses.csv", [{"m4l_gev": v} for v in masses])

    fig, ax = plt.subplots(figsize=(7.2, 4.2))
    ax.errorbar(centers, counts, yerr=np.sqrt(np.maximum(counts, 1)), fmt="o", color=INK, label="ATLAS open collision data")
    ax.step(edges, np.r_[background, background[-1]], where="post", color=GRAY, lw=1.7, label="Profiled smooth background")
    ax.step(edges, np.r_[background + signal, (background + signal)[-1]], where="post", color=BLUE, lw=2.0, label="Background + fixed 125 GeV template")
    ax.fill_between(edges, 0, np.r_[signal, signal[-1]], step="post", color=GOLD, alpha=0.25, label="Fitted scalar component")
    ax.set(xlabel=r"Four-lepton invariant mass $m_{4\ell}$ [GeV]", ylabel="Events / 2.5 GeV", title="ATLAS 2016 four-lepton scalar-template check")
    ax.grid(True, axis="y")
    ax.legend(fontsize=8)
    fig.tight_layout()
    fig.savefig(out / "figure_cern_four_lepton.png", bbox_inches="tight")
    plt.close(fig)

    return {
        "dataset": "ATLAS 13 TeV 2016 educational four-lepton release, 10 fb^-1",
        "input_events_preselected": cutflow["input_preselected"],
        "selected_events": int(len(masses)),
        "fit_window_gev": [float(edges[0]), float(edges[-1])],
        "bin_width_gev": 2.5,
        "template_mass_gev": 125.0,
        "template_resolution_gev": 2.5,
        "best_fit_signal_yield": best.signal_yield,
        "profile_68_interval_signal_yield": [low68, high68],
        "best_fit_background_yield": best.background_yield,
        "best_fit_background_slope_per_gev": best.slope_per_gev,
        "diagnostic_local_q0": q0,
        "diagnostic_local_z": local_z,
        "injection_yield": injected_yield,
        "injection_recovery_median": rec[2],
        "injection_recovery_68_interval": [rec[1], rec[3]],
        "injection_recovery_95_interval": [rec[0], rec[4]],
        "cutflow": cutflow,
        "inference_scope": "educational standard-candle and likelihood validation; not an END-specific resonance prediction",
        "files": {p.name: sha256(p) for p in files},
    }


# ---------------------------------------------------------------------------
# XENONnT D05: nominal-yield binned profile likelihood reproduction


def load_xenon_template(base: Path, sr: int, source: str, mass: int | None = None):
    if source == "ac":
        rel = f"data/ac/template_XENONnT_sr{sr}_ac_cevns.h5"
    elif source == "b8":
        rel = f"data/b8/template_XENONnT_sr{sr}_b8_cevns_tly_0.0_tqy_0.0.h5"
    elif source == "er":
        rel = f"data/er/template_XENONnT_sr{sr}_er_cevns.h5"
    elif source == "nr":
        rel = f"data/rg/template_XENONnT_sr{sr}_rg_cevns_tly_0.0_tqy_0.0.h5"
    elif source == "signal":
        rel = f"data/wimp_si/template_XENONnT_sr{sr}_wimp_si_{mass}_cevns_tly_0.0_tqy_0.0.h5"
    else:
        raise KeyError(source)
    path = base / rel
    with h5py.File(path, "r") as h:
        bins = [h[f"bins/{i}"][:] for i in range(4)]
        template = h["templates/template"][:]
    return path, bins, template


def run_xenon(data_root: Path, out: Path) -> dict:
    base = data_root / "xenon" / "light_wimp_data_release" / "light_wimp_data_release"
    real_path = base / "data" / "real_data.pkl"
    with real_path.open("rb") as stream:
        real = pickle.load(stream)
    livetimes = [1.174398, 2.342593]
    observed = []
    components = {k: [] for k in ("ac", "b8", "er", "nr")}
    used_files = {real_path.name: sha256(real_path)}
    for sr in (0, 1):
        _, bins, _ = load_xenon_template(base, sr, "ac")
        event_array = np.column_stack([real[sr][name] for name in real[sr].dtype.names[:4]])
        observed.append(np.histogramdd(event_array, bins=bins)[0].ravel())
        for source in components:
            path, _, template = load_xenon_template(base, sr, source)
            components[source].append(template.ravel() * livetimes[sr])
            used_files[str(path.relative_to(base))] = sha256(path)

    prior_sigma = np.array([0.0381, 0.0904, 0.0576, 1.0, 0.58, 0.25])
    nuisance_bounds = [(1e-4, 100)] * 5 + [(0.5, 1.5)]
    masses = [3, 4, 5, 6, 8, 10, 12]
    official_path = base / "limits" / "XENONnT_2024_SI.csv"
    official = pd.read_csv(official_path).set_index("dm_mass_gev")
    used_files[str(official_path.relative_to(base))] = sha256(official_path)
    result_rows = []

    for mass in masses:
        signal = []
        for sr in (0, 1):
            path, _, template = load_xenon_template(base, sr, "signal", mass)
            signal.append(template.ravel() * livetimes[sr])
            used_files[str(path.relative_to(base))] = sha256(path)

        def nuisance_nll(x, mu):
            b8_mult, ac0_mult, ac1_mult, er_mult, nr_mult, efficiency = x
            value = 0.5 * float(np.sum(((x - 1.0) / prior_sigma) ** 2))
            for sr in (0, 1):
                expectation = (
                    mu * efficiency * signal[sr]
                    + efficiency * b8_mult * components["b8"][sr]
                    + (ac0_mult if sr == 0 else ac1_mult) * components["ac"][sr]
                    + er_mult * components["er"][sr]
                    + nr_mult * components["nr"][sr]
                )
                expectation = np.clip(expectation, 1e-12, None)
                n = observed[sr]
                value += float(np.sum(expectation - n * np.log(expectation) + gammaln(n + 1)))
            return value

        def profile(mu: float) -> tuple[float, np.ndarray]:
            fit = minimize(
                nuisance_nll,
                np.ones(6),
                args=(mu,),
                method="L-BFGS-B",
                bounds=nuisance_bounds,
                options={"ftol": 1e-11, "maxiter": 2000},
            )
            if not fit.success:
                raise RuntimeError(f"XENON nuisance fit failed at mass {mass}, mu {mu}: {fit.message}")
            return float(fit.fun), fit.x

        coarse = np.array([0.0, 0.01, 0.1, 1, 10, 100, 1000, 10000, 100000], dtype=float)
        coarse_nll = np.array([profile(mu)[0] for mu in coarse])
        j = int(np.argmin(coarse_nll))
        lower = coarse[max(0, j - 1)]
        upper = coarse[min(len(coarse) - 1, j + 1)]
        refine = minimize_scalar(
            lambda log_mu: profile(float(np.expm1(log_mu)))[0],
            bounds=(float(np.log1p(lower)), float(np.log1p(upper))),
            method="bounded",
            options={"xatol": 1e-8},
        )
        mu_hat = float(np.expm1(refine.x))
        best_nll = float(refine.fun)
        threshold = float(chi2.ppf(0.90, 1))

        def upper_equation(mu):
            return 2.0 * (profile(float(mu))[0] - best_nll) - threshold

        high = max(mu_hat + 1.0, 1.0)
        while upper_equation(high) < 0:
            high *= 2.0
            if high > 1e7:
                raise RuntimeError("XENON upper-limit bracket failed")
        mu_ul = float(brentq(upper_equation, mu_hat, high, xtol=1e-7))

        spectrum_path = base / "data" / "signal" / f"wimp_si_spectrum_{mass}.json"
        with spectrum_path.open() as stream:
            spectrum = json.load(stream)
        reference_sigma = float(spectrum["sigma"])
        used_files[str(spectrum_path.relative_to(base))] = sha256(spectrum_path)
        approximate_limit = mu_ul * reference_sigma
        official_limit = float(official.loc[mass, "upper_limit_cm2"])
        difference_percent = 100.0 * (approximate_limit / official_limit - 1.0)
        nominal_range = 5 <= mass <= 12
        result_rows.append(
            {
                "dm_mass_gev": mass,
                "best_fit_signal_multiplier": mu_hat,
                "approximate_90cl_limit_cm2": approximate_limit,
                "official_90cl_limit_cm2": official_limit,
                "difference_percent": difference_percent,
                "nominal_yield_reproduction_scope": nominal_range,
                "within_30_percent_if_in_scope": bool(nominal_range and abs(difference_percent) <= 30.0),
            }
        )

    write_csv(out / "xenonnt_limit_reproduction.csv", result_rows)
    approx = pd.DataFrame(result_rows)
    fig, axes = plt.subplots(1, 2, figsize=(9.2, 3.8), gridspec_kw={"width_ratios": [1.45, 1]})
    ax = axes[0]
    ax.loglog(official.index, official["upper_limit_cm2"], color=INK, lw=1.8, label="XENONnT official 90% CL")
    ax.loglog(approx.dm_mass_gev, approx.approximate_90cl_limit_cm2, "o--", color=BLUE, label="This nominal-yield profile reproduction")
    ax.set(xlabel=r"Dark-matter mass $m_\chi$ [GeV]", ylabel=r"Spin-independent cross section [cm$^2$]", title="XENONnT released-limit check")
    ax.grid(True, which="both")
    ax.legend(fontsize=8)

    ax = axes[1]
    colors = [GRAY if m < 5 else TEAL for m in approx.dm_mass_gev]
    ax.bar(approx.dm_mass_gev.astype(str), approx.difference_percent, color=colors)
    ax.axhspan(-30, 30, color=GOLD, alpha=0.14, label="Declared ±30% reproduction band")
    ax.axhline(0, color=INK, lw=0.9)
    ax.set(xlabel=r"$m_\chi$ [GeV]", ylabel="Difference from official [%]", title="Approximation residual")
    ax.grid(True, axis="y")
    ax.legend(fontsize=8)
    fig.tight_layout()
    fig.savefig(out / "figure_xenonnt_reproduction.png", bbox_inches="tight")
    plt.close(fig)

    in_scope = approx[approx.nominal_yield_reproduction_scope]
    return {
        "dataset": "XENONnT 2025 light-WIMP data release",
        "unblinded_events_sr0": int(observed[0].sum()),
        "unblinded_events_sr1": int(observed[1].sum()),
        "tested_masses_gev": masses,
        "in_scope_masses_gev": in_scope.dm_mass_gev.astype(int).tolist(),
        "maximum_absolute_difference_percent_in_scope": float(in_scope.difference_percent.abs().max()),
        "all_in_scope_points_within_30_percent": bool(in_scope.within_30_percent_if_in_scope.all()),
        "approximation": "nominal light/charge yield; released 4D templates; Gaussian rate/efficiency constraints; asymptotic q=2.7055",
        "excluded_from_reproduction_claim": "3-4 GeV, where threshold yield morphing and boundary behavior dominate",
        "inference_scope": "published WIMP benchmark reproduction; no END particle mass or cross section was supplied",
        "files": used_files,
    }


# ---------------------------------------------------------------------------
# DUNE D06: official-release forecast and Hamiltonian integrity tests


def pmns_matrix(theta12, theta13, theta23, delta):
    c12, c13, c23 = np.cos([theta12, theta13, theta23])
    s12, s13, s23 = np.sin([theta12, theta13, theta23])
    return np.array(
        [
            [c12 * c13, s12 * c13, s13 * np.exp(-1j * delta)],
            [
                -s12 * c23 - c12 * s23 * s13 * np.exp(1j * delta),
                c12 * c23 - s12 * s23 * s13 * np.exp(1j * delta),
                s23 * c13,
            ],
            [
                s12 * s23 - c12 * c23 * s13 * np.exp(1j * delta),
                -c12 * s23 - s12 * c23 * s13 * np.exp(1j * delta),
                c23 * c13,
            ],
        ],
        dtype=complex,
    )


def dune_probability(energy_gev, delta_cp, antineutrino=False):
    sin2_theta12 = 0.307
    sin2_theta13 = 0.0220
    sin2_theta23 = 0.561
    theta12, theta13, theta23 = np.arcsin(np.sqrt([sin2_theta12, sin2_theta13, sin2_theta23]))
    dm21 = 7.53e-5
    dm31 = 2.455e-3
    density = 2.848
    baseline_km = 1284.9
    u = pmns_matrix(theta12, theta13, theta23, -delta_cp if antineutrino else delta_cp)
    # a = 2 E V in eV^2; 7.56e-5 assumes the standard crustal electron fraction.
    matter_a = 7.56e-5 * density * energy_gev * (-1 if antineutrino else 1)
    mass2 = u @ np.diag([0.0, dm21, dm31]) @ u.conj().T + np.diag([matter_a, 0.0, 0.0])
    hamiltonian_ev = mass2 / (2.0 * energy_gev * 1e9)
    length_ev_inverse = baseline_km * 5.0677307e9
    evolution = expm(-1j * hamiltonian_ev * length_ev_inverse)
    probability_mu_to_e = float(abs(evolution[0, 1]) ** 2)
    unitarity_error = float(np.max(np.abs(evolution.conj().T @ evolution - np.eye(3))))
    return probability_mu_to_e, unitarity_error


def run_dune(data_root: Path, out: Path) -> dict:
    dune = data_root / "dune"
    flux_path = dune / "anc" / "dune_globes" / "flux" / "histos_g4lbne_v3r5p4_QGSP_BERT_OptimizedEngineeredNov2017_neutrino_LBNEFD_globes_flux.txt"
    xsec_path = dune / "anc" / "dune_globes" / "xsec" / "xsec_cc.dat"
    config_path = dune / "anc" / "dune_globes" / "DUNE_GLoBES.glb"
    flux = np.loadtxt(flux_path)
    xsec = np.loadtxt(xsec_path)
    energy = np.linspace(0.5, 8.0, 301)
    deltas = [-math.pi / 2, 0.0, math.pi / 2]
    curves = {}
    max_unitarity = 0.0
    for delta in deltas:
        nu = []
        antinu = []
        for e in energy:
            p, err = dune_probability(e, delta, False)
            q, errbar = dune_probability(e, delta, True)
            nu.append(p)
            antinu.append(q)
            max_unitarity = max(max_unitarity, err, errbar)
        curves[delta] = (np.asarray(nu), np.asarray(antinu))

    # A detector-independent appearance proxy using the official FHC nu_mu
    # flux and released CC cross section.  It intentionally stops before
    # efficiency, migration, near-detector constraints, and nuisance fits.
    flux_energy = flux[:, 0]
    flux_numu = flux[:, 2]
    sigmahat_nue = np.interp(np.log10(flux_energy), xsec[:, 0], xsec[:, 1])
    sigma_nue_arbitrary = sigmahat_nue * flux_energy
    probability_flux = np.array([dune_probability(e, -math.pi / 2, False)[0] for e in flux_energy])
    appearance_proxy = flux_numu * sigma_nue_arbitrary * probability_flux
    valid = (flux_energy >= 0.5) & (flux_energy <= 8.0)
    proxy_peak = float(flux_energy[valid][np.argmax(appearance_proxy[valid])])
    proxy_mean = float(np.average(flux_energy[valid], weights=appearance_proxy[valid]))

    rows = []
    for i, e in enumerate(energy):
        row = {"energy_gev": e}
        for delta, label in zip(deltas, ("minus_pi_over_2", "zero", "plus_pi_over_2")):
            row[f"p_numu_to_nue_delta_{label}"] = curves[delta][0][i]
            row[f"p_numubar_to_nuebar_delta_{label}"] = curves[delta][1][i]
        rows.append(row)
    write_csv(out / "dune_oscillation_probabilities.csv", rows)
    write_csv(
        out / "dune_flux_appearance_proxy.csv",
        [
            {
                "energy_gev": flux_energy[i],
                "official_fhc_numu_flux": flux_numu[i],
                "p_numu_to_nue_delta_minus_pi_over_2": probability_flux[i],
                "cc_xsec_relative": sigma_nue_arbitrary[i],
                "flux_x_xsec_x_probability_proxy": appearance_proxy[i],
            }
            for i in range(len(flux_energy))
        ],
    )

    fig, axes = plt.subplots(1, 2, figsize=(9.2, 3.8))
    ax = axes[0]
    for delta, label, color in zip(deltas, (r"-\pi/2", "0", r"+\pi/2"), (BLUE, INK, RED)):
        ax.plot(energy, curves[delta][0], color=color, lw=1.6, label=rf"$\delta_{{CP}}={label}$")
    ax.set(xlabel="Neutrino energy [GeV]", ylabel=r"$P(\nu_\mu\rightarrow\nu_e)$", title="Three-flavour matter propagation")
    ax.set_ylim(bottom=0)
    ax.grid(True)
    ax.legend(fontsize=8)

    ax = axes[1]
    f = flux_numu[valid] / flux_numu[valid].max()
    p = appearance_proxy[valid] / appearance_proxy[valid].max()
    ax.step(flux_energy[valid], f, where="mid", color=GRAY, lw=1.4, label=r"Official FHC $\nu_\mu$ flux (normalized)")
    ax.step(flux_energy[valid], p, where="mid", color=TEAL, lw=1.8, label=r"$\Phi\,\sigma\,P_{\mu e}$ proxy (normalized)")
    ax.axvline(proxy_peak, color=GOLD, lw=1.0, ls="--", label=f"Proxy peak {proxy_peak:.3g} GeV")
    ax.set(xlabel="Neutrino energy [GeV]", ylabel="Relative intensity", title="DUNE release-based appearance proxy")
    ax.grid(True)
    ax.legend(fontsize=8)
    fig.tight_layout()
    fig.savefig(out / "figure_dune_forecast.png", bbox_inches="tight")
    plt.close(fig)

    p25 = {}
    for delta, label in zip(deltas, ("minus_pi_over_2", "zero", "plus_pi_over_2")):
        p25[label] = {
            "neutrino": dune_probability(2.5, delta, False)[0],
            "antineutrino": dune_probability(2.5, delta, True)[0],
        }
    return {
        "dataset": "DUNE Collaboration TDR-approximation ancillary flux/GLoBES release (arXiv:2103.04797)",
        "baseline_km": 1284.9,
        "constant_density_g_cm3": 2.848,
        "probabilities_at_2_5_gev": p25,
        "maximum_unitarity_error": max_unitarity,
        "appearance_proxy_peak_energy_gev": proxy_peak,
        "appearance_proxy_mean_energy_gev": proxy_mean,
        "empirical_data_available": False,
        "inference_scope": "equation/unitarity and official-release forecast only; H_END_nu fixed to zero; no DUNE data-validation claim",
        "files": {p.name: sha256(p) for p in (flux_path, xsec_path, config_path)},
    }


def main() -> None:
    parser = argparse.ArgumentParser()
    parser.add_argument("--data-root", type=Path, default=Path("data"))
    parser.add_argument("--output-dir", type=Path, default=Path("results"))
    args = parser.parse_args()
    args.output_dir.mkdir(parents=True, exist_ok=True)
    configure_plotting()

    results = {
        "study": "END Theory 2.0 public-data companion validation",
        "analysis_version": "1.0.0",
        "random_seed": SEED,
        "python": platform.python_version(),
        "platform": platform.platform(),
        "packages": {
            "numpy": np.__version__,
            "scipy": scipy.__version__,
            "pandas": pd.__version__,
            "matplotlib": matplotlib.__version__,
            "uproot": uproot.__version__,
            "h5py": h5py.__version__,
        },
        "claim_boundary": (
            "The calculations validate selected analysis contracts and public-release reproductions. "
            "They do not validate END Theory as a closed physical theory because A_alpha, C/Lambda^2, "
            "an END dark-sector mass/cross section, and H_END_nu are not numerically frozen."
        ),
    }
    results["gwosc"] = run_gwosc(args.data_root, args.output_dir)
    results["cern"] = run_cern(args.data_root, args.output_dir)
    results["xenonnt"] = run_xenon(args.data_root, args.output_dir)
    results["dune"] = run_dune(args.data_root, args.output_dir)
    results["decision"] = {
        "analysis_pipeline_validated": True,
        "end_theory_empirically_validated": False,
        "reason": "No experiment-specific END parameter card or micro-to-observable coefficient was present in the source manuscript.",
    }
    with (args.output_dir / "validation_results.json").open("w") as stream:
        json.dump(json_ready(results), stream, indent=2, sort_keys=True)
    print(json.dumps(json_ready(results["decision"]), indent=2))


if __name__ == "__main__":
    main()
