/* ======================================================================== *\
!
! *
! * 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: Pixels with valid calibration by the F-Factor-Method
// 16: Mean conversion factor obtained with the Blind Pixel method
// 17: Error on the mean conversion factor obtained with the Blind Pixel method
// 18: Overall F-Factor of the readout obtained with the Blind Pixel method
// 19: Error on Overall F-Factor of the readout obtained with the Blind Pixel method
// 20: Pixels with valid calibration by the Blind Pixel-Method
// 21: Mean conversion factor obtained with the PIN Diode method
// 22: Error on the mean conversion factor obtained with the PIN Diode method
// 23: Overall F-Factor of the readout obtained with the PIN Diode method
// 24: Error on Overall F-Factor of the readout obtained with the PIN Diode method
// 25: Pixels with valid calibration by the PIN Diode-Method
//
// Localized defects:
// ==================
//
// 26: Excluded Pixels
// 27: Number of probable pickup events in the Hi Gain 
// 28: Number of probable pickup events in the Lo Gain
//
// Other classifications of pixels:
// ================================
//
// 29: Pixels with saturated Hi-Gain
//
// Used Pedestals:
// ===============
//
// 30: Mean Pedestal over the entire range of signal extraction
// 31: Error on the Mean Pedestal over the entire range of signal extraction
// 32: Pedestal RMS over the entire range of signal extraction
// 33: Error on the Pedestal RMS over the entire range of signal extraction
//
// Calculated absolute arrival times (very low precision!):
// ========================================================
//
// 34: Absolute Arrival time of the signal
// 35: RMS of the Absolute Arrival time of the signal
//
/////////////////////////////////////////////////////////////////////////////
#include "MCalibrationChargeCam.h"
#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 "MBadPixelsCam.h"
#include "MBadPixelsPix.h"

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

ClassImp(MCalibrationChargeCam);

using namespace std;

const Float_t MCalibrationChargeCam::fgAverageQE                = 0.18;     
const Float_t MCalibrationChargeCam::fgAverageQEErr             = 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);
    fAverageAreas     = new TClonesArray("MCalibrationChargePix",1);
    fAverageSectors   = new TClonesArray("MCalibrationChargePix",1);

    fFluxPhesInnerPixel       = 0.;
    fFluxPhesInnerPixelVar    = 0.;
    fFluxPhesOuterPixel       = 0.;
    fFluxPhesOuterPixelVar    = 0.;
    
    CLRBIT(fFlags,kBlindPixelMethodValid);
    CLRBIT(fFlags,kFFactorMethodValid);
    CLRBIT(fFlags,kPINDiodeMethodValid);
  
    SetAverageQE();
    SetConvFFactorRelErrLimit();
    SetPheFFactorRelErrLimit();
}

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

  
  if (fOffsets)
    delete fOffsets;
  if (fSlopes)
    delete fSlopes;
  if (fOffvsSlope)
    delete fOffvsSlope;

}


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

  fFluxPhesInnerPixel       = 0.;
  fFluxPhesInnerPixelVar    = 0.;
  fFluxPhesOuterPixel       = 0.;
  fFluxPhesOuterPixelVar    = 0.;

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

  MCalibrationCam::Clear();

  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::GetFluxPhesInnerPixelErr()   const
{
  if (fFluxPhesInnerPixelVar <= 0.)
    return -1.;
  return TMath::Sqrt(fFluxPhesInnerPixelVar);
}

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

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

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

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->GetMean()    << " +- " << pix->GetSigma() 
		<< "   Reduced Sigma: "        << pix->GetRSigma() 
		<< "   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;

  TIter Next5(fAverageAreas);
  while ((pix=(MCalibrationChargePix*)Next5()))
  {
    *fLog << all << "Average Area:" 
          << "   Ped.  Rms: "            << pix->GetPedRms()        << " +- " << pix->GetPedRmsErr() 
          << "   Mean signal: "          << pix->GetMean()    << " +- " << pix->GetMeanErr() 
          << "   Sigma signal: "         << pix->GetSigma()    << " +- "<< pix->GetSigmaErr() 
          << "   Reduced Sigma: "        << pix->GetRSigma() 
          << "   Nr Phe's: "             << pix->GetPheFFactorMethod() 
          << endl;
  }

  TIter Next6(fAverageSectors);
  while ((pix=(MCalibrationChargePix*)Next5()))
  {
    *fLog << all << "Average Sector:" 
          << "   Ped.  Rms: "            << pix->GetPedRms()        << " +- " << pix->GetPedRmsErr() 
          << "   Mean signal: "          << pix->GetMean()    << " +- " << pix->GetMeanErr() 
          << "   Sigma signal: "         << pix->GetSigma()    << " +- "<< pix->GetSigmaErr() 
          << "   Reduced Sigma: "        << pix->GetRSigma() 
          << "   Nr Phe's: "             << pix->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 
//
// 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: Pixels with valid calibration by the F-Factor-Method
// 16: Mean conversion factor obtained with the Blind Pixel method
// 17: Error on the mean conversion factor obtained with the Blind Pixel method
// 18: Overall F-Factor of the readout obtained with the Blind Pixel method
// 19: Error on Overall F-Factor of the readout obtained with the Blind Pixel method
// 20: Pixels with valid calibration by the Blind Pixel-Method
// 21: Mean conversion factor obtained with the PIN Diode method
// 22: Error on the mean conversion factor obtained with the PIN Diode method
// 23: Overall F-Factor of the readout obtained with the PIN Diode method
// 24: Error on Overall F-Factor of the readout obtained with the PIN Diode method
// 25: Pixels with valid calibration by the PIN Diode-Method
//
// Localized defects:
// ==================
//
// 26: Excluded Pixels
// 27: Number of probable pickup events in the Hi Gain 
// 28: Number of probable pickup events in the Lo Gain
//
// Other classifications of pixels:
// ================================
//
// 29: Pixels with saturated Hi-Gain
//
// Used Pedestals:
// ===============
//
// 30: Mean Pedestal over the entire range of signal extraction
// 31: Error on the Mean Pedestal over the entire range of signal extraction
// 32: Pedestal RMS over the entire range of signal extraction
// 33: Error on the Pedestal RMS over the entire range of signal extraction
//
// Calculated absolute arrival times (very low precision!):
// ========================================================
//
// 34: Absolute Arrival time of the signal
// 35: 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;

 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*this)[idx];

  switch (type)
    {
    case 0:
      if (pix.IsExcluded())
        return kFALSE;
      val = pix.GetMean();
      break;
    case 1:
      if (pix.IsExcluded())
        return kFALSE;
      val = pix.GetMeanErr();
      break;
    case 2:
      if (pix.IsExcluded())
        return kFALSE;
      val = pix.GetSigma();
      break;
    case 3:
      if (pix.IsExcluded())
        return kFALSE;
      val = pix.GetSigmaErr();
      break;
    case 4:
      if (pix.IsExcluded())
        return kFALSE;
      val = pix.GetProb();
      break;
    case 5:
      if (pix.IsExcluded())
        return kFALSE;
      if (pix.GetRSigma() == -1.)
	  return kFALSE;
      val = pix.GetRSigma();
      break;
    case 6:
      if (pix.IsExcluded())
        return kFALSE;    
      if (pix.GetRSigma() == -1.)
	  return kFALSE;
      val = pix.GetRSigmaErr();
      break;
    case 7:
      if (pix.IsExcluded())
        return kFALSE;
      if (pix.GetRSigma() == -1.)
        return kFALSE;
      if (pix.GetMean() == 0.)
        return kFALSE;
      val = pix.GetRSigma() / pix.GetMean();
      break;
    case 8:
      if (pix.IsExcluded())
        return kFALSE;
      if (pix.GetRSigma() <= 0.)
	  return kFALSE;
      if (pix.GetMean() <= 0.)
        return kFALSE;
      if (pix.GetRSigmaErr() <= 0.)
	  return kFALSE;
      if (pix.GetMeanErr() <= 0.)
        return kFALSE;
      // relative error Rsigma square
      val =    pix.GetRSigmaErr()* pix.GetRSigmaErr() 
            / (pix.GetRSigma()   * pix.GetRSigma()   );
      // relative error  square
      val +=   pix.GetMeanErr() * pix.GetMeanErr()
            / (pix.GetMean()    * pix.GetMean()   );
      // calculate relative error out of squares
      val  =   TMath::Sqrt(val) ;
      // multiply with value to get absolute error
      val  *=  pix.GetRSigma() / pix.GetMean();
      break;
    case 9:
      if (pix.IsExcluded() || !pix.IsFFactorMethodValid())
        return kFALSE;
      val = pix.GetPheFFactorMethod();
      break;
    case 10:
      if (pix.IsExcluded() || !pix.IsFFactorMethodValid())
        return kFALSE;
      val = pix.GetPheFFactorMethodErr();
      break;
    case 11:
      if (pix.IsExcluded() || !pix.IsFFactorMethodValid())
        return kFALSE;
      val = pix.GetMeanConversionFFactorMethod();
      break;
    case 12:
      if (pix.IsExcluded() || !pix.IsFFactorMethodValid())
        return kFALSE;
      val = pix.GetConversionFFactorMethodErr();
      break;
    case 13:
      if (pix.IsExcluded() || !pix.IsFFactorMethodValid())
        return kFALSE;
      val = pix.GetTotalFFactorFFactorMethod();
      break;
    case 14:
      if (pix.IsExcluded() || !pix.IsFFactorMethodValid())
        return kFALSE;
      val = pix.GetTotalFFactorFFactorMethodErr();
      break;
    case 15:
      if (pix.IsExcluded())
        return kFALSE;
      if (pix.IsFFactorMethodValid())
        val = 1;
      else
        return kFALSE;
      break;
    case 16:
      if (pix.IsExcluded() || !pix.IsBlindPixelMethodValid())
        return kFALSE;
      val = pix.GetMeanConversionBlindPixelMethod();
      break;
    case 17:
      if (pix.IsExcluded() || !pix.IsBlindPixelMethodValid())
        return kFALSE;
      val = pix.GetConversionBlindPixelMethodErr();
      break;
    case 18:
      if (pix.IsExcluded() || !pix.IsBlindPixelMethodValid())
        return kFALSE;
      val = pix.GetTotalFFactorBlindPixelMethod();
      break;
    case 19:
      if (pix.IsExcluded() || !pix.IsBlindPixelMethodValid())
        return kFALSE;
      val = pix.GetTotalFFactorBlindPixelMethodErr();
      break;
    case 20:
      if (pix.IsExcluded())
        return kFALSE;
      if (pix.IsBlindPixelMethodValid())
        val = 1;
      else
        return kFALSE;
      break;
    case 21:
      if (pix.IsExcluded() || !pix.IsPINDiodeMethodValid())
        return kFALSE;
      val = pix.GetMeanConversionPINDiodeMethod();
      break;
    case 22:
      if (pix.IsExcluded() || !pix.IsPINDiodeMethodValid())
        return kFALSE;
      val = pix.GetConversionPINDiodeMethodErr();
      break;
    case 23:
      if (pix.IsExcluded() || !pix.IsPINDiodeMethodValid())
        return kFALSE;
      val = pix.GetTotalFFactorPINDiodeMethod();
      break;
    case 24:
      if (pix.IsExcluded() || !pix.IsPINDiodeMethodValid())
        return kFALSE;
      val = pix.GetTotalFFactorPINDiodeMethodErr();
      break;
    case 25:
      if (pix.IsExcluded())
        return kFALSE;
      if (pix.IsPINDiodeMethodValid())
        val = 1;
      else
        return kFALSE;
      break;
    case 26:
      if (pix.IsExcluded())
        val = 1.;
      else
        return kFALSE;
      break;
    case 27:
      if (pix.IsExcluded())
        return kFALSE;
      val = pix.GetHiGainNumPickup();
      break;
    case 28:
      if (pix.IsExcluded())
        return kFALSE;
      val = pix.GetLoGainNumPickup();
      break;
    case 29:
      if (pix.IsExcluded())
        return kFALSE;
      val = pix.IsHiGainSaturation();
      break;
    case 30:
      if (pix.IsExcluded())
        return kFALSE;
      val = pix.GetPed();
      break;
    case 31:
      if (pix.IsExcluded())
        return kFALSE;
      val = pix.GetPedErr();
      break;
    case 32:
      if (pix.IsExcluded())
        return kFALSE;
      val = pix.GetPedRms();
      break;
    case 33:
      if (pix.IsExcluded())
        return kFALSE;
      val = pix.GetPedErr()/2.;
      break;
    case 34:
      if (pix.IsExcluded())
        return kFALSE;
      val = pix.GetAbsTimeMean();
      break;
    case 35:
      if (pix.IsExcluded())
        return kFALSE;
      val = pix.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
{
  MCalibrationChargePix &pix = (MCalibrationChargePix&)(*this)[idx];
  pix.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::CalcFluxPhotonsFFactorMethod(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
    {
      fFluxPhesInnerPixel    = sumphesinner/sumweightsinner;
      fFluxPhesInnerPixelVar = (1./sumweightsinner)*fFluxPhesInnerPixel*fFluxPhesInnerPixel;
    }

  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
    {
      fFluxPhesOuterPixel    = sumphesouter/sumweightsouter;
      fFluxPhesOuterPixelVar = (1./sumweightsouter)*fFluxPhesOuterPixel*fFluxPhesOuterPixel;
    }

  Float_t meanFluxPhotonsRelVar  = fFluxPhesInnerPixelVar
                                / (fFluxPhesInnerPixel * fFluxPhesInnerPixel);

  fFluxPhotonsInnerPixel    =  fFluxPhesInnerPixel/fAverageQE;
  fFluxPhotonsInnerPixelVar =  (meanFluxPhotonsRelVar + avQERelVar)
                                 * fFluxPhotonsInnerPixel * fFluxPhotonsInnerPixel;

  fFluxPhotonsOuterPixel    = 4. *fFluxPhotonsInnerPixel;
  fFluxPhotonsOuterPixelVar = 16.*fFluxPhotonsInnerPixelVar;  
  
  *fLog << inf << " Mean number of photo-electrons from inner pixels (F-Factor Method): "
        << fFluxPhesInnerPixel << " +- " << GetFluxPhesInnerPixelErr() << endl;

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

  //
  // Here starts the second loop discarting pixels out of the range:
  //
  const Float_t innervar = (Float_t)validinner*fPheFFactorRelVarLimit*fFluxPhesInnerPixelVar;
  const Float_t outervar = (Float_t)validouter*fPheFFactorRelVarLimit*fFluxPhesOuterPixelVar;
  
  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 = fFluxPhesInnerPixel - innererr;
  const Float_t upperpheinnerlimit = fFluxPhesInnerPixel + innererr;

  const Float_t lowerpheouterlimit = fFluxPhesOuterPixel - outererr;
  const Float_t upperpheouterlimit = fFluxPhesOuterPixel + 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].SetUncalibrated(MBadPixelsPix::kDeviatingNumPhes);
                  bad[idx].SetUnsuitable(  MBadPixelsPix::kUnreliableRun);
                  continue;
                }

              sumweightsinner += weight;
              sumphesinner    += weight*nphe;
            }
          else
            {
              // 
              // now the outers
              //
              if (nphe < lowerpheouterlimit || nphe > upperpheouterlimit)
                {
                  bad[idx].SetUncalibrated(MBadPixelsPix::kDeviatingNumPhes);
                  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
    {
      fFluxPhesInnerPixel    = sumphesinner/sumweightsinner;
      fFluxPhesInnerPixelVar = (1./sumweightsinner)*fFluxPhesInnerPixel*fFluxPhesInnerPixel;

    }

  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
    {
      fFluxPhesOuterPixel    = sumphesouter/sumweightsouter;
      fFluxPhesOuterPixelVar = (1./sumweightsouter)*fFluxPhesOuterPixel*fFluxPhesOuterPixel;
    }

  meanFluxPhotonsRelVar = fFluxPhesInnerPixelVar
                       / (fFluxPhesInnerPixel    * fFluxPhesInnerPixel);

  fFluxPhotonsInnerPixel    =  fFluxPhesInnerPixel/fAverageQE;
  fFluxPhotonsInnerPixelVar =  (meanFluxPhotonsRelVar + avQERelVar)
                                 * fFluxPhotonsInnerPixel * fFluxPhotonsInnerPixel;

  fFluxPhotonsOuterPixel    = 4. *fFluxPhotonsInnerPixel;
  fFluxPhotonsOuterPixelVar = 16.*fFluxPhotonsInnerPixelVar;  

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

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

  return kTRUE;
}

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

  const Float_t meanphotRelVar  = fFluxPhotonsInnerPixelVar
                               /( fFluxPhotonsInnerPixel    * fFluxPhotonsInnerPixel );

  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               = fFluxPhotonsInnerPixel / pix->GetMean(); 
      else
        conv               = fFluxPhotonsOuterPixel / pix->GetMean(); 

      if (conv <= 0.)
        {
          pix->SetFFactorMethodValid(kFALSE);
          continue;
        }
      
      const Float_t chargeRelVar       =   pix->GetMeanErr() * pix->GetMeanErr()
                                       / ( pix->GetMean()    * pix->GetMean());
      const Float_t rsigmaRelVar =   pix->GetRSigmaErr() *  pix->GetRSigmaErr()
                                        / (pix->GetRSigma()    * pix->GetRSigma()) ;

      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->GetRSigma()/pix->GetMean())
                                    *TMath::Sqrt(fFluxPhotonsInnerPixel);
      
      //
      // Calculate the error of the Total F-Factor of the camera ( in photons )
      //
      const Float_t totalFFactorVar = rsigmaRelVar + 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.GetFluxInsidePlexiglass();
  Float_t fluxerr = blindpix.GetFluxInsidePlexiglassErr();

  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->GetMean();
      const Float_t area      = geom[idx].GetA();
      const Float_t chargeerr = pix->GetMeanErr();          
      
      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.GetFluxOutsidePlexiglass();
  Float_t fluxerr = pindiode.GetFluxOutsidePlexiglassErr();

  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->GetMean();
      const Float_t area      = geom[idx].GetA();
      const Float_t chargeerr = pix->GetMeanErr();          
      
      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)
{

  MCalibrationChargePix &pix = (MCalibrationChargePix&)(*this)[ipx];

  mean  = pix.GetMeanConversionBlindPixelMethod();
  err   = pix.GetConversionBlindPixelMethodErr();
  sigma = pix.GetSigmaConversionBlindPixelMethod();

  return kTRUE;
}


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

  MCalibrationChargePix &pix = (MCalibrationChargePix&)(*this)[ipx];

  Float_t conv = pix.GetMeanConversionFFactorMethod();

  if (conv < 0.)
    return kFALSE;

  mean  = conv;
  err   = pix.GetConversionFFactorMethodErr();
  sigma = pix.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)
{

  MCalibrationChargePix &pix = (MCalibrationChargePix&)(*this)[ipx];

  mean  = pix.GetMeanConversionPINDiodeMethod();
  err   = pix.GetConversionPINDiodeMethodErr();
  sigma = pix.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();
}

*/
