fix(ngmix): emit true half-light radii in r50 columns; dedupe PSF size columns

src/shapepipe/modules/make_cat_package/make_cat.py

@@ -449,7 +449,7 @@ def _save_ngmix_data(self, ngmix_cat_path, moments=False):
                     )
                     self._add2dict(f"NGMIX{m}_ELL_PSFo_{key}", g_psf, idx)
 
-                    t_psfo = ncf_data["T_psfo_ngmix"][ind[0]]
+                    t_psfo = ncf_data["Tpsf"][ind[0]]
                     self._add2dict(f"NGMIX{m}_T_PSFo_{key}", t_psfo, idx)
 
                 self._add2dict(

src/shapepipe/modules/ngmix_package/ngmix.py

@@ -19,6 +19,11 @@
 
 from shapepipe.pipeline import file_io
 
+# Gaussian half-light radius per unit sigma: r50 = sqrt(2 ln 2) * sigma
+# = 1.17741 * sigma.  With the ngmix area parameter T = 2 sigma^2 this
+# gives r50 = SIGMA_TO_R50 * sqrt(T / 2) = sqrt(ln 2 * T).
+SIGMA_TO_R50 = np.sqrt(2.0 * np.log(2.0))
+
 
 def get_prior(pixel_scale, rng, T_range=None, F_range=None):
     """Build ngmix joint prior for a 6-parameter galaxy model.
@@ -342,22 +347,20 @@ def compile_results(self, results):
             'ntry_fit',
             'g1_psfo_ngmix',
             'g2_psfo_ngmix',
-            'T_psfo_ngmix',
-            'T_err_psfo_ngmix',
-            'r50_psfo_ngmix',
             'g1_err_psfo_ngmix',
             'g2_err_psfo_ngmix',
-            'r50_err_psfo_ngmix',
             'g1',
             'g1_err',
             'g2',
             'g2_err',
             'T',
             'T_err',
             'Tpsf',
+            'Tpsf_err',
             'r50',
             'r50_err',
             'r50psf',
+            'r50psf_err',
             'g1_psf',
             'g2_psf',
             'flux',
@@ -402,24 +405,28 @@ def compile_results(self, results):
                 output_dict[name]["g2_err_psfo_ngmix"].append(
                     results[idx]["g_err_PSFo"][1]
                 )
-                output_dict[name]["T_psfo_ngmix"].append(results[idx]["T_PSFo"])
-                output_dict[name]["T_err_psfo_ngmix"].append(
-                    results[idx]["T_err_PSFo"]
-                )
-                output_dict[name]['r50_psfo_ngmix'].append(
-                    results[idx]['r50_PSFo']
-                )
-                output_dict[name]['r50_err_psfo_ngmix'].append(
-                    results[idx]['r50_err_PSFo']
-                )
                 output_dict[name]["T"].append(results[idx][name]["T"])
                 output_dict[name]["T_err"].append(results[idx][name]["T_err"])
                 output_dict[name]["Tpsf"].append(results[idx]["T_PSFo"])
+                output_dict[name]["Tpsf_err"].append(results[idx]["T_err_PSFo"])
                 output_dict[name]["g1_psf"].append(results[idx]["g_PSFo"][0])
                 output_dict[name]["g2_psf"].append(results[idx]["g_PSFo"][1])
-                output_dict[name]['r50'].append(results[idx][name]['pars'][4])
-                output_dict[name]['r50_err'].append(results[idx][name]['pars_err'][4])
+
+                # Galaxy half-light radius from the fitted area T = 2 sigma^2:
+                # r50 = sqrt(ln 2 * T), with d r50 / d T = r50 / (2 T)
+                T_gal = results[idx][name]["T"]
+                T_gal_err = results[idx][name]["T_err"]
+                if T_gal > 0:
+                    r50_gal = SIGMA_TO_R50 * np.sqrt(T_gal / 2)
+                    r50_gal_err = r50_gal * T_gal_err / (2 * T_gal)
+                else:
+                    r50_gal = r50_gal_err = np.nan
+                output_dict[name]['r50'].append(r50_gal)
+                output_dict[name]['r50_err'].append(r50_gal_err)
                 output_dict[name]['r50psf'].append(results[idx]["r50_PSFo"])
+                output_dict[name]['r50psf_err'].append(
+                    results[idx]["r50_err_PSFo"]
+                )
                 output_dict[name]["g1"].append(results[idx][name]["g"][0])
                 output_dict[name]["g2"].append(results[idx][name]["g"][1])
                 output_dict[name]["g1_err"].append(
@@ -632,14 +639,17 @@ def process(self):
                 1 for k in ['noshear', '1p', '1m', '2p', '2m']
                 if res.get(k, {}).get('flags', 0) != 0
             )
-            r50_psfo = np.sqrt(max(psf_res['T_psf'], 0) / 2)
+            # PSF half-light radius r50 = sqrt(2 ln 2) * sigma with
+            # sigma = sqrt(T / 2); error from d sigma / d T = 1 / (4 sigma)
+            sigma_psfo = np.sqrt(max(psf_res['T_psf'], 0) / 2)
             res['g_PSFo'] = psf_res['g_psf']
             res['g_err_PSFo'] = psf_res['g_psf_err']
             res['T_PSFo'] = psf_res['T_psf']
             res['T_err_PSFo'] = psf_res['T_psf_err']
-            res['r50_PSFo'] = r50_psfo
+            res['r50_PSFo'] = SIGMA_TO_R50 * sigma_psfo
             res['r50_err_PSFo'] = (
-                psf_res['T_psf_err'] / (2 * r50_psfo) if r50_psfo > 0 else np.nan
+                SIGMA_TO_R50 * psf_res['T_psf_err'] / (4 * sigma_psfo)
+                if sigma_psfo > 0 else np.nan
             )
             res['mcal_flags'] = 0
             final_res.append(res)
@@ -863,7 +873,7 @@ def get_noise(gal, weight, guess, pixel_scale, thresh=1.2):
         Weight image
     guess : list
         Gaussian parameters fot the window function
-        ``[x0, y0, g1, g2, r50, flux]``
+        ``[x0, y0, g1, g2, T, flux]``
     pixel_scale : float
         Pixel scale of the galaxy image
     thresh : float, optional

src/shapepipe/tests/test_ngmix.py

@@ -135,6 +135,110 @@ def test_metacal_is_reproducible_with_fixed_seed():
     npt.assert_array_equal(_metacal_noshear_g(42), _metacal_noshear_g(42))
 
 
+def _fake_metacal_result(T, T_err, T_psf, T_psf_err):
+    """Build one minimal metacal result as produced by the process loop.
+
+    The PSF size entries mirror the conversion done at the source:
+    ``r50_PSFo = sqrt(2 ln 2) * sigma`` with ``sigma = sqrt(T_psf / 2)``.
+    """
+    from shapepipe.modules.ngmix_package.ngmix import SIGMA_TO_R50
+
+    sigma_psf = np.sqrt(max(T_psf, 0) / 2)
+    per_type = {
+        "nfev": 1,
+        "g": [0.01, -0.02],
+        "g_cov": np.diag([1e-4, 1e-4]),
+        "T": T,
+        "T_err": T_err,
+        "flux": 100.0,
+        "flux_err": 1.0,
+        "s2n": 50.0,
+        "flags": 0,
+    }
+    res = {
+        "obj_id": 1,
+        "n_epoch_model": 1,
+        "moments_fail": 0,
+        "g_PSFo": [0.001, 0.002],
+        "g_err_PSFo": [1e-5, 1e-5],
+        "T_PSFo": T_psf,
+        "T_err_PSFo": T_psf_err,
+        "r50_PSFo": SIGMA_TO_R50 * sigma_psf,
+        "r50_err_PSFo": (
+            SIGMA_TO_R50 * T_psf_err / (4 * sigma_psf)
+            if sigma_psf > 0
+            else np.nan
+        ),
+        "mcal_flags": 0,
+    }
+    res.update(
+        {name: dict(per_type) for name in ("1m", "1p", "2m", "2p", "noshear")}
+    )
+    return res
+
+
+def test_compile_results_size_columns_are_half_light_radii():
+    """Every r50 column is a true half-light radius, on both sides.
+
+    Galaxy ``r50 = sqrt(ln 2 * T)`` (not the raw area ``pars[4]``), PSF
+    ``r50psf = sqrt(2 ln 2) * sigma_psf`` (not bare sigma), and the
+    redundant ``*_psfo_ngmix`` size duplicates are gone.
+    """
+    from shapepipe.modules.ngmix_package.ngmix import Ngmix, SIGMA_TO_R50
+
+    T, T_err = 0.18, 0.02
+    T_psf, T_psf_err = 0.09, 0.001
+
+    inst = object.__new__(Ngmix)
+    inst._zero_point = 30.0
+    out = inst.compile_results([_fake_metacal_result(T, T_err, T_psf, T_psf_err)])
+
+    noshear = out["noshear"]
+
+    # galaxy: r50 = sqrt(ln2 * T) = 1.17741 * sqrt(T / 2), error dT * r50/(2T)
+    r50_expected = np.sqrt(np.log(2) * T)
+    npt.assert_allclose(noshear["r50"], [r50_expected])
+    npt.assert_allclose(noshear["r50_err"], [r50_expected * T_err / (2 * T)])
+
+    # PSF: r50psf = sqrt(2 ln2) * sigma, sigma = sqrt(T_psf / 2)
+    sigma_psf = np.sqrt(T_psf / 2)
+    npt.assert_allclose(noshear["r50psf"], [SIGMA_TO_R50 * sigma_psf])
+    npt.assert_allclose(
+        noshear["r50psf_err"], [SIGMA_TO_R50 * T_psf_err / (4 * sigma_psf)]
+    )
+
+    # galaxy/PSF r50 are now commensurable: same convention on both sides
+    npt.assert_allclose(
+        np.array(noshear["r50"]) / np.array(noshear["r50psf"]),
+        [np.sqrt(T / T_psf)],
+    )
+
+    # areas pass through untouched, with the err asymmetry fixed
+    npt.assert_allclose(noshear["Tpsf"], [T_psf])
+    npt.assert_allclose(noshear["Tpsf_err"], [T_psf_err])
+
+    # the *_psfo_ngmix size duplicates are retired
+    for retired in (
+        "T_psfo_ngmix",
+        "T_err_psfo_ngmix",
+        "r50_psfo_ngmix",
+        "r50_err_psfo_ngmix",
+    ):
+        assert retired not in noshear
+
+
+def test_compile_results_nonpositive_T_gives_nan_r50():
+    """A non-positive fitted area cannot yield a real half-light radius."""
+    from shapepipe.modules.ngmix_package.ngmix import Ngmix
+
+    inst = object.__new__(Ngmix)
+    inst._zero_point = 30.0
+    out = inst.compile_results([_fake_metacal_result(-0.05, 0.02, 0.09, 0.001)])
+
+    assert np.isnan(out["noshear"]["r50"]).all()
+    assert np.isnan(out["noshear"]["r50_err"]).all()
+
+
 def test_weight_map_recovers_injected_inverse_variance():
     """A supplied inverse-variance map must not be renormalized twice."""
     from shapepipe.modules.ngmix_package.ngmix import prepare_ngmix_weights