%
% To run this script, save this file as 'alphasrch_sherpa.sl' in
% the directory that you plan to run it. (Or choose your own favorite
% filename.)  Start up sherpa, then type at the sherpa prompt:
%
%    sherpa> () = evalfile("alphasrch_sherpa.sl")
%
% Followed by (assuming the spectrum is in "pilesim.pi", and you want
% to fit that spectrum between 0.5 and 3 keV)
%
%    sherpa> results = alphasrch("pilesim.pi",0.5,3.0)
%
% "results" will be a structure (with the results, of course), which
% will also be written to the file: 'alphasrch_pilesim.pi_sherpa.dat'
%
% Note:  We are *not* subtracting any background in this fit!
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This reads the data, sets the model, then fits with several methods
%
define alphasrch(filename,lowe,highe)
{
   () = load_dataset(1,filename);
   () = set_ignore(1,,lowe);
   () = set_ignore(1,highe,);
%
% Below, the first are the actual model names, the second is what
% names they are assigned to for the purposes of this fit
%
   variable mods = ["xsphabs","xspowerlaw","xsbbody"];
   variable mod_names = ["phabs1","pl1","bb1"];
%
   variable i=0;
   loop(length(mods))
   {
     ()= create_model(mods[i],mod_names[i]);
     i++;
   }
   () = set_source_expr(1,"phabs1*(bb1+pl1)");
%
% Initialize model parameters, including step sizes for searches
%
   variable parms = ["phabs1.1","pl1.1","pl1.2","bb1.1","bb1.2"];
   variable mins =  [0.0,        0.0,    0.0,    0.1,    0.];
   variable maxs =  [2.0,        4.0,    1.e-4,  0.3,    1.e-4];
   variable vals =  [0.1,        2.0,    1.e-13, 0.192,  1.0e-5];
   variable steps = [0.01,       0.1,    3.e-14, 0.003,  1.e-7];
%
   i=0;
   loop(length(parms))
   {
      () = set_par(parms[i],"min",mins[i]);
      () = set_par(parms[i],"max",maxs[i]);
      () = set_par(parms[i],"value",vals[i]);
      () = set_par(parms[i],"delta",steps[i]);
      i++;
   }
   () = set_frozen("pl1.1");
%
   () = sherpa_eval("statistic chi dvar");
%
% Pileup is a change to the "fitting kernel", so it isn't handled
% exactly the same as other models (i.e., create_model alone doesn't
% invoke it, hence the sherpa_eval below)
%
   () = create_model("jdpileup","jdp");
   () = sherpa_eval("pileup=jdp");
%
% The pileup parameters, however, can be set like other sherpa models.
%
   () = set_par("jdp.5","value",3.2);
%
   variable stats, gp;
   variable r = struct{ xi, nh, bbnorm, kt, plnorm, alpha };
   r.alpha = [100:0:-5];
   r.alpha = r.alpha/100.;
   r.xi = @r.alpha;
   r.nh = @r.alpha;
   r.bbnorm = @r.alpha;
   r.kt = @r.alpha;
   r.plnorm = @r.alpha;
%
% We try several methods a few times each, since pile-up fits
% can be fussier than many models (due to many local minima
% in chi^2 space)
%
   variable method="method "+["simplex","levenberg-marquardt"];
%
% Open the file that we are writing to
%
   variable fp = fopen("alphasrch_"+filename+"_sherpa.dat","a");
%
   variable m;
   i=0;
   loop(length(r.alpha))
   {
      () = set_par("jdp.1","value",r.alpha[i]);
      () = set_frozen("jdp.1");
%
      foreach(method)
      {
         m=();
         () = sherpa_eval(m);
         loop(3)
         {
            () = run_fit;
         } 
      }
%
% Pileup fits can be frought with local minima. In this script, at
% large alpha (i.e., most piled photons are read as real events), very
% little power law is needed.  At small alpha (i.e., most piled
% photons are *not* read as real events), a power law is needed to
% make up the (probably many false) events at high energy.  Thus,
% there is a large dynamic range in power law normalization as we
% sweep from alpha = 1 to 0.  To help Sherpa catch this behaviour,
% we change the "step size" for the power law normalization as we
% sweep through alpha
%
      steps[2] = steps[2]*3.;
      () = set_par(parms[2],"delta",steps[2]);
%
      stats = get_delchi();
      r.xi[i] = sum( (stats)^2 );
%
      gp = get_par("phabs1.nH");
      r.nh[i] = gp.value;
%
      gp = get_par("bb1.norm");
      r.bbnorm[i]  = gp.value;
%
      gp = get_par("bb1.kT");
      r.kt[i] = gp.value;
%
      gp = get_par("pl1.norm");
      r.plnorm[i]  = gp.value;
%
      ()=fprintf(fp,"%11.2e  %11.2e  %11.2e  %11.2e  %11.2e  %12.3e \n",
         r.alpha[i], r.nh[i], r.kt[i], r.bbnorm[i], r.plnorm[i], r.xi[i]);
      i++;
    }
    () = fprintf(fp,"Alpha//Nh//kT//BB_norm//PL_norm//Xi\n");
    fclose(fp);
%
    return r;
}