/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analyzing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!   Author(s): Thomas Bretz,   02/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
!   Author(s): Hendrik Bartko, 01/2004 <mailto:hbartko@mppmu.mpg.de>
!   Author(s): Markus Gaug,    09/2004 <mailto:markus@ifae.es> 
!
!   Copyright: MAGIC Software Development, 2002-2004
!
!
\* ======================================================================== */
//////////////////////////////////////////////////////////////////////////////
//
//   MExtractSlidingWindow
//
//   Extracts the signal from a sliding window of size fHiGainWindowSize and 
//   fLoGainWindowSize, respectively. The signal is the one which maximizes 
//   the clock-noise and pedestal-corrected integral contents.
//
//   The amplitude-weighted arrival time is calculated from the window with 
//   the highest integral using the following formula:
//
//   t = sum(s(i) * i) / sum(i) 
// 
//   where i denotes the FADC slice index and s(i) the clock-noise and 
///  pedestal-corrected FADC value at slice index i. The sum runs over the 
//   extraction window size.
//
//   Call: SetRange(higainfirst, higainlast, logainfirst, logainlast) 
//         to modify the ranges in which the window is allowed to move. 
//
//   Defaults are: 
// 
//   fHiGainFirst =  fgHiGainFirst =  0 
//   fHiGainLast  =  fgHiGainLast  =  14
//   fLoGainFirst =  fgLoGainFirst =  2 
//   fLoGainLast  =  fgLoGainLast  =  14
//
//  Call: SetWindowSize(windowhigain, windowlogain) 
//        to modify the sliding window widths. Windows have to be an even number. 
//        Odd numbers are possible since the clock-noise is corrected for.
//
//  Defaults are:
//
//  fHiGainWindowSize = 6
//  fLoGainWindowSize = 6
//
//////////////////////////////////////////////////////////////////////////////
#include "MExtractTimeAndChargeSlidingWindow.h"

#include "MPedestalPix.h"

#include "MLog.h"
#include "MLogManip.h"

ClassImp(MExtractTimeAndChargeSlidingWindow);

using namespace std;

const Byte_t  MExtractTimeAndChargeSlidingWindow::fgHiGainFirst  = 2;
const Byte_t  MExtractTimeAndChargeSlidingWindow::fgHiGainLast   = 16;
const Byte_t  MExtractTimeAndChargeSlidingWindow::fgLoGainFirst  = 2;
const Byte_t  MExtractTimeAndChargeSlidingWindow::fgLoGainLast   = 14;
const Byte_t  MExtractTimeAndChargeSlidingWindow::fgHiGainWindowSize = 6;
const Byte_t  MExtractTimeAndChargeSlidingWindow::fgLoGainWindowSize = 6;
// --------------------------------------------------------------------------
//
// Default constructor.
//
// Calls: 
// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast)
// 
// Initializes:
// - fWindowSizeHiGain to fgHiGainWindowSize
// - fWindowSizeLoGain to fgLoGainWindowSize
//
MExtractTimeAndChargeSlidingWindow::MExtractTimeAndChargeSlidingWindow(const char *name, const char *title) 
    : fWindowSizeHiGain(fgHiGainWindowSize), 
      fWindowSizeLoGain(fgLoGainWindowSize),
      fHiGainSignal(NULL), fLoGainSignal(NULL)
{
  
  fName  = name  ? name  : "MExtractTimeAndChargeSlidingWindow";
  fTitle = title ? title : "Calculate arrival times and charges using a sliding window";

  SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
}


// --------------------------------------------------------------------------
//
// Destructor: Deletes the arrays
//
MExtractTimeAndChargeSlidingWindow::~MExtractTimeAndChargeSlidingWindow()
{
  
  if (fHiGainSignal)
    delete [] fHiGainSignal;
  if (fLoGainSignal)
    delete [] fLoGainSignal;
  
}

//-------------------------------------------------------------------
//
// Set the ranges
//
// Calls:
// - MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
// - SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain);
//
void MExtractTimeAndChargeSlidingWindow::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast)
{

  MExtractor::SetRange(hifirst, hilast, lofirst, lolast);

  //
  // Redo the checks if the window is still inside the ranges
  //
  SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain);
  
}

// -----------------------------------------------------------------------------------------
//
// Checks:
// - if a window is bigger than the one defined by the ranges, set it to the available range
// - if a window is smaller than 2, set it to 2
// 
// Sets:
// - fNumHiGainSamples to: (Float_t)fWindowSizeHiGain
// - fNumLoGainSamples to: (Float_t)fWindowSizeLoGain
// - fSqrtHiGainSamples to: TMath::Sqrt(fNumHiGainSamples)
// - fSqrtLoGainSamples to: TMath::Sqrt(fNumLoGainSamples)  
//  
void MExtractTimeAndChargeSlidingWindow::SetWindowSize(Byte_t windowh, Byte_t windowl)
{
  
  const Byte_t availhirange = (fHiGainLast-fHiGainFirst+1) & ~1;

  if (fWindowSizeHiGain > availhirange)
    {
      *fLog << warn << GetDescriptor() 
            << Form("%s%2i%s%2i%s%2i%s",": Hi Gain window size: ",(int)fWindowSizeHiGain,
                    " is bigger than available range: [",(int)fHiGainFirst,",",(int)fHiGainLast,"]") << endl;
      *fLog << warn << GetDescriptor() 
            << ": Will set window size to: " << (int)availhirange << endl;
      fWindowSizeHiGain = availhirange;
    }
  
  if (fWindowSizeHiGain<2) 
    {
      fWindowSizeHiGain = 2;
      *fLog << warn << GetDescriptor() << ": High Gain window size too small, set to two samples" << endl;
    }
  
  if (fLoGainLast != 0 && fWindowSizeLoGain != 0)
    {
      const Byte_t availlorange = (fLoGainLast-fLoGainFirst+1) & ~1;
      
      if (fWindowSizeLoGain > availlorange)
        {
          *fLog << warn << GetDescriptor() 
                << Form("%s%2i%s%2i%s%2i%s",": Lo Gain window size: ",(int)fWindowSizeLoGain,
                        " is bigger than available range: [",(int)fLoGainFirst,",",(int)fLoGainLast,"]") << endl;
          *fLog << warn << GetDescriptor() 
                << ": Will set window size to: " << (int)availlorange << endl;
          fWindowSizeLoGain = availlorange;
        }
      
      if (fWindowSizeLoGain<2) 
        {
          fWindowSizeLoGain = 2;
          *fLog << warn << GetDescriptor() << ": Low Gain window size too small set to two samples" << endl;
        }
    }
  
  fNumHiGainSamples = (Float_t)fWindowSizeHiGain;
  fNumLoGainSamples = fLoGainLast ? (Float_t)fWindowSizeLoGain : 0.;
  
  fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
  fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);

}

// --------------------------------------------------------------------------
//
// ReInit
//
// Calls:
// - MExtractTimeAndCharge::ReInit(pList);
// - Deletes all arrays, if not NULL
// - Creates new arrays according to the extraction range
//
Bool_t MExtractTimeAndChargeSlidingWindow::ReInit(MParList *pList)
{

  if (!MExtractTimeAndCharge::ReInit(pList))
    return kFALSE;

  if (fHiGainSignal)
    delete [] fHiGainSignal;
  if (fLoGainSignal)
    delete [] fLoGainSignal;
  
  Int_t range = fHiGainLast - fHiGainFirst + 1 + fHiLoLast;

  fHiGainSignal = new Float_t[range];
  memset(fHiGainSignal,0,range*sizeof(Float_t));

  range = fLoGainLast - fLoGainFirst + 1;

  fLoGainSignal = new Float_t[range];
  memset(fLoGainSignal,0,range*sizeof(Float_t));

  return kTRUE;
}

// --------------------------------------------------------------------------
//
// Calculates the arrival time for each pixel 
//
void MExtractTimeAndChargeSlidingWindow::FindTimeAndChargeHiGain(Byte_t *first, Byte_t *logain, Float_t &sum, Float_t &dsum, 
                                                          Float_t &time, Float_t &dtime, 
                                                          Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
{
  
  Int_t range = fHiGainLast - fHiGainFirst + 1;
  const Byte_t *end = first + range;
  Byte_t       *p  = first;
  
  Float_t max = 0;  // highest integral content of all windows
  Int_t count = 0;  // counter to recognize the AB-flag

  const Float_t pedes  = ped.GetPedestal();
  const Float_t ABoffs = ped.GetPedestalABoffset();

  Float_t PedMean[2];
  PedMean[0] = pedes + ABoffs;
  PedMean[1] = pedes - ABoffs;

  //
  // Check for saturation in all other slices
  //
  while (p<first+fWindowSizeHiGain)
    {
      const Int_t ids      = fHiGainFirst + count;
      const Float_t signal = (Float_t)*p - PedMean[(ids+abflag) & 0x1];      
      sum                 += signal;
      fHiGainSignal[count] = signal;
      
      if (*p++ >= fSaturationLimit)
        sat++;

      count++;
    }

  //
  // Check for saturation in all other slices
  //
  while (p<end)
    if (*p++ >= fSaturationLimit)
      sat++;
  
  if (IsNoiseCalculation())
    return;
   
  //
  // Calculate the i-th sum as
  //    sum_i+1 = sum_i + slice[i+8] - slice[i]
  // This is fast and accurate (because we are using int's)
  //
  count          = 0;
  max            = sum;
  Int_t idx      =  0;  // idx of the first slice of the maximum window
  
  for (p=first; p+fWindowSizeHiGain<end; p++)
    {
      
      const Int_t ids      = fHiGainFirst + count + fWindowSizeHiGain;
      const Float_t signal = (Float_t)*(p+fWindowSizeHiGain) - PedMean[(ids+abflag) & 0x1];      
      sum                 += signal - fHiGainSignal[count];
      fHiGainSignal[count + fWindowSizeHiGain] = signal;
      
      if (sum>max)
        {
          max   = sum;
          idx   = count+1;
        }
      count++;
    }
  
  if (fHiLoLast != 0)
    {
      
      // 
      // overlap bins
      // 
      Byte_t *l = logain;
      
      while (p < end && l < logain+fHiLoLast)
        {
          
          const Int_t   ids    = fHiGainFirst + count + fWindowSizeHiGain;
          const Float_t signal = (Float_t)*l - PedMean[(ids+abflag) & 0x1];          
          sum                 += signal - fHiGainSignal[count];
          fHiGainSignal[count + fWindowSizeHiGain] = signal;
          count++;
          if (*l++ >= fSaturationLimit)
            sat++;
          
          if (sum>max)
            {
              max   = sum;
              idx   = count+1;
            }
          p++;
        }
      
      if (fHiLoLast > fWindowSizeHiGain)
        {
          while (l < logain + fHiLoLast)
            {
              const Int_t   ids    = fHiGainFirst + count + fWindowSizeHiGain;
              const Float_t signal = (Float_t)*l - PedMean[(ids+abflag) & 0x1];          
              sum                 += signal - fHiGainSignal[count];
              fHiGainSignal[count+fWindowSizeHiGain] = signal;
              count++;
              if (*(l+fWindowSizeHiGain) >= fSaturationLimit)
                sat++;
              
              if (sum>max)
                {
                  max   = sum;
                  idx   = count+1;
                }
              l++;
            } /* while (l < logain + fHiLoLast) */
        } /* if (fHiLoLast > fWindowSizeHiGain) */
    } /* if (fHiLoLast != 0) */

  //
  // now calculate the time for the maximum window
  //
  Float_t timesignalsum  = 0.;
  Int_t   timesquaredsum = 0;
  
  for (Int_t i=idx; i<idx+fWindowSizeHiGain; i++)
    {
      timesignalsum   += fHiGainSignal[i]*i;
      timesquaredsum  += i*i;
    }
  
  sum   = max;

  time  = sum != 0 ? timesignalsum / max  + Float_t(fHiGainFirst) : 1.;
  dtime = sum != 0 ? ped.GetPedestalRms() / max * sqrt(timesquaredsum - fWindowSizeHiGain*time) : 1.;

}


// --------------------------------------------------------------------------
//
// Calculates the arrival time for each pixel 
//
void MExtractTimeAndChargeSlidingWindow::FindTimeAndChargeLoGain(Byte_t *first, Float_t &sum, Float_t &dsum, 
                                                          Float_t &time, Float_t &dtime, 
                                                          Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag) 
{
  
  Int_t range = fLoGainLast - fLoGainFirst + 1;
  const Byte_t *end = first + range;
  Byte_t       *p  = first;
  
  Float_t max = 0;  // highest integral content of all windows
  Int_t count = 0;  // counter to recognize the AB-flag

  Float_t pedes        = ped.GetPedestal();
  const Float_t ABoffs = ped.GetPedestalABoffset();

  Float_t PedMean[2];
  PedMean[0] = pedes + ABoffs;
  PedMean[1] = pedes - ABoffs;

  //
  // Check for saturation in all other slices
  //
  while (p<first+fWindowSizeLoGain)
    {
      const Int_t ids      = fLoGainFirst + count;
      const Float_t signal = (Float_t)*p - PedMean[(ids+abflag) & 0x1];      
      sum                 += signal;
      fLoGainSignal[count] = signal;
      
      if (*p++ >= fSaturationLimit)
        sat++;

      count++;
    }

  //
  // Check for saturation in all other slices
  //
  while (p<end)
    if (*p++ >= fSaturationLimit)
      sat++;
  
  if (IsNoiseCalculation())
    return;

  //
  // Calculate the i-th sum as
  //    sum_i+1 = sum_i + slice[i+8] - slice[i]
  // This is fast and accurate (because we are using int's)
  //
  count          = 0;
  max            = sum;
  Int_t idx      =  0;  // idx of the first slice of the maximum window
  
  for (p=first; p+fWindowSizeLoGain<end; p++)
    {

      const Int_t ids      = fLoGainFirst + count + fWindowSizeLoGain;
      const Float_t signal = (Float_t)*(p+fWindowSizeLoGain) - PedMean[(ids+abflag) & 0x1];      
      sum                 += signal - fLoGainSignal[count];
      fLoGainSignal[count + fWindowSizeLoGain] = signal;

      if (sum>max)
	{
	  max   = sum;
          idx   = count+1;
	}
      count++;
    }
  
  //
  // now calculate the time for the maximum window
  //
  Float_t timesignalsum  = 0;
  Int_t   timesquaredsum = 0;
  
  for (Int_t i=idx; i<idx+fWindowSizeLoGain; i++)
    {
      timesignalsum   += fLoGainSignal[i]*i;
      timesquaredsum  += i*i;
    }
  
  sum   = max;

  time  = sum != 0 ? timesignalsum / max  + Float_t(fLoGainFirst) : 1.;
  dtime = sum != 0 ? ped.GetPedestalRms() / max * sqrt(timesquaredsum - fWindowSizeLoGain*time) : 1.;
}

// --------------------------------------------------------------------------
//
// In addition to the resources of the base-class MExtractor:
//   MJPedestal.MExtractor.WindowSizeHiGain: 6
//   MJPedestal.MExtractor.WindowSizeLoGain: 6
//
Int_t MExtractTimeAndChargeSlidingWindow::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{
  
  Byte_t hw = fWindowSizeHiGain;
  Byte_t lw = fWindowSizeLoGain;
  
  Bool_t rc = kFALSE;
  
  if (IsEnvDefined(env, prefix, "HiGainWindowSize", print))
    {
      hw = GetEnvValue(env, prefix, "HiGainWindowSize", hw);
      rc = kTRUE;
    }
  if (IsEnvDefined(env, prefix, "LoGainWindowSize", print))
    {
      lw = GetEnvValue(env, prefix, "LoGainWindowSize", lw);
      rc = kTRUE;
    }
  
  if (rc)
    SetWindowSize(hw, lw);
  
  return MExtractTime::ReadEnv(env, prefix, print) ? kTRUE : rc;

}