#!/usr/bin/env python3
"""
verify_P108.py -- Addendum 108: non-perturbative resonance correction.

This verifier checks the numerical chain in 108_Addendum_ResonanceCorrection.tex:
the P107 light-quark pQCD baseline, the residual by subtraction from the PDG
dispersive value, the rho mass and f_V^2 estimates, the resonance table sums,
and the appendix alpha(MZ) running chain.

The arithmetic mostly reproduces. The flagged issues are internal/status
problems: the residual is defined using the PDG dispersive value rather than
independently derived; the printed GS resonance integral is 0.00365, not the
claimed residual 0.002383; the omega charge-factor prose has an arithmetic/sign
error; and the paper's own open questions leave SU(3)-breaking and higher
resonance contributions unresolved.
"""

from __future__ import annotations

import math
import sys
from pathlib import Path

PASS = FAIL = 0
_N = 0


def check(n, desc, cond):
    global PASS, FAIL
    ok = bool(cond); PASS += ok; FAIL += (not ok)
    print(f"  [{'PASS' if ok else 'FAIL'}] {n:>2}. {desc}")
    return ok


class Verifier:
    """Target-style output adapter; check/record numerics are unchanged."""

    def __init__(self, name):
        print(name)

    def check(self, 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 self.record(label, ok, computed, claimed,
                           err_detail + (f"; {detail}" if detail else ""))

    def record(self, label, ok, computed="", claimed="", detail=""):
        global _N
        ok = bool(ok)
        desc, extra = label, detail
        if not ok and detail:
            desc, extra = f"{label} -- {detail}", ""
        _N += 1
        check(_N, desc, ok)
        if computed != "" or claimed != "":
            print(f"       computed: {computed}")
            print(f"       claimed : {claimed}")
        if extra:
            print(f"       {extra}")
        return ok

    def summary(self):
        print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
        return 1 if FAIL else 0


v = Verifier("P108 -- Resonance Correction")

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

PI = math.pi
MU0 = 4 * PI**3 + PI**2 + PI
ALPHA = 1.0 / MU0
MZ = 91.1876
MCONF = 0.220
ALPHA_S = 0.118
DA_UDS_EXP = 0.019073


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


da_uds_pqcd = (
    ALPHA
    / (3.0 * PI)
    * 2.0
    * (2.0 * math.log(MZ / MCONF) - 5.0 / 3.0)
    * (1.0 + ALPHA_S / PI)
)
da_residual = DA_UDS_EXP - da_uds_pqcd

v.check("P107 uds pQCD baseline", da_uds_pqcd, 0.016690, rel=4e-5)
v.check("rounded residual 0.01907 - 0.01669", 0.01907 - 0.01669, 0.00238, rel=1e-12)
v.check("residual using printed 0.019073 input", da_residual, 0.002383, rel=2e-4)

mrho_toe_gev = 2.0 * math.sqrt(PI) * MCONF
mrho_pdg_gev = 0.77526
v.check("m_rho = 2 sqrt(pi) m_conf in GeV", mrho_toe_gev, 0.780, rel=2e-4)
v.check("m_rho in MeV", 1000.0 * mrho_toe_gev, 779.9, rel=3e-5)
v.check("m_rho residual percent", pct(mrho_toe_gev, mrho_pdg_gev), 0.60, rel=1e-2)

fv2 = 8.0 * PI
frho2_exp = (4.0 * PI * ALPHA**2 / 3.0) * (0.77526 / 7.04e-6)
v.check("f_V^2 = 8 pi", fv2, 25.13, rel=2e-4)
v.check("f_rho^2 from Gamma(rho -> e+e-)", frho2_exp, 24.56, rel=2e-4)
v.check("f_rho^2 residual percent", pct(frho2_exp, fv2), -2.3, rel=2e-2)

c_rho = 1.0 / math.sqrt(2.0)
c_omega = 1.0 / (3.0 * math.sqrt(2.0))
c_phi = 1.0 / 3.0
v.check("C_rho", c_rho, 1.0 / math.sqrt(2.0), rel=1e-12)
v.check("C_omega", c_omega, 1.0 / (3.0 * math.sqrt(2.0)), rel=1e-12)
v.check("C_phi", c_phi, 1.0 / 3.0, rel=1e-12)
v.check(
    "omega prose charge expression",
    (2.0 / 3.0 + 1.0 / 3.0) / math.sqrt(2.0),
    c_omega,
    rel=1e-12,
    detail="Expected fail: with the plus sign printed in the prose this is 1/sqrt(2), not 1/(3sqrt(2)); the correct quark-charge sum uses +2/3 + (-1/3).",
)
v.check("omega charge with physical d-quark sign", (2.0 / 3.0 - 1.0 / 3.0) / math.sqrt(2.0), c_omega, rel=1e-12)

rho_nra = 3.73e-3
omega_nra = 2.81e-4
phi_nra = 4.98e-4
rho_gs = 2.87e-3
omega_gs = 2.80e-4
phi_gs = 4.93e-4
nra_sum = rho_nra + omega_nra + phi_nra
gs_sum = rho_gs + omega_gs + phi_gs

v.check("NRA resonance table sum", nra_sum, 4.51e-3, rel=5e-4)
v.check("GS resonance table sum", gs_sum, 3.65e-3, rel=2e-3)
v.check("rho GS below NRA percent", 100.0 * (1.0 - rho_gs / rho_nra), 23.0, rel=4e-3)
v.record(
    "omega and phi GS estimates are within about 1 percent of NRA",
    max(abs(omega_gs / omega_nra - 1.0), abs(phi_gs / phi_nra - 1.0)) <= 0.011,
    computed=f"omega={100*(omega_gs/omega_nra-1):+.3f}%, phi={100*(phi_gs/phi_nra-1):+.3f}%",
    claimed="within 1%",
    detail="The rounded phi row is just over 1.0%, so this check allows rounding slack.",
)
v.check(
    "GS integral equals claimed residual",
    3.65e-3,
    da_residual,
    rel=2e-4,
    detail="Expected fail: the printed GS integral is about 1.53 times the residual and cannot be the same contribution without an extra subtraction/double-counting argument.",
)
v.check(
    "GS integral / residual ratio",
    3.65e-3 / da_residual,
    1.0,
    rel=2e-4,
    detail="Expected fail: the resonance integral is larger than the residual it is meant to explain.",
)

components = {
    "lep": 0.031497,
    "uds_pqcd": 0.016690,
    "res": 0.002383,
    "cb": 0.008563,
    "top": -0.000072,
}
delta_total = sum(components.values())
alpha0_inv_printed = 137.036304
alpha_mz_inv = alpha0_inv_printed * (1.0 - delta_total)
v.check("appendix total Delta alpha", delta_total, 0.059061, rel=1e-12)
v.check("appendix alpha(MZ)^-1", alpha_mz_inv, 128.943, rel=2e-6)
v.check("appendix residual in sigma units", abs(alpha_mz_inv - 128.944) / 0.010, 0.1, rel=3e-1)
v.check("q=0 inverse-alpha gap propagated to MZ", 0.000305 * (1.0 - delta_total), 0.000287, rel=3e-3)

v.record(
    "residual is independently derived rather than defined from PDG data",
    "Dauds^{\\rm exp} - Dauds^{\\rm pQCD}" not in TEX and "measured $e^+e^-$ cross-section" not in TEX,
    computed="the paper defines da_res = Dauds_exp - Dauds_pQCD and later says the faithful low-energy evaluation requires measured e+e- data",
    claimed="This paper derives that residual from first principles",
    detail="Expected status fail: the displayed residual value is anchored by PDG input, while the geometric pieces identify a source but do not reproduce 0.002383 from only J3(O) data.",
)
v.record(
    "complete parameter-free light-quark account is closed",
    "higher resonances" not in TEX and "open question" not in TEX.lower(),
    computed="open questions remain for SU(3) breaking and broad resonances above 2 GeV",
    claimed="complete, parameter-free account of the light-quark contribution",
    detail="Expected status fail: the paper's own summary leaves material pieces unresolved.",
)
v.record(
    "no experimental input beyond Dauds_exp and mconf enters the displayed resonance derivation",
    "Gamma(\\rho\\to e^+e^-)" not in TEX and "Gounaris-Sakurai" not in TEX,
    computed="the numerical coupling test uses PDG rho leptonic width, and the GS integral is a phenomenological/experimental spectral model",
    claimed="No experimental input beyond the PDG central values Dauds_exp and mconf enters the derivation",
    detail="Expected status fail: these may be cross-checks, but they are still displayed inputs to the numerical account.",
)

sys.exit(v.summary())