Changeset 3889 for trunk/MagicSoft/Mars/mpedestal

Timestamp:

04/29/04 15:48:38 (21 years ago)

Author:

gaug

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.cc (modified) (17 diffs)

trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.cc

@@ -26 +26 @@
 !
 \* ======================================================================== */
-
 /////////////////////////////////////////////////////////////////////////////
 //
@@ -70 +69 @@
 // (and moreover asymmetric) and that they are also correlated.)
 // 
+//  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  =  14
+//   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:
@@ -80 +99 @@
 /////////////////////////////////////////////////////////////////////////////
 #include "MPedCalcPedRun.h"
+#include "MExtractor.h"
 
 #include "MParList.h"
@@ -93 +113 @@
 #include "MPedestalCam.h"
 
-#include "MExtractedSignalPix.h"
-#include "MExtractedSignalCam.h"
-
 #include "MGeomPix.h"
 #include "MGeomCam.h"
@@ -105 +122 @@
 using namespace std;
 
-// --------------------------------------------------------------------------
-//
-// default constructor
+const Byte_t MPedCalcPedRun::fgHiGainFirst      = 0;
+const Byte_t MPedCalcPedRun::fgHiGainLast       = 14;
+const Byte_t MPedCalcPedRun::fgLoGainFirst      = 0;
+const Byte_t MPedCalcPedRun::fgLoGainLast       = 14;
+const Byte_t MPedCalcPedRun::fgHiGainWindowSize = 14;
+const Byte_t MPedCalcPedRun::fgLoGainWindowSize = 0;
+// --------------------------------------------------------------------------
+//
+// Default constructor: 
+//
+// Sets:
+// - fWindowSizeHiGain to fgHiGainWindowSize
+// - fWindowSizeLoGain to fgLoGainWindowSize
+//
+// Calls: 
+// - AddToBranchList("fHiGainPixId");
+// - AddToBranchList("fHiGainFadcSamples");
+// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast)
+// - Clear()
 //
 MPedCalcPedRun::MPedCalcPedRun(const char *name, const char *title)
-    : fRawEvt(NULL), fPedestals(NULL)
-{
-    fName  = name  ? name  : "MPedCalcPedRun";
-    fTitle = title ? title : "Task to calculate pedestals from pedestal runs raw data";
-
-    AddToBranchList("fHiGainPixId");
-    AddToBranchList("fHiGainFadcSamples");
-
-    fNumHiGainSamples = 0;
-    Clear();
-}
-
+    : fWindowSizeHiGain(fgHiGainWindowSize), 
+      fWindowSizeLoGain(fgLoGainWindowSize)
+{
+  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();
+}
+
+// --------------------------------------------------------------------------
+//
+// Sets:
+// - fNumSamplesTot to 0
+// - fRawEvt to NULL
+// - fRunHeader to NULL
+// - fPedestals to NULL
+//
 void MPedCalcPedRun::Clear(const Option_t *o)
 {
@@ -128 +170 @@
 
   fRawEvt    = NULL;
+  fRunHeader = NULL;
   fPedestals = NULL;
 }
 
-
+// --------------------------------------------------------------------------
+//
+// 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;
+    
+  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);
+
+}
+
+  
+    
 // --------------------------------------------------------------------------
 //
@@ -137 +261 @@
 //
 //  - MRawEvtData
+//  - MRawRunHeader
+//  - MGeomCam
 // 
 // The following output containers are also searched and created if
@@ -155 +281 @@
     }
   
+  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)
@@ -165 +298 @@
   if (!fPedestals)
     return kFALSE;
-  
+
   return kTRUE;
 }
@@ -171 +304 @@
 // --------------------------------------------------------------------------
 //
-// The ReInit searches for the following input containers:
-//  - MRawRunHeader
-//
-// It also initializes the data arrays fSumx and fSumx2 
-// (only for the first read file)
-// 
+// 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.
+//
 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->IsMonteCarloRun())
-        return kTRUE;
-
-    Int_t npixels  = fPedestals->GetSize();
-    Int_t areas    = fPedestals->GetAverageAreas();
-    Int_t sectors  = fPedestals->GetAverageSectors();
-
-    if (fSumx.GetSize()==0)
-    {
-        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();
-        fSumx2.Reset();
-    }
-
-    // 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 lastdesired   = (Int_t)fHiGainLast;
+  Int_t 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 FADC Window [",
+                    (int)fHiGainFirst,",",lastdesired,
+                    "] ranges out of the available limits: [0,",lastavailable,"].") << endl;
+      *fLog << GetDescriptor() << ": Will reduce the upper edge to " << (int)(fHiGainLast - diff) << endl;
+      SetRange(fHiGainFirst, fHiGainLast-diff, fLoGainFirst, fLoGainLast);
+    }
+  
+  lastdesired   = (Int_t)(fLoGainLast);
+  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);
+    }
+
+  Int_t npixels  = fPedestals->GetSize();
+  Int_t areas    = fPedestals->GetAverageAreas();
+  Int_t sectors  = fPedestals->GetAverageSectors();
+  
+  if (fSumx.GetSize()==0)
+    {
+      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();
+      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;
+  
 }
 // --------------------------------------------------------------------------
@@ -235 +407 @@
       const UInt_t sector = (*fGeom)[idx].GetSector();      
 
-      Byte_t *ptr = pixel.GetHiGainSamples();
-      const Byte_t *end = ptr + fNumHiGainSamples;
+      Byte_t *ptr = pixel.GetHiGainSamples() + fHiGainFirst;
+      Byte_t *end = ptr + fWindowSizeHiGain;
       
       UInt_t sum = 0;
@@ -247 +419 @@
         }
       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;
@@ -254 +441 @@
       // 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.));
+      //        const Float_t higainped = msum/fNumHiGainSlices;
+      //        const Float_t higainrms = TMath::Sqrt((msqr-msum*msum/fNumHiGainSlices)/(fNumHiGainSlices-1.));
       //        (*fPedestals)[idx].Set(higainped, higainrms);
       
@@ -276 +463 @@
   
   fPedestals->SetReadyToSave();
-  fNumSamplesTot += fNumHiGainSamples;
-  
+  fNumSamplesTot += fWindowSizeHiGain + fWindowSizeLoGain;
   
   return kTRUE;
@@ -319 +505 @@
       Float_t higainVar = (sum2-sum*sum/nevts)/(nevts-1.);
       // 2. Scale the variance to the number of slices:
-      higainVar /= (Float_t)fNumHiGainSamples;
+      higainVar /= (Float_t)(fWindowSizeHiGain+fWindowSizeLoGain);
       // 3. Calculate the RMS from the Variance:
       const Float_t higainrms = TMath::Sqrt(higainVar);
@@ -344 +530 @@
       Float_t higainVar = (sum2-sum*sum/aevts)/(aevts-1.);
       // 2. Scale the variance to the number of slices:
-      higainVar /= (Float_t)fNumHiGainSamples;
+      higainVar /= fWindowSizeHiGain+fWindowSizeLoGain;
       // 3. Calculate the RMS from the Variance:
       Float_t higainrms = TMath::Sqrt(higainVar);
@@ -371 +557 @@
       Float_t higainVar = (sum2-sum*sum/sevts)/(sevts-1.);
       // 2. Scale the variance to the number of slices:
-      higainVar /= (Float_t)fNumHiGainSamples;
+      higainVar /= fWindowSizeHiGain+fWindowSizeLoGain;
       // 3. Calculate the RMS from the Variance:
       Float_t higainrms = TMath::Sqrt(higainVar);