- bug fixes

- more user verbosity
- tune for usage with TRD PID

diff --git a/STAT/TKDInterpolatorBase.cxx b/STAT/TKDInterpolatorBase.cxx
index 3ae38700b0ac490367e95822150faf1287133b15..2afcf08b24df288abfd98504d7fa7e871773ba16 100644 (file)
--- a/STAT/TKDInterpolatorBase.cxx
+++ b/STAT/TKDInterpolatorBase.cxx
@@ -2,6 +2,7 @@
 #include "TKDNodeInfo.h"
 #include "TKDTree.h"
 
+#include "TRandom.h"
 #include "TClonesArray.h"
 #include "TLinearFitter.h"
 #include "TTree.h"
@@ -54,6 +55,10 @@ void TKDInterpolatorBase::Build(Int_t n)
 
   if(fTNodes) delete fTNodes;
   fNTNodes = n;
+  // check granularity
+  if(Int_t((1.+fAlpha)*fLambda) > fNTNodes){
+    Warning("TKDInterpolatorBase::Build()", Form("Minimum number of points [%d] needed for interpolation exceeds number of evaluation points [%d]. Please increase granularity.", Int_t((1.+fAlpha)*fLambda), fNTNodes));
+  }
   fTNodes = new TClonesArray("TKDNodeInfo", fNTNodes);
   for(int in=0; in<fNTNodes; in++) new ((*fTNodes)[in]) TKDNodeInfo(fNSize);
 }
@@ -125,6 +130,7 @@ Double_t TKDInterpolatorBase::Eval(const Double_t *point, Double_t &result, Doub
   for(int idim=0; idim<fNSize; idim++) pointF[idim] = (Float_t)point[idim];
   Int_t nodeIndex = GetNodeIndex(pointF);
   if(nodeIndex<0){
+    Error("TKDInterpolatorBase::Eval()", "Can not retrive node for data point.");
     result = 0.;
     error = 1.E10;
     return 0.;
@@ -140,8 +146,9 @@ Double_t TKDInterpolatorBase::Eval(const Double_t *point, Double_t &result, Doub
       fRefPoints[id] = new Float_t[fNTNodes];
       for(int in=0; in<fNTNodes; in++) fRefPoints[id][in] = ((TKDNodeInfo*)(*fTNodes)[in])->Data()[id];
     }
-    fKDhelper = new TKDTreeIF(fNTNodes, fNSize, 30, fRefPoints);
-    fKDhelper->MakeBoundaries();
+    fKDhelper = new TKDTreeIF(fNTNodes, fNSize, kNhelper, fRefPoints);
+    fKDhelper->Build();
+    fKDhelper->MakeBoundariesExact();
   }
   if(!fFitter) fFitter = new TLinearFitter(fLambda, Form("hyp%d", fLambda-1));
   
@@ -152,17 +159,32 @@ Double_t TKDInterpolatorBase::Eval(const Double_t *point, Double_t &result, Doub
   // prepare workers
 
   Int_t ipar,    // local looping variable
-        npoints = Int_t((1.+fAlpha)*fLambda); // number of data points used for interpolation
-  Int_t *index = new Int_t[2*npoints];  // indexes of NN 
-  Float_t *dist = new Float_t[2*npoints], // distances of NN
+        npoints_new = Int_t((1.+fAlpha)*fLambda),
+        npoints(0); // number of data points used for interpolation
+  Int_t *index = new Int_t[2*npoints_new];  // indexes of NN 
+  Float_t *dist = new Float_t[2*npoints_new], // distances of NN
           d,     // NN normalized distance
           w0,    // work
           w;     // tri-cubic weight function
 
+  Bool_t kDOWN = kFALSE;
   do{
+    if(npoints == npoints_new){
+      Error("TKDInterpolatorBase::Eval()", Form("Interpolation failed and number of interpolation points (%d) Can not be increased further.", npoints));
+      result = 0.;
+      error = 1.E10;
+      return 0.;
+    } else npoints = npoints_new;
+    if(npoints > fNTNodes){
+      Warning("TKDInterpolatorBase::Eval()", Form("The number of interpolation points requested (%d) exceeds number of PDF values (%d). Downscale.", npoints, fNTNodes));
+      npoints = fNTNodes;
+      kDOWN = kTRUE;
+    }
+
     // find nearest neighbors
     for(int idim=0; idim<fNSize; idim++) pointF[idim] = (Float_t)point[idim];
     fKDhelper->FindNearestNeighbors(pointF, npoints+1, index, dist);
+
     // add points to fitter
     fFitter->ClearPoints();
     TKDNodeInfo *tnode = 0x0;
@@ -192,12 +214,12 @@ Double_t TKDInterpolatorBase::Eval(const Double_t *point, Double_t &result, Doub
         w0 = (1. - d*d*d); w = w0*w0*w0;
       } else w = 1.;
       
-//                     printf("x[");
-//                     for(int idim=0; idim<fLambda-1; idim++) printf("%f ", fBuffer[idim]);
-//                     printf("]  v[%f +- %f] (%f, %f)\n", tnode->Val()[0], tnode->Val()[1]/w, tnode->Val()[1], w);
+//       printf("%2d x[", index[in]);
+//       for(int idim=0; idim<fLambda-1; idim++) printf("%f ", fBuffer[idim]);
+//       printf("]  v[%f +- %f] (%f, %f)\n", tnode->Val()[0], tnode->Val()[1]/w, tnode->Val()[1], w);
       fFitter->AddPoint(fBuffer, tnode->Val()[0], tnode->Val()[1]/w);
     }
-    npoints += 4;
+    npoints_new = npoints+ (kDOWN ? 0 : kdN);
   } while(fFitter->Eval());
   delete [] index;
   delete [] dist;
@@ -256,7 +278,7 @@ void TKDInterpolatorBase::DrawBins(UInt_t ax1, UInt_t ax2, Float_t ax1min, Float
   for(int inode = 0; inode < fNTNodes; inode++){
     box = &boxArray[inode];
     box->SetFillStyle(3002);
-    box->SetFillColor(50+inode/*Int_t(gRandom->Uniform()*50.)*/);
+    box->SetFillColor(50+Int_t(gRandom->Uniform()*50.));
     
     bounds = &(((TKDNodeInfo*)(*fTNodes)[inode])->Data()[fNSize]);
     box->DrawBox(bounds[2*ax1]+kBorder, bounds[2*ax2]+kBorder, bounds[2*ax1+1]-kBorder, bounds[2*ax2+1]-kBorder);
@@ -275,6 +297,25 @@ void TKDInterpolatorBase::DrawBins(UInt_t ax1, UInt_t ax2, Float_t ax1min, Float
   return;
 }
 
+//_________________________________________________________________
+void TKDInterpolatorBase::SetAlpha(Float_t a)
+{
+  if(a<0.5){ 
+    Warning("TKDInterpolatorBase::SetAlpha()", "The scale parameter has to be larger than 0.5");
+    fAlpha = 0.5;
+    return;
+  }
+  // check value
+  if(Int_t((a+1.)*fLambda) > fNTNodes){
+    fAlpha = TMath::Max(0.5, Float_t(fNTNodes)/fLambda-1.);
+    Warning("TKDInterpolatorBase::SetAlpha()", Form("Interpolation neighborhood  exceeds number of evaluation points. Downscale alpha to %f", fAlpha));
+    printf("n[%d] nodes[%d]\n", Int_t((fAlpha+1.)*fLambda), fNTNodes);
+    return;
+  }
+  fAlpha = a;
+  return;
+}
+
 //__________________________________________________________________
 void TKDInterpolatorBase::SetInterpolationMethod(Bool_t on)
 {

diff --git a/STAT/TKDInterpolatorBase.h b/STAT/TKDInterpolatorBase.h
index 2674cf2c72feaee1a0de48da035ce779dd8d36db..4bd228e79b3a162dc67df827d5c63374c6a5743a 100644 (file)
--- a/STAT/TKDInterpolatorBase.h
+++ b/STAT/TKDInterpolatorBase.h
@@ -26,6 +26,11 @@ class TKDNodeInfo;
 class TKDInterpolatorBase
 {
 public:        
+  enum EKDInterpolatorBase {
+    kdN = 4       // increase in the number of PDF if fit failled
+   ,kNhelper = 30 // bucket size in helper kdTree
+  };
+
   TKDInterpolatorBase(Int_t size = 0);
   virtual ~TKDInterpolatorBase();
 
@@ -41,7 +46,7 @@ public:
   Bool_t     GetWeights() const {return fStatus&4;}
 
   void       DrawBins(UInt_t ax1 = 0, UInt_t ax2 = 1, Float_t ax1min=-1., Float_t ax1max=1., Float_t ax2min=-1., Float_t ax2max=1.);
-  void       SetAlpha(Float_t a){if(a>0.) fAlpha = a;}
+  void       SetAlpha(Float_t a);
   void       SetInterpolationMethod(Bool_t on = kTRUE);
   void       SetStore(Bool_t on = kTRUE);
   void       SetWeights(Bool_t on = kTRUE);
@@ -62,7 +67,7 @@ protected:
   UChar_t       fStatus;      // status of the interpolator
   UChar_t       fLambda;      // number of parameters in polynom
   Short_t                  fDepth;       //! depth of the KD Tree structure used
-  Float_t       fAlpha;       // alpha parameter
+  Float_t       fAlpha;       // parameter controlling the size of the region to interpolate n = (1+alpha)*lambda
   Float_t       **fRefPoints; //! temporary storage of COG data
   Double_t         *fBuffer;     //! working space [2*fLambda]
   TKDTree<Int_t, Float_t> *fKDhelper;      //! kNN finder

diff --git a/STAT/TKDPDF.cxx b/STAT/TKDPDF.cxx
index 65d3ad268c658e8be02125ffaa9305520cfe866d..9967bd790b87c9e64b11cb716fc9e970bc3263f5 100644 (file)
--- a/STAT/TKDPDF.cxx
+++ b/STAT/TKDPDF.cxx
@@ -71,7 +71,6 @@ TKDPDF::TKDPDF(TTree *t, const Char_t *var, const Char_t *cut, UInt_t bsize, Lon
     v = t->GetV1();
     for(int ip=0; ip<fNPoints; ip++) fData[idim][ip] = (Float_t)v[ip];
   }
-  TKDTreeIF::Build();
   Build();
 }
 
@@ -87,6 +86,7 @@ void TKDPDF::Build(Int_t)
 //  - estimation points 
 //  - corresponding PDF values
 
+  TKDTreeIF::Build();
   fNTNodes = fNPoints/fBucketSize + ((fNPoints%fBucketSize)?1:0);/*TKDTreeIF::GetNTNodes();*/
   if(!fBoundaries) MakeBoundaries();
   fLambda = 1 + fNDim + (fNDim*(fNDim+1)>>1);
@@ -122,7 +122,14 @@ void TKDPDF::Build(Int_t)
   node = (TKDNodeInfo*)(*fTNodes)[inode];
   node->Val()[0] =  Float_t(counts)/fNPoints;
   bounds = GetBoundary(tnode);
-  for(int idim=0; idim<fNDim; idim++) node->Val()[0] /= (bounds[2*idim+1] - bounds[2*idim]);
+  for(int idim=0; idim<fNDim; idim++){ 
+    Float_t dx = bounds[2*idim+1]-bounds[2*idim];
+    if(dx < 1.e-30){
+      Warning("TKDPDF::Build()", Form("Terminal bucket index[%d] too narrow on the %d dimension.", inode, idim));
+      continue;
+    }
+    node->Val()[0] /= (bounds[2*idim+1] - bounds[2*idim]);
+  }
   node->Val()[1] =  node->Val()[0]/TMath::Sqrt(float(counts));
 
   // loop points in this terminal node

diff --git a/STAT/TKDPDF.h b/STAT/TKDPDF.h
index 8ec6a384ce0a47ed72931bf87ea3413cb5e36b2d..257098e980a9421defd5aa4a45398a7095d34193 100644 (file)
--- a/STAT/TKDPDF.h
+++ b/STAT/TKDPDF.h
@@ -20,39 +20,39 @@ class TLinearFitter;
 class TKDPDF : public TKDTreeIF, public TKDInterpolatorBase
 {
 public:
-       TKDPDF();
-       TKDPDF(TTree *t, const Char_t *var, const Char_t *cut = 0, UInt_t bsize = 100, Long64_t nentries = 1000000000, Long64_t firstentry = 0);
-       TKDPDF(Int_t npoints, Int_t ndim, UInt_t bsize, Float_t **data);
-       ~TKDPDF();
+  TKDPDF();
+  TKDPDF(TTree *t, const Char_t *var, const Char_t *cut = 0, UInt_t bsize = 100, Long64_t nentries = 1000000000, Long64_t firstentry = 0);
+  TKDPDF(Int_t npoints, Int_t ndim, UInt_t bsize, Float_t **data);
+  ~TKDPDF();
 
-       inline Bool_t     GetDataPoint(Int_t n, Float_t *p) const;
-       inline Int_t      GetNodeIndex(const Float_t *p);
-                                       void       DrawNode(Int_t tnode, UInt_t ax1=0, UInt_t ax2=1);
+  inline Bool_t GetDataPoint(Int_t n, Float_t *p) const;
+  inline Int_t  GetNodeIndex(const Float_t *p);
+  void          DrawNode(Int_t tnode, UInt_t ax1=0, UInt_t ax2=1);
 
 private:
-       TKDPDF(const TKDPDF &);
-       TKDPDF& operator=(const TKDPDF &);
-               void       Build(Int_t ndim = 0);
+  TKDPDF(const TKDPDF &);
+  TKDPDF& operator=(const TKDPDF &);
+  void     Build(Int_t ndim = 0);
 
-                                       
-       ClassDef(TKDPDF, 1)   // data interpolator based on KD tree
+          
+  ClassDef(TKDPDF, 1)   // data interpolator based on KD tree
 };
 
 
 //__________________________________________________________________
 Bool_t TKDPDF::GetDataPoint(Int_t n, Float_t *p) const
 {
-       if(n < 0 || n >= fNPoints) return kFALSE;
-       if(!fData) return kFALSE;
-               
-       for(int i=0; i<fNDim; i++) p[i] = fData[i][n];
-       return kTRUE;
+  if(n < 0 || n >= fNPoints) return kFALSE;
+  if(!fData) return kFALSE;
+    
+  for(int i=0; i<fNDim; i++) p[i] = fData[i][n];
+  return kTRUE;
 }
 
 //__________________________________________________________________
 Int_t  TKDPDF::GetNodeIndex(const Float_t *p)
 {
-       return FindNode(p) - fNNodes;
+  return FindNode(p) - fNNodes;
 }
 
 #endif