/* ======================================================================== *\
!
! *
! * 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;
}