Plotting Model Components (xfig example)

Plotting Model Components

Even in ISIS, plotting individual model components is somewhat tricky. Once upon a time, there was a useful description to be found [[1]] on pp. 39-43. Unfortunately, since then, the fancy plot commands changed but not the presentation. I will try here to give you a short overview of how to do this at the moment.

Loading the data

load_data("ridge_epoch_4_2.pha");                                     % load your spectra
flux_corr(1);                                                         % the usual setup things
group(1; min_sn=5, min_chan=1, bounds=3, unit="keV");
ignore_en(1, ,3);
ignore_en(1,20, );

Starting the plot

open_plot;                                                            % open a plot window...
multiplot([5,2]);                                                     % ...to get the two panels for spectrum and residuals
mpane(1);                                                             % start plotting in the upper panel
xlog; ylog;                                                           % if that's the usual way you want to have your spectra

load_par("2011-03-09_15:00:21_best.par");                             % load the best-fit parameters for the whole model
() = eval_counts;                                                     % evaluate the model

plot_data(1; dsym=-4, dcol=1, decol=1, mcol=15, xrange={2.9,20.5}, yrange={0.005,1.5}, no_reset=1);
                                                                      % the command to plot your spectrum. note that the
                                                                      % "no_reset" command now has to be used as a qualifier,
                                                                      % different from Mike's presentation. also, the commands
                                                                      % for setting the ranges have changed (from "xrng" to "xrange").

write_plot("write_data_allin");                                       % dump the spectrum and model into an ASCII file. for a data
                                                                      % plot it will end in .dat

and now a single component of the model:

ff = get_fit_fun();                                                   % save the complete model
                                                                      
fit_fun("gauss(1)");                                                  % now only the component you are interested in
                                                                      % **ATTENTION**: next time, reload the best-fit parameters before
                                                                      % adding another model component: "fit_fun(ff)"!!!
() = eval_counts;                                                     % evaluate the model (= the single component)
plot_data(1; dsym=-4, dcol=1, decol=1, mcol=4, oplt=1, no_reset=1);   % plot the component in a different color
write_plot("write_data_gauss");                                       % again, dump the plot into another ASCII file

Of course the colors are unimportant for the moment, as you will have to do them again when doing the plot in xfig, but now the example is complete if you don't have the time for xfig and just want a quick way to plot.

The second panel

mpane(2);                                                             % move to the lower panel
ylin;                                                                 % for residuals, a linear y scale

load_par("2011-03-09_15:00:21_best.par");                             % again, reload the best-fit parameters
() = eval_counts;                                                     % evaluate the model

plot_residuals(1; res=1, rsym=0, rcol=15, recol=1, xrange={2.9,20.5}, yrange={-3,20}, no_reset=1);
                                                                      % plot only the residuals. qualifier are the
                                                                      % same as for "plot_data"
write_plot("write_res_allin");                                        % dump the plot into an ASCII file. for residual plots
                                                                      % it will end in .res

Now the residuals you get when the component is missing in the model (in this case: an iron line complex consisting of three gaussian lines):

set_par([5,8,11],0);                                                  % set the norm of your component to zero
() = eval_counts;                                                     % evaluate the model
plot_residuals(1; res=1, rsym=0, rcol=4, recol=1, oplt=1, no_reset=1);
                                                                      % plot the residuals
write_plot("write_res_gauss.dat");                                    % dump the plot into an ASCII file
close_plot;                                                           % close the plot window when finished

Now you have all the files you need to make a nice xfig plot.

Initializing the plot

require("xfig");                                                      % you need this to do the nice plots

variable p = Struct_Type[2];                                          % this way, you can do a multiplot
p[0] = xfig_plot_new(12.6,6.4);                                       % the size of the upper panel
p[1] = xfig_plot_new(12.6,2.6);                                       % size of the lower panel. proportions are the same as in ISIS plots.

p[0].world (2.9,20.5,0.005,1.5; xlog=1, ylog=1);                      % ranges and scales of the upper panel
p[0].x1axis(;ticlabels=0);                                            % no ticlabels between the two panels
p[0].y2axis(;ticlabels=0);                                            % no ticlabels at the right side of the plot
p[0].ylabel("Counts/s/keV");                                          % ylabel of the upper panel

p[1].world (2.9,20.5,-5,41; xlog=1);                                  % ranges and scales of the lower panel
p[1].ylabel(`$\chi$`);                                                % ylabel of the lower panel
p[1].xlabel("Energy [keV]");                                          % xlabel of the whole plot

Read and plot the ASCII files: upper panel

variable lo, hi, y, e, m;                                             % five variables for the five columns of the .dat file
(lo, hi, y, e, m) = readcol ("write_data_allin_0.dat",1,2,3,4,5);     % read all five columns: bin_lo, bin_hi, value, error, model
variable h = [lo, hi[-1]];                                            % add the upper limit of the last bin to the lower limits of the bins
variable x = (lo+hi)/2;                                               % an array with the central points of each bin

p[0].hplot (h,m; color="blue4");                                      % make a histogram plot of the model
                                                                      % see "p[0].hplot(;help)" for help: 
                                                                      % if length(h)=length(m)+1, the last element is taken as the upper
                                                                      % limit of the last bin. if both arguments have the same length,
                                                                      % the last bin will have the same size as the one before.
                                                                      % **ATTENTION** at the moment, hplot can only plot a histogram with
                                                                      % errorbars when all arrays have the same length!!!
                                                                      
p[0].plot (x,y,e; sym="point");                                       % a normal plot to get data points and errorbars

now the model component:

variable lo, hi, y, e, m;                                             % clear the variables the previous variables
(lo, hi, y, e, m) = readcol ("write_data_gauss_0.dat",1,2,3,4,5);     % read the ASCII file with the component
variable h = [lo, hi[-1]];                                            % add the last bin_hi for histogram plotting

p[0].hplot (h,m; color="green4");                                     % plot a histogram of the model component

Read and plot the ASCII files: lower panel

variable lo, hi, y, el, eu;                                           % five variables for the columns of the .res file
(lo, hi, y, el, eu) = readcol ("write_res_allin_0.res",1,2,3,4,5);    % read all five columns: bin_lo, bin_hi, value, err_lo, err_hi
variable h = [lo, hi[-1]];                                            % add the last bin_hi for histogram plotting
variable x = (lo+hi)/2;                                               % the center of each bin

p[1].hplot (h,y; color="blue4");                                      % plot the residuals histogram of the whole model

p[1].plot (x,y, ,{y-el,eu-y}; sym="point", color="blue4", size=-1, eb_factor=0);
                                                                      % optional: add the errorbars, if you did the plot with an
                                                                      % individually sized last bin (i.e., added the last bin_hi, see above)

now the residuals without the component:

variable lo, hi, y, el, eu;                                           % clear the variables
(lo, hi, y, el, eu) = readcol ("write_res_gauss_0.res",1,2,3,4,5);    % read the .res file
variable h = [lo, hi[-1]];                                            % add the last bin_hi for histogram plotting
variable x = (lo+hi)/2;                                               % center of each bin

p[1].hplot (h,y; color="green4");                                     % plot the residuals histogram without the component\

p[1].plot (x,y, ,{y-el,eu-y}; sym="point", color="green4", size=-1, eb_factor=0);
                                                                      % optional: add the errorbars

Render the plot

xfig_multiplot (p).render ("plot_components.png");

The Result

To remove the annoying vertical lines at the first and last bin of the histogram, you can use the (undocumented) qualifiers y_last and y_first of hplot. If you set them to the value of the first or the last bin respectively, the vertical lines vanish:

 p[0].hplot (h,m; color="blue4", y_first=m[0], y_last=m[-1]);