Chandra X-Ray Observatory
	(CXC)
Skip to the navigation links
Last modified: 5 Jul 2017

URL: http://cxc.harvard.edu/ciao/threads/radial_profile/

Obtain and Fit a Radial Profile

CIAO 4.9 Science Threads


Overview

Synopsis:

The surface brightness flux is determined by finding the net counts in a stack of concentric annuli and then dividing by the respective areas. A specified analytic model may be fit to the resultant histogram. This information can be used, for instance, to provide evidence for extended emission and calculate the hardness ratio thereof.

Purpose:

To produce radial profiles of an HRC or ACIS imaging observation.

Related Links:

Last Update: 5 Jul 2017 - Added admonishment on obtaing the surface brightness profile in flux units.


Contents


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.

[Thumbnail image: The dialog box has fields for setting the center, inner and outer radius, and number of annuli for the region.]

[Version: full-size]

[Print media version: The dialog box has fields for setting the center, inner and outer radius, and number of annuli for the region.]

Figure 1: Annulus editing window

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

[Thumbnail image: The annuli are displayed in green on the data in ds9.]

[Version: full-size]

[Print media version: The annuli are displayed in green on the data in ds9.]

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".

[Thumbnail image: The annulus is displayed in green on the data is ds9.]

[Version: full-size]

[Print media version: The annulus is displayed in green on the data is ds9.]

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(4072,4246,10,15)
annulus(4072,4246,15,20)
annulus(4072,4246,20,25)
.
. (etc.)
.
annulus(4072,4246,190,195)
annulus(4072,4246,195,200)

and the background annulus like this:

unix% more annuli_bgd.reg
# Region file format: CIAO version 1.0
annulus(4070,4250,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.

Having saved the region in ds9:

unix% more contam.reg 
# Region file format: CIAO version 1.0
circle(4246.5,4090.5,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.

[Thumbnail image: The filtered event file is displayed in ds9.]

[Version: full-size]

[Print media version: The filtered event file is displayed in ds9.]

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.

[TIP]
Profile in Flux Units

Suppose you are interested in a radial profile in units of flux (photons/cm2/pix2/s or photons/cm2/arcsec2/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).

Finally, we want to add a column that defines the midpoint of the annular regions (rmid):

unix% punlearn dmtcalc
unix% pset dmtcalc infile=1838_rprofile.fits
unix% pset dmtcalc outfile=1838_rprofile_rmid.fits
unix% pset dmtcalc expression="rmid=0.5*(R[0]+R[1])"
unix% dmtcalc
Input file (1838_rprofile.fits): 
Output file (1838_rprofile_rmid.fits):
expression(s) to evaluate (rmid=0.5*(R[0]+R[1])):

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

The new column has been created in 1838_rprofile_rmid.fits:

unix% dmlist 1838_rprofile_rmid.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 using ChIPS:

unix% chips

chips> make_figure("1838_rprofile_rmid.fits[cols rmid,sur_bri,sur_bri_err]")
chips> log_scale(XY_AXIS)

which produces Figure 6.

Quit ChIPS before continuing:

chips> quit

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_rmid.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()
Statistic value = 233.766 at function evaluation 88
Data points = 38
Degrees of freedom = 35
Probability [Q-value] = 3.08377e-31
Reduced statistic  = 6.67902
   src.r0         121.817     
   src.beta       3.94547     
   src.ampl       4.47076     

sherpa> plot_fit()
sherpa> log_scale(XY_AXIS)

sherpa> limits(X_AXIS, 10, 200)
sherpa> limits(Y_AXIS, 0.0001, 10)

which produces Figure 7.

sherpa> quit

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".



Parameters for /home/username/cxcds_param/dmextract.par


#--------------------------------------------------------------------
#
# DMEXTRACT -- extract columns or counts from an event list
#
#--------------------------------------------------------------------
        infile = acis_1838_evt2.fits[bin sky=@annuli.reg] Input event file 
       outfile = 1838_rprofile.fits    Enter output file name
          (bkg = acis_1838_evt2.fits[bin sky=@annuli_bgd.reg]) Background region file or fixed background
        (error = gaussian)        Method for error determination(poisson|gaussian|<variance file>)
     (bkgerror = gaussian)        Method for background error determination(poisson|gaussian|<variance file>)
      (bkgnorm = 1.0)             Background normalization
          (exp = )                Exposure map image file
       (bkgexp = )                Background exposure map image file
      (sys_err = 0)               Fixed systematic error value for SYS_ERR keyword
          (opt = generic)         Output file type: pha1 
     (defaults = ${ASCDS_CALIB}/cxo.mdb -> /soft/ciao/data/cxo.mdb) Instrument defaults file
         (wmap = )                WMAP filter/binning (e.g. det=8 or default)
      (clobber = no)              OK to overwrite existing output file(s)?
      (verbose = 0)               Verbosity level
         (mode = ql)              
    


Parameters for /home/username/cxcds_param/dmtcalc.par


        infile = 1838_rprofile.fits    Input file
       outfile = 1838_rprofile_rmid.fits Output file
    expression = rmid=0.5*(R[0]+R[1]) expression(s) to evaluate
      (clobber = no)              Clobber output file if it exists?
      (verbose = 0)               Debug level
         (mode = ql)              
    

History

04 Jan 2005 updated for CIAO 3.2: version numbers
20 Dec 2005 updated for CIAO 3.3: default value of dmextract error and bkgerror parameters is "gaussian"
01 Dec 2006 updated for CIAO 3.4: ChIPS and Sherpa versions
22 Jan 2008 updated for CIAO 4.0: updated ChIPS and Sherpa syntax; kernel parameter removed from dmtcalc; filename and contam.reg file updated for reprocessed data (version N002 event file)
09 Feb 2009 updated for CIAO 4.1: images are inline; Python and S-Lang syntax for ChIPS and Sherpa sections
03 Apr 2009 new notes on 2D blurring on images
08 Feb 2010 updated for CIAO 4.2: changes to the ds9 region file format menu; ChIPS and Sherpa version
19 Jul 2010 the S-Lang syntax has been removed from this thread as it is not supported in CIAO 4.2 Sherpa v2.
31 Aug 2010 updated load_data command to use Data Model "cols" syntax
13 Jan 2011 updated for CIAO 4.3: add opt=generic to dmextract command to avoid a warning message
07 Mar 2012 reviewed for CIAO 4.4: changed load_data syntax to work in CIAO 4.4
03 Dec 2012 Review for CIAO 4.5; file name versions.
10 Dec 2013 Review for CIAO 4.6; no changes.
22 Dec 2014 Review for CIAO 4.7; added additional see also info.
05 Jul 2017 Added admonishment on obtaing the surface brightness profile in flux units.


Last modified: 5 Jul 2017
Smithsonian Institute Smithsonian Institute

The Chandra X-Ray Center (CXC) is operated for NASA by the Smithsonian Astrophysical Observatory. 60 Garden Street, Cambridge, MA 02138 USA.   Email:   cxchelp@head.cfa.harvard.edu Smithsonian Institution, Copyright © 1998-2017. All rights reserved.