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

#include "MLog.h"
#include "MLogManip.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::gkAverageQE                = 0.18;     
const Float_t MCalibrationChargePix::gkAverageQEErr             = 0.02;  
const Float_t MCalibrationChargePix::gkConversionHiLo           = 10.;
const Float_t MCalibrationChargePix::gkConversionHiLoErr        = 2.5;

const Float_t MCalibrationChargePix::fgChargeLimit              = 3.;
const Float_t MCalibrationChargePix::fgChargeErrLimit           = 0.;
const Float_t MCalibrationChargePix::fgChargeRelErrLimit        = 1.;
const Float_t MCalibrationChargePix::fgTimeLimit                = 1.5;
const Float_t MCalibrationChargePix::fgTimeErrLimit             = 3.;
const Float_t MCalibrationChargePix::fgConvFFactorRelErrLimit   = 0.25;
// --------------------------------------------------------------------------
//
// Default Constructor: 
//
MCalibrationChargePix::MCalibrationChargePix(const char *name, const char *title)
    : fPixId(-1),
      fFlags(0)
{

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

  Clear();

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

  SetAverageQE();
  SetChargeLimit();
  SetChargeErrLimit();  
  
  SetChargeRelErrLimit();
  SetTimeLimit();
  SetTimeErrLimit();
  SetConvFFactorRelErrLimit();
}

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

  SetHiGainSaturation       ( kFALSE );
  SetLoGainSaturation       ( kFALSE );
  SetHiGainFitted           ( kFALSE );
  SetLoGainFitted           ( kFALSE );
  SetHiGainOscillating      ( kFALSE );
  SetLoGainOscillating      ( kFALSE );
  SetExcluded               ( kFALSE );
  SetChargeValid            ( kFALSE );
  SetBlindPixelMethodValid  ( kFALSE );
  SetFFactorMethodValid     ( kFALSE );
  SetPINDiodeMethodValid    ( kFALSE );
  SetCombinedMethodValid    ( kFALSE );

  SetMeanTimeInFirstBin     ( kFALSE );
  SetMeanTimeInLastBin      ( kFALSE );

  fHiGainMeanCharge                 =  -1.;
  fHiGainMeanChargeErr              =  -1.;
  fHiGainSigmaCharge                =  -1.;
  fHiGainSigmaChargeErr             =  -1.;
  fHiGainChargeProb                 =  -1.;

  fLoGainMeanCharge                 =  -1.;
  fLoGainMeanChargeErr              =  -1.;
  fLoGainSigmaCharge                =  -1.;
  fLoGainSigmaChargeErr             =  -1.;
  fLoGainChargeProb                 =  -1.;

  fRSigmaCharge                     =  -1.;
  fRSigmaChargeErr                  =  -1.;
  
  fHiGainNumPickup                  =  -1;
  fLoGainNumPickup                  =  -1;

  fNumLoGainSamples                 =  -1.;

  fPed                              =  -1.;
  fPedRms                           =  -1.;
  fPedErr                           =  -1.;

  fLoGainPedRms                     =  -1.;
  fLoGainPedRmsErr                  =  -1.;

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

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

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

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

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

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

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

}


void MCalibrationChargePix::DefinePixId(Int_t i)
{
  fPixId = i;
}


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

  fPed       = ped;    
  fPedRms    = pedrms;
  fPedErr    = pederr;
}
  

void  MCalibrationChargePix::SetMeanCharge( const Float_t f ) 
{
    if (IsHiGainSaturation())
	fLoGainMeanCharge = f;
    else
	fHiGainMeanCharge = f;
}

void  MCalibrationChargePix::SetMeanChargeErr( const Float_t f ) 
{
    if (IsHiGainSaturation())
	fLoGainMeanChargeErr = f;
    else
	fHiGainMeanChargeErr = f;

}

void  MCalibrationChargePix::SetSigmaCharge( const Float_t f )  
{
    if (IsHiGainSaturation())
	fLoGainSigmaCharge = f;
    else
	fHiGainSigmaCharge = f;
}


void  MCalibrationChargePix::SetSigmaChargeErr( const Float_t f )   
{
    if (IsHiGainSaturation())
	fLoGainSigmaChargeErr = f;
    else
	fHiGainSigmaChargeErr = f;

}



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

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


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

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

// --------------------------------------------------------------------------
//
// Set the Hi Gain Saturation Bit from outside
//
void MCalibrationChargePix::SetHiGainSaturation(Bool_t b)
{
    b ?  SETBIT(fFlags, kHiGainSaturation) : CLRBIT(fFlags, kHiGainSaturation); 
}

// --------------------------------------------------------------------------
//
// Set the Lo Gain Saturation Bit from outside
//
void MCalibrationChargePix::SetLoGainSaturation(Bool_t b)
{
    b ?  SETBIT(fFlags, kLoGainSaturation) : CLRBIT(fFlags, kLoGainSaturation); 
}

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

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

// --------------------------------------------------------------------------
//
// Set the Fitted Bit from outside 
//
void MCalibrationChargePix::SetHiGainFitted(Bool_t b )
{ 
  b ?  SETBIT(fFlags, kHiGainFitted) : CLRBIT(fFlags, kHiGainFitted); 
}
    
// --------------------------------------------------------------------------
//
// Set the Fitted Bit from outside 
//
void MCalibrationChargePix::SetLoGainFitted(const Bool_t b )
{ 
  b ?  SETBIT(fFlags, kLoGainFitted) : CLRBIT(fFlags, kLoGainFitted); 
}
    
// --------------------------------------------------------------------------
//
// Set the Ocillating Bit from outside 
//
void MCalibrationChargePix::SetHiGainOscillating(const Bool_t b )
{ 
  b ?  SETBIT(fFlags, kHiGainOscillating) : CLRBIT(fFlags, kHiGainOscillating); 
}
    
// --------------------------------------------------------------------------
//
// Set the Ocillating Bit from outside 
//
void MCalibrationChargePix::SetLoGainOscillating(const Bool_t b )
{ 
  b ?  SETBIT(fFlags, kLoGainOscillating) : CLRBIT(fFlags, kLoGainOscillating); 
}
    
// --------------------------------------------------------------------------
//
// Set the Excluded Bit from outside 
//
void MCalibrationChargePix::SetBlindPixelMethodValid(const Bool_t b )
{ 
  b ?  SETBIT(fFlags, kBlindPixelMethodValid) : CLRBIT(fFlags, kBlindPixelMethodValid); 
}    

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

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

// --------------------------------------------------------------------------
//
// Set the Excluded Bit from outside 
//
void MCalibrationChargePix::SetCombinedMethodValid(const Bool_t b )  
{ 
  b ?  SETBIT(fFlags, kCombinedMethodValid) : CLRBIT(fFlags, kCombinedMethodValid); 
}

void MCalibrationChargePix::SetAbsTimeBordersHiGain(const Byte_t f, const Byte_t l)
{

  fTimeFirstHiGain = f;
  fTimeLastHiGain  = l;
  
}

void MCalibrationChargePix::SetAbsTimeBordersLoGain(const Byte_t f, const Byte_t l)
{

  fTimeFirstLoGain = f;
  fTimeLastLoGain  = l;
  
}

void MCalibrationChargePix::SetMeanTimeInFirstBin(const Bool_t b)
{
  b ? SETBIT(fFlags,kMeanTimeInFirstBin) : CLRBIT(fFlags,kMeanTimeInFirstBin);
}

void MCalibrationChargePix::SetMeanTimeInLastBin(const Bool_t b)
{
  b ? SETBIT(fFlags,kMeanTimeInLastBin) : CLRBIT(fFlags,kMeanTimeInLastBin);
}

Float_t  MCalibrationChargePix::GetMeanCharge()   const 
{
    return IsHiGainSaturation() ? fLoGainMeanCharge : fHiGainMeanCharge ;
}

Float_t  MCalibrationChargePix::GetMeanChargeErr()   const 
{
    return IsHiGainSaturation() ? fLoGainMeanChargeErr : fHiGainMeanChargeErr ;
}

Float_t  MCalibrationChargePix::GetChargeProb()   const 
{
    return IsHiGainSaturation() ? fLoGainChargeProb : fHiGainChargeProb ;
}

Float_t  MCalibrationChargePix::GetSigmaCharge()   const 
{
    return IsHiGainSaturation() ? fLoGainSigmaCharge : fHiGainSigmaCharge ;
}

Float_t  MCalibrationChargePix::GetSigmaChargeErr()   const 
{
    return IsHiGainSaturation() ? fLoGainSigmaChargeErr : fHiGainSigmaChargeErr ;
}

Bool_t  MCalibrationChargePix::IsFitted()  const 
{
  return  IsHiGainSaturation() ? IsLoGainFitted() :  IsHiGainFitted();
}

Bool_t  MCalibrationChargePix::IsOscillating()  const 
{
  return  IsHiGainSaturation() ? IsLoGainOscillating() :  IsHiGainOscillating();
}

Bool_t MCalibrationChargePix::IsMeanTimeInFirstBin() const
{
  return TESTBIT(fFlags,kMeanTimeInFirstBin);
}

Bool_t MCalibrationChargePix::IsMeanTimeInLastBin() const
{
  return TESTBIT(fFlags,kMeanTimeInLastBin);
}



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

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

Bool_t MCalibrationChargePix::IsLoGainSaturation()    const
{ 
   return TESTBIT(fFlags,kLoGainSaturation);  
}

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

Bool_t MCalibrationChargePix::IsHiGainFitted()        const 
{ 
  return TESTBIT(fFlags, kHiGainFitted);    
}

Bool_t MCalibrationChargePix::IsLoGainFitted()        const 
{ 
  return TESTBIT(fFlags, kLoGainFitted);    
}

Bool_t MCalibrationChargePix::IsHiGainOscillating()    const
{ 
  return TESTBIT(fFlags, kHiGainOscillating);
}

Bool_t MCalibrationChargePix::IsLoGainOscillating()    const
{ 
  return TESTBIT(fFlags, kLoGainOscillating);
}

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

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

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

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


//
// The check return kTRUE if:
//
// 1) Pixel has a fitted charge greater than fChargeLimit*PedRMS
// 2) Pixel has a fit error greater than fChargeErrLimit
// 3) Pixel has a fitted charge greater its fChargeRelErrLimit times its charge error
// 4) Pixel has a charge sigma bigger than its Pedestal RMS
// 
Bool_t MCalibrationChargePix::CheckChargeValidity()
{
 
  if (!IsFitted())
    return kFALSE;

  if (GetMeanCharge() < fChargeLimit*GetPedRms())
    {
      *fLog << warn << "WARNING: Fitted Charge is smaller than "
            << fChargeLimit << " Pedestal RMS in Pixel  " << fPixId << endl;
      return kFALSE;
    }
  
  if (GetMeanChargeErr() < fChargeErrLimit) 
    {
      *fLog << warn << "WARNING: Error of Fitted Charge is smaller than "
            << fChargeErrLimit << " in Pixel  " << fPixId << endl;
      return kFALSE;
    }
      
   if (GetMeanCharge() < fChargeRelErrLimit*GetMeanChargeErr()) 
    {
      *fLog << warn << "WARNING: Fitted Charge is smaller than "
            << fChargeRelErrLimit << "* its error in Pixel  " << fPixId << endl;
      return kFALSE;
    }
      
  if (GetSigmaCharge() < GetPedRms())
    {
	*fLog << warn << "WARNING: Sigma of Fitted Charge smaller than Pedestal RMS in Pixel  " << fPixId << endl;
	return kFALSE;
    }
 
  SetChargeValid();
  return kTRUE;
}

//
// The check returns kTRUE if:
//
// The mean arrival time is at least 1.0 slices from the used edge slices 
//
Bool_t MCalibrationChargePix::CheckTimeValidity()
{

    const Float_t loweredge     = IsHiGainSaturation() ? fTimeFirstHiGain   : fTimeFirstLoGain;
    const Float_t upperedge     = IsHiGainSaturation() ? fTimeLastHiGain    : fTimeLastLoGain;

    if ( fAbsTimeMean < loweredge+1.)
    {
	*fLog << warn << "WARNING: Mean ArrivalTime in first extraction bin of the Pixel " << fPixId << endl;
	SetMeanTimeInFirstBin();
	return kFALSE;
    }
    
    if ( fAbsTimeMean > upperedge-1.)
    {
	*fLog << warn << "WARNING: Mean ArrivalTime in last extraction bin of the Pixel " << fPixId << endl;
	SetMeanTimeInLastBin();
	return kFALSE;
    }

  return kTRUE;
}

void MCalibrationChargePix::CalcLoGainPed()
{

    Float_t pedRmsSquare                    =       fPedRms * fPedRms;
    Float_t pedRmsSquareErrSquare           =       fPedErr * fPedErr * pedRmsSquare; // fPedRmsErr = fPedErr/2.
    Float_t pedRmsSquareErr;
    
    //
    // 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          =    fElectronicPedRms    * fElectronicPedRms;
    const Float_t elecRmsSquareErrSquare = 4.*fElectronicPedRmsErr * fElectronicPedRmsErr * elecRmsSquare;
    
    Float_t nsbSquare             =  pedRmsSquare          - elecRmsSquare;
    Float_t nsbSquareRelErrSquare = (pedRmsSquareErrSquare + elecRmsSquareErrSquare)
                                	/ (nsbSquare * nsbSquare) ;
    
    if (nsbSquare < 0.)
        nsbSquare = 0.;
    
    //
    // Now, we divide the NSB by the conversion factor and 
    // add it quadratically to the electronic noise
    //
    const Float_t conversionSquare             =    fConversionHiLo    * fConversionHiLo;
    const Float_t conversionSquareRelErrSquare = 4.*fConversionHiLoErr * fConversionHiLoErr / conversionSquare;
    
    const Float_t convertedNsbSquare          =  nsbSquare             / conversionSquare;
    const Float_t convertedNsbSquareErrSquare = (nsbSquareRelErrSquare + conversionSquareRelErrSquare)
                                        	* convertedNsbSquare * convertedNsbSquare;
    
    pedRmsSquare     = convertedNsbSquare           + elecRmsSquare;
    pedRmsSquareErr  = TMath::Sqrt( convertedNsbSquareErrSquare  + elecRmsSquareErrSquare );
    
    //
    // Correct also for the conversion to Hi-Gain:
    //
    fLoGainPedRms    = TMath::Sqrt(pedRmsSquare)               * fConversionHiLo;
    fLoGainPedRmsErr = 0.5 * pedRmsSquareErr /  fLoGainPedRms  * fConversionHiLo;
}

Bool_t MCalibrationChargePix::CalcReducedSigma()
{

  const Float_t sigmaSquare               =    GetSigmaCharge()   * GetSigmaCharge();
  const Float_t sigmaSquareErrSquare      = 4.*GetSigmaChargeErr()* GetSigmaChargeErr() * sigmaSquare;

  Float_t pedRmsSquare;         
  Float_t pedRmsSquareErrSquare;
  
  if (IsHiGainSaturation())
  {
      CalcLoGainPed();

      pedRmsSquare             =      fLoGainPedRms    * fLoGainPedRms;                                          
      pedRmsSquareErrSquare    =  4.* fLoGainPedRmsErr * fLoGainPedRmsErr * pedRmsSquare;
    }   /* if (HiGainSaturation) */
  else
  {

      pedRmsSquare           =       fPedRms * fPedRms;                                          
      pedRmsSquareErrSquare  =       fPedErr * fPedErr * pedRmsSquare; // fPedRmsErr = fPedErr/2.
  }
  //
  // Calculate the reduced sigmas
  //
  const Float_t rsigmachargesquare = sigmaSquare - pedRmsSquare;

  if (rsigmachargesquare <= 0.)
    {
      *fLog << warn 
            << "WARNING: Cannot calculate the reduced sigma: smaller than 0 in pixel " 
            << fPixId << endl;
      return kFALSE;
    }
  
  fRSigmaCharge    = TMath::Sqrt(rsigmachargesquare);
  fRSigmaChargeErr = TMath::Sqrt(sigmaSquareErrSquare + pedRmsSquareErrSquare) / 2. / fRSigmaCharge;

  return kTRUE;
}

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

  if (!IsChargeValid())
      return kFALSE;

  if (fRSigmaCharge < 0.)
      return kFALSE;

  //
  // Square all variables in order to avoid applications of square root
  //
  // First the relative error squares
  //
  const Float_t chargeSquare              =     GetMeanCharge()    * GetMeanCharge();
  const Float_t chargeSquareRelErrSquare  = 4.* GetMeanChargeErr() * GetMeanChargeErr() / chargeSquare;

  const Float_t chargeRelErrSquare        =     GetMeanChargeErr() * GetMeanChargeErr() 
                                             / (GetMeanCharge()    * GetMeanCharge());
  
  const Float_t ffactorsquare             =    gkFFactor    * gkFFactor;
  const Float_t ffactorsquareRelErrSquare = 4.*gkFFactorErr * gkFFactorErr / ffactorsquare;

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

  const Float_t avQEFFactor               = TMath::Sqrt( ( 1. - fAverageQE ) / fAverageQE  );
  const Float_t avQEFFactorErr            = 1./ ( 2. * avQEFFactor ) * fAverageQEErr 
                                            / ( fAverageQE * fAverageQE );

  const Float_t avQEFFactorRelErrSquare  = avQEFFactorErr  * avQEFFactorErr 
                                           / ( avQEFFactor * avQEFFactor) ;

  const Float_t rsigmaSquare              =     fRSigmaCharge    * fRSigmaCharge;
  const Float_t rsigmaSquareRelErrSquare  = 4.* fRSigmaChargeErr * fRSigmaChargeErr / rsigmaSquare;


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

  if (fPheFFactorMethod < 0.)
      return kFALSE;

  //
  // Calculate the Error of Nphe
  //
  const Float_t pheFFactorRelErrSquare =  ffactorsquareRelErrSquare
                                         + chargeSquareRelErrSquare
                                         + rsigmaSquareRelErrSquare;
  fPheFFactorMethodErr                 =  TMath::Sqrt(pheFFactorRelErrSquare) * fPheFFactorMethod;

  //
  // Calculate the conversion factor between PHOTONS and FADC counts
  //
  // Nphot = Nphe / avQE
  // conv  = Nphot / FADC counts
  //
  fMeanConversionFFactorMethod         =  fPheFFactorMethod / fAverageQE / GetMeanCharge();


  //
  // Calculate the error of the mean conversion factor between PHOTONS and FADC counts
  //
  const Float_t convRelErrSquare     =  ( pheFFactorRelErrSquare + chargeRelErrSquare + avQERelErrSquare);

  if (convRelErrSquare < 0.)
      return kFALSE;


  const Float_t convRelErr           =  TMath::Sqrt(convRelErrSquare);
  fConversionFFactorMethodErr        = convRelErr * fMeanConversionFFactorMethod;
  
  if (convRelErr < fConvFFactorRelErrLimit)
    SetFFactorMethodValid();

  //
  // Calculate the Total F-Factor of the camera (in photons)
  //
  fTotalFFactorFFactorMethod =  (fRSigmaCharge/GetMeanCharge())*TMath::Sqrt(fPheFFactorMethod);
  fTotalFFactorFFactorMethod *=  avQEFFactor;
      
  //
  // Calculate the error of the Total F-Factor of the camera ( in photons )
  //
  const Float_t rsigmaChargeRelErrSquare =    fRSigmaChargeErr * fRSigmaChargeErr 
	                                   / (fRSigmaCharge    * fRSigmaCharge) ;
      
  fTotalFFactorErrFFactorMethod = TMath::Sqrt(  rsigmaChargeRelErrSquare
						+ chargeRelErrSquare
						+ pheFFactorRelErrSquare
						+ avQEFFactorRelErrSquare  );
  fTotalFFactorErrFFactorMethod *= fTotalFFactorFFactorMethod;

  //
  // Calculate the sigma of the conversion from FADC counts to photons
  //
  fSigmaConversionFFactorMethod = GetTotalFFactorFFactorMethod()*TMath::Sqrt(fMeanConversionFFactorMethod);

  return kTRUE;
}

void MCalibrationChargePix::ApplyLoGainConversion()
{
  
  const Float_t chargeRelErrSquare     =    fLoGainMeanChargeErr * fLoGainMeanChargeErr
                                         /( fLoGainMeanCharge    * fLoGainMeanCharge  );
  const Float_t sigmaRelErrSquare      =    fLoGainSigmaChargeErr * fLoGainSigmaChargeErr
                                         /( fLoGainSigmaCharge    * fLoGainSigmaCharge );
  const Float_t conversionRelErrSquare =    fConversionHiLoErr * fConversionHiLoErr 
                                         /( fConversionHiLo    * fConversionHiLo  );
  
  fLoGainMeanCharge         *= fConversionHiLo;
  fLoGainMeanChargeErr       = TMath::Sqrt(chargeRelErrSquare + conversionRelErrSquare) * fLoGainMeanCharge;
  
  fLoGainSigmaCharge    *= fConversionHiLo;
  fLoGainSigmaChargeErr =  TMath::Sqrt(sigmaRelErrSquare + conversionRelErrSquare) * fLoGainSigmaCharge;
  
  fElectronicPedRms     = gkElectronicPedRms    * TMath::Sqrt(fNumLoGainSamples);
  fElectronicPedRmsErr  = gkElectronicPedRmsErr * TMath::Sqrt(fNumLoGainSamples);

}

