Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2070) +++ trunk/MagicSoft/Mars/Changelog (revision 2071) @@ -1,3 +1,20 @@ -*-*- END OF LINE -*-*- + + + 2003/05/06: Thomas Bretz + + * manalysis/MRan*, mhist/MHRan*: + - moved to new directory mranforest + + * mranforest, mranforest/Makefile, mranforest/RanForestLinkDef.h, + mranforest/RanForestIncl.h: + - added + + * manalysis/Makefile, manalysis/AnalysisLinkDef.h: + - removed MRan* + + * mhist/Makefile, mhist/HistLinkDef.h: + - removed MHRan* + Index: trunk/MagicSoft/Mars/Makefile =================================================================== --- trunk/MagicSoft/Mars/Makefile (revision 2070) +++ trunk/MagicSoft/Mars/Makefile (revision 2071) @@ -48,10 +48,9 @@ mraw \ mgui \ + mranforest \ mdata \ mhistmc \ mfilter \ mtools - - LIBRARIES = $(SUBDIRS:=.a) Index: trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 2070) +++ trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 2071) @@ -22,10 +22,4 @@ #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+; @@ -57,14 +51,2 @@ #endif - - - - - - - - - - - - Index: trunk/MagicSoft/Mars/manalysis/MRanForest.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MRanForest.cc (revision 2070) +++ (revision ) @@ -1,413 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 "MHMatrix.h" -#include "MRanTree.h" - -#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"; - - fForest=new TObjArray(); - fForest->SetOwner(kTRUE); -} - -// -------------------------------------------------------------------------- -// -// Destructor. -// -MRanForest::~MRanForest() -{ - delete fForest; -} - -void MRanForest::SetNumTrees(Int_t n) -{ - //at least 1 tree - fNumTrees=TMath::Max(n,1); - fTreeHad.Set(fNumTrees); - fTreeHad.Reset(); -} - -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; -} - -Double_t MRanForest::CalcHadroness(TVector &event) -{ - Double_t hadroness=0; - Int_t ntree=0; - MRanTree *tree; - - TIter forest(fForest); - forest.Reset(); - - while ((tree=(MRanTree*)forest.Next())) - { - fTreeHad[ntree]=tree->TreeHad(event); - hadroness+=fTreeHad[ntree]; - ntree++; - } - return hadroness/Double_t(ntree); -} - -Bool_t MRanForest::AddTree(MRanTree *rantree=NULL) -{ - if (rantree) - fRanTree=rantree; - if (!fRanTree) - return kFALSE; - - fForest->Add((MRanTree*)fRanTree->Clone()); - - 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(); - - 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;nRndm()); - - classpopw[fHadTrue[k]]+=fWeight[k]; - winbag[k]+=fWeight[k]; - jinbag[k]=1; - } - - // modifying sorted-data array for in-bag data: - // - // In bagging procedure ca. 2/3 of all elements in the original - // training sample are used to build the in-bag data - datsortinbag=fDataSort; - - ModifyDataSort(datsortinbag,ninbag,jinbag); - - // growing single tree - fRanTree->GrowTree(hadrons,gammas,fNumData,fNumDim,fHadTrue,datsortinbag, - fDataRang,classpopw,jinbag,winbag,fWeight); - - // 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 2070) +++ (revision ) @@ -1,105 +1,0 @@ -#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 - -class MHMatrix; -class MRanTree; -class TVector; - -class MRanForest : public MParContainer -{ -private: - Int_t fNumTrees; - Int_t fTreeNo; - - 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; - -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); - ~MRanForest(); - - // initialize forest - void SetPriors(Float_t prior_had, Float_t prior_gam); - void SetNumTrees(Int_t n); - - // tree growing - //void SetupForest(); - Bool_t SetupGrow(MHMatrix *mhad,MHMatrix *mgam); - Bool_t GrowForest(); - void SetCurTree(MRanTree *rantree) { fRanTree=rantree; } - Bool_t AddTree(MRanTree *rantree); - - // getter methods - TObjArray *GetForest() { return fForest; } - MRanTree *GetCurTree() { return fRanTree; } - 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); } - - // 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 2070) +++ (revision ) @@ -1,150 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 "MHMatrix.h" // must be before MLogManip.h -#include "MDataArray.h" - -#include "MLog.h" -#include "MLogManip.h" - -#include "MParList.h" - -#include "MRanTree.h" -#include "MRanForest.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), fHadronnessName("MHadronness"), 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; - } - - /*if(!fRanForest->GetCurTree()) - { - *fLog << err << dbginf << "MRanForest does not contain trees... aborting." << endl; - return kFALSE; - }*/ - - fData = fRanTree->GetRules(); - - 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(fHadronnessName, "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 2070) +++ (revision ) @@ -1,41 +1,0 @@ -#ifndef MARS_MRanForestCalc -#define MARS_MRanForestCalc - -#ifndef MARS_MTask -#include "MTask.h" -#endif - -class MParList; -class MHadronness; -class MDataArray; -class MRanTree; -class MRanForest; - -class MRanForestCalc : public MTask -{ -private: - Int_t fNum; // number of trees used to compute hadronness - - TString fHadronnessName; // Name of container storing 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 SetHadronnessName(const TString name) { fHadronnessName = name; } - TString GetHadronnessName() const { return fHadronnessName; } - - 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 2070) +++ (revision ) @@ -1,111 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 "MLog.h" -#include "MLogManip.h" - -#include "MParList.h" - -#include "MRanForest.h" - -ClassImp(MRanForestFill); - -static const TString gsDefName = "MRanForestFill"; -static const TString gsDefTitle = "Tree Classification Loop"; - -// -------------------------------------------------------------------------- -// -// -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; - } - - fNum=0; - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// -Bool_t MRanForestFill::Process() -{ - fNum++; - if(!(fRanForest->AddTree(fRanTree))) - return kFALSE; - - return fNumSetNumTrees(fNum); - - return kTRUE; -} - Index: trunk/MagicSoft/Mars/manalysis/MRanForestFill.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MRanForestFill.h (revision 2070) +++ (revision ) @@ -1,35 +1,0 @@ -#ifndef MARS_MRanForestFill -#define MARS_MRanForestFill - -#ifndef MARS_MTask -#include "MTask.h" -#endif - -class MRanTree; -class MRanForest; -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 2070) +++ (revision ) @@ -1,136 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 "MHMatrix.h" // must be before MLogManip.h - -#include "MLog.h" -#include "MLogManip.h" - -#include "MParList.h" - -#include "MRanTree.h" -#include "MRanForest.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->SetCurTree(fRanTree); - fRanForest->SetNumTrees(fNumTrees); - - 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; -} Index: trunk/MagicSoft/Mars/manalysis/MRanForestGrow.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MRanForestGrow.h (revision 2070) +++ (revision ) @@ -1,48 +1,0 @@ -#ifndef MARS_MRanForestGrow -#define MARS_MRanForestGrow - -///////////////////////////////////////////////////////////////////////////// -// // -// MRanForestGrow // -// // -// Task to grow (train) a random forest // -// // -///////////////////////////////////////////////////////////////////////////// - -#ifndef MARS_MTask -#include "MTask.h" -#endif - -class MHMatrix; -class MParList; -class MRanForest; -class MRanTree; - -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){ fNumTrees=n;} - void SetNumTry(Int_t n) { fNumTry=n; } - void SetNdSize(Int_t n) { fNdSize=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 2070) +++ (revision ) @@ -1,549 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 - -#include -#include -#include - -#include "MDataArray.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MRanTree); - -// -------------------------------------------------------------------------- -// -// Default constructor. -// -MRanTree::MRanTree(const char *name, const char *title):fNdSize(0), fNumTry(3), fData(NULL) -{ - - fName = name ? name : "MRanTree"; - fTitle = title ? title : "Storage container for structure of a single tree"; -} - -void MRanTree::SetNdSize(Int_t n) -{ - // threshold nodesize of terminal nodes, i.e. the training data is splitted - // until there is only pure date in the subsets(=terminal nodes) or the - // subset size is LE n - - fNdSize=TMath::Max(1,n);//at least 1 event per node -} - -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); -} - -void MRanTree::GrowTree(TMatrix &mhad,TMatrix &mgam,Int_t numdata, Int_t numdim,TArrayI &hadtrue, - TArrayI &datasort,TArrayI &datarang,TArrayF &tclasspop,TArrayI &jinbag, - TArrayF &winbag,TArrayF &weight) -{ - // arrays have to be initialized with generous size, so number of total nodes (nrnodes) - // is estimated for worst case - Int_t nrnodes=2*numdata+1; - - // number of events in bootstrap sample - Int_t ninbag=0; - for (Int_t n=0;nRndm()); - - // Gini index = rrn/rrd+rln/rld - rrn=pno; - rrd=pdo; - rln=0; - rld=0; - wl.Reset(); - - for (Int_t j=0;j<2;j++) - { - wr[j]=tclasspop[j]; - } - - critvar=-1.0e20; - - for(Int_t nsp=ndstart;nsp<=ndend-1;nsp++) - { - nc=datasort[mvar*numdata+nsp]; - - u=winbag[nc]; - k=hadtrue[nc]; - - rln=rln+u*(2*wl[k]+u); - rrn=rrn+u*(-2*wr[k]+u); - rld=rld+u; - rrd=rrd-u; - - wl[k]=wl[k]+u; - wr[k]=wr[k]-u; - - if (datarang[mvar*numdata+nc]1.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); - - if(jstat==1) { - nodestatus[kbuild]=-1; - continue; - }else{ - fBestVar[kbuild]=msplit; - fGiniDec[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; -} - -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 -#endif - -#ifndef ROOT_TArrayF -#include -#endif - -class TMatrix; -class TVector; -class TRandom; -class MDataArray; - -class MRanTree : public MParContainer -{ -private: - Int_t fNdSize; - Int_t fNumTry; - - Int_t fNumNodes; - Int_t fNumEndNodes; - MDataArray *fData; - - TArrayI fBestVar; - TArrayI fTreeMap1; - TArrayI fTreeMap2; - TArrayF fBestSplit; - - TArrayF fGiniDec; - -public: - MRanTree(const char *name=NULL, const char *title=NULL); - - void SetNdSize(Int_t n); - void SetNumTry(Int_t n); - void SetRules(MDataArray *rules); - - MDataArray *GetRules() { return fData;} - - Int_t GetNdSize() const { return fNdSize; } - Int_t GetNumTry() const { return fNumTry; } - Int_t GetNumNodes() const { return fNumNodes; } - Int_t GetNumEndNodes() const { return fNumEndNodes; } - - Int_t GetBestVar(Int_t i) const { return fBestVar.At(i); } - Int_t GetTreeMap1(Int_t i) const { return fTreeMap1.At(i); } - Int_t GetTreeMap2(Int_t i) const { return fTreeMap2.At(i); } - Int_t GetNodeClass(Int_t i) const { return fBestVar.At(i)+2; } - Int_t GetNodeStatus(Int_t i) const { return TMath::Sign(1,fBestVar.At(i));} - Float_t GetBestSplit(Int_t i)const { return fBestSplit.At(i); } - - Float_t GetGiniDec(Int_t i) const { return fGiniDec.At(i); } - - // functions used in tree growing process - void GrowTree(TMatrix &mhad,TMatrix &mgam,Int_t numdata, Int_t numdim,TArrayI &hadtrue, - TArrayI &datasort,TArrayI &datarang,TArrayF &tclasspop,TArrayI &jinbag, - TArrayF &winbag,TArrayF &weight); - - Int_t FindBestSplit(TArrayI &datasort,TArrayI &datarang,TArrayI &hadtrue,Int_t mdim, - Int_t numdata,Int_t ndstart,Int_t ndend,TArrayF &tclasspop, - Int_t &msplit,Float_t &decsplit,Int_t &nbest,TArrayI &ncase, - TArrayI &jinbag,TArrayI &iv,TArrayF &winbag,TArrayF &wr, - TArrayF &wc,TArrayF &wl,Int_t kbuild); - - void 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); - - void BuildTree(TArrayI &datasort,TArrayI &datarang,TArrayI &hadtrue,Int_t mdim, - Int_t numdata,TArrayI &bestsplit,TArrayI &bestsplitnext, - TArrayI &nodepop,TArrayI &nodestart,TArrayF &tclasspop, - Int_t nrnodes,TArrayI &idmove,TArrayI &ncase,TArrayI &parent, - TArrayI &jinbag,TArrayI &iv,TArrayF &winbag,TArrayF &wr,TArrayF &wc, - TArrayF &wl,Int_t ninbag); - - Double_t TreeHad(TVector &event); - Double_t TreeHad(); - - Bool_t AsciiWrite(ostream &out) const; - - ClassDef(MRanTree, 1) // Storage container for tree structure -}; - -#endif Index: trunk/MagicSoft/Mars/manalysis/Makefile =================================================================== --- trunk/MagicSoft/Mars/manalysis/Makefile (revision 2070) +++ trunk/MagicSoft/Mars/manalysis/Makefile (revision 2071) @@ -41,9 +41,4 @@ MCompProbCalc.cc \ MMultiDimDistCalc.cc \ - MRanTree.cc \ - MRanForest.cc \ - MRanForestGrow.cc \ - MRanForestCalc.cc \ - MRanForestFill.cc \ MCerPhotPix.cc \ MCerPhotEvt.cc \ Index: trunk/MagicSoft/Mars/mhist/MHRanForest.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHRanForest.cc (revision 2070) +++ (revision ) @@ -1,186 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 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); - fGraphSigma->SetMarkerStyle(kFullDotSmall); -} - -// -------------------------------------------------------------------------- -// -// 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, const Stat_t w) -{ - 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/=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 histogram. (For the Meaning see class description) -// -void MHRanForest::Draw(Option_t *) -{ - if (fGraphSigma->GetN()==0) - return; - - TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); - pad->SetBorderMode(0); - - AppendPad(""); - - fGraphSigma->Draw("ALP"); - pad->Modified(); - pad->Update(); - - TH1 *h=fGraphSigma->GetHistogram(); - if (!h) - return; - - h->GetXaxis()->SetRangeUser(0, 1); - h->SetXTitle("No.of Trees"); - h->SetYTitle("\\sigma of est.hadronness"); - - pad->Modified(); - pad->Update(); -} Index: trunk/MagicSoft/Mars/mhist/MHRanForest.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHRanForest.h (revision 2070) +++ (revision ) @@ -1,43 +1,0 @@ -#ifndef MARS_MHRanForest -#define MARS_MHRanForest - -#ifndef MARS_MH -#include "MH.h" -#endif - -#ifndef ROOT_TArrayF -#include -#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 *GetGraphSigma() const { return fGraphSigma; } - - Bool_t SetupFill(const MParList *plist); - Bool_t Fill(const MParContainer *par, const Stat_t w=1); - Bool_t Finalize(); - - void Draw(Option_t *opt=""); - - ClassDef(MHRanForest, 1) // Histogram showing variance of estimated Hadronness -}; - -#endif Index: trunk/MagicSoft/Mars/mhist/MHRanForestGini.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHRanForestGini.cc (revision 2070) +++ (revision ) @@ -1,169 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 "MRanTree.h" -#include "MRanForest.h" - -#include "MLog.h" -#include "MLogManip.h" - -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); - fGraphGini->SetMarkerStyle(kFullDotSmall); -} - -// -------------------------------------------------------------------------- -// -// Delete the histograms. -// -MHRanForestGini::~MHRanForestGini() -{ - delete fGraphGini; -} - -// -------------------------------------------------------------------------- -// -// Setup Filling of the histograms. It needs: -// MMcEvt and MRanForest -// -Bool_t MHRanForestGini::SetupFill(const MParList *plist) -{ - 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, const Stat_t w) -{ - for (Int_t i=0;iGetNumDim();i++) - fGini[i]+=fRanForest->GetCurTree()->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]/=fRanForest->GetNumData(); - - Stat_t ip = i+1.; - 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 histogram. (For the Meaning see class description) -// -void MHRanForestGini::Draw(Option_t *) -{ - if (fGraphGini->GetN()==0) - return; - - TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); - pad->SetBorderMode(0); - - AppendPad(""); - - fGraphGini->Draw("ALP"); - pad->Modified(); - pad->Update(); - - TH1 *h = fGraphGini->GetHistogram(); - if (!h) - return; - - h->GetXaxis()->SetRangeUser(0, 1); - h->SetXTitle("No.of RF-input parameter"); - h->SetYTitle("Mean decrease in Gini-index [a.u.]"); - - pad->Modified(); - pad->Update(); -} Index: trunk/MagicSoft/Mars/mhist/MHRanForestGini.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHRanForestGini.h (revision 2070) +++ (revision ) @@ -1,41 +1,0 @@ -#ifndef MARS_MHRanForestGini -#define MARS_MHRanForestGini - -#ifndef MARS_MH -#include "MH.h" -#endif - -#ifndef ROOT_TArrayF -#include -#endif - -class TH1D; -class TGraph; -class MParList; -class MRanForest; -class MRanTree; - -class MHRanForestGini : public MH -{ -private: - MRanForest *fRanForest; //! - - TArrayF fGini; //! - TGraph *fGraphGini; //-> - -public: - MHRanForestGini(Int_t nbins=100, const char *name=NULL, const char *title=NULL); - ~MHRanForestGini(); - - TGraph *GetGraphGini() const { return fGraphGini; } - - Bool_t SetupFill(const MParList *plist); - Bool_t Fill(const MParContainer *par, const Stat_t w=1); - Bool_t Finalize(); - - void Draw(Option_t *opt=""); - - ClassDef(MHRanForestGini, 1) -}; - -#endif Index: trunk/MagicSoft/Mars/mranforest/MHRanForest.cc =================================================================== --- trunk/MagicSoft/Mars/mranforest/MHRanForest.cc (revision 2071) +++ trunk/MagicSoft/Mars/mranforest/MHRanForest.cc (revision 2071) @@ -0,0 +1,186 @@ +/* ======================================================================== *\ +! +! * +! * 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 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); + fGraphSigma->SetMarkerStyle(kFullDotSmall); +} + +// -------------------------------------------------------------------------- +// +// 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, const Stat_t w) +{ + 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/=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 histogram. (For the Meaning see class description) +// +void MHRanForest::Draw(Option_t *) +{ + if (fGraphSigma->GetN()==0) + return; + + TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); + pad->SetBorderMode(0); + + AppendPad(""); + + fGraphSigma->Draw("ALP"); + pad->Modified(); + pad->Update(); + + TH1 *h=fGraphSigma->GetHistogram(); + if (!h) + return; + + h->GetXaxis()->SetRangeUser(0, 1); + h->SetXTitle("No.of Trees"); + h->SetYTitle("\\sigma of est.hadronness"); + + pad->Modified(); + pad->Update(); +} Index: trunk/MagicSoft/Mars/mranforest/MHRanForest.h =================================================================== --- trunk/MagicSoft/Mars/mranforest/MHRanForest.h (revision 2071) +++ trunk/MagicSoft/Mars/mranforest/MHRanForest.h (revision 2071) @@ -0,0 +1,43 @@ +#ifndef MARS_MHRanForest +#define MARS_MHRanForest + +#ifndef MARS_MH +#include "MH.h" +#endif + +#ifndef ROOT_TArrayF +#include +#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 *GetGraphSigma() const { return fGraphSigma; } + + Bool_t SetupFill(const MParList *plist); + Bool_t Fill(const MParContainer *par, const Stat_t w=1); + Bool_t Finalize(); + + void Draw(Option_t *opt=""); + + ClassDef(MHRanForest, 1) // Histogram showing variance of estimated Hadronness +}; + +#endif Index: trunk/MagicSoft/Mars/mranforest/MHRanForestGini.cc =================================================================== --- trunk/MagicSoft/Mars/mranforest/MHRanForestGini.cc (revision 2071) +++ trunk/MagicSoft/Mars/mranforest/MHRanForestGini.cc (revision 2071) @@ -0,0 +1,169 @@ +/* ======================================================================== *\ +! +! * +! * 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 "MRanTree.h" +#include "MRanForest.h" + +#include "MLog.h" +#include "MLogManip.h" + +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); + fGraphGini->SetMarkerStyle(kFullDotSmall); +} + +// -------------------------------------------------------------------------- +// +// Delete the histograms. +// +MHRanForestGini::~MHRanForestGini() +{ + delete fGraphGini; +} + +// -------------------------------------------------------------------------- +// +// Setup Filling of the histograms. It needs: +// MMcEvt and MRanForest +// +Bool_t MHRanForestGini::SetupFill(const MParList *plist) +{ + 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, const Stat_t w) +{ + for (Int_t i=0;iGetNumDim();i++) + fGini[i]+=fRanForest->GetCurTree()->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]/=fRanForest->GetNumData(); + + Stat_t ip = i+1.; + 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 histogram. (For the Meaning see class description) +// +void MHRanForestGini::Draw(Option_t *) +{ + if (fGraphGini->GetN()==0) + return; + + TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); + pad->SetBorderMode(0); + + AppendPad(""); + + fGraphGini->Draw("ALP"); + pad->Modified(); + pad->Update(); + + TH1 *h = fGraphGini->GetHistogram(); + if (!h) + return; + + h->GetXaxis()->SetRangeUser(0, 1); + h->SetXTitle("No.of RF-input parameter"); + h->SetYTitle("Mean decrease in Gini-index [a.u.]"); + + pad->Modified(); + pad->Update(); +} Index: trunk/MagicSoft/Mars/mranforest/MHRanForestGini.h =================================================================== --- trunk/MagicSoft/Mars/mranforest/MHRanForestGini.h (revision 2071) +++ trunk/MagicSoft/Mars/mranforest/MHRanForestGini.h (revision 2071) @@ -0,0 +1,41 @@ +#ifndef MARS_MHRanForestGini +#define MARS_MHRanForestGini + +#ifndef MARS_MH +#include "MH.h" +#endif + +#ifndef ROOT_TArrayF +#include +#endif + +class TH1D; +class TGraph; +class MParList; +class MRanForest; +class MRanTree; + +class MHRanForestGini : public MH +{ +private: + MRanForest *fRanForest; //! + + TArrayF fGini; //! + TGraph *fGraphGini; //-> + +public: + MHRanForestGini(Int_t nbins=100, const char *name=NULL, const char *title=NULL); + ~MHRanForestGini(); + + TGraph *GetGraphGini() const { return fGraphGini; } + + Bool_t SetupFill(const MParList *plist); + Bool_t Fill(const MParContainer *par, const Stat_t w=1); + Bool_t Finalize(); + + void Draw(Option_t *opt=""); + + ClassDef(MHRanForestGini, 1) +}; + +#endif Index: trunk/MagicSoft/Mars/mranforest/MRanForest.cc =================================================================== --- trunk/MagicSoft/Mars/mranforest/MRanForest.cc (revision 2071) +++ trunk/MagicSoft/Mars/mranforest/MRanForest.cc (revision 2071) @@ -0,0 +1,413 @@ +/* ======================================================================== *\ +! +! * +! * 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 "MHMatrix.h" +#include "MRanTree.h" + +#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"; + + fForest=new TObjArray(); + fForest->SetOwner(kTRUE); +} + +// -------------------------------------------------------------------------- +// +// Destructor. +// +MRanForest::~MRanForest() +{ + delete fForest; +} + +void MRanForest::SetNumTrees(Int_t n) +{ + //at least 1 tree + fNumTrees=TMath::Max(n,1); + fTreeHad.Set(fNumTrees); + fTreeHad.Reset(); +} + +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; +} + +Double_t MRanForest::CalcHadroness(TVector &event) +{ + Double_t hadroness=0; + Int_t ntree=0; + MRanTree *tree; + + TIter forest(fForest); + forest.Reset(); + + while ((tree=(MRanTree*)forest.Next())) + { + fTreeHad[ntree]=tree->TreeHad(event); + hadroness+=fTreeHad[ntree]; + ntree++; + } + return hadroness/Double_t(ntree); +} + +Bool_t MRanForest::AddTree(MRanTree *rantree=NULL) +{ + if (rantree) + fRanTree=rantree; + if (!fRanTree) + return kFALSE; + + fForest->Add((MRanTree*)fRanTree->Clone()); + + 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(); + + 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;nRndm()); + + classpopw[fHadTrue[k]]+=fWeight[k]; + winbag[k]+=fWeight[k]; + jinbag[k]=1; + } + + // modifying sorted-data array for in-bag data: + // + // In bagging procedure ca. 2/3 of all elements in the original + // training sample are used to build the in-bag data + datsortinbag=fDataSort; + + ModifyDataSort(datsortinbag,ninbag,jinbag); + + // growing single tree + fRanTree->GrowTree(hadrons,gammas,fNumData,fNumDim,fHadTrue,datsortinbag, + fDataRang,classpopw,jinbag,winbag,fWeight); + + // 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/mranforest/MRanForest.h =================================================================== --- trunk/MagicSoft/Mars/mranforest/MRanForest.h (revision 2071) +++ trunk/MagicSoft/Mars/mranforest/MRanForest.h (revision 2071) @@ -0,0 +1,105 @@ +#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 + +class MHMatrix; +class MRanTree; +class TVector; + +class MRanForest : public MParContainer +{ +private: + Int_t fNumTrees; + Int_t fTreeNo; + + 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; + +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); + ~MRanForest(); + + // initialize forest + void SetPriors(Float_t prior_had, Float_t prior_gam); + void SetNumTrees(Int_t n); + + // tree growing + //void SetupForest(); + Bool_t SetupGrow(MHMatrix *mhad,MHMatrix *mgam); + Bool_t GrowForest(); + void SetCurTree(MRanTree *rantree) { fRanTree=rantree; } + Bool_t AddTree(MRanTree *rantree); + + // getter methods + TObjArray *GetForest() { return fForest; } + MRanTree *GetCurTree() { return fRanTree; } + 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); } + + // 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/mranforest/MRanForestCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mranforest/MRanForestCalc.cc (revision 2071) +++ trunk/MagicSoft/Mars/mranforest/MRanForestCalc.cc (revision 2071) @@ -0,0 +1,150 @@ +/* ======================================================================== *\ +! +! * +! * 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 "MHMatrix.h" // must be before MLogManip.h +#include "MDataArray.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MRanTree.h" +#include "MRanForest.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), fHadronnessName("MHadronness"), 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; + } + + /*if(!fRanForest->GetCurTree()) + { + *fLog << err << dbginf << "MRanForest does not contain trees... aborting." << endl; + return kFALSE; + }*/ + + fData = fRanTree->GetRules(); + + 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(fHadronnessName, "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/mranforest/MRanForestCalc.h =================================================================== --- trunk/MagicSoft/Mars/mranforest/MRanForestCalc.h (revision 2071) +++ trunk/MagicSoft/Mars/mranforest/MRanForestCalc.h (revision 2071) @@ -0,0 +1,41 @@ +#ifndef MARS_MRanForestCalc +#define MARS_MRanForestCalc + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class MParList; +class MHadronness; +class MDataArray; +class MRanTree; +class MRanForest; + +class MRanForestCalc : public MTask +{ +private: + Int_t fNum; // number of trees used to compute hadronness + + TString fHadronnessName; // Name of container storing 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 SetHadronnessName(const TString name) { fHadronnessName = name; } + TString GetHadronnessName() const { return fHadronnessName; } + + 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/mranforest/MRanForestFill.cc =================================================================== --- trunk/MagicSoft/Mars/mranforest/MRanForestFill.cc (revision 2071) +++ trunk/MagicSoft/Mars/mranforest/MRanForestFill.cc (revision 2071) @@ -0,0 +1,111 @@ +/* ======================================================================== *\ +! +! * +! * 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 "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MRanForest.h" + +ClassImp(MRanForestFill); + +static const TString gsDefName = "MRanForestFill"; +static const TString gsDefTitle = "Tree Classification Loop"; + +// -------------------------------------------------------------------------- +// +// +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; + } + + fNum=0; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// +Bool_t MRanForestFill::Process() +{ + fNum++; + if(!(fRanForest->AddTree(fRanTree))) + return kFALSE; + + return fNumSetNumTrees(fNum); + + return kTRUE; +} + Index: trunk/MagicSoft/Mars/mranforest/MRanForestFill.h =================================================================== --- trunk/MagicSoft/Mars/mranforest/MRanForestFill.h (revision 2071) +++ trunk/MagicSoft/Mars/mranforest/MRanForestFill.h (revision 2071) @@ -0,0 +1,35 @@ +#ifndef MARS_MRanForestFill +#define MARS_MRanForestFill + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class MRanTree; +class MRanForest; +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/mranforest/MRanForestGrow.cc =================================================================== --- trunk/MagicSoft/Mars/mranforest/MRanForestGrow.cc (revision 2071) +++ trunk/MagicSoft/Mars/mranforest/MRanForestGrow.cc (revision 2071) @@ -0,0 +1,136 @@ +/* ======================================================================== *\ +! +! * +! * 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 "MHMatrix.h" // must be before MLogManip.h + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MRanTree.h" +#include "MRanForest.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->SetCurTree(fRanTree); + fRanForest->SetNumTrees(fNumTrees); + + 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; +} Index: trunk/MagicSoft/Mars/mranforest/MRanForestGrow.h =================================================================== --- trunk/MagicSoft/Mars/mranforest/MRanForestGrow.h (revision 2071) +++ trunk/MagicSoft/Mars/mranforest/MRanForestGrow.h (revision 2071) @@ -0,0 +1,48 @@ +#ifndef MARS_MRanForestGrow +#define MARS_MRanForestGrow + +///////////////////////////////////////////////////////////////////////////// +// // +// MRanForestGrow // +// // +// Task to grow (train) a random forest // +// // +///////////////////////////////////////////////////////////////////////////// + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class MHMatrix; +class MParList; +class MRanForest; +class MRanTree; + +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){ fNumTrees=n;} + void SetNumTry(Int_t n) { fNumTry=n; } + void SetNdSize(Int_t n) { fNdSize=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/mranforest/MRanTree.cc =================================================================== --- trunk/MagicSoft/Mars/mranforest/MRanTree.cc (revision 2071) +++ trunk/MagicSoft/Mars/mranforest/MRanTree.cc (revision 2071) @@ -0,0 +1,549 @@ +/* ======================================================================== *\ +! +! * +! * 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 + +#include +#include +#include + +#include "MDataArray.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MRanTree); + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MRanTree::MRanTree(const char *name, const char *title):fNdSize(0), fNumTry(3), fData(NULL) +{ + + fName = name ? name : "MRanTree"; + fTitle = title ? title : "Storage container for structure of a single tree"; +} + +void MRanTree::SetNdSize(Int_t n) +{ + // threshold nodesize of terminal nodes, i.e. the training data is splitted + // until there is only pure date in the subsets(=terminal nodes) or the + // subset size is LE n + + fNdSize=TMath::Max(1,n);//at least 1 event per node +} + +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); +} + +void MRanTree::GrowTree(TMatrix &mhad,TMatrix &mgam,Int_t numdata, Int_t numdim,TArrayI &hadtrue, + TArrayI &datasort,TArrayI &datarang,TArrayF &tclasspop,TArrayI &jinbag, + TArrayF &winbag,TArrayF &weight) +{ + // arrays have to be initialized with generous size, so number of total nodes (nrnodes) + // is estimated for worst case + Int_t nrnodes=2*numdata+1; + + // number of events in bootstrap sample + Int_t ninbag=0; + for (Int_t n=0;nRndm()); + + // Gini index = rrn/rrd+rln/rld + rrn=pno; + rrd=pdo; + rln=0; + rld=0; + wl.Reset(); + + for (Int_t j=0;j<2;j++) + { + wr[j]=tclasspop[j]; + } + + critvar=-1.0e20; + + for(Int_t nsp=ndstart;nsp<=ndend-1;nsp++) + { + nc=datasort[mvar*numdata+nsp]; + + u=winbag[nc]; + k=hadtrue[nc]; + + rln=rln+u*(2*wl[k]+u); + rrn=rrn+u*(-2*wr[k]+u); + rld=rld+u; + rrd=rrd-u; + + wl[k]=wl[k]+u; + wr[k]=wr[k]-u; + + if (datarang[mvar*numdata+nc]1.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); + + if(jstat==1) { + nodestatus[kbuild]=-1; + continue; + }else{ + fBestVar[kbuild]=msplit; + fGiniDec[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; +} + +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 +#endif + +#ifndef ROOT_TArrayF +#include +#endif + +class TMatrix; +class TVector; +class TRandom; +class MDataArray; + +class MRanTree : public MParContainer +{ +private: + Int_t fNdSize; + Int_t fNumTry; + + Int_t fNumNodes; + Int_t fNumEndNodes; + MDataArray *fData; + + TArrayI fBestVar; + TArrayI fTreeMap1; + TArrayI fTreeMap2; + TArrayF fBestSplit; + + TArrayF fGiniDec; + +public: + MRanTree(const char *name=NULL, const char *title=NULL); + + void SetNdSize(Int_t n); + void SetNumTry(Int_t n); + void SetRules(MDataArray *rules); + + MDataArray *GetRules() { return fData;} + + Int_t GetNdSize() const { return fNdSize; } + Int_t GetNumTry() const { return fNumTry; } + Int_t GetNumNodes() const { return fNumNodes; } + Int_t GetNumEndNodes() const { return fNumEndNodes; } + + Int_t GetBestVar(Int_t i) const { return fBestVar.At(i); } + Int_t GetTreeMap1(Int_t i) const { return fTreeMap1.At(i); } + Int_t GetTreeMap2(Int_t i) const { return fTreeMap2.At(i); } + Int_t GetNodeClass(Int_t i) const { return fBestVar.At(i)+2; } + Int_t GetNodeStatus(Int_t i) const { return TMath::Sign(1,fBestVar.At(i));} + Float_t GetBestSplit(Int_t i)const { return fBestSplit.At(i); } + + Float_t GetGiniDec(Int_t i) const { return fGiniDec.At(i); } + + // functions used in tree growing process + void GrowTree(TMatrix &mhad,TMatrix &mgam,Int_t numdata, Int_t numdim,TArrayI &hadtrue, + TArrayI &datasort,TArrayI &datarang,TArrayF &tclasspop,TArrayI &jinbag, + TArrayF &winbag,TArrayF &weight); + + Int_t FindBestSplit(TArrayI &datasort,TArrayI &datarang,TArrayI &hadtrue,Int_t mdim, + Int_t numdata,Int_t ndstart,Int_t ndend,TArrayF &tclasspop, + Int_t &msplit,Float_t &decsplit,Int_t &nbest,TArrayI &ncase, + TArrayI &jinbag,TArrayI &iv,TArrayF &winbag,TArrayF &wr, + TArrayF &wc,TArrayF &wl,Int_t kbuild); + + void 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); + + void BuildTree(TArrayI &datasort,TArrayI &datarang,TArrayI &hadtrue,Int_t mdim, + Int_t numdata,TArrayI &bestsplit,TArrayI &bestsplitnext, + TArrayI &nodepop,TArrayI &nodestart,TArrayF &tclasspop, + Int_t nrnodes,TArrayI &idmove,TArrayI &ncase,TArrayI &parent, + TArrayI &jinbag,TArrayI &iv,TArrayF &winbag,TArrayF &wr,TArrayF &wc, + TArrayF &wl,Int_t ninbag); + + Double_t TreeHad(TVector &event); + Double_t TreeHad(); + + Bool_t AsciiWrite(ostream &out) const; + + ClassDef(MRanTree, 1) // Storage container for tree structure +}; + +#endif Index: trunk/MagicSoft/Mars/mranforest/Makefile =================================================================== --- trunk/MagicSoft/Mars/mranforest/Makefile (revision 2071) +++ trunk/MagicSoft/Mars/mranforest/Makefile (revision 2071) @@ -0,0 +1,53 @@ +################################################################## +# +# makefile +# +# for the MARS software +# +################################################################## +include ../Makefile.conf.$(OSTYPE) +include ../Makefile.conf.general + +# +# Handling name of the Root Dictionary Files +# +CINT = RanForest + +# +# Library name to creatre +# +LIB = mranforest.a + +# +# connect the include files defined in the config.mk file +# +INCLUDES = -I. -I../mbase -I../mhist -I../mdata -I../manalysis -I../mmc +# MParContainer MH MDataArray MHadronness MMcEvt + +#------------------------------------------------------------------------------ + +.SUFFIXES: .c .cc .cxx .h .hxx .o + +SRCFILES = MRanTree.cc \ + MRanForest.cc \ + MRanForestGrow.cc \ + MRanForestCalc.cc \ + MRanForestFill.cc \ + MHRanForest.cc \ + MHRanForestGini.cc \ + +SRCS = $(SRCFILES) +HEADERS = $(SRCFILES:.cc=.h) +OBJS = $(SRCFILES:.cc=.o) + +############################################################ + +all: $(LIB) + +include ../Makefile.rules + +clean: rmcint rmobjs rmcore rmlib + +mrproper: clean rmbak + +# @endcode Index: trunk/MagicSoft/Mars/mranforest/RanForestIncl.h =================================================================== --- trunk/MagicSoft/Mars/mranforest/RanForestIncl.h (revision 2071) +++ trunk/MagicSoft/Mars/mranforest/RanForestIncl.h (revision 2071) @@ -0,0 +1,3 @@ +#ifndef __CINT__ + +#endif // __CINT__ Index: trunk/MagicSoft/Mars/mranforest/RanForestLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mranforest/RanForestLinkDef.h (revision 2071) +++ trunk/MagicSoft/Mars/mranforest/RanForestLinkDef.h (revision 2071) @@ -0,0 +1,16 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#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 MHRanForest+; +#pragma link C++ class MHRanForestGini+; + +#endif