Changeset 4716 for trunk/MagicSoft/Mars/manalysis

Timestamp:

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

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/manalysis

Files:

• MCameraData.cc (modified) (6 diffs)
• MCameraData.h (modified) (2 diffs)

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