Radial and elliptical profiles of Image Data
Sherpa Threads (CIAO 4.8 Sherpa v1)
Overview
Synopsis:
This thread shows how you can use the sherpa_contrib package—part of the CIAO contributed scripts—to plot radial or elliptical profiles of your imaging fits. These can be used to visualize the fit results, and help you interpret the results—e.g. to see how well the model represents the data or to let you fine tune model parameters before a fit to try and reduce the fit time. The thread is based on the Fitting Imaging Data thread, which should be read before this thread.
Last Update: 10 Dec 2015  reviewed for CIAO 4.8, no content change.
Contents
 Create and Read Fits Image Data
 Initial fit to the data: Gaussian plus constant
 Adding a Gaussian to model the core emission
 Using elliptical models
 Overlays and data access
 Conversion from CIAO 3.4
 Scripting It
 Summary
 History

Images
 Figure 1: Initial radial profile of the data
 Figure 2: Initial radial profile of the model
 Figure 3: A Gaussian plus constant fit to the data
 Figure 4: Adding a second Gaussian component
 Figure 5: Residuals (in counts) about the fit
 Figure 6: Residuals (normalized by the error) about the fit
 Figure 7: Elliptical profile
 Figure 8: Overriding the ellipticity
 Figure 9: The core emission
Create and Read Fits Image Data
Download the sample data: 1838 (ACISS, G21.50.9)
unix% download_chandra_obsid 1838
This repeats the data creation step of the imaging thread, which should be read for an explanation of the steps.
unix% punlearn dmcopy unix% dmcopy \ "acisf01838_repro_evt2.fits[sky=box(4059.25,4235.75,521.5,431.5)][bin x=::2,y=::2]" \ image2.fits
After starting Sherpa, we load in both the radialprofile code and the utility routines, then the data (the from .. import line only has to be called once per session):
sherpa> from sherpa_contrib.all import * sherpa> load_image("image2.fits")
An error message of
ImportError: No module named sherpa_contrib
means that the CIAO contributed scripts package has not been installed or is out of date.
We then set up the coordinate system we use for the fit, along with the statistic and optimization method:
sherpa> set_coord("physical") sherpa> set_stat("cash") sherpa> set_method("simplex")
We can review the current choices for these settings using the following commands or with the show_method and show_stat commands:
sherpa> get_coord() 'physical' sherpa> get_stat_name() 'cash' sherpa> get_method_name() 'neldermead'
We finish this section by restricting the region to be fitted:
Initial fit to the data: Gaussian plus constant
We start with a single Gaussian component to represent the bulk of the emission, and restrict the range of the parameters to values appropriate to this dataset. We use values somewhat more restrictive than used in the original thread, the component name src rather than g1, and take advantage of the guess() routine to set sensible limits on the parameter values:
sherpa> set_source(gauss2d.src) sherpa> print(src) gauss2d.src Param Type Value Min Max Units       src.fwhm thawed 10 1.17549e38 3.40282e+38 src.xpos thawed 0 3.40282e+38 3.40282e+38 src.ypos thawed 0 3.40282e+38 3.40282e+38 src.ellip frozen 0 0 0.999 src.theta frozen 0 6.28319 6.28319 radians src.ampl thawed 1 3.40282e+38 3.40282e+38 sherpa> guess(src) sherpa> print(src) gauss2d.src Param Type Value Min Max Units       src.fwhm thawed 10 1.17549e38 3.40282e+38 src.xpos thawed 4073.5 3965.5 4179.5 src.ypos thawed 4246.5 4142.5 4356.5 src.ellip frozen 0 0 0.999 src.theta frozen 0 6.28319 6.28319 radians src.ampl thawed 280 0.28 280000 sherpa> set_par(src.fwhm, min=0.1, max=300, val=20) sherpa> set_par(src.ampl, min=0.1, max=1000, val=20)
The last two sets of lines could also have been written:
sherpa> src.fwhm = 20 sherpa> src.fwhm.min = 0.1 sherpa> src.fwhm.max = 300 sherpa> src.ampl = 20 sherpa> src.ampl.min = 0.1 sherpa> src.ampl.max = 1000
The current source model definition may be displayed:
sherpa> show_model() Model: 1 gauss2d.src Param Type Value Min Max Units       src.fwhm thawed 20 0.1 300 src.xpos thawed 4073.5 3965.5 4179.5 src.ypos thawed 4246.5 4142.5 4356.5 src.ellip frozen 0 0 0.999 src.theta frozen 0 6.28319 6.28319 radians src.ampl thawed 20 0.1 1000
We now create our first radial plot of the data using the prof_data() command:
sherpa> prof_data() sherpa> log_scale() sherpa> limits(X_AXIS, 0.5, AUTO)
The other commands are from ChIPS and change the plot to use logarithmicscaling for both axes and to change the lowerlimit of the Xaxis to be 0.5. The resulting plot is shown in Figure 1:
[Version: PNG, postscript, encapsulated postscript, PDF]
Figure 1: Initial radial profile of the data
We decide to include a background component before fitting, and so add in a constant model to the source expression. The value of the background model was estimated from Figure 1; this suggests a level of 0.05 counts per physical pixel. However, since the data has been binned by 2 in each direction (the cdelt field is 2), then we need to use a value of 2 * 2 * 0.05, as shown below:
sherpa> print(get_data().sky) physical crval = [ 3798.5, 4019.5] crpix = [ 0.5, 0.5] cdelt = [ 2., 2.] sherpa> set_source(src + const2d.bgnd) sherpa> bgnd.c0 = 0.2 sherpa> bgnd.c0.max = 100
The prof_fit() command will produce a radial profile of the fit overlain on the data (similar to how plot_fit() works). Before making this call we change the plot preferences so that both axes are drawn with logscaling. We use the optional label argument of prof_fit to turn off display of the labels that show the profile center. The result is Figure 2.
sherpa> get_data_prof_prefs()["xlog"] = True sherpa> get_data_prof_prefs()["ylog"] = True sherpa> prof_fit(label=False) sherpa> limits(X_AXIS, 0.5, AUTO)
[Version: PNG, postscript, encapsulated postscript, PDF]
Figure 2: Initial radial profile of the model
It is obvious, from Figure 2, that the width of the Gaussian is too small. A value roughly twice the current size would be a better fit, so we change the src.fwhm value before our first fit.
sherpa> src.fwhm = src.fwhm.val * 2 sherpa> fit() Dataset = 1 Method = neldermead Statistic = cash Initial fit statistic = 31194.7 Final fit statistic = 49847.1 at function evaluation 499 Data points = 9171 Degrees of freedom = 9166 Change in statistic = 18652.4 src.fwhm 57.7556 src.xpos 4072.99 src.ypos 4247.45 src.ampl 23.6745 bgnd.c0 0.272061
Note that we had to say "src.fwhm.val * 2" rather than "src.fwhm * 2" in the assignment operator. If we did not include the trailing ".val" to get at the numeric value of the parameter, Sherpa would have tried to create a parameter link to itself, which would have failed with the following error message:
ParameterError: requested parameter link creates a cyclic reference
We now use the prof_fit_resid() and image_fit() commands to see how well the model fits the data. The output of the prof_fit_resid() command is shown in Figure 3.
sherpa> image_fit() sherpa> prof_fit_resid() sherpa> print_window("fit")
The print_window() call creates a postscript file called fit.ps; see the help file for this command to find out what other formats are available.
[Version: PNG, postscript, encapsulated postscript, PDF]
Figure 3: A Gaussian plus constant fit to the data
The yaxis of the top plot in Figure 3 is drawn with a logarithmicscale since we earlier set the data plot preference for ylog to True. The xaxis has remained with a linearscale because the preference for the residual plot has overridden the data setting. We change the xaxis preferences for both styles of residual plot (resid and delchi) so that future such plots (e.g. Figure 5) will have their xaxis drawn with a logscale:
sherpa> get_resid_prof_prefs()["xlog"] = True sherpa> get_delchi_prof_prefs()["xlog"] = True
Adding a Gaussian to model the core emission
We now add a second Gaussian component to represent the core emission.
sherpa> set_source(bgnd + src + gauss2d.core) sherpa> guess(core) sherpa> set_par(core.fwhm, 10, 0.1, 100) sherpa> set_par(core.ampl, 100, 0.1, 1000) sherpa> prof_fit(model=src) sherpa> limits(X_AXIS, 0.5, AUTO)
The model profile now looks like Figure 4. Since the model expression now contains two components with xpos and ypos parameters we have to say which one should be used to define the profile center; in this case we chose the src component by adding model=src to the argumemt list of prof_fit(). If no model had been specified the call would have failed with the error message:
ValueError: Multiple xpos parameters in source expression: ((const2d.bgnd + gauss2d.src) + gauss2d.core)
[Version: PNG, postscript, encapsulated postscript, PDF]
Figure 4: Adding a second Gaussian component
We now fit the whole model. In the original thread only the core component, g2, is fit at this stage, but the results are similar. We then plot the fit results togther with the residuals in two windows: Figure 5, where the residuals are measured as (datamodel) and Figure 6, where they are measured as (datamodel)/error.
sherpa> fit() Dataset = 1 Method = neldermead Statistic = cash Initial fit statistic = 54106.3 Final fit statistic = 55816.9 at function evaluation 1264 Data points = 9171 Degrees of freedom = 9162 Change in statistic = 1710.61 bgnd.c0 0.200077 src.fwhm 64.144 src.xpos 4073.2 src.ypos 4247.87 src.ampl 17.3732 core.fwhm 6.57941 core.xpos 4073.38 core.ypos 4245.66 core.ampl 227.385 sherpa> prof_fit_resid(model=src, group_counts=200) sherpa> limits(X_AXIS, 0.5, AUTO) sherpa> add_window() sherpa> prof_fit_delchi(model=src, group_counts=200) sherpa> limits(X_AXIS, 0.5, AUTO)
[Version: PNG, postscript, encapsulated postscript, PDF]
Figure 5: Residuals (in counts) about the fit
[Version: PNG, postscript, encapsulated postscript, PDF]
Figure 6: Residuals (normalized by the error) about the fit
Using elliptical models
Up until now the models have been circularlysymmetric (i.e. the ellipticity of the components has been zero). We now see what happens if we let the core emission be elliptical (the clear() call is just to remove the two ChIPS windows created in the previous section):
sherpa> clear() sherpa> thaw(core) sherpa> core.ellip = 0.1 sherpa> fit() Dataset = 1 Method = neldermead Statistic = cash Initial fit statistic = 55865.9 Final fit statistic = 56000.3 at function evaluation 1173 Data points = 9171 Degrees of freedom = 9160 Change in statistic = 134.453 bgnd.c0 0.198276 src.fwhm 64.3286 src.xpos 4073.22 src.ypos 4247.89 src.ampl 17.2214 core.fwhm 8.15547 core.xpos 4073.33 core.ypos 4245.63 core.ellip 0.331718 core.theta 0.780865 core.ampl 227.648 sherpa> prof_fit(model=core, group_counts=200) sherpa> limits(X_AXIS, 0.5, AUTO)
Since the range of the theta parameter of the models now defaults to 2π to 2π, rather than having a lower limit of 0, the theta parameter often does not need to be changed before starting a fit.
In earlier versions of CIAO (<4.6) this thread suggested setting
core.theta = 1
before the fit; changing this value before a fit can often be a sensible choice, as with any other parameter, to avoid getting trapped in local minimums.
The resulting plot is shown in Figure 7.
[Version: PNG, postscript, encapsulated postscript, PDF]
Figure 7: Elliptical profile
You can override any of the component values used to define the profile: in the following we use the position of the core component but ignore the model's ellipticity, and use circular annuli instead. The result is shown in Figure 8.
sherpa> prof_fit(model=core, group_counts=200, ellip=0) sherpa> limits(X_AXIS, 0.5, AUTO)
Overlays and data access
In this section we shall highlight some of the other features of the profiles package:
 plot preferences
 plot overlays
 data access
Plot preferences
The individual plots—such as data, model, resid, and delchi—have preference settings. The current settings can be retrieved with the get_<type>_prof_prefs() command: for example
sherpa> print(get_data_prof_prefs()) {'linethickness': None, 'symbolcolor': None, 'symbolfill': False, 'xlog': True, 'ylog': True, 'symbolangle': None, 'errthickness': None, 'fillcolor': None, 'linecolor': None, 'errstyle': 'line', 'linestyle': 0, 'symbolstyle': 4, 'errcolor': None, 'symbolsize': 3, 'fillstyle': None, 'fillopacity': None, 'yerrorbars': True}
These values can be changed, and new ones added, either by saying
sherpa> get_data_prof_prefs()["symbolcolor"] = "green"
or by assigning a variable to the return value of the call and then changing the variable:
sherpa> p = get_data_prof_prefs() sherpa> p["symbolcolor"] = "green" sherpa> p["symbolfill"] = True sherpa> p["errcolor"] = "green"
A full list of the preferences can be found in the following ahelp pages:
 get_data_prof_prefs()
 get_model_prof_prefs()
 get_source_prof_prefs()
 get_resid_prof_prefs()
 get_delchi_prof_prefs()
Overlaying data
The overplot argument can be used to overlay data on an existing plot. The commands
sherpa> prof_data() sherpa> prof_model(overplot=True)
produce a plot similar to plot_fit(). This capability extends to plot_fit(), so you could say
sherpa> prof_fit() sherpa> prof_fit(2, overplot=True)
to overplot the data and model from dataset 2 onto the plot from the default dataset. In this case it can be hard to distinguish the data and model for the two cases so it might be better to say
sherpa> prof_fit() sherpa> prof_data(2, overplot=True) sherpa> set_histogram(["symbol.color", "blue", "err.color", "blue"]) sherpa> prof_model(2, overplot=True) sherpa> set_histogram(["line.color", "blue"])
which plots the data and model for the second dataset in blue.
Accessing the plot data
For the basic plot types—e.g. data, model, source, resid, delchi, and fit—there are corresponding get_<type>_prof() calls which return objects which contain the data used to create the plots. These routines can be called without having to call the corresponding prof_<type>() routine first.
As an example, consider adding the radial profile of just the core component to Figure 8. One way to accomplish this would be to say
sherpa> clear() sherpa> prof_fit(model=core, group_counts=200, ellip=0) sherpa> limits(X_AXIS, 0.5, AUTO) sherpa> dall = get_data_prof(model=core, group_counts=200, ellip=0) sherpa> set_source(core) sherpa> dcore = get_model_prof(rlo=dall.xlo, rhi=dall.xhi, ellip=0) sherpa> set_source(bgnd + src + core) sherpa> print (dcore) xlo = [ 0., 2., 4., ..., 88., 96., 102.] xhi = [ 2., 4., 6., ..., 96., 102., 110.] y = [ 4.7068e+001, 2.6747e+001, 1.1766e+001, ... 8.3337e163, 2.5899e184] xlabel = Radius (physical pixel) ylabel = Counts per physical pixel title = Model of image2.fits histo_prefs = {'linethickness': 2, 'symbolcolor': None, 'symbolfill': None, 'xlog': False, 'ylog': False, 'symbolangle': None, 'errthickness': None, 'fillcolor': None, 'linecolor': 'red', 'errstyle': None, 'linestyle': 1, 'symbolstyle': 0, 'errcolor': None, 'fillstyle': None, 'fillopacity': None, 'yerrorbars': False, 'symbolsize': None} sherpa> add_histogram(dcore.xlo, dcore.xhi, dcore.y, ["line.color", "orange"]) sherpa> limits(Y_AXIS, 0.03, AUTO) sherpa> add_label(10, 1, r"\color{orange}Core", ["size", 16])
which creates Figure 9. The get_data_prof() call was made to get the bin edges used for the profile (the dall.xlo and dall.xhi arrays), so that these could be given directly to get_model_prof(). In this particular case this step was not necessary, since the same bins would have been calculated if we had just said,
dcore = get_model_prof(model=core, group_counts=200, ellip=0);
but it would have been required if we were trying to match the radial binning calculated from a different dataset.
[Version: PNG, postscript, encapsulated postscript, PDF]
Figure 9: The core emission
Alternatively, Figure 9 could have been created by saying:
sherpa> clear() sherpa> prof_fit(model=core, group_counts=200, ellip=0) sherpa> set_source(core) sherpa> prof_model(group_counts=200, ellip=0, overplot=True) sherpa> set_source(bgnd + src + core) sherpa> limits(X_AXIS, 0.5, AUTO) sherpa> limits(Y_AXIS, 0.03, AUTO) sherpa> set_histogram(["line.color", "orange"]) sherpa> add_label(10, 1, r"\color{orange}Core", ["size", 16])
Conversion from CIAO 3.4
The routines described in this thread are a replacement to the CIAO 3.4 plot_rprof() and related routines from the sherpa_plotfns.sl package. The changes are:

The routines can be run from the Python version of Sherpa.

The routines have been renamed to better match the existing Sherpa plot functionality:
CIAO 3.4 name CIAO 4.6 name plot_rprof prof_fit_resid plot_rprofr prof_fit_resid plot_rprofd prof_fit_delchi plot_eprof prof_fit_resid plot_eprofr prof_fit_resid plot_eprofd prof_fit_delchi There are no longer separate routines for circular and elliptical profiles. Profiles can now be produced just for the data, model, or residuals, using the prof_xxx() series of commands. In CIAO 3.4 it was only possible to create profiles of the fit and residuals.

The routines are simpler to use in that you can often get away with giving no arguments and still get a sensible plot. However, there are more ways to change the plot—e.g. what values are used to define the profile (center and ellipticity), or what bin ranges to use—than there are in the CIAO 3.4 versions.

It is now possible to change the plot preferences—using the get_xxx_prof_prefs() calls—and access the profile values—using the get_xxx_prof() calls.
Scripting It
The commands used to run this thread may be saved in a text file such as fit.py, which can then be executed as a script: execfile("fit.py"). Alternatively, the Sherpa session can be saved to a binary file with the save command (restored with the restore command), or to an editable ASCII file with save_all.
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.8 syntax to you.
Summary
This thread showed you how to use the plotting routines from the sherpa_contrib.profiles package. These routines allow you to view a radial plot of your twodimensional data, model, fit, or residuals; this allows you to inspect the fit parameters so that you start closer to the bestfit solution or to visually inspect the fit results. Since they radially average the results you should always also look at the image residuals too, using the image_fit command.
Both sets of routines can be loaded in one go using the sherpa_contrib module.
The ahelp pages for the individual commands provide more information:
 prof_data()
 prof_model()
 prof_source()
 prof_resid()
 prof_delchi()
 prof_fit()
 prof_fit_resid()
 prof_fit_delchi()
 get_data_prof_prefs()
 get_model_prof_prefs()
 get_source_prof_prefs()
 get_resid_prof_prefs()
 get_delchi_prof_prefs()
 get_data_prof()
 get_model_prof()
 get_source_prof()
 get_resid_prof()
 get_delchi_prof()
 get_fit_prof()
History
15 Apr 2009  new version for CIAO 4.1; based on the imagefitting thread 
24 Jul 2009  replaced use of set_param_limits_from_image() by guess() 
06 Jan 2010  updated for CIAO 4.2 
13 Jul 2010  updated for CIAO 4.2 Sherpa v2: removal of SLang version of thread. 
30 Jan 2012  reviewed for CIAO 4.4. 
05 Dec 2013  reviewed for CIAO 4.6: noted that the minimum value of the theta parameter is now 2π rather than 0. 
11 Dec 2014  reviewed for CIAO 4.7: updated input file and fit results 
10 Dec 2015  reviewed for CIAO 4.8, no content change. 