/* ======================================================================== *\
!
! *
! * 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 "MLog.h"
#include "MLogManip.h"

ClassImp(MCalibrationPix);

using namespace std;

const Float_t MCalibrationPix::gkElectronicPedRms         = 1.5;
const Float_t MCalibrationPix::gkElectronicPedRmsErr      = 0.3;
const Float_t MCalibrationPix::gkFFactor                  = 1.15;
const Float_t MCalibrationPix::gkFFactorErr               = 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::gkConvFFactorRelErrLimit   = 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::gkConversionHiLoErr        = 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";

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

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

  Clear();

  //
  // At the moment, we don't have a database, yet, 
  // so we get it from the configuration file
  //
  SetConversionHiLo();
  SetConversionHiLoErr();

  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.;
  fChargeErr                        =  -1.;
  fSigmaCharge                      =  -1.;
  fSigmaChargeErr                   =  -1.;
  fRSigmaCharge                     =  -1.;
  fRSigmaChargeErr                  =  -1.;

  fChargeProb                       =  -1.;
  fPed                              =  -1.;
  fPedRms                           =  -1.;
  fPedRmsErr                        =  -1.;

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

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

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

  fPheFFactorMethod                 =  -1.;
  fPheFFactorMethodErr              =  -1.;

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

  fConversionFFactorMethodErr       =  -1.;
  fConversionBlindPixelMethodErr    =  -1.;
  fConversionPINDiodeMethodErr      =  -1.;
  fConversionCombinedMethodErr      =  -1.;

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

  fTotalFFactorFFactorMethod        =  -1.;
  fTotalFFactorBlindPixelMethod     =  -1.;
  fTotalFFactorPINDiodeMethod       =  -1.;
  fTotalFFactorCombinedMethod       =  -1.;

}


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;
  fConversionFFactorMethodErr   = 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;
  fConversionCombinedMethodErr   = 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;
  fConversionBlindPixelMethodErr   = 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 ;
  fConversionPINDiodeMethodErr   = 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;
  
}

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()         const
{ 
  return TESTBIT(fFlags, kOscillating);
}

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

Bool_t MCalibrationPix::IsFFactorMethodValid()   const
{ 
  return TESTBIT(fFlags, kFFactorMethodValid);     
}

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

Bool_t MCalibrationPix::IsCombinedMethodValid()  const
{ 
  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 eans and RMS of the histogram!!
  //
  fCharge         = fHist->GetChargeMean();
  fChargeErr      = fHist->GetChargeMeanErr(); 
  fSigmaCharge    = fHist->GetChargeSigma();
  fSigmaChargeErr = 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.)
   || (fChargeErr   < 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.* fChargeErr * fChargeErr / chargeSquare;

  const Float_t chargeRelErrSquare        =     fChargeErr * fChargeErr 
                                             / (fCharge    * fCharge);
  
  const Float_t ffactorsquare             =    gkFFactor    * gkFFactor;
  const Float_t ffactorsquareRelErrSquare = 4.*gkFFactorErr * gkFFactorErr / ffactorsquare;

  const Float_t avQERelErrSquare          =    fAverageQEErr * fAverageQEErr / fAverageQE / fAverageQE; 

  const Float_t avQEFFactor               = TMath::Sqrt( ( 1. - fAverageQE ) / fAverageQE  );
  const Float_t avQEFFactorErr            = 1./ ( 2. * avQEFFactor ) * fAverageQEErr 
                                            / ( fAverageQE * fAverageQE );
  const Float_t avQEFFactorRelErrSquare  = avQEFFactorErr  * avQEFFactorErr 
                                           / ( avQEFFactor * avQEFFactor) ;
  //
  // Now the absolute error squares
  //
  const Float_t sigmaSquare               =    fSigmaCharge   * fSigmaCharge;
  const Float_t sigmaSquareErrSquare      = 4.*fSigmaChargeErr* fSigmaChargeErr * sigmaSquare;
  
  Float_t pedRmsSquare                    =       fPedRms * fPedRms;
  Float_t pedRmsSquareErrSquare           = 4.*fPedRmsErr * fPedRmsErr * 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.*gkElectronicPedRmsErr * gkElectronicPedRmsErr * 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.*fConversionHiLoErr * fConversionHiLoErr / 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);
      return kFALSE;
    }
  
  const Float_t rSigmaSquareRelErrSquare = (sigmaSquareErrSquare + pedRmsSquareErrSquare)
                                          / (rsigmachargesquare * rsigmachargesquare) ;
  
  fRSigmaCharge    = TMath::Sqrt(rsigmachargesquare);
  fRSigmaChargeErr = 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;
  //
  // Calculate the number of photons from the F-Factor method
  // FIXME: This is a preliminary solution, the qe shall be 
  // calibrated itself!
  //
  fPheFFactorMethod /= fAverageQE;

  const Float_t pheFFactorRelErrSquare =  ffactorsquareRelErrSquare
                                        + chargeSquareRelErrSquare
                                        + rSigmaSquareRelErrSquare
                                        + avQERelErrSquare;
  
  fPheFFactorMethodErr                 =  TMath::Sqrt(pheFFactorRelErrSquare) * fPheFFactorMethod;
  
  fMeanConversionFFactorMethod         =  fPheFFactorMethod / fCharge ;
  fConversionFFactorMethodErr          =  ( pheFFactorRelErrSquare + chargeRelErrSquare )
                                        * fMeanConversionFFactorMethod * fMeanConversionFFactorMethod;

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

  fSigmaConversionFFactorMethod = GetTotalFFactorFFactorMethod()*TMath::Sqrt(fMeanConversionFFactorMethod);
  
  //
  // Calculate the Total F-Factor of the camera ( in photons )
  //
  if (fPheFFactorMethod > 0)
  {      
      fTotalFFactorFFactorMethod =  (fRSigmaCharge/fCharge)*TMath::Sqrt(fPheFFactorMethod);
      fTotalFFactorFFactorMethod *=  avQEFFactor;
  }

  //
  // Calculate the error of the Total F-Factor of the camera ( in photons )
  //
  const Float_t rSigmaChargeRelErrSquare =    fRSigmaChargeErr * fRSigmaChargeErr 
                                           / (fRSigmaCharge    * fRSigmaCharge) ;

  fTotalFFactorErrFFactorMethod = TMath::Sqrt(  rSigmaChargeRelErrSquare
						+ chargeRelErrSquare
						+ pheFFactorRelErrSquare
						+ avQEFFactorRelErrSquare  );

  fTotalFFactorErrFFactorMethod *= fTotalFFactorFFactorMethod;

  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 (fChargeErr < gkChargeErrLimit) 
    {
      *fLog << warn << "WARNING: Err of Fitted Charge is smaller than "
            << gkChargeErrLimit << " in Pixel " << fPixId << endl;
      return kFALSE;
    }
      
  if (fCharge < gkChargeRelErrLimit*fChargeErr) 
    {
      *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     =    fChargeErr * fChargeErr
                                         /( fCharge    * fCharge  );
  const Float_t sigmaRelErrSquare      =    fSigmaChargeErr * fSigmaChargeErr
                                         /( fSigmaCharge    * fSigmaCharge );
  const Float_t conversionRelErrSquare =    fConversionHiLoErr * fConversionHiLoErr 
                                         /( fConversionHiLo    * fConversionHiLo  );
  
  fCharge         *= fConversionHiLo;
  fChargeErr       = TMath::Sqrt(chargeRelErrSquare + conversionRelErrSquare) * fCharge;
  
  fSigmaCharge    *= fConversionHiLo;
  fSigmaChargeErr =  TMath::Sqrt(sigmaRelErrSquare + conversionRelErrSquare) * fSigmaCharge;
  
}
