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

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
// MCalibrationPix                                                         //
//                                                                         //
// This is the storage container to hold informations about the pedestal   //
// (offset) value of one Pixel (PMT).                                      //
//                                                                         //
// The following values are initialized to meaningful values:
//
// - The Electronic Rms to 1.5 per FADC slice
// - The uncertainty about the Electronic RMS to 0.3 per slice
// - The F-Factor is assumed to have been measured in Munich to 1.13 - 1.17. 
//   with the Munich definition of the F-Factor, thus: 
//   F = Sigma(Out)/Mean(Out) * Mean(In)/Sigma(In)
//   Mean F-Factor  = 1.15
//   Error F-Factor = 0.02
//
// - Average QE: (email David Paneque, 14.2.04):
//
//  The conversion factor that comes purely from QE folded to a Cherenkov
//  spectrum has to be multiplied by:
//  * Plexiglass window -->> 0.96 X 0.96
//  * PMT photoelectron collection efficiency -->> 0.9
//  * Light guides efficiency -->> 0.94
//
//  Concerning the light guides effiency estimation... Daniel Ferenc 
//  is preparing some work (simulations) to estimate it. Yet so far, he has 
//  been busy with other stuff, and this work is still UNfinished.
//
//  The estimation I did comes from:
//  1) Reflectivity of light guide walls is 85 % (aluminum)
//  2) At ZERO degree light incidence, 37% of the light hits such walls 
//    (0.15X37%= 5.6% of light lost)
//  3) When increasing the light incidence angle, more and more light hits 
//     the walls.
//
//  However, the loses due to larger amount of photons hitting the walls is more 
//  or less counteracted by the fact that more and more photon trajectories cross 
//  the PMT photocathode twice, increasing the effective sensitivity of the PMT.
//
//  Jurgen Gebauer did some quick measurements about this issue. I attach a 
//  plot. You can see that the angular dependence is (more or less) in agreement 
//  with a CosTheta function (below 20-25 degrees),
//  which is the variation of teh entrance window cross section. So, in 
//  first approximation, no loses when increasing light incidence angle; 
//  and therefore, the factor 0.94.
//
//  So, summarizing... I would propose the following conversion factors 
//  (while working with CT1 cal box) in order to get the final number of photons 
//  from the detected measured size in ADC counts.
// 
//  Nph = ADC * FmethodConversionFactor * ConvPhe-PhFactor
// 
//  FmethodConversionFactor ; measured for individual pmts
// 
//  ConvPhe-PhFactor = 0.98 * 0.23 * 0.90 * 0.94 * 0.96 * 0.96 = 0.18
// 
//  I would not apply any smearing of this factor (which we have in nature), 
//  since we might be applying it to PMTs in the totally wrong direction.
// 
// 
/////////////////////////////////////////////////////////////////////////////
#include "MCalibrationPix.h"
#include "MCalibrationConfig.h"

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

ClassImp(MCalibrationPix);

using namespace std;

const Float_t MCalibrationPix::gkElectronicPedRms         = 1.5;
const Float_t MCalibrationPix::gkErrElectronicPedRms      = 0.3;
const Float_t MCalibrationPix::gkFFactor                  = 1.15;
const Float_t MCalibrationPix::gkFFactorError             = 0.02;
const Float_t MCalibrationPix::gkChargeLimit              = 3.;
const Float_t MCalibrationPix::gkChargeErrLimit           = 0.;
const Float_t MCalibrationPix::gkChargeRelErrLimit        = 1.;
const Float_t MCalibrationPix::gkTimeLimit                = 1.5;
const Float_t MCalibrationPix::gkTimeErrLimit             = 3.;
const Float_t MCalibrationPix::gkConvFFactorRelErrorLimit = 0.1;

const Float_t MCalibrationPix::gkAverageQE                = 0.18;     
const Float_t MCalibrationPix::gkAverageQEErr             = 0.02;  
const Float_t MCalibrationPix::gkConversionHiLo           = 10.;
const Float_t MCalibrationPix::gkConversionHiLoError      = 2.5;
// --------------------------------------------------------------------------
//
// Default Constructor: 
//
MCalibrationPix::MCalibrationPix(const char *name, const char *title)
    : fPixId(-1),
      fFlags(0)
{

  fName  = name  ? name  : "MCalibrationPixel";
  fTitle = title ? title : "Container of the MHCalibrationPixels and the fit results";

  //
  // At the moment, we don't have a database, yet, 
  // so we get it from the configuration file
  //
  fConversionHiLo      = gkConversionHiLo;
  fConversionHiLoError = gkConversionHiLoError;  

  fHist = new MHCalibrationPixel("MHCalibrationPixel","Calibration Histograms Pixel ");

  if (!fHist)
    *fLog << warn << dbginf << " Could not create MHCalibrationPixel " << endl;

  Clear();

  SetAverageQE();

}

MCalibrationPix::~MCalibrationPix() 
{
  delete fHist;
}


// ------------------------------------------------------------------------
//
// Invalidate values
//
void MCalibrationPix::Clear(Option_t *o)
{

  fHist->Reset();

  CLRBIT(fFlags, kHiGainSaturation);
  CLRBIT(fFlags, kExcluded);
  CLRBIT(fFlags, kExcludeQualityCheck);
  CLRBIT(fFlags, kChargeValid);
  CLRBIT(fFlags, kFitted);
  CLRBIT(fFlags, kOscillating);
  CLRBIT(fFlags, kBlindPixelMethodValid);
  CLRBIT(fFlags, kFFactorMethodValid);
  CLRBIT(fFlags, kPINDiodeMethodValid);
  CLRBIT(fFlags, kCombinedMethodValid);

  fCharge                           =  -1.;
  fErrCharge                        =  -1.;
  fSigmaCharge                      =  -1.;
  fErrSigmaCharge                   =  -1.;
  fRSigmaCharge                     =  -1.;
  fErrRSigmaCharge                  =  -1.;

  fChargeProb                       =  -1.;
  fPed                              =  -1.;
  fPedRms                           =  -1.;
  fErrPedRms                        =   0.;

  fNumHiGainSamples                 =  -1.;
  fNumLoGainSamples                 =  -1.;

  fTimeFirstHiGain                  =   0 ;
  fTimeLastHiGain                   =   0 ;
  fTimeFirstLoGain                  =   0 ;
  fTimeLastLoGain                   =   0 ;

  fAbsTimeMean                      =  -1.;
  fAbsTimeRms                       =  -1.;

  fPheFFactorMethod                 =  -1.;
  fPheFFactorMethodError            =  -1.;

  fMeanConversionFFactorMethod      =  -1.;
  fMeanConversionBlindPixelMethod   =  -1.;
  fMeanConversionPINDiodeMethod     =  -1.;
  fMeanConversionCombinedMethod     =  -1.;

  fErrorConversionFFactorMethod     =  -1.;
  fErrorConversionBlindPixelMethod  =  -1.;
  fErrorConversionPINDiodeMethod    =  -1.;
  fErrorConversionCombinedMethod    =  -1.;

  fSigmaConversionFFactorMethod     =  -1.;
  fSigmaConversionBlindPixelMethod  =  -1.;
  fSigmaConversionPINDiodeMethod    =  -1.;
  fSigmaConversionCombinedMethod    =  -1.;

  fFactorCalculated = kFALSE;

}


void MCalibrationPix::DefinePixId(Int_t i)
{
  
  fPixId = i;
  fHist->ChangeHistId(i);
  
}


// --------------------------------------------------------------------------
//
// Set the pedestals from outside
//
void MCalibrationPix::SetPedestal(const Float_t ped, const Float_t pedrms, 
                                  const Float_t higainsamp, const Float_t logainsamp )
{

  fPed    = ped;    
  fPedRms = pedrms;
  
  fNumHiGainSamples = higainsamp;
  fNumLoGainSamples = logainsamp;

}

// --------------------------------------------------------------------------
//
// Set the conversion factors from outside (only for MC)
//
void MCalibrationPix::SetConversionFFactorMethod(Float_t c, Float_t err, Float_t sig)
{
  fMeanConversionFFactorMethod = c;
  fErrorConversionFFactorMethod = err;
  fSigmaConversionFFactorMethod = sig;
}

// --------------------------------------------------------------------------
//
// Set the conversion factors from outside (only for MC)
//
void MCalibrationPix::SetConversionCombinedMethod(Float_t c, Float_t err, Float_t sig)
{
  fMeanConversionCombinedMethod = c;
  fErrorConversionCombinedMethod = err;
  fSigmaConversionCombinedMethod = sig;
}


// --------------------------------------------------------------------------
//
// Set the conversion factors from outside (only for MC)
//
void MCalibrationPix::SetConversionBlindPixelMethod(Float_t c, Float_t err, Float_t sig)
{
  fMeanConversionBlindPixelMethod  = c;
  fErrorConversionBlindPixelMethod = err;
  fSigmaConversionBlindPixelMethod = sig;
}

// --------------------------------------------------------------------------
//
// Set the conversion factors from outside (only for MC)
//
void MCalibrationPix::SetConversionPINDiodeMethod(Float_t c, Float_t err, Float_t sig)
{
  fMeanConversionPINDiodeMethod = c ;
  fErrorConversionPINDiodeMethod = err;
  fSigmaConversionPINDiodeMethod = sig;
}

// --------------------------------------------------------------------------
//
// Set the Hi Gain Saturation Bit from outside (only for MC)
//
void MCalibrationPix::SetHiGainSaturation(Bool_t b)
{

  if (b) 
    { 
      SETBIT(fFlags, kHiGainSaturation); 
      fHist->SetUseLoGain(1); 
    }
  else    
    {
      CLRBIT(fFlags, kHiGainSaturation); 
      fHist->SetUseLoGain(0); 
    }
}

// --------------------------------------------------------------------------
//
// Set the Excluded Bit from outside 
//
void MCalibrationPix::SetExcluded(Bool_t b )
{ 
  b ?  SETBIT(fFlags, kExcluded) : CLRBIT(fFlags, kExcluded); 
}


// --------------------------------------------------------------------------
//
// Set the Excluded Bit from outside 
//
void MCalibrationPix::SetExcludeQualityCheck(Bool_t b )
{ 
  b ?  SETBIT(fFlags, kExcludeQualityCheck) : CLRBIT(fFlags, kExcludeQualityCheck); 
}

// --------------------------------------------------------------------------
//
// Set the Excluded Bit from outside 
//
void MCalibrationPix::SetChargeValid(Bool_t b )    
{ 
  b ?  SETBIT(fFlags, kChargeValid) : CLRBIT(fFlags, kChargeValid); 
}

// --------------------------------------------------------------------------
//
// Set the Excluded Bit from outside 
//
void MCalibrationPix::SetFitted(Bool_t b )
{ 
  b ?  SETBIT(fFlags, kFitted) : CLRBIT(fFlags, kFitted); 
}
    
// --------------------------------------------------------------------------
//
// Set the Excluded Bit from outside 
//
void MCalibrationPix::SetOscillating(Bool_t b )
{ 
  b ?  SETBIT(fFlags, kOscillating) : CLRBIT(fFlags, kOscillating); 
}
    
// --------------------------------------------------------------------------
//
// Set the Excluded Bit from outside 
//
void MCalibrationPix::SetBlindPixelMethodValid(Bool_t b )
{ 
  b ?  SETBIT(fFlags, kBlindPixelMethodValid) : CLRBIT(fFlags, kBlindPixelMethodValid); 
}    

// --------------------------------------------------------------------------
//
// Set the Excluded Bit from outside 
//
void MCalibrationPix::SetFFactorMethodValid(Bool_t b )
{ 
  b ?  SETBIT(fFlags, kFFactorMethodValid) : CLRBIT(fFlags, kFFactorMethodValid); 
}    

// --------------------------------------------------------------------------
//
// Set the Excluded Bit from outside 
//
void MCalibrationPix::SetPINDiodeMethodValid(Bool_t b )  
{ 
  b ?  SETBIT(fFlags, kPINDiodeMethodValid) : CLRBIT(fFlags, kPINDiodeMethodValid); 
}

void MCalibrationPix::SetAbsTimeBordersHiGain(Byte_t f, Byte_t l)
{

  fTimeFirstHiGain = f;
  fTimeLastHiGain  = l;
  
}

void MCalibrationPix::SetAbsTimeBordersLoGain(Byte_t f, Byte_t l)
{

  fTimeFirstLoGain = f;
  fTimeLastLoGain  = l;
  
}

Float_t MCalibrationPix::GetPheFFactorMethod() 
{
  
  if (!fFactorCalculated)
    CalcFFactorMethod();
  
  return fPheFFactorMethod;
  
}

Float_t MCalibrationPix::GetPheFFactorMethodError() 
{
  
  if (!fFactorCalculated)
    CalcFFactorMethod();
  
  return fPheFFactorMethodError;
  
}


Float_t MCalibrationPix::GetTotalFFactorFFactorMethod() 
{
  if (!fFactorCalculated)
    CalcFFactorMethod();
  
  if (fPheFFactorMethod > 0)
    return (fRSigmaCharge/fCharge)*TMath::Sqrt(fPheFFactorMethod);
  else
    return -1.;
}

Float_t MCalibrationPix::GetTotalFFactorErrorFFactorMethod() 
{

  if (!fFactorCalculated)
   CalcFFactorMethod();

  const Float_t rsigmaChargeRelErrSquare = fErrRSigmaCharge * fErrRSigmaCharge 
                                           / (fRSigmaCharge * fRSigmaCharge) ;
  const Float_t rChargeRelErrSquare      = fErrCharge * fErrCharge 
                                           / (fCharge * fCharge) ;
  const Float_t rPheRelErrSquare         = fPheFFactorMethodError * fPheFFactorMethodError 
                                           / (fPheFFactorMethod * fPheFFactorMethod) ;

  return TMath::Sqrt(rsigmaChargeRelErrSquare+rChargeRelErrSquare+rPheRelErrSquare);
}


Float_t MCalibrationPix::GetTotalFFactorBlindPixelMethod() 
{
  return 1.;
}

Float_t MCalibrationPix::GetTotalFFactorErrorBlindPixelMethod() 
{

  return 1.;
}

Float_t MCalibrationPix::GetTotalFFactorPINDiodeMethod() 
{
  return 1.;
}

Float_t MCalibrationPix::GetTotalFFactorErrorPINDiodeMethod() 
{

  return 1.;
}

Float_t MCalibrationPix::GetTotalFFactorCombinedMethod() 
{
  return 1.;
}

Float_t MCalibrationPix::GetTotalFFactorErrorCombinedMethod() 
{

  return 1.;
}

//
// FIXME: This is a preliminary solution, the qe shall be 
// calibrated itself!
//
Float_t MCalibrationPix::GetMeanConversionFFactorMethod()
{

  if (!fFactorCalculated)
    CalcFFactorMethod();

  return fMeanConversionFFactorMethod/fAverageQE;
  
}

Float_t MCalibrationPix::GetErrorConversionFFactorMethod()
{

  if (!fFactorCalculated)
    CalcFFactorMethod();
  
  Float_t var = fErrorConversionFFactorMethod*fErrorConversionFFactorMethod
             / (fMeanConversionFFactorMethod * fMeanConversionFFactorMethod);

  var += fAverageQEErr * fAverageQEErr 
       / (fAverageQE * fAverageQE);

  if (var > 0)
    return -1.;

  var = TMath::Sqrt(var);

  return var*GetMeanConversionFFactorMethod();
  
}

Float_t MCalibrationPix::GetSigmaConversionFFactorMethod()
{

  if (!fFactorCalculated)
    CalcFFactorMethod();
  
  return fSigmaConversionFFactorMethod;
  
}


Bool_t MCalibrationPix::IsExcluded() const
 { 
   return TESTBIT(fFlags,kExcluded);  
 }

Bool_t MCalibrationPix::IsExcludeQualityCheck() const
 { 
   return TESTBIT(fFlags,kExcludeQualityCheck);  
 }

Bool_t MCalibrationPix::IsHiGainSaturation() const
 { 
   return TESTBIT(fFlags,kHiGainSaturation);  
 }

Bool_t MCalibrationPix::IsChargeValid() const 
{
  return TESTBIT(fFlags, kChargeValid);  
}

Bool_t MCalibrationPix::IsFitted() const 
{ 
  return TESTBIT(fFlags, kFitted);    
}

Bool_t MCalibrationPix::IsOscillating()
{ 

  if (TESTBIT(fFlags, kOscillating))
    return kTRUE;

  if (fHist->CheckOscillations())
    {
      SETBIT(fFlags,kOscillating);
      return kTRUE;
    }
  
  return kFALSE;
}

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

Bool_t MCalibrationPix::IsFFactorMethodValid() 
{ 

  if (!fFactorCalculated)
    CalcFFactorMethod();

  return TESTBIT(fFlags, kFFactorMethodValid);     
}

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

Bool_t MCalibrationPix::IsCombinedMethodValid() 
{ 
  return TESTBIT(fFlags, kCombinedMethodValid);    
}


// --------------------------------------------------------------------------
//
// 1) Return if the charge distribution is already succesfully fitted  
//    or if the histogram is empty
// 2) Set a lower Fit range according to 1.5 Pedestal RMS in order to avoid 
//    possible remaining cosmics to spoil the fit.
// 3) Decide if the LoGain Histogram is fitted or the HiGain Histogram
// 4) Fit the histograms with a Gaussian
// 5) In case of failure set the bit kFitted to false
// 6) Retrieve the results and store them in this class
// 7) Calculate the number of photo-electrons after the F-Factor method
// 8) Calculate the errors of the F-Factor method
//
// The fits are declared valid (fFitValid = kTRUE), if:
//
// 1) Pixel has a fitted charge greater than 3*PedRMS
// 2) Pixel has a fit error greater than 0.
// 3) Pixel has a fit Probability greater than 0.0001 
// 4) Pixel has a charge sigma bigger than its Pedestal RMS
// 5) If FitTimes is used, 
//    the mean arrival time is at least 1.0 slices from the used edge slices 
// (this stage is only performed in the times fit)
//
// If the histogram is empty, all values are set to -1.
//
// The conversion factor after the F-Factor method is declared valid, if:
//
// 1) fFitValid is kTRUE
// 2) Conversion Factor is bigger than 0.
// 3) The error of the conversion factor is smaller than 10%
//
Bool_t MCalibrationPix::FitCharge() 
{

  //
  // 1) Return if the charge distribution is already succesfully fitted  
  //    or if the histogram is empty
  //
  if (fHist->IsChargeFitOK() || fHist->IsEmpty())
    return kTRUE;

  //
  // 2) Set a lower Fit range according to 1.5 Pedestal RMS in order to avoid 
  //    possible remaining cosmics to spoil the fit.
  //
  //  if (fPed && fPedRms)
  //    fHist->SetLowerFitRange(1.5*fPedRms);
  //  else
  //    *fLog << warn << "WARNING: Cannot set lower fit range: Pedestals not available" << endl;

  //
  // 3) Decide if the LoGain Histogram is fitted or the HiGain Histogram
  //
  if (fHist->UseLoGain())
    SetHiGainSaturation();

  //
  // 4) Fit the Lo Gain histograms with a Gaussian
  //
  if (fHist->FitCharge())
    SETBIT(fFlags,kFitted);
  else
    {
      *fLog << warn << "WARNING: Could not fit charges of pixel " << fPixId << endl;
      //          
      // 5) In case of failure set the bit kFitted to false
      //
      CLRBIT(fFlags,kFitted);
    }

  //
  // 6) Retrieve the results and store them in this class
  //    If fFitted is false, we get the means and RMS of the histogram!!
  //
  fCharge         = fHist->GetChargeMean();
  fErrCharge      = fHist->GetChargeMeanErr(); 
  fSigmaCharge    = fHist->GetChargeSigma();
  fErrSigmaCharge = fHist->GetChargeSigmaErr();
  fChargeProb     = fHist->GetChargeProb();


  fAbsTimeMean    = fHist->GetAbsTimeMean();
  fAbsTimeMeanErr = fHist->GetAbsTimeMeanErr();
  fAbsTimeRms     = fHist->GetAbsTimeRms();  

  if (CheckTimeFitValidity())
    SETBIT(fFlags,kTimeFitValid);
  else
    CLRBIT(fFlags,kTimeFitValid);

  //
  // Calculate the conversion factors 
  //
  if (IsHiGainSaturation())
    ApplyLoGainConversion();

  if (CheckChargeValidity())
    SETBIT(fFlags,kChargeValid);
  else
    {
      CLRBIT(fFlags,kChargeValid);
      return kFALSE;
    }
  
  return kTRUE;
  
}

//
// Calculate the number of photo-electrons after the F-Factor method
// Calculate the errors of the F-Factor method
//
Bool_t MCalibrationPix::CalcFFactorMethod()
{

  if ( (fCharge    == -1.)
   || (fErrCharge   < 0.)
   || (fSigmaCharge < 0.)
   || (fPedRms      < 0.) )
    {
      *fLog << warn << GetDescriptor() << "Cannot calculate the FFactor Method! " 
            << "Some of the needed parameters are not available ";
      CLRBIT(fFlags,kFFactorMethodValid);
      return kFALSE;
    }

  //
  // Square all variables in order to avoid applications of square root
  //
  // First the relative error squares
  //
  const Float_t chargeSquare              =       fCharge*   fCharge;
  const Float_t chargeSquareRelErrSquare  = 4.*fErrCharge*fErrCharge / chargeSquare;
  
  const Float_t ffactorsquare             =    gkFFactor      * gkFFactor;
  const Float_t ffactorsquareRelErrSquare = 4.*gkFFactorError * gkFFactorError / ffactorsquare;
  //
  // Now the absolute error squares
  //
  const Float_t sigmaSquare               =       fSigmaCharge*fSigmaCharge;
  const Float_t sigmaSquareErrSquare      = 4.*fErrSigmaCharge*fErrSigmaCharge * sigmaSquare;
  
  Float_t pedRmsSquare                    =       fPedRms*   fPedRms;
  Float_t pedRmsSquareErrSquare           = 4.*fErrPedRms*fErrPedRms * pedRmsSquare;
  
  if (!IsHiGainSaturation())
    {  /* HiGain */

      pedRmsSquare          *= fNumHiGainSamples;
      pedRmsSquareErrSquare *= fNumHiGainSamples*fNumHiGainSamples; 
    } 
  else
    {  /* LoGain */
      
      //
      // We do not know the Lo Gain Pedestal RMS, so we have to retrieve it 
      // from the HI GAIN (all calculation per slice up to now):  
      //
      // We extract the pure NSB contribution:
      //
      const Float_t elecRmsSquare          =    gkElectronicPedRms   *gkElectronicPedRms;
      const Float_t elecRmsSquareErrSquare = 4.*gkErrElectronicPedRms*gkErrElectronicPedRms * elecRmsSquare;
  
      Float_t nsbSquare             =  pedRmsSquare          - elecRmsSquare;
      Float_t nsbSquareRelErrSquare = (pedRmsSquareErrSquare + elecRmsSquareErrSquare)
                                    / (nsbSquare * nsbSquare) ;
      
      if (nsbSquare < 0.)
        nsbSquare = 0.;

      //
      // Now, we divide the NSB by the conversion factor and 
      // add it quadratically to the electronic noise
      //
      const Float_t conversionSquare             =    fConversionHiLo     *fConversionHiLo;
      const Float_t conversionSquareRelErrSquare = 4.*fConversionHiLoError*fConversionHiLoError/conversionSquare;
      
      const Float_t convertedNsbSquare          =  nsbSquare             / conversionSquare;
      const Float_t convertedNsbSquareErrSquare = (nsbSquareRelErrSquare + conversionSquareRelErrSquare)
                                                   * convertedNsbSquare * convertedNsbSquare;

      pedRmsSquare           = convertedNsbSquare           + elecRmsSquare;
      pedRmsSquareErrSquare  = convertedNsbSquareErrSquare  + elecRmsSquareErrSquare;
      
      //
      // Now, correct for the number of used FADC slices in the LoGain:
      //
      pedRmsSquare          *= fNumLoGainSamples;
      pedRmsSquareErrSquare *= fNumLoGainSamples*fNumLoGainSamples; 
      //
      // Correct also for the conversion to Hi-Gain:
      //
      pedRmsSquare          *= fConversionHiLo*fConversionHiLo;
      pedRmsSquareErrSquare *= fConversionHiLo*fConversionHiLo*fConversionHiLo*fConversionHiLo;

    }   /* if (HiGainSaturation) */
  
  //
  // Calculate the reduced sigmas
  //
  const Float_t rsigmachargesquare = sigmaSquare - pedRmsSquare;
  if (rsigmachargesquare <= 0.)
    {
      *fLog << warn 
            << "WARNING: Cannot apply F-Factor calibration: Reduced Sigma smaller than 0 in pixel " 
            << fPixId << endl;
      CLRBIT(fFlags,kFFactorMethodValid);
      fFactorCalculated = kTRUE;
      return kFALSE;
    }
  
  const Float_t rSigmaSquareRelErrSquare = (sigmaSquareErrSquare + pedRmsSquareErrSquare)
                                          / (rsigmachargesquare * rsigmachargesquare) ;
  
  fRSigmaCharge    = TMath::Sqrt(rsigmachargesquare);
  fErrRSigmaCharge = TMath::Sqrt(sigmaSquareErrSquare + pedRmsSquareErrSquare);


  //
  // Calculate the number of phe's from the F-Factor method
  // (independent on Hi Gain or Lo Gain)
  //
  fPheFFactorMethod = ffactorsquare * chargeSquare / rsigmachargesquare;

  const Float_t pheFFactorRelErrSquare =  ffactorsquareRelErrSquare
                                       + chargeSquareRelErrSquare
                                       + rSigmaSquareRelErrSquare ;
  
  fPheFFactorMethodError        =  TMath::Sqrt(pheFFactorRelErrSquare) * fPheFFactorMethod;
  
  const Float_t chargeRelErrSquare = fErrCharge*fErrCharge / (fCharge * fCharge);
  
  fMeanConversionFFactorMethod     =  fPheFFactorMethod / fCharge ;
  fErrorConversionFFactorMethod    =  ( pheFFactorRelErrSquare + chargeRelErrSquare )
                                   * fMeanConversionFFactorMethod * fMeanConversionFFactorMethod;

  const Float_t convrelerror = fErrorConversionFFactorMethod/fMeanConversionFFactorMethod;
  
  if ( (fMeanConversionFFactorMethod > 0.) && (convrelerror < gkConvFFactorRelErrorLimit))
    SETBIT(fFlags,kFFactorMethodValid);

  fFactorCalculated = kTRUE;

  fSigmaConversionFFactorMethod = GetTotalFFactorFFactorMethod()*TMath::Sqrt(fMeanConversionFFactorMethod);

  return kTRUE;
}


//
// The check returns kTRUE if:
//
// 0) Pixel has BIT fitted set: 
//    This means:
//    a)  No result is a nan
//    b)  The NDF is not smaller than fNDFLimit (5)
//    c)  The Probability is greater than gkProbLimit (default 0.001 == 99.9%)
// 1) Pixel has a fitted charge greater than 3*PedRMS
// 2) Pixel has a fit error greater than 0. 
// 3) Pixel has a fitted charge greater its charge error
// 4) Pixel has a fit Probability greater than 0.0001 
// 5) Pixel has a charge sigma bigger than its Pedestal RMS
// 
Bool_t MCalibrationPix::CheckChargeValidity()
{

  if (!IsFitted())
    return kFALSE;

  if (IsExcludeQualityCheck())
    return kTRUE;

  Float_t pedestal;

  if (!IsHiGainSaturation())  /* higain */
    pedestal = GetPedRms()*TMath::Sqrt(fNumHiGainSamples);
  else                         /*  logain */
    pedestal = GetPedRms()*TMath::Sqrt(fNumLoGainSamples);
      

  if (fCharge < gkChargeLimit*pedestal)
    {
      *fLog << warn << "WARNING: Fitted Charge is smaller than "
            << gkChargeLimit << " Pedestal RMS in Pixel " << fPixId << endl;
      return kFALSE;
    }
  
  if (fErrCharge < gkChargeErrLimit) 
    {
      *fLog << warn << "WARNING: Error of Fitted Charge is smaller than "
            << gkChargeErrLimit << " in Pixel " << fPixId << endl;
      return kFALSE;
    }
      
  if (fCharge < gkChargeRelErrLimit*fErrCharge) 
    {
      *fLog << warn << "WARNING: Fitted Charge is smaller than "
            << gkChargeRelErrLimit << "* its error in Pixel " << fPixId << endl;
      return kFALSE;
    }
      
  if (!fHist->IsChargeFitOK()) 
    {
      *fLog << warn << "WARNING: Probability of Fitted Charge too low in Pixel " 
            << fPixId << endl;
      return kFALSE;
    }

  if (fSigmaCharge < pedestal)
    {
      *fLog << warn << "WARNING: Sigma of Fitted Charge smaller than Pedestal RMS in Pixel " 
            << fPixId << endl;
      return kFALSE;
    }
  return kTRUE;
}

//
// The check return kTRUE if:
//
// 0) No value is nan
// 1) Pixel has a fitted rel. time smaller than 3*FADC slices
// 2) Pixel has a fit error greater than 0. 
// 4) Pixel has a fit Probability greater than 0.001 
// 5) The absolute arrival time is at least 1.0 slices from the used edge slices 
//
Bool_t MCalibrationPix::CheckTimeFitValidity()
{

  
  if (IsExcludeQualityCheck())
    return kTRUE;

  if (IsHiGainSaturation())
    {

      if (fAbsTimeMean < (Float_t)fTimeFirstLoGain+1)
        {
          *fLog << warn 
                << "WARNING: Some absolute times smaller than limit in Pixel " 
                << fPixId << " time: " << fAbsTimeMean 
                << " Limit: " << (Float_t)fTimeFirstLoGain+1. << endl;
          return kFALSE;
        }

      if (fAbsTimeMean > (Float_t)fTimeLastLoGain-1)
        {
          *fLog << warn 
                << "WARNING: Some absolute times bigger than limit in Pixel " 
                << fPixId << " time: " << fAbsTimeMean 
                << " Limit: " << (Float_t)fTimeLastLoGain-1. << endl;
          return kFALSE;
        }

    }
  else
    {

      if (fAbsTimeMean < (Float_t)fTimeFirstHiGain+1.)
        {
          *fLog << warn 
                << "WARNING: Some absolute times smaller than limit in Pixel " 
                << fPixId << " time: " << fAbsTimeMean 
                << " Limit: " << (Float_t)fTimeFirstHiGain+1. << endl;
          //          return kFALSE;
        }

      if (fAbsTimeMean > (Float_t)fTimeLastHiGain-1.)
        {
          *fLog << warn 
                << "WARNING: Some absolute times bigger than limit in Pixel " 
                << fPixId << " time: " << fAbsTimeMean 
                << " Limit: " << (Float_t)fTimeLastHiGain-1. << endl;
          //          return kFALSE;
        }

    }



  return kTRUE;
}


void MCalibrationPix::CheckOscillations()
{
  fHist->CheckOscillations();
}

void MCalibrationPix::ApplyLoGainConversion()
{
  
  const Float_t chargeRelErrSquare     =   fErrCharge*fErrCharge
                                         /(   fCharge *  fCharge);
  const Float_t sigmaRelErrSquare      =   fErrSigmaCharge*fErrSigmaCharge
                                         /(   fSigmaCharge *  fSigmaCharge);
  const Float_t conversionRelErrSquare =   fConversionHiLoError*fConversionHiLoError 
                                         /(fConversionHiLo    * fConversionHiLo);
  
  fCharge         *= fConversionHiLo;
  fErrCharge       = TMath::Sqrt(chargeRelErrSquare + conversionRelErrSquare) * fCharge;
  
  fSigmaCharge    *= fConversionHiLo;
  fErrSigmaCharge =  TMath::Sqrt(sigmaRelErrSquare + conversionRelErrSquare) * fSigmaCharge;
  
}
