diff --git a/src/sp_validation/cosmo_val/core.py b/src/sp_validation/cosmo_val/core.py
index 8159a1a0..f3cca848 100644
--- a/src/sp_validation/cosmo_val/core.py
+++ b/src/sp_validation/cosmo_val/core.py
@@ -471,8 +471,10 @@ def ensure_version_exists(ver):
         # B-mode results storage for summarize_bmodes()
         self._pure_eb_results = {}
         self._cosebis_results = {}
-        # The sealed ξ± parts calculate_2pcf_version made, by (version, grid).
+        # The sealed ξ± parts calculate_2pcf_version made, by (version, grid):
+        # the ("all", "all") pair, and every tomographic bin pair.
         self.xi_parts = {}
+        self.xi_tomo_parts = {}
 
     def _output_path(self, *parts):
         """Absolute path under the catalog config's output directory.
@@ -484,13 +486,18 @@ 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.
+    def sacc_nz(self, version, tomography=False):
+        """``nz`` mapping ``{i: (z, nz_i)}`` for the SACC writers.
 
-        The SACC parts carry the ("all", "all") pair only, so the whole-survey
-        n(z) of :meth:`get_redshift` is bin 0.
+        Non-tomographic parts carry the ("all", "all") pair only, so the
+        whole-survey n(z) of :meth:`get_redshift` is bin 0. With
+        ``tomography``, bin ``i`` is column ``i + 1`` of the ``redshift_path``
+        file (column 0 is z), one per tomographic bin.
         """
-        return {0: tuple(self.get_redshift(version))}
+        if not tomography:
+            return {0: tuple(self.get_redshift(version))}
+        z, dndz = self.read_redshift_distribution(version, is_tomography=True)
+        return {i: (z, dndz[:, i]) for i in range(dndz.shape[1])}
 
     def sacc_metadata(self, version):
         """Provenance metadata stored on every SACC part for ``version``."""
@@ -601,17 +608,19 @@ def blind(self, version):
     def _refuse_blinded_tomography(self, version):
         """Refuse to measure tomographic shear signal on a blinded catalogue.
 
-        Only the ``("all", "all")`` pair is born as a SACC part, whose tracer
-        carries the n(z) the blind's theory shifts it by. The tomographic ξ±,
-        ⟨M_ap²⟩ and pseudo-Cℓ are TreeCorr objects and FITS tables with no
-        per-bin tracer, so on a blinded catalogue they would hold true values.
+        A statistic is blinded by sealing its SACC part, whose tracers carry
+        the n(z) the blind's theory shifts it by. ξ± (and ⟨M_ap²⟩, derived from
+        it) is born as such a part for every tomographic bin pair; the
+        tomographic pseudo-Cℓ are still FITS tables with no per-bin tracer, so
+        on a blinded catalogue they would hold true values.
         """
         blind = self.blind(version)
         if blind.name != blinding.NONE.name:
             raise blinding.BlindingError(
                 f"{version} is blinded under {blind.name}, and blinding covers "
-                "its non-tomographic pair only; tomographic shear statistics are "
-                "measured on public (`blind: none`) catalogues"
+                "its tomographic ξ± and non-tomographic statistics only; other "
+                "tomographic shear statistics are measured on public "
+                "(`blind: none`) catalogues"
             )
 
     def basename(
diff --git a/src/sp_validation/cosmo_val/real_space.py b/src/sp_validation/cosmo_val/real_space.py
index a486e8b9..8fde703f 100644
--- a/src/sp_validation/cosmo_val/real_space.py
+++ b/src/sp_validation/cosmo_val/real_space.py
@@ -14,9 +14,7 @@
 from sp_validation.statistics import jackknife_patch_centers
 
 from .. import sacc_io
-from .sacc_writers import xi_to_sacc
-
-ALL_PAIR = "tomo_bin_all_tomo_bin_all"
+from .sacc_writers import xi_pairs_to_sacc
 
 
 class RealSpaceMixin:
@@ -32,23 +30,28 @@ def calculate_2pcf_version(
     ):
         """ξ± of one catalogue version on one binning, for each bin pair.
 
-        The non-tomographic ``("all", "all")`` pair is born as its SACC part,
-        sealed under the version's blind (:func:`sp_validation.sacc_io.seal`)
-        before it is kept in ``self.xi_parts[ver, grid]``, returned or written,
-        so on a blinded catalogue its true values never leave this method. A
-        part already held for ``(ver, grid)`` is reused when no ``npatch``,
-        binning or ``out`` is asked for, so the figure rules draw the parts they
-        are handed. Tomographic pairs are measured only on a public catalogue
-        (:meth:`_refuse_blinded_tomography`). Use :meth:`calculate_2pcf` to run
-        over every version in ``self.versions``.
+        The measured pairs are born as one SACC part — the ``("all", "all")``
+        pair, or every tomographic bin pair with one ``source_i`` tracer and
+        n(z) per bin — sealed under the version's blind
+        (:func:`sp_validation.sacc_io.seal`) before it is kept in
+        ``self.xi_parts[ver, grid]`` (``self.xi_tomo_parts`` for tomography),
+        returned or written, so on a blinded catalogue its true values never
+        leave this method. A part already held for ``(ver, grid)`` is reused
+        when no ``npatch``, binning or ``out`` is asked for, so the figure rules
+        draw the parts they are handed. Use :meth:`calculate_2pcf` to run over
+        every version in ``self.versions``.
+
+        Tomographic bin ids map to SACC bins in ascending order, so bin ``i``
+        takes column ``i + 1`` of the ``redshift_path`` file, and pairs go in
+        as ``(0, 0), (0, 1), …, (n-1, n-1)`` — OneCovariance's bin-pair order.
 
         Parameters:
             ver (str): The catalogue version to measure.
             npatch (int, optional): Jackknife patches; the instance's
                 ``npatch`` by default. With patches the part carries the
-                jackknife covariance, without them the shot-noise diagonal.
-                Seeded patch centres come from the full catalogue and are
-                shared by every bin pair.
+                jackknife covariance (across bin pairs too), without them the
+                shot-noise diagonal. Seeded patch centres come from the full
+                catalogue and are shared by every bin pair.
             compute_tomography (bool, optional): Measure the tomographic bin
                 pairs instead of ``("all", "all")``.
             grid (str): The grid tag the part's points carry.
@@ -57,39 +60,32 @@ def calculate_2pcf_version(
                 e.g. ``min_sep=1``.
 
         Returns:
-            dict: ``"tomo_bin_{b1}_tomo_bin_{b2}"`` → ξ± shaped like a TreeCorr
-            ``GGCorrelation``: the part's
-            :func:`~sp_validation.sacc_io.xi_correlation` view for
-            ``("all", "all")``, the TreeCorr object for a tomographic pair.
+            dict: ``"tomo_bin_{b1}_tomo_bin_{b2}"`` → the part's
+            :func:`~sp_validation.sacc_io.xi_correlation` view of that pair,
+            shaped like a TreeCorr ``GGCorrelation``.
         """
         reuse = npatch is None and out is None and not treecorr_config
         npatch = int(npatch or self.npatch)
         blind = self.blind(ver)
+        parts = self.xi_tomo_parts if compute_tomography else self.xi_parts
 
         if compute_tomography:
-            self._refuse_blinded_tomography(ver)
             tomo_bin_ids, tomo_bin_pairs = self._get_tomo_bins(ver)
             if tomo_bin_ids is None or tomo_bin_pairs is None:
                 raise ValueError(f"Version {ver} does not have tomography information.")
+            index = {b: i for i, b in enumerate(sorted(tomo_bin_ids))}
             self.print_magenta(
                 f"Computing tomographic ξ± for {ver} with {len(tomo_bin_pairs)} bins."
             )
         else:
             self.print_magenta(f"Computing non-tomographic ξ± for {ver}.")
             tomo_bin_pairs = [("all", "all")]
+            index = {"all": 0}
 
-        ggs = {f"tomo_bin_{b1}_tomo_bin_{b2}": None for b1, b2 in tomo_bin_pairs}
-        held = self.xi_parts.get((ver, grid)) if reuse else None
-        if ALL_PAIR in ggs and held is not None:
-            sacc_io.expect(held, blind)
-            ggs[ALL_PAIR] = sacc_io.xi_correlation(held)
-        to_compute = [
-            (b1, b2)
-            for b1, b2 in tomo_bin_pairs
-            if ggs[f"tomo_bin_{b1}_tomo_bin_{b2}"] is None
-        ]
-
-        if to_compute:
+        part = parts.get((ver, grid)) if reuse else None
+        if part is not None:
+            sacc_io.expect(part, blind)
+        else:
             jackknife = npatch > 1
             gg_config = {
                 **self._binning(**treecorr_config),
@@ -98,7 +94,8 @@ def calculate_2pcf_version(
             cols = self._shear_columns(ver, compute_tomography)
             patch_centers = self._patch_centers(cols, npatch)
 
-            for bin1, bin2 in to_compute:
+            measured = []
+            for bin1, bin2 in tomo_bin_pairs:
                 gg = treecorr.GGCorrelation(gg_config)
 
                 cat_gal1 = self._bin_catalog(cols, bin1, npatch, patch_centers)
@@ -109,31 +106,60 @@ def calculate_2pcf_version(
                 )
 
                 gg.process(cat_gal1, cat2=cat_gal2)
-
-                if (bin1, bin2) != ("all", "all"):
-                    ggs[f"tomo_bin_{bin1}_tomo_bin_{bin2}"] = gg
-                    continue
-                s = xi_to_sacc(
-                    self.sacc_nz(ver),
-                    {**self.sacc_metadata(ver), "npatch": npatch},
-                    gg.meanr,
-                    gg.xip,
-                    gg.xim,
-                    grid=grid,
-                    theta_nom=gg.rnom,
-                    npairs=gg.npairs,
-                    weight=gg.weight,
-                    covariance=gg.cov if jackknife else None,
-                    variances=(
-                        None if jackknife else np.concatenate([gg.varxip, gg.varxim])
-                    ),
+                measured.append(((index[bin1], index[bin2]), gg))
+
+            nz = self.sacc_nz(ver, tomography=compute_tomography)
+            if len(nz) != len(index):
+                raise ValueError(
+                    f"{ver}: the n(z) file has {len(nz)} distribution column(s) "
+                    f"but the catalogue has {len(index)} tomographic bin(s); "
+                    "every bin needs its own n(z)"
                 )
-                if out:
-                    part = sacc_io.save(s, out, blind=blind)
-                else:
-                    part = sacc_io.seal(s, blind)
-                self.xi_parts[ver, grid] = part
-                ggs[ALL_PAIR] = sacc_io.xi_correlation(part)
+            correlations = [gg for _, gg in measured]
+            covariance = variances = None
+            if not jackknife:
+                variances = np.concatenate(
+                    [np.concatenate([gg.varxip, gg.varxim]) for gg in correlations]
+                )
+            elif len(correlations) == 1:
+                covariance = correlations[0].cov
+            else:
+                # [ξ+; ξ−] per correlation, in list order: the part's pair-major
+                # order, with the cross-pair covariance a per-pair gg.cov drops.
+                covariance = treecorr.estimate_multi_cov(correlations, "jackknife")
+            s = xi_pairs_to_sacc(
+                nz,
+                {**self.sacc_metadata(ver), "npatch": npatch},
+                [
+                    (
+                        bins,
+                        {
+                            "theta": gg.meanr,
+                            "xip": gg.xip,
+                            "xim": gg.xim,
+                            "theta_nom": gg.rnom,
+                            "npairs": gg.npairs,
+                            "weight": gg.weight,
+                        },
+                    )
+                    for bins, gg in measured
+                ],
+                grid=grid,
+                covariance=covariance,
+                variances=variances,
+            )
+            if out:
+                part = sacc_io.save(s, out, blind=blind)
+            else:
+                part = sacc_io.seal(s, blind)
+            parts[ver, grid] = part
+
+        ggs = {
+            f"tomo_bin_{b1}_tomo_bin_{b2}": sacc_io.xi_correlation(
+                part, bins=(index[b1], index[b2]), grid=grid
+            )
+            for b1, b2 in tomo_bin_pairs
+        }
 
         self.print_done(f"Done 2PCF for {ver}.")
 
@@ -241,9 +267,9 @@ def calculate_aperture_mass_dispersion(
 
         Both are linear in ξ± (TreeCorr's ``calculateMapSq`` sum, Schneider et
         al. 2002 filter), so they and their covariance T·C·Tᵀ come from the ξ±
-        of :meth:`calculate_2pcf_version` on a fine grid: for ``("all",
-        "all")`` its sealed part, concealed on a blinded catalogue, whose shift,
-        the same in every patch, leaves C unchanged.
+        of :meth:`calculate_2pcf_version` on a fine grid: its sealed part, for
+        every bin pair concealed on a blinded catalogue, whose shift, the same
+        in every patch, leaves C unchanged.
         """
         theta_map = np.geomspace(theta_min * 5, theta_max / 2, nbins_map)
         self._map2 = {"theta_map": theta_map}
diff --git a/src/sp_validation/cosmo_val/sacc_writers.py b/src/sp_validation/cosmo_val/sacc_writers.py
index b10bfe81..f6c86044 100644
--- a/src/sp_validation/cosmo_val/sacc_writers.py
+++ b/src/sp_validation/cosmo_val/sacc_writers.py
@@ -7,8 +7,11 @@
 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.
+ξ± is written over any set of bin pairs (:func:`xi_pairs_to_sacc`), and
+:func:`onecov_xi_covariance` puts OneCovariance's ξ± matrix in a part's
+order. The other statistics are single-bin today (``bins=(0, 0)``); their
+interface is tomography-native so a future round supplies real bin pairs
+unchanged.
 """
 
 import numpy as np
@@ -27,6 +30,39 @@
 BIN = (0, 0)
 
 
+def xi_pairs_to_sacc(nz, metadata, pairs, *, grid, variances=None, covariance=None):
+    """One ξ± part over one or more bin pairs, on a named angular grid.
+
+    ``pairs`` is a sequence of ``(bins, xi)``: ``bins`` a source bin pair
+    ``(i, j)`` and ``xi`` a mapping with ``theta``, ``xip``, ``xim`` and,
+    optionally, ``theta_nom``, ``npairs``, ``weight``. Each pair goes in as
+    ``[ξ+; ξ−]`` in the order given, so the data vector is pair-major, and the
+    covariance is tied to that order by position: ``covariance`` is the dense
+    block in it, ``variances`` the matching diagonal. At most one may be given.
+    ``nz`` must hold an n(z) for every bin the pairs name.
+    """
+    s = sio.new_sacc(nz, metadata)
+    for bins, xi in pairs:
+        sio.add_xi(
+            s,
+            bins,
+            xi["theta"],
+            xi["xip"],
+            xi["xim"],
+            grid=grid,
+            theta_nom=xi.get("theta_nom"),
+            npairs=xi.get("npairs"),
+            weight=xi.get("weight"),
+        )
+    if covariance is not None and variances is not None:
+        raise ValueError("give the ξ± writer 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 xi_to_sacc(
     nz,
     metadata,
@@ -43,29 +79,74 @@ def xi_to_sacc(
 ):
     """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.
+    The single-pair form of :func:`xi_pairs_to_sacc`. 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,
+    xi = {
+        "theta": theta,
+        "xip": xip,
+        "xim": xim,
+        "theta_nom": theta_nom,
+        "npairs": npairs,
+        "weight": weight,
+    }
+    return xi_pairs_to_sacc(
+        nz, metadata, [(BIN, xi)], grid=grid, variances=variances, covariance=covariance
     )
-    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 onecov_xi_covariance(s, cov):
+    """OneCovariance's real-space ξ± matrix ``cov`` in the insertion order of ``s``.
+
+    OneCovariance writes ξ± type-major: every ξ+ point, then every ξ− point;
+    within each, bin pairs ``(i, j)`` with ``i <= j`` and ``i`` outermost, then
+    θ ascending. A ξ± part is pair-major (:func:`xi_pairs_to_sacc`), so this
+    finds each point of ``s`` in OneCovariance's order and permutes ``cov``
+    to match. The two orders agree for a single bin pair. OneCovariance's bin
+    ``k`` is column ``k`` of the n(z) file, which must be the part's
+    ``source_k``.
+
+    ``s`` must hold ξ± on one grid only, with the same θ points for every bin
+    pair, as a part does.
+    """
+    cov = np.asarray(cov, float)
+    types = {sio.XI_PLUS: 0, sio.XI_MINUS: 1}
+    if any(dp.data_type not in types for dp in s.data):
+        raise ValueError("onecov_xi_covariance takes a part holding ξ± only")
+    grids = {dp.tags.get("grid") for dp in s.data}
+    if len(grids) != 1:
+        raise ValueError(f"the part's ξ± lies on several grids {sorted(grids)}")
+    if cov.shape != (len(s.data), len(s.data)):
+        raise ValueError(
+            f"OneCovariance matrix is {cov.shape}, but the part has "
+            f"{len(s.data)} ξ± points; check the binning and tomography match"
+        )
+
+    def bin_index(tracer):
+        return int(tracer.removeprefix("source_"))
+
+    pairs = sorted({tuple(sorted(map(bin_index, dp.tracers))) for dp in s.data})
+    thetas = {}
+    for dp in s.data:
+        key = (dp.data_type, tuple(sorted(map(bin_index, dp.tracers))))
+        thetas.setdefault(key, []).append(dp.tags["theta"])
+    ntheta = {len(t) for t in thetas.values()}
+    if len(ntheta) != 1 or len(thetas) != 2 * len(pairs):
+        raise ValueError("every bin pair must carry ξ+ and ξ− on the same θ points")
+    (ntheta,) = ntheta
+
+    order = np.empty(len(s.data), dtype=int)
+    for k, dp in enumerate(s.data):
+        pair = tuple(sorted(map(bin_index, dp.tracers)))
+        theta_rank = sorted(thetas[dp.data_type, pair]).index(dp.tags["theta"])
+        order[k] = (
+            types[dp.data_type] * len(pairs) * ntheta
+            + pairs.index(pair) * ntheta
+            + theta_rank
+        )
+    return cov[np.ix_(order, order)]
 
 
 def pseudo_cl_to_sacc(nz, metadata, ell_eff, cl_all, wsp, covariance=None, nside=None):
diff --git a/src/sp_validation/sacc_io.py b/src/sp_validation/sacc_io.py
index c4bdc5c7..71c552a9 100644
--- a/src/sp_validation/sacc_io.py
+++ b/src/sp_validation/sacc_io.py
@@ -1533,6 +1533,32 @@ def read_nz(path):
     return table[:, 0], table[:, 1:]
 
 
+def read_onecov_matrix(path):
+    """Read a OneCovariance ``style = matrix`` covariance (the workflow's ``.dat``).
+
+    OneCovariance writes the dense matrix with ``np.savetxt`` under a ``#``
+    header; for real-space ξ± its rows are type-major (all ξ+, then all ξ−),
+    which :func:`sp_validation.cosmo_val.sacc_writers.onecov_xi_covariance`
+    reorders to a part's insertion order.
+
+    Parameters
+    ----------
+    path : str or pathlib.Path
+        The matrix file.
+
+    Returns
+    -------
+    numpy.ndarray
+        The square, symmetric covariance matrix.
+    """
+    cov = np.atleast_2d(np.loadtxt(str(path)))
+    if cov.shape[0] != cov.shape[1]:
+        raise ValueError(f"{path}: covariance is not square, shape {cov.shape}")
+    if not np.allclose(cov, cov.T, rtol=1e-6, atol=0):
+        raise ValueError(f"{path}: covariance is not symmetric")
+    return cov
+
+
 def covariance_blocks(cov_list, selectors, *, gaussian=True):
     """Reshape a OneCovariance ``covariance_list`` table into SACC cov blocks.
 
diff --git a/src/sp_validation/tests/regression/test_tomo_plot_cosebis_ignores_tomography.py b/src/sp_validation/tests/regression/test_tomo_plot_cosebis_ignores_tomography.py
index 645a0bb6..fcf1bea2 100644
--- a/src/sp_validation/tests/regression/test_tomo_plot_cosebis_ignores_tomography.py
+++ b/src/sp_validation/tests/regression/test_tomo_plot_cosebis_ignores_tomography.py
@@ -45,12 +45,17 @@ def _make_cv(tmp_path):
     fits.BinTableHDU(t).writeto(cat)
     out = tmp_path / "out"
     out.mkdir()
+    # One n(z) column per bin: the tomographic ξ± part gives each bin a tracer.
+    z = np.linspace(0.05, 2.0, 20)
+    nz_path = tmp_path / "nz_tomo.txt"
+    np.savetxt(nz_path, np.column_stack([z, np.exp(-((z - 0.5) ** 2)), z * 0 + 1]))
     cfg = {
         VER: {
             "blind": "none",
             "subdir": str(tmp_path),
             "shear": {
                 "path": str(cat),
+                "redshift_path": str(nz_path),
                 "e1_col": "e1",
                 "e2_col": "e2",
                 "w_col": "w",
diff --git a/src/sp_validation/tests/test_cosmo_val.py b/src/sp_validation/tests/test_cosmo_val.py
index 9d4ca719..ecc70cb1 100644
--- a/src/sp_validation/tests/test_cosmo_val.py
+++ b/src/sp_validation/tests/test_cosmo_val.py
@@ -592,6 +592,7 @@ def test_cross_tomographic_pair_shares_patch_centres(self, tmp_path, monkeypatch
         params, version = self._write_synthetic_catalogs(
             tmp_path, n_gal=400, with_tomography=True
         )
+        self._use_tomographic_nz(params, version, tmp_path)
         cv = CosmologyValidation(versions=[version], npatch=4, **params)
         cols = cv._shear_columns(version, compute_tomography=True)
         full_catalog = treecorr.Catalog(
@@ -629,9 +630,96 @@ def recording_bin_catalog(cols, bin_id, npatch, patch_centers=None):
             cross_pair_catalogs[1][2], expected_centers, rtol=0, atol=1e-14
         )
 
-    def test_a_blinded_catalogue_refuses_tomographic_signal(self, tmp_path):
-        """Blinding shifts the ("all", "all") parts only, so a blinded
-        catalogue's tomographic ξ± and pseudo-Cℓ are refused before measuring."""
+    @staticmethod
+    def _use_tomographic_nz(params, version, tmp_path, nbin=2):
+        """Point ``version`` at an n(z) file with one column per tomographic bin."""
+        z = np.linspace(0.05, 3.0, 31)
+        columns = [np.exp(-(((z - 0.4 - 0.5 * i) / 0.3) ** 2)) for i in range(nbin)]
+        nz_path = tmp_path / "nz_tomo.txt"
+        np.savetxt(nz_path, np.column_stack([z, *columns]))
+        config = yaml.safe_load(open(params["catalog_config"]))
+        config[version]["shear"]["redshift_path"] = str(nz_path)
+        with open(params["catalog_config"], "w") as f:
+            yaml.dump(config, f, sort_keys=False)
+        return z, columns
+
+    def test_tomographic_xi_is_one_part_with_a_tracer_per_bin(self, tmp_path):
+        """Every bin pair goes into one sealed part: source_i carries column i + 1
+        of the n(z) file, pairs are pair-major in OneCovariance's order, and the
+        jackknife covariance spans the pairs."""
+        params, version = self._write_synthetic_catalogs(
+            tmp_path, n_gal=600, with_tomography=True
+        )
+        z, columns = self._use_tomographic_nz(params, version, tmp_path)
+        cv = CosmologyValidation(versions=[version], npatch=4, **params)
+
+        ggs = cv.calculate_2pcf_version(version, compute_tomography=True)
+
+        part = cv.xi_tomo_parts[version, "reporting"]
+        assert sacc_io.stamp(part) == "none"
+        assert sorted(ggs) == [
+            "tomo_bin_1_tomo_bin_1",
+            "tomo_bin_1_tomo_bin_2",
+            "tomo_bin_2_tomo_bin_2",
+        ]
+        for i, column in enumerate(columns):
+            z_read, nz_read = sacc_io.get_nz(part, i)
+            np.testing.assert_allclose(z_read, z)
+            np.testing.assert_allclose(nz_read, column)
+        nbins = cv.treecorr_config["nbins"]
+        pairs = [dp.tracers for dp in part.data[:: 2 * nbins]]
+        assert pairs == [
+            ("source_0", "source_0"),
+            ("source_0", "source_1"),
+            ("source_1", "source_1"),
+        ]
+        np.testing.assert_allclose(
+            ggs["tomo_bin_1_tomo_bin_2"].xip,
+            sacc_io.get_xi(part, (0, 1), grid="reporting")[1],
+        )
+        cov = part.covariance.dense
+        assert cov.shape == (6 * nbins, 6 * nbins)
+        assert np.any(cov[: 2 * nbins, 2 * nbins : 4 * nbins] != 0)
+
+    def test_tomographic_xi_needs_an_nz_per_bin(self, tmp_path):
+        """A single-column n(z) cannot describe two bins; the part is refused."""
+        params, version = self._write_synthetic_catalogs(
+            tmp_path, n_gal=400, with_tomography=True
+        )
+        cv = CosmologyValidation(versions=[version], npatch=1, **params)
+        with pytest.raises(ValueError, match="n\\(z\\)"):
+            cv.calculate_2pcf_version(version, compute_tomography=True)
+
+    def test_a_blinded_catalogue_conceals_tomographic_xi(self, tmp_path):
+        """On a blinded catalogue every tomographic pair is sealed under the
+        blind, so the ξ± handed back differs from the public measurement."""
+        pytest.importorskip("smokescreen")
+        from sp_validation import blinding
+
+        params, version = self._write_synthetic_catalogs(
+            tmp_path, n_gal=400, with_tomography=True
+        )
+        self._use_tomographic_nz(params, version, tmp_path)
+        public = CosmologyValidation(
+            versions=[version], npatch=1, **params
+        ).calculate_2pcf_version(version, compute_tomography=True)
+
+        config = yaml.safe_load(open(params["catalog_config"]))
+        config["paths"]["blinds"] = str(tmp_path / "blinds")
+        config[version]["blind"] = "toy"
+        with open(params["catalog_config"], "w") as f:
+            yaml.dump(config, f, sort_keys=False)
+        blinding.init("toy", config)
+        cv = CosmologyValidation(versions=[version], npatch=1, **params)
+        blinded = cv.calculate_2pcf_version(version, compute_tomography=True)
+
+        assert sacc_io.stamp(cv.xi_tomo_parts[version, "reporting"]) == "toy"
+        for key in public:
+            assert not np.allclose(blinded[key].xip, public[key].xip)
+
+    def test_a_blinded_catalogue_refuses_tomographic_pseudo_cl(self, tmp_path):
+        """Blinding covers the tomographic ξ± parts but not the pseudo-Cℓ FITS
+        tables, so a blinded catalogue's tomographic pseudo-Cℓ are refused."""
         from sp_validation import blinding
 
         params, version = self._write_synthetic_catalogs(
@@ -645,8 +733,6 @@ def test_a_blinded_catalogue_refuses_tomographic_signal(self, tmp_path):
         blinding.init("toy", config)
         cv = CosmologyValidation(versions=[version], npatch=1, **params)
 
-        with pytest.raises(blinding.BlindingError, match="tomographic"):
-            cv.calculate_2pcf_version(version, compute_tomography=True)
         with pytest.raises(blinding.BlindingError, match="tomographic"):
             cv.calculate_pseudo_cl(compute_tomography=True)
         # Direct producers must refuse before opening even a missing catalogue.
@@ -918,6 +1004,11 @@ def test_sacc_nz_is_the_summed_tomographic_nz(self, tmp_path):
         np.testing.assert_allclose(z_read, z)
         np.testing.assert_allclose(nz_read, nz_bins.sum(axis=0))
 
+        tomographic = cv.sacc_nz(version, tomography=True)
+        assert sorted(tomographic) == [0, 1, 2]
+        for i, nz_bin in enumerate(nz_bins):
+            np.testing.assert_allclose(tomographic[i][1], nz_bin)
+
     @pytest.mark.parametrize("pol_factor", [True, False, 0, 2])
     def test_pol_factor_must_be_plus_or_minus_one(self, tmp_path, pol_factor):
         """pol_factor is a ±1 e2 multiplier; a bool would pass `in (-1, 1)`."""
diff --git a/src/sp_validation/tests/test_glass_mock.py b/src/sp_validation/tests/test_glass_mock.py
index 54014403..1ce1f580 100644
--- a/src/sp_validation/tests/test_glass_mock.py
+++ b/src/sp_validation/tests/test_glass_mock.py
@@ -40,7 +40,24 @@
 
 camb = pytest.importorskip("camb")
 
-HAVE_GLASS = importlib.util.find_spec("glass") is not None
+
+def _have(module):
+    """True if ``module`` is importable. False if it or a parent is missing."""
+    try:
+        return importlib.util.find_spec(module) is not None
+    except ModuleNotFoundError:
+        # find_spec imports the parent package, which raises if it is absent.
+        return False
+
+
+HAVE_GLASS = _have("glass")
+
+# The map path needs more than GLASS itself: build_shells and
+# Cosmology_from_camb import ``cosmology.compat.camb``, which ships with the
+# ``cosmology`` API package rather than with GLASS. Gate on it separately so a
+# missing dependency reports as a skip instead of masquerading as the (real,
+# but different) glass/cosmology API incompatibility tracked in the xfail below.
+HAVE_COSMOLOGY_CAMB = _have("cosmology.compat.camb")
 
 REFERENCE = Path(__file__).parent / "data" / "glass_mock_camb_reference.npz"
 
@@ -131,6 +148,10 @@ def test_config_change_breaks_reference():
 
 
 @pytest.mark.skipif(not HAVE_GLASS, reason="GLASS not installed in this image")
+@pytest.mark.skipif(
+    not HAVE_COSMOLOGY_CAMB,
+    reason="cosmology.compat.camb not installed in this image",
+)
 @pytest.mark.xfail(
     reason=(
         "glass_mock builds its cosmology with cosmology.compat.camb, which an image "
diff --git a/src/sp_validation/tests/test_sacc_io.py b/src/sp_validation/tests/test_sacc_io.py
index d2279969..096dc3ce 100644
--- a/src/sp_validation/tests/test_sacc_io.py
+++ b/src/sp_validation/tests/test_sacc_io.py
@@ -1074,3 +1074,18 @@ def test_merge_partial_files():
     assert cl_bb is not None and cl_eb is None
     n, E, B = sio.get_cosebis(merged, (0, 0))
     assert B is None
+
+
+def test_read_onecov_matrix_skips_the_header(tmp_path):
+    """OneCovariance's `style = matrix` file reads back as the dense matrix."""
+    cov = np.array([[2.0, 0.5], [0.5, 3.0]])
+    path = tmp_path / "cov.dat"
+    np.savetxt(path, cov, fmt="%.6e", header="OneCovariance covariance matrix")
+    np.testing.assert_allclose(sio.read_onecov_matrix(path), cov)
+
+
+def test_read_onecov_matrix_refuses_a_non_square_file(tmp_path):
+    path = tmp_path / "cov.dat"
+    np.savetxt(path, np.ones((2, 3)))
+    with pytest.raises(ValueError, match="not square"):
+        sio.read_onecov_matrix(path)
diff --git a/src/sp_validation/tests/test_sacc_writers.py b/src/sp_validation/tests/test_sacc_writers.py
index a3eb6618..92460f1c 100644
--- a/src/sp_validation/tests/test_sacc_writers.py
+++ b/src/sp_validation/tests/test_sacc_writers.py
@@ -400,3 +400,67 @@ def test_xi_part_rejects_two_covariances():
             covariance=_spd(12, 42),
             variances=np.ones(12),
         )
+
+
+# --------------------------------------------------------------------------- #
+# Tomographic ξ± and the OneCovariance covariance
+# --------------------------------------------------------------------------- #
+def _tomo_xi_part(nbin=2, ntheta=4):
+    """A pair-major ξ± part over every pair of ``nbin`` bins, values labelled."""
+    theta = _theta(ntheta)
+    pairs = [(i, j) for i in range(nbin) for j in range(i, nbin)]
+    measured = [
+        (
+            bins,
+            {
+                "theta": theta,
+                "xip": 1e-5 * (10 * k + np.arange(ntheta)),
+                "xim": -1e-5 * (10 * k + np.arange(ntheta)),
+            },
+        )
+        for k, bins in enumerate(pairs)
+    ]
+    nz = {i: _nz(seed=i) for i in range(nbin)}
+    return sw.xi_pairs_to_sacc(nz, META, measured, grid="reporting"), pairs, theta
+
+
+def test_xi_pairs_to_sacc_writes_every_pair_pair_major(tmp_path):
+    """Each pair goes in as [ξ+; ξ−], pairs in the order given, one tracer per bin."""
+    s, pairs, theta = _tomo_xi_part(nbin=3)
+    s2 = _roundtrip(s, tmp_path, "xitomo")
+    assert {sio.source_name(i) for i in range(3)} <= set(s2.tracers)
+    assert len(s2.mean) == len(pairs) * 2 * len(theta)
+    for k, bins in enumerate(pairs):
+        th, p, m = sio.get_xi(s2, bins, grid="reporting")
+        np.testing.assert_array_equal(th, theta)
+        np.testing.assert_allclose(p, 1e-5 * (10 * k + np.arange(len(theta))))
+        block = slice(k * 2 * len(theta), (k + 1) * 2 * len(theta))
+        np.testing.assert_allclose(s2.mean[block], np.concatenate([p, m]))
+
+
+def test_onecov_xi_covariance_reorders_type_major_to_pair_major():
+    """OneCovariance's [all ξ+; all ξ−] matrix lands on the part's pair-major points."""
+    s, pairs, theta = _tomo_xi_part(nbin=2, ntheta=4)
+    part_order = [(t, b, n) for b in pairs for t in (0, 1) for n in range(len(theta))]
+    onecov_order = sorted(part_order)  # type, then pair, then θ
+    expected = _spd(len(part_order), seed=7)
+    position = {label: k for k, label in enumerate(part_order)}
+    take = [position[label] for label in onecov_order]
+    onecov = expected[np.ix_(take, take)]
+
+    np.testing.assert_array_equal(sw.onecov_xi_covariance(s, onecov), expected)
+
+
+def test_onecov_xi_covariance_is_the_identity_for_one_bin():
+    """For the single (0, 0) pair OneCovariance's order already is the part's."""
+    theta = _theta()
+    s = sw.xi_to_sacc({0: _nz()}, META, theta, theta * 0, theta * 0, grid="reporting")
+    cov = _spd(2 * len(theta), seed=3)
+    np.testing.assert_array_equal(sw.onecov_xi_covariance(s, cov), cov)
+
+
+def test_onecov_xi_covariance_refuses_a_mismatched_matrix():
+    """A matrix from another binning or tomography is refused, not truncated."""
+    s, _, _ = _tomo_xi_part(nbin=2, ntheta=4)
+    with pytest.raises(ValueError, match="binning and tomography"):
+        sw.onecov_xi_covariance(s, np.eye(len(s.mean) - 2))
diff --git a/workflow/rules/cosmo_val.smk b/workflow/rules/cosmo_val.smk
index b416eac1..26cefb0d 100644
--- a/workflow/rules/cosmo_val.smk
+++ b/workflow/rules/cosmo_val.smk
@@ -308,6 +308,7 @@ rule cv_additive_bias:
         "../scripts/cv_additive_bias.py"
 
 
+# LG: DEPRACATED since we do not plot unblinded data vectors.
 rule cv_plot_2pcf:
     """n_pairs / xi± overlay across versions, from the reporting parts."""
     input:
@@ -321,7 +322,7 @@ rule cv_plot_2pcf:
     script:
         "../scripts/cv_plot_2pcf.py"
 
-
+# LG: DEPRACATED since we do not plot unblinded data vectors.
 rule cv_ratio_xi_sys_xi:
     """Ratio of PSF systematics (xi_psf_sys) to the cosmic-shear signal (xi+)."""
     input:
diff --git a/workflow/rules/twopoint.smk b/workflow/rules/twopoint.smk
index 5d6450e6..9f627b2d 100644
--- a/workflow/rules/twopoint.smk
+++ b/workflow/rules/twopoint.smk
@@ -17,15 +17,17 @@ rule xi:
     """TreeCorr ξ±(θ) for one version on one angular grid, as its sealed SACC part.
 
     One rule for every grid: outputs are named by their binning, so a request
-    binds the wildcards and the grid label resolves from them.
+    binds the wildcards and the grid label resolves from them. A `_tomo`
+    suffix measures every tomographic bin pair into the one part.
     """
     input:
         catalog=get_shear_catalog,
     output:
-        sacc=str(COSMO_VAL / "{version}_xi_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.sacc"),
+        sacc=str(COSMO_VAL / "{version}_xi_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}{tomo_suffix}.sacc"),
     threads: 24
     params:
         ver="{version}",
+        compute_tomography=lambda w: w.tomo_suffix == "_tomo",
         min_sep="{min_sep}",
         max_sep="{max_sep}",
         nbins="{nbins}",
@@ -44,6 +46,38 @@ rule xi:
         "../scripts/run_2pcf.py"
 
 
+def xi_onecov_covariance(w):
+    """The OneCovariance matrix on a ξ± request's binning and tomography:
+    Gaussian on the integration grid, non-Gaussian otherwise, as cosmo_val's
+    covariance inputs are."""
+    gaussian = "g" if grid_of(XI_GRIDS, w) == "integration" else "ng"
+    return covariance_path(
+        w.version, gaussian, w.min_sep, w.max_sep, w.nbins,
+        DEFAULT_MASK_SUFFIX, w.tomo_suffix,
+    )
+
+
+rule xi_onecov:
+    """A ξ± part with the OneCovariance matrix of its binning as its covariance.
+
+    The matrix replaces the estimate the part was born with (jackknife or shot
+    noise), reordered to the part's points; the result is a derivation of the
+    part, so it keeps the part's blind.
+    """
+    input:
+        xi=str(COSMO_VAL / "{version}_xi_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}{tomo_suffix}.sacc"),
+        cov=xi_onecov_covariance,
+    output:
+        sacc=str(COSMO_VAL / "{version}_xi_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}{tomo_suffix}_onecov.sacc"),
+    params:
+        blind=lambda w: blind_of(w.version),
+    resources:
+        mem_mb=4000,
+        runtime=10,
+    script:
+        "../scripts/xi_onecov_sacc.py"
+
+
 rule rho_tau_stats:
     output:
         rho_stats=str(COSMO_VAL / "rho_tau_stats/rho_stats_{version}_tomo_bin_all_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.fits"),
diff --git a/workflow/scripts/assemble_sacc.py b/workflow/scripts/assemble_sacc.py
index 0139acdd..1f1a0968 100644
--- a/workflow/scripts/assemble_sacc.py
+++ b/workflow/scripts/assemble_sacc.py
@@ -19,7 +19,10 @@
 
 from sp_validation import sacc_io
 from sp_validation.blinding import Blind
-from sp_validation.cosmo_val.sacc_writers import assemble_analysis_sacc
+from sp_validation.cosmo_val.sacc_writers import (
+    assemble_analysis_sacc,
+    onecov_xi_covariance,
+)
 
 # NaMaster iNKA covariance FITS: per-spectrum HDU names, in SACC insertion order.
 _CL_HDUS = ("COVAR_EE_EE", "COVAR_BB_BB", "COVAR_EB_EB")
@@ -74,7 +77,7 @@ def _attach_cov(part, name, xi_cov, pseudo_cl_cov):
             f"--{_INJECTED[name].replace('_', '-')}, which was not supplied"
         )
     block = (
-        np.loadtxt(supplied)
+        onecov_xi_covariance(part, sacc_io.read_onecov_matrix(supplied))
         if name == "xi_reporting"
         else _pseudo_cl_cov_block(supplied)
     )
diff --git a/workflow/scripts/cv_plot_2pcf.py b/workflow/scripts/cv_plot_2pcf.py
index 5bb3ac32..90c09836 100644
--- a/workflow/scripts/cv_plot_2pcf.py
+++ b/workflow/scripts/cv_plot_2pcf.py
@@ -8,6 +8,7 @@
 
 from cv_runner import _unbuffer_streams, make_cv, take_reporting_parts, touch_sentinels
 
+# LG : This is now DEPRECATED, NOT TO BE USED, since we are no longer plotting unblinded data vectors.
 _unbuffer_streams()
 cv = make_cv(snakemake)
 take_reporting_parts(cv, snakemake.input["xi"])
diff --git a/workflow/scripts/cv_ratio_xi_sys_xi.py b/workflow/scripts/cv_ratio_xi_sys_xi.py
index 9e1e4396..9709445d 100644
--- a/workflow/scripts/cv_ratio_xi_sys_xi.py
+++ b/workflow/scripts/cv_ratio_xi_sys_xi.py
@@ -9,6 +9,7 @@
 
 from cv_runner import _unbuffer_streams, make_cv, take_reporting_parts, verify_outputs
 
+# LG : This is now DEPRECATED, NOT TO BE USED, since we are no longer plotting unblinded data vectors.
 _unbuffer_streams()
 cv = make_cv(snakemake)
 take_reporting_parts(cv, snakemake.input["xi"])
diff --git a/workflow/scripts/run_2pcf.py b/workflow/scripts/run_2pcf.py
index 980eaf34..9625c092 100644
--- a/workflow/scripts/run_2pcf.py
+++ b/workflow/scripts/run_2pcf.py
@@ -10,12 +10,14 @@
         --ver SP_v1.4.6.3_leak_corr \
         --min-sep 1.0 --max-sep 250.0 --nbins 20 --npatch 1 \
         --cat-config /path/to/cosmo_val/cat_config.yaml \
-        --out <output_dir>
+        --out <output_dir> [--tomo]
 
 The measurement is binning-agnostic: the reporting and the fine integration
 grids are the same compute with different ``--min-sep/--max-sep/--nbins``. The
-non-tomographic ``("all", "all")`` ξ± is born as a SACC part, named by its
-binning, tagged with its ``--grid`` and sealed under the catalogue's blind by
+ξ± — the non-tomographic ``("all", "all")`` pair, or with ``--tomo`` every
+tomographic bin pair with one n(z) tracer per bin — is born as a SACC part,
+named by its binning (``_tomo`` appended for tomography), tagged with its
+``--grid`` and sealed under the catalogue's blind by
 ``CosmologyValidation.calculate_2pcf_version``. The part carries the covariance
 the measurement estimated: the dense jackknife covariance when it had patches,
 the shot-noise ``varxip``/``varxim`` diagonal when it had none.
@@ -40,6 +42,7 @@ def run_2pcf(
     output_dir,
     sacc_out=None,
     grid="reporting",
+    compute_tomography=False,
 ):
     """Measure ξ±(θ) for ``ver`` and write its sealed SACC part.
 
@@ -50,6 +53,8 @@ def run_2pcf(
     ``cat_config['paths']['output']``. ``sacc_out`` is the exact destination for
     the part (the Snakemake-declared output); it defaults to a binning-derived
     name under the resolved output directory for the CLI path.
+    ``compute_tomography`` measures every tomographic bin pair instead of
+    ``("all", "all")``.
 
     Returns
     -------
@@ -63,18 +68,21 @@ def run_2pcf(
     )
     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",
+        f"{ver}_xi_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}"
+        f"{'_tomo' if compute_tomography else ''}.sacc",
     )
     cv.calculate_2pcf_version(
         ver,
         npatch=npatch,
+        compute_tomography=compute_tomography,
         grid=grid,
         out=out_path,
         min_sep=min_sep,
         max_sep=max_sep,
         nbins=nbins,
     )
-    part = cv.xi_parts[ver, grid]
+    parts = cv.xi_tomo_parts if compute_tomography else cv.xi_parts
+    part = parts[ver, grid]
     print(f"Wrote {grid} ξ± SACC part: {out_path}")
     return part
 
@@ -91,6 +99,7 @@ def _from_snakemake(smk):
         output_dir=p["output_dir"],
         grid=p.get("grid", "reporting"),
         sacc_out=smk.output["sacc"],
+        compute_tomography=bool(p.get("compute_tomography", False)),
     )
 
 
@@ -120,6 +129,11 @@ def _from_cli(argv=None):
     ap.add_argument(
         "--grid", default="reporting", help="SACC grid tag for the measured points"
     )
+    ap.add_argument(
+        "--tomo",
+        action="store_true",
+        help="Measure every tomographic bin pair instead of the whole sample",
+    )
     a = ap.parse_args(argv)
     run_2pcf(
         ver=a.ver,
@@ -130,6 +144,7 @@ def _from_cli(argv=None):
         cat_config=a.cat_config,
         output_dir=a.out,
         grid=a.grid,
+        compute_tomography=a.tomo,
     )
 
 
diff --git a/workflow/scripts/xi_onecov_sacc.py b/workflow/scripts/xi_onecov_sacc.py
new file mode 100644
index 00000000..878ebed4
--- /dev/null
+++ b/workflow/scripts/xi_onecov_sacc.py
@@ -0,0 +1,74 @@
+"""Give a ξ± SACC part the OneCovariance matrix of its binning as its covariance.
+
+Dual-mode. Under Snakemake (``script:`` directive) the injected ``snakemake``
+object supplies the paths; as a standalone CLI the same step runs from flags:
+
+    python xi_onecov_sacc.py \
+        --xi <version>_xi_..._npatch=1_tomo.sacc \
+        --cov <covariance_dir>/<base>.dat \
+        --blind none \
+        --out <version>_xi_..._npatch=1_tomo_onecov.sacc
+
+OneCovariance writes ξ± type-major; the part is pair-major, so the matrix is
+reordered to the part's points
+(:func:`sp_validation.cosmo_val.sacc_writers.onecov_xi_covariance`) before it
+replaces the part's own estimate. The output is a derivation of the part and
+carries its blind stamp, which must be the catalogue's declared blind.
+"""
+
+import argparse
+
+from sp_validation import sacc_io
+from sp_validation.blinding import Blind
+from sp_validation.cosmo_val.sacc_writers import onecov_xi_covariance
+
+
+def xi_onecov_sacc(xi_path, cov_path, out_path, *, blind):
+    """Write the part at ``xi_path`` with the matrix at ``cov_path`` to ``out_path``.
+
+    Returns
+    -------
+    sacc.Sacc
+        The data set as written.
+    """
+    part = sacc_io.load(xi_path, blind=blind)
+    s = part.copy()
+    s.add_covariance(
+        onecov_xi_covariance(part, sacc_io.read_onecov_matrix(cov_path)),
+        overwrite=True,
+    )
+    out = sacc_io.save(s, out_path, derived_from=[part], blind=Blind(blind))
+    print(f"Wrote ξ± with the OneCovariance covariance {cov_path}: {out_path}")
+    return out
+
+
+def _from_snakemake(smk):
+    xi_onecov_sacc(
+        smk.input["xi"],
+        smk.input["cov"],
+        smk.output["sacc"],
+        blind=smk.params["blind"],
+    )
+
+
+def _from_cli(argv=None):
+    ap = argparse.ArgumentParser(
+        description="Attach a OneCovariance ξ± matrix to a ξ± SACC part."
+    )
+    ap.add_argument("--xi", required=True, help="ξ± SACC part")
+    ap.add_argument("--cov", required=True, help="OneCovariance matrix (.dat)")
+    ap.add_argument(
+        "--blind", required=True, help="The blind the catalogue is declared under"
+    )
+    ap.add_argument("--out", required=True, help="Output SACC path")
+    a = ap.parse_args(argv)
+    xi_onecov_sacc(a.xi, a.cov, a.out, blind=a.blind)
+
+
+if __name__ == "__main__":
+    try:
+        snakemake  # noqa: F821 — injected by Snakemake's script: directive
+    except NameError:
+        _from_cli()
+    else:
+        _from_snakemake(snakemake)  # noqa: F821