Index: trunk/MagicSoft/Mars/mhist/MHMatrix.cc
===================================================================
--- trunk/MagicSoft/Mars/mhist/MHMatrix.cc	(revision 1581)
+++ trunk/MagicSoft/Mars/mhist/MHMatrix.cc	(revision 1582)
@@ -45,7 +45,4 @@
 #include "MHMatrix.h"
 
-#include <math.h>   // abs for Alphas
-#include <stdlib.h> // abs for some Linux machines
-
 #include <fstream.h>
 
@@ -384,7 +381,7 @@
 // represented by the covariance metrix m.
 //  - If n<0 the kernel method is applied and
-//    sum(epx(-d)/n is returned.
-//  - For n>=0 the n nearest neighbors are summed and
-//    sqrt(sum(d))/n is returned.
+//    -log(sum(epx(-d/h))/n) is returned.
+//  - For n>0 the n nearest neighbors are summed and
+//    sqrt(sum(d)/n) is returned.
 //  - if n==0 all distances are summed
 //
@@ -395,5 +392,5 @@
         TVector d = evt;
         d *= m;
-        return sqrt(d*evt); 
+        return TMath::Sqrt(d*evt);
     }
 
@@ -430,18 +427,19 @@
     TMath::Sort(dists.GetSize(), dists.GetArray(), idx.GetArray(), kFALSE);
 
-    const Int_t n = abs(num)<rows ? abs(num) : rows;
+    const Int_t n = TMath::Abs(num)<rows ? TMath::Abs(num) : rows;
 
     if (num<0)
     {
         //
-        // Kernel function sum
+        // Kernel function sum (window size h set according to literature)
         //
-        const Double_t h = 1;
+        const Double_t h    = TMath::Power(rows, -1./(cols+4));
+        const Double_t hwin = h*h*2;
 
         Double_t res = 0;
         for (int i=0; i<n; i++)
-            res += exp(-dists[idx[i]]/h);
-
-        return -log(res/n);
+            res += TMath::Exp(-dists[idx[i]]/hwin);
+
+        return -TMath::Log(res/n);
     }
     else
@@ -454,5 +452,5 @@
             res += dists[idx[i]];
 
-        return sqrt(res/n);
+        return TMath::Sqrt(res/n);
     }
 }