diff --git a/docs/pair-count-noise.md b/docs/pair-count-noise.md
new file mode 100644
index 00000000..e52b07f2
--- /dev/null
+++ b/docs/pair-count-noise.md
@@ -0,0 +1,75 @@
+# Pair-count shape noise
+
+`cosmo_val` workflows (whole-sample and tomographic) use measured effective pairs for OneCovariance's real-space Gaussian shape-noise term by default.
+Set `covariance.shape_noise: analytic` to recover the uniform-area approximation.
+Other workflows retain their analytic default.
+This changes neither blinding nor the mixed, sample-variance, non-Gaussian or SSC terms.
+
+## Cache and normalization
+
+`covariance_pair_counts` reads only the configured RA, Dec, weight and (for tomography) bin-ID columns of the same selected shear catalogue as `xi`.
+It counts on each requested covariance grid, using TreeCorr NN correlations with weights `w` and `w²`, spherical coordinates, logarithmic bins and `min_top=6`, matching `xi`.
+The cache is independent of cosmology and contains no shapes, signals or blind records.
+The reader projects fields directly from the lazy catalogue chunks, including compound HDF5 datasets, rather than materialising rows containing shape fields.
+It is a separate rule so covariance-only reruns need no shape measurement or access to a blind.
+Its `threads`, memory and runtime can be overridden just like `xi` (large catalogues may need more than the default 40 GB).
+For the full v2 catalogue, use `--set-resources covariance_pair_counts:mem_mb=128000 covariance_pair_counts:runtime=1440`, alongside the run's existing resource overrides; the two NN passes can peak above the ξ rule's memory usage.
+
+For unordered auto pairs in bin b, the NPZ stores:
+
+- `weight`: W = Σ w_i w_j.
+- `weight_squared`: Q = Σ (w_i w_j)².
+- `npairs`: the raw unordered count, plus exact bin edges and nominal centres.
+- `ordered_effective_pairs`: 2 W²/Q.
+
+With independent isotropic shapes of constant per-component dispersion s, the pure-noise variance of either ξ component is `2 s⁴ Q/W²`.
+The total two-component RMS is σ_e = √2 s, so the same variance is `σ_e⁴/(2 N_eff)`, where `N_eff = W²/Q` is **unordered**.
+Mixing the two dispersion conventions gives a factor-of-four error.
+
+## OneCovariance interface
+
+Verified against OneCovariance `c8755d6de43aadaa1cd4fe132ecf1c1d534e61f8`, the source used by the v2 run with the `ddd58fe0` container.
+The external checkout is selected by `tools.onecov_executable`; the image supplies its dependencies, not the source.
+`onecov/cov_input.py` reads `[tabulated inputs files]` keys `npair_directory` and `npair_mm_file` (first column θ in arcmin, last column counts).
+`onecov/cov_setup.py:get_npair` accepts a matching grid without rebinning and otherwise falls back to an unsafe-for-effective-counts sum.
+The workflow therefore writes **nominal geometric bin centres**, not measured mean separations, and counts separately on every covariance grid.
+Empty or non-finite pair bins fail instead of silently restoring analytic noise.
+
+OneCovariance's analytic auto count is `A_arcmin² n_eff² π(θ_hi²−θ_lo²)`, which counts **ordered** pairs.
+At `onecov/cov_theta_space.py` lines 1288–1297, the two Kronecker contractions and the `/0.5` yield `4 s⁴/N_ordered` for a single source bin.
+The two-column `.txt` thus supplies `N_ordered = 2 W²/Q`, giving exactly `2 s⁴ Q/W²`.
+No area rescaling or post-run whole-covariance scaling is applied.
+
+The normalization test checks the uniform-field expectation analytically, including non-unit equal weights and both dispersion conventions.
+A four-object positions-only catalogue checks TreeCorr's unordered convention against explicit pair enumeration.
+HDF5 tests forbid whole-dataset or shape-column reads for single, chunked and joined layouts.
+DAG tests verify the cache-to-ini dependency, analytic opt-out, dispersion boundary and tomographic file order and footprint inputs.
+
+## Tomographic pairs
+
+Each positive tomographic bin pair is counted separately, including cross correlations between disjoint bins.
+`shear.tomo_bin_col` names the bin-ID field; zero and negative IDs are excluded, matching ξ± selection.
+One-based indices in filenames refer to the sorted positive catalogue IDs, in the same order as the n(z) columns and `cov_th.tomo` survey parameters.
+The workflow writes `{version}_pairs_minsep=..._maxsep=..._nbins=..._tomo_bin=i-j.npz` and `.txt`, ordered `(1,1), (1,2), …, (2,2), …` in `npair_mm_file`.
+Each cache records the bin indices and the actual catalogue IDs.
+The area's and mask's source remains `cov_th.tomo` as for analytic noise.
+Geometry references must also have identical bin IDs when tomography is enabled.
+
+For an auto pair the table contains **2W²/Q**, with `ordered_effective_pairs` in the cache.
+For a cross pair it contains **W²/Q**, with `effective_cross_pairs` in the cache: TreeCorr counts each cross pair once, and only one of OneCovariance's two Kronecker contractions survives.
+Both yield the independent-shape variance `2 s_i² s_j² Q/W²`.
+Passing 2W²/Q for a cross would underestimate its noise by a factor of two.
+Pure shape noise has no covariance between different unordered bin pairs or between ξ+ and ξ−.
+
+The opt-in test `test_pair_counts_onecov.py` executes the external checkout's actual file reader and real-space noise kernel with unequal bin dispersions and explicit weighted-pair answers.
+Only the Fourier cosmological terms are stubbed to zero; it does not run the full cosmological integration.
+Run it inside the dependency container with `ONECOVARIANCE_PATH` set to the workflow's external source checkout.
+The tomographic probe is verified against candide checkout `70265be3a4bffcfded7b9251a3c4cbe8159d793f`.
+
+## Existing products
+
+Existing covariance files do not encode the noise prescription in their filenames.
+When migrating a completed run, use a fresh output tree or explicitly rerun `covariance_ini` and its descendants; do not interpret an old cached matrix as measured-pair noise.
+Do not sum fine-grid effective counts to obtain coarse-grid counts: sum W and Q separately before forming W²/Q, or count directly on that grid as this rule does.
+Pair noise includes the observed spatial selection and weight variation; it is not an independent mask-only simulation.
+Mask-dependent mixed and sample-variance corrections remain separate validation work.
diff --git a/papers/cosmo_val/config/config.yaml b/papers/cosmo_val/config/config.yaml
index 99f553c9..07e2b3ad 100644
--- a/papers/cosmo_val/config/config.yaml
+++ b/papers/cosmo_val/config/config.yaml
@@ -114,6 +114,9 @@ var_method: jackknife
 
 covariance:
   default_masked: true
+  # Real-space shape noise, including tomography: measured W²/Q includes pair loss.
+  # Set analytic to restore OneCovariance's uniform-area approximation.
+  shape_noise: pair_counts
 
 cl:
   n_ell_bins: 32
diff --git a/src/sp_validation/pair_counts.py b/src/sp_validation/pair_counts.py
new file mode 100644
index 00000000..4f84d48e
--- /dev/null
+++ b/src/sp_validation/pair_counts.py
@@ -0,0 +1,147 @@
+"""Non-cosmological weighted pair counts for real-space shape noise.
+
+TreeCorr auto counts are unordered. OneCovariance's xi_pm noise is
+4*s**4/N_ordered, where s is the per-component dispersion, so supply
+N_ordered = 2*W**2/Q for autos. Cross bins use W**2/Q and one
+Kronecker contraction. No shape columns are accessed.
+"""
+
+from pathlib import Path
+
+import numpy as np
+import treecorr
+
+from sp_validation.io import Catalogue
+
+
+def ordered_effective_pairs(weight, weight_squared):
+    """Convert unordered W=sum(w_i*w_j), Q=sum((w_i*w_j)**2) to 2 W²/Q.
+
+    Empty/non-finite bins cannot define a covariance; fail rather than
+    silently falling back to analytic counts.
+    """
+    weight = np.asarray(weight, dtype=float)
+    weight_squared = np.asarray(weight_squared, dtype=float)
+    if (
+        weight.shape != weight_squared.shape
+        or not np.all(np.isfinite(weight))
+        or not np.all(np.isfinite(weight_squared))
+        or np.any(weight <= 0)
+        or np.any(weight_squared <= 0)
+    ):
+        raise ValueError("Pair counts require matching, finite, positive W and Q")
+    return 2 * weight**2 / weight_squared
+
+
+def measure_pair_counts(
+    catalogue,
+    output,
+    *,
+    min_sep,
+    max_sep,
+    nbins,
+    ra_col="RA",
+    dec_col="Dec",
+    w_col="w",
+    num_threads=1,
+    tomo_bin_col=None,
+    bin_pair=None,
+):
+    """Cache W, Q and raw counts using positions and weights only.
+
+    Uses the xi rule's log binning, spherical metric and min_top=6.
+    The adjacent two-column .txt has nominal log-bin centres in arcmin,
+    never mean pair separations: OneCovariance must not rebin W²/Q.
+    ``bin_pair`` contains one-based indices into the sorted positive IDs in
+    ``tomo_bin_col``, matching the n(z) column order used by OneCovariance.
+    """
+    # Catalogue.read materialises whole compound HDF5 chunks before selecting
+    # columns. Index the public lazy chunk tables field-by-field instead, so
+    # shapes are never loaded, for either FITS or any supported HDF5 layout.
+    columns = (ra_col, dec_col, w_col) + ((tomo_bin_col,) if bin_pair else ())
+    with Catalogue(str(catalogue)) as cat:
+        data = np.empty(len(cat), dtype=cat.dtype(columns))
+        start = 0
+        for table in cat.chunks.values():
+            end = start + len(table)
+            for name in columns:
+                data[name][start:end] = table[name]
+            start = end
+    if bin_pair:
+        ids = np.unique(data[tomo_bin_col])
+        ids = ids[ids > 0]
+        # OneCovariance indexes the sorted positive bins from one, as n(z) does.
+        selected_ids = ids[np.asarray(bin_pair) - 1]
+        selections = [data[tomo_bin_col] == value for value in selected_ids]
+    else:
+        selections = [slice(None)] * 2
+    auto = bin_pair is None or bin_pair[0] == bin_pair[1]
+    weights = np.asarray(data[w_col], dtype=float)
+    if (
+        not np.all(np.isfinite(weights))
+        or np.any(weights < 0)
+        or not np.all(np.isfinite(data[ra_col]))
+        or not np.all(np.isfinite(data[dec_col]))
+    ):
+        raise ValueError("Positions and non-negative weights must be finite")
+    config = dict(
+        min_sep=float(min_sep),
+        max_sep=float(max_sep),
+        nbins=int(nbins),
+        sep_units="arcmin",
+        min_top=6,
+        num_threads=int(num_threads),
+    )
+    constant_weight = (
+        float(weights[0]) if len(weights) and np.all(weights == weights[0]) else None
+    )
+    counts = []
+    for power in (1, 2):
+        cats = [
+            treecorr.Catalog(
+                ra=data[ra_col][mask],
+                dec=data[dec_col][mask],
+                w=weights[mask] ** power,
+                ra_units="deg",
+                dec_units="deg",
+            )
+            for mask in selections[: 1 if auto else 2]
+        ]
+        nn = treecorr.NNCorrelation(config)
+        nn.process(*cats)
+        counts.append(nn.weight.copy())
+        print(f"Completed NN weights^{power}", flush=True)
+        if constant_weight is not None:
+            # For w_i=c, W=c² N and Q=c⁴ N=c² W on the same geometry.
+            # The second traversal carries no additional information.
+            counts.append(counts[0] * constant_weight**2)
+            print("Derived squared-weight counts for constant weights", flush=True)
+            break
+    W, Q = counts
+    # One contraction survives for disjoint cross bins, two for an auto.
+    # OneCovariance therefore needs W²/Q for crosses and 2 W²/Q for autos.
+    effective = ordered_effective_pairs(W, Q) / (1 if auto else 2)
+    edges = np.geomspace(float(min_sep), float(max_sep), int(nbins) + 1)
+    theta = np.sqrt(edges[:-1] * edges[1:])
+    output = Path(output)
+    output.parent.mkdir(parents=True, exist_ok=True)
+    np.savez(
+        output,
+        **(
+            {"bin_pair": np.asarray(bin_pair), "bin_ids": selected_ids}
+            if bin_pair
+            else {}
+        ),
+        theta=theta,
+        edges=edges,
+        npairs=nn.npairs,
+        weight=W,
+        weight_squared=Q,
+        **{("ordered_effective_pairs" if auto else "effective_cross_pairs"): effective},
+    )
+    np.savetxt(
+        output.with_suffix(".txt"),
+        np.column_stack((theta, effective)),
+        header="theta_arcmin "
+        + ("N_ordered_effective" if auto else "N_cross_effective"),
+    )
diff --git a/src/sp_validation/tests/test_pair_counts.py b/src/sp_validation/tests/test_pair_counts.py
new file mode 100644
index 00000000..96776295
--- /dev/null
+++ b/src/sp_validation/tests/test_pair_counts.py
@@ -0,0 +1,194 @@
+"""Weighted-pair normalization, selection, and shape-free I/O."""
+
+import h5py
+import numpy as np
+import pytest
+from astropy.io import fits
+
+from sp_validation.pair_counts import measure_pair_counts, ordered_effective_pairs
+
+
+def test_uniform_limit_and_dispersion_convention():
+    # Exact expectation for a uniform Poisson field, far from boundaries.
+    # All unordered pairs have the same product weight, which must cancel.
+    edges = np.geomspace(1, 83, 31)
+    area, density, s, pair_weight = 2000 * 3600, 5.0, 0.33, 7.0**2
+    unordered = area * density**2 * np.pi * np.diff(edges**2) / 2
+    W = unordered * pair_weight
+    Q = unordered * pair_weight**2
+    ordered = ordered_effective_pairs(W, Q)
+    np.testing.assert_allclose(ordered, 2 * unordered)
+    np.testing.assert_allclose(4 * s**4 / ordered, 2 * s**4 * Q / W**2)
+    sigma_total = np.sqrt(2) * s
+    np.testing.assert_allclose(sigma_total**4 / (2 * unordered), 4 * s**4 / ordered)
+
+
+@pytest.mark.parametrize(
+    "W,Q",
+    [
+        ([0], [1]),
+        ([1], [0]),
+        ([1], [-1]),
+        ([np.nan], [1]),
+        ([1], [np.inf]),
+        ([1, 2], [1]),
+    ],
+)
+def test_invalid_counts_fail_closed(W, Q):
+    with pytest.raises(ValueError):
+        ordered_effective_pairs(W, Q)
+
+
+def test_positions_weights_only_and_unordered_treecorr(tmp_path):
+    # There are deliberately no shape or blind columns/records in this FITS.
+    ra = np.array([0.0, 0.01, 0.03, 0.07])
+    weights = np.array([1.0, 2.0, 3.0, 4.0])
+    cat = tmp_path / "positions.fits"
+    fits.BinTableHDU.from_columns(
+        [
+            fits.Column(name="RA", format="D", array=ra),
+            fits.Column(name="Dec", format="D", array=np.zeros(4)),
+            fits.Column(name="w", format="D", array=weights),
+        ]
+    ).writeto(cat)
+    out = tmp_path / "pairs.npz"
+    measure_pair_counts(cat, out, min_sep=0.1, max_sep=10.0, nbins=1)
+    products = np.array(
+        [weights[i] * weights[j] for i in range(4) for j in range(i + 1, 4)]
+    )
+    with np.load(out) as z:
+        np.testing.assert_allclose(z["weight"], products.sum())
+        np.testing.assert_allclose(z["weight_squared"], (products**2).sum())
+        np.testing.assert_allclose(z["npairs"], 6)
+        np.testing.assert_allclose(z["theta"], [1.0])
+        np.testing.assert_allclose(
+            np.loadtxt(out.with_suffix(".txt"))[1],
+            2 * products.sum() ** 2 / (products**2).sum(),
+        )
+
+
+@pytest.mark.parametrize("weight", [1.0, 2.0, 0.7])
+def test_constant_weights_use_one_tree_traversal(tmp_path, monkeypatch, weight):
+    import treecorr
+
+    calls = []
+    process = treecorr.NNCorrelation.process
+
+    def record(self, *args, **kwargs):
+        calls.append(1)
+        return process(self, *args, **kwargs)
+
+    monkeypatch.setattr(treecorr.NNCorrelation, "process", record)
+    cat = tmp_path / "constant.fits"
+    fits.BinTableHDU.from_columns(
+        [
+            fits.Column(name="RA", format="D", array=[0.0, 0.01, 0.03, 0.07]),
+            fits.Column(name="Dec", format="D", array=np.zeros(4)),
+            fits.Column(name="w", format="D", array=np.full(4, weight)),
+        ]
+    ).writeto(cat)
+    out = tmp_path / "pairs.npz"
+    measure_pair_counts(cat, out, min_sep=0.1, max_sep=10.0, nbins=1)
+    assert len(calls) == 1
+    with np.load(out) as z:
+        np.testing.assert_allclose(z["weight"], 6 * weight**2)
+        np.testing.assert_allclose(z["weight_squared"], 6 * weight**4)
+        np.testing.assert_allclose(z["ordered_effective_pairs"], 12)
+
+
+@pytest.mark.parametrize("layout", ["single", "chunked", "joined"])
+@pytest.mark.parametrize("tomographic", [False, True])
+def test_hdf5_never_materialises_shapes(tmp_path, monkeypatch, layout, tomographic):
+    cat = tmp_path / "catalogue.fits"  # Detection is by content, not extension.
+    data = np.zeros(4, dtype=[(n, "f8") for n in ("RA", "Dec", "w", "e1", "e2", "bin")])
+    data["RA"] = [0.0, 0.01, 0.03, 0.07]
+    data["w"] = [1.0, 2.0, 3.0, 4.0]
+    data["bin"] = [1, 1, 2, 2]
+    allowed = {"RA", "Dec", "w"} | ({"bin"} if tomographic else set())
+    with h5py.File(cat, "w") as f:
+        if layout == "single":
+            f.create_dataset("data", data=data)
+        elif layout == "chunked":
+            group = f.create_group("tiles")
+            group.create_dataset("0001", data=data[:2])
+            group.create_dataset("0002", data=data[2:])
+            f.attrs["n_tiles"] = 2
+        else:
+            f.create_dataset("data", data=data[["RA", "Dec", "e1", "e2", "bin"]])
+            f.create_dataset("data_ext", data=data[["w"]])
+    original = h5py.Dataset.__getitem__
+    reads = []
+
+    def only_pair_fields(dataset, key, *args, **kwargs):
+        # Grammar may inspect zero-row slices to discover field dtypes.
+        if isinstance(key, slice) and key.start == key.stop == 0:
+            return original(dataset, key, *args, **kwargs)
+        # h5py field access recurses with a projected new_dtype and key ().
+        projected = kwargs.get("new_dtype")
+        if projected is not None:
+            assert set(projected.names) <= allowed
+            return original(dataset, key, *args, **kwargs)
+        # Unprojected whole-dataset or shape-field reads remain forbidden.
+        assert isinstance(key, str) and key in allowed, key
+        reads.append(key)
+        return original(dataset, key, *args, **kwargs)
+
+    monkeypatch.setattr(h5py.Dataset, "__getitem__", only_pair_fields)
+    out = tmp_path / "pairs.npz"
+    measure_pair_counts(
+        cat,
+        out,
+        min_sep=0.1,
+        max_sep=10.0,
+        nbins=1,
+        **({"tomo_bin_col": "bin", "bin_pair": (1, 2)} if tomographic else {}),
+    )
+    assert set(reads) == allowed
+    with np.load(out) as z:
+        np.testing.assert_allclose(z["weight"], 21.0 if tomographic else 35.0)
+        np.testing.assert_allclose(z["weight_squared"], 125.0 if tomographic else 273.0)
+
+
+@pytest.mark.parametrize(
+    "pair, indices, factor",
+    [
+        ((1, 1), [(1, 2)], 2),
+        ((1, 2), [(1, 3), (1, 4), (2, 3), (2, 4)], 1),
+        ((2, 2), [(3, 4)], 2),
+    ],
+)
+def test_tomographic_pairs_known_answer(tmp_path, pair, indices, factor):
+    # Negative and zero IDs are excluded; sorted positive IDs needn't be consecutive.
+    weights = np.array([100.0, 1.0, 2.0, 3.0, 4.0, 100.0])
+    cat = tmp_path / "tomo.fits"
+    fits.BinTableHDU.from_columns(
+        [
+            fits.Column(name="RA", format="D", array=np.arange(6) * 0.01),
+            fits.Column(name="Dec", format="D", array=np.zeros(6)),
+            fits.Column(name="w", format="D", array=weights),
+            fits.Column(name="TOM_BIN_ID", format="J", array=[0, 2, 2, 4, 4, -1]),
+        ]
+    ).writeto(cat)
+    output = tmp_path / "counts.npz"
+    measure_pair_counts(
+        cat,
+        output,
+        min_sep=0.1,
+        max_sep=10,
+        nbins=1,
+        tomo_bin_col="TOM_BIN_ID",
+        bin_pair=pair,
+    )
+    products = np.array([weights[i] * weights[j] for i, j in indices])
+    with np.load(output) as z:
+        np.testing.assert_allclose(z["weight"], products.sum())
+        np.testing.assert_allclose(z["weight_squared"], (products**2).sum())
+        np.testing.assert_allclose(z["npairs"], len(indices))
+        np.testing.assert_allclose(
+            z[
+                "ordered_effective_pairs"
+                if pair[0] == pair[1]
+                else "effective_cross_pairs"
+            ],
+            factor * products.sum() ** 2 / (products**2).sum(),
+        )
diff --git a/src/sp_validation/tests/test_pair_counts_onecov.py b/src/sp_validation/tests/test_pair_counts_onecov.py
new file mode 100644
index 00000000..1ea6a427
--- /dev/null
+++ b/src/sp_validation/tests/test_pair_counts_onecov.py
@@ -0,0 +1,102 @@
+"""Opt-in probe of the external OneCovariance table reader and noise kernel.
+
+Run with ONECOVARIANCE_PATH pointing to the checkout used by the workflow.
+Only the cosmological Fourier terms are stubbed; the real table reader and
+real-space noise contractions execute from that checkout.
+"""
+
+import os
+
+import numpy as np
+import pytest
+from astropy.io import fits
+
+from sp_validation.pair_counts import measure_pair_counts
+
+
+@pytest.mark.skipif(
+    not os.environ.get("ONECOVARIANCE_PATH"),
+    reason="Set ONECOVARIANCE_PATH to the workflow's external checkout",
+)
+def test_onecov_tomographic_noise_kernel(tmp_path, monkeypatch):
+    monkeypatch.syspath_prepend(os.environ["ONECOVARIANCE_PATH"])
+    from onecov.cov_input import FileInput
+    from onecov.cov_theta_space import CovTHETASpace
+
+    cat = tmp_path / "positions.fits"
+    weights = np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0])
+    ids = np.array([1, 1, 1, 2, 2, 2])
+    fits.BinTableHDU.from_columns(
+        [
+            fits.Column(name="RA", format="D", array=np.arange(6) * 0.01),
+            fits.Column(name="Dec", format="D", array=np.zeros(6)),
+            fits.Column(name="w", format="D", array=weights),
+            fits.Column(name="bin", format="J", array=ids),
+        ]
+    ).writeto(cat)
+    pairs = [(1, 1), (1, 2), (2, 2)]
+    files, expected = [], {}
+    s = np.array([0.23, 0.37])
+    for i, j in pairs:
+        out = tmp_path / f"pairs_{i}_{j}.npz"
+        measure_pair_counts(
+            cat,
+            out,
+            min_sep=0.1,
+            max_sep=10,
+            nbins=1,
+            tomo_bin_col="bin",
+            bin_pair=(i, j),
+        )
+        files.append(out.with_suffix(".txt").name)
+        products = np.array(
+            [
+                weights[a] * weights[b]
+                for a in range(6)
+                for b in range(6)
+                if ids[a] == i and ids[b] == j and (i != j or a < b)
+            ]
+        )
+        expected[i - 1, j - 1] = (
+            2
+            * s[i - 1] ** 2
+            * s[j - 1] ** 2
+            * (products**2).sum()
+            / products.sum() ** 2
+        )
+
+    reader = FileInput({})
+    reader.clustering, reader.ggl, reader.cosmicshear = False, False, True
+    reader.n_tomo_lens = 2
+    reader._FileInput__get_npair_tabs(
+        {
+            "tabulated inputs files": {
+                "npair_directory": str(tmp_path),
+                "npair_mm_file": ", ".join(files),
+            }
+        }
+    )
+    assert reader.npair_mm.shape == (1, 2, 2, 1)
+
+    cov = CovTHETASpace.__new__(CovTHETASpace)
+    cov.cov_dict = {"gauss": True, "split_gauss": True, "sn_only": True}
+    cov.csmf, cov.gg, cov.gm, cov.mm = False, False, False, True
+    cov.n_tomo_lens = 2
+    cov.sample_dim = 1
+    cov.xi_pp, cov.xi_pm, cov.xi_mm = True, True, True
+    cov.theta_bins_lensing = np.array([1.0])
+    cov.ellrange = np.array([10.0])
+    cov.npair_mm = reader.npair_mm.transpose(0, 3, 1, 2)
+    zero = np.zeros((1, 1, 1, 1, 2, 2, 2, 2))
+    cov.covELL_gaussian = lambda *args: [zero] * 18
+    parts = cov._CovTHETASpace__covTHETA_split_gaussian(
+        {}, {"ellipticity_dispersion": s}, False
+    )
+    for offset in (23, 29):
+        noise = parts[offset][0, 0, 0, 0]
+        for (i, j), value in expected.items():
+            np.testing.assert_allclose(noise[i, j, i, j], value, rtol=1e-12)
+            for k, ell in expected:
+                if (i, j) != (k, ell):
+                    assert noise[i, j, k, ell] == 0
+    assert parts[26] == 0  # Pure-noise xi+ / xi- cross covariance vanishes.
diff --git a/workflow/README.md b/workflow/README.md
index fa84a7b4..f1f41d14 100644
--- a/workflow/README.md
+++ b/workflow/README.md
@@ -134,6 +134,31 @@ snakemake -n --profile workflow/profiles/default -s workflow/Snakefile \
     --configfile papers/bmodes/config/config.yaml
 ```
 
+#### Geometry-only pair counts for paired mocks
+
+The default measures weighted pair counts independently for each catalogue.
+A paired mock run can share this non-cosmological geometry calculation after
+verifying **all** selected positions and weights are identical, for example by
+hashing every RA, Dec and weight value in row order (and every bin ID for tomography).
+Declare the reference explicitly in the run config:
+
+```yaml
+covariance:
+  shape_noise: pair_counts
+  pair_counts_reference:
+    GLASS1_blinded: GLASS1_clean
+    GLASS1_psfleak: GLASS1_clean
+    GLASS1_addB: GLASS1_clean
+```
+
+Each consumer depends directly on the clean entry's pair-count table on the same
+angular grid; the workflow does not copy tables into another producer's outputs.
+The reference must be a configured run version and cannot chain through another
+reference.
+This option does not share shear vectors, rho/tau, fitted coefficients or
+harmonic covariances, and it does not authorize sharing between different
+positions, weights or selections.
+
 ### Which `sp_validation` a rule imports: the launched checkout
 
 The image is the frozen *dependency stack*; the `sp_validation` that runs is
diff --git a/workflow/common.py b/workflow/common.py
index 3d98f78b..b2c9f9a1 100644
--- a/workflow/common.py
+++ b/workflow/common.py
@@ -332,6 +332,21 @@ def catalogue_entry(version):
     return CATALOG_CONFIG[catalogue_name(version)]
 
 
+def pair_count_reference(config, version):
+    """Resolve an explicit geometry-only reference without chaining aliases.
+
+    The run must verify identical positions, weights and selection before
+    declaring a reference. The default measures each catalogue independently.
+    """
+    references = config.get("covariance", {}).get("pair_counts_reference", {})
+    reference = references.get(version, version)
+    if reference not in config["versions"]:
+        raise ValueError(f"Pair-count reference {reference!r} is not a run version")
+    if references.get(reference, reference) != reference:
+        raise ValueError("Pair-count references must point directly to a producer")
+    return reference
+
+
 def redshift_path(version):
     """The n(z) file of ``version``: its catalogue entry's ``shear.redshift_path``,
     as written (as ``CosmologyValidation.get_redshift`` reads it)."""
diff --git a/workflow/rules/covariance.smk b/workflow/rules/covariance.smk
index af8f188d..daa1cd0d 100644
--- a/workflow/rules/covariance.smk
+++ b/workflow/rules/covariance.smk
@@ -114,6 +114,35 @@ def _onecov_tool(key, default=None):
     return lookup
 
 
+def _pair_noise_enabled(w):
+    """Whether this covariance uses measured pair counts on its own grid."""
+    default = "pair_counts" if "cosmo_val" in config else "analytic"
+    mode = config.get("covariance", {}).get("shape_noise", default)
+    if mode not in ("pair_counts", "analytic"):
+        raise ValueError("covariance.shape_noise must be pair_counts or analytic")
+    return mode == "pair_counts"
+
+
+def _pair_count_path(w, suffix=".txt", pair=None):
+    # This explicit reference is for catalogues with verified identical
+    # positions, weights and selection. No shear or rho/tau product is shared.
+    reference = pair_count_reference(config, w.version)
+    return str(COSMO_VAL / (
+        f"{reference}_pairs_minsep={w.min_sep}_maxsep={w.max_sep}"
+        f"_nbins={w.nbins}"
+        + (f"_tomo_bin={pair[0]}-{pair[1]}" if pair else "") + suffix
+    ))
+
+
+def _pair_count_paths(w):
+    if not w.tomo_suffix:
+        return [_pair_count_path(w)]
+    # The ordered n(z) columns and survey parameters define the bin order.
+    n = len(get_onecov_cat_params(w.version, w.tomo_suffix)[1].split(","))
+    return [_pair_count_path(w, pair=(i, j))
+            for i in range(1, n + 1) for j in range(i, n + 1)]
+
+
 # Number of cores to use to execute OneCovariance
 ONECOV_THREADS = 8
 
@@ -183,9 +212,55 @@ rule generate_glass_mock_rhotau_samples:
         "../scripts/generate_glass_mock_rhotau_samples.py"
 
 
+rule covariance_pair_counts:
+    """Non-cosmological cache: positions and weights only, no blind access."""
+    input:
+        catalog=get_shear_catalog,
+    output:
+        cache=str(COSMO_VAL / "{version}_pairs_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}.npz"),
+        tabulated=str(COSMO_VAL / "{version}_pairs_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}.txt"),
+    params:
+        ra_col=lambda w: catalogue_entry(w.version)["shear"].get("ra_col", "RA"),
+        dec_col=lambda w: catalogue_entry(w.version)["shear"].get("dec_col", "Dec"),
+        w_col=lambda w: catalogue_entry(w.version)["shear"]["w_col"],
+    threads: 24
+    resources:
+        mem_mb=40000,
+        disk_mb=20000,
+        runtime=600,
+    script:
+        "../scripts/measure_pair_counts.py"
+
+
+rule covariance_pair_counts_tomo:
+    """Non-cosmological cache: positions and weights only, no blind access."""
+    input:
+        catalog=get_shear_catalog,
+    output:
+        cache=str(COSMO_VAL / "{version}_pairs_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_tomo_bin={bin1}-{bin2}.npz"),
+        tabulated=str(COSMO_VAL / "{version}_pairs_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_tomo_bin={bin1}-{bin2}.txt"),
+    params:
+        ra_col=lambda w: catalogue_entry(w.version)["shear"].get("ra_col", "RA"),
+        dec_col=lambda w: catalogue_entry(w.version)["shear"].get("dec_col", "Dec"),
+        w_col=lambda w: catalogue_entry(w.version)["shear"]["w_col"],
+        tomo_bin_col=lambda w: catalogue_entry(w.version)["shear"]["tomo_bin_col"],
+        bin_pair=lambda w: (int(w.bin1), int(w.bin2)),
+    wildcard_constraints:
+        bin1=r"[1-9][0-9]*",
+        bin2=r"[1-9][0-9]*",
+    threads: 24
+    resources:
+        mem_mb=40000,
+        disk_mb=20000,
+        runtime=600,
+    script:
+        "../scripts/measure_pair_counts.py"
+
+
 rule covariance_ini:
     input:
         nz_file=lambda w: onecov_nz_path(w.version, w.tomo_suffix),
+        npair=lambda w: _pair_count_paths(w) if _pair_noise_enabled(w) else [],
         mask=lambda w: (
             [] if w.mask_suffix != "_masked"
             else [get_footprint_path(w.version, w.tomo_suffix)]
@@ -220,6 +295,12 @@ rule covariance_ini:
         nz_dir=lambda w: os.path.dirname(onecov_nz_path(w.version, w.tomo_suffix)),
         nz_file=lambda w: os.path.basename(onecov_nz_path(w.version, w.tomo_suffix)),
         num_cores=ONECOV_THREADS,
+        pair_section=lambda w: (
+            "[tabulated inputs files]\n"
+            f"npair_directory = {COSMO_VAL}\n"
+            "npair_mm_file = " + ", ".join(os.path.basename(p) for p in _pair_count_paths(w))
+            if _pair_noise_enabled(w) else ""
+        ),
     threads: 1
     shell:
         """
@@ -263,6 +344,8 @@ ellipticity_dispersion = {params.sigma_e_param}
 survey_area_lensing_in_deg2 = {params.area}
 n_eff_lensing = {params.n_e}
 
+{params.pair_section}
+
 [redshift]
 zlens_directory = {params.nz_dir}
 zlens_file = {params.nz_file}
diff --git a/workflow/scripts/measure_pair_counts.py b/workflow/scripts/measure_pair_counts.py
new file mode 100644
index 00000000..51d8ed70
--- /dev/null
+++ b/workflow/scripts/measure_pair_counts.py
@@ -0,0 +1,17 @@
+"""Write a non-cosmological pair-count cache on the requested xi grid."""
+
+from sp_validation.pair_counts import measure_pair_counts
+
+measure_pair_counts(
+    snakemake.input.catalog,
+    snakemake.output.cache,
+    min_sep=float(snakemake.wildcards.min_sep),
+    max_sep=float(snakemake.wildcards.max_sep),
+    nbins=int(snakemake.wildcards.nbins),
+    ra_col=snakemake.params.ra_col,
+    dec_col=snakemake.params.dec_col,
+    w_col=snakemake.params.w_col,
+    num_threads=snakemake.threads,
+    tomo_bin_col=snakemake.params.get("tomo_bin_col"),
+    bin_pair=snakemake.params.get("bin_pair"),
+)
diff --git a/workflow/tests/conftest.py b/workflow/tests/conftest.py
index af4e4258..8cd96d55 100644
--- a/workflow/tests/conftest.py
+++ b/workflow/tests/conftest.py
@@ -112,10 +112,8 @@ def snakemake(self, *args, config=(), cwd=None, env=None, timeout=300):
         environment. The workflow profile is off: it enables apptainer, which CI
         runners lack.
         """
-        cmd = [
-            sys.executable,
-            "-m",
-            "snakemake",
+        launcher = os.environ.get("SPV_TEST_SNAKEMAKE")
+        cmd = ([launcher] if launcher else [sys.executable, "-m", "snakemake"]) + [
             "--cores",
             "1",
             "--workflow-profile",
diff --git a/workflow/tests/test_onecov_inputs.py b/workflow/tests/test_onecov_inputs.py
index ff3ae207..276135ac 100644
--- a/workflow/tests/test_onecov_inputs.py
+++ b/workflow/tests/test_onecov_inputs.py
@@ -8,7 +8,7 @@
 
 # Bin integrals 0.5 and 0.625: not unit-normalised.
 NZ_ROWS = "0.0 0.0 0.0\n0.5 1.0 0.5\n1.0 0.0 1.5\n"
-UNMASKED = ['covariance={"default_masked":false}']
+UNMASKED = ['covariance={"default_masked":false,"shape_noise":"analytic"}']
 
 
 def _ini(toy, tomo_suffix=""):
diff --git a/workflow/tests/test_pair_count_reference.py b/workflow/tests/test_pair_count_reference.py
new file mode 100644
index 00000000..89f37825
--- /dev/null
+++ b/workflow/tests/test_pair_count_reference.py
@@ -0,0 +1,31 @@
+"""Geometry-only pair-count sharing is explicit, direct and opt-in."""
+
+import pytest
+
+
+def test_pair_count_reference_is_opt_in(toy):
+    config = {"versions": ["clean", "blinded"]}
+    assert toy.common.pair_count_reference(config, "blinded") == "blinded"
+    config["covariance"] = {"pair_counts_reference": {"blinded": "clean"}}
+    assert toy.common.pair_count_reference(config, "blinded") == "clean"
+    assert toy.common.pair_count_reference(config, "clean") == "clean"
+
+
+def test_pair_count_reference_rejects_missing_producer(toy):
+    config = {
+        "versions": ["blinded"],
+        "covariance": {"pair_counts_reference": {"blinded": "missing"}},
+    }
+    with pytest.raises(ValueError, match="not a run version"):
+        toy.common.pair_count_reference(config, "blinded")
+
+
+def test_pair_count_reference_rejects_chains(toy):
+    config = {
+        "versions": ["clean", "blinded", "injected"],
+        "covariance": {
+            "pair_counts_reference": {"blinded": "injected", "injected": "clean"}
+        },
+    }
+    with pytest.raises(ValueError, match="directly to a producer"):
+        toy.common.pair_count_reference(config, "blinded")
diff --git a/workflow/tests/test_pair_noise.py b/workflow/tests/test_pair_noise.py
new file mode 100644
index 00000000..b3ae3157
--- /dev/null
+++ b/workflow/tests/test_pair_noise.py
@@ -0,0 +1,121 @@
+"""The covariance depends on a shape-free cache by default in cosmo_val."""
+
+from contextlib import contextmanager
+from pathlib import Path
+
+import yaml
+from conftest import VERSIONS, parse_jobs
+
+
+@contextmanager
+def _implicit_noise_setting(toy):
+    # Snakemake recursively merges CLI config; omitting a nested CLI key does
+    # not remove the explicit setting in the paper template.
+    path = toy.rundir / "config" / "config.yaml"
+    original = path.read_text()
+    config = yaml.safe_load(original)
+    config["covariance"].pop("shape_noise", None)
+    path.write_text(yaml.safe_dump(config))
+    try:
+        yield
+    finally:
+        path.write_text(original)
+
+
+def _ini(toy, tomo_suffix=""):
+    # Unmasked target: no site-specific footprint needed by this test.
+    (toy.root / "data" / "nz_SP_v0.1_A.txt").touch()
+    grid = toy.common.xi_grids(toy.config, toy.config["fiducial"])["integration"]
+    return toy.common.covariance_path(
+        VERSIONS[0],
+        "g",
+        grid["min_sep"],
+        grid["max_sep"],
+        grid["nbins"],
+        "",
+        tomo_suffix,
+        suffix=".ini",
+    )
+
+
+def test_pair_count_dependency_and_ini(toy):
+    target = _ini(toy)
+    with _implicit_noise_setting(toy):
+        result = toy.snakemake("-n", target)
+    assert result.returncode == 0, result.stdout
+    jobs = parse_jobs(result.stdout)
+    (pair,) = [j for j in jobs if j.rule == "covariance_pair_counts"]
+    (ini,) = [j for j in jobs if j.rule == "covariance_ini"]
+    assert pair.input == [str(toy.root / "data" / f"{VERSIONS[0]}.fits")]
+    assert any(p.endswith(".txt") and p in ini.input for p in pair.output)
+    assert not any(".sacc" in p or "blind" in p for p in pair.input + pair.output)
+    # Stand in for measured counts, then execute only the host-side ini writer.
+    for p in pair.output:
+        Path(p).parent.mkdir(parents=True, exist_ok=True)
+        Path(p).touch()
+    result = toy.snakemake(target, "--allowed-rules", "covariance_ini")
+    assert result.returncode == 0, result.stdout
+    text = Path(target).read_text()
+    assert "[tabulated inputs files]\nnpair_directory = " in text
+    assert "npair_mm_file = " in text
+    assert "\\n" not in text
+    entries = yaml.safe_load((toy.root / "cosmo_val" / "cat_config.yaml").read_text())
+    sigma_e = float(entries[VERSIONS[0]]["cov_th"]["sigma_e"])
+    # cov_th.sigma_e is per component, as OneCovariance takes it.
+    assert f"ellipticity_dispersion = {sigma_e!r}" in text
+
+
+def test_analytic_flag_removes_pair_dependency(toy):
+    result = toy.snakemake(
+        "-n",
+        "-F",
+        _ini(toy),
+        config=['covariance={"shape_noise":"analytic","default_masked":false}'],
+    )
+    assert result.returncode == 0, result.stdout
+    jobs = parse_jobs(result.stdout)
+    assert not any(j.rule == "covariance_pair_counts" for j in jobs)
+    (ini,) = [j for j in jobs if j.rule == "covariance_ini"]
+    assert not any("_pairs_" in p for p in ini.input)
+
+
+def test_tomographic_counts_and_file_order(toy):
+    from test_onecov_inputs import _two_bin_catalogue
+
+    with _two_bin_catalogue(toy):
+        target = _ini(toy, "_tomo")
+        result = toy.snakemake("-n", target)
+        assert result.returncode == 0, result.stdout
+        jobs = parse_jobs(result.stdout)
+        pairs = [j for j in jobs if j.rule == "covariance_pair_counts_tomo"]
+        assert len(pairs) == 3
+        assert {(j.wildcards["bin1"], j.wildcards["bin2"]) for j in pairs} == {
+            ("1", "1"),
+            ("1", "2"),
+            ("2", "2"),
+        }
+        for job in pairs:
+            for path in job.output:
+                Path(path).parent.mkdir(parents=True, exist_ok=True)
+                Path(path).touch()
+        result = toy.snakemake(target, "--allowed-rules", "covariance_ini")
+        assert result.returncode == 0, result.stdout
+        text = Path(target).read_text()
+        line = next(
+            line for line in text.splitlines() if line.startswith("npair_mm_file")
+        )
+        names = line.split(" = ")[1].split(", ")
+        assert [name.split("_tomo_bin=")[1] for name in names] == [
+            "1-1.txt",
+            "1-2.txt",
+            "2-2.txt",
+        ]
+        assert "survey_area_lensing_in_deg2 = 7.5" in text
+        result = toy.snakemake(
+            "-n", "-F", target, config=['covariance={"shape_noise":"analytic"}']
+        )
+        assert result.returncode == 0, result.stdout
+        assert not any(
+            j.rule.startswith("covariance_pair_counts")
+            for j in parse_jobs(result.stdout)
+        )