diff --git a/papers/bmodes/rules/claims.smk b/papers/bmodes/rules/claims.smk
index 43962627..3846b9bc 100644
--- a/papers/bmodes/rules/claims.smk
+++ b/papers/bmodes/rules/claims.smk
@@ -355,25 +355,21 @@ rule pure_eb_covariance:
 
 
 rule calculate_pure_eb_ptes:
-    """Calculate PTE matrices for Pure E/B mode scale cut robustness.
+    """PTE matrices for pure E/B-mode scale-cut robustness.
 
-    Per-blind: Uses blind-specific integration covariance for PTE calculation.
-    The pure_eb_data vectors are identical across blinds; only covariance differs.
-
-    In practice, BB covariance is blind-independent (validated by
-    bb_covariance_blind_independence), so downstream consumers (config_space_pte_matrices)
-    only request blind A. The per-blind wildcard is retained for the blind independence test.
+    Nothing here varies with the blind: the data vectors come from the blind-A
+    gather and the PTEs are Hartlap-debiased by the MC draw count, not by a
+    per-blind covariance. The wildcard survives as the filename slot the
+    consumer (config_space_pte_matrices) reads, and only blind A is ever built.
     """
     input:
         pure_eb_data="results/paper_plots/intermediate/{version}_A_pure_eb_semianalytic.npz",
-        cov_integration=lambda w: _cov_integration_path(w.version, w.blind),
     output:
         "results/paper_plots/intermediate/{version}_{blind}_pure_eb_ptes.npz",
     wildcard_constraints:
         blind=r"[ABC]",
     params:
         version="{version}",
-        npatch=FIDUCIAL["npatch"],
         n_samples=config["covariance"]["n_samples"],
     resources:
         mem_mb=16000,
diff --git a/papers/bmodes/scripts/calculate_pure_eb_ptes.py b/papers/bmodes/scripts/calculate_pure_eb_ptes.py
index f8bc709d..3a623669 100644
--- a/papers/bmodes/scripts/calculate_pure_eb_ptes.py
+++ b/papers/bmodes/scripts/calculate_pure_eb_ptes.py
@@ -4,14 +4,13 @@
 pure-E/B ``semianalytic.npz`` (data vectors + MC covariance), evaluates the
 ξ_+^B / ξ_-^B / joint ξ_tot^B χ² PTE matrices over the scale-cut grid via
 ``sp_validation.b_modes.calculate_eb_statistics`` (Hartlap-corrected inverse
-MC covariance), and writes the PTE matrices to
+MC covariance, debiased by the draw count), and writes the PTE matrices to
 ``{out}/{version}_{blind}_pure_eb_ptes.npz``.
 
     python calculate_pure_eb_ptes.py \
         --version SP_v1.4.6.3_leak_corr --blind A \
         --pure-eb-data <..._pure_eb_semianalytic.npz> \
-        --cov-integration <cov ..._processed.txt> \
-        --npatch 1 --n-samples 2000 --out <output_dir>
+        --n-samples 2000 --out <output_dir>
 """
 
 import argparse
@@ -22,31 +21,21 @@
 from sp_validation.b_modes import calculate_eb_statistics
 
 
-class FakeGG:
-    """Minimal GGCorrelation-like object for calculate_eb_statistics."""
-
-    def __init__(self, nbins, npatch):
-        self.nbins = nbins
-        self.npatch1 = npatch
-        self.npatch2 = npatch
-
-
 def calculate_ptes(
     version,
     blind,
     pure_eb_data,
-    cov_integration,
-    npatch,
     n_samples,
     output_dir,
 ):
     dataset = np.load(pure_eb_data)
 
     theta = dataset["theta"]
-    nbins = len(theta)
 
     results = {
-        "gg": FakeGG(nbins, int(npatch)),
+        "theta": theta,
+        # The MC draws are the realisations behind this covariance.
+        "n_eff": int(n_samples),
         "xip_E": dataset["xip_E"],
         "xim_E": dataset["xim_E"],
         "xip_B": dataset["xip_B"],
@@ -57,11 +46,7 @@ def calculate_ptes(
     }
 
     print(f"Calculating PTE matrices for {version}...")
-    results = calculate_eb_statistics(
-        results,
-        cov_path_int=cov_integration,
-        n_samples=int(n_samples),
-    )
+    results = calculate_eb_statistics(results)
 
     pte_matrices = results["pte_matrices"]
     output_data = {
@@ -83,12 +68,6 @@ def _from_cli(argv=None):
     ap.add_argument("--version", required=True)
     ap.add_argument("--blind", default="A")
     ap.add_argument("--pure-eb-data", required=True, help="Gathered semianalytic .npz")
-    ap.add_argument(
-        "--cov-integration",
-        default=None,
-        help="Integration-grid covariance _processed.txt (optional)",
-    )
-    ap.add_argument("--npatch", type=int, default=1)
     ap.add_argument("--n-samples", type=int, default=2000)
     ap.add_argument("--out", required=True, help="Output directory (lc {output})")
     a = ap.parse_args(argv)
@@ -96,8 +75,6 @@ def _from_cli(argv=None):
         version=a.version,
         blind=a.blind,
         pure_eb_data=a.pure_eb_data,
-        cov_integration=a.cov_integration,
-        npatch=a.npatch,
         n_samples=a.n_samples,
         output_dir=a.out,
     )
diff --git a/papers/bmodes/scripts/precompute_pure_eb_chunk.py b/papers/bmodes/scripts/precompute_pure_eb_chunk.py
index 363a3192..99153096 100644
--- a/papers/bmodes/scripts/precompute_pure_eb_chunk.py
+++ b/papers/bmodes/scripts/precompute_pure_eb_chunk.py
@@ -127,15 +127,15 @@ def compute_chunk(
     row_sums = np.array(binning_matrix.sum(axis=1)).flatten()
     binning_matrix = sparse.diags(1 / row_sums) @ binning_matrix
 
-    mean_int = np.concatenate(
-        get_theo_xi(
-            theta=theta_int,
-            z=z_dist[:, 0],
-            nz=z_dist[:, 1],
-            backend="ccl",
-            cosmo=cosmo_cov,
-        )
-    )
+    # One n(z) gives one tracer pair: get_theo_xi's single (xi+, xi-) entry.
+    (xi_pm,) = get_theo_xi(
+        theta=theta_int,
+        z=z_dist[:, 0],
+        nz=z_dist[:, 1],
+        backend="ccl",
+        cosmo=cosmo_cov,
+    ).values()
+    mean_int = np.concatenate(xi_pm)
 
     rng = np.random.default_rng(seed=42 + chunk_id)
 
diff --git a/papers/bmodes/scripts/run_pure_eb_ptes_sweep.sh b/papers/bmodes/scripts/run_pure_eb_ptes_sweep.sh
index 0ea3a862..0e7ecdf5 100644
--- a/papers/bmodes/scripts/run_pure_eb_ptes_sweep.sh
+++ b/papers/bmodes/scripts/run_pure_eb_ptes_sweep.sh
@@ -4,30 +4,27 @@
 # Loops the non-fiducial version list (resolved via sweep_versions.py) and runs
 # calculate_pure_eb_ptes.py once per version — the same χ² PTE-matrix compute the
 # fiducial pure_eb_pte_per_cut recipe calls. Each version's gathered semi-analytic
-# NPZ is read from the pure_eb_sweep output dir and its Gaussian integration
-# covariance from the cov_sweep output dir (both by absolute path — lc does not
-# wire cross-output deps, so the driver reads upstream sweeps directly; run
-# pure_eb_sweep + cov_sweep before this). Each version emits the canonical
-# ``{ver}_{blind}_pure_eb_ptes.npz`` — the exact name config_space_pte_matrices.py
-# reconstructs from --pte-intermediate-dir — straight into --out. Serial over
-# versions; each version is fast (206-pair grid, ~seconds).
+# NPZ (data vectors + MC covariance) is read from the pure_eb_sweep output dir by
+# absolute path — lc does not wire cross-output deps, so the driver reads the
+# upstream sweep directly; run pure_eb_sweep before this. Each version emits the
+# canonical ``{ver}_{blind}_pure_eb_ptes.npz`` — the exact name
+# config_space_pte_matrices.py reconstructs from --pte-intermediate-dir — straight
+# into --out. Serial over versions; each version is fast (206-pair grid, ~seconds).
 #
 # Usage:
 #   run_pure_eb_ptes_sweep.sh --config <config.yaml> --cat-config <cat_config.yaml> \
 #     --pure-eb-sweep-dir <pure_eb/.../pure_eb_sweep> \
-#     --cov-sweep-dir <covariance/.../cov_sweep> \
 #     --out <output_dir> [--blind A] [--versions "v1 v2 ..."]
 set -euo pipefail
 
 . "$(dirname "${BASH_SOURCE[0]}")/container_env.sh"
 
-CONFIG=""; CATCONFIG=""; PUREEBSWEEP=""; COVSWEEP=""; OUT=""; BLIND="A"; VERSIONS=""
+CONFIG=""; CATCONFIG=""; PUREEBSWEEP=""; OUT=""; BLIND="A"; VERSIONS=""
 while [ $# -gt 0 ]; do
   case "$1" in
     --config) CONFIG="$2"; shift 2;;
     --cat-config) CATCONFIG="$2"; shift 2;;
     --pure-eb-sweep-dir) PUREEBSWEEP="$2"; shift 2;;
-    --cov-sweep-dir) COVSWEEP="$2"; shift 2;;
     --out) OUT="$2"; shift 2;;
     --blind) BLIND="$2"; shift 2;;
     --versions) VERSIONS="$2"; shift 2;;
@@ -41,16 +38,11 @@ VERSIONS=$(sweep_versions "$CONFIG")
 
 for ver in $VERSIONS; do
   pureeb="$PUREEBSWEEP/${ver}_${BLIND}_pure_eb_semianalytic.npz"
-  covbase="covariance_${ver}_${BLIND}_g_minsep=0.5_maxsep=300.0_nbins=1000_masked"
-  covint="$COVSWEEP/$covbase/${covbase}_processed.txt"
-  for f in "$pureeb" "$covint"; do
-    [ -f "$f" ] || { echo "MISSING upstream input for $ver: $f" >&2; exit 1; }
-  done
+  [ -f "$pureeb" ] || { echo "MISSING upstream input for $ver: $pureeb" >&2; exit 1; }
   echo "[pure_eb_ptes_sweep] $ver"
   spv_python "$PSCRIPTS/calculate_pure_eb_ptes.py" \
     --version "$ver" --blind "$BLIND" \
-    --pure-eb-data "$pureeb" --cov-integration "$covint" \
-    --npatch 1 --n-samples 2000 --out "$OUT"
+    --pure-eb-data "$pureeb" --n-samples 2000 --out "$OUT"
   echo "[pure_eb_ptes_sweep] $ver -> $OUT/${ver}_${BLIND}_pure_eb_ptes.npz"
 done
 echo "[pure_eb_ptes_sweep] done -> $OUT"
diff --git a/src/sp_validation/b_modes.py b/src/sp_validation/b_modes.py
index 13f4c2e9..7a3cd484 100644
--- a/src/sp_validation/b_modes.py
+++ b/src/sp_validation/b_modes.py
@@ -74,23 +74,34 @@ def scale_cut_to_bins(gg, min_scale=None, max_scale=None):
     stop_bin : int
         Last included bin index + 1 (for slicing notation)
     """
-    nbins = len(gg.meanr)
+    return bins_from_edges(gg.left_edges, gg.right_edges, min_scale, max_scale)
 
-    if min_scale is not None:
-        # Conservative: exclude bins whose left edge is below min_scale
-        start_bin = np.searchsorted(gg.left_edges, min_scale, side="left")
-    else:
-        start_bin = 0
-
-    if max_scale is not None:
-        # Conservative: exclude bins whose right edge is above max_scale
-        stop_bin = np.searchsorted(gg.right_edges, max_scale, side="right")
-    else:
-        stop_bin = nbins
 
+def bins_from_edges(left_edges, right_edges, min_scale=None, max_scale=None):
+    """:func:`scale_cut_to_bins` on bare bin edges."""
+    start_bin = (
+        np.searchsorted(left_edges, min_scale, side="left")
+        if min_scale is not None
+        else 0
+    )
+    stop_bin = (
+        np.searchsorted(right_edges, max_scale, side="right")
+        if max_scale is not None
+        else len(right_edges)
+    )
     return start_bin, stop_bin
 
 
+def log_bin_edges(min_sep, max_sep, nbins):
+    """TreeCorr ``Log`` bin edges — the grid a binning defines.
+
+    SACC ξ± parts store bin centres, not edges, so a consumer working from a
+    part reconstructs the edges from the binning it was measured on.
+    """
+    edges = np.geomspace(float(min_sep), float(max_sep), int(nbins) + 1)
+    return edges[:-1], edges[1:]
+
+
 def correlation_from_covariance(covariance):
     """
     Convert covariance matrix to correlation matrix.
@@ -164,75 +175,40 @@ def pure_EB(corrs):
             parallel=True,
         )
 
-    # Initialize results dictionary with basic E/B mode data
-    results = {"gg": gg, "gg_int": gg_int}
+    # The results dict is self-describing: the grids it was measured on travel
+    # with the modes, so every consumer downstream works from values alone.
+    results = {
+        "theta": gg.meanr,
+        "left_edges": gg.left_edges,
+        "right_edges": gg.right_edges,
+        "xip": gg.xip,
+        "xim": gg.xim,
+        "var_xip": gg.varxip,
+        "var_xim": gg.varxim,
+        "theta_int": gg_int.meanr,
+        "xip_int": gg_int.xip,
+        "xim_int": gg_int.xim,
+        "n_eff": n_samples if cov_path_int is not None else gg.npatch1,
+    }
     results.update(dict(zip(_EB_KEYS, pure_EB([gg, gg_int]))))
 
     if cov_path_int is not None:
-        # Use semi-analytical covariance propagation
-        print("Computing semi-analytical covariance for pure E/B modes")
-        if z_dist is None:
-            raise ValueError("z_dist must be provided for semi-analytical covariance")
-        if cosmo_cov is None:
+        if z_dist is None or cosmo_cov is None:
             raise ValueError(
-                "cosmo_cov must be provided for semi-analytical covariance"
-            )
-
-        # Load covariance matrix
-        cov_int = np.loadtxt(cov_path_int)
-
-        # Set up integration binning and pre-compute binning matrix
-        nbins_int, theta_int = len(gg_int.meanr), gg_int.meanr
-        reporting_bin_edges = np.concatenate([gg.left_edges, [gg.right_edges[-1]]])
-        bin_indices = np.digitize(theta_int, reporting_bin_edges) - 1
-
-        valid_mask = (bin_indices >= 0) & (bin_indices < len(gg.meanr))
-        row_indices, col_indices = (bin_indices[valid_mask], np.where(valid_mask)[0])
-
-        binning_matrix = sparse.csr_matrix(
-            (np.ones(len(row_indices)), (row_indices, col_indices)),
-            shape=(len(gg.meanr), nbins_int),
-        )
-        row_sums = np.array(binning_matrix.sum(axis=1)).flatten()
-        binning_matrix = sparse.diags(1 / row_sums) @ binning_matrix
-
-        # Generate theoretical xi+/xi- predictions and sample
-        mean_int = np.concatenate(
-            get_theo_xi(
-                theta=theta_int,
-                z=z_dist[:, 0],
-                nz=z_dist[:, 1],
-                backend="ccl",
-                cosmo=cosmo_cov,
+                "semi-analytical covariance needs both z_dist and cosmo_cov"
             )
+        cov, eb_samples = pure_eb_covariance_mc(
+            theta=gg.meanr,
+            left_edges=gg.left_edges,
+            right_edges=gg.right_edges,
+            theta_int=gg_int.meanr,
+            cov_int=np.loadtxt(cov_path_int),
+            z=z_dist[:, 0],
+            nz=z_dist[:, 1],
+            cosmo=cosmo_cov,
+            n_samples=n_samples,
         )
-
-        samples_int = np.random.multivariate_normal(mean_int, cov_int, size=n_samples)
-        samples_int_xip = samples_int[:, :nbins_int]
-        samples_int_xim = samples_int[:, nbins_int:]
-        samples_rep_xip = (binning_matrix @ samples_int_xip.T).T
-        samples_rep_xim = (binning_matrix @ samples_int_xim.T).T
-
-        transformed_samples = [
-            np.concatenate(
-                get_pure_EB_modes(
-                    theta=gg.meanr,
-                    theta_int=gg_int.meanr,
-                    xip=samples_rep_xip[i],
-                    xim=samples_rep_xim[i],
-                    xip_int=samples_int_xip[i],
-                    xim_int=samples_int_xim[i],
-                    tmin=min_sep,
-                    tmax=max_sep,
-                    parallel=True,
-                )
-            )
-            for i in tqdm.tqdm(range(n_samples), desc="MC samples")
-        ]
-
-        # Store semi-analytical covariance results
-        eb_samples = np.array(transformed_samples)
-        results.update({"cov": np.cov(eb_samples.T), "eb_samples": eb_samples})
+        results.update({"cov": cov, "eb_samples": eb_samples})
     else:
         # Use existing treecorr covariance estimation
         results["cov"] = treecorr.estimate_multi_cov(
@@ -255,6 +231,114 @@ def pure_EB(corrs):
     return results
 
 
+def pure_eb_from_xi(
+    theta_report, xip_report, xim_report, theta_int, xip_int, xim_int, tmin, tmax
+):
+    """Pure-E/B correlation functions from ξ± arrays through the pipeline kernel.
+
+    The values-only seam of :func:`calculate_pure_eb_correlation`, for callers
+    holding ξ± arrays rather than TreeCorr correlations.
+
+    ``tmin``/``tmax`` are the reporting correlation's TreeCorr *bin edges*
+    (``gg.left_edges[0]`` / ``gg.right_edges[-1]``). The reporting grid must be a
+    strict sub-range of the integration grid: a reporting point on the
+    integration boundary has no interior support and comes back NaN.
+
+    Returns
+    -------
+    dict
+        Keyed by ``_EB_KEYS`` (xip_E, xim_E, xip_B, xim_B, xip_amb, xim_amb).
+    """
+    from cosmo_numba.B_modes.schneider2022 import get_pure_EB_modes
+
+    modes = get_pure_EB_modes(
+        theta=np.asarray(theta_report),
+        xip=np.asarray(xip_report),
+        xim=np.asarray(xim_report),
+        theta_int=np.asarray(theta_int),
+        xip_int=np.asarray(xip_int),
+        xim_int=np.asarray(xim_int),
+        tmin=tmin,
+        tmax=tmax,
+        parallel=True,
+    )
+    return dict(zip(_EB_KEYS, (np.asarray(m) for m in modes)))
+
+
+def pure_eb_covariance_mc(
+    *,
+    theta,
+    left_edges,
+    right_edges,
+    theta_int,
+    cov_int,
+    z,
+    nz,
+    cosmo,
+    n_samples=1000,
+):
+    """Pure-E/B covariance by Monte Carlo through the same kernel as the modes.
+
+    ξ± draws come from ``cov_int``, a ξ± covariance on the integration grid,
+    around the theory mean for ``(z, nz)`` under ``cosmo``; each draw is binned
+    down to the reporting grid and pushed through ``get_pure_EB_modes``. The
+    covariance of the transformed draws is the result, so it depends on the
+    covariance model and the grids, never on the measured data vector.
+
+    Returns ``(cov, eb_samples)`` — the covariance in ``_EB_KEYS`` order and
+    the draws behind it.
+    """
+    from cosmo_numba.B_modes.schneider2022 import get_pure_EB_modes
+
+    theta, theta_int = np.asarray(theta), np.asarray(theta_int)
+    nbins_int = len(theta_int)
+
+    # Each reporting bin averages the integration bins that fall inside it.
+    reporting_bin_edges = np.concatenate([left_edges, [right_edges[-1]]])
+    bin_indices = np.digitize(theta_int, reporting_bin_edges) - 1
+    valid_mask = (bin_indices >= 0) & (bin_indices < len(theta))
+    row_indices, col_indices = (bin_indices[valid_mask], np.where(valid_mask)[0])
+    binning_matrix = sparse.csr_matrix(
+        (np.ones(len(row_indices)), (row_indices, col_indices)),
+        shape=(len(theta), nbins_int),
+    )
+    row_sums = np.array(binning_matrix.sum(axis=1)).flatten()
+    binning_matrix = sparse.diags(1 / row_sums) @ binning_matrix
+
+    # One n(z) gives one tracer pair: get_theo_xi's single (xi+, xi-) entry.
+    (xi_pm,) = get_theo_xi(
+        theta=theta_int, z=z, nz=nz, backend="ccl", cosmo=cosmo
+    ).values()
+    mean_int = np.concatenate(xi_pm)
+    samples_int = np.random.multivariate_normal(mean_int, cov_int, size=n_samples)
+    samples_int_xip, samples_int_xim = (
+        samples_int[:, :nbins_int],
+        samples_int[:, nbins_int:],
+    )
+    samples_rep_xip = (binning_matrix @ samples_int_xip.T).T
+    samples_rep_xim = (binning_matrix @ samples_int_xim.T).T
+
+    eb_samples = np.array(
+        [
+            np.concatenate(
+                get_pure_EB_modes(
+                    theta=theta,
+                    theta_int=theta_int,
+                    xip=samples_rep_xip[i],
+                    xim=samples_rep_xim[i],
+                    xip_int=samples_int_xip[i],
+                    xim_int=samples_int_xim[i],
+                    tmin=left_edges[0],
+                    tmax=right_edges[-1],
+                    parallel=True,
+                )
+            )
+            for i in tqdm.tqdm(range(n_samples), desc="MC samples")
+        ]
+    )
+    return np.cov(eb_samples.T), eb_samples
+
+
 def calculate_cosebis(gg, nmodes=10, scale_cuts=None, cov_path=None):
     """
     Calculate COSEBIs modes from a correlation function for multiple scale cuts.
@@ -279,23 +363,52 @@ def calculate_cosebis(gg, nmodes=10, scale_cuts=None, cov_path=None):
         Each results dictionary contains 'En', 'Bn', 'cov', 'chi2_E', 'chi2_B',
         'pte_B', 'scale_cut', and 'mask' entries.
     """
-    from cosmo_numba.B_modes.cosebis import COSEBIS
+    cov_xipm = np.loadtxt(cov_path) if cov_path is not None else gg.cov
+    return cosebis_scan_from_xi(
+        gg.meanr,
+        gg.xip,
+        gg.xim,
+        cov_xipm,
+        gg.left_edges,
+        gg.right_edges,
+        nmodes=nmodes,
+        scale_cuts=scale_cuts,
+        # A theory covariance has no jackknife realisations to debias.
+        npatch=None if cov_path is not None else gg.npatch1,
+    )
 
-    # Default to full range if no scale cuts provided
-    if scale_cuts is None:
-        scale_cuts = [(gg.left_edges[0], gg.right_edges[-1])]
 
-    # Pre-compute values that don't change across scale cuts
-    nbins = len(gg.meanr)
+def cosebis_scan_from_xi(
+    theta,
+    xip,
+    xim,
+    cov_xipm,
+    left_edges,
+    right_edges,
+    *,
+    nmodes=10,
+    scale_cuts=None,
+    npatch=None,
+):
+    """COSEBIs over a set of scale cuts, from ξ± arrays and their covariance.
 
-    # Load covariance matrix and calculate Hartlap factor once
-    if cov_path is not None:
-        print(f"Loading theoretical covariance from {cov_path}")
-        cov_xipm = np.loadtxt(cov_path)
-        hartlap_factor = 1  # Not defined for analytic covariances
-    else:
-        cov_xipm = gg.cov
-        hartlap_factor = (gg.npatch1 - 2 * nmodes - 2) / (gg.npatch1 - 1)
+    The values-and-covariance seam of :func:`calculate_cosebis`, for callers
+    holding a ξ± data vector rather than a TreeCorr ``GGCorrelation``. The
+    COSEBIs covariance is the ξ± covariance carried through the same linear
+    kernel as the modes, so no estimator re-run is involved. ``npatch`` is the
+    jackknife realisation count behind ``cov_xipm``, which sets the Hartlap
+    debiasing; leave it ``None`` for a theory covariance, which needs none.
+
+    Returns the ``{scale_cut: result}`` mapping :func:`calculate_cosebis` returns.
+    """
+    from cosmo_numba.B_modes.cosebis import COSEBIS
+
+    theta, xip, xim = (np.asarray(a) for a in (theta, xip, xim))
+    cov_xipm = np.asarray(cov_xipm)
+    if scale_cuts is None:
+        scale_cuts = [(left_edges[0], right_edges[-1])]
+    nbins = len(theta)
+    hartlap_factor = 1 if npatch is None else (npatch - 2 * nmodes - 2) / (npatch - 1)
 
     all_results = {}
 
@@ -303,11 +416,12 @@ def calculate_cosebis(gg, nmodes=10, scale_cuts=None, cov_path=None):
     for scale_cut in tqdm.tqdm(scale_cuts, desc="COSEBIs scale cuts"):
         min_theta, max_theta = scale_cut
 
-        # Apply scale cuts using scale_cut_to_bins for consistency
-        start_bin, stop_bin = scale_cut_to_bins(gg, min_theta, max_theta)
+        start_bin, stop_bin = bins_from_edges(
+            left_edges, right_edges, min_theta, max_theta
+        )
         inds = np.arange(start_bin, stop_bin)
 
-        theta_cut, xip_cut, xim_cut = [arr[inds] for arr in [gg.meanr, gg.xip, gg.xim]]
+        theta_cut, xip_cut, xim_cut = [arr[inds] for arr in [theta, xip, xim]]
 
         # Calculate COSEBIs E/B modes using actual theta range (per Axel's recommendation)
         # Use precision=120 (vs default 80) to avoid sympy root convergence failures
@@ -354,11 +468,7 @@ def calculate_cosebis(gg, nmodes=10, scale_cuts=None, cov_path=None):
     return all_results
 
 
-def calculate_eb_statistics(
-    results,
-    cov_path_int=None,
-    n_samples=1000,
-):
+def calculate_eb_statistics(results):
     """
     Calculate E/B mode statistics using 2D PTE analysis for all scale cut combinations.
 
@@ -369,23 +479,17 @@ def calculate_eb_statistics(
     Parameters
     ----------
     results : dict
-        Dictionary containing pure E/B mode results from calculate_pure_eb_correlation
-    cov_path_int : str, optional
-        Path to integration covariance matrix for semi-analytical calculation
-    n_samples : int, optional
-        Number of Monte Carlo samples used for semi-analytical covariance
-    min_bins : int, optional
-        Minimum number of bins required for valid PTE calculation
+        Pure E/B results: the six mode arrays, the ``cov`` block, the reporting
+        ``theta``, and ``n_eff`` — the realisation count behind the covariance
+        (jackknife patches or MC draws), which sets the Hartlap debiasing
 
     Returns
     -------
     dict
         Updated results dictionary with PTE matrices and statistics
     """
-    gg = results["gg"]
-    nbins = gg.nbins
-    npatch = gg.npatch1
-    n_eff = n_samples if cov_path_int is not None else npatch
+    nbins = len(results["theta"])
+    n_eff = results["n_eff"]
 
     # Extract covariance blocks and standard deviations
     cov = results["cov"]
@@ -449,16 +553,16 @@ def calculate_eb_statistics(
     return results
 
 
-def plot_integration_vs_reporting(gg, gg_int, output_path, version):
+def plot_integration_vs_reporting(results, output_path, version):
     """
     Plot integration vs reporting scale comparison.
 
     Parameters
     ----------
-    gg : treecorr.GGCorrelation
-        Reporting scale correlation function
-    gg_int : treecorr.GGCorrelation
-        Integration scale correlation function
+    results : dict
+        Pure E/B results carrying both grids (``theta``/``xip``/``xim`` and the
+        ``theta_int``/``xip_int``/``xim_int`` counterparts), plus the reporting
+        ``var_xip``/``var_xim`` the error bars use
     output_path : str
         Output file path for the plot
     version : str
@@ -468,26 +572,24 @@ def plot_integration_vs_reporting(gg, gg_int, output_path, version):
 
     # Configure plot data for both xi+ and xi- in a consolidated loop
     plot_configs = [
-        ("+", "xip", "varxip", r"$\theta \xi_+(\theta) \times 10^4$"),
-        ("-", "xim", "varxim", r"$\theta \xi_-(\theta) \times 10^4$"),
+        ("+", "xip", r"$\theta \xi_+(\theta) \times 10^4$"),
+        ("-", "xim", r"$\theta \xi_-(\theta) \times 10^4$"),
     ]
 
     data_configs = [
-        (gg_int, "k.", 3, 0.3, "Integration"),
-        (gg, ".", 12, 1, "Reporting"),
+        ("_int", "k.", 3, 0.3, "Integration"),
+        ("", ".", 12, 1, "Reporting"),
     ]
 
-    for ax_idx, (xi_label, xi_attr, var_attr, ylabel) in enumerate(plot_configs):
-        for data, fmt, ms, alpha, label_type in data_configs:
-            xi_val = getattr(data, xi_attr)
-            yerr = (
-                data.meanr * np.sqrt(getattr(data, var_attr)) / 1e-4
-                if hasattr(data, var_attr) and label_type == "Reporting"
-                else None
-            )
+    for ax_idx, (xi_label, xi_attr, ylabel) in enumerate(plot_configs):
+        for suffix, fmt, ms, alpha, label_type in data_configs:
+            theta = results[f"theta{suffix}"]
+            xi_val = results[f"{xi_attr}{suffix}"]
+            var = results.get(f"var_{xi_attr}{suffix}")
+            yerr = theta * np.sqrt(var) / 1e-4 if var is not None else None
             axs[ax_idx].errorbar(
-                data.meanr,
-                data.meanr * xi_val / 1e-4,
+                theta,
+                theta * xi_val / 1e-4,
                 yerr=yerr,
                 fmt=fmt,
                 ms=ms,
@@ -496,8 +598,8 @@ def plot_integration_vs_reporting(gg, gg_int, output_path, version):
                 ls="" if label_type == "Reporting" else None,
                 label=(
                     rf"$\xi_{{{xi_label}}}$, {label_type}: "
-                    rf"${data.min_sep} < \theta < {data.max_sep}$, "
-                    rf"{data.nbins} bins"
+                    rf"${theta[0]:.2g} < \theta < {theta[-1]:.4g}$, "
+                    rf"{len(theta)} bins"
                 ),
             )
         axs[ax_idx].set(
@@ -513,7 +615,7 @@ def plot_integration_vs_reporting(gg, gg_int, output_path, version):
     plt.savefig(output_path, dpi=300, bbox_inches="tight")
 
 
-def _get_pte_from_scale_cut(pte_matrix, gg, scale_cut):
+def _get_pte_from_scale_cut(pte_matrix, edges, scale_cut):
     """
     Extract PTE value from matrix based on scale cut range using conservative logic.
 
@@ -521,8 +623,8 @@ def _get_pte_from_scale_cut(pte_matrix, gg, scale_cut):
     ----------
     pte_matrix : numpy.ndarray
         2D PTE matrix
-    gg : treecorr.GGCorrelation
-        Correlation function object with bin edges
+    edges : tuple of numpy.ndarray
+        ``(left_edges, right_edges)`` of the grid the matrix is indexed on
     scale_cut : tuple
         (min_scale, max_scale) angular range for scale cut
 
@@ -531,7 +633,8 @@ def _get_pte_from_scale_cut(pte_matrix, gg, scale_cut):
     float
         PTE value for the given scale cut, or full-range PTE if scale_cut is None
     """
-    nbins = len(gg.meanr)
+    left_edges, right_edges = edges
+    nbins = len(left_edges)
 
     if scale_cut is None:
         # Return full-range PTE (first row, last column)
@@ -539,8 +642,7 @@ def _get_pte_from_scale_cut(pte_matrix, gg, scale_cut):
 
     min_scale, max_scale = scale_cut
 
-    # Use conservative scale_cut_to_bins helper
-    start_bin, stop_bin = scale_cut_to_bins(gg, min_scale, max_scale)
+    start_bin, stop_bin = bins_from_edges(left_edges, right_edges, min_scale, max_scale)
 
     # Ensure valid range, otherwise fallback to full range
     if stop_bin <= start_bin or start_bin >= nbins or stop_bin <= 0:
@@ -571,19 +673,20 @@ def plot_pure_eb_correlations(
     fiducial_xim_scale_cut : tuple, optional
         (min_scale, max_scale) for xi- fiducial analysis, shown as gray regions
     """
-    gg = results["gg"]
-    nbins = gg.nbins
+    theta = results["theta"]
+    edges = (results["left_edges"], results["right_edges"])
+    nbins = len(theta)
     cov = results["cov"]
 
     # Calculate combined PTE using off-diagonal covariance blocks
     # Get scale cuts for both xi+ and xi-
     if fiducial_xip_scale_cut is not None:
-        xip_start_bin, xip_stop_bin = scale_cut_to_bins(gg, *fiducial_xip_scale_cut)
+        xip_start_bin, xip_stop_bin = bins_from_edges(*edges, *fiducial_xip_scale_cut)
     else:
         xip_start_bin, xip_stop_bin = 0, nbins
 
     if fiducial_xim_scale_cut is not None:
-        xim_start_bin, xim_stop_bin = scale_cut_to_bins(gg, *fiducial_xim_scale_cut)
+        xim_start_bin, xim_stop_bin = bins_from_edges(*edges, *fiducial_xim_scale_cut)
     else:
         xim_start_bin, xim_stop_bin = 0, nbins
 
@@ -618,7 +721,7 @@ def plot_pure_eb_correlations(
     if "eb_samples" in results:  # Semi-analytical case
         n_eff = results["eb_samples"].shape[0]
     else:  # Jackknife case
-        n_eff = gg.npatch1
+        n_eff = results["n_eff"]
 
     hartlap_factor = (n_eff - nbins_eff - 2) / (n_eff - 1)
     chi2_combined = hartlap_factor * (
@@ -628,10 +731,10 @@ def plot_pure_eb_correlations(
 
     # Extract PTE values for fiducial scale cuts (or full range)
     xip_B_pte = _get_pte_from_scale_cut(
-        results["pte_matrices"]["xip_B"], gg, fiducial_xip_scale_cut
+        results["pte_matrices"]["xip_B"], edges, fiducial_xip_scale_cut
     )
     xim_B_pte = _get_pte_from_scale_cut(
-        results["pte_matrices"]["xim_B"], gg, fiducial_xim_scale_cut
+        results["pte_matrices"]["xim_B"], edges, fiducial_xim_scale_cut
     )
 
     fig, axs = plt.subplots(1, 2, figsize=(14, 6), sharex=True, sharey=True)
@@ -643,14 +746,14 @@ def plot_pure_eb_correlations(
         (
             "xip",
             "+",
-            "varxip",
+            "var_xip",
             r"$\xi_{+}=\xi_{+}^{E}+\xi_{+}^{B}+\xi_{+}^{\mathrm{amb}}$",
             xip_B_pte,
         ),
         (
             "xim",
             "-",
-            "varxim",
+            "var_xim",
             r"$\xi_{-}=\xi_{-}^{E}-\xi_{-}^{B}+\xi_{-}^{\mathrm{amb}}$",
             xim_B_pte,
         ),
@@ -660,11 +763,11 @@ def plot_pure_eb_correlations(
         plot_configs
     ):
         # Plot main correlation function
-        xi_val = getattr(gg, xi_type)
+        xi_val = results[xi_type]
         axs[ax_idx].errorbar(
-            gg.meanr,
-            gg.meanr * xi_val / scale_factor,
-            yerr=gg.meanr * np.sqrt(getattr(gg, var_attr)) / scale_factor,
+            theta,
+            theta * xi_val / scale_factor,
+            yerr=theta * np.sqrt(results[var_attr]) / scale_factor,
             fmt="k.",
             capsize=3,
             label=main_label,
@@ -690,9 +793,9 @@ def plot_pure_eb_correlations(
 
         for key, color, alpha, label in plot_data:
             axs[ax_idx].errorbar(
-                gg.meanr,
-                gg.meanr * results[key] / scale_factor,
-                yerr=gg.meanr * results[f"std_{key}"] / scale_factor,
+                theta,
+                theta * results[key] / scale_factor,
+                yerr=theta * results[f"std_{key}"] / scale_factor,
                 color=color,
                 ls="",
                 marker=".",
@@ -722,12 +825,12 @@ def plot_pure_eb_correlations(
             xlim = original_xlims[ax_idx]
 
             # Use conservative scale_cut_to_bins helper for consistency
-            start_bin, stop_bin = scale_cut_to_bins(gg, min_scale, max_scale)
+            start_bin, stop_bin = bins_from_edges(*edges, min_scale, max_scale)
 
             # Show excluded regions based on bin edges used in PTE calculation
             # Lower exclusion: bins 0 to start_bin-1 are excluded
             if start_bin > 0:
-                lower_exclusion_edge = gg.right_edges[start_bin - 1]
+                lower_exclusion_edge = edges[1][start_bin - 1]
                 axs[ax_idx].axvspan(
                     xlim[0],
                     lower_exclusion_edge,
@@ -737,8 +840,8 @@ def plot_pure_eb_correlations(
                 )
 
             # Upper exclusion: bins stop_bin to end are excluded
-            if stop_bin < len(gg.left_edges):
-                upper_exclusion_edge = gg.left_edges[stop_bin]
+            if stop_bin < len(edges[0]):
+                upper_exclusion_edge = edges[0][stop_bin]
                 axs[ax_idx].axvspan(
                     upper_exclusion_edge,
                     xlim[1],
@@ -760,7 +863,7 @@ def plot_pure_eb_correlations(
 
 
 def plot_cosebis_scale_cut_heatmap(
-    cosebis_results, gg, version, output_path, fiducial_scale_cut=None
+    cosebis_results, edges, version, output_path, fiducial_scale_cut=None
 ):
     """
     Create 2D heatmaps showing how COSEBIs statistics vary across different scale cuts.
@@ -769,8 +872,8 @@ def plot_cosebis_scale_cut_heatmap(
     ----------
     cosebis_results : dict
         Dictionary with scale cut tuples as keys, containing 'chi2_E' and 'pte_B' values
-    gg : treecorr.GGCorrelation
-        Correlation function object for bin edges
+    edges : tuple of numpy.ndarray
+        ``(left_edges, right_edges)`` of the grid the scale cuts index
     version : str
         Version string for main title
     output_path : str
@@ -778,7 +881,8 @@ def plot_cosebis_scale_cut_heatmap(
     fiducial_scale_cut : tuple, optional
         (min_scale, max_scale) for cross-hatching
     """
-    nbins = gg.nbins
+    left_edges, right_edges = edges
+    nbins = len(left_edges)
 
     # Initialize matrices
     snrs = [np.sqrt(result["chi2_E"]) for result in cosebis_results.values()]
@@ -794,9 +898,9 @@ def plot_cosebis_scale_cut_heatmap(
         pte_matrix[row, column] = pte
 
     # Fill matrices from COSEBIs results
-    for i in range(len(gg.left_edges)):
-        for j in range(i, len(gg.right_edges)):
-            scale_cut = (gg.left_edges[i], gg.right_edges[j])
+    for i in range(len(left_edges)):
+        for j in range(i, len(right_edges)):
+            scale_cut = (left_edges[i], right_edges[j])
             result = cosebis_results.get(scale_cut)
 
             if result is not None:
@@ -857,7 +961,9 @@ def plot_cosebis_scale_cut_heatmap(
     # Add fiducial scale cut cross-hatching if provided
     if fiducial_scale_cut is not None:
         min_scale, max_scale = fiducial_scale_cut
-        start_bin, stop_bin = scale_cut_to_bins(gg, min_scale, max_scale)
+        start_bin, stop_bin = bins_from_edges(
+            left_edges, right_edges, min_scale, max_scale
+        )
 
         if stop_bin > start_bin and start_bin < nbins and stop_bin > 0:
             # Add cross-hatching at the fiducial scale cut matrix element
@@ -884,10 +990,10 @@ def plot_cosebis_scale_cut_heatmap(
     # Set angular scale ticks
     tick_indices = np.arange(0, nbins)
     x_tick_labels = [
-        f"{gg.left_edges[i]:.1f}" for i in tick_indices if i < len(gg.left_edges)
+        f"{left_edges[i]:.1f}" for i in tick_indices if i < len(left_edges)
     ]
     y_tick_labels = [
-        f"{gg.right_edges[i]:.1f}" for i in tick_indices if i < len(gg.right_edges)
+        f"{right_edges[i]:.1f}" for i in tick_indices if i < len(right_edges)
     ]
     x_tick_positions = tick_indices + 0.5
     y_tick_positions = tick_indices + 0.5
@@ -939,8 +1045,9 @@ def plot_pte_2d_heatmaps(
     fiducial_xim_scale_cut : tuple, optional
         (min_scale, max_scale) for xi- fiducial analysis, shown as cross-hatched
     """
-    gg = results["gg"]
-    nbins = gg.nbins
+    theta = results["theta"]
+    edges = (results["left_edges"], results["right_edges"])
+    nbins = len(theta)
     pte_xip_B = results["pte_matrices"]["xip_B"]
     pte_xim_B = results["pte_matrices"]["xim_B"]
 
@@ -1002,7 +1109,7 @@ def plot_pte_2d_heatmaps(
     for ax_idx, fiducial_scale_cut in enumerate(fiducial_scale_cuts):
         if fiducial_scale_cut is not None:
             min_scale, max_scale = fiducial_scale_cut
-            start_bin, stop_bin = scale_cut_to_bins(gg, min_scale, max_scale)
+            start_bin, stop_bin = bins_from_edges(*edges, min_scale, max_scale)
 
             if stop_bin > start_bin and start_bin < nbins and stop_bin > 0:
                 rect_x = start_bin
@@ -1026,12 +1133,8 @@ def plot_pte_2d_heatmaps(
 
     # Set angular scale ticks
     tick_indices = np.arange(0, nbins)
-    x_tick_labels = [
-        f"{gg.left_edges[i]:.1f}" for i in tick_indices if i < len(gg.left_edges)
-    ]
-    y_tick_labels = [
-        f"{gg.right_edges[i]:.1f}" for i in tick_indices if i < len(gg.right_edges)
-    ]
+    x_tick_labels = [f"{edges[0][i]:.1f}" for i in tick_indices if i < len(edges[0])]
+    y_tick_labels = [f"{edges[1][i]:.1f}" for i in tick_indices if i < len(edges[1])]
     x_tick_positions = tick_indices + 0.5
     y_tick_positions = tick_indices + 0.5
 
@@ -1171,10 +1274,8 @@ def save_pure_eb_results(results, output_path):
     output_path : str
         Output .npz file path
     """
-    gg = results["gg"]
-
     # Data vectors and covariance
-    save_dict = {"theta": gg.meanr, "cov": results["cov"]}
+    save_dict = {"theta": results["theta"], "cov": results["cov"]}
     for key in _EB_KEYS:
         save_dict[key] = results[key]
 
@@ -1183,7 +1284,7 @@ def save_pure_eb_results(results, output_path):
         save_dict[f"pte_matrices_{key}"] = matrix
 
     # Metadata
-    save_dict["npatch"] = np.array(gg.npatch1)
+    save_dict["n_eff"] = np.array(results["n_eff"])
     if "eb_samples" in results:
         save_dict["var_method"] = np.array("semi-analytic")
         save_dict["n_samples"] = np.array(results["eb_samples"].shape[0])
diff --git a/src/sp_validation/cosmo_val/core.py b/src/sp_validation/cosmo_val/core.py
index 48251bbb..1d9fa196 100644
--- a/src/sp_validation/cosmo_val/core.py
+++ b/src/sp_validation/cosmo_val/core.py
@@ -563,18 +563,18 @@ def summarize_bmodes(self, fiducial_scale_cut=(12, 83), versions=None):
             # Pure E/B PTEs from stored results
             if ver in self._pure_eb_results:
                 res = self._pure_eb_results[ver]
-                gg = res["gg"]
+                edges = (res["left_edges"], res["right_edges"])
                 try:
                     for stat in ("xip_B", "xim_B", "combined"):
                         row[stat] = _get_pte_from_scale_cut(
-                            res["pte_matrices"][stat], gg, fiducial_scale_cut
+                            res["pte_matrices"][stat], edges, fiducial_scale_cut
                         )
                 except (KeyError, RuntimeError):
                     pass
                 cov_methods.add(
                     "semi-analytic"
                     if "eb_samples" in res
-                    else f"jackknife ({gg.npatch1} patches)"
+                    else f"jackknife ({res['n_eff']} patches)"
                 )
 
             # COSEBIs PTE from stored results
diff --git a/src/sp_validation/cosmo_val/cosebis.py b/src/sp_validation/cosmo_val/cosebis.py
index aa4657ba..312aadcd 100644
--- a/src/sp_validation/cosmo_val/cosebis.py
+++ b/src/sp_validation/cosmo_val/cosebis.py
@@ -271,7 +271,7 @@ def plot_cosebis(
 
             plot_cosebis_scale_cut_heatmap(
                 results,
-                gg_temp,
+                (gg_temp.left_edges, gg_temp.right_edges),
                 version,
                 out_stub + "_scalecut_ptes.png",
                 fiducial_scale_cut=fiducial_scale_cut,
diff --git a/src/sp_validation/cosmo_val/pure_eb.py b/src/sp_validation/cosmo_val/pure_eb.py
index 71f5d76a..90157456 100644
--- a/src/sp_validation/cosmo_val/pure_eb.py
+++ b/src/sp_validation/cosmo_val/pure_eb.py
@@ -88,9 +88,13 @@ def calculate_pure_eb(
             - "xip_amb": Ambiguity mode for xi+.
             - "xim_amb": Ambiguity mode for xi-.
             - "cov": Covariance matrix for the pure E/B modes.
-            - "gg": The two-point correlation function object for the reporting binning.
-            - "gg_int": The two-point correlation function object for the
-              integration binning.
+            - "theta", "left_edges", "right_edges": Reporting-grid bin centres
+              and edges.
+            - "xip", "xim", "var_xip", "var_xim": Reporting-grid xi+/xi- and
+              their variances.
+            - "theta_int", "xip_int", "xim_int": Integration-grid xi+/xi-.
+            - "n_eff": Realisation count behind "cov" (jackknife patches or
+              MC draws), which sets the Hartlap debiasing.
             - "eb_samples": (only when using semi-analytical covariance) Semi-analytic
               EB samples used for covariance calculation. Shape: (n_samples, 6*nbins)
 
@@ -269,19 +273,13 @@ def plot_pure_eb(
             )
 
             # Calculate E/B statistics for all bin combinations
-            version_results = calculate_eb_statistics(
-                version_results,
-                cov_path_int=cov_path_int,
-                n_samples=n_samples,
-                **kwargs,
-            )
-
-            # Generate all plots using specialized plotting functions
-            gg, gg_int = version_results["gg"], version_results["gg_int"]
+            version_results = calculate_eb_statistics(version_results, **kwargs)
 
             # Integration vs Reporting comparison plot
             plot_integration_vs_reporting(
-                gg, gg_int, out_stub + "_integration_vs_reporting.png", version
+                version_results,
+                out_stub + "_integration_vs_reporting.png",
+                version,
             )
 
             # E/B/Ambiguous correlation functions plot
diff --git a/src/sp_validation/tests/test_b_modes.py b/src/sp_validation/tests/test_b_modes.py
index e8e05c39..67c64376 100644
--- a/src/sp_validation/tests/test_b_modes.py
+++ b/src/sp_validation/tests/test_b_modes.py
@@ -50,6 +50,7 @@ def _eb_inputs():
 
     nbins=4, npatch=50 so the Hartlap factor (n_eff - nbins_eff - 2)/(n_eff-1)
     is well-defined and strictly positive for every scale-cut combination.
+    n_eff is the jackknife patch count, as it is for a jackknife covariance.
     The covariance is built SPD via A @ A.T + I; the B-mode vectors are O(1)
     so the chi-squared (and hence PTE) lands in a meaningful range rather than
     being saturated at 1.0.
@@ -60,8 +61,13 @@ def _eb_inputs():
     cov = A @ A.T + np.eye(6 * nbins)
     xip_B = rng.standard_normal(nbins)
     xim_B = rng.standard_normal(nbins)
-    gg = types.SimpleNamespace(nbins=nbins, npatch1=npatch)
-    return {"gg": gg, "cov": cov, "xip_B": xip_B, "xim_B": xim_B}, nbins
+    return {
+        "theta": np.geomspace(1.0, 100.0, nbins),
+        "n_eff": npatch,
+        "cov": cov,
+        "xip_B": xip_B,
+        "xim_B": xim_B,
+    }, nbins
 
 
 # ---------------------------------------------------------------------------
@@ -229,8 +235,8 @@ def test_calculate_eb_statistics_pte_matrices():
     """Pin representative PTE-matrix entries from the full 2D E/B analysis.
 
     Inputs are fixed (seed 12345, nbins=4, npatch=50, SPD cov = A@A.T + I,
-    O(1) B-mode vectors). With cov_path_int=None the Hartlap correction uses
-    n_eff = npatch = 50. For each of xip_B, xim_B and combined we pin the
+    O(1) B-mode vectors). The Hartlap correction uses n_eff = 50, the patch
+    count behind a jackknife covariance. For each of xip_B, xim_B and combined we pin the
     full-range entry [0, nbins-1] (start=0, stop=nbins) and an interior entry
     [0, 2] (start=0, stop=3). These chi2->sf PTE values are deterministic
     functions of the seeded input.
@@ -240,7 +246,7 @@ def test_calculate_eb_statistics_pte_matrices():
     arithmetic shifts them past tolerance.
     """
     results, nbins = _eb_inputs()
-    out = b_modes.calculate_eb_statistics(results, cov_path_int=None)
+    out = b_modes.calculate_eb_statistics(results)
     pm = out["pte_matrices"]
 
     # Full-range entries [0, nbins-1].
@@ -271,13 +277,13 @@ def test_calculate_eb_statistics_has_teeth():
     combined 0.99999 -> 0.0031.
     """
     results, nbins = _eb_inputs()
-    out = b_modes.calculate_eb_statistics(results, cov_path_int=None)
+    out = b_modes.calculate_eb_statistics(results)
     pm = out["pte_matrices"]
 
     loud, _ = _eb_inputs()
     loud["xip_B"] = loud["xip_B"] * 10.0
     loud["xim_B"] = loud["xim_B"] * 10.0
-    out_loud = b_modes.calculate_eb_statistics(loud, cov_path_int=None)
+    out_loud = b_modes.calculate_eb_statistics(loud)
     pm_loud = out_loud["pte_matrices"]
 
     for key in ("xip_B", "xim_B", "combined"):
@@ -285,3 +291,193 @@ def test_calculate_eb_statistics_has_teeth():
         loud_pte = pm_loud[key][0, nbins - 1]
         assert loud_pte < quiet_pte
         assert loud_pte < 0.05  # louder B-modes are clearly rejected
+
+
+# ---------------------------------------------------------------------------
+# 6. Grid edges and the COSEBIs covariance seam
+# ---------------------------------------------------------------------------
+
+
+def test_log_bin_edges_matches_the_grid_stub():
+    """Edges reconstructed from a binning are the ones TreeCorr would report.
+
+    A part stores bin centres only, so a consumer rebuilds the edges from the
+    binning it was measured on; the two must agree bin for bin.
+    """
+    left, right = b_modes.log_bin_edges(1.0, 100.0, _NBINS_GRID)
+    gg = _grid_gg()
+    npt.assert_allclose(left, gg.left_edges)
+    npt.assert_allclose(right, gg.right_edges)
+    # ...and they index scale cuts identically.
+    assert b_modes.bins_from_edges(left, right, 2.0, 50.0) == (2, 8)
+
+
+def test_cosebis_scan_propagates_the_supplied_covariance(monkeypatch):
+    """The COSEBIs covariance is the ξ± covariance through the same kernel.
+
+    The kernel is stubbed, so what is pinned is the seam: which ξ± covariance
+    sub-block is handed to the transform (the scale cut's, in [ξ+; ξ−] order)
+    and that Hartlap uses the supplied npatch.
+    """
+    nbins, nmodes = _NBINS_GRID, 3
+    theta = np.geomspace(1.2, 90.0, nbins)
+    cov_xipm = np.diag(np.arange(1.0, 2 * nbins + 1))
+    seen = {}
+
+    class _StubCOSEBIS:
+        def __init__(self, **kwargs):
+            seen["init"] = kwargs
+
+        def cosebis_from_xipm(self, theta_cut, xip_cut, xim_cut, parallel=True):
+            seen["n_theta"] = len(theta_cut)
+            return np.ones(nmodes), np.full(nmodes, 2.0)
+
+        def cosebis_covariance_from_xipm_covariance(self, theta_cut, cov_cut):
+            seen["cov_cut"] = cov_cut
+            return np.eye(2 * nmodes)
+
+    module = types.ModuleType("cosmo_numba.B_modes.cosebis")
+    module.COSEBIS = _StubCOSEBIS
+    monkeypatch.setitem(
+        __import__("sys").modules, "cosmo_numba.B_modes.cosebis", module
+    )
+
+    left, right = b_modes.log_bin_edges(1.0, 100.0, nbins)
+    results = b_modes.cosebis_scan_from_xi(
+        theta,
+        np.arange(nbins) * 1e-5,
+        np.arange(nbins) * 2e-5,
+        cov_xipm,
+        left,
+        right,
+        nmodes=nmodes,
+        scale_cuts=[(2.0, 50.0)],
+        npatch=100,
+    )
+
+    (result,) = results.values()
+    # The cut is bins 2..8, so the covariance sub-block is those rows/cols in
+    # both the ξ+ and the ξ− half.
+    inds = np.concatenate([np.arange(2, 8), np.arange(2, 8) + nbins])
+    npt.assert_array_equal(seen["cov_cut"], cov_xipm[np.ix_(inds, inds)])
+    assert seen["n_theta"] == 6
+    npt.assert_allclose(result["hartlap_factor"], (100 - 2 * nmodes - 2) / 99)
+    # χ² carries the Hartlap factor: modes are 1, cov is the identity.
+    npt.assert_allclose(result["chi2_E"], nmodes * result["hartlap_factor"])
+
+
+def test_cosebis_scan_theory_covariance_skips_hartlap(monkeypatch):
+    """A theory covariance has no realisations to debias, so Hartlap is 1."""
+    nbins, nmodes = _NBINS_GRID, 2
+
+    class _StubCOSEBIS:
+        def __init__(self, **kwargs):
+            pass
+
+        def cosebis_from_xipm(self, theta_cut, xip_cut, xim_cut, parallel=True):
+            return np.ones(nmodes), np.ones(nmodes)
+
+        def cosebis_covariance_from_xipm_covariance(self, theta_cut, cov_cut):
+            return np.eye(2 * nmodes)
+
+    module = types.ModuleType("cosmo_numba.B_modes.cosebis")
+    module.COSEBIS = _StubCOSEBIS
+    monkeypatch.setitem(
+        __import__("sys").modules, "cosmo_numba.B_modes.cosebis", module
+    )
+
+    left, right = b_modes.log_bin_edges(1.0, 100.0, nbins)
+    (result,) = b_modes.cosebis_scan_from_xi(
+        np.geomspace(1.2, 90.0, nbins),
+        np.zeros(nbins),
+        np.zeros(nbins),
+        np.eye(2 * nbins),
+        left,
+        right,
+        nmodes=nmodes,
+        npatch=None,
+    ).values()
+    assert result["hartlap_factor"] == 1
+
+
+def test_pure_eb_npz_carries_what_the_summary_reads(tmp_path):
+    """The .npz keys cv_summarize_bmodes reads are the ones the writer emits.
+
+    The two live in different rules, so the contract between them — the PTE
+    matrices under ``pte_matrices_{stat}`` and the realisation count under
+    ``n_eff`` — is pinned here rather than discovered on a cluster run.
+    """
+    results, nbins = _eb_inputs()
+    results.update(
+        {key: np.zeros(nbins) for key in b_modes._EB_KEYS if key not in results}
+    )
+    results = b_modes.calculate_eb_statistics(results)
+
+    out = tmp_path / "pure_eb_data.npz"
+    b_modes.save_pure_eb_results(results, str(out))
+    saved = np.load(out)
+
+    for stat in ("xip_B", "xim_B", "combined"):
+        assert f"pte_matrices_{stat}" in saved
+        assert saved[f"pte_matrices_{stat}"].shape == (nbins, nbins)
+    assert saved["n_eff"] == results["n_eff"]
+    npt.assert_allclose(saved["theta"], results["theta"])
+    for key in b_modes._EB_KEYS:
+        assert key in saved
+
+    # The summary reads the fiducial cut out of those matrices through the same
+    # helper the plots use, so a valid cut must resolve to a finite PTE.
+    edges = b_modes.log_bin_edges(1.0, 100.0, nbins)
+    pte = b_modes._get_pte_from_scale_cut(
+        saved["pte_matrices_xip_B"], edges, (1.0, 100.0)
+    )
+    assert np.isfinite(pte)
+
+
+def test_pure_eb_covariance_mc_draws_around_the_theory_mean(monkeypatch):
+    """The MC draws centre on cs_util's theory ξ±, binned to the reporting grid.
+
+    ``get_theo_xi`` returns ``{pair: (ξ+, ξ−)}``; one n(z) is one pair. With a
+    zero covariance every draw is the theory mean, so a stub kernel that echoes
+    its reporting-grid ξ± pins the unpack, the [ξ+; ξ−] order and the binning.
+    """
+    nrep = 4
+    left, right = b_modes.log_bin_edges(2.0, 50.0, nrep)
+    theta = np.sqrt(left * right)
+    theta_int = np.geomspace(1.0, 100.0, 40)
+    xip_th, xim_th = theta_int.copy(), 2 * theta_int
+
+    monkeypatch.setattr(
+        b_modes, "get_theo_xi", lambda **kw: {"W0xW0": (xip_th, xim_th)}
+    )
+
+    def _echo(theta, theta_int, xip, xim, xip_int, xim_int, tmin, tmax, parallel):
+        zeros = np.zeros_like(xip)
+        return xip, xim, zeros, zeros, zeros, zeros
+
+    module = types.ModuleType("cosmo_numba.B_modes.schneider2022")
+    module.get_pure_EB_modes = _echo
+    monkeypatch.setitem(
+        __import__("sys").modules, "cosmo_numba.B_modes.schneider2022", module
+    )
+
+    cov, eb_samples = b_modes.pure_eb_covariance_mc(
+        theta=theta,
+        left_edges=left,
+        right_edges=right,
+        theta_int=theta_int,
+        cov_int=np.zeros((2 * len(theta_int), 2 * len(theta_int))),
+        z=np.linspace(0.01, 2.0, 50),
+        nz=np.ones(50),
+        cosmo=None,
+        n_samples=3,
+    )
+
+    inside = [(theta_int >= lo) & (theta_int < hi) for lo, hi in zip(left, right)]
+    expected_xip = np.array([xip_th[m].mean() for m in inside])
+    expected_xim = np.array([xim_th[m].mean() for m in inside])
+    assert eb_samples.shape == (3, 6 * nrep)
+    for draw in eb_samples:
+        npt.assert_allclose(draw[:nrep], expected_xip)
+        npt.assert_allclose(draw[nrep : 2 * nrep], expected_xim)
+    npt.assert_allclose(cov, 0.0, atol=1e-20)
diff --git a/src/sp_validation/tests/test_cosmo_val.py b/src/sp_validation/tests/test_cosmo_val.py
index f50992d4..f25ed522 100644
--- a/src/sp_validation/tests/test_cosmo_val.py
+++ b/src/sp_validation/tests/test_cosmo_val.py
@@ -697,4 +697,4 @@ def test_calculate_pure_eb_runs_on_synthetic_catalog(self, tmp_path):
         # shape is pinned, not the values.
         cov = np.asarray(results["cov"])
         assert cov.shape == (6 * nbins, 6 * nbins)
-        assert results["gg"].npatch1 == npatch
+        assert results["n_eff"] == npatch