diff --git a/README.md b/README.md
index e14a2071..36d5589e 100644
--- a/README.md
+++ b/README.md
@@ -87,6 +87,25 @@ docker run --rm -it ghcr.io/cosmostat/sp_validation:develop python -c "import sp
 
 We do not currently build images for Apple Silicon/arm64; however the amd64 images should work on these systems, albeit with reduced performance.
 
+## Local Installation
+
+Requires Python ≥ 3.12. With [uv](https://docs.astral.sh/uv/):
+
+```bash
+uv venv --python 3.12
+uv pip install -e '.[test]'
+```
+
+To also install the data-vector blinding stack (Smokescreen + firecrown, PRD
+[#241](https://github.com/CosmoStat/sp_validation/issues/241)), pass the
+dependency-override file — firecrown is not pip-resolvable without it (see
+`uv-overrides.txt` for why):
+
+```bash
+uv pip install --overrides uv-overrides.txt -e '.[test,blinding]'
+python scripts/patch_firecrown.py  # make pip-installed firecrown importable without NumCosmo
+```
+
 
 
 ## Flow chart
diff --git a/papers/bmodes/rules/claims.smk b/papers/bmodes/rules/claims.smk
index 43962627..388de007 100644
--- a/papers/bmodes/rules/claims.smk
+++ b/papers/bmodes/rules/claims.smk
@@ -86,10 +86,10 @@ def _xi_reporting_path(version):
 
 
 def _xi_integration_path(version):
-    """Path to fine-binned 2PCF integration file."""
+    """Path to fine-binned 2PCF integration file. Unpatched: values only, no covariance."""
     return (
         f"{COSMO_VAL_OUTPUT}/{version}_xi_minsep={FIDUCIAL['min_sep_int']}"
-        f"_maxsep={FIDUCIAL['max_sep_int']}_nbins={FIDUCIAL['nbins_int']}_npatch={FIDUCIAL['npatch']}.txt"
+        f"_maxsep={FIDUCIAL['max_sep_int']}_nbins={FIDUCIAL['nbins_int']}_npatch=1.txt"
     )
 
 
@@ -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/run_xi_sweep.py b/papers/bmodes/scripts/run_xi_sweep.py
index b02283a5..43233474 100644
--- a/papers/bmodes/scripts/run_xi_sweep.py
+++ b/papers/bmodes/scripts/run_xi_sweep.py
@@ -2,8 +2,8 @@
 
 Loops the [non-fiducial version list](sweep_versions.nonfiducial_versions) and
 runs the same ``run_2pcf.run_2pcf`` compute the fiducial two_point recipes call,
-once per version, writing every version's ξ± text dump (+ ξ+/ξ- FITS) into one
-lc ``{output}`` dir under run_2pcf's native, already-canonical name
+once per version, writing every version's ξ± text dump into one lc ``{output}``
+dir under run_2pcf's native, already-canonical name
 ``{ver}_xi_minsep={min}_maxsep={max}_nbins={nbins}_npatch={npatch}.txt`` — the
 exact pattern ``cosebis_version_comparison._xi_integration`` reconstructs.
 
@@ -76,7 +76,7 @@ def _from_cli(argv=None):
                 ver=ver,
                 cat_config=a.cat_config,
                 output_dir=a.out,
-                save_fits=True,
+                grid=grid,
                 **GRIDS[grid],
             )
 
diff --git a/papers/cosmo_val/config/config.yaml b/papers/cosmo_val/config/config.yaml
index 98409371..cbbe33cf 100644
--- a/papers/cosmo_val/config/config.yaml
+++ b/papers/cosmo_val/config/config.yaml
@@ -58,11 +58,11 @@ cosmo_val:
         kmax: 20
         kmax_extrapolate: 500
 
-  # Pure E/B-mode decomposition (config space)
-  pure_eb:
-    min_sep_int: 0.08
-    max_sep_int: 300
-    nbins_int: 1000
+  # The fine ξ± grid the B-mode integrals run over.
+  integration:
+    min_sep: 0.08
+    max_sep: 300
+    nbins: 1000
 
   # COSEBIs decomposition (config space, fine integration binning)
   cosebis:
diff --git a/pyproject.toml b/pyproject.toml
index c2a9aa2c..d278d3f7 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -79,7 +79,8 @@ dependencies = [
     "pymaster",
     "regions",
     "reproject",
-    "sacc>=0.12",
+    # Floor matches the resolved lock; needs sacc's 2.x rewrite (BlockDiagonalCovariance).
+    "sacc>=2.4,<3",
     # scipy 1.18 ported FITPACK from Fortran to C, changing the return shape of
     # RectBivariateSpline(scalar, scalar, grid=False) from 0-d `array(x)` to
     # shape-(1,) `array([x])`. camb's BBN Y_He predictor (bbn.py) wraps the
@@ -152,8 +153,8 @@ glass = [
 # Snakemake workflow and cross-validation runners are available.
 workflow = [
     "snakemake",
-    # run_2pcf_highres.py drives the MPI convergence run; the container ships
-    # OpenMPI (/opt/ompi) so mpi4py builds against it.
+    # Optional MPI runners (the container ships OpenMPI at /opt/ompi, so mpi4py
+    # builds against it).
     "mpi4py",
     # NOTE: workflow/scripts/cv_*.py also import `cv_runner`, which is not
     # published or resolvable (no public repo found) — left undeclared pending
diff --git a/src/sp_validation/b_modes.py b/src/sp_validation/b_modes.py
index 13f4c2e9..c4b9a2b7 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,112 @@ 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
+
+    mean_int = np.concatenate(
+        get_theo_xi(theta=theta_int, z=z, nz=nz, backend="ccl", cosmo=cosmo)
+    )
+    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 +361,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 +414,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 +466,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 +477,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 +551,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 +570,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 +596,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 +613,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 +621,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 +631,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 +640,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 +671,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 +719,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 +729,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 +744,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 +761,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 +791,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 +823,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 +838,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 +861,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 +870,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 +879,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 +896,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 +959,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 +988,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 +1043,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 +1107,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 +1131,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 +1272,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 +1282,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..bff3581a 100644
--- a/src/sp_validation/cosmo_val/core.py
+++ b/src/sp_validation/cosmo_val/core.py
@@ -15,6 +15,7 @@
     find_conservative_scale_cut_key,
 )
 from ..statistics import chi2_and_pte
+from ..version import __version__
 from .catalog_characterization import CatalogCharacterizationMixin
 from .cosebis import CosebisMixin
 from .pseudo_cl import PseudoClMixin
@@ -22,8 +23,42 @@
 from .pure_eb import PureEBMixin
 from .real_space import RealSpaceMixin
 
-
 # %%
+BMODE_COLUMNS = {
+    "xip_B": r"xi+B",
+    "xim_B": r"xi-B",
+    "combined": "Combined",
+    "COSEBIS": "COSEBIS",
+    "C_l_BB": "C_l^BB",
+}
+
+
+def print_bmode_summary(summary, fiducial_scale_cut, cov_methods=()):
+    """Print the B-mode PTE table for ``{version: {statistic: pte}}``.
+
+    Statistics absent from a row print as ``--``.
+    """
+    sc_label = f"[{fiducial_scale_cut[0]}-{fiducial_scale_cut[1]} arcmin]"
+    sep = "\u2500" * 70
+    header = f"{'Version':<28s}" + "".join(
+        f"{label:>10s}" for label in BMODE_COLUMNS.values()
+    )
+
+    print(f"\nB-mode summary {sc_label}")
+    print(sep)
+    print(header)
+    print(sep)
+    for ver, row in summary.items():
+        cells = "".join(
+            f"{row[s]:>10.4f}" if s in row else f"{'--':>10s}" for s in BMODE_COLUMNS
+        )
+        print(f"{ver:<28s}{cells}")
+    print(sep)
+    if cov_methods:
+        print(f"Covariance: {', '.join(sorted(cov_methods))}")
+    print()
+
+
 class CosmologyValidation(
     CosebisMixin,
     PureEBMixin,
@@ -380,6 +415,21 @@ def _output_path(self, *parts):
         """
         return os.path.abspath(os.path.join(self.cc["paths"]["output"], *parts))
 
+    def sacc_nz(self, version):
+        """Single-bin ``nz`` mapping ``{0: (z, nz)}`` for the SACC writers.
+
+        The round is single-bin, so the whole survey n(z) is bin 0.
+        """
+        return {0: tuple(self.get_redshift(version))}
+
+    def sacc_metadata(self, version):
+        """Provenance metadata stored on every SACC part for ``version``."""
+        return {
+            "catalogue_version": version,
+            "sp_validation_version": __version__,
+            "npatch": self.npatch,
+        }
+
     def get_redshift(self, version):
         """Load redshift distribution for a catalog version.
 
@@ -563,18 +613,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
@@ -602,37 +652,5 @@ def summarize_bmodes(self, fiducial_scale_cut=(12, 83), versions=None):
 
             summary[ver] = row
 
-        # Print summary table
-        col_labels = {
-            "xip_B": r"xi+B",
-            "xim_B": r"xi-B",
-            "combined": "Combined",
-            "COSEBIS": "COSEBIS",
-            "C_l_BB": "C_l^BB",
-        }
-        stats_order = list(col_labels)
-
-        sc_label = f"[{fiducial_scale_cut[0]}-{fiducial_scale_cut[1]} arcmin]"
-        sep = "\u2500" * 70
-        header = f"{'Version':<28s}" + "".join(
-            f"{label:>10s}" for label in col_labels.values()
-        )
-
-        print(f"\nB-mode summary {sc_label}")
-        print(sep)
-        print(header)
-        print(sep)
-
-        for ver in versions:
-            row = summary[ver]
-            cells = "".join(
-                f"{row[s]:>10.4f}" if s in row else f"{'--':>10s}" for s in stats_order
-            )
-            print(f"{ver:<28s}{cells}")
-
-        print(sep)
-        if cov_methods:
-            print(f"Covariance: {', '.join(sorted(cov_methods))}")
-        print()
-
+        print_bmode_summary(summary, fiducial_scale_cut, cov_methods)
         return summary
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/pseudo_cl.py b/src/sp_validation/cosmo_val/pseudo_cl.py
index 514049a6..cf74efc6 100644
--- a/src/sp_validation/cosmo_val/pseudo_cl.py
+++ b/src/sp_validation/cosmo_val/pseudo_cl.py
@@ -17,6 +17,7 @@
 from astropy.io import fits
 from cs_util.cosmo import get_theo_c_ell
 
+from .. import sacc_io
 from ..pseudo_cl import (
     apply_random_rotation,
     get_n_gal_map,
@@ -26,6 +27,47 @@
 )
 from ..rho_tau import get_params_rho_tau
 from ..statistics import chi2_and_pte, cov_from_one_covariance
+from .sacc_writers import BIN as SACC_BIN
+from .sacc_writers import pseudo_cl_to_sacc
+
+
+def plot_pseudo_cl_spectrum(datasets, spectrum, output_path):
+    """Two-panel ℓC_ℓ / C_ℓ figure for one spectrum across catalogue versions.
+
+    ``datasets`` maps a version to ``{"ell", "cl", "cov", "style"}``, where
+    ``style`` carries the ``marker`` and ``colour`` the version is drawn with.
+    """
+    fig, ax = plt.subplots(nrows=2, ncols=1, figsize=(8, 8))
+    minor_ticks = [i * 10 for i in range(1, 10)] + [i * 100 for i in range(1, 21)]
+
+    for panel, scaled in ((ax[0], True), (ax[1], False)):
+        for version, data in datasets.items():
+            ell, cl = np.asarray(data["ell"]), np.asarray(data["cl"])
+            err = np.sqrt(np.diag(np.asarray(data["cov"])))
+            style = data.get("style", {})
+            panel.errorbar(
+                ell,
+                ell * cl if scaled else cl,
+                yerr=ell * err if scaled else err,
+                fmt=style.get("marker", "."),
+                color=style.get("colour"),
+                label=f"{version} {spectrum}",
+                capsize=2 if scaled else None,
+            )
+        panel.set_ylabel(r"$\ell C_\ell$" if scaled else r"$C_\ell$")
+        panel.set_xlim(ell.min() - 10, ell.max() + 100)
+        panel.set_xscale("squareroot")
+        panel.set_xticks(np.array([100, 400, 900, 1600]))
+        panel.minorticks_on()
+        panel.tick_params(axis="x", which="minor", length=2, width=0.8)
+        panel.xaxis.set_ticks(minor_ticks, minor=True)
+
+    ax[1].set_xlabel(r"$\ell$")
+    ax[1].set_yscale("log")
+    plt.suptitle(f"Pseudo-Cl {spectrum} (Gaussian covariance)")
+    plt.legend()
+    plt.savefig(output_path)
+    plt.close(fig)
 
 
 class PseudoClMixin:
@@ -451,14 +493,23 @@ def calculate_pseudo_cl_g_ng_cov(self, gaussian_part="iNKA"):
             f"Done Gaussian and Non-Gaussian covariance of the Pseudo-Cl's using {gaussian_part} for the Gaussian part"
         )
 
-    def calculate_pseudo_cl(self):
+    def calculate_pseudo_cl(self, out_path=None):
         """
         Compute the pseudo-Cl of given catalogs.
+
+        ``out_path`` is the exact destination the part is born at (one version
+        only); ``None`` defaults each part to ``pseudo_cl_{ver}.sacc``.
         """
         self.print_start("Computing pseudo-Cl's")
 
         nside = self.nside
 
+        if out_path is not None and len(self.versions) != 1:
+            raise ValueError(
+                "calculate_pseudo_cl(out_path=...) writes one part to one path, "
+                f"but {len(self.versions)} versions are configured; call per version"
+            )
+
         try:
             self._pseudo_cls
         except AttributeError:
@@ -468,20 +519,32 @@ def calculate_pseudo_cl(self):
 
             self._pseudo_cls[ver] = {}
 
-            out_path = self._output_path(f"pseudo_cl_{ver}.fits")
-            if os.path.exists(out_path):
-                self.print_done(f"Skipping Pseudo-Cl's calculation, {out_path} exists")
-                cl_shear = fits.getdata(out_path)
-                self._pseudo_cls[ver]["pseudo_cl"] = cl_shear
+            ver_out_path = out_path or self._output_path(f"pseudo_cl_{ver}.sacc")
+            if os.path.exists(ver_out_path):
+                self.print_done(
+                    f"Skipping Pseudo-Cl's calculation, {ver_out_path} exists"
+                )
+                self._pseudo_cls[ver]["pseudo_cl"] = self._load_pseudo_cl_sacc(
+                    ver_out_path
+                )
             elif self.cell_method == "map":
-                self.calculate_pseudo_cl_map(ver, nside, out_path)
+                self.calculate_pseudo_cl_map(ver, nside, ver_out_path)
             elif self.cell_method == "catalog":
-                self.calculate_pseudo_cl_catalog(ver, out_path)
+                self.calculate_pseudo_cl_catalog(ver, ver_out_path)
             else:
                 raise ValueError(f"Unknown cell method: {self.cell_method}")
 
         self.print_done("Done pseudo-Cl's")
 
+    @staticmethod
+    def _load_pseudo_cl_sacc(out_path):
+        """Read a pseudo-Cl SACC part into the ELL/EE/EB/BB dict."""
+        # Readback of a part this producer just wrote — a legitimate pre-blind
+        # consumer, so the fail-closed load is opted out of.
+        s = sacc_io.load(out_path, allow_unblinded=True)
+        ell, ee, bb, eb, _window = sacc_io.get_pseudo_cl(s, SACC_BIN)
+        return {"ELL": ell, "EE": ee, "EB": eb, "BB": bb}
+
     def calculate_pseudo_cl_map(self, ver, nside, out_path):
         params = get_params_rho_tau(self.cc[ver], survey=ver)
 
@@ -547,10 +610,9 @@ def calculate_pseudo_cl_map(self, ver, nside, out_path):
         cl_shear = cl_shear - cl_noise
 
         self.print_cyan("Saving pseudo-Cl's...")
-        self.save_pseudo_cl(ell_eff, cl_shear, out_path)
+        self.pseudo_cl_to_sacc_part(ver, out_path, ell_eff, cl_shear, wsp)
 
-        cl_shear = fits.getdata(out_path)
-        self._pseudo_cls[ver]["pseudo_cl"] = cl_shear
+        self._pseudo_cls[ver]["pseudo_cl"] = self._load_pseudo_cl_sacc(out_path)
 
     def calculate_pseudo_cl_catalog(self, ver, out_path):
         params = get_params_rho_tau(self.cc[ver], survey=ver)
@@ -563,10 +625,9 @@ def calculate_pseudo_cl_catalog(self, ver, out_path):
         )
 
         self.print_cyan("Saving pseudo-Cl's...")
-        self.save_pseudo_cl(ell_eff, cl_shear, out_path)
+        self.pseudo_cl_to_sacc_part(ver, out_path, ell_eff, cl_shear, wsp)
 
-        cl_shear = fits.getdata(out_path)
-        self._pseudo_cls[ver]["pseudo_cl"] = cl_shear
+        self._pseudo_cls[ver]["pseudo_cl"] = self._load_pseudo_cl_sacc(out_path)
 
     def get_n_gal_map(self, params, nside, cat_gal):
         """Weighted galaxy number-density map (thin wrapper -> primitive)."""
@@ -655,223 +716,50 @@ def apply_random_rotation(self, e1, e2, rng=None):
         """
         return apply_random_rotation(e1, e2, rng)
 
-    def save_pseudo_cl(self, ell_eff, pseudo_cl, out_path):
-        """
-        Save pseudo-Cl's to a FITS file.
+    def pseudo_cl_to_sacc_part(self, version, out_path, ell_eff, cl_all, wsp):
+        """Write the pseudo-Cl SACC part (EE/BB/EB + shared bandpower window).
 
-        Parameters
-        ----------
-        pseudo_cl : np.array
-            Pseudo-Cl's to save.
-        out_path : str
-            Path to save the pseudo-Cl's to.
+        ``cl_all`` is NaMaster's decoupled ``(4, nbp)`` array (EE, EB, BE, BB);
+        the writer takes the shared bandpower window from ``wsp``. No covariance
+        is attached here.
         """
-        # Create columns of the fits file
-        col1 = fits.Column(name="ELL", format="D", array=ell_eff)
-        col2 = fits.Column(name="EE", format="D", array=pseudo_cl[0])
-        col3 = fits.Column(name="EB", format="D", array=pseudo_cl[1])
-        col4 = fits.Column(name="BB", format="D", array=pseudo_cl[3])
-        coldefs = fits.ColDefs([col1, col2, col3, col4])
-        cell_hdu = fits.BinTableHDU.from_columns(coldefs, name="PSEUDO_CELL")
-
-        cell_hdu.writeto(out_path, overwrite=True)
+        s = pseudo_cl_to_sacc(
+            self.sacc_nz(version),
+            self.sacc_metadata(version),
+            ell_eff,
+            cl_all,
+            wsp,
+        )
+        sacc_io.save(s, out_path, type="data")
 
     def plot_pseudo_cl(self):
-        """
-        Plot pseudo-Cl's for given catalogs.
-        """
+        """Plot the EE/EB/BB pseudo-Cl spectra for every version."""
         self.print_cyan("Plotting pseudo-Cl's")
 
-        # Plotting EE
-        out_path = self._output_path("cell_ee.png")
-        fig, ax = plt.subplots(nrows=2, ncols=1, figsize=(8, 8))
-
-        for ver in self.versions:
-            ell = self.pseudo_cls[ver]["pseudo_cl"]["ELL"]
-            cov = self.pseudo_cls[ver]["cov"]["COVAR_EE_EE"].data
-            ax[0].errorbar(
-                ell,
-                ell * self.pseudo_cls[ver]["pseudo_cl"]["EE"],
-                yerr=ell * np.sqrt(np.diag(cov)),
-                fmt=self.cc[ver]["marker"],
-                label=ver + " EE",
-                color=self.cc[ver]["colour"],
-                capsize=2,
-            )
-
-        ax[0].set_ylabel(r"$\ell C_\ell$")
-
-        ax[0].set_xlim(ell.min() - 10, ell.max() + 100)
-        ax[0].set_xscale("squareroot")
-        ax[0].set_xticks(np.array([100, 400, 900, 1600]))
-        ax[0].minorticks_on()
-        ax[0].tick_params(axis="x", which="minor", length=2, width=0.8)
-        minor_ticks = [i * 10 for i in range(1, 10)] + [i * 100 for i in range(1, 21)]
-        ax[0].xaxis.set_ticks(minor_ticks, minor=True)
-
-        for ver in self.versions:
-            ell = self.pseudo_cls[ver]["pseudo_cl"]["ELL"]
-            cov = self.pseudo_cls[ver]["cov"]["COVAR_EE_EE"].data
-            ax[1].errorbar(
-                ell,
-                self.pseudo_cls[ver]["pseudo_cl"]["EE"],
-                yerr=np.sqrt(np.diag(cov)),
-                fmt=self.cc[ver]["marker"],
-                label=ver + " EE",
-                color=self.cc[ver]["colour"],
-            )
-
-        ax[1].set_xlabel(r"$\ell$")
-        ax[1].set_ylabel(r"$C_\ell$")
-
-        ax[1].set_xlim(ell.min() - 10, ell.max() + 100)
-        ax[1].set_xscale("squareroot")
-        ax[1].set_yscale("log")
-        ax[1].set_xticks(np.array([100, 400, 900, 1600]))
-        ax[1].minorticks_on()
-        ax[1].tick_params(axis="x", which="minor", length=2, width=0.8)
-        minor_ticks = [i * 10 for i in range(1, 10)] + [i * 100 for i in range(1, 21)]
-        ax[1].xaxis.set_ticks(minor_ticks, minor=True)
-
-        plt.suptitle("Pseudo-Cl EE (Gaussian covariance)")
-        plt.legend()
-        plt.savefig(out_path)
-
-        # Plotting EB
-        out_path = self._output_path("cell_eb.png")
-
-        fig, ax = plt.subplots(nrows=2, ncols=1, figsize=(8, 8))
-
-        for ver in self.versions:
-            ell = self.pseudo_cls[ver]["pseudo_cl"]["ELL"]
-            cov = self.pseudo_cls[ver]["cov"]["COVAR_EB_EB"].data
-            ax[0].errorbar(
-                ell,
-                ell * self.pseudo_cls[ver]["pseudo_cl"]["EB"],
-                yerr=ell * np.sqrt(np.diag(cov)),
-                fmt=self.cc[ver]["marker"],
-                label=ver + " EB",
-                color=self.cc[ver]["colour"],
-                capsize=2,
+        for spectrum in ("EE", "EB", "BB"):
+            datasets = {
+                ver: {
+                    "ell": self.pseudo_cls[ver]["pseudo_cl"]["ELL"],
+                    "cl": self.pseudo_cls[ver]["pseudo_cl"][spectrum],
+                    "cov": self.pseudo_cls[ver]["cov"][
+                        f"COVAR_{spectrum}_{spectrum}"
+                    ].data,
+                    "style": {
+                        "marker": self.cc[ver]["marker"],
+                        "colour": self.cc[ver]["colour"],
+                    },
+                }
+                for ver in self.versions
+            }
+            plot_pseudo_cl_spectrum(
+                datasets, spectrum, self._output_path(f"cell_{spectrum.lower()}.png")
             )
 
-        ax[0].axhline(0, color="black", linestyle="--")
-        ax[0].set_ylabel(r"$\ell C_\ell$")
-
-        ax[0].set_xlim(ell.min() - 10, ell.max() + 100)
-        ax[0].set_xscale("squareroot")
-        ax[0].set_xticks(np.array([100, 400, 900, 1600]))
-        ax[0].minorticks_on()
-        ax[0].tick_params(axis="x", which="minor", length=2, width=0.8)
-        minor_ticks = [i * 10 for i in range(1, 10)] + [i * 100 for i in range(1, 21)]
-        ax[0].xaxis.set_ticks(minor_ticks, minor=True)
-
-        for ver in self.versions:
-            ell = self.pseudo_cls[ver]["pseudo_cl"]["ELL"]
-            cov = self.pseudo_cls[ver]["cov"]["COVAR_EB_EB"].data
-            ax[1].errorbar(
-                ell,
-                self.pseudo_cls[ver]["pseudo_cl"]["EB"],
-                yerr=np.sqrt(np.diag(cov)),
-                fmt=self.cc[ver]["marker"],
-                label=ver + " EB",
-                color=self.cc[ver]["colour"],
-            )
-
-        ax[1].set_xlabel(r"$\ell$")
-        ax[1].set_ylabel(r"$C_\ell$")
-
-        ax[1].set_xlim(ell.min() - 10, ell.max() + 100)
-        ax[1].set_xscale("squareroot")
-        ax[1].set_yscale("log")
-        ax[1].set_xticks(np.array([100, 400, 900, 1600]))
-        ax[1].minorticks_on()
-        ax[1].tick_params(axis="x", which="minor", length=2, width=0.8)
-        minor_ticks = [i * 10 for i in range(1, 10)] + [i * 100 for i in range(1, 21)]
-        ax[1].xaxis.set_ticks(minor_ticks, minor=True)
-
-        plt.suptitle("Pseudo-Cl EB (Gaussian covariance)")
-        plt.legend()
-        plt.savefig(out_path)
-
-        # Plotting BB
-        out_path = self._output_path("cell_bb.png")
-
-        fig, ax = plt.subplots(nrows=2, ncols=1, figsize=(8, 8))
-
-        for ver in self.versions:
-            ell = self.pseudo_cls[ver]["pseudo_cl"]["ELL"]
-            cov = self.pseudo_cls[ver]["cov"]["COVAR_BB_BB"].data
-            ax[0].errorbar(
-                ell,
-                ell * self.pseudo_cls[ver]["pseudo_cl"]["BB"],
-                yerr=ell * np.sqrt(np.diag(cov)),
-                fmt=self.cc[ver]["marker"],
-                label=ver + " BB",
-                color=self.cc[ver]["colour"],
-                capsize=2,
-            )
-
-        ax[0].axhline(0, color="black", linestyle="--")
-        ax[0].set_ylabel(r"$\ell C_\ell$")
-
-        ax[0].set_xlim(ell.min() - 10, ell.max() + 100)
-        ax[0].set_xscale("squareroot")
-        ax[0].set_xticks(np.array([100, 400, 900, 1600]))
-        ax[0].minorticks_on()
-        ax[0].tick_params(axis="x", which="minor", length=2, width=0.8)
-        minor_ticks = [i * 10 for i in range(1, 10)] + [i * 100 for i in range(1, 21)]
-        ax[0].xaxis.set_ticks(minor_ticks, minor=True)
-
-        for ver in self.versions:
-            ell = self.pseudo_cls[ver]["pseudo_cl"]["ELL"]
-            cov = self.pseudo_cls[ver]["cov"]["COVAR_BB_BB"].data
-            ax[1].errorbar(
-                ell,
-                self.pseudo_cls[ver]["pseudo_cl"]["BB"],
-                yerr=np.sqrt(np.diag(cov)),
-                fmt=self.cc[ver]["marker"],
-                label=ver + " BB",
-                color=self.cc[ver]["colour"],
-            )
-
-        ax[1].set_xlabel(r"$\ell$")
-        ax[1].set_ylabel(r"$C_\ell$")
-
-        ax[1].set_xlim(ell.min() - 10, ell.max() + 100)
-        ax[1].set_xscale("squareroot")
-        ax[1].set_yscale("log")
-        ax[1].set_xticks(np.array([100, 400, 900, 1600]))
-        ax[1].minorticks_on()
-        ax[1].tick_params(axis="x", which="minor", length=2, width=0.8)
-        minor_ticks = [i * 10 for i in range(1, 10)] + [i * 100 for i in range(1, 21)]
-        ax[1].xaxis.set_ticks(minor_ticks, minor=True)
-
-        plt.suptitle("Pseudo-Cl BB (Gaussian covariance)")
-        plt.legend()
-        plt.savefig(out_path)
-
-        # Print C_l^BB PTE for each version and save BB data
-        print("\nC_l^BB PTE summary:")
         for ver in self.versions:
             cl_bb = self.pseudo_cls[ver]["pseudo_cl"]["BB"]
             cov_bb = self.pseudo_cls[ver]["cov"]["COVAR_BB_BB"].data
             chi2_bb, _, pte_bb = chi2_and_pte(cl_bb, cov_bb)
-            chi2_bb = float(chi2_bb)
             print(
                 f"  {ver}: C_l^BB PTE = {pte_bb:.4f} "
-                f"(chi2/dof = {chi2_bb:.1f}/{len(cl_bb)})"
-            )
-
-            # Save BB data + covariance to .npz
-            ell = self.pseudo_cls[ver]["pseudo_cl"]["ELL"]
-            bb_out = self._output_path(f"{ver}_cell_bb_data.npz")
-            np.savez(
-                bb_out,
-                ell=ell,
-                cl_bb=cl_bb,
-                cov_bb=cov_bb,
-                chi2_bb=np.array(chi2_bb),
-                pte_bb=np.array(pte_bb),
+                f"(chi2/dof = {float(chi2_bb):.1f}/{len(cl_bb)})"
             )
-            print(f"  Saved BB data to {bb_out}")
diff --git a/src/sp_validation/cosmo_val/psf_systematics.py b/src/sp_validation/cosmo_val/psf_systematics.py
index 9f8a247e..b20d1085 100644
--- a/src/sp_validation/cosmo_val/psf_systematics.py
+++ b/src/sp_validation/cosmo_val/psf_systematics.py
@@ -16,10 +16,12 @@
 from shear_psf_leakage.rho_tau_stat import PSFErrorFit
 from uncertainties import ufloat
 
+from .. import sacc_io
 from ..rho_tau import (
     get_rho_tau_w_cov,
     get_samples,
 )
+from .sacc_writers import rho_tau_to_sacc
 
 
 class PSFSystematicsMixin:
@@ -41,11 +43,42 @@ def calculate_rho_tau_stats(self):
                 cov_rho=self.compute_cov_rho,
                 npatch=self.npatch,
             )
+            self.rho_tau_to_sacc_part(
+                ver, out_dir, base, rho_stat_handler, tau_stat_handler
+            )
         self.print_done("Rho stats finished")
 
         self._rho_stat_handler = rho_stat_handler
         self._tau_stat_handler = tau_stat_handler
 
+    def rho_tau_to_sacc_part(
+        self, version, out_dir, base, rho_stat_handler, tau_stat_handler
+    ):
+        """Write the ρ/τ SACC part for one version.
+
+        ρ_0…ρ_5 autos and τ_0/τ_2/τ_5 leakage from the handler tables. The
+        ``CovTauTh`` theory covariance ``cov_tau_{base}_th.npy`` is passed as
+        ``tau_cov_th`` when it exists; without it the τ block falls back —
+        loudly — to a variance diagonal.
+        """
+        tau_cov_path = os.path.join(out_dir, f"cov_tau_{base}_th.npy")
+        tau_cov_th = np.load(tau_cov_path) if os.path.exists(tau_cov_path) else None
+        if tau_cov_th is None:
+            self.print_magenta(
+                f"No τ theory covariance at {tau_cov_path}; writing ρ/τ SACC part "
+                "with a diagonal placeholder covariance (τ inference block is a "
+                "variance diagonal, not CovTauTh)."
+            )
+        s = rho_tau_to_sacc(
+            self.sacc_nz(version),
+            self.sacc_metadata(version),
+            rho_stat_handler.rho_stats,
+            tau_stat_handler.tau_stats,
+            tau_cov_th=tau_cov_th,
+        )
+        out_path = os.path.join(out_dir, f"rho_tau_{base}.sacc")
+        sacc_io.save(s, out_path, type="data")
+
     @property
     def rho_stat_handler(self):
         if not hasattr(self, "_rho_stat_handler"):
diff --git a/src/sp_validation/cosmo_val/pure_eb.py b/src/sp_validation/cosmo_val/pure_eb.py
index 71f5d76a..934eaa42 100644
--- a/src/sp_validation/cosmo_val/pure_eb.py
+++ b/src/sp_validation/cosmo_val/pure_eb.py
@@ -269,19 +269,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/cosmo_val/real_space.py b/src/sp_validation/cosmo_val/real_space.py
index 76d05d0a..1521ecc7 100644
--- a/src/sp_validation/cosmo_val/real_space.py
+++ b/src/sp_validation/cosmo_val/real_space.py
@@ -12,12 +12,11 @@
 import matplotlib.ticker as mticker
 import numpy as np
 import treecorr
-from astropy.io import fits
 from cs_util import plots as cs_plots
 
 
 class RealSpaceMixin:
-    def calculate_2pcf(self, ver, npatch=None, save_fits=False, **treecorr_config):
+    def calculate_2pcf(self, ver, npatch=None, **treecorr_config):
         """
         Calculate the two-point correlation function (2PCF) ξ± for a given catalog
         version with TreeCorr.
@@ -34,9 +33,6 @@ def calculate_2pcf(self, ver, npatch=None, save_fits=False, **treecorr_config):
             npatch (int, optional): The number of patches to use for the calculation.
             Defaults to the instance's `npatch` attribute.
 
-            save_fits (bool, optional): Whether to save the ξ± results to FITS files.
-            Defaults to False.
-
             **treecorr_config: Additional TreeCorr configuration parameters that will
             override the instance's default `treecorr_config`. For example, `min_sep=1`.
 
@@ -49,8 +45,7 @@ def calculate_2pcf(self, ver, npatch=None, save_fits=False, **treecorr_config):
               calculation is skipped, and the results are loaded from the file.
             - If a patch file for the given configuration does not exist, it is
               created during the process.
-            - FITS files for ξ+ and ξ− are saved with additional metadata in their
-              headers if `save_fits` is True.
+            - The ``.txt`` TreeCorr dump is the only raw byproduct written here.
         """
 
         self.print_magenta(f"Computing {ver} ξ±")
@@ -99,75 +94,13 @@ def calculate_2pcf(self, ver, npatch=None, save_fits=False, **treecorr_config):
 
             # Process the catalog & write the correlation functions
             gg.process(cat_gal)
-            gg.write(out_fname, write_patch_results=True, write_cov=True)
-
-        # Save xi_p and xi_m results to fits file
-        # (moved outside so it runs even if txt exists)
-        if save_fits:
-            lst = np.arange(1, treecorr_config["nbins"] + 1)
-
-            col1 = fits.Column(name="BIN1", format="K", array=np.ones(len(lst)))
-            col2 = fits.Column(name="BIN2", format="K", array=np.ones(len(lst)))
-            col3 = fits.Column(name="ANGBIN", format="K", array=lst)
-            col4 = fits.Column(name="VALUE", format="D", array=gg.xip)
-            col5 = fits.Column(name="ANG", format="D", unit="arcmin", array=gg.meanr)
-            coldefs = fits.ColDefs([col1, col2, col3, col4, col5])
-            xiplus_hdu = fits.BinTableHDU.from_columns(coldefs, name="XI_PLUS")
-
-            col4 = fits.Column(name="VALUE", format="D", array=gg.xim)
-            coldefs = fits.ColDefs([col1, col2, col3, col4, col5])
-            ximinus_hdu = fits.BinTableHDU.from_columns(coldefs, name="XI_MINUS")
-
-            # append xi_plus header info
-            xiplus_dict = {
-                "2PTDATA": "T",
-                "QUANT1": "G+R",
-                "QUANT2": "G+R",
-                "KERNEL_1": "NZ_SOURCE",
-                "KERNEL_2": "NZ_SOURCE",
-                "WINDOWS": "SAMPLE",
-            }
-            for key in xiplus_dict:
-                xiplus_hdu.header[key] = xiplus_dict[key]
-
-                col1 = fits.Column(name="BIN1", format="K", array=np.ones(len(lst)))
-                col2 = fits.Column(name="BIN2", format="K", array=np.ones(len(lst)))
-                col3 = fits.Column(name="ANGBIN", format="K", array=lst)
-                col4 = fits.Column(name="VALUE", format="D", array=gg.xip)
-                col5 = fits.Column(name="ANG", format="D", unit="arcmin", array=gg.rnom)
-                coldefs = fits.ColDefs([col1, col2, col3, col4, col5])
-                xiplus_hdu = fits.BinTableHDU.from_columns(coldefs, name="XI_PLUS")
-
-                col4 = fits.Column(name="VALUE", format="D", array=gg.xim)
-                coldefs = fits.ColDefs([col1, col2, col3, col4, col5])
-                ximinus_hdu = fits.BinTableHDU.from_columns(coldefs, name="XI_MINUS")
-
-                # append xi_plus header info
-                xiplus_dict = {
-                    "2PTDATA": "T",
-                    "QUANT1": "G+R",
-                    "QUANT2": "G+R",
-                    "KERNEL_1": "NZ_SOURCE",
-                    "KERNEL_2": "NZ_SOURCE",
-                    "WINDOWS": "SAMPLE",
-                }
-                for key in xiplus_dict:
-                    xiplus_hdu.header[key] = xiplus_dict[key]
-            # Use same naming format as txt output
-            fits_base = out_fname.replace(".txt", "").replace("_xi_", "_")
-            xiplus_hdu.writeto(
-                f"{fits_base.replace(ver, f'xi_plus_{ver}')}.fits",
-                overwrite=True,
-            )
-
-            # append xi_minus header info
-            ximinus_dict = {**xiplus_dict, "QUANT1": "G-R", "QUANT2": "G-R"}
-            for key in ximinus_dict:
-                ximinus_hdu.header[key] = ximinus_dict[key]
-            ximinus_hdu.writeto(
-                f"{fits_base.replace(ver, f'xi_minus_{ver}')}.fits",
-                overwrite=True,
-            )
+            # Never write_patch_results: a per-patch ξ± realisation is an
+            # unblinded data vector, and nothing downstream reads one — the
+            # covariance a consumer needs is the matrix, which the SACC part
+            # carries. The .txt keeps the matrix only where there are patches to
+            # estimate it from; at npatch=1 var_method is "shot" and it would add
+            # nothing over the varxip/varxim columns.
+            gg.write(out_fname, write_patch_results=False, write_cov=int(npatch) > 1)
 
         # Add correlation object to class
         if not hasattr(self, "cat_ggs"):
diff --git a/src/sp_validation/cosmo_val/sacc_writers.py b/src/sp_validation/cosmo_val/sacc_writers.py
new file mode 100644
index 00000000..8b7f57b7
--- /dev/null
+++ b/src/sp_validation/cosmo_val/sacc_writers.py
@@ -0,0 +1,237 @@
+"""Born-as-SACC writers for the cosmo_val data products.
+
+A thin, pure layer between the ``cosmo_val`` mixins (which compute statistics as
+TreeCorr / NaMaster / b_modes arrays) and :mod:`sp_validation.sacc_io` (which
+knows the file layout). Each ``*_to_sacc`` function turns one already-computed
+statistic into a single-statistic SACC — a *part* — carrying that statistic's
+own covariance as its one block. :func:`assemble_analysis_sacc` rebuilds the
+single ``{version}.sacc`` analysis file from these parts.
+
+Everything here is single-bin today (``bins=(0, 0)``); the interface is
+tomography-native so a future round supplies real bin pairs unchanged.
+"""
+
+import numpy as np
+import sacc
+
+from .. import sacc_io as sio
+from ..pseudo_cl import bandpower_window_from_workspace
+
+# Statistics carried in the analysis file, and their custom-type k indices.
+RHO_K = range(6)  # ρ_0 … ρ_5
+TAU_K = (0, 2, 5)  # τ_0, τ_2, τ_5
+
+# NaMaster spin-2 × spin-2 decoupled-spectrum row order (EE, EB, BE, BB).
+_NMT_EE, _NMT_EB, _NMT_BB = 0, 1, 3
+
+BIN = (0, 0)
+
+
+def xi_to_sacc(
+    nz,
+    metadata,
+    theta,
+    xip,
+    xim,
+    *,
+    grid,
+    theta_nom=None,
+    npairs=None,
+    weight=None,
+    variances=None,
+    covariance=None,
+):
+    """One ξ± part (``bins=(0, 0)``) on a named angular grid.
+
+    The grid's covariance comes in one of two shapes: ``covariance``, the dense
+    ``[ξ+; ξ−]``-ordered block (a jackknife estimate), or ``variances``, the
+    concatenated ``[varxip; varxim]`` diagonal. At most one may be given.
+    """
+    s = sio.new_sacc(nz, metadata)
+    sio.add_xi(
+        s,
+        BIN,
+        theta,
+        xip,
+        xim,
+        grid=grid,
+        theta_nom=theta_nom,
+        npairs=npairs,
+        weight=weight,
+    )
+    if covariance is not None and variances is not None:
+        raise ValueError("give xi_to_sacc a dense covariance or variances, not both")
+    if covariance is not None:
+        s.add_covariance(np.asarray(covariance))
+    elif variances is not None:
+        sio.add_diagonal_covariance(s, np.asarray(variances))
+    return s
+
+
+def pseudo_cl_to_sacc(nz, metadata, ell_eff, cl_all, wsp, covariance=None):
+    """One pseudo-Cℓ part: EE/BB/EB with the shared bandpower window.
+
+    ``cl_all`` is NaMaster's decoupled ``(4, nbp)`` array (EE, EB, BE, BB); the
+    window comes from :func:`bandpower_window_from_workspace`. ``covariance``,
+    when given, is the dense ``[EE; BB; EB]``-ordered block matching insertion.
+    """
+    window_ells, window_weights = bandpower_window_from_workspace(wsp)
+    s = sio.new_sacc(nz, metadata)
+    sio.add_pseudo_cl(
+        s,
+        BIN,
+        ell_eff,
+        cl_all[_NMT_EE],
+        cl_all[_NMT_BB],
+        cl_all[_NMT_EB],
+        window_ells=window_ells,
+        window_weights=window_weights,
+    )
+    if covariance is not None:
+        s.add_covariance(np.asarray(covariance))
+    return s
+
+
+def cosebis_to_sacc(nz, metadata, result, scale_cut):
+    """One COSEBIs part at the fiducial scale cut.
+
+    ``result`` is a single scale-cut result dict from
+    ``b_modes.calculate_cosebis`` — ``{"En", "Bn", "cov", ...}`` — where ``cov``
+    is the ``[En; Bn]``-ordered COSEBIs covariance.
+    """
+    s = sio.new_sacc(nz, metadata)
+    sio.add_cosebis(s, BIN, result["En"], scale_cut, Bn=result["Bn"])
+    s.add_covariance(np.asarray(result["cov"]))
+    return s
+
+
+def pure_eb_to_sacc(nz, metadata, theta, eb, covariance=None):
+    """One pure-E/B part: the six ``sacc_io.PURE_KEYS`` blocks.
+
+    ``eb`` is a mapping with the six keys (``xip_E`` … ``xim_amb``); each array
+    is sampled at ``theta``. ``covariance``, when given, is the dense block in
+    ``PURE_KEYS`` order (matching ``b_modes._EB_KEYS`` and the insertion order).
+    """
+    s = sio.new_sacc(nz, metadata)
+    sio.add_pure_eb(s, BIN, theta, **{key: eb[key] for key in sio.PURE_KEYS})
+    if covariance is not None:
+        s.add_covariance(np.asarray(covariance))
+    return s
+
+
+def rho_tau_to_sacc(nz, metadata, rho_stats, tau_stats, tau_cov_th=None):
+    """One ρ/τ part: ρ_0…ρ_5 autos and τ_0/τ_2/τ_5 leakage.
+
+    ``rho_stats`` / ``tau_stats`` are the ``shear_psf_leakage`` handler tables
+    (columns ``theta``, ``rho_{k}_p``, ``varrho_{k}_p``, … and the τ analogue).
+    ρ carries a ``varrho`` diagonal; τ carries a ``vartau`` diagonal, with
+    ``tau_cov_th`` — a ``(3·nbin, 3·nbin)`` k-major matrix over the τ-plus points
+    only — scattered into the τ-plus rows/columns when given. ``tau_cov_th=None``
+    leaves the τ block fully diagonal.
+    """
+    s = sio.new_sacc(nz, metadata)
+    theta_rho = np.asarray(rho_stats["theta"])
+    for k in RHO_K:
+        sio.add_rho(
+            s,
+            k,
+            theta_rho,
+            np.asarray(rho_stats[f"rho_{k}_p"]),
+            np.asarray(rho_stats[f"rho_{k}_m"]),
+        )
+    theta_tau = np.asarray(tau_stats["theta"])
+    for k in TAU_K:
+        sio.add_tau(
+            s,
+            BIN,
+            k,
+            theta_tau,
+            np.asarray(tau_stats[f"tau_{k}_p"]),
+            np.asarray(tau_stats[f"tau_{k}_m"]),
+        )
+    nbin = len(theta_tau)
+    rho_var = np.concatenate(
+        [
+            np.concatenate([rho_stats[f"varrho_{k}_p"], rho_stats[f"varrho_{k}_m"]])
+            for k in RHO_K
+        ]
+    )
+    tau_var = np.concatenate(
+        [
+            np.concatenate([tau_stats[f"vartau_{k}_p"], tau_stats[f"vartau_{k}_m"]])
+            for k in TAU_K
+        ]
+    )
+    if tau_cov_th is None:
+        s.add_covariance(np.concatenate([rho_var, tau_var]))
+        return s
+    tau_cov_th = np.asarray(tau_cov_th)
+    n_plus = len(TAU_K) * nbin
+    if tau_cov_th.shape != (n_plus, n_plus):
+        raise ValueError(
+            f"tau_cov_th shape {tau_cov_th.shape} does not match the "
+            f"{n_plus} τ-plus points ({len(TAU_K)} indices × {nbin} bins) — "
+            "CovTauTh.build_cov returns one (plus-folded) component per τ index"
+        )
+    n_rho, n_tau = len(rho_var), len(tau_var)
+    tau_block = np.diag(tau_var)
+    # τ-plus local positions in the τ block, k-major (per-k layout is [+; −]).
+    plus = np.concatenate(
+        [np.arange(2 * i * nbin, 2 * i * nbin + nbin) for i in range(len(TAU_K))]
+    )
+    tau_block[np.ix_(plus, plus)] = tau_cov_th
+    full = np.zeros((n_rho + n_tau, n_rho + n_tau))
+    full[:n_rho, :n_rho] = np.diag(rho_var)
+    full[n_rho:, n_rho:] = tau_block
+    s.add_covariance(full)
+    return s
+
+
+# --------------------------------------------------------------------------- #
+# Analysis-file assembly
+# --------------------------------------------------------------------------- #
+def _copy_data_points(dst, src):
+    """Append every data point of ``src`` into ``dst`` (tags preserved)."""
+    for dp in src.data:
+        dst.add_data_point(dp.data_type, dp.tracers, dp.value, **dp.tags)
+
+
+def assemble_analysis_sacc(parts):
+    """Rebuild the single ``{version}.sacc`` analysis file from parts.
+
+    Each part is a single-statistic Sacc (from a ``*_to_sacc`` writer, loaded
+    from disk) carrying its own covariance = its block. Tracers and metadata are
+    seeded from ``parts[0]`` (every part describes the same catalogue version).
+    Data points are re-added in the order the parts are given, which must be the
+    canonical order (ξ± reporting, pseudo-Cℓ, COSEBIs, pure-E/B, ρ/τ), and the
+    per-part blocks become one ``BlockDiagonalCovariance``. Insertion order and
+    block order therefore agree by construction — validated by
+    :func:`sp_validation.sacc_io.assemble_covariance`.
+
+    Parameters
+    ----------
+    parts : sequence of sacc.Sacc
+        Single-statistic parts, each with a covariance, in canonical order.
+
+    Returns
+    -------
+    sacc.Sacc
+        The analysis Sacc with a ``BlockDiagonalCovariance`` covering every point.
+    """
+    s = sacc.Sacc()
+    s.tracers.update(parts[0].tracers)
+    s.metadata.update(parts[0].metadata)
+    blocks = []
+    cursor = 0
+    for part in parts:
+        if part.covariance is None:
+            raise ValueError(
+                "every analysis part must carry its own covariance block; "
+                f"a part with data types {sorted(set(dp.data_type for dp in part.data))} "
+                "has none"
+            )
+        n = len(part.mean)
+        _copy_data_points(s, part)
+        blocks.append((np.arange(cursor, cursor + n), part.covariance.dense))
+        cursor += n
+    return sio.assemble_covariance(s, blocks)
diff --git a/src/sp_validation/pseudo_cl.py b/src/sp_validation/pseudo_cl.py
index c9355ec9..34cbc68d 100644
--- a/src/sp_validation/pseudo_cl.py
+++ b/src/sp_validation/pseudo_cl.py
@@ -280,3 +280,32 @@ def get_pseudo_cls_catalog(
     cl_all = wsp.decouple_cell(cl_coupled)
 
     return ell_eff, cl_all, wsp
+
+
+# NaMaster spin-2 × spin-2 spectrum order: EE, EB, BE, BB.
+_NMT_EE = 0
+
+
+def bandpower_window_from_workspace(wsp):
+    """Extract the bandpower window matrix ``W`` for a spin-2×spin-2 workspace.
+
+    NaMaster's ``get_bandpower_windows()`` returns a four-index array
+    ``(n_cl_out, n_bpw, n_cl_in, n_ell)`` describing how each output bandpower
+    is built from the input multipoles across the EE/EB/BE/BB spectra. SACC's
+    ``BandpowerWindow`` model (one window per bandpower, shared across the
+    stored spectra) needs the per-spectrum *decoupling* window, i.e. the
+    diagonal EE←EE block (equal to BB←BB and EB←EB, verified identical).
+
+    Returns
+    -------
+    window_ells : np.ndarray
+        Multipoles the window spans, ``arange(n_ell)`` — the ``ell`` axis of
+        ``compute_coupled_cell``.
+    window_weights : np.ndarray
+        ``W`` of shape ``(n_ell, n_bpw)`` — one column per bandpower, the layout
+        :func:`sp_validation.sacc_io.add_pseudo_cl` expects.
+    """
+    bpw = wsp.get_bandpower_windows()  # (n_cl_out, n_bpw, n_cl_in, n_ell)
+    diagonal = bpw[_NMT_EE, :, _NMT_EE, :]  # (n_bpw, n_ell)
+    window_ells = np.arange(diagonal.shape[1], dtype=float)
+    return window_ells, diagonal.T
diff --git a/src/sp_validation/tests/test_assemble_sacc.py b/src/sp_validation/tests/test_assemble_sacc.py
new file mode 100644
index 00000000..f17d5468
--- /dev/null
+++ b/src/sp_validation/tests/test_assemble_sacc.py
@@ -0,0 +1,300 @@
+"""Integration tests for the ``assemble_sacc.py`` workflow script.
+
+The pure assembler (``sacc_writers.assemble_analysis_sacc``) is covered in
+``test_sacc_writers.py``. This file exercises the *script seam* the DAG uses:
+``assemble_sacc.assemble_sacc`` loads per-statistic ``.sacc`` part *files* in
+CANONICAL order, injects the born-cov-less ξ± / pseudo-Cℓ blocks from the real
+CosmoCov ``.txt`` and NaMaster covariance FITS, and writes one
+``{version}.sacc`` whose points and covariance blocks land in canonical order.
+
+The script lives under ``workflow/scripts`` (off the package path); it is loaded
+by file path exactly as the lightcone/ASTRA CLI path imports it.
+"""
+
+import importlib.util
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+from sp_validation import sacc_io as sio
+from sp_validation.cosmo_val import sacc_writers as sw
+
+
+def _load_assemble_module():
+    """Import ``workflow/scripts/assemble_sacc.py`` by file path."""
+    repo_root = next(
+        p for p in Path(__file__).resolve().parents if (p / "pyproject.toml").exists()
+    )
+    path = repo_root / "workflow" / "scripts" / "assemble_sacc.py"
+    spec = importlib.util.spec_from_file_location("assemble_sacc", path)
+    module = importlib.util.module_from_spec(spec)
+    spec.loader.exec_module(module)
+    return module
+
+
+asm = _load_assemble_module()
+
+
+def _nz(seed=0, n=40):
+    rng = np.random.default_rng(seed)
+    return np.linspace(0.01, 2.0, n), rng.uniform(0.1, 1.0, n)
+
+
+def _spd(n, seed):
+    a = np.random.default_rng(seed).normal(size=(n, n))
+    return a @ a.T + n * np.eye(n)
+
+
+def _theta(n=6):
+    return np.geomspace(1.0, 100.0, n)
+
+
+META = {"catalogue_version": "vSYNTH", "npatch": 1}
+
+
+def _xi_cov_txt(tmp_path, n=12, seed=21):
+    """A CosmoCov-format ξ± covariance: the dense ``_processed.txt`` matrix.
+
+    Same writer/format ``covariance_process`` emits and ``--xi-cov`` reads, at
+    the synthetic parts' 12-point ([ξ+; ξ−] over 6 θ) size. Returns
+    ``(path, matrix)``.
+    """
+    cov = _spd(n, seed)
+    path = tmp_path / "xi_cov_processed.txt"
+    np.savetxt(str(path), cov)
+    return str(path), cov
+
+
+def _pseudo_cl_cov_fits(tmp_path, n=3):
+    """A NaMaster covariance FITS: one HDU per spectrum. Returns (path, blocks)."""
+    from astropy.io import fits
+
+    blocks = {"EE": _spd(n, 31), "BB": _spd(n, 32), "EB": _spd(n, 33)}
+    path = tmp_path / "pseudo_cl_cov.fits"
+    fits.HDUList(
+        [fits.PrimaryHDU()]
+        + [fits.ImageHDU(block, name=f"COVAR_{k}_{k}") for k, block in blocks.items()]
+    ).writeto(str(path))
+    return str(path), blocks
+
+
+def _write_parts(tmp_path, *, with_pseudo_cl=True, cov_less=("xi_reporting",)):
+    """Write per-statistic parts to disk; return the ``{name: path}`` mapping.
+
+    Parts named in ``cov_less`` are written without a covariance (mimicking the
+    born-cov-less ξ± reporting / pseudo-Cℓ parts); the rest carry their own block.
+    """
+    nz = {0: _nz()}
+    theta = _theta()
+    ell = np.array([30.0, 60.0, 90.0])
+
+    class _Wsp:
+        def get_bandpower_windows(self):
+            w = np.zeros((4, 3, 4, 20))
+            for out in range(4):
+                for b in range(3):
+                    w[out, b, out, b * 6 : b * 6 + 6] = 1.0
+            return w
+
+    xi = sw.xi_to_sacc(
+        nz, META, theta, np.arange(6) * 1e-5, np.arange(6) * 2e-5, grid="reporting"
+    )
+    if "xi_reporting" not in cov_less:
+        xi.add_covariance(_spd(len(xi.mean), 1))
+
+    cl_all = np.vstack(
+        [np.arange(3) * 1e-9, np.arange(3) * 2e-9, np.zeros(3), np.arange(3) * 3e-9]
+    )
+    cl = sw.pseudo_cl_to_sacc(
+        nz,
+        META,
+        ell,
+        cl_all,
+        _Wsp(),
+        covariance=None if "pseudo_cl" in cov_less else _spd(9, 2),
+    )
+
+    co = sw.cosebis_to_sacc(
+        nz,
+        META,
+        {
+            "En": np.arange(1, 6) * 1e-6,
+            "Bn": np.arange(1, 6) * 1e-7,
+            "cov": _spd(10, 3),
+        },
+        (1.0, 100.0),
+    )
+
+    eb_arrays = {
+        key: np.arange(6) * (i + 1) * 1e-6 for i, key in enumerate(sio.PURE_KEYS)
+    }
+    eb = sw.pure_eb_to_sacc(nz, META, theta, eb_arrays, covariance=_spd(36, 4))
+
+    rho = {"theta": theta}
+    tau = {"theta": theta}
+    rng = np.random.default_rng(5)
+    for k in sw.RHO_K:
+        for suffix in ("p", "m"):
+            rho[f"rho_{k}_{suffix}"] = rng.normal(size=6) * 1e-6
+            rho[f"varrho_{k}_{suffix}"] = rng.uniform(1e-14, 1e-13, 6)
+    for k in sw.TAU_K:
+        for suffix in ("p", "m"):
+            tau[f"tau_{k}_{suffix}"] = rng.normal(size=6) * 1e-6
+            tau[f"vartau_{k}_{suffix}"] = rng.uniform(1e-14, 1e-13, 6)
+    rt = sw.rho_tau_to_sacc(nz, META, rho, tau)
+
+    parts = {
+        "xi_reporting": xi,
+        "pseudo_cl": cl,
+        "cosebis": co,
+        "pure_eb": eb,
+        "rho_tau": rt,
+    }
+    if not with_pseudo_cl:
+        parts.pop("pseudo_cl")
+
+    paths = {}
+    for name, part in parts.items():
+        p = tmp_path / f"{name}.sacc"
+        sio.save(part, str(p), type="mock")
+        paths[name] = str(p)
+    return paths
+
+
+def test_assemble_sacc_canonical_order(tmp_path):
+    """Every point is covered and the blocks land in canonical order
+    (ξ±, pseudo-Cℓ, COSEBIs, pure-E/B, ρ, τ)."""
+    paths = _write_parts(tmp_path, cov_less=("xi_reporting",))
+    cov_path, xi_cov = _xi_cov_txt(tmp_path)
+    cl_cov_path, _blocks = _pseudo_cl_cov_fits(tmp_path)
+    out = tmp_path / "vSYNTH.sacc"
+    s = asm.assemble_sacc(
+        "vSYNTH", paths, str(out), xi_cov=cov_path, pseudo_cl_cov=cl_cov_path
+    )
+    assert out.exists()
+    assert type(s.covariance).__name__ == "BlockDiagonalCovariance"
+    assert s.covariance.dense.shape == (len(s.mean), len(s.mean))
+
+    # Canonical insertion order: the first data types are ξ+ then ξ−.
+    types_in_order = [dp.data_type for dp in s.data]
+    assert types_in_order[0] == sio.XI_PLUS
+    assert sio.XI_MINUS in types_in_order
+    # ξ appears before pseudo-Cℓ before COSEBIs before pure-E/B before ρ/τ.
+    first = {t: types_in_order.index(t) for t in set(types_in_order)}
+    assert first[sio.XI_PLUS] < first[sio.CL_EE] < first[sio.COSEBI_EE]
+    assert first[sio.COSEBI_EE] < first[sio.PURE_TYPES["xip_E"]]
+    assert first[sio.PURE_TYPES["xip_E"]] < first[sio.RHO_PLUS.format(k=0)]
+    assert first[sio.RHO_PLUS.format(k=0)] < first[sio.TAU_PLUS.format(k=0)]
+
+    # The ξ± block is the injected CosmoCov matrix on its own points.
+    tr = ("source_0", "source_0")
+    xi_idx = np.concatenate([s.indices(sio.XI_PLUS, tr), s.indices(sio.XI_MINUS, tr)])
+    dense = s.covariance.dense
+    assert np.allclose(dense[np.ix_(xi_idx, xi_idx)], xi_cov)
+    # ...and it does not bleed into the neighbouring COSEBIs block (cross zero).
+    co_idx = np.concatenate(
+        [s.indices(sio.COSEBI_EE, tr), s.indices(sio.COSEBI_BB, tr)]
+    )
+    assert np.allclose(dense[np.ix_(xi_idx, co_idx)], 0.0)
+
+
+def test_injected_xi_covariance_replaces_the_parts_own(tmp_path):
+    """The analytic ξ± covariance wins over the estimate the part was born with.
+
+    The reporting part carries the jackknife it was measured with — useful as a
+    diagnostic, but the analysis file takes the CosmoCov block.
+    """
+    paths = _write_parts(tmp_path, cov_less=())  # ξ± born with its own jackknife
+    cov_path, xi_cov = _xi_cov_txt(tmp_path)
+    cl_cov_path, _blocks = _pseudo_cl_cov_fits(tmp_path)
+
+    born = sio.load(paths["xi_reporting"], allow_unblinded=True).covariance.dense
+    assert not np.allclose(born, xi_cov)  # the two are distinguishable
+
+    out = tmp_path / "vSYNTH.sacc"
+    s = asm.assemble_sacc(
+        "vSYNTH", paths, str(out), xi_cov=cov_path, pseudo_cl_cov=cl_cov_path
+    )
+    tr = ("source_0", "source_0")
+    xi_idx = np.concatenate([s.indices(sio.XI_PLUS, tr), s.indices(sio.XI_MINUS, tr)])
+    assert np.allclose(s.covariance.dense[np.ix_(xi_idx, xi_idx)], xi_cov)
+
+
+def test_assemble_sacc_injects_pseudo_cl_covariance(tmp_path):
+    """The NaMaster cov FITS (COVAR_EE_EE/BB_BB/EB_EB) → block-diagonal pseudo-Cℓ
+    block, beside the injected CosmoCov ξ± block (the live default)."""
+    paths = _write_parts(tmp_path, cov_less=("xi_reporting", "pseudo_cl"))
+    cov_path, xi_cov = _xi_cov_txt(tmp_path)
+    # pseudo-Cℓ part is 3 ell × {EE, BB, EB} = 9 points; per-spectrum 3×3 blocks.
+    cov_fits, blocks = _pseudo_cl_cov_fits(tmp_path)
+    ee, bb, eb = blocks["EE"], blocks["BB"], blocks["EB"]
+
+    out = tmp_path / "vSYNTH.sacc"
+    s = asm.assemble_sacc(
+        "vSYNTH", paths, str(out), xi_cov=cov_path, pseudo_cl_cov=cov_fits
+    )
+    tr = ("source_0", "source_0")
+    cl_idx = np.concatenate(
+        [s.indices(sio.CL_EE, tr), s.indices(sio.CL_BB, tr), s.indices(sio.CL_EB, tr)]
+    )
+    dense = s.covariance.dense
+    expected = np.zeros((9, 9))
+    expected[0:3, 0:3], expected[3:6, 3:6], expected[6:9, 6:9] = ee, bb, eb
+    assert np.allclose(dense[np.ix_(cl_idx, cl_idx)], expected)
+    # ξ± carries its own CosmoCov block; the two don't bleed into each other.
+    xi_idx = np.concatenate([s.indices(sio.XI_PLUS, tr), s.indices(sio.XI_MINUS, tr)])
+    assert np.allclose(dense[np.ix_(xi_idx, xi_idx)], xi_cov)
+    assert np.allclose(dense[np.ix_(xi_idx, cl_idx)], 0.0)
+
+
+def test_missing_injected_covariance_raises(tmp_path):
+    """A statistic whose covariance is external cannot fall back to its own."""
+    paths = _write_parts(tmp_path, cov_less=())  # every part born with a block
+    out = tmp_path / "vSYNTH.sacc"
+    with pytest.raises(ValueError, match="takes its analysis covariance from"):
+        asm.assemble_sacc("vSYNTH", paths, str(out))
+
+
+def test_assemble_sacc_respects_pseudo_cl_toggle(tmp_path):
+    """With pseudo_cl absent, assembly still succeeds and omits the Cℓ points."""
+    paths = _write_parts(tmp_path, with_pseudo_cl=False, cov_less=("xi_reporting",))
+    assert "pseudo_cl" not in paths
+    cov_path, _xi_cov = _xi_cov_txt(tmp_path)
+    out = tmp_path / "vSYNTH.sacc"
+    s = asm.assemble_sacc("vSYNTH", paths, str(out), xi_cov=cov_path)
+    tr = ("source_0", "source_0")
+    assert len(s.indices(sio.CL_EE, tr)) == 0
+    # Round-trips as a valid BlockDiagonalCovariance over the remaining points.
+    s2 = sio.load(str(out))
+    assert type(s2.covariance).__name__ == "BlockDiagonalCovariance"
+    assert s2.covariance.dense.shape == (len(s2.mean), len(s2.mean))
+
+
+def test_assemble_sacc_expected_part_missing_raises(tmp_path):
+    """A typo'd input keyword drops a part from part_paths; the expected list
+    catches it rather than silently omitting the statistic."""
+    paths = _write_parts(tmp_path, cov_less=("xi_reporting",))
+    # Simulate a rule-input typo: cosebis wired under the wrong key.
+    paths["cosebi"] = paths.pop("cosebis")
+    cov_path, _xi_cov = _xi_cov_txt(tmp_path)
+    out = tmp_path / "vSYNTH.sacc"
+    with pytest.raises(ValueError, match="expected parts \\['cosebis'\\] missing"):
+        asm.assemble_sacc(
+            "vSYNTH",
+            paths,
+            str(out),
+            expected=["xi_reporting", "pseudo_cl", "cosebis", "pure_eb", "rho_tau"],
+            xi_cov=cov_path,
+        )
+
+
+def test_assemble_sacc_expected_rejects_unknown_name(tmp_path):
+    """A typo in the expected list itself is rejected (not a valid statistic)."""
+    paths = _write_parts(tmp_path, cov_less=("xi_reporting",))
+    cov_path, _xi_cov = _xi_cov_txt(tmp_path)
+    out = tmp_path / "vSYNTH.sacc"
+    with pytest.raises(ValueError, match="not assemblable statistics"):
+        asm.assemble_sacc(
+            "vSYNTH", paths, str(out), expected=["cosebi"], xi_cov=cov_path
+        )
diff --git a/src/sp_validation/tests/test_b_modes.py b/src/sp_validation/tests/test_b_modes.py
index e8e05c39..f51913fc 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,144 @@ 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)
diff --git a/src/sp_validation/tests/test_bmodes_workflow_dry_run.py b/src/sp_validation/tests/test_bmodes_workflow_dry_run.py
index 68c09981..3b7eecae 100644
--- a/src/sp_validation/tests/test_bmodes_workflow_dry_run.py
+++ b/src/sp_validation/tests/test_bmodes_workflow_dry_run.py
@@ -1,11 +1,15 @@
-"""Back-pressure guard #2: the B-modes Snakemake workflow dry-runs.
+"""Back-pressure guard #2: the paper Snakemake workflows dry-run.
 
-The reorg is allowed to change the rule graph; this guard only asserts that
-Snakemake can still parse the workflow and construct a dry run.
+The reorg is allowed to change the rule graph; these guards only assert that
+Snakemake can still parse each composed workflow and construct a dry run. One
+guard covers papers/bmodes (config space, no cosmo_val block); a second covers
+papers/cosmo_val, whose config DOES carry a cosmo_val block — so it is the only
+one that includes cosmo_val.smk and hence the born-as-SACC + assemble rules.
 """
 
 import os
 import subprocess
+import sys
 from pathlib import Path
 
 import pytest
@@ -26,31 +30,59 @@ def _repo_root() -> Path:
     raise RuntimeError("could not locate repo root (no pyproject.toml above test)")
 
 
-@requires_candide_data
-def test_bmodes_workflow_dry_runs():
-    """The paper B-mode workflow must still parse and dry-run cleanly."""
-    workflow_dir = _repo_root() / "papers/bmodes"
-    # PYTHONUNBUFFERED satisfies the Snakefile's `envvars:` declaration without
-    # depending on the invoking shell's environment.
+def _dry_run(workflow_dir, targets, *extra_snakemake_args):
+    """Construct a dry run of the paper workflow at ``workflow_dir``.
+
+    PYTHONUNBUFFERED satisfies the Snakefile's ``envvars:`` declaration. A dry
+    run never dispatches jobs, so any inherited SNAKEMAKE_PROFILE is dropped
+    rather than requiring its executor plugin. snakemake is invoked through
+    sys.executable, since a bare python3.12 may resolve off PATH to an
+    interpreter without it.
+    """
     env = os.environ | {"PYTHONNOUSERSITE": "1", "PYTHONUNBUFFERED": "1"}
-    result = subprocess.run(
+    env.pop("SNAKEMAKE_PROFILE", None)
+    return subprocess.run(
         [
-            "python3.12",
+            sys.executable,
             "-m",
             "snakemake",
-            "all_tapestry",
+            *targets,
             "--dry-run",
             "--cores",
             "1",
             "--configfile",
             "config/config.yaml",
+            *extra_snakemake_args,
         ],
         cwd=workflow_dir,
         env=env,
         text=True,
         stdout=subprocess.PIPE,
         stderr=subprocess.STDOUT,
-        timeout=60,
+        timeout=120,
         check=False,
     )
+
+
+@requires_candide_data
+def test_bmodes_workflow_dry_runs():
+    """The paper B-mode workflow must still parse and dry-run cleanly."""
+    result = _dry_run(_repo_root() / "papers/bmodes", ["all_tapestry"])
+    assert result.returncode == 0, result.stdout
+
+
+@requires_candide_data
+def test_cosmo_val_workflow_assemble_dry_runs():
+    """The cosmo_val workflow (the only one including cosmo_val.smk) resolves the
+    born-as-SACC + assemble DAG, and assemble pulls the tagged pseudo-Cl + cov."""
+    version = "SP_v1.4.6.3_leak_corr"
+    result = _dry_run(_repo_root() / "papers/cosmo_val", ["assemble_sacc_all"])
     assert result.returncode == 0, result.stdout
+    # assemble_sacc must pull the tagged pseudo-Cl part + its NaMaster
+    # covariance (not the untagged cv_pseudo_cl diagnostic), plus every part.
+    out = result.stdout
+    assert "rule assemble_sacc:" in out, out
+    assert f"pseudo_cl_{version}_blind=A_powspace_nbins=32.sacc" in out, out
+    assert f"pseudo_cl_cov_{version}_blind=A_powspace_nbins=32.fits" in out, out
+    for part in ("_xi_minsep=", "_cosebis.sacc", "_pure_eb.sacc", "rho_tau_"):
+        assert part in out, f"missing {part} part in assemble DAG:\n{out}"
diff --git a/src/sp_validation/tests/test_cli_seams.py b/src/sp_validation/tests/test_cli_seams.py
new file mode 100644
index 00000000..ce601d25
--- /dev/null
+++ b/src/sp_validation/tests/test_cli_seams.py
@@ -0,0 +1,57 @@
+"""Smoke tests for workflow CLI seams — cheap guards against signature rot.
+
+The compute these scripts drive is cluster-only, so a removed or renamed kwarg
+would only TypeError at invocation. Each test binds the exact call one seam
+makes against the current signature (``inspect.signature(...).bind(...)``) — no
+compute, no data — so the drift fails here instead of on the cluster.
+"""
+
+import importlib.util
+import inspect
+from pathlib import Path
+
+import pytest
+
+
+def _repo_root() -> Path:
+    for parent in Path(__file__).resolve().parents:
+        if (parent / "pyproject.toml").exists():
+            return parent
+    raise RuntimeError("could not locate repo root (no pyproject.toml above test)")
+
+
+def _load(path, name):
+    spec = importlib.util.spec_from_file_location(name, path)
+    module = importlib.util.module_from_spec(spec)
+    spec.loader.exec_module(module)
+    return module
+
+
+def test_run_xi_sweep_run_2pcf_call_binds():
+    """The kwargs run_xi_sweep passes to run_2pcf must bind to its signature."""
+    root = _repo_root()
+    run_2pcf_mod = _load(root / "workflow/scripts/run_2pcf.py", "run_2pcf_seam")
+    sig = inspect.signature(run_2pcf_mod.run_2pcf)
+    # Exactly the keyword set run_xi_sweep._from_cli passes.
+    sig.bind(
+        ver="V",
+        cat_config="/cfg.yaml",
+        output_dir="/out",
+        grid="integration",
+        min_sep=0.5,
+        max_sep=300.0,
+        nbins=1000,
+        npatch=1,
+    )
+    # And the removed kwarg must NOT bind (guards against a silent re-add).
+    with pytest.raises(TypeError):
+        sig.bind(
+            ver="V",
+            cat_config="/cfg.yaml",
+            output_dir="/out",
+            save_fits=True,
+            min_sep=1.0,
+            max_sep=250.0,
+            nbins=20,
+            npatch=1,
+        )
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
diff --git a/src/sp_validation/tests/test_pseudo_cl.py b/src/sp_validation/tests/test_pseudo_cl.py
index 45a2366d..6fa15d0e 100644
--- a/src/sp_validation/tests/test_pseudo_cl.py
+++ b/src/sp_validation/tests/test_pseudo_cl.py
@@ -57,7 +57,9 @@
 import pytest
 import yaml
 
+from sp_validation import sacc_io
 from sp_validation.cosmo_val import CosmologyValidation
+from sp_validation.cosmo_val.sacc_writers import BIN as SACC_BIN
 from sp_validation.rho_tau import get_params_rho_tau
 
 # These tests need the full harmonic-space stack (pymaster/NaMaster + healpy),
@@ -112,6 +114,7 @@ def _write_synthetic_config(tmp_path):
 
     shear_cfg = {
         "path": "shear.fits",
+        "redshift_path": str(nz_dir / "dndz_SP_A.txt"),
         "w_col": "w",
         "e1_col": "e1",
         "e2_col": "e2",
@@ -508,24 +511,20 @@ def test_apply_random_rotation_reproducible_with_seed(cv, cat_and_params):
 # calculate_pseudo_cl_catalog -- deterministic end-to-end catalog path
 # ===========================================================================
 def test_calculate_pseudo_cl_catalog_end_to_end(cv, tmp_path):
-    """End-to-end catalog path: FITS round-trip of ell + EE/EB/BB.
+    """End-to-end catalog path: SACC round-trip of ell + EE/EB/BB.
 
     The catalog method has no random noise debiasing, so it is reproducible to
-    the same ~2e-12 catalog-path float noise. save_pseudo_cl stores ELL/EE/EB/BB
-    (it drops the BE row); we pin the round-tripped table.
+    the same ~2e-12 catalog-path float noise; we pin the round-tripped spectra.
     """
     ver = cv._test_version
     cv._pseudo_cls = {ver: {}}
-    out_path = cv._output_path(f"pseudo_cl_cat_{ver}.fits")
+    out_path = cv._output_path(f"pseudo_cl_{ver}.sacc")
     cv.calculate_pseudo_cl_catalog(ver, out_path)
 
     assert os.path.exists(out_path)
-    d = fits.getdata(out_path)
-    # FITS gives big-endian f8; normalize for value comparison.
-    ell = np.asarray(d["ELL"], dtype=np.float64)
-    ee = np.asarray(d["EE"], dtype=np.float64)
-    eb = np.asarray(d["EB"], dtype=np.float64)
-    bb = np.asarray(d["BB"], dtype=np.float64)
+    s = sacc_io.load(out_path, allow_unblinded=True)
+    ell, ee, bb, eb, window = sacc_io.get_pseudo_cl(s, SACC_BIN)
+    assert window is not None  # the shared BandpowerWindow rides the part
 
     npt.assert_allclose(
         ell,
@@ -590,3 +589,25 @@ def test_calculate_pseudo_cl_catalog_end_to_end(cv, tmp_path):
     params = get_params_rho_tau(cv.cc[ver], survey=ver)
     _, cl_prim, _ = cv.get_pseudo_cls_catalog(catalog=cat_gal, params=params)
     npt.assert_allclose(ee, cl_prim[0], rtol=RTOL_CAT, atol=ATOL_CAT)
+
+
+def test_calculate_pseudo_cl_out_path_born_at_declared_name(cv):
+    """calculate_pseudo_cl(out_path=...) writes to the given path, never the
+    untagged native name — so the tagged and diagnostic rules stay disjoint."""
+    ver = cv._test_version
+    cv._pseudo_cls = {}
+    tagged = cv._output_path(f"pseudo_cl_{ver}_blind=A_powspace_nbins=32.sacc")
+    native = cv._output_path(f"pseudo_cl_{ver}.sacc")
+
+    cv.calculate_pseudo_cl(out_path=tagged)
+
+    assert os.path.exists(tagged)
+    assert not os.path.exists(native)  # no undeclared native basename touched
+
+
+def test_calculate_pseudo_cl_out_path_rejects_multiversion(cv):
+    """out_path targets one part; a multi-version instance must fail loudly
+    rather than write every version to the same path."""
+    cv.versions = [cv._test_version, "SecondVersion"]
+    with pytest.raises(ValueError, match="one part to one path"):
+        cv.calculate_pseudo_cl(out_path=cv._output_path("pseudo_cl_x.sacc"))
diff --git a/src/sp_validation/tests/test_sacc_writers.py b/src/sp_validation/tests/test_sacc_writers.py
new file mode 100644
index 00000000..59ea333d
--- /dev/null
+++ b/src/sp_validation/tests/test_sacc_writers.py
@@ -0,0 +1,351 @@
+"""Tests for :mod:`sp_validation.cosmo_val.sacc_writers`.
+
+Synthetic and fast: each ``*_to_sacc`` writer is exercised with in-memory
+arrays, round-tripped through ``tmp_path``, and checked against the SACC layout
+contract (data types, tags, ordering, covariance alignment). The analysis-file
+assembler is verified to produce a single ``BlockDiagonalCovariance`` covering every
+point with each per-statistic block correctly placed. One real small-nside
+NaMaster round-trip proves the pseudo-Cℓ window survives the writer path.
+"""
+
+import numpy as np
+import pytest
+
+from sp_validation import sacc_io as sio
+from sp_validation.cosmo_val import sacc_writers as sw
+
+
+def _nz(seed=0, n=40):
+    rng = np.random.default_rng(seed)
+    return np.linspace(0.01, 2.0, n), rng.uniform(0.1, 1.0, n)
+
+
+def _spd(n, seed):
+    a = np.random.default_rng(seed).normal(size=(n, n))
+    return a @ a.T + n * np.eye(n)
+
+
+def _theta(n=6):
+    return np.geomspace(1.0, 100.0, n)
+
+
+def _roundtrip(s, tmp_path, name):
+    p = tmp_path / f"{name}.sacc"
+    sio.save(s, str(p), type="mock")
+    return sio.load(str(p))
+
+
+META = {"catalogue_version": "vSYNTH", "npatch": 1}
+
+
+# --------------------------------------------------------------------------- #
+# Per-writer parts
+# --------------------------------------------------------------------------- #
+def test_xi_to_sacc_reporting(tmp_path):
+    theta = _theta()
+    xip, xim = np.arange(6) * 1e-5, np.arange(6) * 2e-5
+    s = sw.xi_to_sacc(
+        {0: _nz()}, META, theta, xip, xim, grid="reporting", theta_nom=theta * 1.01
+    )
+    s2 = _roundtrip(s, tmp_path, "xic")
+    th, p, m = sio.get_xi(s2, (0, 0), grid="reporting")
+    assert np.array_equal(th, theta)
+    assert np.array_equal(p, xip) and np.array_equal(m, xim)
+    assert s2.covariance is None  # reporting part has no cov until assembly
+
+
+def test_xi_to_sacc_integration_diagonal(tmp_path):
+    theta = np.geomspace(0.5, 300.0, 30)
+    xip, xim = np.arange(30) * 1e-5, np.arange(30) * 2e-5
+    varxip, varxim = np.arange(1, 31) * 1e-12, np.arange(1, 31) * 2e-12
+    s = sw.xi_to_sacc(
+        {0: _nz()},
+        META,
+        theta,
+        xip,
+        xim,
+        grid="integration",
+        variances=np.concatenate([varxip, varxim]),
+    )
+    assert type(s.covariance).__name__ == "DiagonalCovariance"
+    s2 = _roundtrip(s, tmp_path, "xif")
+    th, p, _ = sio.get_xi(s2, (0, 0), grid="integration")
+    assert np.array_equal(th, theta) and np.array_equal(p, xip)
+    assert np.array_equal(
+        np.diag(s2.covariance.dense), np.concatenate([varxip, varxim])
+    )
+
+
+def test_pseudo_cl_to_sacc_window_and_rows(tmp_path):
+    ell = np.array([30.0, 60.0, 90.0, 120.0])
+    nbp = len(ell)
+    # NaMaster (4, nbp): EE, EB, BE, BB.
+    cl_all = np.vstack(
+        [
+            np.arange(nbp) * 1e-9,
+            np.arange(nbp) * 2e-9,
+            np.zeros(nbp),
+            np.arange(nbp) * 3e-9,
+        ]
+    )
+
+    class _Wsp:
+        """Stand-in workspace: (n_cl_out, nbp, n_cl_in, nell) window array."""
+
+        def __init__(self, nbp, nell):
+            w = np.zeros((4, nbp, 4, nell))
+            col = np.zeros((nbp, nell))
+            for b in range(nbp):
+                col[b, b * 3 : b * 3 + 3] = 1.0
+            for out in range(4):
+                w[out, :, out, :] = col
+            self._w = w
+
+        def get_bandpower_windows(self):
+            return self._w
+
+    s = sw.pseudo_cl_to_sacc({0: _nz()}, META, ell, cl_all, _Wsp(nbp, 24))
+    s2 = _roundtrip(s, tmp_path, "cl")
+    ell_r, ee, bb, eb, window = sio.get_pseudo_cl(s2, (0, 0))
+    assert np.array_equal(ell_r, ell)
+    assert np.array_equal(ee, cl_all[0])  # EE row
+    assert np.array_equal(bb, cl_all[3])  # BB row (index 3, not 2=BE)
+    assert np.array_equal(eb, cl_all[1])  # EB row
+    assert window.weight.shape == (24, nbp)
+
+
+def test_pseudo_cl_to_sacc_real_namaster(tmp_path):
+    """A real small-nside NaMaster workspace's window survives the writer."""
+    pytest.importorskip("pymaster")
+    from sp_validation.pseudo_cl import get_pseudo_cls_map
+
+    nside = 32
+    mask = np.ones(12 * nside**2)
+    rng = np.random.default_rng(0)
+    shear = (
+        rng.normal(size=12 * nside**2) + 1j * rng.normal(size=12 * nside**2)
+    ) * 1e-2
+    ell_eff, cl_all, wsp = get_pseudo_cls_map(shear, mask, nside, "linear", ell_step=8)
+    s = sw.pseudo_cl_to_sacc({0: _nz()}, META, ell_eff, cl_all, wsp)
+    s2 = _roundtrip(s, tmp_path, "clreal")
+    ell_r, ee, bb, eb, window = sio.get_pseudo_cl(s2, (0, 0))
+    assert np.array_equal(ell_r, ell_eff)
+    assert np.array_equal(ee, cl_all[0]) and np.array_equal(bb, cl_all[3])
+    # window columns correspond to the bandpowers, one per ell_eff
+    assert window.weight.shape[1] == len(ell_eff)
+
+
+def test_cosebis_to_sacc(tmp_path):
+    En, Bn = np.arange(1, 11) * 1e-6, np.arange(1, 11) * 1e-7
+    result = {"En": En, "Bn": Bn, "cov": _spd(20, 7)}
+    s = sw.cosebis_to_sacc({0: _nz()}, META, result, (1.0, 100.0))
+    s2 = _roundtrip(s, tmp_path, "co")
+    n, E, B = sio.get_cosebis(s2, (0, 0))
+    assert np.array_equal(n, np.arange(1, 11))
+    assert np.array_equal(E, En) and np.array_equal(B, Bn)
+    assert type(s2.covariance).__name__ == "FullCovariance"
+    assert np.array_equal(s2.covariance.dense, result["cov"])
+
+
+def test_pure_eb_to_sacc(tmp_path):
+    theta = _theta()
+    eb = {key: np.arange(6) * (i + 1) * 1e-6 for i, key in enumerate(sio.PURE_KEYS)}
+    cov = _spd(6 * len(theta), 9)
+    s = sw.pure_eb_to_sacc({0: _nz()}, META, theta, eb, covariance=cov)
+    s2 = _roundtrip(s, tmp_path, "eb")
+    th, back = sio.get_pure_eb(s2, (0, 0))
+    assert np.array_equal(th, theta)
+    for key in sio.PURE_KEYS:
+        assert np.array_equal(back[key], eb[key])
+    assert np.array_equal(s2.covariance.dense, cov)
+
+
+def _rho_tau_tables(nth=6, seed=0):
+    rng = np.random.default_rng(seed)
+    theta = _theta(nth)
+    rho = {"theta": theta}
+    for k in sw.RHO_K:
+        for suffix in ("p", "m"):
+            rho[f"rho_{k}_{suffix}"] = rng.normal(size=nth) * 1e-6
+            rho[f"varrho_{k}_{suffix}"] = rng.uniform(1e-14, 1e-13, nth)
+    tau = {"theta": theta}
+    for k in sw.TAU_K:
+        for suffix in ("p", "m"):
+            tau[f"tau_{k}_{suffix}"] = rng.normal(size=nth) * 1e-6
+            tau[f"vartau_{k}_{suffix}"] = rng.uniform(1e-14, 1e-13, nth)
+    return rho, tau, theta
+
+
+def test_rho_tau_to_sacc_diagonal(tmp_path):
+    rho, tau, theta = _rho_tau_tables()
+    s = sw.rho_tau_to_sacc({0: _nz()}, META, rho, tau)
+    s2 = _roundtrip(s, tmp_path, "rt")
+    for k in sw.RHO_K:
+        th, p, m = sio.get_rho(s2, k)
+        assert np.array_equal(th, theta)
+        assert np.array_equal(p, rho[f"rho_{k}_p"])
+        assert np.array_equal(m, rho[f"rho_{k}_m"])
+    for k in sw.TAU_K:
+        th, p, m = sio.get_tau(s2, (0, 0), k)
+        assert np.array_equal(p, tau[f"tau_{k}_p"])
+    assert type(s2.covariance).__name__ == "DiagonalCovariance"
+
+
+def test_rho_tau_to_sacc_tau_theory_block(tmp_path):
+    """The (3·nbin) plus-only CovTauTh block scatters into the τ-plus rows/cols;
+    τ-minus keeps a vartau diagonal, and cross plus↔minus stays zero."""
+    rho, tau, theta = _rho_tau_tables()
+    nbin = len(theta)
+    n_plus = len(sw.TAU_K) * nbin  # τ-plus points (k-major, one component per k)
+    tau_cov_th = _spd(n_plus, 11)
+    s = sw.rho_tau_to_sacc({0: _nz()}, META, rho, tau, tau_cov_th=tau_cov_th)
+    assert type(s.covariance).__name__ == "FullCovariance"
+    tr = ("source_0", sio.PSF_TRACER)
+    tau_plus = np.concatenate(
+        [s.indices(sio.TAU_PLUS.format(k=k), tr) for k in sw.TAU_K]
+    )
+    tau_minus = np.concatenate(
+        [s.indices(sio.TAU_MINUS.format(k=k), tr) for k in sw.TAU_K]
+    )
+    s2 = _roundtrip(s, tmp_path, "rttau")
+    dense = s2.covariance.dense
+    # τ-plus sub-block equals the supplied theory covariance (scatter is correct).
+    assert np.allclose(dense[np.ix_(tau_plus, tau_plus)], tau_cov_th)
+    # τ-minus is diagonal from vartau; plus↔minus cross is zero.
+    tau_minus_var = np.concatenate([np.asarray(tau[f"vartau_{k}_m"]) for k in sw.TAU_K])
+    assert np.allclose(np.diag(dense[np.ix_(tau_minus, tau_minus)]), tau_minus_var)
+    assert np.allclose(dense[np.ix_(tau_plus, tau_minus)], 0.0)
+
+
+def test_rho_tau_to_sacc_tau_cov_shape_mismatch():
+    rho, tau, _ = _rho_tau_tables()
+    with pytest.raises(ValueError, match="tau_cov_th shape"):
+        sw.rho_tau_to_sacc({0: _nz()}, META, rho, tau, tau_cov_th=_spd(3, 1))
+
+
+# --------------------------------------------------------------------------- #
+# Analysis-file assembly
+# --------------------------------------------------------------------------- #
+def _make_parts(nz):
+    theta = _theta()
+    ell = np.array([30.0, 60.0, 90.0])
+
+    class _Wsp:
+        def get_bandpower_windows(self):
+            w = np.zeros((4, 3, 4, 20))
+            for out in range(4):
+                for b in range(3):
+                    w[out, b, out, b * 6 : b * 6 + 6] = 1.0
+            return w
+
+    xi = sw.xi_to_sacc(
+        nz, META, theta, np.arange(6) * 1e-5, np.arange(6) * 2e-5, grid="reporting"
+    )
+    xi.add_covariance(_spd(len(xi.mean), 1))
+    cl_all = np.vstack(
+        [np.arange(3) * 1e-9, np.arange(3) * 2e-9, np.zeros(3), np.arange(3) * 3e-9]
+    )
+    cl = sw.pseudo_cl_to_sacc(nz, META, ell, cl_all, _Wsp(), covariance=_spd(9, 2))
+    co = sw.cosebis_to_sacc(
+        nz,
+        META,
+        {
+            "En": np.arange(1, 6) * 1e-6,
+            "Bn": np.arange(1, 6) * 1e-7,
+            "cov": _spd(10, 3),
+        },
+        (1.0, 100.0),
+    )
+    return [xi, cl, co]
+
+
+def test_assemble_analysis_sacc_block_diagonal_covariance(tmp_path):
+    nz = {0: _nz()}
+    parts = _make_parts(nz)
+    s = sw.assemble_analysis_sacc(parts)
+    assert type(s.covariance).__name__ == "BlockDiagonalCovariance"
+    assert s.covariance.dense.shape == (len(s.mean), len(s.mean))
+    # every point covered; blocks placed and cross-blocks zero
+    tr = ("source_0", "source_0")
+    xi_idx = np.concatenate([s.indices(sio.XI_PLUS, tr), s.indices(sio.XI_MINUS, tr)])
+    cl_idx = np.concatenate(
+        [s.indices(sio.CL_EE, tr), s.indices(sio.CL_BB, tr), s.indices(sio.CL_EB, tr)]
+    )
+    co_idx = np.concatenate(
+        [s.indices(sio.COSEBI_EE, tr), s.indices(sio.COSEBI_BB, tr)]
+    )
+    assert len(xi_idx) + len(cl_idx) + len(co_idx) == len(s.mean)
+    dense = s.covariance.dense
+    assert np.array_equal(dense[np.ix_(xi_idx, xi_idx)], parts[0].covariance.dense)
+    assert np.array_equal(dense[np.ix_(cl_idx, cl_idx)], parts[1].covariance.dense)
+    assert np.array_equal(dense[np.ix_(co_idx, co_idx)], parts[2].covariance.dense)
+    assert np.array_equal(
+        dense[np.ix_(xi_idx, cl_idx)], np.zeros((len(xi_idx), len(cl_idx)))
+    )
+    # round-trips
+    s2 = _roundtrip(s, tmp_path, "analysis")
+    assert type(s2.covariance).__name__ == "BlockDiagonalCovariance"
+    assert np.allclose(s2.covariance.dense, s.covariance.dense)
+
+
+def test_assemble_analysis_sacc_requires_covariance():
+    nz = {0: _nz()}
+    parts = _make_parts(nz)
+    parts.append(
+        sw.xi_to_sacc(
+            nz,
+            META,
+            _theta(),
+            np.arange(6) * 1e-5,
+            np.arange(6) * 2e-5,
+            grid="reporting",
+        )
+    )  # no covariance
+    with pytest.raises(ValueError, match="own covariance block"):
+        sw.assemble_analysis_sacc(parts)
+
+
+def test_assemble_from_reloaded_parts(tmp_path):
+    """Parts written to disk then reloaded assemble identically (the DAG path)."""
+    nz = {0: _nz()}
+    parts = _make_parts(nz)
+    reloaded = []
+    for i, part in enumerate(parts):
+        sio.save(part, str(tmp_path / f"part{i}.sacc"), type="mock")
+        reloaded.append(sio.load(str(tmp_path / f"part{i}.sacc")))
+    s = sw.assemble_analysis_sacc(reloaded)
+    assert type(s.covariance).__name__ == "BlockDiagonalCovariance"
+    assert s.covariance.dense.shape == (len(s.mean), len(s.mean))
+
+
+def test_xi_part_carries_a_dense_covariance(tmp_path):
+    """A grid measured with patches puts its jackknife block in the part."""
+    theta = _theta()
+    cov = _spd(2 * len(theta), 41)
+    s = sw.xi_to_sacc(
+        {0: _nz()},
+        META,
+        theta,
+        np.arange(6) * 1e-5,
+        np.arange(6) * 2e-5,
+        grid="cosebis",
+        covariance=cov,
+    )
+    s2 = _roundtrip(s, tmp_path, "xi_cov")
+    assert np.allclose(s2.covariance.dense, cov)
+
+
+def test_xi_part_rejects_two_covariances():
+    """Dense block and variances are alternatives, not a merge."""
+    theta = _theta()
+    with pytest.raises(ValueError, match="not both"):
+        sw.xi_to_sacc(
+            {0: _nz()},
+            META,
+            theta,
+            np.zeros(6),
+            np.zeros(6),
+            grid="reporting",
+            covariance=_spd(12, 42),
+            variances=np.ones(12),
+        )
diff --git a/src/sp_validation/tests/test_xi_grids.py b/src/sp_validation/tests/test_xi_grids.py
new file mode 100644
index 00000000..c35a2e3f
--- /dev/null
+++ b/src/sp_validation/tests/test_xi_grids.py
@@ -0,0 +1,109 @@
+"""Tests for the ξ± grid table in ``workflow/common.py``.
+
+The table names the files the ``xi`` rule writes and the ones every consumer
+asks for, so producer and consumer agree only if the tag is built from
+canonical values. These tests pin that canonicalisation and the grid lookup.
+"""
+
+import importlib.util
+from pathlib import Path
+
+import pytest
+
+pytestmark = pytest.mark.fast
+
+
+def _load_common():
+    root = next(
+        p for p in Path(__file__).resolve().parents if (p / "pyproject.toml").exists()
+    )
+    path = root / "workflow" / "common.py"
+    spec = importlib.util.spec_from_file_location("wf_common_grids", path)
+    module = importlib.util.module_from_spec(spec)
+    spec.loader.exec_module(module)
+    return module
+
+
+common = _load_common()
+
+# A cosmo_val block as YAML delivers it: integer-valued separations stay ints.
+CONFIG = {
+    "cosmo_val": {
+        "theta_min": 1.0,
+        "theta_max": 250.0,
+        "nbins": 20,
+        "npatch": 100,
+        "integration": {"min_sep": 0.08, "max_sep": 300, "nbins": 1000},
+        "cosebis": {
+            "min_sep_int": 0.9,
+            "max_sep_int": 300,
+            "nbins_int": 1000,
+            "npatch": 100,
+        },
+    }
+}
+FIDUCIAL = {
+    "min_sep": 1.0,
+    "max_sep": 250.0,
+    "nbins": 20,
+    "npatch": 1,
+    "min_sep_int": 0.5,
+    "max_sep_int": 300,
+    "nbins_int": 1000,
+}
+
+
+def test_tag_is_built_from_canonical_values():
+    """An integer YAML separation still names the file as a float.
+
+    run_2pcf coerces separations with float() before TreeCorr writes, so a
+    `max_sep: 300` that reached the tag as "300" would have the consumer ask
+    for a path the producer never writes.
+    """
+    grids = common.xi_grids(CONFIG, FIDUCIAL)
+    assert common.grid_binning(grids["integration"]).endswith(
+        "minsep=0.08_maxsep=300.0_nbins=1000_npatch=1"
+    )
+    assert (
+        common.grid_binning(grids["cosebis"])
+        == "minsep=0.9_maxsep=300.0_nbins=1000_npatch=100"
+    )
+    # Counts stay integers, so no "nbins=1000.0" creeps into a name.
+    assert "nbins=1000_" in common.grid_binning(grids["cosebis"])
+
+
+def test_grid_lookup_round_trips_through_the_tag():
+    """Every grid's own binning resolves back to that grid.
+
+    This is the producer/consumer contract: the rule resolves a job's grid from
+    the wildcards its filename bound.
+    """
+    grids = common.xi_grids(CONFIG, FIDUCIAL)
+    for name, grid in grids.items():
+        binning = {key: grid[key] for key in common.XI_KEYS}
+        assert common.grid_of(grids, binning) == name
+        # Wildcards arrive as strings; the comparison is numeric.
+        assert common.grid_of(grids, {k: str(v) for k, v in binning.items()}) == name
+
+
+def test_covariance_mode_follows_the_patches():
+    """Patched grids get a jackknife block, unpatched ones none."""
+    grids = common.xi_grids(CONFIG, FIDUCIAL)
+    assert grids["reporting"]["cov"] == "jackknife"
+    assert grids["cosebis"]["cov"] == "jackknife"
+    assert grids["integration"]["cov"] == "none"
+
+
+def test_unnamed_binning_is_a_reporting_measurement():
+    """The paper's convergence-check binning belongs to no named grid."""
+    grids = common.xi_grids(CONFIG, FIDUCIAL)
+    stray = {"min_sep": 1.0, "max_sep": 250.0, "nbins": 10000, "npatch": 1}
+    assert common.grid_of(grids, stray) == "reporting"
+
+
+def test_workflow_without_cosmo_val_falls_back_to_fiducial():
+    """papers/bmodes carries no cosmo_val block; its grids come from FIDUCIAL."""
+    grids = common.xi_grids({}, FIDUCIAL)
+    assert grids["reporting"]["npatch"] == 1
+    assert grids["integration"]["min_sep"] == 0.5
+    assert "cosebis" not in grids
diff --git a/workflow/common.py b/workflow/common.py
index 3df3413b..ca6901d9 100644
--- a/workflow/common.py
+++ b/workflow/common.py
@@ -5,6 +5,13 @@
 import re
 from pathlib import Path
 
+# The running checkout, for rules that shell out to a script directly rather
+# than through Snakemake's `script:` directive. Anchored on this module's own
+# location, not workflow.basedir — under `module` composition basedir reflects
+# the composing paper, not the running checkout.
+REPO_ROOT = Path(os.path.realpath(__file__)).parents[1]
+WORKFLOW_SCRIPTS = os.path.join(os.path.dirname(os.path.realpath(__file__)), "scripts")
+
 # Output roots are env-overridable so a reproduction run can write into a
 # fresh tree without clobbering (or silently reusing) prior products.
 COSMO_VAL = Path(
@@ -18,6 +25,9 @@
     )
 )
 CAT_CONFIG = "/n17data/cdaley/unions/code/sp_validation/cosmo_val/cat_config.yaml"
+# "blind" is the glass-mock A/B/C realisation convention, NOT Smokescreen
+# blinding (a separate axis: the concealed=True SACC stamp). The name is baked
+# into on-disk filenames we do not own (e.g. nz_{version}_{A|B|C}.txt).
 BLINDS = ["A", "B", "C"]
 BLOCK_PAIRS = [("++", "1"), ("--", "2"), ("+-", "3")]
 
@@ -156,22 +166,122 @@ def covariance_path(
     return str(COSMO_INFERENCE / f"data/covariance/{base}/{base}{suffix}")
 
 
+def base_version(version):
+    """Strip the derived-catalogue suffixes to the base catalogue version.
+
+    The `_leak_corr` / `_ecut{N}` variants share their parent's n(z) and
+    `cov_th` survey parameters, so lookups keyed on either must strip both.
+    """
+    return re.sub(r"_ecut\d+", "", re.sub(r"_leak_corr$", "", version))
+
+
 def build_redshift_path(version, blind):
     """Construct n(z) filepath for given catalog version and blind."""
-    base_version = re.sub(r"_leak_corr$", "", version)
-    base_version = re.sub(r"_ecut\d+", "", base_version)
-    if "v1.4.11" in base_version:
-        base_version = "SP_v1.4.6"
-    version_dir = base_version.replace("SP_", "")
+    base = base_version(version)
+    if "v1.4.11" in base:
+        base = "SP_v1.4.6"
+    version_dir = base.replace("SP_", "")
+    return f"/n17data/sguerrini/UNIONS/WL/nz/{version_dir}/nz_{base}_{blind}.txt"
+
+
+# ---------------------------------------------------------------------------
+# ξ± angular grids
+# ---------------------------------------------------------------------------
+# A grid is a binning plus how its covariance is estimated: (min_sep, max_sep,
+# nbins, npatch, cov). `reporting` is the analysis grid, `integration` the fine
+# one the B-mode integrals run over, `cosebis` the fine patched grid COSEBIs
+# propagates its covariance from. cov is "jackknife" (dense, from the patches),
+# "diagonal" (TreeCorr varxip/varxim) or "none".
+XI_KEYS = (
+    "min_sep",
+    "max_sep",
+    "nbins",
+    "npatch",
+)  # the binning; cov is not in the name
+
+
+def xi_grids(config, fiducial):
+    """The named ξ± grids of a workflow, canonicalised.
+
+    Workflows carrying no cosmo_val block (e.g. papers/bmodes) fall back to
+    ``fiducial``. Values are coerced here — separations to float, counts to int
+    — so the tag the table stamps into a filename is the one the measurement
+    writes: the separations pass through float() on the way to TreeCorr, so a
+    YAML ``300`` must become ``300.0`` before it names a file, or producer and
+    consumer ask for different paths.
+    """
+    cv = config.get("cosmo_val", {})
+    grids = {
+        "reporting": (
+            {
+                "min_sep": cv["theta_min"],
+                "max_sep": cv["theta_max"],
+                "nbins": cv["nbins"],
+                "npatch": cv["npatch"],
+            }
+            if cv
+            else {k: fiducial[k] for k in XI_KEYS}
+        ),
+        "integration": dict(
+            cv.get("integration")
+            or {
+                "min_sep": fiducial["min_sep_int"],
+                "max_sep": fiducial["max_sep_int"],
+                "nbins": fiducial["nbins_int"],
+            }
+        ),
+    }
+    grids["integration"].setdefault("npatch", 1)
+    cb = cv.get("cosebis")
+    if cb:
+        grids["cosebis"] = {
+            "min_sep": cb["min_sep_int"],
+            "max_sep": cb["max_sep_int"],
+            "nbins": cb["nbins_int"],
+            "npatch": cb["npatch"],
+        }
+    for grid in grids.values():
+        for key in ("min_sep", "max_sep"):
+            grid[key] = float(grid[key])
+        for key in ("nbins", "npatch"):
+            grid[key] = int(grid[key])
+        # A jackknife estimate needs patches; at npatch=1 TreeCorr's var_method
+        # is "shot" and the diagonal is all it can offer.
+        grid.setdefault("cov", "jackknife" if grid["npatch"] > 1 else "none")
+    return grids
+
+
+def grid_binning(grid):
+    """The `minsep=..._maxsep=..._nbins=..._npatch=...` tag of one grid."""
     return (
-        f"/n17data/sguerrini/UNIONS/WL/nz/{version_dir}/nz_{base_version}_{blind}.txt"
+        f"minsep={grid['min_sep']}_maxsep={grid['max_sep']}"
+        f"_nbins={grid['nbins']}_npatch={grid['npatch']}"
     )
 
 
+def grid_of(grids, binning):
+    """Name of the grid a binning belongs to, compared numerically.
+
+    A "300" wildcard matches a 300.0 grid value. Binnings matching no named
+    grid (e.g. papers/bmodes' nbins=10000 convergence check) are reporting-style
+    measurements.
+    """
+    key = tuple(float(binning[k]) for k in XI_KEYS)
+    for name, grid in grids.items():
+        if tuple(float(grid[k]) for k in XI_KEYS) == key:
+            return name
+    return "reporting"
+
+
+def pseudo_cl_tag(config):
+    """Fiducial harmonic-binning tag stamped into pseudo-Cl filenames."""
+    fiducial = config["harmonic"]["fiducial"]
+    return f"blind={fiducial['blind']}_{fiducial['binning']}_nbins={fiducial['nbins']}"
+
+
 def get_shear_catalog(wildcards):
     """Resolve shear catalog path from config for a given version."""
-    base_version = wildcards.version.replace("_leak_corr", "")
-    cat_config = CATALOG_CONFIG[base_version]
+    cat_config = CATALOG_CONFIG[wildcards.version.replace("_leak_corr", "")]
     shear_path = cat_config["shear"]["path"]
     if shear_path.startswith("/"):
         return shear_path
diff --git a/workflow/rules/cosmo_val.smk b/workflow/rules/cosmo_val.smk
index 704b4fb0..87138bc0 100644
--- a/workflow/rules/cosmo_val.smk
+++ b/workflow/rules/cosmo_val.smk
@@ -3,29 +3,33 @@
 # The original cosmo_val/run_cosmo_val.py was one linear driver that built a
 # single in-memory `cv` (CosmologyValidation) and called ~13 cv.<method>()
 # diagnostics in sequence, linked only by lazy properties on that object. Here
-# each diagnostic is a rule, and the rules are linked by the *real* data
-# products each method writes under COSMO_VAL (= cosmo_val/output):
+# each diagnostic is a rule, and the rules are linked by the SACC parts and
+# products they write under COSMO_VAL (= cosmo_val/output):
 #
-#   rho/tau FITS  ──┬─→ rho/tau plots
-#                   ├─→ rho_tau_fits (PSF-error MCMC)
-#                   └─────────────────────────────┐
-#   additive bias ──→ xi (2pcf) ──┬─→ 2pcf plot   │
-#                                 ├─→ ratio_xi_sys_xi ←┘ (also needs xi_psf_sys)
-#                                 ├─→ pure_eb (npz) ─┐
-#                                 └─→ cosebis (npz) ─┤
-#   pseudo_cl FITS ──────────────────────────────────┼─→ summarize_bmodes
-#                                                     ┘
+#   catalogue ──→ xi (one job per grid: reporting, integration, cosebis)
+#                   │
+#                   ├─ reporting part ──┬─→ pure_eb (part, npz, figures)
+#                   ├─ integration part ┘        │
+#                   ├─ cosebis part ─────→ cosebis (part, npz, figures)
+#                   └─ reporting .txt ──→ 2pcf plot, ratio_xi_sys_xi
+#   catalogue ──→ pseudo_cl (part) ──┬─→ pseudo-Cl figures
+#   CosmoCov ──→ covariance ─────────┤
+#   rho/tau (part + FITS) ───────────┼─→ summarize_bmodes (reads the products)
+#                                    └─→ assemble_sacc ──→ {version}.sacc
 #
-# Granularity decision: methods that write durable data products
-# (calculate_rho_tau_stats, calculate_2pcf, calculate_pseudo_cl, plot_pure_eb,
-# plot_cosebis) own a compute rule keyed on those files. Methods that only
-# emit figures, or whose figure paths derive from internal handler state
-# (rho/tau plots, rho_tau_fits, objectwise leakage, 2pcf overlay), declare a
-# sentinel under COSMO_VAL/snakemake_sentinels so they stay DAG-trackable.
-# Lazy cv state that the original code never persists (c1/c2, xi_psf_sys) is
-# either materialized to a small JSON (additive bias) or recomputed in the one
-# rule that needs it (xi_psf_sys in ratio_xi_sys_xi) — recompute is cheap next
-# to the science it depends on. See workflow/scripts/cv_runner.py.
+# The B-mode rules are ingests: pure_eb, cosebis, the pseudo-Cl figures and the
+# summary all work from the parts and the covariance inputs, never from a
+# catalogue, so a blinded part keeps everything downstream blinded. The
+# analytic covariances (CosmoCov ξ±, NaMaster pseudo-Cℓ) are what assembly puts
+# in the terminal file, replacing the estimates a part was born with.
+#
+# Methods that only emit figures, or whose figure paths derive from internal
+# handler state (rho/tau plots, rho_tau_fits, objectwise leakage, 2pcf
+# overlay), declare a sentinel under COSMO_VAL/snakemake_sentinels so they stay
+# DAG-trackable. Lazy cv state the original code never persists (c1/c2,
+# xi_psf_sys) is either materialized to a small JSON (additive bias) or
+# recomputed in the one rule that needs it — recompute is cheap next to the
+# science it depends on. See workflow/scripts/cv_runner.py.
 
 CV = config["cosmo_val"]
 CV_VERSIONS = config["versions"]
@@ -52,13 +56,7 @@ def cv_xi_txt(version):
     Mirrors the out_fname f-string in cosmo_val.calculate_2pcf:
     {ver}_xi_minsep=..._maxsep=..._nbins=..._npatch=...txt
     """
-    return str(
-        COSMO_VAL
-        / (
-            f"{version}_xi_minsep={CV['theta_min']}_maxsep={CV['theta_max']}"
-            f"_nbins={CV['nbins']}_npatch={CV['npatch']}.txt"
-        )
-    )
+    return str(COSMO_VAL / f"{version}_xi_{xi_binning('reporting')}.txt")
 
 
 def cv_rho_stats(version):
@@ -73,39 +71,140 @@ def cv_tau_stats(version):
     )
 
 
-def cv_pure_eb_npz(version):
-    eb = CV["pure_eb"]
+def _pure_eb_stub(version):
+    """Shared stem of the pure-E/B diagnostic products (npz + figures)."""
+    eb = CV["integration"]
     return str(
         COSMO_VAL
         / (
             f"{version}_eb_minsep={CV['theta_min']}_maxsep={CV['theta_max']}"
-            f"_nbins={CV['nbins']}_minsepint={eb['min_sep_int']}"
-            f"_maxsepint={eb['max_sep_int']}_nbinsint={eb['nbins_int']}"
-            f"_npatch={CV['npatch']}_varmethod=jackknife_data.npz"
+            f"_nbins={CV['nbins']}_minsepint={eb['min_sep']}"
+            f"_maxsepint={eb['max_sep']}_nbinsint={eb['nbins']}"
+            f"_npatch={CV['npatch']}_varmethod=semi-analytic"
         )
     )
 
 
-def cv_cosebis_npz(version):
+def cv_pure_eb_npz(version):
+    """Pure-E/B data vectors + covariance .npz."""
+    return _pure_eb_stub(version) + "_data.npz"
+
+
+def cv_pure_eb_figures(version):
+    """The pure-E/B companion figures, by output key."""
+    stub = _pure_eb_stub(version)
+    return {
+        "figure_integration_vs_reporting": f"{stub}_integration_vs_reporting.png",
+        "figure_xis": f"{stub}_xis.png",
+        "figure_ptes": f"{stub}_ptes.png",
+        "figure_covariance": f"{stub}_covariance.png",
+    }
+
+
+def cv_xi_cov_integration(version):
+    """CosmoCov gaussian ξ± covariance on the integration grid.
+
+    The covariance model the pure-E/B Monte Carlo draws from; gaussian because
+    the draws only need the scatter a Gaussian field would give.
+    """
+    integ = CV["integration"]
+    return covariance_path(
+        version,
+        FIDUCIAL["blind"],
+        gaussian="g",
+        min_sep=integ["min_sep"],
+        max_sep=integ["max_sep"],
+        nbins=integ["nbins"],
+        mask_suffix=DEFAULT_MASK_SUFFIX,
+    )
+
+
+def _cosebis_stub(version):
+    """Shared stem of the COSEBIs diagnostic products (npz + figures).
+
+    varmethod names where the covariance came from, and these products are the
+    propagated one — which also keeps them clear of the paths plot_cosebis
+    builds for its own byproducts, so nothing overwrites a declared output.
+    """
     cb = CV["cosebis"]
     fsc = CV["fiducial_scale_cut"]
-    # Mirror calculate/plot_cosebis out_stub (cosmo_val.py): a distinct schema
-    # from pure_eb — _cosebis_ prefix, integration nbins, plus _nmodes= and
-    # _scalecut= segments. Must match save_cosebis_results exactly or
-    # verify_outputs raises and cv_summarize_bmodes deadlocks on this input.
     return str(
         COSMO_VAL
         / (
             f"{version}_cosebis_minsep={cb['min_sep_int']}"
             f"_maxsep={cb['max_sep_int']}_nbins={cb['nbins_int']}"
-            f"_npatch={cb['npatch']}_varmethod=jackknife_nmodes={cb['nmodes']}"
-            f"_scalecut={fsc[0]}-{fsc[1]}_data.npz"
+            f"_npatch={cb['npatch']}_varmethod=propagated_nmodes={cb['nmodes']}"
+            f"_scalecut={fsc[0]}-{fsc[1]}"
         )
     )
 
 
-def cv_pseudo_cl_fits(version):
-    return str(COSMO_VAL / f"pseudo_cl_{version}.fits")
+def cv_cosebis_npz(version):
+    """COSEBIs multi-cut diagnostic .npz (the PTE scan)."""
+    return _cosebis_stub(version) + "_data.npz"
+
+
+def cv_cosebis_figures(version):
+    """The COSEBIs companion figures, by output key."""
+    stub = _cosebis_stub(version)
+    return {
+        "figure_modes": f"{stub}_cosebis.png",
+        "figure_covariance": f"{stub}_covariance.png",
+        "figure_scalecut_ptes": f"{stub}_scalecut_ptes.png",
+    }
+
+
+_PSEUDO_CL_TAG = pseudo_cl_tag(config)
+
+
+def cv_pseudo_cl_analysis_sacc(version):
+    """Tagged pseudo-Cl SACC part: the harmonic block of the analysis file."""
+    return str(COSMO_VAL / f"pseudo_cl_{version}_{_PSEUDO_CL_TAG}.sacc")
+
+
+def cv_pseudo_cl_cov(version):
+    """NaMaster pseudo-Cl covariance FITS (COVAR_EE_EE/BB_BB/EB_EB extensions)."""
+    return str(COSMO_VAL / f"pseudo_cl_cov_{version}_{_PSEUDO_CL_TAG}.fits")
+
+
+def cv_xi_cov(version):
+    """CosmoCov-processed ξ± covariance, on the reporting grid's own binning."""
+    return covariance_path(
+        version,
+        FIDUCIAL["blind"],
+        gaussian="ng",
+        min_sep=CV["theta_min"],
+        max_sep=CV["theta_max"],
+        nbins=CV["nbins"],
+        mask_suffix=DEFAULT_MASK_SUFFIX,
+    )
+
+
+def cv_cosebis_sacc(version):
+    """COSEBIs SACC part, at the fiducial scale cut."""
+    return str(COSMO_VAL / f"{version}_cosebis.sacc")
+
+
+def cv_pure_eb_sacc(version):
+    """Pure-E/B SACC part."""
+    return str(COSMO_VAL / f"{version}_pure_eb.sacc")
+
+
+def cv_rho_tau_sacc(version):
+    """ρ/τ SACC part."""
+    return str(
+        COSMO_VAL / "rho_tau_stats" / f"rho_tau_{cv_basename(version, CV_FIDUCIAL)}.sacc"
+    )
+
+
+def cv_xi_sacc(version, grid):
+    """ξ± SACC part for a version on a named grid, named by that grid's binning."""
+    return str(COSMO_VAL / f"{version}_xi_{xi_binning(grid)}.sacc")
+
+
+def cv_analysis_sacc(version):
+    """Terminal assembled analysis file {version}.sacc."""
+    return str(COSMO_VAL / f"{version}.sacc")
 
 
 # Common params block shared by every cosmo_val rule: the cv constructor kwargs
@@ -274,18 +373,32 @@ rule cv_ratio_xi_sys_xi:
 # Harmonic-space pseudo-Cl
 # ---------------------------------------------------------------------------
 
-rule cv_pseudo_cl:
-    """Pseudo-Cl E/B spectra for all versions (NaMaster)."""
+def cv_pseudo_cl_figures():
+    """The pseudo-Cl figures, by output key (one per spectrum, all versions)."""
+    return {
+        f"figure_{name}": str(COSMO_VAL / f"cell_{name}.png")
+        for name in ("ee", "eb", "bb")
+    }
+
+
+rule cv_plot_pseudo_cl:
+    """The EE/EB/BB pseudo-Cl figures, from the analysis parts."""
+    input:
+        pseudo_cl=[cv_pseudo_cl_analysis_sacc(v) for v in CV_VERSIONS],
+        pseudo_cl_cov=[cv_pseudo_cl_cov(v) for v in CV_VERSIONS],
     output:
-        pseudo_cl=[cv_pseudo_cl_fits(v) for v in CV_VERSIONS],
+        **cv_pseudo_cl_figures(),
     params:
-        **cv_params(),
-    threads: 12
+        versions=CV_VERSIONS,
+        # Style is per catalogue, so the derived variants take their parent's.
+        markers=[CATALOG_CONFIG[base_version(v)]["marker"] for v in CV_VERSIONS],
+        colours=[CATALOG_CONFIG[base_version(v)]["colour"] for v in CV_VERSIONS],
+        rundir=CV_RUNDIR,
     resources:
-        mem_mb=32000,
-        runtime=180,
+        mem_mb=8000,
+        runtime=20,
     script:
-        "../scripts/cv_pseudo_cl.py"
+        "../scripts/cv_plot_pseudo_cl.py"
 
 
 # ---------------------------------------------------------------------------
@@ -293,18 +406,27 @@ rule cv_pseudo_cl:
 # ---------------------------------------------------------------------------
 
 rule cv_pure_eb:
-    """Pure E/B-mode decomposition for one version (config-space)."""
+    """Pure E/B-mode decomposition for one version, from its ξ± parts.
+
+    The modes come from the two parts; the covariance is Monte Carlo from the
+    integration-grid covariance model, so no patched estimator run is involved.
+    """
     input:
-        xi=lambda w: cv_xi_txt(w.version),
+        xi_reporting=lambda w: cv_xi_sacc(w.version, "reporting"),
+        xi_integration=lambda w: cv_xi_sacc(w.version, "integration"),
+        cov_integration=lambda w: cv_xi_cov_integration(w.version),
     output:
         npz=cv_pure_eb_npz("{version}"),
+        sacc=cv_pure_eb_sacc("{version}"),
+        **cv_pure_eb_figures("{version}"),
     params:
         version="{version}",
-        min_sep_int=CV["pure_eb"]["min_sep_int"],
-        max_sep_int=CV["pure_eb"]["max_sep_int"],
-        nbins_int=CV["pure_eb"]["nbins_int"],
+        min_sep=CV["theta_min"],
+        max_sep=CV["theta_max"],
+        nbins=CV["nbins"],
+        n_samples=CV.get("n_mc_samples", 1000),
+        cosmo_params=CV["cosmo_params"],
         fiducial_scale_cut=CV["fiducial_scale_cut"],
-        cv_init=lambda w: cv_init_params(config, version_list=[w.version]),
         rundir=CV_RUNDIR,
     threads: 24
     resources:
@@ -315,21 +437,25 @@ rule cv_pure_eb:
 
 
 rule cv_cosebis:
-    """COSEBIs E/B decomposition for one version (config-space, fine binning)."""
+    """COSEBIs E/B decomposition for one version, from its ξ± part.
+
+    Values, covariance and PTEs all come from the part: the COSEBIs covariance
+    is the part's ξ± covariance through the same kernel as the modes.
+    """
     input:
-        xi=lambda w: cv_xi_txt(w.version),
+        xi=lambda w: cv_xi_sacc(w.version, "cosebis"),
     output:
         npz=cv_cosebis_npz("{version}"),
+        sacc=cv_cosebis_sacc("{version}"),
+        **cv_cosebis_figures("{version}"),
     params:
         version="{version}",
-        min_sep_int=CV["cosebis"]["min_sep_int"],
-        max_sep_int=CV["cosebis"]["max_sep_int"],
-        nbins_int=CV["cosebis"]["nbins_int"],
-        npatch=CV["cosebis"]["npatch"],
+        min_sep=CV["cosebis"]["min_sep_int"],
+        max_sep=CV["cosebis"]["max_sep_int"],
+        nbins=CV["cosebis"]["nbins_int"],
         nmodes=CV["cosebis"]["nmodes"],
         scale_cuts=CV["cosebis"]["scale_cuts"],
         fiducial_scale_cut=CV["fiducial_scale_cut"],
-        cv_init=lambda w: cv_init_params(config, version_list=[w.version]),
         rundir=CV_RUNDIR,
     threads: 24
     resources:
@@ -345,32 +471,86 @@ rule cv_summarize_bmodes:
         pure_eb=[cv_pure_eb_npz(v) for v in CV_VERSIONS],
         cosebis=[cv_cosebis_npz(v) for v in CV_VERSIONS],
         pseudo_cl=(
-            [cv_pseudo_cl_fits(v) for v in CV_VERSIONS]
+            [cv_pseudo_cl_analysis_sacc(v) for v in CV_VERSIONS]
+            if CV.get("include_pseudo_cl", False) else []
+        ),
+        pseudo_cl_cov=(
+            [cv_pseudo_cl_cov(v) for v in CV_VERSIONS]
             if CV.get("include_pseudo_cl", False) else []
         ),
     output:
         summary_json=str(COSMO_VAL / "bmode_summary.json"),
     params:
+        versions=CV_VERSIONS,
         fiducial_scale_cut=CV["fiducial_scale_cut"],
-        pure_eb_min_sep_int=CV["pure_eb"]["min_sep_int"],
-        pure_eb_max_sep_int=CV["pure_eb"]["max_sep_int"],
-        pure_eb_nbins_int=CV["pure_eb"]["nbins_int"],
-        cosebis_min_sep_int=CV["cosebis"]["min_sep_int"],
-        cosebis_max_sep_int=CV["cosebis"]["max_sep_int"],
-        cosebis_nbins_int=CV["cosebis"]["nbins_int"],
-        cosebis_npatch=CV["cosebis"]["npatch"],
-        cosebis_nmodes=CV["cosebis"]["nmodes"],
-        cosebis_scale_cuts=CV["cosebis"]["scale_cuts"],
+        min_sep=CV["theta_min"],
+        max_sep=CV["theta_max"],
+        nbins=CV["nbins"],
         include_pseudo_cl=CV.get("include_pseudo_cl", False),
-        **cv_params(),
-    threads: 24
+        rundir=CV_RUNDIR,
     resources:
-        mem_mb=48000,
-        runtime=600,
+        mem_mb=8000,
+        runtime=20,
     script:
         "../scripts/cv_summarize_bmodes.py"
 
 
+# ---------------------------------------------------------------------------
+# Terminal analysis file: assemble the per-statistic SACC parts into {version}.sacc
+# ---------------------------------------------------------------------------
+# The terminal file carries the analysis vector only. The integration-grid ξ± is
+# deliberately not gathered: it stays a per-part intermediate. The two blocks
+# born without a covariance (ξ± reporting, pseudo-Cℓ) get theirs injected from
+# the covariance inputs below.
+
+
+def cv_assemble_inputs(version):
+    """The per-statistic SACC parts + covariance inputs assemble_sacc consumes.
+
+    Each part's filename carries enough to bind its producing rule's wildcards.
+    """
+    parts = dict(
+        xi_reporting=cv_xi_sacc(version, "reporting"),
+        xi_cov=cv_xi_cov(version),
+        cosebis=cv_cosebis_sacc(version),
+        pure_eb=cv_pure_eb_sacc(version),
+        rho_tau=cv_rho_tau_sacc(version),
+    )
+    if CV.get("include_pseudo_cl", False):
+        parts["pseudo_cl"] = cv_pseudo_cl_analysis_sacc(version)
+        parts["pseudo_cl_cov"] = cv_pseudo_cl_cov(version)
+    return parts
+
+
+rule assemble_sacc:
+    """Assemble the terminal {version}.sacc from the per-statistic SACC parts."""
+    input:
+        unpack(lambda w: cv_assemble_inputs(w.version)),
+    output:
+        sacc=cv_analysis_sacc("{version}"),
+    params:
+        version="{version}",
+        type=CV.get("type", "data"),
+        # The statistics this rule wired, so a typo'd input keyword cannot
+        # silently drop one.
+        expected=lambda w: [
+            k
+            for k in cv_assemble_inputs(w.version)
+            if k not in ("xi_cov", "pseudo_cl_cov")
+        ],
+    resources:
+        mem_mb=8000,
+        runtime=20,
+    script:
+        "../scripts/assemble_sacc.py"
+
+
+rule assemble_sacc_all:
+    """Assemble the analysis SACC file for every version."""
+    input:
+        [cv_analysis_sacc(v) for v in CV_VERSIONS],
+
+
 # ---------------------------------------------------------------------------
 # Aggregate target: the whole validation suite
 # ---------------------------------------------------------------------------
@@ -392,3 +572,6 @@ rule cosmo_val_all:
         str(COSMO_VAL / "ratio_xi_sys_xi.png"),
         # B-modes
         str(COSMO_VAL / "bmode_summary.json"),
+        list(cv_pseudo_cl_figures().values()) if CV.get("include_pseudo_cl", False) else [],
+        # Terminal analysis file: the assembled {version}.sacc per version
+        [cv_analysis_sacc(v) for v in CV_VERSIONS],
diff --git a/workflow/rules/covariance.smk b/workflow/rules/covariance.smk
index 49e40201..97feedbe 100644
--- a/workflow/rules/covariance.smk
+++ b/workflow/rules/covariance.smk
@@ -253,7 +253,7 @@ rule covariance_process:
     threads: 1
     shell:
         """
-        python /n17data/cdaley/unions/pure_eb/code/sp_validation/cosmo_inference/scripts/cosmocov_process.py {input} {params.output_stub}
+        python {REPO_ROOT}/cosmo_inference/scripts/cosmocov_process.py {input} {params.output_stub}
         """
 
 
diff --git a/workflow/rules/twopoint.smk b/workflow/rules/twopoint.smk
index 22c09db2..a3ba7e74 100644
--- a/workflow/rules/twopoint.smk
+++ b/workflow/rules/twopoint.smk
@@ -1,13 +1,34 @@
 # Two-point data-vector rules: xi, rho/tau, and pseudo-Cl products.
 
+# ---------------------------------------------------------------------------
+# ξ± angular grids
+# ---------------------------------------------------------------------------
+# The table itself lives in common.py, where it can be built from a plain config
+# dict and tested; these are the workflow's bindings to it.
+XI_GRIDS = xi_grids(config, FIDUCIAL)
+
+
+def xi_binning(grid):
+    """The `minsep=..._maxsep=..._nbins=..._npatch=...` tag of a named grid."""
+    return grid_binning(XI_GRIDS[grid])
+
+
+def xi_grid_of(wildcards):
+    """Grid label for the binning a job was requested with."""
+    return grid_of(XI_GRIDS, {key: getattr(wildcards, key) for key in XI_KEYS})
+
 
 rule xi:
+    """TreeCorr ξ±(θ) for one version on one angular grid.
+
+    One rule for every grid: outputs are named by their binning, so a request
+    binds the wildcards and `xi_grid_of` resolves the grid label from them.
+    """
     input:
         catalog=get_shear_catalog,
     output:
-        str(COSMO_VAL / "{version}_xi_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.txt"),
-        str(COSMO_VAL / "xi_plus_{version}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.fits"),
-        str(COSMO_VAL / "xi_minus_{version}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.fits"),
+        txt=str(COSMO_VAL / "{version}_xi_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.txt"),
+        sacc=str(COSMO_VAL / "{version}_xi_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.sacc"),
     threads: 24
     params:
         ver="{version}",
@@ -15,39 +36,19 @@ rule xi:
         max_sep="{max_sep}",
         nbins="{nbins}",
         npatch="{npatch}",
-        fits=False,
+        cat_config=CAT_CONFIG,
+        grid=lambda w: xi_grid_of(w),
+        cov=lambda w: XI_GRIDS[xi_grid_of(w)]["cov"],
     resources:
-        mem_mb=30000,
+        # The fine integration grid needs more memory and wall time than the
+        # ~20-bin reporting one; scale on nbins rather than splitting the rule.
+        mem_mb=lambda w: 40000 if int(w.nbins) > 100 else 30000,
         disk_mb=20000,
-        runtime=360,
+        runtime=lambda w: 600 if int(w.nbins) > 100 else 360,
     script:
         "../scripts/run_2pcf.py"
 
 
-rule xi_highres:
-    """High-resolution xi for COSEBIS integration."""
-    container: None
-    output:
-        txt=str(COSMO_VAL / f"{FIDUCIAL['version']}_xi_minsep={FIDUCIAL['min_sep_int']}_maxsep={FIDUCIAL['max_sep_int']}_nbins=10000_npatch=1.txt"),
-        xi_plus=str(COSMO_VAL / f"xi_plus_{FIDUCIAL['version']}_minsep={FIDUCIAL['min_sep_int']}_maxsep={FIDUCIAL['max_sep_int']}_nbins=10000_npatch=1.fits"),
-        xi_minus=str(COSMO_VAL / f"xi_minus_{FIDUCIAL['version']}_minsep={FIDUCIAL['min_sep_int']}_maxsep={FIDUCIAL['max_sep_int']}_nbins=10000_npatch=1.fits"),
-    resources:
-        tasks=30,
-        cpus_per_task=12,
-        nodes=6,
-        mem_mb_per_cpu=2000,
-        runtime=2880,
-        slurm_extra="'--exclude=n17,n09,n36 --partition=pscomp'",
-        mpi="/softs/openmpi/5.0.5-slurm-CentOS8/bin/mpiexec",
-    shell:
-        "{resources.mpi} -n {resources.tasks} "
-        "apptainer exec "
-        "--bind /home,/n09data,/n17data,/n23data1,/softs "
-        "--env LD_LIBRARY_PATH=/softs/openmpi/5.0.5-slurm-CentOS8/lib "
-        "/n17data/cdaley/containers/containers "
-        "python /automnt/n17data/cdaley/unions/pure_eb/code/sp_validation/workflow/scripts/run_2pcf_highres.py"
-
-
 rule run_cosmo_val:
     """Full CosmoVal diagnostic suite."""
     output:
@@ -70,6 +71,8 @@ rule rho_tau_stats:
     output:
         rho_stats=str(COSMO_VAL / "rho_tau_stats/rho_stats_{version}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.fits"),
         tau_stats=str(COSMO_VAL / "rho_tau_stats/tau_stats_{version}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.fits"),
+        # Born-as-SACC ρ/τ part, written alongside the FITS.
+        rho_tau=str(COSMO_VAL / "rho_tau_stats/rho_tau_{version}_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.sacc"),
     threads: 48
     params:
         ver="{version}",
@@ -92,9 +95,9 @@ wildcard_constraints:
 
 
 rule pseudo_cl:
-    """Generate pseudo-Cl data vector with configurable binning."""
+    """Generate pseudo-Cl data vector (born as SACC) with configurable binning."""
     output:
-        pseudo_cl=str(COSMO_VAL / "pseudo_cl_{version}_blind={blind}_{binning}_nbins={nbins}.fits"),
+        pseudo_cl=str(COSMO_VAL / "pseudo_cl_{version}_blind={blind}_{binning}_nbins={nbins}.sacc"),
     wildcard_constraints:
         blind="[ABC]",
     params:
@@ -146,7 +149,7 @@ rule pseudo_cl_all:
     """Generate pseudo-Cls for all versions."""
     input:
         expand(
-            str(COSMO_VAL / "pseudo_cl_{version}_blind=A_powspace_nbins=32.fits"),
+            str(COSMO_VAL / "pseudo_cl_{version}_blind=A_powspace_nbins=32.sacc"),
             version=PSEUDO_CL_VERSIONS,
         ),
 
@@ -164,7 +167,7 @@ rule pseudo_cl_fine_all:
     """Generate fine pseudo-Cls for COSEBIS."""
     input:
         expand(
-            str(COSMO_VAL / "pseudo_cl_{version}_blind={blind}_linear_nbins=2040.fits"),
+            str(COSMO_VAL / "pseudo_cl_{version}_blind={blind}_linear_nbins=2040.sacc"),
             version=config["versions"],
             blind=BLINDS,
         ),
diff --git a/workflow/scripts/assemble_sacc.py b/workflow/scripts/assemble_sacc.py
new file mode 100644
index 00000000..50d1d62e
--- /dev/null
+++ b/workflow/scripts/assemble_sacc.py
@@ -0,0 +1,203 @@
+"""Assemble the terminal ``{version}.sacc`` analysis file from per-statistic parts.
+
+Dual-mode: under Snakemake (``script:``) the injected ``snakemake`` object
+supplies the inputs; as a standalone CLI the same assembly runs from flags.
+
+Each part is a single-statistic SACC; they load in CANONICAL order and are
+rebuilt into one Sacc with a single ``BlockDiagonalCovariance``.
+
+Every part must carry a covariance block. ξ± reporting and pseudo-Cℓ take
+theirs from the analytic inputs — the CosmoCov ``.txt`` (``--xi-cov``) and the
+NaMaster covariance FITS (``--pseudo-cl-cov``) — which replace any estimate the
+part was born with; the pseudo-Cℓ cross-spectrum blocks (EE↔BB, …) are dropped,
+matching what the B-mode PTE reads today.
+"""
+
+import argparse
+
+import numpy as np
+
+from sp_validation import sacc_io
+from sp_validation.cosmo_val.sacc_writers import assemble_analysis_sacc
+
+# NaMaster iNKA covariance FITS: per-spectrum HDU names, in SACC insertion order.
+_CL_HDUS = ("COVAR_EE_EE", "COVAR_BB_BB", "COVAR_EB_EB")
+
+# Canonical part order — the order points are inserted in, which must match the
+# covariance block order. Missing parts are simply skipped.
+CANONICAL = ("xi_reporting", "pseudo_cl", "cosebis", "pure_eb", "rho_tau")
+
+
+def _pseudo_cl_cov_block(cov_fits):
+    """Block-diagonal ``[EE; BB; EB]`` from the NaMaster iNKA covariance FITS."""
+    from astropy.io import fits
+
+    with fits.open(cov_fits) as hdul:
+        missing = [name for name in _CL_HDUS if name not in {h.name for h in hdul}]
+        if missing:
+            raise ValueError(f"{cov_fits} lacks the pseudo-Cℓ cov HDUs {missing}")
+        blocks = [np.asarray(hdul[name].data, float) for name in _CL_HDUS]
+    n = blocks[0].shape[0]
+    full = np.zeros((3 * n, 3 * n))
+    for i, block in enumerate(blocks):
+        full[i * n : (i + 1) * n, i * n : (i + 1) * n] = block
+    return full
+
+
+# The statistics whose analysis covariance is external, and the input each one
+# takes it from. A part of one of these types may be born with an estimate of
+# its own — the ξ± reporting part carries the jackknife it was measured with —
+# but the analysis file takes the external one, always.
+_INJECTED = {"xi_reporting": "xi_cov", "pseudo_cl": "pseudo_cl_cov"}
+
+
+def _attach_cov(part, name, xi_cov, pseudo_cl_cov):
+    """Give ``part`` (mutated in place) the covariance the analysis file uses.
+
+    For the two statistics with an external covariance the supplied block
+    replaces whatever the part was born with, loudly; every other part keeps
+    its own. Raises if the block a part needs was not supplied.
+    """
+    if name not in _INJECTED:
+        if part.covariance is None:
+            raise ValueError(
+                f"the {name!r} part carries no covariance and none is injected "
+                "for it; its writer must attach one"
+            )
+        return part
+
+    supplied = xi_cov if name == "xi_reporting" else pseudo_cl_cov
+    if supplied is None:
+        raise ValueError(
+            f"the {name!r} part takes its analysis covariance from "
+            f"--{_INJECTED[name].replace('_', '-')}, which was not supplied"
+        )
+    block = (
+        np.loadtxt(supplied)
+        if name == "xi_reporting"
+        else _pseudo_cl_cov_block(supplied)
+    )
+    if part.covariance is not None:
+        print(
+            f"{name}: replacing the part's own covariance with {supplied} "
+            "(the analysis covariance)"
+        )
+    part.add_covariance(block, overwrite=True)
+    return part
+
+
+def assemble_sacc(
+    version,
+    part_paths,
+    out_path,
+    *,
+    expected=None,
+    xi_cov=None,
+    pseudo_cl_cov=None,
+    allow_unblinded=False,
+):
+    """Assemble ``{version}.sacc`` from the per-statistic ``part_paths`` mapping.
+
+    Parameters
+    ----------
+    version : str
+        Catalogue version, for error messages.
+    part_paths : dict
+        ``{statistic: path}`` with statistic in :data:`CANONICAL`. Only present
+        statistics are assembled; order is forced to canonical.
+    expected : sequence of str, optional
+        Statistics that must be present, from the caller's config toggles. A
+        typo'd input keyword would otherwise silently drop a statistic.
+    xi_cov, pseudo_cl_cov
+        Covariance sourcing — see the module docstring.
+    allow_unblinded : bool, optional
+        Passed to :func:`sacc_io.load` for every part; ``True`` only for mocks.
+    """
+    if expected is not None:
+        unknown = [name for name in expected if name not in CANONICAL]
+        if unknown:
+            raise ValueError(
+                f"expected parts {unknown} are not assemblable statistics; "
+                f"valid names are {CANONICAL}"
+            )
+        missing = [name for name in expected if not part_paths.get(name)]
+        if missing:
+            raise ValueError(
+                f"expected parts {missing} missing from part_paths for {version} "
+                f"(got {sorted(part_paths)}); a required statistic would be "
+                "silently dropped from the terminal analysis file"
+            )
+    parts = []
+    for name in CANONICAL:
+        path = part_paths.get(name)
+        if path is None:
+            continue
+        part = sacc_io.load(path, allow_unblinded=allow_unblinded)
+        parts.append(_attach_cov(part, name, xi_cov, pseudo_cl_cov))
+    if not parts:
+        raise ValueError(f"no parts found for {version}: {part_paths}")
+    s = assemble_analysis_sacc(parts)
+    sacc_io.save(s, out_path, type=s.metadata["type"])
+    print(f"Assembled {len(parts)} parts -> {out_path}")
+    return s
+
+
+def _from_snakemake(smk):
+    p = smk.params
+    inp = smk.input
+    part_paths = {
+        name: getattr(inp, name)
+        for name in CANONICAL
+        if hasattr(inp, name) and getattr(inp, name)
+    }
+    assemble_sacc(
+        version=p["version"],
+        part_paths=part_paths,
+        out_path=str(smk.output[0]),
+        expected=list(p["expected"]),
+        xi_cov=getattr(inp, "xi_cov", None),
+        pseudo_cl_cov=getattr(inp, "pseudo_cl_cov", None),
+        allow_unblinded=(p.get("type", "data") == "mock"),
+    )
+
+
+def _from_cli(argv=None):
+    ap = argparse.ArgumentParser(
+        description="Assemble the terminal {version}.sacc from per-statistic parts."
+    )
+    ap.add_argument("--version", required=True, help="Catalogue version")
+    ap.add_argument("--out", required=True, help="Output {version}.sacc path")
+    ap.add_argument(
+        "--type",
+        choices=("data", "mock"),
+        default="data",
+        help="Run type. 'mock' reads parts freely; 'data' fails closed on "
+        "unblinded parts (only concealed/blinded parts load).",
+    )
+    for name in CANONICAL:
+        ap.add_argument(
+            f"--{name.replace('_', '-')}", default=None, help=f"{name} part"
+        )
+    ap.add_argument("--xi-cov", default=None, help="CosmoCov ξ covariance .txt")
+    ap.add_argument(
+        "--pseudo-cl-cov", default=None, help="NaMaster pseudo-Cℓ covariance FITS"
+    )
+    a = ap.parse_args(argv)
+    part_paths = {name: getattr(a, name) for name in CANONICAL if getattr(a, name)}
+    assemble_sacc(
+        version=a.version,
+        part_paths=part_paths,
+        out_path=a.out,
+        xi_cov=a.xi_cov,
+        pseudo_cl_cov=a.pseudo_cl_cov,
+        allow_unblinded=(a.type == "mock"),
+    )
+
+
+if __name__ == "__main__":
+    try:
+        snakemake  # noqa: F821 — injected by Snakemake's script: directive
+    except NameError:
+        _from_cli()
+    else:
+        _from_snakemake(snakemake)  # noqa: F821
diff --git a/workflow/scripts/cv_cosebis.py b/workflow/scripts/cv_cosebis.py
index 182cda99..170341a3 100644
--- a/workflow/scripts/cv_cosebis.py
+++ b/workflow/scripts/cv_cosebis.py
@@ -1,26 +1,71 @@
 """Rule cv_cosebis: COSEBIs E/B decomposition for one version.
 
-Compute + plot rule (per version). plot_cosebis calls calculate_cosebis over a
-fine integration binning (the 2000-bin TreeCorr is the dominant cost) and
-evaluates the configured scale cuts. Writes the {version}_eb_..._data.npz
-COSEBIs data product (declared output) plus figures, and the per-version
-COSEBIs PTE that cv_summarize_bmodes collects.
+A consumer of the ξ± part alone — values, covariance, PTEs and figures all
+derive from it, so nothing here touches a catalogue. The part's ξ± covariance
+goes through the same linear kernel as the modes to give the COSEBIs
+covariance; its ``npatch`` metadata sets the Hartlap debiasing.
 """
 
-from cv_runner import _unbuffer_streams, make_cv, verify_outputs
+from cv_runner import _unbuffer_streams, verify_outputs
 from snakemake.script import snakemake
 
+from sp_validation import sacc_io
+from sp_validation.b_modes import (
+    cosebis_scan_from_xi,
+    find_conservative_scale_cut_key,
+    log_bin_edges,
+    plot_cosebis_covariance_matrix,
+    plot_cosebis_modes,
+    plot_cosebis_scale_cut_heatmap,
+    save_cosebis_results,
+)
+from sp_validation.cosmo_val.sacc_writers import cosebis_to_sacc
+
 _unbuffer_streams()
-cv = make_cv(snakemake)
 p = snakemake.params
-cv.plot_cosebis(
-    version=p["version"],
-    min_sep_int=p["min_sep_int"],
-    max_sep_int=p["max_sep_int"],
-    nbins_int=p["nbins_int"],
-    npatch=p["npatch"],
+version = p["version"]
+fiducial_scale_cut = tuple(p["fiducial_scale_cut"])
+
+part = sacc_io.load(snakemake.input["xi"])
+theta, xip, xim = sacc_io.get_xi(part, (0, 0), grid="cosebis")
+edges = log_bin_edges(p["min_sep"], p["max_sep"], p["nbins"])
+
+results = cosebis_scan_from_xi(
+    theta,
+    xip,
+    xim,
+    part.covariance.dense,
+    *edges,
     nmodes=p["nmodes"],
     scale_cuts=[tuple(sc) for sc in p["scale_cuts"]],
-    fiducial_scale_cut=tuple(p["fiducial_scale_cut"]),
+    npatch=part.metadata["npatch"],
+)
+
+fiducial_key = find_conservative_scale_cut_key(results, fiducial_scale_cut)
+fiducial = results[fiducial_key]
+
+plot_cosebis_modes(
+    fiducial,
+    version,
+    snakemake.output["figure_modes"],
+    fiducial_scale_cut=fiducial_scale_cut,
+)
+plot_cosebis_covariance_matrix(
+    fiducial, version, "jackknife", snakemake.output["figure_covariance"]
+)
+plot_cosebis_scale_cut_heatmap(
+    results,
+    edges,
+    version,
+    snakemake.output["figure_scalecut_ptes"],
+    fiducial_scale_cut=fiducial_scale_cut,
 )
+
+save_cosebis_results(results, snakemake.output["npz"], fiducial_scale_cut)
+
+# The part inherits the ξ± part's provenance; `type` is re-stamped on save.
+metadata = {k: v for k, v in part.metadata.items() if k != "type"}
+s = cosebis_to_sacc({0: sacc_io.get_nz(part, 0)}, metadata, fiducial, fiducial_key)
+sacc_io.save(s, snakemake.output["sacc"], type="data")
+
 verify_outputs(snakemake)
diff --git a/workflow/scripts/cv_plot_pseudo_cl.py b/workflow/scripts/cv_plot_pseudo_cl.py
new file mode 100644
index 00000000..7d54ef79
--- /dev/null
+++ b/workflow/scripts/cv_plot_pseudo_cl.py
@@ -0,0 +1,40 @@
+"""Rule cv_plot_pseudo_cl: the EE/EB/BB pseudo-Cl figures.
+
+Plot-only, and an ingest like the other B-mode rules: the spectra come from the
+analysis pseudo-Cl parts and their NaMaster covariances, the same pair the
+summary and the terminal file are built from, so the figures cannot show
+something the data products do not.
+"""
+
+import numpy as np
+from astropy.io import fits
+from cv_runner import _unbuffer_streams, verify_outputs
+from snakemake.script import snakemake
+
+from sp_validation import sacc_io
+from sp_validation.cosmo_val.pseudo_cl import plot_pseudo_cl_spectrum
+
+_unbuffer_streams()
+p = snakemake.params
+
+spectra = {}
+for i, version in enumerate(p["versions"]):
+    part = sacc_io.load(snakemake.input["pseudo_cl"][i])
+    ell, ee, bb, eb, _window = sacc_io.get_pseudo_cl(part, (0, 0))
+    with fits.open(snakemake.input["pseudo_cl_cov"][i]) as hdul:
+        covs = {
+            name: np.asarray(hdul[f"COVAR_{name}_{name}"].data, float)
+            for name in ("EE", "EB", "BB")
+        }
+    for name, cl in (("EE", ee), ("EB", eb), ("BB", bb)):
+        spectra.setdefault(name, {})[version] = {
+            "ell": ell,
+            "cl": cl,
+            "cov": covs[name],
+            "style": {"marker": p["markers"][i], "colour": p["colours"][i]},
+        }
+
+for name, datasets in spectra.items():
+    plot_pseudo_cl_spectrum(datasets, name, snakemake.output[f"figure_{name.lower()}"])
+
+verify_outputs(snakemake)
diff --git a/workflow/scripts/cv_pseudo_cl.py b/workflow/scripts/cv_pseudo_cl.py
deleted file mode 100644
index cf04e8e8..00000000
--- a/workflow/scripts/cv_pseudo_cl.py
+++ /dev/null
@@ -1,14 +0,0 @@
-"""Rule cv_pseudo_cl: harmonic-space pseudo-Cl B-mode spectra.
-
-plot_pseudo_cl triggers calculate_pseudo_cl, which writes pseudo_cl_{version}.fits
-for every version (the BB spectrum cv_summarize_bmodes reads) and the cell_ee.png
-figure. The per-version FITS files are the declared outputs.
-"""
-
-from cv_runner import _unbuffer_streams, make_cv, verify_outputs
-from snakemake.script import snakemake
-
-_unbuffer_streams()
-cv = make_cv(snakemake)
-cv.plot_pseudo_cl()
-verify_outputs(snakemake)
diff --git a/workflow/scripts/cv_pure_eb.py b/workflow/scripts/cv_pure_eb.py
index d15a763f..d2acdbfe 100644
--- a/workflow/scripts/cv_pure_eb.py
+++ b/workflow/scripts/cv_pure_eb.py
@@ -1,25 +1,114 @@
 """Rule cv_pure_eb: pure E/B-mode decomposition for one version.
 
-Compute + plot rule (per version). plot_pure_eb calls calculate_pure_eb, which
-runs two TreeCorr correlations (reporting + integration binning); the reporting
-binning reuses the cv_2pcf data vector via calculate_2pcf's skip-if-exists
-path. Writes the {version}_eb_..._data.npz data product (declared output) plus
-companion figures, and the per-version E/B PTEs that cv_summarize_bmodes
-collects.
+A consumer of the two ξ± parts plus one covariance file — nothing here touches
+a catalogue. The modes come from the reporting and integration parts through
+the pipeline kernel; the covariance is Monte Carlo through that same kernel,
+drawn from the CosmoCov integration-grid ξ± covariance around a theory mean, so
+it depends on the covariance model and the grids rather than on the measured
+vector. A jackknife of the transformed modes would need per-patch realisations,
+which are never persisted.
 """
 
-from cv_runner import _unbuffer_streams, make_cv, verify_outputs
+import numpy as np
+from cs_util.cosmo import get_cosmo
+from cv_runner import _unbuffer_streams, verify_outputs
 from snakemake.script import snakemake
 
+from sp_validation import sacc_io
+from sp_validation.b_modes import (
+    calculate_eb_statistics,
+    log_bin_edges,
+    plot_eb_covariance_matrix,
+    plot_integration_vs_reporting,
+    plot_pte_2d_heatmaps,
+    plot_pure_eb_correlations,
+    pure_eb_covariance_mc,
+    pure_eb_from_xi,
+    save_pure_eb_results,
+)
+from sp_validation.cosmo_val.sacc_writers import pure_eb_to_sacc
+
 _unbuffer_streams()
-cv = make_cv(snakemake)
 p = snakemake.params
-cv.plot_pure_eb(
-    versions=[p["version"]],
-    min_sep_int=p["min_sep_int"],
-    max_sep_int=p["max_sep_int"],
-    nbins_int=p["nbins_int"],
-    fiducial_xip_scale_cut=tuple(p["fiducial_scale_cut"]),
-    fiducial_xim_scale_cut=tuple(p["fiducial_scale_cut"]),
+version = p["version"]
+fiducial_scale_cut = tuple(p["fiducial_scale_cut"])
+
+reporting = sacc_io.load(snakemake.input["xi_reporting"])
+integration = sacc_io.load(snakemake.input["xi_integration"])
+theta, xip, xim = sacc_io.get_xi(reporting, (0, 0), grid="reporting")
+theta_int, xip_int, xim_int = sacc_io.get_xi(integration, (0, 0), grid="integration")
+left_edges, right_edges = log_bin_edges(p["min_sep"], p["max_sep"], p["nbins"])
+
+# The reporting grid must sit strictly inside the integration grid: a reporting
+# point on the boundary has no interior support and comes back NaN.
+modes = pure_eb_from_xi(
+    theta, xip, xim, theta_int, xip_int, xim_int, left_edges[0], right_edges[-1]
+)
+
+z, nz = sacc_io.get_nz(reporting, 0)
+cov, eb_samples = pure_eb_covariance_mc(
+    theta=theta,
+    left_edges=left_edges,
+    right_edges=right_edges,
+    theta_int=theta_int,
+    cov_int=np.loadtxt(snakemake.input["cov_integration"]),
+    z=z,
+    nz=nz,
+    cosmo=get_cosmo(**p["cosmo_params"]),
+    n_samples=p["n_samples"],
+)
+
+variances = reporting.covariance.dense.diagonal()
+results = {
+    "theta": theta,
+    "left_edges": left_edges,
+    "right_edges": right_edges,
+    "xip": xip,
+    "xim": xim,
+    "var_xip": variances[: len(theta)],
+    "var_xim": variances[len(theta) :],
+    "theta_int": theta_int,
+    "xip_int": xip_int,
+    "xim_int": xim_int,
+    "n_eff": p["n_samples"],
+    "cov": cov,
+    "eb_samples": eb_samples,
+    **modes,
+}
+results = calculate_eb_statistics(results)
+
+plot_integration_vs_reporting(
+    results, snakemake.output["figure_integration_vs_reporting"], version
+)
+plot_pure_eb_correlations(
+    results,
+    snakemake.output["figure_xis"],
+    version,
+    fiducial_xip_scale_cut=fiducial_scale_cut,
+    fiducial_xim_scale_cut=fiducial_scale_cut,
+)
+plot_pte_2d_heatmaps(
+    results,
+    version,
+    snakemake.output["figure_ptes"],
+    fiducial_xip_scale_cut=fiducial_scale_cut,
+    fiducial_xim_scale_cut=fiducial_scale_cut,
+)
+plot_eb_covariance_matrix(
+    cov, "semi-analytic", snakemake.output["figure_covariance"], version
 )
+
+save_pure_eb_results(results, snakemake.output["npz"])
+
+# The part inherits the ξ± part's provenance; `type` is re-stamped on save.
+metadata = {k: v for k, v in reporting.metadata.items() if k != "type"}
+s = pure_eb_to_sacc(
+    {0: (z, nz)},
+    metadata,
+    theta,
+    {key: results[key] for key in sacc_io.PURE_KEYS},
+    covariance=cov,
+)
+sacc_io.save(s, snakemake.output["sacc"], type="data")
+
 verify_outputs(snakemake)
diff --git a/workflow/scripts/cv_summarize_bmodes.py b/workflow/scripts/cv_summarize_bmodes.py
index 90df5999..1200ed4d 100644
--- a/workflow/scripts/cv_summarize_bmodes.py
+++ b/workflow/scripts/cv_summarize_bmodes.py
@@ -1,56 +1,64 @@
 """Rule cv_summarize_bmodes: collect B-mode PTEs across all statistics.
 
-The terminal diagnostic. summarize_bmodes reads the in-memory
-_pure_eb_results / _cosebis_results / _pseudo_cls dicts, which are populated by
-plot_pure_eb / plot_cosebis / plot_pseudo_cl. The per-version E/B and COSEBIs
-npz products and the pseudo-Cl FITS are declared as inputs (so the DAG forces
-those rules first), but the summary still needs the live result objects (it
-reads each version's TreeCorr `gg`, which the npz cannot hold). So this rule
-re-runs the three B-mode methods in-process: they reload the existing 2pcf /
-data-vector files via their skip-if-exists paths and recompute only the cheap
-PTE statistics, exactly as the original linear driver did on its shared cv.
-
-Writes the summary table to bmode_summary.txt (declared output) — the original
-driver only printed it.
+The terminal diagnostic, and a reader of what the three B-mode rules already
+wrote: the pure-E/B PTE matrices and the COSEBIs B-mode PTE from their .npz
+products, and the pseudo-Cℓ BB spectrum from its SACC part against the NaMaster
+covariance. Nothing is recomputed and no catalogue is touched, so the summary
+cannot disagree with the products it summarises.
 """
 
 import json
 
-from cv_runner import _unbuffer_streams, make_cv, verify_outputs
+import numpy as np
+from cv_runner import _unbuffer_streams, verify_outputs
 from snakemake.script import snakemake
 
+from sp_validation import sacc_io
+from sp_validation.b_modes import _get_pte_from_scale_cut, log_bin_edges
+from sp_validation.cosmo_val.core import print_bmode_summary
+from sp_validation.statistics import chi2_and_pte
+
 _unbuffer_streams()
-cv = make_cv(snakemake)
 p = snakemake.params
 fiducial_scale_cut = tuple(p["fiducial_scale_cut"])
+edges = log_bin_edges(p["min_sep"], p["max_sep"], p["nbins"])
+
+summary = {}
+cov_methods = set()
+
+for i, version in enumerate(p["versions"]):
+    row = {}
+
+    pure_eb = np.load(snakemake.input["pure_eb"][i])
+    for stat in ("xip_B", "xim_B", "combined"):
+        try:
+            row[stat] = _get_pte_from_scale_cut(
+                pure_eb[f"pte_matrices_{stat}"], edges, fiducial_scale_cut
+            )
+        except (KeyError, RuntimeError):
+            pass
+    cov_methods.add(f"pure-E/B: semi-analytic ({int(pure_eb['n_eff'])} draws)")
+
+    # The COSEBIs .npz is written at the fiducial cut, so its PTE is the one
+    # this table wants.
+    cosebis = np.load(snakemake.input["cosebis"][i])
+    row["COSEBIS"] = float(cosebis["pte_B"])
+    cov_methods.add("COSEBIs: propagated from the ξ± covariance")
+
+    if p["include_pseudo_cl"]:
+        from astropy.io import fits
+
+        part = sacc_io.load(snakemake.input["pseudo_cl"][i])
+        _ell, _ee, bb, _eb, _window = sacc_io.get_pseudo_cl(part, (0, 0))
+        with fits.open(snakemake.input["pseudo_cl_cov"][i]) as hdul:
+            cov_bb = np.asarray(hdul["COVAR_BB_BB"].data, float)
+        _chi2, _red, row["C_l_BB"] = chi2_and_pte(bb, cov_bb)
+        cov_methods.add("pseudo-Cℓ: Gaussian (NaMaster)")
+
+    summary[version] = row
+
+print_bmode_summary(summary, fiducial_scale_cut, cov_methods)
 
-# Repopulate the in-memory B-mode result dicts from existing data products.
-cv.plot_pure_eb(
-    min_sep_int=p["pure_eb_min_sep_int"],
-    max_sep_int=p["pure_eb_max_sep_int"],
-    nbins_int=p["pure_eb_nbins_int"],
-    fiducial_xip_scale_cut=fiducial_scale_cut,
-    fiducial_xim_scale_cut=fiducial_scale_cut,
-)
-for version in cv.versions:
-    cv.plot_cosebis(
-        version=version,
-        min_sep_int=p["cosebis_min_sep_int"],
-        max_sep_int=p["cosebis_max_sep_int"],
-        nbins_int=p["cosebis_nbins_int"],
-        npatch=p["cosebis_npatch"],
-        nmodes=p["cosebis_nmodes"],
-        scale_cuts=[tuple(sc) for sc in p["cosebis_scale_cuts"]],
-        fiducial_scale_cut=fiducial_scale_cut,
-    )
-if p.get("include_pseudo_cl", False):
-    cv.plot_pseudo_cl()
-
-summary = cv.summarize_bmodes(fiducial_scale_cut=fiducial_scale_cut)
-
-# summarize_bmodes prints its table and returns {version: {stat: pte}}. Persist
-# the returned dict (the table itself is reproducible from it) so downstream
-# tooling and the all-rule have a real, machine-readable artifact to depend on.
 with open(snakemake.output["summary_json"], "w") as f:
     json.dump(summary, f, indent=2, default=str)
 
diff --git a/workflow/scripts/generate_cosmocov_ini.py b/workflow/scripts/generate_cosmocov_ini.py
new file mode 100644
index 00000000..6996c04f
--- /dev/null
+++ b/workflow/scripts/generate_cosmocov_ini.py
@@ -0,0 +1,149 @@
+"""Generate a CosmoCov ``.ini`` for one (version, blind, grid, flavour, mask).
+
+Cosmology comes from the frozen ``planck18.json`` snapshot, survey parameters
+(area, n_eff, sigma_e) from the catalog config's per-version ``cov_th``, and
+n(z) from ``workflow/common.build_redshift_path``. The footprint mask power
+spectrum is passed explicitly (empty string for the unmasked variant).
+
+    python generate_cosmocov_ini.py \
+        --version SP_v1.4.6.3_leak_corr --blind A \
+        --planck18-json <cosmology_snapshot>/planck18.json \
+        --cat-config <cosmo_val/cat_config.yaml> \
+        --min-sep 0.5 --max-sep 300.0 --nbins 1000 --gaussian g \
+        --mask-cls <mask_cls_footprint_nside_4096_norm.txt> \
+        --out-ini <path/to/covariance.ini>
+"""
+
+import argparse
+import importlib.util
+import json
+import os
+import sys
+
+import yaml
+
+
+def _load_workflow_common():
+    """Load ``workflow/common.py`` (this script also runs outside Snakemake)."""
+    path = os.path.join(
+        os.path.dirname(os.path.dirname(os.path.realpath(__file__))), "common.py"
+    )
+    spec = importlib.util.spec_from_file_location("workflow_common", path)
+    module = importlib.util.module_from_spec(spec)
+    sys.modules[spec.name] = module
+    spec.loader.exec_module(module)
+    return module
+
+
+common = _load_workflow_common()
+
+
+INI_TEMPLATE = """\
+#
+# Cosmological parameters
+#
+Omega_m : {Omega_m}
+Omega_v : {Omega_v}
+sigma_8 : {sigma_8}
+n_spec : {n_s}
+w0 : -1
+wa : 0
+omb : {Omega_b}
+h0 : {h}
+
+
+# Survey and galaxy parameters
+#
+# area in degrees
+# n_gal,lens_n_gal in gals/arcmin^2
+
+area : {area}
+sourcephotoz : multihisto
+lensphotoz : multihisto
+source_tomobins : 1
+lens_tomobins : 1
+sigma_e : {sigma_e}
+source_n_gal : {n_e}
+lens_n_gal : {n_e}
+
+
+shear_REDSHIFT_FILE : {nz}
+clustering_REDSHIFT_FILE : {nz}
+c_footprint_file : {mask}
+
+
+# IA parameters
+IA : 1
+A_ia : 0.0
+eta_ia : 0.0
+
+
+# Covariance parameters
+#
+# tmin,tmax in arcminutes
+tmin : {min_sep}
+tmax : {max_sep}
+ntheta : {nbins}
+ng : {ng}
+cng : {ng}
+
+
+outdir : ./
+filename : cov_tmp
+ss : true
+ls : false
+ll : false
+"""
+
+
+def main(argv=None):
+    ap = argparse.ArgumentParser(description=__doc__.split("\n")[0])
+    ap.add_argument("--version", required=True)
+    ap.add_argument("--blind", default="A")
+    ap.add_argument("--planck18-json", required=True)
+    ap.add_argument("--cat-config", required=True)
+    ap.add_argument("--min-sep", required=True, help="tmin arcmin (string, e.g. 0.5)")
+    ap.add_argument("--max-sep", required=True, help="tmax arcmin (string, e.g. 300.0)")
+    ap.add_argument("--nbins", required=True, help="ntheta (string, e.g. 1000)")
+    ap.add_argument("--gaussian", required=True, choices=["g", "ng"])
+    ap.add_argument(
+        "--mask-cls", default="", help="footprint mask Cl path ('' = unmasked)"
+    )
+    ap.add_argument("--out-ini", required=True)
+    a = ap.parse_args(argv)
+
+    with open(a.planck18_json) as f:
+        cosmo = json.load(f)
+    with open(a.cat_config) as f:
+        cat_config = yaml.safe_load(f)
+
+    cov_th = cat_config[common.base_version(a.version)]["cov_th"]
+
+    ng_value = "1" if a.gaussian == "ng" else "0"
+
+    ini = INI_TEMPLATE.format(
+        Omega_m=cosmo["Omega_m"],
+        Omega_v=cosmo["Omega_v"],
+        sigma_8=cosmo["sigma_8"],
+        n_s=cosmo["n_s"],
+        Omega_b=cosmo["Omega_b"],
+        h=cosmo["h"],
+        area=cov_th["A"],
+        sigma_e=cov_th["sigma_e"],
+        n_e=cov_th["n_e"],
+        nz=common.build_redshift_path(a.version, a.blind),
+        mask=a.mask_cls,
+        min_sep=a.min_sep,
+        max_sep=a.max_sep,
+        nbins=a.nbins,
+        ng=ng_value,
+    )
+
+    os.makedirs(os.path.dirname(os.path.abspath(a.out_ini)), exist_ok=True)
+    with open(a.out_ini, "w") as f:
+        f.write(ini)
+    print(f"Wrote {a.out_ini}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/workflow/scripts/generate_pseudo_cl.py b/workflow/scripts/generate_pseudo_cl.py
index a5c19a56..e25f6ed9 100644
--- a/workflow/scripts/generate_pseudo_cl.py
+++ b/workflow/scripts/generate_pseudo_cl.py
@@ -1,15 +1,11 @@
 """Generate pseudo-Cls (data vector only, no covariance).
 
 Dual-mode. Under Snakemake (``script:`` directive) the injected ``snakemake``
-object supplies the parameters and the native product is renamed to the tagged
-output filename the rule declares; as a standalone CLI (argparse) the same
-compute runs from explicit flags and the primitive's native
-``pseudo_cl_{ver}.fits`` is left in place under ``--out`` (no rename — each
-lc/ASTRA recipe gets its own output directory, so the untagged native name is
-unambiguous and the primitives' skip-if-exists never collides across nbins
-runs). The CLI form is what the lightcone/ASTRA recipe calls, so the
-measurement is driven directly (no nested Snakemake) with lc handling
-orchestration:
+object supplies the parameters; as a standalone CLI (argparse) the same compute
+runs from explicit flags. Either way the part is born at its final path, as
+SACC (EE/BB/EB with a shared bandpower window). The CLI form is what the
+lightcone/ASTRA recipe calls, driving the measurement directly (no nested
+Snakemake) with lc handling orchestration:
 
     python generate_pseudo_cl.py \
         --ver SP_v1.4.6.3_leak_corr \
@@ -28,14 +24,13 @@
 import json
 import os
 
-from astropy.io import fits
-
+from sp_validation import sacc_io
 from sp_validation.cosmo_val import CosmologyValidation
 
 
 def generate_pseudo_cl(
     version: str,
-    output_dir: str,
+    out_path: str,
     cat_config: str,
     nside: int = 1024,
     npatch: int = 1,
@@ -45,16 +40,15 @@ def generate_pseudo_cl(
     nbins: int = None,
     power: float = 0.5,
 ):
-    """Generate a pseudo-Cl data vector into ``output_dir``.
+    """Generate a pseudo-Cl data vector, born as a SACC part at ``out_path``.
 
     Parameters
     ----------
     version : str
         Catalog version (e.g., "SP_v1.4.6_leak_corr")
-    output_dir : str
-        Directory the pseudo-Cl FITS file is written into. The primitive writes
-        its native ``pseudo_cl_{version}.fits`` here; callers that need a tagged
-        filename rename it themselves (see ``_from_snakemake``).
+    out_path : str
+        Exact destination the SACC part is born at — its final (possibly tagged)
+        name. Skip-if-exists keys on it, so no two rules share a basename.
     cat_config : str
         Path to catalog configuration YAML
     nside : int
@@ -76,8 +70,9 @@ def generate_pseudo_cl(
     Returns
     -------
     str
-        Path to the primitive's native ``pseudo_cl_{version}.fits`` product.
+        ``out_path`` (the SACC part written).
     """
+    output_dir = os.path.dirname(out_path)
     os.makedirs(output_dir, exist_ok=True)
 
     blind_str = f" blind={blind}" if blind else ""
@@ -135,26 +130,23 @@ def generate_pseudo_cl(
 
     cv = CosmologyValidation(**cv_kwargs)
 
-    # Calculate pseudo-Cls only (no covariance)
-    cv.calculate_pseudo_cl()
+    # Pseudo-Cls only (no covariance), born directly at the final out_path.
+    cv.calculate_pseudo_cl(out_path=out_path)
 
-    # Report on the native product (renamed by the Snakemake caller, if any)
-    src_cl = os.path.join(output_dir, f"pseudo_cl_{version}.fits")
-    if os.path.exists(src_cl):
-        with fits.open(src_cl) as hdul:
-            data = hdul["PSEUDO_CELL"].data
-            n_ell = len(data["ELL"])
-            print(f"Generated pseudo-Cl with {n_ell} ell bins")
-            print(f"ell range: [{data['ELL'].min():.1f}, {data['ELL'].max():.1f}]")
-    return src_cl
+    if os.path.exists(out_path):
+        # Readback of the part just written — a legitimate pre-blind consumer.
+        s = sacc_io.load(out_path, allow_unblinded=True)
+        ell = sacc_io.get_pseudo_cl(s, (0, 0))[0]
+        print(f"Generated pseudo-Cl with {len(ell)} ell bins")
+        print(f"ell range: [{ell.min():.1f}, {ell.max():.1f}]")
+    return out_path
 
 
 def _from_snakemake(smk):
     p = smk.params
-    output_cl = smk.output.pseudo_cl
-    src_cl = generate_pseudo_cl(
+    generate_pseudo_cl(
         version=p["version"],
-        output_dir=os.path.dirname(output_cl),
+        out_path=smk.output.pseudo_cl,
         cat_config=p["cat_config"],
         nside=int(p["nside"]),
         npatch=int(p["npatch"]),
@@ -164,10 +156,6 @@ def _from_snakemake(smk):
         nbins=int(p["nbins"]),
         power=float(p.get("power", 0.5)),
     )
-    # Snakemake declares a tagged output filename; rename the native product to it.
-    if os.path.exists(src_cl) and src_cl != output_cl:
-        os.rename(src_cl, output_cl)
-        print(f"Saved to: {output_cl}")
 
 
 def _from_cli(argv=None):
@@ -220,9 +208,12 @@ def _from_cli(argv=None):
         with open(a.cosmo_json) as f:
             cosmo_params = json.load(f)
 
+    # lc/ASTRA path: --out is a per-recipe directory, so the untagged name is
+    # unambiguous there.
+    out_path = os.path.join(a.out, f"pseudo_cl_{a.ver}.sacc")
     generate_pseudo_cl(
         version=a.ver,
-        output_dir=a.out,
+        out_path=out_path,
         cat_config=a.cat_config,
         nside=a.nside,
         npatch=a.npatch,
diff --git a/workflow/scripts/run_2pcf.py b/workflow/scripts/run_2pcf.py
index 2e1ccabf..78281c73 100644
--- a/workflow/scripts/run_2pcf.py
+++ b/workflow/scripts/run_2pcf.py
@@ -12,15 +12,26 @@
         --cat-config /path/to/cosmo_val/cat_config.yaml \
         --out <output_dir>
 
-The measurement itself is unchanged — ``CosmologyValidation.calculate_2pcf``
-does the TreeCorr work and writes the ``.txt`` dump plus ξ+/ξ- FITS files into
-``output_dir``. ``output_dir`` is passed explicitly (rather than via the
-``COSMO_VAL`` env hook) so lc can point each run at its own ``{output}`` tree.
+The measurement is binning-agnostic: the reporting and the fine integration
+grids are the same compute with different ``--min-sep/--max-sep/--nbins``.
+``CosmologyValidation.calculate_2pcf`` writes the ``.txt`` dump (a raw
+byproduct); the ξ± data product is born as SACC here, a *part* named by its
+binning and tagged with its ``--grid``. The part carries the covariance its
+grid configures (``--cov``): the dense jackknife estimate from the patches, the
+TreeCorr ``varxip``/``varxim`` diagonal, or none.
+
+``output_dir`` is passed explicitly (rather than via the ``COSMO_VAL`` env hook)
+so lc can point each run at its own ``{output}`` tree.
 """
 
 import argparse
+import os
+
+import numpy as np
 
+from sp_validation import sacc_io
 from sp_validation.cosmo_val import CosmologyValidation
+from sp_validation.cosmo_val.sacc_writers import xi_to_sacc
 
 
 def run_2pcf(
@@ -31,30 +42,68 @@ def run_2pcf(
     npatch,
     cat_config,
     output_dir,
-    save_fits=True,
+    sacc_out=None,
+    grid="reporting",
+    cov="none",
 ):
-    """Measure ξ±(θ) for ``ver`` and write it under ``output_dir``.
+    """Measure ξ±(θ) for ``ver`` and write its reporting SACC part.
 
     Parameters mirror the TreeCorr reporting/integration grids: ``min_sep`` /
     ``max_sep`` in arcmin, ``nbins`` logarithmic bins, ``npatch`` spatial
     patches (1 for the paper fiducial). ``cat_config`` is an absolute path to
     the catalog configuration; ``output_dir`` overrides
-    ``cat_config['paths']['output']`` so products land where lc expects.
+    ``cat_config['paths']['output']`` so the ``.txt`` byproduct lands where lc
+    expects. ``sacc_out`` is the exact destination for the SACC part (the
+    Snakemake-declared output); it defaults to a binning-derived name under
+    the resolved output directory for the CLI path.
+
+    Returns
+    -------
+    treecorr.GGCorrelation
+        The measured correlation object (also the source of the SACC part).
     """
     cv = CosmologyValidation(
         versions=[ver],
         catalog_config=cat_config,
         output_dir=output_dir,
+        # so the SACC provenance metadata stamps the npatch actually measured
+        npatch=npatch,
     )
-    return cv.calculate_2pcf(
+    gg = cv.calculate_2pcf(
         ver=ver,
         npatch=npatch,
-        save_fits=save_fits,
         min_sep=min_sep,
         max_sep=max_sep,
         nbins=nbins,
     )
 
+    if cov == "jackknife" and int(npatch) < 2:
+        raise ValueError(f"cov='jackknife' needs patches; got npatch={npatch}")
+
+    # Born-as-SACC ξ± part. theta = meanr; theta_nom = rnom.
+    s = xi_to_sacc(
+        cv.sacc_nz(ver),
+        cv.sacc_metadata(ver),
+        gg.meanr,
+        gg.xip,
+        gg.xim,
+        grid=grid,
+        theta_nom=gg.rnom,
+        npairs=gg.npairs,
+        weight=gg.weight,
+        covariance=gg.cov if cov == "jackknife" else None,
+        variances=(
+            np.concatenate([gg.varxip, gg.varxim]) if cov == "diagonal" else None
+        ),
+    )
+    out_path = sacc_out or os.path.join(
+        output_dir or cv.cc["paths"]["output"],
+        f"{ver}_xi_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.sacc",
+    )
+    sacc_io.save(s, out_path, type="data")
+    print(f"Wrote {grid} ξ± SACC part: {out_path}")
+    return gg
+
 
 def _from_snakemake(smk):
     p = smk.params
@@ -70,7 +119,11 @@ def _from_snakemake(smk):
         # class defaults (./cat_config.yaml, COSMO_VAL env) otherwise.
         cat_config=p.get("cat_config", "./cat_config.yaml"),
         output_dir=p.get("output_dir", None),
-        save_fits=True,
+        grid=p.get("grid", "reporting"),
+        cov=p.get("cov", "none"),
+        # The SACC part goes exactly where the rule declares it; the .txt
+        # byproduct still lands under the resolved output dir.
+        sacc_out=smk.output["sacc"],
     )
 
 
@@ -97,7 +150,15 @@ def _from_cli(argv=None):
         "--cat-config", required=True, help="Absolute path to cat_config.yaml"
     )
     ap.add_argument("--out", required=True, help="Output directory (lc {output})")
-    ap.add_argument("--no-fits", action="store_true", help="Skip ξ+/ξ- FITS export")
+    ap.add_argument(
+        "--grid", default="reporting", help="SACC grid tag for the measured points"
+    )
+    ap.add_argument(
+        "--cov",
+        default="none",
+        choices=["jackknife", "diagonal", "none"],
+        help="Covariance the part carries",
+    )
     a = ap.parse_args(argv)
     run_2pcf(
         ver=a.ver,
@@ -107,7 +168,8 @@ def _from_cli(argv=None):
         npatch=a.npatch,
         cat_config=a.cat_config,
         output_dir=a.out,
-        save_fits=not a.no_fits,
+        grid=a.grid,
+        cov=a.cov,
     )
 
 
diff --git a/workflow/scripts/run_2pcf_highres.py b/workflow/scripts/run_2pcf_highres.py
deleted file mode 100644
index eb31d2c7..00000000
--- a/workflow/scripts/run_2pcf_highres.py
+++ /dev/null
@@ -1,394 +0,0 @@
-#!/usr/bin/env python
-"""
-High-resolution ξ± measurement for COSEBIS integration.
-
-Computes TreeCorr GGCorrelation with fine angular binning (10,000+ bins)
-required for accurate COSEBIS mode integration. Uses MPI for patch-pair
-distribution across nodes when available; falls back to multi-threaded
-single-process otherwise.
-
-Reference: Asgari et al. 2017 — minimum 10,000 bins for E_7 at 0.5% accuracy.
-
-Usage:
-    # MPI (via Slurm submission script):
-    mpiexec --map-by ppr:1:node python run_2pcf_highres.py \
-        --cat-config /path/to/cosmo_val/cat_config.yaml --out <output_dir>
-
-    # Single-process fallback:
-    python run_2pcf_highres.py \
-        --cat-config /path/to/cosmo_val/cat_config.yaml --out <output_dir>
-"""
-
-import argparse
-import os
-import time
-
-import numpy as np
-import treecorr
-from astropy.io import fits
-
-try:
-    # In-container path: full sp_validation stack available.
-    from sp_validation.cosmo_val import CosmologyValidation
-
-    _HAVE_COSMO_VAL = True
-except ImportError:
-    # Bare-host path (host OpenMPI + host python for the 10k-bin MPI run): the
-    # full sp_validation stack (cs_util.plots -> healpy/healsparse) is not
-    # installed. This measurement only needs the shear catalog path + column
-    # names, which are a pure cat_config.yaml lookup — resolve them standalone.
-    CosmologyValidation = None
-    _HAVE_COSMO_VAL = False
-
-# ---------------------------------------------------------------------------
-# MPI setup (graceful fallback)
-# ---------------------------------------------------------------------------
-try:
-    from mpi4py import MPI
-
-    comm = MPI.COMM_WORLD
-    rank = comm.Get_rank()
-    size = comm.Get_size()
-    USE_MPI = size > 1
-except ImportError:
-    comm = None
-    rank = 0
-    size = 1
-    USE_MPI = False
-
-# ---------------------------------------------------------------------------
-# Configuration
-# ---------------------------------------------------------------------------
-# Shear response (R=1 for all SP catalogs)
-R = 1.0
-
-# Detect threads from Slurm or fall back to OS count
-NUM_THREADS = int(os.environ.get("SLURM_CPUS_PER_TASK", os.cpu_count() or 24))
-
-# The catalog path, ellipticity/weight columns, TreeCorr grid, patch count and
-# output directory are resolved from the CLI in main() (defaults reproduce the
-# historical hardcoded values for a no-arg run). They are declared here as
-# module globals so the rank-aware helpers below resolve them at call time; the
-# catalog path + columns come from cat_config + version exactly as run_2pcf.py
-# resolves them (via CosmologyValidation).
-CAT_PATH = None
-VERSION = None
-E1_COL = None
-E2_COL = None
-W_COL = None
-TMIN = None  # arcmin
-TMAX = None  # arcmin
-NBINS = None
-NPATCH = None
-OUTPUT_DIR = None
-PATCH_FILE = None
-
-
-def parse_args(argv=None):
-    """CLI mirroring run_xi_sweep's signature; defaults reproduce prior behavior."""
-    ap = argparse.ArgumentParser(
-        description="High-resolution TreeCorr ξ± measurement for COSEBIS integration."
-    )
-    ap.add_argument(
-        "--config",
-        default=None,
-        help="Path to bmodes config.yaml (accepted for signature parity with "
-        "run_xi_sweep; not read by this measurement).",
-    )
-    ap.add_argument(
-        "--cat-config", required=True, help="Absolute path to cat_config.yaml"
-    )
-    ap.add_argument(
-        "--version",
-        default="SP_v1.4.6.3_leak_corr",
-        help="Catalog version key in cat_config",
-    )
-    ap.add_argument("--nbins", type=int, default=10000, help="Number of log bins")
-    ap.add_argument("--npatch", type=int, default=50, help="TreeCorr patch count")
-    ap.add_argument(
-        "--min-sep", type=float, default=0.5, help="Min separation [arcmin]"
-    )
-    ap.add_argument(
-        "--max-sep", type=float, default=300.0, help="Max separation [arcmin]"
-    )
-    ap.add_argument("--out", required=True, help="Output directory (lc {output})")
-    return ap.parse_args(argv)
-
-
-def log(msg):
-    """Print with timestamp on rank 0 only."""
-    if rank == 0:
-        print(f"[{time.strftime('%H:%M:%S')}] {msg}", flush=True)
-
-
-def load_catalog():
-    """Load shear catalog and apply mean subtraction."""
-    log(f"Loading catalog: {CAT_PATH}")
-    hdul = fits.open(CAT_PATH, memmap=True)
-    data = hdul[1].data
-
-    ra = np.array(data["ra"], dtype=np.float64)
-    dec = np.array(data["dec"], dtype=np.float64)
-    e1 = np.array(data[E1_COL], dtype=np.float64)
-    e2 = np.array(data[E2_COL], dtype=np.float64)
-    w = np.array(data[W_COL], dtype=np.float64)
-    hdul.close()
-
-    log(f"  {len(ra):,} galaxies loaded")
-
-    # Additive bias: c = <e/R>_w  (R=1 for SP catalogs)
-    c1 = np.average(e1 / R, weights=w)
-    c2 = np.average(e2 / R, weights=w)
-    log(f"  Additive bias: c1={c1:.6e}, c2={c2:.6e}")
-
-    # Calibrated shear: g = (e - c) / R
-    g1 = (e1 - c1) / R
-    g2 = (e2 - c2) / R
-
-    return ra, dec, g1, g2, w
-
-
-def _wait_for_file(path, timeout=300, interval=1.0):
-    """Block until `path` is visible to this node, defeating NFS dir caching.
-
-    On a multi-node run the rank that wrote `path` sees it immediately, but
-    peer nodes can carry a stale negative directory-cache entry past an MPI
-    Barrier. Re-listing the parent directory forces an NFS attribute refresh;
-    poll that until the entry appears (or raise after `timeout` seconds).
-    """
-    parent = os.path.dirname(path) or "."
-    name = os.path.basename(path)
-    waited = 0.0
-    while waited < timeout:
-        try:
-            if name in os.listdir(parent):
-                return
-        except FileNotFoundError:
-            pass
-        time.sleep(interval)
-        waited += interval
-    raise TimeoutError(f"patch-center file not visible after {timeout}s: {path}")
-
-
-def compute_patch_centers(ra, dec):
-    """Compute patch centers from subsampled catalog (rank 0 only)."""
-    if os.path.exists(PATCH_FILE):
-        log(f"Using existing patch centers: {PATCH_FILE}")
-        return
-
-    if rank != 0:
-        return
-
-    log(f"Computing patch centers (npatch={NPATCH}) from 1% subsample...")
-    rng = np.random.default_rng(42)
-    n_sub = max(len(ra) // 100, NPATCH * 100)
-    idx = rng.choice(len(ra), size=n_sub, replace=False)
-
-    cat_sub = treecorr.Catalog(
-        ra=ra[idx],
-        dec=dec[idx],
-        ra_units="degrees",
-        dec_units="degrees",
-        npatch=NPATCH,
-    )
-    cat_sub.write_patch_centers(PATCH_FILE)
-    log(f"  Wrote patch centers to {PATCH_FILE}")
-    del cat_sub
-
-
-def write_xi_fits(gg, prefix, xi_data):
-    """Write ξ+ or ξ- to FITS matching CosmologyValidation format."""
-    out_path = os.path.join(
-        OUTPUT_DIR,
-        f"{prefix}_{VERSION}_minsep={TMIN}_maxsep={TMAX}_nbins={NBINS}_npatch=1.fits",
-    )
-    n = len(xi_data)
-    cols = [
-        fits.Column(name="BIN1", format="K", array=np.ones(n, dtype=int)),
-        fits.Column(name="BIN2", format="K", array=np.ones(n, dtype=int)),
-        fits.Column(name="ANGBIN", format="K", array=np.arange(1, n + 1)),
-        fits.Column(name="VALUE", format="D", array=xi_data),
-        fits.Column(name="ANG", format="D", unit="arcmin", array=gg.meanr),
-    ]
-    ext_name = "XI_PLUS" if "plus" in prefix else "XI_MINUS"
-    hdu = fits.BinTableHDU.from_columns(cols, name=ext_name)
-    for key, val in {
-        "2PTDATA": "T",
-        "QUANT1": "G+R",
-        "QUANT2": "G+R",
-        "KERNEL_1": "NZ_SOURCE",
-        "KERNEL_2": "NZ_SOURCE",
-        "WINDOWS": "SAMPLE",
-    }.items():
-        hdu.header[key] = val
-    hdu.writeto(out_path, overwrite=True)
-    log(f"  Wrote {out_path}")
-
-
-def resolve_shear_config(cat_config_path, version):
-    """Standalone shear-config resolver (bare-host fallback for CosmologyValidation).
-
-    Reproduces exactly the ``cc[version]["shear"]`` fields this measurement reads
-    (path, e1_col, e2_col, w_col), replicating CosmologyValidation's two
-    transforms: (1) subdir-relative path resolution, and (2) the ``_leak_corr``
-    virtual version — deep-copy the base version and swap
-    e1_col/e2_col -> e1_col_corrected/e2_col_corrected. See
-    sp_validation/cosmo_val/core.py.
-    """
-    import copy
-
-    import yaml
-
-    with open(cat_config_path) as fh:
-        cc = yaml.load(fh, Loader=yaml.FullLoader)
-
-    def resolve_paths(ver):
-        subdir = os.fspath(cc[ver]["subdir"])
-        for section in cc[ver].values():
-            if isinstance(section, dict) and "path" in section:
-                p = section["path"]
-                if not os.path.isabs(p):
-                    section["path"] = os.path.join(subdir, p)
-
-    leak_suffix = "_leak_corr"
-    if version in cc:
-        resolve_paths(version)
-    elif version.endswith(leak_suffix):
-        base = version[: -len(leak_suffix)]
-        if base not in cc:
-            raise ValueError(f"Base version '{base}' not in cat_config for '{version}'")
-        resolve_paths(base)
-        base_shear = cc[base]["shear"]
-        if "e1_col_corrected" not in base_shear or "e2_col_corrected" not in base_shear:
-            raise ValueError(
-                f"{base} lacks e1_col_corrected/e2_col_corrected; cannot form {version}"
-            )
-        cc[version] = copy.deepcopy(cc[base])
-        cc[version]["shear"]["e1_col"] = base_shear["e1_col_corrected"]
-        cc[version]["shear"]["e2_col"] = base_shear["e2_col_corrected"]
-        resolve_paths(version)
-    else:
-        raise ValueError(f"Version '{version}' not found in cat_config")
-
-    return cc[version]["shear"]
-
-
-def main():
-    global CAT_PATH, VERSION, E1_COL, E2_COL, W_COL
-    global TMIN, TMAX, NBINS, NPATCH, OUTPUT_DIR, PATCH_FILE
-
-    args = parse_args()
-    VERSION = args.version
-    NBINS = args.nbins
-    NPATCH = args.npatch
-    TMIN = args.min_sep
-    TMAX = args.max_sep
-    OUTPUT_DIR = args.out
-
-    # Resolve catalog path + ellipticity/weight columns from cat_config + version
-    # exactly as run_2pcf.py does (applies the _leak_corr column swap and the
-    # subdir path resolution). In-container this uses CosmologyValidation; bare-
-    # host (10k-bin MPI run) it uses the standalone cat_config resolver, which is
-    # byte-identical for the shear-config fields this measurement reads.
-    if _HAVE_COSMO_VAL:
-        cv = CosmologyValidation(
-            versions=[VERSION], catalog_config=args.cat_config, output_dir=OUTPUT_DIR
-        )
-        shear_cfg = cv.cc[VERSION]["shear"]
-    else:
-        shear_cfg = resolve_shear_config(args.cat_config, VERSION)
-    CAT_PATH = shear_cfg["path"]
-    E1_COL = shear_cfg["e1_col"]
-    E2_COL = shear_cfg["e2_col"]
-    W_COL = shear_cfg["w_col"]
-
-    PATCH_FILE = os.path.join(
-        OUTPUT_DIR,
-        f"patch_centers_{VERSION}_{NPATCH}_{TMIN}_{TMAX}.dat",
-    )
-
-    t0 = time.time()
-
-    log("=" * 60)
-    log("High-resolution ξ± measurement")
-    log(f"  MPI: {'yes' if USE_MPI else 'no'} (ranks={size})")
-    log(f"  Config: {NBINS:,} bins, [{TMIN}, {TMAX}] arcmin")
-    log(f"  Patches: {NPATCH}, Threads/rank: {NUM_THREADS}")
-    log(f"  Version: {VERSION}")
-    log("=" * 60)
-
-    # All ranks load catalog (needed for TreeCorr patch assignment)
-    ra, dec, g1, g2, w = load_catalog()
-
-    # Compute patch centers (rank 0 only; others wait)
-    compute_patch_centers(ra, dec)
-    if USE_MPI:
-        comm.Barrier()
-        # Cross-node visibility: rank 0 wrote PATCH_FILE on its node, but on a
-        # multi-node allocation the other ranks' nodes may not see it yet (NFS
-        # close-to-open + negative-dir caching persists past the Barrier). Poll
-        # with a forced directory refresh until it appears before reading it.
-        _wait_for_file(PATCH_FILE)
-
-    # Create TreeCorr catalog with patch centers
-    log("Creating TreeCorr catalog with patches...")
-    cat = treecorr.Catalog(
-        ra=ra,
-        dec=dec,
-        g1=g1,
-        g2=g2,
-        w=w,
-        ra_units="degrees",
-        dec_units="degrees",
-        patch_centers=PATCH_FILE,
-    )
-    cat.load()
-    cat.get_patches()
-    log(f"  Catalog ready ({cat.nobj:,} objects, {cat.npatch} patches)")
-
-    # Free raw arrays (TreeCorr holds its own copy)
-    del ra, dec, g1, g2, w
-
-    # Compute GG correlation
-    log("Computing GGCorrelation...")
-    gg = treecorr.GGCorrelation(
-        min_sep=TMIN,
-        max_sep=TMAX,
-        nbins=NBINS,
-        sep_units="arcminutes",
-        verbose=2,
-    )
-
-    process_kwargs = {"num_threads": NUM_THREADS}
-    if USE_MPI:
-        process_kwargs["comm"] = comm
-
-    gg.process(cat, **process_kwargs)
-    log(f"  Correlation complete ({time.time() - t0:.0f}s elapsed)")
-
-    # Write output (rank 0 only)
-    if rank == 0:
-        out_txt = os.path.join(
-            OUTPUT_DIR,
-            f"{VERSION}_xi_minsep={TMIN}_maxsep={TMAX}_nbins={NBINS}_npatch=1.txt",
-        )
-        # Write only the main per-bin correlation. The convergence consumer
-        # (cosebis_binning_comparison.py) reads just the per-bin columns
-        # (np.loadtxt max_rows=nbins) and the 1000-bin covariance — the 10k
-        # jackknife cov is used nowhere. write_patch_results/write_cov=True
-        # serialised a 20000x20000 cov + 180 patch blocks (~10 GB) that nothing
-        # reads and also cost the estimate_cov compute; drop both. The patches
-        # still parallelise gg.process; gg.xip/gg.xim (values, FITS) are
-        # unaffected.
-        gg.write(out_txt, write_patch_results=False, write_cov=False)
-        log(f"  Wrote {out_txt}")
-
-        write_xi_fits(gg, "xi_plus", gg.xip)
-        write_xi_fits(gg, "xi_minus", gg.xim)
-
-        elapsed = time.time() - t0
-        log(f"Done! Total time: {elapsed / 3600:.1f}h ({elapsed:.0f}s)")
-
-
-if __name__ == "__main__":
-    main()
diff --git a/workflow/scripts/run_cosmocov_chain.sh b/workflow/scripts/run_cosmocov_chain.sh
new file mode 100644
index 00000000..3a7d3a84
--- /dev/null
+++ b/workflow/scripts/run_cosmocov_chain.sh
@@ -0,0 +1,93 @@
+#!/usr/bin/env bash
+# CosmoCov covariance chain (lc-native, container:none recipe).
+#
+# covariance_ini -> covariance_cosmocov (x3 blocks) -> covariance_cat ->
+# covariance_process. The CosmoCov C++ binary runs on the bare host (module load
+# gcc/intelpython/openmpi); the .ini generation and cosmocov_process steps run
+# inside the sp_validation apptainer container. The 3 shear-shear blocks
+# (++,--,+-) are independent and run in parallel.
+#
+# Usage:
+#   run_cosmocov_chain.sh --version SP_v1.4.6.3_leak_corr --blind A \
+#     --min-sep 0.5 --max-sep 300.0 --nbins 1000 --gaussian g \
+#     --planck18-json <cosmology_snapshot>/planck18.json \
+#     --cat-config <cat_config.yaml> --mask-cls <mask_cls...norm.txt> \
+#     --out <output_dir>
+#
+# The checkout is this script's own (workflow/scripts/../..). Deployment paths
+# come from the environment, with the current candide values as defaults:
+#   SPV_CONTAINER  apptainer image        (default /n17data/cdaley/containers/containers/)
+#   SPV_BIND       apptainer --bind list
+#   COSMOCOV       CosmoCov `cov` binary  (also settable with --cosmocov)
+set -euo pipefail
+
+WT="$(cd "$(dirname "${BASH_SOURCE[0]}")/../.." && pwd)"
+SRC=$WT/src
+CONTAINER=${SPV_CONTAINER:-/n17data/cdaley/containers/containers/}
+BIND=${SPV_BIND:-/home,/scratch,/automnt,/n17data,/n23data1,/n09data}
+COSMOCOV=${COSMOCOV:-/n23data1/n06data/lgoh/scratch/UNIONS/CosmoCov/covs/cov}
+
+VERSION=""; BLIND="A"; MINSEP=""; MAXSEP=""; NBINS=""; GAUSSIAN=""
+PLANCK18=""; CATCONFIG=""; MASKCLS=""; OUT=""
+while [ $# -gt 0 ]; do
+  case "$1" in
+    --version) VERSION="$2"; shift 2;;
+    --blind) BLIND="$2"; shift 2;;
+    --min-sep) MINSEP="$2"; shift 2;;
+    --max-sep) MAXSEP="$2"; shift 2;;
+    --nbins) NBINS="$2"; shift 2;;
+    --gaussian) GAUSSIAN="$2"; shift 2;;
+    --planck18-json) PLANCK18="$2"; shift 2;;
+    --cat-config) CATCONFIG="$2"; shift 2;;
+    --mask-cls) MASKCLS="$2"; shift 2;;
+    --cosmocov) COSMOCOV="$2"; shift 2;;
+    --out) OUT="$2"; shift 2;;
+    *) echo "unknown arg: $1" >&2; exit 2;;
+  esac
+done
+
+mkdir -p "$OUT"
+# Absolutize OUT before the `cd "$OUT"` below (the CosmoCov binary writes its
+# blocks into cwd): every other OUT-relative path would otherwise re-resolve
+# against the new cwd and double-nest. lc templates {output} as a
+# project-relative path, so both relative and absolute --out must work.
+OUT="$(cd "$OUT" && pwd)"
+INI="$OUT/covariance.ini"
+
+echo "[cosmocov] generating .ini"
+apptainer exec --bind "$BIND" --env PYTHONPATH="$SRC" "$CONTAINER" \
+  /usr/local/bin/python "$WT/workflow/scripts/generate_cosmocov_ini.py" \
+  --version "$VERSION" --blind "$BLIND" \
+  --planck18-json "$PLANCK18" --cat-config "$CATCONFIG" \
+  --min-sep "$MINSEP" --max-sep "$MAXSEP" --nbins "$NBINS" --gaussian "$GAUSSIAN" \
+  --mask-cls "$MASKCLS" --out-ini "$INI"
+
+echo "[cosmocov] loading modules + running 3 blocks (parallel)"
+source /etc/profile.d/modules.sh
+module unload gcc 2>/dev/null || true; module load gcc
+module unload intelpython 2>/dev/null || true; module load intelpython/3-2024.1.0
+module load openmpi
+
+cd "$OUT"
+# One CosmoCov invocation per block; see common.py BLOCK_PAIRS.
+for idx in 1 2 3; do
+  ( "$COSMOCOV" "$idx" "$INI" > "$OUT/cosmocov_block_${idx}.log" 2>&1 ) &
+done
+wait
+
+# Concatenate blocks in BLOCK_PAIRS order (++, --, +-), as covariance_cat does.
+CAT="$OUT/covariance.txt"
+: > "$CAT"
+for pm_idx in "++:1" "--:2" "+-:3"; do
+  pm="${pm_idx%%:*}"; idx="${pm_idx##*:}"
+  blk="$OUT/cov_tmp_ssss_${pm}_cov_Ntheta${NBINS}_Ntomo1_${idx}"
+  [ -f "$blk" ] || { echo "MISSING block $blk (see cosmocov_block_${idx}.log)" >&2; exit 1; }
+  cat "$blk" >> "$CAT"
+done
+echo "[cosmocov] concatenated -> $CAT"
+
+echo "[cosmocov] processing (positive-definite check, G/G+NG extract, QA plot)"
+apptainer exec --bind "$BIND" --env PYTHONPATH="$SRC" "$CONTAINER" \
+  /usr/local/bin/python "$WT/cosmo_inference/scripts/cosmocov_process.py" \
+  "$CAT" "$OUT/covariance_processed"
+echo "[cosmocov] done -> $OUT/covariance_processed.txt (+_g.txt, +_plot.pdf)"
diff --git a/workflow/scripts/run_rho_tau.py b/workflow/scripts/run_rho_tau.py
index ea2f35bc..71e3deb7 100644
--- a/workflow/scripts/run_rho_tau.py
+++ b/workflow/scripts/run_rho_tau.py
@@ -48,9 +48,10 @@
 
 cv.calculate_rho_tau_stats()
 
-# Confirm CosmologyValidation produced the requested outputs
+# Confirm CosmologyValidation produced the requested outputs: the rho/tau FITS
+# and the born-as-SACC rho_tau part.
 outputs = snakemake.output  # type: ignore
-for label in ("rho_stats", "tau_stats"):
+for label in ("rho_stats", "tau_stats", "rho_tau"):
     target = Path(outputs[label])
     if not target.exists():
         raise FileNotFoundError(