/* ======================================================================== *\
!
! *
! * 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
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         
// MCalibrationChargePix                                                   
//                                                                         
// Storage container of the calibrated Charge of one pixel.
//                                                                         
// 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 (gkFFactor)     = 1.15
//   Error F-Factor (gkFFactorErr) = 0.02
//
// The following variables are calculated inside this class:
// -  fLoGainPedRmsSquare and fLoGainPedRmsSquareVar (see CalcLoGainPedestal())
// -  fRSigmaSquare and fRSigmaSquareVar             (see CalcReducedSigma()  )
// -  fPheFFactorMethod and fPheFFactorMethodVar     (see CalcFFactor()       )
// -  fMeanConvFADC2Phe and fMeanConvFADC2PheVar     (see CalcConvFFactor()  )
//
// The following variables are set by MHCalibrationChargeCam:
// -  fAbsTimeMean and fAbsTimeRms
// -  all variables in MCalibrationPix
//
// The following variables are set by MCalibrationChargeCalc:
// - fPed, fPedVar and fPedRms                         
// - fMeanConvFADC2Phe
// - fConvFADC2PheVar 
// - fSigmaConvFADC2Phe
// - fTotalFFactorFFactorMethod 
// - fTotalFFactorFFactorMethodVar 
//
// The following variables are not yet implemented:
// - fConversionHiLo and fConversionHiLoVar (now set fixed to 10. +- 2.5)
//
//  Error of all variables are calculated by error-propagation. Note that internally, 
//  all error variables contain Variances in order to save the CPU-intensive square rooting 
// 
//  Low-Gain variables are stored internally unconverted, i.e. directly from the summed 
//  FADC slices extraction results, but can be retrieved converted to High-Gain amplifications 
//  by calls to: GetConvertedMean() or GetConvertedSigma()
// 
// See also: MCalibrationChargeCam, MCalibrationChargeCalc,
//           MHCalibrationChargeCam, MHCalibrationChargePix
//
/////////////////////////////////////////////////////////////////////////////
#include "MCalibrationChargePix.h"

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

#include "MBadPixelsPix.h"

ClassImp(MCalibrationChargePix);

using namespace std;

const Float_t MCalibrationChargePix::gkElectronicPedRmsInner    = 1.5;
const Float_t MCalibrationChargePix::gkElectronicPedRmsOuter    = 1.8;
const Float_t MCalibrationChargePix::gkElectronicPedRmsErr      = 0.35;
const Float_t MCalibrationChargePix::gkFFactor                  = 1.15;
const Float_t MCalibrationChargePix::gkFFactorErr               = 0.02;

const Float_t MCalibrationChargePix::fgConversionHiLo           = 10.;
const Float_t MCalibrationChargePix::fgConversionHiLoErr        = 2.5;
const Float_t MCalibrationChargePix::fgPheFFactorMethodLimit    = 1.;
const Float_t MCalibrationChargePix::fgConvFFactorRelErrLimit   = 0.75;
// --------------------------------------------------------------------------
//
// Default Constructor: 
//
// Sets:
// - fCalibFlags to 0
// - fConversionHiLo to fgConversionHiLo
// - fConversionHiLoVar to square of fgConversionHiLoErr
// - fConvFFactorelErrLimit to fgConvFFactorRelErrLimit*fgConvFFactorelErrLimit
// - fPheFFactorLimit to fgPheFFactorLimit
// 
// Calls:
// - Clear()
//
MCalibrationChargePix::MCalibrationChargePix(const char *name, const char *title)
    : fCalibFlags(0)
{

  fName  = name  ? name  : "MCalibrationChargePix";
  fTitle = title ? title : "Container of the fit results of MHCalibrationChargePixs ";

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

  SetPheFFactorMethodLimit();
  SetConvFFactorRelErrLimit();
  
  Clear();
}

// ------------------------------------------------------------------------
//
// Sets:
// - all flags to kFALSE
// - all variables to -1.
//
// Calls: 
// - MCalibrationPix::Clear()
//
void MCalibrationChargePix::Clear(Option_t *o)
{

  SetFFactorMethodValid     ( kFALSE );

  fRSigmaSquare                     =  -1.;
  fRSigmaSquareVar                  =  -1.;
  
  fPed                              =  -1.;
  fPedRms                           =  -1.;
  fPedVar                           =  -1.;

  fLoGainPedRmsSquare               =  -1.;
  fLoGainPedRmsSquareVar            =  -1.;

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

  fPheFFactorMethod                 =  -1.;
  fPheFFactorMethodVar              =  -1.;

  fMeanConvFADC2Phe                 =  -1.;
  fMeanConvFADC2PheVar              =  -1.;
  fMeanFFactorFADC2Phot             =  -1.;
  fMeanFFactorFADC2PhotVar          =  -1.;  

  MCalibrationPix::Clear();
}


// --------------------------------------------------------------------------
//
// Set F-Factor Method Validity Bit from outside 
//
void MCalibrationChargePix::SetFFactorMethodValid(const Bool_t b )
{ 
  b ?  SETBIT(fCalibFlags, kFFactorMethodValid) : CLRBIT(fCalibFlags, kFFactorMethodValid); 
}    

// --------------------------------------------------------------------------
//
// Set pedestals from outside (done by MCalibrationChargeCalc)
//
void MCalibrationChargePix::SetPedestal(const Float_t ped, const Float_t pedrms, const Float_t pederr)
{

  fPed       = ped;    
  fPedRms    = pedrms;
  fPedVar    = pederr*pederr;
}

// --------------------------------------------------------------------------
//
// Set pedestals from outside (done by MCalibrationChargeCalc)
//
void MCalibrationChargePix::SetPed(const Float_t ped, const Float_t pederr)
{

  fPed       = ped;    
  fPedVar    = pederr*pederr;
}

// --------------------------------------------------------------------------
//
// Set pedestals RMS from outside (done by MHCalibrationChargeCam)
//
void MCalibrationChargePix::SetPedRMS( const Float_t pedrms, const Float_t pedrmserr)
{
  
  fPedRms    = pedrms;
  fPedRmsVar = pedrmserr*pedrmserr;
  
}


// -------------------------------------------------------------------------------
//
// Get the conversion Error Hi-Gain to Low-Gain:
// - If fConversionHiLoVar is smaller than 0 (i.e. has not yet been set), return -1.
//  
Float_t MCalibrationChargePix::GetConversionHiLoErr()  const
{
  if (fConversionHiLoVar < 0.)
    return -1.;

  return TMath::Sqrt(fConversionHiLoVar);
}

// --------------------------------------------------------------------------
//
// Get the relative variance of the conversion factor between higain and logain:
// - If fConversionHiLo    is 0,              return -1.
// - If fConversionHiLoVar is smaller than 0, return -1.
// - Else returns: fConversionHiLoVar / fConversionHiLo^2
//
const Float_t MCalibrationChargePix::GetConversionHiLoRelVar() const 
{

  if (fConversionHiLoVar < 0.)
    return -1.;

  if (fConversionHiLo == 0.)
    return -1.;

  return fConversionHiLoVar / (fConversionHiLo * fConversionHiLo);
}

// --------------------------------------------------------------------------
//
// Get the relative variance of the electronics pedestal RMS
// - If aidx is 0, return rel. variance of gkElectronicPedRmsInner
// - If aidx is 1, return rel. variance of gkElectronicPedRmsOuter
//
const Float_t MCalibrationChargePix::GetElectronicPedRmsRelVar(const Int_t aidx) const 
{

  if (aidx == 0)
    return gkElectronicPedRmsErr * gkElectronicPedRmsErr / gkElectronicPedRmsInner / gkElectronicPedRmsInner;

  if (aidx == 1)
    return gkElectronicPedRmsErr * gkElectronicPedRmsErr / gkElectronicPedRmsOuter / gkElectronicPedRmsOuter;

  return -1.;
}
 

// --------------------------------------------------------------------------
//
// Get the relative variance of the conversion factor between higain and logain:
// - If gkFFactor    is 0,              return -1.
// - If gkFFactorErr is smaller than 0, return -1.
// - Else returns: gkFFactorErr^2 / gkFFactor*^2
//
const Float_t MCalibrationChargePix::GetFFactorRelVar() const 
{

  if (gkFFactorErr < 0.)
    return -1.;

  if (gkFFactor == 0.)
    return -1.;

  return gkFFactorErr * gkFFactorErr / (gkFFactor * gkFFactor);
}
 

// --------------------------------------------------------------------------
//
// Get the pedestals RMS: 
// - Test bit kHiGainSaturation: 
//   If yes, return square root of fLoGainPedRmsSquare (if greater than 0, otherwise -1.), 
//   If no,  return fPedRms
//
Float_t  MCalibrationChargePix::GetPedRms()  const
{

  if (IsHiGainSaturation())
    if (fLoGainPedRmsSquare < 0.)
      return -1.;
    else
      return TMath::Sqrt(fLoGainPedRmsSquare);
  
  return fPedRms;
}

// --------------------------------------------------------------------------
//
// Get the Error of the pedestals RMS: 
// - Test bit kHiGainSaturation: 
//   If yes, return square root of (0.25*fLoGainPedRmsSquareVar/ fLoGainPedRmsSquare) (if greater than 0, otherwise -1.)
//   If no , return square root of (fPedVar) (if greater than 0, otherwise -1.), divided by 2. 
//
Float_t  MCalibrationChargePix::GetPedRmsErr()  const
{
  if (IsHiGainSaturation())
    if (fLoGainPedRmsSquareVar < 0.)
      return -1.;
    else
      return TMath::Sqrt(0.25*fLoGainPedRmsSquareVar/fLoGainPedRmsSquare);
  else
    if (fPedVar < 0.)
      return -1.;
    else
      return TMath::Sqrt(fPedVar)/2.;
}


// --------------------------------------------------------------------------
//
// Get the Low Gain Mean Charge converted to High Gain amplification: 
// Returns fLoGainMean multiplied with fConversionHiLo if IsHiGainSaturation(), 
//         else return fHiGainMean
//
Float_t MCalibrationChargePix::GetConvertedMean()  const 
{

  if (IsHiGainSaturation())
    return fLoGainMean * fConversionHiLo;

  return fHiGainMean;
}

// --------------------------------------------------------------------------
//
// Get the Error of the converted Low Gain Mean: 
//
// Returns -1 if the variable fLoGainMean or fLoGainMeanVar are smaller than 0.
//
// Returns the square root of the quadratic sum of the relative variances of 
// the fLoGainMean and fConversionHiLo, mulitplied with GetConvertedMean()
// in case of HiGain Saturation, 
// else return GetMeanErr()
//
Float_t MCalibrationChargePix::GetConvertedMeanErr()  const
{

  if (IsHiGainSaturation())
    {
      const Float_t logainrelvar = GetLoGainMeanRelVar();
      
      if (logainrelvar < 0.)
        return -1.;

      return TMath::Sqrt(logainrelvar + GetConversionHiLoRelVar()) * GetConvertedMean();
    }
  else
    return GetMeanErr();

}

// --------------------------------------------------------------------------
//
// Get the Low Gain Sigma converted to High Gain amplification: 
// Returns fLoGainSigma multiplied with fConversionHiLo if IsHiGainSaturation()
// else return fHiGainSigma
//
Float_t MCalibrationChargePix::GetConvertedSigma()  const 
{
  
  if (IsHiGainSaturation())
    return fLoGainSigma * fConversionHiLo;
  else
    return fHiGainSigma;
}

// --------------------------------------------------------------------------
//
// Get the Error of the converted Sigma: 
//
// Returns -1 if the variable fLoGainSigma or fLoGainSigmaVar are smaller than 0.
//
// if IsHiGainSaturatio()
// returns the square root of the quadratic sum of the relative variances of 
// the fLoGainSigma and fConversionHiLo, mulitplied with GetConvertedSigma()
// else returns GetSigmaErr()
//
Float_t MCalibrationChargePix::GetConvertedSigmaErr()  const
{

  if (IsHiGainSaturation())
    {
      if (fLoGainSigmaVar < 0.)
        return -1.;
      
      if (fLoGainSigma < 0.)
        return -1.;
      
      const Float_t sigmaRelVar =  fLoGainSigmaVar
                                /( fLoGainSigma * fLoGainSigma );

      return TMath::Sqrt(sigmaRelVar+GetConversionHiLoRelVar()) * GetConvertedSigma();
    }
  else
    return GetSigmaErr();


}

// --------------------------------------------------------------------------
//
// Get the converted reduced Sigma: 
// - If fRSigmaSquare is smaller than 0 (i.e. has not yet been set), return -1.
// - Test bit kHiGainSaturation: 
//   If yes, return square root of fRSigmaSquare, multiplied with fConversionHiLo, 
//   If no , return square root of fRSigmaSquare
//
Float_t MCalibrationChargePix::GetConvertedRSigma()  const
{
  if (fRSigmaSquare < 0)
    return -1;

  const Float_t rsigma = TMath::Sqrt(fRSigmaSquare);
  
  return IsHiGainSaturation() ? rsigma*fConversionHiLo : rsigma ;
} 

// --------------------------------------------------------------------------
//
// Get the error of the converted reduced Sigma: 
// - If fRSigmaSquareVar is smaller than 0 (i.e. has not yet been set), return -1.
// - Calculate the absolute variance of the reduced sigma with the formula:
//   reduced sigma variance = 0.25 * fRSigmaSquareVar / fRSigmaSquare
// - Test bit kHiGainSaturation: 
//   If yes, returns the square root of the quadratic sum of the relative variances of the 
//           reduced sigma and fConversionHiLo, mulitplied with GetRSigma()
//   Else returns the square root of rel. (0.25*fRSigmaSquareVar / fRSigmaSquare)
//
Float_t MCalibrationChargePix::GetConvertedRSigmaErr()  const
{

  if (fRSigmaSquareVar < 0)
    return -1;

  //
  // SigmaSquareVar = 4. * Sigma * Sigma * Var(sigma)
  // ==> Var(sigma) = 0.25 * SigmaSquareVar / (Sigma * Sigma)
  //
  const Float_t rsigmaVar = 0.25 * fRSigmaSquareVar / fRSigmaSquare;

  if (IsHiGainSaturation())
    return TMath::Sqrt(rsigmaVar/fRSigmaSquare + GetConversionHiLoRelVar()) * GetRSigma();
  else
    return TMath::Sqrt(rsigmaVar);

}

// --------------------------------------------------------------------------
//
// Get the reduced Sigma: 
// - If fRSigmaSquare is smaller than 0 (i.e. has not yet been set), return -1.
// - Test bit kHiGainSaturation: 
//   If yes, return square root of fRSigmaSquare, multiplied with fConversionHiLo, 
//   If no , return square root of fRSigmaSquare
//
Float_t MCalibrationChargePix::GetRSigma()  const
{
  if (fRSigmaSquare < 0)
    return -1;

  return TMath::Sqrt(fRSigmaSquare);
  
} 

// --------------------------------------------------------------------------
//
// Get the error of the reduced Sigma: 
// - If fRSigmaSquareVar is smaller than 0 (i.e. has not yet been set), return -1.
// - Calculate the absolute variance of the reduced sigma with the formula:
//   reduced sigma variance = 0.25 * fRSigmaSquareVar / fRSigmaSquare
// - Test bit kHiGainSaturation: 
//   If yes, returns the square root of the quadratic sum of the relative variances of the 
//           reduced sigma and fConversionHiLo, mulitplied with GetRSigma()
//   Else returns the square root of rel. (0.25*fRSigmaSquareVar / fRSigmaSquare)
//
Float_t MCalibrationChargePix::GetRSigmaErr()  const
{

  if (fRSigmaSquareVar < 0)
    return -1;

  //
  // SigmaSquareVar = 4. * Sigma * Sigma * Var(sigma)
  // ==> Var(sigma) = 0.25 * SigmaSquareVar / (Sigma * Sigma)
  //
  return TMath::Sqrt(0.25 * fRSigmaSquareVar / fRSigmaSquare);

}

// --------------------------------------------------------------------------
//
// Get the reduced Sigma per Charge: 
// - If GetRSigma() is smaller or equal 0. (i.e. has not yet been set), return -1.
// - If GetMean() is 0. or -1. (i.e. has not yet been set), return -1.
// - Return GetRSigma() / GetMean() 
//
Float_t MCalibrationChargePix::GetRSigmaPerCharge()  const 
{
  
  const Float_t rsigma = GetRSigma();

  if (rsigma <= 0)
    return -1.;
  

  const Float_t mean   = GetMean();
  
  if (mean == 0. || mean == -1.)
    return -1.;
  
  return rsigma / mean;
} 


// --------------------------------------------------------------------------
//
// Get the error of the reduced Sigma per Charge: 
// - If GetRSigmaRelVar() is smaller or equal 0. (i.e. has not yet been set), return -1.
// - If GetMeanRelVar() is smaller or equal 0. (i.e. has not yet been set), return -1.
// - Return the propagated error of GetRSigmaPerCharge() 
//
Float_t MCalibrationChargePix::GetRSigmaPerChargeErr()  const 
{
  
  const Float_t rsigmarelvar  = GetRSigmaRelVar();

  if (rsigmarelvar <= 0)
    return -1.;
  

  const Float_t meanrelvar   = GetMeanRelVar();
  
  if (meanrelvar <= 0.)
    return -1.;
  
  return TMath::Sqrt(rsigmarelvar + meanrelvar) * GetRSigmaPerCharge();
} 

// --------------------------------------------------------------------------
//
// Get the reduced Sigma Square: 
// - If fRSigmaSquare is smaller than 0 (i.e. has not yet been set), return -1.
// - Test bit kHiGainSaturation: 
//   If yes, return fRSigmaSquare, multiplied with fConversionHiLo^2, 
//   If no , return fRSigmaSquare
//
Float_t MCalibrationChargePix::GetConvertedRSigmaSquare()  const
{
  if (fRSigmaSquare < 0)
    return -1;

  return IsHiGainSaturation() ? fRSigmaSquare*fConversionHiLo*fConversionHiLo : fRSigmaSquare ;
} 

// --------------------------------------------------------------------------
//
// Get the relative variance of the reduced Sigma: 
// - If fRSigmaSquareVar is smaller than 0 (i.e. has not yet been set), return -1.
// - Calculate the relative variance of the reduced sigma squares with the formula:
//   reduced sigma rel. variance = 0.25 * fRSigmaSquareVar / fRSigmaSquare / fRSigmaSquare
// - Test bit kHiGainSaturation: 
//   If yes, returns the sum of the relative variances of the reduced sigma and fConversionHiLo
//   Else returns the relative variance of the reduced sigma
//
Float_t MCalibrationChargePix::GetRSigmaRelVar()  const
{

  if (fRSigmaSquareVar < 0)
    return -1;

  //
  // SigmaSquareVar = 4. * Sigma * Sigma * Var(sigma)
  // ==> Var(sigma) = 0.25 * SigmaSquareVar / (Sigma * Sigma)
  //
  return 0.25 * fRSigmaSquareVar / ( fRSigmaSquare * fRSigmaSquare );

}

// --------------------------------------------------------------------------
//
// Get the error on the number of photo-electrons (F-Factor Method):
// - If fPheFFactorMethodVar is smaller than 0 (i.e. has not yet been set), return -1.
// - Else returns the square root of fPheFFactorMethodVar
//
Float_t MCalibrationChargePix::GetPheFFactorMethodErr()  const
{
  if (fPheFFactorMethodVar < 0.)
    return -1.;
  return TMath::Sqrt(fPheFFactorMethodVar);
}

// --------------------------------------------------------------------------
//
// Get the error on the mean total F-Factor of the signal readout (F-Factor Method):
// - If fMeanFFactorFADC2PhotVar is smaller than 0 (i.e. has not yet been set), return -1.
// - Else returns the square root of fMeanFFactorFADC2PhotVar
//
Float_t MCalibrationChargePix::GetMeanFFactorFADC2PhotErr()  const
{
  if (fMeanFFactorFADC2PhotVar < 0.)
    return -1.;
  return TMath::Sqrt(fMeanFFactorFADC2PhotVar);
}

// --------------------------------------------------------------------------
//
// Get the relative variance on the number of photo-electrons (F-Factor Method):
// - If fPheFFactorMethodVar is smaller than 0 (i.e. has not yet been set), return -1.
// - If fPheFFactorMethod    is 0, return -1.
// - Else returns fPheFFactorMethodVar / fPheFFactorMethod^2
//
Float_t MCalibrationChargePix::GetPheFFactorMethodRelVar()  const
{
  if (fPheFFactorMethodVar < 0.)
    return -1.;
  if (fPheFFactorMethod  == 0.)
    return -1.;

  return fPheFFactorMethodVar / (fPheFFactorMethod * fPheFFactorMethod);
}


// --------------------------------------------------------------------------
//
// Get the error on the mean conversion factor (FFactor  Method):
// - If fMeanConvFADC2PheVar is smaller than 0 (i.e. has not yet been set), return -1.
// - Else returns the square root of fMeanConvFADC2PheVar
//
Float_t MCalibrationChargePix::GetMeanConvFADC2PheErr()  const
{
  if (fMeanConvFADC2PheVar < 0.)
    return -1.;
  return TMath::Sqrt(fMeanConvFADC2PheVar);
}

// --------------------------------------------------------------------------
//
// Test bit kFFactorMethodValid
//
Bool_t MCalibrationChargePix::IsFFactorMethodValid()   const
{ 
  return TESTBIT(fCalibFlags, kFFactorMethodValid);     
}


// ----------------------------------------------------------------------------
// 
// - If fSigma  is smaller than 0 (i.e. has not yet been set), return kFALSE
// - If fPedRms is smaller than 0 (i.e. has not yet been set), return kFALSE
//
// Calculate the reduced sigma of the low-Gain FADC slices:
// - Test bit IsHiGainSaturation() for the Sigma: 
//   If yes, take fLoGainSigma and fLoGainSigmaVar 
//   If no , take fHiGainSigma and fHiGainSigmaVar 
//
// - Test bit IsHiGainSaturation() for the pedRMS: 
//   If yes, take fLoGainPedRmsSquare and fLoGainPedRmsSquareVar
//   If no , take fPedRms and fPedVar
//
// - Calculate the reduced sigma with the formula:
//   fRSigmaSquare = Sigma*Sigma - pedRMS*pedRMS
// 
// - If fRSigmaSquare is smaller than 0, give a warning and return kFALSE
//
// - Calculate the variance of the reduced sigma with the formula:
//   fRSigmaSquareVar = 4.* (sigmaVar*Sigma*Sigma + pedRmsVar*pedRMS*pedRMS)
//
// A back-transformation to the corr. amplification factor of the High-Gain is done 
// in GetRSigma() and GetRSigmaErr()
//
Bool_t MCalibrationChargePix::CalcReducedSigma()
{

  if (GetSigma() < 0.)
    return kFALSE;
  
  if (GetPedRms() < 0.)
    return kFALSE;

  const Float_t sigma           = IsHiGainSaturation() ? fLoGainSigma           : fHiGainSigma   ;
  const Float_t sigmavar        = IsHiGainSaturation() ? fLoGainSigmaVar        : fHiGainSigmaVar;
  const Float_t pedRmsSquare    = IsHiGainSaturation() ? fLoGainPedRmsSquare    : fPedRms*fPedRms;
  const Float_t pedRmsSquareVar = IsHiGainSaturation() ? fLoGainPedRmsSquareVar : 0.25*fPedVar*pedRmsSquare;

  if (IsDebug())
    {
      *fLog << dbginf << "ID: " << GetPixId() 
            << " HiGainSaturation: " << IsHiGainSaturation() 
            << " Sigma: " << sigma 
            << " Var.Sigma: " << sigmavar
            << " PedRmsSquare: " << pedRmsSquare
            << " pedRmsSquareVar: " << pedRmsSquareVar
            << endl;
    }
  
  const Float_t sigmaSquare    =      sigma     * sigma;
  const Float_t sigmaSquareVar = 4. * sigmavar  * sigmaSquare;

  //
  // Calculate the reduced sigmas
  //
  fRSigmaSquare = sigmaSquare - pedRmsSquare;

  if (IsDebug())
    {
      *fLog << dbginf << "ID: " << GetPixId() 
            << " Red.Sigma Square: " << fRSigmaSquare
            << endl;
    }
  
  if (fRSigmaSquare <= 0.)
    {
      if (IsDebug())
        *fLog << warn 
              << "WARNING: Cannot calculate the reduced sigma: smaller than 0 in pixel " 
              << fPixId << endl;
      return kFALSE;
    }

  
  fRSigmaSquareVar = 4. * (sigmaSquareVar + pedRmsSquareVar);

  if (IsDebug())
    {
      *fLog << dbginf << "ID: " << GetPixId() 
            << " Var.Red.Sigma Square: " << fRSigmaSquareVar
            << endl;
    }

  return kTRUE;
}

// ------------------------------------------------------------------
//
// If fRSigmaSquare is smaller than 0 (i.e. has not yet been set),
// return kFALSE
//
// Calculate the number of photo-electrons with the F-Factor method:
// - Test bit IsHiGainSaturation() for the Mean Sum of FADC slices: 
//   If yes, take fLoGainMean and fLoGainMeanVar 
//   If no , take fHiGainMean and fHiGainMeanVar 
//
// - Test bit IsHiGainSaturation() for the pedRMS: 
//   If yes, take fLoGainPedRmsSquare and fLoGainPedRmsSquareVar
//   If no , take fPedRms and fPedVar
//
// - Calculate the number of photo-electrons with the formula:
//   fPheFFactorMethod   = gkFFactor*gkFFactor * Mean * Mean  / fRSigmaSquare
//
// - Calculate the Variance on the photo-electrons with the formula:
//   fPheFFactorMethodVar =  (  4. * gkFFactorErr * gkFFactorErr / ( gkFFactor * gkFFactor ) 
//                            + 4. * Mean Var.   / ( Mean * Mean )
//                            + fRSigmaSquareVar / fRSigmaSquare
//                            ) * fPheFFactorMethod * fPheFFactorMethod
//
// - If fPheFFactorMethod is less than fPheFFactorMethodLimit, 
//   set kFFactorMethodValid to kFALSE and return kFALSE
//
Bool_t MCalibrationChargePix::CalcFFactor()
{

  if (fRSigmaSquare < 0.)
    return kFALSE;
  
  //
  // Square all variables in order to avoid applications of square root
  //
  const Float_t meanSquare          =     GetMean()    * GetMean();
  const Float_t meanSquareRelVar    = 4.* GetMeanRelVar();

  const Float_t ffactorsquare       =     gkFFactor    * gkFFactor;
  const Float_t ffactorsquareRelVar = 4.* GetFFactorRelVar();

  const Float_t rsigmaSquareRelVar  =     fRSigmaSquareVar / fRSigmaSquare / fRSigmaSquare;
  //
  // Calculate the number of phe's from the F-Factor method
  // (independent on Hi Gain or Lo Gain)
  //
  fPheFFactorMethod = ffactorsquare * meanSquare / fRSigmaSquare;

  if (IsDebug())
    {
      *fLog << dbginf << "ID: " << GetPixId() 
            << " F-Factor Square: " << ffactorsquare
            << " Mean Square: " << meanSquare
            << " Red.Sigma Square: " << fRSigmaSquare
            << " Photo-electrons: " << fPheFFactorMethod
            << endl;
    }

  if (fPheFFactorMethod < fPheFFactorMethodLimit)
    return kFALSE;
  
  //
  // Calculate the Error of Nphe
  //
  const Float_t pheRelVar = ffactorsquareRelVar + meanSquareRelVar + rsigmaSquareRelVar;
  fPheFFactorMethodVar    = pheRelVar * fPheFFactorMethod * fPheFFactorMethod;

  if (IsDebug())
    {
      *fLog << dbginf << "ID: " << GetPixId() 
            << " Rel.Var.F-Factor Square: " << ffactorsquareRelVar
            << " Rel.Var. Mean Square: " << meanSquareRelVar
            << " Rel.Var. Red.Sigma Square: " << rsigmaSquareRelVar
            << " Rel.Var. Photo-electrons: " << pheRelVar
            << endl;
    }

  if (fPheFFactorMethodVar < 0. )
    return kFALSE;

  return kTRUE;
}

// ------------------------------------------------------------------
//
// If fPheFFactorMethod is smaller than 0 (i.e. has not yet been set),
// return kFALSE
//
// If GetCovertedMean() is smaller than 0 (i.e. has not yet been set),
// return kFALSE
//
// Calculate fMeanConvFADC2Phe with the following formula:
//
//      fMeanConvFADC2Phe    =  fPheFFactorMethod / GetConvMean();
//
// Calculate the rel. variance of fMeanConvFADC2Phe, taking into account that 
// in the calculation of the number of phe's one mean square has already been used. 
// Now, dividing by another mean, one mean calcels out, one cannot directly propagate
// the errors, but instead havs to take into account this cancellation:
// 
//     convrelvar = ffactorsquareRelVar + GetMeanRelVar() + rsigmaSquareRelVar;
//
// If confrelvar is smaller than 0. or greater than fConvFFactorRelVarLimit, 
// return kFALSE
// 
// Calculate the variance of fMeanConvFADC2Phe with the formula:
//
//    fMeanConvFADC2PheVar =  convrelvar * fMeanConvFADC2Phe * fMeanConvFADC2Phe;
// 
// Set kFFactorMethodValid to kTRUE and 
// return kTRUE
//
Bool_t MCalibrationChargePix::CalcConvFFactor()
{
  
  if (fPheFFactorMethod <= 0.)
    return kFALSE;

  const Float_t convmean = GetConvertedMean();

  if (convmean <= 0.)
    return kFALSE;
  
  fMeanConvFADC2Phe    =  fPheFFactorMethod / convmean;

  if (IsDebug())
    {
      *fLog << dbginf << "ID: " << GetPixId() 
            << " Converted Mean: " << convmean 
            << " Conversion FADC2Phe: " << fMeanConvFADC2Phe
            << endl;
    }

  const Float_t ffactorsquareRelVar = 4.* GetFFactorRelVar();
  const Float_t rsigmaSquareRelVar  = fRSigmaSquareVar / fRSigmaSquare / fRSigmaSquare;
  //
  // In the calculation of the number of phe's one mean square has already been used. 
  // Now, we divide by another mean, so one mean calcels out, we cannot directly propagate
  // the errors, but have to take account of this cancellation:
  // 
  Float_t convrelvar = ffactorsquareRelVar + GetMeanRelVar() + rsigmaSquareRelVar;
  const Float_t limit      = IsHiGainSaturation() ? fConvFFactorRelVarLimit * 4. : fConvFFactorRelVarLimit;

  //
  // Also have to take into account the pixels labelled MBadPixelsPix::kChargeSigmaNotValid which do not 
  // have a fRSigmaSquareVar, calculate their error directly!
  //
  if (fRSigmaSquareVar < 0.)
    convrelvar = GetMeanRelVar() + GetPheFFactorMethodRelVar();

  if (IsDebug())
    {
      *fLog << dbginf << "ID: " << GetPixId() 
            << " Rel.Var.Red.Sigma: " << rsigmaSquareRelVar
            << " Rel.Var.Mean: " << GetMeanRelVar()
            << " Rel.Var.F-Factor: " << ffactorsquareRelVar
            << " Rel.Var.Conversion FADC2Phe: " << convrelvar
            << endl;
    }

  if (convrelvar > limit || convrelvar < 0.)
    {
      *fLog << warn << GetDescriptor() << ":  Conv. F-Factor Method Rel. Var.: " 
            << Form("%4.3f out of limits: [0,%3.2f] in pixel:%4i",convrelvar,limit,fPixId) << endl;
      return kFALSE;
    }
  
  fMeanConvFADC2PheVar =  convrelvar * fMeanConvFADC2Phe * fMeanConvFADC2Phe;
  
  SetFFactorMethodValid(kTRUE);
  return kTRUE;
}

// ----------------------------------------------------------------------------------
//
// If photflux is smaller or equal 0, return kFALSE
//
// Calculate the total F-Factor with the formula:
//   fMeanFFactorFADC2Phot = Sqrt ( fRSigmaSquare ) / GetMean()  * sqrt(nphotons)
//
// Calculate the error of the total F-Factor
//
Bool_t MCalibrationChargePix::CalcMeanFFactor( const Float_t nphotons, const Float_t nphotonsrelvar )
{


  if (IsDebug())
    {
      *fLog << dbginf << "ID: " << GetPixId() 
            << " Number photons: " << nphotons
            << " Rel.Var.Number photons: " << nphotonsrelvar
            << " Red.Sigma Square: " << fRSigmaSquare
            << " Mean: " << GetMean()
            << endl;
    }


  if (nphotons <= 0.)
    {
      *fLog << warn << GetDescriptor() << ": Assumed photon flux is smaller or equal 0." << endl;
      return kFALSE;
    }

  if (nphotonsrelvar < 0.)
    {
      *fLog << warn << GetDescriptor() << ": Assumed photon flux variance is smaller than 0." << endl;
      return kFALSE;
    }

  fMeanFFactorFADC2Phot =  TMath::Sqrt(fRSigmaSquare * nphotons) / GetMean() ;
  
  if (IsDebug())
    {
      *fLog << dbginf << "ID: " << GetPixId() 
            << " F-Factor FADC2Phot: " << fMeanFFactorFADC2Phot
            << endl;
    }

  if (fMeanFFactorFADC2Phot < 0.)
    {
      *fLog << warn << GetDescriptor() << ": F-Factor photons to FADC counts smaller than 0." << endl;
      return kFALSE;
    }
  
  const Float_t ffactorrelvar = 0.25 * fRSigmaSquareVar / ( fRSigmaSquare * fRSigmaSquare) 
                              + GetMeanRelVar() 
                              + 0.25 * nphotonsrelvar;
  
  fMeanFFactorFADC2PhotVar    = ffactorrelvar * fMeanFFactorFADC2Phot * fMeanFFactorFADC2Phot;

  if (IsDebug())
    {
      *fLog << dbginf << "ID: " << GetPixId() 
            << " Rel.Var.Red.Sigma: " << 0.25 * fRSigmaSquareVar / ( fRSigmaSquare * fRSigmaSquare) 
            << " Rel.Var.Mean: " << GetMeanRelVar()
            << " Rel.Var.photons: " << 0.25 * nphotonsrelvar
            << " Rel.Var.F-Factor FADC2Phot: " << ffactorrelvar
            << endl;
    }

  return kTRUE;
}


// ----------------------------------------------------------------------------
// 
// - If fPed    is smaller than 0 (i.e. has not yet been set), return.
// - If fPedVar is smaller than 0 (i.e. has not yet been set), return.
//
// Calculate the electronic pedestal RMS with the formula:
//  - elec. pedestal = gkElectronicPedRms * sqrt(logainsamples)
// 
// Calculate the night sky background ped. RMS contribution ("NSB") in the high-gain 
// from the high gain Pedestal RMS with the formula:
// - HiGain NSB square      = fPedRms * fPedRms - elec.ped.* elec.ped.
// - Var(HiGain NSB square) = fPedVar * fPedRms * fPedRms + 4.*elecPedRmsVar * elec.ped.* elec.ped.
//
// If HiGain NSB square is smaller than 0., set it to zero. (but not the error!)
//
// Convert the NSB ped. RMS contribution to the low-gain with the formula:
// - LoGain NSB square      = - HiGain NSB square / (fConversionHiLo*fConversionHiLo)
// - Var(LoGain NSB square) = ( Var(HiGain NSB square) / (HiGain NSB square * HiGain NSB square)
//                              + GetConversionHiLoRelVar()   
//                            ) * LoGain NSB square * LoGain NSB square
//
// - Low Gain Ped RMS Square       = LoGain NSB square      + elec.ped. square
//   Var (Low Gain Ped RMS Square) = Var(LoGain NSB square) + Var(elec.ped. square)
//
void MCalibrationChargePix::CalcLoGainPedestal(Float_t logainsamples, const Int_t aidx)
{

  if (fPedRms < 0.)
    return;

  if (fPedVar < 0.)
    return;

  const Float_t elecPedRms     = (aidx == 0 ? gkElectronicPedRmsInner : gkElectronicPedRmsOuter ) 
                               * TMath::Sqrt(logainsamples) / fConversionHiLo;
  const Float_t elecPedRmsVar  = ( GetElectronicPedRmsRelVar(aidx) + GetConversionHiLoRelVar() )
                               * elecPedRms * elecPedRms;
  
  Float_t pedRmsSquare      = fPedRms * fPedRms;
  Float_t pedRmsSquareVar   = fPedVar * pedRmsSquare; // fPedRmsErr = fPedErr/2.
  
  //
  // 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    =    elecPedRms    * elecPedRms;
  const Float_t elecRmsSquareVar = 4.*elecPedRmsVar * elecRmsSquare;
  
  if (IsDebug())
    {
      *fLog << dbginf << "ID: " << GetPixId() 
            << " Ped.Rms Square:  " << pedRmsSquare 
            << " Elec.Rms Square: " << elecRmsSquare 
            << " Ped.Rms.Square Var.: " << pedRmsSquareVar 
            << " Elec.Rms Square Var.: " << elecRmsSquareVar 
            << endl;
    }
  

  Float_t higainNsbSquare        =  pedRmsSquare    - elecRmsSquare;
  Float_t higainNsbSquareRelVar  = (pedRmsSquareVar + elecRmsSquareVar);


  if (higainNsbSquare < 0.001)
    higainNsbSquare = 0.;
  else
    higainNsbSquareRelVar  /= (higainNsbSquare * higainNsbSquare) ;
  
  if (IsDebug())
    {
      *fLog << dbginf << "ID: " << GetPixId() 
            << " HiGain NSB Square:  " << higainNsbSquare
            << " Rel.Var.HiGain NSB Square:  " << higainNsbSquareRelVar
            << endl;
    }
  
  //
  // 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 conversionSquareRelVar  = 4.* GetConversionHiLoRelVar();

  const Float_t logainNsbSquare         =   higainNsbSquare       / conversionSquare;
  const Float_t logainNsbSquareVar      = ( higainNsbSquareRelVar + conversionSquareRelVar )
                                	  * logainNsbSquare * logainNsbSquare;
    
  fLoGainPedRmsSquare    = logainNsbSquare    + elecRmsSquare;
  fLoGainPedRmsSquareVar = logainNsbSquareVar + elecRmsSquareVar;

  if (IsDebug())
    {
      *fLog << dbginf << "ID: " << GetPixId() 
            << " LoGain Ped Rms Square:  " << fLoGainPedRmsSquare
            << " Var.Ped Rms Square:  " << fLoGainPedRmsSquareVar
            << endl;
    }
  
  
}