Skip to the navigation links
Search the CIAO website or contact the CXC HelpDesk
Last modified: 31 Aug 2010
Where are the PDFs?
Return to Threads Page: Top | All | Imag

Obtain and Fit a Radial Profile

CIAO 4.2 Science Threads


Overview

Last Update: 31 Aug 2010 - updated load_data command to use Data Model "cols" syntax

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:



Contents


Get Started

Sample ObsID used: 1838 (ACIS-S, G21.5-0.9)

File types needed: evt2

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

unix% dmcopy "acisf01838N002_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

[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).

Figure 2: Source annuli overlaid on data

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

Figure 3: Background annulus overlaid on data

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

Figure 4: Annuli that contain an unwanted point source

[Thumbnail image: The sourceannuli are displayed in blue on the data in ds9.]

[Version: full-size]

[Print media version: The sourceannuli are displayed in blue on the data in ds9.]

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

Figure 5: Event file with source removed

[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% 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.

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
-----------------------------------------
Welcome to ChIPS: CXC's Plotting Package
-----------------------------------------
CIAO 4.2 Tuesday, July 6, 2010

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

Figure 6: Radial profile of the source

[Thumbnail image: The radial profile is plotted on a log-log scale with cross points and white line connecting them.]

[Version: full-size]

[Print media version: The radial profile is plotted on a log-log scale with cross points and white line connecting them.]

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
-----------------------------------------------------
Welcome to Sherpa: CXC's Modeling and Fitting Package
-----------------------------------------------------
CIAO 4.2 Sherpa version 2 Tuesday, July 6, 2010

sherpa> load_data(1, "1838_rprofile_rmid.fits[cols 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

Figure 7: Fit to the radial profile of the source

[Thumbnail image: A log-log plot of x vs y with white datapoints and the fit overlaid in red.]

[Version: full-size]

[Print media version: A log-log plot of x vs y with white datapoints and the fit overlaid in red.]

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



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 = pha1)            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
Return to Threads Page: Top | All | Imag

Where are the PDFs?
Last modified: 31 Aug 2010