source: trunk/MagicSoft/Mars/msignal/MExtractBlindPixel.cc@ 3614

/* ======================================================================== *\
!
! *
! * 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, 02/2004 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//   MExtractBlindPixel
//
//////////////////////////////////////////////////////////////////////////////
#include "MExtractBlindPixel.h"

#include <fstream>

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

#include "MParList.h"

#include "MRawEvtData.h"
#include "MRawEvtPixelIter.h"

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

#include "MExtractedSignalBlindPixel.h"

ClassImp(MExtractBlindPixel);

using namespace std;

const Int_t  MExtractBlindPixel::fgBlindPixelIdx   = 559;
const Byte_t MExtractBlindPixel::fgSaturationLimit = 254;
const Byte_t MExtractBlindPixel::fgFirst =  3;
const Byte_t MExtractBlindPixel::fgLast  = 16;

// --------------------------------------------------------------------------
//
// Default constructor. 
//
MExtractBlindPixel::MExtractBlindPixel(const char *name, const char *title)
    : fSaturationLimit(fgSaturationLimit)
{
  
  fName  = name  ? name  : "MExtractBlindPixel";
  fTitle = title ? title : "Task to extract the signal from the FADC slices";
  
  AddToBranchList("MRawEvtData.*");
  
  SetBlindPixelIdx();
  SetSaturationLimit();
  SetRange();
}

void MExtractBlindPixel::SetRange(Byte_t first, Byte_t last)
{

    fNumSamples = last-first+1;
    fFirst      = first;
    fLast       = last;

    fSqrtSamples = TMath::Sqrt((Float_t)fNumSamples);
}

// --------------------------------------------------------------------------
//
// The PreProcess searches for the following input containers:
//  - MRawEvtData
//
// The following output containers are also searched and created if
// they were not found:
//
//  - MExtractedBlindPixel
//
Int_t MExtractBlindPixel::PreProcess(MParList *pList)
{
    fRawEvt = (MRawEvtData*)pList->FindObject(AddSerialNumber("MRawEvtData"));
    if (!fRawEvt)
    {
        *fLog << err << AddSerialNumber("MRawEvtData") << " not found... aborting." << endl;
        return kFALSE;
    }

    fPedestals = (MPedestalCam*)pList->FindObject(AddSerialNumber("MPedestalCam"));
    if (!fPedestals)
    {
        *fLog << err << AddSerialNumber("MPedestalCam") << " not found... aborting." << endl;
        return kFALSE;
    }

    fBlindPixel = (MExtractedSignalBlindPixel*)pList->FindCreateObj(AddSerialNumber("MExtractedSignalBlindPixel"));
    if (!fBlindPixel)
        return kFALSE;

    fBlindPixel->SetUsedFADCSlices(fFirst, fLast);
    fBlindPixel->SetBlindPixelIdx(fBlindPixelIdx);

    return kTRUE;
}

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

    MRawEvtPixelIter pixel(fRawEvt);

    fBlindPixel->Clear();

    pixel.Jump(fBlindPixelIdx);
 
    const UInt_t nhigain = pixel.GetNumHiGainSamples();

    Byte_t *ptr = pixel.GetHiGainSamples();
    
    //
    // We need a dedicated signal extractor for the blind pixel
    //
    Int_t  diff  = 0;
    UInt_t first = fFirst;
    UInt_t last  = fLast;
    UInt_t sat  = 0;

    if (last > nhigain)
      {
        diff = last - nhigain;
        last = nhigain;
      }
    
    
    Byte_t *start   = ptr + first - 1;
    Byte_t *end     = ptr + last  - 1;
    
    ptr = start;
    
    Int_t sum = 0;
    
    while (ptr<=end)
      {
        sum += *ptr;

        if (*ptr++ >= fSaturationLimit)
          sat++;
      }
    
    if (diff > 0)
      {
        ptr = pixel.GetLoGainSamples();
        end = ptr + diff - 1;
        
        while (ptr<=end)
          {

            sum += *ptr;

            if (*ptr++ >= fSaturationLimit)
              sat++;

          }
      }
    
    fBlindPixel->SetExtractedSignal(sum);
    fBlindPixel->SetNumSaturated(sat);

    if (sat)
      *fLog << warn << "WARNING - saturation occurred in the Blind Pixel " << endl;


    fBlindPixel->SetReadyToSave();

    return kTRUE;
}

Int_t MExtractBlindPixel::PostProcess()
{

    MPedestalPix &pedpix  = (*fPedestals)[fBlindPixelIdx];    

    if (&pedpix)
    {
         fBlindPixel->SetPed      ( pedpix.GetPedestal()   * fNumSamples );
         fBlindPixel->SetPedErr   ( pedpix.GetPedestalRms()* fNumSamples / TMath::Sqrt((Float_t)fPedestals->GetTotalEntries()) );
         fBlindPixel->SetPedRms   ( pedpix.GetPedestalRms()* TMath::Sqrt((Float_t)fNumSamples) );
         fBlindPixel->SetPedRmsErr( fBlindPixel->GetPedErr()/2. );
    }

    return kTRUE;

}

// --------------------------------------------------------------------------
//
// Implementation of SavePrimitive. Used to write the call to a constructor
// to a macro. In the original root implementation it is used to write
// gui elements to a macro-file.
//
void MExtractBlindPixel::StreamPrimitive(ofstream &out) const
{
    out << "   " << ClassName() << " " << GetUniqueName() << "(\"";
    out << "\"" << fName << "\", \"" << fTitle << "\");" << endl;

    if (fSaturationLimit!=fgSaturationLimit)
    {
        out << "   " << GetUniqueName() << ".SetSaturationLimit(";
        out << (int)fSaturationLimit << ");" << endl;
    }

    const Bool_t arg2 = fNumSamples+fFirst-1 != fgLast;
    const Bool_t arg1 = arg2 || fFirst != fgFirst;

    if (!arg1)
        return;

    out << "   " << GetUniqueName() << ".SetRange(";
    out << (int)fFirst;
    if (arg2)
      out << ", " << (int)(fNumSamples+fFirst-1);
    out << ");" << endl;
}