Radial and elliptical profiles of Image Data

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.

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:

The plot shows the radial profile of the data, using the current
model values to define the profile center: in this case
src.xpos and src.ypos.
These values are displayed in the top-right of the plot
as x0 and y0.
The profile shows that a single component is unlikely
to be sufficient to adequately describe the profile.

As no arguments were given to
prof_data(),
the profile was generated using the bin width of
the data (so in this case 2 physical pixels),
and extends from r=0 to the maximum
radius of the data (in this case close to 110
sky pixels).

Note that the errors are calculated using the equation
error = 1 + sqrt(N + 0.75)
where N is the number of counts in a bin (and then
normalized by the bin area).

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:

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 log-scaling. 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.

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:

The top plot shows the data (circles) and best-fit model (red line),
whilst the bottom plot shows the residuals—measured as
(data - model)—of this fit.
Note that the residuals are
not
normalized by the area of each annulus.

The y-axis of the top plot in Figure 3
is drawn with a logarithmic-scale since we earlier set the
data plot preference for ylog to
True. The x-axis has remained with a linear-scale because
the preference for the residual plot has over-ridden the
data setting. We change the x-axis preferences for both
styles of residual plot (resid and delchi)
so that future such plots
(e.g. Figure 5)
will have their x-axis drawn with a log-scale:

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:

The second Gaussian component (core) has been added to try and describe the
emission in the core of the object. The initial parameters we guessed are not ideal
but are close enough to begin the fit.

As there are now two model components with a
center—src and core—we
have to decide which one should be used to define the
profile center. In this case we chose the src component by saying model=src in the
call to prof_fit().

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
(data-model)
and Figure 6, where they are measured as
(data-model)/error.

The addition of the core component produces a much better fit to the
source emission. Note however, that this is a radial profile, so
image_resid()
should still be used to look for
differences that may be hidden by the radial averaging.

The group_counts option was used to change the radial
binning; after calculating the bins using the default
values—e.g. the same as used to create the earlier
plots—the data is grouped together until each new radial group contains at
least 200 counts. Since the source emission is high
in the core regions this only affects the profile at large
radii in this example.

Up until now the models have been circularly-symmetric (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):

Since the model used to define the profile—core—has a
non-zero ellipticity then the annuli used to create this profile are
elliptical. The values used for the ellipticity and position angle
are included in the labels in the top-right of the plot.

You can over-ride 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.

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

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.

The orange histogram shows the radial profile of just the core
component. It was created using circular apertures (i.e.
the model's ellip parameter was over-ridden in the call
to get_model_prof().
The add_label() routine
was used to add a label identifying the profile of the
core component.

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.

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

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
two-dimensional data, model, fit, or residuals;
this allows you to inspect the fit parameters so that
you start closer to the best-fit 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.