Changeset 2629 for trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.cc

Timestamp:

12/10/03 08:24:42 (21 years ago)

Author:

gaug

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.cc (modified) (8 diffs)

trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.cc

@@ -50 +50 @@
 #include "MPedestalCam.h"
 
+#include "MGeomCamMagic.h"
+
 ClassImp(MPedCalcPedRun);
 
@@ -61 +63 @@
     AddToBranchList("fHiGainPixId");
     AddToBranchList("fHiGainFadcSamples");
+    
+    fCounter=0;
 }
 
@@ -75 +79 @@
     if (!fPedestals)
         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);
+    }
 
     return kTRUE;
@@ -110 +125 @@
     while (pixel.Next())
     {
-             Byte_t *ptr = pixel.GetHiGainSamples();
-        const Byte_t *end = ptr + fRawEvt->GetNumHiGainSamples();
+        Byte_t shift=(fNumHiGainSamples/2*2==fNumHiGainSamples) ? 0:1;
+              Byte_t *ptr = pixel.GetHiGainSamples();
+        const Byte_t *end = ptr + fRawEvt->GetNumHiGainSamples()-shift;
 
         const Float_t higainped = CalcHiGainMean(ptr, end);
@@ -118 +134 @@
         const UInt_t pixid = pixel.GetPixelId();
         MPedestalPix &pix = (*fPedestals)[pixid];
-
+        
+        // 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);
-
-    }
+    }
+
+
+    fCounter++;
 
     fPedestals->SetReadyToSave();
@@ -128 +151 @@
 }
 
+Int_t MPedCalcPedRun::PostProcess()
+{
+// Compute pedestals and rms from the whole run
+
+    MRawEvtPixelIter pixel(fRawEvt);
+
+    while (pixel.Next())
+    {
+        const UInt_t pixid = pixel.GetPixelId();
+        MPedestalPix &pix = (*fPedestals)[pixid];
+
+        const Int_t N = fCounter;
+        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);
+
+    }
+    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/fNumHiGainSamples;
+   
+    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;
 }
 
@@ -142 +209 @@
 {
     Float_t rms = 0;
+    Byte_t EvenNumSamples=(fNumHiGainSamples/2*2==fNumHiGainSamples) ? fNumHiGainSamples:fNumHiGainSamples-1;
 
     do
@@ -150 +218 @@
     } while (++ptr != end);
 
-    return sqrt(rms/fNumHiGainSamples);
-}
+    return sqrt(rms/EvenNumSamples);
+}
+
+
+
 /*
 Float_t MPedCalcPedRun::CalcHiGainMeanErr(Float_t higainrms) const