Examples

Basic K-correction Calculation

The way kcorrect works is that the user initializes an object of the Kcorrect class for the specific set of response functions for the photometry they are interested in. This object takes a bit of time to calculate, because it precalculates the projections of the K-correction templates onto the response function on a grid of redshifts.

Given redshifts and photometry in units of maggies (and their inverse variance), the object’s fit_coeffs method fits the templates to the photometry.

These coefficients can then be used to calculate the kcorrections with the kcorrect method.

Below is an example run on a single object, for which redshift is a scalar, and maggies and ivar are 1-dimensional arrays. We can input the arrays as ndarrays or as lists.

import kcorrect.kcorrect

responses = ['sdss_u0', 'sdss_g0', 'sdss_r0', 'sdss_i0', 'sdss_z0']
kc = kcorrect.kcorrect.Kcorrect(responses=responses)

redshift = 0.0826
maggies = [27.26e-9, 73.98e-9, 132.56e-9, 198.52e-9, 238.05e-9]
ivar = [1.528e+18, 5.23e+18, 2.429e+18, 1.042e+18, 0.1653e+18]

# "coeffs" is a [5]-array with coefficients multiplying each template
coeffs = kc.fit_coeffs(redshift=redshift, maggies=maggies, ivar=ivar)

# "k" is a [5]-array with the K-corrections in magnitude units
k = kc.kcorrect(redshift=redshift, coeffs=coeffs)

For multiple objects (and optionally for single objects), redshift is a 1-dimensional array, whereas maggies and ivar are 2-dimensional arrays. Again, these can be input as lists too, as shown below.

import kcorrect.kcorrect

responses = ['sdss_u0', 'sdss_g0', 'sdss_r0', 'sdss_i0']
kc = kcorrect.kcorrect.Kcorrect(responses=responses)

redshift = [0.0826, 0.111]
maggies = [[27.26e-9, 73.98e-9, 132.56e-9, 198.52e-9],
                 [4.022e-9, 29.36e-9, 98.230e-9, 155.63e-9]]
ivar = [[1.528e+18, 5.23e+18, 2.429e+18, 1.042e+18],
              [2.360e+18, 9.87e+18, 3.006e+18, 1.122e+18]]


# "coeffs" is a [2,5]-array coefficients multiplying each template
coeffs = kc.fit_coeffs(redshift=redshift, maggies=maggies, ivar=ivar)

# "k" is a [2,4]-array with the K-corrections in magnitude units
k = kc.kcorrect(redshift=redshift, coeffs=coeffs)

Absolute magnitudes

The absmag method returns absolute magnitudes.

The code returns:

\[m_Q = m_R - {\rm DM}(z) - K_{QR}(z),\]

for those bands for which the observed maggies are positive and have a positive inverse variance. For other bands, the returned absolute magnitude is determined using the reconstructed maggies (see next section).

absmag = kc.absmag(redshift=redshift, maggies=maggies, ivar=ivar, coeffs=coeffs)

If one of the maggies is zero or negative, the corresponding absolute magnitude is returned as -9999. In the code below, the u-band absolute magnitude of the first object should have this value.

maggies = [[-3.26e-9, 73.98e-9, 132.56e-9, 198.52e-9],
           [4.022e-9, 29.36e-9, 98.230e-9, 155.63e-9]]
coeffs = kc.fit_coeffs(redshift=redshift, maggies=maggies, ivar=ivar)
absmag = kc.absmag(redshift=redshift, maggies=maggies, ivar=ivar, coeffs=coeffs)

In these cases, the method can calculate the one-sigma absolute magnitude limit based on the inverse variance, if it is not zero, by setting the limit keyword to True. This returns the one-sigma detection limit for all objects and bands (even the detected ones).

absmag, absmag_limit = kc.absmag(redshift=redshift, maggies=maggies, ivar=ivar, coeffs=coeffs, limit=True)

If one of the input inverse variances is zero or negative, the corresponding absolute magnitude and (if it is requested) the limit are returned as -9999.

ivar = [[0., 0., 2.429e+18, 1.042e+18],
                    [4.022e-9, 29.36e-9, 98.230e-9, 155.63e-9]]
coeffs = kc.fit_coeffs(redshift=redshift, maggies=maggies, ivar=ivar)
absmag, absmag_limit = kc.absmag(redshift=redshift, maggies=maggies, ivar=ivar, coeffs=coeffs, limit=True)

There is also an option to return the reconstructed absolute magnitude from the full SED fit. This will differ somewhat from the absolute magnitude generated through applying the K-correction, of course, because the SED fit is not a perfect fit to the data. It may be a useful quantity in assessing the goodness of fit (though the reconstructed maggies are a more direct comparison to the input data, see below). It also may be useful as a guess of the absolute magnitude for missing bands.

absmag, absmag_reconstruct = kc.absmag(redshift=redshift, maggies=maggies, ivar=ivar, coeffs=coeffs, reconstruct=True)

Reconstructing spectra and fluxes

The model expressed by the coefficients is a full SED. The Kcorrect object has an attribute templates that is an instance of the Template class.

Each template is in units of \({\rm ~erg} {\rm ~s}^{-1} {\rm ~cm}^{-2} {\rm ~A}^{-1} {\rm ~}M_\odot^{-1}\) as observed for a galaxy at 10pc distance. The coefficients are in units of “the solar masses corresponding to a galaxy at 10pc” (see next section).

It is straightforward to use the templates to reconstruct what that SED looks like. If you have found coefficients as in the first example above, the following code will return the SED in \({\rm ~erg} {\rm ~s}^{-1} {\rm ~cm}^{-2} {\rm ~A}^{-1}\). The templates have to be shifted to the observed frame (conserving bolometric flux) to get the correct observed spectrum.

import matplotlib.pyplot as plt
import numpy as np

wave = kc.templates.restframe_wave * (1. + redshift)
spec = coeffs.dot(kc.templates.restframe_flux) / (1. + redshift)

plt.plot(np.log10(wave), np.log10(spec))
plt.xlabel('$\\log_{10} wavelength$')
plt.ylabel('$\\log_{10} flux$ (erg s$^{-1}$ cm$^{-2}$ Ang$^{-1}$')
plt.xlim([3., 4.5])
plt.ylim([-18., -14.])
plt.show()

We can also reconstruct the fluxes from the model (of course we can, because the K-correction determination must do so!). This is useful to do to compare the best fit SED to the observations. In the case here you should find agreement within a few percent between maggies and rmaggies.

rmaggies = kc.reconstruct(redshift=redshift, coeffs=coeffs)

Derived parameters (i.e. stellar mass)

The derived method returns some derived parameters. These parameters are described in the kcorrect paper, but it is important to take them with a grain of salt.

The only one worth taking at all seriously is mremain and mtol, the surviving stellar mass and mass-to-light ratios in the stellar population fit. But this parameter is shown to disagree with other estimates, with a trend of a few tenths of dex across stellar mass (kcorrect declining). While I don’t know if any stellar mass indicator from broad band photometry is great, the one in kcorrect is particularly simple (and also doesn’t come with any error bar).

Importantly, the mass-to-light ratios are in the output bandpasses (and you have to specify band_shift if you want shifted output bandpasses).

Like the absolute magnitudes, the stellar masses use the cosmo attribute of the Kcorrect object, which by default is the Planck18 cosmology from astropy.

derived = kc.derived(redshift=redshift, coeffs=coeffs)

# This has one entry per object
stellar_mass = derived['mtol']

# This has one entry per object per output bandpass
mtol = derived['mtol']

Changing the output responses

As one gets to higher redshift, the K-corrections from a given observed frame bandpass to its rest frame counterpart become more and more dependent on the SED model being correct.

One approach to dealing with this is to define a set of output rest frame bandpasses that are shifted versions of your input bandpasses, where the shift is the typical redshift of your sample. This minimizes the internal error in your sample, at the expense of calculating quantities that are less likely to be comparable to catalogs in the literature. This option can be utilized by just specifying band_shift for the kcorrect, absmag, or derived methods when you use them:

absmag = kc.absmag(redshift=redshift, maggies=maggies, ivar=ivar, coeffs=coeffs, band_shift=0.1)

A second approach is to define a set of output rest frame bandpasses that correspond closely to the effective rest frame wavelength of the observed bandpass for galaxies at the typical observed redshift in your sample. For example, at \(z\sim 0.7\) the observed SDSS \(r\), \(i\), and \(z\) bands are close in wavelength to the rest frame \(U\), \(B\), and \(V\) bands. So if we observed galaxies in the SDSS bands at those redshifts, we could find the K-corrections to \(UBV\). To do so, we have to instantiate a Kcorrect object that can do so using the response_map and response_out arguments:

import kcorrect.kcorrect

responses_in = ['sdss_u0', 'sdss_g0', 'sdss_r0', 'sdss_i0', 'sdss_z0']
responses_out = ['bessell_U', 'bessell_B', 'bessell_V']
responses_map = ['sdss_r0', 'sdss_i0', 'sdss_z0']
kc = kcorrect.kcorrect.Kcorrect(responses=responses_in,
                                responses_out=responses_out,
                                responses_map=responses_map)

# These are the ugriz observations (made up!)
redshift = 0.72
maggies = [27.26e-9, 73.98e-9, 132.56e-9, 198.52e-9, 238.05e-9]
ivar = [1.528e+18, 5.23e+18, 2.429e+18, 1.042e+18, 0.1653e+18]

# "coeffs" is a [5]-array with coefficients multiplying each template
coeffs = kc.fit_coeffs(redshift=redshift, maggies=maggies, ivar=ivar)

# "k" is a [3]-array with the K-corrections in magnitude units,
# from riz to UBV.
k = kc.kcorrect(redshift=redshift, coeffs=coeffs)

# "absmag" is also a [3]-array resuling from applying K-corrections
# and the distance modulus
absmag = kc.absmag(redshift=redshift, maggies=maggies, ivar=ivar, coeffs=coeffs)

Monte Carlo Error Estimates

Here we describe how to get K-correction, absolute magnitudes, and derived quantity errors using a Monte Carlo process built into kcorrect.

The error in the K-correction and/or the resoluting absolute magnitude is a problematic quantity. There is no perfect way to calculate it, because even when the statistical errors are small, there are systematic errors to contend with, having to do with the faithfulness and completeness of the template set, the accuracy of the filter curves, and other issues.

For purely statistical errors, in principle if all of the best-fit coefficients are non-zero, we could use the Hessian at the best-fit position to estimate the errors in the coefficients. We could use a matrix of derivatives (i.e. the Jacobian) between the coefficients and the model fluxes to get errors in the model fluxes, and propagate errors to K-corrections and absolute magnitudes. However, this procedure would not work if any of the non-negative coefficients had been pinned to zero, and would also be rather complicated.

Within kcorrect, we give the user a simple (though expensive) way to calculate errors through a Monte Carlo process. The fit_coeffs method can generate a Monte Carlo sample using the errors in the input maggies (assuming Gaussian errors). It will return the Monte Carlo maggies and their associated best-fit coefficients. Subsequently, calls to the absmag_mc and derived_mc methods will return absolute magnitudes and derived parameters for the Monte Carlo samples. The distribution of these quantities then yields the error distribution around the best-fit value.

The example below shows the use of this technique. It is invoked by providing fit_coeffs with a non-zero input parameter mc, containing the number of Monte Carlos to perform. Then after calling absmag_mc, we can look at the standard deviation of the absolute magnitudes to get an error estimate. In this case, the errors are dominated by the statistical error we input, but try decreasing the u-band inverse variance to get a sense of how the template fitting uncertainty propagates into the errors.

Note that in calculating the standard deviation, we have to account for the possibility that the maggies fluctuate negative, in which case the absolute magnitude is not defined. This is an annoying property of absolute magnitudes, and to properly deal with this you have to perform a more careful analysis of the Monte Carlo distribution. In these circumstances you may also want to look at the absmag_limit ouput of the absmag method.

import kcorrect.kcorrect

responses = ['sdss_u0', 'sdss_g0', 'sdss_r0', 'sdss_i0']
kc = kcorrect.kcorrect.Kcorrect(responses=responses)

redshift = [0.0826, 0.111]
maggies = [[27.26e-9, 73.98e-9, 132.56e-9, 198.52e-9],
           [4.022e-9, 29.36e-9, 98.230e-9, 155.63e-9]]
ivar = [[1.528e+18, 5.23e+18, 2.429e+18, 1.042e+18],
        [2.360e+18, 9.87e+18, 3.006e+18, 1.122e+18]]

# "coeffs" is a [2, 5]-array of coefficients
# "coeffs_mc" is a [2, 5, 100]-array of coefficients for each Monte Carlo
# "maggies_mc" is a [2, 4, 100]-array of the Monte Carloed maggies used
coeffs, coeffs_mc, maggies_mc = kc.fit_coeffs(redshift=redshift, maggies=maggies, ivar=ivar, mc=100)

# "absmag" is a [2, 4]-array of the absolute magnitudes
absmag = kc.absmag(redshift=redshift, maggies=maggies, ivar=ivar, coeffs=coeffs)

# "absmag_mc" is a [2, 4, 100]-array of the Monte Carloed absolute magnitudes
absmag_mc = kc.absmag_mc(redshift=redshift, maggies_mc=maggies_mc, ivar=ivar, coeffs_mc=coeffs_mc)

# The standard deviation over the Monte Carlos is an estimate of errors
absmag_err = absmag_mc.std(axis=2, where=(absmag_mc != -9999.))