Merge pull request #11 from des-science/latest-mdet

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