/* ======================================================================== *\
!
! *
! * 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 MHCalibrationPix-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 
// +- MHCalibrationPix::fPickupLimit (default: 5) sigma (see MHCalibrationPix::RepeatFit())
// In case this does not make the fit valid, the histogram means and RMS's are 
// taken directly (see MHCalibrationPix::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 MHCalibrationPix::fPickupLimit (default: 5) sigmas 
// from the mean are counted as Pickup events (stored in MHCalibrationPix::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 "MHCalibrationChargePix.h"
#include "MHCalibrationPix.h"

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

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

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

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

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

#include "MArrayI.h"
#include "MArrayD.h"

#include <TPad.h>
#include <TVirtualPad.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TF1.h>
#include <TLine.h>
#include <TLatex.h>
#include <TLegend.h>
#include <TGraph.h>

ClassImp(MHCalibrationChargeCam);

using namespace std;

const Int_t   MHCalibrationChargeCam::fgChargeHiGainNbins =  550;
const Axis_t  MHCalibrationChargeCam::fgChargeHiGainFirst = -100.5;
const Axis_t  MHCalibrationChargeCam::fgChargeHiGainLast  = 999.5;
const Int_t   MHCalibrationChargeCam::fgChargeLoGainNbins =  325;
const Axis_t  MHCalibrationChargeCam::fgChargeLoGainFirst = -150.5;
const Axis_t  MHCalibrationChargeCam::fgChargeLoGainLast  =  499.5;
const TString MHCalibrationChargeCam::gsHistName          = "Charge";
const TString MHCalibrationChargeCam::gsHistTitle         = "Signals";
const TString MHCalibrationChargeCam::gsHistXTitle        = "Signal [FADC counts]";
const TString MHCalibrationChargeCam::gsHistYTitle        = "Nr. events";
const TString MHCalibrationChargeCam::gsAbsHistName       = "AbsTime";
const TString MHCalibrationChargeCam::gsAbsHistTitle      = "Abs. Arr. Times";
const TString MHCalibrationChargeCam::gsAbsHistXTitle     = "Time [FADC slices]";
const TString MHCalibrationChargeCam::gsAbsHistYTitle     = "Nr. events";
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.;
// 1Led Green, 1 LED blue, 5 LEDs blue, 10 LEDs blue, 10 LEDs UV, CT1, 5Leds Green
const Float_t MHCalibrationChargeCam::gkHiGainInnerRefLines[7]   = { 245., 323. , 1065., 1467., 180., 211. , 533.5};   
const Float_t MHCalibrationChargeCam::gkHiGainOuterRefLines[7]   = { 217., 307.5, 932. , 1405., 167., 183.5, 405.5};
const Float_t MHCalibrationChargeCam::gkLoGainInnerRefLines[7]   = { 20.8, 28.0 , 121. , 200.2, 16.5, 13.5 , 41.7 };
const Float_t MHCalibrationChargeCam::gkLoGainOuterRefLines[7]   = { 18.9, 26.0 , 108.3, 198. , 14.0, 11.  , 42.  };
// --------------------------------------------------------------------------
//
// Default Constructor. 
//
// Sets:
// - all pointers to NULL
//
// Initializes:
// - fNumHiGainSaturationLimit to fgNumHiGainSaturationLimit
// - fNumLoGainSaturationLimit to fgNumLoGainSaturationLimit
// - fTimeLowerLimit           to fgTimeLowerLimit 
// - fTimeUpperLimit           to fgTimeUpperLimit 
//
// - fNbins to fgChargeHiGainNbins
// - fFirst to fgChargeHiGainFirst
// - fLast  to fgChargeHiGainLast 
//
// - fLoGainNbins to fgChargeLoGainNbins
// - fLoGainFirst to fgChargeLoGainFirst
// - fLoGainLast  to fgChargeLoGainLast 
//
// - fHistName   to gsHistName  
// - fHistTitle  to gsHistTitle 
// - fHistXTitle to gsHistXTitle
// - fHistYTitle to gsHistYTitle
//
// - fAbsHistName   to gsAbsHistName  
// - fAbsHistTitle  to gsAbsHistTitle 
// - fAbsHistXTitle to gsAbsHistXTitle
// - fAbsHistYTitle to gsAbsHistYTitle
//
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();

  SetNbins(fgChargeHiGainNbins);
  SetFirst(fgChargeHiGainFirst);
  SetLast (fgChargeHiGainLast );

  SetLoGainNbins(fgChargeLoGainNbins);
  SetLoGainFirst(fgChargeLoGainFirst);
  SetLoGainLast (fgChargeLoGainLast );

  SetHistName  (gsHistName  .Data());
  SetHistTitle (gsHistTitle .Data());
  SetHistXTitle(gsHistXTitle.Data());
  SetHistYTitle(gsHistYTitle.Data());

  SetAbsHistName  (gsAbsHistName  .Data());
  SetAbsHistTitle (gsAbsHistTitle .Data());
  SetAbsHistXTitle(gsAbsHistXTitle.Data());
  SetAbsHistYTitle(gsAbsHistYTitle.Data());
}

// --------------------------------------------------------------------------
//
// 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
//
// Calls:
// - InitHiGainArrays()
// - InitLoGainArrays()
// 
// Sets:
// - fSumhiarea   to nareas 
// - fSumloarea   to nareas 
// - fTimehiarea  to nareas 
// - fTimeloarea  to nareas
// - fSumhisector to nsectors 
// - fSumlosector to nsectors 
// - fTimehisector to nsectors 
// - fTimelosector to nsectors
// - fSathiarea   to nareas 
// - fSatloarea   to nareas
// - fSathisector to nsectors 
// - fSatlosector to nsectors
//
Bool_t MHCalibrationChargeCam::ReInitHists(MParList *pList)
{

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

  fIntensCam = (MCalibrationIntensityCam*)pList->FindObject(AddSerialNumber("MCalibrationIntensityChargeCam"));
  if (fIntensCam)
    *fLog << inf << "Found MCalibrationIntensityChargeCam ... " << endl;
  else
    {
      fCam = (MCalibrationCam*)pList->FindObject(AddSerialNumber("MCalibrationChargeCam"));
      if (!fCam)
        {
          fCam = (MCalibrationCam*)pList->FindCreateObj(AddSerialNumber("MCalibrationChargeCam"));
          if (!fCam)
            {
              *fLog << err << "Cannot find nor create MCalibrationChargeCam ... abort." << endl;
              return kFALSE;
            }
          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)
    {
      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();

  InitHiGainArrays(npixels,nareas,nsectors);
  InitLoGainArrays(npixels,nareas,nsectors);

  fSumhiarea  .Set(nareas); 
  fSumloarea  .Set(nareas); 
  fTimehiarea .Set(nareas); 
  fTimeloarea .Set(nareas);
  fSumhisector.Set(nsectors); 
  fSumlosector.Set(nsectors); 
  fTimehisector.Set(nsectors); 
  fTimelosector.Set(nsectors);

  fSathiarea  .Set(nareas); 
  fSatloarea  .Set(nareas);
  fSathisector.Set(nsectors); 
  fSatlosector.Set(nsectors);

  return kTRUE;
}

// --------------------------------------------------------------------------
//
// Retrieve:
// - fRunHeader->GetNumSamplesHiGain();
//
// Initializes the High Gain Arrays:
//
// - Expand fHiGainArrays to npixels
// - Expand fAverageHiGainAreas to nareas
// - Expand fAverageHiGainSectors to nsectors
// 
// - For every entry in the expanded arrays: 
//   * Initialize an MHCalibrationPix
//   * Set Binning from  fNbins, fFirst and fLast
//   * Set Binning of Abs Times histogram from  fAbsNbins, fAbsFirst and fAbsLast
//   * Set Histgram names and titles from fHistName and fHistTitle
//   * Set Abs Times Histgram names and titles from fAbsHistName and fAbsHistTitle
//   * Set X-axis and Y-axis titles from fHistXTitle and fHistYTitle
//   * Set X-axis and Y-axis titles of Abs Times Histogram from fAbsHistXTitle and fAbsHistYTitle
//   * Call InitHists
//
//
void MHCalibrationChargeCam::InitHiGainArrays(const Int_t npixels, const Int_t nareas, const Int_t nsectors)
{
  
  TH1F *h;

  const Int_t higainsamples = fRunHeader->GetNumSamplesHiGain();

  if (fHiGainArray->GetEntries()==0)
  {
      fHiGainArray->Expand(npixels);
      for (Int_t i=0; i<npixels; i++)
      {
        (*fHiGainArray)[i] = new MHCalibrationChargePix(Form("%s%s",fHistName.Data(),"HiGainPix"),
                                                  Form("%s%s",fHistTitle.Data()," High Gain Pixel"));

        MHCalibrationChargePix &pix = (MHCalibrationChargePix&)(*this)[i];

        pix.SetNbins(fNbins);
        pix.SetFirst(fFirst);
        pix.SetLast (fLast);

        pix.SetAbsTimeNbins(higainsamples);
        pix.SetAbsTimeFirst(-0.5);
        pix.SetAbsTimeLast(higainsamples-0.5);

        h = pix.GetHGausHist();

        h->SetName (Form("%s%s%s","H",fHistName.Data(),"HiGainPix"));
        h->SetTitle(Form("%s%s",fHistTitle.Data()," High Gain Pixel "));  
        h->SetXTitle(fHistXTitle.Data());
        h->SetYTitle(fHistYTitle.Data());
        
        h = pix.GetHAbsTime();

        h->SetName (Form("%s%s%s","H",fAbsHistName.Data(),"HiGainPix"));
        h->SetTitle(Form("%s%s",fAbsHistTitle.Data()," High Gain Pixel "));
        h->SetXTitle(fAbsHistXTitle.Data());
        h->SetYTitle(fAbsHistYTitle.Data());

        InitHists((*this)[i],(*fBadPixels)[i],i);
      }
  }


  if (fAverageHiGainAreas->GetEntries()==0)
  {
    fAverageHiGainAreas->Expand(nareas);
    
    for (Int_t j=0; j<nareas; j++)
      {
        (*fAverageHiGainAreas)[j] = 
          new MHCalibrationChargePix(Form("%s%s",fHistName.Data(),"HiGainArea"),
                                     Form("%s%s",fHistTitle.Data()," High Gain Area Idx "));
        
        MHCalibrationChargePix &pix = (MHCalibrationChargePix&)GetAverageHiGainArea(j);
        
        pix.SetNbins(fNbins*(Int_t)TMath::Sqrt((Float_t)npixels/nareas));
        pix.SetFirst(fFirst);
        pix.SetLast (fLast);
        
        pix.SetAbsTimeNbins(higainsamples);
        pix.SetAbsTimeFirst(-0.5);
        pix.SetAbsTimeLast(higainsamples-0.5);

        h =  pix.GetHGausHist();
        
        h->SetName (Form("%s%s%s","H",fHistName.Data(),"HiGainArea"));
        h->SetXTitle(fHistXTitle.Data());
        h->SetYTitle(fHistYTitle.Data());
        
        if (fGeom->InheritsFrom("MGeomCamMagic"))
          {
            h->SetTitle(Form("%s%s%s%s",fHistTitle.Data()," averaged on event-by-event basis ",
                             j==0 ? "Inner Pixels " : "Outer Pixels ","High Gain Runs: "));
            pix.InitBins();
            pix.SetEventFrequency(fPulserFrequency);
          }
        else
          {
            h->SetTitle(Form("%s%s",fHistTitle.Data(),
                             " averaged on event-by-event basis High Gain Area Idx "));
            InitHists(pix,fIntensCam ? fIntensCam->GetAverageBadArea(j) : fCam->GetAverageBadArea(j),j);
          }

        h =  pix.GetHAbsTime();
        
        h->SetName (Form("%s%s%s","H",fAbsHistName.Data(),"HiGainArea"));
        h->SetTitle(Form("%s%s",fAbsHistTitle.Data(),
                         " averaged on event-by-event basis High Gain Area Idx "));
        h->SetXTitle(fAbsHistXTitle.Data());
        h->SetYTitle(fAbsHistYTitle.Data());
      }
  }
  
  if (fAverageHiGainSectors->GetEntries()==0)
  {
      fAverageHiGainSectors->Expand(nsectors);

      for (Int_t j=0; j<nsectors; j++)
      {
	  (*fAverageHiGainSectors)[j] = 
            new MHCalibrationChargePix(Form("%s%s",fHistName.Data(),"HiGainSector"),
                                       Form("%s%s",fHistTitle.Data()," High Gain Sector "));

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

          pix.SetNbins(fNbins*(Int_t)TMath::Sqrt((Float_t)npixels/nareas));
          pix.SetFirst(fFirst);
          pix.SetLast (fLast);

          pix.SetAbsTimeNbins(higainsamples);
          pix.SetAbsTimeFirst(-0.5);
          pix.SetAbsTimeLast(higainsamples-0.5);
          
          h =  pix.GetHGausHist();
          
          h->SetName (Form("%s%s%s","H",fHistName.Data(),"HiGainSector"));
          h->SetTitle(Form("%s%s",fHistTitle.Data()," High Gain Sector "));  
          h->SetXTitle(fHistXTitle.Data());
          h->SetYTitle(fHistYTitle.Data());

          h =  pix.GetHAbsTime();
          
          h->SetName (Form("%s%s%s","H",fAbsHistName.Data(),"HiGainSector"));
          h->SetTitle(Form("%s%s",fAbsHistTitle.Data(),
                           " averaged on event-by-event basis High Gain Area Sector "));
          h->SetXTitle(fAbsHistXTitle.Data());
          h->SetYTitle(fAbsHistYTitle.Data());
          
          InitHists(pix,fIntensCam ? fIntensCam->GetAverageBadSector(j) : fCam->GetAverageBadSector(j),j);
      }
  }
}

//--------------------------------------------------------------------------------------
//
// Return, if IsLoGain() is kFALSE 
//
// Retrieve:
// - fRunHeader->GetNumSamplesHiGain();
//
// Initializes the Low Gain Arrays:
//
// - Expand fLoGainArrays to npixels
// - Expand fAverageLoGainAreas to nareas
// - Expand fAverageLoGainSectors to nsectors
// 
// - For every entry in the expanded arrays: 
//   * Initialize an MHCalibrationPix
//   * Set Binning from  fNbins, fFirst and fLast
//   * Set Binning of Abs Times histogram from  fAbsNbins, fAbsFirst and fAbsLast
//   * Set Histgram names and titles from fHistName and fHistTitle
//   * Set Abs Times Histgram names and titles from fAbsHistName and fAbsHistTitle
//   * Set X-axis and Y-axis titles from fHistXTitle and fHistYTitle
//   * Set X-axis and Y-axis titles of Abs Times Histogram from fAbsHistXTitle and fAbsHistYTitle
//   * Call InitHists
//
void MHCalibrationChargeCam::InitLoGainArrays(const Int_t npixels, const Int_t nareas, const Int_t nsectors)
{

  if (!IsLoGain())
    return;

  const Int_t logainsamples = fRunHeader->GetNumSamplesLoGain();

  TH1F *h;

  if (fLoGainArray->GetEntries()==0 )
    {
      fLoGainArray->Expand(npixels);
      
      for (Int_t i=0; i<npixels; i++)
        {
          (*fLoGainArray)[i] = 
            new MHCalibrationChargePix(Form("%s%s",fHistName.Data(),"LoGainPix"),
                                       Form("%s%s",fHistTitle.Data()," Low Gain Pixel"));

          MHCalibrationChargePix &pix = (MHCalibrationChargePix&)(*this)(i);

          pix.SetNbins(fLoGainNbins);
          pix.SetFirst(fLoGainFirst);
          pix.SetLast (fLoGainLast);
          
          pix.SetAbsTimeNbins(logainsamples);
          pix.SetAbsTimeFirst(-0.5);
          pix.SetAbsTimeLast(logainsamples-0.5);
          
          h =  pix.GetHGausHist();
          
          h->SetName (Form("%s%s%s","H",fHistName.Data(),"LoGainPix"));
          h->SetTitle(Form("%s%s",fHistTitle.Data()," Low Gain Pixel "));  
          h->SetXTitle(fHistXTitle.Data());
          h->SetYTitle(fHistYTitle.Data());
          
          h = pix.GetHAbsTime();
          
          h->SetName (Form("%s%s%s","H",fAbsHistName.Data(),"HiGainPix"));
          h->SetTitle(Form("%s%s",fAbsHistTitle.Data()," High Gain Pixel "));
          h->SetXTitle(fAbsHistXTitle.Data());
          h->SetYTitle(fAbsHistYTitle.Data());

          InitHists(pix,(*fBadPixels)[i],i);
      }
  }

  if (fAverageLoGainAreas->GetEntries()==0)
    {
      fAverageLoGainAreas->Expand(nareas);
      
      for (Int_t j=0; j<nareas; j++)
        {
	  (*fAverageLoGainAreas)[j] = 
            new MHCalibrationChargePix(Form("%s%s",fHistName.Data(),"LoGainArea"),
                                       Form("%s%s",fHistTitle.Data()," Low Gain Area Idx "));
        
        MHCalibrationChargePix &pix = (MHCalibrationChargePix&)GetAverageLoGainArea(j);

        pix.SetNbins(fLoGainNbins*(Int_t)TMath::Sqrt((Float_t)npixels/nareas));
        pix.SetFirst(fLoGainFirst);
        pix.SetLast (fLoGainLast);

        pix.SetAbsTimeNbins(logainsamples);
        pix.SetAbsTimeFirst(-0.5);
        pix.SetAbsTimeLast(logainsamples-0.5);

        h =  pix.GetHGausHist();
        
        h->SetName (Form("%s%s%s","H",fHistName.Data(),"LoGainArea"));
        h->SetXTitle(fHistXTitle.Data());
        h->SetYTitle(fHistYTitle.Data());


        if (fGeom->InheritsFrom("MGeomCamMagic"))
          {
            h->SetTitle(Form("%s%s%s%s",fHistTitle.Data()," averaged on event-by-event basis ",
                             j==0 ? "Inner Pixels " : "Outer Pixels ","Low Gain Runs: "));
            pix.InitBins();
            pix.SetEventFrequency(fPulserFrequency);
          }
        else
          {
            h->SetTitle(Form("%s%s",fHistTitle.Data()," averaged on event-by-event basis Low Gain Area Idx ")); 
            InitHists(pix,fIntensCam ? fIntensCam->GetAverageBadArea(j) : fCam->GetAverageBadArea(j),j);
          }

        h =  pix.GetHAbsTime();
        
        h->SetName (Form("%s%s%s","H",fAbsHistName.Data(),"LoGainArea"));
        h->SetTitle(Form("%s%s",fAbsHistTitle.Data(),
                         " averaged on event-by-event basis Low Gain Area Idx "));
        h->SetXTitle(fAbsHistXTitle.Data());
        h->SetYTitle(fAbsHistYTitle.Data());

        }
    }
  

  if (fAverageLoGainSectors->GetEntries()==0 && IsLoGain())
  {
      fAverageLoGainSectors->Expand(nsectors);
  
      for (Int_t j=0; j<nsectors; j++)
      {
	  (*fAverageLoGainSectors)[j] = 
            new MHCalibrationChargePix(Form("%s%s",fHistName.Data(),"LoGainSector"),
                                       Form("%s%s",fHistTitle.Data()," Low Gain Sector "));

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

          pix.SetNbins(fLoGainNbins*(Int_t)TMath::Sqrt((Float_t)npixels/nareas));
          pix.SetFirst(fLoGainFirst);
          pix.SetLast (fLoGainLast);

          pix.SetAbsTimeNbins(logainsamples);
          pix.SetAbsTimeFirst(-0.5);
          pix.SetAbsTimeLast(logainsamples-0.5);
          
          h =  pix.GetHGausHist();
          
          h->SetName (Form("%s%s%s","H",fHistName.Data(),"LoGainSector"));
          h->SetTitle(Form("%s%s",fHistTitle.Data()," Low Gain Sector "));  
          h->SetXTitle(fHistXTitle.Data());
          h->SetYTitle(fHistYTitle.Data());

          h =  pix.GetHAbsTime();
          
          h->SetName (Form("%s%s%s","H",fAbsHistName.Data(),"LoGainSector"));
          h->SetTitle(Form("%s%s",fAbsHistTitle.Data(),
                           " averaged on event-by-event basis Low Gain Area Sector "));
          h->SetXTitle(fAbsHistXTitle.Data());
          h->SetYTitle(fAbsHistYTitle.Data());
          
          //
          // Adapt the range for the case, the intense blue is used:
          // FIXME: this is a nasty workaround, but for the moment necessary 
          // in order to avoid default memory space.
          //
          if (fGeom->InheritsFrom("MGeomCamMagic"))
            {
              if ( fColor == MCalibrationCam::kBLUE)
                {
                  pix.SetFirst(-10.5);
                  pix.SetLast(999.5);
                  pix.SetNbins(3030);
                }
            }
          
          InitHists(pix,fIntensCam ? fIntensCam->GetAverageBadSector(j) : fCam->GetAverageBadSector(j),j);
      }
  }
}

  
// --------------------------------------------------------------------------
//
// 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::AddSaturated()) 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();

  fSumhiarea  .Reset(); 
  fSumloarea  .Reset(); 
  fTimehiarea .Reset(); 
  fTimeloarea .Reset();
  fSumhisector.Reset(); 
  fSumlosector.Reset(); 
  fTimehisector.Reset(); 
  fTimelosector.Reset();

  fSathiarea  .Reset(); 
  fSatloarea  .Reset();
  fSathisector.Reset(); 
  fSatlosector.Reset();

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

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

      if (histhi.IsExcluded())
	continue;

      const MExtractedSignalPix &pix = (*signal)[i];
      
      const Float_t sumhi = pix.GetExtractedSignalHiGain();
      const Int_t   sathi = (Int_t)pix.GetNumHiGainSaturated();

      if (IsOscillations())
        histhi.FillHistAndArray(sumhi);
      else
        histhi.FillHist(sumhi);
        
      histhi.AddSaturated(sathi); 

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

      fSumhiarea[aidx]  += sumhi;
      fSathiarea[aidx]  += sathi;

      fSumhisector[sector]  += sumhi;
      fSathisector[sector]  += sathi;

      if (IsLoGain())
        {
          MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)(*this)(i);
          const Float_t sumlo  = pix.GetExtractedSignalLoGain();
          const Int_t   satlo = (Int_t)pix.GetNumLoGainSaturated();
          
          if (IsOscillations())
            histlo.FillHistAndArray(sumlo);
          else
            histlo.FillHist(sumlo);
          
          histlo.AddSaturated(satlo); 
          
          fSumloarea[aidx]  += sumlo;
          fSatloarea[aidx]  += satlo;
          fSumlosector[sector]  += sumlo;
          fSatlosector[sector]  += satlo;
        }
      
    }

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

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

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

      histhi.FillAbsTime(timehi);

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

      fTimehiarea  [aidx]   += timehi;
      fTimehisector[sector] += timehi;

      if (IsLoGain())
        {
          MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)(*this)(pixid);

          const Float_t timelo = (Float_t)pixel.GetIdxMaxLoGainSample(lofirst);
          histlo.FillAbsTime(timelo);

          fTimeloarea[aidx] += timelo;
          fTimelosector[sector] += timelo;
        }
    }
  
  for (UInt_t j=0; j<nareas; j++)
    {

      const Int_t npix = fAverageAreaNum[j];

      if (npix == 0)
        continue;

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

      if (IsOscillations())
        hipix.FillHistAndArray(fSumhiarea [j]/npix);
      else
        hipix.FillHist(fSumhiarea[j]/npix);
      
      hipix.AddSaturated    ((Float_t)fSathiarea [j]/npix > 0.5 ? 1 : 0); 
      hipix.FillAbsTime     (fTimehiarea[j]/npix);

      if (IsLoGain())
        {
          MHCalibrationChargePix &lopix = (MHCalibrationChargePix&)GetAverageLoGainArea(j);
          if (IsOscillations())
            lopix.FillHistAndArray(fSumloarea [j]/npix);
          else
            lopix.FillHist(fSumloarea [j]/npix);
          lopix.AddSaturated    ((Float_t)fSatloarea [j]/npix > 0.5 ? 1 : 0); 
          lopix.FillAbsTime     (fTimeloarea[j]/npix);
        }
    }

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

      const Int_t npix = fAverageSectorNum[j];

      if (npix == 0)
        continue;

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

      if (IsOscillations())      
        hipix.FillHistAndArray(fSumhisector [j]/npix);
      else
        hipix.FillHist(fSumhisector [j]/npix);

      hipix.AddSaturated    ((Float_t)fSathisector[j]/npix > 0.5 ? 1 : 0); 
      hipix.FillAbsTime     (fTimehisector[j]/npix);

      if (IsLoGain())
        {
          MHCalibrationChargePix &lopix = (MHCalibrationChargePix&)GetAverageLoGainSector(j);

          if (IsOscillations())      
            lopix.FillHistAndArray(fSumlosector [j]/npix);
          else
            lopix.FillHist(fSumlosector [j]/npix);
            
          lopix.AddSaturated    ((Float_t)fSatlosector[j]/npix > 0.5 ? 1 : 0); 
          lopix.FillAbsTime     (fTimelosector[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()
{

  *fLog << endl;

  for (Int_t i=0; i<fHiGainArray->GetSize(); i++)
    {
      
      MHCalibrationChargePix &histhi = (MHCalibrationChargePix&)(*this)[i];
      
      if (histhi.IsExcluded())
        continue;
      
      MCalibrationChargePix  &pix    = fIntensCam 
        ? (MCalibrationChargePix&)(*fIntensCam)[i] 
        : (MCalibrationChargePix&)(*fCam)[i];

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

      MBadPixelsPix          &bad    = (*fBadPixels)[i];

      Stat_t overflow = histhi.GetHGausHist()->GetBinContent(histhi.GetHGausHist()->GetNbinsX()+1);
      if (overflow > 0.1)
        {
          *fLog << warn << GetDescriptor()
                << ": HiGain Histogram Overflow occurred " << overflow 
                << " times in pixel: " << i << " (without saturation!) " << endl;
          bad.SetUncalibrated( MBadPixelsPix::kHiGainOverFlow ); 
        }

      overflow = histhi.GetHGausHist()->GetBinContent(0);
      if (overflow > 0.1)
        {
          *fLog << warn << GetDescriptor()
                << ": HiGain Histogram Underflow occurred " << overflow
                << " times in pixel: " << i << " (without saturation!) " << endl;
          bad.SetUncalibrated( MBadPixelsPix::kHiGainOverFlow ); 
        }
      
      FinalizeAbsTimes(histhi, pix, bad, fFirstHiGain, fLastHiGain);
    }

  if (IsLoGain())
    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 ); 
            if (IsOscillations())
              histlo.CreateFourierSpectrum();
            continue;
          }
        
        Stat_t overflow = histlo.GetHGausHist()->GetBinContent(histlo.GetHGausHist()->GetNbinsX()+1);
        if (overflow > 0.1)
          {
            *fLog << warn << GetDescriptor()
                  << ": LoGain Histogram Overflow occurred " << overflow
                  << " times in pixel: " << i << " (without saturation!) " << endl;
            bad.SetUncalibrated( MBadPixelsPix::kLoGainOverFlow ); 
          }

        overflow = histlo.GetHGausHist()->GetBinContent(0);
        if (overflow > 0.1)
          {
            *fLog << warn << GetDescriptor()
                  << ": LoGain Histogram Underflow occurred " << overflow
                  << " times in pixel: " << i << " (without saturation!) " << endl;
            bad.SetUncalibrated( MBadPixelsPix::kLoGainOverFlow ); 
          }
        
        MCalibrationChargePix &pix = fIntensCam 
          ? (MCalibrationChargePix&)(*fIntensCam)[i]
          : (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    = fIntensCam 
        ? (MCalibrationChargePix&)fIntensCam->GetAverageArea(j)
        : (MCalibrationChargePix&)fCam->GetAverageArea(j);

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

      MBadPixelsPix          &bad    = fIntensCam 
        ? fIntensCam->GetAverageBadArea(j)
        : fCam->GetAverageBadArea(j);

      FinalizeAbsTimes(histhi, pix, bad, fFirstHiGain, fLastHiGain);
   }
  
  if (IsLoGain())
    for (Int_t j=0; j<fAverageLoGainAreas->GetSize(); j++)
      {
        
        MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)GetAverageLoGainArea(j);      

        if (histlo.GetSaturated() > fNumLoGainSaturationLimit*histlo.GetHGausHist()->GetEntries())
          {
            *fLog << warn << "Saturated Lo Gain histogram in area idx: " << j << endl;
            histlo.CreateFourierSpectrum();
            continue;
          }
        
        MCalibrationChargePix  &pix    = fIntensCam
          ? (MCalibrationChargePix&)fIntensCam->GetAverageArea(j) 
          : (MCalibrationChargePix&)fCam->GetAverageArea(j) ;

        if (pix.IsHiGainSaturation())
          {
            MBadPixelsPix          &bad    = fIntensCam
              ? fIntensCam->GetAverageBadArea(j)  
              : fCam->GetAverageBadArea(j);      
            FinalizeAbsTimes(histlo, pix, bad, fFirstLoGain, fLastLoGain);
          }
        
      }
  
  for (Int_t j=0; j<fAverageHiGainSectors->GetSize(); j++)
    {
      
      MHCalibrationChargePix &histhi = (MHCalibrationChargePix&)GetAverageHiGainSector(j);      
      MCalibrationChargePix  &pix    = fIntensCam
        ? (MCalibrationChargePix&)fIntensCam->GetAverageSector(j)  
        : (MCalibrationChargePix&)fCam->GetAverageSector(j);

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

      MBadPixelsPix          &bad    = fIntensCam 
        ? fIntensCam->GetAverageBadSector(j)
        : fCam->GetAverageBadSector(j);      

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

        if (pix.IsHiGainSaturation())
          FinalizeAbsTimes(histlo, pix, bad, fFirstLoGain, fLastLoGain);
      }
  
  //
  // Perform the fitting for the High Gain (done in MHCalibrationCam)
  //
  FitHiGainArrays(fIntensCam ? (MCalibrationCam&)(*fIntensCam->GetCam()) : (MCalibrationCam&)(*fCam),
                  *fBadPixels,
                  MBadPixelsPix::kHiGainNotFitted,
                  MBadPixelsPix::kHiGainOscillating);
  
  //
  // Perform the fitting for the Low Gain (done in MHCalibrationCam)
  //
  if (IsLoGain())
    FitLoGainArrays(fIntensCam ? (MCalibrationCam&)(*fIntensCam->GetCam()) : (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    = fIntensCam ? (*fIntensCam)[i] : (*fCam)[i];

      if (bad.IsUncalibrated( MBadPixelsPix::kHiGainNotFitted ))
        if (!pix.IsHiGainSaturation())
          bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun    );
 
      if (bad.IsUncalibrated( MBadPixelsPix::kLoGainNotFitted ))
        if (pix.IsHiGainSaturation())
          bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun    );
 
      if (bad.IsUncalibrated( MBadPixelsPix::kLoGainSaturation ))
          bad.SetUnsuitable(   MBadPixelsPix::kUnsuitableRun    );

      if (IsOscillations())
        {
          if (bad.IsUncalibrated( MBadPixelsPix::kHiGainOscillating ))
            bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun    );
          
          if (bad.IsUncalibrated( MBadPixelsPix::kLoGainOscillating ))
            if (pix.IsHiGainSaturation())
              bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun    );
        }
    }
}

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

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


// -----------------------------------------------------------------------------
// 
// Default draw:
//
// Displays the averaged areas, both High Gain and Low Gain 
//
// Calls the Draw of the fAverageHiGainAreas and fAverageLoGainAreas objects with options
//
void MHCalibrationChargeCam::Draw(const Option_t *opt)
{

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

  TString option(opt);
  option.ToLower();

  if (!option.Contains("datacheck"))
    {
      MHCalibrationCam::Draw(opt);
      return;
    }

  // 
  // From here on , the datacheck - Draw
  //
  TVirtualPad *pad = gPad ? gPad : MH::MakeDefCanvas(this);  
  pad->SetBorderMode(0);
  pad->Divide(1,nareas);

  //
  // Loop over inner and outer pixels
  //  
  for (Int_t i=0; i<nareas;i++) 
     {
       
       pad->cd(i+1);
       
       MHCalibrationChargePix &hipix = (MHCalibrationChargePix&)GetAverageHiGainArea(i);
       //
       // Ask for Hi-Gain saturation
       //
       if (hipix.GetSaturated() > fNumHiGainSaturationLimit*hipix.GetHGausHist()->GetEntries() && IsLoGain())
         {
           MHCalibrationChargePix &lopix = (MHCalibrationChargePix&)GetAverageLoGainArea(i);
           DrawDataCheckPixel(lopix,i ? gkLoGainOuterRefLines : gkLoGainInnerRefLines);
         }
       else
         DrawDataCheckPixel(hipix,i ? gkHiGainOuterRefLines : gkHiGainInnerRefLines);
    }      
}


// --------------------------------------------------------------------------
//
// Our own clone function is necessary since root 3.01/06 or Mars 0.4
// I don't know the reason. 
//
// Creates new MHCalibrationChargeCam only for the Averaged Areas, 
// the rest has to be retrieved directly, e.g. via: 
//  MHCalibrationChargeCam *cam = MParList::FindObject("MHCalibrationChargeCam");
//  -  cam->GetAverageSector(5).DrawClone();
//  -  (*cam)[100].DrawClone()
//
TObject *MHCalibrationChargeCam::Clone(const char *name) const
{

  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->fAverageHiGainAreas->Expand(navhi);
  //  cam->fAverageHiGainSectors->Expand(nsehi);

  for (int i=0; i<navhi; i++)
    (*cam->fAverageHiGainAreas)  [i] = (*fAverageHiGainAreas)  [i]->Clone();
  //  for (int i=0; i<nsehi; i++)
  //    (*cam->fAverageHiGainSectors)[i] = (*fAverageHiGainSectors)[i]->Clone();

  if (IsLoGain())
    {
      cam->fAverageLoGainAreas->Expand(navlo);
      //      cam->fAverageLoGainSectors->Expand(nselo);
      
      for (int i=0; i<navlo; i++)
        (*cam->fAverageLoGainAreas)  [i] = (*fAverageLoGainAreas)  [i]->Clone();
      //      for (int i=0; i<nselo; i++)
      //        (*cam->fAverageLoGainSectors)[i] = (*fAverageLoGainSectors)[i]->Clone();
    }
  
  cam->fAverageAreaNum         = fAverageAreaNum;
  cam->fAverageAreaSat         = fAverageAreaSat;
  cam->fAverageAreaSigma       = fAverageAreaSigma;      
  cam->fAverageAreaSigmaVar    = fAverageAreaSigmaVar;   
  cam->fAverageAreaRelSigma    = fAverageAreaRelSigma;
  cam->fAverageAreaRelSigmaVar = fAverageAreaRelSigmaVar;   
  cam->fAverageSectorNum       = fAverageSectorNum;      
  cam->fRunNumbers             = fRunNumbers;

  cam->fColor                  = fColor;
  cam->fPulserFrequency        = fPulserFrequency;
  cam->fFlags                  = fFlags;

  return cam;

}

// -----------------------------------------------------------------------------
// 
// Draw the average pixel for the datacheck:
//
// Displays the averaged areas, both High Gain and Low Gain 
//
// Calls the Draw of the fAverageHiGainAreas and fAverageLoGainAreas objects with options
//
void MHCalibrationChargeCam::DrawDataCheckPixel(MHCalibrationChargePix &pix, const Float_t refline[])
{
  
  TVirtualPad *newpad = gPad;
  newpad->Divide(1,2);
  newpad->cd(1);
  
  gPad->SetTicks();
  if (!pix.IsEmpty() && !pix.IsOnlyOverflow() && !pix.IsOnlyUnderflow())
    gPad->SetLogy();

  TH1F *hist = pix.GetHGausHist();

  TH1F *null = new TH1F("Null",hist->GetTitle(),100,
                        pix.GetFirst() > 0. ? pix.GetFirst() : 0.,
                        pix.GetLast() > pix.GetFirst()
                        ? ( pix.GetLast() > 450. ? 450. : pix.GetLast() )
                        : pix.GetFirst()*2.);

  null->SetMaximum(1.1*hist->GetMaximum());
  null->SetDirectory(NULL);
  null->SetBit(kCanDelete);
  null->SetStats(kFALSE);
  //
  // set the labels bigger
  //
  TAxis *xaxe = null->GetXaxis();
  TAxis *yaxe = null->GetYaxis();
  xaxe->CenterTitle();
  yaxe->CenterTitle();    
  xaxe->SetTitleSize(0.07);
  yaxe->SetTitleSize(0.07);    
  xaxe->SetTitleOffset(0.7);
  yaxe->SetTitleOffset(0.55);    
  xaxe->SetLabelSize(0.06);
  yaxe->SetLabelSize(0.06);    
  xaxe->SetTitle(hist->GetXaxis()->GetTitle());
  yaxe->SetTitle(hist->GetYaxis()->GetTitle());  

  null->Draw();
  hist->Draw("same");

  gStyle->SetOptFit();

  TF1 *fit = pix.GetFGausFit();

  if (fit)
  {
    switch ( fColor )
      {
      case MCalibrationCam::kGREEN:
        fit->SetLineColor(kGreen);
        break;
      case MCalibrationCam::kBLUE:
        fit->SetLineColor(kBlue);
        break;
      case MCalibrationCam::kUV:  
        fit->SetLineColor(106);
        break;
      case MCalibrationCam::kCT1: 
        fit->SetLineColor(006);
        break;
      default:
        fit->SetLineColor(kRed);
      }
    fit->Draw("same");
  }

  DisplayRefLines(null,refline);

  newpad->cd(2);
  gPad->SetTicks();

  TH1F *null2 = new TH1F("Null2",hist->GetTitle(),100,0.,pix.GetEvents()->GetSize()/pix.GetEventFrequency());

  null2->SetMinimum(pix.GetMean()-10.*pix.GetSigma());
  null2->SetMaximum(pix.GetMean()+10.*pix.GetSigma());
  null2->SetDirectory(NULL);
  null2->SetBit(kCanDelete);
  null2->SetStats(kFALSE);
  //
  // set the labels bigger
  //
  TAxis *xaxe2 = null2->GetXaxis();
  TAxis *yaxe2 = null2->GetYaxis();
  xaxe2->CenterTitle();
  yaxe2->CenterTitle();    
  xaxe2->SetTitleSize(0.07);
  yaxe2->SetTitleSize(0.07);    
  xaxe2->SetTitleOffset(0.7);
  yaxe2->SetTitleOffset(0.55);    
  xaxe2->SetLabelSize(0.06);
  yaxe2->SetLabelSize(0.06);    

  pix.CreateGraphEvents();
  TGraph *gr = pix.GetGraphEvents();

  xaxe2->SetTitle(gr->GetXaxis()->GetTitle());
  yaxe2->SetTitle(gr->GetYaxis()->GetTitle());  

  null2->Draw();

  pix.DrawEvents("same");
  return;
  
}


void  MHCalibrationChargeCam::DisplayRefLines(const TH1F *hist, const Float_t refline[]) const
{

  TLine *green1 = new TLine(refline[0],0.,refline[0],0.9*hist->GetMaximum());
  green1->SetBit(kCanDelete);
  green1->SetLineColor(kGreen);
  green1->SetLineStyle(2);
  green1->SetLineWidth(3);
  green1->Draw();

  TLine *green5 = new TLine(refline[6],0.,refline[6],0.9*hist->GetMaximum());
  green5->SetBit(kCanDelete);
  green5->SetLineColor(8);
  green5->SetLineStyle(2);
  green5->SetLineWidth(3);
  green5->Draw();

  TLine *blue1   = new TLine(refline[1],0.,refline[1],0.9*hist->GetMaximum());
  blue1->SetBit(kCanDelete);
  blue1->SetLineColor(227);
  blue1->SetLineStyle(2);
  blue1->SetLineWidth(3);
  blue1->Draw();

  TLine *blue5   = new TLine(refline[2],0.,refline[2],0.9*hist->GetMaximum());
  blue5->SetBit(kCanDelete);
  blue5->SetLineColor(68);
  blue5->SetLineStyle(2);
  blue5->SetLineWidth(3);
  blue5->Draw();

  TLine *blue10   = new TLine(refline[3],0.,refline[3],0.9*hist->GetMaximum());
  blue10->SetBit(kCanDelete);
  blue10->SetLineColor(4);
  blue10->SetLineStyle(2);
  blue10->SetLineWidth(3);
  blue10->Draw();

  TLine *uv10    = new TLine(refline[4],0.,refline[4],0.9*hist->GetMaximum());
  uv10->SetBit(kCanDelete);
  uv10->SetLineColor(106);
  uv10->SetLineStyle(2);
  uv10->SetLineWidth(3);
  uv10->Draw();

  TLine *ct1    = new TLine(refline[5],0.,refline[5],0.9*hist->GetMaximum());
  ct1->SetBit(kCanDelete);
  ct1->SetLineColor(6);
  ct1->SetLineStyle(2);
  ct1->SetLineWidth(3);
  ct1->Draw();

  TLegend *leg = new TLegend(0.8,0.35,0.99,0.99);
  leg->SetBit(kCanDelete);
  leg->AddEntry(green1,"1 Led GREEN","l");
  leg->AddEntry(green5,"5 Leds GREEN","l");
  leg->AddEntry(blue1,"1 Led BLUE","l");
  leg->AddEntry(blue5,"5 Leds BLUE","l");
  leg->AddEntry(blue10,"10 Leds BLUE","l");
  leg->AddEntry(uv10,"10 Leds UV","l");
  leg->AddEntry(ct1,"CT1-Pulser","l");

  leg->Draw();
}