-Double_t TKDInterpolator::Eval(const Double_t *point, Double_t &result, Double_t &error)
-{
-// Evaluate PDF for "point". The result is returned in "result" and error in "error". The function returns the chi2 of the fit.
-//
-// Observations:
-//
-// 1. The default method used for interpolation is kCOG.
-// 2. The initial number of neighbors used for the estimation is set to Int(alpha*fLambda) (alpha = 1.5)
-
- Float_t pointF[50]; // local Float_t conversion for "point"
- for(int idim=0; idim<fNDim; idim++) pointF[idim] = (Float_t)point[idim];
- Int_t node = FindNode(pointF) - fNnodes;
- if(fPDFstatus && (fStatus&1) && fPDFstatus[node]) return CookPDF(point, node, result, error);
-
- // Allocate memory
- if(!fBuffer) fBuffer = new Double_t[2*fLambda];
- if(!fKDhelper) fKDhelper = new TKDTreeIF(fNTNodes, fNDim, 30, fRefPoints);
- if(!fFitter) SetIntInterpolation(kFALSE);
-
- // generate parabolic for nD
- //Float_t alpha = Float_t(2*lambda + 1) / fNTNodes; // the bandwidth or smoothing parameter
- //Int_t npoints = Int_t(alpha * fNTNodes);
- //printf("Params : %d NPoints %d\n", lambda, npoints);
- // prepare workers
-
- Int_t *index, // indexes of NN
- ipar, // local looping variable
- npoints = Int_t(1.5*fLambda); // number of data points used for interpolation
- Float_t *dist, // distances of NN
- d, // NN normalized distance
- w0, // work
- w; // tri-cubic weight function
- Double_t sig // bucket error
- = TMath::Sqrt(1./fBucketSize);
- do{
- // find nearest neighbors
- for(int idim=0; idim<fNDim; idim++) pointF[idim] = (Float_t)point[idim];
- if(!fKDhelper->FindNearestNeighbors(pointF, npoints+1, index, dist)){
- Error("Eval()", Form("Failed retriving %d neighbours for point:", npoints));
- for(int idim=0; idim<fNDim; idim++) printf("%f ", point[idim]);
- printf("\n");
- return -1;
- }
- // add points to fitter
- fFitter->ClearPoints();
- for(int in=0; in<npoints; in++){
- if(fStatus&1){ // INT
- for(int idim=0; idim<fNDim; idim++) pointF[idim] = fRefPoints[idim][index[in]];
- Float_t *bounds = GetBoundary(FindNode(pointF));
-
- ipar = 0;
- for(int idim=0; idim<fNDim; idim++){
- fBuffer[ipar++] = .5*(bounds[2*idim] + bounds[2*idim+1]);
- fBuffer[ipar++] = (bounds[2*idim]*bounds[2*idim] + bounds[2*idim] * bounds[2*idim+1] + bounds[2*idim+1] * bounds[2*idim+1])/3.;
- for(int jdim=idim+1; jdim<fNDim; jdim++) fBuffer[ipar++] = (bounds[2*idim] + bounds[2*idim+1]) * (bounds[2*jdim] + bounds[2*jdim+1]) * .25;
- }
- } else { // COG
- for(int idim=0; idim<fNDim; idim++) fBuffer[idim] = fRefPoints[idim][index[in]];
- }
-
- // calculate tri-cubic weighting function
- if(fStatus&4){
- d = dist[in]/ dist[npoints];
- w0 = (1. - d*d*d); w = w0*w0*w0;
- } else w = 1.;
-
- //for(int idim=0; idim<fNDim; idim++) printf("%f ", fBuffer[idim]);
- //printf("\nd[%f] w[%f] sig[%f]\n", d, w, sig);
- fFitter->AddPoint(fBuffer, fRefValues[index[in]], fRefValues[index[in]]*sig/w);
- }
- npoints += 4;
- } while(fFitter->Eval());
-
- // retrive fitter results
- TMatrixD cov(fLambda, fLambda);
- TVectorD par(fLambda);
- fFitter->GetCovarianceMatrix(cov);
- fFitter->GetParameters(par);
- Double_t chi2 = fFitter->GetChisquare()/(npoints - 4 - fLambda);
-
- // store results
- if(fStatus&2 && fStatus&1){
- fPar[node] = par;
- fCov[node] = cov;
- fPDFstatus[node] = 1;
- }
-
- // Build df/dpi|x values
- Double_t *fdfdp = &fBuffer[fLambda];
- ipar = 0;
- fdfdp[ipar++] = 1.;
- for(int idim=0; idim<fNDim; idim++){
- fdfdp[ipar++] = point[idim];
- for(int jdim=idim; jdim<fNDim; jdim++) fdfdp[ipar++] = point[idim]*point[jdim];
- }
-
- // calculate estimation
- result =0.; error = 0.;
- for(int i=0; i<fLambda; i++){
- result += fdfdp[i]*par(i);
- for(int j=0; j<fLambda; j++) error += fdfdp[i]*fdfdp[j]*cov(i,j);
- }
- error = TMath::Sqrt(error);
-
- return chi2;
-}
-
-// //_________________________________________________________________
-// Double_t TKDInterpolator::Eval1(const Double_t *point, Int_t npoints, Double_t &result, Double_t &error)
-// {
-// // Evaluate PDF at k-dimensional position "point". The initial number of
-// // neighbour estimation points is set to "npoints". The default method
-// // used for interpolation is kCOG.
-//
-// // calculate number of parameters in the parabolic expresion
-// Int_t lambda = 1 + fNDim + fNDim*(fNDim+1)/2;
-//
-// if(!fBuffer) fBuffer = new Double_t[lambda-1];
-// if(!fKDhelper) fKDhelper = new TKDTreeIF(GetNTerminalNodes(), fNDim, npoints, fRefPoints);
-//
-// if(!fFitter) fFitter = new TLinearFitter(lambda, Form("hyp%d", fNDim+1));
-// else fFitter->SetFormula(Form("hyp%d", fNDim+1));
-//
-//
-// Float_t pointF[50];
-// for(int idim=0; idim<fNDim; idim++) pointF[idim] = point[idim];
-// Int_t istart = 0;
-// Int_t *index, ipar;
-// Float_t *bounds, *dist, *w = new Float_t[fNDim];
-// Double_t uncertainty = TMath::Sqrt(1./fBucketSize);
-// fFitter->ClearPoints();
-// do{
-// if(!fKDhelper->FindNearestNeighbors(pointF, npoints+1, index, dist)){
-// Error("Eval()", Form("Failed retriving %d neighbours for point:", npoints));
-// for(int idim=0; idim<fNDim; idim++) printf("%f ", point[idim]);
-// printf("\n");
-// return -1;
-// }
-// for(int in=istart; in<npoints; in++){
-// for(int idim=0; idim<fNDim; idim++) w[idim] = fRefPoints[idim][index[in]];
-// bounds = GetBoundary(FindNode(w));
-//
-// ipar = 0;
-// for(int idim=0; idim<fNDim; idim++){
-// fBuffer[ipar++] = .5*(bounds[2*idim] + bounds[2*idim+1]);
-// fBuffer[ipar++] = (bounds[2*idim]*bounds[2*idim] + bounds[2*idim] * bounds[2*idim+1] + bounds[2*idim+1] * bounds[2*idim+1])/3.;
-// for(int jdim=idim+1; jdim<fNDim; jdim++) fBuffer[ipar++] = (bounds[2*idim] + bounds[2*idim+1]) * (bounds[2*jdim] + bounds[2*jdim+1]) * .25;
-// }
-//
-// fFitter->AddPoint(fBuffer, fRefValues[index[in]], fRefValues[index[in]]*uncertainty);
-// }
-// istart = npoints;
-// npoints += 4;
-// } while(fFitter->Eval());
-// delete [] w;
-//
-// // calculate evaluation
-// // fFitter->PrintResults(3);
-// TMatrixD cov(lambda, lambda);
-// TVectorD par(lambda);
-// fFitter->GetCovarianceMatrix(cov);
-// fFitter->GetParameters(par);
-//
-// // Build temporary array to keep values df/dpi|x
-// Double_t f[100];
-// ipar = 0;
-// f[ipar++] = 1.;
-// for(int idim=0; idim<fNDim; idim++){
-// f[ipar++] = point[idim];
-// for(int jdim=idim; jdim<fNDim; jdim++) f[ipar++] = point[idim]*point[jdim];
-// }
-// result =0.; error = 0.;
-// for(int i=0; i<lambda; i++){
-// result += f[i]*par[i];
-// for(int j=0; j<lambda; j++) error += f[i]*f[j]*cov(i,j);
-// }
-// error = TMath::Sqrt(error);
-// Double_t chi2 = fFitter->GetChisquare()/(npoints - 4 - lambda);
-//
-// for(int ipar=0; ipar<lambda; ipar++) printf("%d %8.6e %8.6e\n", ipar, par[ipar], TMath::Sqrt(cov(ipar, ipar)));
-// printf("result %6.3f +- %6.3f [%f]\n", result, error, chi2);
-// return chi2;
-// }
-
-
-//_________________________________________________________________
-void TKDInterpolator::DrawNodes(UInt_t ax1, UInt_t ax2, Int_t depth)