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

/* ======================================================================== *\
!
! *
! * 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
//
//
// ClassVersion 3:
//  + fNumSaturated
//
// Class Version 4:
//  +  Float_t fConversionHiLoSigma;             // Sigma of conversion factor betw. Hi and Lo Gain
//  -  Float_t fMeanConvFADC2PheVar;             // Variance conversion factor (F-factor method)
//  +  Float_t fMeanConvFADC2PheStatVar;         // Variance conversion factor, only stat. error
//  +  Float_t fMeanConvFADC2PheSystVar;         // Variance conversion factor, only syst. error
//  -  Float_t fMeanFFactorFADC2PhotVar;         // Variance mean F-Factor photons (F-factor method)
//  +  Float_t fMeanFFactorFADC2PhotVar;         // Variance mean F-Factor photons, only stat. error
//  -  Float_t fPheFFactorMethod;                // Number Phe's calculated (F-factor method)
//  -  Float_t fPheFFactorMethodVar;             // Variance number of Phe's (F-factor method)
//  +  Float_t fPheFFactorMethod;                // Number Phe's calculated  with F-factor method)
//  +  Float_t fPheFFactorMethodStatVar;         // Variance number of Phe's, only stat. error
//  +  Float_t fPheFFactorMethodSystVar;         // Variance number of Phe's, only syst. error
//
/////////////////////////////////////////////////////////////////////////////
#include "MCalibrationChargePix.h"

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

#include "MBadPixelsPix.h"

ClassImp(MCalibrationChargePix);

using namespace std;

const Float_t MCalibrationChargePix::gkElectronicPedRms    = 1.0;
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      = 0.05;
const Float_t MCalibrationChargePix::fgConversionHiLoSigma    = 2.5;
const Float_t MCalibrationChargePix::fgPheFFactorMethodLimit  = 1.;
const Float_t MCalibrationChargePix::fgConvFFactorRelErrLimit = 0.85;

// --------------------------------------------------------------------------
//
// 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), fNumSaturated(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.;
  fPheFFactorMethodStatVar          =  -1.;
  fPheFFactorMethodSystVar          =  -1.;

  fMeanConvFADC2Phe                 =  -1.;
  fMeanConvFADC2PheStatVar          =  -1.;
  fMeanConvFADC2PheSystVar          =  -1.;
  fMeanFFactorFADC2Phot             =  -1.;
  fMeanFFactorFADC2PhotVar          =  -1.;  

  fNumSaturated                     =   0;

  MCalibrationPix::Clear();
}

// -----------------------------------------------------
//
void MCalibrationChargePix::Copy(TObject& object) const
{
  MCalibrationChargePix &pix = (MCalibrationChargePix&)object;
  
  MCalibrationPix::Copy(pix);

  //
  // Copy the data members
  //
  pix.fAbsTimeMean                = fAbsTimeMean              ;
  pix.fAbsTimeRms                 = fAbsTimeRms               ;
  pix.fCalibFlags                 = fCalibFlags               ;
  pix.fConversionHiLo             = fConversionHiLo           ;
  pix.fConversionHiLoVar          = fConversionHiLoVar        ;
  pix.fConversionHiLoSigma        = fConversionHiLoSigma      ;
  pix.fConvFFactorRelVarLimit     = fConvFFactorRelVarLimit   ;
  pix.fLoGainPedRmsSquare         = fLoGainPedRmsSquare       ;
  pix.fLoGainPedRmsSquareVar      = fLoGainPedRmsSquareVar    ;
  pix.fMeanConvFADC2Phe           = fMeanConvFADC2Phe         ;
  pix.fMeanConvFADC2PheStatVar    = fMeanConvFADC2PheStatVar  ;
  pix.fMeanConvFADC2PheSystVar    = fMeanConvFADC2PheSystVar  ;
  pix.fMeanFFactorFADC2Phot       = fMeanFFactorFADC2Phot     ;
  pix.fMeanFFactorFADC2PhotVar    = fMeanFFactorFADC2PhotVar  ;
  pix.fPed                        = fPed                      ;
  pix.fPedVar                     = fPedVar                   ;
  pix.fPedRms                     = fPedRms                     ;
  pix.fPedRmsVar                  = fPedRmsVar                  ;
  pix.fPheFFactorMethod           = fPheFFactorMethod           ;
  pix.fPheFFactorMethodStatVar    = fPheFFactorMethodStatVar    ;
  pix.fPheFFactorMethodSystVar    = fPheFFactorMethodSystVar    ;
  pix.fPheFFactorMethodLimit      = fPheFFactorMethodLimit      ;
  pix.fRSigmaSquare               = fRSigmaSquare               ;
  pix.fRSigmaSquareVar            = fRSigmaSquareVar            ;
  pix.fNumSaturated               = fNumSaturated               ;
}                                  


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


// Inline Functions:
// -----------------
//
// GetConversionHiLoErr:
//   Get the conversion Error Hi-Gain to Low-Gain:
//   If fConversionHiLoVar is smaller than 0 (i.e. has not yet
//   been set), return -1.
//  
// GetPedRms(): Get the pedestals RMS: Test bit kHiGainSaturation:
//  If yes, return square root of fLoGainPedRmsSquare (if greater than 0,
//  otherwise -1.), If no,  return fPedRms
//
// GetConvertedMean(): Get the Low Gain Mean Charge converted to High Gain
//  amplification: Returns fLoGainMean multiplied with fConversionHiLo if
//  IsHiGainSaturation(), else return fHiGainMean
//
// GetConvertedMeanErr(): 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()
//
// GetConvertedSigma(): Get the Low Gain Sigma converted to High Gain
//  amplification: Returns fLoGainSigma multiplied with fConversionHiLo
//  if IsHiGainSaturation() else return fHiGainSigma
//
// GetConvertedSigmaErr(): 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()
//
// GetConvertedRSigma(): 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
//
// GetRSigma(): Get the reduced Sigma:
//  If fRSigmaSquare is smaller than 0 (i.e. has not yet been set), return -1.
//
// GetConvertedRSigmaSquare(): 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
//
// GetPheFFactorMethodErr(): 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
//
// GetMeanFFactorFADC2PhotErr(): 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
//
// GetPheFFactorMethodRelVar(): 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
//
// GetMeanConvFADC2PheErr(): 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

// --------------------------------------------------------------------------
//
// 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
// - returns rel. variance of gkElectronicPedRms
//
const Float_t MCalibrationChargePix::GetElectronicPedRmsRelVar() const 
{
  return gkElectronicPedRmsErr * gkElectronicPedRmsErr / gkElectronicPedRms / gkElectronicPedRms;
}
 

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

}



// --------------------------------------------------------------------------
//
// 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(const Float_t extractorres)
{

  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 extractorresVar = extractorres * extractorres;

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

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

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

    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  = meanSquareRelVar + rsigmaSquareRelVar;
  fPheFFactorMethodStatVar = pheRelVar * fPheFFactorMethod * fPheFFactorMethod;
  fPheFFactorMethodSystVar = ffactorsquareRelVar * 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 (fPheFFactorMethodStatVar < 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 = GetMeanRelVar() + rsigmaSquareRelVar;
  if (IsHiGainSaturation())
      convrelvar += GetConversionHiLoRelVar();

  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;
    }
  
  fMeanConvFADC2PheStatVar = convrelvar * fMeanConvFADC2Phe  * fMeanConvFADC2Phe;
  fMeanConvFADC2PheSystVar = ffactorsquareRelVar * 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)
{

  if (fPedRms < 0.)
    return;

  if (fPedVar < 0.)
    return;

  const Float_t elecPedRms     = gkElectronicPedRms * TMath::Sqrt(logainsamples);
  const Float_t elecPedRmsVar  = GetElectronicPedRmsRelVar() * elecPedRms * elecPedRms;
  
  Float_t pedRmsSquare         = fPedRms * fPedRms;
  Float_t pedRmsSquareVar      = fPedVar * pedRmsSquare; // fPedRmsErr = fPedErr/2.
  
  //
  // We do not know the Low Gain Pedestal RMS, so we have to retrieve it 
  // from the High Gain:  
  //
  // 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 higainNsbSquareVar  = (pedRmsSquareVar + elecRmsSquareVar);

  if (higainNsbSquare < 0.001)
    higainNsbSquare = 0.;
  
  if (IsDebug())
    {
      *fLog << dbginf << "ID: " << GetPixId() 
            << " HiGain NSB Square:  " << higainNsbSquare
            << " Var.HiGain NSB Square:  " << higainNsbSquareVar
            << 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;
  //
  // Calculation of variance of: c = a/b
  //                    Delta(c)^2 = ( Delta(a)^2 + a^2/b^2*(Delta(b)^2 )  / b^2
  //
  const Float_t logainNsbSquareVar      = ( higainNsbSquareVar 
                                         + conversionSquareRelVar * higainNsbSquare * higainNsbSquare ) 
                                          / conversionSquare / conversionSquare;
    
  fLoGainPedRmsSquare    = logainNsbSquare    + elecRmsSquare;
  fLoGainPedRmsSquareVar = logainNsbSquareVar + elecRmsSquareVar;

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