source: trunk/MagicSoft/Mars/manalysis/MExtractSignal.cc@ 2640

/* ======================================================================== *\
!
! *
! * 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): Markus Gaug  09/2003 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2001
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//                                                                          //
//   MExtractSignal                                                         //
//                                                                          //
//////////////////////////////////////////////////////////////////////////////

#include "MExtractSignal.h"

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

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

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

#include "MParList.h"
#include "MH.h"

#include "MRawRunHeader.h"
#include "MRawEvtData.h"       // MRawEvtData::GetNumPixels
#include "MRawEvtPixelIter.h"

#include "TMath.h"

ClassImp(MExtractSignal);

using namespace std;
// --------------------------------------------------------------------------
//
// Default constructor. 
//
MExtractSignal::MExtractSignal(const char *name, const char *title)
  : fBefore(2), fAfter(4), fSaturationLimit(254), fConversionHiLo(10.)
{

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

    AddToBranchList("MRawEvtData.*");
}

// --------------------------------------------------------------------------
//
// The PreProcess searches for the following input containers:
//  - MRawEvtData
//  - MPedestalCam
//
// The following output containers are also searched and created if
// they were not found:
//
//  - MExtractedSignalCam
//
Int_t MExtractSignal::PreProcess(MParList *pList)
{

    fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
    if (!fRawEvt)
    {
        *fLog << dbginf << "MRawEvtData not found... aborting." << endl;
        return kFALSE;
    }

    const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
    if (!runheader)
        *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl;


    fSignals = (MExtractedSignalCam*)pList->FindCreateObj("MExtractedSignalCam");
    if (!fSignals)
        return kFALSE;

    fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");
    if (!fPedestals)
      {
        *fLog << err << dbginf << "Cannot find MPedestalCam ... aborting" << endl;
        return kFALSE;
      }
    
    return kTRUE;
}


// --------------------------------------------------------------------------
//
// The ReInit searches for the following input containers:
//  - MRawRunHeader
//
Bool_t MExtractSignal::ReInit(MParList *pList )
{
 
    fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
    if (!fRunHeader)
    {
        *fLog << dbginf << "MRawRunHeader not found... aborting." << endl;
        return kFALSE;
    }

    //fNumHiGainSamples =  fRunHeader->GetNumSamplesHiGain();
    //fNumLoGainSamples =  fRunHeader->GetNumSamplesLoGain();
    fNumHiGainSamples =  fAfter+fBefore;
    fNumLoGainSamples =  fAfter+fBefore;

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

    //
    // FIXME: The next statement simply does not work: 
    //        fRawEvt->GetNumPixels() returns always 0
    //

    fSignals->InitSize(577);    
    // fExtractedSignals->InitSize(fRawEvt->GetNumPixels());

    return kTRUE;
}


// --------------------------------------------------------------------------
//
// Calculate the integral of the FADC time slices and store them as a new
// pixel in the MExtractedSignalCam container.
//
Int_t MExtractSignal::Process()
{

    MRawEvtPixelIter pixel(fRawEvt);
    fSignals->Clear();

    while (pixel.Next())
      {

        UShort_t satHi = 0;
        UShort_t satLo = 0;

        Byte_t *ptr   = pixel.GetHiGainSamples();
        Byte_t *max   = ptr + pixel.GetIdxMaxHiGainSample();
        Byte_t *end   = ptr + fRawEvt->GetNumHiGainSamples();
        Byte_t *first;
        Byte_t *last;

        UInt_t  sumHi   = 0;
        UInt_t  sumLo   = 0;

        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();

        if (max-3 < ptr)
          {
            first = ptr+3;
            last  = ptr+3+fAfter+fBefore;
          }
        else if (max+fAfter-1 > end)
          {
            first = end-fAfter-fBefore-1;
            last  = end-1;
          }
        else
          {
            first = max-fBefore;
            last  = max+fAfter;
          }

        for (ptr=first;ptr<last;ptr++)
          {
            sumHi += *ptr;
            if (*ptr >= fSaturationLimit)
              satHi++;
          }

        ptr   = pixel.GetLoGainSamples();
        max   = ptr+pixel.GetIdxMaxLoGainSample();
        end   = ptr+fRawEvt->GetNumLoGainSamples();

        if (max-4 < ptr)
          {
            first = ptr+4;
            last  = ptr+4+fAfter+fBefore;
          }
        else if (max+fAfter > end)
          {
            first = end-fAfter-fBefore;
            last  = end;
          }
        else
          {
            first = max-fBefore;
            last  = max+fAfter;
          }

        for (ptr=first;ptr<last;ptr++)
          {
            sumLo += *ptr;
            if (*ptr >= fSaturationLimit)
              {
                *fLog << err << dbginf 
                      << "Warning: Saturation of Lo Gain reached in pixel: " 
                      << pixid << " " << "   sum = " << sumLo << endl;
                satLo++;
              }
          }

        //
        // FIXME: This is preliminary, we will change to pedestals per slice!!!
        // Assume pedestals per time slice ==> multiply with number of slices
        //

        pix.SetExtractedSignal((float)sumHi - pedes*(float)fNumHiGainSamples,
                               pedrms*fSqrtHiGainSamples,
                               ((float)sumLo - pedes*(float)fNumLoGainSamples)*fConversionHiLo,
                               pedrms*fSqrtLoGainSamples
                               );

        pix.SetGainSaturation((satHi),satHi,satLo);

      } /* while (pixel.Next()) */


    fSignals->SetNumUsedFADCSlices(fNumHiGainSamples);
    fSignals->SetReadyToSave();

    return kTRUE;
}

Int_t MExtractSignal::PostProcess()
{

  return kTRUE;

}