Coverity fixes
[u/mrichter/AliRoot.git] / 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 }
76
77 //_________________________________________________________________
78 TKDPDF::~TKDPDF()
79 {
80 }
81
82 //_________________________________________________________________
83 Bool_t TKDPDF::Build(Int_t)
84 {
85 // Fill interpolator's data array i.e.
86 //  - estimation points 
87 //  - corresponding PDF values
88
89   TKDTreeIF::Build();
90   if(!fBoundaries) MakeBoundaries();
91   fLambda = 1 + fNDim + (fNDim*(fNDim+1)>>1);
92   //printf("after MakeBoundaries() %d\n", memory());
93   
94   // allocate interpolation nodes
95   Int_t fNTNodes = fNPoints/fBucketSize + ((fNPoints%fBucketSize)?1:0);/*TKDTreeIF::GetNTNodes();*/
96   TKDInterpolatorBase::Build(fNTNodes);
97
98   TKDNodeInfo *node = NULL;
99   Float_t *bounds = NULL;
100   Int_t *indexPoints;
101   for(int inode=0, tnode = fNNodes; inode<fNTNodes-1; inode++, tnode++){
102     node = (TKDNodeInfo*)(*fNodes)[inode];
103     node->Val()[0] =  Float_t(fBucketSize)/fNPoints;
104     bounds = GetBoundary(tnode);
105     for(int idim=0; idim<fNDim; idim++) node->Val()[0] /= (bounds[2*idim+1] - bounds[2*idim]);
106     node->Val()[1] =  node->Val()[0]/TMath::Sqrt(float(fBucketSize));
107     
108     indexPoints = GetPointsIndexes(tnode);
109     // loop points in this terminal node
110     for(int idim=0; idim<fNDim; idim++){
111       node->Data()[idim] = 0.;
112       for(int ip = 0; ip<fBucketSize; ip++) node->Data()[idim] += fData[idim][indexPoints[ip]];
113       node->Data()[idim] /= fBucketSize;
114     }
115     memcpy(&(node->Data()[fNDim]), bounds, fNDimm*sizeof(Float_t));
116   }
117
118   // analyze last (incomplete) terminal node
119   Int_t counts = fNPoints%fBucketSize;
120   counts = counts ? counts : fBucketSize;
121   Int_t inode = fNTNodes - 1, tnode = inode + fNNodes;
122   node = (TKDNodeInfo*)(*fNodes)[inode];
123   node->Val()[0] =  Float_t(counts)/fNPoints;
124   bounds = GetBoundary(tnode);
125   for(int idim=0; idim<fNDim; idim++){ 
126     Float_t dx = bounds[2*idim+1]-bounds[2*idim];
127     if(dx < 1.e-30){
128       Warning("TKDPDF::Build()", "Terminal bucket index[%d] too narrow on the %d dimension.", inode, idim);
129       continue;
130     }
131     node->Val()[0] /= (bounds[2*idim+1] - bounds[2*idim]);
132   }
133   node->Val()[1] =  node->Val()[0]/TMath::Sqrt(float(counts));
134
135   // loop points in this terminal node
136   indexPoints = GetPointsIndexes(tnode);
137   for(int idim=0; idim<fNDim; idim++){
138     node->Data()[idim] = 0.;
139     for(int ip = 0; ip<counts; ip++) node->Data()[idim] += fData[idim][indexPoints[ip]];
140     node->Data()[idim] /= counts;
141   }
142   memcpy(&(node->Data()[fNDim]), bounds, fNDimm*sizeof(Float_t));
143
144   delete [] fBoundaries;
145   fBoundaries = NULL;
146
147   return kTRUE;
148 }
149
150
151 //_________________________________________________________________
152 void TKDPDF::DrawNode(Int_t tnode, UInt_t ax1, UInt_t ax2)
153 {
154 // Draw node "node" and the data points within.
155 //
156 // Observation:
157 // This function creates some graphical objects
158 // but don't delete it. Abusing this function may cause memory leaks !
159
160   if(tnode < 0 || tnode >= GetNTNodes()){
161     Warning("DrawNode()", "Terminal node %d outside defined range.", tnode);
162     return;
163   }
164
165   Int_t inode = tnode;
166   tnode += fNNodes;
167   // select zone of interest in the indexes array
168   Int_t *index = GetPointsIndexes(tnode);
169   Int_t nPoints = (tnode == 2*fNNodes) ? fNPoints%fBucketSize : fBucketSize;
170
171   // draw data points
172   TGraph *g = new TGraph(nPoints);
173   g->SetMarkerStyle(7);
174   for(int ip = 0; ip<nPoints; ip++) g->SetPoint(ip, fData[ax1][index[ip]], fData[ax2][index[ip]]);
175
176   // draw estimation point
177   TKDNodeInfo *node = (TKDNodeInfo*)(*fNodes)[inode];
178   TMarker *m=new TMarker(node->Data()[ax1], node->Data()[ax2], 20);
179   m->SetMarkerColor(2);
180   m->SetMarkerSize(1.7);
181   
182   // draw node contour
183   Float_t *bounds = GetBoundary(tnode);
184   TBox *n = new TBox(bounds[2*ax1], bounds[2*ax2], bounds[2*ax1+1], bounds[2*ax2+1]);
185   n->SetFillStyle(0);
186
187   g->Draw("ap");
188   m->Draw();
189   n->Draw();
190   
191   return;
192 }
193