Skip to content

Commit

Permalink
Merge pull request #10 from des-science/more-cuts
Browse files Browse the repository at this point in the history 
ENH add random cut function
  • Loading branch information
@beckermr
beckermr authored Nov 4, 2022
2 parents 25009e9 + 0ccb8c9 commit 808129a
Showing 1 changed file with 150 additions and 78 deletions.
228 changes: 150 additions & 78 deletions des_y6utils/mdet.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,93 +38,58 @@ def make_mdet_cuts(data, version, verbose=False):
raise ValueError("the mdet cut version '%r' is not recognized!" % version)


def _make_mdet_cuts_raw_v12(d, verbose=False):
"""The raw v1 cuts come from extensive analysis over summer 2022. They
reflect a first-pass at a consensus set of cuts.
components/comments
- We use wmom for shears and pgauss for fluxes. The weighted moments
appear to have less noise in terms of precision on the mean shear.
The pgauss fluxes remove the effects of the PSF on the fluxes to first
order and match the aperture for colors across bands.
- 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.
- We use a signal to noise cut of 10 in wmom.
- We use an error dependent cut in the size ratio. This ensures that as the noise
increases we move the size cut higher to eliminate stellar contamination.
- We use "gold-inspired" cuts for crazy colors.
- We cut especially faint objects in each band.
def _make_mdet_cuts(
data,
verbose=False,
min_s2n=10,
n_terr=3,
min_t_ratio=1.2,
max_mfrac=0.1,
max_s2n=np.inf,
mask_name="y6-combined-hleda-gaiafull-des-stars-hsmap16384-nomdet-v2.fits"
):
"""
Returns
-------
msk : np.ndarray of bool
A boolean array with the cuts. To cut the data, use `data[msk]`.
"""
msk = _make_mdet_cuts_raw_v12(
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,
)

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",
"wmom_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["wmom_s2n"] > 10)
& (d["mfrac"] < 0.1)
& (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.9 - 2.8*d["pgauss_T_err"]))
& (
d["wmom_T_ratio"] >= np.maximum(
1.2,
(1.0 + 3.0*d["wmom_T_err"]/d["wmom_psf_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 mdet cuts", np.sum(msk))
print("did mask cuts", np.sum(msk))

return msk


def _make_mdet_cuts_v1(d, verbose=False):

msk = _make_mdet_cuts_raw_v12(d, verbose=verbose)
msk = _make_mdet_cuts_raw_v12(
d,
verbose=verbose,
min_s2n=10,
n_terr=3,
min_t_ratio=1.2,
max_mfrac=0.1,
max_s2n=np.inf,
)

# apply the mask
hmap = _read_hsp_mask("y6-combined-hleda-gaiafull-hsmap16384-nomdet.fits")
Expand All @@ -138,7 +103,15 @@ def _make_mdet_cuts_v1(d, verbose=False):

def _make_mdet_cuts_v2(d, verbose=False):

msk = _make_mdet_cuts_raw_v12(d, verbose=verbose)
msk = _make_mdet_cuts_raw_v12(
d,
verbose=verbose,
min_s2n=10,
n_terr=3,
min_t_ratio=1.2,
max_mfrac=0.1,
max_s2n=np.inf,
)

# apply the mask
hmap = _read_hsp_mask(
Expand Down Expand Up @@ -190,6 +163,11 @@ def _read_hsp_mask(fname):
return healsparse.HealSparseMap.read(mpth)


@lru_cache
def _read_hsp_mask_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)
Expand Down Expand Up @@ -237,3 +215,97 @@ def _download_fname_from_bnl(fpth):
"Could not download mask '%s' from BNL due "
"to wget or curl missing!" % fname,
)


def _make_mdet_cuts_raw_v12(
d,
*,
min_s2n,
n_terr,
min_t_ratio,
max_mfrac,
max_s2n,
verbose=False,
):
"""The raw v1 cuts come from extensive analysis over summer 2022. They
reflect a first-pass at a consensus set of cuts.
components/comments
- We use wmom for shears and pgauss for fluxes. The weighted moments
appear to have less noise in terms of precision on the mean shear.
The pgauss fluxes remove the effects of the PSF on the fluxes to first
order and match the aperture for colors across bands.
- 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.
- We use a signal to noise cut of 10 in wmom.
- We use an error dependent cut in the size ratio. This ensures that as the noise
increases we move the size cut higher to eliminate stellar contamination.
- We use "gold-inspired" cuts for crazy colors.
- We cut especially faint objects in each band.
"""

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",
"wmom_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["wmom_s2n"] > min_s2n)
& (d["wmom_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.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"])
)
)
)
if verbose:
print("did mdet cuts", np.sum(msk))

return msk

0 comments on commit 808129a

Please sign in to comment.