source: trunk/MagicSoft/Mars/msignal/MExtractFixedWindowPeakSearch.cc@ 4342

/* ======================================================================== *\
   !
   ! *
   ! * 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): Abelardo Moralejo,04/2004 <mailto:moralejo@pd.infn.it>
   !              Markus Gaug      ,04/2004 <mailto:markus@ifae.es>
   !   Copyright: MAGIC Software Development, 2000-2004
   !
   !
   \* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//   MExtractFixedWindowPeakSearch
//
//  Calculate the signal integrating fWindowSize time slices, the same for 
//  all pixels. The integrated slices change from event to event, since the
//  pulse positions in the FADC jump between events, but apparently in a 
//  "coherent" fashion. We first loop over all pixels and find the one 
//  which has the highest sum of fPeakSearchWindowSize (default: 4) consecutive 
//  non-saturated high gain slices. The position of the peak in this pixel 
//  determines the integration window position for all pixels of this event. 
//  For the moment we neglect time delays between pixels (which are in most 
//  cases small). The class is based on MExtractSlidingWindow.
//
//  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 =  3 
//   fLoGainLast  =  fgLoGainLast  =  14
//
//  Call: SetWindows(windowhigain, windowlogain,peaksearchwindow) 
//  to modify the sliding window widths. Windows have to be an even number. 
//  In case of odd numbers, the window will be modified.
//
//  Defaults are:
//
//   fHiGainWindowSize     = fgHiGainWindowSize     = 6
//   fLoGainWindowSize     = fgLoGainWindowSize     = 6
//   fPeakSearchWindowSize = fgPeakSearchWindowSize = 4
//   fLowGainPeakShift     = fgLowGainPeakShift     = 0
//
//////////////////////////////////////////////////////////////////////////////
#include "MExtractFixedWindowPeakSearch.h"
#include "MExtractor.h"

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

#include "MRawEvtPixelIter.h"

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

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

ClassImp(MExtractFixedWindowPeakSearch);

using namespace std;

const Byte_t MExtractFixedWindowPeakSearch::fgHiGainFirst          = 0;
const Byte_t MExtractFixedWindowPeakSearch::fgHiGainLast           = 14;
const Byte_t MExtractFixedWindowPeakSearch::fgLoGainFirst          = 3;
const Byte_t MExtractFixedWindowPeakSearch::fgLoGainLast           = 14;
const Byte_t MExtractFixedWindowPeakSearch::fgHiGainWindowSize     = 6;
const Byte_t MExtractFixedWindowPeakSearch::fgLoGainWindowSize     = 6;
const Byte_t MExtractFixedWindowPeakSearch::fgPeakSearchWindowSize = 4;
const Byte_t MExtractFixedWindowPeakSearch::fgOffsetFromWindow     = 1;
const Byte_t MExtractFixedWindowPeakSearch::fgLowGainPeakShift     = 0;
// --------------------------------------------------------------------------
//
// Default constructor. 
//
// Sets:
// - fWindowSizeHiGain to fgWindowSizeHiGain
// - fWindowSizeLoGain to fgWindowSizeLoGain
// - fPeakSearchWindowSize to fgPeakSearchWindowSize
// - fLowGainPeakShift to fgLowGainPeakShift
//
// Calls: 
// - SetOffsetFromWindow()
// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast)
//
MExtractFixedWindowPeakSearch::MExtractFixedWindowPeakSearch(const char *name, const char *title)
    : fWindowSizeHiGain(fgHiGainWindowSize), 
      fWindowSizeLoGain(fgLoGainWindowSize),
      fPeakSearchWindowSize(fgPeakSearchWindowSize),
      fLowGainPeakShift(fgLowGainPeakShift)
{

  fName  = name  ? name  : "MExtractFixedWindowPeakSearch";
  fTitle = title ? title : "Task to extract the signal from the FADC slices";

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

// --------------------------------------------------------------------------
//
// SetRange: 
//
// Calls:
// - MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
// - SetWindows(fWindowSizeHiGain,fWindowSizeLoGain,fPeakSearchWindowSize);
//
void MExtractFixedWindowPeakSearch::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
  //
  SetWindows(fWindowSizeHiGain,fWindowSizeLoGain,fPeakSearchWindowSize);
  
}

// --------------------------------------------------------------------------
//
// Checks:
// - if a window is odd, subtract one
// - 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 MExtractFixedWindowPeakSearch::SetWindows(Byte_t windowh, Byte_t windowl, Byte_t peaksearchwindow)
{

  fWindowSizeHiGain     = windowh & ~1;
  fWindowSizeLoGain     = windowl & ~1;
  fPeakSearchWindowSize = peaksearchwindow & ~1;

  if (fWindowSizeHiGain != windowh)
    *fLog << warn << GetDescriptor() << ": Hi Gain window size has to be even, set to: " 
          << int(fWindowSizeHiGain) << " samples " << endl;
  
  if (fWindowSizeLoGain != windowl)
    *fLog << warn << GetDescriptor() << ": Lo Gain window size has to be even, set to: " 
          << int(fWindowSizeLoGain) << " samples " << endl;
    
  if (fPeakSearchWindowSize != peaksearchwindow)
    *fLog << warn << GetDescriptor() << ": Peak Search window size has to be even, set to: " 
          << int(fPeakSearchWindowSize) << " samples " << endl;

  const Byte_t availhirange = (fHiGainLast-fHiGainFirst+1) & ~1;
  const Byte_t availlorange = (fLoGainLast-fLoGainFirst+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 (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 (fWindowSizeHiGain<2) 
    {
      fWindowSizeHiGain = 2;
      *fLog << warn << GetDescriptor() << ": Hi Gain window size set to two samples" << endl;
    }
  
  if (fWindowSizeLoGain<2) 
    {
      fWindowSizeLoGain = 2;
      *fLog << warn << GetDescriptor() << ": Lo Gain window size set to two samples" << endl;
    }

  if (fPeakSearchWindowSize<2) 
    {
      fPeakSearchWindowSize = 2;
      *fLog << warn << GetDescriptor() 
            << ": Peak Search window size set to two samples" << endl;
    }

  fNumHiGainSamples = (Float_t)fWindowSizeHiGain;
  fNumLoGainSamples = (Float_t)fWindowSizeLoGain;

  fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
  fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
}


// --------------------------------------------------------------------------
//
// The ReInit calls:
// -  MExtractor::ReInit()
// -  fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast+fHiLoLast, fNumHiGainSamples,
//                                fLoGainFirst, fLoGainLast, fNumLoGainSamples);
//
Bool_t MExtractFixedWindowPeakSearch::ReInit(MParList *pList)
{

  MExtractor::ReInit(pList);
  
  fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast+fHiLoLast, fNumHiGainSamples,
                              fLoGainFirst, fLoGainLast, fNumLoGainSamples);

  return kTRUE;
  
}

// --------------------------------------------------------------------------
//
// 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 MExtractFixedWindowPeakSearch::FindPeak(Byte_t *ptr, Byte_t window, Byte_t &startslice, Int_t &max, Int_t &sat) const
{

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

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

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

  while (p<ptr+window)
    {
      sum += *p;
      if (*p++ >= fSaturationLimit)
        sat++;
    }

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

  //
  // Calculate the i-th sum of n as
  //    sum_i+1 = sum_i + slice[i+window] - slice[i]
  // This is fast and accurate (because we are using int's)
  //
  max=sum;
  for (p=ptr; p+window<end; p++)
    {
      sum += *(p+window) - *p;
      if (sum > max)
        {
          max = sum;
          startslice = p-ptr;
        }
    }

  return;
}

// --------------------------------------------------------------------------
//
// FindSignalHiGain:
//
// - Loop from ptr to (ptr+fWindowSizeHiGain)
// - Sum up contents of *ptr
// - If *ptr is greater than fSaturationLimit, raise sat by 1
// 
void MExtractFixedWindowPeakSearch::FindSignalHiGain(Byte_t *ptr, Byte_t *logain, Float_t &sum, Byte_t &sat) const
{

  Byte_t *end = ptr + fWindowSizeHiGain-fHiLoLast;

  Int_t summ = 0;
  //
  // Calculate the sum of the "window" slices starting in ptr
  //
  while (ptr<end)
    {
      summ += *ptr;
      if (*ptr++ >= fSaturationLimit)
        sat++;
    }

  //
  // If part of the "low-Gain" slices are used, 
  // repeat steps one and two for the logain part until fHiLoLast
  //
  Byte_t *p = logain;
  end = logain + fHiLoLast;
  while (p<end)
    {
  
      summ += *p;
      if (*p++ >= fSaturationLimit)
        sat++;
    }

  sum = (Float_t)summ;
  return;
}

// --------------------------------------------------------------------------
//
// FindSignalLoGain:
//
// - Loop from ptr to (ptr+fWindowSizeLoGain)
// - Sum up contents of *ptr
// - If *ptr is greater than fSaturationLimit, raise sat by 1
// 
void MExtractFixedWindowPeakSearch::FindSignalLoGain(Byte_t *ptr, Float_t &sum, Byte_t &sat) const
{
  //
  // Calculate the sum of the "window" slices starting in ptr
  //
  Byte_t *p = ptr;
  Int_t summ = 0;

  while (p<ptr+fWindowSizeLoGain)
    {
      summ += *p;
      if (*p++ >= fSaturationLimit)
        sat++;
    }

  sum = (Float_t)summ;
  return;
}

// --------------------------------------------------------------------------
//
// Process
// First find the pixel with highest sum of fPeakSearchWindowSize slices (default:4)
// "startslice" will mark the slice at which the highest sum begins for that pixel.
// Then define the beginning of the integration window for ALL pixels as the slice
// before that: startslice-fOffsetFromWindow, unless of course startslice-fOffsetFromWindow<=0,
// in which case we start at 0. We will also check that the integration window does not 
// go beyond the FADC limits.
//
Int_t MExtractFixedWindowPeakSearch::Process()
{

  MRawEvtPixelIter pixel(fRawEvt);

  Int_t sat;
  Int_t maxsumhi = -1000000;
  Byte_t startslice;
  Byte_t hiGainFirst = 0;
  Byte_t loGainFirst = 0;
  Byte_t hilolastsave = fHiLoLast;

  while (pixel.Next())
    {
      Int_t sumhi;
      sat = 0;

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

      if (sumhi > maxsumhi && sat == 0 )
        {
          maxsumhi = sumhi;
          if ((startslice-fOffsetFromWindow) > 0)
            hiGainFirst = fHiGainFirst + startslice - fOffsetFromWindow;
          else
            hiGainFirst = fHiGainFirst;
        }
    }


  loGainFirst = ( hiGainFirst+fLowGainPeakShift > fLoGainFirst ) ? 
      hiGainFirst+fLowGainPeakShift : fLoGainFirst;

  // Make sure we will not integrate beyond the hi gain limit:
  if (hiGainFirst+fWindowSizeHiGain > pixel.GetNumHiGainSamples())
    fHiLoLast = hiGainFirst+fWindowSizeHiGain - pixel.GetNumHiGainSamples();
  //    hiGainFirst = pixel.GetNumHiGainSamples()-fWindowSizeHiGain;

  // Make sure we will not integrate beyond the lo gain limit:
  if (loGainFirst+fWindowSizeLoGain > pixel.GetNumLoGainSamples())
    loGainFirst = pixel.GetNumLoGainSamples()-fWindowSizeLoGain;

  pixel.Reset();
  fSignals->Clear();

  sat = 0;
  while (pixel.Next())
    {
      //
      // Find signal in hi- and lo-gain
      //
      Float_t sumhi=0.;
      Byte_t sathi=0;

      FindSignalHiGain(pixel.GetHiGainSamples()+hiGainFirst, pixel.GetLoGainSamples(), sumhi, sathi);

      Float_t sumlo=0.;
      Byte_t satlo=0;
      if (pixel.HasLoGain())
        {
          FindSignalLoGain(pixel.GetLoGainSamples()+loGainFirst, sumlo, satlo);
          if (satlo)
            sat++;
        }

      //
      // Take corresponding pedestal
      //
      const Int_t pixid = pixel.GetPixelId();

      const MPedestalPix  &ped = (*fPedestals)[pixid];
      MExtractedSignalPix &pix = (*fSignals)[pixid];

      const Float_t pedes  = ped.GetPedestal();
      const Float_t pedrms = ped.GetPedestalRms();

      //
      // Set extracted signal with pedestal substracted
      //
      pix.SetExtractedSignal(sumhi - pedes*fNumHiGainSamples, pedrms*fSqrtHiGainSamples,
                             sumlo - pedes*fNumLoGainSamples, pedrms*fSqrtLoGainSamples);

      pix.SetGainSaturation(sathi, sathi, satlo);

      //      pix.SetNumHiGainSlices(fNumHiGainSamples);
      //      pix.SetNumLoGainSlices(fNumLoGainSamples);      
      
    } /* while (pixel.Next()) */


  fHiLoLast = hilolastsave;
  fSignals->SetReadyToSave();

  return kTRUE;
}