added functions to correct for extinction and v6 cuts

des_y6utils/mdet.py

@@ -7,6 +7,52 @@
 import healsparse
 
 
+def add_extinction_correction_columns(fold, fnew):
+    """This function changes and adds columns for extinction corrected magnitudes.
+    pgauss_band_flux_* are the columns that are the corrected fluxes.
+    pgauss_band_flux_*_nodered are the columns that used to be pgauss_band_flux_*.
+
+    Note: This function assumes you're working with fits files, not hdf5.
+
+    Parameters
+    ----------
+    fold : original file path
+    fnew : new file path
+
+    Returns
+    -------
+    None (saves the new file)
+    """
+
+    import fitsio as fio
+    d = fio.read(fold)
+    bands = ['g', 'r', 'i', 'z']
+    fits = fio.FITS(fnew, 'rw')
+
+    # compute dereddened fluxes and
+    # replace the entries of pgauss_flux_band_* with the dereddened fluxes.
+    dustmap = _read_hsp_dustmap_local(
+        '/global/cfs/cdirs/des/y6-shear-catalogs/SFD_dust_4096.hsp'
+    )
+    dered = dustmap.get_values_pos(d["ra"], d["dec"])
+    flux_og = []
+    for ii, b in enumerate(bands):
+        flux = np.copy(d["pgauss_band_flux_" + b])
+        flux_og.append(flux)
+        mag_ = _compute_asinh_dered_mag(d["pgauss_band_flux_" + b], ii, dered)
+        flux_ = _compute_asinh_flux(mag_, ii)
+        d["pgauss_band_flux_" + b] = flux_
+    fits.write(d)
+
+    # make _nodered array with pgauss_band_flux_* entries, and add them to fits.
+    for ii, b in enumerate(bands):
+        fits[-1].insert_column('pgauss_band_flux_' + b + '_nodered', flux_og[ii])
+
+    fits.close()
+
+    return None
+
+
 def make_mdet_cuts(data, version, verbose=False):
     """A function to make the standard metadetection cuts.
 
@@ -40,6 +86,8 @@ def make_mdet_cuts(data, version, verbose=False):
         return _make_mdet_cuts_v4(data, verbose=verbose)
     elif str(version) == "5":
         return _make_mdet_cuts_v5(data, verbose=verbose)
+    elif str(version) == "6":
+        return _make_mdet_cuts_v6(data, verbose=verbose)
     else:
         raise ValueError("the mdet cut version '%r' is not recognized!" % version)
 
@@ -52,7 +100,7 @@ def _make_mdet_cuts_gauss(
     min_t_ratio=0.5,
     max_mfrac=0.1,
     max_s2n=np.inf,
-    mask_name="y6-combined-hleda-gaiafull-des-stars-hsmap16384-nomdet-v3.fits",
+    mask_name="y6-combined-hleda-gaiafull-des-stars-hsmap131k-mdet-v1.hsp",
     max_t=100.0,
 ):
     """
@@ -61,7 +109,7 @@ def _make_mdet_cuts_gauss(
     msk : np.ndarray of bool
         A boolean array with the cuts. To cut the data, use `data[msk]`.
     """
-    msk = _make_mdet_cuts_raw_v345(
+    msk = _make_mdet_cuts_raw_v6(
         data,
         verbose=verbose,
         min_s2n=min_s2n,
@@ -212,8 +260,8 @@ def _make_mdet_cuts_v4(d, verbose=False):
         max_t=np.inf,
     )
 
-    size_sizeerr = (d['gauss_T_ratio']*d['gauss_psf_T']) * d['gauss_T_err']
-    size_s2n = (d['gauss_T_ratio']*d['gauss_psf_T']) / d['gauss_T_err']
+    size_sizeerr = (d['gauss_T_ratio'] * d['gauss_psf_T']) * d['gauss_T_err']
+    size_s2n = (d['gauss_T_ratio'] * d['gauss_psf_T']) / d['gauss_T_err']
     msk_superspreader = ((size_sizeerr > 1) & (size_s2n < 10))
     msk &= ~msk_superspreader
 
@@ -242,8 +290,8 @@ def _make_mdet_cuts_v5(d, verbose=False):
         max_t=np.inf,
     )
 
-    size_sizeerr = (d['gauss_T_ratio']*d['gauss_psf_T']) * d['gauss_T_err']
-    size_s2n = (d['gauss_T_ratio']*d['gauss_psf_T']) / d['gauss_T_err']
+    size_sizeerr = (d['gauss_T_ratio'] * d['gauss_psf_T']) * d['gauss_T_err']
+    size_s2n = (d['gauss_T_ratio'] * d['gauss_psf_T']) / d['gauss_T_err']
     msk_superspreader = ((size_sizeerr > 1) & (size_s2n < 10))
     msk &= ~msk_superspreader
 
@@ -259,6 +307,36 @@ def _make_mdet_cuts_v5(d, verbose=False):
     return msk
 
 
+def _make_mdet_cuts_v6(d, verbose=False):
+
+    msk = _make_mdet_cuts_raw_v6(
+        d,
+        verbose=verbose,
+        min_s2n=10,
+        n_terr=0,
+        min_t_ratio=0.5,
+        max_mfrac=0.1,
+        max_s2n=np.inf,
+        max_t=20.0,
+    )
+
+    size_sizeerr = (d['gauss_T_ratio'] * d['gauss_psf_T']) * d['gauss_T_err']
+    size_s2n = (d['gauss_T_ratio'] * d['gauss_psf_T']) / d['gauss_T_err']
+    msk_superspreader = ((size_sizeerr > 1) & (size_s2n < 10))
+    msk &= ~msk_superspreader
+
+    # apply the mask
+    hmap = _read_hsp_mask(
+        "y6-combined-hleda-gaiafull-des-stars-hsmap131k-mdet-v1.hsp"
+    )
+    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.
 
@@ -291,6 +369,40 @@ def _compute_asinh_mags(flux, i):
     return mag
 
 
+def _compute_asinh_dered_mag(flux, i, map):
+
+    dered_fac = [3.186, 2.140, 1.196, 1.048]
+    mag = _compute_asinh_mags(flux, i)
+    dered_mag = mag - dered_fac[i] * map
+
+    return dered_mag
+
+
+def _compute_asinh_flux(mag, i):
+    """This function does the inverse process of compute_asinh_mag.
+
+    Parameters
+    ----------
+    mag : float or np.ndarray
+        The magnitude.
+    i : int
+        The index of the band in griz (i.e., 0 for g, 1 for r, 2 for i, 3 for z).
+
+    Returns
+    -------
+    flux : float or np.ndarray
+        The asinh flux for the mag.
+    """
+
+    zp = 30.0
+    b_array = np.array([3.27e-12, 4.83e-12, 6.0e-12, 9.0e-12])
+    bscale = np.array(b_array) * 10.**(zp / 2.5)
+    A = (2.5 * np.log10(1.0 / b_array[i]) - mag) * 0.4 * np.log(10)
+    flux = 2 * bscale[i] * np.sinh(A)
+
+    return flux
+
+
 @lru_cache
 def _read_hsp_mask(fname):
     mpth = _get_mask_path(fname)
@@ -302,6 +414,11 @@ def _read_hsp_mask_local(fname):
     return healsparse.HealSparseMap.read(fname)
 
 
+@lru_cache
+def _read_hsp_dustmap_local(fname):
+    return healsparse.HealSparseMap.read(fname)
+
+
 def _get_mask_path(fname):
     # get or make meds dir
     meds_dir = os.environ.get("MEDS_DIR", None)
@@ -431,11 +548,11 @@ def _make_mdet_cuts_raw_v12(
         & (mag_r < 26.5)
         & (mag_i < 26.2)
         & (mag_z < 25.6)
-        & (d["pgauss_T"] < (1.9 - 2.8*d["pgauss_T_err"]))
+        & (d["pgauss_T"] < (1.9 - 2.8 * d["pgauss_T_err"]))
         & (
             d["wmom_T_ratio"] >= np.maximum(
                 min_t_ratio,
-                (1.0 + n_terr*d["wmom_T_err"]/d["wmom_psf_T"])
+                (1.0 + n_terr * d["wmom_T_err"] / d["wmom_psf_T"])
             )
         )
     )
@@ -514,11 +631,98 @@ def _make_mdet_cuts_raw_v345(
         & (mag_r < 26.5)
         & (mag_i < 26.2)
         & (mag_z < 25.6)
-        & (d["pgauss_T"] < (1.2 - 3.1*d["pgauss_T_err"]))
+        & (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_ratio"] * d["gauss_psf_T"]) < max_t)
+    )
+
+    if verbose:
+        print("did mdet cuts", np.sum(msk))
+
+    return msk
+
+
+def _make_mdet_cuts_raw_v6(
+    d,
+    *,
+    min_s2n,
+    n_terr,
+    min_t_ratio,
+    max_mfrac,
+    max_s2n,
+    max_t,
+    verbose=False,
+):
+    """The raw v6 cuts come from analysis in Fall of 2023.
+
+      - We implemented i-band magnitude cut at 24.5. TO-DO: still need to test.
+      - We implemented galaxy size cut at T=20.
+      - We updated a cut in the pgauss T-Terr plane.
+    """
+
+    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 < 24.5)
+        & (mag_z < 25.6)
+        & (d["pgauss_T"] < (1.6 - 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"])
+                (n_terr * d["gauss_T_err"] / d["gauss_psf_T"])
             )
         )
         & ((d["gauss_T_ratio"] * d["gauss_psf_T"]) < max_t)