source: trunk/MagicSoft/Mars/msignal/MExtractSlidingWindow.cc@ 3869

/* ======================================================================== *\
!
! *
! * 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): Thomas Bretz, 02/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
!              Hendrik Bartko, 01/2004 <mailto:hbartko@mppmu.mpg.de>
!              Markus Gaug   , 04/2004 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//   MExtractSlidingWindow
//
//  Extracts the signal from a sliding window of size fHiGainWindowSize and 
//  fLoGainWindowSize. The signal is the one which maximizes the integral 
//  contents. 
//
//////////////////////////////////////////////////////////////////////////////
#include "MExtractSlidingWindow.h"
#include "MExtractor.h"

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

ClassImp(MExtractSlidingWindow);

using namespace std;

const Byte_t MExtractSlidingWindow::fgHiGainFirst      = 3;
const Byte_t MExtractSlidingWindow::fgHiGainLast       = 14;
const Byte_t MExtractSlidingWindow::fgLoGainFirst      = 3;
const Byte_t MExtractSlidingWindow::fgLoGainLast       = 14;
const Byte_t MExtractSlidingWindow::fgHiGainWindowSize = 6;
const Byte_t MExtractSlidingWindow::fgLoGainWindowSize = 6;
// --------------------------------------------------------------------------
//
// Default constructor. 
//
MExtractSlidingWindow::MExtractSlidingWindow(const char *name, const char *title)
{

  fName  = name  ? name  : "MExtractSlidingWindow";
  fTitle = title ? title : "Signal Extractor for a sliding FADC window";

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

}

void MExtractSlidingWindow::SetWindowSize(Byte_t windowh, Byte_t windowl)
{
  
  fWindowSizeHiGain = windowh & ~1;
  fWindowSizeLoGain = windowl & ~1;

  if (fWindowSizeHiGain != windowh)
    *fLog << warn << "MExtractSignal2::SetRange - Hi Gain window size has to be even, set to: " 
          << int(fWindowSizeHiGain) << " samples " << endl;
  
    if (fWindowSizeLoGain != windowl)
      *fLog << warn << "MExtractSignal2::SetRange - Lo Gain window size has to be even, set to: " 
            << int(fWindowSizeLoGain) << " samples " << endl;
    
    if (fWindowSizeHiGain<2) 
      {
        fWindowSizeHiGain = 2;
      *fLog << warn << "MExtractSignal2::SetRange - Hi Gain window size set to two samples" << endl;
      }
    
    if (fWindowSizeLoGain<2) 
    {
      fWindowSizeLoGain = 2;
      *fLog << warn << "MExtractSignal2::SetRange - Lo Gain window size set to two samples" << endl;
    }
    
    fNumHiGainSamples = (Float_t)fWindowSizeHiGain;
    fNumLoGainSamples = (Float_t)fWindowSizeLoGain;
    
    fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
    fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
}



void MExtractSlidingWindow::FindSignalHiGain(Byte_t *ptr, Int_t &max, Byte_t &sat) const
{
    const Byte_t *end = ptr + fHiGainLast - fHiGainFirst + 1;

    Int_t sum=0;

    //
    // Calculate the sum of the first fWindowSize slices
    //
    sat = 0;
    Byte_t *p = ptr;
    while (p<ptr+fWindowSizeHiGain)
    {
        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 as
    //    sum_i+1 = sum_i + slice[i+8] - slice[i]
    // This is fast and accurate (because we are using int's)
    //
    max=sum;
    for (p=ptr; p+fWindowSizeHiGain<end; p++)
    {
        sum += *(p+fWindowSizeHiGain) - *p;
        if (sum>max)
            max = sum;
    }
}


void MExtractSlidingWindow::FindSignalLoGain(Byte_t *ptr, Int_t &max, Byte_t &sat) const
{
    const Byte_t *end = ptr + fLoGainLast - fLoGainFirst + 1;

    Int_t sum=0;

    //
    // Calculate the sum of the first fWindowSize slices
    //
    sat = 0;
    Byte_t *p = ptr;
    while (p<ptr+fWindowSizeLoGain)
    {
        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 as
    //    sum_i+1 = sum_i + slice[i+8] - slice[i]
    // This is fast and accurate (because we are using int's)
    //
    max=sum;
    for (p=ptr; p+fWindowSizeLoGain<end; p++)
    {
        sum += *(p+fWindowSizeLoGain) - *p;
        if (sum>max)
            max = sum;
    }
}