Deproject¶

Deproject is a CIAO Sherpa extension package to facilitate deprojection of two-dimensional circular annular X-ray spectra to recover the three-dimensional source properties. For typical thermal models this would include the radial temperature and density profiles. This basic method has been used extensively for X-ray cluster analysis and is the basis for the XSPEC model projct. The deproject module brings this functionality to Sherpa as a Python module that is straightforward to use and understand.

The module can also be used with the standalone Sherpa release, but as it requires support for XSPEC models - for the thermal plasma emission codes - the documentation will focus on using Sherpa with a CIAO environment.

Overview¶

The deproject module uses specstack to allow for manipulation of a stack of related input datasets and their models. Most of the functions resemble ordinary Sherpa commands (e.g. set_par, set_source, ignore) but operate on a stack of spectra.

The basic physical assumption of deproject is that the extended source emissivity is constant and optically thin within spherical shells whose radii correspond to the annuli used to extract the specta. Given this assumption one constructs a model for each annular spectrum that is a linear volume-weighted combination of shell models. The geometry is illustrated in the figure below (which would be rotated about the line to the observer in three-dimensions):

Model creation¶

It is assumed that prior to starting deproject the user has extracted source and background spectra for each annulus. By convention the annulus numbering starts from the inner radius at 0 and corresponds to the dataset id used within Sherpa. It is not required that the annuli include the center but they must be contiguous between the inner and outer radii.

Given a spectral model M[s] for each shell s, the source model for dataset a (i.e. annulus a) is given by the sum over s >= a of vol_norm[s,a] * M[s] (normalized volume * shell model). The image above shows shell 3 in blue and annulus 2 in red. The intersection of (purple) has a physical volume defined as vol_norm[3,2] * v_sphere where v_sphere is the volume of the sphere enclosing the outer shell (as discussed below).

The bookkeeping required to create all the source models is handled by the deproject module.

Fitting¶

Once the composite source models for each dataset are created the fit analysis can begin. Since the parameter space is typically large the usual procedure is to initally fit using the “onion-peeling” method:

First fit the outside shell model using the outer annulus spectrum

Freeze the model parameters for the outside shell

Fit the next inward annulus / shell and freeze those parameters

Repeat until all datasets have been fit and all shell parameters determined.

From this point the user may choose to do a simultanenous fit of the shell models, possibly freezing some parameters as needed. This process is made manageable with the specstack methods that apply normal Sherpa commands like freeze or set_par to a stack of spectral datasets.

Densities¶

Physical densities (\({\rm cm}^{-3}\)) for each shell can be calculated with deproject assuming the source model is based on a thermal model with the “standard” normalization (from the XSPEC documentation):

\[{\rm norm} = \frac{10^{-14}}{4\pi [D_A (1+z)]^2} \int n_e n_H dV\]

where \(D_A\) is the angular size distance to the source (in cm), \(z\) is the source redshift, and \(n_e\) and \(n_H\) are the electron and Hydrogen densities (in \({\rm cm}^{-3}\)).

Inverting this equation and assuming constant values for \(n_e\) and \(n_H\) gives

\[n_e = \sqrt{\frac{4\pi \, {\rm norm} \, [D_A (1+z)]^2 \, 10^{14}}{\alpha V}}\]

where \(V\) is the volume, \(n_H = \alpha n_e\), and \(\alpha\) is taken to be 1.18 (and this ratio is constant within the source).

Recall that the model components for each volume element (intersection of the annular cylinder a with the spherical shell s) are multiplied by a volume normalization:

\[V_{\rm norm}[s, a] = V[s, a] / V_{\rm sphere}\]

where \(V_{\rm sphere}\) is the volumne of the sphere enclosing the outermost annulus.

With this convention the volume (\(V\)) used above in calculating the electron density for each shell is always \(V_{\rm sphere}\).

Installation¶

As of version 0.2.0, the deproject package can be installed directly from PyPI. This requires CIAO 4.11 or later (as installation with pip in earlier versions of CIAO is not well supported). The module can be used with the standalone release of Sherpa, but it is only useful if Sherpa has been built with XSPEC support which is trickier to achieve than we would like.

Requirements¶

The package uses Astropy and SciPy, for units support and cosmological-distance calculations. It is assumed that Matplotlib is available for plotting (which is included in CIAO 4.11).

Using pip¶

CIAO¶

Installation within CIAO requires:

having sourced the CIAO initialisation script (e.g. ciao.csh or ciao.bash);
and using a constraints file, to avoid updating NumPy in CIAO.

It should be as simple as:

echo "numpy==1.12.1" > constraints.txt
pip install -c constraints.txt 'astropy<3.1' deproject

Note

The constraints are for CIAO 4.11, please adjust accordingly if using a newer version of CIAO (the Astropy restriction is because version 3.1 requires NumPy version 1.13 or later but CIAO 4.11 is shipped with NumPy version 1.12).

Standalone¶

When using the standalone Sherpa intallation, the following should be all that is required:

pip install deproject

Manual installation¶

The source is available on github at https://github.com/sherpa-deproject/deproject, with releases available at https://github.com/sherpa-deproject/deproject/releases.

After downloading the source code (whether from a release or by cloning the repository) and moving into the directory (deproject-<version> or deproject), installation just requires:

python setup.py install

Note

This command should only be run after setting up CIAO or whatever Python environment contains your Standalone Sherpa installation.

Test¶

The source installation includes a basic test suite, which can be run with either

pytest

or

python setup.py test

Example data¶

The example data can be download from either http://cxc.cfa.harvard.edu/contrib/deproject/downloads/m87.tar.gz or from GitHub. The source distribution includes scripts - in the examples directory - that can be used to replicate both the basic M87 example and the follow-on example combining annuli.

As an example (from within an IPython session, such as the Sherpa shell in CIAO):

>>> %run fit_m87.py
...
... a lot of screen output will whizz by
...

The current density estimates can then be displayed with:

>>> dep.density_plot()

the reduced-statistic for the fit to each shell with:

>>> dep.fit_plot('rstat')

and the fit results for the first annulus can be displayed using the Sherpa functions:

>>> set_xlog()
>>> plot_fit_delchi(0)

Changes¶

Version 0.2.0¶

Overview¶

The code has been updated to run with Python 3 and can now be installed from PyPI. Documentation has been moved to Read The Docs.

The deproject module now requires Astropy, which can be installed within CIAO 4.11. The three main areas where Astropy functionality is used are:

the use of Astropy Quantity values (both for arguments to methods and returned values);
Astropy Data Tables are used to return tabular data;
and Cosmology calculations now use the Astropy cosmology module rather than the cosmocalc module.

The deproject_from_xflt() helper function has been introduced, which uses the XFLT0001 to XFLT0005 keywords in the input files to determine the annulus parameters (radii and covering angle). The covering angle (\(\theta\)) can now vary per annulus.

Error values can now be generated using the onion-peeling approach, for the confidence and covariance methods, and the values are returned as an Astropy Table. Parameter values can now be tied together (to combine annuli to try and avoid “ringing”). There is improved support for accessing and plotting values.

Details¶

The code has been re-arranged into the deproject package, which means that you really should say from deproject.deproject import Deproject, but the deproject module re-exports deproject.deproject so that existing scripts still work, and you do not not have to type in the same word multiple times! The package has been updated so that it is available on PyPI.

The scaling between shells (calculated from the intersection between spheres and cylinders) was limited to 5 decimal places, which could cause problems with certain choices of annuli (such as an annulus making no contribution to interior annuli). This restriction has been removed.

Added support for per-annulus theta values (that is, each annulus can have a different opening angle). The radii, theta, and angdist parameters to Deproject all now require values that is an Astropy quantity rather than a dimensionless value.

Added the deproject_from_xflt() helper function, which creates a Deproject instance from PHA files which contain the XSPEC XFLT0001 to XFLT0005 keywords (as used by the projct model), rather than specifying the values from the command line. The routine will error out if the keywords indicate elliptical annuli, and the default is to assume the radii are in arcseconds, but a scaling factor can be given if the radii are in some other units (such as pixels).

Added the guess() method to do an initial fit to each annulus, following the approach suggested in the XSPEC documentation for projct, by just fitting the individual (not de-projected) models to each annulus. This can help speed up the deproject fit - fit() - as well as help avoid the fit getting stuck in a local minimum.

Added covar() and conf() methods that estimate errors - using the covariance and confidence methods respectively - using the onion-skin model (i.e. the errors on the outer annuli are evaluated, then this component is frozen and the errors on the next annulus are evaluated).

The fit(), conf(), and covar() methods now all return Astropy Tables containing the results per annulus. These values can also be retrieved with the get_fit_results(), get_conf_results(), or get_covar_results() methods. A number of columns (radii and density) are returned as Astropy quantities.

The cosmocalc module has been removed and the Astropy cosmology module is used instead. This is only used if the angular-diameter distance to the source is calculated rather than explicitly given. The default cosmology is now set to Planck15.

Values, as a function of radius, can be plotted with a number of new methods: fit_plot(), conf_plot(), and covar_plot() display the last fit results (with the last two including error estimates), and the par_plot() and density_plot() methods show the current values. These support a number of options, including switching between angular and physical distances for the radii.

The get_shells() method has been added to make it easy to see which annuli are combined together, and the get_radii() method to find the radii of the annuli (in a range of units).

Added the tie_par() and untie_par() methods to make it easy to tie (or untie) parameters in neighbouring annuli. The onion-skin approach - used when fitting or running an error analysis - recognizes annuli that are tied together and fits these simultaneously, rather than individually.

The set_source() method can now be called multiple times (previously it would lead to an error).

Added error checking for several routines, such as thaw() when given an unknown parameter name.

Updated to support Python 3.5 and to have better support when the pylab backend is selected. Support for the ChIPS backend is limited. A basic test suite has been added.

Version 0.1.0¶

Initial version.

To Do¶

Use the Python logging module to produce output and allow for a verbosity setting. [Easy]

Move from using obsid to a more-generic label for the multi-data-set case (i.e. when there are multiple data sets in a given annulus)

Allow the theta value to vary between data set in each annulus (for the multi-data-set case)

Support elliptical annuli

Create and use more generalized ModelStack and DataStack classes to allow for general mixing models. [Hard]

Add summary methods (similar to Sherpa’s show series of commands)

Add a specialized version of plot_source_component which takes a general model expression (i.e. without the _annulus suffix) and displays all combos per annulus. A bit tricky to get right. There’s also a plot_model_component.

M87¶

Now we step through in detail the fit_m87.py script in the examples directory to explain each step and illustrate how to use the deproject module. This script should serve as the template for doing your own analysis.

This example uses extracted spectra, response products, and analysis results for the Chandra observation of M87 (obsid 2707). These were kindly provided by Paul Nulsen. Results based on this observation can be found in Forman et al 2005 and via the CXC Archive Obsid 2707 Publications list.

The examples assume this is being run directly from the Sherpa shell, which has already imported the Sherpa module. If you are using IPython or a Jupyter notebook then the following command is needed:

>>> from sherpa.astro.ui import *

Set up¶

The first step is to load in the Deproject class:

>>> from deproject import Deproject

and then set a couple of constants (note that both the angular-diameter distance and thta values must be Astropy quantites):

>>> from astropy import units as u
>>> redshift = 0.004233                     # M87 redshift
>>> angdist = 16 * u.Mpc                    # M87 distance
>>> theta = 75 * u.deg                      # Covering angle of sectors

Next we create an array of the the annular radii in arcsec. The numpy.arange method here returns an array from 30 to 640 in steps of 30. These values were in pixels in the original spectral extraction so we multiply by the pixel size (0.492 arcseconds) to convert to an angle.:

>>> radii = numpy.arange(30., 640., 30) * 0.492 * u.arcesc

The radii parameter must be a list of values that starts with the inner radius of the inner annulus and includes each radius up through the outer radius of the outer annulus. Thus the radii list will be one element longer than the number of annuli.

>>> print(radii)
[ 14.76  29.52  44.28  59.04  73.8   88.56 103.32 118.08 132.84 147.6
 162.36 177.12 191.88 206.64 221.4  236.16 250.92 265.68 280.44 295.2
 309.96] arcsec

Now the key step of creating the Deproject object dep. This object is the interface to the all the deproject methods used for the deprojection analysis.

>>> dep = Deproject(radii, theta=theta, angdist=angdist)

If you are not familiar with object oriented programming, the dep object is just a thingy that stores all the information about the deprojection analysis (e.g. the source redshift, PHA file information and the source model definitions) as object attributes. It also has object methods (i.e. functions) you can call such as dep.get_par(parname) or dep.load_pha(file). The full list of attributes and methods are in the deproject module documentation.

In this particular analysis the spectra were extracted from a 75 degree sector of the annuli, hence theta is set in the object initialization, where the units are set explicitly using the Astropy support for units. Note that this parameter only needs to be set if any of the annuli are sectors rather than the full circle. Since the redshift is not a good distance estimator for M87 we also explicitly set the angular size distance using the angdist parameter.

Load the data¶

Now load the PHA spectral files for each annulus using the Python range function to loop over a sequence ranging from 0 to the last annulus. The load_pha() call is the first example of a deproject method (i.e. function) that mimics a Sherpa function with the same name. In this case dep.load_pha(file, annulus) loads the PHA file using the Sherpa load_pha function but also registers the dataset in the spectral stack:

>>> for annulus in range(len(radii) - 1):
...    dep.load_pha('m87/r%dgrspec.pha' % (annulus + 1), annulus)

Note

The annulus parameter is required in dep.load_pha() to allow multiple data sets per annulus, such as repeated Chandra observations or different XMM instruments.

Create the model¶

With the data loaded we set the source model for each of the spherical shells with the set_source() method. This is one of the more complex bits of deproject. It automatically generates all the model components for each shell and then assigns volume-weighted linear combinations of those components as the source model for each of the annulus spectral datasets:

>>> dep.set_source('xswabs * xsmekal')

The model expression can be any valid Sherpa model expression with the following caveats:

Only the generic model type should be specified in the expression. In typical Sherpa usage one generates the model component name in the model expression, e.g. set_source('xswabs.abs1 * xsmekal.mek1'). This would create model components named abs1 and mek1. In dep.set_source() the model component names are auto-generated as <model_type>_<shell>.

Only one of each model type can be used in the model expression. A source model expression like "xsmekal + gauss1d + gauss1d" would result in an error due to the model component auto-naming.

Next any required parameter values are set and their freeze or thaw status are set.

>>> dep.set_par('xswabs.nh', 0.0255)
>>> dep.freeze('xswabs.nh')

>>> dep.set_par('xsmekal.abundanc', 0.5)
>>> dep.thaw('xsmekal.abundanc')

>>> dep.set_par('xsmekal.redshift', redshift)

As a convenience if any of the model components have a redshift parameter that value will be used as the default redshift for calculating the angular size distance.

Define the data to fit¶

Now the energy range used in the fitting is restricted using the stack version of the Sherpa ignore command. The notice command is also available.

>>> dep.ignore(None, 0.5)
>>> dep.ignore(1.8, 2.2)
>>> dep.ignore(7, None)

Define the optimiser and statistic¶

At this point the model is completely set up and we are ready to do the initial “onion-peeling” fit. As for normal high-signal fitting with binned spectra we issue the commands to set the optimization method and the fit statistic.

In this case we use the Levenberg-Marquardt optimization method with the \(\chi^2\) fit statistic, where the variance is estimated using a similar approach to XSPEC.

>>> set_method("levmar")
>>> set_stat("chi2xspecvar")
>>> dep.subtract()

Note

The chi2xspecvar statistic is used since the values will be compared to results from XSPEC.

Seeding the fit¶

We first try to “guess” a good starting point for the fit (since the default fit parameters, in particular the normalization, are often not close to the expected value). In this case the guess() method implements a scheme suggested in the projct documentation, which fits each annulus with an un-projected model and then corrects the normalizations for the shell/annulus overlaps:

>>> dep.guess()
Projected fit to annulus 19 dataset: 19
Dataset               = 19
...
Change in statistic   = 4.60321e+10
   xsmekal_19.kT   2.748        +/- 0.0551545
   xsmekal_19.Abundanc   0.429754     +/- 0.0350466
   xsmekal_19.norm   0.00147471   +/- 2.45887e-05
Projected fit to annulus 18 dataset: 18
Dataset               = 18
...
Change in statistic   = 4.61945e+10
   xsmekal_18.kT   2.68875      +/- 0.0543974
   xsmekal_18.Abundanc   0.424198     +/- 0.0335045
   xsmekal_18.norm   0.00147443   +/- 2.43454e-05
...
Projected fit to annulus 0 dataset: 0
Dataset               = 0
...
Change in statistic   = 2.15043e+10
   xsmekal_0.kT   1.57041      +/- 0.0121572
   xsmekal_0.Abundanc   1.05816      +/- 0.0373779
   xsmekal_0.norm   0.00152496   +/- 2.93297e-05

As shown in the screen output, the guess routine fits each annulus in turn, from outer to inner. After each fit, the normalization (the xsmekal_*.norm terms) are adjusted by the filling factor of the shell.

Note

The guess step is optional. Parameter values can also be set, either for an individual annulus with the Sherpa set_par function, such as:

set_par('xsmekal_0.kt', 1.5)

or to all annuli with the Deproject method set_par():

dep.set_par('xsmekal.abundanc', 0.3)

The data can be inspected using normal Sherpa commands. The following shows the results of the guess for dataset 0, which corresponds to the inner-most annulus (the datasets attribute can be used to map between annuus and dataset identifier).

>>> set_xlog()
>>> plot_fit(0)

The overall shape of the model looks okay, but the normalization is not (since it has been adjusted by the volume-filling factor of the intersection of the annulus and shell). The gap in the data around 2 keV is because we explicitly excluded this range from the fit.

Note

The Deproject class also provides methods for plotting each annulus in a separate window, such as plot_data() and plot_fit().

Fitting the data¶

The fit() method performs the full onion-skin deprojection (the output is similar to the guess method, with the addition of parameters being frozen as each annulus is processed):

>>> onion = dep.fit()
Fitting annulus 19  dataset: 19
Dataset               = 19
...
Reduced statistic     = 1.07869
Change in statistic   = 1.29468e-08
   xsmekal_19.kT   2.748        +/- 0.0551534
   xsmekal_19.Abundanc   0.429758     +/- 0.0350507
   xsmekal_19.norm   0.249707     +/- 0.00416351
Freezing xswabs_19
Freezing xsmekal_19
Fitting annulus 18  dataset: 18
Dataset               = 18
...
Reduced statistic     = 1.00366
Change in statistic   = 13444.5
   xsmekal_18.kT   2.41033      +/- 0.274986
   xsmekal_18.Abundanc   0.375824     +/- 0.136926
   xsmekal_18.norm   0.0551815    +/- 0.00441495
Freezing xswabs_18
Freezing xsmekal_18
...
Freezing xswabs_1
Freezing xsmekal_1
Fitting annulus 0  dataset: 0
Dataset               = 0
...
Reduced statistic     = 2.72222
Change in statistic   = 12115.6
   xsmekal_0.kT   1.63329      +/- 0.0278325
   xsmekal_0.Abundanc   1.1217       +/- 0.094871
   xsmekal_0.norm   5.7067       +/- 0.262766
Change in statistic   = 13699.6
   xsmekal_0.kT   1.64884      +/- 0.0242336
   xsmekal_0.Abundanc   1.14629      +/- 0.0903398
   xsmekal_0.norm   5.68824      +/- 0.245884
...

The fit to the central annulus (dataset 0) now looks like:

>>> plot_fit_delchi(0)

After the fit process each shell model has an association normalization that can be used to calculate the densities. This is where the source angular diameter distance is used. If the angular diameter distance is not set explicitly in the original dep = Deproject(...) command then it is calculated automatically from the source redshift and an assumed Cosmology, which if not set is taken to be the Planck15 model provided by the astropy.cosmology module.

One can examine the values being used as follows (note that the angdist attribute is an Astropy Quantity):

>>> print("z={:.5f} angdist={}".format(dep.redshift, dep.angdist))
z=0.00423 angdist=16.0 Mpc

New in version 0.2.0 is the behavior of the fit method, which now returns an Astropy table which tabulates the fit results as a function of annulus. This includes the electron density, which can also be retrieved with get_density(). The fit results can also be retrieved with the get_fit_results() method.

>>> print(onion)
annulus rlo_ang rhi_ang ...     xsmekal.norm          density
         arcsec  arcsec ...                           1 / cm3
------- ------- ------- ... ------------------- --------------------
      0   14.76   29.52 ...  5.6882389946381275   0.1100953546292787
      1   29.52   44.28 ...  2.8089409208987233  0.07736622021374819
      2   44.28   59.04 ...   0.814017947154132  0.04164827967805805
      3   59.04    73.8 ...  0.6184339453006981  0.03630168106524076
      4    73.8   88.56 ...  0.2985327157676323 0.025221797991301052
      5   88.56  103.32 ... 0.22395312678845017 0.021845331641349316
    ...     ...     ... ...                 ...                  ...
     13  206.64   221.4 ... 0.07212121560592619 0.012396857131392835
     14   221.4  236.16 ... 0.08384338967492334  0.01336640115325031
     15  236.16  250.92 ... 0.07104455447410102 0.012303975980575187
     16  250.92  265.68 ... 0.08720295254137615 0.013631563529090736
     17  265.68  280.44 ... 0.09192970392746878 0.013996131292837352
     18  280.44   295.2 ... 0.05518150227052414 0.010843683594144967
     19   295.2  309.96 ... 0.24970680803387219 0.023067220584935984
Length = 20 rows

Inspecting the results¶

The electron density can be retrieved directly from the fit results, with the get_density() method:

>>> print(dep.get_density())
[ 0.11009535  0.07736622  0.04164828  0.03630168  0.0252218   0.02184533
  0.01525456  0.01224972  0.01942528  0.01936785  0.01905983  0.01568478
  0.01405426  0.01239686  0.0133664   0.01230398  0.01363156  0.01399613
  0.01084368  0.02306722] 1 / cm3

or plotted using density_plot():

>>> dep.density_plot()

The temperature profile from the deprojection can be plotted using par_plot():

>>> dep.par_plot('xsmekal.kt')

The unphysical temperature oscillations seen here highlights a known issue with this analysis method (e.g. Russell, Sanders, & Fabian 2008).

It can also be instructive to look at various results from the fit, such as the reduced statistic for each annulus (as will be shown below, there are fit_plot(), conf_plot(), and covar_plot() variants):

>>> dep.fit_plot('rstat')

Error analysis¶

Errors can also be calculated, on both the model parameters and the derived densities, with either the conf() or covar() methods. These apply the confidence and covariance error-estimation routines from Sherpa using the same onion-skin model used for the fit, and are new in version 0.2.0. In this example we use the covariance version, since it is generally faster, but confidence is the recommended routine as it is more robust (and calculates asymmetric errors):

>>> errs = dep.covar()

As with the fit method, both conf and covar return the results as an Astropy table. These can also be retrieved with the get_conf_results() or get_covar_results() methods. The columns depend on the command (i.e. fit or error results):

>>> print(errs)
annulus rlo_ang rhi_ang ...        density_lo             density_hi
         arcsec  arcsec ...         1 / cm3                1 / cm3
------- ------- ------- ... ----------------------- ----------------------
      0   14.76   29.52 ...  -0.0023299300992805777  0.0022816336635322065
      1   29.52   44.28 ...  -0.0015878693875514133  0.0015559288091097495
      2   44.28   59.04 ...  -0.0014852395638492444  0.0014340671599212054
      3   59.04    73.8 ...  -0.0011539611188859725    0.00111839214283211
      4    73.8   88.56 ...   -0.001166653616914693   0.001115024462355424
      5   88.56  103.32 ...  -0.0010421595234548324  0.0009946565126391325
    ...     ...     ... ...                     ...                    ...
     13  206.64   221.4 ...  -0.0005679816551559976  0.0005430748019680798
     14   221.4  236.16 ... -0.00048083556526315463 0.00046412880973027357
     15  236.16  250.92 ...  -0.0004951659664768401 0.00047599490797040934
     16  250.92  265.68 ...  -0.0004031089363366082  0.0003915259159180066
     17  265.68  280.44 ...  -0.0003647519140328147 0.00035548459401609986
     18  280.44   295.2 ...                     nan                    nan
     19   295.2  309.96 ... -0.00019648511719753958 0.00019482554589523096
Length = 20 rows

Note

With this set of data, the covariance method failed to calculate errors for the parameters of the last-but one shell, which is indicated by the presence of NaN values in the error columns.

The fit or error results can be plotted as a function of radius with the fit_plot(), conf_plot(), and covar_plot() methods (the latter two include error bars). For example, the following plot mixes these plots with Matplotlib commands to compare the temperature and abundance profiles:

>>> plt.subplot(2, 1, 1)
>>> dep.covar_plot('xsmekal.kt', clearwindow=False)
>>> plt.xlabel('')
>>> plt.subplot(2, 1, 2)
>>> dep.covar_plot('xsmekal.abundanc', clearwindow=False)
>>> plt.subplots_adjust(hspace=0.3)

The derived density profile, along with errors, can also be displayed (the X axis is displayed using angular units rather than as a length):

>>> dep.covar_plot('density', units='arcmin')

Comparing results¶

In the images below the deproject results (blue) are compared with values (gray boxes) from an independent onion-peeling analysis by P. Nulsen using a custom perl script to generate XSPEC model definition and fit commands. These plots were created with the plot_m87.py script in the examples directory. The agreement is good (note that the version of XSPEC used for the two cases does not match):

Combining shells¶

We now look at tie-ing parameters from different annuli (that is, forcing the temperature or abundance of one shell to be the same as a neighbouring shell), which is a technique used to try and reduce the “ringing” that can be seen in deprojected data (e.g. Russell, Sanders, & Fabian 2008 and seen in the temperature profile of the M87 data). Sherpa supports complicated links between model parameters with the link and unlink functions, and the Deproject class provides helper methods - tie_par() and untie_par() - that allows you to easily tie together annuli.

Starting from where the M87 example left off, we have:

The get_shells() method lists the current set of annuli, and show that there are currently no grouped (or tied-together) annuli:

>>> dep.get_shells()
[{'annuli': [19], 'dataids': [19]},
 {'annuli': [18], 'dataids': [18]},
 {'annuli': [17], 'dataids': [17]},
 {'annuli': [16], 'dataids': [16]},
 {'annuli': [15], 'dataids': [15]},
 {'annuli': [14], 'dataids': [14]},
 {'annuli': [13], 'dataids': [13]},
 {'annuli': [12], 'dataids': [12]},
 {'annuli': [11], 'dataids': [11]},
 {'annuli': [10], 'dataids': [10]},
 {'annuli': [9], 'dataids': [9]},
 {'annuli': [8], 'dataids': [8]},
 {'annuli': [7], 'dataids': [7]},
 {'annuli': [6], 'dataids': [6]},
 {'annuli': [5], 'dataids': [5]},
 {'annuli': [4], 'dataids': [4]},
 {'annuli': [3], 'dataids': [3]},
 {'annuli': [2], 'dataids': [2]},
 {'annuli': [1], 'dataids': [1]},
 {'annuli': [0], 'dataids': [0]}]

For this example we are going to see what happens when we tie together the temperature and abundance of annulus 17 and 18 (the inner-most annulus is numbered 0, the outer-most is 19 in this example).

>>> dep.tie_par('xsmekal.kt', 17, 18)
Tying xsmekal_18.kT to xsmekal_17.kT
>>> dep.tie_par('xsmekal.abundanc', 17, 18)
Tying xsmekal_18.Abundanc to xsmekal_17.Abundanc

Looking at the first few elements of the list returned by get_shells shows us that the two annuli have now been grouped together. This means that when a fit is run then these two annuli will be fit simultaneously.

>>> dep.get_shells()[0:5]
[{'annuli': [19], 'dataids': [19]},
 {'annuli': [17, 18], 'dataids': [17, 18]},
 {'annuli': [16], 'dataids': [16]},
 {'annuli': [15], 'dataids': [15]},
 {'annuli': [14], 'dataids': [14]}]

Since the model values for the 18th annulus have been changed, the data needs to be re-fit. The screen output is slightly different, as highlighted below, to reflect this new grouping:

>>> dep.fit()
Note: annuli have been tied together
Fitting annulus 19  dataset: 19
Dataset               = 19
...
Freezing xswabs_19
Freezing xsmekal_19
Fitting annuli [17, 18]  datasets: [17, 18]
Datasets              = 17, 18
Method                = levmar
Statistic             = chi2xspecvar
Initial fit statistic = 469.803
Final fit statistic   = 447.963 at function evaluation 16
Data points           = 439
Degrees of freedom    = 435
Probability [Q-value] = 0.323566
Reduced statistic     = 1.0298
Change in statistic   = 21.8398
   xsmekal_17.kT   2.64505      +/- 0.0986876
   xsmekal_17.Abundanc   0.520943     +/- 0.068338
   xsmekal_17.norm   0.0960039    +/- 0.00370363
   xsmekal_18.norm   0.0516096    +/- 0.00232468
Freezing xswabs_17
Freezing xsmekal_17
Freezing xswabs_18
Freezing xsmekal_18
Fitting annulus 16  dataset: 16
...

The table returned by fit() and get_fit_results() still contains 20 rows, since each annulus has its own row:

>>> tied = dep.get_fit_results()
>>> len(tied)
20

Looking at the outermost annuli, we can see that the temperature and abundance values are the same:

>>> print(tied['xsmekal.kT'][16:])
    xsmekal.kT
------------------
2.5470334673573936
 2.645045455995055
 2.645045455995055
2.7480050952727346
>>> print(tied['xsmekal.Abundanc'][16:])
  xsmekal.Abundanc
-------------------
0.24424232427701525
 0.5209430839965413
 0.5209430839965413
 0.4297577517591726

Note

The field names in the table are case sensitive - that is, you have to say xsmekal.kT and not xsmekal.kt - unlike the Sherpa parameter interface.

The results can be compared to the original, such as in the following plot, which shows that the temperature profile hasn’t changed significantly in the core, but some of the oscillations at large radii have been reduced (the black circles show the temperature values from the original fit):

>>> dep.fit_plot('xsmekal.kt', units='pc')
>>> rmid = (onion['rlo_phys'] + onion['rhi_phys']) / 2
>>> plt.scatter(rmid.to(u.pc), onion['xsmekal.kT', c='k')

_images/m87_temperature_tied_comparison.png

Note

As the X-axis of the plot created by fit_plot() is in pc, the rmid value is expliclty converted to these units before plotting.

The results parameter of fit_plot could have been set to the variable onion to display the results of the previous fit, but there is no easy way to change the color of the points.

More shells could be tied together to try and further reduce the oscillations, or errors calculated, with either conf() or covar().

The parameter ties can be removed with untie_par():

>>> dep.untie_par('xsmekal.kt', 18)
Untying xsmekal_18.kT
>>> dep.untie_par('xsmekal.abundanc', 18)
Untying xsmekal_18.Abundanc
>>> dep.get_shells()[0:4]
[{'annuli': [19], 'dataids': [19]},
 {'annuli': [18], 'dataids': [18]},
 {'annuli': [17], 'dataids': [17]},
 {'annuli': [16], 'dataids': [16]}]

Multiple datasets per annulus (3C186)¶

A second example illustrates the use of deproject for a multi-obsid observation of 3C186. It also shows how to set a background model for fitting with the cstat statistic. The extracted spectral data for this example are not yet publicly available, but were used in Siemiginowska et al. 2010:

The script starts with some setup:

>>> import deproject

>>> radii = ('2.5', '6', '17')
>>> dep = deproject.Deproject(radii=[float(x) for x in radii])

>>> set_method("levmar")
>>> set_stat("cstat")

Now we read in the data as before with dep.load_pha(). The only difference to the M87 example is the additional loop over the obsids. The dep.load_pha() function automatically extracts the obsid (the identifier used to discriminate Chandra observations) from the file header. This is used later in the case of setting a background model.

>>> obsids = (9407, 9774, 9775, 9408)
>>> for ann in range(len(radii)  -1):
...     for obsid in obsids:
...         dep.load_pha('3c186/%d/ellipse%s-%s.pi' % (obsid, radii[ann], radii[ann+1]), annulus=ann)

Create and configure the source model expression as usual:

>>> dep.set_source('xsphabs*xsapec')
>>> dep.ignore(None, 0.5)
>>> dep.ignore(7, None)
>>> dep.freeze("xsphabs.nh")

>>> dep.set_par('xsapec.redshift', 1.06)
>>> dep.set_par('xsphabs.nh', 0.0564)

Set the background model (here we use the IPython %run magic command to evaluate the commands in the file acis-s-bkg.py):

>>> %run acis-s-bkg.py
>>> acis_s_bkg = get_bkg_source()
>>> dep.set_bkg_model(acis_s_bkg)

Fit the projection model:

>>> dep.fit()

The deproject.deproject module¶

Deproject 2-d circular annular spectra to 3-d object properties.

This module implements the “onion-skin” approach popular in X-ray analysis of galaxy clusters and groups to estimate the three-dimensional temperature, metallicity, and density distributions of an optically-thin plasma from the observed (projected) two-dimensional data, arranged in concentric circular annuli.

Copyright:	Smithsonian Astrophysical Observatory (2009, 2019)
Author:	Tom Aldcroft (taldcroft@cfa.harvard.edu), Douglas Burke (dburke@cfa.harvard.edu)

Classes

Deproject(radii[, theta, angdist, cosmology]) Support deprojecting a set of spectra (2-d concentric circular annuli).

Functions

deproject_from_xflt(pat, rscale[, rinner, …]) Set up the projection object from XFLT keywords in the PHA files.

Class Inheritance Diagram¶

The deproject.specstack module¶

Manipulate a stack of spectra in Sherpa.

The methods in the SpecStack class provide a way to automatically apply familiar Sherpa commands such as set_par or freeze or plot_fit to a stack of PHA spectra. This simplifies simultaneous fitting of multiple spectra.

Note that the specstack module is currently distributed within with the deproject package. Specstack is not yet fully documented or tested outside the context of deproject.

Copyright:	Smithsonian Astrophysical Observatory (2009, 2019)
Author:	Tom Aldcroft (taldcroft@cfa.harvard.edu), Douglas Burke (dburke@cfa.harvard.edu)

Classes

SpecStack() Manipulate a stack of spectra in Sherpa.