Index: trunk/MagicSoft/Mars/mhist/MHMatrix.cc
===================================================================
--- trunk/MagicSoft/Mars/mhist/MHMatrix.cc	(revision 1830)
+++ trunk/MagicSoft/Mars/mhist/MHMatrix.cc	(revision 1831)
@@ -591,9 +591,7 @@
 
     TMatrix m(n, fM.GetNcols());
+    TVector vold(fM.GetNcols());
     for (int i=0; i<n; i++)
-    {
-        TVector vold(fM.GetNcols());
         TMatrixRow(m, i) = vold = TMatrixRow(fM, idx[i]);
-    }
 
     fM = m;
@@ -666,4 +664,5 @@
     UInt_t newrow = 0;
 
+    TVector vold(fM.GetNcols());
     while (oldrow<rows)
     {
@@ -671,8 +670,5 @@
 
         if (newrow<=(unsigned int)sum)
-        {
-            TVector vold(fM.GetNcols());
             TMatrixRow(fM, newrow++) = vold = TMatrixRow(fM, oldrow);
-        }
 
         oldrow++;
@@ -743,10 +739,10 @@
     // Calculate normalization factors
     //
-    const int   nbins   = thsh.GetNbinsX();
-    const float frombin = thsh.GetBinLowEdge(1);
-    const float tobin   = thsh.GetBinLowEdge(nbins+1);
-    const float dbin    = thsh.GetBinWidth(1);
-    const Int_t nrows   = fM.GetNrows();
-    const Int_t ncols   = fM.GetNcols();
+    const int    nbins   = thsh.GetNbinsX();
+    const double frombin = thsh.GetBinLowEdge(1);
+    const double tobin   = thsh.GetBinLowEdge(nbins+1);
+    const double dbin    = thsh.GetBinWidth(1);
+    const Int_t  nrows   = fM.GetNrows();
+    const Int_t  ncols   = fM.GetNcols();
 
     //
@@ -757,6 +753,9 @@
         hth.Fill(fM(j, refcolumn-1));
 
+    hth.DrawCopy();
+
     TH1F hthd("thd", "correction factors", nbins, frombin, tobin);
     hthd.Divide((TH1F*)&thsh, &hth, 1, 1);
+
     if (hthd.GetMaximum() <= 0)
     {
@@ -793,19 +792,24 @@
     TArrayF cumulweight(nrows); // keep track for each bin how many events
 
-    TMatrix hindref(fM);
-    hindref -= frombin;
-    hindref *= 1/dbin;
+    //
+    // Project values in reference column into [0,1]
+    //
+    TVector v(fM.GetNrows());
+    v = TMatrixColumn(fM, refcolumn-1);
+    v += -frombin;
+    v *= 1/dbin;
 
     //
     // select events (distribution after renormalization)
     //
-    for (Int_t ir=0; ir<nmaxevts; ir++)
-    {
-        const Int_t indref = (Int_t)hindref(ind[ir], refcolumn);
+    Int_t ir;
+    TVector vold(fM.GetNcols());
+    for (ir=0; ir<nrows; ir++)
+    {
+        const Int_t indref = (Int_t)v(ind[ir]);
 
         cumulweight[indref] += hthd.GetBinContent(indref+1);
         if (cumulweight[indref]<=0.5)
         {
-            TVector vold(fM.GetNcols());
             TMatrixRow(mrest, evtcount2++) = vold = TMatrixRow(fM, ind[ir]);
             continue;
@@ -813,16 +817,12 @@
 
         cumulweight[indref] -= 1.;
-        if (++evtcount1 >= nmaxevts)
+        if (++evtcount1 >=  nmaxevts)
             break;
 
-        TVector vold(fM.GetNcols());
         TMatrixRow(mnewtmp, evtcount1) = vold = TMatrixRow(fM, ind[ir]);
     }
 
-    for (Int_t ir=nmaxevts; ir<nrows; ir++)
-    {
-        TVector vold(fM.GetNcols());
+    for (/*empty*/; ir<nrows; ir++)
         TMatrixRow(mrest, evtcount2++) = vold = TMatrixRow(fM, ind[ir]);
-    }
 
     //
@@ -835,9 +835,6 @@
     fM.ResizeTo(evtcount1, ncols);
     fNumRow = evtcount1;
-    for (Int_t ir=0; ir<evtcount1; ir++)
-    {
-        TVector vold(fM.GetNcols());
+    for (ir=0; ir<evtcount1; ir++)
         TMatrixRow(fM, ir) = vold = TMatrixRow(mnewtmp, ir);
-    }
 
     if (evtcount1 < nmaxevts)
@@ -848,5 +845,5 @@
 
     rest->ResizeTo(evtcount2, ncols);
-    for (Int_t ir=0; ir<evtcount1; ir++)
+    for (ir=0; ir<evtcount1; ir++)
     {
         TVector vold(fM.GetNcols());
@@ -855,59 +852,4 @@
 
     return kTRUE;
-
-    /*
-     TMatrix mnew(evtcount, nconl);
-     for (Int_t ir=0; ir<evtcount; ir++)
-     for (Int_t ic=0; ic<fNcols; ic++)
-     fM(ir,ic)   = mnewtmp(ir,ic);
-
-     //
-     //  test: print new matrix (part)
-     //
-     *fLog << "DefRefMatrix:  Event matrix (output) :" << endl;
-     *fLog << "DefRefMatrix:  Nrows, Ncols = " << mnew.GetNrows();
-     *fLog << " " << mnew.GetNcols() << endl;
-
-     for (Int_t ir=0;ir<10; ir++)
-     {
-     *fLog <<ir <<" ";
-     for (Int_t ic=0; ic<mnew.GetNcols(); ic++)
-     cout<<Mnew(ir,ic)<<" ";
-     *fLog <<endl;
-     }
-     */
-
-    /*
-     //  test print new bin contents
-     *fLog << "MHMatrix::DefRefMatrix;  new histogram: " << endl;
-     for (Int_t j=1; j<=fnbins; j++)
-     {
-     float a = fHthaft->GetBinContent(j);
-     *fLog << j << "  "<< a << "   ";
-     }
-     *fLog <<endl;
-     */
-
-    /*
-     //---------------------------------------------------------
-     // ==== plot four histograms
-     TCanvas *th1 = new TCanvas (fName, fName, 1);
-     th1->Divide(2,2);
-
-     th1->cd(1);
-     ((TH1F*)fHthsh)->DrawCopy();      // target
-
-     th1->cd(2);
-     ((TH1F*)fHth)->DrawCopy();        // real histogram before
-
-     th1->cd(3);
-     ((TH1F*)fHthd)->DrawCopy();       // correction factors
-
-     th1->cd(4);
-     ((TH1F*)fHthaft)->DrawCopy();     // histogram after
-
-     //---------------------------------------------------------
-     */
-    //return kTRUE;
 }
 
@@ -977,14 +919,10 @@
     // select events (distribution after renormalization)
     //
+    TVector vold(fM.GetNcols());
     for (Int_t ir=0; ir<nmaxevts; ir++)
-    {
-        TVector vold(fM.GetNcols());
         TMatrixRow(mnewtmp, evtcount1++) = vold = TMatrixRow(fM, ind[ir]);
-    }
+
     for (Int_t ir=nmaxevts; ir<nrows; ir++)
-    {
-        TVector vold(fM.GetNcols());
         TMatrixRow(mrest, evtcount2++) = vold = TMatrixRow(fM, ind[ir]);
-    }
 
     //