/*!
* @file STools.cpp
* @author Nick Ryder, on behalf of the SoLid collaboration.
* @date 18 Nov 2016
*/

#include "STools.h"

#include <TF1.h>

float PedestalFinder(TH1D *hin){
  TH1D *hNew = (TH1D*)hin->Clone("hNew");
  int maxBin = hNew->GetMaximumBin();
  // temporary fit width (in bins), will receive a final value later on
  int fitWidth = 1;

  while (fitWidth < 50 && maxBin - fitWidth > 0 &&
         hNew->GetBinContent(maxBin - fitWidth) != 0){
    // look at the lower side of the peak for an empty bin,
    // also some security checks to not go out of range
    fitWidth ++;
  }

  double lowSide = hNew->GetXaxis()->GetBinCenter(maxBin - fitWidth);
  double highSide = hNew->GetXaxis()->GetBinCenter(maxBin + fitWidth);
  TF1 *f_g = new TF1("gaus", "gaus", lowSide, highSide);

  // set a seed for the fit
  f_g->SetParameters(
    hNew->GetBinContent(maxBin),            // Gaussian amplitude
    hNew->GetXaxis()->GetBinCenter(maxBin), // gaussian mean
    0.1*(highSide - lowSide));              // gaussian width
  int fitRes = hNew->Fit(f_g,"rq"); // do the fit

  float mu;
  if (fitRes != 0){
    // this is a bad fit, use maximum bin
    std::cout<<"[Warning]: Pedestal fit failed.\n";
    std::cout<<"   Using position of maximum bin instead"<<std::endl;
    mu = hNew->GetXaxis()->GetBinCenter(maxBin);
  } else {
    // good fit, use fit result
    mu = f_g->GetParameter(1);
  }
  delete f_g;
  delete hNew;

  return mu;
}