source: trunk/MagicSoft/Mars/mraw/MRawEvtPixelIter.cc@ 8106

/* ======================================================================== *\
!
! *
! * 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  12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
!   Author(s): Markus Gaus 10/2002 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

///////////////////////////////////////////////////////////////////////////////
//
//  MRawEvtPixelIter
//
//  class to iterate over all pixels of one event.
//  The calling is similar to a root iterator:
//  
//  MRawEvtData *evtdata;            // must be filled with data from somewhere
//  MRawEvtPixelIter pixel(evtdata); // evtdata: ptr to event you want to iterate
//
//  while (pixel.Next())
//  {
//     // here you can access the actual time slices by using
//     //   pixel.GetPixelId();
//     //   pixel.GetHiGainFadcSamples()[i]; // i is the number of the slice
//     //   pixel.HasLoGain();               // check if pixel has
//     //   pixel.GetLoGainFadcSamples()[i]; // i is the number of the slice
//
//     // WARNING: Don't acces more time slices than available.
//     //   Get the numbers by calling: evtdata->GetNum[Lo,Hi]GainSamples()
//     //   This number is constant for one event
//  }
//
///////////////////////////////////////////////////////////////////////////////
#include "MRawEvtPixelIter.h"

#include <TArrayC.h>

#include "MRawEvtData.h"

#include "MArrayS.h"
#include "MArrayB.h"

ClassImp(MRawEvtPixelIter);

using namespace std;

MRawEvtPixelIter::MRawEvtPixelIter(MRawEvtData *dat) : fABFlags(0), fData(dat)
{
    fNumHiGainSamples = dat->GetNumHiGainSamples();
    fNumLoGainSamples = dat->GetNumLoGainSamples();

    Reset();
}

// --------------------------------------------------------------------------
//
// return kTRUE  the lo gains exist for the actual pixel, else return kFALSE
//
Bool_t MRawEvtPixelIter::HasLoGain() const
{
    // We have no lo-gain at all
    if (!fLoGainId)
        return kFALSE;

    // This is to make the first call of this function in Next()
    // work properly! NEVER call this function before Next()
    if (fNumLoGainEntry==0)
        return kTRUE;

    if (fNumLoGainEntry>fData->fLoGainPixId->GetSize())
        return kFALSE;

    Bool_t rc = *fHiGainId!=*fLoGainId;

    if (rc)
        return kFALSE;

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// It steps to the next pixel. If there is no next pixel NULL is returned.
// If a next pixel where found, a pointer to the primary given (constructor)
// data structur is returned.
//
MRawEvtData *MRawEvtPixelIter::Next()
{
    //
    // if we are already at the last entry there is no 'next' entry anymore
    //
    if (fNumHiGainEntry==fData->fHiGainPixId->GetSize())
        return NULL;

    //
    // if we are already at the last entry there is no 'next' entry anymore
    //
    if (HasLoGain())
    {
        //
        // if higainpixid and logainpixid of the actual pixel are
        // identical then we have to move the pointer to the next
        // entry in the lo gains
        //
        fNumLoGainEntry++;
        fLoGainId++;
        fLoGainPos += fNumLoGainSamples;
    }

    //
    // here we have to move the pointer to the next entry in the hi gains
    //
    fNumHiGainEntry++;
    fHiGainId++;
    fHiGainPos += fNumHiGainSamples;

    //
    // return a pointer to the 'source' class if we succeed
    //
    return fData;
}

// --------------------------------------------------------------------------
//
// Reset the iteration. Jump to the first pixel.
//
void MRawEvtPixelIter::Reset()
{
    //
    // set counter to zero
    //
    fNumLoGainEntry = 0;
    fNumHiGainEntry = 0;

    //
    // set pointer to first entry of arrays
    //
    fHiGainId   = fData->fHiGainPixId->GetArray()-1;
    fLoGainId   = fData->fLoGainPixId->GetArray();
    fHiGainPos  = fData->fHiGainFadcSamples->GetArray()-fNumHiGainSamples;
    fLoGainPos  = fData->fLoGainFadcSamples->GetArray()-fNumLoGainSamples;
    fABFlags    = fData->fABFlags->GetSize()==0 ? 0 : fData->fABFlags->GetArray();

    //
    // In case fLoGainPixId.GetSize()=0 some root versions seems to
    // initialize the array with NULL. This makes both cases work.
    //
    if (fLoGainId)
        fLoGainId -= 1;
}

// --------------------------------------------------------------------------
//
// Calls the draw-function of the actual pixel (see MRawEvtData::Draw)
//
void MRawEvtPixelIter::Draw(Option_t *t)
{ 
    fData->Draw(Form("%s%d", t, *fHiGainId));
}

// --------------------------------------------------------------------------
//
// returns the sum of all hi gain fadc samples of the actual pixel
//
ULong_t MRawEvtPixelIter::GetSumHiGainSamples() const
{
    //
    // return the sum of the hi gain samples of the present pixel
    //
    Byte_t *ptr = fHiGainPos;
    const Byte_t *end = fHiGainPos + fNumHiGainSamples;

    ULong_t sum=0;

    do sum += *ptr++;
    while (ptr != end);

    return sum;
}

// --------------------------------------------------------------------------
//
// returns the sum of squares of all hi gain fadc sample of the actual pixel
//
ULong_t MRawEvtPixelIter::GetSumSqrHiGainSamples() const
{
    //
    // return the sum of the squares of the hi gain samples of the present pixel
    //
    Byte_t *ptr = fHiGainPos;
    const Byte_t *end = fHiGainPos + fNumHiGainSamples;

    ULong_t sum=0;

    do sum += (*ptr)*(*ptr);
    while (++ptr != end);

    return sum;
}

// --------------------------------------------------------------------------
//
// Returns the variance (sigma^2) of the HiGainSamples
//
Float_t MRawEvtPixelIter::GetVarHiGainSamples() const
{
    Byte_t *ptr = fHiGainPos;
    const Byte_t *end = fHiGainPos + fNumHiGainSamples;

    ULong_t sum=0;
    ULong_t sqsum=0;

    do {
        sum += *ptr;
        sqsum += (*ptr)*(*ptr);
    } while (++ptr != end);

    return (sqsum-(Float_t)sum*sum/fNumHiGainSamples)/(fNumHiGainSamples-1);
}

// --------------------------------------------------------------------------
//
// Returns the index of the FADC slice the maximum signal in. If the highest
// slices have the same value the last one is returned.
//
Byte_t MRawEvtPixelIter::GetIdxMaxHiGainSample() const
{
    Byte_t *ptr = fHiGainPos+1;
    Byte_t *max = fHiGainPos;
    const Byte_t *end = fHiGainPos + fNumHiGainSamples;

    do if (*ptr>*max) max = ptr;
    while (++ptr != end);

    return max-fHiGainPos;
}

// --------------------------------------------------------------------------
//
// Returns the index of the FADC slice the maximum signal in. If no lo-gains
// are available -1 is returned. If the highest slices have the same value the
// last one is returned.
//
Short_t MRawEvtPixelIter::GetIdxMaxLoGainSample(const Byte_t lofirst) const
{
    if (!HasLoGain())
        return -1; // means: not found

    Byte_t *ptr = fLoGainPos+lofirst+1;
    Byte_t *max = fLoGainPos+lofirst;
    const Byte_t *end = fLoGainPos + fNumLoGainSamples;

    do if (*ptr>*max) max = ptr;
    while (++ptr != end);

    return max-fLoGainPos;
}

// --------------------------------------------------------------------------
//
// Returns the index of the maximum FADC slice from high gain at first. If 
// high gain is saturated it returns the low gain one.
// If no lo-gains are existing and the hi-gains have saturating values
// a negative value (-1) is returned.
// 
Short_t MRawEvtPixelIter::GetIdxMaxHiLoGainSample() const
{
    Byte_t max  = 0;
    Char_t maxi = 0;

    for (int i=fNumHiGainSamples-1; i>=0; i--)
        if (fHiGainPos[i]>max)
        {
            max  = fHiGainPos[i];
            maxi = i;
        }

    return max<0xff ? maxi : GetIdxMaxLoGainSample();
}

// --------------------------------------------------------------------------
//
// Returns the maximum signal of all sliced in the hi gain samples
//
Byte_t MRawEvtPixelIter::GetMaxHiGainSample() const
{
    Byte_t max = 0;

    for (int i=0; i<fNumHiGainSamples; i++)
        if (fHiGainPos[i]>max)
            max = fHiGainPos[i];

    return max;
}

// --------------------------------------------------------------------------
//
// Returns the maximum signal of all sliced in the hi gain samples
//
Byte_t MRawEvtPixelIter::GetMaxLoGainSample() const
{
    Byte_t max = 0;

    for (int i=fNumLoGainSamples-1; i>=0; i--)
        if (fLoGainPos[i]>max)
            max = fLoGainPos[i];

    return max;
}

// --------------------------------------------------------------------------
//
// returns the sum of all lo gain fadc samples of the actual pixel.
// if no lo gain information is available 0 is returned.
//
ULong_t MRawEvtPixelIter::GetSumLoGainSamples() const
{
    //
    // return the sum of the lo gain samples of the present pixel
    //
    if (!HasLoGain())
        return 0;

    Byte_t *ptr = fLoGainPos;
    const Byte_t *end = fLoGainPos + fNumLoGainSamples;

    ULong_t sum=0;

    do sum += *ptr++;
    while (ptr != end);

    return sum;
}

// --------------------------------------------------------------------------
//
// returns the sum of squares of all hi gain fadc sample of the actual pixel
//
ULong_t MRawEvtPixelIter::GetSumSqrLoGainSamples() const
{
    //
    // return the sum of the lo gain samples squares of the present pixel
    //
    if (!HasLoGain())
        return 0;

    Byte_t *ptr = fLoGainPos;
    const Byte_t *end = fLoGainPos + fNumLoGainSamples;

    ULong_t sum=0;

    do sum += (*ptr)*(*ptr);
    while (++ptr != end);

    return sum;
}