Changeset 1828

Timestamp:

03/18/03 18:15:05 (22 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mhist

Files:

• MHMatrix.cc (modified) (10 diffs)
• MHMatrix.h (modified) (1 diff)

trunk/MagicSoft/Mars/mhist/MHMatrix.cc

@@ -688 +688 @@
 //              are not interested in this
 //
-Bool_t MHMatrix::DefRefMatrix(const TH1F &thsh, const Int_t refcolumn,
+Bool_t MHMatrix::DefRefMatrix(const UInt_t refcolumn, const TH1F &thsh,
                               const Int_t nmaxevts, TMatrix *rest)
 {
@@ -694 +694 @@
     {
         *fLog << err << dbginf << "Matrix not initialized" << endl;
+        return kFALSE;
+    }
+
+    if (refcolumn==0)
+    {
+        *fLog << err << dbginf << "Reference column 0 unknown." << endl;
+        return kFALSE;
+    }
+
+    if (thsh.GetMinimum()<0)
+    {
+        *fLog << err << dbginf << "Renormalization not possible: Target Distribution has values < 0" << endl;
+        return kFALSE;
+    }
+
+    if (nmaxevts>fM.GetNrows())
+    {
+        *fLog << err << dbginf << "Number of maximum events exceeds number of events" << endl;
+        return kFALSE;
+    }
+
+    if (nmaxevts<0)
+    {
+        *fLog << err << dbginf << "Number of maximum events < 0" << endl;
         return kFALSE;
     }
@@ -705 +729 @@
     // the variable for which the target distribution is given
     //
-    const Int_t refcol = refcolumn<0 ? -refcolumn - 1 : refcolumn - 1;
 
     //
@@ -722 +745 @@
     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)
-    //
+        hth.Fill(fM(j, refcolumn-1));
+
     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;
-    }
-
+    hthd.Divide((TH1F*)&thsh, &hth, 1, 1);
     if (hthd.GetMaximum() <= 0)
     {
@@ -758 +755 @@
     }
 
+    //
+    // ===== obtain correction factors (normalization factors)
+    //
     hthd.Scale(1/hthd.GetMaximum());
-
-    //
-    // ==== compress matrix fM into Mnew
-    //
 
     //
@@ -781 +777 @@
     Int_t evtcount1 = -1;
     Int_t evtcount2 =  0;
+
     TMatrix mnewtmp(nrows, ncols);
     TMatrix mrest(nrows, ncols);
@@ -793 +790 @@
     // select events (distribution after renormalization)
     //
-    for (Int_t ir=0; ir<nrows; ir++)
-    {
-        const Int_t indref = (Int_t)hindref(ind[ir], refcol);
+    for (Int_t ir=0; ir<nmaxevts; ir++)
+    {
+        const Int_t indref = (Int_t)hindref(ind[ir], refcolumn);
 
         cumulweight[indref] += hthd.GetBinContent(indref+1);
@@ -807 +804 @@
         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]);
+    }
+
+    for (Int_t ir=nmaxevts; ir<nrows; ir++)
+    {
+        TVector vold(fM.GetNcols());
+        TMatrixRow(mrest, evtcount2++) = vold = TMatrixRow(fM, ind[ir]);
     }
 
@@ -827 +823 @@
     // this is the matrix to be used in the g/h separation
     //
-
     fM.ResizeTo(evtcount1, ncols);
     fNumRow = evtcount1;
@@ -907 +902 @@
 }
 
+// ------------------------------------------------------------------------
+//
+// 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)
+//              columns of the primary matrix. Maybe NULL if you
+//              are not interested in this
+//
+Bool_t MHMatrix::DefRefMatrix(const Int_t nmaxevts, TMatrix *rest)
+{
+    if (!fM.IsValid())
+    {
+        *fLog << err << dbginf << "Matrix not initialized" << endl;
+        return kFALSE;
+    }
+
+    const Int_t nrows = fM.GetNrows();
+    const Int_t ncols = fM.GetNcols();
+
+    //
+    // 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 = 0;
+    Int_t evtcount2 = 0;
+
+    TMatrix mnewtmp(nrows, ncols);
+    TMatrix mrest(nrows, ncols);
+
+    //
+    // select events (distribution after renormalization)
+    //
+    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]);
+    }
+
+    //
+    // 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 <<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;
+}
+
 // --------------------------------------------------------------------------
 //

trunk/MagicSoft/Mars/mhist/MHMatrix.h

@@ -99 +99 @@
     Int_t Read(const char *name);
 
-    Bool_t DefRefMatrix(const TH1F &thsh, const Int_t refcolumn=-1,
+    Bool_t DefRefMatrix(const UInt_t refcolumn, const TH1F &thsh,
                         const Int_t nmaxevts=0, TMatrix *mrest=NULL);
+    Bool_t DefRefMatrix(const Int_t nmaxevts=0, TMatrix *mrest=NULL);
 
     ClassDef(MHMatrix, 1) // Multidimensional Matrix to store events