Changeset 1879 for trunk/MagicSoft/Mars/mhist/MHMatrix.cc

Timestamp:

03/28/03 14:24:38 (22 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mhist/MHMatrix.cc (modified) (14 diffs)

trunk/MagicSoft/Mars/mhist/MHMatrix.cc

@@ -250 +250 @@
         return kTRUE;
 
-    TMatrix m(fM);
-
-    fM.ResizeTo(fNumRow, fData->GetNumEntries());
-
-    TVector vold(fM.GetNcols());
-    for (int x=0; x<fM.GetNrows(); x++)
-        TMatrixRow(fM, x) = vold = TMatrixRow(m, x);
+    if (fNumRow != fM.GetNrows())
+    {
+        TMatrix m(fM);
+        CopyCrop(fM, m, fNumRow);
+    }
 
     return kTRUE;
@@ -687 +685 @@
 // ------------------------------------------------------------------------
 //
+//  Used in DefRefMatrix to display the result graphically
+//
+void MHMatrix::DrawDefRefInfo(const TH1 &hth, const TH1 &hthd, const TH1 &thsh, Int_t refcolumn)
+{
+    //
+    // Fill a histogram with the distribution after raduction
+    //
+    TH1F hta;
+    hta.SetName("hta");
+    hta.SetTitle("Distribution after reduction");
+    SetBinning(&hta, &hth);
+
+    for (Int_t i=0; i<fM.GetNrows(); i++)
+        hta.Fill(fM(i, refcolumn));
+
+    TCanvas *th1 = MakeDefCanvas(this);
+    th1->Divide(2,2);
+
+    th1->cd(1);
+    ((TH1&)hth).DrawCopy();   // real histogram before
+
+    th1->cd(2);
+    ((TH1&)hta).DrawCopy();   // histogram after
+
+    th1->cd(3);
+    ((TH1&)hthd).DrawCopy();  // correction factors
+
+    th1->cd(4);
+    ((TH1&)thsh).DrawCopy();  // target
+}
+
+// ------------------------------------------------------------------------
+//
+//  Resizes th etarget matrix to rows*source.GetNcol() and copies
+//  the data from the first (n)rows or the source into the target matrix.
+//
+void MHMatrix::CopyCrop(TMatrix &target, const TMatrix &source, Int_t rows)
+{
+    TVector v(source.GetNcols());
+
+    target.ResizeTo(rows, source.GetNcols());
+    for (Int_t ir=0; ir<rows; ir++)
+        TMatrixRow(target, ir) = v = TMatrixRow(source, ir);
+}
+
+// ------------------------------------------------------------------------
+//
 // Define the reference matrix
-//   refcolumn  number of the column (starting at 1)containing the variable, 
+//   refcolumn  number of the column (starting at 0) containing the variable,
 //              for which a target distribution may be given;
-//              if refcolumn is negative the target distribution will be set 
-//              equal to the real distribution; the events in the reference
-//              matrix will then be simply a random selection of the events 
-//              in the original matrix.
 //   thsh       histogram containing the target distribution of the variable
 //   nmaxevts   maximum number of events in the reference matrix
@@ -709 +750 @@
     }
 
-    if (refcolumn==0)
-    {
-        *fLog << err << dbginf << "Reference column 0 unknown." << endl;
-        return kFALSE;
-    }
-
     if (thsh.GetMinimum()<0)
     {
@@ -737 +772 @@
 
     //
-    // 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
     //
@@ -758 +791 @@
     // set up the real histogram (distribution before)
     //
-    TH1F hth("th", "data at input", nbins, frombin, tobin);
+    TH1F hth("th", "Distribution before reduction", nbins, frombin, tobin);
     for (Int_t j=0; j<nrows; j++)
-        hth.Fill(fM(j, refcolumn-1));
-
-    hth.DrawCopy();
-
-    TH1F hthd("thd", "correction factors", nbins, frombin, tobin);
+        hth.Fill(fM(j, refcolumn));
+
+    TH1F hthd("thd", "Correction factors", nbins, frombin, tobin);
     hthd.Divide((TH1F*)&thsh, &hth, 1, 1);
 
@@ -781 +812 @@
     // 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);
+    GetRandomArrayI(ind);
 
     //
@@ -805 +830 @@
     //
     TVector v(fM.GetNrows());
-    v = TMatrixColumn(fM, refcolumn-1);
+    v = TMatrixColumn(fM, refcolumn);
     v += -frombin;
     v *= 1/dbin;
@@ -842 +867 @@
     // this is the matrix to be used in the g/h separation
     //
-    fM.ResizeTo(evtcount1, ncols);
+    CopyCrop(fM, mnewtmp, evtcount1);
     fNumRow = evtcount1;
-    for (ir=0; ir<evtcount1; ir++)
-        TMatrixRow(fM, ir) = vold = TMatrixRow(mnewtmp, ir);
 
     if (evtcount1 < nmaxevts)
         *fLog << warn << "The reference sample contains less events (" << evtcount1 << ") than required (" << nmaxevts << ")" << endl;
 
+    if (TestBit(kEnableGraphicalOutput))
+        DrawDefRefInfo(hth, hthd, thsh, refcolumn);
+
     if (!rest)
         return kTRUE;
 
-    rest->ResizeTo(evtcount2, ncols);
-    for (ir=0; ir<evtcount1; ir++)
-    {
-        TVector vold(fM.GetNcols());
-        TMatrixRow(*rest, ir) = vold = TMatrixRow(mrest, ir);
-    }
+    CopyCrop(*rest, mrest, evtcount2);
 
     return kTRUE;
@@ -865 +886 @@
 // ------------------------------------------------------------------------
 //
+// Returns a array containing randomly sorted indices
+//
+void MHMatrix::GetRandomArrayI(TArrayI &ind) const
+{
+    const Int_t rows = ind.GetSize();
+
+    TArrayF ranx(rows);
+
+    TRandom3 rnd(0);
+    for (Int_t i=0; i<rows; i++)
+        ranx[i] = rnd.Rndm(i);
+
+    TMath::Sort(rows, ranx.GetArray(), ind.GetArray(), kTRUE);
+}
+
+// ------------------------------------------------------------------------
+//
 // Define the reference matrix
-//   refcolumn  number of the column (starting at 1)containing the variable, 
-//              for which a target distribution may be given;
-//              if refcolumn is negative the target distribution will be set 
-//              equal to the real distribution; the events in the reference
-//              matrix will then be simply a random selection of the events 
-//              in the original matrix.
-//   thsh       histogram containing the target distribution of the variable
 //   nmaxevts   maximum number of events in the reference matrix
 //   rest       a TMatrix conatining the resulting (not choosen)
@@ -878 +909 @@
 //              are not interested in this
 //
+//              the target distribution will be set
+//              equal to the real distribution; the events in the reference
+//              matrix will then be simply a random selection of the events 
+//              in the original matrix.
+//
 Bool_t MHMatrix::DefRefMatrix(Int_t nmaxevts, TMatrix *rest)
 {
@@ -907 +943 @@
     // 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);
+    GetRandomArrayI(ind);
 
     //
@@ -942 +972 @@
     // this is the matrix to be used in the g/h separation
     //
-    fM.ResizeTo(evtcount1, ncols);
+    CopyCrop(fM, mnewtmp, evtcount1);
     fNumRow = evtcount1;
-    for (Int_t ir=0; ir<evtcount1; ir++)
-    {
-        TVector vold(fM.GetNcols());
-        TMatrixRow(fM, ir) = vold = TMatrixRow(mnewtmp, ir);
-    }
 
     if (evtcount1 < nmaxevts)
@@ -956 +981 @@
         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);
-    }
+    CopyCrop(*rest, mrest, evtcount2);
 
     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;
 }
 
@@ -1028 +993 @@
 Int_t MHMatrix::Read(const char *name)
 {
-  Int_t ret = TObject::Read(name);
-  SetName(name);
-
-  return ret;
-}
-
-// --------------------------------------------------------------------------
+    Int_t ret = TObject::Read(name);
+    SetName(name);
+
+    return ret;
+}
+
+// --------------------------------------------------------------------------