source: trunk/MagicSoft/Mars/mcalib/MCalibrationChargeCalc.cc@ 4604

/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Markus Gaug  02/2004 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//   MCalibrationChargeCalc
//
//   Task to calculate the calibration conversion factors and quantum efficiencies 
//   from the fit results to the summed FADC slice distributions delivered by 
//   MCalibrationChargeCam, MCalibrationChargePix, MCalibrationChargeBlindPix and 
//   MCalibrationChargePINDiode, calculated and filled by MHCalibrationChargeCam, 
//   MHCalibrationChargePix, MHCalibrationChargeBlindPix and MHCalibrationChargePINDiode.
//
//   PreProcess(): Initialize pointers to MCalibrationChargeCam, MCalibrationChargeBlindPix
//                 MCalibrationChargePINDiode and  MCalibrationQECam
//
//                 Initialize pulser light wavelength
//               
//   ReInit():     MCalibrationCam::InitSize(NumPixels) is called from MGeomApply (which allocates
//                 memory in a TClonesArray of type MCalibrationChargePix)
//                 Initializes pointer to MBadPixelsCam
//
//   Process():    Nothing to be done, histograms getting filled by MHCalibrationChargeCam
//
//   PostProcess(): - FinalizePedestals()
//                  - FinalizeCharges()
//                  - FinalizeFFactorMethod()
//                  - FinalizeBadPixels()
//                  - FinalizeBlindPixel()
//                  - FinalizePINDiode()
//                  - FinalizeFFactorQECam()
//                  - FinalizeBlindPixelQECam()
//                  - FinalizePINDiodeQECam()
//
//  Input Containers:
//   MCalibrationChargeCam
//   MCalibrationChargeBlindPix
//   MCalibrationChargePINDiode
//   MCalibrationQECam
//   MPedestalCam
//   MBadPixelsCam
//   MGeomCam
//   MTime
//
//  Output Containers:
//   MCalibrationChargeCam
//   MCalibrationChargeBlindPix
//   MCalibrationChargePINDiode
//   MCalibrationQECam
//   MBadPixelsCam
//
//
//  Preliminary description of the calibration in photons (email from 12/02/04)
//
//  Why calibrating in photons:
//  ===========================
//  
//  At the Barcelona meeting in 2002, we decided to calibrate the camera in
//  photons. This for the following reasons:
//  
//  * The physical quantity arriving at the camera are photons. This is
//    the direct physical information from the air shower. The photons
//    have a flux and a spectrum.
//  
//  * The photon fluxes depend mostly on the shower energy (with
//    corrections deriving from the observation conditions), while the photon
//    spectra depend mostly on the observation conditions: zenith angle,
//    quality of the air, also the impact parameter of the shower.
//  
//  * The photomultiplier, in turn, has different response properties
//    (quantum efficiencies) for photons of different colour. (Moreover,
//    different pixels have slightly different quantum efficiencies).
//    The resulting number of photo-electrons is then amplified (linearly)
//    with respect to the photo-electron flux.
//  
//  * In the ideal case, one would like to disentagle the effects
//    of the observation conditions from the primary particle energy (which
//    one likes to measure). To do so, one needs:
//  
//    1) A reliable calibration relating the FADC counts to the photo-electron
//       flux -> This is accomplished with the F-Factor method.
//  
//    2) A reliable calibration of the wavelength-dependent quantum efficiency
//       -> This is accomplished with the combination of the three methods,
//          together with QE-measurements performed by David in order to do
//          the interpolation.
//  
//    3) A reliable calibration of the observation conditions. This means:
//       - Tracing the atmospheric conditions   -> LIDAR
//       - Tracing the observation zenith angle -> Drive System
//
//    4) Some knowlegde about the impact parameter:
//       - This is the only part which cannot be accomplished well with a
//         single telescope. We would thus need to convolute the spectrum
//         over the distribution of impact parameters.
//  
//  
//  How an ideal calibration would look like:
//  =========================================
//  
//  We know from the combined PIN-Diode and Blind-Pixel Method the response of
//  each pixel to well-measured light fluxes in three representative
//  wavelengths (green, blue, UV). We also know the response to these light
//  fluxes in photo-electrons. Thus, we can derive:
//  
//  - conversion factors to photo-electrons
//  - conversion factors to photons in three wavelengths.
//  
//  Together with David's measurements and some MC-simulation, we should be
//  able to derive tables for typical Cherenkov-photon spectra - convoluted
//  with the impact parameters and depending on the athmospheric conditions
//  and the zenith angle (the "outer parameters").
//  
//  From these tables we can create "calibration tables" containing some
//  effective quantum efficiency depending on these outer parameters and which
//  are different for each pixel.
//  
//  In an ideal MCalibrate, one would thus have to convert first the FADC
//  slices to Photo-electrons and then, depending on the outer parameters,
//  look up the effective quantum efficiency and get the mean number of
//  photons which is then used for the further analysis.
//  
//  How the (first) MAGIC calibration should look like:
//  ===================================================
//  
//  For the moment, we have only one reliable calibration method, although
//  with very large systematic errors. This is the F-Factor method. Knowing
//  that the light is uniform over the whole camera (which I would not at all
//  guarantee in the case of the CT1 pulser), one could in principle already
//  perform a relative calibration of the quantum efficiencies in the UV.
//  However, the spread in QE at UV is about 10-15% (according to the plot
//  that Abelardo sent around last time. The spread in photo-electrons is 15%
//  for the inner pixels, but much larger (40%) for the outer ones.
//  
//  I'm not sure if we can already say that we have measured the relative
//  difference in quantum efficiency for the inner pixels and produce a first
//  QE-table for each pixel. To so, I would rather check in other wavelengths
//  (which we can do in about one-two weeks when the optical transmission of
//  the calibration trigger is installed).
//  
//  Thus, for the moment being, I would join Thomas proposal to calibrate in
//  photo-electrons and apply one stupid average quantum efficiency for all
//  pixels. This keeping in mind that we will have much preciser information
//  in about one to two weeks.
//  
//  
//  What MCalibrate should calculate and what should be stored:
//  ===========================================================
//  
//  It is clear that in the end, MCerPhotEvt will store photons.
//  MCalibrationCam stores the conversionfactors to photo-electrons and also
//  some tables of how to apply the conversion to photons, given the outer
//  parameters. This is not yet implemented and not even discussed.
//  
//  To start, I would suggest that we define the "average quantum efficiency"
//  (maybe something like 25+-3%) and apply them equally to all
//  photo-electrons. Later, this average factor can be easily replaced by a
//  pixel-dependent factor and later by a (pixel-dependent) table.
//  
//  
//  
//////////////////////////////////////////////////////////////////////////////
#include "MCalibrationChargeCalc.h"

#include <TSystem.h>
#include <TH1.h>
#include <TF1.h>

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

#include "MParList.h"

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

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

#include "MPedestalCam.h"
#include "MPedestalPix.h"

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

#include "MExtractedSignalCam.h"
#include "MExtractedSignalPix.h"
#include "MExtractedSignalBlindPixel.h"
#include "MExtractedSignalPINDiode.h"

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

#include "MCalibrationQECam.h"
#include "MCalibrationQEPix.h"

#include "MCalibrationCam.h"

ClassImp(MCalibrationChargeCalc);

using namespace std;

const Float_t MCalibrationChargeCalc::fgChargeLimit            = 2.5;
const Float_t MCalibrationChargeCalc::fgChargeErrLimit         = 0.;
const Float_t MCalibrationChargeCalc::fgChargeRelErrLimit      = 1.;
const Float_t MCalibrationChargeCalc::fgLambdaErrLimit         = 0.2;
const Float_t MCalibrationChargeCalc::fgLambdaCheckLimit       = 0.5;
const Float_t MCalibrationChargeCalc::fgPheErrLimit            = 4.5;
const Float_t MCalibrationChargeCalc::fgFFactorErrLimit        = 4.5;
// --------------------------------------------------------------------------
//
// Default constructor. 
//
// Sets the pointer to fQECam and fGeom to NULL
// 
// Calls AddToBranchList for: 
// - MRawEvtData.fHiGainPixId   
// - MRawEvtData.fLoGainPixId   
// - MRawEvtData.fHiGainFadcSamples
// - MRawEvtData.fLoGainFadcSamples
//
// Initializes:
// - fChargeLimit       to fgChargeLimit 
// - fChargeErrLimit    to fgChargeErrLimit   
// - fChargeRelErrLimit to fgChargeRelErrLimit
// - fFFactorErrLimit   to fgFFactorErrLimit
// - fLambdaCheckLimit  to fgLambdaCheckLimit 
// - fLambdaErrLimit    to fgLambdaErrLimit 
// - fPheErrLimit       to fgPheErrLimit  
// - fPulserColor       to MCalibrationCam::kCT1
// - fOutputPath        to "."
// - fOutputFile        to "ChargeCalibStat.txt"
// - flag debug         to kFALSE
//
// Calls:
// - Clear()
//
MCalibrationChargeCalc::MCalibrationChargeCalc(const char *name, const char *title)
    : fQECam(NULL), fGeom(NULL)
{
        
  fName  = name  ? name  : "MCalibrationChargeCalc";
  fTitle = title ? title : "Task to calculate the calibration constants and MCalibrationCam ";
  
  AddToBranchList("MRawEvtData.fHiGainPixId");
  AddToBranchList("MRawEvtData.fLoGainPixId");
  AddToBranchList("MRawEvtData.fHiGainFadcSamples");
  AddToBranchList("MRawEvtData.fLoGainFadcSamples");
  
  SetChargeLimit       ();
  SetChargeErrLimit    ();  
  SetChargeRelErrLimit ();
  SetFFactorErrLimit   ();
  SetLambdaCheckLimit  ();
  SetLambdaErrLimit    ();
  SetPheErrLimit       ();
  SetOutputPath        ();
  SetOutputFile        ();
  SetDebug       ( kFALSE );

  Clear();
  
}

// --------------------------------------------------------------------------
//
// Sets:
// - all variables to 0., 
// - all flags to kFALSE
// - all pointers to NULL
// - the pulser colour to kNONE
// - fBlindPixelFlags to 0
// - fPINDiodeFlags   to 0
//
void MCalibrationChargeCalc::Clear(const Option_t *o)
{
  
    fNumHiGainSamples          = 0.;
    fNumLoGainSamples          = 0.;
    fSqrtHiGainSamples         = 0.;
    fSqrtLoGainSamples         = 0.;
    fNumInnerFFactorMethodUsed = 0;

    fBadPixels  = NULL;
    fCam        = NULL;
    fBlindPixel = NULL;
    fBlindCam   = NULL;
    fPINDiode   = NULL;
    fPedestals  = NULL;

    SetPulserColor ( MCalibrationCam::kNONE );
    
    fBlindPixelFlags.Set(0);
    fPINDiodeFlags  .Set(0);
    fResultFlags    .Set(0);
}


// -----------------------------------------------------------------------------------
//
// The following container are searched for and execution aborted if not in MParList:
//  - MPedestalCam
//
// The following containers are searched and created if they were not found:
//
//  - MCalibrationQECam
//  - MBadPixelsCam
//
Int_t MCalibrationChargeCalc::PreProcess(MParList *pList)
{
  
  // 
  // Containers that have to be there.
  //
  fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");
  if (!fPedestals)
    {
      *fLog << err << "MPedestalCam not found... aborting" << endl;
      return kFALSE;
    }
  
  // 
  // Containers that are created in case that they are not there.
  //
  fQECam = (MCalibrationQECam*)pList->FindCreateObj("MCalibrationQECam");
  if (!fQECam)
    {
      *fLog << err << "Cannot find nor create MCalibrationQECam... aborting" << endl;
      return kFALSE;
    }

  fBadPixels = (MBadPixelsCam*)pList->FindCreateObj("MBadPixelsCam");
  if (!fBadPixels)
    {
      *fLog << err << "Could not find or create MBadPixelsCam ... aborting." << endl;
      return kFALSE;
    }

  
  //
  // Check the pulser colour --> FIXME: this solution is only valid until the arrival of the DM's
  //
  if (fPulserColor == MCalibrationCam::kNONE)
    {
      *fLog << endl;
      *fLog << err << GetDescriptor()
            << ": No Pulser colour has been chosen. Since the installation of the IFAE pulser box,"
            << " you HAVE to provide the LEDs colour, otherwise there is no calibration. " << endl;
      *fLog << "See e.g. the macro calibration.C " << endl;
      return kFALSE;
    }
  
  return kTRUE;
}


// --------------------------------------------------------------------------
//
// Search for the following input containers and abort if not existing:
//  - MGeomCam
//  - MCalibrationChargeCam
// 
// Search for the following input containers and give a warning if not existing:
//  - MCalibrationChargeBlindPix
//  - MCalibrationChargePINDiode
//
// It retrieves the following variables from MCalibrationChargeCam:
//
//  - fNumHiGainSamples 
//  - fNumLoGainSamples 
//
// It defines the PixId of every pixel in:
//
// - MCalibrationChargeCam 
// - MCalibrationQECam
//
// It sets all pixels in excluded which have the flag fBadBixelsPix::IsBad() set in:
// 
// - MCalibrationChargePix
// - MCalibrationQEPix
//
// Sets the pulser colour and tests if it has not changed w.r.t. fPulserColor in:
//
// - MCalibrationChargeCam
// - MCalibrationChargeBlindPix (if existing)
// - MCalibrationChargePINDiode (if existing)
//
Bool_t MCalibrationChargeCalc::ReInit(MParList *pList )
{

  fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
  if (!fGeom)
    {
      *fLog << err << "No MGeomCam found... aborting." << endl;
      return kFALSE;
    }
  
  fCam = (MCalibrationChargeCam*)pList->FindObject("MCalibrationChargeCam");
  if (!fCam)
    {
      *fLog << err << "Cannot find MCalibrationChargeCam... aborting" << endl;
      *fLog << err << "Maybe you forget to call an MFillH for the MHCalibrationChargeCam before..." << endl;
      return kFALSE;
    }

  // 
  // Optional Containers 
  //
  fBlindPixel = (MCalibrationChargeBlindPix*)pList->FindObject("MCalibrationChargeBlindPix");
  if (!fBlindPixel)
    {
      fBlindCam = (MCalibrationChargeBlindCam*)pList->FindObject("MCalibrationChargeBlindCam");
      if (!fBlindCam)
        {
      *fLog << endl;
      *fLog << warn << GetDescriptor() 
            << ": MCalibrationChargeBlindPix nor MCalibrationChargeBlindCam "
            << " found... no Blind Pixel method! " << endl;
        }
    }
  
  fPINDiode = (MCalibrationChargePINDiode*)pList->FindObject("MCalibrationChargePINDiode");
  if (!fPINDiode)
    {
      *fLog << endl;
      *fLog << warn << GetDescriptor() 
            << ": MCalibrationChargePINDiode not found... no PIN Diode method! " << endl;
    }
  

  //
  // Initialize the pulser colours
  //
  if (fCam->GetPulserColor() == MCalibrationCam::kNONE)
    {
      fCam->SetPulserColor( fPulserColor );
      
      if (fBlindPixel)
        fBlindPixel->SetColor( fPulserColor );
      
      if (fBlindCam)
        fBlindCam->SetColor( fPulserColor );
      
      if (fPINDiode)
        fPINDiode->SetColor( fPulserColor );
    }

  if (fPulserColor != fCam->GetPulserColor())
    {
      *fLog << err << GetDescriptor() 
            << ": Pulser colour has changed w.r.t. last file in MCalibrationChargeCam" << endl;
      *fLog << err << "This feature is not yet implemented, sorry ... aborting " << endl;
      return kFALSE;
    }

  if (fBlindPixel)
    if (fPulserColor != fBlindPixel->GetColor())
      {
        *fLog << err << GetDescriptor() 
              << ": Pulser colour has changed w.r.t. last file in MCalibrationChargeBlindPix." << endl;
        *fLog << err << "This feature is not yet implemented, sorry ... aborting " << endl;
        return kFALSE;
      }

  if (fBlindCam)
    if (fPulserColor != fBlindCam->GetColor())
      {
        *fLog << err << GetDescriptor() 
              << ": Pulser colour has changed w.r.t. last file in MCalibrationChargeBlindCam." << endl;
        *fLog << err << "This feature is not yet implemented, sorry ... aborting " << endl;
        return kFALSE;
      }
  
  if (fPINDiode)
    if (fPulserColor != fPINDiode->GetColor())
      {
        *fLog << err << GetDescriptor() 
              << ": Pulser colour has changed w.r.t. last file in MCalibrationChargePINDiode." << endl;
        *fLog << err << "This feature is not yet implemented, sorry ... aborting " << endl;
        return kFALSE;
      }
  

  fNumHiGainSamples  =  fCam->GetNumHiGainFADCSlices();
  fNumLoGainSamples  =  fCam->GetNumLoGainFADCSlices();

  fSqrtHiGainSamples =  TMath::Sqrt(fNumHiGainSamples);
  fSqrtLoGainSamples =  TMath::Sqrt(fNumLoGainSamples);
  
  UInt_t npixels     = fGeom->GetNumPixels();
  
  for (UInt_t i=0; i<npixels; i++)
    {
      
      MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)  [i];
      MCalibrationQEPix     &pqe = (MCalibrationQEPix&)    (*fQECam)[i];
      MBadPixelsPix         &bad = (*fBadPixels)[i];
      
      pix.SetPixId(i);
      pqe.SetPixId(i);
      
      if (bad.IsBad())
        {
          pix.SetExcluded();
          pqe.SetExcluded();
          continue;
        }
      
      if (IsDebug())
        pix.SetDebug();
    }

  return kTRUE;
}

// ----------------------------------------------------------------------------------
//  
// Nothing to be done in Process, but have a look at MHCalibrationChargeCam, instead
// 
Int_t MCalibrationChargeCalc::Process()
{
  return kTRUE;
}

// -----------------------------------------------------------------------
//
// Return if number of executions is null.
//
// First loop over pixels, average areas and sectors, call: 
//  - FinalizePedestals() 
//  - FinalizeCharges()
// for every entry. Count number of valid pixels in loop and return kFALSE
// if there are none (the "Michele check").
//
// Call FinalizeBadPixels()
//
// Call FinalizeFFactorMethod() (second and third loop over pixels and areas)
//
// Call FinalizeBlindPixel() 
// Call FinalizePINDiode()
//
// Call  FinalizeFFactorQECam() (fourth loop over pixels and areas)
// Call  FinalizeBlindPixelQECam() (fifth loop over pixels and areas)
// Call  FinalizePINDiodeQECam() (sixth loop over pixels and areas)
//
// Call  FinalizeUnsuitablePixels()
//  
// Call MParContainer::SetReadyToSave() for fCam, fQECam, fBadPixels and 
//                                          fBlindPixel and fPINDiode if they exist
//
// Print out some statistics
//
Int_t MCalibrationChargeCalc::PostProcess()
{

  if (GetNumExecutions()==0)
    return kFALSE;

  *fLog << endl;

  if (fPINDiode)
    if (!fPINDiode->IsValid())
      {
        *fLog << warn << GetDescriptor() 
              << ": MCalibrationChargePINDiode is declared not valid... no PIN Diode method! " << endl;
        fPINDiode = NULL;
      }

  if (fBlindPixel)
    if (!fBlindPixel->IsValid())
      {
        *fLog << warn << GetDescriptor() 
              << ": MCalibrationChargeBlindPix is declared not valid... no Blind Pixel method! " << endl;
        fBlindPixel = NULL;
      }

  if (fBlindCam)
    if (!fBlindCam->IsValid())
      {
        *fLog << warn << GetDescriptor() 
              << ": MCalibrationChargeBlindCam is declared not valid... no Blind Pixel method! " << endl;
        fBlindCam = NULL;
      }
  
  *fLog << endl;
  //
  // First loop over pixels, call FinalizePedestals and FinalizeCharges
  //
  Int_t   nvalid      = 0;

  for (Int_t pixid=0; pixid<fPedestals->GetSize(); pixid++)
    {

      MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[pixid];
      //
      // Check if the pixel has been excluded from the fits
      //
      if (pix.IsExcluded())
        continue;

      MPedestalPix  &ped = (*fPedestals)[pixid];
      MBadPixelsPix &bad = (*fBadPixels)[pixid];
      const Int_t aidx    = (*fGeom)[pixid].GetAidx();

      FinalizePedestals(ped,pix,aidx);

      if (FinalizeCharges(pix,bad,"pixel  "))
        nvalid++;
    }

  *fLog << endl;  
  //
  // The Michele check ...
  //
  if (nvalid == 0)
  {
      *fLog << err << GetDescriptor() << ": All pixels have non-valid calibration. " 
            << "Did you forget to fill the histograms "
            << "(filling MHCalibrationChargeCam from MExtractedSignalCam using MFillH) ? " << endl;
      *fLog << err << GetDescriptor() << ": Or, maybe, you have used a pedestal run " 
            << "instead of a calibration run " << endl;
      return kFALSE;
  }

  for (UInt_t aidx=0; aidx<fGeom->GetNumAreas(); aidx++)
    {

      const MPedestalPix     &ped = fPedestals->GetAverageArea(aidx);
      MCalibrationChargePix &pix = (MCalibrationChargePix&)fCam->GetAverageArea(aidx);

      FinalizePedestals(ped,pix,aidx);
      FinalizeCharges(pix, fCam->GetAverageBadArea(aidx),"area id");
    }
  
  *fLog << endl;

  for (UInt_t sector=0; sector<fGeom->GetNumSectors(); sector++)
    {

      const MPedestalPix     &ped = fPedestals->GetAverageSector(sector);

      MCalibrationChargePix &pix = (MCalibrationChargePix&)fCam->GetAverageSector(sector);
      FinalizePedestals(ped,pix, 0);
    }
  
  *fLog << endl;

  //
  // Finalize Bad Pixels
  // 
  FinalizeBadPixels();

  // 
  // Finalize F-Factor method
  //
  if (!FinalizeFFactorMethod())
    {
      *fLog << warn << "Could not calculate the photons flux from the F-Factor method " << endl;
      fCam->SetFFactorMethodValid(kFALSE); 
      return kFALSE;
    }
  else
    fCam->SetFFactorMethodValid(kTRUE);
  
  *fLog << endl;
  // 
  // Finalize Blind Pixel
  //
  if (FinalizeBlindPixel())
    fQECam->SetBlindPixelMethodValid(kTRUE);
  else
    fQECam->SetBlindPixelMethodValid(kFALSE); 

  // 
  // Finalize PIN Diode
  //
  if (FinalizePINDiode())
    fQECam->SetPINDiodeMethodValid(kTRUE);
  else
    fQECam->SetPINDiodeMethodValid(kFALSE); 

  //
  // Finalize QE Cam
  //
  FinalizeFFactorQECam();
  FinalizeBlindPixelQECam();
  FinalizePINDiodeQECam();

  //
  // Re-direct the output to an ascii-file from now on:
  //
  MLog asciilog;
  asciilog.SetOutputFile(GetOutputFile(),kTRUE);
  SetLogStream(&asciilog);
  //
  // Finalize calibration statistics
  //
  FinalizeUnsuitablePixels();

  fCam      ->SetReadyToSave();
  fQECam    ->SetReadyToSave();
  fBadPixels->SetReadyToSave();

  if (fBlindPixel)
    fBlindPixel->SetReadyToSave();
  if (fBlindCam)
    fBlindCam->SetReadyToSave();
  if (fPINDiode)
    fPINDiode->SetReadyToSave();
  
  *fLog << inf << endl;
  *fLog << GetDescriptor() << ": Fatal errors statistics:" << endl;  

  PrintUncalibrated(MBadPixelsPix::kChargeIsPedestal,    
                    Form("%s%2.1f%s","Signal less than ",fChargeLimit," Pedestal RMS:                "));
  PrintUncalibrated(MBadPixelsPix::kChargeRelErrNotValid,
                    Form("%s%2.1f%s","Signal Error bigger than ",fChargeRelErrLimit," times Mean Signal:   "));
  PrintUncalibrated(MBadPixelsPix::kChargeSigmaNotValid, 
                    "Signal Sigma smaller than Pedestal RMS:           ");
  PrintUncalibrated(MBadPixelsPix::kLoGainSaturation,    
                    "Pixels with Low Gain Saturation:                  ");
  PrintUncalibrated(MBadPixelsPix::kMeanTimeInFirstBin, 
                    Form("%s%2.1f%s","Mean Abs. Arr. Time in First ",1.," Bin(s):          "));
  PrintUncalibrated(MBadPixelsPix::kMeanTimeInLast2Bins, 
                    Form("%s%2.1f%s","Mean Abs. Arr. Time in Last ",2.," Bin(s):           "));
  PrintUncalibrated(MBadPixelsPix::kDeviatingNumPhes,    
                    "Pixels with deviating number of phes:             ");
  PrintUncalibrated(MBadPixelsPix::kDeviatingFFactor,    
                    "Pixels with deviating F-Factor:                   ");

  *fLog << inf << endl;
  *fLog << GetDescriptor() << ": Unreliable errors statistics:" << endl;  

  PrintUncalibrated(MBadPixelsPix::kHiGainOscillating,   
                    "Pixels with changing Hi Gain signal over time:    ");
  PrintUncalibrated(MBadPixelsPix::kLoGainOscillating,   
                    "Pixels with changing Lo Gain signal over time:    ");
  PrintUncalibrated(MBadPixelsPix::kHiGainNotFitted,     
                    "Pixels with unsuccesful Gauss fit to the Hi Gain: ");
  PrintUncalibrated(MBadPixelsPix::kLoGainNotFitted,     
                    "Pixels with unsuccesful Gauss fit to the Lo Gain: ");

  SetLogStream(&gLog);

  return kTRUE;
}

// ----------------------------------------------------------------------------------
//  
// Retrieves pedestal and pedestal RMS from MPedestalPix 
// Retrieves total entries from MPedestalCam
// Sets pedestal*fNumHiGainSamples and pedestal*fNumLoGainSamples in MCalibrationChargePix
// Sets pedRMS *fSqrtHiGainSamples and pedRMS *fSqrtLoGainSamples in MCalibrationChargePix
//
// If the flag MCalibrationPix::IsHiGainSaturation() is set, call also: 
// - MCalibrationChargePix::CalcLoGainPedestal()
//
void MCalibrationChargeCalc::FinalizePedestals(const MPedestalPix &ped, MCalibrationChargePix &cal, const Int_t aidx)
{
  
  //
  // get the pedestals
  //
  const Float_t pedes  = ped.GetPedestal();
  const Float_t prms   = ped.GetPedestalRms();
  const Float_t num    = TMath::Sqrt((Float_t)fPedestals->GetTotalEntries());


  //
  // set them in the calibration camera
  //
  if (cal.IsHiGainSaturation())
    {
      cal.SetPedestal(pedes* fNumLoGainSamples,
                      prms * fSqrtLoGainSamples,
                      prms * fNumLoGainSamples / num);
      cal.CalcLoGainPedestal((Float_t)fNumLoGainSamples, aidx);
    }
  else
    {
      cal.SetPedestal(pedes* fNumHiGainSamples,
                      prms * fSqrtHiGainSamples,
                      prms * fNumHiGainSamples / num);
    }
  
}

// ----------------------------------------------------------------------------------------------------
//
// Check fit results validity. Bad Pixels flags are set if:
//
// 1) Pixel has a mean smaller than fChargeLimit*PedRMS    ( Flag: MBadPixelsPix::kChargeIsPedestal)
// 2) Pixel has a mean error smaller than fChargeErrLimit  ( Flag: MBadPixelsPix::kChargeErrNotValid)
// 3) Pixel has mean smaller than fChargeRelVarLimit times its mean error 
//                                                         ( Flag: MBadPixelsPix::kChargeRelErrNotValid)
// 4) Pixel has a sigma bigger than its Pedestal RMS       ( Flag: MBadPixelsPix::kChargeSigmaNotValid )
// 
// Further returns if flags: MBadPixelsPix::kUnsuitableRun is set
// 
// Calls MCalibrationChargePix::CalcReducedSigma() and sets flag:  MBadPixelsPix::kChargeIsPedestal 
//       if not succesful.
// 
// Calls MCalibrationChargePix::CalcFFactorMethod() and sets flag: MBadPixelsPix::kDeviatingNumPhes) 
//       if not succesful.
//
Bool_t MCalibrationChargeCalc::FinalizeCharges(MCalibrationChargePix &cal, MBadPixelsPix &bad, const char* what)
{

  if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
    return kFALSE;

  if (cal.GetMean() < fChargeLimit*cal.GetPedRms())
    {
      *fLog << warn << GetDescriptor() 
            << Form(": Fitted Charge: %5.2f is smaller than %2.1f",cal.GetMean(),fChargeLimit)
            << Form("   Pedestal RMS: %5.2f in %s%4i",cal.GetPedRms(),what,cal.GetPixId()) << endl;
      bad.SetUncalibrated( MBadPixelsPix::kChargeIsPedestal);
    }
  
   if (cal.GetMean() < fChargeRelErrLimit*cal.GetMeanErr()) 
    {
      *fLog << warn << GetDescriptor() 
            << Form(": Fitted Charge: %4.2f is smaller than %2.1f",cal.GetMean(),fChargeRelErrLimit) 
            << Form(" times its error: %4.2f in %s%4i",cal.GetMeanErr(),what,cal.GetPixId()) << endl;
      bad.SetUncalibrated( MBadPixelsPix::kChargeRelErrNotValid );
    }

  if (cal.GetSigma() < cal.GetPedRms())
    {
      *fLog << warn << GetDescriptor() 
            << Form(": Sigma of Fitted Charge: %6.2f is smaller than",cal.GetSigma())
            << Form(" Ped. RMS: %5.2f in %s%4i",cal.GetPedRms(),what,cal.GetPixId()) << endl;
      bad.SetUncalibrated( MBadPixelsPix::kChargeSigmaNotValid );
    }

  if (!cal.CalcReducedSigma())
    {
      *fLog << warn << GetDescriptor() 
            << Form(": Could not calculate the reduced sigma in %s:               ",what)
            << Form("         %4i",cal.GetPixId()) 
            << endl;
      bad.SetUncalibrated(MBadPixelsPix::kChargeIsPedestal);
      return kFALSE;
    }

  if (!cal.CalcFFactorMethod())
    {
      *fLog << warn << GetDescriptor() 
            << Form(": Could not calculate the F-Factor in %s:                    ",what)
            << Form("         %4i",cal.GetPixId()) 
            << endl;
      bad.SetUncalibrated(MBadPixelsPix::kDeviatingNumPhes);
      return kFALSE;
    }

  return kTRUE;
}



// -----------------------------------------------------------------------------------
//
// Sets pixel to MBadPixelsPix::kUnsuitableRun, if one of the following flags is set:
// - MBadPixelsPix::kChargeIsPedestal
// - MBadPixelsPix::kChargeRelErrNotValid 
// - MBadPixelsPix::kChargeSigmaNotValid 
// - MBadPixelsPix::kMeanTimeInFirstBin 
// - MBadPixelsPix::kMeanTimeInLast2Bins 
// - MBadPixelsPix::kDeviatingNumPhes
//
// - Call MCalibrationPix::SetExcluded() for the bad pixels
//
void MCalibrationChargeCalc::FinalizeBadPixels()
{
  
  for (Int_t i=0; i<fBadPixels->GetSize(); i++)
    {
      
      MBadPixelsPix    &bad    = (*fBadPixels)[i];
      MCalibrationPix  &pix    = (*fCam)[i];

      if (bad.IsUncalibrated( MBadPixelsPix::kChargeIsPedestal))
        bad.SetUnsuitable(   MBadPixelsPix::kUnsuitableRun   );
 
      if (bad.IsUncalibrated( MBadPixelsPix::kChargeErrNotValid ))
        bad.SetUnsuitable(   MBadPixelsPix::kUnsuitableRun    );

      if (bad.IsUncalibrated( MBadPixelsPix::kChargeRelErrNotValid ))
        bad.SetUnsuitable(   MBadPixelsPix::kUnsuitableRun    );
 
      if (bad.IsUncalibrated( MBadPixelsPix::kChargeSigmaNotValid ))
        bad.SetUnsuitable(   MBadPixelsPix::kUnsuitableRun    );

      if (bad.IsUncalibrated( MBadPixelsPix::kMeanTimeInFirstBin ))
        bad.SetUnsuitable(   MBadPixelsPix::kUnsuitableRun    );

      if (bad.IsUncalibrated( MBadPixelsPix::kMeanTimeInLast2Bins ))
        bad.SetUnsuitable(   MBadPixelsPix::kUnsuitableRun    );

      if (bad.IsUncalibrated( MBadPixelsPix::kDeviatingNumPhes ))
        bad.SetUnsuitable(   MBadPixelsPix::kUnsuitableRun    );

      if (bad.IsUnsuitable(   MBadPixelsPix::kUnsuitableRun    ))
        pix.SetExcluded();

    }
}

// ------------------------------------------------------------------------
//
//
// First loop: Calculate a mean and mean RMS of photo-electrons per area index 
//             Include only pixels which are not MBadPixelsPix::kUnsuitableRun or 
//             MBadPixelsPix::kUnreliableRun (see FinalizeBadPixels()) and set 
//             MCalibrationChargePix::SetFFactorMethodValid(kFALSE) in that case.
//              
// Second loop: Get weighted mean number of photo-electrons and its RMS including 
//              only pixels with flag MCalibrationChargePix::IsFFactorMethodValid() 
//              and further exclude those deviating by more than fPheErrLimit mean 
//              sigmas from the mean (obtained in first loop). Set 
//              MBadPixelsPix::kDeviatingNumPhes if excluded.
// 
//              Set weighted mean and variance of photo-electrons per area index in: 
//              average area pixels of MCalibrationChargeCam (obtained from: 
//              MCalibrationChargeCam::GetAverageArea() )
// 
//              Set weighted mean  and variance of photo-electrons per sector in:
//              average sector pixels of MCalibrationChargeCam (obtained from: 
//              MCalibrationChargeCam::GetAverageSector() )
// 
Bool_t MCalibrationChargeCalc::FinalizeFFactorMethod()
{

  const UInt_t npixels  = fGeom->GetNumPixels();
  const UInt_t nareas   = fGeom->GetNumAreas();
  const UInt_t nsectors = fGeom->GetNumSectors();

  Float_t lowlim      [nareas];
  Float_t upplim      [nareas];
  Double_t areavars   [nareas];
  Double_t areaweights[nareas], sectorweights [nsectors];
  Double_t areaphes   [nareas], sectorphes    [nsectors];
  Int_t   numareavalid[nareas], numsectorvalid[nsectors];

  memset(lowlim        ,0, nareas   * sizeof(Float_t));
  memset(upplim        ,0, nareas   * sizeof(Float_t));
  memset(areaphes      ,0, nareas   * sizeof(Double_t));
  memset(areavars      ,0, nareas   * sizeof(Double_t));
  memset(areaweights   ,0, nareas   * sizeof(Double_t));
  memset(numareavalid  ,0, nareas   * sizeof(Int_t  ));
  memset(sectorweights ,0, nsectors * sizeof(Double_t));
  memset(sectorphes    ,0, nsectors * sizeof(Double_t));
  memset(numsectorvalid,0, nsectors * sizeof(Int_t  ));
  
  //
  // First loop: Get mean number of photo-electrons and the RMS
  //             The loop is only to recognize later pixels with very deviating numbers
  //
  MHCamera camphes(*fGeom,"Camphes","Phes in Camera");

  for (UInt_t i=0; i<npixels; i++)
    {
      
      MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)  [i];
      MBadPixelsPix         &bad = (*fBadPixels)[i];
      
      if (!pix.IsFFactorMethodValid())
        continue;

      if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
        {
          pix.SetFFactorMethodValid(kFALSE);
          continue;
        }

      // if (bad.IsUnsuitable(MBadPixelsPix::kUnreliableRun))
      //        continue;
      
      const Float_t nphe  = pix.GetPheFFactorMethod();
      const Int_t   aidx  = (*fGeom)[i].GetAidx();

      camphes.Fill(i,nphe);
      camphes.SetUsed(i);

      areaphes    [aidx] += nphe;
      areavars    [aidx] += nphe*nphe;
      numareavalid[aidx] ++;
    } 

  for (UInt_t i=0; i<nareas; i++)
    {
      if (numareavalid[i] == 0)
        {
          *fLog << warn << GetDescriptor() << ": No pixels with valid number of photo-electrons found "
                << "in area index: " << i << endl;
          continue;
        }

      if (numareavalid[i] == 1)
        areavars[i] = 0.;
      else if (numareavalid[i] == 0)
        {
          areaphes[i]    = -1.;
          areaweights[i] = -1.;
        }
      else
        {
          areavars[i] = (areavars[i] - areaphes[i]*areaphes[i]/numareavalid[i]) / (numareavalid[i]-1);
          areaphes[i] = areaphes[i] / numareavalid[i];
        }
      
      if (areavars[i] < 0.)
        {
          *fLog << warn << GetDescriptor() << ": No pixels with valid variance of photo-electrons found "
                << "in area index: " << i << endl;
          continue;
        }
      
      lowlim  [i] = areaphes[i] - fPheErrLimit*TMath::Sqrt(areavars[i]);
      upplim  [i] = areaphes[i] + fPheErrLimit*TMath::Sqrt(areavars[i]);

      TArrayI area(1);
      area[0] = i;

      TH1D *hist = camphes.ProjectionS(TArrayI(),area,"_py",100);
      hist->Fit("gaus","Q");
      const Float_t mean  = hist->GetFunction("gaus")->GetParameter(1);
      const Float_t sigma = hist->GetFunction("gaus")->GetParameter(2);
      const Int_t   ndf   = hist->GetFunction("gaus")->GetNDF();

      if (IsDebug())
        camphes.DrawClone();

      if (ndf < 2)
        {
          *fLog << warn << GetDescriptor() << ": Cannot use a Gauss fit to the number of photo-electrons " 
                << "in the camera with area index: " << i << endl;
          *fLog << warn << GetDescriptor() << ": Number of dof.: " << ndf << " is smaller than 2 " << endl;
          *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
          delete hist;
          continue;
        }
      
      const Double_t prob = hist->GetFunction("gaus")->GetProb();

      if (prob < 0.001)
        {
          *fLog << warn << GetDescriptor() << ": Cannot use a Gauss fit to the number of photo-electrons " 
                << "in the camera with area index: " << i << endl;
          *fLog << warn << GetDescriptor() << ": Fit probability " << prob 
                << " is smaller than 0.001 " << endl;
          *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
          delete hist;
          continue;
        }
      
      *fLog << inf << GetDescriptor() << ": Mean number of photo-electrons " 
            << "with area idx " << i << ": " 
            << Form("%7.2f+-%6.2f",mean,sigma) << endl;

      lowlim  [i] = mean  - fPheErrLimit*sigma;
      upplim  [i] = mean  + fPheErrLimit*sigma;

      delete hist;
    }

  *fLog << endl;

  memset(numareavalid,0,nareas*sizeof(Int_t));
  memset(areaphes    ,0,nareas*sizeof(Double_t));
  memset(areavars    ,0,nareas*sizeof(Double_t));

  //
  // Second loop: Get mean number of photo-electrons and its RMS excluding 
  //              pixels deviating by more than fPheErrLimit sigma. 
  //              Set the conversion factor FADC counts to photo-electrons
  // 
  for (UInt_t i=0; i<npixels; i++)
    {
      
      MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];

      if (!pix.IsFFactorMethodValid())
        continue;

      const Float_t nvar  = pix.GetPheFFactorMethodVar();

      if (nvar <= 0.)
        {
          pix.SetFFactorMethodValid(kFALSE);
          continue;
        }
      
      MBadPixelsPix         &bad = (*fBadPixels)[i];

      const Int_t   aidx   = (*fGeom)[i].GetAidx();
      const Int_t   sector = (*fGeom)[i].GetSector();
      const Float_t area   = (*fGeom)[i].GetA();
      const Float_t nphe   = pix.GetPheFFactorMethod();

      if ( nphe < lowlim[aidx] || nphe > upplim[aidx] )
        {
          *fLog << warn << GetDescriptor() << ": Number of phes: " 
                << Form("%7.2f out of %3.1f sigma limit: ",nphe,fPheErrLimit)
                << Form("[%7.2f,%7.2f] pixel%4i",lowlim[aidx],upplim[aidx],i) << endl;
          bad.SetUncalibrated( MBadPixelsPix::kDeviatingNumPhes );
          bad.SetUnsuitable  ( MBadPixelsPix::kUnsuitableRun    );
          pix.SetFFactorMethodValid(kFALSE);
          continue;
        }
      
      areaweights   [aidx]   += nphe*nphe;
      areaphes      [aidx]   += nphe;
      numareavalid  [aidx]   ++;

      if (aidx == 0)
        fNumInnerFFactorMethodUsed++;

      sectorweights [sector] += nphe*nphe/area/area;
      sectorphes    [sector] += nphe/area;
      numsectorvalid[sector] ++;
    } 

  *fLog << endl;

  for (UInt_t aidx=0; aidx<nareas; aidx++)
    {

      MCalibrationChargePix &apix = (MCalibrationChargePix&)fCam->GetAverageArea(aidx);

      if (numareavalid[aidx] == 1)
        areaweights[aidx] = 0.;
      else if (numareavalid[aidx] == 0)
        {
          areaphes[aidx]    = -1.;
          areaweights[aidx] = -1.;
        }
      else
        {
          areaweights[aidx] = (areaweights[aidx] - areaphes[aidx]*areaphes[aidx]/numareavalid[aidx])
                            / (numareavalid[aidx]-1);
          areaphes[aidx]    /=  numareavalid[aidx];
        }
      
      if (areaweights[aidx] < 0. || areaphes[aidx] <= 0.)
        {
          *fLog << warn << GetDescriptor() 
                << ": Mean number phes from area index " << aidx << " could not be calculated: "
                << " Mean: "  << areaphes[aidx] 
                << " Variance: " << areaweights[aidx]    << endl;
          apix.SetFFactorMethodValid(kFALSE);
          continue;
        }
      
      *fLog << inf << "Average total number phes in area idx " << aidx << ": " 
            << Form("%7.2f%s%6.2f",areaphes[aidx]," +- ",TMath::Sqrt(areaweights[aidx])) << endl;
                    
      apix.SetPheFFactorMethod   ( areaphes[aidx] );
      apix.SetPheFFactorMethodVar( areaweights[aidx] / numareavalid[aidx] );      
      apix.SetFFactorMethodValid ( kTRUE );

    }

  *fLog << endl;

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

      if (numsectorvalid[sector] == 1)
        sectorweights[sector] = 0.;
      else if (numsectorvalid[sector] == 0)
        {
          sectorphes[sector]    = -1.;
          sectorweights[sector] = -1.;
        }
      else
        {
          sectorweights[sector] = (sectorweights[sector] 
                                   - sectorphes[sector]*sectorphes[sector]/numsectorvalid[sector]
                                  ) 
                                / (numsectorvalid[sector]-1.);
          sectorphes[sector]     /=  numsectorvalid[sector];
        }
      
      MCalibrationChargePix &spix = (MCalibrationChargePix&)fCam->GetAverageSector(sector);

      if (sectorweights[sector] < 0. || sectorphes[sector] <= 0.)
        {
          *fLog << warn << GetDescriptor() 
                <<": Mean number phes per area for sector " << sector << " could not be calculated: "
                << " Mean: "  << sectorphes[sector] 
                << " Variance: "             << sectorweights[sector]    << endl;
          spix.SetFFactorMethodValid(kFALSE);
          continue;
        }

      *fLog << inf << "Average number phes per area in sector " << sector << ": " 
            << Form("%5.2f+-%4.2f  [phe/mm^2]",sectorphes[sector],TMath::Sqrt(sectorweights[sector])) 
            << endl;

      spix.SetPheFFactorMethod   ( sectorphes[sector] );
      spix.SetPheFFactorMethodVar( sectorweights[sector] / numsectorvalid[sector]);      
      spix.SetFFactorMethodValid ( kTRUE );

    }

  return kTRUE;
}


// ------------------------------------------------------------------------
//
// Returns kFALSE if pointer to MCalibrationChargeBlindPix is NULL
//
// The check returns kFALSE if:
//
// 1) fLambda and fLambdaCheck are separated relatively to each other by more than fLambdaCheckLimit
// 2) BlindPixel has an fLambdaErr greater than fLambdaErrLimit
// 
// Calls:
// - MCalibrationChargeBlindPix::CalcFluxInsidePlexiglass()
//
Bool_t MCalibrationChargeCalc::FinalizeBlindPixel()
{

  if (!fBlindPixel)
    return kFALSE;  

  const Float_t lambda      = fBlindPixel->GetLambda();
  const Float_t lambdaerr   = fBlindPixel->GetLambdaErr();
  const Float_t lambdacheck = fBlindPixel->GetLambdaCheck();

  if (2.*(lambdacheck-lambda)/(lambdacheck+lambda) > fLambdaCheckLimit)
    {
      *fLog << warn << GetDescriptor() 
            << Form("%s%4.2f%s%4.2f%s%4.2f%s",": Lambda: ",lambda," and Lambda-Check: ",
                    lambdacheck," differ by more than ",fLambdaCheckLimit," in the Blind Pixel ")
            << endl;
      return kFALSE;
    }
  
  if (lambdaerr > fLambdaErrLimit) 
    {
      *fLog << warn << GetDescriptor() 
            << Form("%s%4.2f%s%4.2f%s",": Error of Fitted Lambda: ",lambdaerr," is greater than ",
                    fLambdaErrLimit," in Blind Pixel ") << endl;
      return kFALSE;
    }
      
  if (!fBlindPixel->CalcFluxInsidePlexiglass())
    {
      *fLog << warn << "Could not calculate the flux of photons from the Blind Pixel, "
            << "will skip Blind Pixel Calibration " << endl;
      return kFALSE;
    }
  
  return kTRUE;
}

// ------------------------------------------------------------------------
//
// Returns kFALSE if pointer to MCalibrationChargePINDiode is NULL
//
// The check returns kFALSE if:
//
// 1) PINDiode has a fitted charge smaller than fChargeLimit*PedRMS
// 2) PINDiode has a fit error smaller than fChargeErrLimit
// 3) PINDiode has a fitted charge smaller its fChargeRelErrLimit times its charge error
// 4) PINDiode has a charge sigma smaller than its Pedestal RMS
// 
// Calls:
// - MCalibrationChargePINDiode::CalcFluxOutsidePlexiglass()
//
Bool_t MCalibrationChargeCalc::FinalizePINDiode()
{

  if (!fPINDiode)
    return kFALSE;  

  if (fPINDiode->GetMean() < fChargeLimit*fPINDiode->GetPedRms())
    {
      *fLog << warn << GetDescriptor() << ": Fitted Charge is smaller than "
            << fChargeLimit << " Pedestal RMS in PINDiode " << endl;
      return kFALSE;
    }
  
  if (fPINDiode->GetMeanErr() < fChargeErrLimit) 
    {
      *fLog << warn << GetDescriptor() << ": Error of Fitted Charge is smaller than "
            << fChargeErrLimit << " in PINDiode " << endl;
      return kFALSE;
    }
      
  if (fPINDiode->GetMean() < fChargeRelErrLimit*fPINDiode->GetMeanErr()) 
    {
      *fLog << warn << GetDescriptor() << ": Fitted Charge is smaller than "
            << fChargeRelErrLimit << "* its error in PINDiode " << endl;
      return kFALSE;
    }
      
  if (fPINDiode->GetSigma() < fPINDiode->GetPedRms())
    {
      *fLog << warn << GetDescriptor() 
            << ": Sigma of Fitted Charge smaller than Pedestal RMS in PINDiode " << endl;
      return kFALSE;
    }


  if (!fPINDiode->CalcFluxOutsidePlexiglass())
    {
      *fLog << warn << "Could not calculate the flux of photons from the PIN Diode, "
            << "will skip PIN Diode Calibration " << endl;
      return kFALSE;
    }
  
  return kTRUE;
}

// ------------------------------------------------------------------------
//
// Calculate the average number of photons outside the plexiglass with the 
// formula: 
// 
// av.Num.photons(area index) = av.Num.Phes(area index) 
//                            / MCalibrationQEPix::GetDefaultQE(fPulserColor) 
//                            / MCalibrationQEPix::GetPMTCollectionEff()
//                            / MCalibrationQEPix::GetLightGuidesEff(fPulserColor)
//                            / MCalibrationQECam::GetPlexiglassQE()
//
// Calculate the variance on the average number of photons assuming that the error on the 
// Quantum efficiency is reduced by the number of used inner pixels, but the rest of the 
// values keeps it ordinary error since it is systematic.
//
// Loop over pixels: 
//
// - Continue, if not MCalibrationChargePix::IsFFactorMethodValid() and set:
//                    MCalibrationQEPix::SetFFactorMethodValid(kFALSE,fPulserColor)
//
// - Call MCalibrationChargePix::CalcMeanFFactor(av.Num.photons) and set: 
//        MCalibrationQEPix::SetFFactorMethodValid(kFALSE,fPulserColor) if not succesful
//
// - Calculate the quantum efficiency with the formula:
// 
//   QE = ( Num.Phes / av.Num.photons ) * MGeomCam::GetPixRatio()
//
// - Set QE in MCalibrationQEPix::SetQEFFactor ( QE, fPulserColor );
//
// - Set Variance of QE in  MCalibrationQEPix::SetQEFFactorVar ( Variance, fPulserColor );
// - Set bit MCalibrationQEPix::SetFFactorMethodValid(kTRUE,fPulserColor) 
//
// - Call MCalibrationQEPix::UpdateFFactorMethod()
//
void MCalibrationChargeCalc::FinalizeFFactorQECam()
{

  if (fNumInnerFFactorMethodUsed < 2)
    {
      *fLog << warn << GetDescriptor() 
            << ": Could not calculate F-Factor Method: Less than 2 inner pixels valid! " << endl;
      return;
    }
  
  MCalibrationChargePix &avpix = (MCalibrationChargePix&)fCam->GetAverageArea(0);
  MCalibrationQEPix     &qepix = (MCalibrationQEPix&)  fQECam->GetAverageArea(0);  

  const Float_t avphotons   = avpix.GetPheFFactorMethod()       
                           / qepix.GetDefaultQE(fPulserColor) 
                           / qepix.GetPMTCollectionEff()
                           / qepix.GetLightGuidesEff(fPulserColor)
                           / fQECam->GetPlexiglassQE();

  const Float_t avphotrelvar = avpix.GetPheFFactorMethodRelVar() 
                            + qepix.GetDefaultQERelVar(fPulserColor) / fNumInnerFFactorMethodUsed
                            + qepix.GetPMTCollectionEffRelVar()
                            + qepix.GetLightGuidesEffRelVar(fPulserColor)
                            + fQECam->GetPlexiglassQERelVar();

  const UInt_t nareas   = fGeom->GetNumAreas();

  //
  // Set the results in the MCalibrationChargeCam
  //
  fCam->SetNumPhotonsFFactorMethod   (avphotons);
  if (avphotrelvar > 0.)
    fCam->SetNumPhotonsFFactorMethodErr(TMath::Sqrt( avphotrelvar * avphotons * avphotons));  

  Float_t lowlim           [nareas];
  Float_t upplim           [nareas];
  Double_t avffactorphotons [nareas];
  Double_t avffactorphotvar [nareas];
  Int_t   numffactor       [nareas];

  memset(lowlim          ,0, nareas   * sizeof(Float_t));
  memset(upplim          ,0, nareas   * sizeof(Float_t));
  memset(avffactorphotons,0, nareas   * sizeof(Double_t));
  memset(avffactorphotvar,0, nareas   * sizeof(Double_t));
  memset(numffactor      ,0, nareas   * sizeof(Int_t));

  const UInt_t npixels  = fGeom->GetNumPixels();

  MHCamera camffactor(*fGeom,"Camffactor","F-Factor in Camera");

  for (UInt_t i=0; i<npixels; i++)
    {
      
      MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
      MCalibrationQEPix   &qepix = (MCalibrationQEPix&)  (*fQECam)[i];
      MBadPixelsPix         &bad =                   (*fBadPixels)[i];

      if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
        continue;

      const Float_t photons = avphotons / fGeom->GetPixRatio(i);
      const Float_t qe      = pix.GetPheFFactorMethod() / photons ;

      if (!pix.CalcMeanFFactor( photons , avphotrelvar ))
        {
          (*fBadPixels)[i].SetUncalibrated( MBadPixelsPix::kDeviatingNumPhes );
          continue;
        }

      const Float_t qerelvar = avphotrelvar +  pix.GetPheFFactorMethodRelVar();

      qepix.SetQEFFactor    ( qe            , fPulserColor );
      qepix.SetQEFFactorVar ( qerelvar*qe*qe, fPulserColor );      
      qepix.SetFFactorMethodValid(  kTRUE   , fPulserColor );

      if (!qepix.UpdateFFactorMethod())
        *fLog << warn << GetDescriptor() 
              << ": Cannot update Quantum efficiencies with the F-Factor Method" << endl;

      const Int_t aidx = (*fGeom)[i].GetAidx();
      const Float_t ffactor = pix.GetMeanFFactorFADC2Phot();

      camffactor.Fill(i,ffactor);
      camffactor.SetUsed(i);

      avffactorphotons[aidx] += ffactor;
      avffactorphotvar[aidx] += ffactor*ffactor;
      numffactor[aidx]++;
    }

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

      if (numffactor[i] == 0)
        {
          *fLog << warn << GetDescriptor() << ": No pixels with valid total F-Factor found "
                << "in area index: " << i << endl;
          continue;
        }

      avffactorphotvar[i] = (avffactorphotvar[i] - avffactorphotons[i]*avffactorphotons[i]/numffactor[i]) / (numffactor[i]-1.);
      avffactorphotons[i] = avffactorphotons[i] / numffactor[i];

      if (avffactorphotvar[i] < 0.)
        {
          *fLog << warn << GetDescriptor() << ": No pixels with valid variance of total F-Factor found "
                << "in area index: " << i << endl;
          continue;
        }

      lowlim  [i] = 1.1;   // Lowest known F-Factor of a PMT
      upplim  [i] = avffactorphotons[i] + fFFactorErrLimit*TMath::Sqrt(avffactorphotvar[i]);

      TArrayI area(1);
      area[0] = i;

      TH1D *hist = camffactor.ProjectionS(TArrayI(),area,"_py",100);
      hist->Fit("gaus","Q");
      const Float_t mean  = hist->GetFunction("gaus")->GetParameter(1);
      const Float_t sigma = hist->GetFunction("gaus")->GetParameter(2);
      const Int_t   ndf   = hist->GetFunction("gaus")->GetNDF();

      if (IsDebug())
        camffactor.DrawClone();

      if (ndf < 2)
        {
          *fLog << warn << GetDescriptor() << ": Cannot use a Gauss fit to the F-Factor " 
                << "in the camera with area index: " << i << endl;
          *fLog << warn << GetDescriptor() << ": Number of dof.: " << ndf << " is smaller than 2 " << endl;
          *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
          delete hist;
          continue;
        }
      
      const Double_t prob = hist->GetFunction("gaus")->GetProb();

      if (prob < 0.001)
        {
          *fLog << warn << GetDescriptor() << ": Cannot use a Gauss fit to the F-Factor " 
                << "in the camera with area index: " << i << endl;
          *fLog << warn << GetDescriptor() << ": Fit probability " << prob 
                << " is smaller than 0.001 " << endl;
          *fLog << warn << GetDescriptor() << ": Will use the simple mean and rms " << endl;
          delete hist;
          continue;
        }

      *fLog << inf << GetDescriptor() << ": Mean overall F-Factor (photons to FADC counts) " 
            << "with area index: " << i << ": " 
            << Form("%4.2f+-%4.2f",mean,sigma) << endl;

      lowlim  [i] = 1.1;
      upplim  [i] = mean  + fFFactorErrLimit*sigma;

      delete hist;
    }
  
  *fLog << endl;
  
  for (UInt_t i=0; i<npixels; i++)
    {
      
      MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
      MBadPixelsPix         &bad =                   (*fBadPixels)[i];

      if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
        continue;
      
      const Float_t ffactor = pix.GetMeanFFactorFADC2Phot();
      const Int_t   aidx   = (*fGeom)[i].GetAidx();

      if ( ffactor < lowlim[aidx] || ffactor > upplim[aidx] )
        {
          *fLog << warn << GetDescriptor() << ": Overall F-Factor: " 
                << Form("%5.2f out of %3.1f sigma limit: ",ffactor,fFFactorErrLimit)
                << Form("[%5.2f,%5.2f] pixel %4i",lowlim[aidx],upplim[aidx],i) << endl;
          bad.SetUncalibrated( MBadPixelsPix::kDeviatingFFactor );
          bad.SetUnsuitable  ( MBadPixelsPix::kUnsuitableRun    );
        }
    }

  for (UInt_t i=0; i<npixels; i++)
    {
      
      MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
      MCalibrationQEPix   &qepix = (MCalibrationQEPix&)  (*fQECam)[i];
      MBadPixelsPix         &bad =                   (*fBadPixels)[i];
      
      if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
        {
          qepix.SetFFactorMethodValid(kFALSE,fPulserColor);
          pix.SetFFactorMethodValid(kFALSE);
          pix.SetExcluded();
          continue;
        }
    }
}


// ------------------------------------------------------------------------
//
// Loop over pixels: 
//
// - Continue, if not MCalibrationChargeBlindPix::IsFluxInsidePlexiglassAvailable() and set:
//                    MCalibrationQEPix::SetBlindPixelMethodValid(kFALSE,fPulserColor)
//
// - Calculate the quantum efficiency with the formula:
// 
//   QE =  Num.Phes / MCalibrationChargeBlindPix::GetFluxInsidePlexiglass() 
//        / MGeomPix::GetA() * MCalibrationQECam::GetPlexiglassQE()
//
// - Set QE in MCalibrationQEPix::SetQEBlindPixel ( QE, fPulserColor );
// - Set Variance of QE in  MCalibrationQEPix::SetQEBlindPixelVar ( Variance, fPulserColor );
// - Set bit MCalibrationQEPix::SetBlindPixelMethodValid(kTRUE,fPulserColor) 
//
// - Call MCalibrationQEPix::UpdateBlindPixelMethod()
//
void MCalibrationChargeCalc::FinalizeBlindPixelQECam()
{

  const UInt_t npixels  = fGeom->GetNumPixels();

  //
  // Set the results in the MCalibrationChargeCam
  //
  if (fBlindPixel)
    {
      if (fBlindPixel->IsFluxInsidePlexiglassAvailable())
        {

          const Float_t photons = fBlindPixel->GetFluxInsidePlexiglass() * (*fGeom)[0].GetA()
                               / fQECam->GetPlexiglassQE();
          fCam->SetNumPhotonsBlindPixelMethod(photons);

          const Float_t photrelvar = fBlindPixel->GetFluxInsidePlexiglassRelVar() 
                                   + fQECam->GetPlexiglassQERelVar();
          if (photrelvar > 0.)
            fCam->SetNumPhotonsBlindPixelMethodErr(TMath::Sqrt( photrelvar * photons * photons));  
        }
    }
  //
  //  With the knowledge of the overall photon flux, calculate the 
  //  quantum efficiencies after the Blind Pixel and PIN Diode method
  //
  for (UInt_t i=0; i<npixels; i++)
    {
      
      MCalibrationQEPix   &qepix = (MCalibrationQEPix&)  (*fQECam)[i];

      if (!fBlindPixel)
        {
          qepix.SetBlindPixelMethodValid(kFALSE, fPulserColor);
          continue;
        }

      if (!fBlindPixel->IsFluxInsidePlexiglassAvailable())
        {
          qepix.SetBlindPixelMethodValid(kFALSE, fPulserColor);
          continue;
        }
      
      MBadPixelsPix &bad = (*fBadPixels)[i];
      if (bad.IsUnsuitable (MBadPixelsPix::kUnsuitableRun))
        {
          qepix.SetBlindPixelMethodValid(kFALSE, fPulserColor);
          continue;
        }
      
      MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
      MGeomPix              &geo =                        (*fGeom)[i];
      
      const Float_t qe        = pix.GetPheFFactorMethod() 
                             / fBlindPixel->GetFluxInsidePlexiglass() 
                             / geo.GetA() 
                             * fQECam->GetPlexiglassQE();

      const Float_t qerelvar = fBlindPixel->GetFluxInsidePlexiglassRelVar() 
                               + fQECam->GetPlexiglassQERelVar()
                               + pix.GetPheFFactorMethodRelVar();

      qepix.SetQEBlindPixel    ( qe            , fPulserColor );
      qepix.SetQEBlindPixelVar ( qerelvar*qe*qe, fPulserColor );      
      qepix.SetBlindPixelMethodValid(  kTRUE   , fPulserColor );

      if (!qepix.UpdateBlindPixelMethod())
        *fLog << warn << GetDescriptor() 
              << ": Cannot update Quantum efficiencies with the Blind Pixel Method" << endl;
    }
}

// ------------------------------------------------------------------------
//
// Loop over pixels: 
//
// - Continue, if not MCalibrationChargePINDiode::IsFluxOutsidePlexiglassAvailable() and set:
//                    MCalibrationQEPix::SetPINDiodeMethodValid(kFALSE,fPulserColor)
//
// - Calculate the quantum efficiency with the formula:
// 
//   QE =  Num.Phes / MCalibrationChargePINDiode::GetFluxOutsidePlexiglass() / MGeomPix::GetA()
//
// - Set QE in MCalibrationQEPix::SetQEPINDiode ( QE, fPulserColor );
// - Set Variance of QE in  MCalibrationQEPix::SetQEPINDiodeVar ( Variance, fPulserColor );
// - Set bit MCalibrationQEPix::SetPINDiodeMethodValid(kTRUE,fPulserColor) 
//
// - Call MCalibrationQEPix::UpdatePINDiodeMethod()
//
void MCalibrationChargeCalc::FinalizePINDiodeQECam()
{
  
  const UInt_t npixels  = fGeom->GetNumPixels();

  //
  //  With the knowledge of the overall photon flux, calculate the 
  //  quantum efficiencies after the PIN Diode method
  //
  for (UInt_t i=0; i<npixels; i++)
    {
      
      MCalibrationQEPix   &qepix = (MCalibrationQEPix&)  (*fQECam)[i];

      if (!fPINDiode)
        {
          qepix.SetPINDiodeMethodValid(kFALSE, fPulserColor);
          continue;
        }
      
      if (!fPINDiode->IsFluxOutsidePlexiglassAvailable())
        {
          qepix.SetPINDiodeMethodValid(kFALSE, fPulserColor);
          continue;
        }

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

      if (!bad.IsUnsuitable (MBadPixelsPix::kUnsuitableRun))
        {
          qepix.SetPINDiodeMethodValid(kFALSE, fPulserColor);
          continue;
        }
      
      MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i];
      MGeomPix              &geo =                        (*fGeom)[i];
      
      const Float_t qe       =  pix.GetPheFFactorMethod() 
                               / fPINDiode->GetFluxOutsidePlexiglass() 
                               / geo.GetA();

      const Float_t qerelvar = fPINDiode->GetFluxOutsidePlexiglassRelVar() + pix.GetPheFFactorMethodRelVar();

      qepix.SetQEPINDiode    ( qe            , fPulserColor );
      qepix.SetQEPINDiodeVar ( qerelvar*qe*qe, fPulserColor );      
      qepix.SetPINDiodeMethodValid(  kTRUE   , fPulserColor );

      if (!qepix.UpdatePINDiodeMethod())
        *fLog << warn << GetDescriptor() 
              << ": Cannot update Quantum efficiencies with the PIN Diode Method" << endl;
    }
}

// -----------------------------------------------------------------------------------------------
//
// - Print out statistics about BadPixels of type UnsuitableType_t 
// - store numbers of bad pixels of each type in fCam
//
void MCalibrationChargeCalc::FinalizeUnsuitablePixels()
{
  
  *fLog << inf << endl;
  *fLog << GetDescriptor() << ": Charge Calibration status:" << endl;
  *fLog << dec << setfill(' ');

  const Int_t nareas = fGeom->GetNumAreas();

  Int_t counts[nareas];
  memset(counts,0,nareas*sizeof(Int_t));

  for (Int_t i=0; i<fBadPixels->GetSize(); i++)
    {
      MBadPixelsPix &bad = (*fBadPixels)[i];
      if (!bad.IsBad())
        {
          const Int_t aidx = (*fGeom)[i].GetAidx();
          counts[aidx]++;
        }
    }

  if (fGeom->InheritsFrom("MGeomCamMagic"))
    *fLog << " " << setw(7) << "Successfully calibrated Pixels: " 
          << Form("%s%3i%s%3i","Inner: ",counts[0]," Outer: ",counts[1]) << endl;
  
  memset(counts,0,nareas*sizeof(Int_t));

  for (Int_t i=0; i<fBadPixels->GetSize(); i++)
    {
      MBadPixelsPix &bad = (*fBadPixels)[i];
      if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
        {
          const Int_t aidx = (*fGeom)[i].GetAidx();
          counts[aidx]++;
        }
    }

  for (Int_t aidx=0; aidx<nareas; aidx++)
    fCam->SetNumUnsuitable(counts[aidx], aidx);

  if (fGeom->InheritsFrom("MGeomCamMagic"))
    *fLog << " " << setw(7) << "Uncalibrated Pixels:            " 
          << Form("%s%3i%s%3i","Inner: ",counts[0]," Outer: ",counts[1]) << endl;

  memset(counts,0,nareas*sizeof(Int_t));

  for (Int_t i=0; i<fBadPixels->GetSize(); i++)
    {
      MBadPixelsPix &bad = (*fBadPixels)[i];
      if (bad.IsUnsuitable(MBadPixelsPix::kUnreliableRun))
        {
          const Int_t aidx = (*fGeom)[i].GetAidx();
          counts[aidx]++;
        }
    }

  for (Int_t aidx=0; aidx<nareas; aidx++)
    fCam->SetNumUnreliable(counts[aidx], aidx);

  *fLog << " " << setw(7) << "Unreliable Pixels:              "
        << Form("%s%3i%s%3i","Inner: ",counts[0]," Outer: ",counts[1]) << endl;

}

// -----------------------------------------------------------------------------------------------
//
// Print out statistics about BadPixels of type UncalibratedType_t 
// 
void MCalibrationChargeCalc::PrintUncalibrated(MBadPixelsPix::UncalibratedType_t typ, const char *text) const 
{
  
  UInt_t countinner = 0;
  UInt_t countouter = 0;

  for (Int_t i=0; i<fBadPixels->GetSize(); i++)
    {
      MBadPixelsPix &bad = (*fBadPixels)[i];
      if (bad.IsUncalibrated(typ))
        {
          if (fGeom->GetPixRatio(i) == 1.)
            countinner++;
          else
            countouter++;
        }
    }

  *fLog << " " << setw(7) << text  
        << Form("%s%3i%s%3i","Inner: ",countinner," Outer: ",countouter) << endl;
}

// --------------------------------------------------------------------------
//
// Set the path for output file
// 
void MCalibrationChargeCalc::SetOutputPath(TString path)
{
  fOutputPath = path;
  if (fOutputPath.EndsWith("/"))
    fOutputPath = fOutputPath(0, fOutputPath.Length()-1);
}

void MCalibrationChargeCalc::SetOutputFile(TString file)
{ 
  fOutputFile        = file; 
}

// --------------------------------------------------------------------------
//
// Get the output file
// 
const char* MCalibrationChargeCalc::GetOutputFile()
{
  return Form("%s/%s", (const char*)fOutputPath, (const char*)fOutputFile);
}