[u/mrichter/AliRoot.git] / STEER / STEERBase / AliTRDTKDInterpolator.cxx

#include "AliTRDTKDInterpolator.h"

#include "TClonesArray.h"
#include "TLinearFitter.h"
#include "TMath.h"
#include "TRandom.h"

#include "AliLog.h"

#include "iostream"
using namespace std;

ClassImp(AliTRDTKDInterpolator)

//_________________________________________________________________
AliTRDTKDInterpolator::AliTRDTKDInterpolator() :
TKDTreeIF(),
fNDataNodes(0),
fNodes(NULL),
fLambda(0),
fNPointsI(0),
fUseHelperNodes(kFALSE),
fUseWeights(kFALSE),
fPDFMode(kInterpolation),
fStoreCov(kFALSE)
{
  // default constructor
}

//_________________________________________________________________
AliTRDTKDInterpolator::AliTRDTKDInterpolator(Int_t npoints, Int_t ndim, UInt_t bsize, Float_t **data) :
TKDTreeIF(npoints, ndim, bsize, data),
fNDataNodes(0),
fNodes(NULL),
fLambda(1 + ndim + (ndim*(ndim+1)>>1)),
fNPointsI(100),
fUseHelperNodes(kFALSE),
fUseWeights(kFALSE),
fPDFMode(kInterpolation),
fStoreCov(kFALSE)
{
}

//_________________________________________________________________
AliTRDTKDInterpolator::~AliTRDTKDInterpolator()
{
    if(fNodes){
	fNodes->Delete();
	delete fNodes;
	fNodes=NULL;
    }
}
//_________________________________________________________________
AliTRDTKDInterpolator::AliTRDTKDInterpolator(const AliTRDTKDInterpolator &ref):
TKDTreeIF(),
fNDataNodes(ref.fNDataNodes),
fNodes(ref.fNodes),
fLambda(ref.fLambda),
fNPointsI(ref.fNPointsI),
fUseHelperNodes(ref.fUseHelperNodes),
fUseWeights(ref.fUseWeights),
fPDFMode(ref.fPDFMode),
fStoreCov(ref.fStoreCov)
{
    // Copy constructor
    this->Print("");
}

//____________________________________________________________

AliTRDTKDInterpolator &AliTRDTKDInterpolator::operator=(const AliTRDTKDInterpolator &ref){
    //
    // Assignment operator
    //
    if(this == &ref) return *this;

    // Make copy
    TObject::operator=(ref);

    this->Print("");

    return *this;
}

//_________________________________________________________________
Bool_t AliTRDTKDInterpolator::Build()
{
    TKDTreeIF::Build();
    if(!fBoundaries) MakeBoundaries();

    // allocate interpolation nodes
    fNDataNodes = fNPoints/fBucketSize + ((fNPoints%fBucketSize)?1:0);

    if(fNodes){
	Warning("AliTRDTKDInterpolator::Build()", "Data already allocated.");
	fNodes->Delete();
    } else {
	fNodes = new TClonesArray("AliTRDTKDInterpolator::AliTRDTKDNodeInfo", fNDataNodes);
	fNodes->SetOwner();
    }
    for(int in=0; in<fNDataNodes; in++) new ((*fNodes)[in]) AliTRDTKDNodeInfo(fNDim);

    // Set Interpolator nodes

    for(int inode=0, tnode = fNNodes; inode<fNDataNodes-1; inode++, tnode++){
	AliTRDTKDNodeInfo *node =GetNodeInfo(inode);
	memcpy(node->fBounds,GetBoundary(tnode),2*fNDim*sizeof(Float_t));
	node->fVal[0] =  Float_t(fBucketSize)/fNPoints;
	for(int idim=0; idim<fNDim; idim++) node->fVal[0] /= (node->fBounds[2*idim+1] - node->fBounds[2*idim]);
	node->fVal[1] =  node->fVal[0]/TMath::Sqrt(float(fBucketSize));

	Int_t *indexPoints = GetPointsIndexes(tnode);
	// loop points in this terminal node
	for(int idim=0; idim<fNDim; idim++){
	    node->fData[idim] = 0.;
	    for(int ip = 0; ip<fBucketSize; ip++) node->fData[idim] += fData[idim][indexPoints[ip]];
	    node->fData[idim] /= fBucketSize;
	}
    }

    // Analyze last (incomplete) terminal node

    Int_t counts = fNPoints%fBucketSize;
    counts = counts ? counts : fBucketSize;
    Int_t inode = fNDataNodes - 1, tnode = inode + fNNodes;
    AliTRDTKDNodeInfo *ftnode = GetNodeInfo(inode);
    ftnode->fVal[0] =  Float_t(counts)/fNPoints;
    memcpy(ftnode->fBounds,GetBoundary(tnode),2*fNDim*sizeof(Float_t));
 
    for(int idim=0; idim<fNDim; idim++){
	Float_t dx = ftnode->fBounds[2*idim+1]-ftnode->fBounds[2*idim];
	if(dx < 1.e-30){
	    Warning("AliTRDTKDInterpolator::Build()", "Terminal bucket index[%d] too narrow on the %d dimension.", inode, idim);
	    continue;
	}
	ftnode->fVal[0] /= (ftnode->fBounds[2*idim+1] - ftnode->fBounds[2*idim]);
    }
    ftnode->fVal[1] =  ftnode->fVal[0]/TMath::Sqrt(float(counts));

    // loop points in this terminal node
    Int_t *indexPoints = GetPointsIndexes(tnode);
    for(int idim=0; idim<fNDim; idim++){
	ftnode->fData[idim] = 0.;
	for(int ip = 0; ip<counts; ip++) ftnode->fData[idim] += fData[idim][indexPoints[ip]];
	ftnode->fData[idim] /= counts;
    }

    delete [] fBoundaries;
    fBoundaries = NULL;
    // Add Helper Nodes
    if(fUseHelperNodes){BuildBoundaryNodes();}

    if(fNPointsI>GetNTNodes()){fNPointsI=GetNTNodes();}

    BuildInterpolation();

    return kTRUE;
}


//_________________________________________________________________
Bool_t AliTRDTKDInterpolator::Eval(const Double_t *point, Double_t &result, Double_t &error)
{
    AliDebug(3,Form("Eval PDF Mode %d",fPDFMode));
    if((AliLog::GetDebugLevel("",IsA()->GetName()))>0){
	printf("Point [");
	for(int idim=0; idim<fNDim; idim++) printf("%f ", point[idim]);
	printf("] \n");
    }

    Float_t pointF[fNDim]; // local Float_t conversion for "point"
    for(int idim=0; idim<fNDim; idim++) pointF[idim] = (Float_t)point[idim];
    Int_t nodeIndex = GetNodeIndex(pointF);
    if(nodeIndex<0){
	AliError("Can not retrieve node for data point");
	result = 0.;
	error = 1.E10;
	return kFALSE;
    }
    AliTRDTKDNodeInfo *node =GetNodeInfo(nodeIndex);

    if((AliLog::GetDebugLevel("",IsA()->GetName()))>0){
	printf("Node Info: \n");
	node->Print("a");
    }

    return node->CookPDF(point, result, error,fPDFMode);
}

//__________________________________________________________________
void AliTRDTKDInterpolator::Print(const Option_t */*opt*/) const
{
    for(Int_t i=GetNTNodes(); i--;){
	printf("Node %d of %d: \n",i,GetNTNodes());
	GetNodeInfo(i)->Print();
    }

}

//__________________________________________________________________
Int_t AliTRDTKDInterpolator::GetNodeIndex(const Float_t *p)
{
    Int_t inode=FindNode(p)-fNDataNodes+1;
    if(GetNodeInfo(inode)->Has(p)){
	AliDebug(2,Form("Find Node %d",inode));
	return inode;
    }

    // Search extra nodes

    for(inode=fNDataNodes;inode<GetNTNodes();inode++){
	if(GetNodeInfo(inode)->Has(p)){
	    AliDebug(2,Form("Find Extra Node %d",inode));
	    return inode;
	}
    }

    // Search for nearest neighbor
    Float_t dist;
    Float_t closestdist=10000;
    inode=-1;
    for(Int_t ii=0;ii<GetNTNodes();ii++){
	AliTRDTKDNodeInfo *node=GetNodeInfo(ii);
	dist=0;
	for(Int_t idim=0;idim<fNDim;idim++){
	    dist+=TMath::Power((node->fData[idim]-p[idim]),2);
	}
	dist=TMath::Sqrt(dist);
	if(dist<closestdist){closestdist=dist;inode=ii;}
    }
    AliDebug(2,Form("Find Nearest Neighbor Node %d",inode));
    return inode;
}

//_________________________________________________________________
AliTRDTKDInterpolator::AliTRDTKDNodeInfo* AliTRDTKDInterpolator::GetNodeInfo(Int_t inode) const
{
  if(!fNodes || inode >= GetNTNodes()) return NULL;
  return (AliTRDTKDNodeInfo*)(*fNodes)[inode];
}

//_________________________________________________________________
Int_t AliTRDTKDInterpolator::GetNTNodes() const 
{
  return fNodes?fNodes->GetEntriesFast():0;
}

//_________________________________________________________________
Bool_t AliTRDTKDInterpolator::GetRange(Int_t idim,Float_t range[2]) const
{
    if(!fNodes) return kFALSE;
    if(idim<0 || idim>=fNDim){
    range[0]=0.; range[1]=0.;
    return kFALSE;
    }
    range[0]=1.e10; range[1]=-1.e10;
    for(Int_t in=GetNTNodes(); in--; ){
	AliTRDTKDNodeInfo *node = GetNodeInfo(in);

	if(node->fBounds[2*idim]<range[0]) range[0] = node->fBounds[2*idim];
	if(node->fBounds[2*idim+1]>range[1]) range[1] = node->fBounds[2*idim+1];
    }

    return kTRUE;
}

//_________________________________________________________________
TH2Poly *AliTRDTKDInterpolator::Projection(Int_t xdim,Int_t ydim)
{
    Float_t rangex[2],rangey[2];
    GetRange(xdim,rangex);
    GetRange(ydim,rangey);

    TH2Poly* h2 = new TH2Poly("hTKDnodes","hTKDnodes", rangex[0],rangex[1],rangey[0],rangey[1]);
    h2->GetXaxis()->SetTitle(Form("Q_{%d}", xdim));
    h2->GetYaxis()->SetTitle(Form("Q_{%d}", ydim));

    for(Int_t inode=0;inode<GetNTNodes();inode++){

	AliTRDTKDNodeInfo* node=GetNodeInfo(inode);
	h2->AddBin(node->fBounds[2*xdim],node->fBounds[2*ydim],node->fBounds[2*xdim+1],node->fBounds[2*ydim+1]);
	h2->SetBinContent(inode+1, node->fVal[0]);
    }
    return h2;
}

//_________________________________________________________________
void AliTRDTKDInterpolator::BuildInterpolation()
{
    AliInfo("Build Interpolation");

    // Calculate Interpolation

    Double_t buffer[fLambda];

    Float_t **refPoints = new Float_t*[fNDim];
    for(int id=0; id<fNDim; id++){
	refPoints[id] = new Float_t[GetNTNodes()];
	for(int in=0; in<GetNTNodes(); in++) refPoints[id][in] = GetNodeInfo(in)->fData[id];
    }
    TKDTreeIF *KDhelper = new TKDTreeIF(GetNTNodes(), fNDim, 30, refPoints);
    KDhelper->Build();
    KDhelper->MakeBoundariesExact();

    Float_t dist[fNPointsI];
    Int_t ind[fNPointsI];

    TLinearFitter fitter(fLambda, Form("hyp%d", fLambda-1));

    Int_t nodeIndex(0); Float_t param[6], *pp(NULL);
    nodeIndex=GetNTNodes(); pp=&param[0];
    while(nodeIndex--){

	fitter.ClearPoints();

	AliTRDTKDNodeInfo *node = GetNodeInfo(nodeIndex);
	// find nearest neighbors
	KDhelper->FindNearestNeighbors(node->fData,fNPointsI, &ind[0], &dist[0]);

	for(int in=0;in<fNPointsI;in++){
	    AliTRDTKDNodeInfo *nnode = GetNodeInfo(ind[in]);

            Float_t w=1; //weight
	    // calculate tri-cubic weighting function
	    if(fUseWeights){
		Float_t d = dist[in]/dist[fNPointsI-1];
		Float_t w0 = (1. - d*d*d);
		w = w0*w0*w0;
		if(w<1.e-30) continue;
	    }
	    Int_t ipar=0;
	    for(int idim=0; idim<fNDim; idim++){
		buffer[ipar++] = nnode->fData[idim];
		for(int jdim=idim; jdim<fNDim; jdim++) buffer[ipar++] = nnode->fData[idim]*nnode->fData[jdim];
	    }
	    fitter.AddPoint(buffer,nnode->fVal[0], nnode->fVal[1]/w);

	    // Ensure Boundary Condition
	    for(Int_t kdim=0;kdim<fNDim;kdim++){
		if(node->fBounds[2*kdim]==0){
		    Float_t zdata[fNDim];
                    memcpy(&zdata[0],node->fData,fNDim*sizeof(Float_t));
		    zdata[kdim]=0;
		    ipar=0;
		    for(int idim=0; idim<fNDim; idim++){
			buffer[ipar++] = zdata[idim];
			for(int jdim=idim; jdim<fNDim; jdim++) buffer[ipar++] = zdata[idim]*zdata[jdim];
		    }
		    fitter.AddPoint(buffer,0,1);
		}
	    }
	}

	AliDebug(2,Form("Calculate Interpolation for Node %d",nodeIndex));
	fitter.Eval();

	// retrive fitter results
	TMatrixD cov(fLambda, fLambda);
	TVectorD par(fLambda);
	fitter.GetCovarianceMatrix(cov);
	fitter.GetParameters(par);

	// store results
	node->Store(&par,&cov,fStoreCov);
    }

    delete KDhelper;
    for(int id=0; id<fNDim; id++){
	delete refPoints[id];
    }
    delete[] refPoints;
}

//_________________________________________________________________
void AliTRDTKDInterpolator::BuildBoundaryNodes(){

    Int_t nnew=0;

    Float_t treebounds[2*fNDim];
    for(Int_t idim=0;idim<fNDim;idim++){
	GetRange(idim,&treebounds[2*idim]);
    }

    for(int inode=0; inode<GetNTNodes(); inode++){

	AliTRDTKDNodeInfo *node=GetNodeInfo(inode);

	for(Int_t vdim=0;vdim<fNDim;vdim++){

	    // Try expansion to lower and higher values
	    for(Int_t iter=0;iter<2;iter++){
		if(node->fBounds[2*vdim+iter]==treebounds[2*vdim+iter]){

		    // Add new Node
		    new ((*fNodes)[GetNTNodes()]) AliTRDTKDNodeInfo(fNDim);

		    AliTRDTKDNodeInfo *newnode = GetNodeInfo(GetNTNodes()-1);
		    if(iter==0)newnode->fBounds[2*vdim+iter]=0;
                    if(iter==1)newnode->fBounds[2*vdim+iter]=2*treebounds[2*vdim+iter];
		    newnode->fBounds[2*vdim+!iter]=node->fBounds[2*vdim+iter];
		    for(Int_t idim=0;idim<fNDim;idim++){
			if(idim==vdim)continue;
			newnode->fBounds[2*idim]=node->fBounds[2*idim];
			newnode->fBounds[2*idim+1]=node->fBounds[2*idim+1];
		    }
		    newnode->fVal[0]=0;
		    newnode->fVal[1]=Float_t(1)/fNPoints;
		    for(int idim=0; idim<fNDim; idim++){
			newnode->fVal[1] /= (newnode->fBounds[2*idim+1] - newnode->fBounds[2*idim]);
                        newnode->fData[idim]=0.5*(newnode->fBounds[2*idim+1] + newnode->fBounds[2*idim]);
		    }
		    nnew++;
		}
	    }
	}
    }
    AliInfo(Form("%d Boundary Nodes Added \n",nnew));
}

//_________________________________________________________________
AliTRDTKDInterpolator::AliTRDTKDNodeInfo::AliTRDTKDNodeInfo(Int_t ndim):
  TObject()
  ,fNDim(ndim)
  ,fNBounds(2*ndim)
  ,fNPar(1 + ndim + (ndim*(ndim+1)>>1))
  ,fNCov(Int_t((fNPar+1)*fNPar/2))
  ,fData(NULL)
  ,fBounds(NULL)
  ,fPar(NULL)
  ,fCov(NULL)
{
  // Default constructor
    fVal[0] = 0.; fVal[1] = 0.;
    fData=new Float_t[fNDim];
    fBounds=new Float_t[fNBounds];
}

//_________________________________________________________________
AliTRDTKDInterpolator::AliTRDTKDNodeInfo::AliTRDTKDNodeInfo(const AliTRDTKDInterpolator::AliTRDTKDNodeInfo &ref):
TObject(ref)
   ,fNDim(ref.fNDim)
   ,fNBounds(ref.fNBounds)
   ,fNPar(ref.fNPar)
   ,fNCov(ref.fNCov)
   ,fData(NULL)
   ,fBounds(NULL)
   ,fPar(NULL)
   ,fCov(NULL)
{
  // Copy constructor

    if(ref.fData){
	fData = new Float_t[fNDim];
	memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
    }
    if(ref.fBounds){
	fBounds = new Float_t[2*fNDim];
	memcpy(fBounds, ref.fBounds, 2*fNDim*sizeof(Float_t));
    }

    fVal[0] = ref.fVal[0];
    fVal[1] = ref.fVal[1];

    if(ref.fPar){
	fPar=new Double_t[fNPar];
	memcpy(fPar, ref.fPar, fNPar*sizeof(Double_t));
    }
    if(ref.fCov){
	fCov=new Double_t[fNCov];
	memcpy(fCov, ref.fCov, fNCov*sizeof(Double_t));
    }
}

//_________________________________________________________________
AliTRDTKDInterpolator::AliTRDTKDNodeInfo::~AliTRDTKDNodeInfo()
{
  // Destructor
  if(fData) delete [] fData;
  if(fBounds) delete [] fBounds;
  if(fPar) delete [] fPar;
  if(fCov) delete [] fCov;
}

//_________________________________________________________________
AliTRDTKDInterpolator::AliTRDTKDNodeInfo &AliTRDTKDInterpolator::AliTRDTKDNodeInfo::operator=(const AliTRDTKDInterpolator::AliTRDTKDNodeInfo &ref)
{
    //
    // Assignment operator
    //

    if(&ref==this) return *this;
    memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
    memcpy(fBounds, ref.fBounds, 2*fNDim*sizeof(Float_t));

    fVal[0] = ref.fVal[0];
    fVal[1] = ref.fVal[1];

    if(ref.fPar){
	fPar=new Double_t[fNPar];
	memcpy(fPar, ref.fPar, fNPar*sizeof(Double_t));
    }
    if(ref.fCov){
	fCov=new Double_t[fNCov];
	memcpy(fCov, ref.fCov, fNCov*sizeof(Double_t));
    }
    return *this;
}


//_________________________________________________________________
void AliTRDTKDInterpolator::AliTRDTKDNodeInfo::Print(const Option_t *opt) const
{
  // Print the content of the node
  printf("x [");
  for(int idim=0; idim<fNDim; idim++) printf("%f ", fData?fData[idim]:0.);
  printf("] f = [%e +- %e]\n", fVal[0], fVal[1]);

  if(fBounds){
      printf("range [");
      for(int idim=0; idim<fNDim; idim++) printf("(%f %f) ", fBounds[2*idim], fBounds[2*idim+1]);
    printf("]\n");
  }
  if(strcmp(opt, "a")!=0) return;

  if(fPar){
    printf("Fit parameters : \n");
    for(int ip=0; ip<fNPar; ip++) printf("p%d[%e] ", ip, fPar[ip]);
    printf("\n");
  }
  if(!fCov) return;
  for(int ip(0), n(0); ip<fNPar; ip++){
    for(int jp(ip); jp<fNPar; jp++) printf("c(%d %d)[%e] ", ip, jp, fCov[n++]);
    printf("\n");
  }
}

//_________________________________________________________________
void AliTRDTKDInterpolator::AliTRDTKDNodeInfo::Store(TVectorD const *par, TMatrixD const *cov,Bool_t storeCov)
{
    // Store the parameters and the covariance in the node

    AliDebug(2,"Store Node Interpolation Parameters");

     if((AliLog::GetDebugLevel("",IsA()->GetName()))>=10){
	par->Print("");
	cov->Print("");
    }

     if(!fPar){fPar = new Double_t[fNPar];}
     for(int ip=0; ip<fNPar; ip++) fPar[ip] = (*par)[ip];
     if(!cov||!storeCov) return;
     if(!fCov){fCov = new Double_t[fNCov];}
     for(int ip(0), np(0); ip<fNPar; ip++)
	 for(int jp=ip; jp<fNPar; jp++) fCov[np++] = (*cov)(ip,jp);

     if((AliLog::GetDebugLevel("",IsA()->GetName()))>10){this->Print("a");}
}

//_________________________________________________________________
Bool_t AliTRDTKDInterpolator::AliTRDTKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t &error,TRDTKDMode mod) const
{
    // Recalculate the PDF for one node from the results of interpolation (parameters and covariance matrix)

    result =0.; error = 1.;

    if(mod==kNodeVal){
	error=fVal[1];
	result=fVal[0];
	return kTRUE;
    }

    if(!fPar){
	AliDebug(1,"Requested Interpolation Parameters don't exist");
	return kFALSE;
    }

    Double_t fdfdp[fNDim];
    Int_t 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
    for(int i=0; i<fNPar; i++) result += fdfdp[i]*fPar[i];

    if(!fCov){
	AliDebug(3,"Interpolation Error cannot be estimated, Covariance Parameters don't exist");
	return kTRUE;
    }
    // calculate error
    error=0;

    for(int i(0), n(0); i<fNPar; i++){
	error += fdfdp[i]*fdfdp[i]*fCov[n++];
	for(int j(i+1); j<fNPar; j++) error += 2.*fdfdp[i]*fdfdp[j]*fCov[n++];
    }
    if(error>0)error = TMath::Sqrt(error);
    else{error=100;}

    if(mod==kMinError){
	if(error<fVal[1]){
	    return kTRUE;
	}
	if(error==1)AliWarning("Covariance not available, always choose bin content");
	error=fVal[1];
	result=fVal[0];
	return kTRUE;
    }

    // Boundary condition
    if(result<0){
        AliDebug(2,"Using Node Value to ensure Boundary Condition");
	result=fVal[0];
	error=fVal[1];
    }

    AliDebug(1,Form("Cook PDF Result: %e Error: %e",result,error));

    return kTRUE;
}

//_____________________________________________________________________
Bool_t AliTRDTKDInterpolator::AliTRDTKDNodeInfo::Has(const Float_t *p) const
{
  Int_t n(0);
  for(int id=0; id<fNDim; id++)
      if(p[id]>=fBounds[2*id] && p[id]<fBounds[2*id+1]) n++;
  if(n==fNDim) return kTRUE;
  return kFALSE;
}
Commit	Line	Data
db0e2c5f	1	#include "AliTRDTKDInterpolator.h"
	2
	3	#include "TClonesArray.h"
	4	#include "TLinearFitter.h"
	5	#include "TMath.h"
	6	#include "TRandom.h"
	7
f2762b1c	8	#include "AliLog.h"
f2762b1c	9
db0e2c5f	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),
bd58d4b9	24	fPDFMode(kInterpolation),
bd58d4b9	25	fStoreCov(kFALSE)
db0e2c5f	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),
bd58d4b9	39	fPDFMode(kInterpolation),
bd58d4b9	40	fStoreCov(kFALSE)
db0e2c5f	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),
bd58d4b9	62	fPDFMode(ref.fPDFMode),
bd58d4b9	63	fStoreCov(ref.fStoreCov)
db0e2c5f	64	{
db0e2c5f	65	// Copy constructor
f2762b1c	66	this->Print("");
db0e2c5f	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
f2762b1c	80	this->Print("");
f2762b1c	81
db0e2c5f	82	return *this;
	83	}
	84
	85	//_________________________________________________________________
	86	Bool_t AliTRDTKDInterpolator::Build()
	87	{
	88	TKDTreeIF::Build();
	89	if(!fBoundaries) MakeBoundaries();
	90
	91	// allocate interpolation nodes
f2762b1c	92	fNDataNodes = fNPoints/fBucketSize + ((fNPoints%fBucketSize)?1:0);
db0e2c5f	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),2fNDimsizeof(Float_t));
	108	node->fVal[0] = Float_t(fBucketSize)/fNPoints;
	109	for(int idim=0; idim<fNDim; idim++) node->fVal[0] /= (node->fBounds[2idim+1] - node->fBounds[2idim]);
	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),2fNDimsizeof(Float_t));
	129
	130	for(int idim=0; idim<fNDim; idim++){
	131	Float_t dx = ftnode->fBounds[2idim+1]-ftnode->fBounds[2idim];
	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[2idim+1] - ftnode->fBounds[2idim]);
	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	{
f2762b1c	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
db0e2c5f	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){
f2762b1c	175	AliError("Can not retrieve node for data point");
db0e2c5f	176	result = 0.;
	177	error = 1.E10;
	178	return kFALSE;
	179	}
	180	AliTRDTKDNodeInfo *node =GetNodeInfo(nodeIndex);
	181
f2762b1c	182	if((AliLog::GetDebugLevel("",IsA()->GetName()))>0){
	183	printf("Node Info: \n");
	184	node->Print("a");
db0e2c5f	185	}
f2762b1c	186
f2762b1c	187	return node->CookPDF(point, result, error,fPDFMode);
db0e2c5f	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;
f2762b1c	204	if(GetNodeInfo(inode)->Has(p)){
	205	AliDebug(2,Form("Find Node %d",inode));
	206	return inode;
	207	}
db0e2c5f	208
	209	// Search extra nodes
	210
	211	for(inode=fNDataNodes;inode<GetNTNodes();inode++){
f2762b1c	212	if(GetNodeInfo(inode)->Has(p)){
	213	AliDebug(2,Form("Find Extra Node %d",inode));
	214	return inode;
	215	}
db0e2c5f	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	}
f2762b1c	231	AliDebug(2,Form("Find Nearest Neighbor Node %d",inode));
db0e2c5f	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[2idim]<range[0]) range[0] = node->fBounds[2idim];
	261	if(node->fBounds[2idim+1]>range[1]) range[1] = node->fBounds[2idim+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[2xdim],node->fBounds[2ydim],node->fBounds[2xdim+1],node->fBounds[2ydim+1]);
	282	h2->SetBinContent(inode+1, node->fVal[0]);
	283	}
	284	return h2;
	285	}
	286
	287	//_________________________________________________________________
	288	void AliTRDTKDInterpolator::BuildInterpolation()
	289	{
f2762b1c	290	AliInfo("Build Interpolation");
db0e2c5f	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. - ddd);
	328	w = w0w0w0;
	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
f2762b1c	354	AliDebug(2,Form("Calculate Interpolation for Node %d",nodeIndex));
db0e2c5f	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
bd58d4b9	364	node->Store(&par,&cov,fStoreCov);
db0e2c5f	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[2vdim+iter]==treebounds[2vdim+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[2vdim+iter]=2treebounds[2*vdim+iter];
	400	newnode->fBounds[2vdim+!iter]=node->fBounds[2vdim+iter];
	401	for(Int_t idim=0;idim<fNDim;idim++){
	402	if(idim==vdim)continue;
	403	newnode->fBounds[2idim]=node->fBounds[2idim];
	404	newnode->fBounds[2idim+1]=node->fBounds[2idim+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[2idim+1] - newnode->fBounds[2idim]);
	410	newnode->fData[idim]=0.5(newnode->fBounds[2idim+1] + newnode->fBounds[2*idim]);
	411	}
	412	nnew++;
	413	}
	414	}
	415	}
	416	}
f2762b1c	417	AliInfo(Form("%d Boundary Nodes Added \n",nnew));
db0e2c5f	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))
f2762b1c	426	,fNCov(Int_t((fNPar+1)*fNPar/2))
db0e2c5f	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, 2fNDimsizeof(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, 2fNDimsizeof(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[2idim], fBounds[2idim+1]);
521	printf("]\n");
522	}
523	if(strcmp(opt, "a")!=0) return;
524
525	if(fPar){
526	printf("Fit parameters : \n");
f2762b1c	527	for(int ip=0; ip<fNPar; ip++) printf("p%d[%e] ", ip, fPar[ip]);
db0e2c5f	528	printf("\n");
	529	}
	530	if(!fCov) return;
	531	for(int ip(0), n(0); ip<fNPar; ip++){
f2762b1c	532	for(int jp(ip); jp<fNPar; jp++) printf("c(%d %d)[%e] ", ip, jp, fCov[n++]);
db0e2c5f	533	printf("\n");
	534	}
	535	}
	536
	537	//_________________________________________________________________
bd58d4b9	538	void AliTRDTKDInterpolator::AliTRDTKDNodeInfo::Store(TVectorD const par, TMatrixD const cov,Bool_t storeCov)
db0e2c5f	539	{
f2762b1c	540	// Store the parameters and the covariance in the node
db0e2c5f	541
f2762b1c	542	AliDebug(2,"Store Node Interpolation Parameters");
db0e2c5f	543
f2762b1c	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];
bd58d4b9	551	if(!cov\|\|!storeCov) return;
f2762b1c	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");}
db0e2c5f	557	}
	558
	559	//_________________________________________________________________
f2762b1c	560	Bool_t AliTRDTKDInterpolator::AliTRDTKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t &error,TRDTKDMode mod) const
db0e2c5f	561	{
	562	// Recalculate the PDF for one node from the results of interpolation (parameters and covariance matrix)
	563
	564	result =0.; error = 1.;
f2762b1c	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	}
db0e2c5f	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
f2762b1c	588	if(!fCov){
	589	AliDebug(3,"Interpolation Error cannot be estimated, Covariance Parameters don't exist");
	590	return kTRUE;
	591	}
db0e2c5f	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
f2762b1c	602	if(mod==kMinError){
	603	if(error<fVal[1]){
	604	return kTRUE;
	605	}
bd58d4b9	606	if(error==1)AliWarning("Covariance not available, always choose bin content");
f2762b1c	607	error=fVal[1];
	608	result=fVal[0];
	609	return kTRUE;
	610	}
	611
db0e2c5f	612	// Boundary condition
db0e2c5f	613	if(result<0){
f2762b1c	614	AliDebug(2,"Using Node Value to ensure Boundary Condition");
db0e2c5f	615	result=fVal[0];
	616	error=fVal[1];
	617	}
	618
f2762b1c	619	AliDebug(1,Form("Cook PDF Result: %e Error: %e",result,error));
f2762b1c	620
db0e2c5f	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[2id] && p[id]<fBounds[2id+1]) n++;
	630	if(n==fNDim) return kTRUE;
	631	return kFALSE;
	632	}
	633