Coding conventions, removing compilation warnings, consistent usage of const
[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)", Form("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)", Form("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 idim=0; idim<fNDim; idim++) fData[idim] = new Float_t[fNpoints];
69                         kDataOwner = kTRUE;
70                 }
71                 v = t->GetV1();
72                 for(int ip=0; ip<fNpoints; ip++) fData[idim][ip] = (Float_t)v[ip];
73         }
74         TKDTreeIF::Build();
75         Build();
76 }
77
78 //_________________________________________________________________
79 TKDPDF::~TKDPDF()
80 {
81 }
82
83 //_________________________________________________________________
84 void TKDPDF::Build(Int_t)
85 {
86 // Fill interpolator's data array i.e.
87 //  - estimation points 
88 //  - corresponding PDF values
89
90         fNTNodes = fNpoints/fBucketSize + ((fNpoints%fBucketSize)?1:0);/*TKDTreeIF::GetNTNodes();*/
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         TKDInterpolatorBase::Build(fNTNodes);
97
98         TKDNodeInfo *node = 0x0;
99         Float_t *bounds = 0x0;
100         Int_t *indexPoints;
101         for(int inode=0, tnode = fNnodes; inode<fNTNodes-1; inode++, tnode++){
102                 node = (TKDNodeInfo*)(*fTNodes)[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*)(*fTNodes)[inode];
123         node->Val()[0] =  Float_t(counts)/fNpoints;
124         bounds = GetBoundary(tnode);
125         for(int idim=0; idim<fNDim; idim++) node->Val()[0] /= (bounds[2*idim+1] - bounds[2*idim]);
126         node->Val()[1] =  node->Val()[0]/TMath::Sqrt(float(counts));
127
128         // loop points in this terminal node
129         indexPoints = GetPointsIndexes(tnode);
130         for(int idim=0; idim<fNDim; idim++){
131                 node->Data()[idim] = 0.;
132                 for(int ip = 0; ip<counts; ip++) node->Data()[idim] += fData[idim][indexPoints[ip]];
133                 node->Data()[idim] /= counts;
134         }
135         memcpy(&(node->Data()[fNDim]), bounds, fNDimm*sizeof(Float_t));
136
137         delete [] fBoundaries;
138         fBoundaries = 0x0;
139 }
140
141
142 //_________________________________________________________________
143 void TKDPDF::DrawNode(Int_t tnode, UInt_t ax1, UInt_t ax2)
144 {
145 // Draw node "node" and the data points within.
146 //
147 // Observation:
148 // This function creates some graphical objects
149 // but don't delete it. Abusing this function may cause memory leaks !
150
151         if(tnode < 0 || tnode >= fNTNodes){
152                 Warning("DrawNode()", Form("Terminal node %d outside defined range.", tnode));
153                 return;
154         }
155
156         Int_t inode = tnode;
157         tnode += fNnodes;
158         // select zone of interest in the indexes array
159         Int_t *index = GetPointsIndexes(tnode);
160         Int_t nPoints = (tnode == 2*fNnodes) ? fNpoints%fBucketSize : fBucketSize;
161
162         // draw data points
163         TGraph *g = new TGraph(nPoints);
164         g->SetMarkerStyle(7);
165         for(int ip = 0; ip<nPoints; ip++) g->SetPoint(ip, fData[ax1][index[ip]], fData[ax2][index[ip]]);
166
167         // draw estimation point
168         TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[inode];
169         TMarker *m=new TMarker(node->Data()[ax1], node->Data()[ax2], 20);
170         m->SetMarkerColor(2);
171         m->SetMarkerSize(1.7);
172         
173         // draw node contour
174         Float_t *bounds = GetBoundary(tnode);
175         TBox *n = new TBox(bounds[2*ax1], bounds[2*ax2], bounds[2*ax1+1], bounds[2*ax2+1]);
176         n->SetFillStyle(0);
177
178         g->Draw("ap");
179         m->Draw();
180         n->Draw();
181         
182         return;
183 }
184