Changeset 7611 for trunk/MagicSoft/Mars/manalysis

Timestamp:

03/15/06 16:39:37 (19 years ago)

Author:

snruegam

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/manalysis

Files:

• MCameraData.cc (modified) (2 diffs)
• MCameraData.h (modified) (1 diff)

trunk/MagicSoft/Mars/manalysis/MCameraData.cc

@@ -18 +18 @@
 !   Author(s): Thomas Bretz, 10/2003 <mailto:tbretz@astro.uni-wuerzburg.de>
 !   Author(s): Hendrik Bartko, 08/2004 <mailto:hbartko@mppmu.mpg.de>
-!
-!   Copyright: MAGIC Software Development, 2000-2004
+!   Author(s): Stefan Ruegamer, 03/2006 <mailto:snruegam@astro.uni-wuerzburg.de>
+!
+!   Copyright: MAGIC Software Development, 2000-2006
 !
 !
@@ -380 +381 @@
 }
 
+void MCameraData::CalcCleaningArrivalTime(const MSignalCam &evt, const MGeomCam &geom)
+{
+    const Int_t n = geom.GetNumPixels();
+
+    // Reset arrays
+    fData.Set(n);
+    fData.Reset();
+
+    fValidity.Set(n);
+    fValidity.Reset();
+
+    // check validity of noise
+    const Int_t entries = evt.GetNumPixels();
+    TArrayD u(6), w(6);
+    TArrayI ii(6);
+
+    for (int i=0; i<entries; i++)
+    {
+        if (evt[i].IsPixelUnmapped())
+            continue;
+
+        const Double_t t0 = evt[i].GetArrivalTime();
+        const Double_t s0 = evt[i].GetNumPhotons();
+        const Double_t r0 = geom.GetPixRatio(i);
+        const Double_t x0 = geom[i].GetX();
+        const Double_t y0 = geom[i].GetY();
+
+        const Double_t e0 = s0<0.01 ? 100 : 1./s0;
+
+        const Int_t n2 = geom[i].GetNumNeighbors();
+
+        Int_t cnt2=0;
+        for (int j=0; j<n2; j++)
+        {
+            Int_t idx = geom[i].GetNeighbor(j);
+
+            if (evt[idx].IsPixelUnmapped())
+                continue;
+
+            const Double_t tj = evt[idx].GetArrivalTime()-t0;
+            const Double_t sj = evt[idx].GetNumPhotons();
+            const Double_t rj = geom.GetPixRatio(idx)+r0;
+            const Double_t xj = geom[idx].GetX()-x0;
+            const Double_t yj = geom[idx].GetY()-y0;
+
+            const Double_t d2 = xj*xj+yj*yj;
+
+            const Double_t ej = sj<0.01 ? 100 : 1./sj;
+
+            u[cnt2] =  tj*tj/rj;
+            w[cnt2] =  1./(e0+ej)/d2;         // TYPE I+II
+            cnt2++;
+        }
+
+        TMath::Sort(cnt2, u.GetArray(), ii.GetArray(), kFALSE);
+
+        Double_t sumt = 0;
+        Double_t sumw = 0;
+        for (int j=0; j<TMath::Min(cnt2, 3); j++)
+        {
+            sumt += u[ii[j]]*w[ii[j]];
+            sumw += w[ii[j]];
+        }
+
+        const Double_t wuz = TMath::Sqrt(sumt/sumw);
+
+        if (TMath::IsNaN(wuz))
+            *fLog << warn << "WARNING - MCameraData " << sumt << " " << sumw << endl;
+
+        Double_t val = s0<0 || TMath::IsNaN(wuz) || wuz<1e-200 ? 0 : s0*r0/wuz;
+
+        if ((s0>0 && wuz<1e-200) || val>100)
+            val=100;
+
+        fData[i] = TMath::Log10(val+1)*5;
+
+        if (TMath::IsNaN(fData[i]))
+            //     0                1e-6          0              NaN
+            *fLog << warn << "WARNING - MCameraData " << sumt << " " << sumw << " " << wuz << " " << val << endl;
+
+        fValidity[i] = 1;
+    }
+}
+
 // --------------------------------------------------------------------------
 //

trunk/MagicSoft/Mars/manalysis/MCameraData.h

@@ -42 +42 @@
                                  const MGeomCam &geom);
     void CalcCleaningAbsolute(const MSignalCam &evt, const MGeomCam &geom);
+    void CalcCleaningArrivalTime(const MSignalCam &evt, const MGeomCam &geom);
 
     const TArrayD &GetData() const { return fData; }