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

/* ======================================================================== *\
!
! *
! * 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 "MExtractTimeAndCharge.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"

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)
    : 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();
  SetNumEventsDump();
  SetNumAreasDump();
  SetNumSectorsDump();
  
  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;
  fPedestalsIn  = NULL;
  fPedestalsOut = 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... " 
        " raising from " << windows << " to " << flush;
      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;
    }
  
  fPedestalsIn = (MPedestalCam*)pList->FindCreateObj("MPedestalCam", AddSerialNumber(fNamePedestalCamIn.Data()));
  if (!fPedestalsIn)
    {
      *fLog << err << fNamePedestalCamIn.Data() << " could not be found nor created... aborting" << endl;
      return kFALSE;
    }
  
  fPedestalsOut = (MPedestalCam*)pList->FindCreateObj("MPedestalCam", AddSerialNumber(fNamePedestalCamOut.Data()));
  if (!fPedestalsOut)
    {
      *fLog << err << fNamePedestalCamOut.Data() << " could not be found nor created... aborting" << endl;
      return kFALSE;
    }
  
  *fLog << inf;
  
  return kTRUE;
}

// ---------------------------------------------------------------------------------
//
// 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->GetAverageAreas();
    const Int_t sectors  = fPedestalsOut->GetAverageSectors();
    
    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)fExtractor->GetNumHiGainSamples());
    }

  Print();

  return kTRUE;
}

// --------------------------------------------------------------------------
//
//  The following resources are available:
//    ExtractWindowFirst:    15
//    ExtractWindowSize:      6
//    NumEventsDump:        500
//    PedestalUpdate:       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 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.At(pixid);
  const Float_t sum2        = fSumx2.At(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.At(aidx);
  const Float_t sum2        = fAreaSumx2.At(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.At(sector);
  const Float_t sum2        = fSectorSumx2.At(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 incoming MPedestalCam Container:  " << fNamePedestalCamIn.Data() << endl;
    *fLog << "Name of outgoing MPedestalCam Container:  " << fNamePedestalCamOut.Data() << endl;
    *fLog << "Number events for pedestal calculation:   " << fNumEventsDump << endl;
    *fLog << "Number events for av. areas calculation:  " << fNumAreasDump << endl;
    *fLog << "Number events for av. sector calculation: " << fNumSectorsDump << endl;
    *fLog << "Pedestal Update is                        " << (fPedestalUpdate?"on":"off") << endl;
    *fLog << "ExtractWindow from slice " << fExtractWinFirst << " to " << fExtractWinLast << " incl." << endl;
}