Obtain and Fit a Radial Profile

Get Started

Download the sample data: 1838 (ACIS-S, G21.5-0.9)

unix% download_chandra_obsid 1838 evt2

In the following examples, restrict the energy range of the events:

unix% dmcopy "acisf01838N003_evt2.fits[energy=300:8000]" acis_1838_evt2.fits

Creating Radial Profiles

The ability of dmextract to operate on a stack of regions makes it possible to compute radial profiles simply by defining multiple concentric annuli.

1. Creating Multiple Annuli

Display the file:

unix% ds9 acis_1838_evt2.fits &

Select Region → Shape → Annulus and left-click on the image. A singular annular region will appear. To edit the region, make it active (left-click) and select "Get Info..." from the Region menu.

A region editing window (Figure 1) will appear, in which one can adjust the number of annuli and their sizes. Thirty-eight equally-spaced annuli, with minimium and maximum of 10 and 200 pixels respectively, which are located around (but exclude) the core of G21.5-09, are shown in Figure 2.

Figure 1: Annulus editing window

[Version: full-size]

Figure 1: Annulus editing window

The parameters are set to create 38 equally-spaced annuli centered at (4062,4240).

Figure 2: Source annuli overlaid on data

[Version: full-size]

Figure 2: Source annuli overlaid on data

The minimum radius is 10 pixels and the maximum radius is 200 pixels.

Save the annuli:

1. Region → Save Regions... → Save As "annuli.reg".
2. After choosing "OK" in the region filename dialog, a format dialog is opened. Set the format to "CIAO" and the coordinate system to "Physical".

Follow similar steps to create a a background annulus (Figure 3) from 200 to 225 pixels. The background is saved as "annuli_bgd.reg".

Figure 3: Background annulus overlaid on data

[Version: full-size]

Figure 3: Background annulus overlaid on data

The background region has been defined as an annulus with inner radius of 200 pixels and outer radius of 225 pixels.

The source region file looks like this:

unix% more annuli.reg
# Region file format: CIAO version 1.0
annulus(4062,4240,10,15)
annulus(4062,4240,15,20)
annulus(4062,4240,20,25)
.
. (etc.)
.
annulus(4062,4240,190,195)
annulus(4062,4240,195,200)

and the background annulus like this:

unix% more annuli_bgd.reg
# Region file format: CIAO version 1.0
annulus(4062,4240,200,225)

---

2. Removing Contaminating Point Sources

Suppose that the annuli had a maximum radius of 250 pixels in the previous step. The point source circled in green in Figure 4 would then contribute to a few of the radial profiles.

Figure 4: Annuli that contain an unwanted point source

[Version: full-size]

Figure 4: Annuli that contain an unwanted point source

The green circle shows the source that needs to be removed.

Having saved the region in ds9:

unix% more contam.reg
# Region file format: CIAO version 1.0
circle(4235.9495,4084.4343,8)

it is easy to remove this point source before generating the radial profiles:

unix% dmcopy "acis_1838_evt2.fits[exclude sky=region(contam.reg)]" acis_1838_excl_evt2.fits

This command creates a new event file with the point source removed (Figure 5). Use this event file in the rest of the radial profile analysis. This is not an issue in this example, so we continue using acis_1838_evt2.fits.

Figure 5: Event file with source removed

[Version: full-size]

Figure 5: Event file with source removed

The green circle shows where the unwanted source used to be located.

---

3. Run dmextract

It is now possible to run dmextract to extract the radial profiles:

unix% punlearn dmextract
unix% pset dmextract infile="acis_1838_evt2.fits[bin sky=@annuli.reg]"
unix% pset dmextract outfile=1838_rprofile.fits
unix% pset dmextract bkg="acis_1838_evt2.fits[bin sky=@annuli_bgd.reg]"
unix% pset dmextract opt=generic
unix% dmextract
Input event file  (acis_1838_evt2.fits[bin sky=@annuli.reg]):
Enter output file name (1838_rprofile.fits):

The contents of the parameter file may be checked using plist dmextract.

The tool calculates several new columns, the surface brightness (SUR_BRI) and its error (SUR_BRI_ERR) among them:

unix% dmlist 1838_rprofile.fits cols

--------------------------------------------------------------------------------
Columns for Table Block HISTOGRAM
--------------------------------------------------------------------------------

ColNo  Name                 Unit        Type             Range
.
. (output omitted)
.
  20   NET_COUNTS           count        Real8          -Inf:+Inf            Net Counts
  21   NET_ERR              count        Real8          -Inf:+Inf            Error on Net Counts
  22   NET_RATE             count/s      Real8          -Inf:+Inf            Net Count Rate
  23   ERR_RATE             count/s      Real8          -Inf:+Inf            Error Rate
  24   SUR_BRI              count/pixel**2 Real8          -Inf:+Inf            Net Counts per square pixel
  25   SUR_BRI_ERR          count/pixel**2 Real8          -Inf:+Inf            Error on net counts per square pixel
.
.
.

SUR_BRI is calculated as NET_COUNTS/AREA (columns 19 and 7, respectively); SUR_BRI_ERR is NET_ERR/AREA (columns 20 and 7).

Note that since the surface brightness is calculated from the NET_COUNTS column, the background counts are already removed from it: NET_COUNTS = COUNTS - [(BG_COUNTS/BG_AREA) * AREA]. It is therefore not necessary to account for the background separately when fitting this data.

Profile in Flux Units

Suppose you are interested in a radial profile in units of flux (photons/cm²/pix²/s or photons/cm²/arcsec²/s). Exposure-corrected images (e.g. the "flux" image file produced by fluximage and merge_obs) should NOT be used; while the surface brightness profile will be correctly determined, the uncertainties will be incorrect. Instead, in order to appropriately propagate the errors with dmextract, an input image in units of counts (or an input energy-filtered event file) should be used, and in addition an exposure map supplied to provide the counts-to-flux conversion.

unix% fluximage infile=acis_1838_evt2.fits outroot=1838 bands="0.3:8.0:2.5" binsize=1

unix% dmextract infile="1838_0.3-8.0_thresh.img[bin sky=@annuli.reg]" \
? outfile=1838_flux_rprofile.fits \
? bkg="1838_0.3-8.0_thresh.img[bin sky=@annuli_bgd.reg]" \
? exp=1838_0.3-8.0_thresh.expmap \
? bkgexp= 1838_0.3-8.0_thresh.expmap \
? opt=generic

— or extracting from an events file:

unix% dmextract infile="acis_1838_evt2.fits[energy=300:8000][bin sky=@annuli.reg]" \
? outfile=1838_flux_rprofile.fits \
? bkg="acis_1838_evt2.fits[energy=300:8000][bin sky=@annuli_bgd.reg]" \
? exp=1838_0.3-8.0_thresh.expmap \
? bkgexp= 1838_0.3-8.0_thresh.expmap \
? opt=generic

The tool calculates several new columns, including the surface brightness in flux units, (SUR_FLUX) and its error (SUR_FLUX_ERR), where the surface brightness is calculated from the NET_FLUX column with background accounted for: NET_FLUX = FLUX - (AREA / BG_AREA) * BG_COUNTS / (EXPOSURE * MEAN_BG_EXP).

The CEL_FLUX and CEL_FLUX_ERR columns are also provided, with transforms applied to convert values in pixels to ones in arcsec (if the appropriate WCS information is provided in the input file).

New in CIAO 4.11

dmextract in CIAO 4.11 now creates the RMID column. It is the average of the annulus radii.

unix% dmlist 1838_rprofile.fits'[cols R,RMID]' data

--------------------------------------------------------------------------------
Data for Table Block HISTOGRAM
--------------------------------------------------------------------------------

ROW    R[2]                                     RMID

     1                [       10.0        15.0]                12.50
     2                [       15.0        20.0]                17.50
     3                [       20.0        25.0]                22.50
     4                [       25.0        30.0]                27.50
     5                [       30.0        35.0]                32.50
...

Plotting and Fitting

The radial profile can now be plotted in python using matplotlib

unix% python

>>> from pycrates import read_file
>>> import matplotlib.pylab as plt
>>>
>>> tab = read_file("1838_rprofile.fits")
>>> xx = tab.get_column("rmid").values
>>> yy = tab.get_column("sur_bri").values
>>> ye = tab.get_column("sur_bri_err").values
>>>
>>> plt.errorbar(xx,yy,yerr=ye, marker="o")
>>> plt.xscale("log")
>>> plt.yscale("log")
>>> plt.xlabel("R_MID (pixel)")
>>> plt.ylabel("SUR_BRI (photons/cm**2/pixel**2/s)")
>>> plt.title('G21.5-0.9 [Chip S3, T=120, Offsets=-1,0,1]')
>>> plt.show()

which produces Figure 6.

Figure 6: Radial profile of the source

[Version: full-size]

Figure 6: Radial profile of the source

A model can be fit to the measured surface brightness profile using Sherpa. As mentioned before, the background counts are already removed from the surface brightness, so it is not necessary to account for the background separately when fitting the data.

unix% sherpa

sherpa> load_data(1,"1838_rprofile.fits", 3, ["RMID","SUR_BRI","SUR_BRI_ERR"])

sherpa> set_source("beta1d.src")
sherpa> src.r0 = 105
sherpa> src.beta = 4
sherpa> src.ampl = 0.00993448

sherpa> freeze(src.xpos)

sherpa> fit()
Dataset               = 1
Method                = levmar
Statistic             = chi2
Initial fit statistic = 4.24792e+11
Final fit statistic   = 221.594 at function evaluation 30
Data points           = 38
Degrees of freedom    = 35
Probability [Q-value] = 5.65854e-29
Reduced statistic     = 6.33127
Change in statistic   = 4.24792e+11
   src.r0         134.119      +/- 0
   src.beta       4.6542       +/- 0.0351998
   src.ampl       1.51907e-06  +/- 1.70734e-08

sherpa> plot_fit(xlog=True, ylog=True)

which produces Figure 7.

Figure 7: Fit to the radial profile of the source

[Version: full-size]

Figure 7: Fit to the radial profile of the source

Note that the effects of 2D blurring in a 2D image cannot be reproduced by convolving the radial profile of the PSF with a profile of the model. See "Accounting for PSF Effects in 2D Image Fitting".

Aperture Corrections

As is discussed in the Radial Profile Aperture Correction page, there has been an assumption made in this analysis: the data in each of the annuli are statistically independent. This is not entirely true. Due to the PSF, we expect that some counts that are imaged in one annnuli to belong to a neighboring region. We can use the mkrprm script to estimate the fraction of the counts we expect to come from neighboring regions and then use the MatrixModel in sherpa to include it when fitting the radial profile.

First we run mkrprm to create the spatial equivalent of a response matrix.

unix% mkrprm infile=acis_1838_0.3-8_thresh.img region=@annuli.reg outfile=mkrprm.img

The tool may run for several minutes depending on the size of the image and the number of radii. In this example the default psfmethod=map has been used which will approximate the aperture correction my using the PSF map to estimate the PSF fractional contribution in each region from all other regions.

The output from mkrprm is a 2D image. The x-axis is the input region number and the y-axis is the output region number. The pixel values are the fraction of the total flux that the input region contributes to each of the other regions. Values will be between 0 and 1. The output can be visualized with DS9 as shown in Figure 8

Figure 8: Radial Profile Aperture correction matrix

[Version: full-size]

Figure 8: Radial Profile Aperture correction matrix

This image shows the output from the mkrprm tool as displayed using ds9 with the command:

unix% ds9 mkrprm.img -view colorbar yes -cmap cubehelix1 -log -zoom to fit

The output is nearly identically diagonal. That is, the counts in each radial bin are larger independent from neighboring regions. However, there is some spread on each side of the diagonal indicating that a few percent of the counts in each radial bin are attributed to neighboring bins. Given the low amplitudes we expect this to be a small effect.

We will now repeat the fit of the radial profile and include the radial profile aperture correction matrix. To do this we will use the MatrixModel which is part of the sherpa_contrib module:

unix% sherpa

sherpa> load_data(1,"1838_rprofile.fits", 3, ["RMID","SUR_BRI","SUR_BRI_ERR"])
sherpa> x_vals = get_data(1).x

sherpa> from sherpa_contrib.matrix_model import MatrixModel
sherpa> mymatrix = read_file("mkrprm.img").get_image().values
sherpa> mymaxtrix_rsp = MatrixModel(mymatrix, x_vals)

sherpa> set_source(mymaxtrix_rsp("beta1d.src"))

sherpa> src.r0 = 105
sherpa> src.beta = 4
sherpa> src.ampl = 0.00993448

sherpa> freeze(src.xpos)

sherpa> fit()
Dataset               = 1
Method                = levmar
Statistic             = chi2
Initial fit statistic = 4.22307e+11
Final fit statistic   = 220.503 at function evaluation 30
Data points           = 38
Degrees of freedom    = 35
Probability [Q-value] = 9.00736e-29
Reduced statistic     = 6.3001
Change in statistic   = 4.22307e+11
   src.r0         130.225      +/- 0
   src.beta       4.4451       +/- 0.0335157
   src.ampl       1.53312e-06  +/- 1.72717e-08

The result of the fit shows that the width of the 1D Beta model has decreased from 4.65 +/- 0.035 pixels to 4.45 +/- 0.034 pixels. This aligns with our exceptions since the effect of the aperture correction from the PSF is to blur out the data, so by including the response we have removed the blur from the beta model. Also the amplitude of the difference is very small which we expect given the nearly identically diagonal nature of the response. In this specific example it may not have been necessary to include the aperture corrections; however, for some choice of regions it may become significant in the analysis.