Using an Exposure Map in Fitting Image Data

Sherpa Threads (CIAO 4.1)

[Python Syntax]

Overview

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

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.

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

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: 
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.

[Print media version: Contour plot of image data]

Figure 1: Contour plot of the data

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_model, then the initial parameter values are specified:

sherpa> set_model(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.

sherpa> fit()
Dataset               = 1
Method                = levmar
Statistic             = chi2gehrels
Initial fit statistic = 4.88095e+06
Final fit statistic   = 3255.75 at function evaluation 65
Data points           = 6400
Degrees of freedom    = 6393
Probability [Q-value] = 1
Reduced statistic     = 0.509268
Change in statistic   = 4.8777e+06
   b1.r0          12.4625     
   b1.xpos        39.5139     
   b1.ypos        40.8959     
   b1.ellip       0.02592     
   b1.theta       4.72827     
   b1.ampl        1.31312     
   b1.alpha       1.66642

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.

[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.

[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:

sherpa> save("expmap.save")

where expmap.save is the output binary file. The restore() command will restore the session when desired.

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.

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