//
//  Class used to fit lines in the gamma spectrum 
//

#include "gamFit.h"


CgamFit::CgamFit(Int_t scale, Int_t halfWidth)
{
  CHAN_PER_BIN = scale;
  FIT_HW = halfWidth;
  for(int ii=0; ii<4; ii++) {
    error_last[ii] = 0;
    pars_last[ii] = 0;
  }
}

Double_t CgamFit::FindLine(TH1F* gamSpectrum, Double_t gamLineEnergy, Double_t halfW)
{
  gamSpectrum->GetXaxis()->SetRange((gamLineEnergy-halfW)/CHAN_PER_BIN,(gamLineEnergy+halfW)/CHAN_PER_BIN);
  Int_t peakBin = gamSpectrum->GetMaximumBin(); 
  Double_t peakBinE = peakBin*CHAN_PER_BIN;

  return peakBinE;
}

Double_t CgamFit::FitLine(TH1F* gamSpectrum, Double_t gamLineEnergy)
{
  Int_t maxCounts,minCounts,peakBin,firstChannelCont;
  Double_t expectedSigma, peakBinE;

  TF1 *gausMod = new TF1("gausMod","[0]+[1]*exp(-(x-[2])^2/(2*[3]^2))",0,8000);

  Double_t expectedGammaE = gamLineEnergy;
  gamSpectrum->GetXaxis()->SetRange((gamLineEnergy-FIT_HW)/CHAN_PER_BIN,(gamLineEnergy+FIT_HW)/CHAN_PER_BIN);
  peakBin = gamSpectrum->GetMaximumBin(); 
  peakBinE = peakBin*CHAN_PER_BIN;
  maxCounts = gamSpectrum->GetMaximum(); 
  minCounts = gamSpectrum->GetMinimum(0); 
  expectedSigma = this->GetEstSigma(gamSpectrum, peakBin, maxCounts);

  gausMod->SetParameters(minCounts,maxCounts,peakBinE,expectedSigma);
  gausMod->SetParLimits(0,minCounts/2,(minCounts+maxCounts)/2);
  gausMod->SetParLimits(1,0.6*(maxCounts-minCounts),1.2*(maxCounts-minCounts));
  gausMod->SetParLimits(2,peakBinE-expectedSigma,peakBinE+expectedSigma);
  gausMod->SetParLimits(3,expectedSigma/5,expectedSigma*5);
  
  gamSpectrum->Fit("gausMod","MQ+","",peakBinE-FIT_HW,peakBinE+FIT_HW);
  
  for(int ii=0; ii<4; ii++) {
    pars_last[ii] = gausMod->GetParameter(ii);
    error_last[ii] = gausMod->GetParError(ii);
  }

  return pars_last[2];
  
}

Double_t CgamFit::GetEstSigma(TH1F* gamSpectrum, Int_t peakBin, Int_t maxCounts)
{
  Int_t step[2] = {1,1};
  Int_t binCounts;

  binCounts = maxCounts;
  while(binCounts>=maxCounts/2) {
    binCounts = gamSpectrum->GetBinContent(peakBin+step[0]);
    step[0]++;
  }

  binCounts = maxCounts;
  while(binCounts>=maxCounts/2) {
    binCounts = gamSpectrum->GetBinContent(peakBin+step[1]);
    step[1]++;
  }

  return (CHAN_PER_BIN*(step[0]+step[1]-1))/2.35; 

}

Double_t CgamFit::GetParErr(Int_t index)
{
  return error_last[index];
}

Double_t CgamFit::GetPar(Int_t index) 
{
  return pars_last[index];
}

Double_t CgamFit::GetArea()
{
  return (pars_last[3]*pars_last[1])/CHAN_PER_BIN;
}