diff --git a/.gitignore b/.gitignore index b6e4761..8ba67bc 100644 --- a/.gitignore +++ b/.gitignore @@ -127,3 +127,5 @@ dmypy.json # Pyre type checker .pyre/ + +**/.DS_Store diff --git a/des_y6utils/mdet.py b/des_y6utils/mdet.py index c8cf956..55e1a47 100644 --- a/des_y6utils/mdet.py +++ b/des_y6utils/mdet.py @@ -7,6 +7,65 @@ import healsparse +DEFAULT_HSP_MASK = "y6-combined-hleda-gaiafull-des-stars-hsmap131k-mdet-extra-masks-v2.hsp" # noqa: E501 + + +def add_extinction_correction_columns( + data, dust_map_filename="SFD_dust_4096.hsp" +): + """This function adds dereddened columns to the data and copies the existing columns + to columns with '_nodered' appended as a suffix. + + Parameters + ---------- + data : np.ndarray + A structured array of data to be dereddened. + dust_map_filename : str + The path/name of the dust map. Assumed to be in HealSparse format. + + Returns + ------- + new_data : np.ndarray + The new, dereddened data. + """ + + bands = ['g', 'r', 'i', 'z'] + + # compute dereddened fluxes and + # replace the entries of pgauss_flux_band_* with the dereddened fluxes. + if os.path.exists(dust_map_filename): + dustmap = _read_hsp_file_local(dust_map_filename) + else: + dustmap = _read_hsp_file(dust_map_filename) + + dered = dustmap.get_values_pos(data["ra"], data["dec"]) + flux_og = [] + for ii, b in enumerate(bands): + dered_fac = _compute_dered_flux_fac(ii, dered) + flux_og.append(data["pgauss_band_flux_" + b].copy()) + + data["pgauss_band_flux_" + b] *= dered_fac + data["pgauss_band_flux_err_" + b] *= dered_fac + + # make _nodered array with pgauss_band_flux_* entries, and add them to fits. + new_dt = np.dtype( + data.dtype.descr + + [ + ('pgauss_band_flux_g_nodered', 'f4'), + ('pgauss_band_flux_r_nodered', 'f4'), + ('pgauss_band_flux_i_nodered', 'f4'), + ('pgauss_band_flux_z_nodered', 'f4'), + ] + ) + new_data = np.zeros(data.shape, dtype=new_dt) + for col in data.dtype.names: + new_data[col] = data[col] + for ii, b in enumerate(bands): + new_data['pgauss_band_flux_' + b + '_nodered'] = flux_og[ii] + + return new_data + + def make_mdet_cuts(data, version, verbose=False): """A function to make the standard metadetection cuts. @@ -40,6 +99,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 +113,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=DEFAULT_HSP_MASK, max_t=100.0, ): """ @@ -61,7 +122,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, @@ -74,9 +135,9 @@ def _make_mdet_cuts_gauss( # apply the mask if os.path.exists(mask_name): - hmap = _read_hsp_mask_local(mask_name) + hmap = _read_hsp_file_local(mask_name) else: - hmap = _read_hsp_mask(mask_name) + hmap = _read_hsp_file(mask_name) in_footprint = hmap.get_values_pos( data["ra"], data["dec"], valid_mask=True ) @@ -87,6 +148,123 @@ def _make_mdet_cuts_gauss( 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, + ) + + # apply the mask + hmap = _read_hsp_file(DEFAULT_HSP_MASK) + 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 _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.7. + - 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.7) # used to eliminate photo-z outliers motivated by cosmos2020 + # the cut is 24.5 in cosmos mags with + 0.2 mags to correct + # for the pgauss aperture + & (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) + ) + + 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 + + if verbose: + print("did mdet cuts", np.sum(msk)) + + return msk + + def _make_mdet_cuts_wmom( data, verbose=False, @@ -115,9 +293,9 @@ def _make_mdet_cuts_wmom( # apply the mask if os.path.exists(mask_name): - hmap = _read_hsp_mask_local(mask_name) + hmap = _read_hsp_file_local(mask_name) else: - hmap = _read_hsp_mask(mask_name) + hmap = _read_hsp_file(mask_name) in_footprint = hmap.get_values_pos( data["ra"], data["dec"], valid_mask=True ) @@ -141,7 +319,7 @@ def _make_mdet_cuts_v1(d, verbose=False): ) # apply the mask - hmap = _read_hsp_mask("y6-combined-hleda-gaiafull-hsmap16384-nomdet.fits") + hmap = _read_hsp_file("y6-combined-hleda-gaiafull-hsmap16384-nomdet.fits") in_footprint = hmap.get_values_pos(d["ra"], d["dec"], valid_mask=True) msk &= in_footprint if verbose: @@ -163,7 +341,7 @@ def _make_mdet_cuts_v2(d, verbose=False): ) # apply the mask - hmap = _read_hsp_mask( + hmap = _read_hsp_file( "y6-combined-hleda-gaiafull-des-stars-hsmap16384-nomdet-v2.fits" ) in_footprint = hmap.get_values_pos(d["ra"], d["dec"], valid_mask=True) @@ -188,7 +366,7 @@ def _make_mdet_cuts_v3(d, verbose=False): ) # apply the mask - hmap = _read_hsp_mask( + hmap = _read_hsp_file( "y6-combined-hleda-gaiafull-des-stars-hsmap16384-nomdet-v2.fits" ) in_footprint = hmap.get_values_pos(d["ra"], d["dec"], valid_mask=True) @@ -212,13 +390,13 @@ 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 # apply the mask - hmap = _read_hsp_mask( + hmap = _read_hsp_file( "y6-combined-hleda-gaiafull-des-stars-hsmap16384-nomdet-v2.fits" ) in_footprint = hmap.get_values_pos(d["ra"], d["dec"], valid_mask=True) @@ -242,13 +420,13 @@ 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 # apply the mask - hmap = _read_hsp_mask( + hmap = _read_hsp_file( "y6-combined-hleda-gaiafull-des-stars-hsmap16384-nomdet-v3.fits" ) in_footprint = hmap.get_values_pos(d["ra"], d["dec"], valid_mask=True) @@ -291,18 +469,70 @@ def _compute_asinh_mags(flux, i): return mag +def _compute_dered_flux_fac(i, ebv_map_val): + """This function computes the dereddening factor for the flux. + + Eli says we cannot apply the correction to the asinh mag. + + You use like this: + + dered_flux = flux * _compute_dered_flux_fac(i, map) + + Parameters + ---------- + i : int + The index of the band in griz (i.e., 0 for g, 1 for r, 2 for i, 3 for z). + ebv_map_val : float or np.ndarray + The E(B-V) value(s) from the dust map. + + Returns + ------- + dered_flux : float or np.ndarray + The dereddened flux. + """ + + dered_fac = [3.186, 2.140, 1.196, 1.048] + dmag = dered_fac[i] * ebv_map_val + return 10**(dmag / 2.5) + + +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) +def _read_hsp_file(fname): + mpth = _get_hsp_file_path(fname) return healsparse.HealSparseMap.read(mpth) @lru_cache -def _read_hsp_mask_local(fname): +def _read_hsp_file_local(fname): return healsparse.HealSparseMap.read(fname) -def _get_mask_path(fname): +def _get_hsp_file_path(fname): # get or make meds dir meds_dir = os.environ.get("MEDS_DIR", None) if meds_dir is None: @@ -325,25 +555,38 @@ def _download_fname_from_bnl(fpth): wget_res = subprocess.run("which wget", shell=True, capture_output=True) curl_res = subprocess.run("which curl", shell=True, capture_output=True) - bnl = "https://www.cosmo.bnl.gov/www/esheldon/data/y6-healsparse" + base_urls_to_try = [ + "https://www.cosmo.bnl.gov/www/esheldon/data/y6-healsparse", + ] + if wget_res.returncode == 0: - subprocess.run( - "cd %s && wget %s/%s" % ( - fdir, bnl, fname, - ), - shell=True, - check=True, - capture_output=True, - ) + for bnl in base_urls_to_try: + ret = subprocess.run( + "cd %s && wget %s/%s" % ( + fdir, bnl, fname, + ), + shell=True, + capture_output=True, + ) + if ret.returncode == 0: + break + + ret.check_returncode() + elif curl_res.returncode == 0: - subprocess.run( - "cd %s && curl -L %s/%s --output %s" % ( - fdir, bnl, fname, fname, - ), - shell=True, - check=True, - capture_output=True, - ) + for bnl in base_urls_to_try: + ret = subprocess.run( + "cd %s && curl -L %s/%s --output %s" % ( + fdir, bnl, fname, fname, + ), + shell=True, + capture_output=True, + ) + if ret.returncode == 0: + break + + ret.check_returncode() + else: raise RuntimeError( "Could not download mask '%s' from BNL due " @@ -431,11 +674,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 +757,11 @@ 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"]) + (n_terr * d["gauss_T_err"] / d["gauss_psf_T"]) ) ) & ((d["gauss_T_ratio"] * d["gauss_psf_T"]) < max_t) diff --git a/des_y6utils/tests/test_asinh_mag.py b/des_y6utils/tests/test_asinh_mag.py new file mode 100644 index 0000000..22cade2 --- /dev/null +++ b/des_y6utils/tests/test_asinh_mag.py @@ -0,0 +1,23 @@ +from ..mdet import _compute_dered_flux_fac, _compute_asinh_flux, _compute_asinh_mags +import math + + +def test_compute_asinh_flux(): + # this test ensures the function computes magnitudes and fluxes correctly. + flux_input = 10000 + mag_g = _compute_asinh_mags(flux_input, 0) + flux_g = _compute_asinh_flux(mag_g, 0) + + assert math.isclose(flux_input, flux_g, rel_tol=1e-9) + + +def test_compute_dered_flux(): + flux_input = 10000 + dered_flux = flux_input * _compute_dered_flux_fac(0, 2.5 / 3.186) + + assert math.isclose(flux_input*10, dered_flux, rel_tol=1e-9) + + mag_g = _compute_asinh_mags(flux_input, 0) + dered_mag_g = _compute_asinh_mags(dered_flux, 0) + + assert math.isclose(mag_g - 2.5, dered_mag_g, rel_tol=1e-8)