Index: trunk/MagicSoft/Mars/Changelog
===================================================================
--- trunk/MagicSoft/Mars/Changelog	(revision 1825)
+++ trunk/MagicSoft/Mars/Changelog	(revision 1826)
@@ -5,4 +5,8 @@
     * mfileio/MReadTree.cc:
       - fixed a bug in the AddFile function
+
+    * mhist/MHMatrix.[h,cc]:
+      - implemented a request of Th. Hengstebeck: Let DefRefMatrix
+        return the 'unused' events
 
 
Index: trunk/MagicSoft/Mars/mhist/MHMatrix.cc
===================================================================
--- trunk/MagicSoft/Mars/mhist/MHMatrix.cc	(revision 1825)
+++ trunk/MagicSoft/Mars/mhist/MHMatrix.cc	(revision 1826)
@@ -206,6 +206,8 @@
 
     for (int x=0; x<m.GetNcols(); x++)
-        for (int y=0; y<m.GetNrows(); y++)
-            fM(y, x) = m(y, x);
+    {
+        TVector vold(fM.GetNcols());
+        TMatrixColumn(fM, x) = vold = TMatrixColumn(m, x);
+    }
 }
 
@@ -242,6 +244,8 @@
 
     for (int x=0; x<fM.GetNcols(); x++)
-        for (int y=0; y<fM.GetNrows(); y++)
-            fM(y, x) = m(y, x);
+    {
+        TVector vold(fM.GetNcols());
+        TMatrixColumn(fM, x) = vold = TMatrixColumn(m, x);
+    }
 
     return kTRUE;
@@ -366,6 +370,5 @@
         avg /= rows;
 
-        for (int y=0; y<rows; y++)
-            m(y, x) -= avg;
+        TMatrixColumn(m, x) += -avg;
     }
 
@@ -584,8 +587,5 @@
     {
         TVector vold(n);
-        vold = TMatrixRow(fM, idx[i]);
-
-        TMatrixRow rownew(m, i);
-        rownew = vold;
+        TMatrixRow(m, i) = vold = TMatrixRow(fM, idx[i]);
     }
 
@@ -666,26 +666,9 @@
         {
             TVector vold(fM.GetNcols());
-            vold = TMatrixRow(fM, oldrow);
-
-            TMatrixRow rownew(fM, newrow);
-            rownew = vold;
-            newrow++;
+            TMatrixRow(fM, newrow++) = vold = TMatrixRow(fM, oldrow);
         }
 
         oldrow++;
     }
-
-/*
-
-    TMatrix m(n, fM.GetNcols());
-    for (int i=0; i<n; i++)
-    {
-        TVector vold(n);
-        vold = TMatrixRow(fM, idx[i]);
-
-        TMatrixRow rownew(m, i);
-        rownew = vold;
-    }
-*/
 }
 
@@ -701,276 +684,225 @@
 //   thsh       histogram containing the target distribution of the variable
 //   nmaxevts   maximum number of events in the reference matrix
-//
-Bool_t MHMatrix::DefRefMatrix(const Int_t refcolumn, const TH1F *thsh, 
-                              const Int_t nmaxevts)
-{
-   if (!fM.IsValid()) 
-   {
-     *fLog << "MHMatrix::DefRefMatrix; matrix not initialized" << endl;
-     return kFALSE;
-   }
-
-    const TH1F *fHthsh;    // target distribution
-    TH1F *fHthd;     // normalization factors
-    TH1F *fHth;      // distribution before
-    TH1F *fHthaft;   // distribution after renormalization 
-
-    Int_t frefcol;    // number of the column (count from 0) containing 
-                      // the variable for which the target distribution is given
-    Int_t   fnbins;   // histogram parameters for thsh
-    Float_t ffrombin; //
-    Float_t ftobin;   //
-    Float_t fdbin;    //
-
-    TArrayF fnormfac; // array of normalization factors
-    TMatrix fMnew;    // matrix fM having the requested distribution
-                      // and the requested number of rows;
-                      // this is the matrix to be used in the g/h separation
-
-
-   //===================================================================
-   //---------------------------------------------------------
-   // Calculate normalization factors
-   //
-
-   fHthsh = thsh;
-
-   Int_t   fNrows  = fM.GetNrows();
-   Int_t   fNcols  = fM.GetNcols();
-
-   // print characteristics
-   //
-   //Int_t   fColLwb = fM.GetColLwb();
-   //Int_t   fRowLwb = fM.GetRowLwb();
-   //*fLog << "MHMatrix::CalcNormFactors; Event matrix (input): " << endl;
-   //*fLog << "          fNrows, fNcols, fColLwb, fRowLwb = " << fNrows << ",  "
-   //      << fNcols << ",  " << fColLwb << ",  " << fRowLwb << endl;
-
-   
-   //---------------------------------------------------------
-   //
-   // if refcolumn < 0 : select reference events randomly
-   //                    i.e. set the normaliztion factotrs equal to 1.0
-   // refcol is the column number starting at 0; it is >= 0 
-
-   if (refcolumn<0) 
-   {
-     frefcol = -refcolumn - 1;
-   }
-   else
-   {
-     frefcol =  refcolumn - 1;
-   }
-   
-
-   //---------------------------------------------------------
-
-   if (fHthsh == NULL)
-   {
-     *fLog << "MHMatrix::DefRefMatrix; target distribution has to be defined"
-           << endl;
-     return kFALSE;
-   }
-
-   fnbins   = fHthsh->GetNbinsX();
-   ffrombin = fHthsh->GetBinLowEdge(1);
-   ftobin   = fHthsh->GetBinLowEdge(fnbins+1);
-   fdbin    = (ftobin-ffrombin)/fnbins;
-   //*fLog << "CalcNormFactors; Reference column (count from 0) = "
-   //      << frefcol << endl;
-   //*fLog << "CalcNormFactors; Histo params : " << fnbins << ",  " << ffrombin
-   //      <<",  "<< ftobin << endl;
-
-
-   //---------------------------------------------------------
-   // ===== set up the real histogram
-   //
-   fHth = new TH1F("th", "data at input", fnbins, ffrombin, ftobin);
-   for (Int_t j=1; j<=fNrows; j++)  
-   {
-     fHth->Fill(fM(j-1,frefcol));   
-     //*fLog << "j, frefcol, fM(j-1,frefcol) = " << j << ",  " << frefcol
-     //      << ",  " << fM(j-1,frefcol) << endl;
-   }
-
-   //---------------------------------------------------------
-   // ===== obtain correction factors
-   //
-   fHthd = new TH1F ("thd", "correction factors", fnbins, ffrombin, ftobin);
-   
-   Double_t cmax = 0.;
-   Float_t a;
-   Float_t b;
-   Float_t c;
-   for (Int_t j=1; j<=fnbins; j++) 
-   {
-     a = fHth->GetBinContent(j);
-
-     // if refcolumn < 0 set the correction factors equal to 1
-     if ( refcolumn>=0 )
-       b = fHthsh->GetBinContent(j);
-     else
-       b = a;
-
-     if ( a<0.0  || b<0.0 )
+//   rest       a TMatrix conatining the resulting (not choosen)
+//              columns of the primary matrix. Maybe NULL if you
+//              are not interested in this
+//
+Bool_t MHMatrix::DefRefMatrix(const TH1F &thsh, const Int_t refcolumn,
+                              const Int_t nmaxevts, TMatrix *rest)
+{
+    if (!fM.IsValid())
+    {
+        *fLog << err << dbginf << "Matrix not initialized" << endl;
+        return kFALSE;
+    }
+
+    //
+    // if refcolumn < 0 : select reference events randomly
+    //                    i.e. set the normaliztion factotrs equal to 1.0
+    // refcol is the column number starting at 0; it is >= 0
+    //
+    // number of the column (count from 0) containing
+    // the variable for which the target distribution is given
+    //
+    const Int_t refcol = refcolumn<0 ? -refcolumn - 1 : refcolumn - 1;
+
+    //
+    // 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();
+
+    //
+    // set up the real histogram (distribution before)
+    //
+    TH1F hth("th", "data at input", nbins, frombin, tobin);
+    for (Int_t j=0; j<nrows; j++)
+        hth.Fill(fM(j, refcol));
+
+    //
+    // ===== obtain correction factors (normalization factors)
+    //
+    TH1F hthd("thd", "correction factors", nbins, frombin, tobin);
+
+    Double_t cmax = 0.;
+    for (Int_t j=1; j<=nbins; j++)
+    {
+        const float a = hth.GetBinContent(j);
+
+        // if refcolumn < 0 set the correction factors equal to 1
+        const float b = refcolumn>=0 ? thsh.GetBinContent(j) : a;
+
+        if (a<0 || b<0)
+        {
+            *fLog << err << dbginf << "renormalization of input matrix is not possible"   << endl;
+            *fLog << "          a, b = " << a << ",  " << b << endl;
+            return kFALSE;
+        }
+
+        const float c = a==0 ? 0 : b/a;
+
+        hthd.SetBinContent(j, c);
+
+        if (cmax < c)
+            cmax = c;
+    }
+
+    if (hthd.GetMaximum() <= 0)
+    {
+        *fLog << err << dbginf << "Maximum ratio is LE zero" << endl;
+        return kFALSE;
+    }
+
+    hthd.Scale(1/hthd.GetMaximum());
+
+    //
+    // ==== compress matrix fM into Mnew
+    //
+
+    //
+    // get random access
+    //
+    TArrayF ranx(nrows);
+
+    TRandom3 rnd(0);
+    for (Int_t i=0; i<nrows; i++)
+        ranx[i] = rnd.Rndm(i);
+
+    TArrayI ind(nrows);
+    TMath::Sort(nrows, ranx.GetArray(), ind.GetArray(), kTRUE);
+
+    //
+    // define  new matrix
+    //
+    Int_t evtcount1 = -1;
+    Int_t evtcount2 =  0;
+    TMatrix mnewtmp(nrows, ncols);
+    TMatrix mrest(nrows, ncols);
+
+    TArrayF cumulweight(nrows); // keep track for each bin how many events
+
+    TMatrix hindref(fM);
+    hindref -= frombin;
+    hindref *= 1/dbin;
+
+    //
+    // select events (distribution after renormalization)
+    //
+    for (Int_t ir=0; ir<nrows; ir++)
+    {
+        const Int_t indref = (Int_t)hindref(ind[ir], refcol);
+
+        cumulweight[indref] += hthd.GetBinContent(indref+1);
+        if (cumulweight[indref]<=0.5)
+        {
+            TVector vold(fM.GetNcols());
+            TMatrixRow(mrest, evtcount2++) = vold = TMatrixRow(fM, ind[ir]);
+            continue;
+        }
+
+        cumulweight[indref] -= 1.;
+        if (++evtcount1 >= nmaxevts)
+        {
+            for (; ir<nrows; ir++)
+            {
+                TVector vold(fM.GetNcols());
+                TMatrixRow(mrest, evtcount2++) = vold = TMatrixRow(fM, ind[ir]);
+            }
+            break;
+        }
+
+        TVector vold(fM.GetNcols());
+        TMatrixRow(mnewtmp, evtcount1) = vold = TMatrixRow(fM, ind[ir]);
+    }
+
+    //
+    // reduce size
+    //
+    // matrix fM having the requested distribution
+    // and the requested number of rows;
+    // this is the matrix to be used in the g/h separation
+    //
+
+    fM.ResizeTo(evtcount1, ncols);
+    fNumRow = evtcount1;
+    for (Int_t ir=0; ir<evtcount1; ir++)
+    {
+        TVector vold(fM.GetNcols());
+        TMatrixRow(fM, ir) = vold = TMatrixRow(mnewtmp, ir);
+    }
+
+    if (evtcount1 < nmaxevts)
+        *fLog << warn << "The reference sample contains less events (" << evtcount1 << ") than required (" << nmaxevts << ")" << endl;
+
+    if (!rest)
+        return kTRUE;
+
+    rest->ResizeTo(evtcount2, ncols);
+    for (Int_t ir=0; ir<evtcount1; ir++)
+    {
+        TVector vold(fM.GetNcols());
+        TMatrixRow(*rest, ir) = vold = TMatrixRow(mrest, ir);
+    }
+
+    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 << "MHMatrix::DefRefMatrix; renormalization of input matrix is not possible"   << endl;
-       *fLog << "          a, b = " << a << ",  " << b << endl;
-        return kFALSE;
+     *fLog <<ir <<" ";
+     for (Int_t ic=0; ic<mnew.GetNcols(); ic++)
+     cout<<Mnew(ir,ic)<<" ";
+     *fLog <<endl;
      }
-
-     if (a == 0.0)
-       c = 0.0;
-     else
-       c = b/a;
-
-     fHthd->SetBinContent(j, c);
-     if (cmax < c) cmax = c;
-     //*fLog << j << ",  " << a << ",  " << b << ",  " << c << endl;
-   }
-
-   //*fLog << "MHMatrix::CalcNormFactors; maximum ratio between histogram and target "
-   //      << cmax << endl;
-
-   if (cmax <= 0.0)
-   {
-     *fLog << "MHMatrix::CalcNormFactors; maximum ratio is LE zero" << endl;
-     return kFALSE;
-   }
-
-   fnormfac.Set(fnbins);
-   Float_t minbin = 1.0;
-   for (Int_t j=1; j<=fnbins; j++) 
-   {
-     float c = (fHthd->GetBinContent(j)) / cmax;
-     fnormfac[j-1] = c;
-     if (minbin > c  &&  c>0) minbin = c;
-     fHthd->SetBinContent(j, c);
-   }
-   //*fLog << "MHMatrix::CalcNormFactors; maximum weight = 1, minimum non-zero weight = "
-   //      << minbin <<endl;
-
-
-   //---------------------------------------------------------
-   //===================================================================
-
-
-   // ==== compress matrix fM into Mnew
-   //
-   // get random access
-   TArrayF ranx(fNrows);
-   TRandom3 *fRnd = new TRandom3(0);
-
-   for (Int_t i=0;i<fNrows;i++) ranx[i] = fRnd->Rndm(i);
-
-   TArrayI ind(fNrows);
-   TMath::Sort(fNrows, ranx.GetArray(), ind.GetArray(),kTRUE);
-
-   //*fLog << "MHMatrix::DefRefMatrix;  permuting array " <<endl;
-   //for (Int_t i=0; i<10; i++)  cout << ind[i] << " ";
-   //*fLog << endl;
-
-   //---------------------------------------------------------
-   // go through matrix and keep events depending on content of 
-   // reference column
-   //
-   //*fLog << "MHMatrix::DefRefMatrix;  pass at most " << nmaxevts 
-   //      << " events, in random order " << endl;
-
-
-   //---------------------------------------------------------
-   // define  new matrix
-   Int_t evtcount = 0;  
-   TMatrix Mnew_tmp (fNrows,fNcols);
-
-   TArrayF cumul_weight(fNrows);   // keep track for each bin how many events
-
-   for (Int_t ir=0; ir<fNrows; ir++) cumul_weight[ir]=0;
-
-   //--------------------------------
-   // select events
-   //
-   fHthaft = new TH1F ("thaft","data after selection", fnbins, ffrombin, ftobin);
-   for (Int_t ir=0; ir<fNrows; ir++)   
-   {
-     Int_t indref = (Int_t) ((fM(ind[ir],frefcol)-ffrombin)/fdbin);
-     cumul_weight[indref] += fnormfac[indref];
-     if (cumul_weight[indref]>0.5)    
+     */
+
+    /*
+     //  test print new bin contents
+     *fLog << "MHMatrix::DefRefMatrix;  new histogram: " << endl;
+     for (Int_t j=1; j<=fnbins; j++)
      {
-       cumul_weight[indref] -= 1.;
-       if (evtcount++ <  nmaxevts) 
-       {
-	 fHthaft->Fill(fM(ind[ir],frefcol)); 
-	 for (Int_t ic=0; ic<fNcols; ic++)  
-	   Mnew_tmp(evtcount-1,ic) = fM(ind[ir],ic);
-       }
-       else  break;
-     }
-   }
-   evtcount--;
-   //*fLog << "MHMatrix::DefRefMatrix;  passed       " << evtcount 
-   //      << " events" << endl;
-
-   //--------------------------------
-   // reduce size
-
-   TMatrix Mnew(evtcount, fNcols);
-   fM.ResizeTo (evtcount, fNcols);
-   fNumRow = evtcount;
-
-   for (Int_t ir=0;ir<evtcount; ir++) 
-      for (Int_t ic=0;ic<fNcols;ic++)  
-      {
-        Mnew(ir,ic) = Mnew_tmp(ir,ic);
-        fM(ir,ic)   = Mnew_tmp(ir,ic);
-      }
-
-   if (evtcount < nmaxevts)
-   {
-     *fLog << "MHMatrix::DefRefMatrix; Warning : The reference sample contains less events (" << evtcount << ") than required (" << nmaxevts << ")" << endl;
-   }
-
-
-   //---------------------------------------------------------
-   //  test: print new matrix (part)
-   *fLog << "MHMatrix::DefRefMatrix;  Event matrix (output) :" << endl;
-   *fLog << "MHMatrix::DefRefMatrix;  Nrows, Ncols = " << Mnew.GetNrows()
-         << " " << 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;
+     }
+     *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;
 }
 
@@ -989,9 +921,2 @@
 
 // --------------------------------------------------------------------------
-
-
-
-
-
-
-
Index: trunk/MagicSoft/Mars/mhist/MHMatrix.h
===================================================================
--- trunk/MagicSoft/Mars/mhist/MHMatrix.h	(revision 1825)
+++ trunk/MagicSoft/Mars/mhist/MHMatrix.h	(revision 1826)
@@ -99,6 +99,6 @@
     Int_t Read(const char *name);
 
-    Bool_t DefRefMatrix   (const Int_t refcolumn=-1, const TH1F *thsh=NULL, 
-                           const Int_t nmaxevts=0);
+    Bool_t DefRefMatrix(const TH1F &thsh, const Int_t refcolumn=-1,
+                        const Int_t nmaxevts=0, TMatrix *mrest=NULL);
 
     ClassDef(MHMatrix, 1) // Multidimensional Matrix to store events