/**
* @file SPeakFinder.cpp
* @author Dan Saunders, on behalf of the SoLid collaboration.
* @date 17 Feb 2016
*/

#include "SPeakFinder.h"

//==============================================================================

//! Constructor setting up default values.
SPeakFinder::SPeakFinder(SDetector * dtr, SClipboard * cb) :
  ISAlgorithm(dtr, cb, "SPeakFinder"),
  m_nSaveExamples(5),
  m_nExamplesSaved(0),
  m_nPeaksFound(0),
  m_threshold(50),
  m_nWaveformsScanned(0),
  m_integralRange(20),
  m_shortIntegralRange(15),
  m_nSamplesRisingEdge(5)
{
  m_options.push_back(new SOptionDouble("Threshold", &m_threshold));
  m_options.push_back(new SOptionInt("IntegralRange", &m_integralRange));
  m_options.push_back(new SOptionInt("ShortIntegralRange", &m_shortIntegralRange));
  m_options.push_back(new SOptionInt("nSamplesRisingEdge", &m_nSamplesRisingEdge));
}


//==============================================================================

void SPeakFinder::initialize() 
{
  std::cout<<"[Note]: Threshold set to: "<<m_threshold<<" (ADC)"<<std::endl;
  cb()->histosFile()->mkdir((name() + "/Examples").c_str());
}


//==============================================================================

//! Peak finding.
/*! Currently a simple algorithm. Looks for a bin above some threshold, that 
 * should be set far above noise (say 1 PA) and takes a peak as a bin above 
 * its left and right neighbours. Loop is done over channels. Small number of
 * examples are kept for visual inspection.
 */
void SPeakFinder::execute()
{
  for (auto c : (*dtr()->channels())) {
    for (auto w : (*c->waveforms())) findPeaks(w, c);
    timeSortReconObjects(c->peaks());
  }
  timeSortReconObjects(cb()->peaks());
}


//==============================================================================

//! Peak finding for a specific waveform and channel - see execute()
void SPeakFinder::findPeaks(SWaveform * wf, SChannel * chan) 
{
  std::vector<SPeak*> peaks; // Temporary container for peaks found in this waveform.
  auto start = wf->samples()->begin()+1;
  auto end = wf->samples()->end()-1-m_integralRange;
  int iSample=0;
  for (auto iSamp = start; iSamp != end; iSamp++, iSample++) {

    // See if the sample is above its neighbours and above threshold. 
    // Allow cases where two touching samples have the same value.
    if (*iSamp > *(iSamp-1) &&
        *iSamp >= *(iSamp+1) && 
        *iSamp > m_threshold+chan->baseline()) {

      // We have a peak.
      SPeak * p = makePeak(wf, chan, iSample, iSamp);
      cb()->peaks()->push_back(p);
      chan->peaks()->push_back(p);
      peaks.push_back(p);
      m_nPeaksFound++;
    }
  }
  if (m_nExamplesSaved < m_nSaveExamples) saveExample(wf, &peaks);
  m_nWaveformsScanned++;
}

//==============================================================================

//! Creating the peak object from the waveform.
/*! Amplitude and time is currently taken as that of the peak bin. 
 * Integral is taken from the specified range.
 * SPeakFitter can be used to refine this for better estimates.
 */
SPeak * SPeakFinder::makePeak(SWaveform * wf, SChannel * chan, int iSample,
  std::vector<uint16_t>::iterator iSamp) 
{
  uint time = iSample + wf->time();
  double htime = wf->htime() + iSample*chan->sampleWidthNs();
  float amplitude = float(wf->samples()->at(iSample)) - chan->baseline();

  /* Recall integrals to include rising edge. First arg is start position, then 
   * the num samples to integrate over. If the sample doesn't exist, continue.
   */ 
  float integral = wf->integrate(iSample - m_nSamplesRisingEdge, m_integralRange, chan->baseline());
  float shortIntegral = wf->integrate(iSample - m_nSamplesRisingEdge, m_shortIntegralRange, chan->baseline());  
  
  SPeak * p = new SPeak(time, htime, chan, wf);
  p->setAmplitude(amplitude);
  p->setIntegral(integral);
  p->setShortIntegral(shortIntegral);
//  if (amplitude < 2000 && integral > 60000) h_examples.push_back(wf->getHistogram("eg_A" + std::to_string(amplitude) + "_I" + std::to_string(integral)));
  p->setWaveformPos(iSamp); // Since iSamp is an iterator.
  return p;
}


//==============================================================================

void SPeakFinder::finalize() 
{
  std::cout<<"[Note]: Number of peaks found: "<<m_nPeaksFound<<std::endl;
  std::cout<<"[Note]: cf number waveforms scanned: "<<m_nWaveformsScanned<<std::endl;

  cb()->histosFile()->cd((name() + "/Examples").c_str());
  for (auto h : h_examples) h->Write();
}


//==============================================================================

//! Saves the waveform and corresponding peaks in histograms.
void SPeakFinder::saveExample(SWaveform * wf, std::vector<SPeak*> * wf_peaks) 
{
  std::string wfName = "wf_" + std::to_string(m_nExamplesSaved);
  std::string hTitle = "Waveform " + std::to_string(m_nExamplesSaved);
  hTitle += "; DAQ samples (25 ns); Amplitude (ADC)";
  h_examples.push_back(wf->getHistogram(wfName, hTitle));
  std::string peakName = "peaksForWf_" + std::to_string(m_nExamplesSaved);
  std::string peakTitle = "Peaks for waveform " + std::to_string(m_nExamplesSaved);
  peakTitle += "; DAQ samples (25 ns); Amplitude (ADC)";
  TH1D *hPeak = new TH1D(peakName.c_str(), peakTitle.c_str(),
    wf->samples()->size(), 0, wf->samples()->size());

  for (auto p : (*wf_peaks)) 
    hPeak->SetBinContent(p->time()-wf->time()+1, p->amplitude());
  
  h_examples.push_back(hPeak);
  m_nExamplesSaved++;
}


//==============================================================================