STEER/STEERBase/AliTRDTKDInterpolator.cxx

   1 #include "AliTRDTKDInterpolator.h"
   2
   3 #include "TClonesArray.h"
   4 #include "TLinearFitter.h"
   5 #include "TMath.h"
   6 #include "TRandom.h"
   7
   8 #include "AliLog.h"
   9
  10 #include "iostream"
  11 using namespace std;
  12
  13 ClassImp(AliTRDTKDInterpolator)
  14
  15 //_________________________________________________________________
  16 AliTRDTKDInterpolator::AliTRDTKDInterpolator() :
  17 TKDTreeIF(),
  18 fNDataNodes(0),
  19 fNodes(NULL),
  20 fLambda(0),
  21 fNPointsI(0),
  22 fUseHelperNodes(kFALSE),
  23 fUseWeights(kFALSE),
  24 fPDFMode(kInterpolation),
  25 fStoreCov(kFALSE)
  26 {
  27   // default constructor
  28 }
  29
  30 //_________________________________________________________________
  31 AliTRDTKDInterpolator::AliTRDTKDInterpolator(Int_t npoints, Int_t ndim, UInt_t bsize, Float_t **data) :
  32 TKDTreeIF(npoints, ndim, bsize, data),
  33 fNDataNodes(0),
  34 fNodes(NULL),
  35 fLambda(1 + ndim + (ndim*(ndim+1)>>1)),
  36 fNPointsI(100),
  37 fUseHelperNodes(kFALSE),
  38 fUseWeights(kFALSE),
  39 fPDFMode(kInterpolation),
  40 fStoreCov(kFALSE)
  41 {
  42 }
  43
  44 //_________________________________________________________________
  45 AliTRDTKDInterpolator::~AliTRDTKDInterpolator()
  46 {
  47     if(fNodes){
  48         fNodes->Delete();
  49         delete fNodes;
  50         fNodes=NULL;
  51     }
  52 }
  53 //_________________________________________________________________
  54 AliTRDTKDInterpolator::AliTRDTKDInterpolator(const AliTRDTKDInterpolator &ref):
  55 TKDTreeIF(),
  56 fNDataNodes(ref.fNDataNodes),
  57 fNodes(ref.fNodes),
  58 fLambda(ref.fLambda),
  59 fNPointsI(ref.fNPointsI),
  60 fUseHelperNodes(ref.fUseHelperNodes),
  61 fUseWeights(ref.fUseWeights),
  62 fPDFMode(ref.fPDFMode),
  63 fStoreCov(ref.fStoreCov)
  64 {
  65     // Copy constructor
  66     this->Print("");
  67 }
  68
  69 //____________________________________________________________
  70
  71 AliTRDTKDInterpolator &AliTRDTKDInterpolator::operator=(const AliTRDTKDInterpolator &ref){
  72     //
  73     // Assignment operator
  74     //
  75     if(this == &ref) return *this;
  76
  77     // Make copy
  78     TObject::operator=(ref);
  79
  80     this->Print("");
  81
  82     return *this;
  83 }
  84
  85 //_________________________________________________________________
  86 Bool_t AliTRDTKDInterpolator::Build()
  87 {
  88     TKDTreeIF::Build();
  89     if(!fBoundaries) MakeBoundaries();
  90
  91     // allocate interpolation nodes
  92     fNDataNodes = fNPoints/fBucketSize + ((fNPoints%fBucketSize)?1:0);
  93
  94     if(fNodes){
  95         Warning("AliTRDTKDInterpolator::Build()", "Data already allocated.");
  96         fNodes->Delete();
  97     } else {
  98         fNodes = new TClonesArray("AliTRDTKDInterpolator::AliTRDTKDNodeInfo", fNDataNodes);
  99         fNodes->SetOwner();
 100     }
 101     for(int in=0; in<fNDataNodes; in++) new ((*fNodes)[in]) AliTRDTKDNodeInfo(fNDim);
 102
 103     // Set Interpolator nodes
 104
 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));
 111
 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;
 118         }
 119     }
 120
 121     // Analyze last (incomplete) terminal node
 122
 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));
 129
 130     for(int idim=0; idim<fNDim; idim++){
 131         Float_t dx = ftnode->fBounds[2*idim+1]-ftnode->fBounds[2*idim];
 132         if(dx < 1.e-30){
 133             Warning("AliTRDTKDInterpolator::Build()", "Terminal bucket index[%d] too narrow on the %d dimension.", inode, idim);
 134             continue;
 135         }
 136         ftnode->fVal[0] /= (ftnode->fBounds[2*idim+1] - ftnode->fBounds[2*idim]);
 137     }
 138     ftnode->fVal[1] =  ftnode->fVal[0]/TMath::Sqrt(float(counts));
 139
 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;
 146     }
 147
 148     delete [] fBoundaries;
 149     fBoundaries = NULL;
 150     // Add Helper Nodes
 151     if(fUseHelperNodes){BuildBoundaryNodes();}
 152
 153     if(fNPointsI>GetNTNodes()){fNPointsI=GetNTNodes();}
 154
 155     BuildInterpolation();
 156
 157     return kTRUE;
 158 }
 159
 160
 161 //_________________________________________________________________
 162 Bool_t AliTRDTKDInterpolator::Eval(const Double_t *point, Double_t &result, Double_t &error)
 163 {
 164     AliDebug(3,Form("Eval PDF Mode %d",fPDFMode));
 165     if((AliLog::GetDebugLevel("",IsA()->GetName()))>0){
 166         printf("Point [");
 167         for(int idim=0; idim<fNDim; idim++) printf("%f ", point[idim]);
 168         printf("] \n");
 169     }
 170
 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);
 174     if(nodeIndex<0){
 175         AliError("Can not retrieve node for data point");
 176         result = 0.;
 177         error = 1.E10;
 178         return kFALSE;
 179     }
 180     AliTRDTKDNodeInfo *node =GetNodeInfo(nodeIndex);
 181
 182     if((AliLog::GetDebugLevel("",IsA()->GetName()))>0){
 183         printf("Node Info: \n");
 184         node->Print("a");
 185     }
 186
 187     return node->CookPDF(point, result, error,fPDFMode);
 188 }
 189
 190 //__________________________________________________________________
 191 void AliTRDTKDInterpolator::Print(const Option_t */*opt*/) const
 192 {
 193     for(Int_t i=GetNTNodes(); i--;){
 194         printf("Node %d of %d: \n",i,GetNTNodes());
 195         GetNodeInfo(i)->Print();
 196     }
 197
 198 }
 199
 200 //__________________________________________________________________
 201 Int_t AliTRDTKDInterpolator::GetNodeIndex(const Float_t *p)
 202 {
 203     Int_t inode=FindNode(p)-fNDataNodes+1;
 204     if(GetNodeInfo(inode)->Has(p)){
 205         AliDebug(2,Form("Find Node %d",inode));
 206         return inode;
 207     }
 208
 209     // Search extra nodes
 210
 211     for(inode=fNDataNodes;inode<GetNTNodes();inode++){
 212         if(GetNodeInfo(inode)->Has(p)){
 213             AliDebug(2,Form("Find Extra Node %d",inode));
 214             return inode;
 215         }
 216     }
 217
 218     // Search for nearest neighbor
 219     Float_t dist;
 220     Float_t closestdist=10000;
 221     inode=-1;
 222     for(Int_t ii=0;ii<GetNTNodes();ii++){
 223         AliTRDTKDNodeInfo *node=GetNodeInfo(ii);
 224         dist=0;
 225         for(Int_t idim=0;idim<fNDim;idim++){
 226             dist+=TMath::Power((node->fData[idim]-p[idim]),2);
 227         }
 228         dist=TMath::Sqrt(dist);
 229         if(dist<closestdist){closestdist=dist;inode=ii;}
 230     }
 231     AliDebug(2,Form("Find Nearest Neighbor Node %d",inode));
 232     return inode;
 233 }
 234
 235 //_________________________________________________________________
 236 AliTRDTKDInterpolator::AliTRDTKDNodeInfo* AliTRDTKDInterpolator::GetNodeInfo(Int_t inode) const
 237 {
 238   if(!fNodes || inode >= GetNTNodes()) return NULL;
 239   return (AliTRDTKDNodeInfo*)(*fNodes)[inode];
 240 }
 241
 242 //_________________________________________________________________
 243 Int_t AliTRDTKDInterpolator::GetNTNodes() const
 244 {
 245   return fNodes?fNodes->GetEntriesFast():0;
 246 }
 247
 248 //_________________________________________________________________
 249 Bool_t AliTRDTKDInterpolator::GetRange(Int_t idim,Float_t range[2]) const
 250 {
 251     if(!fNodes) return kFALSE;
 252     if(idim<0 || idim>=fNDim){
 253     range[0]=0.; range[1]=0.;
 254     return kFALSE;
 255     }
 256     range[0]=1.e10; range[1]=-1.e10;
 257     for(Int_t in=GetNTNodes(); in--; ){
 258         AliTRDTKDNodeInfo *node = GetNodeInfo(in);
 259
 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];
 262     }
 263
 264     return kTRUE;
 265 }
 266
 267 //_________________________________________________________________
 268 TH2Poly *AliTRDTKDInterpolator::Projection(Int_t xdim,Int_t ydim)
 269 {
 270     Float_t rangex[2],rangey[2];
 271     GetRange(xdim,rangex);
 272     GetRange(ydim,rangey);
 273
 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));
 277
 278     for(Int_t inode=0;inode<GetNTNodes();inode++){
 279
 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]);
 283     }
 284     return h2;
 285 }
 286
 287 //_________________________________________________________________
 288 void AliTRDTKDInterpolator::BuildInterpolation()
 289 {
 290     AliInfo("Build Interpolation");
 291
 292     // Calculate Interpolation
 293
 294     Double_t buffer[fLambda];
 295
 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];
 300     }
 301     TKDTreeIF *KDhelper = new TKDTreeIF(GetNTNodes(), fNDim, 30, refPoints);
 302     KDhelper->Build();
 303     KDhelper->MakeBoundariesExact();
 304
 305     Float_t dist[fNPointsI];
 306     Int_t ind[fNPointsI];
 307
 308     TLinearFitter fitter(fLambda, Form("hyp%d", fLambda-1));
 309
 310     Int_t nodeIndex(0); Float_t param[6], *pp(NULL);
 311     nodeIndex=GetNTNodes(); pp=&param[0];
 312     while(nodeIndex--){
 313
 314         fitter.ClearPoints();
 315
 316         AliTRDTKDNodeInfo *node = GetNodeInfo(nodeIndex);
 317         // find nearest neighbors
 318         KDhelper->FindNearestNeighbors(node->fData,fNPointsI, &ind[0], &dist[0]);
 319
 320         for(int in=0;in<fNPointsI;in++){
 321             AliTRDTKDNodeInfo *nnode = GetNodeInfo(ind[in]);
 322
 323             Float_t w=1; //weight
 324             // calculate tri-cubic weighting function
 325             if(fUseWeights){
 326                 Float_t d = dist[in]/dist[fNPointsI-1];
 327                 Float_t w0 = (1. - d*d*d);
 328                 w = w0*w0*w0;
 329                 if(w<1.e-30) continue;
 330             }
 331             Int_t ipar=0;
 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];
 335             }
 336             fitter.AddPoint(buffer,nnode->fVal[0], nnode->fVal[1]/w);
 337
 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));
 343                     zdata[kdim]=0;
 344                     ipar=0;
 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];
 348                     }
 349                     fitter.AddPoint(buffer,0,1);
 350                 }
 351             }
 352         }
 353
 354         AliDebug(2,Form("Calculate Interpolation for Node %d",nodeIndex));
 355         fitter.Eval();
 356
 357         // retrive fitter results
 358         TMatrixD cov(fLambda, fLambda);
 359         TVectorD par(fLambda);
 360         fitter.GetCovarianceMatrix(cov);
 361         fitter.GetParameters(par);
 362
 363         // store results
 364         node->Store(&par,&cov,fStoreCov);
 365     }
 366
 367     delete KDhelper;
 368     for(int id=0; id<fNDim; id++){
 369         delete refPoints[id];
 370     }
 371     delete[] refPoints;
 372 }
 373
 374 //_________________________________________________________________
 375 void AliTRDTKDInterpolator::BuildBoundaryNodes(){
 376
 377     Int_t nnew=0;
 378
 379     Float_t treebounds[2*fNDim];
 380     for(Int_t idim=0;idim<fNDim;idim++){
 381         GetRange(idim,&treebounds[2*idim]);
 382     }
 383
 384     for(int inode=0; inode<GetNTNodes(); inode++){
 385
 386         AliTRDTKDNodeInfo *node=GetNodeInfo(inode);
 387
 388         for(Int_t vdim=0;vdim<fNDim;vdim++){
 389
 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]){
 393
 394                     // Add new Node
 395                     new ((*fNodes)[GetNTNodes()]) AliTRDTKDNodeInfo(fNDim);
 396
 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];
 405                     }
 406                     newnode->fVal[0]=0;
 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]);
 411                     }
 412                     nnew++;
 413                 }
 414             }
 415         }
 416     }
 417     AliInfo(Form("%d Boundary Nodes Added \n",nnew));
 418 }
 419
 420 //_________________________________________________________________
 421 AliTRDTKDInterpolator::AliTRDTKDNodeInfo::AliTRDTKDNodeInfo(Int_t ndim):
 422   TObject()
 423   ,fNDim(ndim)
 424   ,fNBounds(2*ndim)
 425   ,fNPar(1 + ndim + (ndim*(ndim+1)>>1))
 426   ,fNCov(Int_t((fNPar+1)*fNPar/2))
 427   ,fData(NULL)
 428   ,fBounds(NULL)
 429   ,fPar(NULL)
 430   ,fCov(NULL)
 431 {
 432   // Default constructor
 433     fVal[0] = 0.; fVal[1] = 0.;
 434     fData=new Float_t[fNDim];
 435     fBounds=new Float_t[fNBounds];
 436 }
 437
 438 //_________________________________________________________________
 439 AliTRDTKDInterpolator::AliTRDTKDNodeInfo::AliTRDTKDNodeInfo(const AliTRDTKDInterpolator::AliTRDTKDNodeInfo &ref):
 440 TObject(ref)
 441    ,fNDim(ref.fNDim)
 442    ,fNBounds(ref.fNBounds)
 443    ,fNPar(ref.fNPar)
 444    ,fNCov(ref.fNCov)
 445    ,fData(NULL)
 446    ,fBounds(NULL)
 447    ,fPar(NULL)
 448    ,fCov(NULL)
 449 {
 450   // Copy constructor
 451
 452     if(ref.fData){
 453         fData = new Float_t[fNDim];
 454         memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
 455     }
 456     if(ref.fBounds){
 457         fBounds = new Float_t[2*fNDim];
 458         memcpy(fBounds, ref.fBounds, 2*fNDim*sizeof(Float_t));
 459     }
 460
 461     fVal[0] = ref.fVal[0];
 462     fVal[1] = ref.fVal[1];
 463
 464     if(ref.fPar){
 465         fPar=new Double_t[fNPar];
 466         memcpy(fPar, ref.fPar, fNPar*sizeof(Double_t));
 467     }
 468     if(ref.fCov){
 469         fCov=new Double_t[fNCov];
 470         memcpy(fCov, ref.fCov, fNCov*sizeof(Double_t));
 471     }
 472 }
 473
 474 //_________________________________________________________________
 475 AliTRDTKDInterpolator::AliTRDTKDNodeInfo::~AliTRDTKDNodeInfo()
 476 {
 477   // Destructor
 478   if(fData) delete [] fData;
 479   if(fBounds) delete [] fBounds;
 480   if(fPar) delete [] fPar;
 481   if(fCov) delete [] fCov;
 482 }
 483
 484 //_________________________________________________________________
 485 AliTRDTKDInterpolator::AliTRDTKDNodeInfo &AliTRDTKDInterpolator::AliTRDTKDNodeInfo::operator=(const AliTRDTKDInterpolator::AliTRDTKDNodeInfo &ref)
 486 {
 487     //
 488     // Assignment operator
 489     //
 490
 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));
 494
 495     fVal[0] = ref.fVal[0];
 496     fVal[1] = ref.fVal[1];
 497
 498     if(ref.fPar){
 499         fPar=new Double_t[fNPar];
 500         memcpy(fPar, ref.fPar, fNPar*sizeof(Double_t));
 501     }
 502     if(ref.fCov){
 503         fCov=new Double_t[fNCov];
 504         memcpy(fCov, ref.fCov, fNCov*sizeof(Double_t));
 505     }
 506     return *this;
 507 }
 508
 509
 510 //_________________________________________________________________
 511 void AliTRDTKDInterpolator::AliTRDTKDNodeInfo::Print(const Option_t *opt) const
 512 {
 513   // Print the content of the node
 514   printf("x [");
 515   for(int idim=0; idim<fNDim; idim++) printf("%f ", fData?fData[idim]:0.);
 516   printf("] f = [%e +- %e]\n", fVal[0], fVal[1]);
 517
 518   if(fBounds){
 519       printf("range [");
 520       for(int idim=0; idim<fNDim; idim++) printf("(%f %f) ", fBounds[2*idim], fBounds[2*idim+1]);
 521     printf("]\n");
 522   }
 523   if(strcmp(opt, "a")!=0) return;
 524
 525   if(fPar){
 526     printf("Fit parameters : \n");
 527     for(int ip=0; ip<fNPar; ip++) printf("p%d[%e] ", ip, fPar[ip]);
 528     printf("\n");
 529   }
 530   if(!fCov) return;
 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++]);
 533     printf("\n");
 534   }
 535 }
 536
 537 //_________________________________________________________________
 538 void AliTRDTKDInterpolator::AliTRDTKDNodeInfo::Store(TVectorD const *par, TMatrixD const *cov,Bool_t storeCov)
 539 {
 540     // Store the parameters and the covariance in the node
 541
 542     AliDebug(2,"Store Node Interpolation Parameters");
 543
 544      if((AliLog::GetDebugLevel("",IsA()->GetName()))>=10){
 545         par->Print("");
 546         cov->Print("");
 547     }
 548
 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);
 555
 556      if((AliLog::GetDebugLevel("",IsA()->GetName()))>10){this->Print("a");}
 557 }
 558
 559 //_________________________________________________________________
 560 Bool_t AliTRDTKDInterpolator::AliTRDTKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t &error,TRDTKDMode mod) const
 561 {
 562     // Recalculate the PDF for one node from the results of interpolation (parameters and covariance matrix)
 563
 564     result =0.; error = 1.;
 565
 566     if(mod==kNodeVal){
 567         error=fVal[1];
 568         result=fVal[0];
 569         return kTRUE;
 570     }
 571
 572     if(!fPar){
 573         AliDebug(1,"Requested Interpolation Parameters don't exist");
 574         return kFALSE;
 575     }
 576
 577     Double_t fdfdp[fNDim];
 578     Int_t ipar = 0;
 579     fdfdp[ipar++] = 1.;
 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];
 583     }
 584
 585     // calculate estimation
 586     for(int i=0; i<fNPar; i++) result += fdfdp[i]*fPar[i];
 587
 588     if(!fCov){
 589         AliDebug(3,"Interpolation Error cannot be estimated, Covariance Parameters don't exist");
 590         return kTRUE;
 591     }
 592     // calculate error
 593     error=0;
 594
 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++];
 598     }
 599     if(error>0)error = TMath::Sqrt(error);
 600     else{error=100;}
 601
 602     if(mod==kMinError){
 603         if(error<fVal[1]){
 604             return kTRUE;
 605         }
 606         if(error==1)AliWarning("Covariance not available, always choose bin content");
 607         error=fVal[1];
 608         result=fVal[0];
 609         return kTRUE;
 610     }
 611
 612     // Boundary condition
 613     if(result<0){
 614         AliDebug(2,"Using Node Value to ensure Boundary Condition");
 615         result=fVal[0];
 616         error=fVal[1];
 617     }
 618
 619     AliDebug(1,Form("Cook PDF Result: %e Error: %e",result,error));
 620
 621     return kTRUE;
 622 }
 623
 624 //_____________________________________________________________________
 625 Bool_t AliTRDTKDInterpolator::AliTRDTKDNodeInfo::Has(const Float_t *p) const
 626 {
 627   Int_t n(0);
 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;
 631   return kFALSE;
 632 }
 633