Index: trunk/MagicSoft/Mars/mpedestal/MPedCalcFromLoGain.cc
===================================================================
--- trunk/MagicSoft/Mars/mpedestal/MPedCalcFromLoGain.cc	(revision 4404)
+++ trunk/MagicSoft/Mars/mpedestal/MPedCalcFromLoGain.cc	(revision 4404)
@@ -0,0 +1,633 @@
+/* ======================================================================== *\
+!
+! *
+! * 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): Josep Flix 04/2001 <mailto:jflix@ifae.es>
+!   Author(s): Thomas Bretz 05/2001 <mailto:tbretz@astro.uni-wuerzburg.de>
+!   Author(s): Sebastian Commichau 12/2003 
+!   Author(s): Javier Rico 01/2004 <mailto:jrico@ifae.es>
+!   Author(s): Markus Gaug 01/2004 <mailto:markus@ifae.es>
+!   Author(s): Florian Goebel 06/2004 <mailto:fgoebel@mppmu.mpg.de>
+!
+!   Copyright: MAGIC Software Development, 2000-2004
+!
+!
+\* ======================================================================== */
+/////////////////////////////////////////////////////////////////////////////
+//
+//   MPedCalcLoGain
+//
+//
+//  This task is devide form MPedCalcPedRun, described below. However, It 
+//  calculates the pedstals using the low gain slices, whenever the difference 
+//  between the highest and the lowest slice in the high gain
+//  slices is below a given threshold (SetMaxHiGainVar). In this case the receiver
+//  boards do not switch to lo gain and the so called lo gain slices are actually
+//  high gain slices. 
+//
+//  MPedCalcLoGain also fills the ABoffset in MPedestalPix which allows to correct 
+//  the 150 MHz clock noise.
+//
+//  This task takes a pedestal run file and fills MPedestalCam during
+//  the Process() with the pedestal and rms computed in an event basis.
+//  In the PostProcess() MPedestalCam is finally filled with the pedestal
+//  mean and rms computed in a run basis.
+//  More than one run (file) can be merged
+//
+//  MPedCalcPedRun applies the following formula (1):
+// 
+//  Pedestal per slice = sum(x_i) / n / slices
+//  PedRMS per slice   = Sqrt(  ( sum(x_i^2) - sum(x_i)^2/n ) / n-1 / slices )
+// 
+//  where x_i is the sum of "slices" FADC slices and sum means the sum over all
+//  events. "n" is the number of events, "slices" is the number of summed FADC samples.
+// 
+//  Note that the slice-to-slice fluctuations are not Gaussian, but Poissonian, thus 
+//  asymmetric and they are correlated.
+// 
+//  It is important to know that the Pedestal per slice and PedRMS per slice depend 
+//  on the number of used FADC slices, as seen in the following plots:
+//
+//Begin_Html
+/*
+<img src="images/PedestalStudyInner.gif">
+*/
+//End_Html
+//
+//Begin_Html
+/*
+<img src="images/PedestalStudyOuter.gif">
+*/
+//
+// The plots show the inner and outer pixels, respectivly and have the following meaning:
+// 
+// 1) The calculated mean pedestal per slice (from MPedCalcFromLoGain)
+// 2) The fitted mean pedestal per slice     (from MHPedestalCam)
+// 3) The calculated pedestal RMS per slice  (from MPedCalcFromLoGain)
+// 4) The fitted sigma of the pedestal distribution per slice
+//                                           (from MHPedestalCam)
+// 5) The relative difference between calculation and histogram fit
+//    for the mean
+// 6) The relative difference between calculation and histogram fit
+//    for the sigma or RMS, respectively.
+// 
+// The calculated means do not change significantly except for the case of 2 slices, 
+// however the RMS changes from 5.7 per slice in the case of 2 extracted slices 
+// to 8.3 per slice in the case of 26 extracted slices. This change is very significant.
+// 
+// The plots have been produced on run 20123. You can reproduce them using
+// the macro pedestalstudies.C
+// 
+//  Usage of this class: 
+//  ====================
+// 
+//  Call: SetRange(higainfirst, higainlast, logainfirst, logainlast) 
+//  to modify the ranges in which the window is allowed to move. 
+//  Defaults are: 
+// 
+//   fHiGainFirst =  fgHiGainFirst =  0 
+//   fHiGainLast  =  fgHiGainLast  =  29
+//   fLoGainFirst =  fgLoGainFirst =  0 
+//   fLoGainLast  =  fgLoGainLast  =  14
+//
+//  Call: SetWindowSize(windowhigain, windowlogain) 
+//  to modify the sliding window widths. Windows have to be an even number. 
+//  In case of odd numbers, the window will be modified.
+//
+//  Defaults are:
+//
+//   fHiGainWindowSize = fgHiGainWindowSize = 14
+//   fLoGainWindowSize = fgLoGainWindowSize = 0
+//
+//  Input Containers:
+//   MRawEvtData
+//   MRawRunHeader
+//   MGeomCam
+//
+//  Output Containers:
+//   MPedestalCam
+//
+//  See also: MPedestalCam, MPedestalPix, MHPedestalCam, MExtractor
+//
+/////////////////////////////////////////////////////////////////////////////
+#include "MPedCalcFromLoGain.h"
+#include "MExtractor.h"
+
+#include "MParList.h"
+
+#include "MLog.h"
+#include "MLogManip.h"
+
+#include "MRawRunHeader.h"  
+#include "MRawEvtPixelIter.h"
+#include "MRawEvtData.h"
+
+#include "MPedestalPix.h"
+#include "MPedestalCam.h"
+
+#include "MGeomPix.h"
+#include "MGeomCam.h"
+
+#include "MBadPixelsPix.h"
+#include "MBadPixelsCam.h"
+
+#include "MGeomCamMagic.h"
+
+ClassImp(MPedCalcFromLoGain);
+
+using namespace std;
+
+const Byte_t MPedCalcFromLoGain::fgHiGainFirst      = 0;
+const Byte_t MPedCalcFromLoGain::fgHiGainLast       = 11;
+const Byte_t MPedCalcFromLoGain::fgLoGainFirst      = 1;
+const Byte_t MPedCalcFromLoGain::fgLoGainLast       = 14;
+const Byte_t MPedCalcFromLoGain::fgHiGainWindowSize = 12;
+const Byte_t MPedCalcFromLoGain::fgLoGainWindowSize = 14;
+const Byte_t MPedCalcFromLoGain::fgMaxHiGainVar     = 20;
+
+// --------------------------------------------------------------------------
+//
+// Default constructor: 
+//
+// Sets:
+// - all pointers to NULL
+// - fWindowSizeHiGain to fgHiGainWindowSize
+// - fWindowSizeLoGain to fgLoGainWindowSize
+//
+// Calls: 
+// - AddToBranchList("fHiGainPixId");
+// - AddToBranchList("fHiGainFadcSamples");
+// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast)
+// - Clear()
+//
+MPedCalcFromLoGain::MPedCalcFromLoGain(const char *name, const char *title)
+    : fWindowSizeHiGain(fgHiGainWindowSize), 
+      fWindowSizeLoGain(fgLoGainWindowSize), 
+      fGeom(NULL),fBad(NULL)
+{
+  fName  = name  ? name  : "MPedCalcFromLoGain";
+  fTitle = title ? title : "Task to calculate pedestals from pedestal runs raw data";
+  
+  AddToBranchList("fHiGainPixId");
+  AddToBranchList("fHiGainFadcSamples");
+  
+  SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
+
+  SetMaxHiGainVar(fgMaxHiGainVar);
+  SetPedestalUpdate(kTRUE);
+  Clear();
+}
+
+// --------------------------------------------------------------------------
+//
+// Sets:
+// - fNumSamplesTot to 0
+// - fRawEvt to NULL
+// - fRunHeader to NULL
+// - fPedestals to NULL
+//
+void MPedCalcFromLoGain::Clear(const Option_t *o)
+{
+
+  fRawEvt    = NULL;
+  fRunHeader = NULL;
+  fPedestals = NULL;
+}
+
+// --------------------------------------------------------------------------
+//
+// SetRange: 
+//
+// Calls:
+// - MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
+// - SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain);
+//
+void MPedCalcFromLoGain::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast)
+{
+
+  MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
+
+  //
+  // Redo the checks if the window is still inside the ranges
+  //
+  SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain);
+  
+}
+
+
+// --------------------------------------------------------------------------
+//
+void MPedCalcFromLoGain::SetMaxHiGainVar(Byte_t maxvar) {
+  fMaxHiGainVar = maxvar;
+}
+
+
+// --------------------------------------------------------------------------
+//
+// Checks:
+// - if a window is odd, subtract one
+// - if a window is bigger than the one defined by the ranges, set it to the available range
+// - if a window is smaller than 2, set it to 2
+// 
+// Sets:
+// - fNumHiGainSamples to: (Float_t)fWindowSizeHiGain
+// - fNumLoGainSamples to: (Float_t)fWindowSizeLoGain
+// - fSqrtHiGainSamples to: TMath::Sqrt(fNumHiGainSamples)
+// - fSqrtLoGainSamples to: TMath::Sqrt(fNumLoGainSamples)  
+//  
+void MPedCalcFromLoGain::SetWindowSize(Byte_t windowh, Byte_t windowl)
+{
+  
+  fWindowSizeHiGain = windowh & ~1;
+  fWindowSizeLoGain = windowl & ~1;
+
+  if (fWindowSizeHiGain != windowh)
+    *fLog << warn << GetDescriptor() << ": Hi Gain window size has to be even, set to: " 
+          << int(fWindowSizeHiGain) << " samples " << endl;
+  
+  if (fWindowSizeLoGain != windowl)
+    *fLog << warn << GetDescriptor() << ": Lo Gain window size has to be even, set to: " 
+          << int(fWindowSizeLoGain) << " samples " << endl;
+    
+  const Byte_t availhirange = (fHiGainLast-fHiGainFirst+1) & ~1;
+  const Byte_t availlorange = (fLoGainLast-fLoGainFirst+1) & ~1;
+
+  if (fWindowSizeHiGain > availhirange)
+    {
+      *fLog << warn << GetDescriptor() 
+            << Form("%s%2i%s%2i%s%2i%s",": Hi Gain window size: ",(int)fWindowSizeHiGain,
+                    " is bigger than available range: [",(int)fHiGainFirst,",",(int)fHiGainLast,"]") << endl;
+      *fLog << warn << GetDescriptor() 
+            << ": Will set window size to: " << (int)availhirange << endl;
+      fWindowSizeHiGain = availhirange;
+    }
+  
+  if (fWindowSizeLoGain > availlorange)
+    {
+      *fLog << warn << GetDescriptor() 
+            << Form("%s%2i%s%2i%s%2i%s",": Lo Gain window size: ",(int)fWindowSizeLoGain,
+                    " is bigger than available range: [",(int)fLoGainFirst,",",(int)fLoGainLast,"]") << endl;
+      *fLog << warn << GetDescriptor() 
+            << ": Will set window size to: " << (int)availlorange << endl;
+      fWindowSizeLoGain = availlorange;
+    }
+  
+  
+  fNumHiGainSamples = (Float_t)fWindowSizeHiGain;
+  fNumLoGainSamples = (Float_t)fWindowSizeLoGain;
+  
+  fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
+  fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
+
+}
+
+  
+    
+// --------------------------------------------------------------------------
+//
+// Look for the following input containers:
+//
+//  - MRawEvtData
+//  - MRawRunHeader
+//  - MGeomCam
+//  - MBadPixelsCam
+// 
+// The following output containers are also searched and created if
+// they were not found:
+//
+//  - MPedestalCam
+//
+Int_t MPedCalcFromLoGain::PreProcess( MParList *pList )
+{
+
+  Clear();
+  
+  fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
+  if (!fRawEvt)
+    {
+      *fLog << err << "MRawEvtData not found... aborting." << endl;
+      return kFALSE;
+    }
+  
+  fRunHeader = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader"));
+  if (!fRunHeader)
+    {
+      *fLog << err << AddSerialNumber("MRawRunHeader") << " not found... aborting." << endl;
+      return kFALSE;
+    }
+
+  fGeom   =  (MGeomCam*)pList->FindObject("MGeomCam");
+  if (!fGeom)
+    {
+      *fLog << err << "MGeomCam not found... aborting." << endl;
+      return kFALSE;
+    }
+
+  fPedestals = (MPedestalCam*)pList->FindCreateObj("MPedestalCam");
+  if (!fPedestals)
+    return kFALSE;
+
+  fBad       =  (MBadPixelsCam*)pList->FindObject("MBadPixelsCam");
+
+  return kTRUE;
+}
+
+// --------------------------------------------------------------------------
+//
+// The ReInit searches for:
+// -  MRawRunHeader::GetNumSamplesHiGain()
+// -  MRawRunHeader::GetNumSamplesLoGain()
+//
+// In case that the variables fHiGainLast and fLoGainLast are smaller than 
+// the even part of the number of samples obtained from the run header, a
+// warning is given an the range is set back accordingly. A call to:  
+// - SetRange(fHiGainFirst, fHiGainLast-diff, fLoGainFirst, fLoGainLast) or 
+// - SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast-diff) 
+// is performed in that case. The variable diff means here the difference 
+// between the requested range (fHiGainLast) and the available one. Note that 
+// the functions SetRange() are mostly overloaded and perform more checks, 
+// modifying the ranges again, if necessary.
+//
+// A loop over the MBadPixelsCam is performed and bad pixels are set
+// to MPedestalPix::SetValid(kFALSE);
+//
+Bool_t MPedCalcFromLoGain::ReInit(MParList *pList)
+{
+      
+  Int_t lastdesired   = (Int_t)fLoGainLast;
+  Int_t lastavailable = (Int_t)fRunHeader->GetNumSamplesLoGain()-1;
+  
+  if (lastdesired > lastavailable)
+    {
+      const Int_t diff = lastdesired - lastavailable;
+      *fLog << endl; 
+      *fLog << warn << GetDescriptor()
+            << Form("%s%2i%s%2i%s%2i%s",": Selected Lo Gain FADC Window [",
+                    (int)fLoGainFirst,",",lastdesired,
+                    "] ranges out of the available limits: [0,",lastavailable,"].") << endl;
+      *fLog << GetDescriptor() << ": Will reduce the upper edge to " << (int)(fLoGainLast - diff) << endl;
+      SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast-diff);
+    }
+
+  lastdesired   = (Int_t)fHiGainLast;
+  lastavailable = (Int_t)fRunHeader->GetNumSamplesHiGain()-1;
+  
+  if (lastdesired > lastavailable)
+    {
+      const Int_t diff = lastdesired - lastavailable;
+      *fLog << endl;
+      *fLog << warn << GetDescriptor()
+            << Form("%s%2i%s%2i%s%2i%s",": Selected Hi Gain Range [",
+                    (int)fHiGainFirst,",",lastdesired,
+                    "] ranges out of the available limits: [0,",lastavailable,"].") << endl;
+      *fLog << warn << GetDescriptor() 
+            << Form("%s%2i%s",": Will possibly use ",diff," samples from the Low-Gain for the High-Gain range")
+            << endl;
+      fHiGainLast -= diff;
+      fHiLoLast    = diff;
+    }
+
+  lastdesired   = (Int_t)fHiGainFirst+fWindowSizeHiGain-1;
+  lastavailable = (Int_t)fRunHeader->GetNumSamplesHiGain()-1;
+  
+  if (lastdesired > lastavailable)
+    {
+      const Int_t diff = lastdesired - lastavailable;
+      *fLog << endl;
+      *fLog << warn << GetDescriptor()
+            << Form("%s%2i%s%2i%s",": Selected Hi Gain FADC Window size ",
+                    (int)fWindowSizeHiGain,
+                    " ranges out of the available limits: [0,",lastavailable,"].") << endl;
+      *fLog << warn << GetDescriptor() 
+            << Form("%s%2i%s",": Will use ",diff," samples from the Low-Gain for the High-Gain extraction")
+            << endl;
+
+      if ((Int_t)fWindowSizeHiGain > diff)
+        {
+          fWindowSizeHiGain -= diff;
+          fWindowSizeLoGain += diff;
+        }
+      else
+        {
+          fWindowSizeLoGain += fWindowSizeHiGain;
+          fLoGainFirst       = diff-fWindowSizeHiGain;
+          fWindowSizeHiGain  = 0;
+        }
+    }
+
+
+  Int_t npixels  = fPedestals->GetSize();
+  
+  if (fSumx.GetSize()==0)
+    {
+      fSumx. Set(npixels);
+      fSumx2.Set(npixels);
+      fSumAB0.Set(npixels);
+      fSumAB1.Set(npixels);
+      fNumEventsUsed.Set(npixels);
+      fTotalCounter.Set(npixels);
+      
+      fSumx.Reset();
+      fSumx2.Reset();
+      fSumAB0.Reset();
+      fSumAB1.Reset();
+      fNumEventsUsed.Reset();
+      fTotalCounter.Reset();
+    }
+  
+  if (fWindowSizeHiGain == 0 && fWindowSizeLoGain == 0)
+    {
+      *fLog << err << GetDescriptor() 
+            << ": Number of extracted Slices is 0, abort ... " << endl;
+      return kFALSE;
+    }
+  
+  
+  *fLog << endl;
+  *fLog << inf << GetDescriptor() << ": Taking " << (int)fWindowSizeHiGain
+        << " HiGain FADC samples starting with slice: " << (int)fHiGainFirst << endl;
+  *fLog << inf << GetDescriptor() << ": Taking " << (int)fWindowSizeLoGain
+        << " LoGain FADC samples starting with slice: " << (int)fLoGainFirst << endl;
+  
+  
+  if (fBad)
+    {
+      const Int_t nbads = fBad->GetSize();
+      for (Int_t i=0; i<(nbads>npixels?npixels:nbads);i++)
+        if ((*fBad)[i].IsBad())
+          (*fPedestals)[i].SetValid(kFALSE);
+    }
+  
+  return kTRUE;
+      
+}
+// --------------------------------------------------------------------------
+//
+// Fill the MPedestalCam container with the signal mean and rms for the event.
+// Store the measured signal in arrays fSumx and fSumx2 so that we can 
+// calculate the overall mean and rms in the PostProcess()
+//
+Int_t MPedCalcFromLoGain::Process()
+{
+
+  MRawEvtPixelIter pixel(fRawEvt);
+  
+  //   cout << "fHiGainFirst: " << (Int_t)fHiGainFirst << " fWindowSizeHiGain: " << (Int_t)fWindowSizeHiGain << " fLoGainFirst: " << (Int_t)fLoGainFirst << " fWindowSizeLoGain: " << (Int_t)fWindowSizeLoGain << endl;
+
+
+  while (pixel.Next()) {
+    
+    const UInt_t idx    = pixel.GetPixelId();
+    
+    Byte_t *ptr = pixel.GetHiGainSamples() + fHiGainFirst;
+    Byte_t *end = ptr + fWindowSizeHiGain;
+
+    UInt_t sum = 0;
+    UInt_t sqr = 0;
+    
+    UInt_t max = 0;
+    UInt_t min = 255;
+    
+    // Find the maximum and minimum signal per slice in the high gain window
+    do {
+      if (*ptr > max) {
+	max = *ptr;
+      }
+      if (*ptr < min) {
+	min = *ptr;
+      }
+    } while (++ptr != end);
+    
+    // If the maximum in the high gain window is smaller than 
+    if ((max-min < fMaxHiGainVar) && (max < 255)) {
+      
+      Byte_t *firstSlice = pixel.GetLoGainSamples() + fLoGainFirst;
+      Byte_t *lastSlice  = firstSlice + fWindowSizeLoGain;
+      
+      Byte_t *slice = firstSlice;
+      do {
+	sum += *slice;
+	sqr += *slice * *slice;
+      } while (++slice != lastSlice);
+
+      const Float_t msum = (Float_t)sum;
+      const Float_t sqrsum  = msum*msum;
+      
+      fSumx[idx]  += msum;
+      fSumx2[idx] += sqrsum;
+      fNumEventsUsed[idx]++;
+
+      // Calculate the amplitude of the 150MHz "AB" noise
+
+      Int_t abFlag = (fRunHeader->GetNumSamplesHiGain()
+		       + fLoGainFirst
+		       + pixel.HasABFlag()) & 0x1;
+      for (Int_t islice=0; islice<fWindowSizeLoGain; islice+=2) {
+	Int_t sliceAB0 = islice + abFlag;
+	Int_t sliceAB1 = islice + 1 - abFlag;
+	fSumAB0[idx] += firstSlice[sliceAB0];
+	fSumAB1[idx] += firstSlice[sliceAB1];
+      }
+
+      if (fPedestalUpdate && (fNumEventsUsed[idx] == fNumEventsDump)) {
+
+	const ULong_t n = fNumEventsDump;
+
+	const ULong_t nsamplestot = n*fWindowSizeLoGain;
+	
+	const Float_t sum  = fSumx.At(idx);
+	const Float_t sum2 = fSumx2.At(idx);
+	const Float_t ped  = sum/(nsamplestot);
+    
+	// 1. Calculate the Variance of the sums:
+	Float_t var = (sum2-sum*sum/n)/(n-1.);
+	
+	// 2. Scale the variance to the number of slices:
+	var /= (Float_t)(fWindowSizeLoGain);
+	
+	// 3. Calculate the RMS from the Variance:
+	Float_t rms = var < 0 ? 0. : TMath::Sqrt(var);
+
+	// 4. Calculate the amplitude of the 150MHz "AB" noise
+	Float_t abOffs = (fSumAB0[idx] - fSumAB1[idx]) / nsamplestot;
+
+	(*fPedestals)[idx].Set(ped, rms, abOffs, n);
+	
+	fTotalCounter[idx]++;
+	fNumEventsUsed[idx]=0;
+	fSumx[idx]=0;
+	fSumx2[idx]=0;
+	fSumAB0[idx]=0;
+	fSumAB1[idx]=0;
+      }
+    }
+  }
+  
+  fPedestals->SetReadyToSave();
+  
+  return kTRUE;
+}
+
+// --------------------------------------------------------------------------
+//
+// Compute signal mean and rms in the whole run and store it in MPedestalCam
+//
+Int_t MPedCalcFromLoGain::PostProcess()
+{
+
+  // Compute pedestals and rms from the whole run
+
+  if (!fPedestalUpdate) {
+
+    MRawEvtPixelIter pixel(fRawEvt);
+    
+    while (pixel.Next()) {
+      
+      const Int_t  idx = pixel.GetPixelId();
+      
+      const ULong_t n = fNumEventsUsed[idx];
+      
+      if (n < 2)
+	continue;
+      
+      const ULong_t nsamplestot = n*fWindowSizeLoGain;
+      
+      const Float_t sum  = fSumx.At(idx);
+      const Float_t sum2 = fSumx2.At(idx);
+      const Float_t ped  = sum/(nsamplestot);
+      
+      // 1. Calculate the Variance of the sums:
+      Float_t var = (sum2-sum*sum/n)/(n-1.);
+      
+      // 2. Scale the variance to the number of slices:
+      var /= (Float_t)(fWindowSizeLoGain);
+      
+      // 3. Calculate the RMS from the Variance:
+      Float_t rms = var < 0 ? 0. : TMath::Sqrt(var);
+      
+      // 4. Calculate the amplitude of the 150MHz "AB" noise
+      Float_t abOffs = (fSumAB0[idx] - fSumAB1[idx]) / nsamplestot;
+      
+      (*fPedestals)[idx].Set(ped, rms, abOffs, n);
+      
+      fTotalCounter[idx]++;
+    }
+
+    fPedestals->SetReadyToSave();
+  }
+
+  return kTRUE;
+}
Index: trunk/MagicSoft/Mars/mpedestal/MPedCalcFromLoGain.h
===================================================================
--- trunk/MagicSoft/Mars/mpedestal/MPedCalcFromLoGain.h	(revision 4404)
+++ trunk/MagicSoft/Mars/mpedestal/MPedCalcFromLoGain.h	(revision 4404)
@@ -0,0 +1,69 @@
+#ifndef MARS_MPedCalcFromLoGain
+#define MARS_MPedCalcFromLoGain
+
+#ifndef MARS_MExtractor
+#include "MExtractor.h"
+#endif
+
+#ifndef ROOT_TArrayD
+#include <TArrayD.h>
+#endif
+
+#ifndef ROOT_TArrayI
+#include <TArrayI.h>
+#endif
+
+class MGeomCam;
+class MBadPixelsCam;
+class MPedCalcFromLoGain : public MExtractor
+{
+
+  static const Byte_t fgHiGainFirst;      // First FADC slice Hi-Gain (currently set to: 3) 
+  static const Byte_t fgHiGainLast;       // Last FADC slice Hi-Gain (currently set to: 14) 
+  static const Byte_t fgLoGainFirst;      // First FADC slice Lo-Gain (currently set to: 3) 
+  static const Byte_t fgLoGainLast;       // Last FADC slice Lo-Gain (currently set to: 14) 
+  static const Byte_t fgHiGainWindowSize; // The extraction window Hi-Gain
+  static const Byte_t fgLoGainWindowSize; // The extraction window Lo-Gain
+  static const Byte_t fgMaxHiGainVar;     // The maximum difference between the highest and lowest slice 
+                                          // in the high gain window allowed in order to use low gain 
+                                          // for pedestal calculation
+  Int_t   fNumEventsDump; // Number of event after which MPedestalCam gets updated
+
+  Byte_t  fMaxHiGainVar;
+  Byte_t  fWindowSizeHiGain;             // Number of Hi Gain slices in window
+  Byte_t  fWindowSizeLoGain;             // Number of Lo Gain slices in window  
+
+  Bool_t  fPedestalUpdate;
+  
+  MGeomCam     *fGeom;       // Camera geometry
+  MBadPixelsCam *fBad;       // Bad Pixels
+  
+  TArrayI fNumEventsUsed;   // Number of events used for pedestal calc for each pixel
+  TArrayI fTotalCounter;    // Counter for dumping values to Pedestal Container
+  TArrayD fSumx;            // sum of values
+  TArrayD fSumx2;           // sum of squared values
+  TArrayD fSumAB0;          // sum of ABFlag=0 slices
+  TArrayD fSumAB1;          // sum of ABFlag=1 slices
+  
+  Int_t  PreProcess ( MParList *pList );
+  Bool_t ReInit     ( MParList *pList );
+  Int_t  Process    ();
+  Int_t  PostProcess();
+  
+public:
+
+  MPedCalcFromLoGain(const char *name=NULL, const char *title=NULL);
+  
+  void Clear(const Option_t *o="");
+  void SetRange(Byte_t hifirst=0, Byte_t hilast=0, Byte_t lofirst=0, Byte_t lolast=0);
+  void SetWindowSize(Byte_t windowh=0, Byte_t windowl=0);
+  void SetMaxHiGainVar(Byte_t maxvar=0);
+  void SetDumpEvents(UInt_t dumpevents = 0) {fNumEventsDump = dumpevents;}
+  void SetPedestalUpdate(Bool_t pedupdate)  {fPedestalUpdate = pedupdate;}
+  
+  TArrayI *GetNumEventsUsed() {return &fNumEventsUsed;};
+  
+  ClassDef(MPedCalcFromLoGain, 0)   // Task to calculate pedestals from pedestal runs raw data
+};
+
+#endif
Index: trunk/MagicSoft/Mars/mpedestal/MPedestalCam.cc
===================================================================
--- trunk/MagicSoft/Mars/mpedestal/MPedestalCam.cc	(revision 4384)
+++ trunk/MagicSoft/Mars/mpedestal/MPedestalCam.cc	(revision 4404)
@@ -16,6 +16,7 @@
 !
 !
-!   Author(s): Thomas Bretz  12/2000 <mailto:tbretz@uni-sw.gwdg.de>
-!              Markus Gaug   02/2004 <mailto:markus@ifae.es>
+!   Author(s): Thomas Bretz   12/2000 <mailto:tbretz@uni-sw.gwdg.de>
+!              Markus Gaug    02/2004 <mailto:markus@ifae.es>
+!              Florian Goebel 06/2004 <mailto:fgoebel@mppmu.mpg.de>
 !
 !   Copyright: MAGIC Software Development, 2000-2004
@@ -310,6 +311,5 @@
 }
 
-Bool_t MPedestalCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
-{
+Bool_t MPedestalCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const {
 
   if (GetSize() <= idx)
@@ -319,23 +319,20 @@
     return kFALSE;
 
-  const Float_t ped      = (*this)[idx].GetPedestal();
-  const Float_t rms      = (*this)[idx].GetPedestalRms();
-
-  const Float_t pederr   = rms/TMath::Sqrt((Float_t)fTotalEntries);
-  const Float_t rmserr   = rms/TMath::Sqrt((Float_t)fTotalEntries)/2.;
-
-  switch (type)
-    {
+  switch (type) {
     case 0:
-      val = ped;
+      val = (*this)[idx].GetPedestal();
       break;
     case 1:
-      val = pederr;
+      val = fTotalEntries > 0 ? 
+	  (*this)[idx].GetPedestalRms()/TMath::Sqrt((Float_t)fTotalEntries) 
+	: (*this)[idx].GetPedestalError();
       break;
     case 2:
-      val = rms;
+      val = (*this)[idx].GetPedestalRms();
       break;
     case 3:
-      val = rmserr;
+      val = fTotalEntries > 0 ?
+	  (*this)[idx].GetPedestalRms()/TMath::Sqrt((Float_t)fTotalEntries)/2. 
+	: (*this)[idx].GetPedestalRmsError();
       break;
     default:
Index: trunk/MagicSoft/Mars/mpedestal/MPedestalPix.cc
===================================================================
--- trunk/MagicSoft/Mars/mpedestal/MPedestalPix.cc	(revision 4384)
+++ trunk/MagicSoft/Mars/mpedestal/MPedestalPix.cc	(revision 4404)
@@ -16,6 +16,7 @@
 !
 !
-!   Author(s): Thomas Bretz  12/2000 <mailto:tbretz@uni-sw.gwdg.de>
-!              Markus Gaug   04/2004 <mailto:markus@ifae.es>               
+!   Author(s): Thomas Bretz   12/2000 <mailto:tbretz@uni-sw.gwdg.de>
+!              Markus Gaug    04/2004 <mailto:markus@ifae.es>               
+!              Florian Goebel 06/2004 <mailto:fgoebel@mppmu.mpg.de>
 !
 !   Copyright: MAGIC Software Development, 2000-2004
@@ -31,4 +32,12 @@
 // (offset) value of one Pixel (PMT).                                      //
 //                                                                         //
+// version 2:                                                              //
+// ----------                                                              //
+// added:                                                                  //
+//  fPedestalABoffset   difference between pedestal mean of odd slices and //
+//                      the total pedestal mean (fPedestal)                //
+//  fNumEvents          number of times, the Process was executed          //
+//                      (to estimate the error of pedestal)                //
+//                                                                         //
 /////////////////////////////////////////////////////////////////////////////
 #include "MPedestalPix.h"
@@ -42,5 +51,5 @@
 
 MPedestalPix::MPedestalPix()
-    : fValid(kTRUE)
+  : fValid(kTRUE)
 {
   Clear();
@@ -54,6 +63,8 @@
 void MPedestalPix::Clear(Option_t *o)
 {
-  fPedestal    = -1.;
-  fPedestalRms = -1.;
+  fPedestal         = -1.;
+  fPedestalRms      = -1.;
+  fPedestalABoffset = -1.;
+  fNumEvents        = 0;
 }
 
@@ -61,22 +72,26 @@
 {
 
-  fPedestal    = 0.;
-  fPedestalRms = 0.;
+  fPedestal         = 0.;
+  fPedestalRms      = 0.;
+  fPedestalABoffset = 0.;
+  fNumEvents        = 0;
 }
 
 
-void MPedestalPix::Set(Float_t m, Float_t r)
+void MPedestalPix::Set(Float_t m, Float_t r, Float_t offs, UInt_t n)
 {
-  fPedestal    = m; 
-  fPedestalRms = r; 
+  fPedestal         = m; 
+  fPedestalRms      = r; 
+  fPedestalABoffset = offs;
+  fNumEvents        = n;
 }
 
-Bool_t MPedestalPix::IsValid() const 
+
+Bool_t MPedestalPix::IsValid() const
 {
-
+ 
   if (!fValid)
     return kFALSE;
-
+ 
  return fPedestal>=0||fPedestalRms>=0;
 }
-
Index: trunk/MagicSoft/Mars/mpedestal/MPedestalPix.h
===================================================================
--- trunk/MagicSoft/Mars/mpedestal/MPedestalPix.h	(revision 4384)
+++ trunk/MagicSoft/Mars/mpedestal/MPedestalPix.h	(revision 4404)
@@ -12,5 +12,9 @@
   Float_t fPedestal;     // mean value of pedestal (PMT offset)
   Float_t fPedestalRms;  // root mean square / sigma  / standard deviation of pedestal
-  Bool_t  fValid;        // flag to set pixel valid
+  Float_t fPedestalABoffset; // the difference between odd slice pedestal mean and the 
+                             // total pedestal mean (fPedestal). For even slices pedestal
+                             // use -fPedestalABoffset.
+  UInt_t fNumEvents; // number of times, the Process was executed (to estimate the error of pedestal)
+  Bool_t fValid;   // flag to set pixel valid
   
 public:
@@ -23,8 +27,10 @@
   
     // Setters
-  void SetPedestal(const Float_t f)    { fPedestal = f; }
-  void SetPedestalRms(const Float_t f) { fPedestalRms = f; }
-  
-  void Set(const Float_t m, const Float_t r);
+  void SetPedestal(const Float_t f)         { fPedestal = f; }
+  void SetPedestalRms(const Float_t f)      { fPedestalRms = f; }
+  void SetPedestalABoffset(const Float_t f) { fPedestalABoffset = f; }
+  void SetNumEvents(const UInt_t n)         { fNumEvents = n; }
+
+  void Set(const Float_t m, const Float_t r, const Float_t offs=0, const UInt_t n=0);
   void SetValid(const Bool_t b=kTRUE)    { fValid = b;  }
   
@@ -32,8 +38,11 @@
   Float_t GetPedestal()    const { return fPedestal; }
   Float_t GetPedestalRms() const { return fPedestalRms; }
+  Float_t GetPedestalABoffset() const { return fPedestalABoffset; }
+  Float_t GetPedestalError()    const { return fNumEvents>0 ? fPedestalRms/TMath::Sqrt((Float_t)fNumEvents)   : 0; }
+  Float_t GetPedestalRmsError() const { return fNumEvents>0 ? fPedestalRms/TMath::Sqrt((Float_t)fNumEvents/2) : 0; }
   
   Bool_t IsValid()         const;
 
-  ClassDef(MPedestalPix, 1)	// Storage Container for Pedestal information of one pixel
+  ClassDef(MPedestalPix, 2)	// Storage Container for Pedestal information of one pixel
 };
 
Index: trunk/MagicSoft/Mars/mpedestal/Makefile
===================================================================
--- trunk/MagicSoft/Mars/mpedestal/Makefile	(revision 4384)
+++ trunk/MagicSoft/Mars/mpedestal/Makefile	(revision 4404)
@@ -32,6 +32,7 @@
            MMcPedestalNSBAdd.cc \
            MPedCalcPedRun.cc \
+	   MPedCalcFromLoGain.cc \
            MPedPhotCalc.cc \
- 	   MPedPhotCam.cc \
+           MPedPhotCam.cc \
            MPedPhotPix.cc \
            MPedestalCam.cc \
Index: trunk/MagicSoft/Mars/mpedestal/PedestalLinkDef.h
===================================================================
--- trunk/MagicSoft/Mars/mpedestal/PedestalLinkDef.h	(revision 4384)
+++ trunk/MagicSoft/Mars/mpedestal/PedestalLinkDef.h	(revision 4404)
@@ -13,4 +13,5 @@
 
 #pragma link C++ class MPedCalcPedRun++;
+#pragma link C++ class MPedCalcFromLoGain++;
 #pragma link C++ class MPedestalCam++;
 #pragma link C++ class MPedestalPix++;
