source: trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.cc@ 5502

/* ======================================================================== *\
!
! *
! * 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): Josep Flix 04/2001 <mailto:jflix@ifae.es>
!   Author(s): Thomas Bretz 05/2001 <mailto:tbretz@astro.uni-wuerzburg.de>
!   Author(s): Sebastian Commichau 12/2003 
!   Author(s): Javier Rico 01/2004 <mailto:jrico@ifae.es>
!   Author(s): Markus Gaug 01/2004 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */
/////////////////////////////////////////////////////////////////////////////
//
//   MPedCalcPedRun
//
//  This task takes a pedestal run file and fills MPedestalCam during
//  the Process() with the pedestal and rms computed in an event basis.
//  In the PostProcess() MPedestalCam is finally filled with the pedestal
//  mean and rms computed in a run basis.
//  More than one run (file) can be merged
//
//  MPedCalcPedRun 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 )
// 
//  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">
*/
//
// 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.
// 
// The plots have been produced on run 20123. You can reproduce them using
// the macro pedestalstudies.C
// 
//  Usage of this class: 
//  ====================
// 
//  Call: SetRange(higainfirst, higainlast, logainfirst, logainlast) 
//  to modify the ranges in which the window is allowed to move. 
//  Defaults are: 
// 
//   fHiGainFirst =  fgHiGainFirst =  0 
//   fHiGainLast  =  fgHiGainLast  =  29
//   fLoGainFirst =  fgLoGainFirst =  0 
//   fLoGainLast  =  fgLoGainLast  =  14
//
//  Call: SetWindowSize(windowhigain, windowlogain) 
//  to modify the sliding window widths. Windows have to be an even number. 
//  In case of odd numbers, the window will be modified.
//
//  Defaults are:
//
//   fHiGainWindowSize = fgHiGainWindowSize = 14
//   fLoGainWindowSize = fgLoGainWindowSize = 0
//
//
// ToDo:
//   - Take a setup file (ReadEnv-implementation) as input
//
//
//  Input Containers:
//   MRawEvtData
//   MRawRunHeader
//   MRawEvtHeader
//   MGeomCam
//
//  Output Containers:
//   MPedestalCam
//
//  See also: MPedestalCam, MPedestalPix, MHPedestalCam, MExtractor
//
/////////////////////////////////////////////////////////////////////////////
#include "MPedCalcPedRun.h"
#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"

#include "MStatusDisplay.h"

ClassImp(MPedCalcPedRun);

using namespace std;

const UShort_t MPedCalcPedRun::fgExtractWinFirst       = 0;
const UShort_t MPedCalcPedRun::fgExtractWinSize        = 6;
const UInt_t   MPedCalcPedRun::gkFirstRunWithFinalBits = 45589;
// --------------------------------------------------------------------------
//
// Default constructor: 
//
// Sets:
// - all pointers to NULL
// - fWindowSizeHiGain to fgHiGainWindowSize
// - fWindowSizeLoGain to fgLoGainWindowSize
//
// Calls: 
// - AddToBranchList("fHiGainPixId");
// - AddToBranchList("fHiGainFadcSamples");
// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast)
// - Clear()
//
MPedCalcPedRun::MPedCalcPedRun(const char *name, const char *title)
{
  fName  = name  ? name  : "MPedCalcPedRun";
  fTitle = title ? title : "Task to calculate pedestals from pedestal runs raw data";

  AddToBranchList("fHiGainPixId");
  AddToBranchList("fHiGainFadcSamples");
  
  SetExtractWindow(fgExtractWinFirst, fgExtractWinSize);

  Clear();
  
  fFirstRun   = kTRUE;
  fSkip       = kFALSE;
  fOverlap    = 0;
  fUsedEvents = 0;
}



void MPedCalcPedRun::ResetArrays()
{

  MExtractPedestal::ResetArrays();
  fUsedEvents    = 0;
}

// --------------------------------------------------------------------------
//
// In case that the variables fExtractLast is greater than the number of 
// High-Gain FADC samples obtained from the run header, the flag 
// fOverlap is set to the difference and fExtractLast is set back by the
// same number of slices.
//
// The run type is checked for "Pedestal" (run type: 1) 
// and the variable fSkip is set in that case
//
Bool_t MPedCalcPedRun::ReInit(MParList *pList)
{

  MExtractPedestal::ReInit(pList);

  const UShort_t lastavailable = fRunHeader->GetNumSamplesHiGain()-1;
  
  if (fExtractWinLast > lastavailable)
    {
      const UShort_t diff = fExtractWinLast - lastavailable;
      *fLog << endl;
      *fLog << warn << GetDescriptor()
            << Form("%s%2i%s%2i%s%2i%s",": Selected ExtractWindow [",
                    fExtractWinFirst,",",fExtractWinLast,
                    "] ranges out of available limits: [0,",lastavailable,"].") << endl;
      *fLog << warn << GetDescriptor() 
            << Form("%s%2i%s",": Will use ",diff," samples from the 'Low-Gain' slices for the pedestal extraction")
            << endl;
      fExtractWinLast -= diff;
      fOverlap         = diff;
    }

  const UShort_t runtype = fRunHeader->GetRunType();

  switch (runtype)
    {
    case 1: 
      fFirstRun = kFALSE;
      fSkip     = kFALSE;
      return kTRUE;
      break;
    default:
      fSkip     = kTRUE;
      if (!fFirstRun)
        {
          *fLog << endl;
          *fLog << inf << GetDescriptor() << " : Finalize pedestal calculations..." << flush;
          CallPostProcess();
          Reset();
        }
      return kTRUE;        
      break;
    }
  
  Print();

  return kTRUE;
}

// --------------------------------------------------------------------------
//
// Return kTRUE (without doing anything) in case that the run type is not 
// equal to 1 (pedestal run)
//
// Fill the MPedestalCam container with the signal mean and rms for the event.
// Store the measured signal in arrays fSumx and fSumx2 so that we can 
// calculate the overall mean and rms in the PostProcess()
//
Int_t MPedCalcPedRun::Process()
{

  if (fSkip && !IsPedBitSet())
    return kTRUE;

  /*
  *fLog << err << dec << fEvtHeader->GetTriggerID() << endl;
  for (Int_t i=16; i>= 0; i--)
    *fLog << inf << (fEvtHeader->GetTriggerID() >> i & 1);
  *fLog  << endl;
  *fLog << warn << hex << fEvtHeader->GetTriggerID() << endl;
  */

  fUsedEvents++;

  MRawEvtPixelIter pixel(fRawEvt);
  
  while (pixel.Next())
    {
      const UInt_t idx    = pixel.GetPixelId();
      const UInt_t aidx   = (*fGeom)[idx].GetAidx();
      const UInt_t sector = (*fGeom)[idx].GetSector();      

      Float_t sum = 0.;
      UInt_t  ab0 = 0;
      UInt_t  ab1 = 0;
      
      if (fExtractor)
        CalcExtractor( &pixel, sum, (*fPedestalsIn)[idx]);
      else
        CalcSums( &pixel, sum, ab0, ab1);

      const Float_t msum = (Float_t)sum;
      
      fSumx[idx]          += msum;
      fAreaSumx[aidx]     += msum;
      fSectorSumx[sector] += msum;      

      const Float_t sqrsum  = msum*msum;
      fSumx2[idx]          += sqrsum;
      fAreaSumx2[aidx]     += sqrsum;
      fSectorSumx2[sector] += sqrsum;      

      fSumAB0[idx]        += ab0;
      fSumAB1[idx]        += ab1;

      fAreaSumAB0[aidx]   += ab0;
      fAreaSumAB1[aidx]   += ab1;      
      
      fSectorSumAB0[aidx] += ab0;
      fSectorSumAB1[aidx] += ab1;      
    }
  
  fPedestalsOut->SetReadyToSave();

  return kTRUE;
}



void MPedCalcPedRun::CalcExtractor( MRawEvtPixelIter *pixel, Float_t &sum, MPedestalPix &ped)
{
  
  Byte_t  sat = 0;
  Byte_t *first  = pixel->GetHiGainSamples() + fExtractWinFirst;
  Byte_t *logain = pixel->GetLoGainSamples();
  const Bool_t abflag  = pixel->HasABFlag();
  Float_t dummy;
  fExtractor->FindTimeAndChargeHiGain(first,logain,sum,dummy,dummy,dummy,sat,ped,abflag);
}


void MPedCalcPedRun::CalcSums( MRawEvtPixelIter *pixel, Float_t &sum, UInt_t &ab0, UInt_t &ab1)
{

  Byte_t *higain = pixel->GetHiGainSamples() + fExtractWinFirst;
  Byte_t *ptr = higain;
  Byte_t *end = ptr + fExtractWinSize;

  const Bool_t abflag = pixel->HasABFlag();

  Int_t sumi = 0;
      
  Int_t  cnt = 0;
  do
    {
      sumi += *ptr;
      if (pixel->IsABFlagValid())
        {
          const Int_t abFlag = (fExtractWinFirst + abflag + cnt) & 0x1;
          if (abFlag)
            ab1 += *ptr;
          else
            ab0 += *ptr;
          
          cnt++;
        }
    }
  while (++ptr != end);
  
  if (fOverlap != 0)
    {
      ptr = pixel->GetLoGainSamples();
      end = ptr + fOverlap;
      
      do
        {
          sumi += *ptr;
          if (pixel->IsABFlagValid())
            {
              const Int_t abFlag = (fExtractWinFirst + abflag + cnt) & 0x1;
              if (abFlag)
                ab1 += *ptr;
              else
                ab0 += *ptr;
              
              cnt++;
            }
        }
      while (++ptr != end);
    }

  sum = (Float_t)sumi;

}

// --------------------------------------------------------------------------
//
// Compute signal mean and rms in the whole run and store it in MPedestalCam
//
Int_t MPedCalcPedRun::PostProcess()
{

  if (fUsedEvents == 0)
    return kTRUE;

  MRawEvtPixelIter pixel(fRawEvt);
  
  while (pixel.Next())
    {
      const Int_t  pixid  = pixel.GetPixelId();
      CalcPixResults(fUsedEvents,pixid);
    }

  //
  // Loop over the (two) area indices to get the averaged pedestal per aidx
  //
  for (UInt_t aidx=0; aidx<fAreaValid.GetSize(); aidx++)
    {
      const Int_t   napix = fAreaValid.At(aidx);
      if (napix == 0)
        continue;

      CalcAreaResults(fUsedEvents,napix,aidx);
    }
  
  //
  // Loop over the (six) sector indices to get the averaged pedestal per sector
  //
  for (UInt_t sector=0; sector<fSectorValid.GetSize(); sector++)
    {
      const Int_t   nspix = fSectorValid.At(sector);
      if (nspix == 0)
        continue;
      
      CalcSectorResults(fUsedEvents,nspix,sector);
    }
  
  fPedestalsOut->SetTotalEntries(fUsedEvents*fExtractWinSize);
  fPedestalsOut->SetReadyToSave();

  return kTRUE;
}

// -----------------------------------------------------------------------
//
// Analogue to the MTask::CallPostProcess, but without the manipulation 
// of the bit fIsPreProcessed. Needed in order to call PostProcess multiple 
// times in the intensity calibration.
//
Int_t MPedCalcPedRun::CallPostProcess()
{

  *fLog << all << fName << "... " << flush;
  if (fDisplay)
    fDisplay->SetStatusLine2(*this);
  
  return PostProcess();
}

//-------------------------------------------------------------
// 
// Return if the pedestal bit was set from the calibration trigger box.
// The last but one bit is used for the "pedestal-bit".
//
// This bit is set since run gkFinalizationTriggerBit
//
Bool_t MPedCalcPedRun::IsPedBitSet()
{
  if (fRunHeader->GetRunNumber() < gkFirstRunWithFinalBits)
    return kFALSE;
      
  return (fEvtHeader->GetTriggerID() >> 3 & 1);
}

void MPedCalcPedRun::Print(Option_t *o) const
{

  MExtractPedestal::Print(o);

  *fLog << "Number overlap low-gain slices:           " << fOverlap << endl;
  *fLog << "First run out of sequence:                " << (fFirstRun?"yes":"no") << endl;
  *fLog << "Skip this run:                            " << (fSkip?"yes":"no") << endl;
  *fLog << "Number of used events so far:             " << fUsedEvents << endl;
  *fLog << endl;
}