Dark Energy Survey (DES)

Note

This guide has not been inspected or endorsed by the DES collaboration.

Here, we will show how to cross-correlate BOSS lens galaxies with shape catalogs from DES. We will work with the Y1 data release.

Downloading the Data

DES data can be downloaded from the DES data release website. The following commands should download all the necessary data. The first file downloaded is the main shape catalog, the second one includes the METACALIBRATION photo-z’s used to bin galaxies and the last one includes the MOF Monte-Carlo redshifts we will use to approximate the true redshift distribution \(n(z)\).

wget http://desdr-server.ncsa.illinois.edu/despublic/y1a1_files/shear_catalogs/mcal-y1a1-combined-riz-unblind-v4-matched.fits
wget http://desdr-server.ncsa.illinois.edu/despublic/y1a1_files/photoz_catalogs/y1a1-gold-mof-badregion_BPZ.fits
wget http://desdr-server.ncsa.illinois.edu/despublic/y1a1_files/photoz_catalogs/mcal-y1a1-combined-griz-blind-v3-matched_BPZbase.fits

Preparing the Data

The first step is to put the data into a format easily understandable by dsigma. There are a number of helper functions to make this easy. Additionally, we are going to bin source galaxies by tomographic redshift. (We are using fitsio instead of astropy.io.fits because it allows us to read only certain columns of the very large input table. If memory is not an issue, feel free to use astropy, instead.)

import fitsio
from astropy.table import Table, hstack
from dsigma.helpers import dsigma_table
from dsigma.surveys import des

table_s = []

fname_list = ['mcal-y1a1-combined-riz-unblind-v4-matched.fits',
              'y1a1-gold-mof-badregion_BPZ.fits',
              'mcal-y1a1-combined-griz-blind-v3-matched_BPZbase.fits']
columns_list = [['e1', 'e2', 'R11', 'R12', 'R21', 'R22', 'ra', 'dec',
                 'flags_select', 'flags_select_1p', 'flags_select_1m',
                 'flags_select_2p', 'flags_select_2m'], ['Z_MC'], ['MEAN_Z']]

for fname, columns in zip(fname_list, columns_list):
    table_s.append(Table(fitsio.read(fname, columns=columns), names=columns))

table_s = hstack(table_s)
table_s = dsigma_table(table_s, 'source', survey='DES')
table_s['z_bin'] = des.tomographic_redshift_bin(table_s['z'])

Selection Response

In an intermediate step, we calculate the so-called METACALIBRATION selection response. This factor takes into account how the selection flags of the METACALIBRATION shape measurements used by DES might be biased by shear itself. We need to correct such a bias in order to get unbiased shear and \(\Delta\Sigma\) measurements. See Sheldon & Huff (2017) and McClintock et al. (2018) for details.

import numpy as np

for z_bin in range(4):
    use = table_s['z_bin'] == z_bin
    R_sel = des.selection_response(table_s[use])
    print("Bin {}: R_sel = {:.1f}%".format(
        z_bin + 1, 100 * 0.5 * np.sum(np.diag(R_sel))))
    table_s['R_11'][use] += 0.5 * np.sum(np.diag(R_sel))
    table_s['R_22'][use] += 0.5 * np.sum(np.diag(R_sel))

This will give the following output.

Bin 1: R_sel = 1.2%
Bin 2: R_sel = 1.5%
Bin 3: R_sel = 1.1%
Bin 4: R_sel = 0.9%

You can see that we added this response to the total METACALIBRATION response that also takes into account how the measured shear is affected by intrinsic shear. Ideally, one would calculate the selection response for each radial bin and each specific lens sample (because this affects the source weighting). Additionally, one could also fold in how artificial shear affects the METACALIBRATION redshifts. However, as we can see above, the selection response bias is likely very small and not a strong function of redshift. Thus, we will ignore this complication here (cf. McClintock et al. 2018).

After running this selection response calculation, we are ready to drop all galaxies that are flagged for the unsheared images (and also those galaxies that fall outside the redshift bins).

table_s = table_s[(table_s['flags_select'] == 0) & (table_s['z_bin'] != -1)]

Note on the Estimator

The Dark Energy Survey uses the following estimator for \(\Delta\Sigma\) (excluding random subtraction):

\[\Delta\Sigma = \frac{ \sum_{\mathrm{ls}} w_{\mathrm{sys, l}} w_s \Sigma_{\rm crit}^{-1} (z_l, z_{s, \rm META}) e_t}{ \sum_{\mathrm{ls}} w_{\mathrm{sys, l}} \Sigma_{\rm crit}^{-1} (z_l, z_{s, \rm MOF, MC}) \Sigma_{\rm crit}^{-1} (z_l, z_{s, \rm META}) w_s R_{t,ls}} \, ,\]

where \(z_{s, \rm META}\) is the mean of the \(p(z)\) of the METACALIBRATION photometric redshift and \(z_{s, \rm MOF, MC}\) a Monte-Carlo draw of the \(p(z)\) from the multi-object fitting (MOF) photometric redshift. Finally, \(R_{t,ls}\) denotes the tangential component of the shear response of each individual lens-source pair. Initially, this might look very different from the estimators shown in the earlier section on stacking. However, we can re-arrange the terms to have the following form:

\[\begin{split}\Delta\Sigma =& \frac{ \sum_{\mathrm{ls}} w_{\mathrm{sys, l}} w_{ls}}{ \sum_{\mathrm{ls}} w_{\mathrm{sys, l}} w_{ls} R_{t,ls}} \frac{ \sum_{\mathrm{ls}} w_{\mathrm{sys, l}} w_{ls} R_{t,ls}}{ \sum_{\mathrm{ls}} w_{\mathrm{sys, l}} w_{ls} R_{t,ls} \frac{\Sigma_{\rm crit} (z_l, z_{s, \rm META})}{ \Sigma_{\rm crit} (z_l, z_{s, \rm MOF, MC})}}\\&\frac{ \sum_{\mathrm{ls}} w_{\mathrm{sys, l}} w_{ls} e_t \Sigma_{\rm crit} (z_l, z_{s, \rm META})}{ \sum_{\mathrm{ls}} w_{\mathrm{sys, l}} w_{ls}} \, ,\end{split}\]

where

\[w_{ls} = \frac{w_s}{\Sigma_{\rm crit}^2 (z_l, z_{\rm META})}\]

is the usual weight using the METACALIBRATION redshift. The above equation is much more similar to what dsigma works with. The first of the three fractions corresponds to the correction for the mean shear response. The second term is a response-weighted \(f_{\rm bias}\) factor where the Monte-Carlo draw from the MOF \(p(z)\) is used as the “true” redshift. This implies that we implicitly assume that the Monte-Carlo draws from the MOF \(p(z)\) sample the true redshift distributions \(n(z)\) of source galaxies. Finally, the third time is the usual estimate of the raw excess surface density where we used the METACALIBRATION redshift to calculate weights and critical surface densities.

dsigma implements almost the same estimator with the only difference being how the \(f_{\rm bias}\) term is calculated and applied. First, the DES estimator calculates \(f_{\rm bias}\) using only the actual lens-source pairs in each radial bin and using the actual projected response \(R_{t,ls}\) for the response weighting. On the other hand, dsigma uses all lenses and all sources, regardless of angular separation. Also, because of this, dsigma does not use the projected shear response and we will use the spherically-averaged response \(0.5 (R_{11} + R_{22})\) for each source, instead.

table_c = table_s['z', 'z_true', 'w']
table_c['w_sys'] = 0.5 * (table_s['R_11'] + table_s['R_22'])

This difference should not induce changes as long as neither the tangential response nor the source photometric redshifts are a strong function of lens-source separation. The assumption regarding the tangential response can easily be verified using dsigma in the stacking analysis below. On the other hand, source redshifts, both photometric and intrinsic, will likely change close to lenses. So additional correction may need to be applied, anyway. dsigma uses the lens-source spatial clustering to (optionally) estimate this effect but due to the complex selection function in DES, this estimator is likely biased, too. Overall, \(\Delta\Sigma\) estimates on small scales (\(r_p \lesssim 1 \mathrm{Mpc} / h\)) may currently be biased. The second difference is that the \(f_{\rm bias}\) in DES is applied to the averaged raw \(\Delta\Sigma\) estimate whereas dsigma applies it to each individual raw \(\Delta\Sigma\) estimate as a function of lens redshift, i.e. \(f_{\rm bias} = f_{\rm bias} (z_l)\). This difference should not matter since both estimators can be shown to be unbiased.

Precomputing the Signal

We will now run the computationally expensive pre-computation phase. Here, we first define the lens-source separation cuts. We require that \(z_l + 0.1 < z_s\). Afterwards, we run the actual pre-computation.

from astropy.cosmology import Planck15
from dsigma.precompute import add_maximum_lens_redshift, add_precompute_results

add_maximum_lens_redshift(table_s, dz_min=0.1)
add_maximum_lens_redshift(table_c, dz_min=0.1)

rp_bins = np.logspace(-1, 1.6, 14)
add_precompute_results(table_l, table_s, rp_bins, cosmology=Planck15,
                       comoving=True, table_c=table_c)
add_precompute_results(table_r, table_s, rp_bins, cosmology=Planck15,
                       comoving=True, table_c=table_c)

Stacking the Signal

The total galaxy-galaxy lensing signal can be obtained with the following code. It first filters out all BOSS galaxies for which we couldn’t find any source galaxy nearby. Then we divide it into different jackknife samples that we will later use to estimate uncertainties. Finally, we stack the lensing signal in 4 different BOSS redshift bins and save the data.

We choose to include all the necessary corrections factors. The tensor shear response correction and the photo-z dilution correction are the ones discussed above. Additionally, we perform a random subtraction which is highly recommended but not strictly necessary. Note that we don’t apply a boost correction since this might be biased for DES given our boost estimator.

from dsigma.jackknife import add_continous_fields, jackknife_field_centers
from dsigma.jackknife import add_jackknife_fields, jackknife_resampling
from dsigma.stacking import excess_surface_density

# Drop all lenses that did not have any nearby source.
table_l['n_s_tot'] = np.sum(table_l['sum 1'], axis=1)
table_l = table_l[table_l['n_s_tot'] > 0]

table_r['n_s_tot'] = np.sum(table_r['sum 1'], axis=1)
table_r = table_r[table_r['n_s_tot'] > 0]

centers = compute_jackknife_fields(
    table_l, 100, weights=np.sum(table_l['sum 1'], axis=1))
compute_jackknife_fields(table_r, centers)

z_bins = np.array([0.15, 0.31, 0.43, 0.54, 0.70])

for lens_bin in range(3, len(z_bins) - 1):
    mask_l = ((z_bins[lens_bin] <= table_l['z']) &
              (table_l['z'] < z_bins[lens_bin + 1]))
    mask_r = ((z_bins[lens_bin] <= table_r['z']) &
              (table_r['z'] < z_bins[lens_bin + 1]))

    kwargs = dict(return_table=True, photo_z_dilution_correction=True,
                  matrix_shear_response_correction=True,
                  random_subtraction=True, table_r=table_r[mask_r])

    result = excess_surface_density(table_l[mask_l], **kwargs)
    kwargs['return_table'] = False
    result['ds_err'] = np.sqrt(np.diag(jackknife_resampling(
        excess_surface_density, table_l[mask_l], **kwargs)))

    result.write(f'des_{lens_bin}.csv', overwrite=True)

Acknowledgements

When using the above data and algorithms, please read and follow the acknowledgement section on the DES data release site.