STAT/TKDPDF.cxx

   1 #include "TKDPDF.h"
   2 #include "TKDNodeInfo.h"
   3
   4 #include "TClonesArray.h"
   5 #include "TTree.h"
   6 #include "TH2.h"
   7 #include "TObjArray.h"
   8 #include "TObjString.h"
   9 #include "TBox.h"
  10 #include "TGraph.h"
  11 #include "TMarker.h"
  12
  13 ClassImp(TKDPDF)
  14
  15
  16
  17 //_________________________________________________________________
  18 TKDPDF::TKDPDF() :
  19   TKDTreeIF()
  20   ,TKDInterpolatorBase()
  21 {
  22 // Default constructor. To be used with care since in this case building
  23 // of data structure is completly left to the user responsability.
  24 }
  25
  26 //_________________________________________________________________
  27 TKDPDF::TKDPDF(Int_t npoints, Int_t ndim, UInt_t bsize, Float_t **data) :
  28   TKDTreeIF(npoints, ndim, bsize, data)
  29   ,TKDInterpolatorBase(ndim)
  30 {
  31 // Wrapper constructor for the TKDTree.
  32
  33   Build();
  34 }
  35
  36
  37 //_________________________________________________________________
  38 TKDPDF::TKDPDF(TTree *t, const Char_t *var, const Char_t *cut, UInt_t bsize, Long64_t nentries, Long64_t firstentry) :
  39   TKDTreeIF()
  40   ,TKDInterpolatorBase()
  41 {
  42 // Alocate data from a tree. The variables which have to be analysed are
  43 // defined in the "var" parameter as a colon separated list. The format should
  44 // be identical to that used by TTree::Draw().
  45 //
  46 //
  47
  48   TObjArray *vars = TString(var).Tokenize(":");
  49   fNDim = vars->GetEntriesFast(); fNDimm = 2*fNDim;
  50   fNSize = fNDim;
  51   if(fNDim > 6/*kDimMax*/) Warning("TKDPDF(TTree*, const Char_t, const Char_t, UInt_t)", "Variable number exceed maximum dimension %d. Results are unpredictable.", 6/*kDimMax*/);
  52   fBucketSize = bsize;
  53
  54   Int_t np;
  55   Double_t *v;
  56   for(int idim=0; idim<fNDim; idim++){
  57     if(!(np = t->Draw(((TObjString*)(*vars)[idim])->GetName(), cut, "goff", nentries, firstentry))){
  58       Warning("TKDPDF(TTree*, const Char_t, const Char_t, UInt_t)", "Can not access data for keys %s. Key defined on tree :", ((TObjString*)(*vars)[idim])->GetName());
  59       TIterator *it = (t->GetListOfLeaves())->MakeIterator();
  60       TObject *o;
  61       while((o = (*it)())) printf("\t%s\n", o->GetName());
  62       continue;
  63     }
  64     if(!fNPoints){
  65       fNPoints = np;
  66       //Info("TKDPDF(TTree*, const Char_t, const Char_t, UInt_t)", Form("Allocating %d data points in %d dimensions.", fNpoints, fNDim));
  67       fData = new Float_t*[fNDim];
  68       for(int jdim=fNDim; jdim--;) fData[jdim] = new Float_t[fNPoints];
  69       fDataOwner = kTRUE;
  70     }
  71     v = t->GetV1();
  72     for(int ip=0; ip<fNPoints; ip++) fData[idim][ip] = (Float_t)v[ip];
  73   }
  74   Build();
  75   delete vars;
  76 }
  77
  78 //_________________________________________________________________
  79 TKDPDF::~TKDPDF()
  80 {
  81 }
  82
  83 //_________________________________________________________________
  84 Bool_t TKDPDF::Build(Int_t)
  85 {
  86 // Fill interpolator's data array i.e.
  87 //  - estimation points
  88 //  - corresponding PDF values
  89
  90   TKDTreeIF::Build();
  91   if(!fBoundaries) MakeBoundaries();
  92   fLambda = 1 + fNDim + (fNDim*(fNDim+1)>>1);
  93   //printf("after MakeBoundaries() %d\n", memory());
  94
  95   // allocate interpolation nodes
  96   Int_t fNTNodes = fNPoints/fBucketSize + ((fNPoints%fBucketSize)?1:0);/*TKDTreeIF::GetNTNodes();*/
  97   TKDInterpolatorBase::Build(fNTNodes);
  98
  99   TKDNodeInfo *node = NULL;
 100   Float_t *bounds = NULL;
 101   Int_t *indexPoints;
 102   for(int inode=0, tnode = fNNodes; inode<fNTNodes-1; inode++, tnode++){
 103     node = (TKDNodeInfo*)(*fNodes)[inode];
 104     node->Val()[0] =  Float_t(fBucketSize)/fNPoints;
 105     bounds = GetBoundary(tnode);
 106     for(int idim=0; idim<fNDim; idim++) node->Val()[0] /= (bounds[2*idim+1] - bounds[2*idim]);
 107     node->Val()[1] =  node->Val()[0]/TMath::Sqrt(float(fBucketSize));
 108
 109     indexPoints = GetPointsIndexes(tnode);
 110     // loop points in this terminal node
 111     for(int idim=0; idim<fNDim; idim++){
 112       node->Data()[idim] = 0.;
 113       for(int ip = 0; ip<fBucketSize; ip++) node->Data()[idim] += fData[idim][indexPoints[ip]];
 114       node->Data()[idim] /= fBucketSize;
 115     }
 116     memcpy(&(node->Data()[fNDim]), bounds, fNDimm*sizeof(Float_t));
 117   }
 118
 119   // analyze last (incomplete) terminal node
 120   Int_t counts = fNPoints%fBucketSize;
 121   counts = counts ? counts : fBucketSize;
 122   Int_t inode = fNTNodes - 1, tnode = inode + fNNodes;
 123   node = (TKDNodeInfo*)(*fNodes)[inode];
 124   node->Val()[0] =  Float_t(counts)/fNPoints;
 125   bounds = GetBoundary(tnode);
 126   for(int idim=0; idim<fNDim; idim++){
 127     Float_t dx = bounds[2*idim+1]-bounds[2*idim];
 128     if(dx < 1.e-30){
 129       Warning("TKDPDF::Build()", "Terminal bucket index[%d] too narrow on the %d dimension.", inode, idim);
 130       continue;
 131     }
 132     node->Val()[0] /= (bounds[2*idim+1] - bounds[2*idim]);
 133   }
 134   node->Val()[1] =  node->Val()[0]/TMath::Sqrt(float(counts));
 135
 136   // loop points in this terminal node
 137   indexPoints = GetPointsIndexes(tnode);
 138   for(int idim=0; idim<fNDim; idim++){
 139     node->Data()[idim] = 0.;
 140     for(int ip = 0; ip<counts; ip++) node->Data()[idim] += fData[idim][indexPoints[ip]];
 141     node->Data()[idim] /= counts;
 142   }
 143   memcpy(&(node->Data()[fNDim]), bounds, fNDimm*sizeof(Float_t));
 144
 145   delete [] fBoundaries;
 146   fBoundaries = NULL;
 147
 148   return kTRUE;
 149 }
 150
 151
 152 //_________________________________________________________________
 153 void TKDPDF::DrawNode(Int_t tnode, UInt_t ax1, UInt_t ax2)
 154 {
 155 // Draw node "node" and the data points within.
 156 //
 157 // Observation:
 158 // This function creates some graphical objects
 159 // but don't delete it. Abusing this function may cause memory leaks !
 160
 161   if(tnode < 0 || tnode >= GetNTNodes()){
 162     Warning("DrawNode()", "Terminal node %d outside defined range.", tnode);
 163     return;
 164   }
 165
 166   Int_t inode = tnode;
 167   tnode += fNNodes;
 168   // select zone of interest in the indexes array
 169   Int_t *index = GetPointsIndexes(tnode);
 170   Int_t nPoints = (tnode == 2*fNNodes) ? fNPoints%fBucketSize : fBucketSize;
 171
 172   // draw data points
 173   TGraph *g = new TGraph(nPoints);
 174   g->SetMarkerStyle(7);
 175   for(int ip = 0; ip<nPoints; ip++) g->SetPoint(ip, fData[ax1][index[ip]], fData[ax2][index[ip]]);
 176
 177   // draw estimation point
 178   TKDNodeInfo *node = (TKDNodeInfo*)(*fNodes)[inode];
 179   TMarker *m=new TMarker(node->Data()[ax1], node->Data()[ax2], 20);
 180   m->SetMarkerColor(2);
 181   m->SetMarkerSize(1.7);
 182
 183   // draw node contour
 184   Float_t *bounds = GetBoundary(tnode);
 185   TBox *n = new TBox(bounds[2*ax1], bounds[2*ax2], bounds[2*ax1+1], bounds[2*ax2+1]);
 186   n->SetFillStyle(0);
 187
 188   g->Draw("ap");
 189   m->Draw();
 190   n->Draw();
 191
 192   return;
 193 }
 194