% rg_fit.i
%
% Use forward fold fitting to fit the data in a region.
% The fit function is constructed as the sum of
% wabs(1)*gauss(n) for each bump in the region and
% a constant continuum, wabs(1)*poly(1).
%
%-------------------------------------------------------
% The forward folding here is the fit_counts (Assigned_ARFRMF); case
% of fit options:
%   fit_counts;
%   fit_counts (Assigned_ARFRMF);    % or Assigned_RMFARF
%   fit_counts (Ideal_RMF);
%   fit_counts (Ideal_ARF);
%   fit_counts (Ideal_ARFRMF);       % or Ideal_RMFARF
%   fit_flux (Assigned_RMF);
%   fit_flux (Assigned_ARFRMF);      % ideal ARF is applied anyway
%-------------------------------------------------------

% new variables needed:
variable fun_string, ib, irg, lam_fit_range; 
variable yrangemax, xrangemin, xrangemax;

% assume fitting the counts data w/arf,rmf available.
%%use_counts;

%  index into rgl:
irg=0; while (rgl[irg].id != rgsel) {
  irg=irg+1;}

message("");
message(" --------------------------- rg_fit ---------------------------");
message("");
message("   Data set: "+data_set_name);
message("");
message("   Fit in "+ rgl[irg].desc + " Region ["+string(rgsel)+"]   (z="+string(zpg)+")" );
message("");

lamlow = rgl[irg].lo;
lamhigh = rgl[irg].hi;

% List the not-ignored bumps in this region...
%
%%bpsel = where( (bpl.lam < lamhigh) and (bpl.lam > lamlow)
%%	and (bpl.ignore == 0) );
%
bpsel = where (((1.0+zpg)*array_struct_field (bpl, "lam") < lamhigh)
                 and ((1.0+zpg)*array_struct_field (bpl, "lam") > lamlow)
                 and (array_struct_field (bpl, "ignore") == 0));

message("   Bumps included :");
ib=0; while (ib < length(bpsel)) {
  message("     [" + string(bpl[bpsel[ib]].id) + "]  " +
		bpl[bpsel[ib]].desc +"  " + string(bpl[bpsel[ib]].lam) );
  ib = ib + 1;
  }

message("   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -");

%
%------------- Setup and Show the data --------------------
%

%--------------- ignore all --------------------
ignore([1:12]);

% guess at total counts in line
% notice the range
notice(megm, lamlow, lamhigh);
notice(megp, lamlow, lamhigh);
if (data_n_spectra > 2) {
  notice(hegm, lamlow, lamhigh);
  notice(hegp, lamlow, lamhigh);
  }
% and ignore some:
() = evalfile(data_dir + "/" + data_set_name+"_ignores.i");

xrangemin = lamlow;
xrangemax = lamhigh;

% Y range
ylin;
gds = get_data_counts (megm);
yrangemax = max(gds.value[where(gds.bin_lo > xrangemin and
	gds.bin_hi < xrangemax)]);
gds = get_data_counts (megp);  
yrangemax = max( [ max(gds.value[where(gds.bin_lo > xrangemin and
	gds.bin_hi < xrangemax)]),
	yrangemax] );
if (data_n_spectra > 2) {
  gds = get_data_counts (hegm);  
  if (length(where(gds.bin_lo > xrangemin and gds.bin_hi < xrangemax)) > 0) {
    yrangemax = max( [ max(gds.value[where(gds.bin_lo > xrangemin and
	gds.bin_hi < xrangemax)]),
	yrangemax] );
    }
  gds = get_data_counts (hegp);  
  if (length(where(gds.bin_lo > xrangemin and gds.bin_hi < xrangemax)) > 0) {
    yrangemax = max( [ max(gds.value[where(gds.bin_lo > xrangemin and
	gds.bin_hi < xrangemax)]),
	yrangemax] );
    }
  }

% and measure the flux there...
%%use_flux;
rsout = region_flux(megm, lamlow, lamhigh);
approx_flux = rsout.sum;
rsout = region_flux(megp, lamlow, lamhigh);
approx_flux = (rsout.sum + approx_flux)/2.0;
%%use_counts;
message("");
message(" MEG +/- Ave flux = "+string(approx_flux));
% Case if working at high-energy with HEG spectra
if (data_n_spectra > 2  and  lamhigh < 3.5) {
  %%use_flux;
  rsout = region_flux(hegm, lamlow, lamhigh);
  approx_flux = rsout.sum;
  rsout = region_flux(hegp, lamlow, lamhigh);
  approx_flux = (rsout.sum + approx_flux)/2.0;
  %%use_counts;
  message(" HEG +/- Ave flux = "+string(approx_flux));
  }
message("");

% Scale this flux by the number of bumps in the region
approx_flux = approx_flux/(length(bpsel)+1);

% plot window to use
% - - -
if (psplot == 1) {
  % Postscript plot:
  vps = open_plot(data_set_name+"_rg"+ string(rgsel) +"fit.ps/vcps", 1, 2);
  window(vps);
  cmegm = orange;
  if (data_n_spectra < 4) {cmegm = green;}
  } else {
% - - - OR - - -
  window(winFit);
  erase;
  }
% - - -


plot_unit("Angstrom");
title(data_set_name + ": Fitting in " + rgl[irg].desc + " Region ["+string(rgsel)+"]  " +
	 "(z="+string(zpg)+")");

% Plot the data
% setup the plot 
if (psplot == 1) {
  pane(1);
  } else {
  pane(2);
  erase;
  pane(1);
  erase;
  }

plot_bin_integral;
charsize(1.0);
xlin;
xrange(lamlow, lamhigh);
ylin;
yrange(0.,1.1*yrangemax);
% want error bars on, or not.
errorbars(1);
plot_data_counts(megm,white);
oplot_data_counts(megp,cmegp);
oplot_data_counts(megm,cmegm);
if (data_n_spectra > 2) {
  oplot_data_counts(hegp,chegp);
  oplot_data_counts(hegm,chegm);
  }

% and the groups
() = evalfile("plot_bpl_groups.i");

%
%------------- Now do some fitting  ---------------
%

% Define the function to use:
% One or more "wabs(1)*gauss("+string(ib)+")" plus a final "wabs(1)*poly(1)"
fun_string = "wabs(1)*gauss(1)";
ib=1; while (ib < length(bpsel)) {
  fun_string = fun_string + "+wabs(1)*gauss("+string(ib+1)+")";
  ib=ib+1;
  }
% add poly:
fun_string = fun_string + "+wabs(1)*poly(1)";
% and set it up:
fit_fun(fun_string);

% WABS
set_par(1, data_nh, 1, data_nh, data_nh);
%

% Loop over the "bumps" and set the initial Gaussian parameters:
%  
% Constrain the bump centers to within +/-0.5*( lam_range / (nbumps + 1) )
lam_fit_range = 0.5*((lamhigh-lamlow)/(length(bpsel)+1));
ib=0; while (ib < length(bpsel)) {
  % GAUSSIAN
  % AREA - guess rough number of counts
  set_par(2+3*ib, approx_flux, 0, 0.001*approx_flux, 100.0*approx_flux);
  % CENTER - easy to guess
  lamcent = (1.0+zpg)*bpl[bpsel[ib]].lam;
  set_par(3+3*ib, lamcent, 1, lamcent-lam_fit_range, lamcent+lam_fit_range);
  % SIGMA - use the average between HEG and MEG point-src resolutions
  set_par(4+3*ib, 0.5*(0.012+0.023)/2.35, 1, 0.1*0.004, 3.0*0.021);
  %
  ib = ib + 1;
  }

% POLY
%  Setting continuum to zero:
%%set_par(2+3*length(bpsel), 0.0*approx_flux, 1, 
%%	0.0*approx_flux, 100.0*approx_flux);
%
set_par(2+3*length(bpsel), 0.1*approx_flux, 0, 
	0.0001*approx_flux, 100.0*approx_flux);
set_par(3+3*length(bpsel), 0.0, 1, 0.0,0.0);
set_par(4+3*length(bpsel), 0.0, 1, 0.0,0.0);

% Manual adjustments if desired...
%%  set_par(3, 19.07, 1, lamcent-lam_fit_range, lamcent+lam_fit_range);
%%  set_par(6, 19.02, 1, lamcent-lam_fit_range, lamcent+lam_fit_range);
%%  set_par(4, 0.5*(0.012+0.023)/2.35, 1, 0.004, 3.0*0.021);

% list the function and parameters before the fit
message("");
list_par;

% Do the fit
% and get the returned numbers
% in start.i:
% variable firet = struct {statistic, num_variable_params, num_bins};

%%fit_counts;
fit_counts(Assigned_ARFRMF, &firet);

list_par;

oplot_model_counts(megm, purple);
oplot_model_counts(megp, blue);
if (data_n_spectra > 2) {
  oplot_model_counts(hegm, purple);
  oplot_model_counts(hegp, blue);
  }

% Show the continuum level...
fit_fun("wabs(1)*poly(1)");
eval_counts;
oplot_model_counts(megm, purple);
oplot_model_counts(megp, blue);
if (data_n_spectra > 2) {
  oplot_model_counts(hegm, purple);
  oplot_model_counts(hegp, blue);
  }
fit_fun(fun_string);
eval_counts;

% Save some info to the bump structure...
ib=0; while (ib < length(bpsel)) {
  % fit wavelength
  bpl[bpsel[ib]].lfit = get_par(3+3*ib);
  % fit area
  bpl[bpsel[ib]].afit = get_par(2+3*ib);
  % fit width
  bpl[bpsel[ib]].wfit = 2.35*get_par(4+3*ib);
  % these next ones are the same for all bumps in the fit:
  %   continuum level:
  bpl[bpsel[ib]].cfit = get_par(2+3*length(bpsel));
  %   chi-squared:
  bpl[bpsel[ib]].chifit = firet.statistic /
	(firet.num_bins - firet.num_variable_params);
  %   string to describe the data and region fit:
  bpl[bpsel[ib]].rgfit = data_set_name+": "+rgl[irg].desc;
  ib = ib + 1;
  }
% Now the more time consuming confidence values...
% which can be skipped...
%
% Area
if ( 0 > 0) {
  message("");
  ib=0; while (ib < length(bpsel)) {
    message("Area Gauss("+string(ib+1)+") 1-sigma confidence range ...");
    (conflo, confhi) = conf (2+3*ib, 0, 0.10);
    message("   " +  string(conflo) + " to "+string(confhi) +
	"  or +/- " + string(0.5*(confhi-conflo)) );
    % fit area confidence levels
    bpl[bpsel[ib]].aclo = conflo;
    bpl[bpsel[ib]].achi = confhi;
    ib = ib + 1;
    }
  message("");
  }
% Center
if ( 0 > 0) {
  message("");
  ib=0; while (ib < length(bpsel)) {
    message("Center Gauss("+string(ib+1)+") 1-sigma confidence range ...");
    (conflo, confhi) = conf (3+3*ib, 0, 0.10);
    message("   " +  string(conflo) + " to "+string(confhi) +
	"  or +/- " + string(0.5*(confhi-conflo)) );
    % fit area confidence levels
    bpl[bpsel[ib]].lclo = conflo;
    bpl[bpsel[ib]].lchi = confhi;
    ib = ib + 1;
    }
  message("");
  }
% FWHM
if ( 0 > 0) {
  message("");
  ib=0; while (ib < length(bpsel)) {
    message("FWHM Gauss("+string(ib+1)+") 1-sigma confidence range ...");
    (conflo, confhi) = conf (4+3*ib, 0, 0.10);
    message("   " +  string(2.35*conflo) + " to "+string(2.35*confhi) +
	"  or +/- " + string(2.35*0.5*(confhi-conflo)) );
    % fit FWHM confidence levels
    bpl[bpsel[ib]].wclo = 2.35*conflo;
    bpl[bpsel[ib]].wchi = 2.35*confhi;
    ib = ib + 1;
    }
  message("");
  }
% Continuum
if ( 0 > 0) {
  message("");
    message("Continuum, poly("+string(1)+") 1-sigma confidence range ...");
    (conflo, confhi) = conf (5+3*(length(bpsel)-1), 0, 0.10);
    message("   " +  string(conflo) + " to "+string(confhi) +
	"  or +/- " + string(0.5*(confhi-conflo)) );
    % fit Continuum confidence levels
  ib=0; while (ib < length(bpsel)) {
    bpl[bpsel[ib]].cclo = conflo;
    bpl[bpsel[ib]].cchi = confhi;
    ib = ib + 1;
    }
  message("");
  }

% Plot the residuals!
pane(2);
if (psplot != 1) {
  erase;
  }

errorbars (0);

label ("Wavelength [Angstrom]", "Residuals", "");

gds = get_data_counts (megm);
gdmodel = get_model_counts (megm);

xlin;
xrange(lamlow, lamhigh);
ylin;
% Set yrange of residuals to be the same total
% range as above plot - hence error bars are same size.
% So, use +/- 1.1*yrangemax/2.0 
yrange(-1.1*yrangemax/2.0, 1.1*yrangemax/2.0);
hplot (gds.bin_lo, gds.bin_hi, gds.value - gdmodel.value, cmegm);
oplot([lamlow,lamhigh],[0.0,0.0],white);
ohplot (gds.bin_lo, gds.bin_hi, gds.value - gdmodel.value, cmegm);

gds = get_data_counts (megp);
gdmodel = get_model_counts (megp);
ohplot (gds.bin_lo, gds.bin_hi, gds.value - gdmodel.value, cmegp);

if (data_n_spectra > 2) {
  gds = get_data_counts (hegm);
  gdmodel = get_model_counts (hegm);
  ohplot (gds.bin_lo, gds.bin_hi, gds.value - gdmodel.value, chegm);
  gds = get_data_counts (hegp);
  gdmodel = get_model_counts (hegp);
  ohplot (gds.bin_lo, gds.bin_hi, gds.value - gdmodel.value, chegp);
  }


% when done...
%--------------- notice all --------------------
notice(megm);
notice(megp);
if (data_n_spectra > 2) {
  notice(hegm);
  notice(hegp);
  }
% and ignore some:
() = evalfile(data_dir+"/"+data_set_name+"_ignores.i");

% leave the error bars off
errorbars(0);

% - - - if ps output - - -
if ( psplot == 1) {
  close_plot(vps);
  cmegm = yellow;
  }
% - - -