source: trunk/MagicSoft/Mars/msignal/MExtractor.cc@ 8239

/* ======================================================================== *\
!
! *
! * 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, 04/2004 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */
//////////////////////////////////////////////////////////////////////////////
//
// MExtractor
// ==========
//
//   Base class for the signal extractors, used the functions 
//   FindSignalHiGain() and FindSignalLoGain() to extract the signal and 
//   substract the pedestal value    
//
//   The following variables have to be set by the derived class and 
//   do not have defaults:
//   - fNumHiGainSamples
//   - fNumLoGainSamples
//   - fSqrtHiGainSamples
//   - fSqrtLoGainSamples
//
//   The signal extractor classes can be setup from an environmental
//   setup file. For more information see ReadEnv and MEvtLoop::ReadEnv
//
//
// IMPORTANT: For all classes you derive from MExtractor make sure that:
//    - Print() is correctly implemented
//    - Clone() works
//    - Class Version number != 0 and the I/O works PERFECTLY
//    - only data members which are necessary for the setup (not the ones
//      created in PreProcess and Process) are written
//    - the version number is maintained!
//    - If the flag fNoiseCalculation is set, the signal extractor should 
//      calculate the pure noise contriubtion from a fixed window in time.
//
// The following figure gives and example of possible inheritance trees. 
// An extractor class can inherit from each of the following base classes:
//    - MExtractor
//    - MExtractTime
//    - MExtractTimeAndCharge
//
//Begin_Html
/*
<img src="images/ExtractorClasses.gif">
*/
//End_Html
//
// Class Version 6:
//  +Float_t fResolutionPerPheHiGain; // Extractor-dependent charge resolution per phe for high-gain (see TDAS-0502).
//  +Float_t fResolutionPerPheLoGain; // Extractor-dependent charge resolution per phe for low-gain  (see TDAS-0502).
//
//
// Input Containers:
//   MRawEvtData
//   MRawRunHeader
//   MPedestalCam
//
// Output Containers:
//   MExtractedSignalCam
//
//////////////////////////////////////////////////////////////////////////////
#include "MExtractor.h"

#include <fstream>

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

#include "MParList.h"

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

#include "MPedestalCam.h"
#include "MPedestalPix.h"
//#include "MPedestalSubtractEvt.h"

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

ClassImp(MExtractor);

using namespace std;

const Byte_t  MExtractor::fgSaturationLimit = 245;
const TString MExtractor::fgNamePedestalCam = "MPedestalCam";
const TString MExtractor::fgNameSignalCam   = "MExtractedSignalCam";
const Float_t MExtractor::fgOffsetLoGain    = 1.51;   // 5 ns

// --------------------------------------------------------------------------
//
// Default constructor. 
//
// Set: 
// - all pointers to NULL
// - all variables to 0
// - fSaturationLimit to fgSaturationLimit
// - fNamePedestalCam to fgNamePedestalCam
// - fNameSignalCam to fgNameSignalCam
// - fNoiseCalculation to kFALSE 
//
// Call:
// -  AddToBranchList("MRawEvtData.*")
//
MExtractor::MExtractor(const char *name, const char *title)
    : fResolutionPerPheHiGain(0), fResolutionPerPheLoGain(0),
    fPedestals(NULL), fSignals(NULL), fRawEvt(NULL), fRunHeader(NULL),
    fSignal(NULL),
      /*fHiLoLast(0),*/ fNumHiGainSamples(0), fNumLoGainSamples(0)
{
    fName  = name  ? name  : "MExtractor";
    fTitle = title ? title : "Base class for signal extractors";

    AddToBranchList("MRawEvtData.*");

    SetNamePedestalCam();
    SetNameSignalCam();
    SetOffsetLoGain();
    SetSaturationLimit();
    SetNoiseCalculation(kFALSE);
}

void MExtractor::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast)
{
 
  Clear();

  fHiGainFirst = hifirst;
  fHiGainLast  = hilast;
  
  fLoGainFirst = lofirst;
  fLoGainLast  = lolast;
}

//-----------------------------------------------------------------------
//
// - Set the variable fHiLoLast to 0 (will be initialized later in ReInit()
// - Get the pointers to: 
//     MRawEvtData
//     MRawRunHeader
//     MPedestalCam
//
Int_t MExtractor::PreProcessStd(MParList *pList)
{

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

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

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


    if (fPedestals)
        return kTRUE;

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

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// The PreProcess searches for the following input containers:
//  - MRawEvtData
//  - MRawRunHeader
//  - MPedestalCam
//
// The following output containers are also searched and created if
// they were not found:
//
//  - MExtractedSignalCam
//
Int_t MExtractor::PreProcess(MParList *pList)
{
    fSignals = (MExtractedSignalCam*)pList->FindCreateObj("MExtractedSignalCam",AddSerialNumber(fNameSignalCam));
    if (!fSignals)
        return kFALSE;

    return PreProcessStd(pList);
}

// --------------------------------------------------------------------------
//
// The ReInit searches for:
// -  MRawRunHeader::GetNumSamplesHiGain()
// -  MRawRunHeader::GetNumSamplesLoGain()
//
// In case that the variable fLoGainLast is smaller than 
// the even part of the number of samples obtained from the run header, a
// warning is given an the range is set back accordingly. A call to:  
// - SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast-diff) 
// is performed in that case. The variable diff means here the difference 
// between the requested range (fLoGainLast) and the available one. Note that 
// the functions SetRange() are mostly overloaded and perform more checks, 
// modifying the ranges again, if necessary.
//
// In case that the variable fHiGainLast is smaller than the available range 
// obtained from the run header, a warning is given that a part of the low-gain 
// samples are used for the extraction of the high-gain signal. 
//
Bool_t MExtractor::ReInit(MParList *pList)
{
    const Int_t numl = fRunHeader->GetNumSamplesLoGain();
    const Int_t numh = fRunHeader->GetNumSamplesHiGain();
    const Int_t num  = numh+numl;

    if (fHiGainLast>=num)
    {
        *fLog << err << "ERROR - Last hi-gain slice must not exceed " << num-1 << endl;
        return kFALSE;
    }
    if (fLoGainLast>=num)
    {
        *fLog << err << "ERROR - Last lo-gain slice must not exceed " << num-1 << endl;
        return kFALSE;
    }

/*
  const Int_t logainsamples = fRunHeader->GetNumSamplesLoGain();

  Int_t lastdesired;
  Int_t lastavailable;

  if (logainsamples)
    {

      lastdesired   = (Int_t)(fLoGainLast);
      lastavailable = logainsamples-1;
      
      if (lastavailable < 0)
        *fLog << warn << GetDescriptor() << " - WARNING: Number of available Low-Gain Slices is smaller than or equal zero!" << endl;
      
    if (lastdesired > lastavailable)
      {
        const Int_t diff = lastdesired - lastavailable;
        
        *fLog << endl;
        *fLog << warn << GetDescriptor() << ": Selected Lo Gain FADC Window [";
        *fLog << Form("%2i,%2i", (int)fLoGainFirst, lastdesired);
        *fLog << "] ranges out of the available limits: [0," << Form("%2i", lastavailable) << "]" << endl;
        *fLog << GetDescriptor() << ": Will reduce the upper edge to " << (int)(fLoGainLast - diff) << endl;
        SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast-diff);
      }
    }
  else
    SetRange(fHiGainFirst, fHiGainLast, 0,0);      
  
  const Int_t higainsamples = fRunHeader->GetNumSamplesHiGain();

  if (higainsamples <= 0)
    {
      *fLog << err << GetDescriptor();
      *fLog << " - ERROR: Number of available High-Gain Slices is smaller than or equal zero!" << endl;
      return kFALSE;
    }

  lastdesired   = (Int_t)fHiGainLast;
  lastavailable = higainsamples-1;
  
  if (lastdesired > lastavailable)
    {
      const Int_t diff = lastdesired - lastavailable;
      
      *fLog << endl;
      *fLog << inf << GetDescriptor() << ": Selected Hi Gain FADC Window [";
      *fLog << Form("%2i,%2i", (int)fHiGainFirst,lastdesired);
      *fLog << "] ranges out of the available limits: [0," << Form("%2i", lastavailable) << "]" << endl;
      *fLog << inf << GetDescriptor() << ": Will use ";
      *fLog << Form("%2i", diff) << " samples from the Low-Gain for the High-Gain extraction";
      *fLog << endl;
      
      fHiGainLast -= diff;
      fHiLoLast    = diff;
    }
  
  */
    return kTRUE;
}

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

  MRawEvtPixelIter pixel(fRawEvt);

  while (pixel.Next())
    {
      Float_t sumhi = 0.;
      Byte_t sathi = 0;

      FindSignalHiGain(pixel.GetHiGainSamples()+fHiGainFirst, pixel.GetLoGainSamples(), sumhi, sathi);

      Float_t sumlo  = 0.;
      Byte_t  satlo  = 0;

      if (pixel.HasLoGain())
        FindSignalLoGain(pixel.GetLoGainSamples()+fLoGainFirst, sumlo, satlo);

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

      pix.SetExtractedSignal(sumhi - pedes*fNumHiGainSamples, pedrms*fSqrtHiGainSamples,
                             sumlo - pedes*fNumLoGainSamples, pedrms*fSqrtLoGainSamples);
      
      pix.SetGainSaturation(sathi, satlo);
      
    } /* while (pixel.Next()) */

    fSignals->SetReadyToSave();

    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 MExtractor::StreamPrimitive(ostream &out) const
{
  out << "   " << ClassName() << " " << GetUniqueName() << "(\"";
  out << "\"" << fName << "\", \"" << fTitle << "\");" << endl;
  
  if (fSaturationLimit!=fgSaturationLimit)
  {
      out << "   " << GetUniqueName() << ".SetSaturationLimit(";
      out << (int)fSaturationLimit << ");" << endl;
  }
  
  out << "   " << GetUniqueName() << ".SetRange(";
  out << (int)fHiGainFirst;
  out << ", " << (int)fHiGainLast;
  out << ", " << (int)fLoGainFirst;
  out << ", " << (int)fLoGainLast;
  out << ");" << endl;
}

// --------------------------------------------------------------------------
//
// Read the setup from a TEnv, eg:
//   MJPedestal.MExtractor.HiGainFirst: 5
//   MJPedestal.MExtractor.LoGainFirst: 5
//   MJPedestal.MExtractor.HiGainLast:  10
//   MJPedestal.MExtractor.LoGainLast:  10
//   MJPedestal.MExtractor.SaturationLimit: 88
//
Int_t MExtractor::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{
    Byte_t hf = fHiGainFirst;
    Byte_t lf = fLoGainFirst;
    Byte_t hl = fHiGainLast;
    Byte_t ll = fLoGainLast;

    Bool_t rc = kFALSE;

    if (IsEnvDefined(env, prefix, "HiGainFirst", print))
    {
        hf = GetEnvValue(env, prefix, "HiGainFirst", hf);
        rc = kTRUE;
    }
    if (IsEnvDefined(env, prefix, "LoGainFirst", print))
    {
        lf = GetEnvValue(env, prefix, "LoGainFirst", lf);
        rc = kTRUE;
    }

    if (IsEnvDefined(env, prefix, "HiGainLast", print))
    {
        hl = GetEnvValue(env, prefix, "HiGainLast", hl);
        rc = kTRUE;
    }
    if (IsEnvDefined(env, prefix, "LoGainLast", print))
    {
        ll = GetEnvValue(env, prefix, "LoGainLast", ll);
        rc = kTRUE;
    }

    SetRange(hf, hl, lf, ll);

    if (IsEnvDefined(env, prefix, "OffsetLoGain", print))
    {
        SetOffsetLoGain(GetEnvValue(env, prefix, "OffsetLoGain", fOffsetLoGain));
        rc = kTRUE;
    }

    if (IsEnvDefined(env, prefix, "SaturationLimit", print))
    {
        SetSaturationLimit(GetEnvValue(env, prefix, "SaturationLimit", fSaturationLimit));
        rc = kTRUE;
    }

    if (IsEnvDefined(env, prefix, "NoiseCalculation", print))
    {
        SetNoiseCalculation(GetEnvValue(env, prefix, "NoiseCalculation", fNoiseCalculation));
        rc = kTRUE;
    }

    // Be carefull: Returning kERROR is not forseen in derived classes
    return rc;
}

void MExtractor::Print(Option_t *o) const
{
    if (IsA()==MExtractor::Class())
        *fLog << GetDescriptor() << ":" << endl;

    *fLog << " Hi Gain Range:      " << Form("%2d %2d", fHiGainFirst, fHiGainLast) << endl;
    *fLog << " Lo Gain Range:      " << Form("%2d %2d", fLoGainFirst, fLoGainLast) << endl;
//    *fLog << " Gain Overlap to Lo: " << Form("%2d", fHiLoLast)        << endl;
    *fLog << " Saturation Lim:     " << Form("%3d", fSaturationLimit) << endl;
    *fLog << " Num Samples Hi/Lo:  " << Form("%2.1f %2.1f", fNumHiGainSamples, fNumLoGainSamples) << endl;
    if (fPedestals)
        *fLog << " Pedestals:          " << fPedestals->GetName() << ", " << fPedestals << endl;
}