/* ======================================================================== *\
!
! *
! * 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-2001
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                               
// MCalibrationCam                                               
//                                                               
// Hold the whole Calibration results of the camera:
//                                                               
// 1) MCalibrationCam initializes a TClonesArray whose elements are 
//    pointers to MCalibrationPix 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:
// 
// 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: Pixels with succeeded fit, but apparently wrong results
// 30: Pixels with un-expected behavior in the fourier spectrum (e.g. oscillations)
//
// 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: Error on the Absolute Arrival time of the signal
// 41: RMS of the Absolute Arrival time of the signal
// 42: Error on the RMS of the Absolute Arrival time of the signal
//
/////////////////////////////////////////////////////////////////////////////
#include "MCalibrationCam.h"

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

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

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

#include "MCalibrationPix.h"
#include "MCalibrationConfig.h"
#include "MCalibrationBlindPix.h"
#include "MCalibrationPINDiode.h"

#include "MHCalibrationPixel.h"

ClassImp(MCalibrationCam);

using namespace std;

const Int_t   MCalibrationCam::gkBlindPixelId   =  559;
const Int_t   MCalibrationCam::gkPINDiodeId     = 9999;
const Float_t MCalibrationCam::gkBlindPixelArea = 100;
const Float_t MCalibrationCam::gkPINDiodeArea   = 100;

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

    fPixels     = new TClonesArray("MCalibrationPix",1);
    fBlindPixel = new MCalibrationBlindPix();
    fPINDiode   = new MCalibrationPINDiode();

    Clear();
}

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

  //
  // delete fPixels should delete all Objects stored inside
  // 
  delete fPixels;
  delete fBlindPixel;
  delete fPINDiode;

  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 MCalibrationCam::InitSize(const UInt_t i)
{
  
  //
  // check if we have already initialized to size
  //
  if (CheckBounds(i))
    return;
  
  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 MCalibrationCam::GetSize() const
{
  return fPixels->GetEntriesFast();
}

// --------------------------------------------------------------------------
//
// Check if position i is inside the current bounds of the TClonesArray
//
Bool_t MCalibrationCam::CheckBounds(Int_t i) const 
{
  return i < GetSize();
} 


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

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


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

  fPixels->ForEach(TObject, Clear)();
  fBlindPixel->Clear();
  fPINDiode->Clear();

  fMeanFluxInsidePlexiglass          = -1.;
  fMeanFluxErrInsidePlexiglass       = -1.;
  fMeanFluxOutsidePlexiglass         = -1.;
  fMeanFluxErrOutsidePlexiglass      = -1.;

  fNumExcludedPixels                 = 0;

  CLRBIT(fFlags,kBlindPixelMethodValid);
  CLRBIT(fFlags,kPINDiodeMethodValid);
  CLRBIT(fFlags,kFluxInsidePlexiglassAvailable);
  CLRBIT(fFlags,kFluxOutsidePlexiglassAvailable);  

  return;
}

void MCalibrationCam::SetBlindPixelMethodValid(const Bool_t b)
{

  if (b) 
    SETBIT(fFlags, kBlindPixelMethodValid); 
  else    
    CLRBIT(fFlags, kBlindPixelMethodValid); 
  
}

void MCalibrationCam::SetPINDiodeMethodValid(const Bool_t b)
{

  if (b) 
    SETBIT(fFlags, kPINDiodeMethodValid); 
  else    
    CLRBIT(fFlags, kPINDiodeMethodValid); 
  
  
}

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

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


Bool_t  MCalibrationCam::IsFluxInsidePlexiglassAvailable()   const
{
  return TESTBIT(fFlags,kFluxInsidePlexiglassAvailable);
}

Bool_t  MCalibrationCam::IsFluxOutsidePlexiglassAvailable()   const
{
  return TESTBIT(fFlags,kFluxOutsidePlexiglassAvailable);
}



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

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

  TIter Next(fPixels);
  MCalibrationPix *pix;
  while ((pix=(MCalibrationPix*)Next()))
    {
      
      if (pix->IsChargeValid() && !pix->IsExcluded() && !pix->IsOscillating() && pix->IsFitted()) 
	{

	  *fLog << all << pix->GetPixId() << " Pedestals: " << pix->GetPed() << " +- " 
                << pix->GetPedRms() << " Reduced Charge: " << pix->GetCharge() << " +- " 
                << pix->GetSigmaCharge() << " Reduced Sigma: " << pix->GetRSigmaCharge() 
                << " Nr Phe's: " << pix->GetPheFFactorMethod() << endl;
          id++;
	}
    }
  
  *fLog << all << id << " succesful pixels :-))" << endl;
  id = 0;
  
  *fLog << all << endl;
  *fLog << all << "Pixels with errors:" << endl;
  *fLog << all << endl;
  
  TIter Next2(fPixels);
    while ((pix=(MCalibrationPix*)Next2()))
      {
        
        if (!pix->IsChargeValid() && !pix->IsExcluded() && !pix->IsFitted())
          {

            *fLog << all << pix->GetPixId() << " Pedestals: " << pix->GetPed() << " +- " 
                  << pix->GetPedRms() << " Reduced Charge: " << pix->GetCharge() << " +- " 
                  << pix->GetSigmaCharge() << " Reduced Sigma: " << pix->GetRSigmaCharge() << endl;
            id++;
          }
      }
    *fLog << all << id << " pixels with errors :-((" << endl;

  *fLog << all << endl;
  *fLog << all << "Pixels with oscillations:" << endl;
  *fLog << all << endl;
  
  id = 0;

  TIter Next3(fPixels);
  while ((pix=(MCalibrationPix*)Next3()))
    {
      
      if (pix->IsOscillating() && !pix->IsExcluded())
        {
          
          *fLog << all << pix->GetPixId() << " Pedestals: " << pix->GetPed() << " +- " 
                << pix->GetPedRms() << " Reduced Charge: " << pix->GetCharge() << " +- " 
                << pix->GetSigmaCharge() << " Reduced Sigma: " << pix->GetRSigmaCharge() << endl;
          id++;
        }
    }
  *fLog << all << id << " Oscillating pixels :-((" << endl;
  
    
  *fLog << all << endl;
  *fLog << all << "Excluded pixels:" << endl;
  *fLog << all << endl;
  
  TIter Next4(fPixels);
  while ((pix=(MCalibrationPix*)Next4()))
    if (pix->IsExcluded())
        *fLog << all << pix->GetPixId() << endl;

  *fLog << all << fNumExcludedPixels << " excluded pixels " << endl;
}

// --------------------------------------------------------------------------
//
// Return true if pixel is inside bounds of the TClonesArray fPixels
//
Bool_t MCalibrationCam::IsPixelUsed(Int_t idx) const 
{
  if (!CheckBounds(idx))
    return kFALSE;

  return kTRUE;
}

// --------------------------------------------------------------------------
//
// Return true if pixel has already been fitted once (independent of the result)
//
Bool_t MCalibrationCam::IsPixelFitted(Int_t idx) const 
{

  if (!CheckBounds(idx))
    return kFALSE;

  return (*this)[idx].IsFitted();
}

// --------------------------------------------------------------------------
//
// Sets the user ranges of all histograms such that 
// empty bins at the edges are not used. Additionally, it rebins the 
// histograms such that in total, 50 bins are used.
//
void MCalibrationCam::CutEdges()
{

  fBlindPixel->GetHist()->CutAllEdges();
  fPINDiode->GetHist()->CutAllEdges();

  TIter Next(fPixels);
  MCalibrationPix *pix;
  while ((pix=(MCalibrationPix*)Next()))
    {
      pix->GetHist()->CutAllEdges();
    }

  return;
}
  
// --------------------------------------------------------------------------
//
// 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: Pixels with succeeded fit, but apparently wrong results
// 30: Pixels with un-expected behavior in the fourier spectrum (e.g. oscillations)
//
// 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: Error on the Absolute Arrival time of the signal
// 41: RMS of the Absolute Arrival time of the signal
// 42: Error on the RMS of the Absolute Arrival time of the signal
//
Bool_t MCalibrationCam::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].GetCharge();
      break;
    case 1:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetErrCharge();
      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].GetErrSigmaCharge();
      break;
    case 4:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetChargeProb();
      break;
    case 5:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetRSigmaCharge();
      break;
    case 6:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetErrRSigmaCharge();
      break;
    case 7:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetRSigmaCharge() / (*this)[idx].GetCharge();
      break;
    case 8:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      // relative error RsigmaCharge square
      val =    (*this)[idx].GetErrRSigmaCharge()* (*this)[idx].GetErrRSigmaCharge() 
            / ((*this)[idx].GetRSigmaCharge()   * (*this)[idx].GetRSigmaCharge()   );
      // relative error Charge square
      val +=   (*this)[idx].GetErrCharge() * (*this)[idx].GetErrCharge()
            / ((*this)[idx].GetCharge()    * (*this)[idx].GetCharge()   );
      // calculate relative error out of squares
      val  =   TMath::Sqrt(val) ;
      // multiply with value to get absolute error
      val  *=  (*this)[idx].GetRSigmaCharge() / (*this)[idx].GetCharge();
      break;
    case 9:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetPheFFactorMethod();
      break;
    case 10:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetPheFFactorMethodError();
      break;
    case 11:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetMeanConversionFFactorMethod();
      break;
    case 12:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetErrorConversionFFactorMethod();
      break;
    case 13:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetTotalFFactorFFactorMethod();
      break;
    case 14:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetTotalFFactorErrorFFactorMethod();
      break;
    case 15:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = GetMeanFluxInsidePlexiglass()*area;
      break;
    case 16:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = GetMeanFluxInsidePlexiglass()*area;
      break;
    case 17:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetMeanConversionBlindPixelMethod();
      break;
    case 18:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetErrorConversionBlindPixelMethod();
      break;
    case 19:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetTotalFFactorBlindPixelMethod();
      break;
    case 20:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetTotalFFactorErrorBlindPixelMethod();
      break;
    case 21:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = GetMeanFluxOutsidePlexiglass()*area;
      break;
    case 22:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = GetMeanFluxOutsidePlexiglass()*area;
      break;
    case 23:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetMeanConversionPINDiodeMethod();
      break;
    case 24:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetErrorConversionPINDiodeMethod();
      break;
    case 25:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetTotalFFactorBlindPixelMethod();
      break;
    case 26:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetTotalFFactorErrorBlindPixelMethod();
      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;
      if (!(*this)[idx].IsFitted())
        return kFALSE;
      if (!(*this)[idx].IsChargeValid())
        val = 1;
      else
        return kFALSE;
      break;
    case 30:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      if ((*this)[idx].IsOscillating())
        val = 1;
      else
        return kFALSE;
      break;
    case 31:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      if ((*this)[idx].IsHiGainSaturation())
        val = 1;
      else
        return kFALSE;
      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 = 1.;
      //      val = (*this)[idx].GetPedError();
      break;
    case 37:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetPedRms();
      break;
    case 38:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = 1.;
      //      val = (*this)[idx].GetPedRmsError();
      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].GetAbsTimeMeanErr();
      break;
    case 41:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetAbsTimeRms();
      break;
    case 42:
      if ((*this)[idx].IsExcluded())
        return kFALSE;
      val = (*this)[idx].GetAbsTimeMeanErr()/TMath::Sqrt(2.);
      break;
    default:
      return kFALSE;
    }
  return val!=-1.;
}

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


// --------------------------------------------------------------------------
//
//
//
Bool_t MCalibrationCam::CalcFluxInsidePlexiglass()
{

  if (!fBlindPixel->IsFitOK())
    return kFALSE;
  
  const Float_t mean = fBlindPixel->GetLambda();
  const Float_t merr = fBlindPixel->GetErrLambda();

  //
  // Start calculation of number of photons 
  //
  // The blind pixel has exactly 100 mm^2 area (with negligible error), 
  //
  fMeanFluxInsidePlexiglass    = mean*gkBlindPixelArea;

  // Start calculation of number of photons relative Variance (!!)
  fMeanFluxErrInsidePlexiglass  = merr*merr/mean/mean;
  
  switch (fColor)
    {
    case kECGreen:
      fMeanFluxInsidePlexiglass    /= gkCalibrationBlindPixelQEGreen;   
      fMeanFluxErrInsidePlexiglass += gkCalibrationBlindPixelQEGreenError*gkCalibrationBlindPixelQEGreenError
                                    / gkCalibrationBlindPixelQEGreen  /   gkCalibrationBlindPixelQEGreen;   

      fMeanFluxInsidePlexiglass  *= TMath::Power(10,gkCalibrationBlindPixelAttGreen); // correct for absorption 
      // attenuation has negligible error
      break;
    case kECBlue:
      fMeanFluxInsidePlexiglass    /= gkCalibrationBlindPixelQEBlue;   
      fMeanFluxErrInsidePlexiglass += gkCalibrationBlindPixelQEBlueError*gkCalibrationBlindPixelQEBlueError
                                    / gkCalibrationBlindPixelQEBlue  /   gkCalibrationBlindPixelQEBlue;   

      fMeanFluxInsidePlexiglass *= TMath::Power(10,gkCalibrationBlindPixelAttBlue); // correct for absorption 
      // attenuation has negligible error
      break;
    case kECUV:
      fMeanFluxInsidePlexiglass    /= gkCalibrationBlindPixelQEUV;   
      fMeanFluxErrInsidePlexiglass += gkCalibrationBlindPixelQEUVError*gkCalibrationBlindPixelQEUVError
                                    / gkCalibrationBlindPixelQEUV  /   gkCalibrationBlindPixelQEUV;   

      fMeanFluxInsidePlexiglass *= TMath::Power(10,gkCalibrationBlindPixelAttUV); // correct for absorption 
      // attenuation has negligible error
      break;
    case kECCT1:
    default:
      fMeanFluxInsidePlexiglass    /= gkCalibrationBlindPixelQECT1;   
      fMeanFluxErrInsidePlexiglass += gkCalibrationBlindPixelQECT1Error*gkCalibrationBlindPixelQECT1Error
                                    / gkCalibrationBlindPixelQECT1  /   gkCalibrationBlindPixelQECT1;   

      fMeanFluxInsidePlexiglass *= TMath::Power(10,gkCalibrationBlindPixelAttCT1); // correct for absorption 
      // attenuation has negligible error
      break;
    }

  *fLog << inf << endl;
  *fLog << inf << " Photon flux [ph/mm^2] inside Plexiglass: " 
        << fMeanFluxInsidePlexiglass << endl;

  if (fMeanFluxInsidePlexiglass > 0.)
    SETBIT(fFlags,kFluxInsidePlexiglassAvailable);  
  else
    {
      CLRBIT(fFlags,kFluxInsidePlexiglassAvailable);        
      return kFALSE;
    }

  if (fMeanFluxErrInsidePlexiglass < 0.)
    {
      *fLog << warn << " Relative Variance on Photon flux inside Plexiglass: " 
            << fMeanFluxErrInsidePlexiglass << endl;
      CLRBIT(fFlags,kFluxInsidePlexiglassAvailable);        
      return kFALSE;
    }

  // Finish calculation of errors -> convert from relative variance to absolute error
  fMeanFluxErrInsidePlexiglass = TMath::Sqrt(fMeanFluxErrInsidePlexiglass);
  fMeanFluxErrInsidePlexiglass *= fMeanFluxInsidePlexiglass;

  *fLog << inf << " Error on photon flux [ph/mm^2] inside Plexiglass: " 
        << fMeanFluxErrInsidePlexiglass << endl;
  *fLog << inf << endl;

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

      if(pix->IsChargeValid())
        {
          
          const Int_t idx = pix->GetPixId();

          const Float_t charge    = pix->GetCharge();
          const Float_t area      = (*fGeomCam)[idx].GetA();
          const Float_t chargeerr = pix->GetErrCharge();          

          const Float_t nphot      = fMeanFluxInsidePlexiglass*area;
          const Float_t nphoterr   = fMeanFluxErrInsidePlexiglass*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();
        }
    }
  return kTRUE;
}


Bool_t MCalibrationCam::CalcFluxOutsidePlexiglass()
{

  if (!fPINDiode->IsChargeFitValid())
    return kFALSE;
  
  const Float_t mean = fPINDiode->GetCharge();
  const Float_t merr = fPINDiode->GetErrCharge();

  // Start calculation of number of photons per mm^2 on the camera
  fMeanFluxOutsidePlexiglass  = mean * gkPINDiodeArea;
  // Correct for the distance between camera and PIN Diode and for different areas.
  fMeanFluxOutsidePlexiglass *= gkCalibrationFluxCameravsPINDiode;

  // Start calculation of number of photons relative Variance (!!)
  fMeanFluxErrOutsidePlexiglass  = merr*merr/mean/mean;
  fMeanFluxErrOutsidePlexiglass += gkCalibrationFluxCameravsPINDiodeError*gkCalibrationFluxCameravsPINDiodeError
                                 / gkCalibrationFluxCameravsPINDiode/gkCalibrationFluxCameravsPINDiode;
  
  switch (fColor)
    {
    case kECGreen:
      fMeanFluxOutsidePlexiglass    /= gkCalibrationPINDiodeQEGreen;
      fMeanFluxErrOutsidePlexiglass += gkCalibrationPINDiodeQEGreenError*gkCalibrationPINDiodeQEGreenError
                                     / gkCalibrationPINDiodeQEGreen/gkCalibrationPINDiodeQEGreen;
      break;
    case kECBlue:
      fMeanFluxOutsidePlexiglass    /= gkCalibrationPINDiodeQEBlue;
      fMeanFluxErrOutsidePlexiglass += gkCalibrationPINDiodeQEBlueError*gkCalibrationPINDiodeQEBlueError
                                     / gkCalibrationPINDiodeQEBlue/gkCalibrationPINDiodeQEBlue;
      break; 
    case kECUV:
      fMeanFluxOutsidePlexiglass    /= gkCalibrationPINDiodeQEUV;
      fMeanFluxErrOutsidePlexiglass += gkCalibrationPINDiodeQEUVError*gkCalibrationPINDiodeQEUVError
                                     / gkCalibrationPINDiodeQEUV/gkCalibrationPINDiodeQEUV;
      break;
    case kECCT1:
    default:
      fMeanFluxOutsidePlexiglass    /= gkCalibrationPINDiodeQECT1;
      fMeanFluxErrOutsidePlexiglass += gkCalibrationPINDiodeQECT1Error*gkCalibrationPINDiodeQECT1Error
                                     / gkCalibrationPINDiodeQECT1/gkCalibrationPINDiodeQECT1;
      break;
    }


  *fLog << inf << endl;
  *fLog << inf << " Mean Photon flux [ph/mm^2] outside Plexiglass: " 
        << fMeanFluxOutsidePlexiglass << endl;

  if (fMeanFluxOutsidePlexiglass > 0.)
    SETBIT(fFlags,kFluxOutsidePlexiglassAvailable);  
  else
    {
      CLRBIT(fFlags,kFluxOutsidePlexiglassAvailable);        
      return kFALSE;
    }

  if (fMeanFluxErrOutsidePlexiglass < 0.)
    {
      *fLog << warn << "Relative Variance on Photon flux outside Plexiglass: " 
            << fMeanFluxErrOutsidePlexiglass << endl;
      CLRBIT(fFlags,kFluxOutsidePlexiglassAvailable);        
      return kFALSE;
    }

  // Finish calculation of errors -> convert from relative variance to absolute error
  fMeanFluxErrOutsidePlexiglass = TMath::Sqrt(fMeanFluxErrOutsidePlexiglass);
  fMeanFluxErrOutsidePlexiglass *= fMeanFluxOutsidePlexiglass;

  *fLog << inf << " Error on Photon flux [ph/mm^2] outside Plexiglass: " 
        << fMeanFluxErrOutsidePlexiglass << endl;
  *fLog << inf << endl;

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

      if (pix->IsChargeValid())
        {

          const Int_t idx = pix->GetPixId();

          const Float_t charge    = pix->GetCharge();
          const Float_t area      = (*fGeomCam)[idx].GetA();
          const Float_t chargeerr = pix->GetErrCharge();          

          const Float_t nphot      = fMeanFluxOutsidePlexiglass*area;
          const Float_t nphoterr   = fMeanFluxErrOutsidePlexiglass*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();
          
        }
    }
  return kTRUE;
}



Bool_t MCalibrationCam::GetConversionFactorBlindPixel(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
{
  
  if (ipx < 0 || !IsPixelFitted(ipx))
    return kFALSE;

  if (!IsFluxInsidePlexiglassAvailable())
    if (!CalcFluxInsidePlexiglass())
      return kFALSE;

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

  return kTRUE;
}


Bool_t MCalibrationCam::GetConversionFactorFFactor(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
{
  
  if (ipx < 0 || !IsPixelFitted(ipx))
    return kFALSE;

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

  if (conv < 0.)
    return kFALSE;

  mean  = conv;
  err   = (*this)[ipx].GetErrorConversionFFactorMethod();
  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 
//
// FIXME: The PINDiode is still not working and so is the code 
//
Bool_t MCalibrationCam::GetConversionFactorPINDiode(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
{

  if (ipx < 0 || !IsPixelFitted(ipx))
    return kFALSE;

  if (!IsFluxOutsidePlexiglassAvailable())
    if (!CalcFluxOutsidePlexiglass())
      return kFALSE;

  mean  = (*this)[ipx].GetMeanConversionPINDiodeMethod();
  err   = (*this)[ipx].GetErrorConversionPINDiodeMethod();
  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.
//
// FIXME: The PINDiode is still not working and so is the code 
//
Bool_t MCalibrationCam::GetConversionFactorCombined(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
{

  if (ipx < 0 || !IsPixelFitted(ipx))
    return kFALSE;

  return kFALSE;

}


void MCalibrationCam::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();
}

