Index: trunk/MagicSoft/Mars/Changelog
===================================================================
--- trunk/MagicSoft/Mars/Changelog	(revision 3700)
+++ trunk/MagicSoft/Mars/Changelog	(revision 3701)
@@ -38,4 +38,5 @@
 
    * manalysis/MPedestalCam.[h,cc]
+   * manalysis/MPedCalcPedRun.[h,cc]
    * manalysis/MGeomApply.cc
      - added average pixels in the way like it is done in MCalibrationCam
Index: trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.cc
===================================================================
--- trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.cc	(revision 3700)
+++ trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.cc	(revision 3701)
@@ -96,4 +96,6 @@
 #include "MExtractedSignalCam.h"
 
+#include "MGeomPix.h"
+#include "MGeomCam.h"
 
 #include "MGeomCamMagic.h"
@@ -142,18 +144,26 @@
 Int_t MPedCalcPedRun::PreProcess( MParList *pList )
 {
-    Clear();
-
-    fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
-    if (!fRawEvt)
-    {
-        *fLog << err << "MRawEvtData not found... aborting." << endl;
-        return kFALSE;
-    }
-
-    fPedestals = (MPedestalCam*)pList->FindCreateObj("MPedestalCam");
-    if (!fPedestals)
-        return kFALSE;
-
-    return kTRUE;
+
+  Clear();
+  
+  fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
+  if (!fRawEvt)
+    {
+      *fLog << err << "MRawEvtData 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;
+  
+  return kTRUE;
 }
 
@@ -175,13 +185,23 @@
     }
     else
-        if (runheader->IsMonteCarloRun())
-            return kTRUE;
-
-    Int_t n = fPedestals->GetSize();
+      if (runheader->IsMonteCarloRun())
+        return kTRUE;
+
+    Int_t npixels  = fPedestals->GetSize();
+    Int_t areas    = fPedestals->GetAverageAreas();
+    Int_t sectors  = fPedestals->GetAverageSectors();
 
     if (fSumx.GetSize()==0)
     {
-	fSumx.Set(n);
-	fSumx2.Set(n);
+	fSumx. Set(npixels);
+	fSumx2.Set(npixels);
+
+	fAreaSumx. Set(areas);
+	fAreaSumx2.Set(areas);
+	fAreaValid.Set(areas);
+
+	fSectorSumx. Set(sectors);
+	fSectorSumx2.Set(sectors);
+	fSectorValid.Set(sectors);
 
 	fSumx.Reset();
@@ -204,50 +224,59 @@
 Int_t MPedCalcPedRun::Process()
 {
-    MRawEvtPixelIter pixel(fRawEvt);
-
-    while (pixel.Next())
-    {
-              Byte_t *ptr = pixel.GetHiGainSamples();
-        const Byte_t *end = ptr + fNumHiGainSamples;
-
-	UInt_t sum = 0;
-	UInt_t sqr = 0;
+
+  MRawEvtPixelIter pixel(fRawEvt);
+  
+  while (pixel.Next())
+    {
+
+      const UInt_t idx    = pixel.GetPixelId();
+      const UInt_t aidx   = (*fGeom)[idx].GetAidx();
+      const UInt_t sector = (*fGeom)[idx].GetSector();      
+
+      Byte_t *ptr = pixel.GetHiGainSamples();
+      const Byte_t *end = ptr + fNumHiGainSamples;
+      
+      UInt_t sum = 0;
+      UInt_t sqr = 0;
 	
-	do
-	  {
-	    sum += *ptr;
-	    sqr += *ptr * *ptr;
-	  }
-	while (++ptr != end);
-	
-	const Float_t msum = (Float_t)sum;
-
-	const UInt_t idx = pixel.GetPixelId();
-        //
-        // These three lines have been uncommented by Markus Gaug
-        // If anybody needs them, please contact me!!
-        //
-        //	const Float_t higainped = msum/fNumHiGainSamples;
-        //	const Float_t higainrms = TMath::Sqrt((msqr-msum*msum/fNumHiGainSamples)/(fNumHiGainSamples-1.));
-        //	(*fPedestals)[idx].Set(higainped, higainrms);
-	
-	fSumx[idx]  += msum;
-        //
-        // The old version:
-        //
-	//       const Float_t msqr = (Float_t)sqr;
-        //	fSumx2[idx] += msqr;
-        //
-        // The new version:
-        //
-	fSumx2[idx] += msum*msum;
-
-    }
-    
-    fPedestals->SetReadyToSave();
-    fNumSamplesTot += fNumHiGainSamples;
-
-    
-    return kTRUE;
+      do
+        {
+          sum += *ptr;
+          sqr += *ptr * *ptr;
+        }
+      while (++ptr != end);
+      
+      const Float_t msum = (Float_t)sum;
+      
+      //
+      // These three lines have been uncommented by Markus Gaug
+      // If anybody needs them, please contact me!!
+      //
+      //	const Float_t higainped = msum/fNumHiGainSamples;
+      //	const Float_t higainrms = TMath::Sqrt((msqr-msum*msum/fNumHiGainSamples)/(fNumHiGainSamples-1.));
+      //	(*fPedestals)[idx].Set(higainped, higainrms);
+      
+      fSumx[idx]          += msum;
+      fAreaSumx[aidx]     += msum;
+      fSectorSumx[sector] += msum;      
+      //
+      // The old version:
+      //
+      //       const Float_t msqr = (Float_t)sqr;
+      //	fSumx2[idx] += msqr;
+      //
+      // The new version:
+      //
+      const Float_t sqrsum  = msum*msum;
+      fSumx2[idx]          += sqrsum;
+      fAreaSumx2[aidx]     += sqrsum;
+      fSectorSumx2[sector] += sqrsum;      
+    }
+  
+  fPedestals->SetReadyToSave();
+  fNumSamplesTot += fNumHiGainSamples;
+  
+  
+  return kTRUE;
 }
 
@@ -267,6 +296,12 @@
   while (pixel.Next())
     {
-      const Int_t pixid = pixel.GetPixelId();
-      
+
+      const Int_t  pixid  = pixel.GetPixelId();
+      const UInt_t aidx   = (*fGeom)[pixid].GetAidx();
+      const UInt_t sector = (*fGeom)[pixid].GetSector();      
+      
+      fAreaValid  [aidx]++;
+      fSectorValid[sector]++;
+
       const Float_t sum  = fSumx.At(pixid);
       const Float_t sum2 = fSumx2.At(pixid);
@@ -290,4 +325,58 @@
 
     }
+
+  //
+  // Loop over the (two) area indices to get the averaged pedestal per aidx
+  //
+  for (Int_t aidx=0; aidx<fAreaValid.GetSize(); aidx++)
+    {
+      
+      const Int_t   napix = fAreaValid.At(aidx);
+      const Float_t sum   = fAreaSumx.At(aidx);
+      const Float_t sum2  = fAreaSumx2.At(aidx);
+      const ULong_t an    = napix * n;
+      const ULong_t aevts = napix * nevts;
+      
+      const Float_t higainped = sum/an;
+
+      // 1. Calculate the Variance of the sums:
+      Float_t higainVar = (sum2-sum*sum/aevts)/(aevts-1.);
+      // 2. Scale the variance to the number of slices:
+      higainVar /= (Float_t)fNumHiGainSamples;
+      // 3. Calculate the RMS from the Variance:
+      Float_t higainrms = TMath::Sqrt(higainVar);
+      // 4. Re-scale it with the square root of the number of involved pixels 
+      //    in order to be comparable to the mean of pedRMS of that area
+      higainrms *= TMath::Sqrt((Float_t)napix);
+
+      fPedestals->GetAverageArea(aidx).Set(higainped, higainrms);
+    }
+  
+  //
+  // Loop over the (six) sector indices to get the averaged pedestal per sector
+  //
+  for (Int_t sector=0; sector<fSectorValid.GetSize(); sector++)
+    {
+      
+      const Int_t   nspix = fSectorValid.At(sector);
+      const Float_t sum   = fSectorSumx.At(sector);
+      const Float_t sum2  = fSectorSumx2.At(sector);
+      const ULong_t sn    = nspix * n;
+      const ULong_t sevts = nspix * nevts;
+      
+      const Float_t higainped = sum/sn;
+
+      // 1. Calculate the Variance of the sums:
+      Float_t higainVar = (sum2-sum*sum/sevts)/(sevts-1.);
+      // 2. Scale the variance to the number of slices:
+      higainVar /= (Float_t)fNumHiGainSamples;
+      // 3. Calculate the RMS from the Variance:
+      Float_t higainrms = TMath::Sqrt(higainVar);
+      // 4. Re-scale it with the square root of the number of involved pixels 
+      //    in order to be comparable to the mean of pedRMS of that sector
+      higainrms *= TMath::Sqrt((Float_t)nspix);
+
+      fPedestals->GetAverageSector(sector).Set(higainped, higainrms);
+    }
   
   fPedestals->SetTotalEntries(fNumSamplesTot);
Index: trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.h
===================================================================
--- trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.h	(revision 3700)
+++ trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.h	(revision 3701)
@@ -10,31 +10,43 @@
 #endif
 
+#ifndef ROOT_TArrayI
+#include <TArrayI.h>
+#endif
+
 class MRawEvtData;
 class MPedestalCam;
-
+class MGeomCam;
 class MPedCalcPedRun : public MTask
 {
-    Byte_t fNumHiGainSamples;
-    UInt_t fNumSamplesTot;
 
-    MRawEvtData  *fRawEvt;     // raw event data (time slices)
-    MPedestalCam *fPedestals;  // Pedestals of all pixels in the camera
+  Byte_t fNumHiGainSamples;
+  UInt_t fNumSamplesTot;
+  
+  MRawEvtData  *fRawEvt;     // raw event data (time slices)
+  MPedestalCam *fPedestals;  // Pedestals of all pixels in the camera
+  MGeomCam     *fGeom;       // Camera geometry
+  
+  TArrayF fSumx;         // sum of values
+  TArrayF fSumx2;        // sum of squared values
+  TArrayF fAreaSumx;     // averaged sum of values per area idx
+  TArrayF fAreaSumx2;    // averaged sum of squared values per area idx
+  TArrayI fAreaValid;    // number of valid pixel with area idx  
+  TArrayF fSectorSumx;   // averaged sum of values per sector 
+  TArrayF fSectorSumx2;  // averaged sum of squared values per sector
+  TArrayI fSectorValid;  // number of valid pixel with sector idx  
+  
+  Int_t  PreProcess ( MParList *pList );
+  Bool_t ReInit     ( MParList *pList );
+  Int_t  Process    ();
+  Int_t  PostProcess();
+  
+public:
 
-    TArrayF fSumx;   // sum of values
-    TArrayF fSumx2;  // sum of squared values
-
-    Bool_t ReInit(MParList *pList);
-
-    Int_t PreProcess(MParList *pList);
-    Int_t Process();
-    Int_t PostProcess();
-
-public:
-    MPedCalcPedRun(const char *name=NULL, const char *title=NULL);
-
-    void Clear(const Option_t *o="");
-    void SetNumHiGainSamples(const Byte_t n)      { fNumHiGainSamples = n;   }
-    
-    ClassDef(MPedCalcPedRun, 0)   // Task to calculate pedestals from pedestal runs raw data
+  MPedCalcPedRun(const char *name=NULL, const char *title=NULL);
+  
+  void Clear(const Option_t *o="");
+  void SetNumHiGainSamples(const Byte_t n)      { fNumHiGainSamples = n;   }
+  
+  ClassDef(MPedCalcPedRun, 0)   // Task to calculate pedestals from pedestal runs raw data
 };