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

#include "SPeakMon.h"


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

//! Constructor setting up default values.
SPeakMon::SPeakMon(SDetector * dtr, SClipboard * cb) :
  ISAlgorithm(dtr, cb, "SPeakMon"),
  m_htimeIsolationWindow(300)
{
  m_options.push_back(new SOptionDouble("hTimeIsolationWindow", &m_htimeIsolationWindow));
}


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

void SPeakMon::initialize() 
{
  // Estimate the ratio of amplitude to integral for suitable plot ranges.
  float ampToInt = 6;
  uint n = dtr()->nChannels();

  // Some histo ranges.
  float ampHighBin = 18000;
  float ampLowBin = -2000;
  float ampHighBin_zoom = 1000;
  float ampLowBin_zoom = -200;
  int nZoomBins = 400;
  int nBins_2d = 300;

  h_tBetween = new TH1D("tBetween", 
    "Time between consecutive peaks (all channels); Time (ns); N (/25ns);", 
    400, 0, 10000);

  // Useful to show people the peak amplitude threshold for some plots.
  SPeakFinder * peakFinder = (SPeakFinder*)cb()->algo("SPeakFinder");
  double peakThreshold = peakFinder->threshold();

  h_amplitude_zoom = new TH1D("amplitude_zoom", 
    ("Peak amplitudes, threshold: " + std::to_string(int(peakThreshold)) + " (zoomed); Amplitude (ADC); N;").c_str(), 
    nZoomBins, ampLowBin_zoom, ampHighBin_zoom);
  
  h_amplitudeTimeIsolated_zoom = new TH1D("amplitudeTimeIsolated_zoom", 
    "Peak amplitudes (time isolated and zoomed); Amplitude (ADC); N;", 
    nZoomBins, ampLowBin_zoom, ampHighBin_zoom);
  
  h_shortIntegral_zoom = new TH1D("shortIntegral_zoom", 
    "Peak short integrals (zoomed); Short Integral (ADC Samples); N;", 
    nZoomBins, ampToInt*ampLowBin_zoom, ampToInt*ampHighBin_zoom);
  
  h_shortIntegralTimeIsolated_zoom = new TH1D("shortIntegralTimeIsolated_zoom", 
    "Peak short integrals (time isolated and zoomed); Integral (ADC Samples); N;", 
    nZoomBins, ampToInt*ampLowBin_zoom, ampToInt*ampHighBin_zoom);

  // Per channel stuff.
  h_amplitudes_zoom = initPerChannelPlots1D("amplitude_zoom", 
    "Peak amplitudes (zoomed); Amplitude (ADC); N;", 
    nZoomBins, ampLowBin_zoom, ampHighBin_zoom);
  
  h_amplitudesTimeIsolated_zoom = initPerChannelPlots1D("amplitudeTimeIsolated_zoom", 
    "Peak amplitudes (time isolated and zoomed); Amplitude (ADC); N;", 
    nZoomBins, ampLowBin_zoom, ampHighBin_zoom);
  
  h_shortIntegrals_zoom = initPerChannelPlots1D("shortIntegral_zoom", 
    "Peak short integrals (zoomed); Short Integral (ADC Samples); N;", 
    nZoomBins, ampToInt*ampLowBin_zoom, ampToInt*ampHighBin_zoom);
  
  h_shortIntegralsTimeIsolated_zoom = initPerChannelPlots1D("shortIntegralTimeIsolated_zoom", 
    "Peak short integrals (time isolated and zoomed); Short Integral (ADC Samples); N;", 
    nZoomBins, ampToInt*ampLowBin_zoom, ampToInt*ampHighBin_zoom);
  

  // Integral vs amplitudes plots (including zoomed version with finer binning)
  h_integralVsAmplitude = new TH2D("integralVsAmplitude", 
    "Integral vs. Amplitude; Amplitude (ADC); Integral (ADC Samples);",
    nBins_2d, ampLowBin, ampHighBin, nBins_2d, ampToInt*ampLowBin, ampToInt*ampHighBin);
  
  h_integralVsAmplitude_zoom = new TH2D("integralVsAmplitude_zoom", 
    "Integral vs. Amplitude (zoom); Amplitude (ADC); Integral (ADC Samples);",
    nBins_2d, ampLowBin_zoom, ampHighBin_zoom, nBins_2d, ampToInt*ampLowBin_zoom, ampToInt*ampHighBin_zoom);
  
  // Amplitude Vs channel plots.
  h_amplitudeVsChannel = new TH2D("amplitudeVsChannel", 
    "Peak amplitudes for each channel; Saffron channel number; Amplitude (ADC);", 
    n, 0, n, nBins_2d, ampLowBin, ampHighBin);
  
  h_amplitudeVsChannel_zoom = new TH2D("amplitudeVsChannel_zoom", 
    "Peak amplitudes for each channel (zoomed y range); Saffron channel number; Amplitude (ADC);", 
    n, 0, n, nBins_2d, -ampLowBin_zoom, ampHighBin_zoom);

  // Histo file directories.
  cb()->histosFile()->mkdir((name() + "/amplitudes_zoom/").c_str());
  cb()->histosFile()->mkdir((name() + "/shortIntegrals_zoom/").c_str());
  cb()->histosFile()->mkdir((name() + "/amplitudesTimeIsolated_zoom/").c_str());
  cb()->histosFile()->mkdir((name() + "/shortIntegralsTimeIsolated_zoom/").c_str());
}


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

//! Fill histograms showing peak features.
/*! Care should be taken when editing, as some other code will use these histos,
 * in particular, the time isolation plots are used for online gain finding.
 */
void SPeakMon::execute() 
{
  uint iPeak = 0;
  double htimePrev = 0;
  for (auto p : (*cb()->peaks()) ) {
    h_amplitude_zoom->Fill(p->amplitude());
    h_shortIntegral_zoom->Fill(p->shortIntegral());

    h_integralVsAmplitude->Fill(p->amplitude(), p->integral());
    h_integralVsAmplitude_zoom->Fill(p->amplitude(), p->integral());
    h_amplitudeVsChannel->Fill(p->channel()->ID(), p->amplitude());
    h_amplitudeVsChannel_zoom->Fill(p->channel()->ID(), p->amplitude());
    h_amplitudes_zoom->at(p->channel()->ID())->Fill(p->amplitude());
    h_shortIntegrals_zoom->at(p->channel()->ID())->Fill(p->shortIntegral());

    if (iPeak != 0) h_tBetween->Fill(p->htime() - htimePrev);

    htimePrev = p->htime();
    iPeak++;
  }

  // Time isolation stuff.
  for (auto c : (*dtr()->channels())) {
    auto chanPeaks = c->peaks();
    uint nPeaks = chanPeaks->size();

    // Require at least 3 peaks to isolate just one.
    if (nPeaks < 3) continue;
    for (uint iPeak=1; iPeak<nPeaks-1; iPeak += 1) {
      auto p = chanPeaks->at(iPeak);

      // Get the time between this peak and its neighbours for this channel.
      double delt_backwards = p->htime() - chanPeaks->at(iPeak-1)->htime();
      double delt_forwards = chanPeaks->at(iPeak+1)->htime() - p->htime();

      // Shouldn't use this peak if either of these times are smaller than the isolation window.
      if (std::min(delt_backwards, delt_forwards) < m_htimeIsolationWindow) continue;
      h_amplitudeTimeIsolated_zoom->Fill(p->amplitude());      
      h_amplitudesTimeIsolated_zoom->at(c->ID())->Fill(p->amplitude());
      h_shortIntegralTimeIsolated_zoom->Fill(p->shortIntegral());
      h_shortIntegralsTimeIsolated_zoom->at(c->ID())->Fill(p->shortIntegral());
    }
  }
}


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

void SPeakMon::finalize() 
{
  h_amplitude_zoom->Write();
  h_shortIntegral_zoom->Write();
  h_amplitudeTimeIsolated_zoom->Write();
  h_shortIntegralTimeIsolated_zoom->Write();
  h_integralVsAmplitude->Write();
  h_integralVsAmplitude_zoom->Write();
  h_amplitudeVsChannel->Write();
  h_amplitudeVsChannel_zoom->Write();
  h_tBetween->Write();

  // Per channel stuff.
  cb()->histosFile()->cd((name() + "/amplitudes_zoom/").c_str());
  for (auto h : (*h_amplitudes_zoom)) h->Write();

  cb()->histosFile()->cd((name() + "/amplitudesTimeIsolated_zoom/").c_str());
  for (auto h : (*h_amplitudesTimeIsolated_zoom)) h->Write();

  cb()->histosFile()->cd((name() + "/shortIntegrals_zoom/").c_str());
  for (auto h : (*h_shortIntegrals_zoom)) h->Write();

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


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