Changeset 4609 for trunk/MagicSoft/Mars/mpedestal

Timestamp:

08/13/04 14:59:17 (20 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mpedestal

Files:

• MPedCalcFromLoGain.cc (modified) (12 diffs)
• MPedCalcFromLoGain.h (modified) (1 diff)
• MPedCalcPedRun.cc (modified) (9 diffs)
• MPedCalcPedRun.h (modified) (2 diffs)
• MPedPhotCam.cc (modified) (3 diffs)
• MPedPhotCam.h (modified) (1 diff)
• MPedestalCam.cc (modified) (3 diffs)
• MPedestalCam.h (modified) (2 diffs)
• MPedestalPix.h (modified) (1 diff)

trunk/MagicSoft/Mars/mpedestal/MPedCalcFromLoGain.cc

@@ -185 +185 @@
       fGeom(NULL), fPedContainerName("MPedestalCam")
 {
-  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();
+    fName  = name  ? name  : "MPedCalcFromLoGain";
+    fTitle = title ? title : "Task to calculate pedestals from pedestal runs raw data";
+
+    AddToBranchList("fHiGainPixId");
+    AddToBranchList("fLoGainPixId");
+    AddToBranchList("fHiGainFadcSamples");
+    AddToBranchList("fLoGainFadcSamples");
+
+    SetRange();
+    SetMaxHiGainVar();
+    SetPedestalUpdate(kTRUE);
+
+    Clear();
 }
 
@@ -208 +210 @@
 void MPedCalcFromLoGain::Clear(const Option_t *o)
 {
-  fRawEvt    = NULL;
-  fRunHeader = NULL;
-  fPedestals = NULL;
-
-  // If the size is yet set, set the size
-  if (fSumx.GetSize()>0)
-  {
-      // Reset contents of arrays.
-      fSumx.Reset();
-      fSumx2.Reset();
-      fSumAB0.Reset();
-      fSumAB1.Reset();
-      fNumEventsUsed.Reset();
-      fTotalCounter.Reset();
-  }
+    fRawEvt    = NULL;
+    fRunHeader = NULL;
+    fPedestals = NULL;
+
+    // If the size is yet set, set the size
+    if (fSumx.GetSize()>0)
+    {
+        // Reset contents of arrays.
+        fSumx.Reset();
+        fSumx2.Reset();
+        fSumAB0.Reset();
+        fSumAB1.Reset();
+        fNumEventsUsed.Reset();
+        fTotalCounter.Reset();
+    }
 }
 
@@ -235 +237 @@
 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);
+    MExtractor::SetRange(hifirst, hilast, lofirst, lolast);
+
+    //
+    // Redo the checks if the window is still inside the ranges
+    //
+    SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain);
 }
 
@@ -247 +249 @@
 void MPedCalcFromLoGain::SetMaxHiGainVar(Byte_t maxvar)
 {
-  fMaxHiGainVar = maxvar;
+    fMaxHiGainVar = maxvar;
 }
 
@@ -265 +267 @@
 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;
+    fWindowSizeHiGain = windowh & ~1;
+    fWindowSizeLoGain = windowl & ~1;
   
-  if (fWindowSizeLoGain != windowl)
-    *fLog << warn << GetDescriptor() << ": Lo Gain window size has to be even, set to: " 
-          << int(fWindowSizeLoGain) << " samples " << endl;
+    if (fWindowSizeHiGain != windowh)
+    {
+        *fLog << warn;
+        *fLog << GetDescriptor() << ": HiGain window has to be even, set to: ";
+        *fLog << int(fWindowSizeHiGain) << " 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(": Hi Gain window size: %2i is bigger than available range: [%2i,%2i]",
-                    (int)fWindowSizeHiGain, (int)fHiGainFirst, (int)fHiGainLast) << endl;
-      *fLog << warn << GetDescriptor() 
-            << ": Will set window size to: " << (int)availhirange << endl;
-      fWindowSizeHiGain = availhirange;
-    }
+    if (fWindowSizeLoGain != windowl)
+    {
+        *fLog << warn;
+        *fLog << GetDescriptor() << ": Lo Gain window has to be even, set to: ";
+        *fLog << int(fWindowSizeLoGain) << " samples " << endl;
+    }
+      
+    const Byte_t availhirange = (fHiGainLast-fHiGainFirst+1) & ~1;
+    const Byte_t availlorange = (fLoGainLast-fLoGainFirst+1) & ~1;
   
-  if (fWindowSizeLoGain > availlorange)
-    {
-      *fLog << warn << GetDescriptor() 
-            << Form(": Lo Gain window size: %2i is bigger than available range: [%2i,%2i]",
-                    (int)fWindowSizeLoGain, (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);
+    if (fWindowSizeHiGain > availhirange)
+    {
+        *fLog << warn;
+        *fLog << GetDescriptor() << ": HiGain window " << (int)fWindowSizeHiGain;
+        *fLog << " out of range [" << (int)fHiGainFirst;
+        *fLog << "," << (int)fHiGainLast << "]" << endl;
+        *fLog << "Will set window size to " << (int)availhirange << endl;
+        fWindowSizeHiGain = availhirange;
+    }
+    
+    if (fWindowSizeLoGain > availlorange)
+    {
+        *fLog << warn;
+        *fLog << GetDescriptor() << ": LoGain window " << (int)fWindowSizeLoGain;
+        *fLog << " out of range [" << (int)fLoGainFirst;
+        *fLog << "," << (int)fLoGainLast << "]" << endl;
+        *fLog << "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);
+    */
 }
 
@@ -369 +378 @@
 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(": Selected Lo Gain FADC Window [%2i,%2i] ranges out of the available limits: [0,%2i].",
-                    (int)fLoGainFirst, lastdesired, 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(": Selected Hi Gain Range [%2i,%2i] ranges out of the available limits: [0,%2i].",
-                    (int)fHiGainFirst, lastdesired, lastavailable) << endl;
-      *fLog << warn << GetDescriptor() 
-            << Form(": Will possibly use %2i samples from the Low-Gain for the High-Gain range", diff)
-            << 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(": Selected Hi Gain FADC Window size %2i ranges out of the available limits: [0,%2i].",
-                    (int)fWindowSizeHiGain, lastavailable) << endl;
-      *fLog << warn << GetDescriptor() 
-            << Form(": Will use %2i samples from the Low-Gain for the High-Gain extraction", diff)
-            << endl;
-
-      if ((Int_t)fWindowSizeHiGain > diff)
-        {
-          fWindowSizeHiGain -= diff;
-          fWindowSizeLoGain += diff;
-        }
-      else
-        {
-          fWindowSizeLoGain += fWindowSizeHiGain;
-          fLoGainFirst       = diff-fWindowSizeHiGain;
-          fWindowSizeHiGain  = 0;
-        }
-    }
-
-
-  // If the size is not yet set, set the size
-  if (fSumx.GetSize()==0)
-  {
-      const Int_t npixels = fPedestals->GetSize();
-
-      fSumx. Set(npixels);
-      fSumx2.Set(npixels);
-      fSumAB0.Set(npixels);
-      fSumAB1.Set(npixels);
-      fNumEventsUsed.Set(npixels);
-      fTotalCounter.Set(npixels);
-
-      // Reset contents of arrays.
-      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;
-  
-  return kTRUE;
+    if (fRunHeader->GetNumSamplesHiGain()<1)
+    {
+        *fLog << err << "ERROR - Number of available HiGainSamples<1... abort." << endl;
+        return kFALSE;
+    }
+    if (fRunHeader->GetNumSamplesLoGain()<1)
+    {
+        *fLog << err << "ERROR - Number of available LoGainSamples<1... abort." << endl;
+        return kFALSE;
+    }
+
+    fHiGainFirst = TMath::Max(fHiGainFirst, 0);
+    fHiGainLast  = TMath::Min(fHiGainLast,  fRunHeader->GetNumSamplesHiGain()-1);
+
+    fLoGainFirst = TMath::Max(fLoGainFirst, 0);
+    fLoGainLast  = TMath::Min(fLoGainLast,  fRunHeader->GetNumSamplesLoGain()-1);
+
+    fWindowSizeHiGain = TMath::Min(fWindowSizeHiGain, fHiGainLast-fHiGainFirst+1);
+    fWindowSizeLoGain = TMath::Min(fWindowSizeLoGain, fLoGainLast-fLoGainFirst+1);
+
+    SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast);
+
+    *fLog << inf << endl;
+    *fLog << "Taking " << Form("%2d", (int)fWindowSizeHiGain) << " HiGain from " << (int)fHiGainFirst << endl;
+    *fLog << "Taking " << Form("%2d", (int)fWindowSizeLoGain) << " LoGain from " << (int)fLoGainFirst << endl;
+
+    if (fWindowSizeHiGain==0)
+    {
+        *fLog << err << "ERROR - HiGain windows size == 0... abort." << endl;
+        return kFALSE;
+    }
+    if (fWindowSizeLoGain==0)
+    {
+        *fLog << err << "ERROR - HiGain windows size == 0... abort." << endl;
+        return kFALSE;
+    }
+
+    // If the size is not yet set, set the size
+    if (fSumx.GetSize()==0)
+    {
+        const Int_t npixels = fPedestals->GetSize();
+
+        fSumx. Set(npixels);
+        fSumx2.Set(npixels);
+        fSumAB0.Set(npixels);
+        fSumAB1.Set(npixels);
+        fNumEventsUsed.Set(npixels);
+        fTotalCounter.Set(npixels);
+
+        // Reset contents of arrays.
+        fSumx.Reset();
+        fSumx2.Reset();
+        fSumAB0.Reset();
+        fSumAB1.Reset();
+        fNumEventsUsed.Reset();
+        fTotalCounter.Reset();
+    }
+
+    return kTRUE;
 }
 
@@ -515 +489 @@
         // Find the maximum and minimum signal per slice in the high gain window
         do {
-            if (*ptr > max) {
+            if (*ptr > max)
                 max = *ptr;
-            }
-            if (*ptr < min) {
+            if (*ptr < min)
                 min = *ptr;
-            }
         } while (++ptr != end);
 
@@ -527 +499 @@
             continue;
 
-        Byte_t *firstSlice = pixel.GetLoGainSamples() + fLoGainFirst;
-        Byte_t *lastSlice  = firstSlice + fWindowSizeLoGain;
-
-        Byte_t *slice = firstSlice;
+        ptr = pixel.GetLoGainSamples() + fLoGainFirst;
+        end = ptr + fWindowSizeLoGain;
+
+        Byte_t *firstSlice = ptr;
+
         do {
-            sum += *slice;
-            sqr += *slice * *slice;
-        } while (++slice != lastSlice);
+            sum += *ptr;
+            sqr += *ptr * *ptr;
+        } while (++ptr != end);
 
         const Float_t msum   = (Float_t)sum;
@@ -569 +542 @@
 
     if (fPedestalUpdate)
+    {
+        fPedestals->ReCalc(*fGeom);
         fPedestals->SetReadyToSave();
+    }
 
     return kTRUE;
@@ -581 +557 @@
 {
     // Compute pedestals and rms from the whole run
-    if (fPedestalUpdate)
+    if (fPedestalUpdate || GetNumExecutions()<1)
         return kTRUE;
 
-    *fLog << flush << inf << "Calculating Pedestals..." << flush;
+    *fLog << flush << inf << "Calculating pedestals..." << flush;
+
+    Double_t sum = 0;
 
     const Int_t npix = fGeom->GetNumPixels();
@@ -592 +570 @@
         if (n>1)
             Calc(n, idx);
-    }
-
+        sum += n;
+    }
+
+    *fLog << flush << inf << "Calculating means..." << flush;
+
+    fPedestals->SetTotalEntries((UInt_t)(sum/npix*(fWindowSizeLoGain+fWindowSizeHiGain)));
+    fPedestals->ReCalc(*fGeom);
     fPedestals->SetReadyToSave();
     return kTRUE;
@@ -621 +604 @@
         SetWindowSize(hw, lw);
 
-    Int_t num = fNumEventsDump;
     if (IsEnvDefined(env, prefix, "NumEventsDump", print))
     {
-        num = GetEnvValue(env, prefix, "NumEventsDump", num);
+        SetDumpEvents(GetEnvValue(env, prefix, "NumEventsDump", fNumEventsDump));
         rc = kTRUE;
     }
-    SetDumpEvents(num);
-
-    Byte_t max = fMaxHiGainVar;
+
     if (IsEnvDefined(env, prefix, "MaxHiGainVar", print))
     {
-        max = GetEnvValue(env, prefix, "MaxHiGainVar", max);
+        SetMaxHiGainVar(GetEnvValue(env, prefix, "MaxHiGainVar", fMaxHiGainVar));
         rc = kTRUE;
     }
-    SetMaxHiGainVar(max);
-
-    Bool_t upd = fPedestalUpdate;
+
     if (IsEnvDefined(env, prefix, "PedestalUpdate", print))
     {
-        upd = GetEnvValue(env, prefix, "PedestalUpdate", upd);
+        SetPedestalUpdate(GetEnvValue(env, prefix, "PedestalUpdate", fPedestalUpdate));
         rc = kTRUE;
     }
-    SetPedestalUpdate(upd);
 
     return rc;

trunk/MagicSoft/Mars/mpedestal/MPedCalcFromLoGain.h

@@ -62 +62 @@
 
     // Setter
-    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 SetRange(Byte_t hifirst=fgHiGainFirst, Byte_t hilast=fgHiGainLast, Byte_t lofirst=fgLoGainFirst, Byte_t lolast=fgLoGainLast);
+    void SetWindowSize(Byte_t windowh=fgHiGainWindowSize, Byte_t windowl=fgLoGainWindowSize);
+    void SetMaxHiGainVar(Byte_t maxvar=fgMaxHiGainVar);
     void SetDumpEvents(UInt_t dumpevents = 0) {fNumEventsDump = dumpevents;}
     void SetPedestalUpdate(Bool_t pedupdate)  {fPedestalUpdate = pedupdate;}

trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.cc

@@ -394 +394 @@
 
   const Int_t npixels  = fPedestals->GetSize();
-  const Int_t areas    = fPedestals->GetAverageAreas();
-  const Int_t sectors  = fPedestals->GetAverageSectors();
+//  const Int_t areas    = fPedestals->GetAverageAreas();
+//  const Int_t sectors  = fPedestals->GetAverageSectors();
   
   if (fSumx.GetSize()==0)
@@ -402 +402 @@
       fSumx2.Set(npixels);
       
-      fAreaSumx. Set(areas);
-      fAreaSumx2.Set(areas);
-      fAreaValid.Set(areas);
-      
-      fSectorSumx. Set(sectors);
-      fSectorSumx2.Set(sectors);
-      fSectorValid.Set(sectors);
+//      fAreaSumx. Set(areas);
+//      fAreaSumx2.Set(areas);
+//      fAreaValid.Set(areas);
+      
+//      fSectorSumx. Set(sectors);
+//      fSectorSumx2.Set(sectors);
+//      fSectorValid.Set(sectors);
       
       fSumx.Reset();
@@ -439 +439 @@
 Int_t MPedCalcPedRun::Process()
 {
+    MRawEvtPixelIter pixel(fRawEvt);
+
+    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;
+
+        if (fWindowSizeHiGain != 0)
+        {
+            do
+            {
+                sum += *ptr;
+                sqr += *ptr * *ptr;
+            }
+            while (++ptr != end);
+        }
+
+        if (fWindowSizeLoGain != 0)
+        {
+            ptr = pixel.GetLoGainSamples() + fLoGainFirst;
+            end = ptr + fWindowSizeLoGain;
+
+            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/fNumHiGainSlices;
+        //      const Float_t higainrms = TMath::Sqrt((msqr-msum*msum/fNumHiGainSlices)/(fNumHiGainSlices-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:
+        //
+        const Float_t sqrsum  = msum*msum;
+        fSumx2[idx]          += sqrsum;
+    }
+
+    fNumSamplesTot += fWindowSizeHiGain + fWindowSizeLoGain;
+
+    return kTRUE;
+/*
   MRawEvtPixelIter pixel(fRawEvt);
   
@@ -445 +507 @@
       const UInt_t idx    = pixel.GetPixelId();
       const UInt_t aidx   = (*fGeom)[idx].GetAidx();
-      const UInt_t sector = (*fGeom)[idx].GetSector();      
+      const UInt_t sector = (*fGeom)[idx].GetSector();
 
       Byte_t *ptr = pixel.GetHiGainSamples() + fHiGainFirst;
@@ -490 +552 @@
       fSumx[idx]          += msum;
       fAreaSumx[aidx]     += msum;
-      fSectorSumx[sector] += msum;      
+      fSectorSumx[sector] += msum;
       //
       // The old version:
@@ -502 +564 @@
       fSumx2[idx]          += sqrsum;
       fAreaSumx2[aidx]     += sqrsum;
-      fSectorSumx2[sector] += sqrsum;      
+      fSectorSumx2[sector] += sqrsum;
     }
   
   fPedestals->SetReadyToSave();
   fNumSamplesTot += fWindowSizeHiGain + fWindowSizeLoGain;
-
+*/
   return kTRUE;
 }
@@ -517 +579 @@
 Int_t MPedCalcPedRun::PostProcess()
 {
-  // Compute pedestals and rms from the whole run
-  const ULong_t n     = fNumSamplesTot;
-  const ULong_t nevts = GetNumExecutions();
-
-  MRawEvtPixelIter pixel(fRawEvt);
-  
-  while (pixel.Next())
-    {
-
-      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);
-      const Float_t higainped = sum/n;
-      //
-      // The old version:
-      //
-      //      const Float_t higainrms = TMath::Sqrt((sum2-sum*sum/n)/(n-1.));
-      //
-      // The new version:
-      //
-      // 1. Calculate the Variance of the sums:
-      Float_t higainVar = (sum2-sum*sum/nevts)/(nevts-1.);
-      // 2. Scale the variance to the number of slices:
-      higainVar /= (Float_t)(fWindowSizeHiGain+fWindowSizeLoGain);
-      // 3. Calculate the RMS from the Variance:
-      (*fPedestals)[pixid].Set(higainped, higainVar < 0 ? 0. : TMath::Sqrt(higainVar));
-
-    }
+    const ULong_t nevts = GetNumExecutions();
+    if (nevts<1)
+    {
+        *fLog << err << "ERROR - Not enough events recorded...abort." << endl;
+        return kFALSE;
+    }
+
+    *fLog << flush << inf << "Calculating pedestals..." << flush;
+
+    // Compute pedestals and rms from the whole run
+    const ULong_t n     = fNumSamplesTot;
+    const ULong_t npix  = fGeom->GetNumPixels();
+
+    for (UInt_t pixidx=0; pixidx<npix; pixidx++)
+    {
+        const Float_t sum  = fSumx.At(pixidx);
+        const Float_t sum2 = fSumx2.At(pixidx);
+        const Float_t higainped = sum/n;
+        //
+        // The old version:
+        //
+        //      const Float_t higainrms = TMath::Sqrt((sum2-sum*sum/n)/(n-1.));
+        //
+        // The new version:
+        //
+        // 1. Calculate the Variance of the sums:
+        Float_t higainVar = (sum2-sum*sum/nevts)/(nevts-1.);
+        // 2. Scale the variance to the number of slices:
+        higainVar /= (Float_t)(fWindowSizeHiGain+fWindowSizeLoGain);
+        // 3. Calculate the RMS from the Variance:
+        (*fPedestals)[pixidx].Set(higainped, higainVar < 0 ? 0. : TMath::Sqrt(higainVar), 0, nevts);
+    }
+
+    *fLog << flush << inf << "Calculating means..." << flush;
+
+    fPedestals->SetTotalEntries(fNumSamplesTot);
+    fPedestals->ReCalc(*fGeom);
+    fPedestals->SetReadyToSave();
+
+    return kTRUE;
 
   //
   // Loop over the (two) area indices to get the averaged pedestal per aidx
   //
+  /*
+
+   // Compute pedestals and rms from the whole run
+    const ULong_t n     = fNumSamplesTot;
+    const ULong_t nevts = GetNumExecutions();
+
+    MRawEvtPixelIter pixel(fRawEvt);
+
+    while (pixel.Next())
+    {
+        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);
+        const Float_t higainped = sum/n;
+        //
+        // The old version:
+        //
+        //      const Float_t higainrms = TMath::Sqrt((sum2-sum*sum/n)/(n-1.));
+        //
+        // The new version:
+        //
+        // 1. Calculate the Variance of the sums:
+        Float_t higainVar = (sum2-sum*sum/nevts)/(nevts-1.);
+        // 2. Scale the variance to the number of slices:
+        higainVar /= (Float_t)(fWindowSizeHiGain+fWindowSizeLoGain);
+        // 3. Calculate the RMS from the Variance:
+        (*fPedestals)[pixid].Set(higainped, higainVar < 0 ? 0. : TMath::Sqrt(higainVar), 0, nevts);
+    }
   for (Int_t aidx=0; aidx<fAreaValid.GetSize(); aidx++)
     {
@@ -582 +684 @@
       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++)
     {
@@ -613 +716 @@
       fPedestals->GetAverageSector(sector).Set(higainped, higainrms);
     }
-  
   fPedestals->SetTotalEntries(fNumSamplesTot);
   fPedestals->SetReadyToSave();
 
   return kTRUE;
+    */
 }
 

trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.h

@@ -33 +33 @@
     TArrayD fSumx;             // sum of values
     TArrayD fSumx2;            // sum of squared values
+    /*
     TArrayD fAreaSumx;         // averaged sum of values per area idx
     TArrayD fAreaSumx2;        // averaged sum of squared values per area idx
@@ -39 +40 @@
     TArrayD fSectorSumx2;      // averaged sum of squared values per sector
     TArrayI fSectorValid;      // number of valid pixel with sector idx
+    */
 
     Int_t  PreProcess (MParList *pList);

trunk/MagicSoft/Mars/mpedestal/MPedPhotCam.cc

@@ -229 +229 @@
 void MPedPhotCam::ReCalc(const MGeomCam &geom, MBadPixelsCam *bad)
 {
-    
     const Int_t np = GetSize();
     const Int_t ns = GetNumSectors();
     const Int_t na = GetNumAreas();
-
-    
 
     TArrayI acnt(na);
@@ -243 +240 @@
     TArrayD ssumr(ns);
 
-
     for (int i=0; i<np; i++)
     {
-        
-
-      if (bad && (*bad)[i].IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
-        continue; //was: .IsBad()       
+        if (bad && (*bad)[i].IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
+            continue; //was: .IsBad()
 
         // Create sums for areas and sectors
@@ -268 +262 @@
         ssumr[sect] += ne*rms;
         scnt[sect]  += ne;
-    
-
-    }
-
-    for (int i=0; i<ns; i++){
-      if (scnt[i]>0)  GetSector(i).Set(ssumx[i]/scnt[i], ssumr[i]/scnt[i], scnt[i]);
-      else GetSector(i).Set(-1., -1., 0);
-    }
-
-    for (int i=0; i<na; i++){
-      if (acnt[i]>0) GetArea(i).Set(asumx[i]/acnt[i], asumr[i]/acnt[i], acnt[i]);
-      else  GetArea(i).Set(-1., -1., 0);
-    }
+    }
+
+    for (int i=0; i<ns; i++)
+        if (scnt[i]>0)
+            GetSector(i).Set(ssumx[i]/scnt[i], ssumr[i]/scnt[i], scnt[i]);
+        else
+            GetSector(i).Clear();
+
+    for (int i=0; i<na; i++)
+        if (acnt[i]>0)
+            GetArea(i).Set(asumx[i]/acnt[i], asumr[i]/acnt[i], acnt[i]);
+        else
+            GetArea(i).Clear();
 }
 

trunk/MagicSoft/Mars/mpedestal/MPedPhotCam.h

@@ -51 +51 @@
     void Print(Option_t *o="") const;
 
-    void ReCalc(const MGeomCam &geom, MBadPixelsCam *bad);
+    void ReCalc(const MGeomCam &geom, MBadPixelsCam *bad=NULL);
 
     Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const;

trunk/MagicSoft/Mars/mpedestal/MPedestalCam.cc

@@ -35 +35 @@
 #include "MPedestalPix.h"
 
+#include <TArrayI.h>
+#include <TArrayD.h>
+
 #include <TClonesArray.h>
 
@@ -43 +46 @@
 
 #include "MGeomCam.h"
+#include "MGeomPix.h"
+
+#include "MBadPixelsCam.h"
+#include "MBadPixelsPix.h"
 
 ClassImp(MPedestalCam);
@@ -311 +318 @@
 }
 
-Bool_t MPedestalCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const {
-
-  if (GetSize() <= idx)
-    return kFALSE;
-
-  if (!(*this)[idx].IsValid())
-    return kFALSE;
-
-  switch (type) {
+// --------------------------------------------------------------------------
+//
+// Calculates the avarage pedestal and pedestal rms for all sectors
+// and pixel sizes. The geometry container is used to get the necessary
+// geometry information (sector number, size index) If the bad pixel
+// container is given all pixels which have the flag 'bad' are ignored
+// in the calculation of the sector and size average.
+//
+void MPedestalCam::ReCalc(const MGeomCam &geom, MBadPixelsCam *bad)
+{
+    const Int_t np = GetSize();
+    const Int_t ns = geom.GetNumSectors();
+    const Int_t na = geom.GetNumAreas();
+
+    TArrayI acnt(na);
+    TArrayI scnt(ns);
+    TArrayD asumx(na);
+    TArrayD ssumx(ns);
+    TArrayD asumr(na);
+    TArrayD ssumr(ns);
+
+    for (int i=0; i<np; i++)
+    {
+        if (bad && (*bad)[i].IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
+            continue; //was: .IsBad()
+
+        // Create sums for areas and sectors
+        const UInt_t aidx = geom[i].GetAidx();
+        const UInt_t sect = geom[i].GetSector();
+
+        const MPedestalPix &pix = (*this)[i];
+
+        const UInt_t  ne   = pix.GetNumEvents();
+        const Float_t mean = pix.GetPedestal();
+        const Float_t rms  = pix.GetPedestalRms();
+
+        asumx[aidx] += ne*mean;
+        asumr[aidx] += ne*rms;
+        acnt[aidx]  += ne;
+
+        ssumx[sect] += ne*mean;
+        ssumr[sect] += ne*rms;
+        scnt[sect]  += ne;
+    }
+
+    for (int i=0; i<ns; i++)
+        if (scnt[i]>0)
+            GetAverageSector(i).Set(ssumx[i]/scnt[i], ssumr[i]/scnt[i], 0, scnt[i]);
+        else
+            GetAverageSector(i).Clear();
+
+    for (int i=0; i<na; i++)
+        if (acnt[i]>0)
+            GetAverageArea(i).Set(asumx[i]/acnt[i], asumr[i]/acnt[i], 0, acnt[i]);
+        else
+            GetAverageArea(i).Clear();
+}
+
+Bool_t MPedestalCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
+{
+    if (GetSize() <= idx)
+        return kFALSE;
+
+    if (!(*this)[idx].IsValid())
+        return kFALSE;
+
+    switch (type)
+    {
     case 0:
-      val = (*this)[idx].GetPedestal();
-      break;
+        val = (*this)[idx].GetPedestal();
+        break;
     case 1:
-      val = fTotalEntries > 0 ? 
-          (*this)[idx].GetPedestalRms()/TMath::Sqrt((Float_t)fTotalEntries) 
-        : (*this)[idx].GetPedestalError();
-      break;
+        val = fTotalEntries > 0 ?
+            (*this)[idx].GetPedestalRms()/TMath::Sqrt((Float_t)fTotalEntries)
+            : (*this)[idx].GetPedestalError();
+        break;
     case 2:
-      val = (*this)[idx].GetPedestalRms();
-      break;
+        val = (*this)[idx].GetPedestalRms();
+        break;
     case 3:
-      val = fTotalEntries > 0 ?
-          (*this)[idx].GetPedestalRms()/TMath::Sqrt((Float_t)fTotalEntries)/2. 
-        : (*this)[idx].GetPedestalRmsError();
-      break;
+        val = fTotalEntries > 0 ?
+            (*this)[idx].GetPedestalRms()/TMath::Sqrt((Float_t)fTotalEntries)/2.
+            : (*this)[idx].GetPedestalRmsError();
+        break;
     default:
-      return kFALSE;
+        return kFALSE;
     }
-  return kTRUE;
+    return kTRUE;
 }
 
 void MPedestalCam::DrawPixelContent(Int_t idx) const
 {
-  *fLog << warn << "MPedestalCam::DrawPixelContent - not available." << endl;
-}
+    *fLog << warn << "MPedestalCam::DrawPixelContent - not available." << endl;
+}

trunk/MagicSoft/Mars/mpedestal/MPedestalCam.h

@@ -13 +13 @@
 class MGeomCam;
 class MPedestalPix;
+class MBadPixelsCam;
 
 class MPedestalCam : public MParContainer, public MCamEvent
@@ -51 +52 @@
   void  InitAverageSectors             ( const UInt_t i      );
 
+  void ReCalc(const MGeomCam &geom, MBadPixelsCam *bad=NULL);
+
   void Print(Option_t *o="") const;
   

trunk/MagicSoft/Mars/mpedestal/MPedestalPix.h

@@ -37 +37 @@
     Float_t GetPedestalRmsError() const { return fNumEvents>0 ? fPedestalRms/TMath::Sqrt((Float_t)fNumEvents/2) : 0; }
 
+    UInt_t  GetNumEvents() const { return fNumEvents; }
+
     Bool_t IsValid() const;