Changeset 7142 for trunk/MagicSoft/Mars/mranforest

Timestamp:

06/10/05 13:10:09 (19 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mranforest

Files:

• MRFEnergyEst.cc (modified) (9 diffs)
• MRFEnergyEst.h (modified) (3 diffs)
• MRanTree.cc (modified) (4 diffs)
• MRanTree.h (modified) (2 diffs)

trunk/MagicSoft/Mars/mranforest/MRFEnergyEst.cc

@@ -69 +69 @@
 //
 MRFEnergyEst::MRFEnergyEst(const char *name, const char *title)
-    : fNumTrees(-1), fNumTry(-1), fNdSize(-1), fData(0), fEnergyEst(0),
-    fTestMatrix(0)
+    : fNumTrees(-1), fNumTry(-1), fNdSize(-1), fDebug(kFALSE),
+    fData(0), fEnergyEst(0), fTestMatrix(0)
 {
     fName  = name  ? name  : gsDefName.Data();
@@ -91 +91 @@
     }
 
-    const Int_t ncols = matrixtrain.GetM().GetNcols();
-    const Int_t nrows = matrixtrain.GetM().GetNrows();
+    const TMatrix &m = matrixtrain.GetM();
+
+    const Int_t ncols = m.GetNcols();
+    const Int_t nrows = m.GetNrows();
     if (ncols<=0 || nrows <=0)
     {
@@ -114 +116 @@
     }
 
+    const Int_t nobs = 3; // Number of obsolete columns
+
     MDataArray usedrules;
-    for(Int_t i=0; i<ncols-3; i++) // -3 is important!!!
+    for(Int_t i=0; i<ncols-nobs; i++) // -3 is important!!!
         usedrules.AddEntry((*matrixtrain.GetColumns())[i].GetRule());
-
-    const TMatrix &m = matrixtrain.GetM();
 
     // training of RF for each energy bin
@@ -140 +142 @@
         }
 
-        const Bool_t valid = irow1*irow0>0;
-
-        if (!valid)
+        const Bool_t invalid = irow1==0 || irow0==0;
+
+        if (invalid)
             *fLog << warn << "WARNING - Skipping";
         else
             *fLog << inf << "Training RF for";
 
-        *fLog << " energy bin " << ie << " (" << grid[ie] << ", " << grid[ie+1] << ")" << endl;
-
-        if (!valid)
+        *fLog << " energy bin " << ie << " (" << grid[ie] << ", " << grid[ie+1] << ") " << irow0 << " " << irow1 << endl;
+
+        if (invalid)
             continue;
 
-        gLog.SetNullOutput(kTRUE);
+        if (fDebug)
+            gLog.SetNullOutput(kTRUE);
 
         // last column must be removed (true energy col.)
-        mat1.ResizeTo(irow1, ncols-1);
-        mat0.ResizeTo(irow0, ncols-1);
+        mat1.ResizeTo(irow1, ncols-nobs);
+        mat0.ResizeTo(irow0, ncols-nobs);
 
         // create MHMatrix as input for RF
@@ -183 +186 @@
         evtloop.SetParList(&plist);
 
-        gLog.SetNullOutput(kFALSE);
+        if (fDebug)
+            gLog.SetNullOutput(kFALSE);
 
         if (!evtloop.Eventloop())
             return kFALSE;
 
+        const Double_t E = (log10(grid[ie])+log10(grid[ie+1]))/2;
+
         // save whole forest
         MRanForest *forest=(MRanForest*)plist.FindObject("MRanForest");
-        const TString title = Form("%f", TMath::Sqrt(grid[ie]*grid[ie+1]));
-        //const TString title = Form("%f", (grid[ie]+grid[ie+1])/2);
+        const TString title = Form("%.10f", E);
         forest->SetTitle(title);
         forest->Write(title);
@@ -203 +208 @@
     return kTRUE;
 }
-/*
-Int_t MRFEnergyEst::Test(const MHMatrix &matrixtest)
-{
-    gLog.Separator("MRFEnergyEst - Test");
-
-    const Int_t ncols = matrixtest.GetM().GetNcols();
-    const Int_t nrows = matrixtest.GetM().GetNrows();
-
-    if (ncols<=0 || nrows <=0)
-    {
-        *fLog << err << dbginf << "No. of columns or no. of rows of matrixtrain equal 0 ... aborting." << endl;
-        return kFALSE;
-    }
-
-    if (!ReadForests())
-    {
-        *fLog << err << dbginf << "Reading RFs failed... aborting." << endl;
-        return kFALSE;
-    }
-
-    const Int_t nbins=fEForests.GetSize();
-
-    Double_t elow =  FLT_MAX;
-    Double_t eup  = -FLT_MAX;
-
-    for(Int_t j=0;j<nbins;j++)
-    {
-        elow = TMath::Min(atof(fEForests[j]->GetTitle()), elow);
-        eup  = TMath::Max(atof(fEForests[j]->GetTitle()), eup);
-    }
-
-    TH1F hres("hres", "", 1000, -10, 10);
-    TH2F henergy("henergy", "",100, elow, eup,100, elow, eup);
-
-    const TMatrix &m=matrixtest.GetM();
-    for(Int_t i=0;i<nrows;i++)
-    {
-        Double_t etrue = m(i,ncols-1);
-        Double_t eest  = 0;
-        Double_t hsum  = 0;
-
-        for(Int_t j=0;j<nbins;j++)
-        {
-            const TVector  v = TMatrixFRow_const(m,i);
-
-            const Double_t h = ((MRanForest*)fEForests[j])->CalcHadroness(v);
-            const Double_t e = atof(fEForests[j]->GetTitle());
-
-            hsum += h;
-            eest += h*e;
-        }
-        eest /= hsum;
-        eest  = pow(10., eest);
-
-        //if (etrue>80.)
-        //    hres.Fill((eest-etrue)/etrue);
-
-        henergy.Fill(log10(etrue),log10(eest));
-    }
-
-    if(gStyle)
-        gStyle->SetOptFit(1);
-
-    // show results
-    TCanvas *c=MH::MakeDefCanvas("c","",300,900);
-
-    c->Divide(1,3);
-    c->cd(1);
-    henergy.SetTitle("Estimated vs true energy");
-    henergy.GetXaxis()->SetTitle("log10(E_{true}[GeV])");
-    henergy.GetYaxis()->SetTitle("log10(E_{est}[GeV])");
-    henergy.DrawCopy();
-
-    c->cd(2);
-    TH1F *hptr=(TH1F*)henergy.ProfileX("_px",1,100,"S");
-    hptr->SetTitle("Estimated vs true energy - profile");
-    hptr->GetXaxis()->SetTitle("log10(E_{true}[GeV])");
-    hptr->GetYaxis()->SetTitle("log10(E_{est}[GeV])");
-    hptr->DrawCopy();
-
-    c->cd(3);
-    hres.Fit("gaus");
-    hres.SetTitle("Energy resolution for E_{true}>80Gev");
-    hres.GetXaxis()->SetTitle("(E_{estimated}-E_{true})/E_{true})");
-    hres.GetYaxis()->SetTitle("counts");
-    hres.DrawCopy();
-
-    c->GetPad(1)->SetGrid();
-    c->GetPad(2)->SetGrid();
-    c->GetPad(3)->SetGrid();
-
-    return kTRUE;
-}
-*/
-Int_t MRFEnergyEst::ReadForests(MParList *plist)
+
+Int_t MRFEnergyEst::ReadForests(MParList &plist)
 {
     TFile fileRF(fFileName,"read");
@@ -321 +233 @@
 
         fEForests.Add(forest);
-
-    }
-
-    if (plist)
-    {
-        fData = (MDataArray*)plist->FindCreateObj("MDataArray");
-        fData->Read("rules");
-    }
-
-    fileRF.Close();
+    }
+
+    if (fData->Read("rules")<=0)
+    {
+        *fLog << err << "ERROR - Reading 'rules' from file " << fFileName << endl;
+        return kFALSE;
+    }
 
     return kTRUE;
@@ -341 +250 @@
         return kFALSE;
 
-    if (!ReadForests(plist))
+    fData = (MDataArray*)plist->FindCreateObj("MDataArray");
+    if (!fData)
+        return kFALSE;
+
+    if (!ReadForests(*plist))
     {
         *fLog << err << "Reading RFs failed... aborting." << endl;
         return kFALSE;
     }
+
+    *fLog << inf << "RF read from " << fFileName << endl;
 
     if (fTestMatrix)
         return kTRUE;
 
-    if (!fData)
-    {
-        *fLog << err << "MDataArray not found... aborting." << endl;
-        return kFALSE;
-    }
+    fData->Print();
 
     if (!fData->PreProcess(plist))
@@ -387 +298 @@
 
         hsum += h;
-        eest += h*e;
-    }
-
-    fEnergyEst->SetVal(eest/hsum);
+        eest += e*h;
+    }
+
+    fEnergyEst->SetVal(pow(10, eest/hsum));
     fEnergyEst->SetReadyToSave();
 
     return kTRUE;
 }
+
+// --------------------------------------------------------------------------
+//
+//
+Int_t MRFEnergyEst::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
+{
+    Bool_t rc = kFALSE;
+    if (IsEnvDefined(env, prefix, "FileName", print))
+    {
+        rc = kTRUE;
+        SetFileName(GetEnvValue(env, prefix, "FileName", fFileName));
+    }
+    if (IsEnvDefined(env, prefix, "Debug", print))
+    {
+        rc = kTRUE;
+        SetDebug(GetEnvValue(env, prefix, "Debug", fDebug));
+    }
+    return rc;
+}

trunk/MagicSoft/Mars/mranforest/MRFEnergyEst.h

@@ -23 +23 @@
     Int_t        fNdSize;    // Training parameters
 
+    Bool_t       fDebug;     // Debugging of eventloop while training on/off
+
     TString      fFileName;
     TObjArray    fEForests;
@@ -34 +36 @@
     Int_t Process();
 
-    Int_t ReadForests(MParList *plist=NULL);
+    Int_t ReadForests(MParList &plist);
+
+    Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print);
 
 public:
@@ -41 +45 @@
     void  SetFileName(TString str)     { fFileName = str; }
 
-    void  SetNumTrees(Int_t n=-1)   { fNumTrees = n; }
-    void  SetNdSize(Int_t n=-1)     { fNdSize   = n; }
-    void  SetNumTry(Int_t n=-1)     { fNumTry   = n; }
+    void  SetNumTrees(Int_t n=-1)      { fNumTrees = n; }
+    void  SetNdSize(Int_t n=-1)        { fNdSize   = n; }
+    void  SetNumTry(Int_t n=-1)        { fNumTry   = n; }
 
     void  SetTestMatrix(MHMatrix *m=0) { fTestMatrix=m; }
     void  InitMapping(MHMatrix *m=0)   { fTestMatrix=m; }
+
+    void SetDebug(Bool_t b=kTRUE)      { fDebug = b; }
 
     Int_t Train(const MHMatrix &n, const TArrayD &grid);

trunk/MagicSoft/Mars/mranforest/MRanTree.cc

@@ -447 +447 @@
     // are coded into fBestVar:
     // status of node kt = TMath::Sign(1,fBestVar[kt])
-    // class  of node kt = fBestVar[kt]+2 (class defined by larger
-    //  node population, actually not used)
+    // class  of node kt = fBestVar[kt]+2
+    //  (class defined by larger node population, actually not used)
     // hadronness assigned to node kt = fBestSplit[kt]
 
@@ -463 +463 @@
 }
 
-Double_t MRanTree::TreeHad(const TMatrixRow &event)
+Double_t MRanTree::TreeHad(const TMatrixFRow_const &event)
 {
     Int_t kt=0;
@@ -504 +504 @@
 Bool_t MRanTree::AsciiWrite(ostream &out) const
 {
-    TString str;
+    out.width(5);
+    out << fNumNodes << endl;
+
     Int_t k;
-
-    out.width(5);out<<fNumNodes<<endl;
-
-    for (k=0;k<fNumNodes;k++)
-    {
-        str=Form("%f",GetBestSplit(k));
+    for (k=0; k<fNumNodes; k++)
+    {
+        TString str=Form("%f", GetBestSplit(k));
 
         out.width(5);  out << k;
@@ -521 +520 @@
         out.width(5);  out << GetNodeClass(k);
     }
-    out<<endl;
-
-    return k==fNumNodes;
-}
+    out << endl;
+
+    return kTRUE;
+}

trunk/MagicSoft/Mars/mranforest/MRanTree.h

@@ -16 +16 @@
 class TMatrix;
 class TMatrixRow;
+class TMatrixFRow_const;
 class TVector;
 class TRandom;
@@ -84 +85 @@
 
     Double_t TreeHad(const TVector &event);
-    Double_t TreeHad(const TMatrixRow &event);
+    Double_t TreeHad(const TMatrixFRow_const &event);
     Double_t TreeHad(const TMatrix &m, Int_t ievt);
     Double_t TreeHad();