/* ======================================================================== *\ ! ! * ! * 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 ! ! Copyright: MAGIC Software Development, 2006 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // MJTrainSeparation // //////////////////////////////////////////////////////////////////////////// #include "MJTrainSeparation.h" #include #include #include #include #include #include #include #include #include "MHMatrix.h" #include "MLog.h" #include "MLogManip.h" // tools #include "MMath.h" #include "MDataSet.h" #include "MTFillMatrix.h" #include "MStatusDisplay.h" // eventloop #include "MParList.h" #include "MTaskList.h" #include "MEvtLoop.h" // tasks #include "MReadMarsFile.h" #include "MReadReports.h" #include "MContinue.h" #include "MFillH.h" #include "MSrcPosRndm.h" #include "MHillasCalc.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" // wobble #include "MPointingPos.h" #include "MPointingDevCalc.h" #include "../mastro/MObservatory.h" #include "MSrcPosCalc.h" #include "MSrcPosCorrect.h" #include "../mimage/MHillasCalc.h" ClassImp(MJTrainSeparation); using namespace std; void MJTrainSeparation::DisplayResult(MH3 &h31, MH3 &h32, Float_t ontime) { TH2D &g = (TH2D&)h32.GetHist(); TH2D &h = (TH2D&)h31.GetHist(); h.SetMarkerColor(kRed); g.SetMarkerColor(kGreen); 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 sig1 = MMath::SignificanceLiMa(s+b, b); const Float_t sig2 = s<1 ? 0 : MMath::SignificanceLiMa(s+b, b)*TMath::Log10(s); gr3.SetPoint(y, h.GetYaxis()->GetBinLowEdge(y+2), sig1); gr4.SetPoint(y, h.GetYaxis()->GetBinLowEdge(y+2), sig2); } delete hx; delete gx; if (fDisplay) { TCanvas &c = fDisplay->AddTab("OptCut"); c.SetBorderMode(0); c.Divide(2,2); 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); 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 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"); if (!set.AddFilesOn(chain)) return kFALSE; 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, Double_t &num) const { TChain chain1("Events"); if (!set.AddFilesOff(chain1)) return kFALSE; 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"); if (!set.AddFilesOff(chain)) return kFALSE; chain.Draw("MEffectiveOnTime.fVal", "MEffectiveOnTime.fVal", "goff"); TH1 *h = dynamic_cast(gROOT->FindObject("htemp")); if (!h) { *fLog << err << "ERROR - Weird things are happening (htemp not found)!" << endl; return -1; } const Double_t ontime = h->Integral(); delete h; if (ontime<1) { *fLog << err << "ERROR - Less than 1s of effective observation time found in Train-Data." << endl; return -1; } return num/ontime; } Double_t MJTrainSeparation::GetNumMC(MDataSet &set) const { TChain chain1("Events"); if (!set.AddFilesOn(chain1)) return kFALSE; 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; } return num; } Float_t MJTrainSeparation::AutoTrain(MDataSet &set, Type_t typon, Type_t typof, Float_t flux) { Double_t num, min, max; if (!GetEventsProduced(set, num, min, max)) return -1; *fLog << inf << "Using build-in radius of 300m to calculate collection area!" << endl; // Target spectrum TF1 flx("Flux", "[0]/1000*(x/1000)^(-2.6)", min, max); flx.SetParameter(0, flux); // 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] *fLog << "Source Spectrum: " << flux << " * (E/TeV)^(-2.6) * TeV*m^2*s" << endl; *fLog << "Gamma rate from the source inside the MC production area: " << R << "Hz" << endl; // Number N of events produced (in trainings sample) const Double_t N = num; //[#] *fLog << "Events produced by MC inside the production area: " << TMath::Nint(num) << endl; // This correponds to an observation time T [s] const Double_t T = N/R; //[s] *fLog << "Total time produced by the Monte Carlo: " << T << "s" << endl; // With an average data rate after star of Double_t data=0; const Double_t r = GetDataRate(set, data); //[Hz] Double_t ontime = data/r; *fLog << "Events measured per second effective on time: " << r << "Hz" << endl; *fLog << "Total effective on time: " << ontime << "s" << endl; const Double_t ratio = T/ontime; *fLog << "Ratio of Monte Carlo to data observation time: " << ratio << endl; // 3570.5/43440.2 = 0.082 // this yields a number of n events to be read for training const Double_t n = r*T; //[#] *fLog << "Events to be read from the data sample: " << TMath::Nint(n) << endl; *fLog << "Events available in data sample: " << data << endl; if (r<0) return -1; Double_t nummc = GetNumMC(set); *fLog << "Events available in MC sample: " << nummc << endl; // *fLog << "MC read probability: " << data/n << endl; // more data requested than available => Scale down num MC events Double_t on, off; if (data0 && fNum[typon]Separator("Auto-Training -- Train-Data"); if (AutoTrain(fDataSetTrain, kTrainOn, kTrainOff, fFluxTrain)<0) return kFALSE; fLog->Separator("Auto-Training -- Test-Data"); ontime = AutoTrain(fDataSetTest, kTestOn, kTestOff, fFluxTest); if (ontime<0) return kFALSE; } // --------------------- Setup files -------------------- MReadReports read1;//("Events"); MReadMarsFile read2("Events"); MReadMarsFile read3("Events"); MReadReports read4;//("Events"); //read1.DisableAutoScheme(); read2.DisableAutoScheme(); read3.DisableAutoScheme(); //read4.DisableAutoScheme(); read1.AddTree("Events", "MTime.", MReadReports::kMaster); read4.AddTree("Events", "MTime.", MReadReports::kMaster); read1.AddTree("Drive", MReadReports::kRequired); read4.AddTree("Drive", MReadReports::kRequired); read1.AddTree("Starguider", MReadReports::kRequired); read4.AddTree("Starguider", MReadReports::kRequired); // Setup four reading tasks with the on- and off-data of the two datasets if (!fDataSetTrain.AddFilesOn(read1)) return kFALSE; const Bool_t setrc1 = fDataSetTrain.IsWobbleMode() ? fDataSetTrain.AddFilesOn(read3) : fDataSetTrain.AddFilesOff(read3); const Bool_t setrc2 = fDataSetTest.IsWobbleMode() ? fDataSetTest.AddFilesOn(read2) : fDataSetTest.AddFilesOff(read2); if (!setrc1 || !setrc2) return kFALSE; if (!fDataSetTest.AddFilesOn(read4)) return kFALSE; // ----------------------- Setup RF Matrix ---------------------- MHMatrix train("Train"); train.AddColumns(fRules); if (fEnableWeights[kTrainOn] || fEnableWeights[kTrainOff]) train.AddColumn("MWeight.fVal"); train.AddColumn("MHadronness.fVal"); // ----------------------- Fill Matrix RF ---------------------- // Setup the hadronness container identifying gammas and off-data // and setup a container for the weights MParameterD had("MHadronness"); MParameterD wgt("MWeight"); // Add them to the parameter list MParList plistx; plistx.AddToList(this); // take care of fDisplay! plistx.AddToList(&had); plistx.AddToList(&wgt); // Setup the tool class to fill the matrix MTFillMatrix fill; fill.SetLogStream(fLog); fill.SetDisplay(fDisplay); fill.AddPreCuts(fPreCuts); fill.AddPreCuts(fTrainCuts); // Set classifier for gammas had.SetVal(0); wgt.SetVal(1); // How to get source position from off- and on-data? MPointingPos source("MSourcePos"); MObservatory obs; MSrcPosCalc scalc; MSrcPosCorrect scor; MHillasCalc hcalcw; MHillasCalc hcalcw2; MPointingDevCalc devcalc; scalc.SetMode(MSrcPosCalc::kDefault); // kWobble for off-source hcalcw.SetFlags(MHillasCalc::kCalcHillasSrc); hcalcw2.SetFlags(MHillasCalc::kCalcHillasSrc); hcalcw2.SetNameHillasSrc("MHillasSrcAnti"); hcalcw2.SetNameSrcPosCam("MSrcPosAnti"); if (fDataSetTrain.IsWobbleMode()) { // ******************************************************************* // Possible source position (eg. Wobble Mode) if (fDataSetTrain.HasSource()) { if (!fDataSetTrain.GetSourcePos(source)) return -1; *fLog << all; source.Print("RaDec"); } else *fLog << all << "No source position applied..." << endl; // La Palma Magic1 plistx.AddToList(&obs); plistx.AddToList(&source); TList tlist2; tlist2.Add(&scalc); tlist2.Add(&scor); tlist2.Add(&hcalcw); tlist2.Add(&hcalcw2); devcalc.SetStreamId("Starguider"); tlist2.Add(&devcalc); fill.AddPreTasks(tlist2); // ******************************************************************* } // Setup the tool class to read the gammas and read them fill.SetName("FillGammas"); fill.SetDestMatrix1(&train, fNum[kTrainOn]); fill.SetReader(&read1); fill.AddPreTasks(fPreTasksSet[kTrainOn]); fill.AddPreTasks(fPreTasks); fill.AddPostTasks(fPostTasksSet[kTrainOn]); fill.AddPostTasks(fPostTasks); if (!fill.Process(plistx)) return kFALSE; // Check the number or read events const Int_t numgammastrn = train.GetNumRows(); if (numgammastrn==0) { *fLog << err << "ERROR - No gammas available for training... aborting." << endl; return kFALSE; } // Remove possible post tasks fill.ClearPreTasks(); fill.ClearPostTasks(); // Set classifier for background had.SetVal(1); wgt.SetVal(1); // In case of wobble mode we have to do something special MSrcPosRndm srcrndm; srcrndm.SetDistOfSource(0.4); MHillasCalc hcalc; MHillasCalc hcalc2; hcalc.SetFlags(MHillasCalc::kCalcHillasSrc); hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc); hcalc2.SetNameHillasSrc("MHillasSrcAnti"); hcalc2.SetNameSrcPosCam("MSrcPosAnti"); if (fDataSetTrain.IsWobbleMode()) { scalc.SetMode(MSrcPosCalc::kWobble); // kWobble for off-source fPreTasksSet[kTrainOff].AddFirst(&hcalc2); fPreTasksSet[kTrainOff].AddFirst(&hcalc); //fPreTasksSet[kTrainOff].AddFirst(&srcrndm); fPreTasksSet[kTrainOff].AddFirst(&scor); fPreTasksSet[kTrainOff].AddFirst(&scalc); } // Setup the tool class to read the background and read them fill.SetName("FillBackground"); fill.SetDestMatrix1(&train, fNum[kTrainOff]); fill.SetReader(&read3); fill.AddPreTasks(fPreTasksSet[kTrainOff]); fill.AddPreTasks(fPreTasks); fill.AddPostTasks(fPostTasksSet[kTrainOff]); fill.AddPostTasks(fPostTasks); if (!fill.Process(plistx)) return kFALSE; // Check the number or read events const Int_t numbackgrndtrn = train.GetNumRows()-numgammastrn; if (numbackgrndtrn==0) { *fLog << err << "ERROR - No background available for training... aborting." << endl; return kFALSE; } // ------------------------ Train RF -------------------------- MRanForestCalc rf; rf.SetNumTrees(fNumTrees); rf.SetNdSize(fNdSize); rf.SetNumTry(fNumTry); rf.SetNumObsoleteVariables(1); rf.SetLastDataColumnHasWeights(fEnableWeights[kTrainOn] || fEnableWeights[kTrainOff]); rf.SetDebug(fDebug); rf.SetDisplay(fDisplay); rf.SetLogStream(fLog); rf.SetFileName(out); rf.SetNameOutput("MHadronness"); // Train the random forest either by classification or regression if (fUseRegression) { if (!rf.TrainRegression(train)) // regression return kFALSE; } else { if (!rf.TrainSingleRF(train)) // classification return kFALSE; } // 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; if (fDataSetTrain.IsWobbleMode()) *fLog << " * random source position in a distance of 0.4°" << endl; *fLog << endl; *fLog << "Events used for training:" << endl; *fLog << " * Gammas: " << numgammastrn << endl; *fLog << " * Background: " << numbackgrndtrn << endl; *fLog << endl; *fLog << "Gamma/Background ratio:" << endl; *fLog << " * Requested: " << (float)fNum[kTrainOn]/fNum[kTrainOff] << endl; *fLog << " * Result: " << (float)numgammastrn/numbackgrndtrn << 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; // Chekc if testing is requested if (!fDataSetTest.IsValid()) return kTRUE; // --------------------- Display result ---------------------- fLog->Separator("Test"); clock.Continue(); // 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); // ----- 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"); MHHadronness hist; // ----- Setup tasks ----- MFillH fillh0(&hist, "", "FillHadronness"); MFillH fillh1(&h31); MFillH fillh2(&h32); MFillH fillh4(&h40); 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("Background"); fillh2.SetNameTab("GammasH"); fillh4.SetNameTab("GammasE"); fillh0.SetBit(MFillH::kDoNotDisplay); // ----- Setup filter ----- MFilterList precuts; precuts.AddToList(fPreCuts); precuts.AddToList(fTestCuts); MContinue c0(&precuts); c0.SetName("PreCuts"); c0.SetInverted(); MFEventSelector sel; // FIXME: USING IT (WITH PROB?) in READ will by much faster!!! sel.SetNumSelectEvts(fNum[kTestOff]); MContinue c1(&sel); c1.SetInverted(); // ----- Setup tasklist ----- tlist.AddToList(&read2); if (fDataSetTest.IsWobbleMode()) { tlist.AddToList(&srcrndm); tlist.AddToList(&hcalc); } tlist.AddToList(&c1); tlist.AddToList(fPreTasksSet[kTestOff]); tlist.AddToList(fPreTasks); tlist.AddToList(&c0); tlist.AddToList(&rf); tlist.AddToList(fPostTasksSet[kTestOff]); tlist.AddToList(fPostTasks); tlist.AddToList(&fillh0); tlist.AddToList(&fillh1); // Enable Acceleration tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime); // ----- Run eventloop on background ----- MEvtLoop loop; loop.SetDisplay(fDisplay); loop.SetLogStream(fLog); loop.SetParList(&plist); wgt.SetVal(1); if (!loop.Eventloop()) return kFALSE; /* if (!loop.GetDisplay()) { gLog << warn << "Display closed by user... execution aborted." << endl << endl; return kFALSE; } */ // ----- Setup and run eventloop on gammas ----- sel.SetNumSelectEvts(fNum[kTestOn]); fillh0.ResetBit(MFillH::kDoNotDisplay); // Remove PreTasksOff and PostTasksOff from the list tlist.RemoveFromList(fPreTasksSet[kTestOff]); tlist.RemoveFromList(fPostTasksSet[kTestOff]); // replace the reading task by a new one tlist.Replace(&read4); if (fDataSetTest.IsWobbleMode()) { // ******************************************************************* // Possible source position (eg. Wobble Mode) if (fDataSetTest.HasSource()) { if (!fDataSetTest.GetSourcePos(source)) return -1; *fLog << all; source.Print("RaDec"); } else *fLog << all << "No source position applied..." << endl; // La Palma Magic1 plist.AddToList(&obs); plist.AddToList(&source); // How to get source position from off- and on-data? tlist.AddToListAfter(&scalc, &read4); tlist.AddToListAfter(&scor, &scalc); tlist.AddToListAfter(&hcalcw, &scor); tlist.AddToList(&devcalc, "Starguider"); // ******************************************************************* } // Add the PreTasksOn directly after the reading task tlist.AddToListAfter(fPreTasksSet[kTestOn], &c1); // Add the PostTasksOn after rf tlist.AddToListAfter(fPostTasksSet[kTestOn], &rf); // Replace fillh1 by fillh2 tlist.Replace(&fillh2); // Add fillh4 after the new fillh2 tlist.AddToListAfter(&fillh4, &fillh2); // Enable Acceleration tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime); wgt.SetVal(1); if (!loop.Eventloop()) return kFALSE; // Show what was going on in the testing const Double_t numgammastst = h32.GetHist().GetEntries(); const Double_t numbackgrndtst = 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; if (fDataSetTrain.IsWobbleMode()) *fLog << " * random source position in a distance of 0.4°" << endl; *fLog << endl; *fLog << "Events used for test:" << endl; *fLog << " * Gammas: " << numgammastst << endl; *fLog << " * Background: " << numbackgrndtst << endl; *fLog << endl; *fLog << "Gamma/Background ratio:" << endl; *fLog << " * Requested: " << (float)fNum[kTestOn]/fNum[kTestOff] << endl; *fLog << " * Result: " << (float)numgammastst/numbackgrndtst << endl; *fLog << endl; // Display the result plots 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 the display if (!WriteDisplay(out)) return kFALSE; return kTRUE; }