%% Solution to "Spectral Analysis of High Resolution X-ray Binary Data"
%%    Exercise by M. Nowak
%%    http://cxc.cfa.harvard.edu/xrayschool/analysis/mnowak_example_I.pdf
%%
%% ISIS script by L. Corrales
%% 2015.08.17 - lia@space.mit.edu

%%--------------------------------------------------------------------------------------------
%% BEFORE ANYTHING: Use TGcat (http://tgcat.mit.edu) to download gratings data for ObsId 5461

%%--------------------------------------------------------------------------------------------
%% 1. Load the TGcat data into ISIS

% We can load the MEG and HEG +/- 1 orders by specifying what rows of the PHA file to read in
%   Rows 3 and 4 are HEG -1 and +1; Rows 9 and 10 are MEG -1 and +1

variable HEG1 = load_data("pha2.gz", [3,4]);
variable MEG1 = load_data("pha2.gz", [9,10]);

% Show what we have loaded so far
% Note the "id" column, which shows the integer value used to identify each particular dataset

list_data;

% Notice that the arf and rmf were not automatically loaded.
% They will have to be read in separately and assigned to the proper data set

() = load_arf("heg_-1.arf.gz");
() = load_arf("heg_1.arf.gz");
() = load_arf("meg_-1.arf.gz");
() = load_arf("meg_1.arf.gz");

% Examine the loaded arfs
list_arf;

% Now assign them to the datasets - be sure that the ids match properly!
assign_arf([1:4], [HEG1, MEG1]);
list_data;

% Now do the same for the rmfs

() = load_rmf("heg_-1.rmf.gz");
() = load_rmf("heg_1.rmf.gz");
() = load_rmf("meg_-1.rmf.gz");
() = load_rmf("meg_1.rmf.gz");

list_rmf;
assign_rmf([1:4], [HEG1, MEG1]);
list_data;

%%--------------------------------------------------------------------------------------------
%% 2. Play with plotting the data
%%    Unfortunately, the default ISIS plotting tools are very limited.
%%    That's why we asked you to download the Remeis ISIS scripts
%%    http://www.sternwarte.uni-erlangen.de/isis/

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

%% Here's how to do counts/bin with the default ISIS plot commands
plot_unit("Angstrom");
plot_bin_integral;
xrange(1.5,30);
plot_data_counts( MEG1[0], 2 );  % MEG -1 data (red)
oplot_data_counts( HEG1[0], 4 ); % HEG -1 data (blue)

%% Here's how to do it with the ISIS scripts
fancy_plot_unit("Angstrom");
plot_counts( MEG1[0]; dsym=0, dcol=2, decol=15 );
plot_counts( HEG1[0]; dsym=0, dcol=4, decol=15, oplt=1 );

%% Now for counts/bin/sec
plot_data( MEG1[0]; dsym=0, dcol=2, decol=15 );
plot_data( HEG1[0]; dsym=0, dcol=4, decol=15, oplt=1 );


%%--------------------------------------------------------------------------------------------
%% 3. Play with combining the data

%% **VERY IMPORTANT** -- First, the datasets need to be placed on the same spectral grid

match_dataset_grids( [MEG1, HEG1] );

%% This will improve the fit statistic, and will give correct values for overall reduced chi^2
%%    see:   help combine_datasets
variable cid = combine_datasets( [MEG1, HEG1] );

%% Now group the data so that the *combined* s/n is >=5 and at least two channels per bin
%%    see:   help group
group( [MEG1, HEG1]; min_sn=5, min_chan=2 );

%%--------------------------------------------------------------------------------------------
%% 4. Plot the flux corrected spectrum.  This first requires running fluxcorr

flux_corr( [MEG1, HEG1] );

%% The default ISIS way of plotting flux corrected spectrum
plot_bin_density; ylog;
plot_data_flux( MEG1[0], 2 );
oplot_data_flux( HEG1[0], 4 );
% (notice that the endpoints are background dominated)

%% Let's zoom in on 12-16 Angs spectrum, then plot using the Remeis ISIS scripts way
xrange(12, 16);
plot_unfold( MEG1[0]; dsym=0, dcol=2, decol=15 );
plot_unfold( HEG1[0]; dsym=0, dcol=4, decol=15, oplt=1 );

%% One great feature of the Remeis plotting tools -- you can plot the combined data
%%    by giving it a list of dataset ids
%%    (This is correct, only if the dataset grids have been matched)
plot_unfold( [MEG1, HEG1]; dsym=0, dcol=1, decol=15 );

%%--------------------------------------------------------------------------------------------
%% Fitting fun!  "You should notice an absorption edge in these data, as well as 
%%    several prominent absorption lines".  Let's fit them!
%%--------------------------------------------------------------------------------------------


%%--------------------------------------------------------------------------------------------
%% 5. Restrict the range of the noticed data and use a power law to
%%    fit the local continuum model

xnotice( [MEG1, HEG1], 13.5, 15 ); % Angs
list_data;

fit_fun("powerlaw(1)");
() = fit_counts;

variable pstyle = struct{ dsym=0, dcol=1, decol=15, mcol=2 };
xrange(13.5, 15);
plot_unfold( [MEG1, HEG1], pstyle; res=6 );
% See the plot_unfold help for information on plotting residual options
%    isis> plot_unfold;

%%--------------------------------------------------------------------------------------------
%% 6. Add an edge model to the data

fit_fun("powerlaw(1)*edge(1)");

% Let's check out what the parameters are
list_par;

% We need to give it a preliminary value for the edge position
%    From the plot, it looks like my edge is around 14.3 Angstroms
%    You can convert to keV using the Const_hc parameter
%
% I'm also going to set the min and max possible values to be within a reasonable range
%    see 'help set_par'
set_par("edge(1).edgeE", Const_hc/14.3, 0, Const_hc/14.5, Const_hc/13.8);

% Now fit it!
() = fit_counts;

% ... and plot
plot_unfold( [MEG1, HEG1], pstyle; res=6 );


%%--------------------------------------------------------------------------------------------
%% 7. Add Gaussian to fit the most prominent absorption feature

fit_fun("powerlaw(1)*edge(1) - gaussian(1)");
list_par;

set_par("gaussian(1).LineE", Const_hc/14.6, 0, Const_hc/14.7, Const_hc/14.55);

%% You can freeze parameters according to name
%    I like to freeze the original model before fitting the line
%    You can freeze all of the parameters for a given model using a *
freeze("powerlaw(1).*");
freeze("edge(1).*");
() = fit_counts;

% Then thaw and fit again to get best fit
thaw("powerlaw(1).*");
thaw("edge(1).*");
() = fit_counts;

plot_unfold( [MEG1, HEG1], pstyle; res=6 );

%%--------------------------------------------------------------------------------------------
%% 8. Incorporate the next three additional absorption lines

% I see three prominent lines at 14.3, 14.5, and 14.8
% The latter two are easier to get at, so do them first
fit_fun("powerlaw(1)*edge(1) - gaussian(1) - gaussian(2) - gaussian(3)");

% You can also freeze parameters by parameter id index
%    Use list_par and make note of the "idx" column
list_par;
freeze([1:7]);

% Make the suggest lines similarly thin to the first absorption feature
% The 14.5 Angs line:
set_par("gaussian(2).LineE", Const_hc/14.5, 0, Const_hc/14.55, Const_hc/14.45);
set_par("gaussian(2).Sigma", get_par("gaussian(1).Sigma"));

% The 14.8 Angs line:
set_par("gaussian(3).LineE", Const_hc/14.8, 0, Const_hc/14.90, Const_hc/14.70);
set_par("gaussian(3).Sigma", get_par("gaussian(1).Sigma"));

() = fit_counts;
plot_unfold( [MEG1, HEG1], pstyle; res=6 );

% Looks pretty good, so finally I need to incorporate the line that
% looks like it is part of the edge structure

fit_fun("powerlaw(1)*edge(1) - gaussian(1) - gaussian(2) - gaussian(3) - gaussian(4)");
list_par;

freeze([1:13]);

% For now, I'm going to freeze the line position at 14.3 Angs
set_par("gaussian(4).LineE", Const_hc/14.3, 1, Const_hc/14.33, Const_hc/14.27);
set_par("gaussian(4).Sigma", get_par("gaussian(1).Sigma"));

() = fit_counts;
thaw("gaussian(4).LineE");
() = fit_counts;

plot_unfold( [MEG1, HEG1], pstyle; res=6 );
% Looks pretty good!

% Thaw everything to get final fit
thaw([1:16]);
() = fit_counts;

% Here I've changed the residual plot to just show me chi values
plot_unfold( [MEG1, HEG1], pstyle; res=4 );

%%--------------------------------------------------------------------------------------------
%% 9. Get error bars

% See 'help conf' and 'help conf_loop'

% Get 1-sigma CI (level = 0 means 68% C.I.) on line widths
variable smin, smax;
(smin,smax) = conf_loop("gaussian(*).Sigma", 0, 0);