Index: trunk/MagicSoft/Mars/macros/RanForest.C =================================================================== --- trunk/MagicSoft/Mars/macros/RanForest.C (revision 1859) +++ trunk/MagicSoft/Mars/macros/RanForest.C (revision 1859) @@ -0,0 +1,172 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +void RanForest() +{ + // + // 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", "~/MagicData/data1/CT1Data/ONTrain/*.root"); + read.DisableAutoScheme(); + tlist.AddToList(&read); + + MFParticleId fgamma("MMcEvt", '=', kGAMMA); + tlist.AddToList(&fgamma); + + MFParticleId fhadrons("MMcEvt", '!', kGAMMA); + tlist.AddToList(&fhadrons); + + MHMatrix matrixg("MatrixGammas"); + matrixg.AddColumn("cos(MMcEvt.fTelescopeTheta)"); + matrixg.AddColumn("MSigmabar.fSigmabar"); + matrixg.AddColumn("log10(Hillas.fSize)"); + matrixg.AddColumn("HillasSrc.fDist"); + matrixg.AddColumn("Hillas.fWidth"); + matrixg.AddColumn("Hillas.fLength"); + matrixg.AddColumn("log10(Hillas.fSize)/Hillas.fWidth*Hillas.fLength"); + matrixg.AddColumn("abs(Hillas.fM3Long)"); + matrixg.AddColumn("Hillas.fConc"); + matrixg.AddColumn("Hillas.fConc1"); + + //matrixg.AddColumn("abs(Hillas.fAsym)"); + //matrixg.AddColumn("abs(Hillas.fM3Trans)"); + //matrixg.AddColumn("abs(HillasSrc.fHeadTail)"); + //matrixg.AddColumn("abs(HillasSrc.fAlpha)"); + + plist.AddToList(&read); + + plist.AddToList(&matrixg); + + MHMatrix matrixh("MatrixHadrons"); + matrixh.AddColumns(matrixg.GetColumns()); + plist.AddToList(&matrixh); + + MFillH fillmatg("MatrixGammas"); + fillmatg.SetFilter(&fgamma); + tlist.AddToList(&fillmatg); + + MFillH fillmath("MatrixHadrons"); + fillmath.SetFilter(&fhadrons); + tlist.AddToList(&fillmath); + + // + // Create and setup the eventloop + // + MEvtLoop evtloop; + evtloop.SetParList(&plist); + + // + // 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(50); + rfgrow2.SetNumTry(3); + rfgrow2.SetNdSize(1); + + MWriteRootFile rfwrite2("RF.root"); + rfwrite2.AddContainer("MRanTree","Tree"); //save all trees + MFillH fillh2("MHRanForestGini"); + + tlist2.AddToList(&rfgrow2); + tlist2.AddToList(&rfwrite2); + tlist2.AddToList(&fillh2); + + // + // Execute tree growing (now the eventloop is actually a treeloop) + // + 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", "~/MagicData/data1/CT1Data/ONTest/*.root"); + read3.DisableAutoScheme(); + tlist3.AddToList(&read3); + + MRanForestCalc calc; + tlist3.AddToList(&calc); + + MFillH fillh3a("MHHadronness"); + MFillH fillh3b("MHRanForest"); + + tlist3.AddToList(&fillh3a); + tlist3.AddToList(&fillh3b); + + MProgressBar bar; + evtloop.SetProgressBar(&bar); + + // + // Execute your analysis + // + if (!evtloop.Eventloop()) + return; + + tlist3.PrintStatistics(); + + plist.FindObject("MHRanForest")->DrawClone(); + plist.FindObject("MHHadronness")->DrawClone(); + plist.FindObject("MHHadronness")->Print();//*/ +} Index: trunk/MagicSoft/Mars/macros/RanForest2.C =================================================================== --- trunk/MagicSoft/Mars/macros/RanForest2.C (revision 1859) +++ trunk/MagicSoft/Mars/macros/RanForest2.C (revision 1859) @@ -0,0 +1,110 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +void RanForest2() +{ + // + // 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); + // + // --------------------------------------------------------------- + // --------------------------------------------------------------- + // first event loop: the trees of the random forest are read in + // --------------------------------------------------------------- + // --------------------------------------------------------------- + // + MReadTree read("Tree","RF.root"); + read.DisableAutoScheme(); + + MRanForestFill rffill; + rffill.SetNumTrees(100); + + tlist.AddToList(&read); + tlist.AddToList(&rffill); + + // + // Create and setup the eventloop + // + MEvtLoop evtloop; + evtloop.SetParList(&plist); + + // + // Execute tree reading (now the eventloop is actually a treeloop) + // + if (!evtloop.Eventloop()) + return; + + tlist.PrintStatistics(); + + // --------------------------------------------------------------- + // --------------------------------------------------------------- + // second event loop: the control sample 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 tlist2; + + plist.Replace(&tlist2); + + MReadMarsFile read2("Events", "~/MagicData/data1/CT1Data/ONTest/*.root"); + read2.DisableAutoScheme(); + tlist2.AddToList(&read2); + + MRanForestCalc calc; + tlist2.AddToList(&calc); + + MFillH fillh2a("MHHadronness"); + MFillH fillh2b("MHRanForest"); + + tlist2.AddToList(&fillh2a); + tlist2.AddToList(&fillh2b); + + // + // Execute your analysis + // + MProgressBar bar; + evtloop.SetProgressBar(&bar); + + if (!evtloop.Eventloop()) + return; + + tlist2.PrintStatistics(); + + plist.FindObject("MHRanForest")->DrawClone(); + plist.FindObject("MHHadronness")->DrawClone(); + plist.FindObject("MHHadronness")->Print(); +} Index: trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 1858) +++ trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 1859) @@ -25,4 +25,10 @@ #pragma link C++ class MCompProbCalc+; #pragma link C++ class MMultiDimDistCalc+; + +#pragma link C++ class MRanTree+; +#pragma link C++ class MRanForest+; +#pragma link C++ class MRanForestGrow+; +#pragma link C++ class MRanForestCalc+; +#pragma link C++ class MRanForestFill+; #pragma link C++ class MSrcPosCam+; Index: trunk/MagicSoft/Mars/manalysis/MRanForest.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MRanForest.cc (revision 1859) +++ trunk/MagicSoft/Mars/manalysis/MRanForest.cc (revision 1859) @@ -0,0 +1,417 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MRanForest // +// // +// ParameterContainer for Forest structure // +// // +// A random forest can be grown by calling GrowForest. // +// In advance SetupGrow must be called in order to initialize arrays and // +// do some preprocessing. // +// GrowForest() provides the training data for a single tree (bootstrap // +// aggregate procedure) // +// // +// Essentially two random elements serve to provide a "random" forest, // +// namely bootstrap aggregating (which is done in GrowForest()) and random // +// split selection (which is subject to MRanTree::GrowTree()) // +///////////////////////////////////////////////////////////////////////////// +#include "MRanForest.h" +#include +#include + +#include // gFile + +#include "MLog.h" +#include "MLogManip.h" + + +ClassImp(MRanForest); +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MRanForest::MRanForest(const char *name, const char *title) : fNumTrees(100), fRanTree(NULL),fUsePriors(kFALSE) +{ + fName = name ? name : "MRanForest"; + fTitle = title ? title : "Storage container for Random Forest"; +} + +void MRanForest::SetNumTrees(Int_t n) +{ + //at least 1 tree + fNumTrees=TMath::Max(n,1); + fTreeHad.Set(fNumTrees); + fTreeHad.Reset(); + + return; +} + +void MRanForest::SetPriors(Float_t prior_had, Float_t prior_gam) +{ + Float_t sum; + + sum=prior_gam+prior_had; + + prior_gam/=sum; + prior_had/=sum; + + fPrior[0]=prior_had; + fPrior[1]=prior_gam; + + fUsePriors=kTRUE; + + return; +} + +Double_t MRanForest::CalcHadroness(TVector &event) +{ + Double_t hadroness=0; + Double_t ntree=0; + MRanTree *tree; + TIter forest(fForest); + while ((tree=(MRanTree*)forest.Next())) + { + fTreeHad[ntree]=tree->TreeHad(event); + hadroness+=fTreeHad[ntree]; + ntree++; + } + return hadroness/ntree; +} + +void MRanForest::SetupForest() +{ + fForest=new TObjArray(); + fForest->SetOwner(kTRUE); + + return; +} + +void MRanForest::SetTree(MRanTree *rantree) +{ + // initialize current tree for tree-growing loop in GrowForest + // ndsize + numtry are set in MRanForestGrow!! + fRanTree=rantree; + + return; +} + +Bool_t MRanForest::AddTree() +{ + if (!fRanTree) + return kFALSE; + MRanTree *newtree=(MRanTree*)fRanTree->Clone(); + fForest->Add(newtree); + + return kTRUE; +} + +Bool_t MRanForest::SetupGrow(MHMatrix *mhad,MHMatrix *mgam) +{ + // pointer to training data + fHadrons=mhad; + fGammas=mgam; + + // determine data entries and dimension of Hillas-parameter space + fNumHad=fHadrons->GetM().GetNrows(); + fNumGam=fGammas->GetM().GetNrows(); + fNumDim=fHadrons->GetM().GetNcols(); + + if (fNumDim!=fHadrons->GetM().GetNcols()) return kFALSE; + + fNumData=fNumHad+fNumGam; + + // allocating and initializing arrays + fHadTrue.Set(fNumData); + fHadTrue.Reset(); + fHadEst.Set(fNumData); + + fPrior.Set(2); + fClassPop.Set(2); + fWeight.Set(fNumData); + fNTimesOutBag.Set(fNumData); + fNTimesOutBag.Reset(); + + fGiniDec.Set(fNumDim); + + fDataSort.Set(fNumDim*fNumData); + fDataRang.Set(fNumDim*fNumData); + + // calculating class populations (= no. of gammas and hadrons) + fClassPop.Reset(); + for(Int_t n=0;n=1) ? + Float_t(fNumData)/Float_t(fClassPop[j]):0; + + for(Int_t n=0;nSetRules(fGammas->GetColumns()); + fTreeNo=0; + + return kTRUE; +} + +Bool_t MRanForest::GrowForest() +{ + Int_t ninbag=0; + TArrayI datsortinbag; + TArrayF classpopw; + TArrayI jinbag; + TArrayF winbag; + + jinbag.Set(fNumData); + winbag.Set(fNumData); + classpopw.Set(2); + + TMatrix hadrons=fHadrons->GetM(); + TMatrix gammas=fGammas->GetM(); + + fTreeNo++; + + // initialize running output + if(fTreeNo==1) + { + cout< overestim. of error)"< sampling with replacement: + // + // The integer k is randomly (uniformly) chosen from the set + // {0,1,...,fNumData-1}, which is the set of the index numbers of + // all events in the training sample + + for (Int_t n=0;nGrowTree(hadrons,gammas,fNumData,fNumDim,fHadTrue,datsortinbag, + fDataRang,fGiniDec,classpopw,jinbag,winbag,fWeight,fRand); + + // error-estimates from out-of-bag data (oob data): + // + // For a single tree the events not(!) contained in the bootstrap sample of + // this tree can be used to obtain estimates for the classification error of + // this tree. + // If you take a certain event, it is contained in the oob-data of 1/3 of + // the trees (see comment to ModifyData). This means that the classification error + // determined from oob-data is underestimated, but can still be taken as upper limit. + + TVector event(fNumDim); + for(Int_t ievt=0;ievt0)continue; + event=TMatrixRow(hadrons,ievt); + fHadEst[ievt]+=fRanTree->TreeHad(event); + fNTimesOutBag[ievt]++; + } + for(Int_t ievt=0;ievt0)continue; + event=TMatrixRow(gammas,ievt); + fHadEst[fNumHad+ievt]+=fRanTree->TreeHad(event); + fNTimesOutBag[fNumHad+ievt]++; + } + + Int_t n=0; + fErr=0; + for(Int_t ievt=0;ievtGetM(); + TMatrix gammas=fGammas->GetM(); + + for (Int_t j=0;j=ninbag) break; + } + } + return; +} + +Bool_t MRanForest::AsciiWrite(ostream &out) const +{ + Int_t n=0; + MRanTree *tree; + TIter forest(fForest); + + while ((tree=(MRanTree*)forest.Next())) + { + tree->AsciiWrite(out); + n++; + } + + return n==fNumTrees; +} Index: trunk/MagicSoft/Mars/manalysis/MRanForest.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MRanForest.h (revision 1859) +++ trunk/MagicSoft/Mars/manalysis/MRanForest.h (revision 1859) @@ -0,0 +1,116 @@ +#ifndef MARS_MRanForest +#define MARS_MRanForest + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +#ifndef ROOT_TArrayI +#include +#endif + +#ifndef ROOT_TArrayF +#include +#endif + +#ifndef ROOT_TArrayD +#include +#endif + +#ifndef ROOT_TObjArray +#include +#endif + +#ifndef ROOT_TRandom +#include +#endif + +#ifndef ROOT_TMatrix +#include +#endif + +#ifndef MARS_MHMatrix +#include "MHMatrix.h" +#endif + +#ifndef MARS_MRanTree +#include "MRanTree.h" +#endif + +class MRanForest : public MParContainer +{ +private: + Int_t fNumTrees; + Int_t fTreeNo; + + TRandom fRand; + MRanTree *fRanTree; + TObjArray *fForest; + + // training data + MHMatrix *fGammas; + MHMatrix *fHadrons; + + Int_t fNumGam; + Int_t fNumHad; + Int_t fNumData; + Int_t fNumDim; + + // true and estimated hadronness + TArrayI fHadTrue; + TArrayF fHadEst; + + // data sorted according to parameters + TArrayI fDataSort; + TArrayI fDataRang; + TArrayI fClassPop; + + // weights + Bool_t fUsePriors; + TArrayF fPrior; + TArrayF fWeight; + TArrayI fNTimesOutBag; + + // estimates for classification error of growing forest + TArrayD fTreeHad; + Float_t fErr; + + // decrease in Gini-index + TArrayF fGiniDec; + +protected: + // create and modify (->due to bagging) fDataSort + void CreateDataSort(); + void ModifyDataSort(TArrayI &datsortinbag,Int_t ninbag,TArrayI &jinbag); + +public: + MRanForest(const char *name=NULL, const char *title=NULL); + + // initialize forest + void SetPriors(Float_t prior_had, Float_t prior_gam); + void SetNumTrees(Int_t n); + + // tree growing + void SetupForest(); + void SetTree(MRanTree *rantree); + Bool_t AddTree(); + Bool_t SetupGrow(MHMatrix *mhad,MHMatrix *mgam); + Bool_t GrowForest(); + + // getter methods + TObjArray *GetForest() { return fForest; } + Int_t GetNumTrees() const { return fNumTrees; } + Int_t GetNumData() const { return fNumData; } + Int_t GetNumDim() const { return fNumDim; } + Double_t GetTreeHad(Int_t i) const { return fTreeHad.At(i); } + Float_t GetGiniDec(Int_t i) const { return fGiniDec.At(i); } + + // use forest to calculate hadronness of event + Double_t CalcHadroness(TVector &event); + + Bool_t AsciiWrite(ostream &out) const; + + ClassDef(MRanForest, 1) // Storage container for tree structure +}; + +#endif Index: trunk/MagicSoft/Mars/manalysis/MRanForestCalc.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MRanForestCalc.cc (revision 1859) +++ trunk/MagicSoft/Mars/manalysis/MRanForestCalc.cc (revision 1859) @@ -0,0 +1,143 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MRanForestCalc +// +// Calculates the hadroness of an event. It calculates a mean value of all +// classifications by the trees in a previously grown random forest. +// +// To use only n trees for your calculation use: +// MRanForestCalc::SetUseNumTrees(n); +// +//////////////////////////////////////////////////////////////////////////// +#include "MRanForestCalc.h" + +#include + +#include "MHMatrix.h" // must be before MLogManip.h + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" +#include "MDataChain.h" +#include "MDataArray.h" + +#include "MHadronness.h" + +ClassImp(MRanForestCalc); + +static const TString gsDefName = "MRanForestCalc"; +static const TString gsDefTitle = "Tree Classification Loop 1/2"; +// -------------------------------------------------------------------------- +// +// Setup histograms and the number of distances which are used for +// avaraging in CalcDist +// +MRanForestCalc::MRanForestCalc(const char *name, const char *title) + : fNum(100), fData(NULL) +{ + // + // set the name and title of this object + // + fName = name ? name : gsDefName.Data(); + fTitle = title ? title : gsDefTitle.Data(); +} + +// -------------------------------------------------------------------------- +// +// Delete the data chains +// +MRanForestCalc::~MRanForestCalc() +{ + // delete fData; +} + +// -------------------------------------------------------------------------- +// +// Needs: +// - MatrixGammas [MHMatrix] +// - MatrixHadrons {MHMatrix] +// - MHadronness +// - all data containers used to build the matrixes +// +// The matrix object can be filles using MFillH. And must be of the same +// number of columns (with the same meaning). +// +Bool_t MRanForestCalc::PreProcess(MParList *plist) +{ + fRanForest = (MRanForest*)plist->FindObject("MRanForest"); + if (!fRanForest) + { + *fLog << err << dbginf << "MRanForest not found... aborting." << endl; + return kFALSE; + } + + fRanTree = (MRanTree*)plist->FindObject("MRanTree"); + if (!fRanTree) + { + *fLog << err << dbginf << "MRanTree not found... aborting." << endl; + return kFALSE; + } + + fData = fRanTree->GetData(); + + if (!fData) + { + *fLog << err << dbginf << "Error matrix doesn't contain columns... aborting." << endl; + return kFALSE; + } + + if (!fData->PreProcess(plist)) + { + *fLog << err << dbginf << "PreProcessing of the MDataArray failed for the columns failed... aborting." << endl; + return kFALSE; + } + + fHadroness = (MHadronness*)plist->FindCreateObj("MHadronness"); + if (!fHadroness) + return kFALSE; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// +Bool_t MRanForestCalc::Process() +{ + const Double_t ncols = fData->GetNumEntries(); + TVector event(ncols); + + for (int i=0; iGetNumEntries(); i++) + event(i) = (*fData)(i); + + Double_t hadroness=fRanForest->CalcHadroness(event); + fHadroness->SetHadronness(hadroness); + + return kTRUE; +} + Index: trunk/MagicSoft/Mars/manalysis/MRanForestCalc.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MRanForestCalc.h (revision 1859) +++ trunk/MagicSoft/Mars/manalysis/MRanForestCalc.h (revision 1859) @@ -0,0 +1,38 @@ +#ifndef MARS_MRanForestCalc +#define MARS_MRanForestCalc + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef MARS_MRanForest +#include "MRanForest.h" +#endif + +class MParList; +class MHadronness; +class MDataArray; + +class MRanForestCalc : public MTask +{ +private: + Int_t fNum; // number of trees used to compute hadronness + + MHadronness *fHadroness; //! Output container for calculated hadroness + MDataArray *fData; //! Used to store the MDataChains to get the event values + MRanForest *fRanForest; + MRanTree *fRanTree; + +public: + MRanForestCalc(const char *name=NULL, const char *title=NULL); + ~MRanForestCalc(); + + void SetUseNumTrees(UShort_t n=100) { fNum = n; } + + Bool_t PreProcess(MParList *plist); + Bool_t Process(); + + ClassDef(MRanForestCalc, 0) // Task +}; + +#endif Index: trunk/MagicSoft/Mars/manalysis/MRanForestFill.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MRanForestFill.cc (revision 1859) +++ trunk/MagicSoft/Mars/manalysis/MRanForestFill.cc (revision 1859) @@ -0,0 +1,115 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MRanForestFill +// +// Calculates the hadroness of an event. It calculates a mean value of all +// classifications by the trees in a previously grown random forest. +// +// To use only n trees for your calculation use: +// MRanForestFill::SetUseNumTrees(n); +// +//////////////////////////////////////////////////////////////////////////// +#include "MRanForestFill.h" + +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" +#include "MDataChain.h" +#include "MDataArray.h" + +ClassImp(MRanForestFill); + +static const TString gsDefName = "MRanForestFill"; +static const TString gsDefTitle = "Tree Classification Loop 1/2"; +// -------------------------------------------------------------------------- +// +// +MRanForestFill::MRanForestFill(const char *name, const char *title):fNumTrees(100) +{ + // + // set the name and title of this object + // + fName = name ? name : gsDefName.Data(); + fTitle = title ? title : gsDefTitle.Data(); +} + +// -------------------------------------------------------------------------- +// +// Delete the data chains +// +MRanForestFill::~MRanForestFill() +{ + // delete fData; +} + +// -------------------------------------------------------------------------- +Bool_t MRanForestFill::PreProcess(MParList *plist) +{ + fRanTree = (MRanTree*)plist->FindObject("MRanTree"); + if (!fRanTree) + { + *fLog << err << dbginf << "MRanTree not found... aborting." << endl; + return kFALSE; + } + + fRanForest = (MRanForest*)plist->FindCreateObj("MRanForest"); + if (!fRanForest) + { + *fLog << err << dbginf << "MRanForest not found... aborting." << endl; + return kFALSE; + } + + fRanForest->SetupForest(); + fNum=0; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// +Bool_t MRanForestFill::Process() +{ + + fNum++; + fRanForest->SetTree(fRanTree); + if(!(fRanForest->AddTree())) + return kFALSE; + + return fNumSetNumTrees(fNum); + + return kTRUE; +} + Index: trunk/MagicSoft/Mars/manalysis/MRanForestFill.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MRanForestFill.h (revision 1859) +++ trunk/MagicSoft/Mars/manalysis/MRanForestFill.h (revision 1859) @@ -0,0 +1,41 @@ +#ifndef MARS_MRanForestFill +#define MARS_MRanForestFill + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef MARS_MRanForest +#include "MRanForest.h" +#endif + +#ifndef MARS_MRanTree +#include "MRanTree.h" +#endif + +class MParList; +class MDataArray; + +class MRanForestFill : public MTask +{ +private: + MRanTree *fRanTree; + MRanForest *fRanForest; + MDataArray *fData; + Int_t fNumTrees; + Int_t fNum; + +public: + MRanForestFill(const char *name=NULL, const char *title=NULL); + ~MRanForestFill(); + + void SetNumTrees(UShort_t n=100) { fNumTrees = n; } + + Bool_t PreProcess(MParList *plist); + Bool_t Process(); + Bool_t PostProcess(); + + ClassDef(MRanForestFill, 0) // Task +}; + +#endif Index: trunk/MagicSoft/Mars/manalysis/MRanForestGrow.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MRanForestGrow.cc (revision 1859) +++ trunk/MagicSoft/Mars/manalysis/MRanForestGrow.cc (revision 1859) @@ -0,0 +1,153 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MRanForestGrow // +// // +// Grows a random forest. // +// // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MRanForestGrow.h" + +#include + +#include "MHMatrix.h" // must be before MLogManip.h + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +ClassImp(MRanForestGrow); + +static const TString gsDefName = "MRanForestGrow"; +static const TString gsDefTitle = "Tree Classification Loop 1/2"; +// -------------------------------------------------------------------------- +// +// Setup histograms and the number of distances which are used for +// avaraging in CalcDist +// +MRanForestGrow::MRanForestGrow(const char *name, const char *title) + : fNumTrees(100),fNumTry(3),fNdSize(1) +{ + // + // set the name and title of this object + // + fName = name ? name : gsDefName.Data(); + fTitle = title ? title : gsDefTitle.Data(); + +} + +// -------------------------------------------------------------------------- +// +// Needs: +// - MatrixGammas [MHMatrix] +// - MatrixHadrons {MHMatrix] +// - MHadroness +// - all data containers used to build the matrixes +// +// The matrix object can be filles using MFillH. And must be of the same +// number of columns (with the same meaning). +// +Bool_t MRanForestGrow::PreProcess(MParList *plist) +{ + fMGammas = (MHMatrix*)plist->FindObject("MatrixGammas", "MHMatrix"); + if (!fMGammas) + { + *fLog << err << dbginf << "MatrixGammas [MHMatrix] not found... aborting." << endl; + return kFALSE; + } + + fMHadrons = (MHMatrix*)plist->FindObject("MatrixHadrons", "MHMatrix"); + if (!fMHadrons) + { + *fLog << err << dbginf << "MatrixHadrons [MHMatrix] not found... aborting." << endl; + return kFALSE; + } + + if (fMGammas->GetM().GetNcols() != fMHadrons->GetM().GetNcols()) + { + *fLog << err << dbginf << "Error matrices have different numbers of columns... aborting." << endl; + return kFALSE; + } + + fRanTree = (MRanTree*)plist->FindCreateObj("MRanTree"); + if (!fRanTree) + { + *fLog << err << dbginf << "MRanTree not found... aborting." << endl; + return kFALSE; + } + + fRanForest = (MRanForest*)plist->FindCreateObj("MRanForest"); + if (!fRanForest) + { + *fLog << err << dbginf << "MRanForest not found... aborting." << endl; + return kFALSE; + } + + fRanTree->SetNumTry(fNumTry); + fRanTree->SetNdSize(fNdSize); + + fRanForest->SetNumTrees(fNumTrees); + fRanForest->SetTree(fRanTree); + + return fRanForest->SetupGrow(fMHadrons,fMGammas); +} + +// -------------------------------------------------------------------------- +// +// +Bool_t MRanForestGrow::Process() +{ + Bool_t not_last=fRanForest->GrowForest(); + fRanTree->SetReadyToSave(); + + return not_last; +} +Bool_t MRanForestGrow::PostProcess() +{ + fRanTree->SetReadyToSave(); + fRanForest->SetReadyToSave(); + + return kTRUE; +} +void MRanForestGrow::SetNumTrees(Int_t n) +{ + fNumTrees=n; + return; +} + +void MRanForestGrow::SetNumTry(Int_t n) +{ + fNumTry=n; + return; +} + +void MRanForestGrow::SetNdSize(Int_t n) +{ + fNdSize=n; + return; +} Index: trunk/MagicSoft/Mars/manalysis/MRanForestGrow.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MRanForestGrow.h (revision 1859) +++ trunk/MagicSoft/Mars/manalysis/MRanForestGrow.h (revision 1859) @@ -0,0 +1,52 @@ +#ifndef MARS_MRanForestGrow +#define MARS_MRanForestGrow + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#include "TRandom.h" +#include "TMath.h" + +///////////////////////////////////////////////////////////////////////////// +// // +// MRanForestGrow // +// // +// Task to grow (train) a random forest // +// // +///////////////////////////////////////////////////////////////////////////// + +#include "MRanForest.h" + +class MHMatrix; +class MParList; +class MRanForest; + +class MRanForestGrow : public MTask +{ +private: + MRanTree *fRanTree; + MRanForest *fRanForest; + MHMatrix *fMGammas; //! Gammas describing matrix + MHMatrix *fMHadrons; //! Hadrons (non gammas) describing matrix + + Int_t fNumTrees; + Int_t fNumTry; + Int_t fNdSize; + +public: + MRanForestGrow(const char *name=NULL, const char *title=NULL); + + void SetNumTrees(Int_t n); + void SetNumTry(Int_t n); + void SetNdSize(Int_t n); + + Bool_t PreProcess(MParList *pList); + Bool_t Process(); + Bool_t PostProcess(); + + ClassDef(MRanForestGrow, 0) // Task to grow a random forest +}; + +#endif + Index: trunk/MagicSoft/Mars/manalysis/MRanTree.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MRanTree.cc (revision 1859) +++ trunk/MagicSoft/Mars/manalysis/MRanTree.cc (revision 1859) @@ -0,0 +1,534 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MRanTree // +// // +// ParameterContainer for Tree structure // +// // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MRanTree.h" +#include +ClassImp(MRanTree); + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MRanTree::MRanTree(const char *name, const char *title):fData(NULL) +{ + fName = name ? name : "MRanTree"; + fTitle = title ? title : "Storage container for structure of a single tree and additional information"; +} + +void MRanTree::SetNdSize(Int_t n) +{ + // minimum nodesize of terminal nodes + fNdSize=TMath::Max(1,n);//at least 1 event per node + + return; +} + +void MRanTree::SetNumTry(Int_t n) +{ + // number of trials in random split selection: + // choose at least 1 variable to split in... + fNumTry=TMath::Max(1,n); + // and not more candidates than available + if(fData) + fNumTry=TMath::Min(fData->GetNumEntries(),n); + + return; +} + +void MRanTree::GrowTree(TMatrix &mhad,TMatrix &mgam,Int_t numdata, Int_t numdim,TArrayI &hadtrue,TArrayI &datasort, + TArrayI &datarang,TArrayF &ginidec,TArrayF &tclasspop,TArrayI &jinbag, + TArrayF &winbag,TArrayF &weight,TRandom &rand) +{ + //Int_t nrnodes=2*(numdata/fNdSize)+1; + Int_t nrnodes=2*numdata+1; + Int_t ninbag=0; + for (Int_t n=0;n1.0e-5) + { + crit=(rln/rld)+(rrn/rrd); + if (crit>critvar) + { + nbestvar=nsp; + critvar=crit; + } + } + } + } + + if (critvar>critmax) { + msplit=mvar; + nbest=nbestvar; + critmax=critvar; + } + } + + decsplit=critmax-crit0; + if (critmax<-1.0e10) jstat=1; + + return jstat; +} + +void MRanTree::MoveData(TArrayI &datasort,Int_t mdim,Int_t numdata,Int_t ndstart, + Int_t ndend,TArrayI &idmove,TArrayI &ncase,Int_t msplit, + Int_t nbest,Int_t &ndendl) +{ + // This is the heart of the BuildTree construction. Based on the best split + // the data in the part of datasort corresponding to the current node is moved to the + // left if it belongs to the left child and right if it belongs to the right child-node. + + Int_t nc,k,ih; + TArrayI tdatasort(numdata); + + // compute idmove = indicator of case nos. going left + + for (Int_t nsp=ndstart;nsp<=nbest;nsp++) + { + nc=datasort[msplit*numdata+nsp]; + idmove[nc]=1; + } + for (Int_t nsp=nbest+1;nsp<=ndend;nsp++) + { + nc=datasort[msplit*numdata+nsp]; + idmove[nc]=0; + } + ndendl=nbest; + + // shift case. nos. right and left for numerical variables. + + for(Int_t msh=0;mshncur) break; + if (nodestatus[kbuild]!=2) continue; + + // initialize for next call to FindBestSplit + + ndstart=nodestart[kbuild]; + ndend=ndstart+nodepop[kbuild]-1; + for (Int_t j=0;j<2;j++) + { + tclasspop[j]=classpop[j*nrnodes+kbuild]; + } + + jstat=FindBestSplit(datasort,datarang,hadtrue,mdim,numdata, + ndstart,ndend,tclasspop,msplit,decsplit, + nbest,ncase,jinbag,iv,winbag,wr,wc,wl, + kbuild,rand); + + if(jstat==1) { + nodestatus[kbuild]=-1; + continue; + }else{ + fBestVar[kbuild]=msplit; + ginidec[msplit]+=decsplit; + + bestsplit[kbuild]=datasort[msplit*numdata+nbest]; + bestsplitnext[kbuild]=datasort[msplit*numdata+nbest+1]; + } + + MoveData(datasort,mdim,numdata,ndstart,ndend,idmove,ncase, + msplit,nbest,ndendl); + + // leftnode no.= ncur+1, rightnode no. = ncur+2. + + nodepop[ncur+1]=ndendl-ndstart+1; + nodepop[ncur+2]=ndend-ndendl; + nodestart[ncur+1]=ndstart; + nodestart[ncur+2]=ndendl+1; + + // find class populations in both nodes + + for (Int_t n=ndstart;n<=ndendl;n++) + { + nc=ncase[n]; + Int_t j=hadtrue[nc]; + classpop[j*nrnodes+ncur+1]+=winbag[nc]; + } + + for (Int_t n=ndendl+1;n<=ndend;n++) + { + nc=ncase[n]; + Int_t j=hadtrue[nc]; + classpop[j*nrnodes+ncur+2]+=winbag[nc]; + + } + + // check on nodestatus + + nodestatus[ncur+1]=2; + nodestatus[ncur+2]=2; + if (nodepop[ncur+1]<=fNdSize) nodestatus[ncur+1]=-1; + if (nodepop[ncur+2]<=fNdSize) nodestatus[ncur+2]=-1; + popt1=0; + popt2=0; + for (Int_t j=0;j<2;j++) + { + popt1+=classpop[j*nrnodes+ncur+1]; + popt2+=classpop[j*nrnodes+ncur+2]; + } + + for (Int_t j=0;j<2;j++) + { + if (classpop[j*nrnodes+ncur+1]==popt1) nodestatus[ncur+1]=-1; + if (classpop[j*nrnodes+ncur+2]==popt2) nodestatus[ncur+2]=-1; + } + + fTreeMap1[kbuild]=ncur+1; + fTreeMap2[kbuild]=ncur+2; + parent[ncur+1]=kbuild; + parent[ncur+2]=kbuild; + nodestatus[kbuild]=1; + ncur+=2; + if (ncur>=nrnodes) break; + } + + // determine number of nodes + fNumNodes=nrnodes; + for (Int_t k=nrnodes-1;k>=0;k--) + { + if (nodestatus[k]==0) fNumNodes-=1; + if (nodestatus[k]==2) nodestatus[k]=-1; + } + + fNumEndNodes=0; + for (Int_t kn=0;knpp) + { + // class + status of node kn coded into fBestVar[kn] + fBestVar[kn]=j-2; + pp=classpop[j*nrnodes+kn]; + } + } + fBestSplit[kn] =classpop[1*nrnodes+kn]; + fBestSplit[kn]/=(classpop[0*nrnodes+kn]+classpop[1*nrnodes+kn]); + } + } + return; +} + +void MRanTree::SetRules(MDataArray *rules) +{ + fData=rules; + + return; +} + +Double_t MRanTree::TreeHad(TVector &event) +{ + Int_t kt=0; + // to optimize on storage space node status and node class + // are coded into fBestVar: + // status of node kt = TMath::Sign(1,fBestVar[kt]) + // hadronness assigned to node kt = fBestSplit[kt] + + for (Int_t k=0;kGetNumEntries(); + TVector event(ncols); + + for (int i=0; iGetNumEntries(); i++) + event(i) = (*fData)(i); + + Int_t kt=0; + // to optimize on storage space node status and node class + // are coded into fBestVar: + // status of node kt = TMath::Sign(1,fBestVar[kt]) + // class of node kt = fBestVar[kt]+2 (class defined by larger + // node population, actually not used) + // hadronness assigned to node kt = fBestSplit[kt] + + for (Int_t k=0;k +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MHRanForest +// +// This histogram shows the evolution of the standard deviation +// of est. hadronness as the number of trees increases. +// It helps you to find out how many trees are really needed in g/h-sep. +// +//////////////////////////////////////////////////////////////////////////// +#include "MHRanForest.h" + +#include +#include +#include +#include +#include + +#include "MParList.h" +#include "MBinning.h" +#include "MRanForest.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MMcEvt.hxx" + +ClassImp(MHRanForest); + +// -------------------------------------------------------------------------- +// +// Setup histograms, nbins is the number of bins used for the evaluation. +// The default is 100 bins. +// +MHRanForest::MHRanForest(Int_t nbins, const char *name, const char *title) +{ + // + // set the name and title of this object + // + fName = name ? name : "MHRanForest"; + fTitle = title ? title : "Histogram showing Standard deviation of estimated hadronness"; + + fGraphSigma = new TGraph; + fGraphSigma->SetTitle("Evolution of Standard deviation of estimated hadronness in tree combination"); + fGraphSigma->SetMaximum(1); +} + +// -------------------------------------------------------------------------- +// +// Delete the histograms. +// +MHRanForest::~MHRanForest() +{ + delete fGraphSigma; +} + +// -------------------------------------------------------------------------- +// +// Setup Filling of the histograms. It needs: +// MMcEvt and MRanForest +// +Bool_t MHRanForest::SetupFill(const MParList *plist) +{ + fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt"); + if (!fMcEvt) + { + *fLog << err << dbginf << "MMcEvt not found... aborting." << endl; + return kFALSE; + } + + fRanForest = (MRanForest*)plist->FindObject("MRanForest"); + if (!fRanForest) + { + *fLog << err << dbginf << "MRanForest not found... aborting." << endl; + return kFALSE; + } + + fSigma.Set(fRanForest->GetNumTrees()); + fNumEvent=0; + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// +Bool_t MHRanForest::Fill(const MParContainer *par) +{ + fNumEvent++; + Double_t hest=0; + Double_t htrue=fMcEvt->GetPartId()==kGAMMA ? 0. : 1.; + + Int_t ntrees=fRanForest->GetNumTrees(); + + for (Int_t i=0;iGetTreeHad(j); + + hest/=Double_t(i+1); + fSigma[i]+=(htrue-hest)*(htrue-hest); + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Finalize the histogram: +// calculate standard deviation and set histogram max and min +// +Bool_t MHRanForest::Finalize() +{ + Int_t n = fSigma.GetSize(); + + fGraphSigma->Set(n); + + Stat_t max=0.; + Stat_t min=0.; + for (Int_t i=0; iSetPoint(i,ip,ig); + } + fGraphSigma->SetMaximum(1.05*max); + fGraphSigma->SetMinimum(0.95*min); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Draw clone of histogram +// +TObject *MHRanForest::DrawClone(Option_t *opt) const +{ + if (fGraphSigma->GetN()==0) + return NULL; + + TCanvas &c = *MakeDefCanvas("RanForest", fTitle); + gROOT->SetSelectedPad(NULL); + + gStyle->SetOptStat(10); + TGraph &g = (TGraph&)*fGraphSigma->DrawClone("AL"); + g.SetBit(kCanDelete); + gPad->Modified(); + gPad->Update(); + if (g.GetHistogram()) + { + g.GetXaxis()->SetRangeUser(0, fNumEvent); + g.GetXaxis()->SetTitle("Number of Trees"); + g.GetYaxis()->SetTitle("Standard deviation of estimated hadronness"); + g.SetMarkerStyle(kFullDotlarge); + g.Draw("P"); + + } + gPad->SetGrid(); + + return &c; +} + +// -------------------------------------------------------------------------- +// +// Draw histogram. (For the Meaning see class description) +// +void MHRanForest::Draw(Option_t *) +{ + if (fGraphSigma->GetN()==0) + return; + + if (!gPad) + MakeDefCanvas("RanForest", fTitle); + + gStyle->SetOptStat(10); + fGraphSigma->Draw("ALP"); + gPad->Modified(); + gPad->Update(); + if (fGraphSigma->GetHistogram()) + { + fGraphSigma->GetXaxis()->SetRangeUser(0, 1); + fGraphSigma->GetXaxis()->SetTitle("Number of Trees"); + fGraphSigma->GetYaxis()->SetTitle("Standard deviation of estimated hadronness"); + + fGraphSigma->SetMarkerStyle(kFullDotSmall); + fGraphSigma->Draw("P"); + gPad->Modified(); + gPad->Update(); + } +} +void MHRanForest::Print(Option_t *) const +{ + return; +} Index: trunk/MagicSoft/Mars/mhist/MHRanForest.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHRanForest.h (revision 1859) +++ trunk/MagicSoft/Mars/mhist/MHRanForest.h (revision 1859) @@ -0,0 +1,46 @@ +#ifndef MARS_MHRanForest +#define MARS_MHRanForest + +#ifndef MARS_MH +#include "MH.h" +#endif + +#ifndef TArrayF +#include "TArrayF.h" +#endif + +class TH1D; +class TGraph; +class MParList; +class MMcEvt; +class MRanForest; + +class MHRanForest : public MH +{ +private: + const MMcEvt *fMcEvt; + const MRanForest *fRanForest; + + TArrayF fSigma; + Int_t fNumEvent; + TGraph *fGraphSigma; + +public: + MHRanForest(Int_t nbins=100, const char *name=NULL, const char *title=NULL); + ~MHRanForest(); + + TGraph *GetGrphSigma() const { return fGraphSigma; } + + Bool_t SetupFill(const MParList *plist); + Bool_t Fill(const MParContainer *par); + Bool_t Finalize(); + + void Print(Option_t *option="") const; + + void Draw(Option_t *opt=""); + TObject *DrawClone(Option_t *opt="") const; + + ClassDef(MHRanForest, 1) // Histogram showing variance of estimated Hadronness +}; + +#endif Index: trunk/MagicSoft/Mars/mhist/MHRanForestGini.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHRanForestGini.cc (revision 1859) +++ trunk/MagicSoft/Mars/mhist/MHRanForestGini.cc (revision 1859) @@ -0,0 +1,207 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MHRanForest +// +// This histogram shows a measure of variable importance (mean decrease in +// Gini-index) +// +//////////////////////////////////////////////////////////////////////////// +#include "MHRanForestGini.h" + +#include +#include +#include +#include +#include + +#include "MParList.h" +#include "MBinning.h" +#include "MRanForest.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MMcEvt.hxx" + +ClassImp(MHRanForestGini); + +// -------------------------------------------------------------------------- +// +// Setup histograms, nbins is the number of bins used for the evaluation. +// The default is 100 bins. +// +MHRanForestGini::MHRanForestGini(Int_t nbins, const char *name, const char *title) +{ + // + // set the name and title of this object + // + fName = name ? name : "MHRanForestGini"; + fTitle = title ? title : "Measure of importance of Random Forest-input parameters"; + + fGraphGini = new TGraph; + fGraphGini->SetTitle("Importance of RF-input parameters measured by mean Gini decrease"); + fGraphGini->SetMaximum(1); +} + +// -------------------------------------------------------------------------- +// +// Delete the histograms. +// +MHRanForestGini::~MHRanForestGini() +{ + delete fGraphGini; +} + +// -------------------------------------------------------------------------- +// +// Setup Filling of the histograms. It needs: +// MMcEvt and MRanForest +// +Bool_t MHRanForestGini::SetupFill(const MParList *plist) +{ + fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt"); + if (!fMcEvt) + { + *fLog << err << dbginf << "MMcEvt not found... aborting." << endl; + return kFALSE; + } + + fRanForest = (MRanForest*)plist->FindObject("MRanForest"); + if (!fRanForest) + { + *fLog << err << dbginf << "MRanForest not found... aborting." << endl; + return kFALSE; + } + + fGini.Set(fRanForest->GetNumDim()); + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Fill the RanForest from a MRanForest container into the corresponding +// histogram dependant on the particle id. +// +// +Bool_t MHRanForestGini::Fill(const MParContainer *par) +{ + for (Int_t i=0;iGetNumDim();i++) + fGini[i]+=fRanForest->GetGiniDec(i); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// +Bool_t MHRanForestGini::Finalize() +{ + Int_t n = fGini.GetSize(); + + fGraphGini->Set(n); + + Stat_t max=0.; + Stat_t min=0.; + for (Int_t i=0; iGetNumTrees()); + fGini[i]/=Stat_t(fRanForest->GetNumData()); + Stat_t ig = fGini[i]; + max=TMath::Max(max,ig); + min=TMath::Min(min,ig); + fGraphGini->SetPoint(i,ip,ig); + } + fGraphGini->SetMaximum(1.05*max); + fGraphGini->SetMinimum(0.95*min); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Draw clone of histogram (For the Meaning see class description) +// +TObject *MHRanForestGini::DrawClone(Option_t *opt) const +{ + if (fGraphGini->GetN()==0) + return NULL; + + TCanvas &c = *MakeDefCanvas("RanForestGini", fTitle); + gROOT->SetSelectedPad(NULL); + + gStyle->SetOptStat(10); + TGraph &g = (TGraph&)*fGraphGini->DrawClone("AL"); + g.SetBit(kCanDelete); + gPad->Modified(); + gPad->Update(); + if (g.GetHistogram()) + { + g.GetXaxis()->SetRangeUser(0, fRanForest->GetNumTrees()); + g.GetXaxis()->SetTitle("No. of RF-input parameter parameter"); + g.GetYaxis()->SetTitle("Mean decrease in Gini-index [a.u.]"); + g.SetMarkerStyle(kFullDotlarge); + g.Draw("P"); + + } + gPad->SetGrid(); + + return &c; +} + +// -------------------------------------------------------------------------- +// +// Draw histogram. (For the Meaning see class description) +// +void MHRanForestGini::Draw(Option_t *) +{ + if (fGraphGini->GetN()==0) + return; + + if (!gPad) + MakeDefCanvas("RanForest", fTitle); + + gStyle->SetOptStat(10); + fGraphGini->Draw("ALP"); + gPad->Modified(); + gPad->Update(); + if (fGraphGini->GetHistogram()) + { + fGraphGini->GetXaxis()->SetRangeUser(0, 1); + fGraphGini->GetXaxis()->SetTitle("No. of parameter"); + fGraphGini->GetYaxis()->SetTitle("Mean decrease in Gini-index [a.u.]"); + + fGraphGini->SetMarkerStyle(kFullDotSmall); + fGraphGini->Draw("P"); + gPad->Modified(); + gPad->Update(); + } +} +void MHRanForestGini::Print(Option_t *) const +{ + return; +} Index: trunk/MagicSoft/Mars/mhist/MHRanForestGini.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHRanForestGini.h (revision 1859) +++ trunk/MagicSoft/Mars/mhist/MHRanForestGini.h (revision 1859) @@ -0,0 +1,45 @@ +#ifndef MARS_MHRanForestGini +#define MARS_MHRanForestGini + +#ifndef MARS_MH +#include "MH.h" +#endif + +#ifndef TArrayF +#include "TArrayF.h" +#endif + +class TH1D; +class TGraph; +class MParList; +class MMcEvt; +class MRanForest; + +class MHRanForestGini : public MH +{ +private: + const MMcEvt *fMcEvt; + const MRanForest *fRanForest; + + TArrayF fGini; + TGraph *fGraphGini; + +public: + MHRanForestGini(Int_t nbins=100, const char *name=NULL, const char *title=NULL); + ~MHRanForestGini(); + + TGraph *GetGrphGini() const { return fGraphGini; } + + Bool_t SetupFill(const MParList *plist); + Bool_t Fill(const MParContainer *par); + Bool_t Finalize(); + + void Print(Option_t *option="") const; + + void Draw(Option_t *opt=""); + TObject *DrawClone(Option_t *opt="") const; + + ClassDef(MHRanForestGini, 1) +}; + +#endif Index: trunk/MagicSoft/Mars/mhist/Makefile =================================================================== --- trunk/MagicSoft/Mars/mhist/Makefile (revision 1858) +++ trunk/MagicSoft/Mars/mhist/Makefile (revision 1859) @@ -52,4 +52,6 @@ MHCompProb.cc \ MHHadronness.cc \ + MHRanForest.cc \ + MHRanForestGini.cc \ MHMcEnergy.cc \ MHMcEfficiency.cc \