/* ======================================================================== *\
!
! *
! * 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 derives from MExtractPedestal. 
//  It calculates the pedestals 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">
*/
//End_Html
//
// 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 
//   fCheckWinLast    =  fgCheckWinLast    =  29
//   fExtractWinFirst =  fgExtractWinFirst =  15 
//   fExtractWinSize  =  fgExtractWinSize  =  6
//   fMaxSignalVar    =  fgMaxSignalVar    = 40;
//
//  Call: 
//  SetCheckRange(fCheckWinFirst,fCheckWinLast); 
//  to set the Window in which a signal is searched
//
//  SetExtractWindow(fExtractWinFirst,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 "MExtractPedestal.h"

#include "MExtractTimeAndCharge.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 UShort_t MPedCalcFromLoGain::fgCheckWinFirst   =  0;
const UShort_t MPedCalcFromLoGain::fgCheckWinLast    = 29;
const UShort_t MPedCalcFromLoGain::fgExtractWinFirst = 15;
const UShort_t MPedCalcFromLoGain::fgExtractWinSize  =  6;
const UShort_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)
{

  fName  = name  ? name  : "MPedCalcFromLoGain";
  fTitle = title ? title : "Task to calculate pedestals from lo-gains";
  
  SetCheckRange(fgCheckWinFirst, fgCheckWinLast);
  SetExtractWindow(fgExtractWinFirst, fgExtractWinSize);
  
  SetMaxSignalVar(fgMaxSignalVar);
  
  Clear();
}

void MPedCalcFromLoGain::ResetArrays()
{

  MExtractPedestal::ResetArrays();
  
  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
//
Bool_t MPedCalcFromLoGain::SetCheckRange(UShort_t chfirst, UShort_t chlast)
{

  Bool_t rc = kTRUE;
  
  if (chlast<=chfirst)
    {
      *fLog << warn << GetDescriptor();
      *fLog << " - WARNING: Last slice in SetCheckRange smaller than first slice... set to first+2" << endl;
      chlast = chfirst+1;
      rc = kFALSE;
    }

  fCheckWinFirst = chfirst;
  fCheckWinLast  = chlast;

  return rc;
}

// ---------------------------------------------------------------------------------
//
// 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
//
Bool_t MPedCalcFromLoGain::ReInit(MParList *pList)
{

  const UShort_t hisamples = fRunHeader->GetNumSamplesHiGain();
  const UShort_t losamples = fRunHeader->GetNumSamplesLoGain();

  UShort_t lastavailable   = hisamples+losamples-1;

  if (fExtractor)
    lastavailable = losamples-1;

  // If the size is not yet set, set the size
  if (fSumx.GetSize()==0)
  {
    const Int_t npixels  = fPedestalsOut->GetSize();
    fNumEventsUsed.Set(npixels);
    fTotalCounter.Set(npixels); 
  }
  
  if (fExtractWinLast > lastavailable) //changed to override check
   {
     const UShort_t diff = fExtractWinLast - lastavailable;
     *fLog << warn << GetDescriptor();
     *fLog << " - WARNING: Selected Extract Window ranges out of range...adjusting to last available slice ";
     *fLog << lastavailable << endl;
     
     fExtractWinLast -= diff;
     fExtractWinSize -= diff;
   }
  
   lastavailable = hisamples + losamples -1;
   
   if (fCheckWinLast > lastavailable) //changed to override check
     {
       *fLog << warn << GetDescriptor();
       *fLog << " - WARNING: Last Check Window slice out of range...adjusting to last available slice ";
       *fLog << lastavailable << endl;
       
       fCheckWinLast = lastavailable;
     }

  MExtractPedestal::ReInit(pList);
  
  return kTRUE;
}

// ---------------------------------------------------------------------------------
//
// Returns the pointer to slice "slice".
//
UShort_t MPedCalcFromLoGain::GetSlice(MRawEvtPixelIter *pixel, UInt_t slice)
{
    const UShort_t nh = (Int_t)fRunHeader->GetNumSamplesHiGain();

    Byte_t *ptr;

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

   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();
      const UInt_t aidx   = (*fGeom)[idx].GetAidx();
      const UInt_t sector = (*fGeom)[idx].GetSector();      
      
      UShort_t max = 0;
      UShort_t min = (UShort_t)-1;
      
      // Find the maximum and minimum signal per slice in the high gain window
      for (Int_t slice=fCheckWinFirst; slice<=fCheckWinLast; slice++)
        {
          const UShort_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;
      
      Float_t sum  = 0.;
      UInt_t  sumi = 0;
      
      //extract pedestal
      if (fExtractor)
        CalcExtractor( &pixel, sum, (*fPedestalsIn)[idx] );
      else
        {
          for(Int_t slice=fExtractWinFirst; slice<=fExtractWinLast; slice++)
            sumi += GetSlice(&pixel,slice);
          sum = (Float_t)sumi;
        }
      
      const Float_t sqrsum = sum*sum;
      
      fSumx[idx]           += sum;
      fSumx2[idx]          += sqrsum;
      fAreaSumx[aidx]      += sum;
      fAreaSumx2[aidx]     += sqrsum;
      fSectorSumx[sector]  += sum;      
      fSectorSumx2[sector] += sqrsum;      
      
      fNumEventsUsed[idx]   ++;
      fAreaFilled   [aidx]  ++;
      fSectorFilled [sector]++;
      
      if (!fExtractor)
        {
          //
          // Calculate the amplitude of the 150MHz "AB" noise
          //
          const UShort_t abFlag = (fExtractWinFirst + pixel.HasABFlag()) & 0x1;
          
          for (UShort_t islice=fExtractWinFirst; islice<=fExtractWinLast; islice+=2)
            {
              const UShort_t sliceAB0 = islice + abFlag;
              const UShort_t sliceAB1 = islice - abFlag + 1;
              const UShort_t ab0 = GetSlice(&pixel, sliceAB0);
              const UShort_t ab1 = GetSlice(&pixel, sliceAB1);
              fSumAB0[idx]        += ab0;
              fSumAB1[idx]        += ab1;
              fAreaSumAB0[aidx]   += ab0;
              fAreaSumAB1[aidx]   += ab1;      
              fSectorSumAB0[aidx] += ab0;
              fSectorSumAB1[aidx] += ab1;      
            }
        }
      
      if (!fPedestalUpdate || (UInt_t)fNumEventsUsed[idx]<fNumEventsDump)
        continue;
      
      CalcPixResults(fNumEventsDump, idx);
      fTotalCounter[idx]++;
      
      fNumEventsUsed[idx]=0;
      fSumx[idx]=0;
      fSumx2[idx]=0;
      fSumAB0[idx]=0;
      fSumAB1[idx]=0;
    }
  
  if (!(GetNumExecutions() % fNumAreasDump))
    {
      //
      // 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(fSectorFilled[aidx],napix,aidx);
        }
    }
  
  if (!(GetNumExecutions() % fNumSectorsDump))
    {
      //
      // Loop over the (two) 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(fSectorFilled[sector],nspix,sector);
        }
    }
  
  if (fPedestalUpdate)
    fPedestalsOut->SetReadyToSave();
  
  return kTRUE;
}

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


// --------------------------------------------------------------------------
//
// 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)
        {
          CalcPixResults(n, idx);
          fTotalCounter[idx]++;
        }
    }

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

      CalcAreaResults(fAreaFilled[aidx],napix,aidx);
    }
  
  //
  // Loop over the (six) sector indices to get the averaged pedestal per sector
  //
  for (UInt_t sector=0; sector<fSectorFilled.GetSize(); sector++)
    {
      const Int_t   nspix = fSectorValid.At(sector);
      if (nspix == 0)
        continue;
      
      CalcSectorResults(fSectorFilled[sector],nspix,sector);
    }

  fPedestalsOut->SetReadyToSave();

  return kTRUE;
}


// --------------------------------------------------------------------------
//
//  The following resources are available:
//    FirstCheckWindowSlice:  0
//    LastCheckWindowSlice:  29
//    MaxSignalVar:          40
//
Int_t MPedCalcFromLoGain::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{

  Bool_t rc=kFALSE;
  
  // Find resources for CheckWindow
  Int_t fs = fCheckWinFirst;
  Int_t ls = fCheckWinLast;
  if (IsEnvDefined(env, prefix, "CheckWinFirst", print))
    {
      fs = GetEnvValue(env, prefix, "CheckWinFirst", fs);
      rc = kTRUE;
    }
  if (IsEnvDefined(env, prefix, "CheckWinLast", print))
    {
      ls = GetEnvValue(env, prefix, "CheckWinLast", ls);
      rc = kTRUE;
    }
  
  SetCheckRange(fs,ls);
  
  // find resource for maximum signal variation
  if (IsEnvDefined(env, prefix, "MaxSignalVar", print))
    {
      SetMaxSignalVar(GetEnvValue(env, prefix, "MaxSignalVar", fMaxSignalVar));
      rc = kTRUE;
    }
  
  return MExtractPedestal::ReadEnv(env,prefix,print) ? kTRUE : rc;
}

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

  MExtractPedestal::Print(o);

  *fLog << "CheckWindow   from slice  " << fCheckWinFirst   << " to " << fCheckWinLast << " incl." << endl;
  *fLog << "Max. allowed signal variation:             " << fMaxSignalVar << endl;
  *fLog << endl;
}
