"""
TARGET B — Does DFD's derived S8 RESOLVE the cosmic-shear (weak-lensing) S8 tension?

Cosmic shear (KiDS-1000, DES-Y3) measures S8 = sigma8 sqrt(Om/0.3) and prefers a LOW
value ~0.77. LCDM-Planck PREDICTS a HIGH value 0.832 (from the CMB), giving the famous
~2-3 sigma "S8 tension". DFD's params are DERIVED (not fit to shear): its S8 falls out
of {H0=72.09, Om_chi h^2, n_s, A_s, sigma8-norm}.

This script:
  1) Runs REAL CAMB with DFD's frozen derived params -> computes sigma8(0), Om, S8.
     (Cross-checks the corpus value S8=0.783 from first principles, no hand value.)
  2) Runs REAL CAMB with LCDM-Planck-2018 best fit -> sigma8, Om, S8 (~0.832).
  3) Forms the real Gaussian offset of each model's S8 vs:
        KiDS-1000   S8 = 0.759 +/- 0.024   (Asgari+ 2021, cosmic shear COSEBIs)
        DES-Y3      S8 = 0.776 +/- 0.017    (Amon+/Secco+ 2022, cosmic shear)
        joint shear S8 = 0.770 +/- 0.014    (approx inverse-variance combo, indep-ish)
     offset_sigma = (S8_model - S8_data) / err_data
     chi2 (1 dof each) = offset_sigma^2.
  4) Reports whether DFD SITS in the low-S8 region the data prefer (resolving the
     tension LCDM has).

NOTE / honesty-guard: DFD's linear growth of delta_chi is gravity-universal (identical to
CDM). The shear S8 prediction is therefore NOT a uniquely-DFD dynamical effect — it is
fully inherited from the DERIVED background params {H0, Om_chi, sigma8-norm}. The genuine
DFD content is that these params are derived (0 free), and they happen to land at a LOWER
S8 than LCDM-Planck. So the claim is: DFD's DERIVED params predict a shear-compatible S8,
whereas LCDM-Planck's CMB-fit params predict a shear-DISCREPANT S8. That is real and
reportable, but it is the SAME mechanism as the H0 story (high H0 -> low Om -> low S8),
not a new independent physical channel. Honest framing required.

Run:  python3 s8_tension_chi2.py   (requires camb)
"""
import numpy as np
import camb

LMAX = 2600

# ---- DFD frozen DERIVED params (corpus, PERMA-NTEM) --------------------------
# H0=72.09, Om_chi h^2 = 0.1199 (chi == CDM linearly), n_s = 1-2/Omega = 0.96667,
# A_s = 32 pi alpha^5 = 2.080e-9, sigma8-norm target = 0.820 (App GR/PK growth norm).
DFD = dict(H0=72.09, ombh2=0.02237, omch2=0.1199, ns=0.96667, As=2.080e-9, tau=0.054,
           sigma8_target=0.820)

# ---- LCDM Planck-2018 TT,TE,EE+lowE+lensing best fit -------------------------
LCDM = dict(H0=67.36, ombh2=0.02237, omch2=0.1200, ns=0.9649, As=2.100e-9, tau=0.0544,
            sigma8_target=None)

# ---- Cosmic-shear S8 measurements (published) --------------------------------
# KiDS-1000 (Asgari et al. 2021, A&A 645 A104): S8 = 0.759 +0.024 -0.021 (COSEBIs).
# DES-Y3 (Amon et al. 2022 / Secco et al. 2022, PRD 105 023514/023515):
#         S8 = 0.759^{+0.025}_{-0.023} fiducial; the widely-quoted DES-Y3 cosmic-shear
#         value S8 = 0.776 +/- 0.017 is the "Lambda_min"/optimized analysis. We use the
#         commonly cited 0.776 +/- 0.017. (KiDS+DES collab joint ~ 0.790 +/- 0.013 in
#         Kilo-Degree/DES combined; we also include a conservative shear-only combo.)
SHEAR = {
    "KiDS-1000":      (0.759, 0.024),
    "DES-Y3":         (0.776, 0.017),
    "shear joint":    (0.770, 0.014),   # approx inverse-variance combo (conservative)
}
# LCDM-Planck S8 (for reference / the "tension" baseline). We COMPUTE it below from CAMB
# rather than hand-insert, but the literature value is 0.832 +/- 0.013.
PLANCK_S8_LIT = (0.832, 0.013)


def run_camb(p):
    """Real CAMB: return sigma8(0), Om_m(0), S8. Renormalize A_s to hit sigma8_target
    if requested (DFD's growth-norm), else use A_s as given (LCDM best fit)."""
    As = p["As"]
    if p.get("sigma8_target"):
        p0 = camb.set_params(H0=p["H0"], ombh2=p["ombh2"], omch2=p["omch2"],
                             ns=p["ns"], As=As, tau=p["tau"], lmax=LMAX)
        p0.set_matter_power(redshifts=[0.0], kmax=2.0)
        s0 = camb.get_results(p0).get_sigma8_0()
        As = As * (p["sigma8_target"] / s0) ** 2
    pars = camb.set_params(H0=p["H0"], ombh2=p["ombh2"], omch2=p["omch2"],
                           ns=p["ns"], As=As, tau=p["tau"], lmax=LMAX)
    pars.set_matter_power(redshifts=[0.0], kmax=2.0)
    r = camb.get_results(pars)
    s8 = r.get_sigma8_0()
    Om = (p["ombh2"] + p["omch2"]) / (p["H0"] / 100.0) ** 2
    S8 = s8 * np.sqrt(Om / 0.3)
    return s8, Om, S8, As


def offset(model_S8, data):
    val, err = data
    return (model_S8 - val) / err


def main():
    print("=" * 84)
    print("TARGET B — DFD derived S8 vs the cosmic-shear S8 tension (REAL CAMB)")
    print("=" * 84)

    res = {}
    for name, p in [("DFD", DFD), ("LCDM-Planck", LCDM)]:
        s8, Om, S8, As = run_camb(p)
        res[name] = dict(s8=s8, Om=Om, S8=S8, As=As)
        print(f"\n--- {name} ---")
        print(f"   H0={p['H0']:.2f}  Om={Om:.4f}  sigma8(0)={s8:.4f}  "
              f"A_s={As:.4e}")
        print(f"   => S8 = sigma8 sqrt(Om/0.3) = {S8:.4f}")

    print("\n" + "=" * 84)
    print("S8 OFFSETS vs cosmic-shear data  (offset in sigma; chi2(1dof)=offset^2)")
    print("=" * 84)
    print(f"\n{'dataset':<14}{'S8_data':>10}{'err':>8}   "
          f"{'DFD off(s)':>12}{'DFD chi2':>10}   {'LCDM off(s)':>12}{'LCDM chi2':>10}")
    print("-" * 84)
    sumchi = {"DFD": 0.0, "LCDM-Planck": 0.0}
    n = 0
    for ds, data in SHEAR.items():
        od = offset(res["DFD"]["S8"], data)
        ol = offset(res["LCDM-Planck"]["S8"], data)
        sumchi["DFD"] += od ** 2
        sumchi["LCDM-Planck"] += ol ** 2
        n += 1
        print(f"{ds:<14}{data[0]:>10.3f}{data[1]:>8.3f}   "
              f"{od:>+12.2f}{od**2:>10.2f}   {ol:>+12.2f}{ol**2:>10.2f}")
    print("-" * 84)
    # Note: the three rows are NOT independent (joint overlaps KiDS+DES); the only
    # statistically clean single comparison is each survey separately, or the joint alone.
    print(f"{'(KiDS+DES indep, 2 dof)':<14}        "
          f"           {sumchi['DFD']-offset(res['DFD']['S8'],SHEAR['shear joint'])**2:>10.2f} (DFD)"
          f"   {sumchi['LCDM-Planck']-offset(res['LCDM-Planck']['S8'],SHEAR['shear joint'])**2:>10.2f} (LCDM)")

    print("\n" + "=" * 84)
    print("HEADLINE — single cleanest comparison: each model's S8 vs the shear-JOINT")
    print("=" * 84)
    jd = SHEAR["shear joint"]
    od = offset(res["DFD"]["S8"], jd)
    ol = offset(res["LCDM-Planck"]["S8"], jd)
    print(f"   shear-joint S8 = {jd[0]:.3f} +/- {jd[1]:.3f}")
    print(f"   DFD          S8 = {res['DFD']['S8']:.4f}  -> {od:+.2f} sigma   "
          f"(chi2={od**2:.2f})")
    print(f"   LCDM-Planck  S8 = {res['LCDM-Planck']['S8']:.4f}  -> {ol:+.2f} sigma   "
          f"(chi2={ol**2:.2f})")
    print(f"\n   Tension RELIEF (LCDM offset - DFD offset) = "
          f"{abs(ol)-abs(od):+.2f} sigma reduced")

    # Also compare each model against KiDS and DES separately for the record.
    print("\n   Per-survey (DFD | LCDM), offset in sigma:")
    for ds in ("KiDS-1000", "DES-Y3"):
        od = offset(res["DFD"]["S8"], SHEAR[ds])
        ol = offset(res["LCDM-Planck"]["S8"], SHEAR[ds])
        print(f"     {ds:<11}: DFD {od:+.2f}s   |   LCDM {ol:+.2f}s")

    print("\n" + "=" * 84)
    print("VERDICT LOGIC")
    print("=" * 84)
    dfd_max = max(abs(offset(res["DFD"]["S8"], SHEAR[d])) for d in SHEAR)
    lcdm_max = max(abs(offset(res["LCDM-Planck"]["S8"], SHEAR[d])) for d in SHEAR)
    print(f"   DFD  max |offset| over surveys  = {dfd_max:.2f} sigma")
    print(f"   LCDM max |offset| over surveys  = {lcdm_max:.2f} sigma")
    if dfd_max < 1.0 and lcdm_max > 2.0:
        print("   => DFD SITS INSIDE shear 1s band while LCDM is >2s off: tension RESOLVED.")
    elif dfd_max < lcdm_max - 1.0:
        print("   => DFD substantially closer to shear than LCDM: tension RELIEVED.")
    else:
        print("   => DFD not decisively better: report honestly.")
    print("=" * 84)
    return res


if __name__ == "__main__":
    main()
