Changeset 5487 for trunk/MagicSoft/Mars/mpedestal

Timestamp:

11/27/04 16:36:36 (20 years ago)

Author:

gaug

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mpedestal

Files:

• MPedCalcPedRun.cc (modified) (15 diffs)
• MPedCalcPedRun.h (modified) (2 diffs)

trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.cc

@@ -109 +109 @@
 //   MRawEvtData
 //   MRawRunHeader
+//   MRawEvtHeader
 //   MGeomCam
 //
@@ -118 +119 @@
 /////////////////////////////////////////////////////////////////////////////
 #include "MPedCalcPedRun.h"
-#include "MExtractor.h"
+#include "MExtractPedestal.h"
+
+#include "MExtractTimeAndCharge.h"
 
 #include "MParList.h"
@@ -126 +129 @@
 
 #include "MRawRunHeader.h"  
+#include "MRawEvtHeader.h"  
 #include "MRawEvtPixelIter.h"
 #include "MRawEvtData.h"
@@ -135 +139 @@
 #include "MGeomCam.h"
 
+#include "MStatusDisplay.h"
+
 ClassImp(MPedCalcPedRun);
 
 using namespace std;
 
-const Byte_t MPedCalcPedRun::fgHiGainFirst      = 0;
-const Byte_t MPedCalcPedRun::fgHiGainLast       = 29;
-const Byte_t MPedCalcPedRun::fgLoGainFirst      = 0;
-const Byte_t MPedCalcPedRun::fgLoGainLast       = 14;
-const Byte_t MPedCalcPedRun::fgHiGainWindowSize = 14;
-const Byte_t MPedCalcPedRun::fgLoGainWindowSize = 0;
+const UShort_t MPedCalcPedRun::fgExtractWinFirst       = 0;
+const UShort_t MPedCalcPedRun::fgExtractWinSize        = 6;
+const UInt_t   MPedCalcPedRun::gkFirstRunWithFinalBits = 45589;
 // --------------------------------------------------------------------------
 //
@@ -161 +164 @@
 //
 MPedCalcPedRun::MPedCalcPedRun(const char *name, const char *title)
-    : fWindowSizeHiGain(fgHiGainWindowSize), 
-      fWindowSizeLoGain(fgLoGainWindowSize), 
-      fGeom(NULL)
-{
-    fName  = name  ? name  : "MPedCalcPedRun";
-    fTitle = title ? title : "Task to calculate pedestals from pedestal runs raw data";
-
-    AddToBranchList("fHiGainPixId");
-    AddToBranchList("fHiGainFadcSamples");
-
-    SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
-    Clear();
+{
+  fName  = name  ? name  : "MPedCalcPedRun";
+  fTitle = title ? title : "Task to calculate pedestals from pedestal runs raw data";
+
+  AddToBranchList("fHiGainPixId");
+  AddToBranchList("fHiGainFadcSamples");
+  
+  SetExtractWindow(fgExtractWinFirst, fgExtractWinSize);
+
+  Clear();
+  
+  fFirstRun   = kTRUE;
+  fSkip       = kFALSE;
+  fOverlap    = 0;
+  fUsedEvents = 0;
+}
+
+
+
+void MPedCalcPedRun::ResetArrays()
+{
+
+  MExtractPedestal::ResetArrays();
+  fUsedEvents    = 0;
 }
 
 // --------------------------------------------------------------------------
 //
-// Sets:
-// - fNumSamplesTot to 0
-// - fRawEvt to NULL
-// - fRunHeader to NULL
-// - fPedestals to NULL
-//
-void MPedCalcPedRun::Clear(const Option_t *o)
-{
-    fNumSamplesTot = 0;
-
-    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();
-
-        fAreaSumx. Reset();
-        fAreaSumx2.Reset();
-        fAreaSumAB0.Reset();
-        fAreaSumAB1.Reset();
-        fAreaValid.Reset();
-        
-        fSectorSumx. Reset();
-        fSectorSumx2.Reset();
-        fSectorSumAB0.Reset();
-        fSectorSumAB1.Reset();
-        fSectorValid.Reset();
-    }
-}
-
-// --------------------------------------------------------------------------
-//
-// SetRange: 
-//
-// Calls:
-// - MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
-// - SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain);
-//
-void MPedCalcPedRun::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);
-}
-
-// --------------------------------------------------------------------------
-//
-// 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 MPedCalcPedRun::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;
-    
-  if (fWindowSizeHiGain == 0)
-    {
-      *fLog << warn;
-      *fLog << GetDescriptor() << ": HiGain window currently set to 0, will set it to 2 samples " << endl;
-      fWindowSizeHiGain = 2;
-    }
-  
-  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
-// 
-// The following output containers are also searched and created if
-// they were not found:
-//
-//  - MPedestalCam
-//
-Int_t MPedCalcPedRun::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;
-
-  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.
+// In case that the variables fExtractLast is greater than the number of 
+// High-Gain FADC samples obtained from the run header, the flag 
+// fOverlap is set to the difference and fExtractLast is set back by the
+// same number of slices.
+//
+// The run type is checked for "Pedestal" (run type: 1) 
+// and the variable fSkip is set in that case
 //
 Bool_t MPedCalcPedRun::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;
+  MExtractPedestal::ReInit(pList);
+
+  const UShort_t lastavailable = fRunHeader->GetNumSamplesHiGain()-1;
+  
+  if (fExtractWinLast > lastavailable)
+    {
+      const UShort_t diff = fExtractWinLast - 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;
+            << Form("%s%2i%s%2i%s%2i%s",": Selected ExtractWindow [",
+                    fExtractWinFirst,",",fExtractWinLast,
+                    "] ranges out of 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")
+            << Form("%s%2i%s",": Will use ",diff," samples from the 'Low-Gain' slices for the pedestal extraction")
             << 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)
+      fExtractWinLast -= diff;
+      fOverlap         = diff;
+    }
+
+  const UShort_t runtype = fRunHeader->GetRunType();
+
+  switch (runtype)
+    {
+    case 1: 
+      fFirstRun = kFALSE;
+      fSkip     = kFALSE;
+      return kTRUE;
+      break;
+    default:
+      fSkip     = kTRUE;
+      if (!fFirstRun)
         {
-          fWindowSizeHiGain -= diff;
-          fWindowSizeLoGain += diff;
+          *fLog << endl;
+          *fLog << inf << GetDescriptor() << " : Finalize pedestal calculations..." << flush;
+          CallPostProcess();
+          Reset();
         }
-      else
-        {
-          fWindowSizeLoGain += fWindowSizeHiGain;
-          fLoGainFirst       = diff-fWindowSizeHiGain;
-          fWindowSizeHiGain  = 0;
-        }
-    }
-
-
-  const Int_t npixels  = fPedestals->GetSize();
-  const Int_t areas    = fPedestals->GetAverageAreas();
-  const Int_t sectors  = fPedestals->GetAverageSectors();
-  
-  if (fSumx.GetSize()==0)
-    {
-      fSumx.  Set(npixels);
-      fSumx2. Set(npixels);
-      fSumAB0.Set(npixels);
-      fSumAB1.Set(npixels);
-      
-      fAreaSumx.  Set(areas);
-      fAreaSumx2. Set(areas);
-      fAreaSumAB0.Set(areas);
-      fAreaSumAB1.Set(areas);
-      fAreaValid. Set(areas);
-      
-      fSectorSumx.  Set(sectors);
-      fSectorSumx2. Set(sectors);
-      fSectorSumAB0.Set(sectors);
-      fSectorSumAB1.Set(sectors);
-      fSectorValid. Set(sectors);
-      
-      fSumx. Reset();
-      fSumx2.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;        
+      break;
+    }
+  
+  Print();
 
   return kTRUE;
@@ -468 +250 @@
 
 // --------------------------------------------------------------------------
+//
+// Return kTRUE (without doing anything) in case that the run type is not 
+// equal to 1 (pedestal run)
 //
 // Fill the MPedestalCam container with the signal mean and rms for the event.
@@ -476 +261 @@
 {
 
+  if (fSkip && !IsPedBitSet())
+    return kTRUE;
+
+  /*
+  *fLog << err << dec << fEvtHeader->GetTriggerID() << endl;
+  for (Int_t i=16; i>= 0; i--)
+    *fLog << inf << (fEvtHeader->GetTriggerID() >> i & 1);
+  *fLog  << endl;
+  *fLog << warn << hex << fEvtHeader->GetTriggerID() << endl;
+  */
+
+  fUsedEvents++;
+
   MRawEvtPixelIter pixel(fRawEvt);
   
@@ -484 +282 @@
       const UInt_t sector = (*fGeom)[idx].GetSector();      
 
-      Byte_t *ptr = pixel.GetHiGainSamples() + fHiGainFirst;
-      Byte_t *end = ptr + fWindowSizeHiGain;
-      
-      UInt_t sum = 0;
-      UInt_t sqr = 0;
-      UInt_t ab0 = 0;
-      UInt_t ab1 = 0;
-      Int_t  cnt = 0;
-
-      if (fWindowSizeHiGain != 0)
-        {
-          do
-            {
-              sum += *ptr;
-              sqr += *ptr * *ptr;
-
-              if (pixel.IsABFlagValid())
-                {
-                  const Int_t abFlag = (fHiGainFirst + pixel.HasABFlag() + cnt) & 0x1;
-                  if (abFlag)
-                    ab1 += *ptr;
-                  else
-                    ab0 += *ptr;
-
-                  cnt++;
-                }
-            }
-          while (++ptr != end);
-        }
-
-      cnt = 0;
-      
-      if (fWindowSizeLoGain != 0)
-        {
-          
-          ptr = pixel.GetLoGainSamples() + fLoGainFirst;
-          end = ptr + fWindowSizeLoGain;
-      
-          do
-            {
-              sum += *ptr;
-              sqr += *ptr * *ptr;
-
-              if (pixel.IsABFlagValid())
-                {
-                  const Int_t abFlag = (fLoGainFirst + pixel.GetNumHiGainSamples() + pixel.HasABFlag() + cnt) 
-                                      & 0x1;
-                  if (abFlag)
-                    ab1 += *ptr;
-                  else
-                    ab0 += *ptr;
-
-                  cnt++;
-                }
-
-            }
-          while (++ptr != end);
-          
-        }
-      
+      Float_t sum = 0.;
+      UInt_t  ab0 = 0;
+      UInt_t  ab1 = 0;
+      
+      if (fExtractor)
+        CalcExtractor( &pixel, sum, (*fPedestalsIn)[idx]);
+      else
+        CalcSums( &pixel, sum, ab0, ab1);
+
       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;
@@ -558 +297 @@
       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;
@@ -571 +302 @@
       fSectorSumx2[sector] += sqrsum;      
 
-      //
-      // Now, the sums separated for AB0 and AB1
-      //
       fSumAB0[idx]        += ab0;
       fSumAB1[idx]        += ab1;
@@ -584 +312 @@
     }
   
-  fPedestals->SetReadyToSave();
-  fNumSamplesTot += fWindowSizeHiGain + fWindowSizeLoGain;
+  fPedestalsOut->SetReadyToSave();
 
   return kTRUE;
 }
 
+
+
+void MPedCalcPedRun::CalcExtractor( MRawEvtPixelIter *pixel, Float_t &sum, MPedestalPix &ped)
+{
+  
+  Byte_t  sat = 0;
+  Byte_t *first  = pixel->GetHiGainSamples() + fExtractWinFirst;
+  Byte_t *logain = pixel->GetLoGainSamples();
+  const Bool_t abflag  = pixel->HasABFlag();
+  Float_t dummy;
+  fExtractor->FindTimeAndChargeHiGain(first,logain,sum,dummy,dummy,dummy,sat,ped,abflag);
+}
+
+
+void MPedCalcPedRun::CalcSums( MRawEvtPixelIter *pixel, Float_t &sum, UInt_t &ab0, UInt_t &ab1)
+{
+
+  Byte_t *higain = pixel->GetHiGainSamples() + fExtractWinFirst;
+  Byte_t *ptr = higain;
+  Byte_t *end = ptr + fExtractWinSize;
+
+  const Bool_t abflag = pixel->HasABFlag();
+
+  Int_t sumi = 0;
+      
+  Int_t  cnt = 0;
+  do
+    {
+      sumi += *ptr;
+      if (pixel->IsABFlagValid())
+        {
+          const Int_t abFlag = (fExtractWinFirst + abflag + cnt) & 0x1;
+          if (abFlag)
+            ab1 += *ptr;
+          else
+            ab0 += *ptr;
+          
+          cnt++;
+        }
+    }
+  while (++ptr != end);
+  
+  if (fOverlap != 0)
+    {
+      ptr = pixel->GetLoGainSamples();
+      end = ptr + fOverlap;
+      
+      do
+        {
+          sumi += *ptr;
+          if (pixel->IsABFlagValid())
+            {
+              const Int_t abFlag = (fExtractWinFirst + abflag + cnt) & 0x1;
+              if (abFlag)
+                ab1 += *ptr;
+              else
+                ab0 += *ptr;
+              
+              cnt++;
+            }
+        }
+      while (++ptr != end);
+    }
+
+  sum = (Float_t)sumi;
+
+}
+
 // --------------------------------------------------------------------------
 //
@@ -597 +392 @@
 {
 
-  // Compute pedestals and rms from the whole run
-  const ULong_t n     = fNumSamplesTot;
-  const ULong_t nevts = GetNumExecutions();
+  if (fUsedEvents == 0)
+    return kTRUE;
 
   MRawEvtPixelIter pixel(fRawEvt);
@@ -605 +399 @@
   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:
-      const Float_t rms = higainVar<0 ? 0. : TMath::Sqrt(higainVar);
-      // 4. Calculate the amplitude of the 150MHz "AB" noise
-      const Float_t abOffs = (fSumAB0.At(pixid) - fSumAB1.At(pixid)) / n;
-      
-      (*fPedestals)[pixid].Set(higainped,rms,abOffs);
+      CalcPixResults(fUsedEvents,pixid);
     }
 
@@ -638 +406 @@
   // Loop over the (two) area indices to get the averaged pedestal per aidx
   //
-  for (Int_t aidx=0; aidx<fAreaValid.GetSize(); aidx++)
-    {
-      
+  for (UInt_t aidx=0; aidx<fAreaValid.GetSize(); aidx++)
+    {
       const Int_t   napix = fAreaValid.At(aidx);
-
       if (napix == 0)
         continue;
 
-      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 /= fWindowSizeHiGain+fWindowSizeLoGain;
-      // 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);
-      // 5. Calculate the amplitude of the 150MHz "AB" noise
-      const Float_t abOffs = (fAreaSumAB0.At(aidx) - fAreaSumAB1.At(aidx)) / an;
-
-      fPedestals->GetAverageArea(aidx).Set(higainped, higainrms,abOffs);
+      CalcAreaResults(fUsedEvents,napix,aidx);
     }
   
@@ -671 +418 @@
   // Loop over the (six) sector indices to get the averaged pedestal per sector
   //
-  for (Int_t sector=0; sector<fSectorValid.GetSize(); sector++)
-    {
-      
+  for (UInt_t sector=0; sector<fSectorValid.GetSize(); sector++)
+    {
       const Int_t   nspix = fSectorValid.At(sector);
-
       if (nspix == 0)
         continue;
       
-      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 /= fWindowSizeHiGain+fWindowSizeLoGain;
-      // 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);
-      // 5. Calculate the amplitude of the 150MHz "AB" noise
-      const Float_t abOffs = (fSectorSumAB0.At(sector) - fSectorSumAB1.At(sector)) / sn;
-
-      fPedestals->GetAverageSector(sector).Set(higainped, higainrms, abOffs);
-    }
-  
-  fPedestals->SetTotalEntries(fNumSamplesTot);
-  fPedestals->SetReadyToSave();
+      CalcSectorResults(fUsedEvents,nspix,sector);
+    }
+  
+  fPedestalsOut->SetTotalEntries(fUsedEvents*fExtractWinSize);
+  fPedestalsOut->SetReadyToSave();
 
   return kTRUE;
 }
 
-Int_t MPedCalcPedRun::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
-{
-    if (MExtractor::ReadEnv(env, prefix, print)==kERROR)
-        return kERROR;
-
-    Byte_t hw = fWindowSizeHiGain;
-    Byte_t lw = fWindowSizeLoGain;
-
-    Bool_t rc = kFALSE;
-
-    if (IsEnvDefined(env, prefix, "WindowSizeHiGain", print))
-    {
-        hw = GetEnvValue(env, prefix, "WindowSizeHiGain", hw);
-        rc=kTRUE;
-    }
-
-    if (IsEnvDefined(env, prefix, "WindowSizeLoGain", print))
-    {
-        lw = GetEnvValue(env, prefix, "WindowSizeLoGain", lw);
-        rc=kTRUE;
-    }
-
-    if (rc)
-        SetWindowSize(hw, lw);
-
-    return rc;
-}
+// -----------------------------------------------------------------------
+//
+// Analogue to the MTask::CallPostProcess, but without the manipulation 
+// of the bit fIsPreProcessed. Needed in order to call PostProcess multiple 
+// times in the intensity calibration.
+//
+Int_t MPedCalcPedRun::CallPostProcess()
+{
+
+  *fLog << all << fName << "... " << flush;
+  if (fDisplay)
+    fDisplay->SetStatusLine2(*this);
+  
+  return PostProcess();
+}
+
+//-------------------------------------------------------------
+// 
+// Return if the pedestal bit was set from the calibration trigger box.
+// The last but one bit is used for the "pedestal-bit".
+//
+// This bit is set since run gkFinalizationTriggerBit
+//
+Bool_t MPedCalcPedRun::IsPedBitSet()
+{
+  if (fRunHeader->GetRunNumber() < gkFirstRunWithFinalBits)
+    return kFALSE;
+      
+  return (fEvtHeader->GetTriggerID() >> 3 & 1);
+}
+
+void MPedCalcPedRun::Print(Option_t *o) const
+{
+
+  MExtractPedestal::Print(o);
+
+  *fLog << "Number overlap low-gain slices:           " << fOverlap << endl;
+  *fLog << "First run out of sequence:                " << (fFirstRun?"yes":"no") << endl;
+  *fLog << "Skip this run:                            " << (fSkip?"yes":"no") << endl;
+  *fLog << "Number of used events so far:             " << fUsedEvents << endl;
+  *fLog << endl;
+}

trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.h

@@ -2 +2 @@
 #define MARS_MPedCalcPedRun
 
-#ifndef MARS_MExtractor
-#include "MExtractor.h"
+#ifndef MARS_MExtractPedestal
+#include "MExtractPedestal.h"
 #endif
 
-#ifndef ROOT_TArrayD
-#include <TArrayD.h>
+#ifndef MARS_MArrayD
+#include <MArrayD.h>
 #endif
 
-#ifndef ROOT_TArrayI
-#include <TArrayI.h>
+#ifndef MARS_MArrayI
+#include <MArrayI.h>
 #endif
 
-class MGeomCam;
-class MPedCalcPedRun : public MExtractor
+class MRawEvtPixelIter;
+class MPedestalPix;
+class MPedCalcPedRun : public MExtractPedestal
 {
 
-  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 UShort_t fgExtractWinFirst;  // First FADC slice Hi-Gain (currently set to: 3)
+  static const UShort_t fgExtractWinSize;   // Extraction Size Hi-Gain (currently set to: 14)
+  static const UInt_t   gkFirstRunWithFinalBits; // First Run with pedestal trigger bit at place 3 
   
-  UInt_t  fNumSamplesTot;
-  Byte_t  fWindowSizeHiGain; // Number of Hi Gain slices in window
-  Byte_t  fWindowSizeLoGain; // Number of Lo Gain slices in window
+  UShort_t fOverlap;         // Number of overlapping slices from High-Gain to Low-Gain
+
+  Bool_t  fFirstRun;         // Flag to tell if the first run out of many is used
+  Bool_t  fSkip;             // Flag to tell if the Process has to be skipped 
+  ULong_t fUsedEvents;       // Number of used (not skipped) events
   
-  MGeomCam *fGeom;           // Camera geometry
+  Bool_t IsPedBitSet();
   
-  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
-  TArrayD fAreaSumx;         // averaged sum of values per area idx
-  TArrayD fAreaSumx2;        // averaged sum of squared values per area idx
-  TArrayD fAreaSumAB0;       // averaged sum of ABFlag=0 slices per area idx
-  TArrayD fAreaSumAB1;       // averaged sum of ABFlag=1 slices per area idx
-  TArrayI fAreaValid;        // number of valid pixel with area idx
-  TArrayD fSectorSumx;       // averaged sum of values per sector
-  TArrayD fSectorSumx2;      // averaged sum of squared values per sector
-  TArrayD fSectorSumAB0;     // averaged sum of ABFlag=0 slices per sector
-  TArrayD fSectorSumAB1;     // averaged sum of ABFlag=1 slices 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();
-  Int_t  ReadEnv(const TEnv &env, TString prefix, Bool_t print);
+  
+  void ResetArrays();
+  void CalcSums   ( MRawEvtPixelIter *pixel, Float_t &sum, UInt_t &ab0, UInt_t &ab1);
+  void CalcExtractor( MRawEvtPixelIter *pixel, Float_t &sum, MPedestalPix &ped);
   
 public:
@@ -56 +43 @@
   MPedCalcPedRun(const char *name=NULL, const char *title=NULL);
   
-  void Clear(const Option_t *o="");
+  Int_t CallPostProcess();  
+
+  void Print(Option_t *o="") const;
   
-  void SetRange        ( const Byte_t hifirst=0, const Byte_t hilast=0, 
-                         const Byte_t lofirst=0, const Byte_t lolast=0  );
-  void SetWindowSize   ( const Byte_t windowh=0, const Byte_t windowl=0 );
-  
-  ClassDef(MPedCalcPedRun, 0)   // Task to calculate pedestals from pedestal runs raw data
+  ClassDef(MPedCalcPedRun, 1)   // Task to calculate pedestals from pedestal runs
 };