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