source: trunk/MagicSoft/Mars/mcalib/MCalibrationChargeCam.cc@ 3529

/* ======================================================================== *\
!
! *
! * 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   11/2003 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                               
// MCalibrationChargeCam                                               
//                                                               
// Hold the whole Calibration results of the camera:
//                                                               
// 1) MCalibrationChargeCam initializes a TClonesArray whose elements are 
//    pointers to MCalibrationChargePix Containers
// 2) It initializes a pointer to an MCalibrationBlindPix container
// 3) It initializes a pointer to an MCalibrationPINDiode container
//
//
// The calculated values (types of GetPixelContent) are:
// 
// --------------------------------------------------------------------------
//
// The types are as follows:
// 
// Fitted values:
// ============== 
//
// 0: Fitted Charge
// 1: Error of fitted Charge
// 2: Sigma of fitted Charge
// 3: Error of Sigma of fitted Charge
//
// Useful variables derived from the fit results:
// =============================================
//
// 4: Returned probability of Gauss fit to Charge distribution
// 5: Reduced Sigma of fitted Charge --> sqrt(sigma_Q^2 - PedRMS^2)
// 6: Error Reduced Sigma of fitted Charge 
// 7: Reduced Sigma per Charge 
// 8: Error of Reduced Sigma per Charge 
//
// Results of the different calibration methods:
// =============================================
//
// 9: Number of Photo-electrons obtained with the F-Factor method
// 10: Error on Number of Photo-electrons obtained with the F-Factor method
// 11: Mean conversion factor obtained with the F-Factor method
// 12: Error on the mean conversion factor obtained with the F-Factor method
// 13: Overall F-Factor of the readout obtained with the F-Factor method
// 14: Error on Overall F-Factor of the readout obtained with the F-Factor method
// 15: Number of Photons inside Plexiglass obtained with the Blind Pixel method
// 16: Error on Number of Photons inside Plexiglass obtained with the Blind Pixel method
// 17: Mean conversion factor obtained with the Blind Pixel method
// 18: Error on the mean conversion factor obtained with the Blind Pixel method
// 19: Overall F-Factor of the readout obtained with the Blind Pixel method
// 20: Error on Overall F-Factor of the readout obtained with the Blind Pixel method
// 21: Number of Photons outside Plexiglass obtained with the PIN Diode method
// 22: Error on Number of Photons outside Plexiglass obtained with the PIN Diode method
// 23: Mean conversion factor obtained with the PIN Diode method
// 24: Error on the mean conversion factor obtained with the PIN Diode method
// 25: Overall F-Factor of the readout obtained with the PIN Diode method
// 26: Error on Overall F-Factor of the readout obtained with the PIN Diode method
//
// Localized defects:
// ==================
//
// 27: Excluded Pixels
// 28: Pixels where the fit did not succeed --> results obtained only from the histograms
// 29: Number of probable pickup events in the Hi Gain 
// 30: Number of probable pickup events in the Lo Gain
//
// Other classifications of pixels:
// ================================
//
// 31: Pixels with saturated Hi-Gain
//
// Classification of validity of the calibrations:
// ===============================================
//
// 32: Pixels with valid calibration by the F-Factor-Method
// 33: Pixels with valid calibration by the Blind Pixel-Method
// 34: Pixels with valid calibration by the PIN Diode-Method
//
// Used Pedestals:
// ===============
//
// 35: Mean Pedestal over the entire range of signal extraction
// 36: Error on the Mean Pedestal over the entire range of signal extraction
// 37: Pedestal RMS over the entire range of signal extraction
// 38: Error on the Pedestal RMS over the entire range of signal extraction
//
// Calculated absolute arrival times (very low precision!):
// ========================================================
//
// 39: Absolute Arrival time of the signal
// 40: RMS of the Absolute Arrival time of the signal
//
/////////////////////////////////////////////////////////////////////////////
#include "MCalibrationChargeCam.h"

#include <TH2.h>
#include <TCanvas.h>
#include <TClonesArray.h>

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

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

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

#include "MCalibrationChargePix.h"
#include "MCalibrationChargeBlindPix.h"
#include "MCalibrationChargePINDiode.h"

ClassImp(MCalibrationChargeCam);

using namespace std;

const Float_t MCalibrationChargeCam::gkAverageQE                = 0.18;     
const Float_t MCalibrationChargeCam::gkAverageQEErr             = 0.02;  
const Float_t MCalibrationChargeCam::fgConvFFactorRelErrLimit   = 0.35;
const Float_t MCalibrationChargeCam::fgPheFFactorRelErrLimit    = 5.;
// --------------------------------------------------------------------------
//
// Default constructor. 
//
// Creates a TClonesArray of MCalibrationPix containers, initialized to 1 entry
// Later, a call to MCalibrationChargeCam::InitSize(Int_t size) has to be performed
//
// Creates an MCalibrationBlindPix container 
//
MCalibrationChargeCam::MCalibrationChargeCam(const char *name, const char *title)
    : fOffsets(NULL),
      fSlopes(NULL),
      fOffvsSlope(NULL)
{
    fName  = name  ? name  : "MCalibrationChargeCam";
    fTitle = title ? title : "Storage container for the Calibration Information in the camera";

    fPixels          = new TClonesArray("MCalibrationChargePix",1);
    fAverageInnerPix = new MCalibrationChargePix("AverageInnerPix","Container of the fit results of the camera average inner pixels");
    fAverageOuterPix = new MCalibrationChargePix("AverageOuterPix","Container of the fit results of the camera average outer pixels");
    fAverageInnerBadPix = new MBadPixelsPix("AverageInnerBadPix","Bad Pixel Container of the camera average inner pixels");
    fAverageOuterBadPix = new MBadPixelsPix("AverageOuterBadPix","Bad Pixel Container of the camera average outer pixels");

    Clear();

    SetAverageQE();
    SetConvFFactorRelErrLimit();
    SetPheFFactorRelErrLimit();
}

// --------------------------------------------------------------------------
//
// Delete the TClonesArray of MCalibrationPix containers
// Delete the MCalibrationPINDiode and the MCalibrationBlindPix
//
// Delete the histograms if they exist
//
MCalibrationChargeCam::~MCalibrationChargeCam()
{

  //
  // delete fPixels should delete all Objects stored inside
  // 
  delete fPixels;
  delete fAverageInnerPix;
  delete fAverageOuterPix;
  
  delete fAverageInnerBadPix;
  delete fAverageOuterBadPix;
  
  if (fOffsets)
    delete fOffsets;
  if (fSlopes)
    delete fSlopes;
  if (fOffvsSlope)
    delete fOffvsSlope;

}

// -------------------------------------------------------------------
//
// This function simply allocates memory via the ROOT command:
// (TObject**) TStorage::ReAlloc(fCont, newSize * sizeof(TObject*),
//                                      fSize * sizeof(TObject*));
// newSize corresponds to size in our case
// fSize is the old size (in most cases: 1)
//
void MCalibrationChargeCam::InitSize(const UInt_t i)
{
  
  fPixels->ExpandCreate(i);

}

// --------------------------------------------------------------------------
//
// This function returns the current size of the TClonesArray 
// independently if the MCalibrationPix is filled with values or not.
//
// It is the size of the array fPixels.
//
Int_t MCalibrationChargeCam::GetSize() const
{
  return fPixels->GetEntriesFast();
}


// --------------------------------------------------------------------------
//
// Get i-th pixel (pixel number)
//
MCalibrationChargePix &MCalibrationChargeCam::operator[](UInt_t i)
{
  return *static_cast<MCalibrationChargePix*>(fPixels->UncheckedAt(i));
}

// --------------------------------------------------------------------------
//
// Get i-th pixel (pixel number)
//
const MCalibrationChargePix &MCalibrationChargeCam::operator[](UInt_t i) const
{
  return *static_cast<MCalibrationChargePix*>(fPixels->UncheckedAt(i));
}


// --------------------------------------
//
void MCalibrationChargeCam::Clear(Option_t *o)
{

  fPixels->ForEach(TObject, Clear)();
  fAverageInnerPix->Clear();
  fAverageOuterPix->Clear();

  fAverageInnerBadPix->Clear();
  fAverageOuterBadPix->Clear();

  fNumExcludedPixels            = 0;
  fMeanFluxPhesInnerPixel       = 0.;
  fMeanFluxPhesInnerPixelVar    = 0.;
  fMeanFluxPhesOuterPixel       = 0.;
  fMeanFluxPhesOuterPixelVar    = 0.;

  CLRBIT(fFlags,kBlindPixelMethodValid);
  CLRBIT(fFlags,kFFactorMethodValid);
  CLRBIT(fFlags,kPINDiodeMethodValid);

  return;
}

void MCalibrationChargeCam::SetFFactorMethodValid(const Bool_t b)
{
    b ? SETBIT(fFlags, kFFactorMethodValid) : CLRBIT(fFlags, kFFactorMethodValid); 
}

void MCalibrationChargeCam::SetBlindPixelMethodValid(const Bool_t b)
{
    b ? SETBIT(fFlags, kBlindPixelMethodValid) : CLRBIT(fFlags, kBlindPixelMethodValid); 
}

void MCalibrationChargeCam::SetPINDiodeMethodValid(const Bool_t b)
{
    b ? SETBIT(fFlags, kPINDiodeMethodValid) : CLRBIT(fFlags, kPINDiodeMethodValid); 
}

Float_t MCalibrationChargeCam::GetMeanFluxPhesInnerPixelErr()   const
{
  if (fMeanFluxPhesInnerPixelVar <= 0.)
    return -1.;
  return TMath::Sqrt(fMeanFluxPhesInnerPixelVar);
}

Float_t MCalibrationChargeCam::GetMeanFluxPhesOuterPixelErr()   const
{
  if (fMeanFluxPhesOuterPixelVar <= 0.)
    return -1.;
  return TMath::Sqrt(fMeanFluxPhesOuterPixelVar);
}

Float_t MCalibrationChargeCam::GetMeanFluxPhotonsInnerPixelErr()   const
{
  if (fMeanFluxPhotonsInnerPixelVar <= 0.)
    return -1.;
  return TMath::Sqrt(fMeanFluxPhotonsInnerPixelVar);
}

Float_t MCalibrationChargeCam::GetMeanFluxPhotonsOuterPixelErr()   const
{
  if (fMeanFluxPhotonsOuterPixelVar <= 0.)
    return -1.;
  return TMath::Sqrt(fMeanFluxPhotonsOuterPixelVar);
}

Bool_t  MCalibrationChargeCam::IsBlindPixelMethodValid()   const
{
  return TESTBIT(fFlags,kBlindPixelMethodValid);
}

Bool_t  MCalibrationChargeCam::IsPINDiodeMethodValid() const
{
  return TESTBIT(fFlags,kPINDiodeMethodValid);  
}


// --------------------------------------------------------------------------
//
// Print first the well fitted pixels 
// and then the ones which are not FitValid
//
void MCalibrationChargeCam::Print(Option_t *o) const
{

  *fLog << all << GetDescriptor() << ":" << endl;
  int id = 0;
  
  *fLog << all << "Calibrated pixels:" << endl;
  *fLog << all << endl;

  TIter Next(fPixels);
  MCalibrationChargePix *pix;
  while ((pix=(MCalibrationChargePix*)Next()))
    {
      
      if (!pix->IsExcluded()) 
        {                            

          *fLog << all << "Pix " << pix->GetPixId() 
                << ":  Ped.  Rms: "            << pix->GetPedRms()        << " +- " << pix->GetPedRmsErr() 
                << "   Mean signal: "          << pix->GetMeanCharge()    << " +- " << pix->GetSigmaCharge() 
                << "   Reduced Sigma: "        << pix->GetRSigmaCharge() 
                << "   Nr Phe's: "             << pix->GetPheFFactorMethod() 
                << " Saturated? :" << pix->IsHiGainSaturation() 
                << endl;
          id++;
        }
    }
  
  *fLog << all << id << " pixels" << endl;
  id = 0;
  
   
  *fLog << all << endl;
  *fLog << all << "Excluded pixels:" << endl;
  *fLog << all << endl;
  
  id = 0;

  TIter Next4(fPixels);
  while ((pix=(MCalibrationChargePix*)Next4()))
  {
      if (pix->IsExcluded())
      {
          *fLog << all << pix->GetPixId() << endl;
          id++;
      }
  }
  *fLog << all << id << " Excluded pixels " << endl;
  *fLog << endl;

  *fLog << all << "Average Inner Pix:" 
                << "   Ped.  Rms: "            << fAverageInnerPix->GetPedRms()        << " +- " << fAverageInnerPix->GetPedRmsErr() 
                << "   Mean signal: "          << fAverageInnerPix->GetMeanCharge()    << " +- " << fAverageInnerPix->GetMeanChargeErr() 
                << "   Sigma signal: "         << fAverageInnerPix->GetSigmaCharge()    << " +- "<< fAverageInnerPix->GetSigmaChargeErr() 
                << "   Reduced Sigma: "        << fAverageInnerPix->GetRSigmaCharge() 
                << "   Nr Phe's: "             << fAverageInnerPix->GetPheFFactorMethod() 
                << endl;
  *fLog << all << "Average Outer Pix:" 
                << "   Ped.  Rms: "            << fAverageOuterPix->GetPedRms()        << " +- " << fAverageOuterPix->GetPedRmsErr() 
                << "   Mean signal: "          << fAverageOuterPix->GetMeanCharge()    << " +- " << fAverageOuterPix->GetMeanChargeErr() 
                << "   Sigma signal: "         << fAverageOuterPix->GetSigmaCharge()    << " +- "<< fAverageOuterPix->GetSigmaChargeErr() 
                << "   Reduced Sigma: "        << fAverageOuterPix->GetRSigmaCharge() 
                << "   Nr Phe's: "             << fAverageOuterPix->GetPheFFactorMethod() 
                << endl;

}


// --------------------------------------------------------------------------
//
// The types are as follows:
// 
// Fitted values:
// ============== 
//
// 0: Fitted Charge
// 1: Error of fitted Charge
// 2: Sigma of fitted Charge
// 3: Error of Sigma of fitted Charge
//
// Useful variables derived from the fit results:
// =============================================
//
// 4: Returned probability of Gauss fit to Charge distribution
// 5: Reduced Sigma of fitted Charge --> sqrt(sigma_Q^2 - PedRMS^2)
// 6: Error Reduced Sigma of fitted Charge 
// 7: Reduced Sigma per Charge 
// 8: Error of Reduced Sigma per Charge 
//
// Results of the different calibration methods:
// =============================================
//
// 9: Number of Photo-electrons obtained with the F-Factor method
// 10: Error on Number of Photo-electrons obtained with the F-Factor method
// 11: Mean conversion factor obtained with the F-Factor method
// 12: Error on the mean conversion factor obtained with the F-Factor method
// 13: Overall F-Factor of the readout obtained with the F-Factor method
// 14: Error on Overall F-Factor of the readout obtained with the F-Factor method
// 15: Number of Photons inside Plexiglass obtained with the Blind Pixel method
// 16: Error on Number of Photons inside Plexiglass obtained with the Blind Pixel method
// 17: Mean conversion factor obtained with the Blind Pixel method
// 18: Error on the mean conversion factor obtained with the Blind Pixel method
// 19: Overall F-Factor of the readout obtained with the Blind Pixel method
// 20: Error on Overall F-Factor of the readout obtained with the Blind Pixel method
// 21: Number of Photons outside Plexiglass obtained with the PIN Diode method
// 22: Error on Number of Photons outside Plexiglass obtained with the PIN Diode method
// 23: Mean conversion factor obtained with the PIN Diode method
// 24: Error on the mean conversion factor obtained with the PIN Diode method
// 25: Overall F-Factor of the readout obtained with the PIN Diode method
// 26: Error on Overall F-Factor of the readout obtained with the PIN Diode method
//
// Localized defects:
// ==================
//
// 27: Excluded Pixels
// 28: Pixels where the fit did not succeed --> results obtained only from the histograms
// 29: Number of probable pickup events in the Hi Gain 
// 30: Number of probable pickup events in the Lo Gain
//
// Other classifications of pixels:
// ================================
//
// 31: Pixels with saturated Hi-Gain
//
// Classification of validity of the calibrations:
// ===============================================
//
// 32: Pixels with valid calibration by the F-Factor-Method
// 33: Pixels with valid calibration by the Blind Pixel-Method
// 34: Pixels with valid calibration by the PIN Diode-Method
//
// Used Pedestals:
// ===============
//
// 35: Mean Pedestal over the entire range of signal extraction
// 36: Error on the Mean Pedestal over the entire range of signal extraction
// 37: Pedestal RMS over the entire range of signal extraction
// 38: Error on the Pedestal RMS over the entire range of signal extraction
//
// Calculated absolute arrival times (very low precision!):
// ========================================================
//
// 39: Absolute Arrival time of the signal
// 40: RMS of the Absolute Arrival time of the signal
//
Bool_t MCalibrationChargeCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
{

  if (idx > GetSize())
    return kFALSE;

  Float_t area = cam[idx].GetA();

 if (area == 0)
    return kFALSE;

  switch (type)
    {
    case 0:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetMeanCharge();
      break;
    case 1:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetMeanChargeErr();
      break;
    case 2:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetSigmaCharge();
      break;
    case 3:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetSigmaChargeErr();
      break;
    case 4:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetChargeProb();
      break;
    case 5:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      if ((*this)[idx].GetRSigmaCharge() == -1.)
          return kFALSE;
      val = (*this)[idx].GetRSigmaCharge();
      break;
    case 6:
      if ((*this)[idx].IsExcluded())
        return kFALSE;    
      if ((*this)[idx].GetRSigmaCharge() == -1.)
          return kFALSE;
      val = (*this)[idx].GetRSigmaChargeErr();
      break;
    case 7:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      if ((*this)[idx].GetRSigmaCharge() == -1.)
        return kFALSE;
      if ((*this)[idx].GetMeanCharge() == 0.)
        return kFALSE;
      val = (*this)[idx].GetRSigmaCharge() / (*this)[idx].GetMeanCharge();
      break;
    case 8:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      if ((*this)[idx].GetRSigmaCharge() <= 0.)
          return kFALSE;
      if ((*this)[idx].GetMeanCharge() <= 0.)
        return kFALSE;
      if ((*this)[idx].GetRSigmaChargeErr() <= 0.)
          return kFALSE;
      if ((*this)[idx].GetMeanChargeErr() <= 0.)
        return kFALSE;
      // relative error RsigmaCharge square
      val =    (*this)[idx].GetRSigmaChargeErr()* (*this)[idx].GetRSigmaChargeErr() 
            / ((*this)[idx].GetRSigmaCharge()   * (*this)[idx].GetRSigmaCharge()   );
      // relative error Charge square
      val +=   (*this)[idx].GetMeanChargeErr() * (*this)[idx].GetMeanChargeErr()
            / ((*this)[idx].GetMeanCharge()    * (*this)[idx].GetMeanCharge()   );
      // calculate relative error out of squares
      val  =   TMath::Sqrt(val) ;
      // multiply with value to get absolute error
      val  *=  (*this)[idx].GetRSigmaCharge() / (*this)[idx].GetMeanCharge();
      break;
    case 9:
      if ((*this)[idx].IsExcluded() || !(*this)[idx].IsFFactorMethodValid())
        return kFALSE;
      val = (*this)[idx].GetPheFFactorMethod();
      break;
    case 10:
      if ((*this)[idx].IsExcluded() || !(*this)[idx].IsFFactorMethodValid())
        return kFALSE;
      val = (*this)[idx].GetPheFFactorMethodErr();
      break;
    case 11:
      if ((*this)[idx].IsExcluded() || !(*this)[idx].IsFFactorMethodValid())
        return kFALSE;
      val = (*this)[idx].GetMeanConversionFFactorMethod();
      break;
    case 12:
      if ((*this)[idx].IsExcluded() || !(*this)[idx].IsFFactorMethodValid())
        return kFALSE;
      val = (*this)[idx].GetConversionFFactorMethodErr();
      break;
    case 13:
      if ((*this)[idx].IsExcluded() || !(*this)[idx].IsFFactorMethodValid())
        return kFALSE;
      val = (*this)[idx].GetTotalFFactorFFactorMethod();
      break;
    case 14:
      if ((*this)[idx].IsExcluded() || !(*this)[idx].IsFFactorMethodValid())
        return kFALSE;
      val = (*this)[idx].GetTotalFFactorFFactorMethodErr();
      break;
    case 15:
      if ((*this)[idx].IsExcluded() || !(*this)[idx].IsBlindPixelMethodValid())
        return kFALSE;
      //      val = fBlindPixel->GetMeanFluxInsidePlexiglass()*area;
      val = 1.;
      break;
    case 16:
      if ((*this)[idx].IsExcluded() || !(*this)[idx].IsBlindPixelMethodValid())
        return kFALSE;
      //      val = fBlindPixel->GetMeanFluxErrInsidePlexiglass()*area;
      val = 1.;
      break;
    case 17:
      if ((*this)[idx].IsExcluded() || !(*this)[idx].IsBlindPixelMethodValid())
        return kFALSE;
      val = (*this)[idx].GetMeanConversionBlindPixelMethod();
      break;
    case 18:
      if ((*this)[idx].IsExcluded() || !(*this)[idx].IsBlindPixelMethodValid())
        return kFALSE;
      val = (*this)[idx].GetConversionBlindPixelMethodErr();
      break;
    case 19:
      if ((*this)[idx].IsExcluded() || !(*this)[idx].IsBlindPixelMethodValid())
        return kFALSE;
      val = (*this)[idx].GetTotalFFactorBlindPixelMethod();
      break;
    case 20:
      if ((*this)[idx].IsExcluded() || !(*this)[idx].IsBlindPixelMethodValid())
        return kFALSE;
      val = (*this)[idx].GetTotalFFactorBlindPixelMethodErr();
      break;
    case 21:
      if ((*this)[idx].IsExcluded() || !(*this)[idx].IsPINDiodeMethodValid())
        return kFALSE;
      //      val = fPINDiode->GetMeanFluxOutsidePlexiglass()*area;
      val = 1.;
      break;
    case 22:
      if ((*this)[idx].IsExcluded() || !(*this)[idx].IsPINDiodeMethodValid())
        return kFALSE;
      //      val = fPINDiode->GetMeanFluxErrOutsidePlexiglass()*area;
      val = 1.;
      break;
    case 23:
      if ((*this)[idx].IsExcluded() || !(*this)[idx].IsPINDiodeMethodValid())
        return kFALSE;
      val = (*this)[idx].GetMeanConversionPINDiodeMethod();
      break;
    case 24:
      if ((*this)[idx].IsExcluded() || !(*this)[idx].IsPINDiodeMethodValid())
        return kFALSE;
      val = (*this)[idx].GetConversionPINDiodeMethodErr();
      break;
    case 25:
      if ((*this)[idx].IsExcluded() || !(*this)[idx].IsPINDiodeMethodValid())
        return kFALSE;
      val = (*this)[idx].GetTotalFFactorPINDiodeMethod();
      break;
    case 26:
      if ((*this)[idx].IsExcluded() || !(*this)[idx].IsPINDiodeMethodValid())
        return kFALSE;
      val = (*this)[idx].GetTotalFFactorPINDiodeMethodErr();
      break;
    case 27:
      if ((*this)[idx].IsExcluded())
        val = 1.;
      else
        return kFALSE;
      break;
    case 28:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      if (!(*this)[idx].IsFitted())
        val = 1;
      else
        return kFALSE;
      break;
    case 29:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetHiGainNumPickup();
      break;
    case 30:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetLoGainNumPickup();
      break;
    case 31:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].IsHiGainSaturation();
      break;
    case 32:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      if ((*this)[idx].IsFFactorMethodValid())
        val = 1;
      else
        return kFALSE;
      break;
    case 33:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      if ((*this)[idx].IsBlindPixelMethodValid())
        val = 1;
      else
        return kFALSE;
      break;
    case 34:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      if ((*this)[idx].IsPINDiodeMethodValid())
        val = 1;
      else
        return kFALSE;
      break;
    case 35:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetPed();
      break;
    case 36:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetPedErr();
      break;
    case 37:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetPedRms();
      break;
    case 38:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetPedErr()/2.;
      break;
    case 39:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetAbsTimeMean();
      break;
    case 40:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetAbsTimeRms();
      break;
    default:
      return kFALSE;
    }

  return val!=-1.;

}

// --------------------------------------------------------------------------
//
// What MHCamera needs in order to draw an individual pixel in the camera
//
void MCalibrationChargeCam::DrawPixelContent(Int_t idx) const
{
    if (idx == -1)
        fAverageInnerPix->DrawClone();
    if (idx == -2)
        fAverageOuterPix->DrawClone();

    (*this)[idx].DrawClone();
}

//
// Calculate the weighted mean of the phe's of all inner and outer pixels, respectively.
// Bad pixels are excluded from the calculation. Two loops are performed to exclude pixels 
// which are fPheFFactorRelLimit sigmas from the mean.
//
Bool_t MCalibrationChargeCam::CalcMeanFluxPhotonsFFactorMethod(const MGeomCam &geom, MBadPixelsCam &bad)
{

  const Float_t avQERelVar = fAverageQEVar / (fAverageQE * fAverageQE ); 

  Float_t sumweightsinner = 0.;
  Float_t sumphesinner    = 0.;
  Float_t sumweightsouter = 0.;
  Float_t sumphesouter    = 0.;
  Int_t   validinner      = 0;
  Int_t   validouter      = 0;

  TIter Next(fPixels);
  MCalibrationChargePix *pix;
  while ((pix=(MCalibrationChargePix*)Next()))
    {

      if (!pix->IsFFactorMethodValid())
        continue;

      const Int_t idx = pix->GetPixId();

      if(!bad[idx].IsCalibrationResultOK())
        continue;
      
      const Float_t nphe    = pix->GetPheFFactorMethod();
      const Float_t npheerr = pix->GetPheFFactorMethodErr();
      const Float_t ratio   = geom.GetPixRatio(idx);

      if (npheerr > 0.)
        {
          const Float_t weight = nphe*nphe/npheerr/npheerr;
          // 
          // first the inner pixels:
          //
          if (ratio == 1.)
            {
              sumweightsinner += weight;
              sumphesinner    += weight*nphe;
              validinner++;
            }
          else
            {
              // 
              // now the outers
              //
              sumweightsouter += weight;
              sumphesouter    += weight*nphe;
              validouter++;
            }
        } /* if npheerr != 0 */
    } /* while ((pix=(MCalibrationChargePix*)Next())) */

  if (sumweightsinner <= 0. || sumphesinner <= 0.)
    {
      *fLog << warn << " Mean number of phe's from inner pixels cannot be calculated: "
            << " Sum of weights: " << sumweightsinner 
            << " Sum of weighted phes: " << sumphesinner << endl;
      return kFALSE;
    }
  else
    {
      fMeanFluxPhesInnerPixel    = sumphesinner/sumweightsinner;
      fMeanFluxPhesInnerPixelVar = (1./sumweightsinner)*fMeanFluxPhesInnerPixel*fMeanFluxPhesInnerPixel;
    }

  if (sumweightsouter <= 0. || sumphesouter <= 0.)
    {
      *fLog << warn << " Mean number of phe's from outer pixels cannot be calculated: "
            << " Sum of weights or sum of weighted phes is 0. " << endl;
    }
  else
    {
      fMeanFluxPhesOuterPixel    = sumphesouter/sumweightsouter;
      fMeanFluxPhesOuterPixelVar = (1./sumweightsouter)*fMeanFluxPhesOuterPixel*fMeanFluxPhesOuterPixel;
    }

  Float_t meanFluxPhotonsRelVar  = fMeanFluxPhesInnerPixelVar
                                / (fMeanFluxPhesInnerPixel * fMeanFluxPhesInnerPixel);

  fMeanFluxPhotonsInnerPixel    =  fMeanFluxPhesInnerPixel/fAverageQE;
  fMeanFluxPhotonsInnerPixelVar =  (meanFluxPhotonsRelVar + avQERelVar)
                                 * fMeanFluxPhotonsInnerPixel * fMeanFluxPhotonsInnerPixel;

  fMeanFluxPhotonsOuterPixel    = 4. *fMeanFluxPhotonsInnerPixel;
  fMeanFluxPhotonsOuterPixelVar = 16.*fMeanFluxPhotonsInnerPixelVar;  
  
  *fLog << inf << " Mean number of photo-electrons from inner pixels (F-Factor Method): "
        << fMeanFluxPhesInnerPixel << " +- " << GetMeanFluxPhesInnerPixelErr() << endl;

  *fLog << inf << " Mean number of photons from inner pixels (F-Factor Method): "
        << fMeanFluxPhotonsInnerPixel << " +- " << GetMeanFluxPhotonsInnerPixelErr() << endl;

  //
  // Here starts the second loop discarting pixels out of the range:
  //
  const Float_t innervar = (Float_t)validinner*fPheFFactorRelVarLimit*fMeanFluxPhesInnerPixelVar;
  const Float_t outervar = (Float_t)validouter*fPheFFactorRelVarLimit*fMeanFluxPhesOuterPixelVar;
  
  Float_t innererr;
  Float_t outererr;

  if (innervar > 0.)
    innererr = TMath::Sqrt(innervar);
  else
    {
      *fLog << warn << GetDescriptor() << ": Variance of mean number of photo-electrons for inner pixels "
            << " is smaller than 0. " << endl;
      SetFFactorMethodValid(kFALSE);
      return kFALSE;
    }
  
  if (outervar > 0.)
    outererr = TMath::Sqrt(outervar);
  else
    {
      *fLog << warn << GetDescriptor() << ": Variance of mean number of photo-electrons for outer pixels "
            << " is smaller than 0. " << endl;
      SetFFactorMethodValid(kFALSE);
      return kFALSE;
    }

  const Float_t lowerpheinnerlimit = fMeanFluxPhesInnerPixel - innererr;
  const Float_t upperpheinnerlimit = fMeanFluxPhesInnerPixel + innererr;

  const Float_t lowerpheouterlimit = fMeanFluxPhesOuterPixel - outererr;
  const Float_t upperpheouterlimit = fMeanFluxPhesOuterPixel + outererr;

  sumweightsinner = 0.;
  sumphesinner    = 0.;
  sumweightsouter = 0.;
  sumphesouter    = 0.;
   
  TIter Next2(fPixels);
  MCalibrationChargePix *pix2;
  while ((pix2=(MCalibrationChargePix*)Next2()))
    {

      if (!pix2->IsFFactorMethodValid())
        continue;

      const Int_t idx = pix2->GetPixId();

      if(!bad[idx].IsCalibrationResultOK())
        continue;

      const Float_t nphe    = pix2->GetPheFFactorMethod();
      const Float_t npheerr = pix2->GetPheFFactorMethodErr();
      const Float_t ratio   = geom.GetPixRatio(idx);

      if (npheerr > 0.)
        {
          const Float_t weight = nphe*nphe/npheerr/npheerr;
          // 
          // first the inner pixels:
          //
          if (ratio == 1.)
            {
              
              if (nphe < lowerpheinnerlimit || nphe > upperpheinnerlimit)
                {
                  bad[idx].SetUnsuitable(MBadPixelsPix::kUnreliableRun);
                  continue;
                }

              sumweightsinner += weight;
              sumphesinner    += weight*nphe;
            }
          else
            {
              // 
              // now the outers
              //
              if (nphe < lowerpheouterlimit || nphe > upperpheouterlimit)
                {
                  bad[idx].SetUnsuitable(MBadPixelsPix::kUnreliableRun);
                  continue;
                }

              sumweightsouter += weight;
              sumphesouter    += weight*nphe;
            }
        } /* if npheerr != 0 */
    } /* while ((pix2=(MCalibrationChargePix*)Next2())) */

  if (sumweightsinner <= 0. || sumphesinner <= 0.)
    {
      *fLog << warn << " Mean number of phe's from inner pixels cannot be calculated: "
            << " Sum of weights: " << sumweightsinner 
            << " Sum of weighted phes: " << sumphesinner << endl;
      return kFALSE;
    }
  else
    {
      fMeanFluxPhesInnerPixel    = sumphesinner/sumweightsinner;
      fMeanFluxPhesInnerPixelVar = (1./sumweightsinner)*fMeanFluxPhesInnerPixel*fMeanFluxPhesInnerPixel;

    }

  if (sumweightsouter <= 0. || sumphesouter <= 0.)
    {
      *fLog << warn << " Mean number of phe's from outer pixels cannot be calculated: "
            << " Sum of weights or sum of weighted phes is 0. " << endl;
    }
  else
    {
      fMeanFluxPhesOuterPixel    = sumphesouter/sumweightsouter;
      fMeanFluxPhesOuterPixelVar = (1./sumweightsouter)*fMeanFluxPhesOuterPixel*fMeanFluxPhesOuterPixel;
    }

  meanFluxPhotonsRelVar = fMeanFluxPhesInnerPixelVar
                       / (fMeanFluxPhesInnerPixel    * fMeanFluxPhesInnerPixel);

  fMeanFluxPhotonsInnerPixel    =  fMeanFluxPhesInnerPixel/fAverageQE;
  fMeanFluxPhotonsInnerPixelVar =  (meanFluxPhotonsRelVar + avQERelVar)
                                 * fMeanFluxPhotonsInnerPixel * fMeanFluxPhotonsInnerPixel;

  fMeanFluxPhotonsOuterPixel    = 4. *fMeanFluxPhotonsInnerPixel;
  fMeanFluxPhotonsOuterPixelVar = 16.*fMeanFluxPhotonsInnerPixelVar;  

  *fLog << inf << " Mean number of photo-electrons from inner pixels (F-Factor Method): "
        << fMeanFluxPhesInnerPixel << " +- " << GetMeanFluxPhesInnerPixelErr() << endl;

  *fLog << inf << " Mean number of photons from inner pixels (F-Factor Method): "
        << fMeanFluxPhotonsInnerPixel << " +- " << GetMeanFluxPhotonsInnerPixelErr() << endl;

  return kTRUE;
}

void MCalibrationChargeCam::ApplyFFactorCalibration(const MGeomCam &geom, const MBadPixelsCam &bad)
{

  const Float_t meanphotRelVar  = fMeanFluxPhotonsInnerPixelVar
                               /( fMeanFluxPhotonsInnerPixel    * fMeanFluxPhotonsInnerPixel );

  TIter Next(fPixels);
  MCalibrationChargePix *pix;
  while ((pix=(MCalibrationChargePix*)Next()))
    {

      const Int_t idx = pix->GetPixId();

      if (!(bad[idx].IsCalibrationResultOK()))
        {
          pix->SetFFactorMethodValid(kFALSE);
          continue;
        }

      if (!(pix->IsFFactorMethodValid()))
        continue;

      const Float_t ratio   = geom.GetPixRatio(idx);
      //
      // Calculate the conversion factor between PHOTONS and FADC counts
      //
      // Nphot = Nphe / avQE
      // conv  = Nphot / FADC counts
      //
      Float_t conv;
      
      if (ratio == 1.)
        conv               = fMeanFluxPhotonsInnerPixel / pix->GetMeanCharge(); 
      else
        conv               = fMeanFluxPhotonsOuterPixel / pix->GetMeanCharge(); 

      if (conv <= 0.)
        {
          pix->SetFFactorMethodValid(kFALSE);
          continue;
        }
      
      const Float_t chargeRelVar       =   pix->GetMeanChargeErr() * pix->GetMeanChargeErr()
                                       / ( pix->GetMeanCharge()    * pix->GetMeanCharge());
      const Float_t rsigmaChargeRelVar =   pix->GetRSigmaChargeErr() *  pix->GetRSigmaChargeErr()
                                        / (pix->GetRSigmaCharge()    * pix->GetRSigmaCharge()) ;

      const Float_t convrelvar = meanphotRelVar + chargeRelVar;

      if (convrelvar > fConvFFactorRelVarLimit)
        {
          *fLog << warn << GetDescriptor() << ": Conversion Factor F-Factor Method Rel. Variance: " 
                << convrelvar << " above limit of: " << fConvFFactorRelVarLimit 
                << " in pixel: " << idx << endl;
          pix->SetFFactorMethodValid(kFALSE);
          continue;
        }
      
      //
      // Calculate the Total F-Factor of the camera (in photons)
      //
      const Float_t totalFFactor  =  (pix->GetRSigmaCharge()/pix->GetMeanCharge())
                                    *TMath::Sqrt(fMeanFluxPhotonsInnerPixel);
      
      //
      // Calculate the error of the Total F-Factor of the camera ( in photons )
      //
      const Float_t totalFFactorVar = rsigmaChargeRelVar + chargeRelVar + meanphotRelVar;

      if (convrelvar > 0. && conv > 0.)
        pix->SetConversionFFactorMethod(conv,
                                        TMath::Sqrt(convrelvar)*conv,
                                        totalFFactor*TMath::Sqrt(conv));

      pix->SetTotalFFactorFFactorMethod(   totalFFactor   );

      if (totalFFactorVar > 0.)
        pix->SetTotalFFactorFFactorMethodErr(TMath::Sqrt(totalFFactorVar));      

      pix->SetFFactorMethodValid();
    }
}



void MCalibrationChargeCam::ApplyBlindPixelCalibration(const MGeomCam &geom, const MBadPixelsCam &bad, const MCalibrationChargeBlindPix &blindpix)
{

  Float_t flux    = blindpix.GetMeanFluxInsidePlexiglass();
  Float_t fluxerr = blindpix.GetMeanFluxErrInsidePlexiglass();

  TIter Next(fPixels);
  MCalibrationChargePix *pix;
  while ((pix=(MCalibrationChargePix*)Next()))
    {

      const Int_t idx = pix->GetPixId();

      if (!(bad[idx].IsCalibrationResultOK()))
        {
          pix->SetBlindPixelMethodValid(kFALSE);
          continue;
        }
      
      const Float_t charge    = pix->GetMeanCharge();
      const Float_t area      = geom[idx].GetA();
      const Float_t chargeerr = pix->GetMeanChargeErr();          
      
      const Float_t nphot      = flux    * area;
      const Float_t nphoterr   = fluxerr * area;
      const Float_t conversion = nphot/charge;
      Float_t conversionerr;
      
      conversionerr  = nphoterr/charge 
                     * nphoterr/charge ;
      conversionerr += chargeerr/charge
                     * chargeerr/charge
                     * conversion*conversion;
      conversionerr = TMath::Sqrt(conversionerr);

      const Float_t conversionsigma = 0.;
      
      pix->SetConversionBlindPixelMethod(conversion, conversionerr, conversionsigma);
      
      if (conversionerr/conversion < 0.1) 
        pix->SetBlindPixelMethodValid();
    }
}

void MCalibrationChargeCam::ApplyPINDiodeCalibration(const MGeomCam &geom, const MBadPixelsCam &bad, const MCalibrationChargePINDiode &pindiode)
{

  Float_t flux    = pindiode.GetMeanFluxOutsidePlexiglass();
  Float_t fluxerr = pindiode.GetMeanFluxErrOutsidePlexiglass();

  TIter Next(fPixels);
  MCalibrationChargePix *pix;
  while ((pix=(MCalibrationChargePix*)Next()))
    {

      const Int_t idx = pix->GetPixId();
      
      if (!(bad[idx].IsCalibrationResultOK()))
        {
          pix->SetPINDiodeMethodValid(kFALSE);
          continue;
        }
      
      const Float_t charge    = pix->GetMeanCharge();
      const Float_t area      = geom[idx].GetA();
      const Float_t chargeerr = pix->GetMeanChargeErr();          
      
      const Float_t nphot      = flux    * area;
      const Float_t nphoterr   = fluxerr * area;
      const Float_t conversion = nphot/charge;
      
      Float_t conversionerr;
      
      conversionerr  = nphoterr/charge 
                     * nphoterr/charge ;
      conversionerr += chargeerr/charge
                     * chargeerr/charge
                     * conversion*conversion;
      if (conversionerr > 0.)
        conversionerr = TMath::Sqrt(conversionerr);
      
      const Float_t conversionsigma = 0.;
      
      pix->SetConversionPINDiodeMethod(conversion, conversionerr, conversionsigma);
      
      if (conversionerr/conversion < 0.1) 
        pix->SetPINDiodeMethodValid();
      
    }
}


Bool_t MCalibrationChargeCam::GetConversionFactorBlindPixel(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
{

  mean  = (*this)[ipx].GetMeanConversionBlindPixelMethod();
  err   = (*this)[ipx].GetConversionBlindPixelMethodErr();
  sigma = (*this)[ipx].GetSigmaConversionBlindPixelMethod();

  return kTRUE;
}


Bool_t MCalibrationChargeCam::GetConversionFactorFFactor(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
{

  Float_t conv = (*this)[ipx].GetMeanConversionFFactorMethod();

  if (conv < 0.)
    return kFALSE;

  mean  = conv;
  err   = (*this)[ipx].GetConversionFFactorMethodErr();
  sigma = (*this)[ipx].GetSigmaConversionFFactorMethod();

  return kTRUE;
}


//-----------------------------------------------------------------------------------
//
// Calculates the conversion factor between the integral of FADCs slices 
// (as defined in the signal extractor MExtractSignal.cc)
// and the number of photons reaching the plexiglass for one Inner Pixel 
//
Bool_t MCalibrationChargeCam::GetConversionFactorPINDiode(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
{

  mean  = (*this)[ipx].GetMeanConversionPINDiodeMethod();
  err   = (*this)[ipx].GetConversionPINDiodeMethodErr();
  sigma = (*this)[ipx].GetSigmaConversionPINDiodeMethod();

  return kFALSE;

}

//-----------------------------------------------------------------------------------
//
// Calculates the best combination of the three used methods possible 
// between the integral of FADCs slices 
// (as defined in the signal extractor MExtractSignal.cc)
// and the number of photons reaching one Inner Pixel. 
// The procedure is not yet defined.
//
Bool_t MCalibrationChargeCam::GetConversionFactorCombined(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
{
  return kFALSE;

}

/*
void MCalibrationChargeCam::DrawHiLoFits()
{

  if (!fOffsets)
    fOffsets = new TH1D("pp","Offsets of the HiGain LoGain Fit",100,-600.,400.);
  if (!fSlopes)
    fSlopes  = new TH1D("mm","Slopes of the HiGain LoGain Fit",100,-2.,2.);
  if (!fOffvsSlope)
    fOffvsSlope = new TH2D("aa","Slopes vs Offsets of the HiGain LoGain Fit",100,-600.,400.,100,-2.,2.);

  TIter Next(fPixels);
  MCalibrationPix *pix;
  MHCalibrationPixel *hist;
  while ((pix=(MCalibrationPix*)Next()))
    {
      hist = pix->GetHist();
      hist->FitHiGainvsLoGain();
      fOffsets->Fill(hist->GetOffset(),1.);
      fSlopes->Fill(hist->GetSlope(),1.);
      fOffvsSlope->Fill(hist->GetOffset(),hist->GetSlope(),1.);
    }

   TCanvas *c1 = new TCanvas();

   c1->Divide(1,3);
   c1->cd(1);
   fOffsets->Draw();
   gPad->Modified();
   gPad->Update();

   c1->cd(2);
  fSlopes->Draw();
  gPad->Modified();
  gPad->Update();

  c1->cd(3);
  fOffvsSlope->Draw("col1");
  gPad->Modified();
  gPad->Update();
}

*/