Changeset 4716 for trunk/MagicSoft/Mars

Timestamp:

08/23/04 17:49:53 (20 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (3 diffs)
• NEWS (modified) (1 diff)
• manalysis/MCameraData.cc (modified) (6 diffs)
• manalysis/MCameraData.h (modified) (2 diffs)
• mbase/MMath.cc (modified) (3 diffs)
• mbase/MMath.h (modified) (1 diff)
• mhist/MHEvent.cc (modified) (2 diffs)
• mhist/MHEvent.h (modified) (1 diff)
• mimage/MImgCleanStd.cc (modified) (1 diff)
• mjobs/MSequence.h (modified) (1 diff)
• msignal/MArrivalTime.h (modified) (1 diff)

trunk/MagicSoft/Mars/Changelog

@@ -79 +79 @@
    * mhist/MHEvent.[h,cc]:
      - added new option for island index
+     - added kEvtCleaningData
 
    * mimage/MImgCleanStd.[h,cc]:
@@ -85 +86 @@
      - added ReadEnv
      - updated StreamPrimitive
+     - added new cleaning method (probability cleaning)
 
    * mimage/Makefile:
@@ -142 +144 @@
    * msignal/MArrivalTime.h:
      - removed operator()
+     - added operator[] const
+
+   * manalysis/MCameraData.[h,cc]:
+     - added algorithm for 'Probability cleaning'
+
+   * mbase/MMath.[h,cc]:
+     - added GaussProb
+
+   * mjobs/MSequence.h:
+     - added IsValid
 
 

trunk/MagicSoft/Mars/NEWS

@@ -25 +25 @@
    - implemented new image parameters displaying the number of islands,
      saturated hi-gain and lo-gain pixels (MImagePar, MHImagePar)
+
+   - event display in executable changed to support also calibrated files
+     (done with MJCalibrateSignal)
+
+   - added a cleaning which takes signal height _and_ arrival time into
+     account: probability cleaning (for more details see MImgCleanStd)
 
 

trunk/MagicSoft/Mars/manalysis/MCameraData.cc

@@ -34 +34 @@
 #include "MCameraData.h"
 
+#include "MMath.h"
+
+#include "MLog.h"
+#include "MLogManip.h"
+
+#include "MGeomCam.h"
+#include "MGeomPix.h"
+
+#include "MPedPhotCam.h"
+#include "MPedPhotPix.h"
+
 #include "MCerPhotEvt.h"
 #include "MCerPhotPix.h"
 
-#include "MGeomCam.h"
-#include "MGeomPix.h"
-
-#include "MLog.h"
-#include "MLogManip.h"
-
-#include "MPedPhotCam.h"
-#include "MPedPhotPix.h"
-
 #include "MSigmabar.h"
+#include "MArrivalTime.h"
 
 ClassImp(MCameraData);
@@ -101 +104 @@
     const Int_t n = geom.GetNumPixels();
 
-    fData.Set(n);
-    fData.Reset();
-
-    fValidity.Set(n);
-    fValidity.Reset();
-
-    const Int_t entries = evt.GetNumPixels();
-
-    //
-    // check the number of all pixels against the noise level and
-    // set them to 'unused' state if necessary
-    // 
+    // Reset arrays
+    fData.Set(n);
+    fData.Reset();
+
+    fValidity.Set(n);
+    fValidity.Reset();
+
+    const Int_t entries = evt.GetNumPixels();
+
+    // calculate cleaning levels
     for (Int_t i=0; i<entries; i++)
     {
@@ -149 +150 @@
     const Int_t n = geom.GetNumPixels();
 
-    fData.Set(n);
-    fData.Reset();
-
-    fValidity.Set(n);
-    fValidity.Reset();
-
-    const Int_t entries   = evt.GetNumPixels();
+    // reset arrays
+    fData.Set(n);
+    fData.Reset();
+
+    fValidity.Set(n);
+    fValidity.Reset();
+
+    // check validity of rms with area index 0
     const Float_t anoise0 = cam.GetArea(0).GetRms();
     if (anoise0<=0)
         return;
 
-    //
-    // check the number of all pixels against the noise level and
-    // set them to 'unused' state if necessary
-    // 
+    // calculate cleaning levels
+    const Int_t entries = evt.GetNumPixels();
     for (Int_t i=0; i<entries; i++)
     {
@@ -200 +200 @@
     const Int_t n = geom.GetNumPixels();
 
-    fData.Set(n);
-    fData.Reset();
-
-    fValidity.Set(n);
-    fValidity.Reset();
-
+    // reset arrays
+    fData.Set(n);
+    fData.Reset();
+
+    fValidity.Set(n);
+    fValidity.Reset();
+
+    // check validity of noise
     if (noise<=0)
         return;
 
-    const Int_t entries = evt.GetNumPixels();
-
-    //
-    // check the number of all pixels against the noise level and
-    // set them to 'unused' state if necessary
-    // 
+    // calculate cleaning levels
+    const Int_t entries = evt.GetNumPixels();
     for (Int_t i=0; i<entries; i++)
     {
@@ -244 +242 @@
     const Int_t n = geom.GetNumPixels();
 
-    fData.Set(n);
-    fData.Reset();
-
-    fValidity.Set(n);
-    fValidity.Reset();
-
-    const Int_t   entries = evt.GetNumPixels();
-    const Float_t noise0  = cam.GetArea(0).GetRms();
+    // reset arrays
+    fData.Set(n);
+    fData.Reset();
+
+    fValidity.Set(n);
+    fValidity.Reset();
+
+    // check validity of rms with area index 0
+    const Float_t noise0 = cam.GetArea(0).GetRms();
     if (noise0<=0)
         return;
 
-    //
-    // check the number of all pixels against the noise level and
-    // set them to 'unused' state if necessary
-    // 
+    // calculate cleaning levels
+    const Int_t entries = evt.GetNumPixels();
     for (Int_t i=0; i<entries; i++)
     {
@@ -280 +277 @@
 // --------------------------------------------------------------------------
 //
+// Function to calculate the cleaning level for all pixels in a given event.
+// The level is the probability that the signal is a real signal (taking
+// the signal height and the fluctuation of the background into account)
+// times one minus the probability that the signal is a background
+// fluctuation (calculated from the time spread of arrival times
+// around the pixel with the highest signal)
+//
+// FIXME: Should the check noise<=0 be replaced by MBadPixels?
+//
+void MCameraData::CalcCleaningProbability(const MCerPhotEvt &evt, const MPedPhotCam &pcam,
+                                          const MGeomCam &geom,   const MArrivalTime *tcam)
+{
+    const Int_t n = geom.GetNumPixels();
+
+    // Reset arrays
+    fData.Set(n);
+    fData.Reset();
+
+    fValidity.Set(n);
+    fValidity.Reset();
+
+    // check validity of noise
+    const Float_t anoise0 = pcam.GetArea(0).GetRms();
+    if (anoise0<=0)
+        return;
+
+    const Int_t entries = evt.GetNumPixels();
+
+    // Find pixel with max signal
+    Int_t maxidx = 0;
+    if (tcam)
+    {
+        // Find pixel enty with maximum signal
+        for (Int_t i=0; i<entries; i++)
+        {
+            const Double_t s0 = evt[i].GetNumPhotons()      * geom.GetPixRatio(i);
+            const Double_t s1 = evt[maxidx].GetNumPhotons() * geom.GetPixRatio(maxidx);
+            if (s0>s1)
+                maxidx = i;
+        }
+
+        // Get software index of pixel
+        maxidx = evt[maxidx].GetPixId();
+    }
+
+    const Double_t timemean = tcam ? (*tcam)[maxidx] : 0;
+    const Double_t timerms  = 0.75; //[slices] rms time spread around highest pixel
+
+    // calculate cleaning levels
+    for (Int_t i=0; i<entries; i++)
+    {
+        const MCerPhotPix &spix = evt[i];
+
+        const Int_t   idx = spix.GetPixId();
+        const Float_t rms = pcam[idx].GetRms();
+        if (rms<=0) // fData[idx]=0, fValidity[idx]=0
+            continue;
+
+        fValidity[idx]=1;
+
+        // get signal and arrival time
+        const UInt_t  aidx     = geom[idx].GetAidx();
+        const Float_t ratio    = pcam.GetArea(aidx).GetRms()/anoise0;
+
+        const Double_t signal  = spix.GetNumPhotons() * geom.GetPixRatio(idx) * ratio / rms;
+        const Double_t time    = tcam ? (*tcam)[idx] : 1;
+
+        // if signal<0 the probability is equal 0
+        if (signal<0)
+            continue;
+
+        // Just for convinience for easy readybility
+        const Double_t means   = 0;
+        const Double_t meant   = timemean;
+
+        const Double_t sigmas  = rms;
+        const Double_t sigmat  = timerms;
+
+        // Get probabilities
+        const Double_t psignal = MMath::GaussProb(signal, sigmas, means);
+        const Double_t pbckgnd = MMath::GaussProb(time,   sigmat, meant);
+
+        // Set probability
+        fData[idx]     = psignal*(1-pbckgnd);
+        fValidity[idx] = 1;
+
+        // Make sure, that we don't run into trouble because of
+        // a little numerical uncertanty
+        if (fData[idx]>1)
+            fData[idx]=1;
+        if (fData[idx]<0)
+            fData[idx]=0;
+    }
+}
+
+// --------------------------------------------------------------------------
+//
 // Returns the contents of the pixel.
 //

trunk/MagicSoft/Mars/manalysis/MCameraData.h

@@ -19 +19 @@
 class MCerPhotEvt;
 class MPedPhotCam;
+class MArrivalTime;
 
 class MCameraData : public MParContainer, public MCamEvent
@@ -38 +39 @@
     void CalcCleaningLevel(const MCerPhotEvt &evt, Double_t noise,
                            const MGeomCam &geom);
-
     void CalcCleaningLevel2(const MCerPhotEvt &evt, const MPedPhotCam &fCam,
-                           const MGeomCam &geom);
-
+                            const MGeomCam &geom);
     void CalcCleaningLevelDemocratic(const MCerPhotEvt &evt, const MPedPhotCam &cam,
                                      const MGeomCam &geom);
+    void CalcCleaningProbability(const MCerPhotEvt &evt, const MPedPhotCam &pcam,
+                                 const MGeomCam &geom,   const MArrivalTime *tcam);
 
     const TArrayD &GetData() const { return fData; }

trunk/MagicSoft/Mars/mbase/MMath.cc

@@ -47 +47 @@
     const Double_t f = s+k*k*b;
 
-    return f==0 ? 0 : (s-b)/TMath::Sqrt(f);
+    return f==0 ? 0 : (s-b)/Sqrt(f);
 }
 
@@ -86 +86 @@
         return -1;
 
-    const Double_t l = s*TMath::Log(s/sum*(alpha+1)/alpha);
-    const Double_t m = b*TMath::Log(b/sum*(alpha+1)      );
+    const Double_t l = s*Log(s/sum*(alpha+1)/alpha);
+    const Double_t m = b*Log(b/sum*(alpha+1)      );
 
-    return l+m<0 ? -1 : TMath::Sqrt((l+m)*2);
+    return l+m<0 ? -1 : Sqrt((l+m)*2);
 }
 
@@ -104 +104 @@
         return 0;
 
-    return TMath::Sign(sig, s-b);
+    return Sign(sig, s-b);
 }
+
+// --------------------------------------------------------------------------
+//
+// Returns: 2/(sigma*sqrt(2))*integral[0,x](exp(-(x-mu)^2/(2*sigma^2)))
+//
+Double_t MMath::GaussProb(Double_t x, Double_t sigma, Double_t mean)
+{
+    static const Double_t sqrt2 = Sqrt(2.);
+    return Erf((x-mean)/(sigma*sqrt2));
+}
+

trunk/MagicSoft/Mars/mbase/MMath.h

@@ -9 +9 @@
 {
 public:
+    static Double_t GaussProb(Double_t x, Double_t sigma, Double_t mean=0);
+
     static Double_t Significance(Double_t s, Double_t b);
     static Double_t SignificanceSym(Double_t s, Double_t b);

trunk/MagicSoft/Mars/mhist/MHEvent.cc

@@ -152 +152 @@
         fHist->SetYTitle("L");
         break;
+    case kEvtCleaningData:
+        fHist->SetName("CleanData");
+        fHist->SetYTitle("L");
+        break;
      case kEvtIdxMax:
         fHist->SetName("Max Slice Idx");
@@ -220 +224 @@
             fHist->SetLevels(lvl);
         }
+        break;
+    case kEvtCleaningData:
+        fHist->SetCamContent(*event, 0);
         break;
     case kEvtIdxMax:

trunk/MagicSoft/Mars/mhist/MHEvent.h

@@ -22 +22 @@
         kEvtSignalRaw, kEvtSignalDensity, kEvtPedestal, 
         kEvtPedestalRMS, kEvtRelativeSignal, kEvtCleaningLevels,
+        kEvtCleaningData,
         kEvtIdxMax, kEvtArrTime, kEvtTrigPix, kEvtIslandIndex
     };

trunk/MagicSoft/Mars/mimage/MImgCleanStd.cc

@@ -549 +549 @@
         fData->CalcCleaningLevelDemocratic(*fEvt, *fPed, *fCam);
         break;
-        /*
     case kProbability:
         fData->CalcCleaningProbability(*fEvt, *fPed, *fCam, fTime);
-        break;*/
+        break;
     default:
         break;

trunk/MagicSoft/Mars/mjobs/MSequence.h

@@ -43 +43 @@
     void Print(Option_t *o="") const;
 
+    Bool_t IsValid() const { return fSequence!=(UInt_t)-1; }
+
     void SetupPedRuns(MDirIter &iter, const char *path=0) const;
     void SetupDatRuns(MDirIter &iter, const char *path=0) const;

trunk/MagicSoft/Mars/msignal/MArrivalTime.h

@@ -40 +40 @@
     const TArrayF &GetDataErr() const { return fDataErr; }
 
-    Double_t operator[](int i) { return fData[i]; }
+    Double_t operator[](int i)       { return fData[i]; }
+    Double_t operator[](int i) const { return fData[i]; }
 
     // Do not do such things! It is highly dangerous to use two very similar operators