"""
CALIBRATE THE MEASURING STICK: how well does the BOSS galaxy-density proxy-kappa trace TRUE mass?
Cross-correlate the BOSS proxy-kappa MAP with the real Planck CMB-lensing kappa MAP over the BOSS footprint.
The correlation r quantifies the proxy's signal-vs-noise as a mass tracer = how SHARP the knife is.
HONEST CAVEAT: Planck kappa is CMB lensing (source z~1100), so its kernel only partly overlaps the foreground
(z<0.7) the proxy traces -> r(proxy, CMB) is a LOWER BOUND on r(proxy, true foreground kappa). A real shear
kappa map (foreground kernel) would give the exact calibration; that awaits DECADE/LSST/Roman.
"""
import numpy as np, healpy as hp, glob, os
from astropy.io import fits
from scipy.integrate import quad
import dfd_params as P, background as B

BD="boss_data"; NSIDE=256; SMOOTH=1.0; d2r=np.pi/180
def chi(z): return (P.C_KMS/P.H0)*quad(lambda x:1/B.E(x),0,max(z,1e-4),limit=60)[0]
def load_gal(f):
    d=fits.open(f, ignore_missing_end=True)[1].data
    ra=np.asarray(d['RA'],float); dec=np.asarray(d['DEC'],float); z=np.asarray(d['Z'],float)
    try: w=np.asarray(d['WEIGHT_SYSTOT'],float)*(np.asarray(d['WEIGHT_CP'],float)+np.asarray(d['WEIGHT_NOZ'],float)-1.0)
    except Exception: w=np.ones(len(ra))
    return ra,dec,z,w
def load_rand(f):
    with open(f,'rb') as fh: raw=bytearray(fh.read())
    h=fits.open(f, ignore_missing_end=True); off=int(h[1]._data_offset); nb1=int(h[1].header['NAXIS1']); h.close()
    nrow=(len(raw)-off)//nb1; idx=raw.find(b'NAXIS2'); raw[idx:idx+80]=("NAXIS2  = %20d"%nrow).ljust(80).encode()
    tmp=f+'.pp'; open(tmp,'wb').write(bytes(raw[:off+nrow*nb1])); d=fits.open(tmp,ignore_missing_end=True)[1].data
    return np.asarray(d['RA'],float),np.asarray(d['DEC'],float),np.asarray(d['Z'],float)

# build BOSS proxy-kappa map (galaxy-density convergence, CMB-source-weighted W~chi_l(1+z_l), bias-corrected)
SAMP=[(t,b) for t,b in [("CMASS_N",2.0),("LOWZ_N",1.85),("CMASS_S",2.0),("LOWZ_S",1.85)]
      if os.path.exists(f"{BD}/galaxy_{t}.fits.gz")]
npix=hp.nside2npix(NSIDE); kproxy=np.zeros(npix); foot=np.zeros(npix,bool)
ZB=np.arange(0.10,0.75,0.05)
GAL={t:load_gal(f"{BD}/galaxy_{t}.fits.gz") for t,_ in SAMP}; RND={t:load_rand(f"{BD}/rand_{t}.fits") for t,_ in SAMP}
for i in range(len(ZB)-1):
    zlo,zhi=ZB[i],ZB[i+1]; zc=0.5*(zlo+zhi); Nd=np.zeros(npix); Nr=np.zeros(npix); bz=2.0
    for t,bias in SAMP:
        gra,gdec,gz,gw=GAL[t]; rra,rdec,rz=RND[t]; gs=(gz>=zlo)&(gz<zhi); rs=(rz>=zlo)&(rz<zhi)
        if gs.sum()==0 or rs.sum()==0: continue
        np.add.at(Nd,hp.ang2pix(NSIDE,(90-gdec[gs])*d2r,gra[gs]*d2r),gw[gs])
        np.add.at(Nr,hp.ang2pix(NSIDE,(90-rdec[rs])*d2r,rra[rs]*d2r),1.0); bz=bias
    if Nr.sum()<100: continue
    alpha=Nd.sum()/Nr.sum(); cov=Nr>0; foot|=cov
    dg=np.zeros(npix); dg[cov]=(Nd[cov]-alpha*Nr[cov])/(alpha*Nr[cov])
    kproxy += chi(zc)*(1+zc)*(dg/bz)
kproxy=hp.smoothing(kproxy,fwhm=SMOOTH*d2r)

# Planck CMB-lensing kappa (Galactic -> Equatorial)
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=hp.Rotator(coord=['G','C']).rotate_map_pixel(hp.smoothing(hp.alm2map(klm,NSIDE),fwhm=SMOOTH*d2r))

m=foot&np.isfinite(kproxy)&np.isfinite(kcmb)
x=kproxy[m]-kproxy[m].mean(); y=kcmb[m]-kcmb[m].mean()
r=np.sum(x*y)/np.sqrt(np.sum(x*x)*np.sum(y*y))
print("="*78); print("MEASURING-STICK CALIBRATION — BOSS proxy-kappa x Planck CMB-lensing kappa"); print("="*78)
print(f"\n  BOSS footprint pixels: {m.sum()} (~{100*m.sum()/npix:.0f}% sky)")
print(f"  Pearson r(kappa_proxy_BOSS, kappa_Planck) = {r:+.3f}")
print(f"  (SDSS proxy gave +0.032; higher r = the proper-randoms+bias proxy traces mass better.)")
print(f"\n  INTERPRETATION (how sharp is the knife):")
print(f"  r is a LOWER BOUND on the proxy<->true-foreground-kappa correlation (CMB-lensing kernel mismatch).")
print(f"  The BOSS slope is ATTENUATED by ~this correlation: a true physical DFD slope b_phys shows up in")
print(f"  proxy units as ~ b_phys * (attenuation). With r~{r:.2f}, the 648-SN BOSS proxy test has LOW physical")
print(f"  power -> consistent with the null being uninformative on the physical amplitude. A real FOREGROUND")
print(f"  shear-kappa map (DECADE/LSST/Roman) would raise r and sharpen the knife. This is the honest bound now.")
