Index: trunk/MagicSoft/Mars/mcalib/MCalibrationChargePix.cc
===================================================================
--- trunk/MagicSoft/Mars/mcalib/MCalibrationChargePix.cc	(revision 3648)
+++ trunk/MagicSoft/Mars/mcalib/MCalibrationChargePix.cc	(revision 3650)
@@ -22,10 +22,9 @@
 \* ======================================================================== */
 /////////////////////////////////////////////////////////////////////////////
-//                                                                         //
-// MCalibrationChargePix                                                   //
-//                                                                         //
-// Storage container to hold informations about the calibration values     //
-// values of one Pixel (PMT).                                              //
-//                                                                         //
+//                                                                         
+// MCalibrationChargePix                                                   
+//                                                                         
+// Storage container of the calibrated Charge of one pixel.
+//                                                                         
 // The following values are initialized to meaningful values:
 //
@@ -35,54 +34,38 @@
 //   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.
-// 
+//   Mean F-Factor (gkFFactor)     = 1.15
+//   Error F-Factor (gkFFactorErr) = 0.02
+//
+// The following variables are calculated inside this class:
+// -  fLoGainPedRmsSquare and LoGainPedRmsSquareVar (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
+// - fSigmaConversionFFactorMethodm, 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 
 // 
+// See also: MCalibrationChargeCam, MCalibrationChargeCalc,
+//           MHCalibrationChargeCam, MHCalibrationChargePix
+//
 /////////////////////////////////////////////////////////////////////////////
 #include "MCalibrationChargePix.h"
@@ -109,4 +92,13 @@
 // Default Constructor: 
 //
+// Sets:
+// - fCalibFlags to 0
+// - fConversionHiLo to fgConversionHiLo
+// - fConversionHiLoErr to fgConversionHiLoErr
+// - PheFFactorMethodLimit to fgPheFFactorMethodLimit
+// 
+// Calls:
+// - Clear()
+//
 MCalibrationChargePix::MCalibrationChargePix(const char *name, const char *title)
     : fCalibFlags(0)
@@ -131,5 +123,10 @@
 // ------------------------------------------------------------------------
 //
-// Invalidate values
+// Sets:
+// - all flags to 0
+// - all variables to -1.
+//
+// Calls: 
+// - MCalibrationPix::Clear()
 //
 void MCalibrationChargePix::Clear(Option_t *o)
@@ -141,6 +138,6 @@
   SetCombinedMethodValid    ( kFALSE );
 
-  fRSigma                           =  -1.;
-  fRSigmaVar                        =  -1.;
+  fRSigmaSquare                     =  -1.;
+  fRSigmaSquareVar                  =  -1.;
   
   fPed                              =  -1.;
@@ -148,6 +145,6 @@
   fPedVar                           =  -1.;
 
-  fLoGainPedRms                     =  -1.;
-  fLoGainPedRmsVar                  =  -1.;
+  fLoGainPedRmsSquare               =  -1.;
+  fLoGainPedRmsSquareVar            =  -1.;
 
   fAbsTimeMean                      =  -1.;
@@ -183,19 +180,29 @@
 // --------------------------------------------------------------------------
 //
-// Set the pedestals from outside
-//
-void MCalibrationChargePix::SetPedestal(const Float_t ped, const Float_t pedrms, const Float_t pederr)
-{
-
-  fPed       = ped;    
-  fPedRms    = pedrms;
-  fPedVar    = pederr*pederr;
-}
-  
-
-
-// --------------------------------------------------------------------------
-//
-// Set the conversion factors from outside (only for MC)
+// Set the 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 the 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 the conversion factors F-Factor Method from outside (only for MC)
 //
 void MCalibrationChargePix::SetConversionFFactorMethod(Float_t c, Float_t err, Float_t sig)
@@ -208,28 +215,5 @@
 // --------------------------------------------------------------------------
 //
-// Set the conversion factors 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 the conversion factors 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 the conversion factors from outside (only for MC)
+// Set the conversion factors PIN Diode Method from outside (only for MC)
 //
 void MCalibrationChargePix::SetConversionPINDiodeMethod(Float_t c, Float_t err, Float_t sig)
@@ -242,5 +226,5 @@
 // --------------------------------------------------------------------------
 //
-// Set the Excluded Bit from outside 
+// Set the BlindPixelMethod Validity Bit from outside 
 //
 void MCalibrationChargePix::SetBlindPixelMethodValid(const Bool_t b )
@@ -251,5 +235,14 @@
 // --------------------------------------------------------------------------
 //
-// Set the Excluded Bit from outside 
+// Set the CombinedMethod Validity Bit from outside 
+//
+void MCalibrationChargePix::SetCombinedMethodValid(const Bool_t b )  
+{ 
+  b ?  SETBIT(fCalibFlags, kCombinedMethodValid) : CLRBIT(fCalibFlags, kCombinedMethodValid); 
+}
+
+// --------------------------------------------------------------------------
+//
+// Set the FFactorMethod Validity Bit from outside 
 //
 void MCalibrationChargePix::SetFFactorMethodValid(const Bool_t b )
@@ -260,5 +253,5 @@
 // --------------------------------------------------------------------------
 //
-// Set the Excluded Bit from outside 
+// Set the PINDiodeMethod Validity Bit from outside 
 //
 void MCalibrationChargePix::SetPINDiodeMethodValid(const Bool_t b )  
@@ -269,28 +262,70 @@
 // --------------------------------------------------------------------------
 //
-// Set the Excluded Bit from outside 
-//
-void MCalibrationChargePix::SetCombinedMethodValid(const Bool_t b )  
-{ 
-  b ?  SETBIT(fCalibFlags, kCombinedMethodValid) : CLRBIT(fCalibFlags, kCombinedMethodValid); 
-}
-
-
+// Set the pedestals from outside
+//
+void MCalibrationChargePix::SetPedestal(const Float_t ped, const Float_t pedrms, const Float_t pederr)
+{
+
+  fPed       = ped;    
+  fPedRms    = pedrms;
+  fPedVar    = pederr*pederr;
+}
+  
+
+// --------------------------------------------------------------------------
+//
+// 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
 {
-  return IsHiGainSaturation() ? fLoGainPedRms : fPedRms;
-}
-
+
+  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 kLoGainPedRms, else fPedRms
+//
 Float_t  MCalibrationChargePix::GetPedRmsErr()  const
 {
-  return IsHiGainSaturation() ? TMath::Sqrt(fLoGainPedRmsVar) : TMath::Sqrt(fPedVar)/2.;
-}
-
-Float_t  MCalibrationChargePix::GetPedErr()  const
-{
-  return  TMath::Sqrt(fPedVar);
-}
-
-
+  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: 
+// Returns fLoGainMean multiplied with fConversionHiLo
+//
 Float_t MCalibrationChargePix::GetLoGainMean()  const 
 {
@@ -298,4 +333,10 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// Get the Error of the Low Gain Mean: 
+// Returns the quadratic sum of the relative low Gain Mean error and the 
+// the relative conversion High-to-Low error, mulitplied with GetLoGainMean()
+//
 Float_t MCalibrationChargePix::GetLoGainMeanErr()  const
 {
@@ -304,10 +345,12 @@
                                  /( fLoGainMean * fLoGainMean );
 
-  const Float_t conversionRelVar =  fConversionHiLoVar
-                                 /( fConversionHiLo   * fConversionHiLo   );
-
-  return TMath::Sqrt(chargeRelVar+conversionRelVar) * GetLoGainMean();
-}
-
+  return TMath::Sqrt(chargeRelVar+GetConversionHiLoRelVar()) * GetLoGainMean();
+}
+
+// --------------------------------------------------------------------------
+//
+// Get the Low Gain Sigma: 
+// Returns fLoGainSigma multiplied with fConversionHiLo
+//
 Float_t MCalibrationChargePix::GetLoGainSigma()  const 
 {
@@ -315,4 +358,10 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// Get the Error of the Low Gain Sigma: 
+// Returns the quadratic sum of the relative low Gain Sigma error and the 
+// the relative conversion High-to-Low error, mulitplied with GetLoGainSigma()
+//
 Float_t MCalibrationChargePix::GetLoGainSigmaErr()  const
 {
@@ -321,37 +370,75 @@
                                 /( fLoGainSigma * fLoGainSigma );
 
-  const Float_t conversionRelVar =  fConversionHiLoVar
-                                 /( fConversionHiLo   * fConversionHiLo   );
-
-  return TMath::Sqrt(sigmaRelVar+conversionRelVar) * GetLoGainSigma();
-}
-
+  return TMath::Sqrt(sigmaRelVar+GetConversionHiLoRelVar()) * GetLoGainSigma();
+}
+
+// --------------------------------------------------------------------------
+//
+// 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
 {
-  return IsHiGainSaturation() ? fRSigma*fConversionHiLo : fRSigma ;
+  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:
+//   sigma variance = 0.25 * fRSigmaSquareVar / fRSigmaSquare
+// - Test bit kHiGainSaturation: 
+//   If yes, returns the quadratic sum of the relative reduced Sigma error and the 
+//    the relative conversion High-to-Low error, mulitplied with GetRSigma()
+//   Else returns the square root of rel. (0.25*fRSigmaSquareVar / fRSigmaSquare)
+//
 Float_t MCalibrationChargePix::GetRSigmaErr()  const
 {
- if (IsHiGainSaturation())
-    {
-      const Float_t rsigmaRelVar     =  fRSigmaVar
-                                    /( fRSigma * fRSigma );
-      const Float_t conversionRelVar =  fConversionHiLoVar
-                                     /( fConversionHiLo   * fConversionHiLo   );
-      return TMath::Sqrt(rsigmaRelVar+conversionRelVar) * GetRSigma();
-    }
+
+  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(fRSigmaVar);
-
-}
-
+   return TMath::Sqrt(rsigmaVar);
+
+}
+
+// --------------------------------------------------------------------------
+//
+// Get the conversion Error Hi-Gain to Low-Gain:
+// - If fConversionHiLoVar is smaller than 0 (i.e. has not yet been set), return -1.
+// - Else returns the square root of fConversionHiLoVar
+//
 Float_t MCalibrationChargePix::GetConversionHiLoErr()  const
 {
   if (fConversionHiLoVar < 0.)
     return -1.;
+
   return TMath::Sqrt(fConversionHiLoVar);
 }
 
+// --------------------------------------------------------------------------
+//
+// 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
 {
@@ -361,4 +448,10 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// 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
 {
@@ -368,4 +461,10 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// 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
 {
@@ -375,4 +474,10 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// 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
 {
@@ -382,4 +487,10 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// 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
 {
@@ -389,4 +500,10 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// 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
 {
@@ -396,4 +513,10 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// 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
 {
@@ -403,4 +526,10 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// 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
 {
@@ -410,4 +539,10 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// 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
 {
@@ -417,4 +552,8 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// Test bit kBlindPixelMethodValid
+//
 Bool_t MCalibrationChargePix::IsBlindPixelMethodValid() const 
 { 
@@ -422,4 +561,8 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// Test bit kFFactorMethodValid
+//
 Bool_t MCalibrationChargePix::IsFFactorMethodValid()   const
 { 
@@ -427,4 +570,8 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// Test bit kPINDiodeMethodValid
+//
 Bool_t MCalibrationChargePix::IsPINDiodeMethodValid()  const 
 { 
@@ -432,4 +579,8 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// Test bit kCombinedMethodValid
+//
 Bool_t MCalibrationChargePix::IsCombinedMethodValid()  const
 { 
@@ -437,34 +588,46 @@
 }
 
-//
+// ----------------------------------------------------------------------------
 // 
+// - 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:
+// - 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 fPedRmsVar
+//
+// - Calculate the reduced sigma with the formula:
+//   fRSigmaSquare = Sigma*Sigma - pedRMS*pedRMS
+// 
+// - Calculate the variance of the reduced sigma with the formula:
+//   fRSigmaSquareVar = 4.* (sigmaVar*Sigma*Sigma + pedRmsVar*pedRMS*pedRMS)
 //
 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 sigmaSquare     =      sigma     * sigma;
-  const Float_t sigmaSquareVar  = 4.*  sigmavar  * sigmaSquare;
-
-  Float_t pedRmsSquare ;         
-  Float_t pedRmsSquareVar;
-  
-  if (IsHiGainSaturation())
-    {
-      pedRmsSquare     =      fLoGainPedRms    * fLoGainPedRms;
-      pedRmsSquareVar  =  4.* fLoGainPedRmsVar * pedRmsSquare;
-    }
-  else
-  {
-      pedRmsSquare       =  fPedRms * fPedRms;                                          
-      pedRmsSquareVar    =  fPedVar * pedRmsSquare; // fPedRmsErr = fPedErr/2.
-  }
+  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
   //
-  const Float_t rsigmasquare = sigmaSquare - pedRmsSquare;
-  if (rsigmasquare <= 0.)
+  fRSigmaSquare = sigmaSquare - pedRmsSquare;
+  if (fRSigmaSquare <= 0.)
     {
       *fLog << warn 
@@ -473,20 +636,40 @@
       return kFALSE;
     }
-  
-
-  fRSigma    = TMath::Sqrt(rsigmasquare);
-  fRSigmaVar = 0.25 * (sigmaSquareVar + pedRmsSquareVar) / rsigmasquare;
+
+  fRSigmaSquareVar = 4. * (sigmaSquareVar + pedRmsSquareVar);
 
   return kTRUE;
 }
 
-//
-// Calculate the number of photo-electrons after the F-Factor method
-// Calculate the errors of the F-Factor method
+// ------------------------------------------------------------------
+//
+// 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 fPedRmsVar
+//
+// - 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. * fMeanVar                    / ( mean      * mean      )
+//                            + 4. * fRSigmaVar                  / ( fRSigma   * fRSigma   )
+//                            ) * fPheFFactor * fPheFFactor
+// - 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 (fRSigma < 0.)
+  if (fRSigmaSquare < 0.)
     {
       SetFFactorMethodValid(kFALSE);
@@ -500,14 +683,11 @@
   // Square all variables in order to avoid applications of square root
   //
-  // First the relative error squares
-  //
-  const Float_t meanSquare        =     mean * mean;
-  const Float_t meanSquareRelVar  = 4.* var/ meanSquare;
-
-  const Float_t ffactorsquare       =    gkFFactor    * gkFFactor;
-  const Float_t ffactorsquareRelVar = 4.*gkFFactorErr * gkFFactorErr / ffactorsquare;
-
-  const Float_t rsigmaSquare        =     fRSigma    * fRSigma;
-  const Float_t rsigmaSquareRelVar  = 4.* fRSigmaVar / rsigmaSquare;
+  const Float_t meanSquare          =     mean * mean;
+  const Float_t meanSquareRelVar    = 4.* var / meanSquare;
+
+  const Float_t ffactorsquare       =     gkFFactor    * gkFFactor;
+  const Float_t ffactorsquareRelVar = 4.* gkFFactorErr * gkFFactorErr / ffactorsquare;
+
+  const Float_t rsigmaSquareRelVar  =     fRSigmaSquareVar / fRSigmaSquare;
 
   //
@@ -515,5 +695,5 @@
   // (independent on Hi Gain or Lo Gain)
   //
-  fPheFFactorMethod = ffactorsquare * meanSquare / rsigmaSquare;
+  fPheFFactorMethod = ffactorsquare * meanSquare / fRSigmaSquare;
 
   if (fPheFFactorMethod < fPheFFactorMethodLimit)
@@ -534,6 +714,36 @@
 
 
+// ----------------------------------------------------------------------------
+// 
+// - 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 LONS ped. RMS contribution 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 PedRMS (LONS,lowgain) square is smaller than 0., set it to zero. (but not the error!)
+//
+// Convert the LONS 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);
@@ -552,10 +762,10 @@
   const Float_t elecRmsSquareVar = 4.*elecPedRmsVar * elecRmsSquare;
   
-  Float_t nsbSquare             =  pedRmsSquare    - elecRmsSquare;
-  Float_t nsbSquareRelVar       = (pedRmsSquareVar + elecRmsSquareVar)
-                                 / (nsbSquare * nsbSquare) ;
-  
-  if (nsbSquare < 0.)
-    nsbSquare = 0.;
+  Float_t higainNsbSquare        =  pedRmsSquare    - elecRmsSquare;
+  Float_t higainNsbSquareRelVar  = (pedRmsSquareVar + elecRmsSquareVar)
+                                 / (higainNsbSquare * higainNsbSquare) ;
+  
+  if (higainNsbSquare < 0.)
+    higainNsbSquare = 0.;
   
   //
@@ -563,17 +773,21 @@
   // add it quadratically to the electronic noise
   //
-  const Float_t conversionSquare        =    fConversionHiLo    * fConversionHiLo;
-  const Float_t convertedNsbSquare      =    nsbSquare       / conversionSquare;
-  const Float_t convertedNsbSquareVar   =    nsbSquareRelVar
-                                	    * convertedNsbSquare * convertedNsbSquare;
+  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;
     
-  pedRmsSquare     = convertedNsbSquare    + elecRmsSquare;
-  pedRmsSquareVar  = convertedNsbSquareVar + elecRmsSquareVar;
-  
-  fLoGainPedRms    = TMath::Sqrt(pedRmsSquare);
-  fLoGainPedRmsVar = 0.25 * pedRmsSquareVar /  pedRmsSquare;
-  
+  fLoGainPedRmsSquare    = logainNsbSquare    + elecRmsSquare;
+  fLoGainPedRmsSquareVar = logainNsbSquareVar + elecRmsSquareVar;
 }
  
+const Float_t MCalibrationChargePix::GetConversionHiLoRelVar() const 
+{
+  if (fConversionHiLo == 0.)
+    return 0.;
+
+  return fConversionHiLoVar / (fConversionHiLo * fConversionHiLo);
+}
  
- 
Index: trunk/MagicSoft/Mars/mcalib/MCalibrationPedCam.cc
===================================================================
--- trunk/MagicSoft/Mars/mcalib/MCalibrationPedCam.cc	(revision 3650)
+++ trunk/MagicSoft/Mars/mcalib/MCalibrationPedCam.cc	(revision 3650)
@@ -0,0 +1,139 @@
+/* ======================================================================== *\
+!
+! *
+! * This file is part of MARS, the MAGIC Analysis and Reconstruction
+! * Software. It is distributed to you in the hope that it can be a useful
+! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
+! * It is distributed WITHOUT ANY WARRANTY.
+! *
+! * Permission to use, copy, modify and distribute this software and its
+! * documentation for any purpose is hereby granted without fee,
+! * provided that the above copyright notice appear in all copies and
+! * that both that copyright notice and this permission notice appear
+! * in supporting documentation. It is provided "as is" without express
+! * or implied warranty.
+! *
+!
+!
+!   Author(s): Markus Gaug   11/2003 <mailto:markus@ifae.es>
+!
+!   Copyright: MAGIC Software Development, 2000-2004
+!
+!
+\* ======================================================================== */
+
+/////////////////////////////////////////////////////////////////////////////
+//                                                               
+// MCalibrationPedCam                                               
+//                                                               
+// Hold the pedestal Calibration results obtained from MHPedestalCam of the camera:
+//                                                               
+// 1) MCalibrationPedCam initializes a TClonesArray whose elements are 
+//    pointers to MCalibrationPix Containers
+//
+// The calculated values (types of GetPixelContent) are:
+// 
+// Fitted values:
+// ============== 
+//
+// 0: Fitted Pedestal
+// 1: Error of fitted Pedestal
+// 2: Sigma of fitted Pedestal
+// 3: Error of Sigma of fitted Pedestal
+//
+// Useful variables derived from the fit results:
+// =============================================
+//
+// 4: Returned probability of Gauss fit to Pedestal distribution
+//
+/////////////////////////////////////////////////////////////////////////////
+#include "MCalibrationPedCam.h"
+#include "MCalibrationCam.h"
+
+#include <TClonesArray.h>
+
+#include "MLog.h"
+#include "MLogManip.h"
+
+#include "MGeomCam.h"
+#include "MGeomPix.h"
+
+#include "MCalibrationPix.h"
+
+ClassImp(MCalibrationPedCam);
+
+using namespace std;
+
+// --------------------------------------------------------------------------
+//
+// Default constructor. 
+//
+// Creates a TClonesArray of MCalibrationPix containers, initialized to 1 entry
+// Later, a call to MCalibrationPedCam::InitSize(Int_t size) has to be performed
+//
+MCalibrationPedCam::MCalibrationPedCam(const char *name, const char *title)
+{
+    fName  = name  ? name  : "MCalibrationPedCam";
+    fTitle = title ? title : "Storage container for the Pedestal Calibration Results in the camera";
+
+    fPixels           = new TClonesArray("MCalibrationPix",1);
+    fAverageAreas     = new TClonesArray("MCalibrationPix",1);
+    fAverageSectors   = new TClonesArray("MCalibrationPix",1);
+
+}
+
+// --------------------------------------------------------------------------
+//
+// The calculated values (types of GetPixelContent) are:
+// 
+// Fitted values:
+// ============== 
+//
+// 0: Fitted Pedestal
+// 1: Error of fitted Pedestal
+// 2: Sigma of fitted Pedestal
+// 3: Error of Sigma of fitted Pedestal
+//
+// Useful variables derived from the fit results:
+// =============================================
+//
+// 4: Returned probability of Gauss fit to Pedestal distribution
+//
+Bool_t MCalibrationPedCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
+{
+
+  if (idx > GetSize())
+    return kFALSE;
+
+  Float_t area = cam[idx].GetA();
+
+ if (area == 0)
+    return kFALSE;
+
+ if ((*this)[idx].IsExcluded())
+   return kFALSE;
+
+  switch (type)
+    {
+    case 0:
+      val = (*this)[idx].GetMean();
+      break;
+    case 1:
+      val = (*this)[idx].GetMeanErr();
+      break;
+    case 2:
+      val = (*this)[idx].GetSigma();
+      break;
+    case 3:
+      val = (*this)[idx].GetSigmaErr();
+      break;
+    case 4:
+      val = (*this)[idx].GetProb();
+      break;
+    default:
+      return kFALSE;
+    }
+
+  return val!=-1.;
+  
+}
Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationCam.cc
===================================================================
--- trunk/MagicSoft/Mars/mcalib/MHCalibrationCam.cc	(revision 3648)
+++ trunk/MagicSoft/Mars/mcalib/MHCalibrationCam.cc	(revision 3650)
@@ -703,4 +703,5 @@
   
   hist.Renorm();
+
   //
   // 4) Check for oscillations
@@ -708,4 +709,5 @@
   hist.CreateFourierSpectrum();
   
+
   if (!hist.IsFourierSpectrumOK())
     bad.SetUncalibrated( osctyp );
Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeCam.cc
===================================================================
--- trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeCam.cc	(revision 3648)
+++ trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeCam.cc	(revision 3650)
@@ -216,7 +216,24 @@
 {
 
-  fCam = (MCalibrationCam*)pList->FindCreateObj("MCalibrationChargeCam");
+  const Int_t npixels  = fGeom->GetNumPixels();
+  const Int_t nsectors = fGeom->GetNumSectors();
+  const Int_t nareas   = fGeom->GetNumAreas();
+
+  fCam = (MCalibrationCam*)pList->FindObject("MCalibrationChargeCam");
   if (!fCam)
-      return kFALSE;
+    {
+      fCam = (MCalibrationCam*)pList->FindCreateObj(AddSerialNumber("MCalibrationChargeCam"));
+      if (!fCam)
+        {
+          gLog << err << "Cannot find nor create MCalibrationChargeCam ... abort." << endl;
+          return kFALSE;
+        }
+      else 
+        {
+          fCam->InitSize(npixels);
+          fCam->InitAverageAreas(nareas);
+          fCam->InitAverageSectors(nsectors);
+        }
+    }
 
   MExtractedSignalCam *signal = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
@@ -226,8 +243,4 @@
       return kFALSE;
   }
-
-  const Int_t npixels  = fGeom->GetNumPixels();
-  const Int_t nsectors = fGeom->GetNumSectors();
-  const Int_t nareas   = fGeom->GetNumAreas();
 
   if (fHiGainArray->GetEntries()==0)
Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimeCam.cc
===================================================================
--- trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimeCam.cc	(revision 3648)
+++ trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimeCam.cc	(revision 3650)
@@ -153,8 +153,21 @@
 {
 
-  fCam = (MCalibrationCam*)pList->FindCreateObj("MCalibrationRelTimeCam");
+  fCam = (MCalibrationCam*)pList->FindObject("MCalibrationRelTimeCam");
   if (!fCam)
-      return kFALSE;
-
+    {
+      fCam = (MCalibrationCam*)pList->FindCreateObj(AddSerialNumber("MCalibrationRelTimeCam"));
+      if (!fCam)
+        {
+          gLog << err << "Cannot find nor create MCalibrationRelTimeCam ... abort." << endl;
+          return kFALSE;
+        }
+      else 
+        {
+          fCam->InitSize(npixels);
+          fCam->InitAverageAreas(nareas);
+          fCam->InitAverageSectors(nsectors);
+        }
+    }
+  
   MArrivalTimeCam *signal = (MArrivalTimeCam*)pList->FindObject("MArrivalTimeCam");
   if (!signal)
@@ -376,9 +389,9 @@
 {
 
-  FitHiGainArrays((MCalibrationCam&)(*fCam),*fBadPixels,
+  FitHiGainArrays((*fCam),*fBadPixels,
                   MBadPixelsPix::kRelTimeNotFitted,
                   MBadPixelsPix::kRelTimeOscillating);
 
-  FitLoGainArrays((MCalibrationCam&)(*fCam),*fBadPixels,
+  FitLoGainArrays((*fCam),*fBadPixels,
                   MBadPixelsPix::kRelTimeNotFitted,
                   MBadPixelsPix::kRelTimeOscillating);
