Line data Source code
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 128 : ClassImp(TKDPDF)
14 :
15 :
16 :
17 : //_________________________________________________________________
18 : TKDPDF::TKDPDF() :
19 220 : TKDTreeIF()
20 220 : ,TKDInterpolatorBase()
21 1320 : {
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 440 : }
25 :
26 : //_________________________________________________________________
27 : TKDPDF::TKDPDF(Int_t npoints, Int_t ndim, UInt_t bsize, Float_t **data) :
28 0 : TKDTreeIF(npoints, ndim, bsize, data)
29 0 : ,TKDInterpolatorBase(ndim)
30 0 : {
31 : // Wrapper constructor for the TKDTree.
32 :
33 0 : Build();
34 0 : }
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 0 : TKDTreeIF()
40 0 : ,TKDInterpolatorBase()
41 0 : {
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 0 : TObjArray *vars = TString(var).Tokenize(":");
49 0 : fNDim = vars->GetEntriesFast(); fNDimm = 2*fNDim;
50 0 : fNSize = fNDim;
51 0 : 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 0 : fBucketSize = bsize;
53 :
54 : Int_t np;
55 : Double_t *v;
56 0 : for(int idim=0; idim<fNDim; idim++){
57 0 : if(!(np = t->Draw(((TObjString*)(*vars)[idim])->GetName(), cut, "goff", nentries, firstentry))){
58 0 : 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 0 : TIterator *it = (t->GetListOfLeaves())->MakeIterator();
60 : TObject *o;
61 0 : while((o = (*it)())) printf("\t%s\n", o->GetName());
62 : continue;
63 : }
64 0 : if(!fNPoints){
65 0 : 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 0 : fData = new Float_t*[fNDim];
68 0 : for(int jdim=fNDim; jdim--;) fData[jdim] = new Float_t[fNPoints];
69 0 : fDataOwner = kTRUE;
70 0 : }
71 0 : v = t->GetV1();
72 0 : for(int ip=0; ip<fNPoints; ip++) fData[idim][ip] = (Float_t)v[ip];
73 0 : }
74 0 : Build();
75 0 : delete vars;
76 0 : }
77 :
78 : //_________________________________________________________________
79 0 : TKDPDF::~TKDPDF()
80 0 : {
81 0 : }
82 :
83 : //_________________________________________________________________
84 0 : Bool_t TKDPDF::Build(Int_t)
85 : {
86 : // Fill interpolator's data array i.e.
87 : // - estimation points
88 : // - corresponding PDF values
89 :
90 0 : TKDTreeIF::Build();
91 0 : if(!fBoundaries) MakeBoundaries();
92 0 : fLambda = 1 + fNDim + (fNDim*(fNDim+1)>>1);
93 : //printf("after MakeBoundaries() %d\n", memory());
94 :
95 : // allocate interpolation nodes
96 0 : Int_t fNTNodes = fNPoints/fBucketSize + ((fNPoints%fBucketSize)?1:0);/*TKDTreeIF::GetNTNodes();*/
97 0 : TKDInterpolatorBase::Build(fNTNodes);
98 :
99 : TKDNodeInfo *node = NULL;
100 : Float_t *bounds = NULL;
101 : Int_t *indexPoints;
102 0 : for(int inode=0, tnode = fNNodes; inode<fNTNodes-1; inode++, tnode++){
103 0 : node = (TKDNodeInfo*)(*fNodes)[inode];
104 0 : node->Val()[0] = Float_t(fBucketSize)/fNPoints;
105 0 : bounds = GetBoundary(tnode);
106 0 : for(int idim=0; idim<fNDim; idim++) node->Val()[0] /= (bounds[2*idim+1] - bounds[2*idim]);
107 0 : node->Val()[1] = node->Val()[0]/TMath::Sqrt(float(fBucketSize));
108 :
109 0 : indexPoints = GetPointsIndexes(tnode);
110 : // loop points in this terminal node
111 0 : for(int idim=0; idim<fNDim; idim++){
112 0 : node->Data()[idim] = 0.;
113 0 : for(int ip = 0; ip<fBucketSize; ip++) node->Data()[idim] += fData[idim][indexPoints[ip]];
114 0 : node->Data()[idim] /= fBucketSize;
115 : }
116 0 : memcpy(&(node->Data()[fNDim]), bounds, fNDimm*sizeof(Float_t));
117 : }
118 :
119 : // analyze last (incomplete) terminal node
120 0 : Int_t counts = fNPoints%fBucketSize;
121 0 : counts = counts ? counts : fBucketSize;
122 0 : Int_t inode = fNTNodes - 1, tnode = inode + fNNodes;
123 0 : node = (TKDNodeInfo*)(*fNodes)[inode];
124 0 : node->Val()[0] = Float_t(counts)/fNPoints;
125 0 : bounds = GetBoundary(tnode);
126 0 : for(int idim=0; idim<fNDim; idim++){
127 0 : Float_t dx = bounds[2*idim+1]-bounds[2*idim];
128 0 : if(dx < 1.e-30){
129 0 : Warning("TKDPDF::Build()", "Terminal bucket index[%d] too narrow on the %d dimension.", inode, idim);
130 0 : continue;
131 : }
132 0 : node->Val()[0] /= (bounds[2*idim+1] - bounds[2*idim]);
133 0 : }
134 0 : node->Val()[1] = node->Val()[0]/TMath::Sqrt(float(counts));
135 :
136 : // loop points in this terminal node
137 0 : indexPoints = GetPointsIndexes(tnode);
138 0 : for(int idim=0; idim<fNDim; idim++){
139 0 : node->Data()[idim] = 0.;
140 0 : for(int ip = 0; ip<counts; ip++) node->Data()[idim] += fData[idim][indexPoints[ip]];
141 0 : node->Data()[idim] /= counts;
142 : }
143 0 : memcpy(&(node->Data()[fNDim]), bounds, fNDimm*sizeof(Float_t));
144 :
145 0 : delete [] fBoundaries;
146 0 : fBoundaries = NULL;
147 :
148 0 : 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 0 : if(tnode < 0 || tnode >= GetNTNodes()){
162 0 : Warning("DrawNode()", "Terminal node %d outside defined range.", tnode);
163 0 : return;
164 : }
165 :
166 : Int_t inode = tnode;
167 0 : tnode += fNNodes;
168 : // select zone of interest in the indexes array
169 0 : Int_t *index = GetPointsIndexes(tnode);
170 0 : Int_t nPoints = (tnode == 2*fNNodes) ? fNPoints%fBucketSize : fBucketSize;
171 :
172 : // draw data points
173 0 : TGraph *g = new TGraph(nPoints);
174 0 : g->SetMarkerStyle(7);
175 0 : for(int ip = 0; ip<nPoints; ip++) g->SetPoint(ip, fData[ax1][index[ip]], fData[ax2][index[ip]]);
176 :
177 : // draw estimation point
178 0 : TKDNodeInfo *node = (TKDNodeInfo*)(*fNodes)[inode];
179 0 : TMarker *m=new TMarker(node->Data()[ax1], node->Data()[ax2], 20);
180 0 : m->SetMarkerColor(2);
181 0 : m->SetMarkerSize(1.7);
182 :
183 : // draw node contour
184 0 : Float_t *bounds = GetBoundary(tnode);
185 0 : TBox *n = new TBox(bounds[2*ax1], bounds[2*ax2], bounds[2*ax1+1], bounds[2*ax2+1]);
186 0 : n->SetFillStyle(0);
187 :
188 0 : g->Draw("ap");
189 0 : m->Draw();
190 0 : n->Draw();
191 :
192 : return;
193 0 : }
194 :
|
