diff --git a/src/shapepipe/modules/make_cat_package/make_cat.py b/src/shapepipe/modules/make_cat_package/make_cat.py index 72024b3a4..890e24713 100644 --- a/src/shapepipe/modules/make_cat_package/make_cat.py +++ b/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( diff --git a/src/shapepipe/modules/ngmix_package/ngmix.py b/src/shapepipe/modules/ngmix_package/ngmix.py index 6ab987005..16c24d876 100644 --- a/src/shapepipe/modules/ngmix_package/ngmix.py +++ b/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,12 +347,8 @@ 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', @@ -355,9 +356,11 @@ def compile_results(self, results): '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 diff --git a/src/shapepipe/tests/test_ngmix.py b/src/shapepipe/tests/test_ngmix.py index 819489b9d..a5ce0277e 100644 --- a/src/shapepipe/tests/test_ngmix.py +++ b/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