%; Time-stamp: <2002-06-06 19:43:25 dph> 
%; MIT Directory: ~dph/d1/CXO_Data/Capella/Sl
%; File: R_1103.sl
%; Author: D. Huenemoerder
%; Original version: 2002.05.08
%;====================================================================
% version: 1
%
% purpose: Fit a line profile. Simulate many times to examine
%          statistical deviations from the model profile.
%          Do for single component line, use gaussian+constant model.
%          Do for multicomponent line, and use aped plasma model.
%          In each case, the model is convolved by the rmf*arf.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reprocessed Capella ObsID 1103; see ../Sh/Proc_1103.sh
%   de-afterglowed
%   destreaked
%   re-did acis_process_events, tg*, ... tgextract,  fullgarf.
%
%  NOTE: saturation in core of zero-order, but centroid looks good.
%  NOTE: osip not so good on S1.  Some Line fluxes 20% low.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

_auto_declare = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%% Load Data, responses, and assign %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
% Define the data, arfs, rmfs:

variable f_pha = "./Pha/acisf01103_002N002_pha2.fits";
variable pha_rows = [1,2,3,4];       % HEG -1,1, then MEG -1,1 rows

% ARFs in same order as data, for convenience:
%
variable f_arfs = "./Arf/acisf01103_002N002" +
       ["HEG_-1", "HEG_1", "MEG_-1", "MEG_1"] +
       "_garf.fits" ;

% RMFs: HEG, then MEG:
variable f_rmfs = ["./Ard/acisheg1D1999-07-22rmfN0004.fits",
		   "./Ard/acismeg1D1999-07-22rmfN0004.fits"];

% Load the data +-1st orders, HEG and MEG, and associated arfs, and rmfs.
% Then assign the responses to the data histograms.  Don't save return
% values. 
%

() = load_data(f_pha, pha_rows);

() = load_arf( f_arfs[0] );
() = load_arf( f_arfs[1] );
() = load_arf( f_arfs[2] );
() = load_arf( f_arfs[3] );

() = load_rmf( f_rmfs[0] );
() = load_rmf( f_rmfs[1] );

assign_arf( 1, 1);          % HEG -1
assign_arf( 2, 2);          % HEG +1
assign_arf( 3, 3);          % MEG -1
assign_arf( 4, 4);          % MEG +1

assign_rmf( 1, [1,2] );      % HEG
assign_rmf( 2, [3,4] );      % MEG

flux_corr([1:4]); 

%-


%%
% Load the plasma database for use in modeling/line marking.
% Create a plot window.

%
plasma(aped);

variable p1=open_plot("/xwin"); resize(18, 0.6); erase;

%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% Set up some variables %%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+

Fit_Verbose = -1;         % turn off fit output.

variable n = 24;    %  number of runs to do, 1 w/ real data, 23 fake.

%%%%%%%%%%%%%%%%%%
%      n = 3; %%%% FOR QUICK TESTING.....  **** xxxx XXXX DEBUG
%%%%%%%%%%%%%%%%%

variable v = Float_Type[n],        % vrad or wavelength array
         f = Float_Type[n],        % flux array
	 c = Float_Type[n],        % stat array
	 o = Float_Type[n],        % sigma array
	 s = Struct_Type;          % fit stat structure

define Chisq( z )   % utility to reduce chisq from fit stats.
  {
    return z.statistic / (z.num_bins - z.num_variable_params);
  }


% we will make confidence contours for the observed lines, and overlay
% contours from one of the simulations.

variable dcf,dcv, dcs,  % data contour x,y for flux, velocity (or
                        %         wavelength), sigma
         scf, scv, scs, % simulation contour x,y, sigma
	 dfvmap,  sfvmap, sfsmap;  % data and simulation 2D maps.

%-


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%  Fit a line parametrically %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
%
%%%%%%%  Pick an unblended line  and fit with gaussian + constant
%

ignore([1:4]);
variable xlo=14.165,  xhi=14.235;
xnotice(2, xlo, xhi);
xrange(xlo-.02, xhi+.02);
yrange(0,40);
label("Wavelength [\\A]", "Flux [counts/bin]", "ObsID 1103/Capella/HEG +1");

fit_fun( "gauss(1) + poly(1)" );

set_par( "gauss(1).area", 0.016, 0, 0.0, 0.1 );
set_par( "gauss(1).center", 14.205, 0, 14.18, 14.23 );
set_par( "gauss(1).sigma", 0.004, 0, 0.001, 0.008 );

set_par( "poly(1).a0", 0.01, 0, 0, 0.1 );
set_par( "poly(1).a1", 0, 1 );
set_par( "poly(1).a2", 0, 1 );

() = renorm_counts;

window(p1);    rplot_counts(2);

() = fit_counts( &s );
list_par;

()=printf("Conf.sigma: %e, %e\n", conf("gauss(1).sigma"));


% compute contour map; ad hoc limits.
  dcv = conf_grid("gauss(1).center", 14.201, 14.210, 25);
  dcf = conf_grid("gauss(1).area", 6.e-4, 2.e-3, 25);
  dcs = conf_grid("gauss(1).sigma", 0.0005, 0.01, 25);
  dfvmap = conf_map_counts(dcf, dcv);
  dfsmap = conf_map_counts(dcf, dcs);


rplot_counts(2);
freeze( "poly(1).a0" );     

variable p = get_params;    % Copy the parameters into an array of structs.
                            % This will be used to re-instantiate the same
                            % model before the next fakeit.

variable sp1, sp2, sp3;     % save params vs run;

%  save the fit norm, vrad, and stat
%
f[0] = get_par("gauss(1).area");
v[0] = get_par("gauss(1).center");
o[0] = get_par("gauss(1).sigma");
c[0] = Chisq( s );


%%% fakeit n-1 times... and display in mosaic w/ the data.
%%% Save the stat, flux, and centroid from each iteration.

variable p2=open_plot("/xwin",6,4);resize(8*2.54, 8.5/11. ); erase;

xrange(xlo-.02, xhi+.02);
yrange(0,36);
xnotice(2, xlo, xhi);

variable pid=dup_plot("1103_14A_HEG_p1_fakeits.ps/cps",p2);
resize(10*2.54, 8.5/11.);

window(p2);

rplot_counts(2);
window(pid);          % switch to the postscript device...
rplot_counts(2);      % ... and replot it.
window(p2);           % go back to screen.

ignore(2);
assign_arf(2,5);   % assigning to non-existent histogram index creates one.
assign_rmf(1,5);
xnotice(5, xlo, xhi);

variable i;            % loop over fakeits.....................
for (i=0;i<n-1;i++) {
  set_params(p);
  fakeit;
  () = fit_counts( &s );
  rplot_counts(5);   window(pid);   rplot_counts(5);   window(p2);
  f[i+1] = get_par("gauss(1).area");
  v[i+1] = get_par("gauss(1).center");
  o[i+1] = get_par("gauss(1).sigma");
  c[i+1] = Chisq( s );
}


% save results for summary/comparison of param values and ranges;
variable f1 = @f;
variable v1 = @v;
variable o1 = @o;
variable sp1 = @p;

close_plot(pid);
close_plot(p2);

% print out results of all fits (could plot these, or
% do distributions if doing cash statistic and needed distributions
% for confidence limits)

vmessage("Moment of fluxes: \n"); print(moment(f));
vmessage("Moment of center: \n"); print(moment(v));
vmessage("Moment of sigma: \n"); print(moment(o));
vmessage("Moment of stat: \n"); print(moment(c));

% compute contour map for last simulation:
  scv = conf_grid("gauss(1).center", 14.201, 14.21, 25);
  scf = conf_grid("gauss(1).area", 6.e-4, 2.e-3, 25);
  scs = conf_grid("gauss(1).sigma", 0.0005, 0.01, 25);
  sfvmap = conf_map_counts(scf, scv);
  sfsmap = conf_map_counts(scf, scs);

window(p1);
label("Flux", "Wavelength", "Conf Maps, data and sim");
variable line = struct { width, type };
line.width=2; line.type=2;
limits;

plot_conf(dfvmap);  oplot_conf(sfvmap, line);
connect_points(0); pointstyle(4); oplot(f1, v1, 5 ); connect_points(1);

pid=dup_plot("1103_14A_HEG_p1_conf_map_1.ps/vcps",p1); resize(6*2.54,1.0);

plot_conf(dfvmap);  oplot_conf(sfvmap,line);
connect_points(0); pointstyle(4); oplot(f1, v1, 5 ); connect_points(1);

label("Flux", "Sigma", "Conf Maps, data and sim");
plot_conf(dfsmap);  oplot_conf(sfsmap,line);
connect_points(0); pointstyle(4); oplot(f1, o1, 5 ); connect_points(1);

close_plot(pid);

%-

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Fit with APED plasma %%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-
%  These models take longer to evaluate, so depending on the speed of
%  the computer, it may be a few seconds per fit, or much longer.
%
% define the model functions.
%
% aped_fit_models.sl is from dph and defines several plasma functions;
% the one used here is a single temperature plasma, and has parameters
% for temperature, radial velocity, turbulent broadening, and
% normalization.   
%
% Available from:
%
% http://space.mit.edu/CXC/ISIS/contrib/Aped_fit/aped_fit_models.sl
% (and in the isis contrib directory, with a bug)

()=evalfile("aped_fit_models.sl");   % the aped model definitions

% define a single-temperature plasma model
%  (list_par will show the params)

fit_fun("apedx(1)");    % (renamed apedx because of
                        % collision w/ aped variable for database.)

%%%
% the following are the result of an ad hoc fit to Ne X 12A:
%
set_par("apedx(1).norm", 0.018, 0, 0.0001, 0.08);
set_par("apedx(1).temp", 10.^6.8, 1, 1.e6, 3.e7);   % a good enough temperature

% The lines do seem to need more-than-instrumental broadening.  Is
% thermal broadening  significant for Ne X at 12A in the HEG?

use_thermal_profile;                          % turn on broadening, 
set_par("apedx(1).vturb", 0, 1, 0.0, 200);    %  but freeze vturb to 0.


% Since there seems to be a problem with the wavelength scale, we will 
% work around it by letting the radial velocity be free.  (This is OK
% for single orders, but not for joint fits w/ different offsets.)

set_par("apedx(1).vrad", 100, 0, 0, 200);   % vrad

% ignore all but HEG -1 for this test:
%
xlo=12.09;   xhi=12.17;
xnotice(1, xlo, xhi);
ignore([2:5]);

yrange(0,80);

() = renorm_counts;   % adjust the normalization to the data.

list_data;   % inspect lists
list_par;

xrange(xlo-.01, xhi+.01);         % select Ne X + Fe  region
label("Wavelength [\\A]", "Flux [counts/bin]", "ObsID 1103/Capella/HEG -1");
plot_data_counts(1,red);
oplot_model_counts(1,green);

rplot_counts(1);             % residual plot

() = fit_counts;

rplot_counts(1);

% the profile isn't really as bad as it looks. (but the model isn't
%   perfect; residuals aren't quite random);
list_par;


% compute contour map for data fit
  dcv = conf_grid("apedx(1).vrad", 0.0, 150.0, 20);
  dcf = conf_grid("apedx(1).norm", 0.012, 0.023, 20);
  dfvmap = conf_map_counts(dcf, dcv);


%  Create some fake data with the same model, and then fit it.
% do a few times successively.

assign_arf(1,5);   % assigning to non-existent histogram index creates one.
assign_rmf(1,5);
fakeit;                    % uses model to make fake data.
rplot_counts(5);           % look at it.
ignore(1);                 % switch fit to histogram 5:
xnotice(5, xlo, xhi);  

%%%%
% let's do several fakeits and plot in a m*n window...
% First go back to histogram 1 and fit, to put that in the window for
% comparison:

ignore(5); 
xnotice(1, xlo, xhi);

%charsize(1.33);  % a little bigger, but not too much to fall out the box.

p2 = open_plot("/xwin",6,4); resize(8*2.54,8.5/11.); erase;
xrange(xlo, xhi);
yrange(0, 80);

%  set up a ps device.  We will plot to each device.
%  The ps device will inherit from the xwin if the ranges are set first.
%
pid=dup_plot("1103_12A_HEG_m1_fakeits.ps/cps",p2); resize(10.5*2.54, 8.5/11.);
window(p2);


() = fit_counts( &s );   % do the fit, and save the statistics structure.
rplot_counts(1);
window(pid);
rplot_counts(1);
window(p2);

p = get_params;    % Copy the parameters into an array of structs.
                   % This will be used to re-instantiate the same
                   % model before the next fakeit.

%  save the fit norm, vrad, and stat:
%
f[0] = get_par("apedx(1).norm");
v[0] = get_par("apedx(1).vrad");
c[0] = Chisq( s );

ignore(1);
xnotice(5, xlo, xhi);

% Do a series of simulations, each time using the same underlying model.
%
for ( i=0; i<n-1; i++ ) {
  set_params( p );
  fakeit;
  () = fit_counts( &s );
  rplot_counts(5);  window(pid); rplot_counts(5); window(p2);
  f[i+1] = get_par("apedx(1).norm");
  v[i+1] = get_par("apedx(1).vrad");
  c[i+1] = Chisq( s );
}
close_plot(pid);

f2 = @f;
v2 = @v;
sp2 = @p;

vmessage("Moment of fluxes: \n"); print(moment(f));
vmessage("Moment of vrad: \n"); print(moment(v));
vmessage("Moment of stat: \n"); print(moment(c));

% compute contour map for last simulation:
  scv = conf_grid("apedx(1).vrad", 0.0, 150.0, 20);
  scf = conf_grid("apedx(1).norm", 0.012, 0.023, 20);
  sfvmap = conf_map_counts(scf, scv);

window(p1);
label("Flux", "Radial Velocity", "Conf Maps, data and sim");
line.width=2; line.type=2;
limits;

plot_conf(dfvmap);  oplot_conf(sfvmap, line);
connect_points(0); pointstyle(4); oplot(f2, v2, 5 );connect_points(1);

pid=dup_plot("1103_12A_HEG_m1_conf_map_2.ps/vcps",p1); resize(6*2.54,1.0);
plot_conf(dfvmap);  oplot_conf(sfvmap,line);
connect_points(0); pointstyle(4); oplot(f2, v2, 5 );connect_points(1);
close_plot(pid);

%-

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% plot the flux model, to see lines:

window(p1);
ignore(5);
xnotice(1, xlo, xhi);
xrange(xlo, xhi);
yrange(0, 6.e-4);
label("Wavelength [\\A]", "Flux [counts/bin]", "ObsID 1103/Capella/HEG -1");

pid = dup_plot("1103_12A_model_flux_H_m1.ps/cps", p1);  % set up ps
window(p1);                                   % ...but select xwin

set_par("apedx(1).vrad", 0 );  % turn off shift

use_delta_profile;
eval_counts;
plot_model_flux(1,2);
plot_group(brightest(10,where(wl(12.1,12.17))),1);

window(pid);
plot_model_flux(1,2);
plot_group(brightest(10,where(wl(12.1,12.17))),1);
window(p1);

use_thermal_profile;
eval_counts;
oplot_model_flux(1,3);

window(pid);
oplot_model_flux(1,3);

close_plot(pid); window(p1);
close_plot(p2); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Repeat, but without thermal broadening, to see resids:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-

use_delta_profile;
set_params( p );

xlo=12.09;   xhi=12.17;
xnotice(1, xlo, xhi);
ignore([2:5]);
yrange(0,80);

p2 = open_plot("/xwin",6,4); resize(8*2.54,8.5/11.); erase;
xrange(xlo, xhi);
yrange(0, 80);

pid=dup_plot("1103_12A_HEG_m1_fakeits_delta.ps/cps",p2); resize(10.5*2.54, 8.5/11.);
window(p2);

() = fit_counts( &s );   % do the fit, and save the statistics structure.
rplot_counts(1);
window(pid);
rplot_counts(1);
window(p2);

p = get_params;    % Copy the parameters into an array of structs.
                   % This will be used to re-instantiate the same
                   % model before the next fakeit.

%  save the fit norm, vrad, and stat:
%
f[0] = get_par("apedx(1).norm");
v[0] = get_par("apedx(1).vrad");
c[0] = Chisq( s );

% compute contour map for data fit
  dcv = conf_grid("apedx(1).vrad", 0.0, 150.0, 20);
  dcf = conf_grid("apedx(1).norm", 0.012, 0.023, 20);
  dfvmap = conf_map_counts(dcf, dcv);

ignore(1);
xnotice(5, xlo, xhi);

% Do a series of simulations, each time using the same underlying model.
%
for ( i=0; i<n-1; i++ ) {
  set_params( p );
  fakeit;
  () = fit_counts( &s );
  rplot_counts(5);  window(pid); rplot_counts(5); window(p2);
  f[i+1] = get_par("apedx(1).norm");
  v[i+1] = get_par("apedx(1).vrad");
  c[i+1] = Chisq( s );
}
close_plot(pid);

v3 = @v;
f3 = @f;
sp3 = @p;

vmessage("Moment of fluxes: \n"); print(moment(f));
vmessage("Moment of dlambda: \n"); print(moment(v));
vmessage("Moment of stat: \n"); print(moment(c));

% compute contour map for last simulation:
  scv = conf_grid("apedx(1).vrad", 0.0, 150.0, 20);
  scf = conf_grid("apedx(1).norm", 0.012, 0.023, 15);
  sfvmap = conf_map_counts(dcf, dcv);

window(p1);
label("Flux", "Radial Velocity", "Conf Maps, data and sim");
line.width=2; line.type=2;
limits;

plot_conf(dfvmap);  oplot_conf(sfvmap, line);
connect_points(0); pointstyle(4); oplot(f3, v3, 5 ); connect_points(1);
pid=dup_plot("1103_12A_HEG_m1_conf_map_3.ps/vcps",p1); resize(6*2.54,1.0);
plot_conf(dfvmap);  oplot_conf(sfvmap,line);
connect_points(0); pointstyle(4); oplot(f3, v3, 5 ); connect_points(1);
close_plot(pid);

%-

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Fit_Verbose = 0;