%
%  This script loads in the simulated data, sets the pileup kernel,
%  defines the model, and then fits the data.  It steps through different
%  values of alpha, and writes the results to a file (which were then plotted
%  with idl).
%
%  This script is written in S-Lang.  More information about S-Lang can
%  be found at:  www.s-lang.org
%
%  As far as the s-lang scripts go, this is extremely basic.  One can
%  do much fancier things and then run them in ISIS.  As written here,
%  this is designed to be a simple executable file (the default for
%  scripting in Sherpa); however, within ISIS, any of this could be done
%  interactively during the session, allowing one to "program on the fly".
%  Also note below that within ISIS, I could capture many intermediate and
%  final variable values, and pass them back to ISIS for further use.
%  For the most part, that ability within ISIS is not being utilized here.
%
%  To run the script, assuming it has been saved as "pilesim.sl", in isis type:
%       isis> .load pilesim.sl
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  First, we define a few functions.  Normally I have these defined
%  upon start-up in ISIS.  I write them out here for completeness
%
%  ISIS defaults bins to wavelength.  The following allows me to
%  eXclusively NOTICE data bins by keV instead.
%

public define keV_note(a,b,c)
     {
     variable bi, ci ;
     variable conv=12.398424 ;
     if ( b != 0.) bi = conv/b; else return 0;
     if ( c !=0. ) ci = conv/c; else return 0;
     xnotice (a,ci,bi);
     return 1;
     }
%
% Raised on ftools and XSPEC, the following acts like a
% "group min minv" in grppha, starting the grouping at the
% lowest energy bin, and working towards higher energy (as opposed
% to doing it low to high in wavelength space)
%

public define grppha(a,minv)
     {
      rebin_data(a,0);
      variable s = get_data_counts(a);
      variable len = length(s.value);
      variable idx = [0:len-1:1];
      idx[*] = 1;
      variable sn = 1;
      variable ssum = 0.;
      variable i;
      for (i=len-1; i >=0; i--)
       {
        idx[i] = sn*idx[i];
        ssum += s.value[i];
        if (ssum >= minv)
           {
            sn = sn*(-1);
            ssum = 0.;
           }
        }
       rebin_data(a,idx);
       s = get_data_counts(a);
%
% Here I'm choosing the error bars to be square root of N.
% One can readily redefine this to provide any error bars you like.
%
       s.err = sqrt(s.value);
       put_data_counts(a,s);
%
% Here I'm just returning a "1" to indicate successful completion.
% However, I could have easily returned a structure giving information
% (new length, average signal to noise, whatever) about my newly
% rebinned data set.
%
       return 1;
      }
%
% Load the data file in
%

() = load_data("pilesim.pi");

%
% Group to a minimum of 20 counts per bin, and only notice ~0.5-3 keV
%

() = grppha(1,20);
() = keV_note(1,0.52,2.8);

%
% Set the pileup kernel, import XSPEC functions, and define the model
%

set_kernel(1,"pileup");
import("xspec");
fit_fun("phabs(1)*(bbody(1))");

%
% Normally, one can now define the starting values for the fits
% by either reading in a file, or editing them in a (vi, emacs, or
% choose your favorite editor) window.  Here, however, I'll define
% them one parameter at a time. Syntax is: parameter value,
% freeze (1) or thaw (0), min value, max value
%

% N_H
set_par(1,0.2,0,0.05,0.4);

% Bbody Norm
set_par(2,2.e-6,0,1.e-7,1.e-5);

% Bbody Temp
set_par(3,0.2,0,0.05,0.4);

% Pileup parameter alpha, which we freeze
set_par(6,1.,1);

%
% Next we define a function to loop through the values of alpha,
% fit the model, and then write the results to a file
%

public define alphasrch()
   {
   variable info=struct {statistic, num_variable_params, num_bins};
   variable alpha=[1.0:-0.01:-0.05];
   variable xi=@alpha;
   variable nh=@alpha;
   variable bbnorm=@alpha;
   variable kt=@alpha;
   variable i=0;
%
% Pileup fits can be fussy, so we try a couple of methods and
% loop through each a few times.  More clever ways of ensuring
% true minima can be programmed.
%
   foreach(alpha)
   {
% Me just being brain dead with roundoff errors above
      if (alpha[i] < 0.) alpha[i] = 0.;
      set_par(6,alpha[i],1);
      renorm_counts;
      set_fit_method("subplex;maxnfe=1000");
      loop(6)
          {
           fit_counts;
          }
      set_fit_method("marquardt;max_loops=100");
      loop(2)
         {
         fit_counts(&info);
         }
      xi[i]=info.statistic;
      nh[i]=get_par(1);
      bbnorm[i]=get_par(2);
      kt[i]=get_par(3);
       i++;
   }
   variable results=struct{nh, bbnorm, kt, alpha, xi};
   results.nh=@nh;
   results.bbnorm=@bbnorm;
   results.kt=@kt;
   results.alpha=@alpha;
   results.xi=@xi;
%
% Write out the results
%
   variable fp = fopen("pilesim.results","w");
   i=0;
   foreach(alpha)
      {
 ()=fprintf(fp,"%11.2e  %11.2e  %11.2e  %11.2e  %12.3e  \n",
                       alpha[i], nh[i], kt[i], bbnorm[i], xi[i]);
      i++;
      }
   () = fprintf(fp,"Alpha            Nh            kT          BB_norm      Xi \n");
   fclose(fp);
%
% One can return the results as a variable or structure, which can
% then be used in other S-Lang programs within ISIS
%
   return results;
   }

%
% Run the above, and just write the results to a file without doing
% anything further to them.
%

alphasrch();