#!/usr/bin/env python3
"""
verify_P107.py -- Addendum 107: W-mass closure from hadronic running.

This verifier checks the pQCD hadronic-vacuum-polarisation arithmetic in
107_Addendum_WMassClosure.tex.  The displayed delta-alpha contributions,
running-alpha value, and linear W-mass update reproduce.  The flagged issues
are status/proof limitations: the paper says OP3 is completed, while the
non-perturbative resonance spectral function remains open; the light-quark
cutoff is called a UV cutoff in one place even though it is used as an IR
cutoff; and the W-mass update is a linear sensitivity estimate rather than a
fresh full electroweak calculation.
"""

from __future__ import annotations

import math
import sys
from pathlib import Path



PASS = FAIL = 0
_N = 0

def record(label, ok, computed="", claimed="", detail=""):
    """Modern-format check line; behavior-preserving port of verify_common."""
    global PASS, FAIL, _N
    _N += 1
    ok = bool(ok)
    desc = label
    if ok:
        PASS += 1
    else:
        FAIL += 1
        if "Expected" in detail:
            i = detail.find("Expected")
            desc = f"{label} -- {detail[i:]}"
            detail = detail[:i].rstrip().rstrip(";")
    print(f"  [{'PASS' if ok else 'FAIL'}] {_N:>2}. {desc}")
    if computed != "" or claimed != "":
        print(f"        computed: {computed}")
        print(f"        claimed : {claimed}")
    if detail:
        print(f"        {detail}")
    return ok

def check(label, computed, claimed, *, rel=1e-3, abs_tol=None, detail=""):
    if abs_tol is not None:
        ok = abs(computed - claimed) <= abs_tol
        err_detail = f"abs err={abs(computed - claimed):.6g}, tol={abs_tol:.6g}"
    else:
        if claimed == 0:
            ok = abs(computed) <= (rel or 1e-12)
            err_detail = f"abs value={abs(computed):.6g}, tol={rel:.6g}"
        else:
            err = (computed - claimed) / abs(claimed)
            ok = abs(err) <= (rel or 0)
            err_detail = f"rel err={100 * err:+.6g}%, tol={100 * (rel or 0):.6g}%"
    return record(label, ok, computed, claimed, err_detail + (f"; {detail}" if detail else ""))

print("P107 -- W-Mass Closure")

ROOT = Path(__file__).resolve().parents[1]
TEX = (ROOT / "107_Addendum_WMassClosure.tex").read_text()

PI = math.pi
MU0 = 4 * PI**3 + PI**2 + PI
ALPHA_EM = 1.0 / MU0
ALPHA_S_MZ = 0.1186
MZ_MEV = 91188.0
MCONF_MEV = PI * MU0 * 0.511
MC_MEV = 1275.0
MB_MEV = 4170.0
DALPHA_LEPT = 0.03142
DALPHA_HAD_PDG = 0.02764
MW_BASELINE = 79.04
MW_PDG = 80.379
SENS = 41.61


def delta_alpha(charge_factor: float, mass_mev: float) -> float:
    bracket = 2.0 * math.log(MZ_MEV / mass_mev) - 5.0 / 3.0
    return ALPHA_EM / (3.0 * PI) * charge_factor * bracket * (1.0 + ALPHA_S_MZ / PI)


def pct(value: float, target: float) -> float:
    return 100.0 * (value - target) / abs(target)


check("mu0", MU0, 137.036, rel=3e-6)
check("m_conf=pi*mu0*m_e in MeV", MCONF_MEV, 219.99, rel=2e-5)
check("alpha_s NLO factor", 1.0 + ALPHA_S_MZ / PI, 1.0378, rel=6e-5)
check(
    "alpha/(3pi) prefactor",
    ALPHA_EM / (3.0 * PI),
    7.772e-4,
    rel=6e-5,
    detail="Expected fail: exact 1/(137.036*3*pi) is about 7.7427e-4; the final delta-alpha rows use the exact-scale value.",
)
check("uds charge/color factor", 3.0 * (4.0 / 9.0 + 1.0 / 9.0 + 1.0 / 9.0), 2.0, rel=1e-12)

br_uds = 2.0 * math.log(MZ_MEV / 220.0) - 5.0 / 3.0
br_c = 2.0 * math.log(MZ_MEV / MC_MEV) - 5.0 / 3.0
br_b = 2.0 * math.log(MZ_MEV / MB_MEV) - 5.0 / 3.0
check("uds logarithm bracket", br_uds, 10.385, rel=3e-4)
check("charm logarithm bracket", br_c, 6.861, rel=2e-3)
check("bottom logarithm bracket", br_b, 4.505, rel=5e-4)

da_uds = delta_alpha(2.0, 220.0)
da_c = delta_alpha(4.0 / 3.0, MC_MEV)
da_b = delta_alpha(1.0 / 3.0, MB_MEV)
da_total = da_uds + da_c + da_b
check("Delta alpha uds", da_uds, 0.01669, rel=2e-4)
check("Delta alpha c", da_c, 0.00736, rel=5e-4)
check("Delta alpha b", da_b, 0.00121, rel=4e-3)
check("Delta alpha had pQCD total", da_total, 0.02526, rel=1e-4)

res = DALPHA_HAD_PDG - da_total
check("non-perturbative residual delta alpha", res, 0.00238, rel=2e-3)
check("residual percent of total", 100.0 * res / DALPHA_HAD_PDG, 8.6, rel=3e-4)

alpha_mz_inv = MU0 * (1.0 - DALPHA_LEPT - da_total)
check("alpha_EM(MZ)^-1 pQCD", alpha_mz_inv, 129.27, rel=1e-5)
check("missing inverse-alpha amount", MU0 * res, 0.326, rel=6e-4)
check("alpha inverse overshoot vs PDG", alpha_mz_inv - 128.946, 0.324, rel=5e-3)

mw_pqcd = MW_BASELINE + SENS * da_total
mw_full = MW_BASELINE + SENS * DALPHA_HAD_PDG
check("W mass pQCD linear update", mw_pqcd, 80.09, rel=3e-5)
check("W mass pQCD residual percent", pct(mw_pqcd, MW_PDG), -0.36, rel=6e-3)
check("W mass with PDG hadronic running", mw_full, 80.19, rel=2e-5)
check("W mass full-row residual percent", pct(mw_full, MW_PDG), -0.23, rel=3e-2)
check("resonance contribution to W shift", SENS * res, 0.099, rel=1e-3)

record(
    "present addendum fully completes OP3",
    False,
    computed="summary still has a -0.23% row, a final '? -> 0%' arrow, and defines OP3-Deltaalpha as open",
    claimed="present addendum completes OP3",
    detail="Expected status fail: the pQCD layer is closed, but full all-order/non-perturbative closure is explicitly not closed.",
)
record(
    "non-perturbative hadronic spectral function is derived here",
    "OP3-$\\Delta\\alpha$ open problem" not in TEX,
    computed="paper defines deriving R(s) in the resonance region as OP3-Deltaalpha open problem",
    claimed="resonance contribution closes the gap",
    detail="Expected status fail.",
)
record(
    "m_conf is consistently described as an IR cutoff",
    "UV cutoff for light quarks" not in TEX,
    computed="setup line calls m_conf a UV cutoff for light quarks",
    claimed="m_conf used as the infrared cutoff in the theorem",
    detail="Expected wording/consistency fail.",
)
record(
    "constant alpha_s(MZ) pQCD factor is an all-scale derived dispersion calculation",
    False,
    computed="formula uses a single (1+alpha_s(MZ)/pi) factor over all logarithmic thresholds",
    claimed="hadronic vacuum polarisation computed from TOE alpha_s and quark masses",
    detail="Expected proof-status fail: this is a leading pQCD estimate, not the non-perturbative dispersion integral.",
)
record(
    "W mass is recomputed from the full electroweak relation",
    False,
    computed="paper applies dMW/d(Delta alpha_had)=41.61 GeV linearly to the P104 baseline",
    claimed="complete Delta r / W-mass closure",
    detail="Expected method-status fail: the arithmetic is reproducible as a linear sensitivity update.",
)

print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
sys.exit(0 if FAIL == 0 else 1)
