#include <RAT/TimingPeaks.hh>
using namespace RAT;
using namespace RAT::Classifiers;

#include <TH1D.h>
using namespace ROOT;

#include <math.h>
using namespace std;

void
TimingPeaks::Initialise( const std::string& )
{
  fPeakLevel = 0.05;
  fPeakBinWidth = 1;
}

DS::ClassifierResult
TimingPeaks::GetClassification()
{
  fClassifierResult.Reset();
  if( fPMTData.empty() )
    return fClassifierResult;

  // First produce a histogram of the hit times with 5ns bins
  const int nBins = 100; const double startBin = 0.0; const double endBin = 500.0;
  TH1D timeHist( "time", "time", nBins, startBin, endBin );
  timeHist.SetDirectory(0); // Memory manage myself

  for( vector<FitterPMT>::const_iterator iPMT = fPMTData.begin(); iPMT != fPMTData.end(); ++iPMT )
	timeHist.Fill( iPMT->GetTime() );
  //timeHist->Smooth( 1 );

  const double histMax = timeHist.GetMaximum();
  timeHist.Scale( 1.0 / histMax );

  TH1D derivative( "deriv", "deriv", nBins, startBin, endBin );
  derivative.SetDirectory(0); // Memory manage myself
  // Find simple derivative
  for( int iLoop = 3; iLoop <= nBins - 2; iLoop++ )
    {
      double value = timeHist.GetBinContent( iLoop + 1 ) - timeHist.GetBinContent( iLoop - 1 );
      value /= 2.0; // Two bins
	  derivative.SetBinContent( iLoop, value );
    }

  // Find and count peaks
  int numPeaks = 0;
  for( int iLoop = 1; iLoop <= nBins; iLoop += 1 )
    {
	  // Only a peak if it is fPeakBinWidth wide of bins greater than fPeakLevel
	  if( derivative.GetBinContent( iLoop ) > fPeakLevel ) // Start of a Peak
		{
		  int iPeak = iLoop + 1;
		  while( derivative.GetBinContent( iPeak ) > fPeakLevel && iPeak <= nBins )
			iPeak++;
		  // Peak has now ended
		  if( iPeak - iLoop > fPeakBinWidth )
			numPeaks++;
		  iLoop = iPeak + fPeakBinWidth; // Move fPeakBinWidth bin ahead and renew search
		}
    }

  fClassifierResult.SetClassification( "timingPeaks", numPeaks );
  fClassifierResult.SetValid( true );
  return fClassifierResult;
}