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

/* ======================================================================== *\
!
! *
! * 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, 2006
!
!
\* ======================================================================== */

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

#include <TF1.h>
#include <TH2.h>
#include <TChain.h>
#include <TGraph.h>
#include <TVirtualPad.h>

#include "MHMatrix.h"

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

// tools
#include "MMath.h"
#include "MDataSet.h"
#include "MTFillMatrix.h"
#include "MChisqEval.h"
#include "MStatusDisplay.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 "MFEventSelector.h"
#include "MFilterList.h"

ClassImp(MJTrainSeparation);

using namespace std;

void MJTrainSeparation::DisplayResult(MH3 &h31, MH3 &h32)
{
    TH2 &g = (TH2&)h32.GetHist();
    TH2 &h = (TH2&)h31.GetHist();

    h.SetMarkerColor(kRed);
    g.SetMarkerColor(kGreen);

    const Int_t nx = h.GetNbinsX();
    const Int_t ny = h.GetNbinsY();

    gROOT->SetSelectedPad(NULL);

    TGraph gr1;
    TGraph gr2;
    for (int x=0; x<nx; x++)
    {
        TH1 *hx = h.ProjectionY("H_py", x+1, x+1);
        TH1 *gx = g.ProjectionY("G_py", x+1, x+1);

        Double_t max1 = -1;
        Double_t max2 = -1;
        Int_t maxy1 = 0;
        Int_t maxy2 = 0;
        for (int y=0; y<ny; y++)
        {
            const Float_t s = gx->Integral(1, y+1);
            const Float_t b = hx->Integral(1, y+1);
            const Float_t sig1 = MMath::SignificanceLiMa(s+b, b);
            const Float_t sig2 = s<1 ? 0 : MMath::SignificanceLiMa(s+b, b)*TMath::Log10(s);
            if (sig1>max1)
            {
                maxy1 = y;
                max1 = sig1;
            }
            if (sig2>max2)
            {
                maxy2 = y;
                max2 = sig2;
            }
        }

        gr1.SetPoint(x, h.GetXaxis()->GetBinCenter(x+1), h.GetYaxis()->GetBinCenter(maxy1+1));
        gr2.SetPoint(x, h.GetXaxis()->GetBinCenter(x+1), h.GetYaxis()->GetBinCenter(maxy2+1));

        delete hx;
        delete gx;
    }

    fDisplay->AddTab("OptCut");
    gPad->SetLogx();
    h.DrawCopy();
    g.DrawCopy("same");
    gr1.SetMarkerStyle(kFullDotMedium);
    gr1.DrawClone("LP")->SetBit(kCanDelete);
    gr2.SetLineColor(kBlue);
    gr2.SetMarkerStyle(kFullDotMedium);
    gr2.DrawClone("LP")->SetBit(kCanDelete);
}

/*
Bool_t  MJSpectrum::InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const
{
    fLog->Separator("Initialize energy weighting");

    if (!CheckEnv(w))
    {
        *fLog << err << "ERROR - Reading resources for MMcSpectrumWeight failed." << endl;
        return kFALSE;
    }

    TChain chain("RunHeaders");
    set.AddFilesOn(chain);

    MMcCorsikaRunHeader *h=0;
    chain.SetBranchAddress("MMcCorsikaRunHeader.", &h);
    chain.GetEntry(1);

    if (!h)
    {
        *fLog << err << "ERROR - Couldn't read MMcCorsikaRunHeader from DataSet." << endl;
        return kFALSE;
    }

    if (!w.Set(*h))
    {
        *fLog << err << "ERROR - Initializing MMcSpectrumWeight failed." << endl;
        return kFALSE;
    }

    w.Print();
    return kTRUE;
}

Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &weight) const
{
    // Some debug output
    fLog->Separator("Compiling original MC distribution");

    weight.SetNameMcEvt("MMcEvtBasic");
    const TString w(weight.GetFormulaWeights());
    weight.SetNameMcEvt();

    *fLog << inf << "Using weights: " << w << endl;
    *fLog << "Please stand by, this may take a while..." << flush;

    if (fDisplay)
        fDisplay->SetStatusLine1("Compiling MC distribution...");

    // Create chain
    TChain chain("OriginalMC");
    set.AddFilesOn(chain);

    // Prepare histogram
    h.Reset();

    // Fill histogram from chain
    h.SetDirectory(gROOT);
    if (h.InheritsFrom(TH2::Class()))
    {
        h.SetNameTitle("ThetaEMC", "Event-Distribution vs Theta and Energy for MC (produced)");
        h.SetXTitle("\\Theta [\\circ]");
        h.SetYTitle("E [GeV]");
        h.SetZTitle("Counts");
        chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", w, "goff");
    }
    else
    {
        h.SetNameTitle("ThetaMC", "Event-Distribution vs Theta for MC (produced)");
        h.SetXTitle("\\Theta [\\circ]");
        h.SetYTitle("Counts");
        chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", w, "goff");
    }
    h.SetDirectory(0);

    *fLog << "done." << endl;
    if (fDisplay)
        fDisplay->SetStatusLine2("done.");

    if (h.GetEntries()>0)
        return kTRUE;

    *fLog << err << "ERROR - Histogram with original MC distribution empty..." << endl;

    return h.GetEntries()>0;
}
*/

Bool_t MJTrainSeparation::GetEventsProduced(MDataSet &set, Double_t &num, Double_t &min, Double_t &max) const
{
    TChain chain("OriginalMC");
    set.AddFilesOn(chain);

    min = chain.GetMinimum("MMcEvtBasic.fEnergy");
    max = chain.GetMaximum("MMcEvtBasic.fEnergy");

    num = chain.GetEntries();

    if (num<100)
        *fLog << err << "ERROR - Less than 100 entries in OriginalMC-Tree of MC-Train-Data found." << endl;

    return num>=100;
}

Double_t MJTrainSeparation::GetDataRate(MDataSet &set) const
{
    TChain chain1("Events");
    set.AddFilesOff(chain1);

    const Double_t num = chain1.GetEntries();
    if (num<100)
    {
        *fLog << err << "ERROR - Less than 100 entries in Events-Tree of Train-Data found." << endl;
        return -1;
    }

    TChain chain("EffectiveOnTime");
    set.AddFilesOff(chain);

    TH1F h;
    h.SetDirectory(gROOT);
    h.SetNameTitle("OnTime", "Effective on-time");
    chain.Draw("MEffectiveOnTime.fVal>>OnTime", "", "goff");
    h.SetDirectory(0);

    if (h.Integral()<1)
    {
        *fLog << err << "ERROR - Less than 1s of effective observation time found in Train-Data." << endl;
        return -1;
    }

    *fLog << inf << "Found " << num << " events in " << h.Integral();
    *fLog << "s (" << num/h.Integral() << "Hz)" << endl;

    return num/h.Integral();
}

/*
 Scale:


 TF1 fold("old", "x^(-2.6)", emin, emax);
 TF1 fnew("new", "x^(-4.0)", emin, emax);

 TF1 q("q", "new/old", emin, emax);

 Double_t scale = 1./q.GetMaximum(emin, emax);

 // Anzahl produzierter Events vor MFEnergySlope:
 Double_t nold = fold.Integral(emin, emax);

 // Anzahl produzierter Events nach MFEnergySlope:
 Double_t nnew = fnew.Integral(emin, emax)*scale;

 class MFSpectrum : MMcSpectrumWeight
 {
 Double_t fScale;
 Bool_t   fResult;

 MFSpectrum::MFSpectrum(const char *name, const char *title)
 {
    fName  = name  ? name  : "MMcSpectrumWeight";
    fTitle = title ? title : "Task to calculate weights to change the energy spectrum"; 

    Init(fName, fTitle);

 }

 Int_t PreProcess(MParList *pList)
 {
     Int_t rc = MFSpectrumWeight::PreProcess(pList);
     if (rc!=kTRUE)
        return rc;

     fScale = fEval->GetMaximum(fEnergyMin, fEnergyMax);

     return kTRUE;
 }

 Int_t Process()
 {
    const Double_t e = fMcEvt->GetEnergy();

    Double_t prob = fFunc->Eval(e)/fScale;

    const Float_t Nexp = fN0 * pow(energy,fMcSlope-fNewSlope);
    const Float_t Nrnd = ;

    fResult = Nexp >= gRandom->Uniform();
 }

 }


 */

Bool_t MJTrainSeparation::AutoTrain()
{
    Double_t num, min, max;
    if (!GetEventsProduced(fDataSetTrain, num, min, max))
        return kFALSE;

    *fLog << inf << "Using build-in radius of 300m to calculate collection area!" << endl;

    // Target spectrum
    TF1 flx("Flux", "[0]*(x/1000)^(-2.6)", min, max);
    flx.SetParameter(0, 1e-5);

    // Number n0 of events this spectrum would produce per s and m^2
    const Double_t n0 = flx.Integral(min, max);    //[#]

    // Area produced in MC
    const Double_t A = TMath::Pi()*300*300;        //[m²]

    // Rate R of events this spectrum would produce per s
    const Double_t R = n0*A;                       //[Hz]

    // Number N of events produced (in trainings sample)
    const Double_t N = num;                        //[#]

    // This correponds to an observation time T [s]
    const Double_t T = N/R;                        //[s]

    // With an average data rate after star of
    const Double_t r = GetDataRate(fDataSetTrain); //[Hz]

    // this yields a number of n events to be read for training
    const Double_t n = r*T;                        //[#]

    if (r<0)
        return kFALSE;

    *fLog << "Calculated a total Monte Carlo observation time of " << T << "s" << endl;
    *fLog << "For a data rate of " << r << "Hz this corresponds to " << n << " data events." << endl;

    fNumTrainOn  = (UInt_t)-1;
    fNumTrainOff = TMath::Nint(n);

    /*
     An event rate dependent selection?
     ----------------------------------
     Total average data rate:      R
     Goal number of events:        N
     Number of data events:        N0
     Rate assigned to single evt:  r

     Selection probability: N/N0 * r/R

     f := N/N0 * r

     MF f("f * MEventRate.fRate < rand");
     */


    return kTRUE;
}

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;
    }

    // ----------------------- Auto Train? ----------------------

    if (fAutoTrain)
        if (!AutoTrain())
            return kFALSE;

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

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

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

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

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

    MParameterD had("MHadronness");

    MParList plistx;
    plistx.AddToList(&had);
    plistx.AddToList(this);

    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(&read1);
    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.SetNumTrees(fNumTrees);
    rf.SetNdSize(fNdSize);
    rf.SetNumTry(fNumTry);
    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;

    //fDisplay = rf.GetDisplay();

    // --------------------- Display result ----------------------
    gLog.Separator("Test");

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

    MMcEvt mcevt;
    plist.AddToList(&mcevt);

    // ----- Setup histograms -----
    MBinning binsy(100, 0 , 1,      "BinningMH3Y", "lin");
    MBinning binsx( 50, 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");

    MHHadronness hist;

    // ----- Setup tasks -----
    MFillH fillh0(&hist, "", "FillHadronness");
    MFillH fillh1(&h31);
    MFillH fillh2(&h32);
    fillh1.SetNameTab("Background");
    fillh2.SetNameTab("Gammas");
    fillh0.SetBit(MFillH::kDoNotDisplay);

    // ----- Setup filter -----
    MFilterList precuts;
    precuts.AddToList(fPreCuts);
    precuts.AddToList(fTestCuts);

    MContinue c0(&precuts);
    c0.SetName("PreCuts");
    c0.SetInverted();

    MFEventSelector sel;
    sel.SetNumSelectEvts(fNumTestOn);

    MContinue c1(&sel);
    c1.SetInverted();

    // ----- Setup tasklist -----
    tlist.AddToList(&read2);
    tlist.AddToList(&c0);
    tlist.AddToList(&c1);
    tlist.AddToList(&rf);
    tlist.AddToList(&fillh0);
    tlist.AddToList(&fillh1);

    // ----- Run eventloop on gammas -----
    MEvtLoop loop;
    loop.SetDisplay(fDisplay);
    loop.SetLogStream(fLog);
    loop.SetParList(&plist);

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

    // ----- Setup and run eventloop on background -----
    sel.SetNumSelectEvts(fNumTestOff);
    fillh0.ResetBit(MFillH::kDoNotDisplay);

    tlist.Replace(&read4);
    tlist.Replace(&fillh2);

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

    DisplayResult(h31, h32);

    if (!WriteDisplay(out))
        return kFALSE;

    return kTRUE;
}