Example: Deprojection, gas density profiles and gas masses

This thread shows how to use PyProffit to extract gas density profiles, gas masses, integrated count rates and luminosities from imaging data. This code implements the method presented in Eckert et al. (2020) to apply deprojection and PSF deconvolution.

We start by loading the data and extracting a profile…

import numpy as np
import pyproffit
import matplotlib.pyplot as plt

import os

# Change this to the proper directory containing the test data
os.chdir('../../validation/')

Let’s load the data into a Data structure…

dat=pyproffit.Data(imglink='b_37.fits.gz',explink='expose_mask_37.fits.gz',
                   bkglink='back_37.fits.gz')

WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' / Equatorial system reference
the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]

Now we extract a profile in circular annuli from the surface brightness peak…

prof=pyproffit.Profile(dat,center_choice='peak',maxrad=45.,binsize=20.)

prof.SBprofile()

prof.Plot()

Determining X-ray peak
Coordinates of surface-brightness peak: 272.0 281.0
Corresponding FK5 coordinates:  55.72147733144434 -53.628226297404545

<Figure size 432x288 with 0 Axes>

PSF modeling

To correct for PSF smearing, we need to create a PSF mixing matrix which describes the leaking of photons from each annulus to the others. The shape of the PSF needs to be known a priori. The definition of the PSF can be done either by providing an input image or through an analytical function.

In the following example we model the ROSAT/PSPC PSF as a King function with a core radius of 25 arcsec

def fking(x):
    r0=25./60. # arcmin
    alpha=1.5 # King slope
    return np.power(1. + (x/r0) ** 2, -alpha)

Now we pass this function to the PSF method of the Profile class. To create the mixing matrix, the method creates normalized images of each annulus individually, convolves them with the PSF using FFT and computes the fraction of events in the FFT-convolved image that fall within each annulus. See Eckert et al. (2020) for more details.

prof.PSF(psffunc = fking)

Let’s inspect the PSF mixing matrix. Each row describes the fraction originating from that annulus which are in fact recorded in any other row

fig = plt.figure(figsize=(20,20))
plt.imshow(np.log10(prof.psfmat), aspect='auto')

<matplotlib.image.AxesImage at 0x7f8ab85e3310>

Deprojection

For deprojection, PyProffit implements both the standard onion-peeling deprojection similar to that of plain Proffit and a new method based on multiscale decomposition of the observed profile. The new method is suitable in the low count-rate regime, provides on-the-fly propagation of the background value and PSF deconvolution.

The two methods can be accessed with the Deproject class. For the extraction of density profiles the class requires the following input:

• profile=prof: a Profile object containing the data
• z=redshift: the redshift of the source
• cf=factor: the conversion between count rate and emission measure

The conversion factor can be computed using the calc_emissivity function, which goes through XSPEC to simulate an absorbed APEC model with an emission measure of 1 and retrieve the corresponding count rate. XSPEC needs to be accessible in the PATH for this command to work properly.

Note: Since we have corrected from vignetting using the exposure map, the count rates are all rescaled to equivalent on-axis count rates. Therefore the on-axis effective area of the telescope must be used to convert correctly count rates into emission measure.

Note: In some cases (in particular Chandra), the exposure maps are usually provided in units of cm\(^2\) s instead of just s, such that the unit of the surface brightness measured by PyProffit is phot/cm\(^2\)/s instead of count/s. In this case, the type=”cr” option should be used.

z_a3158 = 0.059 # Source redshift, here 0.059 for the test cluster A3158
kt_a3158 = 5.0 # Plasma temperature; if a soft band is used the profile is mildly dependent on it
nh_a3158 = 0.0118 # Source NH in 1e22 cm**(-2) unit
rsp = 'pspcb_gain2_256.rsp' # Response file, here ROSAT/PSPC in RSP format
elow = 0.42 # Lower energy boundary of the image
ehigh = 2.01 # Upper energy boundary of the image

cf = prof.Emissivity(z=z_a3158,
                    kt=kt_a3158,
                    nh=nh_a3158,
                    rmf=rsp,
                    elow=elow,
                    ehigh=ehigh)

print(cf)

43.53

We are now ready to declare the Deproject object. Note that in case the redshift and conversion factor are not known, it is still possible to run the PSF deconvolution and profile reconstruction, however the gas density profile and gas mass cannot be computed.

depr = pyproffit.Deproject(z=z_a3158, cf=cf, profile=prof)

Let’s start with the multiscale decomposition method. It can be launched with the Multiscale method of the Deproject class. The parameters of the method are the following:

• backend=’pymc3’: choose whether the optimization will be performed with the PyMC3 or the Stan backend
• nmcmc=1000: number of points in Hamiltonian Monte Carlo chain
• tune=500: number of tuning steps to set up the NUTS sampler
• bkglim=rad: radius beyond which it is assumed that the source is 0 (i.e. background only)
• samplefile=file.dat: output file where the HMC samples will be stored

The sampling time with HMC will depend on a number of factors, including the number of bins in the profile, the number of points in the output chain, and the bkglim value.

depr.Multiscale(nmcmc=1000, tune=1000, bkglim=30.)

Running MCMC...

100.00% [331/331 00:00<00:00 logp = -964.13, ||grad|| = 0.021414]

/home/deckert/.local/lib/python3.8/site-packages/pyproffit/deproject.py:628: FutureWarning: In v4.0, pm.sample will return an `arviz.InferenceData` object instead of a `MultiTrace` by default. You can pass return_inferencedata=True or return_inferencedata=False to be safe and silence this warning.
  trace = pm.sample(nmcmc, tune=tune, start=tm)
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [bkg, coefs]

100.00% [8000/8000 02:50<00:00 Sampling 4 chains, 3,374 divergences]

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 171 seconds.
There were 738 divergences after tuning. Increase `target_accept` or reparameterize.
There were 841 divergences after tuning. Increase `target_accept` or reparameterize.
There were 861 divergences after tuning. Increase `target_accept` or reparameterize.
There were 934 divergences after tuning. Increase `target_accept` or reparameterize.
The rhat statistic is larger than 1.05 for some parameters. This indicates slight problems during sampling.
The estimated number of effective samples is smaller than 200 for some parameters.

Done.
 Total computing time is:  2.9870954672495524  minutes

The quality of the restruction can be viewed through the PlotSB method of the Deproject class

depr.PlotSB()

<Figure size 432x288 with 0 Axes>

Here the total model (PSF convolved source + background) is shown in orange. The reconstructed source profile (PSF deconvolved) is shown in blue, and the fitted sky background is in black. The residuals (bottom panel) allow the user to assess the quality of the reconstruction.

Count Rates and luminosities

The count rates can be computed easily from the reconstructed surface brightness profile within any user given apertures. This is done through the CountRate method, which integrates the PSF deconvolved model over the area. The posterior distribution of the count rate can be displayed as well.

cr, cr_lo, cr_hi = depr.CountRate(0., 30.)

Reconstructed count rate: 2.01158 (1.97752 , 2.04536)

<Figure size 432x288 with 0 Axes>

And the luminosity can be obtained similarly through the Luminosity method. Similarly to the gas density, the luminosity requires the emissivity conversion to be calculated

lum, lum_lo, lum_hi = depr.Luminosity(0., 30.)

Reconstructed luminosity: 2.01001e+44 (1.97598e+44 , 2.04376e+44)

<Figure size 432x288 with 0 Axes>

Gas density profile

Once a Multiscale reconstruction has been performed, and if the source redshift and the emission measure conversion factor have been provided, it is straightforward to measure the gas density profile of the source. This is done through the Density method of the Deproject class. The gas density profile can then be plotted through the PlotDensity method

depr.Density()

depr.PlotDensity()

<Figure size 432x288 with 0 Axes>

Onion Peeling deprojection

If instead of the multiscale approach one wishes to compute the deprojected profile using the classical onion peeling approach, in which case the projection kernel is directly inverted, the Deproject class contains the OnionPeeling method.

Note that in this case the background is not reconstructed on-the-fly, thus this method should be used directly on background subtracted profiles. Here we provide an example of the use of the OnionPeeling method. First, let us fit the surface brightness profile beyond 30 arcmin with a constant,

mod = pyproffit.Model(pyproffit.Const)

fitconst = pyproffit.Fitter(model=mod, profile=prof, fitlow=30., fithigh=40., bkg=-3.5)

fitconst.Migrad()

┌──────────────────────────────────┬──────────────────────────────────────┐
│ FCN = 47.4                       │              Nfcn = 22               │
│ EDM = 5.94e-06 (Goal: 0.0002)    │                                      │
├───────────────┬──────────────────┼──────────────────────────────────────┤
│ Valid Minimum │ Valid Parameters │        No Parameters at limit        │
├───────────────┴──────────────────┼──────────────────────────────────────┤
│ Below EDM threshold (goal x 10)  │           Below call limit           │
├───────────────┬──────────────────┼───────────┬─────────────┬────────────┤
│  Covariance   │     Hesse ok     │ Accurate  │  Pos. def.  │ Not forced │
└───────────────┴──────────────────┴───────────┴─────────────┴────────────┘
┌───┬──────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐
│   │ Name │   Value   │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit-  │ Limit+  │ Fixed │
├───┼──────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤
│ 0 │ bkg  │  -3.670   │   0.018   │            │            │         │         │       │
└───┴──────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘
┌─────┬──────────┐
│     │      bkg │
├─────┼──────────┤
│ bkg │ 0.000322 │
└─────┴──────────┘
Best fit chi-squared: 47.3985 for 135 bins and 134 d.o.f
Reduced chi-squared: 0.35372

prof.Plot(model=mod, axes=[30., 40., 5e-5, 1e-3])

<Figure size 432x288 with 0 Axes>

Now we define a new Profile object with a logarithmic binning, from which we will subtract the fitted background

prof2 = pyproffit.Profile(dat, center_choice='peak', binsize=30, maxrad=30., binning='log')

prof2.SBprofile()

prof2.Backsub(fitconst)

Determining X-ray peak
Coordinates of surface-brightness peak: 272.0 281.0
Corresponding FK5 coordinates:  55.72147733144434 -53.628226297404545

prof2.Plot()

<Figure size 432x288 with 0 Axes>

Now we are ready to define a new Deproject object and apply the OnionPeeling method,

depr_op = pyproffit.Deproject(profile=prof2, cf=cf, z=z_a3158)

depr_op.OnionPeeling()

depr_op.PlotDensity(xunit='arcmin')

<Figure size 432x288 with 0 Axes>

Clearly we want to know how the two methods compare. The plot_multi_methods function allows the user to easily compare the results of several density profile reconstructions

outfig = pyproffit.plot_multi_methods(deps = (depr, depr_op),
                            profs = (prof, prof2),
                            labels = ('Multiscale', 'Onion Peeling'))

Showing 2 density profiles

Clearly the two methods are consistent, but the Multiscale approach is much less noisy. In the central regions we can easily notice the effect of the PSF deconvolution in the Multiscale case; in the OnionPeeling case no PSF deconvolution can be applied.

If instead of PyMC3 one wishes to use the (usually more computationally efficient) Stan backend, it is easy to set it up when calling Multiscale by using the backend=’stan’ option. The results of PyMC3 and Stan are usually indistinguishable.

Gas masses

The gas mass profile is the integral of the gas density profile over the volume. The Deproject class contains two methods to compute the gas mass profile and the posterior distribution of \(M_{gas}\) evaluated at a specific radius. PlotMgas allows the user to view the total reconstructed \(M_{gas}\) profile, wherease Mgas computes \(M_{gas}\) at any user given radius (in kpc) and plot the posterior distribution of this value.

depr.PlotMgas()

Now let’s say for instance that we want to compute \(M_{gas}\) at \(R_{500}=1123 \pm 50\) kpc. The Mgas method evaluates the gas mass at the provided overdensity radius. The uncertainty in the overdensity radius can be propagated to the posterior \(M_{gas}\) distribution by randomizing the radius out to which the profile is integrated

mg_r500, mg_lo, mg_hi, rho = depr.Mgas(radius = 1123., radius_err=50.)

<Figure size 432x288 with 0 Axes>