"""
FREEZEOUT-WEIGHTING TEST (the referee's lemma). The gate Pi_b=R/(1+R) is forced as the trace-bearing fraction of the
photon-baryon acoustic oscillator's inertia (enthalpy (4/3)rho_g + rho_b; only baryons carry the trace). OPEN:
at WHICH epoch is it evaluated? Candidates:
  (a) freezeout/visibility surface  -> Pi_b(z_d) ~ instant value where the feature imprints   [the referee's claim]
  (b) full sound-horizon integral   -> <Pi_b>_{r_s}  (averaged over the whole pre-drag propagation)
  (c) late-time-screen weighted     -> <Pi_b>_{psi}  (where the optical screen actually lives)
The DATA require an effective gate = A_data/(2 sqrt(alpha)). If ONLY the freezeout/visibility surface matches and
the full integral is REJECTED, that is strong, physically-motivated support for the freezeout lemma (the acoustic
feature imprints at last scattering -> standard visibility weighting). Honesty-guard: report which candidates fail.
"""
import numpy as np
from scipy.integrate import quad
import dfd_params as P, background as B

ZSTAR, ZDRAG = P.Z_STAR, P.Z_DRAG
ALPHA = getattr(P,'ALPHA',1/137.035999084)
A_DATA = 0.0686                       # amplitude that pins Planck theta* (four_case_decomposition.py, Case 2)
gate_data = A_DATA/(2*np.sqrt(ALPHA))

def R(z):  return (3*P.OMEGA_B_H2)/(4*P.OMEGA_G_H2)/(1+z)
def Pi(z): return R(z)/(1+R(z))
def cs(z): return 1/np.sqrt(3*(1+R(z)))
def w_rs(z): return P.C_KMS*cs(z)/B.Hofz(z)      # sound-horizon integrand (propagation weight)

# (a) freezeout/visibility surface. Visibility g(z) sharply peaked at last scattering; model Gaussian sigma~80.
def g_vis(z, zc, sig=80.0): return np.exp(-0.5*((z-zc)/sig)**2)
def gate_vis(zc):
    num=quad(lambda z: Pi(z)*g_vis(z,zc), zc-400, zc+400, limit=300)[0]
    den=quad(lambda z: g_vis(z,zc),       zc-400, zc+400, limit=300)[0]
    return num/den
# (b) full sound-horizon integral
num=quad(lambda z: Pi(z)*w_rs(z), ZSTAR,1e6,limit=400)[0]; den=quad(w_rs,ZSTAR,1e6,limit=400)[0]
gate_int = num/den
# (c) late-time-screen weighted (psi ~ Omega_Lambda(z), lives at low z where R is huge -> Pi->1)
OmL=P.OMEGA_LAMBDA
def psi_w(z): return OmL/(OmL+P.OMEGA_M*(1+z)**3)
num2=quad(lambda z: Pi(z)*psi_w(z), 0, 10, limit=300)[0]; den2=quad(psi_w,0,10,limit=300)[0]
gate_psi = num2/den2

print("="*82)
print("FREEZEOUT-WEIGHTING: which epoch sets the trace gate Pi_b = R/(1+R)?")
print("="*82)
print(f"\n  DATA-REQUIRED effective gate  = A_data/(2 sqrt alpha) = {A_DATA}/{2*np.sqrt(ALPHA):.4f} = {gate_data:.4f}")
print(f"\n  CANDIDATE WEIGHTINGS:")
print(f"    (a) freezeout/visibility @ z_drag={ZDRAG:.0f}   Pi_b = {gate_vis(ZDRAG):.4f}   "
      f"[instant Pi_b(z_d)={Pi(ZDRAG):.4f}]   <-- the referee's freezeout claim")
print(f"    (a') visibility @ z_star={ZSTAR:.0f}          Pi_b = {gate_vis(ZSTAR):.4f}   "
      f"[instant Pi_b(z*)={Pi(ZSTAR):.4f}]")
print(f"    (b) full sound-horizon integral           <Pi_b> = {gate_int:.4f}   (pre-drag propagation average)")
print(f"    (c) late-time-screen weighted             <Pi_b> = {gate_psi:.4f}   (where psi actually lives, R huge)")
print(f"\n  MATCH TO DATA ({gate_data:.3f}):")
for lab,val in [("(a) freezeout/visibility z_d",gate_vis(ZDRAG)),("(a') visibility z*",gate_vis(ZSTAR)),
                ("(b) sound-horizon integral",gate_int),("(c) late-time-screen",gate_psi)]:
    dev = (val-gate_data)/gate_data*100
    verdict = "MATCHES" if abs(dev)<8 else "REJECTED"
    print(f"    {lab:30s} {val:.4f}   {dev:+6.1f}%   {verdict}")

print("\n" + "-"*82)
print("READ: if (a)/(a') [freezeout surface] MATCH and (b)/(c) are REJECTED, the data discriminate among")
print("weightings and SELECT the freezeout/visibility surface -- the physically-motivated answer (the acoustic")
print("feature imprints at last scattering, weighted by the visibility function, NOT averaged over propagation).")
print("This UPGRADES the freezeout lemma from 'data-backed' to 'data-SELECTED among candidates + physically")
print("motivated by standard visibility weighting'. HONEST RESIDUAL: a from-the-action proof that the trace")
print("coupling's effect on the observable is visibility-weighted (not the c_s integral) is still the open map.")
