Changeset 3476 for trunk/MagicSoft/Mars/mbadpixels

Timestamp:

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

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mbadpixels

Files:

• MBadPixelsCalc.cc (modified) (4 diffs)
• MBadPixelsCalc.h (modified) (1 diff)
• MBadPixelsCam.cc (modified) (2 diffs)
• MBadPixelsPix.h (modified) (2 diffs)
• MBadPixelsTreat.cc (modified) (8 diffs)
• MBadPixelsTreat.h (modified) (1 diff)
• MMcBadPixelsSet.cc (modified) (1 diff)

trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCalc.cc

@@ -56 +56 @@
 #include "MBadPixelsCalc.h"
 
+#include <TArrayD.h>
+
 #include "MLog.h"
 #include "MLogManip.h"
@@ -62 +64 @@
 
 #include "MGeomCam.h"
+#include "MGeomPix.h"
+
 #include "MSigmabar.h"
 
@@ -150 +154 @@
 
         if (pixPedRms*nratio > fPedestalLevel * meanPedRMS || pixPedRms == 0)
-            (*fBadPixels)[i].SetUnsuitableEvt();
-    }
-}
+            (*fBadPixels)[i].SetUnsuitable(MBadPixelsPix::kUnsuitableEvt);
+    }
+}
+
+// --------------------------------------------------------------------------
+// Check the pedestal Rms of the pixels: compute with 2 iterations the mean 
+// for inner and outer pixels. Set as blind the pixels with too small or 
+// too high pedestal Rms with respect to the mean.
+// 
+Bool_t MBadPixelsCalc::CheckPedestalRms() const
+{
+    const Int_t entries = fPedPhotCam->GetSize();
+
+    const Int_t na = fGeomCam->GetNumAreas();
+
+    TArrayD meanrms(na);
+    TArrayI npix(na);
+
+    for (Int_t i=0; i<entries; i++)
+    {
+        const Double_t rms = (*fPedPhotCam)[i].GetRms();
+
+        if (rms<=0 || rms>=200*fGeomCam->GetPixRatioSqrt(i))
+            continue;
+
+        const Byte_t aidx = (*fGeomCam)[i].GetAidx();
+
+        meanrms[aidx] += rms;
+        npix[aidx]++;
+    }
+
+    //if no pixel has a minimum signal, return
+    for (int i=0; i<na; i++)
+    {
+        if (npix[i]==0 || meanrms[i]==0)
+        {
+            //fErrors[1]++;          //no valid Pedestals Rms
+            return kFALSE;
+        }
+
+        meanrms[i] /= npix[i];
+        npix[i]=0;
+    }
+
+    TArrayD meanrms2(na);
+    for (Int_t i=0; i<entries; i++)
+    {
+        const Double_t rms = (*fPedPhotCam)[i].GetRms();
+        const Byte_t  aidx = (*fGeomCam)[i].GetAidx();
+
+        //Calculate the corrected means:
+
+        if (rms<=0.5*meanrms[aidx] || rms>=1.5*meanrms[aidx])
+            continue;
+
+        meanrms2[aidx] += rms;
+        npix[aidx]++;
+    }
+
+    //if no pixel has a minimum signal, return
+    for (int i=0; i<na; i++)
+    {
+        if (npix[i]==0 || meanrms2[i]==0)
+        {
+            //fErrors[1]++;          //no valid Pedestals Rms
+            return kFALSE;
+        }
+
+        meanrms2[i] /= npix[i];
+    }
+
+    Int_t bads = 0;
+
+    //Blind the Bad Pixels
+    for (Int_t i=0; i<entries; i++)
+    {
+        const Double_t rms = (*fPedPhotCam)[i].GetRms();
+        const Byte_t  aidx = (*fGeomCam)[i].GetAidx();
+
+        if (rms>meanrms2[aidx]/3 && rms<=meanrms2[aidx]*3)
+            continue;
+
+        (*fBadPixels)[i].SetUnsuitable(MBadPixelsPix::kUnsuitableEvt);
+        bads++;
+    }
+
+    // Check if the number of pixels to blind is > 60% of total number of pixels
+    //
+    if (bads>0.6*entries)
+    {
+        //fErrors[2]++;
+        return kFALSE;
+    }
+
+    return kTRUE;
+}
+
+
 // --------------------------------------------------------------------------
 //
@@ -159 +258 @@
 {
     if (fPedestalLevel>0)
-        CheckPedestalRMS();
+        CheckPedestalRms();
 
     return kTRUE;

trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCalc.h

@@ -23 +23 @@
 
     void CheckPedestalRMS() const;
+    Bool_t CheckPedestalRms() const;
 
     Int_t PreProcess(MParList *pList);

trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCam.cc

@@ -280 +280 @@
             InitSize(idx+1);
 
-        (*this)[idx].SetUnsuitableRun();
+        (*this)[idx].SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
     }
 }
@@ -295 +295 @@
 
     for (int i=0; i<GetSize(); i++)
-        if ((*this)[i].IsUnsuitableRun())
+        if ((*this)[i].IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
             fout << " " << i;
 

trunk/MagicSoft/Mars/mbadpixels/MBadPixelsPix.h

@@ -52 +52 @@
 
     // Setter
-    void   SetUnsuitableRun() { fInfo[0] |=  kUnsuitableRun; }
-    void   SetSuitableRun()   { fInfo[0] &= ~kUnsuitableRun; }
-
-    void   SetUnsuitableEvt() { fInfo[0] |=  kUnsuitableEvt; }
-    void   SetSuitableEvt()   { fInfo[0] &= ~kUnsuitableEvt; }
-
-    void   SetUnreliableRun() { fInfo[0] |=  kUnreliableRun; }
-    void   SetReliableRun()   { fInfo[0] &= ~kUnreliableRun; }
+    void SetUnsuitable(UnsuitableType_t typ) { fInfo[0] |= typ; }
 
     // Calibration
@@ -78 +71 @@
 
     // Getter
-    Bool_t IsUnsuitableRun() const { return   fInfo[0]&kUnsuitableRun; }
-    Bool_t IsSuitableRun() const   { return !(fInfo[0]&kUnsuitableRun); }
-
-    Bool_t IsUnsuitableEvt() const { return   fInfo[0]&kUnsuitableEvt; }
-    Bool_t IsSuitableEvt() const   { return !(fInfo[0]&kUnsuitableEvt); }
-
-    Bool_t IsUnreliableRun() const { return   fInfo[0]&kUnreliableRun; }
-    Bool_t IsReliableRun() const   { return !(fInfo[0]&kUnreliableRun); }
+    Bool_t IsUnsuitable(UnsuitableType_t typ) const { return fInfo[0]&typ; }
 
     Bool_t IsOK() const  { return fInfo[0]==0; }

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;
 }

trunk/MagicSoft/Mars/mbadpixels/MBadPixelsTreat.h

@@ -27 +27 @@
     };
 
-    static Double_t Sqr(Double_t x) { return x*x; }
+    static Double_t Pow2(Double_t x) { return x*x; }
+
+    void InterpolateSignal() const;
+    void InterpolatePedestals() const;
 
     void  Interpolate() const;

trunk/MagicSoft/Mars/mbadpixels/MMcBadPixelsSet.cc

@@ -131 +131 @@
     // Case for Crab Nebula FOV
     //
-    (*fBadPixels)[400].SetUnsuitableRun();
-    (*fBadPixels)[401].SetUnsuitableRun();
-    (*fBadPixels)[402].SetUnsuitableRun();
-    (*fBadPixels)[437].SetUnsuitableRun();
-    (*fBadPixels)[438].SetUnsuitableRun();
-    (*fBadPixels)[439].SetUnsuitableRun();
+    (*fBadPixels)[400].SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
+    (*fBadPixels)[401].SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
+    (*fBadPixels)[402].SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
+    (*fBadPixels)[437].SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
+    (*fBadPixels)[438].SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
+    (*fBadPixels)[439].SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
 
     *fLog << inf;