source: trunk/MagicSoft/Mars/mtemp/mpadova/macros/RanForestPD.C@ 6668

/* ======================================================================== *\
!
! *
! * 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 3/2003 <mailto:hengsteb@alwa02.physik.uni-siegen.de>!
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

#include "MMcEvt.hxx"

void RanForestPD(Int_t minsize=0)
{
    //
    // Create a empty Parameter List and an empty Task List
    // The tasklist is identified in the eventloop by its name
    //
    MParList  plist;

    MTaskList tlist;
    plist.AddToList(&tlist);

    MReadMarsFile  read("Events", "star_gamma_train.root");
    Float_t numgammas = read.GetEntries();

    read.AddFile("star_proton_train.root");
    read.AddFile("star_helium_train.root");

    Float_t numhadrons = read.GetEntries() - numgammas;

    // Fraction of gammas to be used in training:
    // (less than total sample, to speed up execution)
    //
    Float_t gamma_frac = 2.*(numhadrons / numgammas);


    read.DisableAutoScheme();

    tlist.AddToList(&read);

     MFParticleId fhadrons("MMcEvt", '!', MMcEvt::kGAMMA);
    tlist.AddToList(&fhadrons);

    MHMatrix matrixg("MatrixGammas");

    //    matrixg.AddColumn("floor(0.5+(cos(MMcEvt.fTelescopeTheta)/0.01))");
    //    matrixg.AddColumn("MMcEvt.fTelescopePhi");
    //    matrixg.AddColumn("MSigmabar.fSigmabar");

    matrixg.AddColumn("MHillas.fSize");
    matrixg.AddColumn("MHillasSrc.fDist");

    //    matrixg.AddColumn("MHillas.fWidth");
    //    matrixg.AddColumn("MHillas.fLength");

    // Scaled Width and Length:
    //
    matrixg.AddColumn("MHillas.fWidth/(-13.1618+(10.0492*log10(MHillas.fSize)))");   
    matrixg.AddColumn("MHillas.fLength/(-57.3784+(32.6131*log10(MHillas.fSize)))");

    //    matrixg.AddColumn("(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))");
    matrixg.AddColumn("MNewImagePar.fConc");

    //    matrixg.AddColumn("MNewImagePar.fConc1");
    //    matrixg.AddColumn("MNewImagePar.fAsym");
    //    matrixg.AddColumn("MHillasExt.fM3Trans");
    //    matrixg.AddColumn("abs(MHillasSrc.fHeadTail)");
    //    matrixg.AddColumn("MNewImagePar.fLeakage1");

    //
    // Third moment, with sign referred to source position:
    //
    matrixg.AddColumn("MHillasExt.fM3Long*MHillasSrc.fCosDeltaAlpha/abs(MHillasSrc.fCosDeltaAlpha)");

    //#    matrixg.AddColumn("MHillasSrc.fAlpha");

    //
    // SIZE cut:
    //
    TString dummy("(MHillas.fSize<");
    dummy += minsize;
    dummy += ")";

    // AM, useful FOR High Energy only:
    //    dummy += " || (MImagePar.fNumSatPixelsLG>0)";

    MContinue SizeCut(dummy);
    tlist.AddToList(&SizeCut);

    // AM TEST!!!, no muons in train
    //    dummy = "(MMcEvt.fMuonCphFraction>0.1)";
    //    MContinue MuonsCut(dummy);
    //    tlist.AddToList(&MuonsCut);

    plist.AddToList(&matrixg);

    MHMatrix matrixh("MatrixHadrons");
    matrixh.AddColumns(matrixg.GetColumns());
    plist.AddToList(&matrixh);

    MFillH fillmath("MatrixHadrons");
    fillmath.SetFilter(&fhadrons);
    tlist.AddToList(&fillmath);

    MContinue skiphad(&fhadrons);
    tlist.AddToList(&skiphad);

    MFEventSelector reduce_training_gammas;
    reduce_training_gammas.SetSelectionRatio(gamma_frac);
    tlist.AddToList(&reduce_training_gammas);

    MFillH fillmatg("MatrixGammas");
    fillmatg.SetFilter(&reduce_training_gammas);
    tlist.AddToList(&fillmatg);

    //
    // Create and setup the eventloop
    //
    MEvtLoop evtloop;
    evtloop.SetParList(&plist);

    MProgressBar bar;
    evtloop.SetProgressBar(&bar);
    //
    // Execute your analysis
    //
    if (!evtloop.Eventloop())
        return;

    tlist.PrintStatistics();

    // ---------------------------------------------------------------
    // ---------------------------------------------------------------
    // second event loop: the trees of the random forest are grown,
    // the event loop is now actually a tree loop (loop of training
    // process)
    // ---------------------------------------------------------------
    // ---------------------------------------------------------------

    MTaskList tlist2;
    plist.Replace(&tlist2);

    MRanForestGrow rfgrow2;
    rfgrow2.SetNumTrees(100);
    rfgrow2.SetNumTry(3);
    rfgrow2.SetNdSize(10);

    MWriteRootFile rfwrite2("RF.root");
    rfwrite2.AddContainer("MRanTree","Tree");       //save all trees
    MFillH fillh2("MHRanForestGini");

    tlist2.AddToList(&rfgrow2);
    tlist2.AddToList(&rfwrite2);
    tlist2.AddToList(&fillh2);

    // gRandom is accessed from MRanForest (-> bootstrap aggregating)
    // and MRanTree (-> random split selection) and should be initialized
    // here if you want to set a certain random number generator
    if(gRandom)
        delete gRandom;
    //    gRandom = new TRandom3(0); // Takes seed from the computer clock
    gRandom = new TRandom3(); // Uses always  same seed
    //
    // Execute tree growing (now the eventloop is actually a treeloop)
    //

    evtloop.SetProgressBar(&bar);

    if (!evtloop.Eventloop())
        return;

    tlist2.PrintStatistics();

    plist.FindObject("MHRanForestGini")->DrawClone();

    // ---------------------------------------------------------------
    // ---------------------------------------------------------------
    // third event loop: the control sample (star2.root) is processed
    // through the previously grown random forest,
    //
    // the histograms MHHadronness (quality of g/h-separation) and
    // MHRanForest are created and displayed.
    // MHRanForest shows the convergence of the classification error
    // as function of the number of grown (and combined) trees
    // and tells the user how many trees are actually needed in later
    // classification tasks.
    // ---------------------------------------------------------------
    // ---------------------------------------------------------------

    MTaskList tlist3;

    plist.Replace(&tlist3);

    MReadMarsFile  read3("Events", "star_gamma_test.root");

    read3.AddFile("star_proton_test.root");
    read3.AddFile("star_helium_test.root");

    read3.DisableAutoScheme();
    tlist3.AddToList(&read3);

    tlist3.AddToList(&SizeCut);

    MRanForestCalc calc;
    tlist3.AddToList(&calc);

    // parameter containers you want to save
    //

    TString outfname = "star_S";
    outfname += minsize;
    outfname += "_all.root";
    MWriteRootFile write(outfname, "recreate");

    write.AddContainer("MMcEvt",         "Events");
    write.AddContainer("MHillas",        "Events");
    write.AddContainer("MHillasExt",     "Events");
    write.AddContainer("MImagePar",      "Events");
    write.AddContainer("MNewImagePar",   "Events");
    write.AddContainer("MHillasSrc",     "Events");
    write.AddContainer("MHadronness",    "Events");

    write.AddContainer("MRawRunHeader", "RunHeaders");
    write.AddContainer("MSrcPosCam",    "RunHeaders");
    write.AddContainer("MMcRunHeader",  "RunHeaders");

    /*
    write.AddContainer("MMcTrig",       "Events");
    write.AddContainer("MRawEvtData",   "Events");
    write.AddContainer("MRawEvtHeader", "Events");
    */


    MFillH fillh3a("MHHadronness");
    MFillH fillh3b("MHRanForest");

    tlist3.AddToList(&fillh3a);
    tlist3.AddToList(&fillh3b);


    tlist3.AddToList(&write);

    evtloop.SetProgressBar(&bar);

    //
    // Execute your analysis
    //
    if (!evtloop.Eventloop())
        return;

    tlist3.PrintStatistics();

    plist.FindObject("MHRanForest")->DrawClone();

    TCanvas *had = new TCanvas;
    had->cd();
    plist.FindObject("MHHadronness")->Draw();

    TString gifname = "had_S";
    gifname += minsize;
    gifname += ".gif";

    had->SaveAs(gifname);
    had->Close();

    plist.FindObject("MHHadronness")->DrawClone();
    plist.FindObject("MHHadronness")->Print();//*/

    return;
}