source: trunk/MagicSoft/Mars/mpedestal/MPedCalcFromLoGain.cc@ 5357

/* ======================================================================== *\
!
! *
! * 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>
!   Author(s): Florian Goebel 06/2004 <mailto:fgoebel@mppmu.mpg.de>
!   Author(s): Nepomuk Otte 10/2004 <mailto:otte@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
// 
// 
//
//   MPedCalcLoGain
//
//
//  This task is devide form MPedCalcPedRun, described below. However, It 
//  calculates the pedstals using the low gain slices, whenever the difference 
//  between the highest and the lowest slice in the high gain
//  slices is below a given threshold (SetMaxHiGainVar). In this case the receiver
//  boards do not switch to lo gain and the so called lo gain slices are actually
//  high gain slices. 
//
//  MPedCalcLoGain also fills the ABoffset in MPedestalPix which allows to correct 
//  the 150 MHz clock noise.
//
//  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 MPedCalcFromLoGain)
// 2) The fitted mean pedestal per slice     (from MHPedestalCam)
// 3) The calculated pedestal RMS per slice  (from MPedCalcFromLoGain)
// 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: 
//  ====================
//  
// 
//   fCheckWinFirst =  fgCheckWinFirst =  0 
//   fHiGainLast  =  fgCheckWinLast  =  29
//   fExtractWinFirst =  fgExtractWinFirst =  0 
//   fExtractWinSize  =  fgExtractWinSize  =  6
//   fMaxSignalVar   =   fgMaxSignalVar  = 40;
//
//  Call: 
//  SetCheckRange(fCheckWinFirst,fCheckWinLast); 
//  to set the Window in which a signal is searched
//
//  SetExtractWindow(fExtractWindFirst,fExtractWinSize);
//  to set the Window from which a signal is extracted
//
//  SetMaxSignalVar(fMaxSignalVar);
//  set the maximum allowed difference between maximum and minimal signal in CheckWindow  
//
//   Variables:
//   fgCheckWinFirst;      First FADC slice to check for signal (currently set to: 0)
//   fgCheckWinLast:       Last FADC slice to check for signal (currently set to: 29)
//   fgExtractWinFirst:    First FADC slice to be used for pedestal extraction (currently set to: 15)
//   fgExtractWinSize:     Window size in slices used for pedestal extraction (currently set to: 6)
//   fgMaxSignalVar:       The maximum difference between the highest and lowest slice
//                         in the check window allowed in order to use event
//
//  Input Containers:
//   MRawEvtData
//   MRawRunHeader
//   MGeomCam
//
//  Output Containers:
//   MPedestalCam
//
//  See also: MPedestalCam, MPedestalPix, MHPedestalCam, MExtractor
//
/////////////////////////////////////////////////////////////////////////////
#include "MPedCalcFromLoGain.h"
#include "MExtractor.h"

#include "MParList.h"

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

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

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

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

ClassImp(MPedCalcFromLoGain);

using namespace std;

const Int_t MPedCalcFromLoGain::fgCheckWinFirst      = 0;
const Int_t MPedCalcFromLoGain::fgCheckWinLast       = 29;
const Int_t MPedCalcFromLoGain::fgExtractWinFirst    = 15;
const Int_t MPedCalcFromLoGain::fgExtractWinSize     = 6;
const Int_t MPedCalcFromLoGain::fgMaxSignalVar       = 40;

// --------------------------------------------------------------------------
//
// Default constructor: 
//
// Sets:
// - all pointers to NULL
//
// Calls: 
// - AddToBranchList("fHiGainPixId");
// - AddToBranchList("fHiGainFadcSamples");
// - SetCheckRange(fgCheckWinFirst, fgCheckWinLast, fgExtractWinFirst, fgExtractWinSize)
// - Clear()
//
MPedCalcFromLoGain::MPedCalcFromLoGain(const char *name, const char *title)
    : fGeom(NULL), fPedContainerName("MPedestalCam")
{

  fName  = name  ? name  : "MPedCalcFromLoGain";
  fTitle = title ? title : "Task to calculate pedestals from pedestal runs raw data";
  
  AddToBranchList("fHiGainPixId");
  AddToBranchList("fHiGainFadcSamples");
  
  SetCheckRange(fgCheckWinFirst, fgCheckWinLast);
  SetExtractWindow(fgExtractWinFirst, fgExtractWinSize);

  SetMaxSignalVar(fgMaxSignalVar);
  SetPedestalUpdate(kTRUE);
  Clear();

}

// --------------------------------------------------------------------------
//
// Sets:
// - fRawEvt to NULL
// - fRunHeader to NULL
// - fPedestals to NULL
//
// Resets:
// - fSumx
// - fSumx2
// - fSumAB0
// - fSumAB1
// - fNumEventsUsed
// - fTotalCounter
//
void MPedCalcFromLoGain::Clear(const Option_t *o)
{

  fRawEvt    = NULL;
  fRunHeader = NULL;
  fPedestals = NULL;

  // If the size is yet set, set the size
  if (fSumx.GetSize()>0)
  {
      // Reset contents of arrays.
      fSumx.Reset();
      fSumx2.Reset();
      fSumAB0.Reset();
      fSumAB1.Reset();
      fNumEventsUsed.Reset();
      fTotalCounter.Reset();
  }
}

// --------------------------------------------------------------------------
//
// SetCheckRange: 
// 
// Exits, if the first argument is smaller than 0
// Exits, if the the last argument is smaller than the first
//
void MPedCalcFromLoGain::SetCheckRange(Int_t chfirst, Int_t chlast)
{

  if(chfirst<0)
    {
      *fLog << warn << GetDescriptor() 
            << ": First slice in window to check for Signal <0, adjust:" << endl;
      exit(-1);
    }

  if(chlast<=chfirst)
    {
      *fLog << warn << GetDescriptor() 
            << ": Last slice in Check window smaller than first slice in window, adjust:" << endl;
      exit(-1);
    }
  
  fCheckWinFirst = chfirst;
  fCheckWinLast = chlast;

}

// --------------------------------------------------------------------------
//
// Exits:
// - if a window is odd
//
void MPedCalcFromLoGain::SetExtractWindow(Int_t windowf, Int_t windows)
{
 
  if(windowf<0)
    {
      *fLog << warn << GetDescriptor() 
            << Form(": First slice in Extract window has to be >0, adjust:")<< endl;
      exit(-1);
    }
 
  Int_t odd  = windows & 0x1;
 
  if (odd||(windows==0))
    {
      *fLog << warn << GetDescriptor() << ": Extract window size has to be even and larger 0, adjust!"<< endl;
      exit(-1);
    }
  
  fExtractWinSize  = windows;
  fExtractWinFirst = windowf;
  fExtractWinLast  = fExtractWinFirst+fExtractWinSize-1;
  
/*
    const Byte_t availhirange = (fHiGainLast-fHiGainFirst+1) & ~1;
    const Byte_t availlorange = (fLoGainLast-fLoGainFirst+1) & ~1;
  
    if (fWindowSizeHiGain > availhirange)
    {
        *fLog << warn;
        *fLog << GetDescriptor() << ": HiGain window " << (int)fWindowSizeHiGain;
        *fLog << " out of range [" << (int)fHiGainFirst;
        *fLog << "," << (int)fHiGainLast << "]" << endl;
        *fLog << "Will set window size to " << (int)availhirange << endl;
        fWindowSizeHiGain = availhirange;
    }
    
    if (fWindowSizeLoGain > availlorange)
    {
        *fLog << warn;
        *fLog << GetDescriptor() << ": LoGain window " << (int)fWindowSizeLoGain;
        *fLog << " out of range [" << (int)fLoGainFirst;
        *fLog << "," << (int)fLoGainLast << "]" << endl;
        *fLog << "Will set window size to " << (int)availlorange << endl;
        fWindowSizeLoGain = availlorange;
    }
*/
}

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

    Clear();

    fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
    if (!fRawEvt)
    {
        *fLog << err << "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;
    }

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


    fPedestals = (MPedestalCam*)pList->FindCreateObj("MPedestalCam", AddSerialNumber(fPedContainerName));
    if (!fPedestals)
        return kFALSE;

    
    *fLog << inf << "Parameters used for pedestal calculation:"<<endl;
    *fLog << inf <<"CheckWindow: First slice: "<<fCheckWinFirst<<"; Last slice: "<<fCheckWinLast<<endl;
    *fLog << inf <<"ExtractWindow: First slice: "<<fExtractWinFirst<<"; Last slice: "<<fExtractWinLast<<endl;
    *fLog << inf <<"Max allowed signal variation: "<<fMaxSignalVar<<endl;
    return kTRUE;
}

// ---------------------------------------------------------------------------------
//
// Checks:
// - if the number of available slices 
//   (fRunHeader->GetNumSamplesHiGain()+(Int_t)fRunHeader->GetNumSamplesLoGain()-1)
//   is smaller than the number of used slices
//   (fExtractWinSize+ fExtractWinFirst-1) or 
//    fCheckWinLast
//
// Sets the size (from MPedestalCam::GetSize() ) and resets the following arrays:
//  - fSumx
//  - fSumx2
//  - fSumAB0
//  - fSumAB1
//  - fNumEventsUsed
//  - fTotalCounter
//
Bool_t MPedCalcFromLoGain::ReInit(MParList *pList)
{

  Int_t lastavailableslice = (Int_t)fRunHeader->GetNumSamplesHiGain()
                            +(Int_t)fRunHeader->GetNumSamplesLoGain()-1;
  Int_t lastextractslice   =  fExtractWinSize+ fExtractWinFirst - 1;

   if ( lastextractslice  > lastavailableslice ) //changed to override check
    {
      *fLog << endl; 
      *fLog << warn << GetDescriptor()
            << Form(": Selected Extract Window ranges out of the available limits adjust. Last available slice is %2i",
                    lastavailableslice) << endl;
      return kFALSE;
    }

  if ( fCheckWinLast  > lastavailableslice ) //changed to override check
    {
      *fLog << endl; 
      *fLog << warn << GetDescriptor()
            << Form(": Last Check Window slice is out of the available limits adjust. Last available slice is %2i",
                    lastavailableslice) << endl;
      return kFALSE;
    }


 
  // If the size is not yet set, set the size
  if (fSumx.GetSize()==0)
  {
      const Int_t npixels = fPedestals->GetSize();

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

      // Reset contents of arrays.
      fSumx.Reset();
      fSumx2.Reset();
      fSumAB0.Reset();
      fSumAB1.Reset();
      fNumEventsUsed.Reset();
      fTotalCounter.Reset();
  }
  
  return kTRUE;
}

// ---------------------------------------------------------------------------------
//
// 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;
//
// Sets fTotalCounter up by 1.
//
void MPedCalcFromLoGain::Calc(ULong_t n, UInt_t idx)
{

  const ULong_t nsamplestot = n*fExtractWinSize;
  
  const Float_t sum  = fSumx.At(idx);
  const Float_t sum2 = fSumx2.At(idx);
  const Float_t ped  = sum/nsamplestot;
  
  // 1. Calculate the Variance of the sums:
  Float_t var = (sum2-sum*sum/n)/(n-1.);

  // 2. Scale the variance to the number of slices:
  var /= (Float_t)(fExtractWinSize);
  
  // 3. Calculate the RMS from the Variance:
  const Float_t rms = var<0 ? 0 : TMath::Sqrt(var);
  
  // 4. Calculate the amplitude of the 150MHz "AB" noise
  const Float_t abOffs = (fSumAB0[idx] - fSumAB1[idx]) / nsamplestot;
  
  (*fPedestals)[idx].Set(ped, rms, abOffs, n);
  
  fTotalCounter[idx]++;
}

// ---------------------------------------------------------------------------------
//
// Returns the pointer to slice "slice".
//
UInt_t MPedCalcFromLoGain::GetSlice(MRawEvtPixelIter *pixel, UInt_t slice)
{

 const UInt_t nh = (Int_t)fRunHeader->GetNumSamplesHiGain();

 Byte_t *ptr;

 if(slice>=nh)
   {
     ptr = pixel->GetLoGainSamples();
     ptr += slice - nh;
   }
 else
   {
     ptr = pixel->GetHiGainSamples();
     ptr += slice;
   }

   return *ptr;
}

// --------------------------------------------------------------------------
//
// 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 MPedCalcFromLoGain::Process()
{
    MRawEvtPixelIter pixel(fRawEvt);

    while (pixel.Next())
    {

         const UInt_t idx = pixel.GetPixelId();

        Int_t max = 0;
        Int_t min = 1025;
    
        // Find the maximum and minimum signal per slice in the high gain window
     
        for(Int_t slice=fCheckWinFirst; slice<=fCheckWinLast; slice++)
          {
    
          Int_t svalue = GetSlice(&pixel,slice);
            if (svalue > max) {
                max = svalue;
            }
            if (svalue < min) {
                min = svalue;
            }
        }

        // If the maximum in the high gain window is smaller than
        if (max-min>=fMaxSignalVar || max>=250)
            continue;

        UInt_t sum = 0;
        UInt_t sqr = 0;

        //extract pedestal
        for(Int_t slice=fExtractWinFirst; slice<=fExtractWinLast; slice++)
          {    
          UInt_t svalue = GetSlice(&pixel,slice);
          sum += svalue;
          sqr += svalue*svalue;
        }

        const Float_t msum   = (Float_t)sum;
        const Float_t sqrsum = msum*msum;

        fSumx[idx]  += msum;
        fSumx2[idx] += sqrsum;
        fNumEventsUsed[idx]++;

        // Calculate the amplitude of the 150MHz "AB" noise

        Int_t abFlag = (fExtractWinFirst + pixel.HasABFlag()) & 0x1;

        //cout << " MPedCalcFromLoGain: idx: " << idx << " abFlag: " << abFlag << endl;

        for (Int_t islice=fExtractWinFirst; islice<=fExtractWinLast; islice+=2)
        {
            Int_t sliceAB0 = islice + abFlag;
            Int_t sliceAB1 = islice - abFlag + 1;
            fSumAB0[idx] += GetSlice(&pixel, sliceAB0);
            fSumAB1[idx] += GetSlice(&pixel, sliceAB1);
        }

        if (!fPedestalUpdate || fNumEventsUsed[idx]<fNumEventsDump)
          continue;

        Calc(fNumEventsDump, idx);

        fNumEventsUsed[idx]=0;
        fSumx[idx]=0;
        fSumx2[idx]=0;
        fSumAB0[idx]=0;
        fSumAB1[idx]=0;
    }

    if (fPedestalUpdate)
      fPedestals->SetReadyToSave();

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Compute signal mean and rms in the whole run and store it in MPedestalCam
//
Int_t MPedCalcFromLoGain::PostProcess()
{
    // Compute pedestals and rms from the whole run
    if (fPedestalUpdate)
      return kTRUE;

    *fLog << flush << inf << "Calculating Pedestals..." << flush;

    const Int_t npix = fGeom->GetNumPixels();
    for (Int_t idx=0; idx<npix; idx++)
      {
        const ULong_t n = fNumEventsUsed[idx];
        if (n>1)
          Calc(n, idx);
      }

    fPedestals->SetReadyToSave();
    return kTRUE;
}


Int_t MPedCalcFromLoGain::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{
    if (MExtractor::ReadEnv(env, prefix, print)==kERROR)
        return kERROR;

    Bool_t rc=kFALSE;

    Int_t lw = fExtractWinSize;
    Int_t wf = fExtractWinFirst;

    if (IsEnvDefined(env, prefix, "ExtractWindowSize", print))
    {
        lw = GetEnvValue(env, prefix, "ExtractWindowSize", lw);
        wf = GetEnvValue(env, prefix, "ExtractWindowFirst", wf);
        rc = kTRUE;

    }
    SetExtractWindow(wf,lw);

    Int_t num = fNumEventsDump;
    if (IsEnvDefined(env, prefix, "NumEventsDump", print))
    {
        num = GetEnvValue(env, prefix, "NumEventsDump", num);
        rc = kTRUE;
    }
    SetNumEventsDump(num);

    Byte_t max = fMaxSignalVar;
    if (IsEnvDefined(env, prefix, "MaxSignalVar", print))
    {
        max = GetEnvValue(env, prefix, "MaxSignalVar", max);
        rc = kTRUE;
    }
    SetMaxSignalVar(max);

    Bool_t upd = fPedestalUpdate;
    if (IsEnvDefined(env, prefix, "PedestalUpdate", print))
    {
        upd = GetEnvValue(env, prefix, "PedestalUpdate", upd);
        rc = kTRUE;
    }
    SetPedestalUpdate(upd);

    return rc;
}