1 #include "AliTRDTKDInterpolator.h"
3 #include "TClonesArray.h"
4 #include "TLinearFitter.h"
11 ClassImp(AliTRDTKDInterpolator)
13 //_________________________________________________________________
14 AliTRDTKDInterpolator::AliTRDTKDInterpolator() :
20 fUseHelperNodes(kFALSE),
22 fPDFMode(kInterpolation)
24 // default constructor
27 //_________________________________________________________________
28 AliTRDTKDInterpolator::AliTRDTKDInterpolator(Int_t npoints, Int_t ndim, UInt_t bsize, Float_t **data) :
29 TKDTreeIF(npoints, ndim, bsize, data),
32 fLambda(1 + ndim + (ndim*(ndim+1)>>1)),
34 fUseHelperNodes(kFALSE),
36 fPDFMode(kInterpolation)
40 //_________________________________________________________________
41 AliTRDTKDInterpolator::~AliTRDTKDInterpolator()
49 //_________________________________________________________________
50 AliTRDTKDInterpolator::AliTRDTKDInterpolator(const AliTRDTKDInterpolator &ref):
52 fNDataNodes(ref.fNDataNodes),
55 fNPointsI(ref.fNPointsI),
56 fUseHelperNodes(ref.fUseHelperNodes),
57 fUseWeights(ref.fUseWeights),
58 fPDFMode(ref.fPDFMode)
63 //____________________________________________________________
65 AliTRDTKDInterpolator &AliTRDTKDInterpolator::operator=(const AliTRDTKDInterpolator &ref){
67 // Assignment operator
69 if(this == &ref) return *this;
72 TObject::operator=(ref);
77 //_________________________________________________________________
78 Bool_t AliTRDTKDInterpolator::Build()
81 if(!fBoundaries) MakeBoundaries();
83 // allocate interpolation nodes
84 fNDataNodes = fNPoints/fBucketSize + ((fNPoints%fBucketSize)?1:0);/*TKDTreeIF::GetNTNodes();*/
87 Warning("AliTRDTKDInterpolator::Build()", "Data already allocated.");
90 fNodes = new TClonesArray("AliTRDTKDInterpolator::AliTRDTKDNodeInfo", fNDataNodes);
93 for(int in=0; in<fNDataNodes; in++) new ((*fNodes)[in]) AliTRDTKDNodeInfo(fNDim);
95 // Set Interpolator nodes
97 for(int inode=0, tnode = fNNodes; inode<fNDataNodes-1; inode++, tnode++){
98 AliTRDTKDNodeInfo *node =GetNodeInfo(inode);
99 memcpy(node->fBounds,GetBoundary(tnode),2*fNDim*sizeof(Float_t));
100 node->fVal[0] = Float_t(fBucketSize)/fNPoints;
101 for(int idim=0; idim<fNDim; idim++) node->fVal[0] /= (node->fBounds[2*idim+1] - node->fBounds[2*idim]);
102 node->fVal[1] = node->fVal[0]/TMath::Sqrt(float(fBucketSize));
104 Int_t *indexPoints = GetPointsIndexes(tnode);
105 // loop points in this terminal node
106 for(int idim=0; idim<fNDim; idim++){
107 node->fData[idim] = 0.;
108 for(int ip = 0; ip<fBucketSize; ip++) node->fData[idim] += fData[idim][indexPoints[ip]];
109 node->fData[idim] /= fBucketSize;
113 // Analyze last (incomplete) terminal node
115 Int_t counts = fNPoints%fBucketSize;
116 counts = counts ? counts : fBucketSize;
117 Int_t inode = fNDataNodes - 1, tnode = inode + fNNodes;
118 AliTRDTKDNodeInfo *ftnode = GetNodeInfo(inode);
119 ftnode->fVal[0] = Float_t(counts)/fNPoints;
120 memcpy(ftnode->fBounds,GetBoundary(tnode),2*fNDim*sizeof(Float_t));
122 for(int idim=0; idim<fNDim; idim++){
123 Float_t dx = ftnode->fBounds[2*idim+1]-ftnode->fBounds[2*idim];
125 Warning("AliTRDTKDInterpolator::Build()", "Terminal bucket index[%d] too narrow on the %d dimension.", inode, idim);
128 ftnode->fVal[0] /= (ftnode->fBounds[2*idim+1] - ftnode->fBounds[2*idim]);
130 ftnode->fVal[1] = ftnode->fVal[0]/TMath::Sqrt(float(counts));
132 // loop points in this terminal node
133 Int_t *indexPoints = GetPointsIndexes(tnode);
134 for(int idim=0; idim<fNDim; idim++){
135 ftnode->fData[idim] = 0.;
136 for(int ip = 0; ip<counts; ip++) ftnode->fData[idim] += fData[idim][indexPoints[ip]];
137 ftnode->fData[idim] /= counts;
140 delete [] fBoundaries;
143 if(fUseHelperNodes){BuildBoundaryNodes();}
145 if(fNPointsI>GetNTNodes()){fNPointsI=GetNTNodes();}
147 BuildInterpolation();
153 //_________________________________________________________________
154 Bool_t AliTRDTKDInterpolator::Eval(const Double_t *point, Double_t &result, Double_t &error)
156 Float_t pointF[fNDim]; // local Float_t conversion for "point"
157 for(int idim=0; idim<fNDim; idim++) pointF[idim] = (Float_t)point[idim];
158 Int_t nodeIndex = GetNodeIndex(pointF);
160 Error("AliTRDTKDInterpolator::Eval()", "Can not retrive node for data point.");
165 AliTRDTKDNodeInfo *node =GetNodeInfo(nodeIndex);
170 return node->CookPDF(point, result, error);
172 node->CookPDF(point, result, error);
173 if(error<node->fVal[1]){
177 result=node->fVal[0];
181 result=node->fVal[0];
188 //__________________________________________________________________
189 void AliTRDTKDInterpolator::Print(const Option_t */*opt*/) const
191 for(Int_t i=GetNTNodes(); i--;){
192 printf("Node %d of %d: \n",i,GetNTNodes());
193 GetNodeInfo(i)->Print();
198 //__________________________________________________________________
199 Int_t AliTRDTKDInterpolator::GetNodeIndex(const Float_t *p)
201 Int_t inode=FindNode(p)-fNDataNodes+1;
202 if(GetNodeInfo(inode)->Has(p))return inode;
204 // Search extra nodes
206 for(inode=fNDataNodes;inode<GetNTNodes();inode++){
207 if(GetNodeInfo(inode)->Has(p)){return inode;}
210 // Search for nearest neighbor
212 Float_t closestdist=10000;
214 for(Int_t ii=0;ii<GetNTNodes();ii++){
215 AliTRDTKDNodeInfo *node=GetNodeInfo(ii);
217 for(Int_t idim=0;idim<fNDim;idim++){
218 dist+=TMath::Power((node->fData[idim]-p[idim]),2);
220 dist=TMath::Sqrt(dist);
221 if(dist<closestdist){closestdist=dist;inode=ii;}
226 //_________________________________________________________________
227 AliTRDTKDInterpolator::AliTRDTKDNodeInfo* AliTRDTKDInterpolator::GetNodeInfo(Int_t inode) const
229 if(!fNodes || inode >= GetNTNodes()) return NULL;
230 return (AliTRDTKDNodeInfo*)(*fNodes)[inode];
233 //_________________________________________________________________
234 Int_t AliTRDTKDInterpolator::GetNTNodes() const
236 return fNodes?fNodes->GetEntriesFast():0;
239 //_________________________________________________________________
240 Bool_t AliTRDTKDInterpolator::GetRange(Int_t idim,Float_t range[2]) const
242 if(!fNodes) return kFALSE;
243 if(idim<0 || idim>=fNDim){
244 range[0]=0.; range[1]=0.;
247 range[0]=1.e10; range[1]=-1.e10;
248 for(Int_t in=GetNTNodes(); in--; ){
249 AliTRDTKDNodeInfo *node = GetNodeInfo(in);
251 if(node->fBounds[2*idim]<range[0]) range[0] = node->fBounds[2*idim];
252 if(node->fBounds[2*idim+1]>range[1]) range[1] = node->fBounds[2*idim+1];
258 //_________________________________________________________________
259 TH2Poly *AliTRDTKDInterpolator::Projection(Int_t xdim,Int_t ydim)
261 Float_t rangex[2],rangey[2];
262 GetRange(xdim,rangex);
263 GetRange(ydim,rangey);
265 TH2Poly* h2 = new TH2Poly("hTKDnodes","hTKDnodes", rangex[0],rangex[1],rangey[0],rangey[1]);
266 h2->GetXaxis()->SetTitle(Form("Q_{%d}", xdim));
267 h2->GetYaxis()->SetTitle(Form("Q_{%d}", ydim));
269 for(Int_t inode=0;inode<GetNTNodes();inode++){
271 AliTRDTKDNodeInfo* node=GetNodeInfo(inode);
272 h2->AddBin(node->fBounds[2*xdim],node->fBounds[2*ydim],node->fBounds[2*xdim+1],node->fBounds[2*ydim+1]);
273 h2->SetBinContent(inode+1, node->fVal[0]);
278 //_________________________________________________________________
279 void AliTRDTKDInterpolator::BuildInterpolation()
282 // Calculate Interpolation
284 Double_t buffer[fLambda];
286 Float_t **refPoints = new Float_t*[fNDim];
287 for(int id=0; id<fNDim; id++){
288 refPoints[id] = new Float_t[GetNTNodes()];
289 for(int in=0; in<GetNTNodes(); in++) refPoints[id][in] = GetNodeInfo(in)->fData[id];
291 TKDTreeIF *KDhelper = new TKDTreeIF(GetNTNodes(), fNDim, 30, refPoints);
293 KDhelper->MakeBoundariesExact();
295 Float_t dist[fNPointsI];
296 Int_t ind[fNPointsI];
298 TLinearFitter fitter(fLambda, Form("hyp%d", fLambda-1));
300 Int_t nodeIndex(0); Float_t param[6], *pp(NULL);
301 nodeIndex=GetNTNodes(); pp=¶m[0];
304 fitter.ClearPoints();
306 AliTRDTKDNodeInfo *node = GetNodeInfo(nodeIndex);
307 // find nearest neighbors
308 KDhelper->FindNearestNeighbors(node->fData,fNPointsI, &ind[0], &dist[0]);
310 for(int in=0;in<fNPointsI;in++){
311 AliTRDTKDNodeInfo *nnode = GetNodeInfo(ind[in]);
313 Float_t w=1; //weight
314 // calculate tri-cubic weighting function
316 Float_t d = dist[in]/dist[fNPointsI-1];
317 Float_t w0 = (1. - d*d*d);
319 if(w<1.e-30) continue;
322 for(int idim=0; idim<fNDim; idim++){
323 buffer[ipar++] = nnode->fData[idim];
324 for(int jdim=idim; jdim<fNDim; jdim++) buffer[ipar++] = nnode->fData[idim]*nnode->fData[jdim];
326 fitter.AddPoint(buffer,nnode->fVal[0], nnode->fVal[1]/w);
328 // Ensure Boundary Condition
329 for(Int_t kdim=0;kdim<fNDim;kdim++){
330 if(node->fBounds[2*kdim]==0){
331 Float_t zdata[fNDim];
332 memcpy(&zdata[0],node->fData,fNDim*sizeof(Float_t));
335 for(int idim=0; idim<fNDim; idim++){
336 buffer[ipar++] = zdata[idim];
337 for(int jdim=idim; jdim<fNDim; jdim++) buffer[ipar++] = zdata[idim]*zdata[jdim];
339 fitter.AddPoint(buffer,0,1);
346 // retrive fitter results
347 TMatrixD cov(fLambda, fLambda);
348 TVectorD par(fLambda);
349 fitter.GetCovarianceMatrix(cov);
350 fitter.GetParameters(par);
353 node->Store(&par,&cov);
357 for(int id=0; id<fNDim; id++){
358 delete refPoints[id];
363 //_________________________________________________________________
364 void AliTRDTKDInterpolator::BuildBoundaryNodes(){
368 Float_t treebounds[2*fNDim];
369 for(Int_t idim=0;idim<fNDim;idim++){
370 GetRange(idim,&treebounds[2*idim]);
373 for(int inode=0; inode<GetNTNodes(); inode++){
375 AliTRDTKDNodeInfo *node=GetNodeInfo(inode);
377 for(Int_t vdim=0;vdim<fNDim;vdim++){
379 // Try expansion to lower and higher values
380 for(Int_t iter=0;iter<2;iter++){
381 if(node->fBounds[2*vdim+iter]==treebounds[2*vdim+iter]){
384 new ((*fNodes)[GetNTNodes()]) AliTRDTKDNodeInfo(fNDim);
386 AliTRDTKDNodeInfo *newnode = GetNodeInfo(GetNTNodes()-1);
387 if(iter==0)newnode->fBounds[2*vdim+iter]=0;
388 if(iter==1)newnode->fBounds[2*vdim+iter]=2*treebounds[2*vdim+iter];
389 newnode->fBounds[2*vdim+!iter]=node->fBounds[2*vdim+iter];
390 for(Int_t idim=0;idim<fNDim;idim++){
391 if(idim==vdim)continue;
392 newnode->fBounds[2*idim]=node->fBounds[2*idim];
393 newnode->fBounds[2*idim+1]=node->fBounds[2*idim+1];
396 newnode->fVal[1]=Float_t(1)/fNPoints;
397 for(int idim=0; idim<fNDim; idim++){
398 newnode->fVal[1] /= (newnode->fBounds[2*idim+1] - newnode->fBounds[2*idim]);
399 newnode->fData[idim]=0.5*(newnode->fBounds[2*idim+1] + newnode->fBounds[2*idim]);
406 printf("%d Boundary Nodes Added \n",nnew);
409 //_________________________________________________________________
410 AliTRDTKDInterpolator::AliTRDTKDNodeInfo::AliTRDTKDNodeInfo(Int_t ndim):
414 ,fNPar(1 + ndim + (ndim*(ndim+1)>>1))
421 // Default constructor
422 fVal[0] = 0.; fVal[1] = 0.;
423 fData=new Float_t[fNDim];
424 fBounds=new Float_t[fNBounds];
427 //_________________________________________________________________
428 AliTRDTKDInterpolator::AliTRDTKDNodeInfo::AliTRDTKDNodeInfo(const AliTRDTKDInterpolator::AliTRDTKDNodeInfo &ref):
431 ,fNBounds(ref.fNBounds)
442 fData = new Float_t[fNDim];
443 memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
446 fBounds = new Float_t[2*fNDim];
447 memcpy(fBounds, ref.fBounds, 2*fNDim*sizeof(Float_t));
450 fVal[0] = ref.fVal[0];
451 fVal[1] = ref.fVal[1];
454 fPar=new Double_t[fNPar];
455 memcpy(fPar, ref.fPar, fNPar*sizeof(Double_t));
458 fCov=new Double_t[fNCov];
459 memcpy(fCov, ref.fCov, fNCov*sizeof(Double_t));
463 //_________________________________________________________________
464 AliTRDTKDInterpolator::AliTRDTKDNodeInfo::~AliTRDTKDNodeInfo()
467 if(fData) delete [] fData;
468 if(fBounds) delete [] fBounds;
469 if(fPar) delete [] fPar;
470 if(fCov) delete [] fCov;
473 //_________________________________________________________________
474 AliTRDTKDInterpolator::AliTRDTKDNodeInfo &AliTRDTKDInterpolator::AliTRDTKDNodeInfo::operator=(const AliTRDTKDInterpolator::AliTRDTKDNodeInfo &ref)
477 // Assignment operator
480 if(&ref==this) return *this;
481 memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
482 memcpy(fBounds, ref.fBounds, 2*fNDim*sizeof(Float_t));
484 fVal[0] = ref.fVal[0];
485 fVal[1] = ref.fVal[1];
488 fPar=new Double_t[fNPar];
489 memcpy(fPar, ref.fPar, fNPar*sizeof(Double_t));
492 fCov=new Double_t[fNCov];
493 memcpy(fCov, ref.fCov, fNCov*sizeof(Double_t));
499 //_________________________________________________________________
500 void AliTRDTKDInterpolator::AliTRDTKDNodeInfo::Print(const Option_t *opt) const
502 // Print the content of the node
504 for(int idim=0; idim<fNDim; idim++) printf("%f ", fData?fData[idim]:0.);
505 printf("] f = [%e +- %e]\n", fVal[0], fVal[1]);
509 for(int idim=0; idim<fNDim; idim++) printf("(%f %f) ", fBounds[2*idim], fBounds[2*idim+1]);
512 if(strcmp(opt, "a")!=0) return;
515 printf("Fit parameters : \n");
516 for(int ip=0; ip<fNPar; ip++) printf("p%d[%f] ", ip, fPar[ip]);
520 for(int ip(0), n(0); ip<fNPar; ip++){
521 for(int jp(ip); jp<fNPar; jp++) printf("c(%d %d)[%f] ", ip, jp, fCov[n++]);
526 //_________________________________________________________________
527 void AliTRDTKDInterpolator::AliTRDTKDNodeInfo::Store(TVectorD const *par, TMatrixD const *cov)
529 // Store the parameters and the covariance in the node
531 if(!fPar){fPar = new Double_t[fNPar];}
532 for(int ip=0; ip<fNPar; ip++) fPar[ip] = (*par)[ip];
535 if(!fCov){fCov = new Double_t[fNCov];}
536 for(int ip(0), np(0); ip<fNPar; ip++)
537 for(int jp=ip; jp<fNPar; jp++) fCov[np++] = (*cov)(ip,jp);
540 //_________________________________________________________________
541 Bool_t AliTRDTKDInterpolator::AliTRDTKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t &error) const
543 // Recalculate the PDF for one node from the results of interpolation (parameters and covariance matrix)
545 result =0.; error = 1.;
546 if(!fPar) return kFALSE;
548 Double_t fdfdp[fNDim];
551 for(int idim=0; idim<fNDim; idim++){
552 fdfdp[ipar++] = point[idim];
553 for(int jdim=idim; jdim<fNDim; jdim++) fdfdp[ipar++] = point[idim]*point[jdim];
556 // calculate estimation
557 for(int i=0; i<fNPar; i++) result += fdfdp[i]*fPar[i];
559 if(!fCov)return kTRUE;
563 for(int i(0), n(0); i<fNPar; i++){
564 error += fdfdp[i]*fdfdp[i]*fCov[n++];
565 for(int j(i+1); j<fNPar; j++) error += 2.*fdfdp[i]*fdfdp[j]*fCov[n++];
567 if(error>0)error = TMath::Sqrt(error);
570 // Boundary condition
579 //_____________________________________________________________________
580 Bool_t AliTRDTKDInterpolator::AliTRDTKDNodeInfo::Has(const Float_t *p) const
583 for(int id=0; id<fNDim; id++)
584 if(p[id]>=fBounds[2*id] && p[id]<fBounds[2*id+1]) n++;
585 if(n==fNDim) return kTRUE;