source: trunk/MagicSoft/Mars/mranforest/MRFEnergyEst.cc@ 7143

/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Thomas Hengstebeck 2/2005 <mailto:hengsteb@physik.hu-berlin.de>
!
!   Copyright: MAGIC Software Development, 2000-2005
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  MRFEnergyEst
//
//
////////////////////////////////////////////////////////////////////////////
#include "MRFEnergyEst.h"

#include <TFile.h>
#include <TList.h>

#include <TH1.h>
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TVector.h>

#include "MHMatrix.h"

#include "MLog.h"
#include "MLogManip.h"

#include "MParList.h"
#include "MTaskList.h"
#include "MEvtLoop.h"

#include "MRanTree.h"
#include "MRanForest.h"
#include "MRanForestGrow.h"

#include "MData.h"
#include "MParameters.h"

ClassImp(MRFEnergyEst);

using namespace std;

static const TString gsDefName  = "MRFEnergyEst";
static const TString gsDefTitle = "RF for energy estimation";

// --------------------------------------------------------------------------
//
//  Default constructor. Set the name and title of this object
//
MRFEnergyEst::MRFEnergyEst(const char *name, const char *title)
    : fNumTrees(-1), fNumTry(-1), fNdSize(-1), fDebug(kFALSE),
    fData(0), fEnergyEst(0), fTestMatrix(0)
{
    fName  = name  ? name  : gsDefName.Data();
    fTitle = title ? title : gsDefTitle.Data();
}

// --------------------------------------------------------------------------
//
// Train the RF with the goven MHMatrix. The last column must contain the
// True energy.
//
Int_t MRFEnergyEst::Train(const MHMatrix &matrixtrain, const TArrayD &grid)
{
    gLog.Separator("MRFEnergyEst - Train");

    if (!matrixtrain.GetColumns())
    {
        *fLog << err << "ERROR - MHMatrix does not contain rules... abort." << endl;
        return kFALSE;
    }

    const TMatrix &m = matrixtrain.GetM();

    const Int_t ncols = m.GetNcols();
    const Int_t nrows = m.GetNrows();
    if (ncols<=0 || nrows <=0)
    {
        *fLog << err << "ERROR - No. of columns or no. of rows of matrixtrain equal 0 ... abort." << endl;
        return kFALSE;
    }

    const Int_t nbins = grid.GetSize()-1;
    if (nbins<=0)
    {
        *fLog << err << "ERROR - Energy grid not vaild... abort." << endl;
        return kFALSE;
    }

    // rules (= combination of image par) to be used for energy estimation
    TFile fileRF(fFileName, "recreate");
    if (!fileRF.IsOpen())
    {
        *fLog << err << "ERROR - File to store RFs could not be opened... abort." << endl;
        return kFALSE;
    }

    const Int_t nobs = 3; // Number of obsolete columns

    MDataArray usedrules;
    for(Int_t i=0; i<ncols-nobs; i++) // -3 is important!!!
        usedrules.AddEntry((*matrixtrain.GetColumns())[i].GetRule());

    // training of RF for each energy bin
    for (Int_t ie=0; ie<nbins; ie++)
    {
        TMatrix mat1(nrows, ncols);
        TMatrix mat0(nrows, ncols);

        // prepare matrix for current energy bin
        Int_t irow1=0;
        Int_t irow0=0;

        for (Int_t j=0; j<nrows; j++)
        {
            const Double_t energy = m(j,ncols-1);

            if (energy>grid[ie] && energy<=grid[ie+1])
                TMatrixFRow(mat1, irow1++) = TMatrixFRow_const(m,j);
            else
                TMatrixFRow(mat0, irow0++) = TMatrixFRow_const(m,j);
        }

        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] << ") " << irow0 << " " << irow1 << endl;

        if (invalid)
            continue;

        if (fDebug)
            gLog.SetNullOutput(kTRUE);

        // last column must be removed (true energy col.)
        mat1.ResizeTo(irow1, ncols-nobs);
        mat0.ResizeTo(irow0, ncols-nobs);

        // create MHMatrix as input for RF
        MHMatrix matrix1(mat1, "MatrixHadrons");
        MHMatrix matrix0(mat0, "MatrixGammas");

        matrix1.AddColumns(&usedrules);
        matrix0.AddColumns(&usedrules);

        // training of RF
        MTaskList tlist;
        MParList plist;
        plist.AddToList(&tlist);
        plist.AddToList(&matrix0);
        plist.AddToList(&matrix1);

        MRanForestGrow rfgrow;
        rfgrow.SetNumTrees(fNumTrees); // number of trees
        rfgrow.SetNumTry(fNumTry);     // number of trials in random split selection
        rfgrow.SetNdSize(fNdSize);     // limit for nodesize

        tlist.AddToList(&rfgrow);
    
        MEvtLoop evtloop;
        evtloop.SetDisplay(fDisplay);
        evtloop.SetParList(&plist);

        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("%.10f", E);
        forest->SetTitle(title);
        forest->Write(title);
    }

    // save rules
    usedrules.Write("rules");

    fileRF.Close();

    return kTRUE;
}

Int_t MRFEnergyEst::ReadForests(MParList &plist)
{
    TFile fileRF(fFileName,"read");
    if (!fileRF.IsOpen())
    {
        *fLog << err << dbginf << "File containing RFs could not be opened... aborting." << endl;
        return kFALSE;
    }

    fEForests.Delete();

    TIter Next(fileRF.GetListOfKeys());
    TObject *o=0;
    while ((o=Next()))
    {
        MRanForest *forest;
        fileRF.GetObject(o->GetName(), forest);
        if (!forest)
            continue;

        forest->SetTitle(o->GetTitle());
        forest->SetBit(kCanDelete);

        fEForests.Add(forest);
    }

    if (fData->Read("rules")<=0)
    {
        *fLog << err << "ERROR - Reading 'rules' from file " << fFileName << endl;
        return kFALSE;
    }

    return kTRUE;
}

Int_t MRFEnergyEst::PreProcess(MParList *plist)
{
    fEnergyEst = (MParameterD*)plist->FindCreateObj("MParameterD", "MEnergyEst");
    if (!fEnergyEst)
        return kFALSE;

    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;

    fData->Print();

    if (!fData->PreProcess(plist))
    {
        *fLog << err << "PreProcessing of the MDataArray failed... aborting." << endl;
        return kFALSE;
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
//
Int_t MRFEnergyEst::Process()
{
    TVector event;
    if (fTestMatrix)
        *fTestMatrix >> event;
    else
        *fData >> event;

    Double_t eest = 0;
    Double_t hsum = 0;

    TIter Next(&fEForests);
    MRanForest *rf = 0;
    while ((rf=(MRanForest*)Next()))
    {
        const Double_t h = rf->CalcHadroness(event);
        const Double_t e = atof(rf->GetTitle());

        hsum += h;
        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;
}