/ ciao4.17 / ahelp / mkrprm.html

Jump to: Description · Examples · Parameters · Algorithm: psfmethod · Including the Response Matrix when fitting with Sherpa · Changes in the scripts 4.17.0 (Decemeber 2024) release · About Contributed Software · Bugs · See Also

AHELP for CIAO 4.17

mkrprm

Context: Tools::Response

Synopsis

Compute the aperture correction for radial profiles

Syntax

mkrprm  infile regions outfile [psfmethod] [ecf] [function] [energy]
[flux] [random_seed] [marx_root] [tmpdir] [clobber] [verbose]

Description

'mkrprm' will approximate the contribution of the PSF scattering in each region from all regions. Phrased slightly differently: it will compute an approximate aperture correction for each region.

It is generally assumed that when creating radial profiles that the data in each radial bin are independent. For example if extracting counts in a series of annulii: the counts in each annular ring are assumed to be independent of all other annuli. Typically to "ensure" that this assumption is valid users will try to select regions/annuli that are (much) larger than the size of the Point Spread Function (PSF).

This assumption of independence is not always valid due to the variable size of the Chandra PSF. Depending on how far off-axis the data are imaged, it may be the case that there is some fraction of the counts in one annular ring that is scattered from the emission imaged in a different region.

This is the spatial equivalent of the standard spectral response/redistribution matrix, a.k.a. the "RMF". The output is in fact a matrix that provides an approximation of the fraction of the aperture convolved with PSF that is collected in each bin from all bins.

Details of the PSF convolution are discussed in the Algorithms section below.

Examples

Example 1

dmextract "img.fits[bin sky=@ciao.reg]" radial.prof op=generic
mkrprm img.fits @ciao.reg matrix.img"[opt kernel=text/dtf]"

We start with a simple example with only five (5) annulii. The stack of regions looks like:

$ cat ciao.reg
annulus(4100,4123,0,2)
annulus(4100,4123,2,4)
annulus(4100,4123,4,6)
annulus(4100,4123,6,8)
annulus(4100,4123,8,10)

First we extract the radial profile in the usual way using dmextract. Then we compute the radial profile response matrix using mkrprm. The inputs are the same input image and stack of regions. For this example, we have specified that the output file, matrix.img, should be an ASCII file using the "[opt kernel=text/dtf]" datamodel option. This is only done is this example so that we can more easily display the values in this help file. In this example the output looks like the following. Note: we have manually aligned the columns to be easier to read.

$ cat matrix.img
...
END

8.0967010351261e-01 1.9030467774856e-01 2.5218738826281e-05 0.0                 0.0
6.8721133631426e-02 8.0291679360828e-01 1.2835295705841e-01 9.1157018863392e-06 0.0
5.1225563240883e-06 7.2198538345354e-02 8.2203443497993e-01 1.0575165397976e-01 1.0250138635154e-05
0.0                 3.9067293798596e-06 8.0572688746482e-02 8.1328640656448e-01 1.0613308357164e-01
0.0                 0.0                 5.4667406054155e-06 7.4293158500145e-02 8.3350501980292e-01

The first row corresponds to the first region in the stack. We see that 80.967% of the PSF convolved with first annulus is imaged in the first annuls. The 2nd column represents the 2nd annulus, so we see that 19.03% of the first annulus convolved with the PSF is actually imaged in the second annulus. And then there is 0.00252% of the PSF convolved with the first annuls imaged in the third annulus. The sum of all the values in the first row equals 1; that is 100% of the PSF is imaged in those three annulii.

And so on for the second row: 6.87% of the second annulus is imaged in the 1st, 80.29% of the 2nd annulus is imaged in the 2nd annulus, 12.84% of the 2nd annulus is imaged in the 3rd, and 0.00091% of the 2nd annulus is imaged in the 4th. Again, the values sum to 1.0.

The matrix is nearly diagonal as we expect. Most of the PSF contribution originating in any region is actually imaged in that region. It is worth noting that the last two rows do not sum to 1.0. The PSF scatter in the outer regions extends beyond those regions so that the total aperture correction is less than 1.0 for those regions.

Example 2

mkrprm img.fits @ciao.reg matrix.img
ds9 matrix.img -zoom to fit

The same example as before but this time the output file will be a FITS image that can be display with ds9 or with matplotlib's "imshow()" command.

Parameters

name	type	ftype	def	min	max	units	reqd	stacks
infile	file	input					yes
regions	string						yes	yes
outfile	file	output					yes
psfmethod	string		map
ecf	real		0.393	0	1
function	string		gaus
energy	real					keV
flux	real					photon/cm**2/sec
random_seed	integer		-1
marx_root	file
tmpdir	file		${ASCDS_WORK_PATH}
clobber	boolean		no
verbose	integer		1	0	5

Detailed Parameter Descriptions

Parameter=infile (file required filetype=input)

Input image file.

The image used to generate the radial profile. If users used an event file, then it must be binned into an image. The image can be filtered (cropped) to just include the area around the extraction region; this will make the script run faster. If cropped, then it is suggest to be slightly larger extraction region to allow for correct PSF normalization in the outer most regions.

Parameter=regions (string required stacks=yes)

The stack of extraction regions used to generate the radial profile.

The stack of extraction regions should be in CIAO region form.

Using dmextract's annulus

The dmextract tool has a special syntax for a stack of annulus shaped regions that is not supported by the rest of CIAO.

$ dmextract img.fits"[bin sky=annulus(4100,4123,0:10:2)]" rad.prof op=generic

Luckily the dmextract output contains the region information; it just requires a little special syntax to use it. In this case users can use the CIAO stack "igrid()" syntax to access each individual row in the dmextract output file

$ mkrprm img.fits "rad.prof[component=igrid(1:5:1)]" rad.matrix

Since each row in the dmextract output file has a unique value in the COMPONENT column, we can use that as a filter. We can see the result of using "igrid" using the stk_build command.

$ stk_build "rad.prof[component=igrid(1:5:1)]" out=stdout
rad.prof[component=1]
rad.prof[component=2]
rad.prof[component=3]
rad.prof[component=4]
rad.prof[component=5]

Users could also choose to filter on the #row virtual column.

ds9 format regions: annulus, panda, epanda, bpanda

Certain ds9 regions support multiple radii and angles. These include the annulus, panda, epanda, and bpanda shapes. CIAO does not recongize these shapes so they must be converted into the equivalent CIAO syntax. This can be done using the convert_ds9_region_to_ciao_stack script:

$ convert_ds9_region_to_ciao_stack ds9.reg ciao.lis

Note that the region must be saved in physical coordinates.

Parameter=outfile (file required filetype=output)

The output file name.

The output is the aperture response matrix. It provides the fraction of the region convolved with the PSF from one aperture (eg annulus) is present in all other bin.

The output is stored as a 2D image. The X-axis is the input radius bin number. The Y-axis is the output radius bin number. Bin are counted 1 to N based on the number of input regions.

Parameter=psfmethod (string default=map)

The method used to convolve with the PSF.

Details of the PSF map method are described in the Algorithms section below.

Currently users can only select psfmethod=map.

There is a method that uses MARX to simulate the PSF convolution with each region; however there is a bug in MARX when using images as input that currently makes this choice unusable.

Parameter=ecf (real default=0.393 min=0 max=1)

The encircled counts fraction (ECF).

The enclosed counts fraction of the PSF to use when smoothing with the psfmethod=map method. The default value equal to 0.393 represents a 1 sigma value for a 2D Gaussian.

Parameter=function (string default=gaus)

The shape of the approximate PSF convolution function

The Chandra PSF cannot be expressed analytically. These various shape functions can be used to approximate the morphology of the PSF.

List of available PSF convolution functions

gaus : Gaussian
tophat : A flat, circular shape
lor : A Lorentzian
hemisphere : half of a circular ball
cone : A circular base with linear sides

The Chandra PSF becomes less circular farther away from the optical axis and the aperture correction approximation becomes less accurate.

Parameter=energy (real units=keV)

Energy of the PSF to lookup in keV.

The energy value is used with all psfmethod options.

Parameter=flux (real units=photon/cm**2/sec)

The flux for the marx method. Not currently used.

Parameter=random_seed (integer default=-1)

Initial random seed for simulations.

The random_seed is set to initialize the simulation. With random_seed equal to the default value of -1, the script will use a random value for the initial random_seed.

Parameter=marx_root (file)

The directory where MARX was installed.

The marx executable must be $marx_root/bin/marx.

Parameter=tmpdir (file default=${ASCDS_WORK_PATH})

Directory for temporary files

Parameter=clobber (boolean default=no)

Overwrite output files if they already exist?

Parameter=verbose (integer default=1 min=0 max=5)

Amount of tool chatter level.

Algorithm: psfmethod

psfmethod=map

The psfmethod=map method follows the PSF Size Image Smoothing thread.

First the PSF map is created using the mkpsfmap tool at the requested input energy and ecf values. The units are set to "logical" pixel size which is required for the dmimgadapt step.

Next an image is created for each input region; it is the same size as the infile image. Pixels inside the region are set to 1 divided by the region area (specifically number of pixels). Pixels outside the region are set to 0.

Then dmimgadapt is used to adaptively smooth each region map with PSF map using the specified functional shape.

Since the input image was normalized to sum to 1.0; the output from dmimgadapt gives the fraction of the region with the adaptively smoothed convolution function. We then just iterated through all the region to determine the fraction of the PSF imaged in each region.

psfmethod=marx

This option is not currently available.

Including the Response Matrix when fitting with Sherpa

The traditional spectral RMF is specific to fitting energy distributions and cannot be used for fitting radial profiles.

To use the radial response matrix requires a generic matrix multiplication model. The MatrixModel is available in the sherpa_contrib package. Below is an example of fitting a 1D Gaussian to a radial profile that includes convolution with the radial profile response matrix:

import sherpa.astro.ui as ui
from sherpa_contrib.matrix_model import MatrixModel
from pycrates import read_file

ui.load_data(1,"radial_profile.fits",3,["RMID","SUR_BRI","SUR_BRI_ERR"])

response_matrix = read_file("rp.matrix").get_image().values*1.0
x_vals = ui.get_data().x
mymatrix = MatrixModel(response_matrix, x_vals, name="rprm")

gg2 = ui.gauss1d("gg2")
ui.set_source(mymatrix(gg2))
ui.fit()

The matrix values are loaded using pycrates' "read_file" routine and along with the full set of X-coordinates (x_vals) are used to instantiate the MatrixModel: mymatrix. The matrix is then applied to Gaussian model, gg2, in the "set_source(mymatrix(gg))" command and fit(). The MatrixModel is the type of convolution model akin to the traditional spectral RMF model.

Changes in the scripts 4.17.0 (Decemeber 2024) release

This script is new for CIAO 4.17.

About Contributed Software

This script is not an official part of the CIAO release but is made available as "contributed" software via the CIAO scripts page . Please see this page for installation instructions.

Bugs

There are no known bugs for this tool.