"""
RETRACTED HEADLINE (2026-07 tension hunt, 3/3-verified). This script's original claim -- that the
frozen A = 2*sqrt(alpha)*R_d/(1+R_d) = 0.0667 works as well as the fit -- was a cross-pipeline
calibration transplant: this toy pipeline is miscalibrated in absolute theta* (-0.32%, i.e. -10.9
sigma vs Planck at Planck's own parameters, where no DFD physics exists) and was compared against a
CAMB-convention measurement. At the consistent convention the frozen A=0.0667 gives theta* at
+12.2 sigma and plik-lite dchi2 = +202/+172 (per-bin variants). A is ONE FITTED amplitude:
A* = 0.0908 (-26.6% off the 2*sqrt(alpha)-gate value); the fitted headline (+1.0 sigma,
dchi2 = +23.3/+8.0) stands. The screen-factor RATIO response remains trustworthy (toy vs CAMB
ratios agree to 0.001%). Numeric logic below is untouched; only this equivalence claim is retracted.

Also tests the referee's LITERAL multiplicative prescription D_M^opt=e^{Dpsi(z)}D_M vs the
integrated-LOS kernel, because the SHAPE/picture interacts with whether theta* is pinned.
"""
import numpy as np
from scipy.integrate import cumulative_trapezoid
import dfd_params as P, background as B

C, H0 = P.C_KMS, P.H0
ALPHA = P.ALPHA
R_d = 0.640
A_DERIVED = 2*np.sqrt(ALPHA)*R_d/(1+R_d)     # = 0.0667, FROZEN, no fit
A_FIT = 0.0685                                # the earlier best fit (for reference only)
THP, THE = 1.04109, 0.00030

zg = np.linspace(0.0, P.Z_STAR, 120000)
cH = C/(H0*np.array([B.E(z) for z in zg]))
DM_raw = cumulative_trapezoid(cH, zg, initial=0.0)
rs_star = B.sound_horizon(P.Z_STAR); rd = B.r_drag()
def cH_at(z): return C/(H0*B.E(z))
OL = 1 - P.OMEGA_M - P.OMEGA_R_H2/(H0/100)**2
def shape(z): return OL/(OL + P.OMEGA_M*(1+z)**3)   # Omega_Lambda(z), the DE-fraction shape

def run(A, mode):
    dps = A*shape(zg)
    if mode == 'integrated':                  # transverse = integral of local index (LOS-correct)
        DM_opt = cumulative_trapezoid(np.exp(dps)*cH, zg, initial=0.0)
    else:                                     # multiplicative (the referee's literal D_M^opt=e^{Dpsi}D_M)
        DM_opt = np.exp(dps)*DM_raw
    theta = 100*rs_star/DM_opt[-1]
    bao = 0.0
    for b in B.BAO:
        z = b["z"]; DMo = np.interp(z, zg, DM_opt)
        if mode == 'integrated':
            DHo = np.exp(A*shape(z))*cH_at(z)                       # local radial
        else:
            dpsidz = A*(shape(z+1e-4)-shape(z-1e-4))/2e-4
            DHo = np.exp(A*shape(z))*(cH_at(z)+DMo*dpsidz)          # the referee derivative-corrected
        pred = DMo/rd if b["kind"]=="DM" else DHo/rd if b["kind"]=="DH" else (z*DMo**2*DHo)**(1/3)/rd
        bao += ((pred-b["val"])/b["err"])**2
    return theta, bao                                              # report theta* and PURE BAO separately

print("="*78)
print("FROZEN-AMPLITUDE TEST: A_DFD = 2*sqrt(alpha)*R_d/(1+R_d), NO refit")
print("="*78)
print(f"\n  2*sqrt(alpha)={2*np.sqrt(ALPHA):.5f}  R_d/(1+R_d)={R_d/(1+R_d):.5f}  ->  A_DERIVED={A_DERIVED:.5f}")
print(f"  (fitted reference A_FIT={A_FIT:.4f}; the derived value is {100*(A_DERIVED-A_FIT)/A_FIT:+.1f}% from it)\n")
print(f"  {'picture':14s} {'A':>8s} {'100*theta*':>11s} {'theta* dev':>11s} {'BAO chi2/dof':>13s}")
for mode in ('integrated','multiplicative'):
    for A,lab in ((A_DERIVED,'derived'),(A_FIT,'fit')):
        th,c = run(A,mode)
        nb = len(B.BAO)
        print(f"  {mode:14s} {A:8.5f} {th:11.5f} {(th-THP)/THE:+9.1f}σ {c/nb:13.2f}   [{lab}]")

print("\n" + "-"*78)
th_d,c_d = run(A_DERIVED,'integrated')
print(f"  RESULT (integrated LOS kernel, the physically-correct picture):")
print(f"    FROZEN derived A=0.0667 gives 100*theta* = {th_d:.5f} = {(th_d-THP)/THE:+.1f}sigma from Planck,")
print(f"    BAO chi^2/dof = {c_d/len(B.BAO):.2f}. No refit. Compare fitted A=0.0685: theta* perfect, BAO 0.78.")
print(f"  The multiplicative-literal picture does NOT pin theta* (Omega_L shape -> Dpsi(z*)~0), confirming")
print(f"  theta* is fixed by the INTEGRATED low-z screen, not a local kick at z*.")
print(f"  Intrinsic TT/TE/EE (heights/ratios/pol) are A-INDEPENDENT (gate_spectrum.py: <0.3%/1.5%); only")
print(f"  theta* depends on A. So the frozen derived amplitude is the WHOLE test, and it lands at {(th_d-THP)/THE:+.1f}sigma.")
