"""
HALO-LOS AUDIT (the referee's two checks): (1) show the FULL x-distribution behind <1/mu> (Jensen gap, tail,
outer-cutoff sensitivity); (2) define EXACTLY what is predicted vs which target. No prose, one table.
Predicted object = galaxy-weighted screen EXCESS <Dpsi>_gal = kappa_gal * <1/mu>  (anisotropic, galaxy-selected
sightlines vs mean) -> compares to the BAO DIFFERENTIAL (~0.030), NOT the absolute Dpsi(0.5)~0.039.
"""
import numpy as np
from scipy.integrate import quad
import camb, dfd_params as P

G=6.674e-11; Msun=1.989e30; Mpc=3.0857e22; c_si=2.998e8; H0=P.H0*1e3/Mpc
a0=2*np.sqrt(P.ALPHA)*c_si*H0
z=0.5; rho_crit_z=3*(H0*np.sqrt(P.OMEGA_M*(1+z)**3+P.OMEGA_LAMBDA))**2/(8*np.pi*G)
rho_m_z=P.OMEGA_M*(1+z)**3*3*H0**2/(8*np.pi*G)
pars=camb.CAMBparams(); pars.set_cosmology(H0=P.H0,ombh2=P.OMEGA_B_H2,omch2=P.OMEGA_CHI_H2,mnu=0.06,tau=0.054)
pars.InitPower.set_params(As=P.A_S,ns=P.N_S); pars.set_matter_power(redshifts=[z],kmax=80.0)
camb.get_results(pars); PK=camb.get_matter_power_interpolator(pars,nonlinear=False,kmax=80,zmax=1)
def sigma_M(Mm):
    R=(3*Mm*Msun/(4*np.pi*rho_m_z))**(1/3)/Mpc
    def ig(k):
        kR=k*R; W=3*(np.sin(kR)-kR*np.cos(kR))/kR**3 if kR>1e-3 else 1.0
        return PK.P(z,k)*W**2*k**2/(2*np.pi**2)
    return np.sqrt(quad(ig,1e-4,80,limit=200)[0])
def st_f(nu): A,a,p=0.3222,0.707,0.3; return A*np.sqrt(2*a/np.pi)*(1+(1/(a*nu**2))**p)*nu*np.exp(-a*nu**2/2)
def Rvir(M): return (3*M/(4*np.pi*200*rho_crit_z))**(1/3)
def conc(Mm): return 5.0*(Mm/1e14)**(-0.1)

def halo_xw(Mm, rout_mult):                      # return (x_array, weight_array) for one halo, profile-weighted
    M=Mm*Msun; Rv=Rvir(M)*rout_mult; cc=conc(Mm); rs=Rvir(M)/cc; mc=np.log(1+cc)-cc/(1+cc)
    rr=np.logspace(np.log10(rs*1e-2),np.log10(Rv),400); xn=rr/rs
    Mr=M*(np.log(1+xn)-xn/(1+xn))/mc; xacc=G*Mr/rr**2/a0
    rho=1.0/(xn*(1+xn)**2); w=rho*rr**3            # rho r^2 dr in dlnr
    return xacc, w

masses=np.logspace(11,15.3,26); sig=np.array([sigma_M(m) for m in masses]); nu=1.686/sig
dndlnM=(rho_m_z/(masses*Msun))*st_f(nu)*np.abs(np.gradient(np.log(sig),np.log(masses)))
Ngal=np.where(masses>1e12,(masses/1e12)**0.9,0.0)

def population(rout_mult, halo_w):                # halo_w: 'mass' or 'gal'
    X=[]; W=[]
    for i,Mm in enumerate(masses):
        xacc,w=halo_xw(Mm,rout_mult)
        hw = dndlnM[i]*masses[i] if halo_w=='mass' else dndlnM[i]*Ngal[i]
        X.append(xacc); W.append(w*hw)
    X=np.concatenate(X); W=np.concatenate(W); W/=W.sum()
    o=np.argsort(X); Xs=X[o]; Ws=W[o]; cw=np.cumsum(Ws)
    med=Xs[np.searchsorted(cw,0.5)]
    mean=np.sum(X*W); invmu=np.sum(((1+X)/X)*W); frac_lt01=np.sum(W[X<0.1])
    return dict(mean=mean,med=med,Q=invmu,frac01=frac_lt01)

KGAL=0.006; BAO_DIFF=0.030; ABS_DPSI05=0.039
print("="*84); print("HALO-LOS AUDIT"); print("="*84)
print(f"a0={a0:.3e}.  mu(x)=x/(1+x).  kappa_gal baseline={KGAL}.")
print("\nJENSEN / TAIL AUDIT (why Q>1/mu(<x>)):  1/mu(<x>) vs <1/mu>, and the low-x tail")
print(f"{'weighting':10s} {'<x>':>7s} {'med x':>7s} {'1/mu(<x>)':>10s} {'<1/mu>':>8s} {'frac x<0.1':>11s}")
for w in ('mass','gal'):
    r=population(1.0,w)
    print(f"{w:10s} {r['mean']:7.3f} {r['med']:7.3f} {(1+r['mean'])/r['mean']:10.2f} {r['Q']:8.2f} {r['frac01']:11.2f}")
print("\nOUTER-CUTOFF SENSITIVITY (is the tail an artifact of how far the halo extends?) [mass-weighted]")
print(f"{'r_out':>8s} {'<x>':>7s} {'<1/mu>':>8s} {'pred excess=kgal*Q':>19s}")
for rom in (0.5,1.0,1.5,2.0):
    r=population(rom,'mass'); print(f"{rom:7.1f}Rv {r['mean']:7.3f} {r['Q']:8.2f} {KGAL*r['Q']:19.4f}")

print("\n" + "="*84); print("THE TABLE (mass-weighted, r_out=1.0 Rvir)"); print("="*84)
r=population(1.0,'mass'); rg=population(1.0,'gal')
rows=[("halo mass range","1e11 - 2e15 Msun"),
      ("weighting used","mass (trace~rho)  |  galaxy (HOD)"),
      ("<x>",f"{r['mean']:.3f}  |  {rg['mean']:.3f}"),
      ("median x",f"{r['med']:.3f}  |  {rg['med']:.3f}"),
      ("<1/mu>",f"{r['Q']:.2f}  |  {rg['Q']:.2f}"),
      ("predicted EXCESS = kgal*<1/mu>",f"{KGAL*r['Q']:.4f}  |  {KGAL*rg['Q']:.4f}"),
      ("target BAO differential",f"{BAO_DIFF:.3f}"),
      ("predicted absolute Dpsi(0.5)","N/A (model gives the EXCESS, not the monopole)"),
      ("target absolute Dpsi(0.5)",f"{ABS_DPSI05:.3f}  (different object; do NOT compare to excess)")]
for k,v in rows: print(f"  {k:34s} {v}")
print("\nNOTE: the predicted object is the galaxy-weighted EXCESS -> compare to BAO DIFFERENTIAL 0.030 ONLY.")
print("The Jensen gap (Q>>1/mu(<x>)) is driven by the low-x halo-outskirt tail; cutoff sensitivity above")
print("shows how much of Q is robust vs an artifact of the assumed outer radius.")
