Index: trunk/MagicSoft/Mars/mimage/MHillasExt.cc
===================================================================
--- trunk/MagicSoft/Mars/mimage/MHillasExt.cc	(revision 7784)
+++ trunk/MagicSoft/Mars/mimage/MHillasExt.cc	(revision 8619)
@@ -20,5 +20,5 @@
 !   Author(s): Wolfgang Wittek 06/2002 <mailto:wittek@mppmu.mpg.de>
 !
-!   Copyright: MAGIC Software Development, 2000-2004
+!   Copyright: MAGIC Software Development, 2000-2007
 !
 !
@@ -51,6 +51,12 @@
 // fMaxDist  added. Distance between center and most distant used pixel
 //
-//
-// WARNING: Before you can use fAsym, fM3Long and fM3Trans you must
+// Version 4:
+// ----------
+// fMaxDist      removed.
+// fSlopeLong    added
+// fSlopeTrans   added
+//
+//
+// WARNING: Before you can use fAsym, fSlope* and fM3* you must
 //          multiply by the sign of MHillasSrc::fCosDeltaAlpha
 //
@@ -66,4 +72,5 @@
 #include <TMarker.h>
 #include <TVector2.h>
+#include <TMinuit.h>
 #include <TVirtualPad.h>
 
@@ -97,17 +104,13 @@
 // -------------------------------------------------------------------------
 //
-// Reset values to:
-//    fAsym    =  0;
-//    fM3Long  =  0;
-//    fM3Trans =  0;
-//    fMaxDist =  0;
+// Reset all values to 0
 //
 void MHillasExt::Reset()
 {
-    fAsym    =  0;
-    fM3Long  =  0;
-    fM3Trans =  0;
-
-    fMaxDist =  0;
+    fAsym       =  0;
+    fM3Long     =  0;
+    fM3Trans    =  0;
+    fSlopeLong  =  0;
+    fSlopeTrans =  0;
 }
 
@@ -128,14 +131,23 @@
     //  the complex matrix multiplication and sum is evaluated correctly.
     //
-    Double_t m3x = 0;
-    Double_t m3y = 0;
-
     Int_t maxpixid  = 0;
     Float_t maxpix  = 0;
-    Float_t maxdist = 0;
-
-    //    MSignalPix *pix = 0;
-    //    TIter Next(evt);
-    //    while ((pix=(MSignalPix*)Next()))
+
+    // Variables to caluclate time slope
+    // Formula: http://mo.mathematik.uni-stuttgart.de/inhalt/erlaeuterung/erlaeuterung300/
+    UInt_t cnt = 0;
+
+    Double_t sumx   = 0;
+    Double_t sumy   = 0;
+    Double_t sumt   = 0;
+    Double_t sumxy  = 0;
+    Double_t sumxt  = 0;
+    Double_t sumyt  = 0;
+    Double_t sumx2  = 0;
+    Double_t sumy2  = 0;
+    Double_t sumx3  = 0;
+    Double_t sumy3  = 0;
+    Double_t sumx2y = 0;
+    Double_t sumxy2 = 0;
 
     const UInt_t npix = evt.GetNumPixels();
@@ -149,21 +161,33 @@
             continue;
 
-        //const Int_t pixid = pix->GetPixId();
-
-        const MGeomPix &gpix = geom[i/*pixid*/];
-        const Double_t dx = gpix.GetX() - hil.GetMeanX();      // [mm]
-        const Double_t dy = gpix.GetY() - hil.GetMeanY();      // [mm]
-
-        Double_t nphot = pix.GetNumPhotons();                  // [1]
-
-        const Double_t dzx =  hil.GetCosDelta()*dx + hil.GetSinDelta()*dy; // [mm]
-        const Double_t dzy = -hil.GetSinDelta()*dx + hil.GetCosDelta()*dy; // [mm]
-
-        const Double_t dist = dx*dx+dy*dy;
-        if (dist>TMath::Abs(maxdist))
-            maxdist = dzx<0 ? -dist : dist;                    // [mm^2]
-
-        m3x += nphot * dzx*dzx*dzx;                            // [mm^3]
-        m3y += nphot * dzy*dzy*dzy;                            // [mm^3]
+        const MGeomPix &gpix = geom[i];
+
+        const Double_t x = gpix.GetX();
+        const Double_t y = gpix.GetY();
+        const Double_t t = pix.GetArrivalTime();
+
+        // --- time slope ----
+        sumx  += x;
+        sumy  += y;
+        sumt  += t;
+
+        sumx2 += x*x;
+        sumy2 += y*y;
+        sumxy += x*y;
+        sumxt += x*t;
+        sumyt += y*t;
+
+        // --- 3rd moment ---
+        const Double_t dx = x - hil.GetMeanX();      // [mm]
+        const Double_t dy = y - hil.GetMeanY();      // [mm]
+
+        Double_t nphot = pix.GetNumPhotons();        // [1]
+
+        sumx3  += nphot * dx*dx*dx;
+        sumy3  += nphot * dy*dy*dy;
+        sumx2y += nphot * dx*dx*dy;
+        sumx2y += nphot * dx*dy*dy;
+
+        cnt++;
 
         //
@@ -171,35 +195,49 @@
         // must take pixel size into account
         //
-        nphot *= geom.GetPixRatio(i/*pixid*/);
-
+        nphot *= geom.GetPixRatio(i);
+
+        // --- max pixel ---
         if (nphot>maxpix)
         {
-            maxpix   = nphot;                                  // [1]
-            maxpixid = i;//pixid;
+            maxpix   = nphot;                        // [1]
+            maxpixid = i;
         }
     }
 
+    const Double_t c = hil.GetCosDelta();
+    const Double_t s = hil.GetSinDelta();
+
+    //
+    // Time slope
+    //
+    const Double_t dxt  =  c*sumxt + s*sumyt;
+    const Double_t dyt  = -s*sumxt + c*sumyt;
+
+    const Double_t dx   =  c*sumx  + s*sumy;
+    const Double_t dy   = -s*sumx  + c*sumy;
+
+    const Double_t dx2  =  c*c*sumx2 + 2*c*s*sumxy + s*s*sumy2;
+    const Double_t dy2  =  s*s*sumx2 - 2*c*s*sumxy + c*c*sumy2;
+
+    const Double_t detx =  cnt*dx2 - dx*dx;
+    const Double_t dety =  cnt*dy2 - dy*dy;
+
+    fSlopeLong  = detx==0 ? 0 : (cnt*dxt - sumt*dx)/detx;
+    fSlopeTrans = dety==0 ? 0 : (cnt*dyt - sumt*dy)/dety;
+
+    //
+    // Third moments along axes get normalized
+    //
+    const Double_t m3l = c*c*c*sumx3 + 3*(s*c*c*sumx2y + c*s*s*sumxy2) + s*s*s*sumy3;
+    const Double_t m3t = c*c*c*sumy3 + 3*(s*s*c*sumx2y - s*c*c*sumxy2) - s*s*s*sumx3;
+
+    fM3Long  = MMath::Sqrt3(m3l/hil.GetSize());      // [mm]
+    fM3Trans = MMath::Sqrt3(m3t/hil.GetSize());      // [mm]
+
+    //
+    // Asymmetry
+    //
     const MGeomPix &maxp = geom[maxpixid];
-
-    //
-    // Asymmetry
-    //
-    fAsym = (hil.GetMeanX()-maxp.GetX())*hil.GetCosDelta() +
-            (hil.GetMeanY()-maxp.GetY())*hil.GetSinDelta();            // [mm]
-
-    //
-    // Third moments along axes get normalized
-    //
-    m3x /= hil.GetSize();
-    m3y /= hil.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]
-
-    //
-    // Distance between max signal and COG
-    //
-    const Double_t md = TMath::Sqrt(TMath::Abs(maxdist));
-    fMaxDist = maxdist<0 ? -md : md;                           // [mm]
+    fAsym = (hil.GetMeanX()-maxp.GetX())*c + (hil.GetMeanY()-maxp.GetY())*s;            // [mm]
 
     SetReadyToSave();
@@ -215,11 +253,12 @@
 void MHillasExt::Set(const TArrayF &arr)
 {
-    if (arr.GetSize() != 4)
+    if (arr.GetSize() != 5)
         return;
 
-    fAsym    = arr.At(0);
-    fM3Long  = arr.At(1);
-    fM3Trans = arr.At(2);
-    fMaxDist = arr.At(3);
+    fAsym       = arr.At(0);
+    fM3Long     = arr.At(1);
+    fM3Trans    = arr.At(2);
+    fSlopeLong  = arr.At(3);
+    fSlopeTrans = arr.At(4);
 }
 
@@ -228,12 +267,17 @@
 // Print contents of MHillasExt to *fLog.
 //
-void MHillasExt::Print(Option_t *) const
-{
+void MHillasExt::Print(Option_t *o) const
+{
+    const Bool_t showtrans = TString(o).Contains("trans");
+
     *fLog << all;
     *fLog << GetDescriptor() << endl;
     *fLog << " - Asymmetry      [mm]  = " << fAsym    << endl;
     *fLog << " - 3.Moment Long  [mm]  = " << fM3Long  << endl;
-    *fLog << " - 3.Moment Trans [mm]  = " << fM3Trans << endl;
-    *fLog << " - Max.Dist       [mm]  = " << fMaxDist << endl;
+    if (showtrans)
+        *fLog << " - 3.Moment Trans [mm]  = " << fM3Trans << endl;
+    *fLog << " - Slope Long     [mm]  = " << fSlopeLong  << endl;
+    if (showtrans)
+        *fLog << " - Slope Trans    [mm]  = " << fSlopeTrans << endl;
 }
 
@@ -249,6 +293,7 @@
     *fLog << " - Asymmetry      [deg] = " << fAsym   *geom.GetConvMm2Deg() << endl;
     *fLog << " - 3.Moment Long  [deg] = " << fM3Long *geom.GetConvMm2Deg() << endl;
-    *fLog << " - 3.Moment Trans [deg] = " << fM3Trans*geom.GetConvMm2Deg() << endl;
-    *fLog << " - Max.Dist       [deg] = " << fMaxDist*geom.GetConvMm2Deg() << endl;
+    *fLog << " - 3.Moment Trans [mm]  = " << fM3Trans*geom.GetConvMm2Deg() << endl;
+    *fLog << " - Slope Long     [mm]  = " << fSlopeLong*geom.GetConvMm2Deg()  << endl;
+    *fLog << " - Slope Trans    [mm]  = " << fSlopeTrans*geom.GetConvMm2Deg() << endl;
 }
 
@@ -272,5 +317,5 @@
 
     TMarker m;
-    m.SetMarkerColor(15);
+    m.SetMarkerColor(kRed);
     m.SetMarkerStyle(kFullDotLarge);
     m.PaintMarker(v.X(), v.Y());
