Index: trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.cc
===================================================================
--- trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.cc	(revision 2818)
+++ trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.cc	(revision 2821)
@@ -18,6 +18,8 @@
 !   Author(s): Josep Flix 04/2001 <mailto:jflix@ifae.es>
 !   Author(s): Thomas Bretz 05/2001 <mailto:tbretz@uni-sw.gwdg.de>
+!   Author(s): Sebastian Commichau 12/2003 
+!   Author(s): Javier Rico 01/2004 
 !
-!   Copyright: MAGIC Software Development, 2000-2001
+!   Copyright: MAGIC Software Development, 2000-2004
 !
 !
@@ -28,4 +30,10 @@
 //   MPedCalcPedRun                                                        //
 //                                                                         //
+//  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                                 //
+//                                                                         //
 //  Input Containers:                                                      //
 //   MRawEvtData                                                           //
@@ -33,4 +41,5 @@
 //  Output Containers:                                                     //
 //   MPedestalCam                                                          //
+//                                                                         //
 //                                                                         //
 /////////////////////////////////////////////////////////////////////////////
@@ -70,5 +79,5 @@
     if (!fRawEvt)
     {
-        *fLog << dbginf << "MRawEvtData not found... aborting." << endl;
+        *fLog << err << "MRawEvtData not found... aborting." << endl;
         return kFALSE;
     }
@@ -78,69 +87,75 @@
         return kFALSE;
 
-    MGeomCamMagic magiccam;
-
-    fSumx.Set(magiccam.GetNumPixels());
-    fSumx2.Set(magiccam.GetNumPixels());
-    
-    for(UInt_t i=0;i<magiccam.GetNumPixels();i++)
-    {
-	fSumx.AddAt(0,i);
-	fSumx2.AddAt(0,i);
-    }
+    fNumSamplesTot=0;
 
     return kTRUE;
 }
 
-Bool_t MPedCalcPedRun::ReInit(MParList *pList )   
+Bool_t MPedCalcPedRun::ReInit(MParList *pList)
 {
+    const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
+    if (!runheader)
+    {
+        *fLog << warn << dbginf;
+        *fLog << "Warning - cannot check file type, MRawRunHeader not found." << endl;
+    }
+    else
+        if (runheader->GetRunType() == kRTMonteCarlo)
+            return kTRUE;
 
-    fRunheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
-    if (!fRunheader)
-	{
-        *fLog << warn << dbginf << 
-        	"Warning - cannot check file type, MRawRunHeader not found." << endl;
-	}
-    else
-        if (fRunheader->GetRunType() == kRTMonteCarlo)
-        {
-            return kTRUE;
-        }
+    fNumPixels = fPedestals->GetSize();
 
-    fNumHiGainSamples =  fRunheader->GetNumSamplesHiGain();
+    if(fSumx.GetSize()==0)
+      {
+	fSumx.Set(fNumPixels);
+	fSumx2.Set(fNumPixels);
 
-    fPedestals->InitSize(fRunheader->GetNumPixel());
+	memset(fSumx.GetArray(),  0, sizeof(Float_t)*fNumPixels);
+	memset(fSumx2.GetArray(), 0, sizeof(Float_t)*fNumPixels);
+      }
+
+    // Calculate an even number for the hi gain samples to avoid
+    // biases due to the fluctuation in pedestal from one slice to
+    // the other one
+    fNumHiGainSamples = runheader->GetNumSamplesHiGain() & ~1;
 
     return kTRUE;
-
 }
-
 
 Int_t MPedCalcPedRun::Process()
 {
-
     MRawEvtPixelIter pixel(fRawEvt);
 
     while (pixel.Next())
     {
-	Byte_t shift=(fNumHiGainSamples/2*2==fNumHiGainSamples) ? 0:1;
               Byte_t *ptr = pixel.GetHiGainSamples();
-        const Byte_t *end = ptr + fRawEvt->GetNumHiGainSamples()-shift;
+        const Byte_t *end = ptr + fNumHiGainSamples;
 
-        const Float_t higainped = CalcHiGainMean(ptr, end);
-        const Float_t higainrms = CalcHiGainRms(ptr, end, higainped);
-
-        const UInt_t pixid = pixel.GetPixelId();
-        MPedestalPix &pix = (*fPedestals)[pixid];
+	UInt_t sum = 0;
+	UInt_t sqr = 0;
 	
-	// cumulate the sum of pedestals and of pedestal squares
-	fSumx.AddAt(higainped+fSumx.At(pixid),pixid);
-	fSumx2.AddAt(GetSumx2(ptr, end)+fSumx2.At(pixid),pixid);
-
-	// set the value of the pedestal and rms computed from the processed event
-        pix.Set(higainped, higainrms);
+	do
+	  {
+	    sum += *ptr;
+	    sqr += *ptr * *ptr;
+	  }
+	while (++ptr != end);
+	
+	const Float_t msum = (Float_t)sum;
+	const Float_t msqr = (Float_t)sqr;
+	
+	const Float_t higainped = msum/fNumHiGainSamples;
+	const Float_t higainrms = TMath::Sqrt((msqr-msum*msum/fNumHiGainSamples)/(fNumHiGainSamples-1.));
+	
+	const UInt_t idx = pixel.GetPixelId();
+	(*fPedestals)[idx].Set(higainped, higainrms);
+	
+	fSumx[idx]  += msum;
+	fSumx2[idx] += sqr;
     }
-
+    
     fPedestals->SetReadyToSave();
-
+    fNumSamplesTot+=fNumHiGainSamples;
+    
     return kTRUE;
 }
@@ -148,84 +163,24 @@
 Int_t MPedCalcPedRun::PostProcess()
 {
-// Compute pedestals and rms from the whole run
-
-    MRawEvtPixelIter pixel(fRawEvt);
-
-    while (pixel.Next())
+  // Compute pedestals and rms from the whole run
+  const Int_t n  = fNumSamplesTot;
+  
+  MRawEvtPixelIter pixel(fRawEvt);
+  
+  while (pixel.Next())
     {
-        const UInt_t pixid = pixel.GetPixelId();
-        MPedestalPix &pix = (*fPedestals)[pixid];
-
-	const Int_t N = GetNumExecutions();
-	const Float_t sum = fSumx.At(pixid);
-	const Float_t sum2 = fSumx2.At(pixid);
-	const Float_t higainped = sum/N;
-        const Float_t higainrms = sqrt(1./(N-1.)*(sum2-sum*sum/N));
-        pix.Set(higainped,higainrms);
-
+      const UInt_t pixid = pixel.GetPixelId();
+      
+      const Float_t sum  = fSumx.At(pixid);
+      const Float_t sum2 = fSumx2.At(pixid);
+      
+      const Float_t higainped = sum/n;
+      const Float_t higainrms = TMath::Sqrt((sum2-sum*sum/n)/(n-1.));
+      
+      (*fPedestals)[pixid].Set(higainped, higainrms);
     }
-    return kTRUE;
-
+  
+  return kTRUE;
 }
 
 
-Float_t MPedCalcPedRun::CalcHiGainMean(Byte_t *ptr, const Byte_t *end) const
-{
-    Int_t sum=0;
-    Byte_t EvenNumSamples=(fNumHiGainSamples/2*2==fNumHiGainSamples) 
-                        ? fNumHiGainSamples : fNumHiGainSamples-1;
-
-    do sum += *ptr;
-    while (++ptr != end);
-   
-    return (Float_t)sum/EvenNumSamples;
-}
-
-
-
-Float_t MPedCalcPedRun::GetSumx2(Byte_t *ptr, const Byte_t *end) const
-{
-    Float_t square = 0;
-
-    // Take an even number of time slices to avoid biases due to A/B effect
-    Byte_t EvenNumSamples=(fNumHiGainSamples/2*2==fNumHiGainSamples) ? fNumHiGainSamples:fNumHiGainSamples-1;
-
-    do
-    {
-        const Float_t val = (Float_t)(*ptr);
-
-        square += val*val;
-    } while (++ptr != end);
-
-    return square/EvenNumSamples;
-}
-
-
-Float_t MPedCalcPedRun::CalcHiGainRms(Byte_t *ptr, const Byte_t *end, Float_t higainped) const
-{
-    Float_t rms = 0;
-    Byte_t EvenNumSamples=(fNumHiGainSamples/2*2==fNumHiGainSamples) ? fNumHiGainSamples:fNumHiGainSamples-1;
-
-    do
-    {
-        const Float_t diff = (Float_t)(*ptr)-higainped;
-
-        rms += diff*diff;
-    } while (++ptr != end);
-
-    return sqrt(rms/(EvenNumSamples-1));
-}
-
-
-
-/*
-Float_t MPedCalcPedRun::CalcHiGainMeanErr(Float_t higainrms) const
-{
-    return higainrms/sqrt((float)fNumHiGainSamples);
-}
-
-Float_t MPedCalcPedRun::CalcHiGainRmsErr(Float_t higainrms) const
-{
-    return higainrms/sqrt(2.*fNumHiGainSamples);
-}
-*/
Index: trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.h
===================================================================
--- trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.h	(revision 2818)
+++ trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.h	(revision 2821)
@@ -16,5 +16,4 @@
 #include <TArrayF.h>
 
-class MRawRunHeader;
 class MRawEvtData;
 class MPedestalCam;
@@ -22,7 +21,8 @@
 class MPedCalcPedRun : public MTask
 {
-    Byte_t fNumHiGainSamples;
+    Byte_t   fNumHiGainSamples;
+    UShort_t fNumPixels;
+    ULong_t   fNumSamplesTot;
 
-    MRawRunHeader *fRunheader; // raw event run header
     MRawEvtData  *fRawEvt;     // raw event data (time slices)
     MPedestalCam *fPedestals;  // Pedestals of all pixels in the camera
@@ -30,10 +30,4 @@
     TArrayF fSumx;   // sum of values
     TArrayF fSumx2;  // sum of squared values
-
-    Float_t CalcHiGainMean(Byte_t *ptr, const Byte_t *end) const;
-    Float_t CalcHiGainRms(Byte_t *ptr, const Byte_t *end, Float_t higainped) const;
-    Float_t GetSumx2(Byte_t* ptr, const Byte_t* end) const; 
-   //Float_t CalcHiGainMeanErr(Float_t higainrms) const;
-    //Float_t CalcHiGainRmsErr(Float_t higainrms) const;
 
     Bool_t ReInit(MParList *pList);
@@ -44,5 +38,4 @@
 
 public:
-
     MPedCalcPedRun(const char *name=NULL, const char *title=NULL);
 
