Mike's fancy plots (xfig example) - Revision history

Obst: Obst moved page Mike's fancy plots to Mike's fancy plots (xfig example) without leaving a redirect

2018-04-18T10:31:36Z

Obst moved page Mike's fancy plots to Mike's fancy plots (xfig example) without leaving a redirect

Obst: /* Finally combining all panels */

2018-04-18T10:25:10Z

Finally combining all panels

Obst: /* Mike's fancy plots */

2018-04-18T10:23:27Z

Mike's fancy plots

Obst: /* Mike's fancy plots */

2018-04-18T09:54:24Z

Mike's fancy plots

Obst: Created page with "--- //example added by Manfred// --- = Mike's fancy plots = A passionate ISIS user (hereafter referred to as "Mike") has compiled a set of http://space.mit.edu/home/mnowak/..."

2018-04-18T09:50:23Z

Created page with "--- //example added by Manfred// --- = Mike's fancy plots = A passionate ISIS user (hereafter referred to as "Mike") has compiled a set of http://space.mit.edu/home/mnowak/..."

New page

--- //example added by Manfred// ---

= Mike's fancy plots =
A passionate ISIS user (hereafter referred to as "Mike")
has compiled a set of [[http://space.mit.edu/home/mnowak/isis_vs_xspec/plots.html|powerful plot functions]]
for spectra in ISIS, which rely on its PGPLOT module.
To my knowledge, these functions do currently (2011 April) not yet offer an interface
to extract the data directly within ISIS.
However, there is a ''write_plot'' function that //"creates ASCII files
with the data from the last plot made using//
''plot_counts'', ''plot_data'', ''plot_residuals'', //or// ''plot_unfold''."
The numbers that appear in any of these fancy plots can hence be dumped to a file,
subsequently read from S-Lang, and plotted using [[isis:slxfig|SLxfig]].

Ideally, one should at some point modify Mike's fancy plot functions,
creating a public interface for the data.
But for the time being, using ''write_plot'' is a simple workaround,
and it's actually quite useful to have the data, which shall be displayed in your plots,
stored in plain ASCII files -- independent of the full set of spectra and response matrices etc.

=== If get an error message ===
fpd is undefined
(or the like) from ''write_plot'', you are most probably using an old version of Mike's scripts (0.X),
which is not supported any more,
and should therefore seriously consider to change to version 1.x -- even though it is not 100% backward compatible.

= An example =
This example is based on Fig. 2.79 of a document I was preparing recently...

As usually, it is recommended to store all files
needed to generate //a single plot// for your thesis, talk, publication, etc.
in //a single subdirectory//.
This advice also holds for the script that loads the data:

== Loading the data ==
You certainly have already a script that loads your data, don't you?
It's nevertheless reasonable to copy this file to the directory with the files needed for the plot,
since you might modify the load data script in your analysis script,
but will once have to find out what (the heck) you did (a long time before) to create (exactly) this plot.
% loaddata.sl

require("isisscripts");

variable dir = getcwd(); % store current working directory
()=chdir("/path/to/your/data"); % go to your analysis directory
variable pn = load_data("src_sd.pha"); % load the data (obviously)
()=chdir(dir); % go back to the previous directory

group(pn; min_sn=20, bounds=0.6, unit="keV"); % do all the setup...
ignore_en(pn, , 0.6);
ignore_en(pn, 5, );

Even though ''loaddata.sl'' is here only needed to get the data (next section),
it might be useful to keep it separate.

== Getting the data / model / residuals ==
I want to compare three models, for which I copied the parameter files
''powerlaw.par'', ''powerlaw_diskbb.par'', and ''powerlaw_diskbb_gainshift.par''
to the directory storing all information how the plot was created.

()=evalfile("loaddata.sl");

variable m;
foreach m (["powerlaw"+["", "_diskbb"+["", "_gainshift"]]]) % a Horner scheme for file names? ;-)
{
message(m);
load_par(m+".par"); % load the model
()=eval_counts(); % evaluate the model
plot_data(1; res=1); % call Mike's fancy plot function <==
write_plot(m); % dump data into ASCII file <==
}

The statements marked by <code> <== </code> create the ASCII data files
''powerlaw_0.dat'', ''powerlaw_diskbb_0.dat'', and ''powerlaw_diskbb_gainshift_0.dat''

in the following format
# Plot Type : plot_data
#
# X-axis : \frEnergy (keV)
# Y-Axis : \frCounts s\u-1\d keV\u-1\d
# Residual : \fr\gx
#
# Index- 1, Model Plotted- yes, Data File: src_sd.pha/part=0/order=0, Scale: 1.0000e+00
#
# Background included in data and model- no
#
# XAXIS_COLS : 1 set (X 2 columns - bin_lo, bin_hi)
# DATA_COLS : 1 set (x 2 columns - Data, Error)
# MODEL_COLS : 1 set (x 1 column - Value)
# RESIDUAL_COLS: 1 set (x 3 columns - Mean, Mean-1Sigma, Mean+1Sigma)

{| class="wikitable"
|-
! X_LO
! X_HI
! DATA_VALUE
! DATA_ERROR
! MODEL_VALUE
! RES_VALUE
! RES_VAL-ERR
! RES_VAL+ERR
|-
| ...
| ...
| ...
| ...
| ...
| ...
| ...
| ...
|}

== Reading and plotting the data ==
Now I can go on with these ASCII data files only, even on my own laptop.

==== Some global definitions ====
require("isisscripts");

variable file = ["powerlaw"+["", "_diskbb"+["", "_gainshift"]]];
variable chi2red = [`7.56`, `2.18`, `1.09`];

variable W=14, H=6, % width and height of first panel
Emin=.55, Emax=5.5, % common xrange of all panels
Rmin=-11, Rmax=9, % common yrange of the residual panels
Maj=[.6, 1, 1.5, 2, 3, 4, 5], % major tic marks (energy axis)
Min=[[.5:2:.1], [2:6:.2]]; % minor tic marks (energy axis)
xfig_new_color("col0", 0xC00000);
xfig_new_color("col1", 0x0000E0);
xfig_new_color("col2", 0x00A000);

==== Data panel ====
edit: alternatively to the example below using readcol, you can also use

Structure str = read_data_from_write_plot(String_Type file)

...with readcol:

% reading four columns of the first ASCII data file (X_LO, X_HI, DATA_VALUE, DATA_ERROR)
variable lo, hi, val, err; (lo, hi, val, err) = readcol(file[0]+"_0.dat", 1, 2, 3, 4);

variable X = xfig_plot_new(W, H);
X.world(Emin, Emax, min_max(val/1e3; logpad=0.05); xlog, ylog);
X.xaxis(; major=Maj, minor=Min, ticlabels2=0);
X.ylabel(`Count Rate [$10^3\rm\,s^{-1}\,keV^{-1}$]`);

% plotting the data
X.plot((lo+hi)/2., val/1e3, err/1e3; sym="point", size=-1, eb_factor=0, depth=100);

% I add some labels using the world0 coordinate system (in which the plot window extends from 0 to 1).
% I use `` as string delimiters for LaTeX code and avoid to have to escape every backslash (as in "\\").
X.xylabel(.05, .3, `\tt TBnew\:*\: powerlaw`, -.5, 0; world0, color="col0");
X.xylabel(.05, .2, `\tt TBnew\:*\:(powerlaw\,+\,diskbb)`, -.5, 0; world0, color="col1");
X.xylabel(.05, .1, `\tt TBnew\:*\:(powerlaw\,+\,diskbb) $\otimes$ gainshift`, -.5, 0; world0, color="col2");
X.xylabel(.01, 1, `a)`, -.5, .8; world0);

==== Loop over the models, creating the panels with the residuals ====
variable R = Struct_Type[3];
variable i;
_for i (2, 0, -1)
{
% read (MODEL_VALUE, RES_VALUE) from the corresponding ASCII data file:
variable m, r; (m, r) = readcol(file[i]+"_0.dat", 5, 6);
% add model to the data panel
X.plot((lo+hi)/2., m/1e3; color="col$i"$);

% create residual panel
R[i] = xfig_plot_new(W, H/3);
R[i].world(Emin, Emax, Rmin, Rmax; xlog);
R[i].yaxis(; major=[-10:10:5], minor=[-20:10], ticlabels2=0);
R[i].xaxis(; major=Maj, minor=Min);
R[i].ylabel(`$\chi$`; color="col$i"$);
R[i].plot((lo+hi)/2., r, 1; sym="point", size=-1, eb_factor=0, color="col$i"$);
R[i].plot([0,1], [0,0]; line=1, color="#808080", depth=100, world01);
R[i].xylabel(.98, 0, `\footnotesize$\chi^2_{\rm red}=`+chi2red[i]+"$", .5, -1; world0, color="col$i"$);
R[i].xylabel(.01, 1, sprintf("%c)", 'b'+i), -.5, .8; world0);
}

==== Panel with the flux model and its components ====
variable M = xfig_plot_new(W, H*.75);
M.world(Emin, Emax, 0.03, 1.2; xlog, ylog);
M.xaxis(; major=Maj, minor=Min, ticlabels2=0);
M.xlabel(`Energy [keV]`);
M.ylabel(`Flux \footnotesize[s$^{-1}$\,cm$^{-2}$\,keV$^{-1}$]`);
M.xylabel(.01, 1, `e)`, -.5, .8; world0);

% evaluate the model on a fine grid
load_par(file[-1]+".par");
variable Lo, Hi; (Lo, Hi) = log_grid(.6, 5, 1000);
variable y = eval_fun_keV(Lo, Hi)/(Hi-Lo); % flux in photons/s/cm^2/keV
variable e = sqrt(Lo*Hi); % Since the grid is logarithmic, I tend to use the geometric mean energy.
M.plot(e, y; xlog, ylog, color="col2", depth=90, width=3);

% disk flux
variable coldisk = "#908000";
variable ff = get_fit_fun();
fit_fun(strreplace(ff, "powerlaw(1)+", ""));
y = eval_fun_keV(Lo, Hi)/(Hi-Lo);
M.plot(e, y; color=coldisk, line=3, width=3);

% powerlaw flux
variable colpow = "#006080";
fit_fun("powerlaw(1)*tbnew_simple(1)");
y = eval_fun_keV(Lo, Hi)/(Hi-Lo);
M.plot(e, y; color=colpow, line=1, width=3);

% some labels
M.xylabel(.98, .85, `\small$N_{\tt powerlaw} = 1.98\pm0.04$`, .5, 0; world0, color=colpow);
M.xylabel(.98, .7, `\small$\Gamma = 1.774\pm0.016$`, .5, 0; world0, color=colpow);
M.xylabel(.98, .3, `\small$N_{\tt diskbb} = (1.4^{+0.3}_{-0.2})\times10^{5}$`, .5, 0; world0, color=coldisk);
M.xylabel(.98, .15, `\small$kT_{\rm in} = (0.227\pm0.006)\rm\,keV$`, .5, 0; world0, color=coldisk);
M.xylabel(.02, .15, `\small$N_{\rm H}=(7.82\pm0.15)\times10^{21}\rm\,cm^{-2}$`, -.5, 0; world0, color="col2");

==== Finally combining all panels ====
simply using [[isis:slxfig:multiplot|xfig_multiplot]]:
xfig_multiplot(X, R, M; xlabel=`Energy [keV]`).render("EPIC-pn_spectrum.pdf");

=== The result ===
[[File:fancy_plot.png]]

Conclusion: An EPIC-pn burst mode spectrum may require a ''gainshift'' correction!

[[Category:SLxfig]]

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

@@ Line 199: / Line 199: @@
 ==== Finally combining all panels ====
-simply using [[isis:slxfig:multiplot|xfig_multiplot]]:
+simply using [[Compounds and multi-panel plots (xfig example)|xfig_multiplot]]:
    xfig_multiplot(X, R, M; xlabel=`Energy [keV]`).render("EPIC-pn_spectrum.pdf");

@@ Line 2: / Line 2: @@
 = Mike's fancy plots =
-A passionate ISIS user (hereafter referred to as "Mike")
+A passionate ISIS user (hereafter referred to as "Mike") has compiled a set of [http://space.mit.edu/home/mnowak/isis_vs_xspec/plots.html powerful plot functions]
-has compiled a set of [http://space.mit.edu/home/mnowak/isis_vs_xspec/plots.html powerful plot functions]
+for spectra in ISIS, which rely on its PGPLOT module. To my knowledge, these functions do currently (2011 April) not yet offer an interface
-for spectra in ISIS, which rely on its PGPLOT module.
+to extract the data directly within ISIS. However, there is a <code>write_plot</code> function that ''"creates ASCII files with the data from the last plot made using'' <code>plot_counts</code>, <code>plot_data</code>, <code>plot_residuals</code>, ''or'' <code>plot_unfold</code>." The numbers that appear in any of these fancy plots can hence be dumped to a file, subsequently read from S-Lang, and plotted using [[Plotting with the S-Lang Xfig module|SLxfig]].
-Ideally, one should at some point modify Mike's fancy plot functions,
+Ideally, one should at some point modify Mike's fancy plot functions, creating a public interface for the data. But for the time being, using <code>write_plot</code> is a simple workaround, and it's actually quite useful to have the data, which shall be displayed in your plots, stored in plain ASCII files -- independent of the full set of spectra and response matrices etc.
 === If get an error message ===
    fpd is undefined
-(or the like) from ''write_plot'', you are most probably using an old version of Mike's scripts (0.X),
+(or the like) from <code>write_plot</code>, you are most probably using an old version of Mike's scripts (0.X), which is not supported any more, and should therefore seriously consider to change to version 1.x -- even though it is not 100% backward compatible.
 = An example =