Changeset 1460 for trunk/MagicSoft/Mars/manalysis/MHillasExt.cc

Timestamp:

07/31/02 11:27:21 (22 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/manalysis/MHillasExt.cc (modified) (7 diffs)

trunk/MagicSoft/Mars/manalysis/MHillasExt.cc

@@ -72 +72 @@
     MHillas::Reset();
 
-    fConc    = 0;
-    fConc1   = 0; 
-    fAsym    = 0;
-    fM3Long  = 0;
-    fM3Trans = 0;
+    fConc    = -1;
+    fConc1   = -1;
+    fAsym    =  0;
+    fM3Long  =  0;
+    fM3Trans =  0;
 }
 
@@ -107 +107 @@
     // --------------------------------------------
     //
-    //   loop to get third moments along ellipse axes and two max pixels
-    //
-    float m3x = 0;
-    float m3y = 0;
-
-    float maxpix1 = 0;                                               // [#phot]
-    float maxpix2 = 0;                                               // [#phot]
-
-    float maxpixx = 0;                                               // [mm]
-    float maxpixy = 0;                                               // [mm]
+    //  loop to get third moments along ellipse axes and two max pixels
+    //
+    //  For the moments double precision is used to make sure, that
+    //  the complex matrix multiplication and sum is evaluated correctly.
+    //
+    Double_t m3x = 0;
+    Double_t m3y = 0;
+
+    Float_t maxpix1 = 0;                                               // [#phot]
+    Float_t maxpix2 = 0;                                               // [#phot]
+
+    Int_t maxpixid = 0;
 
     const UInt_t npixevt = evt.GetNumPixels();
+
+    const Float_t A0 = geom[0].GetA();
 
     for (UInt_t i=0; i<npixevt; i++)
@@ -127 +131 @@
 
         const MGeomPix &gpix = geom[pix.GetPixId()];
-        const float dx = gpix.GetX() - GetMeanX();                   // [mm]
-        const float dy = gpix.GetY() - GetMeanY();                   // [mm]
-
-        const float nphot = pix.GetNumPhotons();                     // [1]
-
-        const float dzx =  GetCosDelta()*dx + GetSinDelta()*dy;      // [mm]
-        const float dzy = -GetSinDelta()*dx + GetCosDelta()*dy;      // [mm]
+        const Double_t dx = gpix.GetX() - GetMeanX();                // [mm]
+        const Double_t dy = gpix.GetY() - GetMeanY();                // [mm]
+
+        Double_t nphot = pix.GetNumPhotons();                        // [1]
+
+        const Double_t dzx =  GetCosDelta()*dx + GetSinDelta()*dy;   // [mm]
+        const Double_t dzy = -GetSinDelta()*dx + GetCosDelta()*dy;   // [mm]
 
         m3x += nphot * dzx*dzx*dzx;                                  // [mm^3]
         m3y += nphot * dzy*dzy*dzy;                                  // [mm^3]
 
+        /*
+         //
+         // count number of photons in pixels at the edge of the camera
+         //
+         if (gpix.IsInOutermostRing())
+            edgepix1 += nphot;
+         if (gpix.IsInOuterRing())
+            edgepix2 += nphot;
+         */
+
+        //
+        // Now we are working on absolute values of nphot, which
+        // must take pixel size into account
+        //
+        const Double_t r = A0/gpix.GetA();
+        nphot *= r;
+
         if (nphot>maxpix1)
         {
-            maxpix2 = maxpix1;
-            maxpix1 = nphot;                                         // [1]
-            maxpixx = dx + GetMeanX();                               // [mm]
-            maxpixy = dy + GetMeanY();                               // [mm]
+            maxpix2  = maxpix1;
+            maxpix1  = nphot;                                        // [1]
+            maxpixid = pix.GetPixId();
             continue;                                                // [1]
         }
@@ -149 +169 @@
         if (nphot>maxpix2)
             maxpix2 = nphot;                                         // [1]
+
+        /*
+         //
+         // get sums for calculating fAsymna
+         // the outer pixels are 4 times as big (in area)
+         // as the inner pixels !
+         //
+         const Double_t dummy = pow(nphot, na)/r;
+
+         sna +=     dummy;
+         xna += dzx*dummy;
+
+         sna1 += sna/nphot;
+         xna1 += xna/nphot;
+
+         //
+         // forward-backward asymmetry
+         //
+         fb += dzx<0 ? -nphot: nphot;
+         */
     }
 
-    fAsym  = (GetMeanX()-maxpixx)*GetCosDelta() + (GetMeanY()-maxpixy)*GetSinDelta(); // [mm]
+    const MGeomPix &maxpix = geom[maxpixid];
+
+    fAsym  =
+        (GetMeanX()*2-maxpix.GetX())*GetCosDelta() +
+        (GetMeanY()*2-maxpix.GetY())*GetSinDelta();                  // [mm]
+
     fConc  = (maxpix1+maxpix2)/GetSize();                            // [ratio]
     fConc1 = maxpix1/GetSize();                                      // [ratio]
 
+    /*
+
+     //
+     // power na for calculating fAsymna;
+     // the value 1.5 was suggested by Thomas Schweizer
+     //
+     Double_t na = 1.5;
+
+     fLeakage1 = edgepix1 / GetSize();
+     fLeakage2 = edgepix2 / GetSize();
+     fAsym0    =       fb / GetSize();
+
+     fAsymna   = na * (sna*xna1 - sna1*xna) / (sna*sna);
+     */
+
     //
     // Third moments along axes get normalized
@@ -161 +221 @@
     m3y /= GetSize();
 
-    fM3Long  = m3x<0 ? -pow(-m3x, 1./3) : pow(m3x, 1./3); // [mm]
-    fM3Trans = m3y<0 ? -pow(-m3y, 1./3) : pow(m3y, 1./3); // [mm]
+    fM3Long  = m3x<0 ? -pow(-m3x, 1./3) : pow(m3x, 1./3);          // [mm]
+    fM3Trans = m3y<0 ? -pow(-m3y, 1./3) : pow(m3y, 1./3);          // [mm]
 
     SetReadyToSave();
@@ -169 +229 @@
 }
 
+/*
 // -------------------------------------------------------------------------
 //
@@ -181 +242 @@
     fin >> fM3Trans;
 }
-
+*/
 // -------------------------------------------------------------------------
 /*