Index: trunk/Mars/mjtrain/MJTrainCuts.cc =================================================================== --- trunk/Mars/mjtrain/MJTrainCuts.cc (revision 9870) +++ trunk/Mars/mjtrain/MJTrainCuts.cc (revision 9870) @@ -0,0 +1,1033 @@ +/* ======================================================================== *\ +! +! * +! * 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 08/2010 +! +! Copyright: MAGIC Software Development, 2010 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MJTrainCuts +// ========= +// +// This class is meant as a tool to understand better what a trained +// random forest is doing in the multi-dimensional phase space. +// Consequently, it can also be used to deduce good one or two dimensional +// cuts from the results by mimicing the behaviour of the random forest. +// +// +// Usage +// ----- +// +// The instance is created by its default constructor +// +// MJTrainCuts opt; +// +// In a first step a random forest must be trained and in a second step its +// performance can be evaluated with an independent test sample. The used +// samples are defined by two MDataSet objects, one for the on-data (e.g. +// gammas) and the other one for the off-data (e.g. protons). SequencesOn +// and SequencesOff are used for testing and training respectively. +// +// MDataSet seton ("myondata.txt"); +// MDataSet setoff("myoffdata.txt"); +// +// // If you want to use all available events +// opt.SetDataSetOn(seton); +// opt.SetDataSetOff(setoff); +// +// // Try to select 10000 and 30000 events for training and testing resp. +// // opt.SetDataSetOn(seton, 10000, 30000); +// // opt.SetDataSetOff(setoff, 10000, 30000); +// +// Note that by using several data set in one file (see MDataSet) you can +// have everything in a single file. +// +// The variables which are used for training are now setup as usual +// +// Int_t p1 = opt.AddParameter("MHillas.fSize"); +// Int_t p2 = opt.AddParameter("MHillas.GetArea*MGeomCam.fConvMm2Deg^2"); +// Int_t p3 = opt.AddParameter("MHillasSrc.fDist*MGeomCam.fConvMm2Deg"); +// +// In addition you can now setup a binning for the display of each train +// parameter as follows (for details see MBinning) +// +// opt.AddBinning(p1, MBinning(40, 10, 10000, "", "log")); +// opt.AddBinning(p2, MBinning(50, 0, 0.25)); +// opt.AddBinning(p3, MBinning(50, 0, 2.5)); +// +// Since with increasing number of variables the possibly combinations +// increase to fast you have to define which plots you are interested in, +// for example: +// +// opt.AddHist(p3); // A 1D plot dist +// opt.AddHist(p1, p2); // A 2D plot area vs. size +// opt.AddHist(p3, p2); // A 2D plot dist vs. size +// +// Also 3D plots are avaiable but they are most probably difficult to +// interprete. +// +// In addition to this you have the usual user interface, i.e. that +// - PreCuts +// - TrainCuts +// - TestCuts +// - PreTasks +// - PostTasks +// - TestTasks +// are available. For details see MJOptimizeBase +// +// void EnableRegression() / void EnableClassification() +// Defines whether to use the random forest's regression of +// classification method. Classification is the default. +// +// +// The produced plots +// ------------------ +// +// The tab with the plots filled will always look like this: +// +// +--------+--------+ +// |1 |2 | +// +--------+--------| +// |3 |4 | +// +--------+ | +// |5 | | +// +--------+--------+ +// +// Pad1 and Pad2 contain the weighted event distribution of the test-sample. +// +// Pad2 and Pad5 conatin a profile of the hadronness distribution of the +// test-sample. +// +// Pad4 contains a profile of the hadronness distribution of on-data and +// off-data together of the test-data. +// +// If the profiles for on-data and off-data are identical the displayed +// hadronness is obviously independant of the other (non shown) trainings +// variables. Therefore the difference between the two plots show how much +// the variables are correlated. The same is true if the prfiles in +// pad3 and pad5 don't differe from the profile in pad4. +// +// In the most simple case - the random forest is only trained with the +// variables displayed - all three plots should be identical (apart from the +// difference in the distrubution of the three sets). +// +// The plot in pad4 can now be used to deduce a good classical cut in the +// displayed variables. +// +//////////////////////////////////////////////////////////////////////////// +#include "MJTrainCuts.h" + +#include +#include +#include +#include +#include + +#include "MHMatrix.h" + +#include "MLog.h" +#include "MLogManip.h" + +// tools +#include "MMath.h" +#include "MBinning.h" +#include "MTFillMatrix.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" + +// container +#include "MParameters.h" + +// histograms +#include "MHn.h" +#include "MHHadronness.h" + +// filter +#include "MFEventSelector.h" +#include "MFilterList.h" + +using namespace std; + +class HistSet1D : public TObject +{ +protected: + UInt_t fNx; + + virtual void AddHist(MHn &h, const char *rx, const char *, const char *) const + { + h.AddHist(rx); + } + virtual void AddProf(MHn &h, const char *rx, const char *, const char *) const + { + h.AddHist(rx, "MHadronness.fVal", MH3::kProfileSpread); + } + virtual void SetupName(MHn &h, const char *name) const + { + h.InitName(Form("%s%d;%d", name, fNx, fNx)); + } + void SetupHist(MHn &h, const char *name, const char *title) const + { + SetupName(h, name); + + h.SetAutoRange(kFALSE, kFALSE, kFALSE); + h.InitTitle(title); + h.SetDrawOption("colz"); + } + + void CreateHist(MHn &h, const char *rx, const char *ry=0, const char *rz=0) const + { + h.SetLayout(MHn::kComplex); + h.SetBit(MHn::kDoNotReset); + + AddHist(h, rx, ry, rz); + SetupHist(h, "DistOn", "Distribution of on-data"); + h.SetWeight("Type.fVal"); + + AddHist(h, rx, ry, rz); + SetupHist(h, "DistOff", "Distribution of off-data"); + h.SetWeight("1-Type.fVal"); + + AddProf(h, rx, ry, rz); + SetupHist(h, "HadOn", "Hadronness profile for on-data"); + h.SetWeight("Type.fVal"); + + AddProf(h, rx, ry, rz); + SetupHist(h, "Had", "Hadronness profile for all events"); + + AddProf(h, rx, ry, rz); + SetupHist(h, "HadOff", "Hadronness profile for off-data"); + h.SetWeight("1-Type.fVal"); + } + + +public: + HistSet1D(UInt_t nx) : fNx(nx) { } + + virtual MHn *GetHistN(const TList &rules) const + { + if (!rules.At(fNx)) + return 0; + + MHn *h = new MHn(Form("%d", fNx)); + + CreateHist(*h, rules.At(fNx)->GetName()); + + return h; + } + + virtual Bool_t CheckBinning(const TObjArray &binnings) const + { + return binnings.FindObject(Form("Binning%d", fNx)); + } +}; + +class HistSet2D : public HistSet1D +{ +protected: + UInt_t fNy; + + void AddHist(MHn &h, const char *rx, const char *ry, const char *) const + { + h.AddHist(rx, ry); + } + void AddProf(MHn &h, const char *rx, const char *ry, const char *) const + { + h.AddHist(rx, ry, "MHadronness.fVal", MH3::kProfileSpread); + } + void SetupName(MHn &h, const char *name) const + { + h.InitName(Form("%s%d:%d;%d;%d", name, fNx, fNy, fNx, fNy)); + } + +public: + HistSet2D(UInt_t nx, UInt_t ny) : HistSet1D(nx), fNy(ny) { } + + MHn *GetHistN(const TList &rules) const + { + if (!rules.At(fNx) || !rules.At(fNy)) + return 0; + + MHn *h = new MHn(Form("%d:%d", fNx, fNy)); + + CreateHist(*h, rules.At(fNx)->GetName(), rules.At(fNy)->GetName()); + + return h; + } + Bool_t CheckBinning(const TObjArray &binnings) const + { + return HistSet1D::CheckBinning(binnings) && binnings.FindObject(Form("Binning%d", fNy)); + } +}; + + +class HistSet3D : public HistSet2D +{ +private: + UInt_t fNz; + + void AddHist(MHn &h, const char *rx, const char *ry, const char *rz) const + { + h.AddHist(rx, ry, rz); + } + void AddProf(MHn &h, const char *rx, const char *ry, const char *rz) const + { + h.AddHist(rx, ry, rz, "MHadronness.fVal"); + } + void SetupName(MHn &h, const char *name) const + { + h.InitName(Form("%s%d:%d:%d;%d;%d;%d", name, fNx, fNy, fNz, fNx, fNy, fNz)); + } + +public: + HistSet3D(UInt_t nx, UInt_t ny, UInt_t nz) : HistSet2D(nx, ny), fNz(nz) { } + + MHn *GetHistN(const TList &rules) const + { + if (!rules.At(fNx) || !rules.At(fNy) || !rules.At(fNz)) + return 0; + + MHn *h = new MHn(Form("%d:%d:%d", fNx, fNy, fNz)); + + CreateHist(*h, + rules.At(fNx)->GetName(), + rules.At(fNy)->GetName(), + rules.At(fNy)->GetName()); + + return h; + } + Bool_t CheckBinning(const TObjArray &binnings) const + { + return HistSet2D::CheckBinning(binnings) && binnings.FindObject(Form("Binning%d", fNz)); + } +}; + +// --------------------------------------------------------------------------------------- + +void MJTrainCuts::AddHist(UInt_t nx) +{ + fHists.Add(new HistSet1D(nx)); +} + +void MJTrainCuts::AddHist(UInt_t nx, UInt_t ny) +{ + fHists.Add(new HistSet2D(nx, ny)); +} + +void MJTrainCuts::AddHist(UInt_t nx, UInt_t ny, UInt_t nz) +{ + fHists.Add(new HistSet3D(nx, ny, nz)); +} + +void MJTrainCuts::AddBinning(UInt_t n, const MBinning &bins) +{ + const char *name = Form("Binning%d", n); + + TObject *o = fBinnings.FindObject(name); + if (o) + { + delete fBinnings.Remove(o); + *fLog << warn << "WARNING - Binning for parameter " << n << " (" << name << ") already exists... replaced." << endl; + } + + // FIXME: Check for existence + fBinnings.Add(new MBinning(bins, name, bins.GetTitle())); +} + +/* +void MJTrainCuts::AddBinning(const MBinning &bins) +{ + // FIXME: Check for existence + fBinnings.Add(new MBinning(bins, bins.GetName(), bins.GetTitle())); +} +*/ + +// --------------------------------------------------------------------------------------- + +// -------------------------------------------------------------------------- +// +void MJTrainCuts::DisplayResult(MH3 &h31, MH3 &h32, Float_t ontime) +{ + TH2D &g = (TH2D&)h32.GetHist(); + TH2D &h = (TH2D&)h31.GetHist(); + + h.SetMarkerColor(kRed); + g.SetMarkerColor(kBlue); + + TH2D res1(g); + TH2D res2(g); + + h.SetTitle("Hadronness-Distribution vs. Size"); + res1.SetTitle("Significance Li/Ma"); + res1.SetXTitle("Size [phe]"); + res1.SetYTitle("Hadronness"); + res2.SetTitle("Significance-Distribution"); + res2.SetXTitle("Size-Cut [phe]"); + res2.SetYTitle("Hadronness-Cut"); + res1.SetContour(50); + res2.SetContour(50); + + const Int_t nx = h.GetNbinsX(); + const Int_t ny = h.GetNbinsY(); + + gROOT->SetSelectedPad(NULL); + + + Double_t Stot = 0; + Double_t Btot = 0; + + Double_t max2 = -1; + + TGraph gr1; + TGraph gr2; + for (int x=nx-1; x>=0; x--) + { + TH1 *hx = h.ProjectionY("H_py", x+1, x+1); + TH1 *gx = g.ProjectionY("G_py", x+1, x+1); + + Double_t S = 0; + Double_t B = 0; + + Double_t max1 = -1; + Int_t maxy1 = 0; + Int_t maxy2 = 0; + for (int y=ny-1; y>=0; 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 = MMath::SignificanceLiMa(s+Stot+b+Btot, b+Btot)*TMath::Log10(s+Stot+1); + if (sig1>max1) + { + maxy1 = y; + max1 = sig1; + } + if (sig2>max2) + { + maxy2 = y; + max2 = sig2; + + S=s; + B=b; + } + + res1.SetBinContent(x+1, y+1, sig1); + } + + Stot += S; + Btot += B; + + 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; + } + + //cout << "--> " << MMath::SignificanceLiMa(Stot+Btot, Btot) << " "; + //cout << Stot << " " << Btot << endl; + + + Int_t mx1=0; + Int_t my1=0; + Int_t mx2=0; + Int_t my2=0; + Int_t s1=0; + Int_t b1=0; + Int_t s2=0; + Int_t b2=0; + Double_t sig1=-1; + Double_t sig2=-1; + for (int x=0; xIntegral(1, y+1); + const Float_t b = hx->Integral(1, y+1); + const Float_t sig = MMath::SignificanceLiMa(s+b, b); + res2.SetBinContent(x+1, y+1, sig); + + // Search for top-rightmost maximum + if (sig>=sig1) + { + mx1=x+1; + my1=y+1; + s1 = TMath::Nint(s); + b1 = TMath::Nint(b); + sig1=sig; + } + if (TMath::Log10(s)*sig>=sig2) + { + mx2=x+1; + my2=y+1; + s2 = TMath::Nint(s); + b2 = TMath::Nint(b); + sig2=TMath::Log10(s)*sig; + } + } + delete hx; + delete gx; + } + + TGraph gr3; + TGraph gr4; + gr4.SetTitle("Significance Li/Ma vs. Hadronness-cut"); + + TH1 *hx = h.ProjectionY("H_py"); + TH1 *gx = g.ProjectionY("G_py"); + for (int y=0; yIntegral(1, y+1); + const Float_t b = hx->Integral(1, y+1); + const Float_t sg1 = MMath::SignificanceLiMa(s+b, b); + const Float_t sg2 = s<1 ? 0 : MMath::SignificanceLiMa(s+b, b)*TMath::Log10(s); + + gr3.SetPoint(y, h.GetYaxis()->GetBinLowEdge(y+2), sg1); + gr4.SetPoint(y, h.GetYaxis()->GetBinLowEdge(y+2), sg2); + } + delete hx; + delete gx; + + if (fDisplay) + { + TCanvas &c = fDisplay->AddTab("OptCut"); + c.SetBorderMode(0); + c.Divide(2,2); + + gROOT->SetSelectedPad(0); + c.cd(1); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + gPad->SetLogx(); + gPad->SetGridx(); + gPad->SetGridy(); + h.DrawCopy(); + g.DrawCopy("same"); + gr1.SetMarkerStyle(kFullDotMedium); + gr1.DrawClone("LP")->SetBit(kCanDelete); + gr2.SetLineColor(kBlue); + gr2.SetMarkerStyle(kFullDotMedium); + gr2.DrawClone("LP")->SetBit(kCanDelete); + + gROOT->SetSelectedPad(0); + c.cd(3); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + gPad->SetGridx(); + gPad->SetGridy(); + gr4.SetMinimum(0); + gr4.SetMarkerStyle(kFullDotMedium); + gr4.DrawClone("ALP")->SetBit(kCanDelete); + gr3.SetLineColor(kBlue); + gr3.SetMarkerStyle(kFullDotMedium); + gr3.DrawClone("LP")->SetBit(kCanDelete); + + c.cd(2); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + gPad->SetLogx(); + gPad->SetGridx(); + gPad->SetGridy(); + gPad->AddExec("color", "gStyle->SetPalette(1, 0);"); + res1.SetMaximum(7); + res1.DrawCopy("colz"); + + c.cd(4); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + gPad->SetLogx(); + gPad->SetGridx(); + gPad->SetGridy(); + gPad->AddExec("color", "gStyle->SetPalette(1, 0);"); + res2.SetMaximum(res2.GetMaximum()*1.05); + res2.DrawCopy("colz"); + + // Int_t mx, my, mz; + // res2.GetMaximumBin(mx, my, mz); + + TMarker m; + m.SetMarkerStyle(kStar); + m.DrawMarker(res2.GetXaxis()->GetBinCenter(mx1), res2.GetYaxis()->GetBinCenter(my1)); + m.SetMarkerStyle(kPlus); + m.DrawMarker(res2.GetXaxis()->GetBinCenter(mx2), res2.GetYaxis()->GetBinCenter(my2)); + } + + if (ontime>0) + *fLog << all << "Observation Time: " << TMath::Nint(ontime/60) << "min" << endl; + *fLog << "Maximum Significance: " << Form("%.1f", sig1); + if (ontime>0) + *fLog << Form(" [%.1f/sqrt(h)]", sig1/TMath::Sqrt(ontime/3600)); + *fLog << endl; + + *fLog << "Significance: S=" << Form("%.1f", sig1) << " E=" << s1 << " B=" << b1 << " h<"; + *fLog << Form("%.2f", res2.GetYaxis()->GetBinCenter(my1)) << " s>"; + *fLog << Form("%3d", TMath::Nint(res2.GetXaxis()->GetBinCenter(mx1))) << endl; + *fLog << "Significance*LogE: S=" << Form("%.1f", sig2/TMath::Log10(s2)) << " E=" << s2 << " B=" << b2 << " h<"; + *fLog << Form("%.2f", res2.GetYaxis()->GetBinCenter(my2)) << " s>"; + *fLog << Form("%3d", TMath::Nint(res2.GetXaxis()->GetBinCenter(mx2))) << endl; + *fLog << endl; +} + +// -------------------------------------------------------------------------- +// +Bool_t MJTrainCuts::Process(const char *out) +{ + // =========================== Consistency checks ================================== + if (!fDataSetOn.IsValid()) + { + *fLog << err << "ERROR - DataSet for on-data invalid!" << endl; + return kFALSE; + } + if (!fDataSetOff.IsValid()) + { + *fLog << err << "ERROR - DataSet for off-data invalid!" << endl; + return kFALSE; + } + + if (fDataSetOn.IsWobbleMode()!=fDataSetOff.IsWobbleMode()) + { + *fLog << err << "ERROR - On- and Off-DataSet have different observation modes!" << endl; + return kFALSE; + } + + if (fDataSetOn.IsMonteCarlo()!=fDataSetOff.IsMonteCarlo()) + { + *fLog << err << "ERROR - On- and Off-DataSet have different monte carlo modes!" << endl; + return kFALSE; + } + + if (!HasWritePermission(out)) + return kFALSE; + + // Check if needed binning exists + TIter NextH(&fHists); + TObject *o = 0; + while ((o=NextH())) + { + const HistSet1D *hs = static_cast(o); + if (hs->CheckBinning(fBinnings)) + continue; + + *fLog << err << "ERROR - Not all needed binnning exist." << endl; + return kFALSE; + } + + // =========================== Preparation ================================== + + if (fDisplay) + fDisplay->SetTitle(out); + + TStopwatch clock; + clock.Start(); + + // ------------------ Setup reading -------------------- + MReadMarsFile read1("Events"); + MReadMarsFile read2("Events"); + MReadMarsFile read3("Events"); + MReadMarsFile read4("Events"); + read1.DisableAutoScheme(); + read2.DisableAutoScheme(); + read3.DisableAutoScheme(); + read4.DisableAutoScheme(); + + // Setup four reading tasks with the on- and off-data of the two datasets + // Training -- On + if (!fDataSetOn.AddFilesOn(read1)) + return kFALSE; + // Testing -- On + if (!fDataSetOn.AddFilesOff(read4)) + return kFALSE; + // Training -- Off + if (!fDataSetOff.AddFilesOn(read3)) + return kFALSE; + // Testing -- Off + if (!fDataSetOff.AddFilesOff(read2)) + return kFALSE; + + // =============================================================================== + // ====================== Training ========================= + // =============================================================================== + + // ---------------- Setup RF Matrix ---------------- + MHMatrix train("Train"); + train.AddColumns(fRules); +// if (fEnableWeights[kTrainOn] || fEnableWeights[kTrainOff]) +// train.AddColumn("MWeight.fVal"); + train.AddColumn("MHadronness.fVal"); + + // ----------------- Prepare filling Matrix RF ------------------ + + // Setup the hadronness container identifying gammas and off-data + // and setup a container for the weights + MParameterD had("MHadronness"); + MParameterD wgt("MWeight"); + MParameterD typ("Type"); + + // Add them to the parameter list + MParList plistx; + plistx.AddToList(this); // take care of fDisplay! + plistx.AddToList(&had); + plistx.AddToList(&wgt); + plistx.AddToList(&typ); + + // Setup the tool class to fill the matrix + MTFillMatrix fill; + fill.SetLogStream(fLog); + fill.SetDisplay(fDisplay); + fill.AddPreCuts(fPreCuts); + fill.AddPreCuts(fTrainCuts); + + // ----------------- Fill on data into matrix ------------------ + + // Setup the tool class to read the gammas and read them + fill.SetName("FillOn"); + fill.SetDestMatrix1(&train, fNum[kTrainOn]); + fill.SetReader(&read1); +// fill.AddPreTasks(fPreTasksSet[kTrainOn]); + fill.AddPreTasks(fPreTasks); +// fill.AddPostTasks(fPostTasksSet[kTrainOn]); + fill.AddPostTasks(fPostTasks); + + // Set classifier for gammas + had.SetVal(0); + wgt.SetVal(1); + typ.SetVal(0); + + // Fill matrix + if (!fill.Process(plistx)) + return kFALSE; + + // Check the number or read events + const Int_t numontrn = train.GetNumRows(); + if (numontrn==0) + { + *fLog << err << "ERROR - No on-data events available for training... aborting." << endl; + return kFALSE; + } + + // Remove possible post tasks + fill.ClearPreTasks(); + fill.ClearPostTasks(); + + // ----------------- Fill off data into matrix ------------------ + + // In case of wobble mode we have to do something special + // Setup the tool class to read the background and read them + fill.SetName("FillOff"); + fill.SetDestMatrix1(&train, fNum[kTrainOff]); + fill.SetReader(&read3); +// fill.AddPreTasks(fPreTasksSet[kTrainOff]); + fill.AddPreTasks(fPreTasks); +// fill.AddPostTasks(fPostTasksSet[kTrainOff]); + fill.AddPostTasks(fPostTasks); + + // Set classifier for background + had.SetVal(1); + wgt.SetVal(1); + typ.SetVal(1); + + // Fiull matrix + if (!fill.Process(plistx)) + return kFALSE; + + // Check the number or read events + const Int_t numofftrn = train.GetNumRows()-numontrn; + if (numofftrn==0) + { + *fLog << err << "ERROR - No off-data available for training... aborting." << endl; + return kFALSE; + } + + // ------------------------ Train RF -------------------------- + + MRanForestCalc rf("TrainSeparation", fTitle); + rf.SetNumTrees(fNumTrees); + rf.SetNdSize(fNdSize); + rf.SetNumTry(fNumTry); + rf.SetNumObsoleteVariables(1); +// rf.SetLastDataColumnHasWeights(fEnableWeights[kTrainOn] || fEnableWeights[kTrainOff]); + rf.SetDebug(fDebug>1); + rf.SetDisplay(fDisplay); + rf.SetLogStream(fLog); + rf.SetFileName(out); + rf.SetNameOutput("MHadronness"); + + // Train the random forest either by classification or regression + if (!rf.Train(train, fUseRegression)) + return kFALSE; + + // ----------------- Print result of training ------------------ + + // Output information about what was going on so far. + *fLog << all; + fLog->Separator("The forest was trained with..."); + + *fLog << "Training method:" << endl; + *fLog << " * " << (fUseRegression?"regression":"classification") << endl; + /* + if (fEnableWeights[kTrainOn]) + *fLog << " * weights for on-data" << endl; + if (fEnableWeights[kTrainOff]) + *fLog << " * weights for off-data" << endl; + */ + *fLog << endl; + *fLog << "Events used for training:" << endl; + *fLog << " * Gammas: " << numontrn << endl; + *fLog << " * Background: " << numofftrn << endl; + *fLog << endl; + *fLog << "Gamma/Background ratio:" << endl; + *fLog << " * Requested: " << (float)fNum[kTrainOn]/fNum[kTrainOff] << endl; + *fLog << " * Result: " << (float)numontrn/numofftrn << endl; + *fLog << endl; + *fLog << "Run-Time: " << Form("%.1f", clock.RealTime()/60) << "min (CPU: "; + *fLog << Form("%.1f", clock.CpuTime()/60) << "min)" << endl; + *fLog << endl; + *fLog << "Output file name: " << out << endl; + + // =============================================================================== + // ====================== Testing ========================= + // =============================================================================== + fLog->Separator("Test"); + + clock.Continue(); + + // ---------------------- Prepare eventloop off-data --------------------- + + // Setup parlist and tasklist for testing + MParList plist; + MTaskList tlist; + plist.AddToList(this); // Take care of display + plist.AddToList(&tlist); + +// MMcEvt mcevt; +// plist.AddToList(&mcevt); + + plist.AddToList(&wgt); + plist.AddToList(&typ); + + // ----- Setup histograms ----- + MBinning binsy(50, 0 , 1, "BinningMH3Y", "lin"); + MBinning binsx(40, 10, 100000, "BinningMH3X", "log"); + + plist.AddToList(&binsx); + plist.AddToList(&binsy); + + MH3 h31("MHillas.fSize", "MHadronness.fVal"); + MH3 h32("MHillas.fSize", "MHadronness.fVal"); + MH3 h40("MMcEvt.fEnergy", "MHadronness.fVal"); + h31.SetTitle("Background probability vs. Size:Size [phe]:Hadronness h"); + h32.SetTitle("Background probability vs. Size:Size [phe]:Hadronness h"); + h40.SetTitle("Background probability vs. Energy:Energy [GeV]:Hadronness h"); + + plist.AddToList(&fBinnings); + + MHHadronness hist; + + // ----- Setup tasks ----- + MFillH fillh0(&hist, "", "FillHadronness"); + MFillH fillh1(&h31, "", "FillHadVsSize"); + MFillH fillh2(&h32, "", "FillHadVsSize"); + MFillH fillh4(&h40, "", "FillHadVsEnergy"); + fillh0.SetWeight("MWeight"); + fillh1.SetWeight("MWeight"); + fillh2.SetWeight("MWeight"); + fillh4.SetWeight("MWeight"); + fillh1.SetDrawOption("colz profy"); + fillh2.SetDrawOption("colz profy"); + fillh4.SetDrawOption("colz profy"); + fillh1.SetNameTab("HadSzOff"); + fillh2.SetNameTab("HadSzOn"); + fillh4.SetNameTab("HadEnOn"); + fillh0.SetBit(MFillH::kDoNotDisplay); + + // ----- Setup filter ----- + MFilterList precuts; + precuts.AddToList(fPreCuts); + precuts.AddToList(fTestCuts); + + MContinue cont0(&precuts); + cont0.SetName("PreCuts"); + cont0.SetInverted(); + + MFEventSelector sel; // FIXME: USING IT (WITH PROB?) in READ will by much faster!!! + sel.SetNumSelectEvts(fNum[kTestOff]); + + MContinue contsel(&sel); + contsel.SetInverted(); + + // ----- Setup tasklist ----- + tlist.AddToList(&read2); // Reading task + tlist.AddToList(&contsel); // event selector +// tlist.AddToList(fPreTasksSet[kTestOff]); + tlist.AddToList(fPreTasks); // list of pre tasks + tlist.AddToList(&cont0); // list of pre cuts and test cuts + tlist.AddToList(&rf); // evaluate random forest +// tlist.AddToList(fPostTasksSet[kTestOff]); + tlist.AddToList(fPostTasks); // list of post tasks + tlist.AddToList(&fillh1); // Fill HadSzOff + + TList autodel; + autodel.SetOwner(); + + TPRegexp regexp("([0-9]:*)+"); + + NextH.Reset(); + while ((o=NextH())) + { + HistSet1D *hs = static_cast(o); + + // FIXME: Move to beginning of function + // Check if needed binning exists + if (!hs->CheckBinning(fBinnings)) + return kFALSE; + + MHn *hist = hs->GetHistN(fRules); + MFillH *fill = new MFillH(hist, "", Form("Fill%s", hist->GetName())); + + fill->SetWeight("MWeight"); + fill->SetDrawOption("colz"); + fill->SetNameTab(hist->GetName()); + fill->SetBit(MFillH::kDoNotDisplay); + + tlist.AddToList(fill); + + autodel.Add(fill); + autodel.Add(hist); + } + tlist.AddToList(&fillh0); // Fill MHHadronness (not displayed in first loop) + tlist.AddToList(&fTestTasks); // list of test tasks + + // Enable Acceleration + tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime); + + // ---------------------- Run eventloop on background --------------------- + MEvtLoop loop; + loop.SetDisplay(fDisplay); + loop.SetLogStream(fLog); + loop.SetParList(&plist); + //if (!SetupEnv(loop)) + // return kFALSE; + + wgt.SetVal(1); + typ.SetVal(0); + if (!loop.Eventloop()) + return kFALSE; + + // ---------------------- Prepare eventloop on-data --------------------- + + sel.SetNumSelectEvts(fNum[kTestOn]); // set number of target events + + fillh0.ResetBit(MFillH::kDoNotDisplay); // Switch on display MHHadronness + + TIter NextF(&autodel); + while ((o=NextF())) + { + MFillH *fill = dynamic_cast(o); + if (fill) + fill->ResetBit(MFillH::kDoNotDisplay); + } + + // Remove PreTasksOff and PostTasksOff from the list +// tlist.RemoveFromList(fPreTasksSet[kTestOff]); +// tlist.RemoveFromList(fPostTasksSet[kTestOff]); + + tlist.Replace(&read4); // replace reading off-data by on-data + + // Add the PreTasksOn directly after the reading task +// tlist.AddToListAfter(fPreTasksSet[kTestOn], &c1); + + // Add the PostTasksOn after rf +// tlist.AddToListAfter(fPostTasksSet[kTestOn], &rf); + + tlist.Replace(&fillh2); // Fill HadSzOn instead of HadSzOff + tlist.AddToListAfter(&fillh4, &fillh0); // Filling of HadEnOn + + // Enable Acceleration + tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime); + + // ---------------------- Run eventloop on-data --------------------- + + wgt.SetVal(1); + typ.SetVal(1); + if (!loop.Eventloop()) + return kFALSE; + + // ---------------------- Print/Display result --------------------- + + // Show what was going on in the testing + const Double_t numontst = h32.GetHist().GetEntries(); + const Double_t numofftst = h31.GetHist().GetEntries(); + + *fLog << all; + fLog->Separator("The forest was tested with..."); + *fLog << "Test method:" << endl; + *fLog << " * Random Forest: " << out << endl; + /* + if (fEnableWeights[kTestOn]) + *fLog << " * weights for on-data" << endl; + if (fEnableWeights[kTestOff]) + *fLog << " * weights for off-data" << endl; + */ + *fLog << endl; + *fLog << "Events used for test:" << endl; + *fLog << " * Gammas: " << numontst << endl; + *fLog << " * Background: " << numofftst << endl; + *fLog << endl; + *fLog << "Gamma/Background ratio:" << endl; + *fLog << " * Requested: " << (float)fNum[kTestOn]/fNum[kTestOff] << endl; + *fLog << " * Result: " << (float)numontst/numofftst << endl; + *fLog << endl; + + // Display the result plots + DisplayResult(h31, h32, -1); + //DisplayResult(h31, h32, ontime); + + *fLog << "Total Run-Time: " << Form("%.1f", clock.RealTime()/60) << "min (CPU: "; + *fLog << Form("%.1f", clock.CpuTime()/60) << "min)" << endl; + fLog->Separator(); + + // ----------------- Write result ------------------ + + fDataSetOn.SetName("DataSetOn"); + fDataSetOff.SetName("DataSetOff"); + + // Write the display + TObjArray arr; + arr.Add(const_cast(&fDataSetOn)); + arr.Add(const_cast(&fDataSetOff)); + if (fDisplay) + arr.Add(fDisplay); + + SetPathOut(out); + return WriteContainer(arr, 0, "UPDATE"); +} + Index: trunk/Mars/mjtrain/MJTrainCuts.h =================================================================== --- trunk/Mars/mjtrain/MJTrainCuts.h (revision 9870) +++ trunk/Mars/mjtrain/MJTrainCuts.h (revision 9870) @@ -0,0 +1,117 @@ +#ifndef MARS_MJTrainCuts +#define MARS_MJTrainCuts + +#ifndef MARS_MJTrainRanForest +#include "MJTrainRanForest.h" +#endif + +#ifndef ROOT_TObjArray +#include +#endif + +#ifndef MARS_MDataSet +#include "MDataSet.h" +#endif + +class MH3; +class MBinning; + +class MJTrainCuts : public MJTrainRanForest +{ +public: + enum Type_t { kTrainOn, kTrainOff, kTestOn, kTestOff }; + +private: + MDataSet fDataSetOn; + MDataSet fDataSetOff; + + UInt_t fNum[4]; +/* + TList fPreTasksSet[4]; + TList fPostTasksSet[4]; + + Bool_t fEnableWeights[4]; +*/ + Bool_t fUseRegression; + + TList fHists; + TObjArray fBinnings; + + // Result + void DisplayResult(MH3 &h31, MH3 &h32, Float_t ontime); + +public: + MJTrainCuts() : fUseRegression(kFALSE) + { + for (int i=0; i<4; i++) + { + //fEnableWeights[i]=kFALSE; + fNum[i] = (UInt_t)-1; + } + + fHists.SetOwner(); + fBinnings.SetOwner(); + } + + void SetDataSetOn(const MDataSet &ds, UInt_t ntrain=(UInt_t)-1, UInt_t ntest=(UInt_t)-1) + { + ds.Copy(fDataSetOn); + + fDataSetOn.SetNumAnalysis(1); + + fNum[kTestOn] = ntrain; + fNum[kTrainOn] = ntest; + } + void SetDataSetOff(const MDataSet &ds, UInt_t ntrain=(UInt_t)-1, UInt_t ntest=(UInt_t)-1) + { + ds.Copy(fDataSetOff); + + fDataSetOff.SetNumAnalysis(2); + + fNum[kTestOff] = ntrain; + fNum[kTrainOff] = ntest; + } + + // Add Histogram + void AddHist(UInt_t nx); + void AddHist(UInt_t nx, UInt_t ny); + void AddHist(UInt_t nx, UInt_t ny, UInt_t nz); + + void AddBinning(UInt_t n, const MBinning &bins); + //void AddBinning(const MBinning &bins); + + // Standard user interface +/* + void SetWeights(Type_t typ, const char *rule) { if (fEnableWeights[typ]) return; fEnableWeights[typ]=kTRUE; AddPostTask(typ, rule); } + void SetWeights(Type_t typ, MTask *t) { if (fEnableWeights[typ]) return; fEnableWeights[typ]=kTRUE; AddPostTask(typ, t); } + + void AddPreTaskOn(MTask *t) { AddPreTask(kTrainOn, t); AddPreTask(kTestOn, t); } + void AddPreTaskOn(const char *rule, const char *name="MWeight") { AddPreTask(kTrainOn, rule, name); AddPreTask(kTestOn, rule, name); } + void AddPreTaskOff(MTask *t) { AddPreTask(kTrainOff, t); AddPreTask(kTestOff, t); } + void AddPreTaskOff(const char *rule, const char *name="MWeight") { AddPreTask(kTrainOff, rule, name); AddPreTask(kTestOff, rule, name); } + + void AddPostTaskOn(MTask *t) { AddPostTask(kTrainOn, t); AddPostTask(kTestOn, t); } + void AddPostTaskOn(const char *rule, const char *name="MWeight") { AddPostTask(kTrainOn, rule, name); AddPostTask(kTestOn, rule, name); } + void AddPostTaskOff(MTask *t) { AddPostTask(kTrainOff, t); AddPostTask(kTestOff, t); } + void AddPostTaskOff(const char *rule, const char *name="MWeight") { AddPostTask(kTrainOff, rule, name); AddPostTask(kTestOff, rule, name); } + + void AddPreTask(Type_t typ, MTask *t) { Add(fPreTasksSet[typ], t); } + void AddPreTask(Type_t typ, const char *rule, const char *name="MWeight") { AddPar(fPreTasksSet[typ], rule, name); } + void AddPostTask(Type_t typ, MTask *t) { Add(fPostTasksSet[typ], t); } + void AddPostTask(Type_t typ, const char *rule, const char *name="MWeight") { AddPar(fPostTasksSet[typ], rule, name); } + + void SetWeightsOn(const char *rule) { SetWeights(kTrainOn, rule); SetWeights(kTestOn, rule); } + void SetWeightsOn(MTask *t) { SetWeights(kTrainOn, t); SetWeights(kTestOn, t); } + void SetWeightsOff(const char *rule) { SetWeights(kTrainOff, rule); SetWeights(kTestOff, rule); } + void SetWeightsOff(MTask *t) { SetWeights(kTrainOff, t); SetWeights(kTestOff, t); } +*/ + void EnableRegression(Bool_t b=kTRUE) { fUseRegression = b; } + void EnableClassification(Bool_t b=kTRUE) { fUseRegression = !b; } + + // Main function to start processing + Bool_t Process(const char *out); + + ClassDef(MJTrainCuts, 0)//Class to help finding cuts using the random forest +}; + +#endif Index: trunk/Mars/mjtrain/Makefile =================================================================== --- trunk/Mars/mjtrain/Makefile (revision 9869) +++ trunk/Mars/mjtrain/Makefile (revision 9870) @@ -28,5 +28,6 @@ MJTrainDisp.cc \ MJTrainImpact.cc \ - MJTrainSeparation.cc + MJTrainSeparation.cc \ + MJTrainCuts.cc ############################################################ Index: trunk/Mars/mjtrain/TrainLinkDef.h =================================================================== --- trunk/Mars/mjtrain/TrainLinkDef.h (revision 9869) +++ trunk/Mars/mjtrain/TrainLinkDef.h (revision 9870) @@ -10,4 +10,5 @@ #pragma link C++ class MJTrainImpact+; #pragma link C++ class MJTrainSeparation+; +#pragma link C++ class MJTrainCuts+; #endif