source: trunk/MagicSoft/Mars/mjtrain/MJTrainSeparation.cc@ 7413

/* ======================================================================== *\
!
! *
! * 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 Bretz 11/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2005
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MJTrainSeparation
//
////////////////////////////////////////////////////////////////////////////
#include "MJTrainSeparation.h"

#include "MHMatrix.h"

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

// tools
#include "MDataSet.h"
#include "MTFillMatrix.h"
#include "MChisqEval.h"

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

// tasks
#include "MReadMarsFile.h"
#include "MContinue.h"
#include "MFillH.h"
#include "MRanForestCalc.h"
#include "MParameterCalc.h"

// container
#include "MMcEvt.hxx"
#include "MParameters.h"

// histograms
#include "MBinning.h"
#include "MH3.h"
#include "MHHadronness.h"

// filter
#include "MF.h"
#include "MFilterList.h"

ClassImp(MJTrainSeparation);

using namespace std;

Bool_t MJTrainSeparation::Train(const char *out)
{
    if (!fDataSetTrain.IsValid())
    {
        *fLog << err << "ERROR - DataSet for training invalid!" << endl;
        return kFALSE;
    }
    if (!fDataSetTest.IsValid())
    {
        *fLog << err << "ERROR - DataSet for testing invalid!" << endl;
        return kFALSE;
    }

    // --------------------- Setup files --------------------
    MReadMarsFile read("Events");
    MReadMarsFile read2("Events");
    MReadMarsFile read3("Events");
    read.DisableAutoScheme();
    read2.DisableAutoScheme();
    read3.DisableAutoScheme();

    fDataSetTrain.AddFilesOn(read);
    fDataSetTrain.AddFilesOff(read3);

    fDataSetTest.AddFilesOff(read2);
    fDataSetTest.AddFilesOn(read2);

    read2.VetoBranch("MMcTrig");
    read2.VetoBranch("MTime");
    read2.VetoBranch("MMcRunHeader");
    read2.VetoBranch("MMcTrigHeader");
    read2.VetoBranch("MMcFadcHeader");
    read2.VetoBranch("MMcCorsikaRunHeader");
    read2.VetoBranch("MMcConfigRunHeader");


    // ----------------------- Setup RF ----------------------
    MHMatrix train("Train");
    train.AddColumns(fRules);
    train.AddColumn("MHadronness.fVal");

    // ----------------------- Fill Matrix RF ----------------------

    MParameterD had("MHadronness");

    MParList plistx;
    plistx.AddToList(&had);

    MTFillMatrix fill;
    fill.SetLogStream(fLog);
    fill.SetDisplay(fDisplay);
    fill.AddPreCuts(fPreCuts);
    fill.AddPreCuts(fTrainCuts);

    // Set classifier for gammas
    had.SetVal(0);
    fill.SetName("FillGammas");
    fill.SetDestMatrix1(&train, fNumTrainOn);
    fill.SetReader(&read);
    if (!fill.Process(plistx))
        return kFALSE;

    // Set classifier for hadrons
    had.SetVal(1);
    fill.SetName("FillBackground");
    fill.SetDestMatrix1(&train, fNumTrainOff);
    fill.SetReader(&read3);
    if (!fill.Process(plistx))
        return kFALSE;

    // ------------------------ Train RF --------------------------

    MRanForestCalc rf;
    rf.SetNumObsoleteVariables(1);
    rf.SetDebug(fDebug);
    rf.SetDisplay(fDisplay);
    rf.SetLogStream(fLog);
    rf.SetFileName(out);
    rf.SetNameOutput("MHadronness");

    //MBinning b(2, -0.5, 1.5, "BinningHadronness", "lin");

    //if (!rf.TrainMultiRF(train, b.GetEdgesD()))    // classification
    //    return;

    //if (!rf.TrainSingleRF(train, b.GetEdgesD()))   // classification
    //    return;

    if (!rf.TrainSingleRF(train))                  // regression
        return kFALSE;

    // --------------------- Display result ----------------------

    gLog.Separator("Test");

    MParList  plist;
    MTaskList tlist;
    plist.AddToList(this);
    plist.AddToList(&tlist);

    MMcEvt mcevt;
    plist.AddToList(&mcevt);

    MBinning binsy(100, 0 , 1,      "BinningMH3Y", "lin");
    MBinning binsx(100, 10, 100000, "BinningMH3X", "log");

    plist.AddToList(&binsx);
    plist.AddToList(&binsy);

    MH3 h31("MHillas.fSize", "MHadronness.fVal");
    MH3 h32("MHillas.fSize", "MHadronness.fVal");
    h31.SetTitle("Background probability vs. Size:Size [phe]:Hadronness");
    h32.SetTitle("Background probability vs. Size:Size [phe]:Hadronness");

    MF f("MRawRunHeader.fRunType>255");
    MFilterList list;
    list.SetInverted();
    list.AddToList(&f);

    MHHadronness hist;

    MFillH fillh(&hist, "", "FillHadronness");
    MFillH fillh1(&h31, "", "FillGammas");
    MFillH fillh2(&h32, "", "FillBackground");
    fillh1.SetNameTab("Gammas");
    fillh2.SetNameTab("Background");

    fillh1.SetFilter(&f);
    fillh2.SetFilter(&list);

    MFilterList precuts;
    precuts.AddToList(fPreCuts);
    precuts.AddToList(fTestCuts);

    MContinue c0(&precuts);
    c0.SetInverted();

    tlist.AddToList(&read2);
    tlist.AddToList(&c0);
    tlist.AddToList(&rf);
    tlist.AddToList(&fillh);
    tlist.AddToList(&list);
    tlist.AddToList(&fillh1);
    tlist.AddToList(&fillh2);

    MEvtLoop loop;
    loop.SetDisplay(fDisplay);
    loop.SetLogStream(fLog);
    loop.SetParList(&plist);

    if (!loop.Eventloop())
        return kFALSE;

    if (!WriteDisplay(out))
        return kFALSE;

    return kTRUE;
}