Using an Exposure Map in Fitting Image Data

Sherpa Threads (CIAO 4.5 Sherpa v1)

Overview

Synopsis:

This thread shows how to use an exposure map when fitting 2-D spatial data. The exposure map file is input to Sherpa as a file-based exposure map model via the load_table_model function.

Last Update: 13 Dec 2012 - reviewed for CIAO 4.5 (no changes)

Getting Started

Please follow the Sherpa Getting Started thread.

Reading and Plotting 2-D FITS Data

We are using 2-D spatial data from the FITS datafile img.fits. This data set is input to Sherpa with the load_image command:

sherpa> load_image("img.fits")

sherpa> show_data()
Data Set: 1
Filter: 
{header information omitted for brevity}
name      = img.fits
x0        = Float64[6400]
x1        = Float64[6400]
y         = Float64[6400]
shape     = (80, 80)
staterror = None
syserror  = None
sky       = physical
 crval    = [ 3944.  3920.]
 crpix    = [ 0.5  0.5]
 cdelt    = [ 5.  5.]
eqpos     = world
 crval    = [ 40.0117  59.9967]
 crpix    = [ 4096.5  4096.5]
 cdelt    = [-0.0001  0.0001]
 crota    = 0
 epoch    = 2000
 equinox  = 2000
coord     = logical

The data set may be viewed as a contour plot (contour_data) or an image (image_data). Here we show the contour plot method, creating a PostScript file of the output as well:

sherpa> contour_data()
sherpa> print_window("contour_plot")

This creates Figure 1.

Figure 1: Contour plot of the data

[Print media version: Contour plot of image data]

Figure 1: Contour plot of the data

Contour plot of image data set created with the contour_data function.

Setting the Exposure Map

We define a file-based exposure map model by loading an exposure map file with the load_table_model command:

sherpa> load_table_model("emap", "expmap.fits")

To display the status of the model emap, use the print() command; notice that Sherpa identifies the exposure map model as "tablemodel.emap":

sherpa> print(emap)
tablemodel.emap
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   emap.ampl    thawed            1 -3.40282e+38  3.40282e+38

Defining and Fitting the Source

One can now define a model to be used as a source model. After viewing Figure 1, the BETA2D model is found to be a promising candidate for the source. The BETA2D model is defined for the source with set_source, then the initial parameter values are specified:

sherpa> set_source(beta2d.b1*emap) 

sherpa> show_model()

Model: 1
(beta2d.b1 * tablemodel.emap)
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   b1.r0        thawed           10  1.17549e-38  3.40282e+38           
   b1.xpos      thawed            0 -3.40282e+38  3.40282e+38           
   b1.ypos      thawed            0 -3.40282e+38  3.40282e+38           
   b1.ellip     frozen            0            0        0.999           
   b1.theta     frozen            0            0      6.28319    radians
   b1.ampl      thawed            1 -3.40282e+38  3.40282e+38           
   b1.alpha     thawed            1          -10           10           
   emap.ampl    thawed            1 -3.40282e+38  3.40282e+38      

sherpa> b1.r0 = 30
sherpa> b1.xpos = 40
sherpa> b1.ypos = 40
sherpa> b1.ellip = 0.3
sherpa> b1.theta = 5
sherpa> b1.ampl = 3.0
sherpa> b1.alpha = 1.5

sherpa> thaw(b1.ellip, b1.theta)
 
sherpa> freeze(emap.ampl)

sherpa> show_model()
Model: 1
(beta2d.b1 * tablemodel.emap)
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   b1.r0        thawed           30  1.17549e-38  3.40282e+38           
   b1.xpos      thawed           40 -3.40282e+38  3.40282e+38           
   b1.ypos      thawed           40 -3.40282e+38  3.40282e+38           
   b1.ellip     thawed          0.3            0        0.999           
   b1.theta     thawed            5            0      6.28319    radians
   b1.ampl      thawed            3 -3.40282e+38  3.40282e+38           
   b1.alpha     thawed          1.5          -10           10           
   emap.ampl    frozen            1 -3.40282e+38  3.40282e+38

Next, we fit the model to the data. By default, Sherpa uses Chi2Gehrels as the fit statistic and LevMar as the fit optimization method. We decide to change the optimization method for this fit to Nelder-Mead:

sherpa> set_method("neldermead")
sherpa> fit()
Dataset               = 1
Method                = neldermead
Statistic             = chi2gehrels
Initial fit statistic = 4.88095e+06
Final fit statistic   = 3255.75 at function evaluation 1595
Data points           = 6400
Degrees of freedom    = 6393
Probability [Q-value] = 1
Reduced statistic     = 0.509268
Change in statistic   = 4.8777e+06
   b1.r0          12.4631     
   b1.xpos        39.5138     
   b1.ypos        40.8959     
   b1.ellip       0.0259077   
   b1.theta       4.72822     
   b1.ampl        1.31309     
   b1.alpha       1.66649

To display the fit and residuals of the plot, we use image_fit and image_resid; the former opens in DS9 a display of the the data image, model image, and fit image, and the latter the (data - model) fit residuals image.

sherpa> image_fit()

sherpa> image_resid()

These commands create Figure 2 and Figure 3.

Figure 2: Data, model, and fit images

[DS9 display of data image, model image, and fit image]

Figure 2: Data, model, and fit images

DS9 display of data image, model image, and fit image.

Figure 3: Fit residuals image

[Print media version: DS9 display of fit residuals image]

Figure 3: Fit residuals image

DS9 display of the (data-model) fit residuals image.

Saving a Sherpa Session

To save the Sherpa session in order to return to the analysis at a later point:

sherpa> save("expmap.save")

sherpa> save_all("expmap.ascii")

where the save function records all the information about the current session to the binary file expmap.save, and the save_all function records the session settings to an editable ASCII file.

To restore the session that was saved to the binary file expmap.save or ASCII file expmap.ascii:

sherpa> restore("expmap.save")

sherpa> execfile("expmap.ascii")

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)

(Note that restoring a Sherpa session from such a file could be problematic since it may include syntax errors, unwanted fitting trials, etcetera.)

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.5 syntax to you.

Summary

This thread is complete, so we can exit the Sherpa session:

sherpa> quit

History

14 Jan 2005	reviewed for CIAO 3.2: no changes
21 Dec 2005	reviewed for CIAO 3.3: no changes
01 Dec 2006	reviewed for CIAO 3.4: no changes
07 Dec 2008	updated for Sherpa 4.1
29 Apr 2009	new `script` command is available with CIAO 4.1.2
07 Jan 2010	updated for CIAO 4.2
13 Jul 2010	updated for CIAO 4.2 Sherpa v2: removal of S-Lang version of thread.
30 Jan 2012	reviewed for CIAO 4.4 (no changes)
13 Dec 2012	reviewed for CIAO 4.5 (no changes)