Mike's fancy plots (xfig example)

From Remeis-Wiki
Jump to navigation Jump to search

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

Mike's fancy plots

A passionate ISIS user (hereafter referred to as "Mike") has compiled a set of 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 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 <== 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)
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 xfig_multiplot: 

 xfig_multiplot(X, R, M; xlabel=`Energy [keV]`).render("EPIC-pn_spectrum.pdf");

The result

Fancy plot.png

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

Retrieved from "https://www.sternwarte.uni-erlangen.de/wiki/index.php?title=Mike%27s_fancy_plots_(xfig_example)&oldid=1203"