Calculating photometry for Roman/Rubin 2023 Diffsky galaxies¶
This notebook illustrates how to calculate photometry predictions from the physical parameters of each diffsky galaxy. There’s demo code showing how to compute exact photometry for a single galaxy at a time, and also a demo of calculating approximate photometry for large galaxy populations at once.
Download mock galaxy data¶
First we’ll download a very small dataset that stores a downsampling of data from a single healpixel of the roman_rubin_2023 mock.
[1]:
! curl https://portal.nersc.gov/project/hacc/aphearin/lsstdesc_diffsky_data/roman_rubin_2023_z_0_1_cutout_9043.testdata.hdf5 > diffsky.testdata.hdf5
- more-to-come:
% Total % Received % Xferd Average Speed Time Time Time Current Dload Upload Total Spent Left Speed
- 0 0 0 0 0 0 0 0 –:–:– –:–:– –:–:– 0
</pre>
- 0 0 0 0 0 0 0 0 –:–:– –:–:– –:–:– 0
end{sphinxVerbatim}
0 0 0 0 0 0 0 0 –:–:– –:–:– –:–:– 0
- more-to-come:
- 0 0 0 0 0 0 0 0 –:–:– –:–:– –:–:– 0
</pre>
- 0 0 0 0 0 0 0 0 –:–:– –:–:– –:–:– 0
end{sphinxVerbatim}
0 0 0 0 0 0 0 0 –:–:– –:–:– –:–:– 0
- 100 2134k 100 2134k 0 0 2376k 0 –:–:– –:–:– –:–:– 2377k
</pre>
- 100 2134k 100 2134k 0 0 2376k 0 –:–:– –:–:– –:–:– 2377k
end{sphinxVerbatim}
100 2134k 100 2134k 0 0 2376k 0 –:–:– –:–:– –:–:– 2377k
Load the Diffsky data from the test healpixel¶
The next cell directly reads the hdf5 file storing the mock data. This test file is formatted in the same way as the healpixels distributed on NERSC: the data is separated by the simulation snapshot, and there is a metadata column storing additional information. For demonstration purposes, we’ll just load galaxies directly from the hdf5 file with the load_diffsky_healpixel convenience function. The returned mock stores the full collection of data from the snapshots, concatenated and stored
as a flat ndarray for each column.
DESC users working at NERSC may instead wish to use the GCR.
[2]:
from lsstdesc_diffsky.io_utils import load_healpixel
fn = "diffsky.testdata.hdf5"
patlist = ('LSST', )
mock, metadata = load_healpixel(fn, patlist)
Download and inspect template SEDs for SSPs¶
In this next cell we’ll download the template SEDs of the simple stellar populations used to compute the galaxy SEDs.
Note: The demos below are the same single-metallicity kernels used to generate the roman_rubin_2023 mock, and so in the following cell we download single-metallicitiy SSP SEDs. Multi-metallicity kernels and SSPs are also included in the lsstdesc-diffsky library, so take care to use the single-metallicity versions as shown below.
[3]:
! curl https://portal.nersc.gov/project/hacc/aphearin/DSPS_data/ssp_data_fsps_v3.2_age.h5 > dsps_ssp_data_singlemet.h5
- more-to-come:
% Total % Received % Xferd Average Speed Time Time Time Current Dload Upload Total Spent Left Speed
- 0 0 0 0 0 0 0 0 –:–:– –:–:– –:–:– 0
</pre>
- 0 0 0 0 0 0 0 0 –:–:– –:–:– –:–:– 0
end{sphinxVerbatim}
0 0 0 0 0 0 0 0 –:–:– –:–:– –:–:– 0
- more-to-come:
- 0 0 0 0 0 0 0 0 –:–:– –:–:– –:–:– 0
</pre>
- 0 0 0 0 0 0 0 0 –:–:– –:–:– –:–:– 0
end{sphinxVerbatim}
0 0 0 0 0 0 0 0 –:–:– –:–:– –:–:– 0
- 100 5060k 100 5060k 0 0 5778k 0 –:–:– –:–:– –:–:– 5776k
</pre>
- 100 5060k 100 5060k 0 0 5778k 0 –:–:– –:–:– –:–:– 5776k
end{sphinxVerbatim}
100 5060k 100 5060k 0 0 5778k 0 –:–:– –:–:– –:–:– 5776k
The ssp_data quantity stores all the information needed from the Simple Stellar Population templates to compute our galaxy SEDs. Note that these kernels are imported a custom version of dsps located within lsstdesc_diffsky:
[4]:
from lsstdesc_diffsky.legacy.roman_rubin_2023.dsps.data_loaders.load_ssp_data import load_ssp_templates_singlemet
ssp_data = load_ssp_templates_singlemet(fn='dsps_ssp_data_singlemet.h5')
print(ssp_data._fields)
print('ssp_lg_age_gyr.shape = {}'.format(ssp_data.ssp_lg_age_gyr.shape))
print('ssp_wave.shape = {}'.format(ssp_data.ssp_wave.shape))
print('ssp_flux.shape = {}'.format(ssp_data.ssp_flux.shape))
('ssp_lg_age_gyr', 'ssp_wave', 'ssp_flux')
ssp_lg_age_gyr.shape = (107,)
ssp_wave.shape = (5994,)
ssp_flux.shape = (107, 5994)
Load transmission curves¶
The dsps library ships with a few transmission curves as a convenience for getting started. Up-to-date transmission curves for a wide range of instruments can be found from publicly available sources such as the kcorrect library and elsewhere.
Note on units: Wherever you get your transmission curves, double-check that your wavelengths are in angstroms, as these are the \(\lambda\) units used in the dsps library we’ll use to calculate photometry.
In this next cell, we’ll just use dsps to generate on-the-fly a few transmission curves that roughly approximate LSST-band photometry in ugrizy.
[5]:
from dsps.data_loaders.retrieve_fake_fsps_data import load_fake_filter_transmission_curves
wave, *trans_curves = load_fake_filter_transmission_curves()
Interpolate transmission curves to a common wavelength¶
We will vectorize our photometry computations across the set of filters, so with these calculations it’s always necessary to interpolate the collection of filters to be defined by arrays of the same length. In this demo this is trivial since the fake transmission curve generator already returns all transmission curves to be defined on the same wavelength grid, but a real transmission curve is typically defined on its own specialized grid in wavelength. The code below handles this with the
interpolate_filter_trans_curves function:
[6]:
from lsstdesc_diffsky.photometry.precompute_ssp_tables import interpolate_filter_trans_curves
wave_filters = [wave for x in trans_curves]
trans_filters = [x for x in trans_curves]
rest_filter_waves, rest_filter_trans = interpolate_filter_trans_curves(wave_filters, trans_filters)
obs_filter_waves, obs_filter_trans = interpolate_filter_trans_curves(wave_filters, trans_filters)
Retrieve Diffsky parameters of each galaxy¶
Each individual diffsky galaxy has its own parameters controlling its assembly history and SED. The load_diffsky_params function interprets the columns of the mock that store these parameters, and returns a collection of arrays that are formatted and shaped in the form expected by the function used to compute the SED of the disk, bulge, and knots.
[7]:
from lsstdesc_diffsky.io_utils import load_diffsky_params
diffsky_params = load_diffsky_params(mock)
Retrieve DiffskyPop parameters for the roman_rubin_2023 mock¶
The DiffskyPop model has a number of parameters controlling the probabilistic relationships of the galaxy–halo connection. The next cell retrieves the values of these parameters used to generate the roman_rubin_2023 mock. The returned quantity is a NamedTuple with a field name for the parameters of each DiffskyPop model ingredient.
[8]:
from lsstdesc_diffsky.param_data import read_diffskypop_params
diffskypop_params = read_diffskypop_params('roman_rubin_2023')
Calculating exact photometry for an individual galaxy¶
The code below shows how to compute the photometry of an individual diffsky galaxy. We’ll see the corresponding kernel for calculating photometry of populations of diffsky galaxies. After that, we’ll demo a much faster and more memory efficient kernel for calculating approximate photometry.
Note: The demos below are the same single-metallicity kernels used to generate the roman_rubin_2023 mock. Multi-metallicity kernels are also included in the lsstdesc-diffsky library, so take care to import from modules with _singlemet in the name as shown below.
[9]:
from lsstdesc_diffsky.defaults import OUTER_RIM_COSMO_PARAMS
[10]:
from lsstdesc_diffsky.photometry.photometry_kernels_singlemet import calc_photometry_singlegal
igal = 0
args = (
mock['redshift'][igal],
diffsky_params.mah_params[igal],
diffsky_params.ms_params[igal],
diffsky_params.q_params[igal],
ssp_data,
diffskypop_params,
rest_filter_waves,
rest_filter_trans,
obs_filter_waves,
obs_filter_trans,
OUTER_RIM_COSMO_PARAMS)
_res = calc_photometry_singlegal(*args)
rest_mags, obs_mags, rest_mags_nodust, obs_mags_nodust = _res
Calculating exact photometry for a galaxy population¶
[11]:
from lsstdesc_diffsky.photometry.photometry_kernels_singlemet import calc_photometry_galpop
args = (
mock['redshift'],
diffsky_params.mah_params,
diffsky_params.ms_params,
diffsky_params.q_params,
ssp_data,
diffskypop_params,
rest_filter_waves,
rest_filter_trans,
obs_filter_waves,
obs_filter_trans,
OUTER_RIM_COSMO_PARAMS)
_res = calc_photometry_galpop(*args)
rest_mags, obs_mags, rest_mags_nodust, obs_mags_nodust = _res
Calculating approximate photometry for a galaxy population¶
There are a couple of approximations we can make in calculating the photometry that offer orders-of-magnitude improvement in memory efficiency and runtime. These approximations give a statistically unbiased magnitude with typical scatter of \(\sigma(\delta{\rm mag})\approx0.05\), and so calculating approximate photometry can be much more practical when studying large galaxy populations. The details behind the approximations are described in these notes.
Precompute photometry of SSP template SEDs¶
In the next cell we compute the restframe magnitudes of our SSPs through each filter, and in the following cell we do the same calculation for apparent magnitudes of our SSPs on a grid of redshift.
[12]:
from lsstdesc_diffsky.photometry.precompute_ssp_tables import precompute_ssp_obsmags_on_z_table_singlemet
import numpy as np
z_table = np.linspace(0.01, mock['redshift'].max()+0.05, 51)
ssp_obsmag_table = precompute_ssp_obsmags_on_z_table_singlemet(
ssp_data.ssp_wave, ssp_data.ssp_flux,
obs_filter_waves, obs_filter_trans, z_table,
*OUTER_RIM_COSMO_PARAMS[:-1])
[13]:
from lsstdesc_diffsky.photometry.precompute_ssp_tables import precompute_ssp_restmags_singlemet
ssp_restmag_table = precompute_ssp_restmags_singlemet(
ssp_data.ssp_wave, ssp_data.ssp_flux, rest_filter_waves, rest_filter_trans)
Compute galaxy photometry¶
[14]:
from lsstdesc_diffsky.photometry.photometry_lc_interp_singlemet import get_diffsky_sed_info_singlemet
from dsps.cosmology.flat_wcdm import age_at_z0
from jax import random as jran
ran_key = jran.PRNGKey(0)
t0 = age_at_z0(*OUTER_RIM_COSMO_PARAMS[:-1])
gal_t_table = np.linspace(0.1, t0, 100)
ssp_z_table = np.linspace(mock['redshift'].min()/2, mock['redshift'].max()+0.1, 51)
args = (ran_key,
mock['redshift'],
diffsky_params.mah_params,
diffsky_params.ms_params,
diffsky_params.q_params,
ssp_z_table,
ssp_restmag_table,
ssp_obsmag_table,
ssp_data,
gal_t_table,
rest_filter_waves,
rest_filter_trans,
obs_filter_waves,
obs_filter_trans,
diffskypop_params,
OUTER_RIM_COSMO_PARAMS)
sed_info = get_diffsky_sed_info_singlemet(*args)
Interpreting the results¶
The returned sed_info is a namedtuple that stores the calculated photometry along with numerous other quantities computed along the way. * Columns gal_obsmags_dust stores the apparent magnitudes of each galaxy on each input observer-frame filter * Columns gal_restmags_dust stores the absolute magnitude in the restframe of the galaxy * The _dust and _nodust quantities refer to whether or not dust attenuation in the emitting galaxy is included in the magnitude computation
(in particular, neither quantity takes any account of dust in the Milky Way)
[15]:
print(sed_info._fields)
('gal_ssp_weights', 'gal_frac_trans_obs', 'gal_frac_trans_rest', 'gal_att_curve_params', 'gal_frac_unobs', 'gal_fburst', 'gal_burstshape_params', 'gal_frac_bulge_t_obs', 'gal_fbulge_params', 'gal_fknot', 'gal_obsmags_nodust', 'gal_restmags_nodust', 'gal_obsmags_dust', 'gal_restmags_dust')
[16]:
print(sed_info.gal_obsmags_dust.shape)
(46, 6)
[17]:
from matplotlib import pyplot as plt
fig, ax = plt.subplots(1, 1)
ylim = ax.set_ylim(-0.2, 1.25)
__=ax.scatter(
mock['redshift'],
sed_info.gal_restmags_dust[:, 1]-sed_info.gal_restmags_dust[:, 2], s=2)
xlabel = ax.set_xlabel('redshift')
ylabel = ax.set_ylabel('g-r')
[18]:
iband = 1
fig, ax = plt.subplots(1, 1)
xlim = ax.set_xlim(-0.05, 0.05)
__=ax.hist(rest_mags[:, iband] - sed_info.gal_restmags_dust[:, iband])
xlabel = ax.set_xlabel(r'$\delta m_{g}$')
[19]:
assert np.allclose(rest_mags, sed_info.gal_restmags_dust, atol=0.1)
Now clean up the temporary files¶
[20]:
! rm dsps_ssp_data_singlemet.h5
! rm diffsky.testdata.hdf5