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

Timestamp:

02/19/03 13:47:29 (22 years ago)

Author:

wittek

Message:

*** empty log message ***

File:

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

trunk/MagicSoft/Mars/mhist/MHMatrix.cc

@@ -16 +16 @@
 !
 !
-!   Author(s): Thomas Bretz  2002 <mailto:tbretz@astro.uni-wuerzburg.de>
+!   Author(s): Thomas Bretz   2002 <mailto:tbretz@astro.uni-wuerzburg.de>
+!              Rudy Boeck     2003 <mailto:
+!              Wolfgang Wittek2003 <mailto:wittek@mppmu.mpg.de>
 !
-!   Copyright: MAGIC Software Development, 2000-2002
+!   Copyright: MAGIC Software Development, 2000-2003
 !
 !
@@ -52 +54 @@
 #include <TArrayI.h>
 
+#include <TH1.h>
+#include <TCanvas.h>
+#include <TRandom3.h>
+
 #include "MLog.h"
 #include "MLogManip.h"
@@ -134 +140 @@
 }
 
+// --------------------------------------------------------------------------
+//
 void MHMatrix::AddColumns(MDataArray *matrix)
 {
@@ -341 +349 @@
 }
 
+// --------------------------------------------------------------------------
+//
 const TMatrix *MHMatrix::InvertPosDef()
 {
@@ -541 +551 @@
 }
 
+// --------------------------------------------------------------------------
+//
 const TArrayI MHMatrix::GetIndexOfSortedColumn(Int_t ncol, Bool_t desc) const
 {
@@ -558 +570 @@
 }
 
+// --------------------------------------------------------------------------
+//
 void MHMatrix::SortMatrixByColumn(Int_t ncol, Bool_t desc)
 {
@@ -577 +591 @@
 }
 
+// --------------------------------------------------------------------------
+//
 Bool_t MHMatrix::Fill(MParList *plist, MTask *read)
 {
@@ -617 +633 @@
 // --------------------------------------------------------------------------
 //
-// Return a comma seperated list of all data members used in the matrix.
-// This is mainly used in MTask::AddToBranchList
 //
 void MHMatrix::ReduceNumberOfRows(UInt_t numrows, const TString opt)
@@ -674 +688 @@
 }
 
+// ------------------------------------------------------------------------
+//
+// Define the reference matrix
+//   refcolumn  number of the column 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
+//
+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();
+   Int_t   fColLwb = fM.GetColLwb();
+   Int_t   fRowLwb = fM.GetRowLwb();
+
+   // print characteristics
+   //
+   //*fLog << "MHMatrix::CalcNormFactors; Event matrix (input): " << endl;
+   //*fLog << "          fNrows, fNcols, fColLwb, fRowLwb = " << fNrows << ",  "
+   //      << fNcols << ",  " << fColLwb << ",  " << fRowLwb << endl;
+
+   
+   //---------------------------------------------------------
+   //
+   // if refcol < 0 : select reference events randomly
+   //                 i.e. set the normaliztion factotrs equal to 1.0
+
+   if (refcolumn<0) 
+   {
+     frefcol = -refcolumn;
+   }
+   else
+   {
+     frefcol = refcolumn;
+   }
+   
+
+   //---------------------------------------------------------
+
+   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;
+   }
+
+   // if refcolumn<0 set target distribution equal to the real distribution
+   // in order to obtain normalization factors = 1.0
+   //if (refcolumn<0) 
+   //{
+   //  for (Int_t j=1; j<=fnbins; j++) 
+   //  {
+   //    float cont = fHth-> GetBinContent(j);
+   //    fHthsh->SetBinContent(j, cont);
+   //  }
+   //}
+
+   //---------------------------------------------------------
+   // ===== 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.0 )
+       b = fHthsh->GetBinContent(j);
+     else
+       b = a;
+
+     if ( a<0.0  || b<0.0 )
+     {
+       *fLog << "MHMatrix::DefRefMatrix; renormalization of input matrix is not possible"   << endl;
+       *fLog << "          a, b = " << a << ",  " << b << endl;
+        return kFALSE;
+     }
+
+     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= (fM(ind[ir],frefcol)-ffrombin)/fdbin;
+     cumul_weight[indref] += fnormfac[indref];
+     if (cumul_weight[indref]>0.5)    
+     {
+       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<13;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);
+     if (a>0) *fLog << j << "  "<< a << "   ";
+   }
+   *fLog <<endl;
+
+   //---------------------------------------------------------
+   // ==== plot four histograms 
+   TCanvas *th1 = new TCanvas (fName, fName, 1);   
+   th1->Divide(2,2);
+
+   th1->cd(1);
+   ((TH1F*)fHthsh)->Draw();      // target
+
+   th1->cd(2);
+   ((TH1F*)fHth)->Draw();        // real histogram before
+
+   th1->cd(3);
+   ((TH1F*)fHthd) -> Draw();     // correction factors
+
+   th1->cd(4);
+   ((TH1F*)fHthaft) -> Draw();   // histogram after
+
+   //---------------------------------------------------------
+   // --- write onto output file
+   //
+   //TFile *outfile = new TFile(fileNameout, "RECREATE", "");
+   //Mnew.Write(fileNameout);
+   //outfile->Close();
+
+   return kTRUE;
+}
+
 // --------------------------------------------------------------------------
 //
@@ -691 +1005 @@
 
 
+
+
+
+