Changeset 3476 for trunk/MagicSoft/Mars/mbadpixels/MBadPixelsTreat.cc

Timestamp:

03/11/04 16:32:11 (21 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mbadpixels/MBadPixelsTreat.cc (modified) (8 diffs)

trunk/MagicSoft/Mars/mbadpixels/MBadPixelsTreat.cc

@@ -36 +36 @@
 //  Input Containers:
 //   MCerPhotEvt
+//   MPedPhotCam
 //   MBadPixelsCam
+//   [MGeomCam]
 //
 //  Output Containers:
@@ -45 +47 @@
 
 #include <fstream>
+
+#include <TArrayD.h>
 
 #include "MLog.h"
@@ -114 +118 @@
     if (!fGeomCam && TESTBIT(fFlags, kUseInterpolation))
     {
-        *fLog << err << "No camera geometry available... can't use interpolation." << endl;
+        *fLog << err << "MGeomCam not found... can't use interpolation." << endl;
         return kFALSE;
     }
 
     return kTRUE;
+}
+
+// --------------------------------------------------------------------------
+//
+//  Replaces each pixel (signal, signal error, pedestal, pedestal rms)
+//  by the average of its surrounding pixels.
+//  If TESTBIT(fFlags, kUseCentralPixel) is set the central pixel is also
+//  included.
+//
+void MBadPixelsTreat::InterpolateSignal() const
+{
+    const UShort_t entries = fGeomCam->GetNumPixels();
+
+    //
+    // Create arrays
+    //
+    TArrayD nphot(entries);
+    TArrayD perr(entries);
+ 
+    //
+    // Loop over all pixels
+    //
+    for (UShort_t i=0; i<entries; i++)
+    {
+        //
+        // Check whether pixel with idx i is blind
+        //
+        if (!(*fBadPixels)[i].IsBad())
+            continue;
+
+        //
+        // Get a pointer to this pixel. If it is not yet existing
+        // create a new entry for this pixel in MCerPhotEvt
+        //
+        MCerPhotPix *pix = fEvt->GetPixById(i);
+        if (!pix)
+            pix = fEvt->AddPixel(i, 0, 0);
+
+        //
+        // Get the corresponding geometry and pedestal
+        //
+        const MGeomPix &gpix = (*fGeomCam)[i];
+
+        // Do Not-Use-Central-Pixel
+        const Bool_t nucp = !TESTBIT(fFlags, kUseCentralPixel);
+
+        Int_t num = nucp ? 0 : 1;
+
+        nphot[i]  = nucp ? 0 : pix->GetNumPhotons();
+        perr[i]   = nucp ? 0 : Pow2(pix->GetErrorPhot());
+
+        //
+        // The values are rescaled to the small pixels area for the right comparison
+        //
+        const Double_t ratio = fGeomCam->GetPixRatio(i);
+
+        nphot[i] *= ratio;
+        perr[i]  *= ratio;
+
+        //
+        // Loop over all its neighbors
+        //
+        const Int_t n = gpix.GetNumNeighbors();
+        for (int j=0; j<n; j++)
+        {
+            const UShort_t nidx = gpix.GetNeighbor(j);
+
+            //
+            // Do not use blind neighbors
+            //
+            if (!(*fBadPixels)[i].IsBad())
+                continue;
+
+            //
+            // Check whether the neighbor has a signal stored
+            //
+            const MCerPhotPix *evtpix = fEvt->GetPixById(nidx);
+            if (!evtpix)
+                continue;
+
+            //
+            // Get the geometry for the neighbor
+            //
+            const Double_t nratio = fGeomCam->GetPixRatio(nidx);
+
+            //
+            //The error is calculated as the quadratic sum of the errors
+            //
+            nphot[i] += nratio*evtpix->GetNumPhotons();
+            perr[i]  += nratio* Pow2(evtpix->GetErrorPhot());
+
+            num++;
+        }
+
+        //
+        // Now the mean is calculated and the values rescaled back
+        // to the pixel area
+        //
+        nphot[i] /= (num*ratio);
+        perr[i]   = TMath::Sqrt(perr[i]/(num*ratio));
+ 
+        // Check if there are enough neighbors to calculate the mean
+        // If not, unmap the pixel. The maximum number of blind neighbors
+        // should be 2
+        if (num<3)
+        {
+            pix->SetPixelUnmapped();
+            continue;
+        }
+
+        pix->Set(nphot[i], perr[i]);
+    }
+}
+
+// --------------------------------------------------------------------------
+//
+void MBadPixelsTreat::InterpolatePedestals() const
+{
+    const Int_t entries = fPedPhot->GetSize();
+
+    TArrayD ped(entries);
+    TArrayD rms(entries);
+
+    //
+    // Loop over all pixels
+    //
+    for (UShort_t i=0; i<entries; i++)
+    {
+        //
+        // Check whether pixel with idx i is blind
+        //
+        if (!(*fBadPixels)[i].IsBad())
+            continue;
+
+        //
+        // Get the corresponding geometry and pedestal
+        //
+        const MGeomPix    &gpix = (*fGeomCam)[i];
+        const MPedPhotPix &ppix = (*fPedPhot)[i];
+
+        // Do Not-Use-Central-Pixel
+        const Bool_t nucp = !TESTBIT(fFlags, kUseCentralPixel);
+
+        Int_t num = nucp ? 0 : 1;
+
+        ped[i] = nucp ? 0 : ppix.GetMean();
+        rms[i] = nucp ? 0 : Pow2(ppix.GetRms());
+
+        //
+        // The values are rescaled to the small pixels area for the right comparison
+        //
+        const Double_t ratio = fGeomCam->GetPixRatio(i);
+
+        ped[i] *= ratio;
+        rms[i] *= ratio;
+
+        //
+        // Loop over all its neighbors
+        //
+        const Int_t n = gpix.GetNumNeighbors();
+        for (int j=0; j<n; j++)
+        {
+            const UShort_t nidx = gpix.GetNeighbor(j);
+
+            //
+            // Do not use blind neighbors
+            //
+            if (!(*fBadPixels)[i].IsBad())
+                continue;
+
+            //
+            // Get the geometry for the neighbor
+            //
+            const Double_t    nratio = fGeomCam->GetPixRatio(nidx);
+            const MPedPhotPix &nppix = (*fPedPhot)[nidx];
+
+            //
+            //The error is calculated as the quadratic sum of the errors
+            //
+            ped[i] += (nratio*nppix.GetMean());
+            rms[i] += nratio*Pow2(nppix.GetRms());
+
+            num++;
+        }
+
+        // Check if there are enough neighbors to calculate the mean
+        // If not, unmap the pixel. The minimum number of good neighbors
+        // should be 3
+        if (num < 3)
+        {
+            MCerPhotPix *pix =fEvt->GetPixById(i);
+            if (!pix)
+                pix = fEvt->AddPixel(i, 0, 0);
+            pix->SetPixelUnmapped();
+            continue;
+        }
+
+        //
+        // Now the mean is calculated and the values rescaled back
+        // to the pixel area
+        //
+        ped[i] /=  (num*ratio);
+        rms[i]  = TMath::Sqrt(rms[i]/(num*ratio));
+
+        (*fPedPhot)[i].Set(ped[i], rms[i]);
+    }
 }
 
@@ -177 +387 @@
         nphot[i]  = nucp ? 0 : pix->GetNumPhotons();
         ped[i]    = nucp ? 0 : ppix.GetMean();
-        perr[i]   = nucp ? 0 : Sqr(pix->GetErrorPhot());
-        pedrms[i] = nucp ? 0 : Sqr(ppix.GetRms());
+        perr[i]   = nucp ? 0 : Pow2(pix->GetErrorPhot());
+        pedrms[i] = nucp ? 0 : Pow2(ppix.GetRms());
 
         //
@@ -187 +397 @@
         nphot[i]  *= ratio;
         ped[i]    *= ratio;
-        perr[i]   *= Sqr(ratio);
-        pedrms[i] *= Sqr(pedrms[i]);
+        perr[i]   *= Pow2(ratio);
+        pedrms[i] *= Pow2(pedrms[i]);
 
         //
@@ -222 +432 @@
             ped[i]    += nratio*nppix.GetMean();
             nphot[i]  += nratio*evtpix->GetNumPhotons();
-            perr[i]   += Sqr(nratio*evtpix->GetErrorPhot());
-            pedrms[i] += Sqr(nratio*nppix.GetRms());
+            perr[i]   += Pow2(nratio*evtpix->GetErrorPhot());
+            pedrms[i] += Pow2(nratio*nppix.GetRms());
 
             num++;
@@ -294 +504 @@
 Int_t MBadPixelsTreat::Process()
 {
+    /*
+    if (TESTBIT(fFlags, kCheckPedestalRms))
+    {
+        // if the number of blind pixels is too high, do not interpolate
+       if (CheckPedestalRms()==kFALSE)
+           return kTRUE;
+
+       if (TESTBIT(fFlags, kUseInterpolation))
+           InterpolatePedestals();
+    }
+    */
+
     if (TESTBIT(fFlags, kUseInterpolation) && fGeomCam)
         Interpolate();
@@ -299 +521 @@
         Unmap();
 
+
+    //fErrors[0]++;
+
     return kTRUE;
 }