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

/* ======================================================================== *\
!
! *
! * 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
//
//  Extracts the signal from a fixed window in a given range.
//
//  Call: SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast) 
//  to modify the ranges. The "lo-gain" ranges are used for the NSB rejection 
//  whereas the high-gain ranges for blind pixel signal extraction. "High-gain" 
//  ranges can extend to the slices stored as "low-gain" in MRawEvtPixelIter
//  Defaults are: 
// 
//   fHiGainFirst =  fgHiGainFirst =  12 
//   fHiGainLast  =  fgHiGainLast  =  16
//   fLoGainFirst =  fgLoGainFirst =  0 
//   fLoGainLast  =  fgLoGainLast  =  10
//
//////////////////////////////////////////////////////////////////////////////
#include "MExtractBlindPixel.h"

#include <fstream>

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

#include "MParList.h"

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

#include "MExtractedSignalBlindPixel.h"

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

ClassImp(MExtractBlindPixel);

using namespace std;

const Int_t  MExtractBlindPixel::fgBlindPixelIdx  = 559;
const Int_t  MExtractBlindPixel::fgNSBFilterLimit = 800;
const Byte_t MExtractBlindPixel::fgHiGainFirst    =  12;
const Byte_t MExtractBlindPixel::fgHiGainLast     =  29;
const Byte_t MExtractBlindPixel::fgLoGainFirst    =  0;
const Byte_t MExtractBlindPixel::fgLoGainLast     =  10;
// --------------------------------------------------------------------------
//
// Default constructor. 
//
// Initializes:
// - fBlindPixelIdx to fgBlindPixelIdx
// - fNSBFilterLimit to fgNSBFilterLimit
//
// Calls:
// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
//
MExtractBlindPixel::MExtractBlindPixel(const char *name, const char *title)
{
  
  fName  = name  ? name  : "MExtractBlindPixel";
  fTitle = title ? title : "Task to extract the signal from the FADC slices";
  
  AddToBranchList("MRawEvtData.*");
  
  SetBlindPixelIdx();
  SetNSBFilterLimit();
  SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);

}

void MExtractBlindPixel::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast)
{

  MExtractor::SetRange(hifirst,hilast,lofirst,lolast);

  fNumHiGainSamples = (Float_t)(fHiGainLast-fHiGainFirst+1);
  fNumLoGainSamples = (Float_t)(fLoGainLast-fLoGainFirst+1);  

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

// --------------------------------------------------------------------------
//
// 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)
{

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

  fBlindPixel->SetBlindPixelIdx(fBlindPixelIdx);
  fBlindPixel->SetUsedFADCSlices(fHiGainFirst, fHiGainLast);
  
  MPedestalPix &pedpix  = (*fPedestals)[fBlindPixelIdx];    
  
  if (&pedpix)
    {
      fBlindPixel->SetPed      ( pedpix.GetPedestal()   * fNumLoGainSamples );
      fBlindPixel->SetPedErr   ( pedpix.GetPedestalRms()* fNumLoGainSamples 
                                 / TMath::Sqrt((Float_t)fPedestals->GetTotalEntries()) );
      fBlindPixel->SetPedRms   ( pedpix.GetPedestalRms()* TMath::Sqrt((Float_t)fNumLoGainSamples) );
      fBlindPixel->SetPedRmsErr( fBlindPixel->GetPedErr()/2. );
    }
  
  return kTRUE;
}

// --------------------------------------------------------------------------
//
// FindSignalHiGain:
//
// - Loop from ptr to (ptr+fHiGainLast-fHiGainFirst)
// - Sum up contents of *ptr
// - If *ptr is greater than fSaturationLimit, raise sat by 1
// - If fHiLoLast is set, loop from logain to (logain+fHiLoLast)
// - Add contents of *logain to sum
// 
void MExtractBlindPixel::FindSignalHiGain(Byte_t *ptr, Byte_t *logain, Int_t &sum, Byte_t &sat) const
{

  Byte_t *end = ptr + fHiGainLast - fHiGainFirst + 1;
  
  while (ptr<end)
    {
      sum += *ptr;
      
      if (*ptr++ >= fSaturationLimit)
        sat++;
    }

  if (fHiLoLast == 0)
    return;
  
  end = logain + fHiLoLast;
  while (logain<end)
    {
      sum += *logain;
      
      if (*logain++ >= fSaturationLimit)
            sat++;
    }

}

// --------------------------------------------------------------------------
//
// FindSignalFilter:
//
// - Loop from ptr to (ptr+fLoGainLast-fLoGainFirst)
// - Sum up contents of *ptr
// - If *ptr is greater than fSaturationLimit, raise sat by 1
// 
void MExtractBlindPixel::FindSignalFilter(Byte_t *ptr, Int_t &sum, Byte_t &sat) const
{

  Byte_t *end = ptr + fLoGainLast - fLoGainFirst + 1;
  
  while (ptr<end)
    {
      sum += *ptr;
      
      if (*ptr++ >= fSaturationLimit)
        sat++;
    }
}

// --------------------------------------------------------------------------
//
// 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);
  
  Int_t sum   = 0;
  Byte_t sat  = 0;

  FindSignalFilter(pixel.GetHiGainSamples()+fLoGainFirst, sum, sat);

  if (sum > fNSBFilterLimit)
    {
      sum = -1;
      fBlindPixel->SetExtractedSignal(sum);
      fBlindPixel->SetNumSaturated(sat);
      fBlindPixel->SetReadyToSave();
      return kTRUE;
    }

  sum = 0;
  sat = 0;

  FindSignalHiGain(pixel.GetHiGainSamples()+fLoGainFirst, pixel.GetLoGainSamples(), sum, sat);
  
  fBlindPixel->SetExtractedSignal(sum);
  fBlindPixel->SetNumSaturated(sat);
  fBlindPixel->SetReadyToSave();
  
  return kTRUE;
}