diff --git a/des_y6utils/mdet.py b/des_y6utils/mdet.py index 94a24ca..9386368 100644 --- a/des_y6utils/mdet.py +++ b/des_y6utils/mdet.py @@ -34,11 +34,56 @@ def make_mdet_cuts(data, version, verbose=False): return _make_mdet_cuts_v1(data, verbose=verbose) elif str(version) == "2": return _make_mdet_cuts_v2(data, verbose=verbose) + elif str(version) == "3": + return _make_mdet_cuts_v3(data, verbose=verbose) else: raise ValueError("the mdet cut version '%r' is not recognized!" % version) -def _make_mdet_cuts( +def _make_mdet_cuts_gauss( + data, + verbose=False, + min_s2n=10, + n_terr=0, + min_t_ratio=0.5, + max_mfrac=0.1, + max_s2n=np.inf, + mask_name="y6-combined-hleda-gaiafull-des-stars-hsmap16384-nomdet-v2.fits", + max_t=100.0, +): + """ + Returns + ------- + msk : np.ndarray of bool + A boolean array with the cuts. To cut the data, use `data[msk]`. + """ + msk = _make_mdet_cuts_raw_v3( + data, + verbose=verbose, + min_s2n=min_s2n, + n_terr=n_terr, + min_t_ratio=min_t_ratio, + max_mfrac=max_mfrac, + max_s2n=max_s2n, + max_t=max_t, + ) + + # apply the mask + if os.path.exists(mask_name): + hmap = _read_hsp_mask_local(mask_name) + else: + hmap = _read_hsp_mask(mask_name) + in_footprint = hmap.get_values_pos( + data["ra"], data["dec"], valid_mask=True + ) + msk &= in_footprint + if verbose: + print("did mask cuts", np.sum(msk)) + + return msk + + +def _make_mdet_cuts_wmom( data, verbose=False, min_s2n=10, @@ -125,6 +170,31 @@ def _make_mdet_cuts_v2(d, verbose=False): return msk +def _make_mdet_cuts_v3(d, verbose=False): + + msk = _make_mdet_cuts_raw_v3( + d, + verbose=verbose, + min_s2n=10, + n_terr=0, + min_t_ratio=0.5, + max_mfrac=0.1, + max_s2n=np.inf, + max_t=np.inf, + ) + + # apply the mask + hmap = _read_hsp_mask( + "y6-combined-hleda-gaiafull-des-stars-hsmap16384-nomdet-v2.fits" + ) + in_footprint = hmap.get_values_pos(d["ra"], d["dec"], valid_mask=True) + msk &= in_footprint + if verbose: + print("did mask cuts", np.sum(msk)) + + return msk + + def _compute_asinh_mags(flux, i): """This function and coefficients are from from Eli. Ask him. @@ -309,3 +379,87 @@ def _make_mdet_cuts_raw_v12( print("did mdet cuts", np.sum(msk)) return msk + + +def _make_mdet_cuts_raw_v3( + d, + *, + min_s2n, + n_terr, + min_t_ratio, + max_mfrac, + max_s2n, + max_t, + verbose=False, +): + """The raw v3 cuts come from analysis in Fall of 2022. + + - We use a cut in the pgauss T-Terr plane. This cut removes junk detections + near the wings of stars. For this cut we require pgauss_T_flags == 0 as well. + """ + + mag_g = _compute_asinh_mags(d["pgauss_band_flux_g"], 0) + mag_r = _compute_asinh_mags(d["pgauss_band_flux_r"], 1) + mag_i = _compute_asinh_mags(d["pgauss_band_flux_i"], 2) + mag_z = _compute_asinh_mags(d["pgauss_band_flux_z"], 3) + gmr = mag_g - mag_r + rmi = mag_r - mag_i + imz = mag_i - mag_z + + msk = np.ones(d.shape[0]).astype(bool) + + potential_flag_columns = [ + "psfrec_flags", + "gauss_flags", + "pgauss_T_flags", + "pgauss_band_flux_flags_g", + "pgauss_band_flux_flags_r", + "pgauss_band_flux_flags_i", + "pgauss_band_flux_flags_z", + "mask_flags", + ] + for col in potential_flag_columns: + if col in d.dtype.names: + msk &= (d[col] == 0) + if verbose: + print("did cut " + col, np.sum(msk)) + + if "shear_bands" in d.dtype.names: + msk &= (d["shear_bands"] == "123") + if verbose: + print("did cut shear_bands", np.sum(msk)) + + if "pgauss_s2n" in d.dtype.names: + msk &= (d["pgauss_s2n"] > 5) + if verbose: + print("did cut pgauss_s2n", np.sum(msk)) + + # now do the rest + msk &= ( + (d["gauss_s2n"] > min_s2n) + & (d["gauss_s2n"] < max_s2n) + & (d["mfrac"] < max_mfrac) + & (np.abs(gmr) < 5) + & (np.abs(rmi) < 5) + & (np.abs(imz) < 5) + & np.isfinite(mag_g) + & np.isfinite(mag_r) + & np.isfinite(mag_i) + & np.isfinite(mag_z) + & (mag_g < 26.5) + & (mag_r < 26.5) + & (mag_i < 26.2) + & (mag_z < 25.6) + & (d["pgauss_T"] < (1.2 - 3.1*d["pgauss_T_err"])) + & ( + d["gauss_T_ratio"] >= np.maximum( + min_t_ratio, + (n_terr*d["gauss_T_err"]/d["gauss_psf_T"]) + ) + ) + & (d["gauss_T"] < max_t) + ) + if verbose: + print("did mdet cuts", np.sum(msk)) + + return msk