/* ======================================================================== *\
!
! *
! * 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): Markus Gaug   02/2004 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */
/////////////////////////////////////////////////////////////////////////////
//                                                               
// MHCalibrationChargeCam                                               
//
// Fills the extracted signals of MExtractedSignalCam into the MHGausEvents-classes 
// MHCalibrationChargeHiGainPix and MHCalibrationChargeLoGainPix for every:
//
// - Pixel, stored in the TObjArray's MHCalibrationCam::fHiGainArray and 
//   MHCalibrationCam::fLoGainArray
//
// - Average pixel per AREA index (e.g. inner and outer for the MAGIC camera), 
//   stored in the TObjArray's MHCalibrationCam::fAverageHiGainAreas and 
//   MHCalibrationCam::fAverageLoGainAreas
//
// - Average pixel per camera SECTOR (e.g. sectors 1-6 for the MAGIC camera), 
//   stored in the TObjArray's MHCalibrationCam::fAverageHiGainSectors and 
//   MHCalibrationCam::fAverageLoGainSectors
// 
// Every signal is taken from MExtractedSignalCam and filled into a histogram and 
// an array, in order to perform a Fourier analysis (see MHGausEvents). 
// The signals are moreover averaged on an event-by-event basis and written into 
// the corresponding average pixels.
//
// Additionally, the (FADC slice) position of the maximum is stored in an Absolute 
// Arrival Time histogram. This histogram serves for a rough cross-check if the 
// signal does not lie at or outside the edges of the extraction window. 
//
// The Charge histograms are fitted to a Gaussian, mean and sigma with its errors 
// and the fit probability are extracted. If none of these values are NaN's and 
// if the probability is bigger than MHGausEvents::fProbLimit (default: 0.5%), 
// the fit is declared valid.
// Otherwise, the fit is repeated within ranges of the previous mean 
// +- MHGausEvents::fPickupLimit (default: 5) sigma (see MHGausEvents::RepeatFit())
// In case this does not make the fit valid, the histogram means and RMS's are 
// taken directly (see MHGausEvents::BypassFit()) and the following flags are set:
// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainNotFitted ) or  
// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainNotFitted ) and 
// - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnreliableRun   ) 
// 
// Outliers of more than MHGausEvents::fPickupLimit (default: 5) sigmas 
// from the mean are counted as Pickup events (stored in MHGausEvents::fPickup) 
//
// Unless more than fNumHiGainSaturationLimit (default: 1%) of the overall FADC 
// slices show saturation, the following flag is set:
// - MCalibrationChargePix::SetHiGainSaturation();
// In that case, the calibration constants are derived from the low-gain results.
//
// If more than fNumLoGainSaturationLimit (default: 1%) of the overall 
// low-gain FADC slices saturate, the following flags are set:
// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainSaturation ) and
// - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnsuitableRun    )
// 
// The class also fills arrays with the signal vs. event number, creates a fourier 
// spectrum and investigates if the projected fourier components follow an exponential 
// distribution. In case that the probability of the exponential fit is less than 
// MHGausEvents::fProbLimit (default: 0.5%), the following flags are set:
// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainOscillating ) or
// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainOscillating ) and
// - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnreliableRun     )
// 
// This same procedure is performed for the average pixels.
//
// The following results are written into MCalibrationChargeCam:
//
// - MCalibrationPix::SetHiGainSaturation() 
// - MCalibrationPix::SetHiGainMean()
// - MCalibrationPix::SetHiGainMeanErr()
// - MCalibrationPix::SetHiGainSigma()
// - MCalibrationPix::SetHiGainSigmaErr()
// - MCalibrationPix::SetHiGainProb()
// - MCalibrationPix::SetHiGainNumPickup()
//
// - MCalibrationPix::SetLoGainMean()
// - MCalibrationPix::SetLoGainMeanErr()
// - MCalibrationPix::SetLoGainSigma()
// - MCalibrationPix::SetLoGainSigmaErr()
// - MCalibrationPix::SetLoGainProb()
// - MCalibrationPix::SetLoGainNumPickup()
//
// - MCalibrationChargePix::SetAbsTimeMean()
// - MCalibrationChargePix::SetAbsTimeRms()
//
// For all averaged areas, the fitted sigma is multiplied with the square root of 
// the number involved pixels in order to be able to compare it to the average of 
// sigmas in the camera.
//
/////////////////////////////////////////////////////////////////////////////
#include "MHCalibrationChargeCam.h"
#include "MHCalibrationCam.h"

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

#include "MParList.h"

#include "MHCalibrationChargeHiGainPix.h"
#include "MHCalibrationChargeLoGainPix.h"
#include "MHCalibrationChargePix.h"

#include "MCalibrationChargeCam.h"
#include "MCalibrationChargePix.h"

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

#include "MHGausEvents.h"

#include "MBadPixelsCam.h"
#include "MBadPixelsPix.h"

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

#include "MExtractedSignalCam.h"
#include "MExtractedSignalPix.h"

ClassImp(MHCalibrationChargeCam);

using namespace std;

const Float_t MHCalibrationChargeCam::fgNumHiGainSaturationLimit = 0.01;
const Float_t MHCalibrationChargeCam::fgNumLoGainSaturationLimit = 0.005;
const Float_t MHCalibrationChargeCam::fgTimeLowerLimit           = 1.;
const Float_t MHCalibrationChargeCam::fgTimeUpperLimit           = 2.;
// --------------------------------------------------------------------------
//
// Default Constructor. 
//
// Sets:
// - all pointers to NULL
//
// Initializes:
// - fNumHiGainSaturationLimit to fgNumHiGainSaturationLimit
// - fNumLoGainSaturationLimit to fgNumLoGainSaturationLimit
// - fTimeLowerLimit           to fgTimeLowerLimit 
// - fTimeUpperLimit           to fgTimeUpperLimit 
//
MHCalibrationChargeCam::MHCalibrationChargeCam(const char *name, const char *title)
    : fRawEvt(NULL)
{
  fName  = name  ? name  : "MHCalibrationChargeCam";
  fTitle = title ? title : "Class to fill the calibration histograms ";
  
  SetNumHiGainSaturationLimit(fgNumHiGainSaturationLimit);
  SetNumLoGainSaturationLimit(fgNumLoGainSaturationLimit);
  SetTimeLowerLimit();
  SetTimeUpperLimit();
}

// --------------------------------------------------------------------------
//
// Gets the pointers to:
// - MRawEvtData
//
Bool_t MHCalibrationChargeCam::SetupHists(const MParList *pList)
{

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

  return kTRUE;
}

// --------------------------------------------------------------------------
//
// Gets or creates the pointers to:
// - MExtractedSignalCam
// - MCalibrationChargeCam
// - MBadPixelsCam
//
// Initializes the number of used FADC slices from MExtractedSignalCam 
// into MCalibrationChargeCam and test for changes in that variable
//
// Initializes, if empty to MGeomCam::GetNumPixels():
// - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray
//
// Initializes, if empty to MGeomCam::GetNumAreas() for:
// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas
//
// Initializes, if empty to MGeomCam::GetNumSectors() for:
// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors
//
// Calls MHCalibrationCam::InitHists() for every entry in:
// - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray
// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas
// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors
//
// Sets Titles and Names for the Charge Histograms:
// - MHCalibrationCam::fAverageHiGainAreas
// - MHCalibrationCam::fAverageHiGainSectors
// 
// Sets number of bins to MHCalibrationCam::fAverageNbins for:
// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas
// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors
// 
Bool_t MHCalibrationChargeCam::ReInitHists(MParList *pList)
{

  MExtractedSignalCam *signal = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
  if (!signal)
  {
      *fLog << err << "MExtractedSignalCam not found... abort." << endl;
      return kFALSE;
  }

  fCam = (MCalibrationCam*)pList->FindObject("MCalibrationChargeCam");
  if (!fCam)
    {
      fCam = (MCalibrationCam*)pList->FindCreateObj(AddSerialNumber("MCalibrationChargeCam"));
      if (!fCam)
        {
          gLog << err << "Cannot find nor create MCalibrationChargeCam ... abort." << endl;
          return kFALSE;
        }
      else 
        fCam->Init(*fGeom);
    }

  fFirstHiGain = signal->GetFirstUsedSliceHiGain();
  fLastHiGain  = signal->GetLastUsedSliceHiGain();
  fFirstLoGain = signal->GetFirstUsedSliceLoGain();
  fLastLoGain  = signal->GetLastUsedSliceLoGain();

  const Float_t numhigain = signal->GetNumUsedHiGainFADCSlices();
  const Float_t numlogain = signal->GetNumUsedLoGainFADCSlices();  

  if (fCam->GetNumHiGainFADCSlices() == 0.)
    fCam->SetNumHiGainFADCSlices ( numhigain );
  else if (fCam->GetNumHiGainFADCSlices() != numhigain)
    {
      *fLog << err << GetDescriptor() 
            << ": Number of High Gain FADC extraction slices has changed, abort..." << endl;
      return kFALSE;
    }

  if (fCam->GetNumLoGainFADCSlices() == 0.)
    fCam->SetNumLoGainFADCSlices ( numlogain );
  else if (fCam->GetNumLoGainFADCSlices() != numlogain)
    {
      *fLog << err << GetDescriptor() 
            << ": Number of Low Gain FADC extraction slices has changes, abort..." << endl;
      return kFALSE;
    }

  const Int_t npixels  = fGeom->GetNumPixels();
  const Int_t nsectors = fGeom->GetNumSectors();
  const Int_t nareas   = fGeom->GetNumAreas();

  if (fHiGainArray->GetEntries()==0)
  {
      fHiGainArray->Expand(npixels);
      for (Int_t i=0; i<npixels; i++)
      {
	  (*fHiGainArray)[i] = new MHCalibrationChargeHiGainPix;
          InitHists((*this)[i],(*fBadPixels)[i],i);
      }
  }

  if (fLoGainArray->GetEntries()==0)
  {
      fLoGainArray->Expand(npixels);
      
      for (Int_t i=0; i<npixels; i++)
      {
	  (*fLoGainArray)[i] = new MHCalibrationChargeLoGainPix;
          InitHists((*this)(i),(*fBadPixels)[i],i);
      }
      
  }

  if (fAverageHiGainAreas->GetEntries()==0)
  {
    fAverageHiGainAreas->Expand(nareas);
    
    for (Int_t j=0; j<nareas; j++)
      {
        (*fAverageHiGainAreas)[j] = 
          new MHCalibrationChargeHiGainPix("AverageHiGainArea",
                                           "Average HiGain FADC sums area idx ");

        MHCalibrationChargePix &hist = (MHCalibrationChargePix&)GetAverageHiGainArea(j);

        hist.GetHGausHist()->SetTitle("Summed FADC slices average HiGain Area Idx ");
        hist.SetNbins(fAverageNbins);
        hist.GetHAbsTime()->SetTitle("Absolute Arrival Time average HiGain Area Idx ");

        InitHists(hist,fCam->GetAverageBadArea(j),j);

      }
  }


  if (fAverageLoGainAreas->GetEntries()==0)
    {
      fAverageLoGainAreas->Expand(nareas);
      
      for (Int_t j=0; j<nareas; j++)
        {
	  (*fAverageLoGainAreas)[j] = 
            new MHCalibrationChargeLoGainPix("AverageLoGainArea",
                                             "Average LoGain FADC sums of pixel area idx ");

        MHCalibrationChargePix &hist = (MHCalibrationChargePix&)GetAverageLoGainArea(j);

        hist.GetHGausHist()->SetTitle("Summed FADC slices average LoGain Area Idx ");
        hist.SetNbins(fAverageNbins);
        hist.GetHAbsTime()->SetTitle("Absolute Arrival Time average LoGain Area Idx ");

        InitHists(hist,fCam->GetAverageBadArea(j),j);

        }
    }
  
  if (fAverageHiGainSectors->GetEntries()==0)
  {
      fAverageHiGainSectors->Expand(nsectors);

      for (Int_t j=0; j<nsectors; j++)
      {
	  (*fAverageHiGainSectors)[j] = 
            new MHCalibrationChargeHiGainPix("AverageHiGainSector",
                                             "Average HiGain FADC sums of pixel sector ");

          MHCalibrationChargePix &hist = (MHCalibrationChargePix&)GetAverageHiGainSector(j);

          hist.GetHGausHist()->SetTitle("Summed FADC slices average HiGain Sector ");
          hist.SetNbins(fAverageNbins);
          hist.GetHAbsTime()->SetTitle("Absolute Arrival Time average HiGain Sector ");

          InitHists(hist,fCam->GetAverageBadSector(j),j);

      }
  }

  if (fAverageLoGainSectors->GetEntries()==0)
  {
      fAverageLoGainSectors->Expand(nsectors);
  
      for (Int_t j=0; j<nsectors; j++)
      {
	  (*fAverageLoGainSectors)[j] = 
            new MHCalibrationChargeLoGainPix("AverageLoGainSector",
                                             "Average LoGain FADC sums of pixel sector ");

          MHCalibrationChargePix &hist = (MHCalibrationChargePix&)GetAverageLoGainSector(j);

          hist.GetHGausHist()->SetTitle("Summed FADC slices average LoGain Sector ");
          hist.SetNbins(fAverageNbins);
          hist.GetHAbsTime()->SetTitle("Absolute Arrival Time average LoGain Sector ");

          InitHists(hist,fCam->GetAverageBadSector(j),j);
          
      }
  }

  return kTRUE;
}

  
// --------------------------------------------------------------------------
//
// Retrieves from MExtractedSignalCam:
// - first used LoGain FADC slice
//
// Retrieves from MGeomCam:
// - number of pixels
// - number of pixel areas
// - number of sectors
//
// For all TObjArray's (including the averaged ones), the following steps are performed: 
//
// 1) Fill Charges histograms (MHGausEvents::FillHistAndArray()) with:
// - MExtractedSignalPix::GetExtractedSignalHiGain();
// - MExtractedSignalPix::GetExtractedSignalLoGain();
//
// 2) Set number of saturated slices (MHCalibrationChargePix::SetSaturated()) with:
// - MExtractedSignalPix::GetNumHiGainSaturated();
// - MExtractedSignalPix::GetNumLoGainSaturated();
//
// 3) Fill AbsTime histograms (MHCalibrationChargePix::FillAbsTime()) with:
// - MRawEvtPixelIter::GetIdxMaxHiGainSample();       
// - MRawEvtPixelIter::GetIdxMaxLoGainSample(first slice);
//
Bool_t MHCalibrationChargeCam::FillHists(const MParContainer *par, const Stat_t w)
{

  MExtractedSignalCam *signal = (MExtractedSignalCam*)par;
  if (!signal)
    {
      *fLog << err << "No argument in MExtractedSignalCam::Fill... abort." << endl;
      return kFALSE;
    }
  
  const UInt_t npixels  = fGeom->GetNumPixels();
  const UInt_t nareas   = fGeom->GetNumAreas();
  const UInt_t nsectors = fGeom->GetNumSectors();
  const UInt_t lofirst  = signal->GetFirstUsedSliceLoGain();

  Float_t sumhiarea  [nareas],   sumloarea  [nareas],   timehiarea  [nareas],   timeloarea  [nareas];
  Float_t sumhisector[nsectors], sumlosector[nsectors], timehisector[nsectors], timelosector[nsectors];
  Int_t   sathiarea  [nareas],   satloarea  [nareas];
  Int_t   sathisector[nsectors], satlosector[nsectors];

  memset(sumhiarea,   0, nareas * sizeof(Float_t));
  memset(sumloarea,   0, nareas * sizeof(Float_t));
  memset(timehiarea,  0, nareas * sizeof(Float_t));
  memset(timeloarea,  0, nareas * sizeof(Float_t));
  memset(sathiarea,   0, nareas * sizeof(Int_t  ));
  memset(satloarea,   0, nareas * sizeof(Int_t  ));
  memset(sumhisector, 0, nsectors*sizeof(Float_t));
  memset(sumlosector, 0, nsectors*sizeof(Float_t));
  memset(timehisector,0, nsectors*sizeof(Float_t));
  memset(timelosector,0, nsectors*sizeof(Float_t));
  memset(sathisector, 0, nsectors*sizeof(Int_t  ));
  memset(satlosector, 0, nsectors*sizeof(Int_t  ));

  for (UInt_t i=0; i<npixels; i++)
    {

      MHCalibrationChargePix &histhi = (MHCalibrationChargePix&)(*this)[i];
      MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)(*this)(i);

      if (histhi.IsExcluded())
	continue;

      const MExtractedSignalPix &pix = (*signal)[i];
      
      const Float_t sumhi  = pix.GetExtractedSignalHiGain();
      const Float_t sumlo  = pix.GetExtractedSignalLoGain();
      
      histhi.FillHistAndArray(sumhi);
      histlo.FillHistAndArray(sumlo);

      const Int_t sathi = (Int_t)pix.GetNumHiGainSaturated();
      const Int_t satlo = (Int_t)pix.GetNumLoGainSaturated();

      histhi.SetSaturated(sathi); 
      histlo.SetSaturated(satlo); 

      const Int_t aidx   = (*fGeom)[i].GetAidx();
      const Int_t sector = (*fGeom)[i].GetSector();

      sumhiarea[aidx]  += sumhi;
      sumloarea[aidx]  += sumlo;
      sathiarea[aidx]  += sathi;
      satloarea[aidx]  += satlo;

      sumhisector[sector]  += sumhi;
      sumlosector[sector]  += sumlo;
      sathisector[sector]  += sathi;
      satlosector[sector]  += satlo;
    }

  MRawEvtPixelIter pixel(fRawEvt);
  while (pixel.Next())
    {
      
      const UInt_t pixid = pixel.GetPixelId();

      MHCalibrationChargePix &histhi = (MHCalibrationChargePix&)(*this)[pixid];
      MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)(*this)(pixid);

      if (histhi.IsExcluded())
         continue;
      
      const Float_t timehi = (Float_t)pixel.GetIdxMaxHiGainSample();
      const Float_t timelo = (Float_t)pixel.GetIdxMaxLoGainSample(lofirst);

      histhi.FillAbsTime(timehi);
      histlo.FillAbsTime(timelo);

      const Int_t aidx   = (*fGeom)[pixid].GetAidx();
      const Int_t sector = (*fGeom)[pixid].GetSector();

      timehiarea[aidx] += timehi;
      timeloarea[aidx] += timelo;

      timehisector[sector] += timehi;
      timelosector[sector] += timelo;
    }
  
  for (UInt_t j=0; j<nareas; j++)
    {

      const Int_t npix = fAverageAreaNum[j];

      MHCalibrationChargePix &hipix = (MHCalibrationChargePix&)GetAverageHiGainArea(j);
      MHCalibrationChargePix &lopix = (MHCalibrationChargePix&)GetAverageLoGainArea(j);

      hipix.FillHistAndArray(sumhiarea[j]/npix);
      lopix.FillHistAndArray(sumloarea[j]/npix);

      hipix.SetSaturated((Float_t)sathiarea[j]/npix); 
      lopix.SetSaturated((Float_t)satloarea[j]/npix); 

      hipix.FillAbsTime(timehiarea[j]/npix);
      lopix.FillAbsTime(timeloarea[j]/npix);

    }

  for (UInt_t j=0; j<nsectors; j++)
    {

      const Int_t npix = fAverageSectorNum[j];

      MHCalibrationChargePix &hipix = (MHCalibrationChargePix&)GetAverageHiGainSector(j);
      MHCalibrationChargePix &lopix = (MHCalibrationChargePix&)GetAverageLoGainSector(j);

      hipix.FillHistAndArray(sumhisector[j]/npix);
      lopix.FillHistAndArray(sumlosector[j]/npix);

      hipix.SetSaturated((Float_t)sathisector[j]/npix); 
      lopix.SetSaturated((Float_t)satlosector[j]/npix); 

      hipix.FillAbsTime(timehisector[j]/npix);
      lopix.FillAbsTime(timelosector[j]/npix);

    }

  return kTRUE;
}

// --------------------------------------------------------------------------
//
// For all TObjArray's (including the averaged ones), the following steps are performed: 
//
// 1) Returns if the pixel is excluded.
// 2) Tests saturation. In case yes, set the flag: MCalibrationPix::SetHiGainSaturation()
//    or the flag: MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainSaturated )
// 3) Store the absolute arrival times in the MCalibrationChargePix's. If flag 
//    MCalibrationPix::IsHiGainSaturation() is set, the Low-Gain arrival times are stored, 
//    otherwise the Hi-Gain ones.
// 4) Calls to MHCalibrationCam::FitHiGainArrays() and MCalibrationCam::FitLoGainArrays() 
//    with the flags:
//    - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainNotFitted )
//    - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainNotFitted )
//    - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainOscillating )
//    - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainOscillating )
//
Bool_t MHCalibrationChargeCam::FinalizeHists()
{

  for (Int_t i=0; i<fHiGainArray->GetSize(); i++)
    {
      
      MHCalibrationChargePix &histhi = (MHCalibrationChargePix&)(*this)[i];
      MCalibrationChargePix  &pix    = (MCalibrationChargePix&)(*fCam)[i];
      MBadPixelsPix          &bad    = (*fBadPixels)[i];
      
      if (histhi.IsExcluded())
        continue;
      
      if (histhi.GetSaturated() > fNumHiGainSaturationLimit*histhi.GetHGausHist()->GetEntries())
        {
          pix.SetHiGainSaturation();
          histhi.CreateFourierSpectrum();
          continue;
        }

      FinalizeAbsTimes(histhi, pix, bad, fFirstHiGain, fLastHiGain);
    }

  for (Int_t i=0; i<fLoGainArray->GetSize(); i++)
    {
      
      MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)(*this)(i);
      MBadPixelsPix          &bad    = (*fBadPixels)[i];

      if (histlo.IsExcluded())
        continue;
      
      if (histlo.GetSaturated() > fNumLoGainSaturationLimit*histlo.GetHGausHist()->GetEntries())
        {
          *fLog << warn << "Saturated Lo Gain histogram in pixel: " << i << endl;
          bad.SetUncalibrated( MBadPixelsPix::kLoGainSaturation ); 
          histlo.CreateFourierSpectrum();
          continue;
        }
  
      MCalibrationChargePix &pix    = (MCalibrationChargePix&)(*fCam)[i];
      
      if (pix.IsHiGainSaturation())
        FinalizeAbsTimes(histlo, pix, bad, fFirstLoGain, fLastLoGain);
    }

  for (Int_t j=0; j<fAverageHiGainAreas->GetSize(); j++)
    {
      
      MHCalibrationChargePix &histhi = (MHCalibrationChargePix&)GetAverageHiGainArea(j);      
      MCalibrationChargePix  &pix    = (MCalibrationChargePix&)fCam->GetAverageArea(j);
      MBadPixelsPix          &bad    = fCam->GetAverageBadArea(j);
      
      if (histhi.GetSaturated() > fNumHiGainSaturationLimit*histhi.GetHGausHist()->GetEntries())
        {
          pix.SetHiGainSaturation();
          histhi.CreateFourierSpectrum();
          continue;
        }

      FinalizeAbsTimes(histhi, pix, bad, fFirstHiGain, fLastHiGain);
   }
  
  for (Int_t j=0; j<fAverageLoGainAreas->GetSize(); j++)
    {
      
      MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)GetAverageLoGainArea(j);      
      MCalibrationChargePix  &pix    = (MCalibrationChargePix&)fCam->GetAverageArea(j);
      MBadPixelsPix          &bad    = fCam->GetAverageBadArea(j);      

      if (histlo.GetSaturated() > fNumLoGainSaturationLimit*histlo.GetHGausHist()->GetEntries())
        {
          *fLog << warn << "Saturated Lo Gain histogram in area idx: " << j << endl;
          histlo.CreateFourierSpectrum();
          continue;
        }

      if (pix.IsHiGainSaturation())
        FinalizeAbsTimes(histlo, pix, bad, fFirstLoGain, fLastLoGain);
    }

  for (Int_t j=0; j<fAverageHiGainSectors->GetSize(); j++)
    {
      
      MHCalibrationChargePix &histhi = (MHCalibrationChargePix&)GetAverageHiGainSector(j);      
      MCalibrationChargePix  &pix    = (MCalibrationChargePix&)fCam->GetAverageSector(j);
      MBadPixelsPix          &bad    = fCam->GetAverageBadSector(j);      

      if (histhi.GetSaturated() > fNumHiGainSaturationLimit*histhi.GetHGausHist()->GetEntries())
        {
          pix.SetHiGainSaturation();
          histhi.CreateFourierSpectrum();
          continue;
        }

      FinalizeAbsTimes(histhi, pix, bad, fFirstHiGain, fLastHiGain);
    }
  
  for (Int_t j=0; j<fAverageLoGainSectors->GetSize(); j++)
    {
      
      MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)GetAverageLoGainSector(j);      
      MCalibrationChargePix  &pix    = (MCalibrationChargePix&)fCam->GetAverageSector(j);
      MBadPixelsPix          &bad    = fCam->GetAverageBadSector(j);        
      
      if (histlo.GetSaturated() > fNumLoGainSaturationLimit*histlo.GetHGausHist()->GetEntries())
        {
          *fLog << warn << "Saturated Lo Gain histogram in sector: " << j << endl;
          bad.SetUncalibrated( MBadPixelsPix::kLoGainSaturation ); 
          histlo.CreateFourierSpectrum();
          continue;
        }

      if (pix.IsHiGainSaturation())
        FinalizeAbsTimes(histlo, pix, bad, fFirstLoGain, fLastLoGain);
    }
  
  //
  // Perform the fitting for the High Gain (done in MHCalibrationCam)
  //
  FitHiGainArrays((MCalibrationCam&)(*fCam),(*fBadPixels),
                  MBadPixelsPix::kHiGainNotFitted,
                  MBadPixelsPix::kHiGainOscillating);
  //
  // Perform the fitting for the Low Gain (done in MHCalibrationCam)
  //
  FitLoGainArrays((MCalibrationCam&)(*fCam),(*fBadPixels),
                  MBadPixelsPix::kLoGainNotFitted,
                  MBadPixelsPix::kLoGainOscillating);
      
  return kTRUE;
}

// --------------------------------------------------------------------------------
//
// Fill the absolute time results into MCalibrationChargePix
//
// Check absolute time validity:
// - Mean arrival time is at least fTimeLowerLimit slices from the lower edge 
// - Mean arrival time is at least fUpperLimit     slices from the upper edge
//
void MHCalibrationChargeCam::FinalizeAbsTimes(MHCalibrationChargePix &hist, MCalibrationChargePix &pix, MBadPixelsPix &bad, 
                                              Byte_t first, Byte_t last)
{
  
  const Float_t mean = hist.GetAbsTimeMean();
  const Float_t rms  = hist.GetAbsTimeRms();

  pix.SetAbsTimeMean ( mean );
  pix.SetAbsTimeRms  ( rms  );
  
  const Float_t lowerlimit = (Float_t)first + fTimeLowerLimit;
  const Float_t upperlimit = (Float_t)last  + fTimeUpperLimit;  

  if ( mean < lowerlimit)
    {
      *fLog << warn << GetDescriptor() 
            << Form("%s%3.1f%s%2.1f%s%4i",": Mean ArrivalTime: ",mean," smaller than ",fTimeLowerLimit,
                    " FADC slices from lower edge in pixel ",hist.GetPixId()) << endl;
      bad.SetUncalibrated( MBadPixelsPix::kMeanTimeInFirstBin );
    }
  
  if ( mean  > upperlimit )
    {
      *fLog << warn << GetDescriptor() 
            << Form("%s%3.1f%s%2.1f%s%4i",": Mean ArrivalTime: ",mean," greater than ",fTimeUpperLimit,
                    " FADC slices from upper edge in pixel ",hist.GetPixId()) << endl;
      bad.SetUncalibrated( MBadPixelsPix::kMeanTimeInLast2Bins );
    }
}

// --------------------------------------------------------------------------
//
// Sets all pixels to MBadPixelsPix::kUnsuitableRun, if following flags are set:
// - MBadPixelsPix::kLoGainSaturation
//
// Sets all pixels to MBadPixelsPix::kUnreliableRun, if following flags are set:
// - if MBadPixelsPix::kHiGainNotFitted   and !MCalibrationPix::IsHiGainSaturation()
// - if MBadPixelsPix::kHiGainOscillating and !MCalibrationPix::IsHiGainSaturation()
// - if MBadPixelsPix::kLoGainNotFitted   and  MCalibrationPix::IsLoGainSaturation()
// - if MBadPixelsPix::kLoGainOscillating and  MCalibrationPix::IsLoGainSaturation()
//
void MHCalibrationChargeCam::FinalizeBadPixels()
{
      
  for (Int_t i=0; i<fBadPixels->GetSize(); i++)
    {
      
      MBadPixelsPix    &bad    = (*fBadPixels)[i];
      MCalibrationPix  &pix    = (*fCam)[i];

      if (bad.IsUncalibrated( MBadPixelsPix::kHiGainNotFitted ))
        if (!pix.IsHiGainSaturation())
          bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun    );
 
      if (bad.IsUncalibrated( MBadPixelsPix::kHiGainOscillating ))
        bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun    );

      if (bad.IsUncalibrated( MBadPixelsPix::kLoGainNotFitted ))
        if (pix.IsHiGainSaturation())
          bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun    );
 
      if (bad.IsUncalibrated( MBadPixelsPix::kLoGainOscillating ))
        if (pix.IsHiGainSaturation())
          bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun    );

      if (bad.IsUncalibrated( MBadPixelsPix::kLoGainSaturation ))
          bad.SetUnsuitable(   MBadPixelsPix::kUnsuitableRun    );
    }
}

// --------------------------------------------------------------------------
//
// Dummy, needed by MCamEvent
//
Bool_t MHCalibrationChargeCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
{
  return kTRUE;
}

// --------------------------------------------------------------------------
//
// Calls MHGausEvents::DrawClone() for pixel idx
//
void MHCalibrationChargeCam::DrawPixelContent(Int_t idx) const
{
  (*this)[idx].DrawClone();
}