% fit_line.i
%
% do .load fit_setup.i once before running this

% some new variables
variable approxcounts, fxtitle, data_yunit;
use_counts;
data_yunit = "counts";

%--------------- ignore all --------------------
ignore(megm);
ignore(megp);
ignore(hegm);
ignore(hegp);
% if simulation present:
ignore(simm);
ignore(simp);
% and default labels off:
LABEL_BY_DEFAULT=0;

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

% Select the data and wavelength to fit
% My favorite lines:
%  18.972   ( 0.6535 keV)
%  15.012   ( 0.826 keV)
%  12.1315  ( 1.022 keV)
%   9.169   ( 1.352 keV)  Mg
%
% Data set selection:
%%fdata = hegm;
%%cfdata = chegm;
fperiod = p_heg ;
%
% Wavelength selection:
lamcent = 15.0;

message("");
message(" --------------------------- fit_line ---------------------------");
message("");
message("   Data set: "+data_set_name+"; hist_index = "+string(fdata));
message("");
message("   Fit around: "+string(lamcent)+" A   ("+string(hc/lamcent)+" keV)");
message("");

% parameters of the fitting
lamrange = 0.75*2.0*0.3;
corefwhm = 2.0*0.0180;
wingcont = 0.0;

% the lamrange and corefwhm above are set for MEG spectra,
% scale them to other grating types:
range_scale = fperiod/p_meg;
lamrange = lamrange * range_scale;
corefwhm = corefwhm * range_scale;

% calculate the limits...
lamlow = lamcent - lamrange/2. ;
lamhigh = lamcent + lamrange/2. ;
corelow = lamcent -  2.0*corefwhm;
corehigh = lamcent +  2.0*corefwhm;

% guess at total counts in line
% notice the core range
notice(fdata, corelow, corehigh);

% and see how many counts are there...
rsout = region_sum(fdata, corelow, corehigh);
approxcounts = rsout.sum;

% plot window to use
window(winFit);
pane(1);
erase;
pane(2);
erase;

%plot_unit("Angstrom");

ftitle = "Fit to CORE of "+string(lamcent)+" A LRF,  "+data_set_name;
fxtitle = data_xtitle;

% Plot the data
% setup the plot 
pane(1,1);
charsize(1.2);
plot_bin_integral;
xrange(corelow, corehigh);
% Use log scale, one size fits all!
%ylog;
%yrange(-0.5, 5.) ; 
ylin;
plot_data(fdata);
yrange; 
% want error bars on
errorbars(1);
erase;
plot_data(fdata);
color(white);
label(fxtitle, data_yunit+"/bin", ftitle);
% and the groups
%%() = evalfile("plot_all_groups.i");

%
%------------- Now do some fitting to the CORE ---------------
%

% ignore all of this data set
ignore(fdata);

% Notice the region we want to fit
% Fit the core only
% limited core
notice(fdata, lamcent-0.7*corefwhm, lamcent+0.7*corefwhm);
% the full core
% notice(fdata, corelow, corehigh);

% Define the function to use
fit_fun("gauss(1)+poly(1)");

% GAUSSIAN
% Guess the parameters:
% AREA - guess rough number of counts
set_par(1, approxcounts, 0, 0.1*approxcounts, 10.0*approxcounts);
% CENTER - easy to guess
set_par(2, lamcent, 0, lamlow, lamhigh);
% FWHM or SIGMA - use corefwhm as guess
set_par(3, corefwhm/2.35, 0, 0.2*corefwhm/2.35, 5.*corefwhm/2.35);

% POLY
% use constant part only
% Fix it at 0 as default
set_par(4, 0.0*approxcounts, 1, 0.0, 10.0*approxcounts);
%%set_par(4, approxcounts, 0.001*approxcounts, 100.0*approxcounts, 0);
set_par(5, 0.0, 1, 0.0,0.0);
set_par(6, 0.0, 1, 0.0,0.0);


% list the parameters before the fit
%%list_par;
%%eval_counts(Ideal_ARFRMF);
%%oplot_model(fdata);

% Do the CORE fit
fit_counts(Ideal_ARFRMF);

% find the error on the gaussian width
% 0 is for 1-sigma
(conflo,confhi) = conf(3,1,0.05);
%%fit_counts;
%%(conflo,confhi) = conf(3,1,0.05);

% Show the best fit parameters
message("");
message("  ROI CORE counts: "+string(approxcounts));
message("");
message("  Parameters for the CORE fit:");
message("");
list_par;
message("");
%%message("  E/dE for the CORE fit = "+string(get_par(2)/(2.35*get_par(3))) );
%%message("   ( "+string(get_par(2)/(2.35*confhi)) + 
%%	" to "+string(get_par(2)/(2.35*conflo))+", 90% confidence )");
% this message added for use in zo_1d_analysis:
message("  FWHM of the CORE fit = " + 
	string(24.0*(2.35*get_par(3))/use_angpix) + " um" );
message("   ( "+string(24.0*(2.35*conflo)/use_angpix) + 
	" to "+string(24.0*(2.35*confhi)/use_angpix)+", 90% confidence )");
message("");

% add the model
% undo the data color - arg!
unset_data_color(fdata);
oplot_model(fdata, cmodel);
set_data_color(fdata, cfdata);

% save the core's fit center for use below
variable core_fit_center;
core_fit_center = get_par(2);

%
%----------- Now do some fitting to the WINGS -------------
%

% setup the plot 
window(winFit);
pane(1,2);
charsize(1.2);
plot_bin_integral;
xrange(lamlow,lamhigh);
% Use log scale, one size fits all!
ylog;
yrange(0.1^0.5, 10.0^5.) ; 

ftitle = "Fit to WINGS of LRF,  "+data_set_name;
fxtitle = data_xtitle;

% Plot the data
notice(fdata, lamlow,lamhigh);
% want error bars on
errorbars(1);
erase;
plot_data(fdata);
color(white);
label(fxtitle, "Log10( "+data_yunit+"/bin )", ftitle);
% and the groups
%%() = evalfile("plot_all_groups.i");

% ignore all of this data set
ignore(fdata);

% Notice the two wing regions... and get total counts in them
notice(fdata, lamlow, lamcent-2.0*corefwhm);
rsout = region_sum(fdata, lamlow, lamcent-2.0*corefwhm);
approxcounts = rsout.sum;
appwingcounts = approxcounts;
notice(fdata, lamcent+2.0*corefwhm, lamhigh);
rsout = region_sum(fdata, lamcent+2.0*corefwhm, lamhigh);
approxcounts = rsout.sum;
appwingcounts = appwingcounts + approxcounts;

fit_fun("lorentz(1)");

% Guess the parameters:
% AREA - guess fraction of number of wing counts
set_par(1, 3.0*appwingcounts, 0, 0.3*appwingcounts, 30.0*appwingcounts);
% CENTER - freeze at CORE fit value
set_par(2, core_fit_center, 1, lamlow, lamhigh);
% FWHM or SIGMA - use corefwhm as guess
set_par(3, corefwhm, 0, 0.2*corefwhm, 5.*corefwhm);
% list the parameters before the fit
%%list_par;

% Do the WINGS fit
fit_counts;
% if fit doesn't work do this:
% eval_fun;

% Show the best fit parameters
message("");
message("  ROI WING counts: "+string(appwingcounts));
message("");
message("  Parameters for the WINGS fit:");
message("");
list_par;
message("");

% and the model
unset_data_color(fdata);
oplot_model(fdata, cmodel);
set_data_color(fdata,cfdata);

% when done...
%--------------- notice all --------------------
notice(1);
notice(2);
notice(3);
notice(4);
% and labels on
LABEL_BY_DEFAULT=0;

% show the emission model FYI
pane(1,2);
%%ohplot(mod_binlo, mod_binhi, mod_em/1.E42, yellow);

% leave the error bars off
errorbars(0);