source: trunk/MagicSoft/Mars/mpedestal/MPedCalcFromData.cc@ 5219

/* ======================================================================== *\
!
! *
! * 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 expressed
! * or implied warranty.
! *
!
!
!   Author(s): Josep Flix 06/2004 <mailto:jflix@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
//   MPedCalcFromData                                                      //
//                                                                         //
//  Input Containers:                                                      //
//   MRawEvtData                                                           //
//                                                                         //
//  Output Containers:                                                     //
//   MPedestalCam                                                          //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////

#include "MPedCalcFromData.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"

ClassImp(MPedCalcFromData);

using namespace std;

MPedCalcFromData::MPedCalcFromData(const char *name, const char *title) 
{
    fName  = name  ? name  : "MPedCalcFromData";
    fTitle = title ? title : "Task to calculate pedestals from data runs";

    AddToBranchList("fHiGainPixId");
    AddToBranchList("fHiGainFadcSamples");

    SetRange(fLoGainFirst, fLoGainLast);
    Clear();
}

void MPedCalcFromData::Clear(const Option_t *o)
{

  fRawEvt    = NULL;
  fRunHeader = NULL;

}

void MPedCalcFromData::SetLoRange(Byte_t lofirst, Byte_t lolast)
{
    fLoGainFirst = lofirst;
    fLoGainLast = lolast;
    
    fWindowSizeLoGain = lolast - lofirst;

}

Int_t MPedCalcFromData::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;
    }
    
    fPedestals = (MPedestalCam*)pList->FindCreateObj("MPedestalCam");
    if (!fPedestals)
        return kFALSE;
    
    return kTRUE;
}


Bool_t MPedCalcFromData::ReInit(MParList *pList)
{
  
  Int_t npixels  = fPedestals->GetSize();
  
  if (fSumx.GetSize()==0)
    {
      fSumx. Set(npixels);
      fSumx2.Set(npixels);
      fEvtCounter.Set(npixels);
      fTotalCounter.Set(npixels);

      fEvtCounter.Reset();
      fTotalCounter.Reset();
      fSumx.Reset();
      fSumx2.Reset();
    }
  
  return kTRUE;
      
}

Int_t MPedCalcFromData::Process()
{

  MRawEvtPixelIter pixel(fRawEvt);
  
  while (pixel.Next())
    {

      const UInt_t idx    = pixel.GetPixelId();
      
      if ( (UInt_t)pixel.GetMaxHiGainSample() < fHiGainThreshold ) {
          
          fEvtCounter[idx]++;

          Byte_t *ptr = pixel.GetLoGainSamples() + fLoGainFirst;
          Byte_t *end = ptr + fWindowSizeLoGain;
          
          UInt_t sum = 0;
          UInt_t sqr = 0;
          
          if (fWindowSizeLoGain != 0)
          {
              do
              {
                  sum += *ptr;
                  sqr += *ptr * *ptr;
              }
              while (++ptr != end);
          }
          
          const Float_t msum = (Float_t)sum;
          fSumx[idx]          += msum;
          
          const Float_t sqrsum  = msum*msum;
          fSumx2[idx]          += sqrsum;

          if ((UInt_t)fEvtCounter[idx] == fDump){
              
              // Compute pedestals and rms from the sample
              const ULong_t n     = fWindowSizeLoGain*fDump;

              const Float_t sum  = fSumx.At(idx);
              const Float_t sum2 = fSumx2.At(idx);
              const Float_t higainped = sum/n;        

              // 1. Calculate the Variance of the sums:
              Float_t higainVar = (sum2-sum*sum/fDump)/(fDump-1.);
              // 2. Scale the variance to the number of slices:
              higainVar /= (Float_t)(fWindowSizeLoGain);
              // 3. Calculate the RMS from the Variance:
              (*fPedestals)[idx].Set(higainped, higainVar < 0 ? 0. : TMath::Sqrt(higainVar));

              fTotalCounter[idx]++;
              fEvtCounter[idx]=0;
              fSumx[idx]=0;
              fSumx2[idx]=0;
          };
      }
      
    };
  
  fPedestals->SetReadyToSave();
  return kTRUE;
}