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

#include "SOnlineNeutronID.h"

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

//! Constructor setting up default values.
SOnlineNeutronID::SOnlineNeutronID(SDetector * dtr, SClipboard * cb) :
  ISAlgorithm(dtr, cb, "SOnlineNeutronID"),
  m_totSampleWindow(40), // As in SM1.
  m_totThreshold(50),
  m_nPeaksHtimeWindow(1000)
{
  m_options.push_back(new SOptionInt("TotSampleWindow", &m_totSampleWindow));
  m_options.push_back(new SOptionInt("TotThreshold", &m_totThreshold));
  m_options.push_back(new SOptionDouble("nPeaksHtimeWindow", &m_nPeaksHtimeWindow));
}


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

void SOnlineNeutronID::initialize() 
{}


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

//! Online neutron ID.
/*! Current implimentation is very simple. Loop over all events, and draw histos
 * of usual nID parameters (tot, nPeaks, I/A). The parameters are summed across 
 * all channels involved in the event (no selections placed of the channels yet).
 * Since nPeaks, tot, and IonA are not yet found for peaks, this needs to be done 
 * manually. Things will likely move around after the starterkit.
 */ 
void SOnlineNeutronID::execute() 
{
  for (auto event : (*cb()->events())) {
    findIonA(event);
    findTot(event);
    findNpeaks(event);
  }
}


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

void SOnlineNeutronID::finalize() 
{}


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

//! For each event, assign a value of IonA. Take as the sum from all channels.
void SOnlineNeutronID::findIonA(SEvent * event) 
{
  double IonA = 0;
  for (auto p : (*event->peaks())) IonA += p->integral()/p->amplitude();
  event->setNID_IonA(IonA);
}


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

//! For each event, assign a value of tot.
/*! Find number of samples above some threshold for a preset window (in units of 
 * sampled). Need to loop over the waveforms in the vicinity of the event for 
 * m_totSampleWindow samples, and see if over threshold. Sum over channels.
 * Use the iterator provided by each peak. Need to be careful not to access 
 * something that doesn't exist.
 */
void SOnlineNeutronID::findTot(SEvent * event) 
{
  // Has to be updated for NZS data.
  unsigned int tot = 0;
  for (auto p : (*event->peaks())) {
    uint thres = uint(m_totThreshold) + uint(p->channel()->baseline());
    // Loop over the samples at the peak.
    auto upperBound = std::min(p->waveform()->samples()->end(), p->waveformPos()+m_totSampleWindow);
    for (auto iSamp = p->waveformPos(); 
         iSamp != upperBound;
         iSamp++) {
      if ((*iSamp) > thres) tot++;
    }
  }
  event->setNID_tot(tot);
}


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

//! For each event, assign a value of nPeaks.
/*! Need to look at the channels forming the event, and look for other peaks. 
 * Since Saffron2 has already done peak finding at a low threshold, we'll
 * count those peaks from the forming channels in the htime window.
 */
void SOnlineNeutronID::findNpeaks(SEvent * event) 
{
  unsigned int nPeaks = 0;  
  for (auto c : (*event->channels())) {
    // Loop over peaks on this channel.
    for (auto p : (*c->peaks())) {
      double delt = event->htime() - p->htime();
      if (delt >= 0 && delt < m_nPeaksHtimeWindow) nPeaks++;
    }
  }
  event->setNID_nPeaks(nPeaks);
}


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