% bp_fit_counts.i
%
% Fit directly to the counts histogram without
% using arf or rmf.
%
% - - - 
% %fdata = megp:
% ibb=5; while(ibb < 43) {bpsel=ibb; ibb=ibb+1;
%        ()=evalfile("bp_fit_counts.i"); }
% megp_counts = array_struct_field(bpl,"afit");
% - - - 
% new variables needed:
variable fun_string, ib, lam_fit_range; 
variable found_better_fit=0;
variable how_many_times = 0;
variable c0,c1,c2,c3;

%----------- Select the histogram to fit ------------
% values set outside of this routine.
%%fdata = megm;
%%cfdata = cmegm;

% Find the bump index...
%  index into bpl:
ii=0; while (bpl[ii].id != bpsel) {
  ii=ii+1;}

message("");
message(" --------------------------- bp_fit_counts_fixfwhm ----------------");
message("");
message("   Data set: "+data_set_name+" [index="+string(fdata)+"]");
message("");
message("   Fit the "+ bpl[ii].desc + " Bump Counts ["+string(bpsel)+"]   (z="+string(zpg)+")" );
message("");

lamlow = (1.0+zpg)*bpl[ii].lam - 0.5*rgbprange;
lamhigh = (1.0+zpg)*bpl[ii].lam + 0.5*rgbprange;

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

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

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

% guess at total counts in line
% notice the range
notice(fdata, lamlow, lamhigh);
% measure the counts there...
rsout = region_counts(fdata, lamlow, lamhigh);
approx_flux = rsout.sum;
message("");
message(" ROI counts = "+string(approx_flux));
message("");

% plot window to use
% - - -
if (psplot == 1) {
  % Postscript plot:
  vps = open_plot(data_set_name+"_i"+string(fdata)+"_bp" +
	 string(bpsel) +"_fixfwhm.ps/vcps", 1, 2);
  window(vps);
  pane(1);
  % replace yellow with orange, if there:
  if (cfdata == yellow) { cfdata = orange;}
  } else {
% - - - OR - - -
  window(winFit);
  pane(2);
  erase;
  pane(1);
  erase;
  }
% - - -

% Move plotting to after the fitting...

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

% Define the function to use:
% One or more "gauss("+string(ib)+")" plus a final "poly(1)"
fun_string = "gauss(1)";

% add poly:
fun_string = fun_string + "+poly(1)";
% and set it up:
fit_fun(fun_string);

% Set the initial Gaussian parameters:
%  
% Wide-ish range...
lam_fit_range = (lamhigh-lamlow)/2.0;
ib=0;
  % GAUSSIAN
  % AREA - guess rough number of counts
  set_par(1+3*ib, approx_flux, 0, 0.001*approx_flux, 100.0*approx_flux);
  % CENTER - easy to guess
  lamcent = (1.0+zpg)*bpl[ii].lam;
  set_par(2+3*ib, lamcent, 0, lamcent-lam_fit_range, lamcent+lam_fit_range);
  % SIGMA - !!!
  % Use the value from the polynomial fit to MARX simulation
  % from http://space.mit.edu/HETG/technotes/ede_0102/ede_0102.html :
  if (fdata == megm  or  fdata == megp) {
    c0 = 0.018851868;
    c1 = -0.00033946575;
    c2 = 2.9781558e-05;
    c3 = -4.1097023e-07;
  }  else  {
    c0 = 0.010019662;
    c1 = -0.00014101982;
    c2 = 1.4591413e-05;
    c3 = -1.9632994e-07;
  }
  set_par(3+3*ib, (c0 + c1*lamcent + c2*lamcent^2 + c3*lamcent^3)/2.35, 
	1, 0.002, 0.7*rgbprange/2.35);

% POLY
set_par(1+3*1, 0.1*approx_flux, 0, 
	0.0001*approx_flux, 100.0*approx_flux);
set_par(2+3*1, 0.0, 1, 0.0,0.0);
set_par(3+3*1, 0.0, 1, 0.0,0.0);

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

% flag to do the fitting...
found_better_fit = 1;
how_many_times = 0;

while (found_better_fit > 0) {
% - - - - top of fitting loop - - - -

% Do the fit
% with no responses used:
fit_counts(Ideal_ARFRMF,&firet);
list_par;

how_many_times = how_many_times + 1;

% show the fit
% blue for preliminary...
%%oplot_model_counts(fdata, blue);

% Save the parameters
ib=0;
  % fit area
  bpl[ii].afit = get_par(1+3*ib);
  % fit wavelength
  bpl[ii].lfit = get_par(2+3*ib);
  % fit FWHM
  bpl[ii].wfit = 2.35*get_par(3+3*ib);
  % identify this fitting...
  bpl[ii].rgfit = data_set_name+" [index="+string(fdata)+"]: bp_fit_counts";
  % Chi-squared of fit
  bpl[ii].chifit = firet.statistic / 
	(firet.num_bins - firet.num_variable_params);
  % Continuum value
  bpl[ii].cfit = get_par(4);

% OK, flag done fitting unless conf's give better value...
found_better_fit = 0;

% re-do it if the fit wavelength is too far from the expected value
% (given as the average of lamlow and lamhigh ...)
% in order to better center the data around the bump...
if ( abs(bpl[ii].lfit - (lamlow+lamhigh)/2.0) > 0.05*rgbprange )
	{found_better_fit = 1;}

% Now the more time consuming confidence values...
% which can be skipped...
if ( 1 > 0) {
  message("");
  ib=0; while (ib < 1) {
    %
    message("Area of Gauss("+string(ib+1)+") 1-sigma confidence range ...");
    (conflo, confhi) = conf (1+3*ib, 0, 0.05);
    message("   " +  string(conflo) + " to "+string(confhi) +
	"  or +/- " + string(0.5*(confhi-conflo)) );
    % fit area confidence levels
    bpl[ii].aclo = conflo;
    bpl[ii].achi = confhi;
    if (conflo == confhi) {found_better_fit = 1;}
    %
    message("Lambda of Gauss("+string(ib+1)+") 1-sigma confidence range ...");
    (conflo, confhi) = conf (2+3*ib, 0, 0.05);
    message("   " +  string(conflo) + " to "+string(confhi) +
	"  or +/- " + string(0.5*(confhi-conflo)) );
    % fit area confidence levels
    bpl[ii].lclo = conflo;
    bpl[ii].lchi = confhi;
    if (conflo == confhi) {found_better_fit = 1;}
    %
%    message("FWHM of Gauss("+string(ib+1)+") 1-sigma confidence range ...");
%    (conflo, confhi) = conf (3+3*ib, 0, 0.05);
%    message("   " +  string(2.35*conflo) + " to "+string(2.35*confhi) +
%	"  or +/- " + string(2.35*0.5*(confhi-conflo)) );
%    % fit area confidence levels
    bpl[ii].wclo = bpl[ii].wfit+0.0001;
    bpl[ii].wchi = bpl[ii].wfit-0.0001;
%    if (conflo == confhi) {found_better_fit = 1;}
    %
    message("Continuum, poly("+string(1)+") 1-sigma confidence range ...");
    (conflo, confhi) = conf (4, 0, 0.05);
    message("   " +  string(conflo) + " to "+string(confhi) +
	"  or +/- " + string(0.5*(confhi-conflo)) );
    % continuum confidence levels
    bpl[ii].cclo = conflo;
    bpl[ii].cchi = confhi;
    if (conflo == confhi) {found_better_fit = 1;}
    %
    ib = ib + 1;
    }
  }

% if too many times just get out of here
if (how_many_times > 5) {
  found_better_fit = 0;
  message(" * ");
  message(" * Too many times - stop this!");
  message(" * ");
}
if (found_better_fit > 0) {
  message(" * ");
  message(" * Oh, man... Going to refit...");
  message(" * ");
  % re-select the data around current wavelength...
  lamlow = bpl[ii].lfit - 0.5*rgbprange;
  lamhigh = bpl[ii].lfit + 0.5*rgbprange;
  xnotice(fdata, lamlow, lamhigh);
}

% - - - - -  end of fitting loop - - - - 
}

% Plot the data
plot_unit("Angstrom");

title(data_set_name+" [index="+string(fdata)+"] : Fitting the " + bpl[ii].desc + " Bump ["+string(bpsel)+"]  (z="+string(zpg)+")" );


plot_bin_integral;
charsize(1.0);
xlin;
xrange(lamlow, lamhigh);
ylin;
yrange(0.,);
% want error bars on, or not.
errorbars(1);
plot_data_counts(fdata,white);
oplot_data_counts(fdata,cfdata);

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

% now plot the model in the color I like...
oplot_model_counts(fdata, purple);

% Plot the residuals!
pane(2);

errorbars (0);

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

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

% This region...
gdsel = where(gds.bin_lo > lamlow and gds.bin_hi < lamhigh);

xlin;
xrange(lamlow, lamhigh);
ylin;
% Set the y range to mimic the (upper) data plot in scale:
yrange( -0.5*max(abs(gds.value[gdsel])), 
	0.5*max(abs(gds.value[gdsel])) );
hplot (gds.bin_lo, gds.bin_hi, gds.value - gdmodel.value, cfdata);
oplot([lamlow,lamhigh],[0.0,0.0],white);
ohplot (gds.bin_lo, gds.bin_hi, gds.value - gdmodel.value, cfdata);

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

% when done...
%--------------- notice all --------------------
notice(fdata);

% leave the error bars off
errorbars(0);