%
% file: analysis_lines.txt
%


%------------------
% get the model basics
%
% FH's model modified for HETG
%%  .source e0102_fh2_model
.source e0102_cal_model

% Look at the spatial components:
s3d_proj_comps;

% Get the line image array(s):
.load e0102_line_image

% Just need this one for the 'fh2 model:
E0102_All_v3d = e0102_line_image([0.4,5.5], NULL, 1.0, 0);
% the 'cal model uses all of these line-image arrays:
E0102_LowE_v3d  = e0102_line_image([0.4,0.5], NULL, 1.0, 0);
E0102_OVII_v3d  = e0102_line_image([0.5,0.6], [0.4,0.5], 3.0, 0);
E0102_OVIII_v3d = e0102_line_image([0.62,0.70], [0.73,0.76], 2.0, 0);
E0102_FeXVII_v3d    = e0102_line_image([0.75,0.85], [0.85,0.87], 1.5, 0);
E0102_NeIX_v3d  = e0102_line_image([0.87,0.96], [0.97,0.99], 2.0, 0);
E0102_NeX_v3d   = e0102_line_image([0.98,1.15], [1.15,1.25], 2.0, 0);
E0102_MgSiS_v3d = e0102_line_image([1.25,2.6], [2.7,5.5], 5.0, 0);

% evaluate new model and look at it

s3d_update;
s3d_list_model;
s3d_proj_comps;


%------------------
% get the data:
.source e0102_data


% On fast Linux box use:
e2d_meta[1].oversim=25.0;
e2d_meta[2].oversim=25.0;
e2d_meta[3].oversim=25.0;
e2d_meta[4].oversim=25.0;
e2d_meta[5].oversim=25.0;
e2d_meta[6].oversim=25.0;
e2d_meta[7].oversim=25.0;
e2d_meta[8].oversim=25.0;
e2d_meta[10].oversim=9.0;
e2d_meta[11].oversim=9.0;
e2d_meta[12].oversim=9.0;
e2d_meta[13].oversim=9.0;

% ignore the non-dispersed data
e2d_ignore(10);
e2d_ignore(11);
e2d_ignore(12);
e2d_ignore(13);

e2d_list_data;


% evaluate the model
e2d_pars2fom;



% - - - Global Images for Comparison - - -

variable iloop, ithis;
_for iloop (1,8,1)
  {
  e2d_meta[iloop].scisigma=0.05;
  e2d_view[iloop].chigreen=1.5;
  e2d_view[iloop].nsmoo=2;
  };

% Nominal dispersed binning and views for global FOM computation:
variable npixbin=3.0;
e2d_set_view (1, [ 4.60 , 23.0 , npixbin*0.0111] , [ -28.0 , 26.0 , npixbin*0.492 ] );
e2d_set_view (2, [ 4.60 , 23.0 , npixbin*0.0111] , [ -28.0 , 26.0 , npixbin*0.492 ] );
e2d_set_view (3, [ 4.60 , 23.0 , npixbin*0.0111] , [ -28.0 , 26.0 , npixbin*0.492 ] );
e2d_set_view (4, [ 4.60 , 23.0 , npixbin*0.0111] , [ -28.0 , 26.0 , npixbin*0.492 ] );
e2d_set_view (5, [ 4.60 , 16.7 , npixbin*0.0111/2.] , [ -28.0 , 26.0 , npixbin*0.492 ] );
e2d_set_view (6, [ 4.60 , 16.7 , npixbin*0.0111/2.] , [ -28.0 , 26.0 , npixbin*0.492 ] );
e2d_set_view (7, [ 4.60 , 16.7 , npixbin*0.0111/2.] , [ -28.0 , 26.0 , npixbin*0.492 ] );
e2d_set_view (8, [ 4.60 , 16.7 , npixbin*0.0111/2.] , [ -28.0 , 26.0 , npixbin*0.492 ] );

e2d_list_fom;


% Look at Data/Model 2D residuals:
e2d_view_resid(7);


% Print number of events in the data sets in the above ranges:
_for ithis (1,8,1)
{
  message( "Dataset "+string(ithis)+" counts in range = "+
  string(sum(e2d_data[ithis].dhist)) );
}

%Dataset 1 counts in range = 35572
%Dataset 2 counts in range = 26626
%Dataset 3 counts in range = 28619    bkg: 680
%Dataset 4 counts in range = 21509    bkg: 679
%Dataset 5 counts in range = 10189
%Dataset 6 counts in range = 7566
%Dataset 7 counts in range = 7978     bkg: 963
%Dataset 8 counts in range = 5982     bkg: 540

% Counts in 0.5 - 2.0 keV:
e2d_set_view (3, [ Const_hc/2.0 , Const_hc/0.5 , npixbin*0.0111] , [ -28.0 , 26.0 , npixbin*0.492 ] );
e2d_view_data;   %  28416   bkg:  657  2.31%
e2d_set_view (4, [ Const_hc/2.0 , Const_hc/0.5 , npixbin*0.0111] , [ -28.0 , 26.0 , npixbin*0.492 ] );
e2d_view_data;   %  21188   bkg:  591  2.79%
e2d_set_view (7, [ Const_hc/2.0 , Const_hc/0.5 , npixbin*0.0111] , [ -28.0 , 26.0 , npixbin*0.492 ] );
e2d_view_data;   %   8009   bkg:  925 11.55%
e2d_set_view (8, [ Const_hc/2.0 , Const_hc/0.5 , npixbin*0.0111] , [ -28.0 , 26.0 , npixbin*0.492 ] );
e2d_view_data;   %   5732   bkg:  434  7.57%



% - - - Global comparison - - -

% Set parameters to determine the Chi^2 etc.

variable iloop, ithis;
_for iloop (1,8,1)
  {
  e2d_meta[iloop].scisigma=0.000001;
  e2d_view[iloop].chigreen=1.5;
  e2d_view[iloop].nsmoo=1;
  };

% Nominal dispersed binning and views for global FOM computation:
variable npixbin=40.0;
e2d_set_view (1, [ 4.60 , 23.0 , npixbin*0.0111] , [ -28.0 , 26.0 , npixbin*0.492 ] );
e2d_set_view (2, [ 4.60 , 23.0 , npixbin*0.0111] , [ -28.0 , 26.0 , npixbin*0.492 ] );
e2d_set_view (3, [ 4.60 , 23.0 , npixbin*0.0111] , [ -28.0 , 26.0 , npixbin*0.492 ] );
e2d_set_view (4, [ 4.60 , 23.0 , npixbin*0.0111] , [ -28.0 , 26.0 , npixbin*0.492 ] );
e2d_set_view (5, [ 4.60 , 16.7 , npixbin*0.0111/2.] , [ -28.0 , 26.0 , npixbin*0.492 ] );
e2d_set_view (6, [ 4.60 , 16.7 , npixbin*0.0111/2.] , [ -28.0 , 26.0 , npixbin*0.492 ] );
e2d_set_view (7, [ 4.60 , 16.7 , npixbin*0.0111/2.] , [ -28.0 , 26.0 , npixbin*0.492 ] );
e2d_set_view (8, [ 4.60 , 16.7 , npixbin*0.0111/2.] , [ -28.0 , 26.0 , npixbin*0.492 ] );

e2d_list_fom;

% v1p9 as is:
%      o Data-Model 1,  chi^2 = 798.081  from 126 cells;   chi^2/cell = 6.33398;  sum(D)/sum(M) = 0.961844  +/- 0.00509977
%      o Data-Model 2,  chi^2 = 674.744  from 125 cells;   chi^2/cell = 5.39795;  sum(D)/sum(M) = 0.985968  +/- 0.00604241

%      o Data-Model 3,  chi^2 = 366.436  from 125 cells;   chi^2/cell = 2.93149;  sum(D)/sum(M) = 1.03036  +/- 0.00609063
%      o Data-Model 4,  chi^2 = 343.361  from 121 cells;   chi^2/cell = 2.8377;  sum(D)/sum(M) = 1.02185  +/- 0.00696747

%      o Data-Model 5,  chi^2 = 626.315  from 165 cells;   chi^2/cell = 3.79585;  sum(D)/sum(M) = 1.18565  +/- 0.011746
%      o Data-Model 6,  chi^2 = 499.951  from 154 cells;   chi^2/cell = 3.24644;  sum(D)/sum(M) = 1.15107  +/- 0.0132333
%      o Data-Model 7,  chi^2 = 368.696  from 163 cells;   chi^2/cell = 2.26194;  sum(D)/sum(M) = 1.18387  +/- 0.0132543
%      o Data-Model 8,  chi^2 = 310.86  from 151 cells;   chi^2/cell = 2.05867;  sum(D)/sum(M) = 1.15  +/- 0.0148687
%
% v1p9 with FIT adjustments:
%      o Data-Model 1,  chi^2 = 828.394  from 126 cells;   chi^2/cell = 6.57456;  sum(D)/sum(M) = 0.947842  +/- 0.00502553
%      o Data-Model 2,  chi^2 = 698.708  from 125 cells;   chi^2/cell = 5.58967;  sum(D)/sum(M) = 0.980451  +/- 0.0060086

%      o Data-Model 3,  chi^2 = 351.315  from 125 cells;   chi^2/cell = 2.81052;  sum(D)/sum(M) = 1.01877  +/- 0.0060221
%      o Data-Model 4,  chi^2 = 348.805  from 121 cells;   chi^2/cell = 2.88269;  sum(D)/sum(M) = 1.01737  +/- 0.00693693

%      o Data-Model 5,  chi^2 = 632.944  from 165 cells;   chi^2/cell = 3.83602;  sum(D)/sum(M) = 1.18205  +/- 0.0117104
%      o Data-Model 6,  chi^2 = 498.496  from 154 cells;   chi^2/cell = 3.23699;  sum(D)/sum(M) = 1.151  +/- 0.0132325
%      o Data-Model 7,  chi^2 = 360.772  from 163 cells;   chi^2/cell = 2.21332;  sum(D)/sum(M) = 1.17864  +/- 0.0131957
%      o Data-Model 8,  chi^2 = 306.265  from 151 cells;   chi^2/cell = 2.02825;  sum(D)/sum(M) = 1.1511  +/- 0.014883


% 1d plots, residuals and chi^2 :

PSwin=open_plot("e0102_megm_"+fit_str+"1d.ps/CPS");
_pgsvp(0.20, 0.85, 0.20, 0.85);
e2d_plot_1d(3,1);
close_plot(PSwin);
PSwin=open_plot("e0102_megm_"+fit_str+"1drat.ps/CPS",1,3);
_pgsvp(0.20, 0.85, 0.20, 0.85);
charsize(3.0);
set_frame_line_width(2) ;
e2d_plot_1drat(3,1);
close_plot(PSwin);

PSwin=open_plot("e0102_megp_"+fit_str+"1d.ps/CPS");
_pgsvp(0.20, 0.85, 0.20, 0.85);
e2d_plot_1d(4,1);
close_plot(PSwin);
PSwin=open_plot("e0102_megp_"+fit_str+"1drat.ps/CPS",1,3);
_pgsvp(0.20, 0.85, 0.20, 0.85);
charsize(3.0);
set_frame_line_width(2) ;
e2d_plot_1drat(4,1);
close_plot(PSwin);


% Look at the K-S results for each data set:

%% .load e2d_plot_ksWGT
%% _for iloop (1,8,1) {e2d_plot_ksWGT(iloop,1);};

% v1p9 as is:
%     o Data-Model 3,  K-S on TG_LAM: Max |diff| is 0.0150856  at 16.2582 Angstrom
%                             ks_test2 statistic of 0.0150856  gives p-value of 7.37855e-06
%      o D/M ratio = 1.03102
%      o Data counts = 28619
%      o Data-Model 4,  K-S on TG_LAM: Max |diff| is 0.0257091  at 9.27127 Angstrom
%                             ks_test2 statistic of 0.0257091  gives p-value of 2.73848e-12
%      o D/M ratio = 1.02156
%      o Data counts = 21509

% FIT version:
%      o Data-Model 3,  K-S on TG_LAM: Max |diff| is 0.0113296  at 8.16428 Angstrom
%                             ks_test2 statistic of 0.0113296  gives p-value of 0.00171853
%      o D/M ratio = 1.01877
%      o Data counts = 28619
%      o Data-Model 4,  K-S on TG_LAM: Max |diff| is 0.0254588  at 9.27127 Angstrom
%                             ks_test2 statistic of 0.0254588  gives p-value of 4.6293e-12
%      o D/M ratio = 1.01737
%      o Data counts = 21509



%
% - - -


.source fill_regions.sl

.source plot_regions.sl

%
message(" *** Now do:  .source write_regions.sl  ");
message("     and do isis> !cp load_v1p5_regions_EMPTY.sl "+
  "load_v1p9_regions.sl ");
message("     and    isis> ! cat >> load_v1p9_regions.sl ");
message("");



% - - -
%  Merge every thing into one MEG and one HEG
% MEGs:
e2d_merge_data(0, [1,2,3,4], [0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0]);
% HEGs:
e2d_merge_data(9, 4+[1,2,3,4], [0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0]);
%
% Merge the models too:
e2d_data[0].mxs = [ 
  e2d_data[1].mxs,e2d_data[2].mxs,e2d_data[3].mxs,e2d_data[4].mxs ];
e2d_data[0].mys = [
  e2d_data[1].mys,e2d_data[2].mys,e2d_data[3].mys,e2d_data[4].mys ];
e2d_meta[0].modvalid=1;

e2d_data[9].mxs = [ 
  e2d_data[5].mxs,e2d_data[6].mxs,e2d_data[7].mxs,e2d_data[8].mxs ];
e2d_data[9].mys = [
  e2d_data[5].mys,e2d_data[6].mys,e2d_data[7].mys,e2d_data[8].mys ];
e2d_meta[9].modvalid=1;

iloop=0;
e2d_meta[iloop].scisigma=0.05;
e2d_view[iloop].chigreen=1.5;
e2d_view[iloop].nsmoo=2;
iloop=9;
e2d_meta[iloop].scisigma=0.05;
e2d_view[iloop].chigreen=1.5;
e2d_view[iloop].nsmoo=2;

% Nominal dispersed binning and views for global FOM computation:
variable npixbin=3.0;
e2d_set_view (0, [ 4.60 , 23.0 , npixbin*0.0111] , [ -28.0 , 26.0 , npixbin*0.492 ] );
e2d_set_view (9, [ 4.60 , 16.7 , npixbin*0.0111/2.] , [ -28.0 , 26.0 , npixbin*0.492 ] );

e2d_view_resid(0);
e2d_view_resid(9);



% - - -
% Look at the ARFs compared with where the lines are
%   include ACIS-S3 too:
e2d_setup_arf(14,"aciss_aimpt_cy10.arf");

variable PSwin = open_plot("compare_arfs.ps/CPS");

%% variable grat_clrs = [6,8,4,5];

xrange(0.3,2.0); xlog;
yrange(0.1,1.e3); ylog;
label("Energy (keV)","ARF value (cm2 x counts/photon)","");
hplot(e2d_meta[14].arfelo,s3d_hi_from_lo(e2d_meta[14].arfelo),e2d_meta[14].arfval,1);
ohplot(s3d_spec[10].es, s3d_hi_from_lo(s3d_spec[10].es),
  300.0*(s3d_spec[10].fs+s3d_spec[9].fs),15);
ohplot(e2d_meta[3].arfelo,s3d_hi_from_lo(e2d_meta[3].arfelo),e2d_meta[3].arfval,6);
ohplot(e2d_meta[4].arfelo,s3d_hi_from_lo(e2d_meta[4].arfelo),e2d_meta[4].arfval,8);
ohplot(e2d_meta[7].arfelo,s3d_hi_from_lo(e2d_meta[7].arfelo),e2d_meta[7].arfval,4);
ohplot(e2d_meta[8].arfelo,s3d_hi_from_lo(e2d_meta[8].arfelo),e2d_meta[8].arfval,5);
color(15);
xylabel(1.4,12.0,"model flux");
color(6);
xylabel(0.39,0.8,"MEG -1");
color(8);
xylabel(0.45,0.15,"MEG +1");
color(4);
xylabel(0.60,0.15,"HEG -1");
color(5);
xylabel(0.80,0.15,"HEG +1");
color(1);
xylabel(0.35,200.0,"ACIS-S3 Cy10");

close_plot(PSwin);