Changeset 3701 for trunk/MagicSoft/Mars/manalysis

Timestamp:

04/09/04 19:17:17 (21 years ago)

Author:

gaug

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/manalysis

Files:

• MPedCalcPedRun.cc (modified) (6 diffs)
• MPedCalcPedRun.h (modified) (1 diff)

trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.cc

@@ -96 +96 @@
 #include "MExtractedSignalCam.h"
 
+#include "MGeomPix.h"
+#include "MGeomCam.h"
 
 #include "MGeomCamMagic.h"
@@ -142 +144 @@
 Int_t MPedCalcPedRun::PreProcess( MParList *pList )
 {
-    Clear();
-
-    fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
-    if (!fRawEvt)
-    {
-        *fLog << err << "MRawEvtData not found... aborting." << endl;
-        return kFALSE;
-    }
-
-    fPedestals = (MPedestalCam*)pList->FindCreateObj("MPedestalCam");
-    if (!fPedestals)
-        return kFALSE;
-
-    return kTRUE;
+
+  Clear();
+  
+  fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
+  if (!fRawEvt)
+    {
+      *fLog << err << "MRawEvtData 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;
 }
 
@@ -175 +185 @@
     }
     else
-        if (runheader->IsMonteCarloRun())
-            return kTRUE;
-
-    Int_t n = fPedestals->GetSize();
+      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(n);
-        fSumx2.Set(n);
+        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();
@@ -204 +224 @@
 Int_t MPedCalcPedRun::Process()
 {
-    MRawEvtPixelIter pixel(fRawEvt);
-
-    while (pixel.Next())
-    {
-              Byte_t *ptr = pixel.GetHiGainSamples();
-        const Byte_t *end = ptr + fNumHiGainSamples;
-
-        UInt_t sum = 0;
-        UInt_t sqr = 0;
+
+  MRawEvtPixelIter pixel(fRawEvt);
+  
+  while (pixel.Next())
+    {
+
+      const UInt_t idx    = pixel.GetPixelId();
+      const UInt_t aidx   = (*fGeom)[idx].GetAidx();
+      const UInt_t sector = (*fGeom)[idx].GetSector();      
+
+      Byte_t *ptr = pixel.GetHiGainSamples();
+      const Byte_t *end = ptr + fNumHiGainSamples;
+      
+      UInt_t sum = 0;
+      UInt_t sqr = 0;
         
-        do
-          {
-            sum += *ptr;
-            sqr += *ptr * *ptr;
-          }
-        while (++ptr != end);
-        
-        const Float_t msum = (Float_t)sum;
-
-        const UInt_t idx = pixel.GetPixelId();
-        //
-        // These three lines have been uncommented by Markus Gaug
-        // 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.));
-        //      (*fPedestals)[idx].Set(higainped, higainrms);
-        
-        fSumx[idx]  += msum;
-        //
-        // The old version:
-        //
-        //       const Float_t msqr = (Float_t)sqr;
-        //      fSumx2[idx] += msqr;
-        //
-        // The new version:
-        //
-        fSumx2[idx] += msum*msum;
-
-    }
-    
-    fPedestals->SetReadyToSave();
-    fNumSamplesTot += fNumHiGainSamples;
-
-    
-    return kTRUE;
+      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/fNumHiGainSamples;
+      //        const Float_t higainrms = TMath::Sqrt((msqr-msum*msum/fNumHiGainSamples)/(fNumHiGainSamples-1.));
+      //        (*fPedestals)[idx].Set(higainped, higainrms);
+      
+      fSumx[idx]          += msum;
+      fAreaSumx[aidx]     += msum;
+      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;
+      fAreaSumx2[aidx]     += sqrsum;
+      fSectorSumx2[sector] += sqrsum;      
+    }
+  
+  fPedestals->SetReadyToSave();
+  fNumSamplesTot += fNumHiGainSamples;
+  
+  
+  return kTRUE;
 }
 
@@ -267 +296 @@
   while (pixel.Next())
     {
-      const Int_t pixid = pixel.GetPixelId();
-      
+
+      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);
@@ -290 +325 @@
 
     }
+
+  //
+  // Loop over the (two) area indices to get the averaged pedestal per aidx
+  //
+  for (Int_t aidx=0; aidx<fAreaValid.GetSize(); aidx++)
+    {
+      
+      const Int_t   napix = fAreaValid.At(aidx);
+      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 /= (Float_t)fNumHiGainSamples;
+      // 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);
+
+      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++)
+    {
+      
+      const Int_t   nspix = fSectorValid.At(sector);
+      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 /= (Float_t)fNumHiGainSamples;
+      // 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);
+
+      fPedestals->GetAverageSector(sector).Set(higainped, higainrms);
+    }
   
   fPedestals->SetTotalEntries(fNumSamplesTot);

trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.h

@@ -10 +10 @@
 #endif
 
+#ifndef ROOT_TArrayI
+#include <TArrayI.h>
+#endif
+
 class MRawEvtData;
 class MPedestalCam;
-
+class MGeomCam;
 class MPedCalcPedRun : public MTask
 {
-    Byte_t fNumHiGainSamples;
-    UInt_t fNumSamplesTot;
 
-    MRawEvtData  *fRawEvt;     // raw event data (time slices)
-    MPedestalCam *fPedestals;  // Pedestals of all pixels in the camera
+  Byte_t fNumHiGainSamples;
+  UInt_t fNumSamplesTot;
+  
+  MRawEvtData  *fRawEvt;     // raw event data (time slices)
+  MPedestalCam *fPedestals;  // Pedestals of all pixels in the camera
+  MGeomCam     *fGeom;       // Camera geometry
+  
+  TArrayF fSumx;         // sum of values
+  TArrayF fSumx2;        // sum of squared values
+  TArrayF fAreaSumx;     // averaged sum of values per area idx
+  TArrayF fAreaSumx2;    // averaged sum of squared values per area idx
+  TArrayI fAreaValid;    // number of valid pixel with area idx  
+  TArrayF fSectorSumx;   // averaged sum of values per sector 
+  TArrayF fSectorSumx2;  // averaged sum of squared values 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();
+  
+public:
 
-    TArrayF fSumx;   // sum of values
-    TArrayF fSumx2;  // sum of squared values
-
-    Bool_t ReInit(MParList *pList);
-
-    Int_t PreProcess(MParList *pList);
-    Int_t Process();
-    Int_t PostProcess();
-
-public:
-    MPedCalcPedRun(const char *name=NULL, const char *title=NULL);
-
-    void Clear(const Option_t *o="");
-    void SetNumHiGainSamples(const Byte_t n)      { fNumHiGainSamples = n;   }
-    
-    ClassDef(MPedCalcPedRun, 0)   // Task to calculate pedestals from pedestal runs raw data
+  MPedCalcPedRun(const char *name=NULL, const char *title=NULL);
+  
+  void Clear(const Option_t *o="");
+  void SetNumHiGainSamples(const Byte_t n)      { fNumHiGainSamples = n;   }
+  
+  ClassDef(MPedCalcPedRun, 0)   // Task to calculate pedestals from pedestal runs raw data
 };