Changeset 2502

Timestamp:

11/10/03 16:20:53 (21 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• manalysis/MBlindPixelCalc.cc (modified) (2 diffs)
• manalysis/MBlindPixelCalc.h (modified) (1 diff)
• manalysis/MCerPhotEvt.h (modified) (2 diffs)
• mdatacheck (deleted)
• meventdisp (deleted)
• mmain/MGEvtDisplay.cc (added)
• mmain/MGEvtDisplay.h (added)

trunk/MagicSoft/Mars/Changelog

@@ -116 +116 @@
    * mmain/MCameraDisplay.cc:
      - minor changes
-     
+
+   * manalysis/MBlindPixelCalc.[h,cc]:
+     - corrected interpolation of all values (thanks to Nadia)
+     - fixed algorithm for pixels not existing yet
+
+   * manalysis/MCerPhotEvt.h:
+     - added return value to AddPixel
 
 

trunk/MagicSoft/Mars/manalysis/MBlindPixelCalc.cc

@@ -165 +165 @@
 // --------------------------------------------------------------------------
 //
-//  Replaces each pixel by the average of its surrounding pixels.
+//  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.
 //
+//  NT: Informations about the interpolated pedestals are added. 
+//      When the option Interpolated is called, the blind pixel with the new
+//      values of signal and fluttuation is included in the calculation of
+//      the Image Parameters.
+//
 void MBlindPixelCalc::Interpolate() const
 {
-    const UShort_t entries = fEvt->GetNumPixels();
-
-    Double_t *nphot = new Double_t[entries];
-    Double_t *perr  = new Double_t[entries];
-    Double_t *ped   = new Double_t[entries];
-
-    //
-    // remove the pixels in fPixelsIdx if they are set to be used,
-    // (set them to 'unused' state)
+    const UShort_t entries = fGeomCam->GetNumPixels();
+
+    //
+    // Create arrays
+    //
+    Double_t *nphot  = new Double_t[entries];
+    Double_t *perr   = new Double_t[entries];
+    Double_t *ped    = new Double_t[entries];
+    Double_t *pedrms = new Double_t[entries];
+ 
+    //
+    // Loop over all pixels
     //
     for (UShort_t i=0; i<entries; i++)
     {
-        // FIXME: Loop over blinds instead of all pixels!!!
-        MCerPhotPix &pix = (*fEvt)[i];
-
-        const Int_t idx = pix.GetPixId();
-
-        if (!fPixels->IsBlind(idx))
+        //
+        // Check whether pixel with idx i is blind
+        //
+        if (!fPixels->IsBlind(i))
             continue;
 
-        const MGeomPix &gpix = (*fGeomCam)[idx];
-
-        Int_t num = TESTBIT(fFlags, kUseCentralPixel) ? 1 : 0;
-
-        nphot[i] = TESTBIT(fFlags, kUseCentralPixel) ? pix.GetNumPhotons() : 0;
-        perr[i]  = TESTBIT(fFlags, kUseCentralPixel) ? (*fPed)[idx].GetPedestalRms() : 0;
-        ped[i]   = TESTBIT(fFlags, kUseCentralPixel) ? (*fPed)[idx].GetPedestal() : 0;
-
-        nphot[i] *= fGeomCam->GetPixRatio(idx);
-        // FIXME? perr
-        // FIXME? ped
-
+        //
+        // Get a pointer to this pixel. If it is not yet existing
+        // create a new entry for this pixel in MCerPhotEvt
+        //
+        const MCerPhotPix *pix = fEvt->GetPixById(i);
+        if (!pix)
+            pix = fEvt->AddPixel(i, 0, 0);
+
+        //
+        // Get the corresponding geometry and pedestal
+        //
+        const MGeomPix &gpix = (*fGeomCam)[i];
+        const MPedestalPix &ppix = (*fPed)[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();
+        ped[i]    = nucp ? 0 : ppix.GetPedestal();
+        perr[i]   = nucp ? 0 : Sqr(pix->GetErrorPhot());
+        pedrms[i] = nucp ? 0 : Sqr(ppix.GetPedestalRms());
+
+        //
+        // The values are rescaled to the small pixels area for the right comparison
+        //
+        const Double_t ratio = fGeomCam->GetPixRatio(i);
+
+        nphot[i]  *= ratio;
+        ped[i]    *= ratio;
+        perr[i]   *= Sqr(ratio);
+        pedrms[i] *= Sqr(pedrms[i]);
+
+        //
+        // Loop over all its neighbors
+        //
         const Int_t n = gpix.GetNumNeighbors();
         for (int j=0; j<n; j++)
@@ -208 +240 @@
             const UShort_t nidx = gpix.GetNeighbor(j);
 
+            //
+            // Do not use blind neighbors
+            //
             if (fPixels->IsBlind(nidx))
                 continue;
 
+            //
+            // Check whether the neighbor has a signal stored
+            //
             const MCerPhotPix *evtpix = fEvt->GetPixById(nidx);
             if (!evtpix)
                 continue;
 
-            nphot[i] += evtpix->GetNumPhotons()*fGeomCam->GetPixRatio(nidx);
-            perr[i]  += (*fPed)[nidx].GetPedestalRms();
-            ped[i]   += (*fPed)[nidx].GetPedestal();
-            // FIXME? perr
-            // FIXME? ped
+            //
+            // Get the geometry for the neighbor
+            //
+            const Double_t nratio = fGeomCam->GetPixRatio(nidx);
+            MPedestalPix &nppix = (*fPed)[nidx];
+
+            //
+            //The error is calculated as the quadratic sum of the errors
+            //
+            ped[i]    += nratio*nppix.GetPedestal();
+            nphot[i]  += nratio*evtpix->GetNumPhotons();
+            perr[i]   += Sqr(nratio*evtpix->GetErrorPhot());
+            pedrms[i] += Sqr(nratio*nppix.GetPedestalRms());
 
             num++;
         }
 
-        nphot[i] /= num*fGeomCam->GetPixRatio(idx);
-        perr[i]  /= num/*FIXME:???*/;
-        ped[i]   /= num/*FIXME:???*/;
-    }
-
-    if (TESTBIT(fFlags, kUseInterpolation) && fGeomCam)
-        for (UShort_t i=0; i<entries; i++)
-        {
-            MCerPhotPix &pix = (*fEvt)[i];
-
-            if (fPixels->IsBlind(pix.GetPixId()))
-            {
-                pix.Set(nphot[i], -1);
-                (*fPed)[pix.GetPixId()].Set(ped[i], perr[i]);
-            }
-        }
-
+        //
+        // Now the mean is calculated and the values rescaled back to the pixel area
+        //
+        nphot[i] /= num*ratio;
+        ped[i]   /= num*ratio;
+        perr[i]   = TMath::Sqrt(perr[i]/num)*ratio;
+        pedrms[i] = TMath::Sqrt(pedrms[i]/num)*ratio;
+
+    }
+
+    //
+    // Now the new pixel values are calculated and can be replaced in
+    // the corresponding containers
+    //
+    for (UShort_t i=0; i<entries; i++)
+    {
+        //
+        // Do not use blind neighbors
+        //
+        if (!fPixels->IsBlind(i))
+            continue;
+
+        //
+        // It must exist, we have created it in the loop before.
+        //
+        fEvt->GetPixById(i)->Set(nphot[i], perr[i]);
+        (*fPed)[i].Set(ped[i], pedrms[i]);
+    }
+
+    //
+    // Delete the intermediat arrays
+    //
     delete nphot;
     delete perr;
     delete ped;
+    delete pedrms;
 }
 

trunk/MagicSoft/Mars/manalysis/MBlindPixelCalc.h

@@ -35 +35 @@
     };
 
+    static Double_t Sqr(Double_t x) { return x*x; }
+
     void Interpolate() const;
     void Unmap() const;

trunk/MagicSoft/Mars/manalysis/MCerPhotEvt.h

@@ -32 +32 @@
     //void   InitSize(UInt_t num) { fPixels->Expand(num); }
 
-    void   AddPixel(Int_t idx, Float_t nph, Float_t er)
+    MCerPhotPix *AddPixel(Int_t idx, Float_t nph=0, Float_t er=0)
     {
         //
@@ -47 +47 @@
 
         fLut[idx] = fNumPixels;
-        new ((*fPixels)[fNumPixels++]) MCerPhotPix(idx, nph, er);
+        return new ((*fPixels)[fNumPixels++]) MCerPhotPix(idx, nph, er);
     }