Accounting for PSF Effects in 2D Image Fitting

Sherpa Threads (CIAO 4.1)

[Python Syntax]

Overview

Last Update: 29 Apr 2009 - new script command is available with CIAO 4.1.2

Synopsis:

In this thread, we fit 2-D image data using a Point-Spread-Function (PSF) image as the convolution kernel.

Proceed to the HTML or hardcopy (PDF: A4 | letter) version of the thread.

Getting Started
Reading & Filtering Image Data
Defining the Instrument Model with a PSF Image
- Instrument model parameters
Defining a Multi-Component Source Model Expression
Modifying Statistic Setting
Fitting
Examining Fit Results
Scripting It
History
Images

Getting Started

Please follow the "Obtaining data used in Sherpa threads" thread.

Reading & Filtering Image Data

In this thread, we fit 2-D image data from the FITS datafile center_box_0.25pix.fits. This is the image created from an event file by binning over a region to 0.25 ACIS pixel size with dmcopy; for example:

unix% dmcopy "event.fits[sky=region(center_box.reg)][bin x=::0.25,y=::0.25]" \
      center_box_0.25pix.fits

This dataset is input into Sherpa with the load_data command:

sherpa> load_data("center_box_0.25pix.fits")

Now the dataset may be displayed:

sherpa> image_data()

The input data image looks like this .

Image data may be filtered in Sherpa interactively with DS9 and if scripted, if region parameters are already known.

DS9 regions may be used to filter the data. After the data have been displayed, go to the Region box in DS9 and choose a designed region shape. In this example we use the box shape. When the box is displayed, the size of the box can be changed with the cursor or within the Region Info box which is displayed from the Region box with the Get Info... button. After the desired region size is set, the region can be used to filter the data as follows:

sherpa> image_getregion()
'rotbox(88.16875,75.8625,70.416667,68.508334,0);'
sherpa> notice2d("rotbox(88.16875,75.8625,70.416667,68.508334,0)")

In this case, the dataset is using physical/image coordinates, so we have to make sure that DS9 is set to interact with Sherpa in physical coordinates by setting Region-->File-->File Coordinate System-->Physical in DS9.

sherpa> show_all() 
Data Set: 1
Filter: Box(88.1688,75.8625,70.4167,68.5083)
name      = center_box_0.25pix.fits
x0        = Float64[24843]
x1        = Float64[24843]
y         = Float64[24843]
shape     = (147, 169)
staterror = None
syserror  = None
sky       = physical
 crval    = [ 4075.    4039.25]
 crpix    = [ 0.5  0.5]
 cdelt    = [ 0.25  0.25]
eqpos     = world
 crval    = [ 248.6211   70.531 ]
 crpix    = [ 4096.5  4096.5]
 cdelt    = [-0.0001  0.0001]
 crota    = 0
 epoch    = 2000
 equinox  = 2000
coord     = logical

The same filter can be set on the command-line with the region definition if it is already known. Note that Sherpa requires the Image coordinates (not the Physical coordinates) to be used in the region definitions:

sherpa> notice2d("box(88.16875,75.8625,70.416667,68.508334)")

Defining the Instrument Model with a PSF Image

An instrument model, named psf0 may be established using an image of a PSF:

sherpa> load_psf("psf0","psf_f1_norm_0.25pix.fits")

Note that the PSF is automatically renormalized to 1. Renormalization is done by summing over all image pixels, regardless of the setting of xsize and ysize. To learn more about 2D PSF instrument models, please read the ahelp file by typing:

sherpa> ahelp "set_psf"

The PSF image file was created using the CIAO tool mkpsf (see the CIAO Create a PSF thread. Note that in recent CIAO versions the best PSF image can be created with The Chandra Ray Tracer (ChaRT)). The file center_box_0.25pix.fits was used as input to mkpsf in order to match the binning of the resulting PSF image file (psf_f1_norm_0.25pix.fits). Sherpa requires that the binning of the PSF image file match the binning of the input data image.

To view the PSF image , load it into DS9 (from outside Sherpa).

Instrument model parameters

size:

The PSF image provided may be much larger than the instrument PSF size (in number of pixels), and therefore larger than needed. In order to speed the convolution process, a sub-image kernel is specified. In this example, the input PSF image file has 256x256 pixels, but the portion that will be used for the convolution is only the center 32x32 pixels sub-image.

sherpa> set_psf(psf0)
sherpa> psf0.size = [32,32]
sherpa> psf0.center = [128,129]
sherpa> print get_psf()
psfmodel.psf0
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   psf0.kernel  frozen psf_f1_norm_0.25pix.fits
   psf0.size    frozen     (32, 32)     (32, 32)     (32, 32)
   psf0.center  frozen   (128, 129)   (128, 129)   (128, 129)
   psf0.radial  frozen            0            0            1           
   psf0.norm    frozen            1            0            1           
sherpa> image_psf()
PSF frac: 0.976006417088

This sub-image contains only a part of the input file defined by the size and center parameters. Note that the sub-image will be empty if the PSF centroid is located outside the sub-image. When the image_psf command is issued by Sherpa, information about the PSF fraction included in the sub-image is printed.

In some cases, a larger sub-image will be needed: when the source is located off-axis, or when including the PSF wings is important for the analysis. We may update the parameters to increase the PSF fraction to about 99%

sherpa> set_psf(psf0)
sherpa> psf0.size = [72,72]
sherpa> print get_psf()
psfmodel.psf0
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   psf0.kernel  frozen psf_f1_norm_0.25pix.fits
   psf0.size    frozen     (72, 72)     (72, 72)     (72, 72)
   psf0.center  frozen   (128, 129)   (128, 129)   (128, 129)
   psf0.radial  frozen            0            0            1           
   psf0.norm    frozen            1            0            1           
sherpa> image_psf()
PSF frac: 0.991270796001

The instrument PSF sub-images look like this . The left panel shows the smaller sub-image of 32x32 pixel. The right panel shows an expanded to 72x72 pixels sub-image. Notice that the larger sub-image contains a significant amount of the PSF wings structure.

center:
The instrument PSF (kernel) centroid must always be at the center of the extracted sub-image. The parameter center=[center0,center1] determines the location of the original PSF image from which the sub-image is centered on.

Defining a Multi-Component Source Model Expression

Now we will define a source model expression, using a two-dimensional Gaussian function called gauss2d, and a constant called const2d which will be convolved by the PSF model using FFT. We will also guess reasonable parameter values to begin the fitting.

sherpa> set_source(const2d.c1 + gauss2d.g1)
sherpa> c1.c0 = 1.0
sherpa> g1.fwhm = 6
sherpa> g1.xpos = 90
sherpa> g1.ypos = 80
sherpa> g1.ampl = 100

sherpa> show_source()
Model: 1
(const2d.c1 + gauss2d.g1)
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   c1.c0        frozen            1            0  3.40282e+38
   g1.fwhm      thawed            6  1.17549e-38  3.40282e+38
   g1.xpos      thawed           90 -3.40282e+38  3.40282e+38
   g1.ypos      thawed           80 -3.40282e+38  3.40282e+38
   g1.ellip     frozen            0            0        0.999
   g1.theta     frozen            0            0      6.28319    radians
   g1.ampl      thawed          100 -3.40282e+38  3.40282e+38

The model component c1 is interpreted as the constant background contribution to the data. We can also see how the source model is convolved with the PSF, represented as "(kernel)," by using show_model.

sherpa> show_model()
Model: 1
(kernel) o ((const2d.c1 + gauss2d.g1))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   c1.c0        thawed            1            0  3.40282e+38           
   g1.fwhm      thawed            6  1.17549e-38  3.40282e+38           
   g1.xpos      thawed           90 -3.40282e+38  3.40282e+38           
   g1.ypos      thawed           80 -3.40282e+38  3.40282e+38           
   g1.ellip     frozen            0            0        0.999           
   g1.theta     frozen            0            0      6.28319    radians
   g1.ampl      thawed          100 -3.40282e+38  3.40282e+38           
   psf0.kernel  frozen psf_f1_norm_0.25pix.fits
   psf0.size    frozen     (72, 72)     (72, 72)     (72, 72)
   psf0.center  frozen   (128, 129)   (128, 129)   (128, 129)
   psf0.radial  frozen            0            0            1           
   psf0.norm    frozen            1            0            1

Modifying Statistic Setting

Since the data being fit has low counts, we will change the statistics to cash:

sherpa> set_stat("cash")

Note that truncation is turned on when the Cash statistic is used. This setting prevents negative model-predicted data values from affecting the convergence process.

Further details about the Cash statistic method are available by typing:

sherpa> ahelp "cash"

Fitting

The data is first fit assuming a constant background (i.e. c1 is frozen). We also want to change our optimization method, since the default Sherpa method is Levenberg-Marquardt (levmar), which is not robust with our choice of statistics. We update the optimization method to the Nelder-Mead Simplex method.

sherpa> c1.c0 = 1.0
sherpa> freeze(c1.c0)
sherpa> set_method("neldermead")
sherpa> set_method_opt("iquad",0)
sherpa> set_method_opt("finalsimplex",0)
sherpa> fit()
Dataset               = 1
Method                = neldermead
Statistic             = cash
Initial fit statistic = 7648.86
Final fit statistic   = 4066.78 at function evaluation 305
Data points           = 4899
Degrees of freedom    = 4895
Change in statistic   = 3582.08
   g1.fwhm        2.80117     
   g1.xpos        88.661      
   g1.ypos        77.2271     
   g1.ampl        166.649

The fit is run again with the background component thawed.

sherpa> thaw(c1.c0)
sherpa> fit() 
Dataset               = 1
Method                = neldermead
Statistic             = cash
Initial fit statistic = 4066.78
Final fit statistic   = -4440.06 at function evaluation 200
Data points           = 4899
Degrees of freedom    = 4894
Change in statistic   = 8506.84
   c1.c0          0.00282397  
   g1.fwhm        3.1077      
   g1.xpos        88.5928     
   g1.ypos        77.2558     
   g1.ampl        167.44

Now, we will get the confidence level of the fit, and in doing so, recalculate the best-fit parameters using the projection method.

sherpa> projection()
Confidence limits search for c1.c0 finished (1 of 5)
Confidence limits search for g1.fwhm finished (2 of 5)
Confidence limits search for g1.xpos finished (3 of 5)
Confidence limits search for g1.ypos finished (4 of 5)
WARNING: New minimum statistic found while computing confidence limits
WARNING: New best-fit parameters:
Method                = neldermead
Statistic             = cash
Initial fit statistic = -4441.01
Final fit statistic   = -4457.71 at function evaluation 259
Data points           = 4899
Degrees of freedom    = 4894
Change in statistic   = 16.7061
   c1.c0          0.0027721   
   g1.fwhm        3.73769     
   g1.xpos        88.5796     
   g1.ypos        77.241      
   g1.ampl        113.139     
Confidence limits search for c1.c0 finished (1 of 5)
Confidence limits search for g1.fwhm finished (2 of 5)
Confidence limits search for g1.xpos finished (3 of 5)
Confidence limits search for g1.ypos finished (4 of 5)
Confidence limits search for g1.ampl finished (5 of 5)
Dataset               = 1
Confidence Method     = projection
Fitting Method        = neldermead
Statistic             = cash
projection 1-sigma (68.2689%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   c1.c0           0.0027721  -0.00124881   0.00147312
   g1.fwhm           3.73769    -0.151937     0.151556
   g1.xpos           88.5796   -0.0688209    0.0687766
   g1.ypos            77.241   -0.0682646    0.0676121
   g1.ampl           113.139      -9.0049      10.1892

Examining Fit Results

The fit results may be examined with the image_fit command :

sherpa> image_fit()

This command displays the data (upper left), best fit model (upper right), and residuals (lower left) in three DS9 frames, as shown in this image . This montage of images can be used to see how well the model describes the data.

The residuals may also be examined using a contour plot:

sherpa> contour_resid()
sherpa> limits(X_AXIS, 55, 120)
sherpa> limits(Y_AXIS, 45, 110)
sherpa> get_contour_zrange()
            [-12.241783285130133, 12.760612823792094]

The resulting contour plot looks like this , and get_contour_zrange() provides minimum and maximum values of the plot.

The image model and residual data (in counts) may be written to external FITS image files with the following series of commands:

sherpa> image_crate = pack_image()
sherpa> set_imagevals(image_crate,get_model_image().y)
sherpa> write_file(image_crate,"2dpsf_model_image.fits")

sherpa> set_imagevals(image_crate, get_resid_image().y)
sherpa> write_file(image_crate, "2dpsf_resid_cnts.fits")

This residuals file "2dpsf_resid_cnts.fits" may then be used in further analysis, such as smoothing the residuals with aconvolve.

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.

Scripting It

The file fit.py is a Python script which performs the primary commands used above; it can be executed by typing execfile("fit.py") on the Sherpa command line.

The Sherpa script command may be used to save everything typed on the command line in a Sherpa session:

sherpa> script(filename="sherpa.log", clobber=False)

The CXC is committed to helping Sherpa users transition to new syntax as smoothly as possible. If you have existing Sherpa scripts or save files, submit them to us via the CXC Helpdesk and we will provide the CIAO/Sherpa 4.1 syntax to you.

History

14 Dec 2004	reviewed for CIAO 3.2: no changes
21 Dec 2005	reviewed for CIAO 3.3: minor changes to fit results
29 Jun 2006	added `() = evalfile("sherpa_plotfns.sl");` command in Examining Fit Results section
01 Dec 2006	reviewed for CIAO 3.4: no changes
03 Dec 2008	reviewed for CIAO 4.1: updated for Sherpa 4.1
03 Apr 2009	added note on 2D blurring
24 Apr 2009	information on PSF models is in the "set_psf" ahelp file (previously was "ahelp psf2d")
29 Apr 2009	new `script` command is available with CIAO 4.1.2