Tomographic covariances use the binned sample's own footprint and area

src/sp_validation/cosmo_val/catalog_characterization.py

@@ -145,6 +145,21 @@ def area(self):
             self.calculate_area()
         return self._area
 
+    def tomographic_footprint(self, ver):
+        """Footprint map of ``ver``'s tomographic (binned) sample.
+
+        ``cov_th.tomo.mask`` when the entry has one, since bins can need data
+        the whole sample lacks on part of the sky; otherwise the entry's mask.
+        """
+        tomo = self.cc[ver].get("cov_th", {}).get("tomo", {})
+        return tomo.get("mask", self.cc[ver].get("mask"))
+
+    def tomographic_area(self, ver):
+        """Area in deg² of ``ver``'s tomographic footprint."""
+        if "mask" in self.cc[ver].get("cov_th", {}).get("tomo", {}):
+            return self._area_from_mask(self.tomographic_footprint(ver))
+        return self.area[ver]
+
     @property
     def n_eff_gal(self):
         if not hasattr(self, "_n_eff_gal"):
@@ -224,18 +239,20 @@ def calculate_n_eff_gal(self, tomography=False):
             else:
                 tomo_bin_ids, tomo_bin_pairs = ["all"], [("all", "all")]
 
+            # Each density is per unit area of the footprint its sample covers.
+            area = self.tomographic_area(ver) if tomography else self.area[ver]
             with self.results[ver].temporarily_read_data():
                 w = self._read_shear_cols(ver, "w_col")
                 for tomo_bin_id in tomo_bin_ids:
                     if tomo_bin_id == "all":
                         self._n_eff_gal[ver][f"tomo_bin_{tomo_bin_id}"] = n_eff_density(
-                            w, self.area[ver]
+                            w, area
                         )
                     else:
                         tomo_bin = self._read_shear_cols(ver, "tomo_bin_col")
                         mask = tomo_bin == tomo_bin_id
                         self._n_eff_gal[ver][f"tomo_bin_{tomo_bin_id}"] = n_eff_density(
-                            w[mask], self.area[ver]
+                            w[mask], area
                         )
                     print(
                         f"n_eff_gal for tomo_bin_{tomo_bin_id} = {self._n_eff_gal[ver][f'tomo_bin_{tomo_bin_id}']:.2f} gal./arcmin^-2"

src/sp_validation/cosmo_val/pseudo_cl.py

@@ -778,8 +778,16 @@ def calculate_pseudo_cl_onecovariance(self, tomography=False):
                 self.print_done(f"Skipping OneCovariance calculation, {out_dir} exists")
                 self._load_onecovariance_cov(out_dir, ver, tomography)
             else:
-                mask_path = self.cc[ver]["mask"]
+                mask_path = (
+                    self.tomographic_footprint(ver)
+                    if tomography
+                    else self.cc[ver]["mask"]
+                )
                 if not os.path.exists(mask_path):
+                    if mask_path != self.cc[ver]["mask"]:
+                        raise FileNotFoundError(
+                            f"{ver}: tomographic footprint {mask_path} not found"
+                        )
                     print("Mask file does not exist")
                     print("Computing the mask from the binned catalog and saving...")
                     mask = self._get_binned_catalog_mask(ver)
@@ -1126,7 +1134,8 @@ def _modify_onecov_config(
         mask_folder = os.path.dirname(os.path.abspath(mask_path))
         config["survey specs"]["mask_directory"] = mask_folder
         config["survey specs"]["mask_file_lensing"] = mask_base
-        config["survey specs"]["survey_area_lensing_in_deg2"] = str(self.area[ver])
+        area = self.tomographic_area(ver) if tomography else self.area[ver]
+        config["survey specs"]["survey_area_lensing_in_deg2"] = str(area)
 
         # Update ellipticity dispersion and effective number density
         # Account for tomography if needed

workflow/rules/covariance.smk

@@ -4,10 +4,15 @@
 # DEFAULT_MASK_SUFFIX defined in Snakefile
 
 
-def get_footprint_path(version):
-    """Healpix FITS footprint map of ``version``'s catalogue entry, which
-    OneCovariance reads as ``mask_file_lensing``."""
-    return catalogue_entry(version)["mask"]
+def get_footprint_path(version, tomo_suffix=""):
+    """Healpix FITS footprint map OneCovariance reads as ``mask_file_lensing``.
+
+    A tomographic covariance reads ``cov_th.tomo.mask`` when the entry has one:
+    the binned sample can cover less sky than the whole catalogue (bins need
+    photometry the whole sample does not). Otherwise the entry's ``mask``."""
+    entry = catalogue_entry(version)
+    tomo = entry.get("cov_th", {}).get("tomo", {}) if tomo_suffix else {}
+    return tomo.get("mask", entry["mask"])
 
 
 # ---------------------------------------------------------------------------
@@ -20,13 +25,20 @@ def get_onecov_cat_params(version, tomo_suffix=""):
     catalogue entry, as comma-separated strings with one entry per bin.
 
     A tomographic covariance reads ``cov_th.tomo`` (per-bin ``sigma_e`` and
-    ``n_e``) when the entry has one; a non-tomographic covariance reads the
-    whole-sample ``cov_th.sigma_e`` and ``cov_th.n_e``. A tomographic covariance
+    ``n_e``, and the binned sample's area ``A`` if it differs from the whole
+    sample's) when the entry has one; a non-tomographic covariance reads the
+    whole-sample ``cov_th.sigma_e``, ``cov_th.n_e`` and ``cov_th.A``. Each
+    ``n_e`` must be computed on the area it is paired with. A tomographic covariance
     of an entry that declares ``shear.tomo_bin_col`` refuses single
     whole-sample values, which would stand in for every bin."""
     entry = catalogue_entry(version)
     cov_th = entry["cov_th"]
     survey = cov_th.get("tomo", cov_th) if tomo_suffix else cov_th
+    if survey is not cov_th and ("A" in survey) != ("mask" in survey):
+        raise ValueError(
+            f"{version}: cov_th.tomo sets A and mask together: the binned "
+            "sample's area and the footprint it is measured on"
+        )
 
     def per_bin(value):
         values = value if isinstance(value, (list, tuple)) else str(value).split(",")
@@ -39,7 +51,7 @@ def get_onecov_cat_params(version, tomo_suffix=""):
             f"{version}: a tomographic covariance needs per-bin n_e and sigma_e "
             "in cov_th.tomo"
         )
-    return ", ".join(sigma_e), ", ".join(n_e), cov_th["A"]
+    return ", ".join(sigma_e), ", ".join(n_e), survey.get("A", cov_th["A"])
 
 
 def onecov_nz_path(version, tomo_suffix=""):
@@ -174,7 +186,10 @@ rule generate_glass_mock_rhotau_samples:
 rule covariance_ini:
     input:
         nz_file=lambda w: onecov_nz_path(w.version, w.tomo_suffix),
-        mask=lambda w: [] if w.mask_suffix != "_masked" else [get_footprint_path(w.version)],
+        mask=lambda w: (
+            [] if w.mask_suffix != "_masked"
+            else [get_footprint_path(w.version, w.tomo_suffix)]
+        ),
     output:
         cov_output(".ini"),
     params:
@@ -193,12 +208,12 @@ rule covariance_ini:
         area=lambda w: get_onecov_cat_params(w.version, w.tomo_suffix)[2],
         # OneCovariance takes directory and file separately.
         mask_dir=lambda w: (
-            os.path.dirname(get_footprint_path(w.version))
+            os.path.dirname(get_footprint_path(w.version, w.tomo_suffix))
             if w.mask_suffix == "_masked"
             else ""
         ),
         mask_file=lambda w: (
-            os.path.basename(get_footprint_path(w.version))
+            os.path.basename(get_footprint_path(w.version, w.tomo_suffix))
             if w.mask_suffix == "_masked"
             else ""
         ),

workflow/tests/test_onecov_inputs.py

@@ -38,7 +38,14 @@ def _two_bin_catalogue(toy, tomo_cov_th=True, nz_rows=NZ_ROWS):
         entries[VERSIONS[0]]["cov_th"]["tomo"] = {
             "n_e": [2.0, 3.0],
             "sigma_e": [0.25, 0.35],
+            "A": 7.5,
+            "mask": str(toy.root / "data" / "tomo_footprint.fits"),
         }
+        if isinstance(tomo_cov_th, dict):  # overrides; None drops the key
+            tomo = entries[VERSIONS[0]]["cov_th"]["tomo"]
+            tomo.update(tomo_cov_th)
+            for key in [k for k, v in tomo_cov_th.items() if v is None]:
+                tomo.pop(key)
     entries[VERSIONS[0]]["shear"]["tomo_bin_col"] = "TOM_BIN_ID"
     cat_config.write_text(yaml.safe_dump(entries))
     try:
@@ -53,8 +60,8 @@ def _two_bin_catalogue(toy, tomo_cov_th=True, nz_rows=NZ_ROWS):
 
 def test_tomographic_catalogue_covariance_inputs(toy):
     """One multi-bin entry serves both covariances: the tomographic ini reads its
-    per-bin n(z) and cov_th.tomo, the non-tomographic one the bin-summed n(z)
-    and the whole-sample cov_th."""
+    per-bin n(z) and cov_th.tomo (with the binned sample's own area), the
+    non-tomographic one the bin-summed n(z) and the whole-sample cov_th."""
     with _two_bin_catalogue(toy) as cov_th:
         texts = {}
         for suffix in ("", "_tomo"):
@@ -72,6 +79,8 @@ def test_tomographic_catalogue_covariance_inputs(toy):
     assert "zlens_file = nz_SP_v0.1_A.txt" in texts["_tomo"]
     assert "ellipticity_dispersion = 0.25, 0.35" in texts["_tomo"]
     assert "n_eff_lensing = 2.0, 3.0" in texts["_tomo"]
+    assert "survey_area_lensing_in_deg2 = 7.5\n" in texts["_tomo"]
+    assert f"survey_area_lensing_in_deg2 = {cov_th['A']}\n" in texts[""]
     assert f"ellipticity_dispersion = {float(cov_th['sigma_e'])!r}\n" in texts[""]
     assert f"n_eff_lensing = {float(cov_th['n_e'])!r}\n" in texts[""]
     summed = toy.cosmo_val / "nz" / f"{VERSIONS[0]}_nz_bin_sum.txt"
@@ -93,6 +102,15 @@ def test_tomographic_covariance_refuses_whole_sample_cov_th(toy):
     assert "per-bin n_e and sigma_e in cov_th.tomo" in result.stdout
 
 
+def test_tomographic_area_needs_its_footprint(toy):
+    """cov_th.tomo.A without cov_th.tomo.mask is refused: the area and the
+    footprint OneCovariance reads would describe different samples."""
+    with _two_bin_catalogue(toy, tomo_cov_th={"mask": None}):
+        result = toy.snakemake("-n", _ini(toy, "_tomo"), config=UNMASKED)
+    assert result.returncode != 0
+    assert "sets A and mask together" in result.stdout
+
+
 def test_unit_normalised_bins_refused(toy):
     """Bins that each integrate to one carry no weights; summing them is refused."""
     with _two_bin_catalogue(toy, nz_rows="0.0 0.0 0.0\n0.5 2.0 2.0\n1.0 0.0 0.0\n"):