/* ======================================================================== *\ ! ! * ! * 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 Hengstebeck 3/2003 ! ! ! Copyright: MAGIC Software Development, 2000-2003 ! ! \* ======================================================================== */ #include "MMcEvt.hxx" void RanForestPD(Int_t minsize=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", "star_gamma_train.root"); Float_t numgammas = read.GetEntries(); read.AddFile("star_proton_train.root"); read.AddFile("star_helium_train.root"); Float_t numhadrons = read.GetEntries() - numgammas; // Fraction of gammas to be used in training: // (less than total sample, to speed up execution) // Float_t gamma_frac = 2.*(numhadrons / numgammas); read.DisableAutoScheme(); tlist.AddToList(&read); MFParticleId fhadrons("MMcEvt", '!', MMcEvt::kGAMMA); tlist.AddToList(&fhadrons); MHMatrix matrixg("MatrixGammas"); // matrixg.AddColumn("floor(0.5+(cos(MMcEvt.fTelescopeTheta)/0.01))"); // matrixg.AddColumn("MMcEvt.fTelescopePhi"); // matrixg.AddColumn("MSigmabar.fSigmabar"); matrixg.AddColumn("MHillas.fSize"); matrixg.AddColumn("MHillasSrc.fDist"); // matrixg.AddColumn("MHillas.fWidth"); // matrixg.AddColumn("MHillas.fLength"); // Scaled Width and Length: // matrixg.AddColumn("MHillas.fWidth/(-13.1618+(10.0492*log10(MHillas.fSize)))"); matrixg.AddColumn("MHillas.fLength/(-57.3784+(32.6131*log10(MHillas.fSize)))"); // matrixg.AddColumn("(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))"); matrixg.AddColumn("MNewImagePar.fConc"); // matrixg.AddColumn("MNewImagePar.fConc1"); // matrixg.AddColumn("MNewImagePar.fAsym"); // matrixg.AddColumn("MHillasExt.fM3Trans"); // matrixg.AddColumn("abs(MHillasSrc.fHeadTail)"); // matrixg.AddColumn("MNewImagePar.fLeakage1"); // // Third moment, with sign referred to source position: // matrixg.AddColumn("MHillasExt.fM3Long*MHillasSrc.fCosDeltaAlpha/abs(MHillasSrc.fCosDeltaAlpha)"); //# matrixg.AddColumn("MHillasSrc.fAlpha"); // // SIZE cut: // TString dummy("(MHillas.fSize<"); dummy += minsize; dummy += ")"; // AM, useful FOR High Energy only: // dummy += " || (MImagePar.fNumSatPixelsLG>0)"; MContinue SizeCut(dummy); tlist.AddToList(&SizeCut); // AM TEST!!!, no muons in train // dummy = "(MMcEvt.fMuonCphFraction>0.1)"; // MContinue MuonsCut(dummy); // tlist.AddToList(&MuonsCut); plist.AddToList(&matrixg); MHMatrix matrixh("MatrixHadrons"); matrixh.AddColumns(matrixg.GetColumns()); plist.AddToList(&matrixh); MFillH fillmath("MatrixHadrons"); fillmath.SetFilter(&fhadrons); 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; 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); 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_gamma_test.root"); read3.AddFile("star_proton_test.root"); read3.AddFile("star_helium_test.root"); read3.DisableAutoScheme(); tlist3.AddToList(&read3); tlist3.AddToList(&SizeCut); MRanForestCalc calc; tlist3.AddToList(&calc); // parameter containers you want to save // TString outfname = "star_S"; outfname += minsize; outfname += "_all.root"; MWriteRootFile write(outfname, "recreate"); write.AddContainer("MMcEvt", "Events"); 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("MRawRunHeader", "RunHeaders"); write.AddContainer("MSrcPosCam", "RunHeaders"); write.AddContainer("MMcRunHeader", "RunHeaders"); /* write.AddContainer("MMcTrig", "Events"); write.AddContainer("MRawEvtData", "Events"); write.AddContainer("MRawEvtHeader", "Events"); */ MFillH fillh3a("MHHadronness"); MFillH fillh3b("MHRanForest"); tlist3.AddToList(&fillh3a); tlist3.AddToList(&fillh3b); tlist3.AddToList(&write); evtloop.SetProgressBar(&bar); // // Execute your analysis // if (!evtloop.Eventloop()) return; tlist3.PrintStatistics(); plist.FindObject("MHRanForest")->DrawClone(); TCanvas *had = new TCanvas; had->cd(); plist.FindObject("MHHadronness")->Draw(); TString gifname = "had_S"; gifname += minsize; gifname += ".gif"; had->SaveAs(gifname); had->Close(); plist.FindObject("MHHadronness")->DrawClone(); plist.FindObject("MHHadronness")->Print();//*/ return; }