diff --git a/des_y6utils/mdet.py b/des_y6utils/mdet.py index c8cf956..cbe39ee 100644 --- a/des_y6utils/mdet.py +++ b/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, @@ -258,6 +306,35 @@ 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 +368,41 @@ 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) @@ -301,6 +413,10 @@ def _read_hsp_mask(fname): 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 @@ -528,3 +644,90 @@ def _make_mdet_cuts_raw_v345( 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"]) + ) + ) + & ((d["gauss_T_ratio"] * d["gauss_psf_T"]) < max_t) + ) + + if verbose: + print("did mdet cuts", np.sum(msk)) + + return msk \ No newline at end of file