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

/* ======================================================================== *\
!
! *
! * 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 from the FADC
//   time slices. The integrated time slices have to be delivered by an 
//   MExtractedSignalCam. The pedestals by an MPedestalCam.
//
//   The output container MCalibrationCam holds one entry of type MCalibrationChargePix 
//   for every pixel. It is filled in the following way:
//
//   ProProcess: Initialize MCalibrationCam
//               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 done by this class, histograms are filled by 
//               MHCalibrationChargeCam
//
//   PostProcess:  Fit results from MHCalibrationChargeCam are retrieved
//                and used for the calculation of the reduced sigma, 
//                the F-Factor method, the blind pixel method (photon flux 
//                                     inside plexiglass) and
//                                     the PINDiode method (photon flux
//                                     outside plexiglass)
//
//                Hi-Gain vs. Lo-Gain Calibration (very memory-intensive)
//                can be skipped with the command:
//                MalibrationCam::SkipHiLoGainCalibration()
//
//  Input Containers:
//   MRawEvtData
//   MPedestalCam
//   MBadPixelsCam
//
//  Output Containers:
//   MCalibrationCam
//   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 "MLog.h"
#include "MLogManip.h"

#include "MParList.h"

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

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

#include "MCalibrationChargeCam.h"
#include "MCalibrationChargePix.h"
#include "MCalibrationChargePINDiode.h"
#include "MCalibrationChargeBlindPix.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"


ClassImp(MCalibrationChargeCalc);

using namespace std;

const Float_t MCalibrationChargeCalc::fgChargeLimit              = 3.;
const Float_t MCalibrationChargeCalc::fgChargeErrLimit           = 0.;
const Float_t MCalibrationChargeCalc::fgChargeRelErrLimit        = 1.;
const Float_t MCalibrationChargeCalc::fgTimeLowerLimit           = 1.;
const Float_t MCalibrationChargeCalc::fgTimeUpperLimit           = 2.;
// --------------------------------------------------------------------------
//
// Default constructor. 
//
MCalibrationChargeCalc::MCalibrationChargeCalc(const char *name, const char *title)
    : fPedestals(NULL), fCam(NULL), fQECam(NULL),
      fRawEvt(NULL), fRunHeader(NULL), fGeom(NULL), 
      fBadPixels(NULL), fEvtTime(NULL),
      fSignals(NULL), fSigBlind(NULL), fSigPIN(NULL),
      fPINDiode(NULL), fBlindPixel(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");

    Clear();

    SetChargeLimit();
    SetChargeErrLimit();  
    
    SetChargeRelErrLimit();
    SetTimeLowerLimit();
    SetTimeUpperLimit();
}

void MCalibrationChargeCalc::Clear(const Option_t *o)
{
  
    SETBIT(fFlags, kUseQualityChecks);
    SETBIT(fFlags, kHiLoGainCalibration);

    fNumHiGainSamples  = 0.;
    fNumLoGainSamples  = 0.;
    fSqrtHiGainSamples = 0.;
    fSqrtLoGainSamples = 0.;
    fConversionHiLo    = 0;
    SkipQualityChecks      ( kFALSE );
    SkipHiLoGainCalibration( kFALSE );    

}


// --------------------------------------------------------------------------
//
// The PreProcess searches for the following input containers:
//  - MRawEvtData
//  - MPedestalCam
//  - MExtractedSignalCam
//  - MExtractedSignalBlindPixel
//  - MExtractedSignalPINDiode
//
// The following output containers are also searched and created if
// they were not found:
//
//  - MCalibrationCam
//  - MCalibrationQECam
//
// The following output containers are only searched, but not created
//
//  - MTime
//
Int_t MCalibrationChargeCalc::PreProcess(MParList *pList)
{
  
  fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
  if (!fRawEvt)
    {
      *fLog << err << "MRawEvtData not found... aborting." << endl;
      return kFALSE;
    }
  
  fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");
  if (!fPedestals)
    {
      *fLog << err << "MPedestalCam not found... aborting" << endl;
      return kFALSE;
    }
  
  fSignals = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
  if (!fSignals)
    {
      *fLog << err << "MExtractedSignalCam not found... aborting" << endl;
      return kFALSE;
    }
  
  fSigBlind = (MExtractedSignalBlindPixel*)pList->FindObject("MExtractedSignalBlindPixel");
  if (!fSigBlind)
    *fLog << warn << "MExtractedSignalBlindPixel not found... no blind pixel method! " << endl;
  else
    {
      fBlindPixel = (MCalibrationChargeBlindPix*)pList->FindCreateObj("MCalibrationChargeBlindPix");
      if (!fBlindPixel)
        {
          *fLog << err << "Cannot find nor create MCalibrationChargeBlindPix... aborting" << endl;
          return kFALSE;
        }
    }
  
  fSigPIN  = (MExtractedSignalPINDiode*)pList->FindObject("MExtractedSignalPINDiode");
  if (!fSigPIN)
    *fLog << warn << "MExtractedSignalPINDiode not found... no PIN Diode method! " << endl;
  else
    {
      fPINDiode = (MCalibrationChargePINDiode*)pList->FindCreateObj("MCalibrationChargePINDiode");
      if (!fPINDiode)
        {
          *fLog << err << "Cannot find nor create MCalibrationChargePINDiode... aborting" << endl;
          return kFALSE;
        }
    }
  
  fCam = (MCalibrationChargeCam*)pList->FindCreateObj("MCalibrationChargeCam");
  if (!fCam)
    {
      *fLog << err << "Cannot find nor create MCalibrationChargeCam... aborting" << endl;
      return kFALSE;
    }
  
  fQECam = (MCalibrationQECam*)pList->FindCreateObj("MCalibrationQECam");
  if (!fQECam)
    {
      *fLog << err << "Cannot find nor create MCalibrationQECam... aborting" << endl;
      return kFALSE;
    }
  
  fEvtTime = (MTime*)pList->FindObject("MTime");
  
  return kTRUE;
}


// --------------------------------------------------------------------------
//
// The ReInit searches for the following input containers:
//  - MRawRunHeader
//  - MGeomCam
//  - MBadPixelsCam
//
// It retrieves the following variables from MExtractedSignalCam:
//
//  fNumHiGainSamples 
//  fNumLoGainSamples 
//
//  fFirstUsedSliceHiGain
//  fLastUsedSliceHiGain
//  fFirstUsedSliceLoGain
//  fLastUsedSliceLoGain
//
// It defines the PixId of every pixel in MCalibrationChargeCam and MCalibrationQECam
// It sets all pixels excluded which have the flag fBadBixelsPix::IsBad() set.
//
Bool_t MCalibrationChargeCalc::ReInit(MParList *pList )
{

  fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
  if (!fRunHeader)
    {
      *fLog << err << "MRawRunHeader not found... aborting." << endl;
      return kFALSE;
    }

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

  fNumHiGainSamples  =  fSignals->GetNumUsedHiGainFADCSlices();
  fNumLoGainSamples  =  fSignals->GetNumUsedLoGainFADCSlices();
  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;
      }

    }

    return kTRUE;
}


Int_t MCalibrationChargeCalc::Process()
{
  return kTRUE;
}

// --------------------------------------------------------------------------
//  
// Finalize pedestals: 
// 
// - Retrieve pedestal and pedestal RMS from MPedestalPix
// - Retrieve total entries from MPedestalCam
// - Sum up pedestal and pedestalRMS for the average pixel
// - Set pedestal*number of used samples in MCalibrationChargePix
// - Set pedestal RMS * sqrt of number of used samples in MCalibrationChargePix
// 
//
void MCalibrationChargeCalc::FinalizePedestals(const MPedestalPix &ped, MCalibrationChargePix &cal,
                                               Float_t &avped, Float_t &avrms)
{
  
  //
  // 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());

  //
  // Calculate the average pedestal
  //
  avped += pedes;
  avrms += prms;
  
  //
  // set them in the calibration camera
  //
  if (cal.IsHiGainSaturation())
    {
      cal.SetPedestal(pedes* fNumLoGainSamples,
                      prms * fSqrtLoGainSamples,
                      prms * fNumLoGainSamples / num);
      cal.CalcLoGainPedestal((Float_t)fNumLoGainSamples);
    }
  else
    {
      cal.SetPedestal(pedes* fNumHiGainSamples,
                      prms * fSqrtHiGainSamples,
                      prms * fNumHiGainSamples / num);
    }
  
}

void MCalibrationChargeCalc::FinalizeAvPedestals(MCalibrationChargePix &cal, 
                                                 Float_t avped, Float_t avrms, Int_t avnum)
{
  
  //
  // set the pedestans in the calibration camera
  //
  if (cal.IsHiGainSaturation())
    {
      cal.SetPedestal(avped/avnum * fNumLoGainSamples,
                      avrms/avnum * fSqrtLoGainSamples,
                      avrms/avnum * fSqrtLoGainSamples/avnum);
      cal.CalcLoGainPedestal((Float_t)fNumLoGainSamples);
    }
  else
    {
      cal.SetPedestal(avped/avnum * fNumHiGainSamples,
                      avrms/avnum * fSqrtHiGainSamples,
                      avrms/avnum * fSqrtHiGainSamples/avnum);
    }
}

// ---------------------------------------------------------------------
//
// Finalize charges per pixel: 
// - Check chage validity
// - Check absolute time validity
// - Calculate the reduced sigma
// - Calculate the number of photo-electrons
//
Bool_t MCalibrationChargeCalc::FinalizeCharges(MCalibrationChargePix &cal, MBadPixelsPix &bad)
{
  
  //
  // The check return kTRUE if:
  //
  // 1) Pixel has a fitted charge greater than fChargeLimit*PedRMS
  // 2) Pixel has a fit error greater than fChargeVarLimit
  // 3) Pixel has a fitted charge greater its fChargeRelVarLimit times its charge error
  // 4) Pixel has a charge sigma bigger than its Pedestal RMS
  // 
  if (cal.GetMean() < fChargeLimit*cal.GetPedRms())
    {
      *fLog << warn << "WARNING: Fitted Charge: " << cal.GetMean() << " is smaller than "
            << fChargeLimit << " Pedestal RMS: " <<  cal.GetPedRms() << " in Pixel  " << cal.GetPixId() << endl;
      bad.SetUncalibrated( MBadPixelsPix::kChargeIsPedestal);
      bad.SetUnsuitable(   MBadPixelsPix::kUnsuitableRun   );
    }
  
  if (cal.GetMeanErr() < fChargeErrLimit) 
    {
      *fLog << warn << "WARNING: Sigma of Fitted Charge: " << cal.GetMeanErr() << " is smaller than "
            << fChargeErrLimit << " in Pixel  " << cal.GetPixId() << endl;
      bad.SetUncalibrated( MBadPixelsPix::kChargeErrNotValid );
      bad.SetUnsuitable(   MBadPixelsPix::kUnsuitableRun     );
    }
      
   if (cal.GetMean() < fChargeRelErrLimit*cal.GetMeanErr()) 
    {
      *fLog << warn << "WARNING: Fitted Charge: " << cal.GetMean() << " is smaller than "
            << fChargeRelErrLimit << "* its error: " << cal.GetMeanErr() << " in Pixel  " << cal.GetPixId() << endl;
      bad.SetUncalibrated( MBadPixelsPix::kChargeRelErrNotValid );
      bad.SetUnsuitable(   MBadPixelsPix::kUnsuitableRun        );
    }

  if (cal.GetSigma() < cal.GetPedRms())
    {
      *fLog << warn << "WARNING: Sigma of Fitted Charge: " << cal.GetSigma() 
            << " smaller than Pedestal RMS: " << cal.GetPedRms() << " in Pixel  " << cal.GetPixId() << endl;
      bad.SetUncalibrated( MBadPixelsPix::kChargeSigmaNotValid );
      bad.SetUnsuitable(   MBadPixelsPix::kUnsuitableRun       );
    }

  //
  // The check returns kTRUE if:
  //
  // The mean arrival time is at least fTimeLowerLimit slices from the lower edge 
  // and fUpperLimit slices from the upper edge
  //
  const Byte_t loweredge  = cal.IsHiGainSaturation() ? fSignals->GetFirstUsedSliceLoGain() 
                                                     : fSignals->GetFirstUsedSliceHiGain();
  const Byte_t upperedge  = cal.IsHiGainSaturation() ? fSignals->GetLastUsedSliceLoGain()
                                                     : fSignals->GetLastUsedSliceHiGain();

  const Float_t lowerlimit = (Float_t)loweredge + fTimeLowerLimit;
  const Float_t upperlimit = (Float_t)upperedge + fTimeUpperLimit;  

  if ( cal.GetAbsTimeMean() < lowerlimit)
    {
      *fLog << warn << "WARNING: Mean ArrivalTime in first " << fTimeLowerLimit 
            << " extraction bin of the Pixel " << cal.GetPixId() << endl;
      *fLog << cal.GetAbsTimeMean() << "   " << lowerlimit << endl;
      bad.SetUncalibrated( MBadPixelsPix::kMeanTimeInFirstBin );
      bad.SetUnsuitable(   MBadPixelsPix::kUnsuitableRun      );
    }
  
  if ( cal.GetAbsTimeMean() > upperlimit )
    {
      *fLog << warn << "WARNING: Mean ArrivalTime in last " << fTimeUpperLimit 
            << " two extraction bins of the Pixel " << cal.GetPixId() << endl;
      *fLog << cal.GetAbsTimeMean() << "   " << upperlimit << endl;
      bad.SetUncalibrated( MBadPixelsPix::kMeanTimeInLast2Bins );
      bad.SetUnsuitable(   MBadPixelsPix::kUnsuitableRun       );
    }
      
  if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
    return kFALSE;

  if (!cal.CalcReducedSigma())
    {
      *fLog << warn << GetDescriptor() 
            << ": Could not calculate reduced sigmas of pixel: " << cal.GetPixId() << endl;
      bad.SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
      return kFALSE;
    }
  
  if (!cal.CalcFFactorMethod())
    {
      *fLog << warn << GetDescriptor() 
            << ": Could not calculate F-Factor of pixel: " << cal.GetPixId() << endl;
      bad.SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
      return kFALSE;
    }
  return kTRUE;
}

// ------------------------------------------------------------------------
//
// Returns kFALSE if pointer to MExtractedSignalPINDiode is NULL
// Returns kFALSE if pointer to MCalibrationChargePINDiode is NULL
//
// The check return kTRUE if:
//
// 1) PINDiode has a fitted charge greater than fChargeLimit*PedRMS
// 2) PINDiode has a fit error greater than fChargeErrLimit
// 3) PINDiode has a fitted charge greater its fChargeRelErrLimit times its charge error
// 4) PINDiode has a charge sigma bigger than its Pedestal RMS
// 5) The mean arrival time is at least fTimeLowerLimit slices from the lower edge 
//    and fUpperLimit slices from the upper edge
// 
Bool_t MCalibrationChargeCalc::FinalizePINDiode()
{

  if (!fSigPIN)
    return kFALSE;
  
  if (!fPINDiode)
    return kFALSE;  

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

  const Byte_t  loweredge  = fSigPIN->GetFirstUsedSlice();
  const Byte_t  upperedge  = fSigPIN->GetLastUsedSlice();
  const Float_t lowerlimit = (Float_t)loweredge + fTimeLowerLimit;
  const Float_t upperlimit = (Float_t)upperedge + fTimeUpperLimit;  

  if (fPINDiode->GetAbsTimeMean() < lowerlimit)
    {
      *fLog << warn << "WARNING: Mean ArrivalTime in first " << fTimeLowerLimit 
            << " extraction bin in PIN Diode " << endl;
      *fLog << fPINDiode->GetAbsTimeMean() << "   " << lowerlimit << endl;
      return kFALSE;
    }
  
  if ( fPINDiode->GetAbsTimeMean() > upperlimit )
    {
      *fLog << warn << "WARNING: Mean ArrivalTime in last " << fTimeUpperLimit 
            << " two extraction bins in PIN Diode " << endl;
      *fLog << fPINDiode->GetAbsTimeMean() << "   " << upperlimit << 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;
}


// -----------------------------------------------------------------------
//
// - Finalize the pedestals
// - Do the quality checks 
// - Calculate the reduced sigma 
// - Calculate the F-Factor Method
//
Int_t MCalibrationChargeCalc::PostProcess()
{

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

  //
  // loop over the pedestal events and check if we have calibration
  //
  Int_t   nvalid      = 0;
  Float_t avinnerped  = 0.;
  Float_t avinnerprms = 0.;
  Int_t   avinnernum  = 0;
  Float_t avouterped  = 0.;
  Float_t avouterprms = 0.;
  Int_t   avouternum  = 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];

      if (fGeom->GetPixRatio(pixid) == 1.)
        {
          FinalizePedestals(ped,pix,avinnerped,avinnerprms);
          avinnernum++;
        }
      else
        {
          FinalizePedestals(ped,pix,avouterped,avouterprms);
          avouternum++;
        }

      if (FinalizeCharges(pix,bad))
          nvalid++;
    }
  
  //
  // 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++)
    {

      FinalizeAvPedestals((MCalibrationChargePix&)fCam->GetAverageArea(aidx), 
                          avinnerped, avinnerprms,avinnernum);
      FinalizeCharges((MCalibrationChargePix&)fCam->GetAverageArea(aidx),
                      fCam->GetAverageBadArea(aidx));
    }
  
  for (UInt_t sector=0; sector<fGeom->GetNumSectors(); sector++)
    {

      FinalizeAvPedestals((MCalibrationChargePix&)fCam->GetAverageSector(sector), 
                          avinnerped, avinnerprms,avinnernum);
      FinalizeCharges((MCalibrationChargePix&)fCam->GetAverageSector(sector),
                      fCam->GetAverageBadSector(sector));
    }
  
  // 
  // F-Factor calibration
  //
  if (fCam->CalcMeanFluxPhotonsFFactorMethod(*fGeom, *fBadPixels))
  {
    fCam->ApplyFFactorCalibration(*fGeom,*fBadPixels);
    fCam->SetFFactorMethodValid(kTRUE);
  }
  else
    {
      *fLog << warn << "Could not calculate the flux of photo-electrons from the F-Factor method, " << endl;
      fCam->SetFFactorMethodValid(kFALSE);
    }
  
  // 
  // Blind Pixel calibration
  //
  if (!fBlindPixel->CheckChargeFitValidity())
  {
      *fLog << warn << "Could not calculate the flux of photons from the Blind Pixel, "
            << "charge fit not valid " << endl;
      fCam->SetBlindPixelMethodValid(kFALSE);
  }
  else
  { 
      if (!fBlindPixel->CalcFluxInsidePlexiglass())
      {
          *fLog << warn << "Could not calculate the flux of photons from the Blind Pixel, "
                << "will skip PIN Diode Calibration " << endl;
          fCam->SetBlindPixelMethodValid(kFALSE);
      }
      else
      {
          fCam->SetBlindPixelMethodValid(kTRUE);
          fCam->ApplyBlindPixelCalibration(*fGeom,*fBadPixels, *fBlindPixel);
      }
  }

  // 
  // PIN Diode calibration
  //
  if (!FinalizePINDiode())
    fCam->SetPINDiodeMethodValid(kFALSE); 
  else
    {
      fCam->SetPINDiodeMethodValid(kTRUE);
      fCam->ApplyPINDiodeCalibration(*fGeom,*fBadPixels, *fPINDiode);
    }

  fCam->SetReadyToSave();
  
  *fLog << inf << endl;
  *fLog << GetDescriptor() << ": Calibration statistics:" << endl;
  *fLog << dec << setfill(' ');

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

  *fLog << " " << setw(7) << "Successfully calibrated Pixels: " 
        << "Inner: " << countinner << " Outer: " << countouter << endl;
  
  PrintUnsuitable(MBadPixelsPix::kUnsuitableRun, "Bad Pixels:                     ");
  PrintUnsuitable(MBadPixelsPix::kUnreliableRun, "Unreliable Pixels:              ");

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

  PrintUncalibrated(MBadPixelsPix::kChargeIsPedestal,    
                    Form("%s%2.1f%s","Signal less than ",fChargeLimit," Pedestal RMS:                "));
  PrintUncalibrated(MBadPixelsPix::kChargeErrNotValid, 
                    Form("%s%2.1f%s","Signal Error smaller than ",fChargeErrLimit,":                    "));
  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 ",fTimeLowerLimit," Bin(s):          "));
  PrintUncalibrated(MBadPixelsPix::kMeanTimeInLast2Bins, 
                    Form("%s%2.1f%s","Mean Abs. Arr. Time in Last ",fTimeUpperLimit," Bin(s):           "));
  PrintUncalibrated(MBadPixelsPix::kHiGainOscillating,   
                    "Pixels with changing Hi Gain signal over time:    ");
  PrintUncalibrated(MBadPixelsPix::kLoGainOscillating,   
                    "Pixels with changing Lo Gain signal over time:    ");
  PrintUncalibrated(MBadPixelsPix::kDeviatingNumPhes,    
                    "Pixels with deviating number of phes:             ");
  PrintUncalibrated(MBadPixelsPix::kHiGainNotFitted,     
                    "Pixels with unsuccesful Gauss fit to the Hi Gain: ");
  PrintUncalibrated(MBadPixelsPix::kLoGainNotFitted,     
                    "Pixels with unsuccesful Gauss fit to the Lo Gain: ");

  return kTRUE;
}

void MCalibrationChargeCalc::PrintUnsuitable(MBadPixelsPix::UnsuitableType_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.IsUnsuitable(typ))
        {
          if (fGeom->GetPixRatio(i) == 1.)
            countinner++;
          else
            countouter++;
        }
    }

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

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;
}