source: trunk/MagicSoft/Mars/mextralgo/MExtralgoDigitalFilter.cc@ 8049

/* ======================================================================== *\
!
! *
! * 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): Thomas Bretz, 08/2006 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2006
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//   MExtralgoDigitalFilter
//
//////////////////////////////////////////////////////////////////////////////
#include "MExtralgoDigitalFilter.h"

using namespace std;

Float_t MExtralgoDigitalFilter::ExtractNoise(Int_t iter, Int_t windowsize) const
{
    Double_t sum, dummy;
    Eval(sum, dummy, windowsize, 0, iter-fWeightsPerBin/2);
    return sum;
}

void MExtralgoDigitalFilter::Extract(Int_t windowsize, Float_t timeshift)
{
    fSignal    =  0; // default is: no pulse found
    fTime      = -1; // default is: out if range (--> random)
    fSignalDev =  0; // default is: valid
    fTimeDev   =  0; // default is: valid

    // FIXME: How to handle saturation?

    Float_t maxamp  = -FLT_MAX;
    Float_t maxtime =  0;
    Int_t   maxp    = -1;

    //
    // Calculate the sum of the first fWindowSize slices
    //
    for (Int_t i=0; i<fNum-windowsize+1; i++)
    {
        Double_t sumamp, sumtime;
        Eval(sumamp, sumtime, windowsize, i);

        if (sumamp>maxamp)
        {
            maxamp  = sumamp;
            maxtime = sumtime;
            maxp    = i;
        }
    }

    // The number of available slices were smaller than the
    // extraction window size of the extractor
    if (maxp<0)
    {
        fSignalDev = -1;  // means: is invalid
        fTimeDev   = -1;  // means: is invalid
        return;
    }

    // For some reason (by chance or because all slices contain 0)
    // maxamp is 0. This means the signal is zero and no arrival
    // time can be extracted (but both informations are valid)
    if (maxamp==0)
        return;

    // This is the time offset from the extraction position
    const Float_t time = maxtime/maxamp;

    // This is used to align the weights to bins between
    // -0.5/fWeightsPerBin and 0.5/fWeightsPerBin instead of
    // 0 and 1./fWeightsPerBin
    const Float_t halfres = 0.5/fWeightsPerBin;

    // move the extractor by an offset number of slices
    // determined by the extracted time
    const Int_t integ = TMath::FloorNint(time+halfres+0.5);
    maxp -= integ;

    // Convert the residual fraction of one slice into an
    // offset position in the extraction weights
    Int_t frac = TMath::FloorNint((time+halfres-integ)*fWeightsPerBin);

    // Align maxp into available range
    if (maxp < 0)
    {
        maxp = 0;
        frac = -fWeightsPerBin/2; // Assume peak at the beginning of the first slice
    }
    if (maxp+windowsize > fNum)
    {
        maxp = fNum-windowsize;
        frac = fWeightsPerBin/2-1; // Assume peak at the end of the last slice
    }

    Double_t sumamp, sumtime;
    Eval(sumamp, sumtime, windowsize, maxp, frac);

    fSignal = sumamp;
    if (sumamp == 0)
        return;

    // here, the first high-gain slice is already included
    // in the fTimeShiftHiGain this shifts the time to the
    // start of the rising edge
    fTime = maxp - Float_t(frac)/fWeightsPerBin + timeshift;

    const Float_t timefineadjust = sumtime/sumamp;

    // FIXME: Is 3 a good value?
    if (TMath::Abs(timefineadjust)*fWeightsPerBin < 3.)
        fTime -= timefineadjust;
}