source: trunk/MagicSoft/Mars/mpedestal/MExtractPedestal.cc@ 6192

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

/////////////////////////////////////////////////////////////////////////////
// 
//   MExtractPedestal
//
//  Pedestal Extractor base class
//
//  Input Containers:
//   MRawEvtData
//   MRawRunHeader
//   MRawEvtHeader
//   MGeomCam
//   MPedestalCam
//
//  Output Containers:
//   MPedestalCam
//
//  This class should be used for pedestal extractors with the following facilities:
//  a) Standardized calculation of AB-noise, mean pedestals and RMS
//  b) Standardized treatment of area- and sector averaged pedestals values
//  c) Possibility to use a signal extractor to be applied on the pedestals
//  d) Possibility to handle two MPedestalCams: one for the signal extractor and 
//     a second to be filled during the pedestal calculating process.  
//
//  ad a): Every calculated value is referred to one FADC slice (e.g. Mean pedestal per slice),
//         RMS per slice. 
//         MExtractPedestal applies the following formula (1):
// 
//         Pedestal per slice = sum(x_i) / n / slices
//         PedRMS per slice   = Sqrt(  ( sum(x_i^2) - sum(x_i)^2/n ) / n-1 / slices )
//         AB-Offset per slice = (sumAB0 - sumAB1) / n / slices 
//
//         where x_i is the sum of "slices" FADC slices and sum means the sum over all
//         events. "n" is the number of events, "slices" is the number of summed FADC samples.
// 
//         Note that the slice-to-slice fluctuations are not Gaussian, but Poissonian, thus 
//         asymmetric and they are correlated.
// 
//         It is important to know that the Pedestal per slice and PedRMS per slice depend 
//         on the number of used FADC slices, as seen in the following plots:
//
//Begin_Html
/*
<img src="images/PedestalStudyInner.gif">
*/
//End_Html
//
//Begin_Html
/*
<img src="images/PedestalStudyOuter.gif">
*/
//End_Html
//
//        The plots show the inner and outer pixels, respectivly and have the following meaning:
// 
//        1) The calculated mean pedestal per slice (from MPedCalcPedRun)
//        2) The fitted mean pedestal per slice     (from MHPedestalCam)
//        3) The calculated pedestal RMS per slice  (from MPedCalcPedRun)
//        4) The fitted sigma of the pedestal distribution per slice
//                                                  (from MHPedestalCam)
//        5) The relative difference between calculation and histogram fit
//           for the mean
//        6) The relative difference between calculation and histogram fit
//           for the sigma or RMS, respectively.
// 
//        The calculated means do not change significantly except for the case of 2 slices, 
//        however the RMS changes from 5.7 per slice in the case of 2 extracted slices 
//        to 8.3 per slice in the case of 26 extracted slices. This change is very significant.
// 
// ad b)  Every calculated value is referred to one FADC slice and one (averaged) pixel,
//        (e.g. Mean Pedestal per area index per slice per pixel, etc. )
//
//         MExtractPedestal applies the following formula (2):
// 
//         Averaged Pedestal per slice = sum(x_i) / n / slices / n_pix
//         PedRMS per slice   = Sqrt(  ( sum(x_i^2) - sum(x_i)^2/n ) / n-1 / slices / n_pix )
//         AB-Offset per slice = (sumAB0 - sumAB1) / n / slices / n_pix
//
//         where x_i is the sum of "slices" FADC slices and sum means the sum over all
//         events and all concerned pixels. 
//         "n" is the number of events, "slices" is the number of summed FADC samples and
//         "n_pix" is the number of pixels belonging to the specific area index or camera sector. 
// 
//         Calculating these averaged event-by-event values is very important to trace coherent
//         fluctuations. An example is given in the following plots:
//
//Begin_Html
/*
<img src="images/PedestalOscillations.gif">
*/
//End_Html
//
//        The plots show the extracted pedestals of the inner pixels (obtained
//        with MHPedestalCam), averaged on an event-by-event basis from 
//        run 13428 with switched off camera LV. 
//        The meaning of the four plots is:
// 
//        1) The distribution of the averaged pedestals
//        2) The averaged pedestals vs. time. 
//           One can see clearly the oscillation pattern
//        3) The fourier transform of the averaged pedestals vs. time. 
//           One can see clearly a peak at a certain frequency
//        4) The projection of the fourier components with the non-exponential
//           (and therefore significant) outlier.
//
//    ad c) Many signal extractors, especially those using a sliding window
//          have biases and their resolutions for zero-signals do not agree 
//          with the pedestal RMS. For the F-Factor method in the calibration
//          and the image cleaning, however, both have to be known and measured. 
//          
//          For this reason, a signal extractor can be handed over to the 
//          pedestal extractor and applied on the pedestal events with the 
//          function SetExtractor(). 
//          The results will get stored in an MPedestalCam. 
//        
//          Note that only extractors deriving from MExtractTimeAndCharge
//          can be used.
//
//   ad d)  The signal extractors themselves need a pedestal to be subtracted 
//          from the FADC slices. 
//          If the user wishes that the pededestals do not get overwritten by 
//          the results from the signal extractor, a different named MPedestalCam
//          can be created with the function: SetNamePedestalOut(). 
//
//  See also: MPedestalCam, MPedestalPix, MPedCalcPedRun, MPedCalcFromLoGain
//
/////////////////////////////////////////////////////////////////////////////
#include "MExtractPedestal.h"

#include "MParList.h"

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

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

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

#include "MGeomPix.h"
#include "MGeomCam.h"

#include "MTaskEnv.h"
#include "MExtractTimeAndCharge.h"

ClassImp(MExtractPedestal);

using namespace std;

const TString MExtractPedestal::fgNamePedestalCam = "MPedestalCam";

// --------------------------------------------------------------------------
//
// Default constructor: 
//
// Sets:
// - all pointers to NULL
//
// Calls: 
// - AddToBranchList("fHiGainPixId");
// - AddToBranchList("fHiGainFadcSamples");
// - Clear()
//
MExtractPedestal::MExtractPedestal(const char *name, const char *title)
    : fRandomCalculation(kTRUE), fGeom(NULL), fPedestalsIn(NULL), fPedestalsOut(NULL), 
      fExtractor(NULL), fExtractWinFirst(0), fExtractWinSize(0)
{
    fName  = name  ? name  : "MExtractPedestal";
    fTitle = title ? title : "Base class to calculate pedestals";

    AddToBranchList("fHiGainPixId");
    AddToBranchList("fLoGainPixId");
    AddToBranchList("fHiGainFadcSamples");
    AddToBranchList("fLoGainFadcSamples");

    SetPedestalUpdate(kTRUE);

    SetNamePedestalCamIn();
    SetNamePedestalCamOut();
    SetNumDump();

    Clear();
}

void MExtractPedestal::ResetArrays()
{
    // Reset contents of arrays.
    fSumx.Reset();
    fSumx2.Reset();
    fSumAB0.Reset();
    fSumAB1.Reset();
    fAreaSumx.Reset();
    fAreaSumx2.Reset();
    fAreaSumAB0.Reset();
    fAreaSumAB1.Reset();
    fAreaFilled.Reset();
    fAreaValid.Reset();
    fSectorSumx.Reset();
    fSectorSumx2.Reset();
    fSectorSumAB0.Reset();
    fSectorSumAB1.Reset();
    fSectorFilled.Reset();
    fSectorValid.Reset();

}

// --------------------------------------------------------------------------
//
// Resets Arrays:
//
// Sets:
// - fRawEvt to NULL
// - fRunHeader to NULL
// - fEvtHeader to NULL
// - fPedestalsIn to NULL
// - fPedestalsOut to NULL
//
void MExtractPedestal::Clear(const Option_t *o)
{
 
  fRawEvt       = NULL;
  fRunHeader    = NULL;
  fEvtHeader    = NULL;

  // If the size is yet set, set the size
  if (fSumx.GetSize()>0)
    ResetArrays();

}

// --------------------------------------------------------------------------
//
// Checks:
// - if a window is odd
// 
Bool_t MExtractPedestal::SetExtractWindow(UShort_t windowf, UShort_t windows)
{

  Bool_t rc = kTRUE;

  const Int_t odd  = windows & 0x1;
  
  if (odd)
  {
      *fLog << warn << GetDescriptor();
      *fLog << " - WARNING: Window size in SetExtraxtWindow has to be even... ";
      *fLog << " raising from " << windows << " to ";
      windows += 1;
      *fLog << windows << "!" << endl;
      rc = kFALSE;
  }

  if (windows==0)
  {
      *fLog << warn << GetDescriptor();
      *fLog << " - WARNING: Window size in SetExtraxtWindow has to be > 0... adjusting to 2!" << endl;
      windows = 2;
      rc = kFALSE;
  }

  fExtractWinSize  = windows;
  fExtractWinFirst = windowf;
  fExtractWinLast  = fExtractWinFirst+fExtractWinSize-1;
  
  return rc;
}

// --------------------------------------------------------------------------
//
// Look for the following input containers:
//
//  - MRawEvtData
//  - MRawRunHeader
//  - MRawEvtHeader
//  - MGeomCam
// 
// The following output containers are also searched and created if
// they were not found:
//
//  - MPedestalCam with the name fPedContainerName
//
Int_t MExtractPedestal::PreProcess(MParList *pList)
{

  Clear();
  
  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;
  }
  
  fEvtHeader = (MRawEvtHeader*)pList->FindObject(AddSerialNumber("MRawEvtHeader"));
  if (!fEvtHeader)
  {
      *fLog << err << AddSerialNumber("MRawEvtHeader") << " not found... aborting." << endl;
      return kFALSE;
  }
  
  fGeom = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
  if (!fGeom)
  {
      *fLog << err << AddSerialNumber("MGeomCam") << " not found... aborting." << endl;
      return kFALSE;
  }

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

  if (!fPedestalsOut)
  {
      fPedestalsOut = (MPedestalCam*)pList->FindCreateObj("MPedestalCam", AddSerialNumber(fNamePedestalCamOut));
      if (!fPedestalsOut)
          return kFALSE;
  }

  *fLog << inf;
  Print();

  // FIMXE: MUST call fExtractor->PreProcess()???
  return kTRUE;
}

Int_t MExtractPedestal::Process()
{
    if (fExtractor)
        fExtractor->SetNoiseCalculation(fRandomCalculation);

    const Int_t rc = Calc();

    if (fExtractor)
        fExtractor->SetNoiseCalculation(kFALSE);

    return rc;
}

// ---------------------------------------------------------------------------------
//
// Sets the size (from MPedestalCam::GetSize() ) and resets the following arrays:
//  - fSumx
//  - fSumx2
//  - fSumAB0
//  - fSumAB1
//  - fAreaSumx
//  - fAreaSumx2
//  - fAreaSumAB0
//  - fAreaSumAB1
//  - fAreaFilled
//  - fAreaValid 
//  - fSectorSumx 
//  - fSectorSumx2
//  - fSectorSumAB0
//  - fSectorSumAB1
//  - fSectorFilled
//  - fSectorValid
//
Bool_t MExtractPedestal::ReInit(MParList *pList)
{
  // If the size is not yet set, set the size
  if (fSumx.GetSize()==0)
  {
      const Int_t npixels  = fPedestalsOut->GetSize();
      const Int_t areas    = fPedestalsOut->GetNumAverageArea();
      const Int_t sectors  = fPedestalsOut->GetNumAverageSector();

      fSumx.  Set(npixels);
      fSumx2. Set(npixels);
      fSumAB0.Set(npixels);
      fSumAB1.Set(npixels);

      fAreaSumx.  Set(areas);
      fAreaSumx2. Set(areas);
      fAreaSumAB0.Set(areas);
      fAreaSumAB1.Set(areas);
      fAreaFilled.Set(areas);
      fAreaValid .Set(areas);

      fSectorSumx.  Set(sectors);
      fSectorSumx2. Set(sectors);
      fSectorSumAB0.Set(sectors);
      fSectorSumAB1.Set(sectors);
      fSectorFilled.Set(sectors);
      fSectorValid .Set(sectors);

      for (Int_t i=0; i<npixels; i++)
      {
          const UInt_t aidx   = (*fGeom)[i].GetAidx();
          const UInt_t sector = (*fGeom)[i].GetSector();

          fAreaValid  [aidx]  ++;
          fSectorValid[sector]++;
      }
  }

  if (fExtractor)
  {
      if (!fExtractor->InitArrays())
          return kFALSE;

      SetExtractWindow(fExtractor->GetHiGainFirst(), (Int_t)TMath::Nint(fExtractor->GetNumHiGainSamples()));
  }

  return kTRUE;
}

Int_t MExtractPedestal::PostProcess()
{
    fPedestalsIn = NULL;
    return kTRUE;
}


// --------------------------------------------------------------------------
//
//  The following resources are available:
//    ExtractWindowFirst:    15
//    ExtractWindowSize:      6
//    NumEventsDump:        500
//    PedestalUpdate:       yes
//    RandomCalculation:    yes
//
Int_t MExtractPedestal::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{
    Bool_t rc=kFALSE;

    // find resource for numeventsdump
    if (IsEnvDefined(env, prefix, "NumEventsDump", print))
    {
        SetNumEventsDump(GetEnvValue(env, prefix, "NumEventsDump", (Int_t)fNumEventsDump));
        rc = kTRUE;
    }

    // find resource for numeventsdump
    if (IsEnvDefined(env, prefix, "NumAreasDump", print))
    {
        SetNumAreasDump(GetEnvValue(env, prefix, "NumAreasDump", (Int_t)fNumAreasDump));
        rc = kTRUE;
    }

    // find resource for numeventsdump
    if (IsEnvDefined(env, prefix, "NumSectorsDump", print))
    {
        SetNumSectorsDump(GetEnvValue(env, prefix, "NumSectorsDump", (Int_t)fNumSectorsDump));
        rc = kTRUE;
    }

    // find resource for pedestal update
    if (IsEnvDefined(env, prefix, "PedestalUpdate", print))
    {
        SetPedestalUpdate(GetEnvValue(env, prefix, "PedestalUpdate", fPedestalUpdate));
        rc = kTRUE;
    }

    // find resource for random calculation
    if (IsEnvDefined(env, prefix, "RandomCalculation", print))
    {
        SetRandomCalculation(GetEnvValue(env, prefix, "RandomCalculation", fRandomCalculation));
        rc = kTRUE;
    }

    // Find resources for ExtractWindow
    Int_t ef = fExtractWinFirst;
    Int_t es = fExtractWinSize;
    if (IsEnvDefined(env, prefix, "ExtractWinFirst", print))
    {
        ef = GetEnvValue(env, prefix, "ExtractWinFirst", ef);
        rc = kTRUE;
    }
    if (IsEnvDefined(env, prefix, "ExtractWinSize", print))
    {
        es = GetEnvValue(env, prefix, "ExtractWinSize", es);
        rc = kTRUE;
    }

    SetExtractWindow(ef,es);

    // find resource for MPedestalCam
    if (IsEnvDefined(env, prefix, "NamePedestalCamIn", print))
    {
        SetNamePedestalCamIn(GetEnvValue(env, prefix, "NamePedestalCamIn", fNamePedestalCamIn));
        rc = kTRUE;
    }

    if (IsEnvDefined(env, prefix, "NamePedestalCamOut", print))
    {
        SetNamePedestalCamOut(GetEnvValue(env, prefix, "NamePedestalCamOut", fNamePedestalCamOut));
        rc = kTRUE;
    }

    return rc;
}

// ---------------------------------------------------------------------------------
//
// Calculates for pixel "idx":
//
// Ped per slice      = sum / n / fExtractWinSize;
// RMS per slice      = sqrt { (sum2 -  sum*sum/n) / (n-1) / fExtractWinSize }
// ABOffset per slice = (fSumAB0[idx] - fSumAB1[idx]) / n / fExtractWinSize;
//
// Stores the results in MPedestalCam[pixid]
//
void MExtractPedestal::CalcPixResults(const UInt_t nevts, const UInt_t pixid)
{
    const Float_t sum  = fSumx[pixid];
    const Float_t sum2 = fSumx2[pixid];

    // 1. Calculate the mean of the sums:
    Float_t ped        = sum/nevts;

    // 2. Calculate the Variance of the sums:
    Float_t var = (sum2-sum*sum/nevts)/(nevts-1.);

    // 3. Calculate the amplitude of the 150MHz "AB" noise
    Float_t abOffs = (fSumAB0[pixid] - fSumAB1[pixid]) / nevts;

    // 4. Scale the mean, variance and AB-noise to the number of slices:
    ped    /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;
    var    /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;
    abOffs /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;

    // 5. Calculate the RMS from the Variance:
    const Float_t rms = var<0 ? 0 : TMath::Sqrt(var);

    (*fPedestalsOut)[pixid].Set(ped, rms, abOffs, nevts);
}

// ---------------------------------------------------------------------------------
//
// Calculates for area idx "aidx" with "napix" valid pixels:
//
// Ped per slice      = sum / nevts / fExtractWinSize / napix;
// RMS per slice      = sqrt { (sum2 -  sum*sum/nevts) / (nevts-1) / fExtractWinSize / napix }
// ABOffset per slice = (fSumAB0[idx] - fSumAB1[idx]) / nevts / fExtractWinSize / napix;
//
// Stores the results in MPedestalCam::GetAverageArea(aidx)
//
void MExtractPedestal::CalcAreaResults(const UInt_t nevts, const UInt_t napix, const UInt_t aidx)
{
    const Float_t sum  = fAreaSumx[aidx];
    const Float_t sum2 = fAreaSumx2[aidx];

    // 1. Calculate the mean of the sums:
    Float_t ped        = sum/nevts;

    // 2. Calculate the Variance of the sums:
    Float_t var = (sum2-sum*sum/nevts)/(nevts-1.);

    // 3. Calculate the amplitude of the 150MHz "AB" noise
    Float_t abOffs = (fAreaSumAB0[aidx] - fAreaSumAB1[aidx]) / nevts;

    // 4. Scale the mean, variance and AB-noise to the number of slices:
    ped    /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;
    var    /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;
    abOffs /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;

    // 5. Scale the mean, variance and AB-noise to the number of pixels:
    ped    /= napix;
    var    /= napix;
    abOffs /= napix;

    // 6. Calculate the RMS from the Variance:
    const Float_t rms = var<0 ? 0 : TMath::Sqrt(var);

    fPedestalsOut->GetAverageArea(aidx).Set(ped, rms, abOffs, nevts);
}

// ---------------------------------------------------------------------------------
//
// Calculates for sector idx "sector" with "nspix" valid pixels:
//
// Ped per slice      = sum / nevts / fExtractWinSize / nspix;
// RMS per slice      = sqrt { (sum2 -  sum*sum/nevts) / (nevts-1) / fExtractWinSize / nspix }
// ABOffset per slice = (fSumAB0[idx] - fSumAB1[idx]) / nevts / fExtractWinSize / nspix;
//
// Stores the results in MPedestalCam::GetAverageSector(sector)
//
void MExtractPedestal::CalcSectorResults(const UInt_t nevts, const UInt_t nspix, const UInt_t sector)
{
    const Float_t sum  = fSectorSumx[sector];
    const Float_t sum2 = fSectorSumx2[sector];

    // 1. Calculate the mean of the sums:
    Float_t ped        = sum/nevts;

    // 2. Calculate the Variance of the sums:
    Float_t var = (sum2-sum*sum/nevts)/(nevts-1.);

    // 3. Calculate the amplitude of the 150MHz "AB" noise
    Float_t abOffs = (fSectorSumAB0[sector] - fSectorSumAB1[sector]) / nevts;

    // 4. Scale the mean, variance and AB-noise to the number of slices:
    ped    /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;
    var    /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;
    abOffs /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;

    // 5. Scale the mean, variance and AB-noise to the number of pixels:
    ped    /= nspix;
    var    /= nspix;
    abOffs /= nspix;

    // 6. Calculate the RMS from the Variance:
    const Float_t rms = var<0 ? 0 : TMath::Sqrt(var);

    fPedestalsOut->GetAverageSector(sector).Set(ped, rms, abOffs, nevts);
}

void MExtractPedestal::Print(Option_t *o) const
{
    *fLog << GetDescriptor() << ":" << endl;
    *fLog << "Name of input  MPedestalCam:  " << (fPedestalsIn?fPedestalsIn->GetName():fNamePedestalCamIn.Data()) << " (" << fPedestalsIn << ")" << endl;
    *fLog << "Name of output MPedestalCam:  " << (fPedestalsOut?fPedestalsOut->GetName():fNamePedestalCamOut.Data()) << " (" << fPedestalsOut << ")" << endl;
    *fLog << "Pedestal Update is            " << (fPedestalUpdate?"on":"off") << endl;
    if (fPedestalUpdate)
    {
        *fLog << "Num evts for pedestal   calc: " << fNumEventsDump << endl;
        *fLog << "Num evts for avg.areas  calc: " << fNumAreasDump << endl;
        *fLog << "Num evts for avg.sector calc: " << fNumSectorsDump << endl;
    }
    if (fExtractor)
    {
        *fLog << "Extractor used:               " << fExtractor->ClassName() << " (";
        *fLog << (fRandomCalculation?"":"non-") << "random)" << endl;
    }
}