1 #include "AliTRDTKDInterpolator.h"
3 #include "TClonesArray.h"
4 #include "TLinearFitter.h"
13 ClassImp(AliTRDTKDInterpolator)
15 //_________________________________________________________________
16 AliTRDTKDInterpolator::AliTRDTKDInterpolator() :
22 fUseHelperNodes(kFALSE),
24 fPDFMode(kInterpolation),
27 // default constructor
30 //_________________________________________________________________
31 AliTRDTKDInterpolator::AliTRDTKDInterpolator(Int_t npoints, Int_t ndim, UInt_t bsize, Float_t **data) :
32 TKDTreeIF(npoints, ndim, bsize, data),
35 fLambda(1 + ndim + (ndim*(ndim+1)>>1)),
37 fUseHelperNodes(kFALSE),
39 fPDFMode(kInterpolation),
44 //_________________________________________________________________
45 AliTRDTKDInterpolator::~AliTRDTKDInterpolator()
53 //_________________________________________________________________
54 AliTRDTKDInterpolator::AliTRDTKDInterpolator(const AliTRDTKDInterpolator &ref):
56 fNDataNodes(ref.fNDataNodes),
59 fNPointsI(ref.fNPointsI),
60 fUseHelperNodes(ref.fUseHelperNodes),
61 fUseWeights(ref.fUseWeights),
62 fPDFMode(ref.fPDFMode),
63 fStoreCov(ref.fStoreCov)
69 //____________________________________________________________
71 AliTRDTKDInterpolator &AliTRDTKDInterpolator::operator=(const AliTRDTKDInterpolator &ref){
73 // Assignment operator
75 if(this == &ref) return *this;
78 TObject::operator=(ref);
85 //_________________________________________________________________
86 Bool_t AliTRDTKDInterpolator::Build()
89 if(!fBoundaries) MakeBoundaries();
91 // allocate interpolation nodes
92 fNDataNodes = fNPoints/fBucketSize + ((fNPoints%fBucketSize)?1:0);
95 Warning("AliTRDTKDInterpolator::Build()", "Data already allocated.");
98 fNodes = new TClonesArray("AliTRDTKDInterpolator::AliTRDTKDNodeInfo", fNDataNodes);
101 for(int in=0; in<fNDataNodes; in++) new ((*fNodes)[in]) AliTRDTKDNodeInfo(fNDim);
103 // Set Interpolator nodes
105 for(int inode=0, tnode = fNNodes; inode<fNDataNodes-1; inode++, tnode++){
106 AliTRDTKDNodeInfo *node =GetNodeInfo(inode);
107 memcpy(node->fBounds,GetBoundary(tnode),2*fNDim*sizeof(Float_t));
108 node->fVal[0] = Float_t(fBucketSize)/fNPoints;
109 for(int idim=0; idim<fNDim; idim++) node->fVal[0] /= (node->fBounds[2*idim+1] - node->fBounds[2*idim]);
110 node->fVal[1] = node->fVal[0]/TMath::Sqrt(float(fBucketSize));
112 Int_t *indexPoints = GetPointsIndexes(tnode);
113 // loop points in this terminal node
114 for(int idim=0; idim<fNDim; idim++){
115 node->fData[idim] = 0.;
116 for(int ip = 0; ip<fBucketSize; ip++) node->fData[idim] += fData[idim][indexPoints[ip]];
117 node->fData[idim] /= fBucketSize;
121 // Analyze last (incomplete) terminal node
123 Int_t counts = fNPoints%fBucketSize;
124 counts = counts ? counts : fBucketSize;
125 Int_t inode = fNDataNodes - 1, tnode = inode + fNNodes;
126 AliTRDTKDNodeInfo *ftnode = GetNodeInfo(inode);
127 ftnode->fVal[0] = Float_t(counts)/fNPoints;
128 memcpy(ftnode->fBounds,GetBoundary(tnode),2*fNDim*sizeof(Float_t));
130 for(int idim=0; idim<fNDim; idim++){
131 Float_t dx = ftnode->fBounds[2*idim+1]-ftnode->fBounds[2*idim];
133 Warning("AliTRDTKDInterpolator::Build()", "Terminal bucket index[%d] too narrow on the %d dimension.", inode, idim);
136 ftnode->fVal[0] /= (ftnode->fBounds[2*idim+1] - ftnode->fBounds[2*idim]);
138 ftnode->fVal[1] = ftnode->fVal[0]/TMath::Sqrt(float(counts));
140 // loop points in this terminal node
141 Int_t *indexPoints = GetPointsIndexes(tnode);
142 for(int idim=0; idim<fNDim; idim++){
143 ftnode->fData[idim] = 0.;
144 for(int ip = 0; ip<counts; ip++) ftnode->fData[idim] += fData[idim][indexPoints[ip]];
145 ftnode->fData[idim] /= counts;
148 delete [] fBoundaries;
151 if(fUseHelperNodes){BuildBoundaryNodes();}
153 if(fNPointsI>GetNTNodes()){fNPointsI=GetNTNodes();}
155 BuildInterpolation();
161 //_________________________________________________________________
162 Bool_t AliTRDTKDInterpolator::Eval(const Double_t *point, Double_t &result, Double_t &error)
164 AliDebug(3,Form("Eval PDF Mode %d",fPDFMode));
165 if((AliLog::GetDebugLevel("",IsA()->GetName()))>0){
167 for(int idim=0; idim<fNDim; idim++) printf("%f ", point[idim]);
171 Float_t pointF[fNDim]; // local Float_t conversion for "point"
172 for(int idim=0; idim<fNDim; idim++) pointF[idim] = (Float_t)point[idim];
173 Int_t nodeIndex = GetNodeIndex(pointF);
175 AliError("Can not retrieve node for data point");
180 AliTRDTKDNodeInfo *node =GetNodeInfo(nodeIndex);
182 if((AliLog::GetDebugLevel("",IsA()->GetName()))>0){
183 printf("Node Info: \n");
187 return node->CookPDF(point, result, error,fPDFMode);
190 //__________________________________________________________________
191 void AliTRDTKDInterpolator::Print(const Option_t */*opt*/) const
193 for(Int_t i=GetNTNodes(); i--;){
194 printf("Node %d of %d: \n",i,GetNTNodes());
195 GetNodeInfo(i)->Print();
200 //__________________________________________________________________
201 Int_t AliTRDTKDInterpolator::GetNodeIndex(const Float_t *p)
203 Int_t inode=FindNode(p)-fNDataNodes+1;
204 if(GetNodeInfo(inode)->Has(p)){
205 AliDebug(2,Form("Find Node %d",inode));
209 // Search extra nodes
211 for(inode=fNDataNodes;inode<GetNTNodes();inode++){
212 if(GetNodeInfo(inode)->Has(p)){
213 AliDebug(2,Form("Find Extra Node %d",inode));
218 // Search for nearest neighbor
220 Float_t closestdist=10000;
222 for(Int_t ii=0;ii<GetNTNodes();ii++){
223 AliTRDTKDNodeInfo *node=GetNodeInfo(ii);
225 for(Int_t idim=0;idim<fNDim;idim++){
226 dist+=TMath::Power((node->fData[idim]-p[idim]),2);
228 dist=TMath::Sqrt(dist);
229 if(dist<closestdist){closestdist=dist;inode=ii;}
231 AliDebug(2,Form("Find Nearest Neighbor Node %d",inode));
235 //_________________________________________________________________
236 AliTRDTKDInterpolator::AliTRDTKDNodeInfo* AliTRDTKDInterpolator::GetNodeInfo(Int_t inode) const
238 if(!fNodes || inode >= GetNTNodes()) return NULL;
239 return (AliTRDTKDNodeInfo*)(*fNodes)[inode];
242 //_________________________________________________________________
243 Int_t AliTRDTKDInterpolator::GetNTNodes() const
245 return fNodes?fNodes->GetEntriesFast():0;
248 //_________________________________________________________________
249 Bool_t AliTRDTKDInterpolator::GetRange(Int_t idim,Float_t range[2]) const
251 if(!fNodes) return kFALSE;
252 if(idim<0 || idim>=fNDim){
253 range[0]=0.; range[1]=0.;
256 range[0]=1.e10; range[1]=-1.e10;
257 for(Int_t in=GetNTNodes(); in--; ){
258 AliTRDTKDNodeInfo *node = GetNodeInfo(in);
260 if(node->fBounds[2*idim]<range[0]) range[0] = node->fBounds[2*idim];
261 if(node->fBounds[2*idim+1]>range[1]) range[1] = node->fBounds[2*idim+1];
267 //_________________________________________________________________
268 TH2Poly *AliTRDTKDInterpolator::Projection(Int_t xdim,Int_t ydim)
270 Float_t rangex[2],rangey[2];
271 GetRange(xdim,rangex);
272 GetRange(ydim,rangey);
274 TH2Poly* h2 = new TH2Poly("hTKDnodes","hTKDnodes", rangex[0],rangex[1],rangey[0],rangey[1]);
275 h2->GetXaxis()->SetTitle(Form("Q_{%d}", xdim));
276 h2->GetYaxis()->SetTitle(Form("Q_{%d}", ydim));
278 for(Int_t inode=0;inode<GetNTNodes();inode++){
280 AliTRDTKDNodeInfo* node=GetNodeInfo(inode);
281 h2->AddBin(node->fBounds[2*xdim],node->fBounds[2*ydim],node->fBounds[2*xdim+1],node->fBounds[2*ydim+1]);
282 h2->SetBinContent(inode+1, node->fVal[0]);
287 //_________________________________________________________________
288 void AliTRDTKDInterpolator::BuildInterpolation()
290 AliInfo("Build Interpolation");
292 // Calculate Interpolation
294 Double_t buffer[fLambda];
296 Float_t **refPoints = new Float_t*[fNDim];
297 for(int id=0; id<fNDim; id++){
298 refPoints[id] = new Float_t[GetNTNodes()];
299 for(int in=0; in<GetNTNodes(); in++) refPoints[id][in] = GetNodeInfo(in)->fData[id];
301 TKDTreeIF *KDhelper = new TKDTreeIF(GetNTNodes(), fNDim, 30, refPoints);
303 KDhelper->MakeBoundariesExact();
305 Float_t dist[fNPointsI];
306 Int_t ind[fNPointsI];
308 TLinearFitter fitter(fLambda, Form("hyp%d", fLambda-1));
310 Int_t nodeIndex(0); Float_t param[6], *pp(NULL);
311 nodeIndex=GetNTNodes(); pp=¶m[0];
314 fitter.ClearPoints();
316 AliTRDTKDNodeInfo *node = GetNodeInfo(nodeIndex);
317 // find nearest neighbors
318 KDhelper->FindNearestNeighbors(node->fData,fNPointsI, &ind[0], &dist[0]);
320 for(int in=0;in<fNPointsI;in++){
321 AliTRDTKDNodeInfo *nnode = GetNodeInfo(ind[in]);
323 Float_t w=1; //weight
324 // calculate tri-cubic weighting function
326 Float_t d = dist[in]/dist[fNPointsI-1];
327 Float_t w0 = (1. - d*d*d);
329 if(w<1.e-30) continue;
332 for(int idim=0; idim<fNDim; idim++){
333 buffer[ipar++] = nnode->fData[idim];
334 for(int jdim=idim; jdim<fNDim; jdim++) buffer[ipar++] = nnode->fData[idim]*nnode->fData[jdim];
336 fitter.AddPoint(buffer,nnode->fVal[0], nnode->fVal[1]/w);
338 // Ensure Boundary Condition
339 for(Int_t kdim=0;kdim<fNDim;kdim++){
340 if(node->fBounds[2*kdim]==0){
341 Float_t zdata[fNDim];
342 memcpy(&zdata[0],node->fData,fNDim*sizeof(Float_t));
345 for(int idim=0; idim<fNDim; idim++){
346 buffer[ipar++] = zdata[idim];
347 for(int jdim=idim; jdim<fNDim; jdim++) buffer[ipar++] = zdata[idim]*zdata[jdim];
349 fitter.AddPoint(buffer,0,1);
354 AliDebug(2,Form("Calculate Interpolation for Node %d",nodeIndex));
357 // retrive fitter results
358 TMatrixD cov(fLambda, fLambda);
359 TVectorD par(fLambda);
360 fitter.GetCovarianceMatrix(cov);
361 fitter.GetParameters(par);
364 node->Store(&par,&cov,fStoreCov);
368 for(int id=0; id<fNDim; id++){
369 delete refPoints[id];
374 //_________________________________________________________________
375 void AliTRDTKDInterpolator::BuildBoundaryNodes(){
379 Float_t treebounds[2*fNDim];
380 for(Int_t idim=0;idim<fNDim;idim++){
381 GetRange(idim,&treebounds[2*idim]);
384 for(int inode=0; inode<GetNTNodes(); inode++){
386 AliTRDTKDNodeInfo *node=GetNodeInfo(inode);
388 for(Int_t vdim=0;vdim<fNDim;vdim++){
390 // Try expansion to lower and higher values
391 for(Int_t iter=0;iter<2;iter++){
392 if(node->fBounds[2*vdim+iter]==treebounds[2*vdim+iter]){
395 new ((*fNodes)[GetNTNodes()]) AliTRDTKDNodeInfo(fNDim);
397 AliTRDTKDNodeInfo *newnode = GetNodeInfo(GetNTNodes()-1);
398 if(iter==0)newnode->fBounds[2*vdim+iter]=0;
399 if(iter==1)newnode->fBounds[2*vdim+iter]=2*treebounds[2*vdim+iter];
400 newnode->fBounds[2*vdim+!iter]=node->fBounds[2*vdim+iter];
401 for(Int_t idim=0;idim<fNDim;idim++){
402 if(idim==vdim)continue;
403 newnode->fBounds[2*idim]=node->fBounds[2*idim];
404 newnode->fBounds[2*idim+1]=node->fBounds[2*idim+1];
407 newnode->fVal[1]=Float_t(1)/fNPoints;
408 for(int idim=0; idim<fNDim; idim++){
409 newnode->fVal[1] /= (newnode->fBounds[2*idim+1] - newnode->fBounds[2*idim]);
410 newnode->fData[idim]=0.5*(newnode->fBounds[2*idim+1] + newnode->fBounds[2*idim]);
417 AliInfo(Form("%d Boundary Nodes Added \n",nnew));
420 //_________________________________________________________________
421 AliTRDTKDInterpolator::AliTRDTKDNodeInfo::AliTRDTKDNodeInfo(Int_t ndim):
425 ,fNPar(1 + ndim + (ndim*(ndim+1)>>1))
426 ,fNCov(Int_t((fNPar+1)*fNPar/2))
432 // Default constructor
433 fVal[0] = 0.; fVal[1] = 0.;
434 fData=new Float_t[fNDim];
435 fBounds=new Float_t[fNBounds];
438 //_________________________________________________________________
439 AliTRDTKDInterpolator::AliTRDTKDNodeInfo::AliTRDTKDNodeInfo(const AliTRDTKDInterpolator::AliTRDTKDNodeInfo &ref):
442 ,fNBounds(ref.fNBounds)
453 fData = new Float_t[fNDim];
454 memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
457 fBounds = new Float_t[2*fNDim];
458 memcpy(fBounds, ref.fBounds, 2*fNDim*sizeof(Float_t));
461 fVal[0] = ref.fVal[0];
462 fVal[1] = ref.fVal[1];
465 fPar=new Double_t[fNPar];
466 memcpy(fPar, ref.fPar, fNPar*sizeof(Double_t));
469 fCov=new Double_t[fNCov];
470 memcpy(fCov, ref.fCov, fNCov*sizeof(Double_t));
474 //_________________________________________________________________
475 AliTRDTKDInterpolator::AliTRDTKDNodeInfo::~AliTRDTKDNodeInfo()
478 if(fData) delete [] fData;
479 if(fBounds) delete [] fBounds;
480 if(fPar) delete [] fPar;
481 if(fCov) delete [] fCov;
484 //_________________________________________________________________
485 AliTRDTKDInterpolator::AliTRDTKDNodeInfo &AliTRDTKDInterpolator::AliTRDTKDNodeInfo::operator=(const AliTRDTKDInterpolator::AliTRDTKDNodeInfo &ref)
488 // Assignment operator
491 if(&ref==this) return *this;
492 memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
493 memcpy(fBounds, ref.fBounds, 2*fNDim*sizeof(Float_t));
495 fVal[0] = ref.fVal[0];
496 fVal[1] = ref.fVal[1];
499 fPar=new Double_t[fNPar];
500 memcpy(fPar, ref.fPar, fNPar*sizeof(Double_t));
503 fCov=new Double_t[fNCov];
504 memcpy(fCov, ref.fCov, fNCov*sizeof(Double_t));
510 //_________________________________________________________________
511 void AliTRDTKDInterpolator::AliTRDTKDNodeInfo::Print(const Option_t *opt) const
513 // Print the content of the node
515 for(int idim=0; idim<fNDim; idim++) printf("%f ", fData?fData[idim]:0.);
516 printf("] f = [%e +- %e]\n", fVal[0], fVal[1]);
520 for(int idim=0; idim<fNDim; idim++) printf("(%f %f) ", fBounds[2*idim], fBounds[2*idim+1]);
523 if(strcmp(opt, "a")!=0) return;
526 printf("Fit parameters : \n");
527 for(int ip=0; ip<fNPar; ip++) printf("p%d[%e] ", ip, fPar[ip]);
531 for(int ip(0), n(0); ip<fNPar; ip++){
532 for(int jp(ip); jp<fNPar; jp++) printf("c(%d %d)[%e] ", ip, jp, fCov[n++]);
537 //_________________________________________________________________
538 void AliTRDTKDInterpolator::AliTRDTKDNodeInfo::Store(TVectorD const *par, TMatrixD const *cov,Bool_t storeCov)
540 // Store the parameters and the covariance in the node
542 AliDebug(2,"Store Node Interpolation Parameters");
544 if((AliLog::GetDebugLevel("",IsA()->GetName()))>=10){
549 if(!fPar){fPar = new Double_t[fNPar];}
550 for(int ip=0; ip<fNPar; ip++) fPar[ip] = (*par)[ip];
551 if(!cov||!storeCov) return;
552 if(!fCov){fCov = new Double_t[fNCov];}
553 for(int ip(0), np(0); ip<fNPar; ip++)
554 for(int jp=ip; jp<fNPar; jp++) fCov[np++] = (*cov)(ip,jp);
556 if((AliLog::GetDebugLevel("",IsA()->GetName()))>10){this->Print("a");}
559 //_________________________________________________________________
560 Bool_t AliTRDTKDInterpolator::AliTRDTKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t &error,TRDTKDMode mod) const
562 // Recalculate the PDF for one node from the results of interpolation (parameters and covariance matrix)
564 result =0.; error = 1.;
573 AliDebug(1,"Requested Interpolation Parameters don't exist");
577 Double_t fdfdp[fNDim];
580 for(int idim=0; idim<fNDim; idim++){
581 fdfdp[ipar++] = point[idim];
582 for(int jdim=idim; jdim<fNDim; jdim++) fdfdp[ipar++] = point[idim]*point[jdim];
585 // calculate estimation
586 for(int i=0; i<fNPar; i++) result += fdfdp[i]*fPar[i];
589 AliDebug(3,"Interpolation Error cannot be estimated, Covariance Parameters don't exist");
595 for(int i(0), n(0); i<fNPar; i++){
596 error += fdfdp[i]*fdfdp[i]*fCov[n++];
597 for(int j(i+1); j<fNPar; j++) error += 2.*fdfdp[i]*fdfdp[j]*fCov[n++];
599 if(error>0)error = TMath::Sqrt(error);
606 if(error==1)AliWarning("Covariance not available, always choose bin content");
612 // Boundary condition
614 AliDebug(2,"Using Node Value to ensure Boundary Condition");
619 AliDebug(1,Form("Cook PDF Result: %e Error: %e",result,error));
624 //_____________________________________________________________________
625 Bool_t AliTRDTKDInterpolator::AliTRDTKDNodeInfo::Has(const Float_t *p) const
628 for(int id=0; id<fNDim; id++)
629 if(p[id]>=fBounds[2*id] && p[id]<fBounds[2*id+1]) n++;
630 if(n==fNDim) return kTRUE;