diff --git a/doc/changes.rst b/doc/changes.rst index 34b74376d..57c00b4bd 100644 --- a/doc/changes.rst +++ b/doc/changes.rst @@ -5,7 +5,7 @@ desitarget Change Log 0.44.1 (unreleased) ------------------- -* No changes yet. +* Add RF files for DR9 and correct threshold for each photometric footprint [`PR #648`]. 0.44.0 (2020-11-12) ------------------- diff --git a/py/desitarget/cmx/cmx_cuts.py b/py/desitarget/cmx/cmx_cuts.py index be5e644ab..ff7bfb487 100644 --- a/py/desitarget/cmx/cmx_cuts.py +++ b/py/desitarget/cmx/cmx_cuts.py @@ -2158,16 +2158,13 @@ def apply_cuts(objects, cmxdir=None, noqso=False): photsys_north, photsys_south, obs_rflux, gflux, rflux, zflux, \ w1flux, w2flux, gfiberflux, rfiberflux, zfiberflux, \ - objtype, release, gfluxivar, rfluxivar, zfluxivar, \ + objtype, release, ra, dec, gfluxivar, rfluxivar, zfluxivar, \ gnobs, rnobs, znobs, gfracflux, rfracflux, zfracflux, \ gfracmasked, rfracmasked, zfracmasked, \ gfracin, rfracin, zfracin, gallmask, rallmask, zallmask, \ - gsnr, rsnr, zsnr, w1snr, w2snr, dchisq, deltaChi2, maskbits, refcat = \ + gsnr, rsnr, zsnr, w1snr, w2snr, dchisq, deltaChi2, maskbits, refcat = \ _prepare_optical_wise(objects) - # ADM in addition, cmx needs ra and dec. - ra, dec = objects["RA"], objects["DEC"] - # ADM Currently only coded for objects with Gaia matches # ADM (e.g. DR6 or above). Fail for earlier Data Releases. if np.any(release < 6000): diff --git a/py/desitarget/cuts.py b/py/desitarget/cuts.py index 6127017d1..c3f2b2c62 100644 --- a/py/desitarget/cuts.py +++ b/py/desitarget/cuts.py @@ -1245,7 +1245,7 @@ def isQSO_colors(gflux=None, rflux=None, zflux=None, w1flux=None, w2flux=None, def isQSO_randomforest(gflux=None, rflux=None, zflux=None, maskbits=None, w1flux=None, w2flux=None, objtype=None, release=None, gnobs=None, rnobs=None, znobs=None, deltaChi2=None, - primary=None, south=True, return_probs=False): + primary=None, ra=None, dec=None, south=True, return_probs=False): """Define QSO targets from a Random Forest. Returns a boolean array. Parameters @@ -1343,6 +1343,8 @@ def isQSO_randomforest(gflux=None, rflux=None, zflux=None, maskbits=None, rf_DR7_HighZ_fileName = pathToRF + '/rf_model_dr7_HighZ.npz' rf_DR8_fileName = pathToRF + '/rf_model_dr8.npz' rf_DR8_HighZ_fileName = pathToRF + '/rf_model_dr8_HighZ.npz' + rf_DR9_fileName = pathToRF + '/rf_model_dr9.npz' + rf_DR9_HighZ_fileName = pathToRF + '/rf_model_dr9_HighZ.npz' tmpReleaseOK = releaseReduced < 5000 if np.any(tmpReleaseOK): @@ -1391,7 +1393,7 @@ def isQSO_randomforest(gflux=None, rflux=None, zflux=None, maskbits=None, qsohiz[colorsReducedIndex[tmpReleaseOK]] = \ (tmp_rf_HighZ_proba >= pcut_HighZ) - tmpReleaseOK = releaseReduced >= 8000 + tmpReleaseOK = (releaseReduced >= 8000) & (releaseReduced < 9000) if np.any(tmpReleaseOK): # rf initialization - colors data duplicated within "myRF" rf = myRF(colorsReduced[tmpReleaseOK], pathToRF, @@ -1420,6 +1422,51 @@ def isQSO_randomforest(gflux=None, rflux=None, zflux=None, maskbits=None, # ADM populate a mask specific to the "HighZ" selection. pqsohiz[colorsReducedIndex[tmpReleaseOK]] = tmp_rf_HighZ_proba + tmpReleaseOK = releaseReduced >= 9000 + if np.any(tmpReleaseOK): + # rf initialization - colors data duplicated within "myRF" + rf = myRF(colorsReduced[tmpReleaseOK], pathToRF, + numberOfTrees=500, version=2) + rf_HighZ = myRF(colorsReduced[tmpReleaseOK], pathToRF, + numberOfTrees=500, version=2) + # rf loading + rf.loadForest(rf_DR9_fileName) + rf_HighZ.loadForest(rf_DR9_HighZ_fileName) + # Compute rf probabilities + tmp_rf_proba = rf.predict_proba() + tmp_rf_HighZ_proba = rf_HighZ.predict_proba() + # Compute optimized proba cut + tmp_r_Reduced = r_Reduced[tmpReleaseOK] + if not south: + #threshold selection for North footprint + pcut = 0.857 - 0.03*np.tanh(tmp_r_Reduced - 20.5) + pcut_HighZ = 0.7 + else: + pcut = np.ones(tmp_rf_proba.size) + pcut_HighZ = np.ones(tmp_rf_HighZ_proba.size) + is_des = (gnobs[preSelection][tmpReleaseOK] > 4) &\ + (rnobs[preSelection][tmpReleaseOK] > 4) &\ + (znobs[preSelection][tmpReleaseOK] > 4) &\ + ((ra[preSelection][tmpReleaseOK] >= 320) | (ra[preSelection][tmpReleaseOK] <= 100)) &\ + (dec[preSelection][tmpReleaseOK] <= 10) + #threshold selection for Des footprint + pcut[is_des] = 0.75 - 0.05*np.tanh(tmp_r_Reduced[is_des] - 20.5) + pcut_HighZ[is_des] = 0.50 + #threshold selection for South footprint + pcut[~is_des] = 0.85 - 0.04*np.tanh(tmp_r_Reduced[~is_des] - 20.5) + pcut_HighZ[~is_des] = 0.65 + + # Add rf proba test result to "qso" mask + qso[colorsReducedIndex[tmpReleaseOK]] = \ + (tmp_rf_proba >= pcut) | (tmp_rf_HighZ_proba >= pcut_HighZ) + # ADM populate a mask specific to the "HighZ" selection. + qsohiz[colorsReducedIndex[tmpReleaseOK]] = \ + (tmp_rf_HighZ_proba >= pcut_HighZ) + # ADM store the probabilities in case they need returned. + pqso[colorsReducedIndex[tmpReleaseOK]] = tmp_rf_proba + # ADM populate a mask specific to the "HighZ" selection. + pqsohiz[colorsReducedIndex[tmpReleaseOK]] = tmp_rf_HighZ_proba + # In case of call for a single object passed to the function with # scalar arguments. Return "numpy.bool_" instead of "~numpy.ndarray". if nbEntries == 1: @@ -1576,6 +1623,9 @@ def _prepare_optical_wise(objects, mask=True): objtype = objects['TYPE'] release = objects['RELEASE'] + ra = objects['RA'] + dec = objects['DEC'] + gfluxivar = objects['FLUX_IVAR_G'] rfluxivar = objects['FLUX_IVAR_R'] zfluxivar = objects['FLUX_IVAR_Z'] @@ -1629,7 +1679,7 @@ def _prepare_optical_wise(objects, mask=True): return (photsys_north, photsys_south, obs_rflux, gflux, rflux, zflux, w1flux, w2flux, gfiberflux, rfiberflux, zfiberflux, - objtype, release, gfluxivar, rfluxivar, zfluxivar, + objtype, release, ra, dec, gfluxivar, rfluxivar, zfluxivar, gnobs, rnobs, znobs, gfracflux, rfracflux, zfracflux, gfracmasked, rfracmasked, zfracmasked, gfracin, rfracin, zfracin, gallmask, rallmask, zallmask, @@ -1757,7 +1807,7 @@ def unextinct_fluxes(objects): def set_target_bits(photsys_north, photsys_south, obs_rflux, gflux, rflux, zflux, w1flux, w2flux, gfiberflux, rfiberflux, zfiberflux, - objtype, release, gfluxivar, rfluxivar, zfluxivar, + objtype, release, ra, dec, gfluxivar, rfluxivar, zfluxivar, gnobs, rnobs, znobs, gfracflux, rfracflux, zfracflux, gfracmasked, rfracmasked, zfracmasked, gfracin, rfracin, zfracin, gallmask, rallmask, zallmask, @@ -1842,6 +1892,8 @@ def set_target_bits(photsys_north, photsys_south, obs_rflux, resolvetargs : :class:`boolean`, optional, defaults to ``True`` If ``True``, if only northern (southern) sources are passed then only apply the northern (southern) cuts to those sources. + ra, dec : :class:`~numpy.ndarray` + The Ra, Dec position of objects Returns ------- @@ -1924,7 +1976,7 @@ def set_target_bits(photsys_north, photsys_south, obs_rflux, primary=primary, zflux=zflux, rflux=rflux, gflux=gflux, w1flux=w1flux, w2flux=w2flux, deltaChi2=deltaChi2, maskbits=maskbits, gnobs=gnobs, rnobs=rnobs, znobs=znobs, - objtype=objtype, release=release, south=south + objtype=objtype, release=release, ra=ra, dec=dec, south=south ) else: raise ValueError('Unknown qso_selection {}; valid options are {}'.format( @@ -2298,11 +2350,11 @@ def apply_cuts(objects, qso_selection='randomforest', gaiamatch=False, # ADM process the Legacy Surveys columns for Target Selection. photsys_north, photsys_south, obs_rflux, gflux, rflux, zflux, \ w1flux, w2flux, gfiberflux, rfiberflux, zfiberflux, \ - objtype, release, gfluxivar, rfluxivar, zfluxivar, \ + objtype, release, ra, dec, gfluxivar, rfluxivar, zfluxivar, \ gnobs, rnobs, znobs, gfracflux, rfracflux, zfracflux, \ gfracmasked, rfracmasked, zfracmasked, \ gfracin, rfracin, zfracin, gallmask, rallmask, zallmask, \ - gsnr, rsnr, zsnr, w1snr, w2snr, dchisq, deltaChi2, maskbits, refcat = \ + gsnr, rsnr, zsnr, w1snr, w2snr, dchisq, deltaChi2, maskbits, refcat = \ _prepare_optical_wise(objects, mask=mask) # Process the Gaia inputs for target selection. @@ -2332,7 +2384,7 @@ def apply_cuts(objects, qso_selection='randomforest', gaiamatch=False, photsys_north, photsys_south, obs_rflux, gflux, rflux, zflux, w1flux, w2flux, gfiberflux, rfiberflux, zfiberflux, - objtype, release, gfluxivar, rfluxivar, zfluxivar, + objtype, release, ra, dec, gfluxivar, rfluxivar, zfluxivar, gnobs, rnobs, znobs, gfracflux, rfracflux, zfracflux, gfracmasked, rfracmasked, zfracmasked, gfracin, rfracin, zfracin, gallmask, rallmask, zallmask, @@ -2341,8 +2393,8 @@ def apply_cuts(objects, qso_selection='randomforest', gaiamatch=False, gaiagmag, gaiabmag, gaiarmag, gaiaaen, gaiadupsource, gaiaparamssolved, gaiabprpfactor, gaiasigma5dmax, galb, tcnames, qso_optical_cuts, qso_selection, - maskbits, Grr, refcat, primary, resolvetargs=resolvetargs - ) + maskbits, Grr, refcat, primary, resolvetargs=resolvetargs, + ) return desi_target, bgs_target, mws_target diff --git a/py/desitarget/data/rf_model_dr9.npz b/py/desitarget/data/rf_model_dr9.npz new file mode 100644 index 000000000..fe766ea4b Binary files /dev/null and b/py/desitarget/data/rf_model_dr9.npz differ diff --git a/py/desitarget/data/rf_model_dr9_HighZ.npz b/py/desitarget/data/rf_model_dr9_HighZ.npz new file mode 100644 index 000000000..037b75d8e Binary files /dev/null and b/py/desitarget/data/rf_model_dr9_HighZ.npz differ diff --git a/py/desitarget/randoms.py b/py/desitarget/randoms.py index c527f8229..1125f98ec 100644 --- a/py/desitarget/randoms.py +++ b/py/desitarget/randoms.py @@ -1131,7 +1131,7 @@ def pixmap(randoms, targets, rand_density, nside=256, gaialoc=None): mbcomb = [] mbstore = [] for mb in [[10, 12, 13], - [1, 10, 12, 13], + [1, 5, 6, 7, 10, 12, 13], [1, 5, 6, 7, 11, 12, 13]]: bitint = np.sum(2**np.array(mb)) mbcomb.append(bitint) diff --git a/py/desitarget/sv1/sv1_cuts.py b/py/desitarget/sv1/sv1_cuts.py index b6bdadd48..7864499fb 100644 --- a/py/desitarget/sv1/sv1_cuts.py +++ b/py/desitarget/sv1/sv1_cuts.py @@ -1429,7 +1429,7 @@ def isMWS_WD(primary=None, gaia=None, galb=None, astrometricexcessnoise=None, def set_target_bits(photsys_north, photsys_south, obs_rflux, gflux, rflux, zflux, w1flux, w2flux, gfiberflux, rfiberflux, zfiberflux, - objtype, release, gfluxivar, rfluxivar, zfluxivar, + objtype, release, ra, dec, gfluxivar, rfluxivar, zfluxivar, gnobs, rnobs, znobs, gfracflux, rfracflux, zfracflux, gfracmasked, rfracmasked, zfracmasked, gfracin, rfracin, zfracin, gallmask, rallmask, zallmask,