Using an Exposure Map in Fitting Image Data
Sherpa Threads (CIAO 4.10 Sherpa v1)
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: 1 Jun 2018 - reviewed for CIAO 4.10, no content change
- Getting Started
- Reading and Plotting 2-D FITS Data
- Setting the Exposure Map
- Defining and Fitting the Source
- Saving a Sherpa Session
- Scripting It
Please follow the Sherpa Getting Started thread.
Reading and Plotting 2-D FITS Data
We are using 2-D spatial data from the FITS data file 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 x1 = Float64 y = Float64 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
sherpa> contour_data() sherpa> print_window("contour_plot")
This creates Figure 1.
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_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 -6.28319 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 -6.28319 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 2294 Data points = 6400 Degrees of freedom = 6393 Probability [Q-value] = 1 Reduced statistic = 0.509268 Change in statistic = 4.8777e+06 b1.r0 12.4624 b1.xpos 39.5139 b1.ypos 40.8959 b1.ellip 0.0259202 b1.theta -4.6965 b1.ampl 1.31312 b1.alpha 1.66641
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()
Figure 2: Data, model, and fit images
Figure 3: Fit residuals image
Saving a Sherpa Session
To save the Sherpa session in order to return to the analysis at a later point:
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> exec(open("expmap.ascii").read())
The file fit.py is a Python script which performs the primary commands used above; it can be executed by typing exec(open("fit.py").read()) 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, et cetera.)
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.10 syntax to you.
This thread is complete, so we can exit the Sherpa session:
|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)|
|05 Dec 2013||reviewed for CIAO 4.6: no changes|
|18 Mar 2015||reviewed for CIAO 4.7, updated model parameter boundaries.|
|10 Dec 2015||reviewed for CIAO 4.8, updated screen output|
|01 Nov 2016||reviewed for CIAO 4.9, no content change|
|01 Jun 2018||reviewed for CIAO 4.10, no content change|