//
// pulse_analyser.cc
//
// Class pulse_analyser provides a common set of functions for
// evaluating the height of pulses recorded in a digitized time
// series from the fADC250 VME module operating in raw mode.
//
// author: richard.t.jones at uconn.edu
// version: august 20, 2014
// notes:
//
// 1) Three distinct ways of finding pulses in the time series
//    data from the fADC250 are supported:
//       a) pulser trigger - assumes the pulses are generated by
//          light from the fast laser diode pulser mounted 
//          inside the fiber testing dark box at UConn. The 
//          pulser generates its own separate trigger signal 
//          that is fed into a discriminator in the VME crate
//          and used to trigger the DAQ. All pulses in this
//          setting are assumed to come at the same time.
//       b) summer trigger - assumes that the DAQ trigger is
//          generated by the VME discriminator watching one or
//          more of the same signals that are being digitized
//          by the fADC250 module. The is similar to case (a)
//          above except that where the pulse occurs in the time
//          series is different.
//       c) scanning trigger - assumes that the DAQ trigger is
//          random (eg. a regular clock) and that the analyzer
//          must search for pulses in the time series and set
//          a unique integration window around each pulse it 
//          finds. A fixed number of samples at the start and
//          end of the pulse are left out of the search, and
//          used to set the pedestal level for the event.
// 2) All of the measured pulse heights are divided by a common
//    normalization factor for the event which is evaluated based
//    on the peak height observed in the laser diode amplitude
//    feedback signal, which is digitized in fADC250 channel 0.
// 3) Basic consistency checks are made on the trace data before
//    a pulse height is returned to make sure the baseline level
//    is stable, otherwise the algorithms return 0 for the pulse
//    height.
// 4) The overall scale is designed to put healthy fiber pulse
//    heights from the laser diode in the range 700 - 1000.
//    There is no other significance to it than that.

#ifndef _PULSE_ANALYZER_
#include <pulse_analyzer.h>
#include <math.h>

#define FIRST_SPLIT_CHANNEL 99
#define LAST_SPLIT_CHANNEL 99

double pulse_analyzer::pheight_trig_pulser(int channel)
{
   int pulse_start;
   if (channel >= FIRST_SPLIT_CHANNEL &&
       channel <= LAST_SPLIT_CHANNEL)
      pulse_start = 105;
   else
      pulse_start = 86;
   pulse_start += 28;
   return pheight_trig_scanning(channel, pulse_start, pulse_start);
}

double pulse_analyzer::pheight_trig_summer(int channel)
{
   int pulse_start;
   if (channel >= FIRST_SPLIT_CHANNEL &&
       channel <= LAST_SPLIT_CHANNEL)
      pulse_start = 105;
   else
      pulse_start = 86;
   return pheight_trig_scanning(channel, pulse_start, pulse_start);
}

double pulse_analyzer::pheight_trig_scanning(int channel)
{
   int pulse_start = 12;
   int pulse_stop = fWindowSize - fPulseWindowWidth - 12;
   return pheight_trig_scanning(channel, pulse_start, pulse_stop);
}

double pulse_analyzer::pheight_trig_scanning(int channel,
                                             int scan_start,
                                             int scan_stop)
{
   double pedestal = pedestal_value(channel);
   if (pedestal <= 0) {
      return -99;
   }

   double pulse_max = 0;
   for (int i0=scan_start; i0 <= scan_stop; ++i0) {
      double pulse_sum = 0;
      for (int i=0; i < fPulseWindowWidth; ++i) {
         if (i+i0 < fWindowSize) {
            pulse_sum += fADC_raw[channel*fWindowSize + i+i0] - pedestal;
         }
      }
      pulse_max = (pulse_sum > pulse_max)? pulse_sum : pulse_max;
   }
   return pulse_max / fStandardPulseNormFactor;
}

double pulse_analyzer::pedestal_value(int channel)
{
   int ped_window[] = {1, 12};
   double mean[2] = {0,0};
   double rms[2] = {0,0};
   int count = 0;
   for (int i=ped_window[0]; i < ped_window[1]; ++i) {
      double h0 = fADC_raw[channel*fWindowSize + i-1];
      double h1 = fADC_raw[channel*fWindowSize + fWindowSize-i];
      mean[0] += h0;
      mean[1] += h1;
      rms[0] += h0*h0;
      rms[1] += h1*h1;
      ++count;
   }
   mean[0] /= count;
   mean[1] /= count;
   rms[0] /= count;
   rms[1] /= count;
   rms[0] = sqrt(rms[0] - mean[0]*mean[0] + 1e-10);
   rms[1] = sqrt(rms[1] - mean[1]*mean[1] + 1e-10);
   if (rms[0] > 0.03 * mean[0] || rms[1] > 0.10 * mean[1]) {
      return 0;
   }
   return mean[0];
}

double pulse_analyzer::pulse_normalization()
{
   double pedsum[2] = {0,0};
   for (int i=fPulseNormPedWindow[0]; i < fPulseNormPedWindow[1]; ++i) {
      pedsum[1] += fADC_raw[fPulseNormChannel*fWindowSize + i];
      pedsum[0]++;
   }
   double sigsum[2] = {0,0};
   for (int i=fPulseNormSigWindow[0]; i < fPulseNormSigWindow[1]; ++i) {
      sigsum[1] += fADC_raw[fPulseNormChannel*fWindowSize + i];
      sigsum[0]++;
   }
   return (sigsum[1] - pedsum[1]*sigsum[0]/pedsum[0])/fPulseNormReference;
}