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

/* ======================================================================== *\
   !
   ! *
   ! * 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, 4/2004 <mailto:moralejo@pd.infn.it>
   !              Markus Gaug      , 4/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 MExtractSignal2.
//
//////////////////////////////////////////////////////////////////////////////
#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::fgHiGainWindowSize     = 6;
const Byte_t MExtractFixedWindowPeakSearch::fgLoGainWindowSize     = 6;
const Byte_t MExtractFixedWindowPeakSearch::fgPeakSearchWindowSize = 4;
const Byte_t MExtractFixedWindowPeakSearch::fgHiGainFirst          = 0;
const Byte_t MExtractFixedWindowPeakSearch::fgHiGainLast           = 14;
const Byte_t MExtractFixedWindowPeakSearch::fgLoGainFirst          = 0;
const Byte_t MExtractFixedWindowPeakSearch::fgLoGainLast           = 14;
// --------------------------------------------------------------------------
//
// Default constructor. 
//
MExtractFixedWindowPeakSearch::MExtractFixedWindowPeakSearch(const char *name, const char *title)
{

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

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

void MExtractFixedWindowPeakSearch::SetWindows(Byte_t windowh, Byte_t windowl, Byte_t peaksearchwindow)
{
  //
  // Set windows to even number of slices due to clock noise (odd/even slice effect).
  //
  fWindowSizeHiGain = windowh & ~1;
  fWindowSizeLoGain = windowl & ~1;
  fPeakSearchWindowSize = peaksearchwindow & ~1;


  if (fWindowSizeHiGain != windowh)
    *fLog << endl << warn << 
      "MExtractFixedWindowPeakSearch::SetWindows - Hi Gain window size has to be even, set to: " 
          << int(fWindowSizeHiGain) << " samples " << endl;
    
  if (fWindowSizeLoGain != windowl)
    *fLog << endl << warn << 
      "MExtractFixedWindowPeakSearch::SetWindows - Lo Gain window size has to be even, set to: " 
          << int(fWindowSizeLoGain) << " samples " << endl;

  if (fPeakSearchWindowSize != peaksearchwindow)
    *fLog << endl << warn << 
      "MExtractFixedWindowPeakSearch::SetWindows - Peak Search window size has to be even, set to: " 
          << int(fPeakSearchWindowSize) << " samples " << endl;


  if (fWindowSizeHiGain<2) 
    {
      fWindowSizeHiGain = 2;
      *fLog << warn 
            << "MExtractFixedWindowPeakSearch::SetWindows - Hi Gain window size set to two samples" << endl;
    }

  if (fWindowSizeLoGain<2) 
    {
      fWindowSizeLoGain = 2;
      *fLog << warn 
            << "MExtractFixedWindowPeakSearch::SetWindows - Lo Gain window size set to two samples" << endl;
    }

  if (fPeakSearchWindowSize<2) 
    {
      fPeakSearchWindowSize = 2;
      *fLog << warn 
            << "MExtractFixedWindowPeakSearch::SetWindows - 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);
}

// --------------------------------------------------------------------------
//
// 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;

  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;
}

// --------------------------------------------------------------------------
//
// FindSignal
// Adds up "window" slices starting in slice to which "ptr" points. The result 
// is stored in "sum", and the number of saturated  slices in "sat".
//
void MExtractFixedWindowPeakSearch::FindSignalHiGain(Byte_t *ptr, Int_t &sum, Byte_t &sat) const
{
  //
  // Calculate the sum of the "window" slices starting in ptr
  //
  sum = 0;
  sat = 0;
  Byte_t *p = ptr;

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

  return;
}

// --------------------------------------------------------------------------
//
// FindSignal
// Adds up "window" slices starting in slice to which "ptr" points. The result 
// is stored in "sum", and the number of saturated  slices in "sat".
//
void MExtractFixedWindowPeakSearch::FindSignalLoGain(Byte_t *ptr, Int_t &sum, Byte_t &sat) const
{
  //
  // Calculate the sum of the "window" slices starting in ptr
  //
  sum = 0;
  sat = 0;
  Byte_t *p = ptr;

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

  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-1 (this is somehow arbitrary), unless of course startslice=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 = 0;
  Byte_t startslice;
  Byte_t hiGainFirst = 0;
  Byte_t loGainFirst = 0;

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

      FindPeak(pixel.GetHiGainSamples(), fPeakSearchWindowSize, startslice, sumhi, sat);
      if (sumhi > maxsumhi && sat == 0 )
        {
          maxsumhi = sumhi;
          if (startslice > 0)
            hiGainFirst = startslice-1;
          else
            hiGainFirst = 0;
        }
    }


  loGainFirst = hiGainFirst;

  // Make sure we will not integrate beyond the hi gain limit:
  if (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
      //
      Int_t sumhi;
      Byte_t sathi;

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

      Int_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);
    } /* while (pixel.Next()) */

  fSignals->SetReadyToSave();

  return kTRUE;
}