%
% file: do_streak_fit.sl
% 
% Get and fit the streak events from an obsid. 10/17/08 - dd
% 
% Usage:  
%   - edit this file to select the desired obsid, then:
%   - hydra> .source do_streak_fit
%
% Note that the default findzo energy range of 0.5 to 4 keV is used.

% - - -
% Add get_streak_events to findzo by doing:
%    ! cat findzo.sl get_streak_events.sl > findzo_plus.sl
% where get_streak_events.sl is just the "define findzo ..." routine
% but end it prematurely with:
%  ...
%   % And draw a first region box to isolate transfer streak events  
%   % without zero order and with energy filtering. 
%   %
%      l=where(xp > (gx-boxw_strc/2) and xp < (gx+boxw_strc/2) and 
%              yp > (gy-boxl_strc/2) and yp < (gy+boxl_strc/2) and
%              x > 0.0            and y    > 0.0         and  
%              hypot(xp-gx,yp-gy)>rad_strc and (en > n_eng_lo and en < n_eng_hi));
%   
%      % offset w.r.t. zo for the returned event locations:
%      return xp[l]-gx, yp[l]-gy;
%   }
% - - -

.load findzo_plus
  
% Set the location of the tgcat dir:
variable tgcat_dir = "/data/tgcat";

variable xc, yc, xs, ys;

variable obsid_str, obs_to_do, iobs=0;

variable grat_str = "m";

% ----- Loop over a bunch...
%% obs_to_do=[702,1921,2741,2739,2708];
%%obs_to_do=[3814,3823,4572,4584,5040,6601,7449,7451];
%%obs_to_do=[1010,1103,1318,2583,3674,5955,6471];
%%obs_to_do=[337,1014,1705,3167,3706,3708,5173,7291];
%%obs_to_do=[105,106,1252,1451];
%%obs_to_do=[706,609,4552,7280];
%%obs_to_do=[1928,4574,15,425,3504,];
%%_for iobs (0,length(obs_to_do)-1,1)
%%{
%%obsid_str = string(obs_to_do[iobs]); 

% - - - - -  Or just do one
% Select the obsid from available ones:
%obsid_str="702";
%obsid_str="1921";
%obsid_str="2741";
%obsid_str="2739";
%obsid_str="2708";
%obsid_str="4420";  % very high counts
%obsid_str="3814";
%obsid_str="3823";
%%%obsid_str="4572";  % skip: double source ~ aligned by roll
%obsid_str="4584";
%obsid_str="5040";
%obsid_str="6601";
%obsid_str="7449";
%obsid_str="7451";
%obsid_str="7485";
%obsid_str="9638";
%obsid_str="9858";
%obsid_str="9076";
%obsid_str="9705";
%obsid_str="9712";
%obsid_str="9703";
%obsid_str="7511";
%obsid_str="8525";
%obsid_str="10659";
%obsid_str="10661";
%obsid_str="10695";
%obsid_str="10857";
obsid_str="10679";
% - - - - - 
% LETG-ACIS for fun:
%grat_str="l";
%obsid_str="8319";
% - - - - -

% clear out previous plots
close_plot; close_plot; close_plot;

% find the center with findzo using MEG (or LEG):
(xc,yc) = findzo(tgcat_dir+"/obs_"+obsid_str+"/evt2",grat_str);

% now get the rotated streak events:
(xs, ys) = get_streak_events(tgcat_dir+"/obs_"+obsid_str+"/evt2",grat_str,4,xc,yc);

% plot them
window(1);
xrange;yrange;
connect_points(0);
plot(xs,ys);
connect_points(1);

% Load the real MEG-1 data:
delete_data(all_data);
if (grat_str == "l")
{
   ()=load_data(tgcat_dir+"/obs_"+obsid_str+"/pha2",4);   
}
else 
{
   ()=load_data(tgcat_dir+"/obs_"+obsid_str+"/pha2",9);   
}

%% plot_data_counts(1);
list_data;


% Get the data counts strucuture: (bin_lo, bin_hi, value, err)
variable dstruct;
dstruct = get_data_counts(1);

% create a fake counts array based on the streak values converted
% to wavelength, centered at 15.0 A.
% The units are sky pixels (0.492" each)
variable streak_pha;
streak_pha = hist1d (xs+15.0, dstruct.bin_lo);
put_data_counts (1, dstruct.bin_lo, dstruct.bin_hi, streak_pha, sqrt(streak_pha));

% Now group, notice, and fit...
if (grat_str == "l")
{
   % LEG
   group_data(1,20);
   xnotice(1,12.0,18.0);   
}
else 
{
   % MEG pha has binning of 0.005 A
   % so group x40 for 0.2 sky-pixel bins
   group_data(1,40);
   % limit the fit to +/- 2.0 pixels:
   xnotice(1,13.0,17.0);   
}

xrange(10.0,20.0);
xlin;
ylog;
()=open_plot("/XWIN");
plot_data_counts(1);

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

% guess the number of event per exposure time:
set_par("gauss(1).area",0.5*length(xs)/get_data_exposure(1), 0);
set_par("gauss(1).center",15.0, 0, 14.0,16.0);
set_par("gauss(1).sigma",(38.0/24.0)/2.35, 0);
% background ~ constant with exposure
set_par("poly(1).a0",0.0005, 0);
set_par("poly(1).a1",0.0, 1);
set_par("poly(1).a2",0.0, 1);

set_fit_method("lmdif");

%% ()=eval_counts;
%% rplot_counts(1);
()=fit_counts;

message("  Parameters for the Streak fit:");
message("");
list_par;

% Get the 90% confidence limits:
variable conflo, confhi;
(conflo,confhi)=conf("gauss(1).sigma",1,0.001);

message("    Obsid "+obsid_str+" :");
message("");
message("  Streak area (counts): "+string(get_par(1)*get_data_exposure(1)));
message("");
message("  FWHM of the Streak fit = " + string(24.0*(2.35*get_par(3))) + " um" );
message("   ( "+string(24.0*(2.35*conflo)) + 
" to "+string(24.0*(2.35*confhi))+", 90% confidence;  +/- "+
string(0.5*24.0*(2.35*(confhi-conflo)))+" )");
message("");

% Save it to .gif
variable GIFwin=open_plot("streak_"+obsid_str+"_plot.gif/GIF");
variable titstr="Obsid "+obsid_str+" :  FWHM of Streak = " + string(24.0*(2.35*get_par(3))) + " um" +
" +/- "+string(0.5*24.0*(2.35*(confhi-conflo)))+" (90%)";
title(titstr);
rplot_counts(1);
close_plot(GIFwin);

title(titstr);
rplot_counts(1);


% Loop over a bunch end
%%}