Obtain and Fit a Radial Profile
![[CXC Logo]](../../imgs/cxc-logo.gif)
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:
- Analysis Guide: HRC Imaging
- Analysis Guide: Extended Sources
Contents
- Get Started
- Creating Radial Profiles
- Plotting and Fitting
- Parameter files:
- History
-
Images
- Figure 1: Annulus editing window
- Figure 2: Source annuli overlaid on data
- Figure 3: Background annulus overlaid on data
- Figure 4: Annuli that contain an unwanted point source
- Figure 5: Event file with source removed
- Figure 6: Radial profile of the source
- Figure 7: Fit to the radial profile of the source
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.
![[Thumbnail image: The dialog box has fields for setting the center, inner and outer radius, and number of annuli for the region.]](edit.thumb300.png)
[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.]](edit.png)
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.]](source.thumb300.png)
[Version: full-size]
![[Print media version: The annuli are displayed in green on the data in ds9.]](source.png)
Figure 2: Source annuli overlaid on data
The minimum radius is 10 pixels and the maximum radius is 200 pixels.
Save the annuli:
- Region → Save Regions... → Save As "annuli.reg".
- 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.]](bg.thumb300.png)
[Version: full-size]
![[Print media version: The annulus is displayed in green on the data is ds9.]](bg.png)
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.
![[Thumbnail image: The sourceannuli are displayed in blue on the data in ds9.]](ptsrc.thumb300.png)
[Version: full-size]
![[Print media version: The sourceannuli are displayed in blue on the data in ds9.]](ptsrc.png)
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.
![[Thumbnail image: The filtered event file is displayed in ds9.]](clean.thumb300.png)
[Version: full-size]
![[Print media version: The filtered event file is displayed in ds9.]](clean.png)
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
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
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 |