/* ======================================================================== *\
!
! *
! * 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 fHiGainArray and fLoGainArray
//
// - average pixel per area index (e.g. inner and outer for the MAGIC camera), 
//   stored in the TObjArray's fAverageHiGainAreas and fAverageLoGainAreas
//
// - average pixel per camera sector (e.g. sectors 1-6 for the MAGIC camera), 
//   stored in the TObjArray's fAverageHiGainSectors and fAverageLoGainSectors
// 
// Every signal is filled into a Charge 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 TObjArray's fAverageHiGainAreas, 
// fAverageLoGainAreas, fAverageHiGainSectors and fAverageLoGainSectors.
//
// 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 fProbLimit (default: 0.5%), the fit is valid.
// Otherwise, the fit is repeated within ranges of the previous mean +- 5 sigma. 
// In case this does not make the fit valid, the histogram means and RMS's are 
// taken directly and the following flags are set:
// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainNotFitted ) or  
// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainNotFitted ) and 
// - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnreliableRun   ) 
// 
// Outliers of more than MHCalibrationChargePix::fPickupLimit (default: 5) sigmas 
// from the mean are counted as PickUp events (stored in MHCalibrationChargePix::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 
// 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:
//
// - MCalibrationChargePix::SetHiGainSaturation(Bool_t b) 
// - MCalibrationChargePix::SetHiGainMeanCharge()
// - MCalibrationChargePix::SetHiGainMeanChargeErr()
// - MCalibrationChargePix::SetHiGainSigmaCharge()
// - MCalibrationChargePix::SetHiGainSigmaChargeErr()
// - MCalibrationChargePix::SetHiGainChargeProb()
// - MCalibrationChargePix::SetHiGainNumPickup()
//
// - MCalibrationChargePix::SetLoGainMeanCharge()
// - MCalibrationChargePix::SetLoGainMeanChargeErr()
// - MCalibrationChargePix::SetLoGainSigmaCharge()
// - MCalibrationChargePix::SetLoGainSigmaChargeErr()
// - MCalibrationChargePix::SetLoGainChargeProb()
// - MCalibrationChargePix::SetLoGainNumPickup()
//
// - MCalibrationChargePix::SetAbsTimeMean()
// - MCalibrationChargePix::SetAbsTimeRms()
//
// The last two items are filled with the Hi-Gain or Lo-Gain result, depending on 
// the saturation status. 
//
// 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 <TVirtualPad.h>
#include <TCanvas.h>
#include <TPad.h>
#include <TText.h>
#include <TPaveText.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 "MBadPixelsCam.h"
#include "MBadPixelsPix.h"

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

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

ClassImp(MHCalibrationChargeCam);

using namespace std;

const Int_t   MHCalibrationChargeCam::fgAverageNbins             = 4000;
const Float_t MHCalibrationChargeCam::fgNumHiGainSaturationLimit = 0.01;
const Float_t MHCalibrationChargeCam::fgNumLoGainSaturationLimit = 0.005;
const Int_t   MHCalibrationChargeCam::fgPulserFrequency          = 500;
// --------------------------------------------------------------------------
//
// Default Constructor. 
//
// Sets:
// - all pointers to NULL
//
// Initializes and sets owner of:
// - fHiGainArray, fLoGainArray
// - fAverageHiGainAreas, fAverageLoGainAreas
// - fAverageHiGainSectors, fAverageLoGainSectors
//
// Initializes:
// - fNumHiGainSaturationLimit to fgNumHiGainSaturationLimit
// - fNumLoGainSaturationLimit to fgNumLoGainSaturationLimit
// - fPulserFrequency to fgPulserFrequency
//
MHCalibrationChargeCam::MHCalibrationChargeCam(const char *name, const char *title)
    : fCam(NULL), fRawEvt(NULL), fGeom(NULL), fBadPixels(NULL)
{
    fName  = name  ? name  : "MHCalibrationChargeCam";
    fTitle = title ? title : "Class to fill the calibration histograms ";

    fHiGainArray = new TObjArray;
    fHiGainArray->SetOwner();
    
    fLoGainArray = new TObjArray;
    fLoGainArray->SetOwner();

    fAverageHiGainAreas = new TObjArray;
    fAverageHiGainAreas->SetOwner();

    fAverageLoGainAreas = new TObjArray;
    fAverageLoGainAreas->SetOwner();

    fAverageHiGainSectors = new TObjArray;
    fAverageHiGainSectors->SetOwner();

    fAverageLoGainSectors = new TObjArray;
    fAverageLoGainSectors->SetOwner();

    SetNumHiGainSaturationLimit();
    SetNumLoGainSaturationLimit();
    SetPulserFrequency();
}

// --------------------------------------------------------------------------
//
// Deletes the TClonesArray of:
// - fHiGainArray, fLoGainArray
// - fAverageHiGainAreas, fAverageLoGainAreas
// - fAverageHiGainSectors, fAverageLoGainSectors
//
MHCalibrationChargeCam::~MHCalibrationChargeCam()
{
  delete fHiGainArray;
  delete fLoGainArray;

  delete fAverageHiGainAreas;
  delete fAverageLoGainAreas;

  delete fAverageHiGainSectors;
  delete fAverageLoGainSectors;
}

// --------------------------------------------------------------------------
//
// Get i-th High Gain pixel (pixel number)
//
MHCalibrationChargeHiGainPix &MHCalibrationChargeCam::operator[](UInt_t i)
{
  return *static_cast<MHCalibrationChargeHiGainPix*>(fHiGainArray->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Get i-th High Gain pixel (pixel number)
//
const MHCalibrationChargeHiGainPix &MHCalibrationChargeCam::operator[](UInt_t i) const
{
  return *static_cast<MHCalibrationChargeHiGainPix*>(fHiGainArray->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Get i-th Low Gain pixel (pixel number)
//
MHCalibrationChargeLoGainPix &MHCalibrationChargeCam::operator()(UInt_t i)
{
  return *static_cast<MHCalibrationChargeLoGainPix*>(fLoGainArray->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Get i-th Low Gain pixel (pixel number)
//
const MHCalibrationChargeLoGainPix &MHCalibrationChargeCam::operator()(UInt_t i) const
{
  return *static_cast<MHCalibrationChargeLoGainPix*>(fLoGainArray->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Get i-th High Gain pixel Area (area number)
//
MHCalibrationChargeHiGainPix &MHCalibrationChargeCam::GetAverageHiGainArea(UInt_t i)
{
  return *static_cast<MHCalibrationChargeHiGainPix*>(fAverageHiGainAreas->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Get i-th High Gain pixel Area (area number)
//
const MHCalibrationChargeHiGainPix &MHCalibrationChargeCam::GetAverageHiGainArea(UInt_t i) const
{
  return *static_cast<MHCalibrationChargeHiGainPix*>(fAverageHiGainAreas->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Get i-th Low Gain pixel Area (area number)
//
MHCalibrationChargeLoGainPix &MHCalibrationChargeCam::GetAverageLoGainArea(UInt_t i)
{
  return *static_cast<MHCalibrationChargeLoGainPix*>(fAverageLoGainAreas->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Get i-th Low Gain pixel Area (area number)
//
const MHCalibrationChargeLoGainPix &MHCalibrationChargeCam::GetAverageLoGainArea(UInt_t i) const
{
  return *static_cast<MHCalibrationChargeLoGainPix*>(fAverageLoGainAreas->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Get i-th High Gain Sector (sector number)
//
MHCalibrationChargeHiGainPix &MHCalibrationChargeCam::GetAverageHiGainSector(UInt_t i)
{
  return *static_cast<MHCalibrationChargeHiGainPix*>(fAverageHiGainSectors->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Get i-th High Gain Sector (sector number)
//
const MHCalibrationChargeHiGainPix &MHCalibrationChargeCam::GetAverageHiGainSector(UInt_t i) const
{
  return *static_cast<MHCalibrationChargeHiGainPix*>(fAverageHiGainSectors->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Get i-th Low Gain Sector (sector number)
//
MHCalibrationChargeLoGainPix &MHCalibrationChargeCam::GetAverageLoGainSector(UInt_t i)
{
  return *static_cast<MHCalibrationChargeLoGainPix*>(fAverageLoGainSectors->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Get i-th Low Gain Sector (sector number)
//
const MHCalibrationChargeLoGainPix &MHCalibrationChargeCam::GetAverageLoGainSector(UInt_t i) const
{
  return *static_cast<MHCalibrationChargeLoGainPix*>(fAverageLoGainSectors->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Our own clone function is necessary since root 3.01/06 or Mars 0.4
// I don't know the reason. 
//
// Creates new MHCalibrationChargeCam
// Deletes the TObjArray's and Clones them individually
// Copies the TArray's
//
TObject *MHCalibrationChargeCam::Clone(const char *) const
{

  const Int_t nhi   = fHiGainArray->GetEntries();
  const Int_t nlo   = fLoGainArray->GetEntries();
  const Int_t navhi = fAverageHiGainAreas->GetEntries();
  const Int_t navlo = fAverageLoGainAreas->GetEntries();
  const Int_t nsehi = fAverageHiGainSectors->GetEntries();
  const Int_t nselo = fAverageLoGainSectors->GetEntries();
  
  //
  // FIXME, this might be done faster and more elegant, by direct copy.
  //
  MHCalibrationChargeCam *cam = new MHCalibrationChargeCam();

  cam->fHiGainArray->Expand(nhi);
  cam->fLoGainArray->Expand(nlo);
  cam->fAverageHiGainAreas->Expand(navhi);
  cam->fAverageLoGainAreas->Expand(navlo);
  cam->fAverageHiGainSectors->Expand(nsehi);
  cam->fAverageLoGainSectors->Expand(nselo);

  for (int i=0; i<nhi; i++)
    {
      delete (*cam->fHiGainArray)[i];
      (*cam->fHiGainArray)[i] = (MHCalibrationChargeHiGainPix*)(*fHiGainArray)[i]->Clone();
    }
  for (int i=0; i<nlo; i++)
    {
      delete (*cam->fLoGainArray)[i];
      (*cam->fLoGainArray)[i] = (MHCalibrationChargeLoGainPix*)(*fLoGainArray)[i]->Clone();
    }
  for (int i=0; i<navhi; i++)
    {
      delete (*cam->fAverageHiGainAreas)[i];
      (*cam->fAverageHiGainAreas)[i] = (MHCalibrationChargeHiGainPix*)(*fAverageHiGainAreas)[i]->Clone();
    }
  for (int i=0; i<navlo; i++)
    {
      delete (*cam->fAverageLoGainAreas)[i];
      (*cam->fAverageLoGainAreas)[i] = (MHCalibrationChargeLoGainPix*)(*fAverageLoGainAreas)[i]->Clone();
    }
  for (int i=0; i<nsehi; i++)
    {
      delete (*cam->fAverageHiGainSectors)[i];
      (*cam->fAverageHiGainSectors)[i] = (MHCalibrationChargeHiGainPix*)(*fAverageHiGainSectors)[i]->Clone();
    }
  for (int i=0; i<nselo; i++)
    {
      delete (*cam->fAverageLoGainSectors)[i];
      (*cam->fAverageLoGainSectors)[i] = (MHCalibrationChargeLoGainPix*)(*fAverageLoGainSectors)[i]->Clone();
    }

  cam->fAverageAreaNum         = fAverageAreaNum;
  cam->fAverageAreaSat         = fAverageAreaSat;
  cam->fAverageAreaSigma       = fAverageAreaSigma;      
  cam->fAverageAreaSigmaErr    = fAverageAreaSigmaErr;   
  cam->fAverageAreaRelSigma    = fAverageAreaRelSigma;
  cam->fAverageAreaRelSigmaErr = fAverageAreaRelSigmaErr;   
  cam->fAverageSectorNum       = fAverageSectorNum;      

  return cam;
}

// --------------------------------------------------------------------------
//
// Gets the pointers to:
// - MRawEvtData
// - MGeomCam
//
// Calls Delete-Function of:
// - fHiGainArray, fLoGainArray
// - fAverageHiGainAreas, fAverageLoGainAreas
// - fAverageHiGainSectors, fAverageLoGainSectors
//
Bool_t MHCalibrationChargeCam::SetupFill(const MParList *pList)
{
  
  fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
  if (!fRawEvt)
  {
      *fLog << err << dbginf << "MRawEvtData not found... aborting." << endl;
      return kFALSE;
  }

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

  fHiGainArray->Delete();
  fLoGainArray->Delete();

  fAverageHiGainAreas->Delete();
  fAverageLoGainAreas->Delete();

  fAverageHiGainSectors->Delete();
  fAverageLoGainSectors->Delete();

  return kTRUE;
}

// --------------------------------------------------------------------------
//
// Gets or creates the pointers to:
// - MBadPixelsCam
// - MCalibrationChargeCam
//
// Searches pointer to:
// - MExtractedSignalCam
//
// Initializes, if empty to MExtractedSignalCam::GetSize() for:
// - fHiGainArray, fLoGainArray
//
// Initializes, if empty to MGeomCam::GetNumAreas() for:
// - fAverageHiGainAreas, fAverageLoGainAreas
//
// Initializes, if empty to MGeomCam::GetNumSectors() for:
// - fAverageHiGainSectors, fAverageLoGainSectors
// 
// Initializes TArrays to MGeomCam::GetNumAreas and MGeomCam::GetNumSectors, respectively
// Fills with number of valid pixels (if !MBadPixelsPix::IsBad()):
// - fAverageAreaNum[area index]
// - fAverageSectorNum[area index]
//
// Calls InitializeHiGainHists() and InitializeLoGainHists() for every entry in:
// - fHiGainArray, fLoGainArray
// - fAverageHiGainAreas, fAverageLoGainAreas
// - fAverageHiGainSectors, fAverageLoGainSectors
//
// Sets Titles and Names for the Charge Histograms and 
// Sets number of bins to fAverageNbins for:
// - fAverageHiGainAreas, fAverageLoGainAreas
// - fAverageHiGainSectors, fAverageLoGainSectors
// 
Bool_t MHCalibrationChargeCam::ReInit(MParList *pList)
{

  fBadPixels = (MBadPixelsCam*)pList->FindCreateObj("MBadPixelsCam");
  if (!fBadPixels)
      return kFALSE;


  fCam = (MCalibrationChargeCam*)pList->FindCreateObj("MCalibrationChargeCam");
  if (!fCam)
      return kFALSE;

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

  const Int_t n        = signal->GetSize();
  const Int_t nsectors = fGeom->GetNumSectors();
  const Int_t nareas   = fGeom->GetNumAreas();

  //
  // The function TArrayF::Set() already sets all entries to 0.
  //
  fAverageAreaNum.        Set(nareas);
  fAverageAreaSat.        Set(nareas);           
  fAverageAreaSigma.      Set(nareas);      
  fAverageAreaSigmaErr.   Set(nareas);   
  fAverageAreaRelSigma.   Set(nareas);   
  fAverageAreaRelSigmaErr.Set(nareas);
  fAverageSectorNum.      Set(nsectors);

  for (Int_t i=0; i<n; i++)
    {
      if ((*fBadPixels)[i].IsBad())
        continue;
      fAverageAreaNum  [(*fGeom)[i].GetAidx()  ]++;
      fAverageSectorNum[(*fGeom)[i].GetSector()]++;
    }

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

 
  if (fLoGainArray->GetEntries()==0)
  {
      fLoGainArray->Expand(n);
      
      for (Int_t i=0; i<n; i++)
      {
	  (*fLoGainArray)[i] = new MHCalibrationChargeLoGainPix;
          InitializeLoGainHists((*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 ");
        InitializeHiGainHists(GetAverageHiGainArea(j),fCam->GetAverageBadArea(j),j);        
        GetAverageHiGainArea(j).GetHGausHist()->SetTitle("Summed FADC slices average HiGain Area Idx ");
        GetAverageHiGainArea(j).GetHAbsTime()->SetTitle("Absolute Arrival Time average HiGain Area Idx ");
        GetAverageHiGainArea(j).SetChargeNbins(fAverageNbins);
      }
  }

  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 ");
          InitializeLoGainHists(GetAverageLoGainArea(j),fCam->GetAverageBadArea(j),j);        
          GetAverageLoGainArea(j).GetHGausHist()->SetTitle("Summed FADC slices average LoGain Area Idx ");
          GetAverageLoGainArea(j).GetHAbsTime()->SetTitle("Absolute Arrival Time average LoGain Area Idx ");
          GetAverageLoGainArea(j).SetChargeNbins(fAverageNbins);
        }
    }
  
  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 ");
          InitializeHiGainHists(GetAverageHiGainSector(j),fCam->GetAverageBadSector(j),j);        
          GetAverageHiGainSector(j).GetHGausHist()->SetTitle("Summed FADC slices average HiGain Sector ");
          GetAverageHiGainSector(j).GetHAbsTime()->SetTitle("Absolute Arrival Time average HiGain Sector ");
          GetAverageHiGainSector(j).SetChargeNbins(fAverageNbins);
      }
  }

  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 ");
          InitializeLoGainHists(GetAverageLoGainSector(j),fCam->GetAverageBadSector(j),j); 
          GetAverageLoGainSector(j).GetHGausHist()->SetTitle("Summed FADC slices average LoGain Sector ");
          GetAverageLoGainSector(j).GetHAbsTime()->SetTitle("Absolute Arrival Time average LoGain Sector ");
          GetAverageLoGainSector(j).SetChargeNbins(fAverageNbins);
      }
  }

  return kTRUE;
}

// -------------------------------------------------------------
//
// If MBadPixelsPix::IsBad():
// - calls MHCalibrationChargeHiGainPix::SetExcluded()
//
// Calls:
// - MHCalibrationChargePix::Init()
// - MHCalibrationChargePix::ChangeHistId(i)
// - MHGausEvents::SetEventFrequency(fPulserFrequency)
// 
void MHCalibrationChargeCam::InitializeHiGainHists(MHCalibrationChargeHiGainPix &hist,MBadPixelsPix &bad, const Int_t i)
{
  
  if (bad.IsBad())
    {
      *fLog << warn << "Excluded pixel: " << i << " from calibration " << endl;
      hist.SetExcluded();
    }
  
  hist.Init();
  hist.ChangeHistId(i);
  hist.SetEventFrequency(fPulserFrequency);
          
}

// -------------------------------------------------------------
//
// If MBadPixelsPix::IsBad():
// - calls MHCalibrationChargeLoGainPix::SetExcluded()
//
// Calls:
// - MHCalibrationChargePix::Init()
// - MHCalibrationChargePix::ChangeHistId(i)
// - MHGausEvents::SetEventFrequency(fPulserFrequency)
// 
void MHCalibrationChargeCam::InitializeLoGainHists(MHCalibrationChargeLoGainPix &hist,MBadPixelsPix &bad, const Int_t i)
{

  if (bad.IsBad())
    hist.SetExcluded();
  
  hist.Init();
  hist.ChangeHistId(i);
  hist.SetEventFrequency(fPulserFrequency);
}

  
// --------------------------------------------------------------------------
//
// Retrieves from MExtractedSignalCam:
// - number of pixels
// - first used LoGain FADC slice
//
// Retrieves from MGeomCam:
// - number of pixel areas
// - number of sectors
//
// For all TObjArray's (including the averaged ones), the following steps are performed: 
//
// 1) Test size and return kFALSE if not matching
//
// 2) Fill Charges histograms (MHGausEvents::FillHistAndArray()) with:
// - MExtractedSignalPix::GetExtractedSignalHiGain();
// - MExtractedSignalPix::GetExtractedSignalLoGain();
//
// 3) Set number of saturated slices (MHCalibrationChargePix::SetSaturated()) with:
// - MExtractedSignalPix::GetNumHiGainSaturated();
// - MExtractedSignalPix::GetNumLoGainSaturated();
//
// 4) Fill AbsTime histograms (MHCalibrationChargePix::FillAbsTime()) with:
// - MRawEvtPixelIter::GetIdxMaxHiGainSample();       
// - MRawEvtPixelIter::GetIdxMaxLoGainSample(first slice);
//
Bool_t MHCalibrationChargeCam::Fill(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 Int_t npixels  = signal->GetSize();
  const Int_t lofirst  = signal->GetFirstUsedSliceLoGain();
  const Int_t nareas   = fGeom->GetNumAreas();
  const Int_t nsectors = fGeom->GetNumSectors();

  if (fHiGainArray->GetEntries() != npixels)
    {
      *fLog << err << "ERROR - Size mismatch in number of pixels ... abort." << endl;
      return kFALSE;
    }

  if (fLoGainArray->GetEntries() != npixels)
    {
      *fLog << err << "ERROR - Size mismatch in number of pixels ... abort." << endl;
      return kFALSE;
    }

  if (fAverageHiGainAreas->GetEntries() != nareas)
    {
      *fLog << err << "ERROR - Size mismatch in number of areas ... abort." << endl;
      return kFALSE;
    }

  if (fAverageLoGainAreas->GetEntries() != nareas)
    {
      *fLog << err << "ERROR - Size mismatch in number of areas ... abort." << endl;
      return kFALSE;
    }

  if (fAverageHiGainSectors->GetEntries() != nsectors)
    {
      *fLog << err << "ERROR - Size mismatch in number of sectors ... abort." << endl;
      return kFALSE;
    }

  if (fAverageLoGainSectors->GetEntries() != nsectors)
    {
      *fLog << err << "ERROR - Size mismatch in number of sectors ... abort." << endl;
      return kFALSE;
    }

  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];

  for (UInt_t j=0; j<nareas; j++)
    {
      sumhiarea  [j] = sumloarea  [j] = timehiarea  [j] =  timeloarea  [j] = 0.;
      sumhisector[j] = sumlosector[j] = timehisector[j] =  timelosector[j] = 0.;
      sathiarea  [j] = satloarea  [j] = 0;
      sathisector[j] = satlosector[j] = 0;
    }
  
  for (Int_t i=0; i<npixels; i++)
    {

      if ((*this)[i].IsExcluded())
	continue;

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

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

      (*this)[i].SetSaturated(sathi); 
      (*this)(i).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();
       if ((*this)[pixid].IsExcluded())
         continue;

      const Float_t timehi = (Float_t)pixel.GetIdxMaxHiGainSample();
      const Float_t timelo = (Float_t)pixel.GetIdxMaxLoGainSample(lofirst);

      (*this)[pixid].FillAbsTime(timehi);
      (*this)(pixid).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];

      MHCalibrationChargeHiGainPix &hipix = GetAverageHiGainArea(j);
      MHCalibrationChargeLoGainPix &lopix = 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];

      MHCalibrationChargeHiGainPix &hipix = GetAverageHiGainSector(j);
      MHCalibrationChargeLoGainPix &lopix = 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) Return if the pixel is excluded
// 
// 2) Call to FinalizeHiGainHists() and FinalizeLoGainHists()
//
// For all averaged areas, the fitted sigma is multiplied with the square root of 
// the number involved pixels
//
Bool_t MHCalibrationChargeCam::Finalize()
{

  for (Int_t i=0; i<fHiGainArray->GetSize(); i++)
    {
      
      if ((*this)[i].IsExcluded())
        continue;
      
      MHCalibrationChargeHiGainPix &histhi = (*this)[i];
      MCalibrationChargePix        &pix    = (*fCam)[i];
      MBadPixelsPix                &bad    = (*fBadPixels)[i];
      
      FinalizeHiGainHists(histhi,pix,bad);
      
    }
  
  for (Int_t i=0; i<fLoGainArray->GetSize(); i++)
    {
      
      if ((*this)(i).IsExcluded())
        continue;
      
      MHCalibrationChargeLoGainPix &histlo = (*this)(i);
      MCalibrationChargePix        &pix    = (*fCam)[i];
      MBadPixelsPix                &bad    = (*fBadPixels)[i];
      
      FinalizeLoGainHists(histlo,pix,bad);
      
    }
  
  for (Int_t j=0; j<fAverageHiGainAreas->GetSize(); j++)
    {
      
      MHCalibrationChargeHiGainPix &histhi = GetAverageHiGainArea(j);      
      MCalibrationChargePix        &pix    = fCam->GetAverageArea(j);
      MBadPixelsPix                &bad    = fCam->GetAverageBadArea(j);        
      
      FinalizeHiGainHists(histhi,pix,bad);
    }
  
  for (Int_t j=0; j<fAverageLoGainAreas->GetSize(); j++)
    {
      
      MHCalibrationChargeLoGainPix &histlo = GetAverageLoGainArea(j);      
      MCalibrationChargePix        &pix    = fCam->GetAverageArea(j);
      MBadPixelsPix                &bad    = fCam->GetAverageBadArea(j);        
      
      FinalizeLoGainHists(histlo,pix,bad);
    }

  for (Int_t j=0; j<fAverageHiGainSectors->GetSize(); j++)
    {
      
      MHCalibrationChargeHiGainPix &histhi = GetAverageHiGainSector(j);      
      MCalibrationChargePix        &pix    = fCam->GetAverageSector(j);
      MBadPixelsPix                &bad    = fCam->GetAverageBadSector(j);        
      
      FinalizeHiGainHists(histhi,pix,bad);
    }
  
  for (Int_t j=0; j<fAverageLoGainSectors->GetSize(); j++)
    {
      
      MHCalibrationChargeLoGainPix &histlo = GetAverageLoGainSector(j);      
      MCalibrationChargePix        &pix    = fCam->GetAverageSector(j);
      MBadPixelsPix                &bad    = fCam->GetAverageBadSector(j);        
      
      FinalizeLoGainHists(histlo,pix,bad);
    }
  
  for (Int_t j=0; j<fGeom->GetNumAreas();j++)
    {
      
      MCalibrationChargePix &pix = fCam->GetAverageArea(j);

      if (pix.IsHiGainSaturation())
        fAverageAreaSat[j]++;

      fAverageAreaSigma[j]    = pix.GetSigmaCharge    () * TMath::Sqrt((Float_t)fAverageAreaNum[j]);
      fAverageAreaSigmaErr[j] = pix.GetSigmaChargeErr () * TMath::Sqrt((Float_t)fAverageAreaNum[j]);

      pix.SetSigmaCharge   (fAverageAreaSigma[j]);
      pix.SetSigmaChargeErr(fAverageAreaSigmaErr[j]);

      fAverageAreaRelSigma[j]   = fAverageAreaSigma[j] / pix.GetMeanCharge();
      
      Float_t relsigmaerr       =  fAverageAreaSigmaErr[j]*fAverageAreaSigmaErr[j] 
                                / (fAverageAreaSigma[j]   *fAverageAreaSigma[j]   );
      relsigmaerr               += pix.GetMeanChargeErr()*pix.GetMeanChargeErr() 
                                / (pix.GetMeanCharge()   *pix.GetMeanCharge()   );
      relsigmaerr               *= fAverageAreaRelSigma[j];
      fAverageAreaRelSigmaErr[j] = TMath::Sqrt(relsigmaerr);

    }

  return kTRUE;
}

// ---------------------------------------------------------------------------
//
// Returns if the histogram is empty and sets the following flag:
// - MBadPixelsPix::SetUnsuitable(MBadPixelsPix::kUnsuitableRun)
//
// Returns if the number of saturated slices (MHCalibrationChargePix::GetSaturated()) is 
// higher than the allowed limit (fNumHiGainSaturationLimit*histogram-entries), creates 
// the fourier spectrum and sets the following flags:
// - MCalibrationChargePix::SetHiGainSaturation()
//
// Fits the histograms with a Gaussian, in case of failure 
// calls MHCalibrationChargePix::RepeatFit(), in case of repeated failure 
// calls MHGausEvents::BypassFit() and sets the following flags:
// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainNotFitted )
// - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnreliableRun   )
// 
// Counts the number of pickup events
//
// Creates the fourier spectrum and tests MHGausEvents::IsFourierSpectrumOK(). 
// In case no, sets the following flags:
// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainOscillating )
// - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnreliableRun     )
//
// Retrieves the results and store them in MCalibrationChargePix
//
void MHCalibrationChargeCam::FinalizeHiGainHists(MHCalibrationChargeHiGainPix &hist, 
                                                 MCalibrationChargePix &pix, 
                                                 MBadPixelsPix &bad)
{

    if (hist.IsEmpty())
      {
	*fLog << warn << "Empty Hi Gain histogram in pixel: " << pix.GetPixId() << endl;
        bad.SetUnsuitable(MBadPixelsPix::kUnsuitableRun);        
	return;
      }
    
    if (hist.GetSaturated() > fNumHiGainSaturationLimit*hist.GetHGausHist()->GetEntries())
    {
      *fLog << warn << "Saturated Hi Gain histogram in pixel: " << pix.GetPixId() << endl;
      pix.SetHiGainSaturation();
      hist.CreateFourierSpectrum();
      return;
    }

    //
    // 2) Fit the Hi Gain histograms with a Gaussian
    //
    if (!hist.FitGaus())
    //
    // 3) In case of failure set the bit Fitted to false and take histogram means and RMS
    //
      if (!hist.RepeatFit())
        {
          hist.BypassFit();
	  *fLog << warn << "Hi Gain could not be fitted in pixel: " << pix.GetPixId() << endl;
          bad.SetUncalibrated( MBadPixelsPix::kHiGainNotFitted );
          bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun   );
        }

    //
    // 4) Check for pickup
    //
    hist.CountPickup();
    
    //
    // 5) Check for oscillations
    // 
    hist.CreateFourierSpectrum();
    
    if (!hist.IsFourierSpectrumOK())
      {
	*fLog << warn << "Oscillating Hi Gain in pixel: " << pix.GetPixId() << endl;
	bad.SetUncalibrated( MBadPixelsPix::kHiGainOscillating );
        bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun     );
      }

    //
    // 6) Retrieve the results and store them in this class
    //
    pix.SetHiGainMeanCharge(     hist.GetMean()      );
    pix.SetHiGainMeanChargeErr(  hist.GetMeanErr()   );
    pix.SetHiGainSigmaCharge(    hist.GetSigma()     );
    pix.SetHiGainSigmaChargeErr( hist.GetSigmaErr()  );
    pix.SetHiGainChargeProb    ( hist.GetProb()      );
    
    pix.SetAbsTimeMean         ( hist.GetAbsTimeMean());
    pix.SetAbsTimeRms          ( hist.GetAbsTimeRms() );
    
    pix.SetHiGainNumPickup     ( hist.GetPickup()     );
    
}

// ---------------------------------------------------------------------------
//
// Returns if the histogram is empty. If yes and the flag MCalibrationChargePix::IsHiGainSaturation()
// is set, sets the following flag:
// - MBadPixelsPix::SetUnsuitable(MBadPixelsPix::kUnsuitableRun)
//
// Returns if the number of saturated slices (MHCalibrationChargePix::GetSaturated()) is 
// higher than the allowed limit (fNumLoGainSaturationLimit*histogram-entries), creates 
// the fourier spectrum and sets the following flags:
// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainSaturation ); 
// - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnsuitableRun    );
//
// Fits the histograms with a Gaussian, in case of failure calls MHCalibrationChargePix::RepeatFit(), 
// in case of repeated failure calls MHGausEvents::BypassFit() and sets the following flags:
// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainNotFitted )
//   If the flag MCalibrationChargePix::IsHiGainSaturation() is set, sets additionally:
// - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnreliableRun   )
// 
// Counts the number of pickup events
//
// Creates the fourier spectrum and tests MHGausEvents::IsFourierSpectrumOK(). 
// In case no, sets the following flags:
// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainOscillating )
//   If the flag MCalibrationChargePix::IsHiGainSaturation() is set, sets additionally:
// - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnreliableRun     )
//
// Retrieves the results and store them in MCalibrationChargePix
//
void MHCalibrationChargeCam::FinalizeLoGainHists(MHCalibrationChargeLoGainPix &hist, 
                                                 MCalibrationChargePix &pix, 
                                                 MBadPixelsPix &bad)
{

    if (hist.IsEmpty())
      {
        if (pix.IsHiGainSaturation())
          bad.SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
        return;
      }

    if (hist.GetSaturated() > fNumLoGainSaturationLimit*hist.GetHGausHist()->GetEntries())
    {
      *fLog << warn << "Saturated Lo Gain histogram in pixel: " << pix.GetPixId() << endl;
      bad.SetUncalibrated( MBadPixelsPix::kLoGainSaturation ); 
      bad.SetUnsuitable(   MBadPixelsPix::kUnsuitableRun    );
      hist.CreateFourierSpectrum();
      return;
    }

    //
    // 2) Fit the Lo Gain histograms with a Gaussian
    //
    if (!hist.FitGaus())
    //
    // 3) In case of failure set the bit kFitted to false and take histogram means and RMS
    //
      if (!hist.RepeatFit())
        {
          hist.BypassFit();
	  bad.SetUncalibrated( MBadPixelsPix::kLoGainNotFitted );
          if (pix.IsHiGainSaturation())
            bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun);
        }

    //
    // 4) Check for pickup
    //
    hist.CountPickup();
    
    //
    // 5) Check for oscillations
    // 
    hist.CreateFourierSpectrum();
    
    if (!hist.IsFourierSpectrumOK())
      {
        bad.SetUncalibrated( MBadPixelsPix::kLoGainOscillating );
        if (pix.IsHiGainSaturation())
          bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun);
      }
    //
    // 6) Retrieve the results and store them in this class
    //
    pix.SetLoGainMeanCharge(     hist.GetMean()     );
    pix.SetLoGainMeanChargeErr(  hist.GetMeanErr()  );
    pix.SetLoGainSigmaCharge(    hist.GetSigma()    );
    pix.SetLoGainSigmaChargeErr( hist.GetSigmaErr() );
    pix.SetLoGainChargeProb    ( hist.GetProb()     );
    
    if (pix.IsHiGainSaturation())
    {
	pix.SetAbsTimeMean     ( hist.GetAbsTimeMean());
	pix.SetAbsTimeRms      ( hist.GetAbsTimeRms() );
    }	    
    
    pix.SetLoGainNumPickup     (  hist.GetPickup()    );

}    

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

// --------------------------------------------------------------------------
//
// What MHCamera needs in order to draw an individual pixel in the camera
//
void MHCalibrationChargeCam::DrawPixelContent(Int_t idx) const
{
  (*this)[idx].DrawClone();
}

// -----------------------------------------------------------------------------
// 
// Default draw:
//
// Displays the averaged areas, both High Gain and Low Gain 
//
// The following options can be chosen:
//
// "": displays the fHGausHist and the fHAbsTime
// "all": executes additionally MHGausEvents::Draw(), with options
//
void MHCalibrationChargeCam::Draw(const Option_t *opt)
{

  const Int_t nareas = fAverageHiGainAreas->GetEntries();
  if (nareas == 0)
    return;

  TVirtualPad *pad = gPad ? gPad : MH::MakeDefCanvas(this);  
  pad->SetBorderMode(0);

  pad->Divide(2,nareas);

  for (Int_t i=0; i<nareas;i++) 
    {
      pad->cd(2*(i+1)-1);
      GetAverageHiGainArea(i).Draw(opt);

      if (!fAverageAreaSat[i])
        DrawAverageSigma(fAverageAreaSat[i], i,
                         fAverageAreaSigma[i],    fAverageAreaSigmaErr[i],
                         fAverageAreaRelSigma[i], fAverageAreaRelSigmaErr[i]);

      pad->cd(2*(i+1));
      GetAverageLoGainArea(i).Draw(opt);
      
      if (fAverageAreaSat[i])
        DrawAverageSigma(fAverageAreaSat[i], i,
                         fAverageAreaSigma[i], fAverageAreaSigmaErr[i],
                         fAverageAreaRelSigma[i], fAverageAreaRelSigmaErr[i]);
    }
}

// -----------------------------------------------------------------------------
// 
// Default draw:
//
// Displays a TPaveText with the re-normalized sigmas of the average area
//
void MHCalibrationChargeCam::DrawAverageSigma(Bool_t sat, Bool_t inner,
                                              Float_t sigma, Float_t sigmaerr,
                                              Float_t relsigma, Float_t relsigmaerr) const 
{
  
  if (sigma != 0)
    {
      
      TPad *newpad = new TPad("newpad","transparent",0,0,1,1);
      newpad->SetFillStyle(4000);
      newpad->Draw();
      newpad->cd();
      
      TPaveText *text = new TPaveText(sat? 0.1 : 0.35,0.7,sat ? 0.4 : 0.7,1.0);
      text->SetTextSize(0.07);
      const TString line1 = Form("%s%s%s",inner ? "Outer" : "Inner",
                                 " Pixels ", sat ? "Low Gain" : "High Gain");
      TText *txt1 = text->AddText(line1.Data());
      const TString line2 = Form("Sigma per Pixel: %2.2f #pm %2.2f",sigma,sigmaerr);
      TText *txt2 = text->AddText(line2.Data());
      const TString line3 = Form("Rel. Sigma per Pixel: %2.2f #pm %2.2f",relsigma,relsigmaerr);
      TText *txt3 = text->AddText(line3.Data());
      text->Draw("");
      
      text->SetBit(kCanDelete);
      txt1->SetBit(kCanDelete);
      txt2->SetBit(kCanDelete);
      txt3->SetBit(kCanDelete);
      newpad->SetBit(kCanDelete);
    }
}

