source: trunk/MagicSoft/Mars/mtemp/mpadova/macros/RanForestDISP.C@ 6396

/* ======================================================================== *\
   !
   ! *
   ! * 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): Abelardo Moralejo 1/2005 <mailto:moralejo@pd.infn.it>!
   !
   !   Copyright: MAGIC Software Development, 2000-2003
   !
   !
   \* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////////
//
// macro RanForestDISP.C
//
// Version of Random Forest macro, intended to run on real data without assuming 
// any particular source position in the camera FOV. So we do not use the DIST
// parameter, but rather (one version of) DISP method. We do not use either the 
// sign of the third longitudinal moment with respect to the center of the camera 
// (since the source can be somewhere else like in the case of wobble mode).
//
//////////////////////////////////////////////////////////////////////////////////


#include "MMcEvt.hxx"

void RanForestDISP(Int_t minsize=1000)
{
    Int_t mincorepixels = 6; // minimum number of core pixels.
//  Int_t mincorepixels = 0;

    TString skipevents;
    //  skipevents = "(MMcEvt.fImpact>15000.)";
    //  skipevents = "(MNewImagePar.fLeakage1>0.001)";
    //  skipevents = "(MHillas.fSize>1000)";

    // Skip calibration events, leaking events and flash events:
    skipevents = "(MNewImagePar.fLeakage1>0.1) || ({1.5+log10(MNewImagePar.fConc)*2.33333}-{6.5-log10(MHillas.fSize)}>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", "/data1/magic/mc_data/root/Period021_0.73_mirror/FixedWindowPeakSearch/WOBBLE_Calib_6slices/tc_3.84_2.56/star_gamma_train_new.root");
    Float_t numgammas = read.GetEntries();

    read.AddFile("star_train_20050110_CrabNebulaW1.root");
    Float_t numhadrons = read.GetEntries() - numgammas;

    // Fraction of gammas to be used in training:
    //    Float_t gamma_frac = 1.5*(numhadrons / numgammas);
    Float_t gamma_frac = 1.;

    // Total hadron number to be processed in training:
    // If you run RF over different data runs, using always the same # of training gammas,
    // you should also use always a fixed number of hadrons for training (although the data 
    // runs may be of different duration). Otherwise the meaning of the output hadronness 
    // would be different for the different subsamples, since its absolute values depend on 
    // the statistics of the training sample. The number set here is the number BEFORE the
    // "a priori cuts" set above through "skipevents".

    Float_t select_n_hadrons = 37000.;  
    Float_t hadron_frac = select_n_hadrons / numhadrons;

    read.DisableAutoScheme();

    tlist.AddToList(&read);

    MFParticleId fgamma("MMcEvt", '=', MMcEvt::kGAMMA);
    tlist.AddToList(&fgamma);

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


    MHMatrix matrixg("MatrixGammas");

    // TEST TEST TEST
//     matrixg.AddColumn("MMcEvt.fImpact");
//     matrixg.AddColumn("MMcEvt.fMuonCphFraction");

    //    matrixg.AddColumn("MNewImagePar.fLeakage1");
    //    matrixg.AddColumn("floor(0.5+(cos(MMcEvt.fTelescopeTheta)/0.01))");

    //    matrixg.AddColumn("MMcEvt.fTelescopePhi");
    //    matrixg.AddColumn("MSigmabar.fSigmabar");


    matrixg.AddColumn("MHillas.fSize");

    // We do not use DIST, since it depends on knowing the source position on camera.
    // In wobble mode the position changes with time. Besides, we intend to do a 2-d map
    // with the DISP method, and so we cannot use DIST in the cuts (since it would amde the
    // camera center a "provoleged" position).

    matrixg.AddColumn("(1.-(MHillas.fWidth/MHillas.fLength))*(0.9776+(0.101062*log10(MHillas.fSize)))");

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


    // Scaling from real data?
    //    matrixg.AddColumn("MHillas.fWidth/(-8.37847+(9.58826*log10(MHillas.fSize)))");   
    //    matrixg.AddColumn("MHillas.fLength/(-40.1409+(29.3044*log10(MHillas.fSize)))");

    // Scaling from MC:
    matrixg.AddColumn("(MHillas.fWidth*3.37e-3)/(-0.028235+(0.03231*log10(MHillas.fSize)))");   
    matrixg.AddColumn("(MHillas.fLength*3.37e-3)/(-0.13527+(0.09876*log10(MHillas.fSize)))");

    matrixg.AddColumn("MNewImagePar.fConc");

    // TEST TEST TEST TEST
    //    matrixg.AddColumn("MConcentration.Conc7");

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


    // We cannot use the 3rd moment with a sign referred to the nominal source position, if we
    // assume we do not know it a priori. We cannot use either the sign as referred to the 
    // reconstructed source position through the DISP method, since that method uses precisely 
    // fM3Long to solve the "head-tail" asimetry and hence all showers would be "gamma-like" from
    // that point of view. We add however the 3rd moment with no sign to the separation parameters.

    matrixg.AddColumn("abs(MHillasExt.fM3Long)");

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

    TString dummy("(MHillas.fSize<");
    dummy += minsize;

    // AM TEST
    dummy += ") || (MNewImagePar.fNumCorePixels<";
    dummy += mincorepixels;

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

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


    if (skipevents != "")
    {
        MContinue SCut(skipevents);
        tlist.AddToList(&SCut);
    }

    plist.AddToList(&matrixg);

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

    MFEventSelector reduce_training_hadrons;
    reduce_training_hadrons.SetSelectionRatio(hadron_frac);
    MFilterList hadfilter;
    hadfilter.AddToList(&reduce_training_hadrons);
    hadfilter.AddToList(&fhadrons);
    tlist.AddToList(&hadfilter);

    MFillH fillmath("MatrixHadrons");
    fillmath.SetFilter(&hadfilter);
    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 = new MProgressBar;
    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);

    //
    //    matrixh.ShuffleRows(9999);
    //    matrixg.ShuffleRows(1111);
    //

    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_20050110_CrabNebulaW1.root");
    //    MReadMarsFile  read3("Events", "/data1/magic/mc_data/root/Period021_0.73_mirror/FixedWindowPeakSearch/WOBBLE_Calib_6slices/tc_3.84_2.56/star_gamma_test.root");

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

    if (skipevents != "")
        tlist3.AddToList(&SCut);

    MRanForestCalc calc;
    tlist3.AddToList(&calc);


    TString outfname = "star_S";
    outfname += minsize;
    outfname += "_20050110_CrabNebulaW1.root";

    //    outfname += "_mcgammas.root";

    // parameter containers you want to save
    //

    MWriteRootFile write(outfname, "recreate");

    //    write.AddContainer("MSigmabar", "Events");

    write.AddContainer("MRawEvtHeader",  "Events");
    write.AddContainer("MMcEvt",         "Events", kFALSE);
    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("MPointingPos",   "Events");
    write.AddContainer("MTime",          "Events", kFALSE);
    write.AddContainer("MConcentration", "Events", kFALSE);

    tlist3.AddToList(&write);

    evtloop.SetProgressBar(bar);

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

    tlist3.PrintStatistics();

    //    bar->DestroyWindow();

    return;
}