Obtain and Fit a Radial Profile
[CIAO 3.4 Science Threads]
OverviewLast Update: 1 Dec 2006 - updated for CIAO 3.4: ChIPS and Sherpa versions 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 in Sherpa. This information can be used, for instance, to provide evidence for extended emission and calculate the hardness ratio thereof. Purpose: To produce radial profiles, then fit a model to them in Sherpa. Read this thread if: you would like to create a radial profile of an HRC or ACIS imaging observation. Related Links:
|
Contents
- Get Started
- Creating Radial Profiles
- Plotting and Fitting
- Parameter files:
- History
- Images
Get Started
Sample ObsID used: 1838 (ACIS-S, G21.5-09)
File types needed: evt2
In the following examples, restrict the energy range of the events:
unix% dmcopy "acisf01838N001_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 ![[Link to Image 1: ds9 region information/edit window]](../../imgs/imageicon.gif){kind=small}
will appear, in which one can adjust the number of
annuli and their sizes. Thirty-eight equal-radii
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 . We also created
a background annulus from
200 to 225 pixels.
Save the annuli:
- Create the annuli
- Region -> File Format-> Ciao
- Region -> File Coordinate System -> Physical
- Region -> Save Regions... -> Save As "annuli.reg"
Follow similar steps to create a file containing the background annulus, here named "annuli_bgd.reg".
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(4245,4094.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 . 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.
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 in Sherpa.
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, version CIAO 3.4 Copyright (C) 1999-2003, Smithsonian Astrophysical Observatory chips> plot "1838_rprofile_rmid.fits[cols rmid,sur_bri,sur_bri_err]" x 1 y 2 yerr 3 chips> log Warning: negative and zero values ignored in log scale
which produces Figure 6 . Exit
ChIPS before continuing:
chips> exit
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 Program
-----------------------------------------------------
Version: CIAO 3.4
Type AHELP SHERPA for overview.
Type EXIT, QUIT, or BYE to leave the program.
Notes:
Temporary files for visualization will be written to the directory:
/tmp
To change this so that these files are not deleted when you exit Sherpa,
edit $ASCDS_WORK_PATH in your 'ciao' setup script.
Abundances set to Anders & Grevesse
sherpa> read data 1 "1838_rprofile_rmid.fits[columns rmid,sur_bri]" FITSBIN
sherpa> read errors 1 "1838_rprofile_rmid.fits[columns rmid,sur_bri_err]" FITSBIN
sherpa> beta1d[sbr1]
sbr1.r0 parameter value [105]
sbr1.beta parameter value [1e-05]
sbr1.xpos parameter value [0]
sbr1.ampl parameter value [0.00993448]
sherpa> sbr1.ampl.max=10
sherpa> show sbr1
beta1d[sbr1] (integrate: off)
Param Type Value Min Max Units
----- ---- ----- --- --- -----
1 r0 thawed 105 1 197.5
2 beta thawed 1e-05 1e-05 10
3 xpos frozen 0 0 197.5
4 ampl thawed 9.9345e-03 9.9345e-05 10
sherpa> source=sbr1
sherpa> fit
LVMQT: V2.0
LVMQT: initial statistic value = 18548.3
LVMQT: final statistic value = 197.351 at iteration 25
sbr1.r0 116.969
sbr1.beta 3.67579
sbr1.ampl 4.50021
sherpa> lplot fit
sherpa> log
Warning: negative and zero values ignored in log scale
sherpa> limits y 0.0001 10
sherpa> limits x 10 200
sherpa> redraw
which produces Figure 7 .
sherpa> exit Goodbye.
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
(kernel = default) Data Model creation/copy kernel
(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 |
