source: trunk/MagicSoft/Mars/mcalib/MCalibrationChargePix.cc@ 3677

/* ======================================================================== *\
!
! *
! * 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 CalcFFactorMethod() )
//
// 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                         
// - fMeanConversionFFactorMethod, fMeanConversionBlindPixelMethod, 
//   fMeanConversionPINDiodeMethod and fMeanConversionCombinedMethod   
// - fConversionFFactorMethodVar, fConversionBlindPixelMethodVar 
//   fConversionPINDiodeMethodVar and fConversionCombinedMethodVar
// - fSigmaConversionFFactorMethod, fSigmaConversionBlindPixelMethod
//   fSigmaConversionPINDiodeMethod and fSigmaConversionCombinedMethod  
// - fTotalFFactorFFactorMethod, fTotalFFactorBlindPixelMethod   
//   fTotalFFactorPINDiodeMethod and fTotalFFactorCombinedMethod     
// - fTotalFFactorFFactorMethodVar, fTotalFFactorBlindPixelMethodVar, 
//   fTotalFFactorPINDiodeMethodVar and fTotalFFactorCombinedMethodVar  
//
// 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: GetConvertedLoGainMean() or GetConvertedLoGainSigma()
// 
// 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::gkElectronicPedRms         = 1.5;
const Float_t MCalibrationChargePix::gkElectronicPedRmsErr      = 0.3;
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    = 5.;
// --------------------------------------------------------------------------
//
// Default Constructor: 
//
// Sets:
// - fCalibFlags to 0
// - fConversionHiLo to fgConversionHiLo
// - fConversionHiLoVar to square of fgConversionHiLoErr
// - fPheFFactorMethodLimit to fgPheFFactorMethodLimit
// 
// 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();
  
  Clear();
}

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

  SetBlindPixelMethodValid  ( kFALSE );
  SetFFactorMethodValid     ( kFALSE );
  SetPINDiodeMethodValid    ( kFALSE );
  SetCombinedMethodValid    ( kFALSE );

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

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

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

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

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

  fConversionFFactorMethodVar       =  -1.;
  fConversionBlindPixelMethodVar    =  -1.;
  fConversionPINDiodeMethodVar      =  -1.;
  fConversionCombinedMethodVar      =  -1.;

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

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

  MCalibrationPix::Clear();
}


// --------------------------------------------------------------------------
//
// Set conversion factors Blind Pixel Method from outside (only for MC)
//
void MCalibrationChargePix::SetConversionBlindPixelMethod(Float_t c, Float_t err, Float_t sig)
{
  fMeanConversionBlindPixelMethod  = c;
  fConversionBlindPixelMethodVar   = err*err;
  fSigmaConversionBlindPixelMethod = sig;
}


// --------------------------------------------------------------------------
//
// Set conversion factors Combined Method from outside (only for MC)
//
void MCalibrationChargePix::SetConversionCombinedMethod(Float_t c, Float_t err, Float_t sig)
{
  fMeanConversionCombinedMethod  = c;
  fConversionCombinedMethodVar   = err*err;
  fSigmaConversionCombinedMethod = sig;
}


// --------------------------------------------------------------------------
//
// Set conversion factors F-Factor Method from outside (only for MC)
//
void MCalibrationChargePix::SetConversionFFactorMethod(Float_t c, Float_t err, Float_t sig)
{
  fMeanConversionFFactorMethod  = c;
  fConversionFFactorMethodVar   = err*err;
  fSigmaConversionFFactorMethod = sig;
}

// --------------------------------------------------------------------------
//
// Set conversion factors PIN Diode Method from outside (only for MC)
//
void MCalibrationChargePix::SetConversionPINDiodeMethod(Float_t c, Float_t err, Float_t sig)
{
  fMeanConversionPINDiodeMethod  = c ;
  fConversionPINDiodeMethodVar   = err*err;
  fSigmaConversionPINDiodeMethod = sig;
}

// --------------------------------------------------------------------------
//
// Set Blind Pixel Method Validity Bit from outside 
//
void MCalibrationChargePix::SetBlindPixelMethodValid(const Bool_t b )
{ 
  b ?  SETBIT(fCalibFlags, kBlindPixelMethodValid) : CLRBIT(fCalibFlags, kBlindPixelMethodValid); 
}    

// --------------------------------------------------------------------------
//
// Set Combined Method Validity Bit from outside 
//
void MCalibrationChargePix::SetCombinedMethodValid(const Bool_t b )  
{ 
  b ?  SETBIT(fCalibFlags, kCombinedMethodValid) : CLRBIT(fCalibFlags, kCombinedMethodValid); 
}

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

// --------------------------------------------------------------------------
//
// Set PIN Diode Method Validity Bit from outside 
//
void MCalibrationChargePix::SetPINDiodeMethodValid(const Bool_t b )  
{ 
  b ?  SETBIT(fCalibFlags, kPINDiodeMethodValid) : CLRBIT(fCalibFlags, kPINDiodeMethodValid); 
}

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

// -------------------------------------------------------------------------------
//
// 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 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 Error of the Mean pedestals: 
// Returns square root of fPedVar
//
Float_t  MCalibrationChargePix::GetPedErr()  const
{
  return  TMath::Sqrt(fPedVar);
}

// --------------------------------------------------------------------------
//
// 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 converted to High Gain amplification: 
// Returns fLoGainMean multiplied with fConversionHiLo
//
Float_t MCalibrationChargePix::GetConvertedLoGainMean()  const 
{
  return fLoGainMean * fConversionHiLo;
}

// --------------------------------------------------------------------------
//
// 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 GetConvertedLoGainMean()
//
Float_t MCalibrationChargePix::GetConvertedLoGainMeanErr()  const
{

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

  return TMath::Sqrt(logainrelvar + GetConversionHiLoRelVar()) * GetConvertedLoGainMean();
}

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

// --------------------------------------------------------------------------
//
// Get the Error of the converted Low Gain Sigma: 
//
// Returns -1 if the variable fLoGainSigma or fLoGainSigmaVar are smaller than 0.
//
// Returns the square root of the quadratic sum of the relative variances of 
// the fLoGainSigma and fConversionHiLo, mulitplied with GetConvertedLoGainSigma()
//
Float_t MCalibrationChargePix::GetConvertedLoGainSigmaErr()  const
{

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

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

  const Float_t sigmaRelVar     =  fLoGainSigmaVar
                                /( fLoGainSigma * fLoGainSigma );

  return TMath::Sqrt(sigmaRelVar+GetConversionHiLoRelVar()) * GetConvertedLoGainSigma();
}

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

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

// --------------------------------------------------------------------------
//
// 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)
  //
  const Float_t rsigmaVar = 0.25 * fRSigmaSquareVar / fRSigmaSquare;

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

}

// --------------------------------------------------------------------------
//
// 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)
  //
  const Float_t rsigmaRelVar = 0.25 * fRSigmaSquareVar / ( fRSigmaSquare * fRSigmaSquare );

  if (IsHiGainSaturation())
    return rsigmaRelVar + GetConversionHiLoRelVar();
  else
    return rsigmaRelVar;
}

// --------------------------------------------------------------------------
//
// 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 conversion factor (Combined Method):
// - If fConversionCombinedMethodVar is smaller than 0 (i.e. has not yet been set), return -1.
// - Else returns the square root of fConversionCombinedMethodVar
//
Float_t MCalibrationChargePix::GetConversionCombinedMethodErr()  const
{
  if (fConversionCombinedMethodVar < 0.)
    return -1.;
  return TMath::Sqrt(fConversionCombinedMethodVar);
}

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

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

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

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

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

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

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

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

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

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

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

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

  const Float_t sigmaSquare    = sigma     * sigma;
  const Float_t sigmaSquareVar = sigmavar  * sigmaSquare;

  //
  // Calculate the reduced sigmas
  //
  fRSigmaSquare = sigmaSquare - pedRmsSquare;
  if (fRSigmaSquare <= 0.)
    {
      *fLog << warn 
            << "WARNING: Cannot calculate the reduced sigma: smaller than 0 in pixel " 
            << fPixId << endl;
      return kFALSE;
    }

  fRSigmaSquareVar = 4. * (sigmaSquareVar + pedRmsSquareVar);

  return kTRUE;
}

// ------------------------------------------------------------------
//
// If fRSigmaSquare is smaller than 0 (i.e. has not yet been set),
// set kFFactorMethodValid to kFALSE and 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
//   else: Set kFFactorMethodValid to kTRUE and return kTRUE
//
Bool_t MCalibrationChargePix::CalcFFactorMethod()
{

  if (fRSigmaSquare < 0.)
    {
      SetFFactorMethodValid(kFALSE);
      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;

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

  if (fPheFFactorMethod < fPheFFactorMethodLimit)
    {
      SetFFactorMethodValid(kFALSE);
      return kFALSE;
    }
  
  //
  // Calculate the Error of Nphe
  //
  fPheFFactorMethodVar =  (ffactorsquareRelVar + meanSquareRelVar + rsigmaSquareRelVar) 
                         * fPheFFactorMethod * fPheFFactorMethod;

  SetFFactorMethodValid(kTRUE);
  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)
{

  if (fPedRms < 0.)
    return;

  if (fPedVar < 0.)
    return;

  const Float_t elecPedRms     = gkElectronicPedRms    * TMath::Sqrt(logainsamples);
  const Float_t elecPedRmsVar  = gkElectronicPedRmsErr * gkElectronicPedRmsErr * logainsamples;
  
  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;
  
  Float_t higainNsbSquare        =  pedRmsSquare    - elecRmsSquare;
  Float_t higainNsbSquareRelVar  = (pedRmsSquareVar + elecRmsSquareVar)
                                 / (higainNsbSquare * higainNsbSquare) ;
  
  if (higainNsbSquare < 0.)
    higainNsbSquare = 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 conversionSquareRelVar  = 4.* GetConversionHiLoRelVar();

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