Plotting Model Components (xfig example) - Revision history

Obst: /* The Result */

2018-04-18T11:07:01Z

The Result

Obst at 11:06, 18 April 2018

2018-04-18T11:06:17Z

Obst: /* Plotting Model Components */

2018-04-18T11:04:53Z

Plotting Model Components

Obst: Obst moved page Plotting Model Components to Plotting Model Components (xfig example) without leaving a redirect

2018-04-18T10:29:33Z

Obst moved page Plotting Model Components to Plotting Model Components (xfig example) without leaving a redirect

Falkner at 16:01, 13 April 2018

2018-04-13T16:01:28Z

Falkner: /* The Result */

2018-04-13T15:58:43Z

The Result

Falkner: Created page with " Category:SLxfig ====== Plotting Model Components ====== Even in ISIS, plotting individual model components is somewhat tricky. Once upon a time, there was a useful des..."

2018-04-13T15:57:38Z

Created page with " Category:SLxfig ====== Plotting Model Components ====== Even in ISIS, plotting individual model components is somewhat tricky. Once upon a time, there was a useful des..."

New page

[[Category:SLxfig]]

====== 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 [[http://www.black-hole.eu/media/summerschool2/X-ray_Spectra_Part_II.pdf|here]] 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.

===== The plot in ISIS =====

=== Loading the data ===
<code>

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, );

</code>

=== Starting the plot ===
<code>
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

</code>

and now a single component of the model:

<code>
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

</code>

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 ===

<code>
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
</code>

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

<code>
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
</code>

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

===== Reading and plotting the data with xfig =====

=== Initializing the plot ===
<code>
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
</code>

=== Read and plot the ASCII files: upper panel ===
<code>
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
</code>

now the model component:

<code>
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
</code>

=== Read and plot the ASCII files: lower panel ===
<code>
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)
</code>

now the residuals without the component:

<code>
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
</code>

=== Render the plot ===
<code>
xfig_multiplot (p).render ("plot_components.png");
</code>

===== The Result =====
{{:isis:slxfig:plot_individual_components.png|}}

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]);

@@ Line 169: / Line 169: @@
 [[File:plot_individual_components.png]]
-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:
+To remove the annoying vertical lines at the first and last bin of the histogram, you can use the (undocumented) qualifiers <code>y_last</code> and <code>y_first</code> of <code>hplot</code>. 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]);

@@ Line 2: / Line 2: @@
 [[Category:SLxfig]]
 = 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 [http://www.black-hole.eu/media/summerschool2/X-ray_Spectra_Part_II.pdf here] 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.
-===== The plot in ISIS =====
+= The plot in ISIS =
-=== Loading the data ===
+== Loading the data ==
 <pre>
@@ Line 20: / Line 19: @@
 </pre>
-=== Starting the plot ===
+== Starting the plot ==
 <pre>
 open_plot;                                                            % open a plot window...
@@ Line 57: / Line 56: @@
 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 ===
+== The second panel ==
 <pre>
@@ Line 86: / Line 85: @@
 Now you have all the files you need to make a nice xfig plot.
-===== Reading and plotting the data with xfig =====
+= Reading and plotting the data with xfig =
-=== Initializing the plot ===
+== Initializing the plot ==
 <pre>
 require("xfig");                                                      % you need this to do the nice plots
@@ Line 106: / Line 105: @@
 </pre>
-=== Read and plot the ASCII files: upper panel ===
+== Read and plot the ASCII files: upper panel ==
 <pre>
 variable lo, hi, y, e, m;                                             % five variables for the five columns of the .dat file
@@ Line 134: / Line 133: @@
 </pre>
-=== Read and plot the ASCII files: lower panel ===
+== Read and plot the ASCII files: lower panel ==
 <pre>
 variable lo, hi, y, el, eu;                                           % five variables for the columns of the .res file
@@ Line 162: / Line 161: @@
 </pre>
-=== Render the plot ===
+== Render the plot ==
 <pre>
 xfig_multiplot (p).render ("plot_components.png");
 </pre>
-===== The Result =====
+= The Result =
 [[File:plot_individual_components.png]]
 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]);

@@ Line 5: / Line 5: @@
 ====== 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 [[http://www.black-hole.eu/media/summerschool2/X-ray_Spectra_Part_II.pdf|here]] 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.
+Even in ISIS, plotting individual model components is somewhat tricky. Once upon a time, there was a useful description to be found [http://www.black-hole.eu/media/summerschool2/X-ray_Spectra_Part_II.pdf here] 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.
 ===== The plot in ISIS =====

← Older revision	Revision as of 10:29, 18 April 2018
(No difference)

@@ Line 11: / Line 11: @@
 === Loading the data ===
-<code>
+<pre>
 load_data("ridge_epoch_4_2.pha");                                     % load your spectra
@@ Line 19: / Line 19: @@
 ignore_en(1,20, );
-</code>
+</pre>
 === Starting the plot ===
-<code>
+<pre>
 open_plot;                                                            % open a plot window...
 multiplot([5,2]);                                                     % ...to get the two panels for spectrum and residuals
@@ Line 40: / Line 40: @@
                                                                        % plot it will end in .dat
-</code>
+</pre>
 and now a single component of the model:
-<code>
+<pre>
 ff = get_fit_fun();                                                   % save the complete model
@@ Line 54: / Line 54: @@
 write_plot("write_data_gauss");                                       % again, dump the plot into another ASCII file
-</code>
+</pre>
 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.
@@ Line 60: / Line 60: @@
 === The second panel ===
-<code>
+<pre>
 mpane(2);                                                             % move to the lower panel
 ylin;                                                                 % for residuals, a linear y scale
@@ Line 72: / Line 72: @@
 write_plot("write_res_allin");                                        % dump the plot into an ASCII file. for residual plots
                                                                        % it will end in .res
-</code>
+</pre>
 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):
-<code>
+<pre>
 set_par([5,8,11],0);                                                  % set the norm of your component to zero
 () = eval_counts;                                                     % evaluate the model
@@ Line 83: / Line 83: @@
 write_plot("write_res_gauss.dat");                                    % dump the plot into an ASCII file
 close_plot;                                                           % close the plot window when finished
-</code>
+</pre>
 Now you have all the files you need to make a nice xfig plot.
@@ Line 90: / Line 90: @@
 === Initializing the plot ===
-<code>
+<pre>
 require("xfig");                                                      % you need this to do the nice plots
@@ Line 105: / Line 105: @@
 p[1].ylabel(`$\chi$`);                                                % ylabel of the lower panel
 p[1].xlabel("Energy [keV]");                                          % xlabel of the whole plot
-</code>
+</pre>
 === Read and plot the ASCII files: upper panel ===
-<code>
+<pre>
 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
@@ Line 123: / Line 123: @@
 p[0].plot (x,y,e; sym="point");                                       % a normal plot to get data points and errorbars
-</code>
+</pre>
 now the model component:
-<code>
+<pre>
 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
@@ Line 133: / Line 133: @@
 p[0].hplot (h,m; color="green4");                                     % plot a histogram of the model component
-</code>
+</pre>
 === Read and plot the ASCII files: lower panel ===
-<code>
+<pre>
 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
@@ Line 147: / Line 147: @@
                                                                        % optional: add the errorbars, if you did the plot with an
                                                                        % individually sized last bin (i.e., added the last bin_hi, see above)
-</code>
+</pre>
 now the residuals without the component:
-<code>
+<pre>
 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
@@ Line 161: / Line 161: @@
 p[1].plot (x,y, ,{y-el,eu-y}; sym="point", color="green4", size=-1, eb_factor=0);
                                                                        % optional: add the errorbars
-</code>
+</pre>
 === Render the plot ===
-<code>
+<pre>
 xfig_multiplot (p).render ("plot_components.png");
-</code>
+</pre>
 ===== The Result =====