Changeset 1762

Timestamp:

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

Author:

wittek

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• manalysis/MHillasCalc.cc (modified) (1 diff)
• manalysis/MHillasSrcCalc.cc (modified) (3 diffs)
• manalysis/MHillasSrcCalc.h (modified) (2 diffs)
• manalysis/MSelBasic.cc (added)
• manalysis/MSelBasic.h (added)
• manalysis/MSelStandard.cc (added)
• manalysis/MSelStandard.h (added)
• mfilter/MFEventSelector.cc (modified) (3 diffs)
• mfilter/MFEventSelector.h (modified) (2 diffs)
• mhist/MHHillasExt.cc (modified) (3 diffs)
• mhist/MHHillasExt.h (modified) (1 diff)
• mhist/MHMatrix.cc (modified) (10 diffs)
• mhist/MHMatrix.h (modified) (4 diffs)

trunk/MagicSoft/Mars/manalysis/MHillasCalc.cc

@@ -84 +84 @@
     }
 
-    fHillas = (MHillas*)pList->FindCreateObj(fHilName);
+    fHillas = (MHillas*)pList->FindCreateObj("MHillas",fHilName);
     if (!fHillas)
         return kFALSE;

trunk/MagicSoft/Mars/manalysis/MHillasSrcCalc.cc

@@ -56 +56 @@
 // The default is "MHillasSrc"
 //
+//MHillasSrcCalc::MHillasSrcCalc(const char *src, const char *hil,
+//                               const char *name, const char *title)
+//    : fHillas(NULL), fSrcPos(NULL), fHillasSrc(NULL)
+//{
+//    fName  = name  ? name  : gsDefName.Data();
+//    fTitle = title ? title : gsDefTitle.Data();
+
+//    fSrcName    = src;
+//    fHillasName = hil;
+//    fHillasInput = "MHillas";
+//}
+// -------------------------------------------------------------------------
+//
 MHillasSrcCalc::MHillasSrcCalc(const char *src, const char *hil,
-                               const char *name, const char *title)
+                               const char *name, const char *title,
+                               const char *hilinput)
     : fHillas(NULL), fSrcPos(NULL), fHillasSrc(NULL)
 {
@@ -65 +79 @@
     fSrcName    = src;
     fHillasName = hil;
+    fHillasInput = hilinput;
 }
 
@@ -71 +86 @@
 Bool_t MHillasSrcCalc::PreProcess(MParList *pList)
 {
-    fHillas = (MHillas*)pList->FindObject("MHillas");
+    fHillas = (MHillas*)pList->FindObject(fHillasInput);
     if (!fHillas)
     {

trunk/MagicSoft/Mars/manalysis/MHillasSrcCalc.h

@@ -19 +19 @@
     TString     fSrcName;
     TString     fHillasName;
+    TString     fHillasInput;
 
     Int_t       fErrors;
@@ -29 +30 @@
 
 public:
+    //    MHillasSrcCalc(const char *src="MSrcPosCam", const char *hil="MHillasSrc",
+    //               const char *name=NULL, const char *title=NULL);
+
     MHillasSrcCalc(const char *src="MSrcPosCam", const char *hil="MHillasSrc",
-                   const char *name=NULL, const char *title=NULL);
+                   const char *name=NULL, const char *title=NULL,
+                   const char *hilinput="MHillas");
 
     ClassDef(MHillasSrcCalc, 1) // task to calculate the source position depandant hillas parameters

trunk/MagicSoft/Mars/mfilter/MFEventSelector.cc

@@ -94 +94 @@
 // the class description above.
 //
-MFEventSelector::MFEventSelector(const char *name, const char *title)
+MFEventSelector::MFEventSelector(const char *name, const char *title,
+                                 const char *read)
 : fNumTotalEvts(-1), fNumSelectEvts(-1), fSelRatio(-1), fNumSelectedEvts(0)
 {
     fName  = name  ? name  : gsDefName.Data();
     fTitle = title ? title : gsDefTitle.Data();
+
+    fRead = read;
 }
 
@@ -127 +130 @@
             return kFALSE;
         }
-        MRead *read = (MRead*)tlist->FindObject("MRead");
+        MRead *read = (MRead*)tlist->FindObject(fRead);
         if (!read)
         {
@@ -134 +137 @@
         }
         fNumTotalEvts = read->GetEntries();
+
+        *fLog << "MFEventSelector::PreProcess; fNumTotalEvts = " 
+              << fNumTotalEvts << endl;
 
         SetBit(kNumTotalFromFile);

trunk/MagicSoft/Mars/mfilter/MFEventSelector.h

@@ -24 +24 @@
 
     Bool_t  fResult;
+    TString fRead;
 
     void StreamPrimitive(ofstream &out) const;
@@ -38 +39 @@
 public:
     // MFEventSelector();
-    MFEventSelector(const char *name=NULL, const char *title=NULL);
+    MFEventSelector(const char *name=NULL, const char *title=NULL,
+                    const char *read="MRead");
     ~MFEventSelector();
 

trunk/MagicSoft/Mars/mhist/MHHillasExt.cc

@@ -57 +57 @@
 // Setup four histograms for Width, Length
 //
-MHHillasExt::MHHillasExt(const char *name, const char *title)
+MHHillasExt::MHHillasExt(const char *name, const char *title,
+                         const char *hil)
     : fMm2Deg(1), fUseMmScale(kTRUE)
 {
@@ -65 +66 @@
     fName  = name  ? name  : "MHHillasExt";
     fTitle = title ? title : "Container for extended Hillas histograms";
+
+    fHilName = hil;
+    //*fLog << "MHHillasExt : fHilName = " << fHilName << endl;
 
     //
@@ -139 +143 @@
 Bool_t MHHillasExt::SetupFill(const MParList *plist)
 {
-    TObject *obj = plist->FindObject("MHillas");
+    TObject *obj = plist->FindObject(fHilName, "MHillas");
     if (!obj)
     {
-        *fLog << err << dbginf << "Sorry 'MHillas' not found in parameter list... aborting." << endl;
+      *fLog << err << dbginf << "Sorry '" << fHilName 
+            << "' not found in parameter list... aborting." << endl;
         return kFALSE;
     }

trunk/MagicSoft/Mars/mhist/MHHillasExt.h

@@ -25 +25 @@
     Bool_t  fUseMmScale;
 
+    TString fHilName;
+
 public:
-    MHHillasExt(const char *name=NULL, const char *title=NULL);
+    //MHHillasExt(const char *name=NULL, const char *title=NULL);
+    MHHillasExt(const char *name=NULL, const char *title=NULL,
+                const char *hil="MHillas");
+
     ~MHHillasExt();
 

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 @@
 
 
+
+
+
+

trunk/MagicSoft/Mars/mhist/MHMatrix.h

@@ -13 +13 @@
 #endif
 
+#include <TArrayF.h>
+
 class TArrayI;
+class TArrayF;
 
+class TH1F;
 class MTask;
 class MParList;
@@ -32 +36 @@
 
     MDataArray *fData;  // List of data members (columns)
+
 
     enum {
@@ -94 +99 @@
     Int_t Read(const char *name);
 
+    Bool_t DefRefMatrix   (const Int_t refcolumn=-1, const TH1F *thsh=NULL, 
+                           const Int_t nmaxevts=0);
+
     ClassDef(MHMatrix, 1) // Multidimensional Matrix to store events
 };
@@ -99 +107 @@
 #endif
 
+