Calculating single-metallicity SEDs for the Roman/Rubin 2023 Diffsky Mock

This notebook illustrates how to calculate SEDs of diffsky galaxies, including the SED of the entire galaxy, and also its decomposition into contributions from the disk, bulge, and star-forming knots.

First we’ll download a very small dataset that stores a downsampling of data from a single healpixel of the roman_rubin_2023 mock.

! 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:

2 2134k 2 57055 0 0 131k 0 0:00:16 –:–:– 0:00:16 131k

</pre>

2 2134k 2 57055 0 0 131k 0 0:00:16 –:–:– 0:00:16 131k

end{sphinxVerbatim}

2 2134k 2 57055 0 0 131k 0 0:00:16 –:–:– 0:00:16 131k

100 2134k 100 2134k 0 0 2798k 0 –:–:– –:–:– –:–:– 2797k

</pre>

100 2134k 100 2134k 0 0 2798k 0 –:–:– –:–:– –:–:– 2797k

end{sphinxVerbatim}

100 2134k 100 2134k 0 0 2798k 0 –:–:– –:–:– –:–:– 2797k

In this next cell we’ll download the template SEDs of the simple stellar populations used to compute the galaxy SEDs.

! 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 5813k 0 –:–:– –:–:– –:–:– 5809k

</pre>

100 5060k 100 5060k 0 0 5813k 0 –:–:– –:–:– –:–:– 5809k

end{sphinxVerbatim}

100 5060k 100 5060k 0 0 5813k 0 –:–:– –:–:– –:–:– 5809k

Retrieve model parameters

The diffsky model has a number of parameters controlling 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.

from lsstdesc_diffsky import read_diffskypop_params
all_diffskypop_params = read_diffskypop_params("roman_rubin_2023")
print(all_diffskypop_params._fields)

('lgfburst_u_params', 'burstshape_u_params', 'lgav_dust_u_params', 'delta_dust_u_params', 'funo_dust_u_params', 'lgmet_params')

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.

from lsstdesc_diffsky.io_utils import load_healpixel

fn = "diffsky.testdata.hdf5"
patlist = ('LSST', )
mock, metadata = load_healpixel(fn, patlist)

print(mock.keys())

odict_keys(['snapnum', 'diffmah_logmp_fit', 'diffmah_mah_logtc', 'diffmah_early_index', 'diffmah_late_index', 'diffstar_lgmcrit', 'diffstar_lgy_at_mcrit', 'diffstar_indx_lo', 'diffstar_indx_hi', 'diffstar_tau_dep', 'diffstar_lg_qt', 'diffstar_qlglgdt', 'diffstar_lg_drop', 'diffstar_lg_rejuv', 'fbulge_tcrit', 'fbulge_early', 'fbulge_late', 'burstshape_lgyr_peak', 'burstshape_lgyr_max', 'fknot', 'fburst', 'redshift', 'LSST_obs_g', 'LSST_obs_g_nodust', 'LSST_obs_i', 'LSST_obs_i_nodust', 'LSST_obs_r', 'LSST_obs_r_nodust', 'LSST_obs_u', 'LSST_obs_u_nodust', 'LSST_obs_y', 'LSST_obs_y_nodust', 'LSST_obs_z', 'LSST_obs_z_nodust', 'LSST_rest_G', 'LSST_rest_G_nodust', 'LSST_rest_I', 'LSST_rest_I_nodust', 'LSST_rest_R', 'LSST_rest_R_nodust', 'LSST_rest_U', 'LSST_rest_U_nodust', 'LSST_rest_Y', 'LSST_rest_Y_nodust', 'LSST_rest_Z', 'LSST_rest_Z_nodust'])

Each diffsky galaxy has its own individual 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.

from lsstdesc_diffsky.io_utils import load_diffsky_params
diffsky_param_data = load_diffsky_params(mock)

Inspect the SSP data

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:

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)

Working with customized lower-resolution SEDs

The SSP templates used in the Roman/Rubin 2023 mocks are high-resolution spectra (\(R\sim6000\) in the optical). You can speed up the galaxy SED computations if you use a lower-resolution version of these spectra. This next cell creates such a lower-res version of the SSP templates in a very crude way for demonstration purposes.

thin_ssp_wave = ssp_data.ssp_wave[::10]
thin_ssp_flux = ssp_data.ssp_flux[:, ::10]

However you may choose to thin the SSP templates, all downstream computations remain unchanged if you just pack your SSPs into an SSPData NamedTuple as shown below.

from lsstdesc_diffsky.legacy.roman_rubin_2023.dsps.data_loaders.defaults import SSPDataSingleMet
thin_ssp_data = SSPDataSingleMet(ssp_data.ssp_lg_age_gyr, thin_ssp_wave, thin_ssp_flux)

Compute the disk, bulge, and knot component SEDs of an individual galaxy

There are separate diffsky functions that can be used to compute component SEDs of an individual object, or of a population of objects at a time. Computing SEDs of a population all at once can be much more efficient than the SEDs one object at a time, but taking advantage of speedups from vectorization comes at the cost of increasing the memory footprint of the calculation.

from lsstdesc_diffsky.sed.disk_bulge_sed_kernels_singlemet import calc_rest_sed_disk_bulge_knot_singlegal
from lsstdesc_diffsky.defaults import OUTER_RIM_COSMO_PARAMS

igal = 0
args = (mock['redshift'][igal],
    diffsky_param_data.mah_params[igal],
    diffsky_param_data.ms_params[igal],
    diffsky_param_data.q_params[igal],
    diffsky_param_data.fbulge_params[igal],
    diffsky_param_data.fknot[igal],
    ssp_data,
    all_diffskypop_params,
    OUTER_RIM_COSMO_PARAMS)

disk_bulge_sed_info = calc_rest_sed_disk_bulge_knot_singlegal(*args)
print(disk_bulge_sed_info._fields)

print(disk_bulge_sed_info.rest_sed_bulge.shape)
print(disk_bulge_sed_info.rest_sed_diffuse_disk.shape)
print(disk_bulge_sed_info.rest_sed_knot.shape)

('rest_sed_bulge', 'rest_sed_diffuse_disk', 'rest_sed_knot', 'rest_sed_bulge_nodust', 'rest_sed_diffuse_disk_nodust', 'rest_sed_knot_nodust', 'mstar_bulge', 'mstar_diffuse_disk', 'mstar_knot', 'mstar_burst', 'mstar_total')
(5994,)
(5994,)
(5994,)

Compute the disk, bulge, and knot SED of an entire galaxy population at once

from lsstdesc_diffsky.sed.disk_bulge_sed_kernels_singlemet import calc_rest_sed_disk_bulge_knot_galpop
from lsstdesc_diffsky.defaults import OUTER_RIM_COSMO_PARAMS

args = (mock['redshift'],
    diffsky_param_data.mah_params,
    diffsky_param_data.ms_params,
    diffsky_param_data.q_params,
    diffsky_param_data.fbulge_params,
    diffsky_param_data.fknot,
    ssp_data,
    all_diffskypop_params,
    OUTER_RIM_COSMO_PARAMS)

disk_bulge_sed_info = calc_rest_sed_disk_bulge_knot_galpop(*args)
print(disk_bulge_sed_info._fields)

print(disk_bulge_sed_info.rest_sed_bulge.shape)
print(disk_bulge_sed_info.rest_sed_diffuse_disk.shape)
print(disk_bulge_sed_info.rest_sed_knot.shape)

('rest_sed_bulge', 'rest_sed_diffuse_disk', 'rest_sed_knot', 'rest_sed_bulge_nodust', 'rest_sed_diffuse_disk_nodust', 'rest_sed_knot_nodust', 'mstar_bulge', 'mstar_diffuse_disk', 'mstar_knot', 'mstar_burst', 'mstar_total')
(46, 5994)
(46, 5994)
(46, 5994)

Plot the disk/bulge/knot SEDs of an example galaxy

The cell below shows how to interpret the named tuple returned by the above computation. The in-panel annotation shows the fraction of stellar mass in the galaxy, which is also returned above.

from matplotlib import pyplot as plt

fig, ax = plt.subplots(1, 1)
xlim = ax.set_xlim(1e3, 1e5)
ylim = ax.set_ylim(5e-11, 2e-5)

__=ax.loglog()

igal = 5

__=ax.plot(
    ssp_data.ssp_wave,
    disk_bulge_sed_info.rest_sed_bulge[igal, :],
    label=r'${\rm bulge}$', color='red')
__=ax.plot(ssp_data.ssp_wave,
           disk_bulge_sed_info.rest_sed_diffuse_disk[igal, :],
           label=r'${\rm diffuse\ disk}$', color='green')
__=ax.plot(ssp_data.ssp_wave,
           disk_bulge_sed_info.rest_sed_knot[igal, :],
           label=r'${\rm star-forming\ knots}$', color='purple')

frac_bulge_igal = disk_bulge_sed_info.mstar_bulge[igal]/disk_bulge_sed_info.mstar_total[igal]
title = ax.set_title('B/T={0:.1f}'.format(frac_bulge_igal))
leg = ax.legend(loc='upper right')
xlabel = ax.set_xlabel(r'$\lambda\ [{\rm Angstrom}]$')
ylabel = ax.set_ylabel(r'${\rm SED\ [L_{\odot}/Hz]}$')

Now clean up the temporary files

! rm dsps_ssp_data.h5

rm: cannot remove 'dsps_ssp_data.h5': No such file or directory

! rm diffsky.testdata.hdf5