Index: trunk/MagicSoft/Mars/mranforest/MRanForestCalc.cc
===================================================================
--- trunk/MagicSoft/Mars/mranforest/MRanForestCalc.cc	(revision 7697)
+++ trunk/MagicSoft/Mars/mranforest/MRanForestCalc.cc	(revision 7710)
@@ -1,3 +1,3 @@
-/* ======================================================================== *\
+/* ================================q======================================== *\
 !
 ! *
@@ -32,4 +32,6 @@
 #include "MRanForestCalc.h"
 
+#include <TF1.h>
+#include <TGraph.h>
 #include <TVector.h>
 
@@ -279,6 +281,5 @@
     }
 
-    // Maybe fEForests[0].fRules yould be used instead?
-
+    // Maybe fEForests[0].fRules could be used instead?
     if (fData->Read("rules")<=0)
     {
@@ -322,6 +323,4 @@
 }
 
-#include <TGraph.h>
-#include <TF1.h>
 Int_t MRanForestCalc::Process()
 {
Index: trunk/MagicSoft/Mars/mranforest/MRanTree.cc
===================================================================
--- trunk/MagicSoft/Mars/mranforest/MRanTree.cc	(revision 7697)
+++ trunk/MagicSoft/Mars/mranforest/MRanTree.cc	(revision 7710)
@@ -41,4 +41,6 @@
 #include "MArrayF.h"
 
+#include "MMath.h"
+
 #include "MLog.h"
 #include "MLogManip.h"
@@ -176,5 +178,6 @@
     Double_t pdo=0;
 
-    for (Int_t j=0; j<nclass; j++)
+    // tclasspop: sum of weights for events in class
+    for (Int_t j=0; j<nclass; j++) // loop over number of classes to classifiy
     {
         pno+=tclasspop[j]*tclasspop[j];
@@ -182,5 +185,5 @@
     }
 
-    const Double_t crit0=pno/pdo;
+    const Double_t crit0=pno/pdo;  // weighted mean of weights
 
     // start main loop through variables to find best split,
@@ -190,5 +193,5 @@
 
     // random split selection, number of trials = fNumTry
-    for (Int_t mt=0; mt<fNumTry; mt++)
+    for (Int_t mt=0; mt<fNumTry; mt++) // we could try ALL variables???
     {
         const Int_t mvar= gRandom->Integer(mdim);
@@ -212,23 +215,46 @@
 
             // do classification, Gini index as split rule
-            rln+=u*( 2*wl[k]+u);
-            rrn+=u*(-2*wr[k]+u);
-
-            rld+=u;
-            rrd-=u;
-
-            wl[k]+=u;
-            wr[k]-=u;
-
+            rln   +=u*(2*wl[k]+u);  // += u*(wl[k]{i-1} + wl[k]{i-1}+u{i})
+            rld   +=u;   // sum of weights left  from cut total
+            wl[k] +=u;   // sum of weights left  from cut for class k
+
+            rrn   -=u*(2*wr[k]-u);  // -= u*(wr[k]{i-1} + wr[k]{i-1}-u{i})
+            //  rr0=0; rr0+=u*2*tclasspop[k]
+            //  rrn = pno - rr0 + rln
+            rrd   -=u;   // sum of weights right from cut total
+            wr[k] -=u;   // sum of weights right from cut for class k
+
+            // REPLACE BY?
+            // rr0   = 0
+            // rr0  += u*2*tclasspop[k]
+            // rrn   = pno - rr0 + rln
+            // rrd   = pdo - rld
+            // wr[k] = tclasspop[k] - wl[k]
+
+            // crit = (rln*(pdo - rld + 1) + pno - rr0) / rld*(pdo - rld)
+
+            /*
+             if (k==background)
+                continue;
+             crit = TMath::Max(MMath::SignificanceLiMa(rld, rld-wl[k]),
+                               MMath::SignificanceLiMa(rrd, rrd-wr[k]))
+             */
+
+            // This condition is in fact a == (> cannot happen at all)
+            // This is because we cannot set the cut between two identical values
+            //if (datarang[mn+datasort[mn+nsp]]>=datarang[mn+datasort[mn+nsp+1]])
             if (datarang[mn+nc]>=datarang[mn+datasort[mn+nsp+1]])
                 continue;
 
-            if (TMath::Min(rrd,rld)<=1.0e-5)
+            // If crit starts to become pretty large do WHAT???
+            if (TMath::Min(rrd,rld)<=1.0e-5) // FIXME: CHECKIT FOR WEIGHTS!
                 continue;
 
             const Double_t crit=(rln/rld)+(rrn/rrd);
 
+            // Search for the highest value of crit
             if (crit<=critvar) continue;
 
+            // store the highest crit value and the corresponding event to cut at
             nbestvar=nsp;
             critvar=crit;
@@ -237,9 +263,12 @@
         if (critvar<=critmax) continue;
 
-        msplit=mvar;
-        nbest=nbestvar;
+        msplit=mvar;      // Variable in which to split
+        nbest=nbestvar;   // event at which the best split was found
         critmax=critvar;
     }
 
+    // crit0 = MMath::SignificanceLiMa(pdo, pdo-tclasspop[0])
+    // mean increase of sensitivity
+    // decsplit = sqrt(critmax/crit0)
     decsplit=critmax-crit0;
 
@@ -299,6 +328,4 @@
         const Int_t mn  = mvar*numdata;
 
-        Double_t rrn=0, rrd=0, rln=0, rld=0;
-
         Double_t esumr =mean;
         Double_t e2sumr=square;
@@ -316,30 +343,4 @@
             const Float_t &u=winbag[nc];
 
-            e2sumr-=u*f*f;
-            esumr -=u*f;
-            wr    -=u;
-
-            //-------------------------------------------
-            // resolution
-            //rrn=(wr*e2sumr-esumr*esumr)*wr;
-            //rrd=(wr-1.)*esumr*esumr;
-
-            // resolution times n
-            //rrn=(wr*e2sumr-esumr*esumr)*wr;
-            //rrd=esumr*esumr;
-
-            // sigma
-            //rrn=(e2sumr-esumr*esumr/wr);
-            //rrd=(wr-1.);
-
-            // sigma times n
-            rrn=(e2sumr-esumr*esumr/wr);
-            rrd=1.;
-
-            // 1./(n*variance)
-            //rrn=1.;
-            //rrd=(e2sumr-esumr*esumr/wr);
-            //-------------------------------------------
-
             e2suml+=u*f*f;
             esuml +=u*f;
@@ -348,22 +349,49 @@
             //-------------------------------------------
             // resolution
-            //rln=(wl*e2suml-esuml*esuml)*wl;
-            //rld=(wl-1.)*esuml*esuml;
+            //const Double_t rln=(wl*e2suml-esuml*esuml)*wl;
+            //const Double_t rld=(wl-1.)*esuml*esuml;
 
             // resolution times n
-            //rln=(wl*e2suml-esuml*esuml)*wl;
-            //rld=esuml*esuml;
+            //const Double_t rln=(wl*e2suml-esuml*esuml)*wl;
+            //const Double_t rld=esuml*esuml;
 
             // sigma
-            //rln=(e2suml-esuml*esuml/wl);
-            //rld=(wl-1.);
+            //const Double_t rln=(e2suml-esuml*esuml/wl);
+            //const Double_t rld=(wl-1.);
 
             // sigma times n
-            rln=(e2suml-esuml*esuml/wl);
-            rld=1.;
+            Double_t rln=(e2suml-esuml*esuml/wl);
+            Double_t rld=1.;
 
             // 1./(n*variance)
-            //rln=1.;
-            //rld=(e2suml-esuml*esuml/wl);
+            //const Double_t rln=1.;
+            //const Double_t rld=(e2suml-esuml*esuml/wl);
+            //-------------------------------------------
+
+            // REPLACE BY??? 
+            e2sumr-=u*f*f;   // e2sumr = square       - e2suml
+            esumr -=u*f;     // esumr  = mean         - esuml
+            wr    -=u;       // wr     = tclasspop[0] - wl
+
+            //-------------------------------------------
+            // resolution
+            //const Double_t rrn=(wr*e2sumr-esumr*esumr)*wr;
+            //const Double_t rrd=(wr-1.)*esumr*esumr;
+
+            // resolution times n
+            //const Double_t rrn=(wr*e2sumr-esumr*esumr)*wr;
+            //const Double_t rrd=esumr*esumr;
+
+            // sigma
+            //const Double_t rrn=(e2sumr-esumr*esumr/wr);
+            //const Double_t rrd=(wr-1.);
+
+            // sigma times n
+            const Double_t rrn=(e2sumr-esumr*esumr/wr);
+            const Double_t rrd=1.;
+
+            // 1./(n*variance)
+            //const Double_t rrn=1.;
+            //const Double_t rrd=(e2sumr-esumr*esumr/wr);
             //-------------------------------------------
 
@@ -538,7 +566,9 @@
               const int j=idclass[nc];
                    
+              // statistics left from cut
               mean[ncur+1]+=hadtrue[nc]*winbag[nc];
               square[ncur+1]+=hadtrue[nc]*hadtrue[nc]*winbag[nc];
 
+              // sum of weights left from cut
               classpop[j*nrnodes+ncur+1]+=winbag[nc];
           }
@@ -549,7 +579,9 @@
               const int j=idclass[nc];
 
+              // statistics right from cut
               mean[ncur+2]  +=hadtrue[nc]*winbag[nc];
               square[ncur+2]+=hadtrue[nc]*hadtrue[nc]*winbag[nc];
 
+              // sum of weights right from cut
               classpop[j*nrnodes+ncur+2]+=winbag[nc];
           }
