How can we calculate a flux in Sherpa¶

This notebook is geared towards people who want to calculate the flux of a source for which they have an X-ray spectrum in PHA format. It is also aimed at people using Sherpa from CIAO - that is, they are used to the Sherpa threads.

The presentation is done using a Python (aka Jupyter) notebook, and is broken down into the following sections:

• How do I ...
• calculate a flux?
• generate erros?
• 1. Introduction
• 1.1 How to load Sherpa
• 1.2 How to create a model
• 1.3 Getting help
• 1.4 What model instances have been created?
• 1.5 The many ways of specifying a model expression
• 1.6 Skip ahead
• 1.7 How do I delete a model instance
• 1.8 Evaluating a model
• 2. Manually calculating the flux
• 3. Setting up the data in Sherpa
• 4. calculate_energy_flux
• 4.1 Absorbed versus unabsorbed fluxes
• 5. How about parameter errors?
• 5.1 sample_photon_flux and sample_energy_flux
• 5.2 correlated=True
• 5.3 A note on clipping
• 5.4 What errors are used?
• 5.5 sample_flux
• 5.6 sample_flux and unabsorbed fluxes
• 5.7 sample_energy_flux versus sample_flux
• 6. XSPEC's cflux model
• 7. Calculating luminosities
• 7.1 The K correction
• 7.2 luminosity distance
• 7.3 Putting it all together

History¶

CIAO 4.16¶

There have been no significant changes to the flux-calculation routines in this release, but there are several changes we do take advantage of:

• use the sherpa.astro.charge_e value for converting between keV and erg rather than a copy of this value (i.e. no more KEV_TO_ERG = 1.60217653e-09)
• the group_counts now defaults to only selecting the noticed data range, meaning we do not have to set the tabStops argument
• the new set_rng call is used to fix the random-number generator (so that repeated runs of this notebook will create the same output)
• CIAO 4.16 includes XSPEC version 12.13.1e and so now has access to the xscglumin convolution model

CIAO 4.14¶

The sample_flux command has seen a number of improvements so that it matches sample_energy_flux.

This release fixed a number of corner cases with filtering PHA data, and as part of that the get_filter output was changed to record the full energy range (previously it showed the mid-points for the start and end ranges: as the energy bins have a finite width, and are often grouped together the mid-point can make it harder to see the chosen range). This does not change the results of this notebook but may make small changes to your analysis.

The default xschatter setting of XSPEC models has been changed from 0 to 10, to match XSPEC, which may create more screen output when a model is evaluated.

CIAO 4.13¶

Version 4.13 contained bug fixes and improvements to the sample_energy_flux and sample_photon_flux functions (the handling of the samples parameter was improved, the model and clip parameters were added, and the return value gained an extra column). It also adds support for the XSPEC convolution models, which simplifies the XSPEC's cflux model section below compared to the CIAO 4.12.1 version.

The notebook integration of Sherpa was also improved, with "rich text" output used for various objects, as shown below.

CIAO 4.12.1¶

Version 4.12.1 - which was released July 14, 2020 - contained bug fixes and improvements to the calc_energy_flux and calc_photon_flux functions (the addition of the model parameter).

How do I ...¶

This section is for those people who have already read this document, or need an answer quickly. It is strongly recomended that you read the remainder of the document to explain the various trade offs in the various approaches.

The examples in this section assume the data and model use the default dataset (namely id=1).

... calculate a flux?¶

The calc_energy_flux function is used to calculate a model flux.

If you have a model expression like xsphabs.gal * powlaw1d.pl and want to calculate the unabsorbed flux of the power-law component, then use

calc_energy_flux(lo, hi, model=pl)

It is important to include the low and high limits for the calculation otherwise the response limits will be used.

If the model argument is not used then the absorbed flux will be returned.

The calc_photon_flux calculates fluxes in photon/cm$^2$/s rather than erg/cm$^2$/s.

... generate errors?¶

The sample_flux function will create n samples of the parameter distribution (estimating errors using the covar routine) and return the flux and error values for both the full model and (if the modelcomponent argument is set), the unabsorbed model. It also returns the flux distribution, but only for the full (absorbed) moodel.

The sample_energy_flux and sample_photon_flux functions have more functionality and should be used if you need something more complex than sample_flux.

First we start with some background material on using Sherpa:

1. Introduction¶

This document is a Jupyter notebook, a way of writing Python code along with documentation that be read or even re-run, allowing you to change variables or calls. This is not a guide on how to use Jupyter notebooks, for which there are many resources which start from the Jupyter home page. CIAO comes with a version of Jupyter, so all you need to do is say either of

jupyter notebook
jupyter lab

and either a web page will open in your browser, or you can use the address printed to the screen. Alternatives are to run Python directly, either with ipython or sherpa (which is just sets up ipython to display the Sherpa prompt and loads in Sherpa).

Note that the lab version requires you to have installed the jupyterlab package with pip or conda.

Other packages that are used by this notebook (you do not need them to calculate fluxes with Sherpa, but they may come in handy) can be installed with either command:

pip install scipy astropy
conda install scipy astropy

and (this must be installed with pip:

pip install corner

1.1 How to load Sherpa¶

In this notebook I am going to be importing the sherpa.astro.ui module without qualification. This is essentially the environment provided by the Sherpa "application" (which is actually just IPython with some minor configuration changes):

As this is a Jupyter notebook, I want any plots to appear within the page, and not as a separate window:

For reproducability, let's check what version of the packages we have installed (the versions you get may differ from these, as they depend on how CIAO was installed and what version of Python was selected):

Although not essential, the Python corner package provides helpful visualization capabilities which I will use in this notebook, and SciPy is used with the AstroPy cosmology routines to calculate a luminosity distance:

I am also going to ensure we have a fixed set of abundance and cross-section values used in this notebook (this is only needed to make the notebook repeatable, since these settings default to the values in your ~/.xspec/Xspec.init file):

These commands should create screen output, but due to the way XSPECs model library works the output may get directed to the terminal running the notebook application rather than the notebook itself. For this particular set of commands the output is:

 Solar Abundance Vector set to grsa:  Grevesse, N. & Sauval, A. J. Space Science Reviews 85, 161 (1998)
 Cross Section Table set to vern:  Verner, Ferland, Korista, and Yakovlev 1996

1.2 How to create a model¶

In Sherpa we have model classes (or types) and instances of those models, which we refer to with a user-supplied name. It is the model instance that we need here, and we can create in several different ways:

• explicitly, with create_model_component
• implicitly, using the "classname.instancename" syntax

Model types that come from the XSPEC model library - which in CIAO 4.16 is version XSPEC 12.13.1e - have names that begin with xs, and Sherpa directly supports additive, multiplicative, and convolution XSPEC models. That is, sources of emission (additive), simple changes to that emission (multiplicative), and models that can change the shape of the emission (convolution). These models are generally specific to X-ray astronomy, and are likely to be the ones you use to describe the emission from your source.

Sherpa also contains a number of simple one-dimensional geometric models which can also be used. These names normally end in 1d. There are also a number of other models which are intended for the optical pass band.

The names of the model types can be found with the list_models function, and the list_model_components function to list the model instances that have been created. As I haven't done anything, there are no model instances yet:

However, there are plenty of model types to chose from:

First I'm going to use create_model_component to create an instance called spl of the Sherpa power-law model type:

Certain Sherpa objects take advantage of the Python notebook dispay capabilities, as shown here with the model object (displayed as a HTML table). You can still display the data as you would at the Sherpa or IPython prompt (with print, as shown below), but this notebook will rely on the notebook output in most cases):

1.3 Getting help¶

You can find more information on the model by using the Python help function on the model instance: that is

help(spl)

Unfortunately you can't do this with the model name - that is help(powlaw1d) does not return anything useful (this is on our list of things to improve).

I am also going to create an instance of the XSPEC powerlaw model, imaginatively called xpl, using the "dotted" syntax approach:

We can see that the parameters have different names (e.g. gamma versus PhoIndex) and that Sherpa's powerlaw model lets you change the reference point (but that by default it is fixed to 1 which matches the XSPEC model).

1.4 What model instances have been created?¶

Since we now have some model instances, we now get a response from list_model_components:

We can retrieve models (if we have accidentally let a variable go out of scope or set it to something else) with the get_model_component function:

Note that as the spl variable still exists in this notebook, changing the gamma parameter of sillyName has also changed the spl variable (since they are actually referencing the same Python object):

This can be very useful, but can occasionally cause confusion, so I suggest trying to keep the Python variable names the same as the instance name you created. In other words, I should have said

spl = get_model_compoenent('spl')

earlier.

1.5 The many ways of specifying a model expression¶

It's rare that we have a situation where a single model type can describe the observed emission from a source, so Sherpa allows you to combine model instances with the normal mathematical operations (normally + and * but - and even / work too). You can even include numeric constants in the expression.

So, if I wanted to model a power law that has been absorbed by the galactic medium (using the XSPEC phabs model) I could say

set_source(xsphabs.gal * powlaw1d.pl)

or

mdl = xsphabs.gal * powlaw1d.pl)
set_source(mdl)

or

gal = create_model_component('xsphabs', 'gal')
pl = create_model_component('powlaw1d', 'pl')
set_source(gal * pl)

There is even a variant where you specify the model combination as a string: set_source("xsphabs.gal * powlaw1d.pl").

Once you have created a model instance, you do not need to specify the model type again, but it is not a problem if you do so. That is, the following two sets of commands are identical

set_source(1, xsphabs.gal * powlaw1d.pl1)
set_source(2, xsphabs.gal * powlaw1d.pl2)

and

set_source(1, xsphabs.gal * powlaw1d.pl1)
set_source(2, gal * powlaw1d.pl2)

1.6 Skip ahead¶

Normally we use models associated with an actual dataset, but in the next few subsections (1.6 to 2) we are going to explore how much you can do in the purely theoretical model space where you have not loaded any observational data. This is particularly useful when you just want to evaluate a model for a simple grid (say 0.1 to 10 keV with a bin width of 10 eV) for a plot, or calculation, and do not want to bother with a PHA file and response.

Skip ahead to section 3 to get back to our normal programming.

1.7 How do I delete a model instance¶

The delete_model_component function will remove the model instance from Sherpa. Fortunately Sherpa will stop you deleting any of the model components we used in set_source, so to illustrate this we have to make a new model component 'bngd' purely so we can delete it:

The delete_model_component call will also remove the variable from the default namespace, which means that the variable bgnd is no-longer defined, and trying to use it will return an error. I don't show that here as it breaks attempts to re-run the notebook!

1.8 Evaluating a model¶

Normally you don't need to evaluate a model directly (that is, calculate the predicted values on a given energy grid), but it can be useful, for instance

• to see how the model varies with parameter values
• for calculating things (such as the flux ;-)

When evaluating the models directly, as I do below, it is simplest to always use an grid of values in keV (there are ways to use wavelengths, with bin values in Angstroms, but it depends on what combination of models are being used). With that out of the way, I create a "typical" (at least for Chandra imaging-mode data) energy grid:

For this section I am going to use the XSPEC apec model, and see how varying the temperature and metallicity of a zero-redshift source changes the emission:

Evaluating the model is as simple as calling it with the grid you want to use. It uses the current parameter values, so the following evaluates the model for

Note that the first time the model was evaluated a message was displayed (Reading APEC data from ...). This is because the XSPEC chatter setting defaults to 10 to match XSPEC. It can be set to 0 to silence most messages if the output is distracting. See the help for get_xschatter and set_xschatter for more information.

These can be displayed (the return value is a NumPy array the same size as the input grid). For XSPEC additive models, such as xsapec, the return value has units of photon/cm$^2$/s.

Normally this is displayed "per keV", wihch we can get by dividing each bin by ehi - elo, but in this case each bin has the same width (0.01 keV), so all this would do would be to change the Y axis values, but not the shape.

No instruments!¶

The above has not used any instrument response (RMF or ARF), since we care here about the spectra before it has interacted with a telescope.

2. Manually calculating the flux¶

Since we can evaluate the model, we can use this to directly calculate the flux from a given model, as I'll discuss in this section. The following sections are all ways that Sherpa provides that automate this process, but I feel it's useful to know how to do it manually.

As mentioned when plotting the models, XSPEC additive models produce results with units of photon/cm$^2$/s. This means that summing up the values over over the desired energy range will give us the flux. As an example, I go back to the XSPEC power-law model (since this has an analytic form we can compare numerical results to):

If we represent the model as $S(E)$ with $$S(E) = A E^{-\gamma}$$ - where $A$ is the normalization and $\gamma$ equal to the PhoIndex parameter - then the photon flux for the energy range $E_1$ to $E_2$ is

$$ \begin{eqnarray} \rm{pflux}(E_1, E_2) & = & A \int_{E_1}^{E_2} E^{-\gamma} dE \\ & = & \frac{A}{1 - \gamma} [E^{1 - \gamma}]^{E_2}_{E_1} \end{eqnarray} $$

as long as $\gamma \ne 1$.

We can compare this analytic expression to the model-calculated version. First I set up a power-law model with a slope of 1.7 and normaliation of $10^{-4}$ photon/cm$^2$/s/keV:

So, the expected photon flux for the range 0.5-7 keV is:

For this example I am going to use the 0.5 to 7.0 keV passband. I could evaluate the model on elo and ehi (which are defined over 0.1 to 11 keV) and then filter the answer to the bins we want, but I've decided to just use the energy range I need. As I may just be re-using this a few times in the notebook, I have decided to make it into a small routine:

So, in this particular case the bin size doesn't matter, as even with only 10 bins the photon flux matches the expected value. This is because the power-law model is rather simple. It becomes more important when calculating the energy flux.

The energy flux is given by $\int E S(E) dE$ and so for the powerlaw model the analytic solution is

$$ \begin{eqnarray} \rm{eflux}(E_1, E_2) & = & A \int_{E_1}^{E_2} E . E^{-\gamma} dE \\ & = & \frac{A}{2 - \gamma} [E^{2 - \gamma}]^{E_2}_{E_1} \end{eqnarray} $$

as long as $\gamma \ne 2$. Note that this gives a result in keV/cm$^2$/s, so we are going to have to convert this to erg/cm$^2$/s using the same constant that Sherpa does.

The expected energy-flux for our power-law model is therefore:

This time the bin size does make a difference when calculating the model flux, since we are approximating the integral by the trapezoidal rule, and weighting each bin with its mid-point:

3. Setting up the data in Sherpa¶

Whilst you can calculate the flux manually, there are a number of routine provided by Sherpa for doing this for you. As they are intended for use as part of a "normal" Sherpa session (or script) they rely on having a PHA dataset loaded.

For this, I'll use a spectrum from release 2 of the Chandra Source Catalog, as I happen to have it on disk and it is well described by an absorbed power-law model.

For this example I want to use chi-square statistics, so I group the data into bins of at least 15 counts with the group_counts call (prior to CIAO 4.16 I had to use the tabStops argument to only apply the grouping over the selected energy range, but this is now the default behavior):

The data can be displayed in a notebook (for PHA data we get a plot of the data and tables summarizing the observation):

For this example I shall use Chi-square statistics (errors set to the square-root of the number of counts in a bin):

and for with an absorbed powerlaw model:

For this example I am going to fix the absorption to use the column density from the colden tool, which requires converting the approximate position of this source (00$^h$ 03$^m$ 39$^s$ +16$^\circ$ 02$^{'}$ 20$^{''}$) to decimal and then running colden, using Python modules provided with the CIAO contributed package:

The absorption models have units of

whereas colden returns values in units of 10$^{20}$ atoms cm$^{-2}$, so let's scale the result before freezing the parameter.

Now I can fit:

The reason for fixing the absorption to the Galactic value is that it is large enough to make a difference in the predicted spectrum (the blue line in the plot below, compared to the unabsorbed model shown by the orange line). If the $n_H$ parameter were allowed to vary the best-fit is much closer to 0.

4. calculate_energy_flux¶

The basic routines for PHA data are calc_photon_flux and calc_energy_flux. You can call them without any arguments, but this is not advised as the energy range is then set by the instrument response. I am going to always use them with a given pass band (i.e. energy range), so that I know what I'm calculating:

These results can be compared to the analytic answer (using the photon flux for demonstration purposes, and noting that the XSPEC and Sherpa power-law models are the same, they just have different parameter names):

4.1 Absorbed versus unabsorbed fluxes¶

The fact that the analytic solution is larger than the calculated value is because both calc_photon_flux and calc_energy_flux use the full model expression - that is an absorbed power-law - which is going to result in a smaller flux than the power-law only (i.e. unabsorbed) component.

The CIAO 4.12.1 release added the model argument to calc_energy_flux which makes it easy to calculate unabsorbed fluxes too. Setting the model argument to pl in this example will calculate the unabsorbed flux.

The photon-flux is now the same as the expected value (at least to the fourth decimal place in this particular example):

It is important to only calculate a flux for an XSPEC additive model component (or a combination of multiplicative * additive). The calc_energy_flux routine will allow you to specify a XSPEC multiplicative model - here this would be model=gal - but the number is meaningless.

5. How about parameter errors?¶

A simple solution for estimating the error on the flux is to just use the fractional error on the normalization (or amplitude) parameter, and apply that to the flux. However, this ignores the influence on any other model parameters on the flux calculation (so in our example the slope of the power law).

Sherpa provides three routines that will simulate the flux distribution based on the parameter errors:

• sample_photon_flux
• sample_energy_flux
• sample_flux

If you require photon fluxes then use sample_photon_flux, otherwise you can use either of sample_energy_flux or sample_flux. The sample_energy_flux routins has more features, but sample_flux may be a bit easier to use.

All three routines use the same technique: create a random set of the free parameters in the model using a Gaussian distribution centered on the best-fit parameter values and the widths from the errors on these parameters, and then use that set of parameters to calculate the flux. This step is repeated a given number of times, and the resulting flux distribution can be used to estimate the flux errors. Note that the Gaussian distribution can include correlations between the parameters, as we show below.

5.1 sample_photon_flux and sample_energy_flux¶

I am going to first look at sample_energy_flux. This will create num samples of the free parameters of the model (where num is an input argument and defaults to 1) assuming a normal distribution about the current parameter values. The default is to assume un-correlated errors, and to calculate these errors with the covar routine. We show how to change these assumptions in sections 5.2 and 5.4. Each sample is then used to calculate a flux, using calc_energy_flux, resulting in a distribution of fluxes.

Although not shown here, sample_photon_flux works the same way as sample_energy_flux, but returns fluxes in units of photon/cm$^2$/s rather than erg/cm$^2$/s.

We can get a visual handle of sample_energy_flux by using the plot_energy_flux function to plot up 1000 samples of the 0.5 to 7 keV photon flux (the number of bins used for the histogram has been reduced from 75 to 30 as it works better for this dataset, and 1000 iterations have been used as it doesn't take too long for this notebook rather than providing any particular level of accuracy). We can take advantage of the model parameter to calculate the fluxes for an unabsorbed component. As expected, the unabsorbed fluxes are on average larger than the absorbed fluxes!

Note that the process involves using random numbers, and so - to ensure the results do not change between runs - I set the random generator to a known value using the (new to CIAO 4.16) set_rng routine coupled with the NumPy random generator call:

The sample_energy_flux routine - which plot_energy_flux uses - returns the actual flux values per iteration (along with the parameter values), as a two-dimensional array (iteration number, and a combination of parameter values, clipping flag, and flux):

We can see that the fluxes - which are the first column - show the same behavior as seen in the plot_energy_flux plot (note that the results will not exactly match because the binning is different and the results shown here are a different realisation than in the plot_energy_flux call):

In the rest of this section I shall concentrate on the unabsorbed, rather than absorbed, fluxes.

Here I plot the ampl and $\gamma$ parameter values used (there's no assumed correlation between these two parameters since the correlated argument to sample_energy_flux was left at its default setting of False), along with a histogram of the flux distribution:

We can compare the power-law normalization (i.e. the ampl parameter) to the calculated flux, to show that they are correlated, but that there is scatter (thanks to $\gamma$ varying).

The last column (here referred to as clip2) is a flag column, and indicates whether any of the parameter samples fell outside of the limits set by the clip parameter and were replaced bu the limit. The default setting for clip is 'hard', which means that only those rows where a parameter falls outside its hard limits will be flagged. For this model and dataset it is unlikely to happen, which we can check by looking if any value in clip2 is not 0 (this column only has values of 0 and 1 but is stored as a floating-point value):

This is the time for the corner module to shine, since it is designed to visualize the correlations for data like the return value from sample_energy_flux. There are a number of customization options for the plot, but let's see what the defaults look like, other than adding in the column labels (thanks to the way that the notebook works, I add a semi-colon so that the return value isn't displayed automatically, as this then leads to the plot being duplicated).

Since the last column is all zeros I remove it from the call with sample2[:, :-1]:

The diagonal elements show the histograms of the columns, and the other terms plot the distribution of the two column values. The binning used in these histograms is designed to show the structure of the data (and is more advanced than the simple scatter plots I showed above).

We can even overlay the distributions we got for the unabsorbed (orange) on top of the absorbed (black) results. We can clearly see the shift in the flux distribution, but the parameter values remain similar:

We can use NumPy routines to extract the calculated energy flux values; in this case I chose to look at the median value along with the one-sigma range from the distribution.

We can use Sherpa's plot_cdf function for displaying the cumulative distribution function:

5.2 correlated=True¶

As shown earlier, there was no assumed correlation between the amplitude and slope parameters when simulating the fluxes. However, we normally find there is some correlation: for this fit we can use the reg_proj function to show that as the slope - i.e. $\gamma$ - increases, so does the amplitude.

If we set the correlated parameter to True when calling sample_energy_flux then the covariance matrix will be used to define parameter correlations.

First let's check whether we really do see a difference in the parameter values (i.e. were the correlations used in the samples?). I have overdrawn the contours from reg_proj to show that the correlated version matches, in shape at least.

So, the first run used uncorrelated samples, which results in an ellipse parallel to the horizontal, whereas the correlated version has the major axis rotated.

Does this make a difference to the flux distribution?

We can see below, where I overlay the cumulative distributions, that the medians are essentially the same$^\dagger$, but that the distribution of fluxes when the parameters are correlated is narrower (and so has a smaller error range).

$^\dagger$ although the exact difference depends on the random generator used.

This can also be seen by using the corner module to overplot the two distributions (the correlated version is in orange):

We can see that the flux distribution is narrower in the correlated version but centered about the same location, and the actual model parameters have similar one-dimensional distributions.

5.3 A note on clipping¶

If a sampled parameter exceeds the clip limit (by default the hard limit, but it can also be the soft limits) then that parameter will be replaced by the value at the limit, and the clip flag for the row (the last column) will be set to 1. This can be useful when a source is faint and so is not well detected, or a parameter is not well constrained (e.g. the n$_H$ parameter of an absorption model), and you want to post-process the fluxes to include or exclude parameters at the limit.

Note that the hard limits of a parameter are provided by the hard_min and hard_max attributes, as shown below for the power-law $\gamma$ and $n_H$ column density:

and the soft limits by the min and max attributes:

Section 5.7 contains an example when the clipping parameter makes a difference.

5.4 What errors are used?¶

The sample_energy_flux and sample_photon_flux routines require an error estimate to generate the parameter samples. The default is to run the covar routine to calculate the covariance matrix, and use that (just the diagonal elements when correlated is False).

The scales parameter can be used to send in either a one- or two-dimensional matrix (size and ordering match the free parameters of the model) to define the errors used in the simulation. A normal (i.e. symmetric) distribution is used in the sampling.

5.5 sample_flux¶

The sample_flux is similar to sample_energy_flux but has several differences in what values are returned. It also automatically clips the flux results to exclude any rows that are are equal to (or outside) the soft limits of the parameters, which can be problematic when calculating fluxes for un-detected sources (where parameter values may hit the soft limits0.

NOTE: Several restrictions and limitations in sample_flux have been lifted - in particular it can now be used with datasets other than the default identifier and the clipping has been changed. See "help(sample_flux)" for more details.

NOTE: In CIAO 4.15, the sample_flux routine is restricted to calculating the energy flux for PHA data sets (setting Xrays=False will result in an error).

Unlike sample_energy_flux, sample_flux returns three items:

• the errors for the "full" model (median and one-sigma limits);
• the errors for the "component" model (median and one-sigma limits);
• and the parameter values used for the simulations along with the fluxes for the "full" model.

It will also display the calculated errors on the screen. In the following I have not set the modelcomponent parameter - which is the same as the model argument for calc_energy_flux and sample_energy_flux - and so the first two arguments will be the same.

Section 5.6 discusses using the modelcomponent to calculate the fluxes for a single component.

The flux argument gives the median and one-sigma upper and lower limits (since I did not change the confidence parameter from its default value of 68):

The third element is similar to the return value of sample_energy_flux - that is the flux, clip flag, and parameter values for each iteration - except that it has an extra value, the statistic value for the sample:

We can use this extra information, for instance by coloring the points with the statistic value (top plot):

As with sample_energy_flux we can change the correlated flag or explicitly set the errors, with the scales parameter. Remember that the third samples argument returned by sample_flux contains - in its first column - the flux using the "full" model, even when the modelcomponent is given.

We again see that the flux distribution calculated using correlated errors is tighter than the un-correlated version (this time I use a histogram/PDF rather than cumulative plot):

As with sample_energy_flux, the corner module can also be used to display and compare the distributions, although it is a bit-more messy to remove the clipped column as it is no-longer the last column:

5.6 sample_flux and unabsorbed fluxes¶

Setting the modelcomponent argument tells sample_flux to calculate two sets of fluxes: the full source expression (original) value, reported first, and the component-only flux (component) second.

Here we can see that the unabsorbed flux (i.e. just pl and not gal * pl) is higher than the absorbed flux, as it should be.

Just to re-iterate, the median and one-sigma limits of the absorbed flux (first element) compared to the unabsorbed flux (second element):

That is

Note that in CIAO 4.16 the flux distribution returned in this case is still for the full model expression, and not the modelcomponent. We can show this by comparing the median of the flux column to the first elements of ures[0] and ures[1] and seeing which one matches:

5.7 sample_energy_flux versus sample_flux¶

Here I compare the flux distribution for the fluxes from sample_energy_flux and sample_flux, using the correlated=True case and for the unabsorbed power-law component (re-running everything since I have decided to chose a softer energy band for this calculation).

Note that the clip parameter for sample_energy_flux is set to soft to more-closely match the behavior of sample_flux (this behavior changed in CIAO 4.14).

In order to discuss the clipping behavior, I have thawed the n$_H$ parameter since the error range is large enough to cover the minimum value for the parameter (which in this case is n$_H$ = 0).

The sample_energy_flux and sample_flux calls are as before, but now show off a warning about hitting the minimum value for the gal.nH parameter:

We can see how the boundary has affected the analysis by comparing the n$_H$ and $\gamma$ parameter values: the cut off at n$_H$=0 is clear, and the two distributions appear similar.

We can view the sample_energy_flux parameters that were clipped (that is, in this case the n$_H$ parameter was reset to 0 from a negative value) by plotting the trace of all the fluxes (as a line) and then overplotting those iterations where the parameters were clipped as circles:

We can see what difference these boundary points can make by comparing the full distribution to only those points that were not clipped:

In this case removing the clipped points leads to larger fluxes, which can bias the calculation of upper-limits.

Prior to CIAO 4.14 the sample_flux command removed points at the soft limit, which meant that it could bias the results. It now matches sample_energy_flux as we can see below, where the limits returned by sample_flux are close to the "only within the soft limits" version of the sample_energy_flux command.

We have to compare the un-clipped data from sample_energy_flux to the sample_flux results to get a match.

We can identify the number of rows that are affected from sample_energy_flux and sample_flux directly (as shown above):

We can use this to look at the statistic values and see where the clipped values lie:

6. XSPEC's cflux model¶

XSPEC provides several "convolution" models for calculating flux and luminosites.

You do not need to use xscflux, as the previous approaches (hopefully) suffice, but I mention it if you want to compare your results to a colleague's XSPEC analysis.

To show how to use the cflux component I need a copy of the source dataset. Using get_data and set_data makes this fast (although I have to be careful not to change settings, such as the noticed energy range, as it will also change the values for dataset 1).

I first want to set the n$_H$ value back to a fixed value, to make it easy to compare to earlier results (and re-fit):

To start with I need to create an instance of the cflux model. Unlike normal models, the experimental interface uses the load_xxx pattern to create a component with the given name:

The parameters are described in the xscflux ahelp page, but briefly:

• Emin and Emax give the energy range for the flux
• lg10Flux is the log (base 10) of the flux of the model components it "contains"

In this case I am going to apply it just to the power-law model, so I can get the unabsorbed flux:

We can verify the model, and see that the convolution component is only applied to the emission model (i.e. power-law):

With this setup, I can fit as normal, but instead of the amplitude the model tells me the flux:

The fit looks believable (as it should since the reduced statistic is $\sim 0.9$):

We can see that the power-law slope is essentially identical to that I got earlier, as is the best-fit statistic (that is, the use of xscflux hasn't significantly changed the fit results):

The flux (in this case it is an energy flux) can be calculated directly from the model:

We can get an error using covar (or conf but that takes longer to run):

I can now compare the sample_flux results from above, namely:

to the cflux results:

There is also the xscpflux model for if you want to calculate the photon flux, and the xsclumin and xscglumin (new in CIAO 4.16) models for calculating the luminosity, although I suggest doing the flux-to-luminosity calculation yourself taking advantage of the AstroPy cosmology module, as discussed below.

7. Calculating luminosities¶

The following section is only relevant for "properly" extragalactic objects, that is those objects for which you need to worry about luminosity distances and k corrections.

7.1 The K correction¶

The K correction is used to account for the fact that the observed and rest-frame pass bands are not the same. For simple models (such as a power-law) it can be calculated analytically, but in general we will not be so lucky. Although it can be calculated manually (using the same techniques as you would use to calculate the flux manually), Sherpa provides the calc_kcorr function. The following section just shows how to call this routine, as there's much-more information provided in the Calculating K-corrections using Sherpa thread, which should be reviewed by the interested reader (particularly for models which contain a redshift parameter, as this needs to be handled carefully when used with calc_kcorr).

For this example I am going to assume the source is at a redshift of 0.4, and that I want the flux in the rest-frame 0.5 to 7 keV:

7.2 luminosity distance¶

So, how do we go about using the AstroPy cosmology module then? For this example I'm going to use the Planck15 cosmology - but there are a number of pre-defined cosmologies I could have picked from (or created one myself):

This step does require that you have installed AstroPy and SciPy into your CIAO installation. For conda users this can be done with

conda install astropy scipy

and for people whose CIAO was installed with ciao-install

pip install astropy scipy

should also work.

NOTE it is possible that installing packages into CIAO could install an incompatible version of a library, so please be careful (and we recommend that you keep a copy of the install files if you used ciao-install so that you can re-create your CIAO installation in case of problems).

With a cosmology, we can now calculate the luminosity distance to our source (using my ficticious redshift of 0.4).

7.3 Putting it all together¶

Taking advantage of the support for units in AstroPy, let's take the median unabsorbed flux calculated by sample_flux and give it its units, then use it to calculate the rest-frame 0.5 to 7 keV luminosity for the source:

There is special code to display the units "nicely" in a Python notebook:

or we can print it out normally:

Unfortunately this is not the friendly version (I tend to use erg/s) so can we simplify this? We can use the decompose method, which returns a value in units of Watts (although split into its fundamental components):

Instead, we can use to and tell it the units we want: