/////////////////////////////////////////////////////////////////////////////
//
// EEezUtils
//
// Author: Jason C. Webb <jwebb@iucf.indiana.edu>
//
// Utility class to hold miscellaneous analysis functions, such as SMD
// fitting routines, etc...
//
// Implemened as a singleton.  Retrieve a pointer to the only instance
// of this class using the call
//
// EEezUtils *util = EEezUtils::instance();
//
/////////////////////////////////////////////////////////////////////////////

#include "EEezUtils.h"

#include "eeSinglePeak.h"
#include "eeDoublePeak.h"

#include <TH1F.h>
#include <TH2F.h>
#include <TF1.h>

#include <iostream>

ClassImp(EEezUtils);

EEezUtils *EEezUtils::sInstance = 0;

/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////

 EEezUtils::EEezUtils() 
{
  // Class Constructor.  
  // 
  // Use EEezUtils *util = EEezUtils::instance() to
  // get a pointer to this class.

}

/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////

 EEezUtils *EEezUtils::instance()
{
  // Returns the single instance of this class.
  sInstance = (sInstance != 0) ? sInstance : new EEezUtils();
  return sInstance;
}

/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////

 TF1 *EEezUtils::fitSmdSingleInit( TH1F *histo, Option_t *opts )
{
  // Initializes fit function for the specifed histogram.

  Int_t verbose = 0;
  if ( TString(opts).Contains("V") ) verbose++;

  //-- For each bin, sum with its neighbors and enter into new
  //-- histogram at that point.  We use this for an initial
  //-- guess at the mean of the distribution.  We also try
  //-- to reject mip signatures by requiring 3/3 strips to be
  //-- nonzero.
  Int_t nbins = histo -> GetNbinsX();
  TH1F *sum3 = (TH1F *)histo -> Clone("sum3");
  for ( Int_t i = 2; i <= nbins-1; i++ ) {
    Float_t sum = 0.;
    Int_t count = 0;
    for ( Int_t j = i-1; j<=i+1; j++ ) if ( histo->GetBinContent(j) > 1.0 ) count++;
    if ( count > 2 )
    for ( Int_t j = i-1; j<=i+1; j++ ) sum += histo -> GetBinContent(j);
    sum3 -> SetBinContent(i,sum);
  }
  //-- Maximum bin is here
  Int_t   max3  = sum3 -> GetMaximumBin();
  Float_t mean3 = sum3 -> GetBinLowEdge(max3) + 0.5;
  if ( verbose ) Info("fitSmdSingleInit","max bin = %i mean = %f", max3, mean3);

  //-- Starting values for mean
  Float_t start_mean = mean3;
  Float_t min_mean   = start_mean - 4.5;
  Float_t max_mean   = start_mean + 4.5;
  Int_t   bin_mean   = sum3 -> FindBin( start_mean );

  //-- Note: for now, just start with a constant for the width
  Float_t start_sigma = 0.84;
  Float_t min_sigma   = 0.65 * start_sigma;
  Float_t max_sigma   = 2.00 * start_sigma;

  //-- Estimate the yield
  Float_t start_yield = 1.5 * sum3 -> GetBinContent( bin_mean );
  Float_t min_yield   = 0.1 * start_yield;
  Float_t max_yield   = 2.5 * start_yield;

  TF1 *fit = new TF1("eeSinglePeak", eeSinglePeak, 0, 288, 5 );

  fit -> SetLineColor(2);
  fit -> SetFillStyle(0);
  fit -> SetParameter( 0, start_yield );
  fit -> SetParameter( 1, start_mean );
  fit -> SetParameter( 2, start_sigma );
  fit -> FixParameter( 3, 0.2 );
  fit -> FixParameter( 4, 3.5 );
  fit -> SetParameter( 3, 0.2 );
  fit -> SetParameter( 4, 3.5 );

  fit -> SetParLimits( 0, min_yield, max_yield );
  fit -> SetParLimits( 1, min_mean,  max_mean );
  fit -> SetParLimits( 2, min_sigma, max_sigma );

  delete sum3; // Clean up after self 

  return fit;

}

//////////////////////////////////////////////////////////////////

 TF1 *EEezUtils::fitSmdSingle( TH1F *histo, Option_t *opts, Float_t min, Float_t max )
{
  // Fit histogram of SMD response to a single shower-shape
  // function.  The fit is stored with the histogram, and
  // a pointer to the fit is returned (actually, a pointer
  // to the first function stored in the histogram, which
  // should always be this one... i.e. no "+" in the opts
  // arguement!).

  //-- Restrict the range of the histogram
  histo -> GetXaxis() -> SetRangeUser(min,max);

  //-- Initialize the fit function
  TF1 *fit = fitSmdSingleInit( histo, opts );

  //-- Perform the fit
  histo -> Fit(fit,TString(opts)+"0QN");

  //-- Restrict the range of the fit to +/- 3 sigma and refit
  Float_t mean = fit -> GetParameter(1);
  Float_t sigma = fit -> GetParameter(2);
  fit -> SetRange(mean-5.0*sigma,mean+5.0*sigma);
  histo -> Fit(fit,opts);

  //-- Unzoom the histogram
  histo -> GetXaxis() -> UnZoom();

  //-- And return with a pointer to the fit.  It looks like,
  //-- as usual, the root documentation is incomplete, and
  //-- the TF1 is copied into the histogram.  So we delete
  //-- our copy of this function and return a pointer to
  //-- the first function in the histogram's function list.

  return fit;

}

 TF1 *EEezUtils::fitSmdDouble( TH1F *histo, Option_t *opts, Float_t min, Float_t max, TF1 *infit )
{
  // Fit histogram of SMD response to a double shower-shape
  // function.  The fit is stored with the histogram, and
  // a pointer to the fit is returned.

  //-- Begin by fitting the histogram to a single shower-
  //-- shape distribution.

    TF1 *fit; 

  if ( infit == NULL ) {
  //-- Initialize parameters automatically
  TF1 *single = fitSmdSingle(histo,opts);

    //-- Compute the residual histogram
  TString name = histo -> GetName();
  TH1F *hresid = (TH1F *)histo -> Clone(name+"_residual");
  assert(single);
  hresid -> Add(single,-1.0);
  
  //-- Fit the residual histogram
  TF1 *fresid = fitSmdSingle(hresid,opts);

  //-- Now initialize the double-peak fit based on the
  //-- two single-peak fits
  fit = new TF1("eeDoublePeak", eeDoublePeak, 0, 288, 10 );
  fit -> SetLineColor(2);
  fit -> SetFillStyle(0);
  for ( Int_t i = 0; i < 5; i++ ) {
    fit -> SetParameter( i, single -> GetParameter(i) );
    Double_t min, max;
    single -> GetParLimits(i,min,max);
    fit -> SetParLimits( i, min, max );
  }
  for ( Int_t i = 0; i < 5; i++ ) {
    fit -> SetParameter( i+5, fresid -> GetParameter(i) );
    Double_t min, max;
    fresid -> GetParLimits(i,min,max);
    fit -> SetParLimits(i+5,min,max);
  }
  }
  else {
      fit = infit; 
  }

  //-- Perform the fit
  histo -> Fit( fit, opts );

  delete fit;
  return (TF1 *)histo -> GetListOfFunctions() -> At(0);

}