Calculating K-corrections using S-Lang and Sherpa

1. Define the energy grid of the model

Since we are not fitting the model to a spectrum, Sherpa does not know what energy grid to calculate the model. We therefore make use of the DATASPACE command to set one up. It is important to ensure that the energy grid covers both the observed and rest-frame energy bands and that the HISTOGRAM option is supplied, since the XSpec models are integrated over each energy bin.

For the rest-frame energy range (0.5-2.0 keV) and redshift range (0-2) considered in this example, we need to ensure that the model is evaluated over at least 0.5 to 4.0 keV. It does not matter if the energy range is larger than that, so we chose 0.01 to 10 keV, with a grid spacing of 0.01 keV:

sherpa> dataspace (0.01:10:0.01) histogram

The Sherpa utilities package contains a simple function (set_dataspace()) that calls the dataspace command and ensures that the histogram mode is set. It is of more use from S-Lang scripts than from the Sherpa command line, but would be used like:

sherpa> set_dataspace(1,0.01,10,0.01)

The first parameter (here 1), gives the source/dataset number for which the energy grid is for, and the remaining parameters are the minimum, maximum, and grid step in keV.

2. Display the model

Here we use Sherpa commands to display the model and then modify the plot using ChIPS commands (since we are in Sherpa we need to issue a REDRAW command to display the final plot):

sherpa> set_log
sherpa> sherpa.multiplot.xfilter_lo = 0.09
sherpa> lp source
sherpa> xlabel "Energy (keV)"
sherpa> ylabel "Spectrum (photon/cm^2/s)"
sherpa> title "XSMEKAL model: kT = 7 keV, Abundance = 0.3"
sherpa> xlabel size 1.4
sherpa> ylabel size 1.4
sherpa> title size 1.4
sherpa> redraw

This image  shows the resulting plot.

The set_log call ensures that both axes are plotted with a logarithmic scale. The sherpa.multiplot.xfilter_lo is used to limit all plots to have a minimum "x" value of 0.08.

Using calc_kcorr()

Now we have defined the source spectrum, we can use the calc_kcorr() S-Lang routine to calculate the k correction. So, if we want the k-correction for an object at redshift 0.4, we just say:

sherpa> calc_kcorr( 1, 0.4, 0.5, 2.0 )
0.838862

This is the k-correction for a flux measured in the observed 0.5-2.0 keV energy band to get the flux in the rest-frame 0.5-2.0 keV band for an object at a redshift of 0.4. If you want to convert to a different rest-frame energy band you just supply the band as extra parameters, as shown here for the rest-frame 0.1 to 2.4 keV energy range.

sherpa> calc_kcorr( 1, 0.4, 0.5, 2.0, 0.1, 2.4 )
1.35431

Instead of a single (ie scalar) value for the redshift, a one-dimensional array can be supplied, containing a range of redshifts. Here we calculate the correction values for the redshift range 0 to 2 with a step size of 0.1, and then plot it:

sherpa> z = [0.0:2.1:0.1]
sherpa> kc = calc_kcorr( 1, z, 0.5, 2.0 )
sherpa> clear
sherpa> () = curve(z,kc)
sherpa> symbol none
sherpa> simpleline
sherpa> xlabel "Redshift"
sherpa> ylabel "k correction"
sherpa> title "XSMEKAL model: kT = 7 keV, Abundance = 0.3"
sherpa> xlabel size 1.4
sherpa> ylabel size 1.4
sherpa> title size 1.4
sherpa> redraw

This image  shows the resulting plot.

How calc_kcorr() works

The k-correction is calculated as the ratio of the flux in the rest-frame energy band to that in the observed band (for instance see Appendix B of Jones et al. 1998, ApJ, 495, 100-114). If we denote the low and high energies of the observed band as elo and ehi respectively, then, for a source at redshift z, we need to calculate the flux between these bands:

1. observed band: elo * (1+z) to ehi * (1_z)
2. rest-frame band: elo to ehi

Use the get_eflux command:

sherpa> erange = [0.5,2.0]
sherpa> z = 0.4
sherpa> rest_flux = get_eflux(1,erange)
sherpa> obs_flux  = get_eflux(1,erange*(1+z))
sherpa> print(rest_flux)
dataset          =  1
range            =  Double_Type[2]
comp             =  NULL
value            =  5.57426e-10
units            =  counts
sherpa> print(obs_flux)
dataset          =  1
range            =  Double_Type[2]
comp             =  NULL
value            =  6.64502e-10
units            =  counts
sherpa> kcorr = rest_flux.value / obs_flux.value
sherpa> kcorr
0.838862

The get_eflux() routine returns a structure whose field value contains the source flux within the requested range - here 0.5-2.0 and 0.7-2.8 (0.5*1.4-2.0*1.4) keV for rest_flux and obs_flux respectively. Dividing the two values gives the k correction for this spectral model, redshift, and energy range.

Note that Sherpa is not able to always correctly guess the units of the models; here we have not set up any instrument model so get_eflux() is not able to work out that the source flux is in units of erg/cm²/s and so defaults to counts.

The calc_kcorr() routine contains additional code to check for errors, deal with an array of redshifts - rather than the scalar value used above - and to allow different rest-frame and observed-frame energy ranges to be used (e.g. to convert an observed flux in the 0.5-2.0 keV band to the rest-frame 0.1-2.4 keV band).

Displaying the models being evaluated

Here we show the source model over the energy ranges corresponding to the observed and rest-frame 0.5-2.0 keV band used here. It demonstrates how to use the Sherpa routines described in ahelp get to read models into S-Lang variables and then display them.

sherpa> # set up models representing the rest & observed energy ranges
sherpa> source 2 = plasma
sherpa> source 3 = plasma
sherpa> dataspace 2 (0.5:2.0:0.01) histogram
sherpa> dataspace 3 (0.7:2.8:0.01) histogram
sherpa> # get the source data (ie the model amplitudes)
sherpa> rest = get_source(2)
sherpa> obs  = get_source(3)
sherpa> rest
Double_Type[150]
sherpa> obs
Double_Type[210]
sherpa> # get the energy axis information
sherpa> rest_axis = get_photon_energy_axes(2)
sherpa> obs_axis  = get_photon_energy_axes(3)
sherpa> print(rest_axis)
axistype         =  Energy
axisunits        =  keV
lo               =  Double_Type[150]
hi               =  Double_Type[150]
mid              =  NULL
sherpa> print(obs_axis)
axistype         =  Energy
axisunits        =  keV
lo               =  Double_Type[210]
hi               =  Double_Type[210]
mid              =  NULL
sherpa> # set up the plot
sherpa> clear
sherpa> chips.curvestyle = _chips->step
sherpa> chips.symbolstyle = _chips->none
sherpa> limits 0.4 3 0.04 0.8
sherpa> log y
sherpa> # display the two energy ranges (with a slight offset in y)
sherpa> () = curve( rest_axis.lo, rest*0.8 )
sherpa> red
sherpa> () = curve( obs_axis.lo, obs*1.2 )
sherpa> green
sherpa> # Tidy up the plot
sherpa> xlabel "Energy (keV)"
sherpa> ylabel Flux
sherpa> title "XSMEKAL model: kT = 7 keV, Abundance = 0.3"
sherpa> xlabel size 1.4
sherpa> ylabel size 1.4
sherpa> title size 1.4
sherpa> redraw

This image  shows the resulting plot: the models were offset by a factor of 20 per cent in flux for display purposes.

1. Accuracy

The accuracy of the k-correction is determined by the width of the energy bins. We shall illustrate this by using the xspowerlaw model, for which an analytic solution for the k-correction exists: for a photon index of gamma, the k-correction is given by k(z) = (1+z)^gamma-2. Below we calculate the corrections for two bin widths - 0.1 and 0.01 keV - and compare them to the analytic solution.

sherpa> # set up a model (XSPEC models always have integrate set to on)
sherpa> erase all
sherpa> paramprompt off
Model parameter prompting is off
sherpa> source 1 = xspowerlaw[p1]
sherpa> p1.phoindx = 3
sherpa> # set up the redshift range and the analytic solutions
sherpa> z = [0:1.1:0.2]
sherpa> kz = (1.0+z)^(3-2)
sherpa> # calculate the k-correction using two different energy grids
sherpa> dataspace (0.1:10:0.1) histogram
sherpa> kz1 = calc_kcorr( 1, z, 0.5, 2.0 )
sherpa> dataspace (0.1:10:0.01) histogram
sherpa> kz2 = calc_kcorr( 1, z, 0.5, 2.0 )
sherpa> # write out the results
sherpa> writeascii( stdout, z, kz, kz1, kz2 )
0       1       1       1
0.2     1.2     1.20158 1.20002
0.4     1.4     1.40296 1.40003
0.6     1.6     1.60421 1.60004
0.8     1.8     1.80538 1.80005
1       2       2.00649 2.00007

The output columns list the redshift, analytic solution, then the values calculated by calc_kcorr() for the two bin widths (0.1 and then 0.01 keV). As can be seen, the smaller grid spacing produces results closer to the true value. The choice of grid spacing depends on the model(s) you use and the desired accuracy: for instance, in this example the error introduced by using a bin width of 0.1 keV is likely to be much smaller than the error on the measured flux of the source or uncertainties in the model parameters used to calculate the K correction.