/* ======================================================================== *\
!
! *
! * 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 analysing 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): Hendrik Bartko, 09/2004 <mailto:hbartko@mppmu.mpg.de> 
!   Author(s): Markus Gaug, 05/2004 <mailto:markus@ifae.es>
!   Author(s): Diego Tescaro, 05/2004 <mailto:tescaro@pd.infn.it>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */
//////////////////////////////////////////////////////////////////////////////
//
//   MExtractTimeAndChargeDigitalFilterPeakSearch
//
//   Hendrik has promised to write more documentation
//
//
//   The following variables have to be set by the derived class and 
//   do not have defaults:
//   - fNumHiGainSamples
//   - fNumLoGainSamples
//   - fSqrtHiGainSamples
//   - fSqrtLoGainSamples
//
// Input Containers:
//   MRawEvtData
//   MRawRunHeader
//   MPedestalCam
//
// Output Containers:
//   MArrivalTimeCam
//   MExtractedSignalCam
//
//////////////////////////////////////////////////////////////////////////////
#include "MExtractTimeAndChargeDigitalFilterPeakSearch.h"

#include <errno.h>
#include <fstream>

#include <TFile.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TString.h>
#include <TMatrix.h>

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

#include "MPedestalPix.h"
#include "MPedestalCam.h"

#include "MRawEvtPixelIter.h"

#include "MExtractedSignalCam.h"
#include "MExtractedSignalPix.h"

#include "MArrivalTimeCam.h"
#include "MArrivalTimePix.h"

ClassImp(MExtractTimeAndChargeDigitalFilterPeakSearch);

using namespace std;

const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgHiGainFirst             =  0;
const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgHiGainLast              = 18;
const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgLoGainFirst             =  3;
const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgLoGainLast              = 14;
const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgOffsetLeftFromPeak      = 1;
const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgOffsetRightFromPeak     = 2;
const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgPeakSearchWindowSize    = 2;
// --------------------------------------------------------------------------
//
// Default constructor. 
//
// Sets:
// - fOffsetLeftFromPeak  to fgOffsetLeftFromPeak
// - fOffsetRightFromPeak to fgOffsetRightFromPeak
//
MExtractTimeAndChargeDigitalFilterPeakSearch::MExtractTimeAndChargeDigitalFilterPeakSearch(const char *name, const char *title) 
{
    fName  = name  ? name  : "MExtractTimeAndChargeDigitalFilterPeakSearch";
    fTitle = title ? title : "Digital Filter";

    SetOffsetLeftFromPeak();
    SetOffsetRightFromPeak();
    SetPeakSearchWindowSize();
}

// --------------------------------------------------------------------------
//
// FindPeak
//
// Finds highest sum of "window" consecutive FADC slices in a pixel, and store
// in "startslice" the first slice of the group which yields the maximum sum.
// In "max" the value of the maximum sum is stored, in "sat" the number of 
// saturated slices.
//
void MExtractTimeAndChargeDigitalFilterPeakSearch::FindPeak(Byte_t *ptr, Byte_t &startslice, Int_t &max, 
							    Int_t &sat, Byte_t &satpos) const
{

  const Byte_t *end = ptr + fHiGainLast - fHiGainFirst + 1;

  sat = 0;
  satpos = 0;
  
  startslice = 0;
  Int_t sum=0;

  //
  // Calculate the sum of the first "fPeakSearchWindowSize" slices
  //
  sat = 0;
  Byte_t *p = ptr;

  while (p<ptr+fPeakSearchWindowSize)
    {
      sum += *p;
      if (*p++ >= fSaturationLimit)
        {
          if (sat == 0)
              satpos = p-ptr;
          sat++;
        }
    }

  //
  // Check for saturation in all other slices
  //
  while (p<end)
    if (*p++ >= fSaturationLimit)
      {
        if (sat == 0)
          satpos = p-ptr;
        sat++;
      }
  
  //
  // Calculate the i-th sum of n as
  //    sum_i+1 = sum_i + slice[i+fPeakSearchWindowSize] - slice[i]
  // This is fast and accurate (because we are using int's)
  //
  max=sum;
  for (p=ptr; p+fPeakSearchWindowSize<end; p++)
    {
      sum += *(p+fPeakSearchWindowSize) - *p;
      if (sum > max)
	{
	  max = sum;
	  startslice = p-ptr+1;
	}
    }

  return;
}


// --------------------------------------------------------------------------
//
// Process
//
// First find the pixel with highest sum of fPeakSearchWindowSize slices (default:2)
// "startslice" will mark the slice at which the highest sum begins for that pixel.
//
// Then define the beginning of the digital filter search window for ALL pixels as the slice
// before that: startslice-fOffsetLeftFromPeak until: startslice+fOffsetRightFromPeak
//
Int_t MExtractTimeAndChargeDigitalFilterPeakSearch::Process()
{

  MRawEvtPixelIter pixel(fRawEvt);

  Int_t   sat;
  Byte_t  satpos;
  ULong_t gsatpos    = 0;

  Int_t   maxsumhi   = -1000000;
  Int_t   numsat     = 0;
  Byte_t startslice;

  Byte_t hiGainFirstsave = fHiGainFirst;
  Byte_t hiGainLastsave  = fHiGainLast;
  Byte_t loGainFirstsave = fLoGainFirst;
  Byte_t loGainLastsave  = fLoGainLast; 

  Byte_t higainfirst     = fHiGainFirst;

  while (pixel.Next())
    {

      Int_t sumhi;
      sat = 0;

      FindPeak(pixel.GetHiGainSamples()+fHiGainFirst+fOffsetLeftFromPeak, startslice, sumhi, sat, satpos);

      if (sumhi > maxsumhi && sat == 0)
	{
	  maxsumhi     = sumhi;
	  higainfirst  = fHiGainFirst + startslice;
	}
      else if (sat)
        {
          numsat++;
          gsatpos += satpos;
        }
    }

  //
  // Check necessary for calibration events
  //
  if (numsat > fSignals->GetSize()*0.9)
    higainfirst = gsatpos/numsat - 1;

  //
  // Shift the start slice to the left:
  //
  if (higainfirst > fHiGainFirst + fOffsetLeftFromPeak)
    fHiGainFirst = higainfirst - fOffsetLeftFromPeak;
  else
    *fLog << warn << " High Gain ranges out of limits to the left!!! " << (Int_t)higainfirst << endl;

  //
  // Shift the last slice to the right:
  //
  if (higainfirst + fOffsetRightFromPeak + fWindowSizeHiGain < hiGainLastsave)
    fHiGainLast  = higainfirst + fOffsetRightFromPeak + fWindowSizeHiGain;
  else 
    *fLog << warn << " High Gain ranges out of limits to the right!!! " << (Int_t)higainfirst << endl;
      

  if ( fHiGainFirst+(Int_t)fOffsetLoGain > fLoGainFirst ) 
    fLoGainFirst = fHiGainFirst + (Int_t)fOffsetLoGain;
  else
    *fLog << inf << " High Gain ranges out of limits to the left!!! " << (Int_t)fHiGainFirst << endl;
  
  //
  // Make sure we will not integrate beyond the lo gain limit:
  //
  if (fLoGainFirst+fWindowSizeLoGain+fOffsetRightFromPeak <= pixel.GetNumLoGainSamples())
    fLoGainLast = fLoGainFirst+fWindowSizeLoGain+fOffsetRightFromPeak;
  else
    *fLog << warn << " Low Gain ranges out of limits to the right!!! " << endl;

  pixel.Reset();

  while (pixel.Next())
    {
      //
      // Find signal in hi- and lo-gain
      //
      Float_t sumhi =0., deltasumhi =0; // Set hi-gain of MExtractedSignalPix valid
      Float_t timehi=0., deltatimehi=0; // Set hi-gain of MArrivalTimePix valid
      Byte_t sathi=0;

      // Initialize fMaxBinContent for the case, it gets not set by the derived class
      fMaxBinContent = fLoGainSwitch + 1; 

      const Int_t pixidx = pixel.GetPixelId();
      const MPedestalPix  &ped = (*fPedestals)[pixidx];
      const Bool_t higainabflag = pixel.HasABFlag();

      FindTimeAndChargeHiGain(pixel.GetHiGainSamples()+fHiGainFirst, pixel.GetLoGainSamples(), 
                              sumhi, deltasumhi, 
                              timehi, deltatimehi, 
                              sathi, ped, higainabflag);

      //
      // Make sure that in cases the time couldn't be correctly determined
      // more meaningfull default values are assigned
      //
      if (timehi<0)
          timehi = -1;
      if (timehi>pixel.GetNumHiGainSamples())
          timehi = pixel.GetNumHiGainSamples();
      
      Float_t sumlo =0., deltasumlo =-1.; // invalidate logain of MExtractedSignalPix
      Float_t timelo=0., deltatimelo=-1;  // invalidate logain of MArrivalTimePix
      Byte_t satlo=0;
      
      //
      // Adapt the low-gain extraction range from the obtained high-gain time
      //
      if (pixel.HasLoGain() && (fMaxBinContent > fLoGainSwitch) )
      {
          deltasumlo  = 0; // make logain of MExtractedSignalPix valid
          deltatimelo = 0; // make logain of MArrivalTimePix valid

          fLoGainFirstSave = fLoGainFirst;
          const Byte_t logainstart = sathi 
            ? sathi + (Int_t)fLoGainStartShift
            : (Byte_t)(timehi + fLoGainStartShift);
        
	  fLoGainFirst = logainstart > fLoGainFirstSave ? logainstart : fLoGainFirstSave;

          if ( fLoGainFirst < fLoGainLast )
            {
              const Bool_t logainabflag = (higainabflag + pixel.GetNumHiGainSamples()) & 0x1;
              FindTimeAndChargeLoGain(pixel.GetLoGainSamples()+fLoGainFirst,
                                      sumlo, deltasumlo,
                                      timelo, deltatimelo,
                                      satlo, ped, logainabflag);
          }
          fLoGainFirst = fLoGainFirstSave;

          // Make sure that in cases the time couldn't be correctly determined
          // more meaningfull default values are assigned
          if (timelo<0)
              timelo = -1;
          if (timelo>pixel.GetNumLoGainSamples())
              timelo = pixel.GetNumLoGainSamples();
      }

      MExtractedSignalPix &pix = (*fSignals)[pixidx];
      MArrivalTimePix     &tix = (*fArrTime)[pixidx];
      pix.SetExtractedSignal(sumhi, deltasumhi,sumlo, deltasumlo);
      pix.SetGainSaturation(sathi, sathi, satlo);

      tix.SetArrivalTime(timehi, deltatimehi, timelo-fOffsetLoGain, deltatimelo);
      tix.SetGainSaturation(sathi, sathi, satlo);
 
    } /* while (pixel.Next()) */

  fArrTime->SetReadyToSave();
  fSignals->SetReadyToSave();

  fHiGainFirst  = hiGainFirstsave;
  fHiGainLast   = hiGainLastsave ; 
  fLoGainFirst  = loGainFirstsave;
  fLoGainLast   = loGainLastsave ; 

  return kTRUE;
}



// --------------------------------------------------------------------------
//
// Read the setup from a TEnv, eg:
//   MJPedestal.MExtractor.WindowSizeHiGain: 6
//   MJPedestal.MExtractor.WindowSizeLoGain: 6
//   MJPedestal.MExtractor.BinningResolutionHiGain: 10
//   MJPedestal.MExtractor.BinningResolutionLoGain: 10
//   MJPedestal.MExtractor.WeightsFile: filename
//
Int_t MExtractTimeAndChargeDigitalFilterPeakSearch::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{

  Bool_t rc = kFALSE;

  if (IsEnvDefined(env, prefix, "OffsetLeftFromPeak", print))
    {
      fOffsetLeftFromPeak = GetEnvValue(env, prefix, fOffsetLeftFromPeak);
      rc = kTRUE;
    }
  
  if (IsEnvDefined(env, prefix, "OffsetRightFromPeak", print))
    {
      fOffsetRightFromPeak = GetEnvValue(env, prefix, fOffsetRightFromPeak);
      rc = kTRUE;
    }
  
  if (IsEnvDefined(env, prefix, "PeakSearchWindowSize", print))
    {
      fPeakSearchWindowSize = GetEnvValue(env, prefix, fPeakSearchWindowSize);
      rc = kTRUE;
    }
  
  return MExtractTimeAndChargeDigitalFilter::ReadEnv(env, prefix, print) ? kTRUE : rc;
}


void MExtractTimeAndChargeDigitalFilterPeakSearch::Print(Option_t *o) const
{
    if (IsA()==Class())
        *fLog << GetDescriptor() << ":" << endl;

    MExtractTimeAndChargeDigitalFilter::Print(o);
    *fLog << " Offset from Peak left:   " << fOffsetLeftFromPeak   << endl;
    *fLog << " Offset from Peak right:  " << fOffsetRightFromPeak  << endl;
    *fLog << " Peak search window size: " << fPeakSearchWindowSize << endl;
}
