"""
EXTERNAL SIGN CHECK (the referee point #4): does my SDSS galaxy-density kappa-proxy have the correct sign
convention? It MUST correlate POSITIVELY with the Planck CMB-lensing kappa (the galaxy x CMB-lensing
cross-correlation is a robustly-measured POSITIVE signal: galaxies trace mass, CMB lensing traces mass).
If r(kappa_proxy_SDSS, kappa_Planck) > 0, then positive proxy = real overdensity/convergence, so a NEGATIVE
SN-residual slope genuinely means the LCDM/standard-lensing sign (magnification). If it came out negative,
my whole sign interpretation would be inverted -- this is the honesty-guard the injection test could not provide.
"""
import numpy as np, healpy as hp, glob, os
from scipy.integrate import quad
import dfd_params as P, background as B

NSIDE=256; d2r=np.pi/180; FD="foreground_data"
def chi(z): return (P.C_KMS/P.H0)*quad(lambda x:1/B.E(x),0,z,limit=60)[0]

# SDSS galaxy-density convergence proxy (CMB-source-weighted: W ~ chi_l(1+z_l))
kproxy=np.zeros(hp.nside2npix(NSIDE)); foot=None
for f in sorted(glob.glob(os.path.join(FD,"sdss_z*.csv"))):
    g=np.loadtxt(f,delimiter=',',ndmin=2)
    if g.size==0: continue
    ra,dec,z=g[:,0],g[:,1],g[:,2]; zc=z.mean()
    cnt=np.bincount(hp.ang2pix(NSIDE,(90-dec)*d2r,ra*d2r),minlength=hp.nside2npix(NSIDE)).astype(float)
    nfg=hp.smoothing(cnt,fwhm=0.5*d2r); nbg=hp.smoothing(cnt,fwhm=10*d2r)
    dg=np.where(nbg>0.05,nfg/np.maximum(nbg,1e-6)-1.0,0.0)
    kproxy += chi(zc)*(1+zc)*dg
    cov=nbg>0.05; foot=cov if foot is None else (foot&cov)

# Planck CMB lensing kappa map (GALACTIC coords) -> rotate to EQUATORIAL to match SDSS (RA/Dec)
klm=hp.read_alm("lensing_map/COM_Lensing_4096_R3.00/MV/dat_klm.fits")
lm=hp.Alm.getlmax(len(klm)); klm=hp.almxfl(klm,(np.arange(lm+1)<=400).astype(float))
kcmb_gal=hp.smoothing(hp.alm2map(klm,NSIDE),fwhm=0.5*d2r)
kcmb=hp.Rotator(coord=['G','C']).rotate_map_pixel(kcmb_gal)   # Galactic -> Celestial/Equatorial
print("  [coordinate fix] rotated Planck kappa Galactic->Equatorial to match SDSS RA/Dec")

# smooth proxy to the same scale, correlate over the SDSS footprint
kproxy_s=hp.smoothing(kproxy,fwhm=0.5*d2r)
m=foot & np.isfinite(kproxy_s) & np.isfinite(kcmb)
x=kproxy_s[m]-kproxy_s[m].mean(); y=kcmb[m]-kcmb[m].mean()
r=np.sum(x*y)/np.sqrt(np.sum(x*x)*np.sum(y*y))
slope=np.sum(x*y)/np.sum(x*x)
print("="*72); print("EXTERNAL SIGN CHECK — SDSS kappa-proxy vs Planck CMB-lensing kappa"); print("="*72)
print(f"\n  footprint pixels: {m.sum()}  (~{100*m.sum()/len(foot):.0f}% sky)")
print(f"  Pearson r(kappa_proxy_SDSS, kappa_Planck) = {r:+.4f}")
print(f"  regression slope = {slope:+.4e}")
print(f"\n  {'VERDICT':8s}: " + ("POSITIVE -> sign convention CORRECT. My proxy's positive = real convergence;"
      if r>0 else "NEGATIVE -> sign convention INVERTED; the SN-residual interpretation must flip!"))
if r>0:
    print("  so the SN-residual NEGATIVE slope genuinely = LCDM/standard-lensing sign (magnification).")
    print("  This is the known galaxy x CMB-lensing positive cross-correlation, recovered -> proxy validated.")
