- extend user interface
authorabercuci <abercuci@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 21 Oct 2009 06:34:52 +0000 (06:34 +0000)
committerabercuci <abercuci@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 21 Oct 2009 06:34:52 +0000 (06:34 +0000)
- protect code against errors
- documentation

STAT/TKDInterpolator.cxx
STAT/TKDInterpolator.h
STAT/TKDInterpolatorBase.cxx
STAT/TKDInterpolatorBase.h
STAT/TKDNodeInfo.cxx
STAT/TKDNodeInfo.h

index 2abe0d1..8c4ddd5 100644 (file)
@@ -1,6 +1,7 @@
 #include "TKDInterpolator.h"
 #include "TKDNodeInfo.h"
 
+#include "TError.h"
 #include "TClonesArray.h"
 
 ClassImp(TKDInterpolator)
@@ -9,7 +10,7 @@ ClassImp(TKDInterpolator)
 
 //_________________________________________________________________
 TKDInterpolator::TKDInterpolator() :
-       TKDInterpolatorBase()
+  TKDInterpolatorBase()
 {
 // Default constructor. To be used with care since in this case building
 // of data structure is completly left to the user responsability.
@@ -17,11 +18,11 @@ TKDInterpolator::TKDInterpolator() :
 
 //_________________________________________________________________
 TKDInterpolator::TKDInterpolator(Int_t ndim, Int_t npoints) :
-       TKDInterpolatorBase(ndim)
+  TKDInterpolatorBase(ndim)
 {
 // Wrapper constructor for the TKDTree.
 
-       if(npoints) TKDInterpolatorBase::Build(npoints);
+  if(npoints) TKDInterpolatorBase::Build(npoints);
 }
 
 
@@ -33,51 +34,52 @@ TKDInterpolator::~TKDInterpolator()
 //_________________________________________________________________
 void TKDInterpolator::AddNode(const TKDNodeInfo &node)
 {
-       if(!fTNodes){
-               printf("W - TKDInterpolator::SetNode() : Node array not defined.\n");
-               return;
-       }
+  if(!fTNodes){
+    Warning("TKDInterpolator::SetNode()", "Node array not defined.");
+    return;
+  }
 
-       new((*fTNodes)[fNTNodes++]) TKDNodeInfo(node);
+  new((*fTNodes)[fNTNodes++]) TKDNodeInfo(node);
 }
 
 //_________________________________________________________________
 void TKDInterpolator::Build(Int_t npoints, Int_t ndim)
 {
-       fNSize = ndim;
-       TKDInterpolatorBase::Build(npoints);
+  fNSize = ndim;
+  TKDInterpolatorBase::Build(npoints);
 }
 
 //_________________________________________________________________
 Int_t TKDInterpolator::GetNodeIndex(const Float_t *p)
 {
-/*     printf("TKDInterpolator::GetNodeIndex() ...\n");
-       printf("Looking for p[");
-       for(int i=0; i<fNSize; i++) printf("%f ", p[i]);
-       printf("] ...\n");*/
-       TKDNodeInfo *node;
-       for(int inode=0; inode<fNTNodes; inode++){
-               node = (TKDNodeInfo*)(*fTNodes)[inode];
-               //node->Print();
-               if(node->Has(p)) return inode;
-       }
-       return -1;
+//     printf("TKDInterpolator::GetNodeIndex() ...\n");
+//   printf("Looking for p[");
+//   for(int i=0; i<fNSize; i++) printf("%f ", p[i]);
+//   printf("] ...\n");
+
+  for(Int_t i=fNTNodes; i--;)
+    if(((TKDNodeInfo*)(*fTNodes)[i])->Has(p)) return i;
+  
+  printf("Point p[");
+  for(int i=0; i<fNSize; i++) printf("%f ", p[i]);
+  printf("] outside range.\n");
+  return -1;
 }
 
 
 //_________________________________________________________________
 Bool_t TKDInterpolator::SetNode(Int_t inode, const TKDNodeInfo &ref)
 {
-       if(!fTNodes){
-               printf("W - TKDInterpolator::SetNode() : Node array not defined.\n");
-               return kFALSE;
-       }
-       if(inode >= fNTNodes){
-               printf("W - TKDInterpolator::SetNode() : Node array defined up to %d.\n", fNTNodes);
-               return kFALSE;
-       }
-       TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[inode];
-       (*node) = ref;
-       return kTRUE;
+  if(!fTNodes){
+    Warning("TKDInterpolator::SetNode()", "Node array not defined.");
+    return kFALSE;
+  }
+  if(inode >= fNTNodes){
+    Warning("TKDInterpolator::SetNode()", Form("Node array defined up to %d.", fNTNodes));
+    return kFALSE;
+  }
+  TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[inode];
+  (*node) = ref;
+  return kTRUE;
 }
 
index 8970a9d..1c3f6a9 100644 (file)
@@ -8,22 +8,22 @@
 class TKDInterpolator : public TKDInterpolatorBase
 {
 public:
-       TKDInterpolator();
-       TKDInterpolator(Int_t ndim, Int_t npoints=0);
-       ~TKDInterpolator();
-       void       AddNode(const TKDNodeInfo &ref);
-       void       Build(Int_t ndim) {TKDInterpolatorBase::Build(ndim);}
-       void       Build(Int_t npoints, Int_t ndim);
-       Int_t      GetNodeIndex(const Float_t *p);
-       Bool_t     SetNode(Int_t i, const TKDNodeInfo &ref);
+  TKDInterpolator();
+  TKDInterpolator(Int_t ndim, Int_t npoints=0);
+  ~TKDInterpolator();
+  void       AddNode(const TKDNodeInfo &ref);
+  void       Build(Int_t ndim) {TKDInterpolatorBase::Build(ndim);}
+  void       Build(Int_t npoints, Int_t ndim);
+  Int_t      GetNodeIndex(const Float_t *p);
+  Bool_t     SetNode(Int_t i, const TKDNodeInfo &ref);
 
 private:
-       TKDInterpolator(const TKDInterpolator &);
-       TKDInterpolator& operator=(const TKDInterpolator &);    
+  TKDInterpolator(const TKDInterpolator &);
+  TKDInterpolator& operator=(const TKDInterpolator &); 
 
 private:
-       
-       ClassDef(TKDInterpolator, 1)   // LOWESS data interpolator
+  
+  ClassDef(TKDInterpolator, 1)   // LOWESS data interpolator
 };
 
 
index df0c193..074c434 100644 (file)
@@ -33,7 +33,7 @@ TKDInterpolatorBase::TKDInterpolatorBase(Int_t dim) :
   ,fNTNodes(0)
   ,fTNodes(0x0)
   ,fTNodesDraw(0x0)
-  ,fStatus(4)
+  ,fStatus(0)
   ,fLambda(1 + dim + (dim*(dim+1)>>1))
   ,fDepth(-1)
   ,fAlpha(.5)
@@ -44,6 +44,7 @@ TKDInterpolatorBase::TKDInterpolatorBase(Int_t dim) :
 {
 // Default constructor. To be used with care since in this case building
 // of data structure is completly left to the user responsability.
+  UseWeights();
 }
 
 //_________________________________________________________________
@@ -102,18 +103,39 @@ TKDNodeInfo* TKDInterpolatorBase::GetNodeInfo(Int_t inode) const
   return (TKDNodeInfo*)(*fTNodes)[inode];
 }
 
+//_________________________________________________________________
+Bool_t TKDInterpolatorBase::GetRange(Int_t ax, Float_t &min, Float_t &max) const
+{
+  if(!fTNodes) return kFALSE;
+  Int_t ndim = ((TKDNodeInfo*)(*fTNodes)[0])->GetDimension();
+  if(ax<0 || ax>=ndim){
+    min=0.; max=0.;
+    return kFALSE;
+  }
+  min=1.e10; max=-1.e10;
+  Float_t axmin, axmax;
+  for(Int_t in=fNTNodes; in--; ){ 
+    TKDNodeInfo *node = (TKDNodeInfo*)((*fTNodes)[in]);
+    node->GetBoundary(ax, axmin, axmax);
+    if(axmin<min) min = axmin;
+    if(axmax>max) max = axmax;
+  }
+  
+  return kTRUE;
+}
 
 //__________________________________________________________________
-void TKDInterpolatorBase::GetStatus()
+void TKDInterpolatorBase::GetStatus(Option_t *opt)
 {
 // Prints the status of the interpolator
 
-  printf("Interpolator Status :\n");
+  printf("Interpolator Status[%d] :\n", fStatus);
   printf("  Dim    : %d [%d]\n", fNSize, fLambda);
-  printf("  Method : %s\n", fStatus&1 ? "INT" : "COG");
-  printf("  Store  : %s\n", fStatus&2 ? "YES" : "NO");
-  printf("  Weights: %s\n", fStatus&4 ? "YES" : "NO");
+  printf("  Method : %s\n", UseCOG() ? "COG" : "INT");
+  printf("  Store  : %s\n", HasStore() ? "YES" : "NO");
+  printf("  Weights: %s\n", UseWeights() ? "YES" : "NO");
   
+  if(strcmp(opt, "all") != 0 ) return;
   printf("fNTNodes %d\n", fNTNodes);        //Number of evaluation data points
   for(int i=0; i<fNTNodes; i++){
     TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[i]; 
@@ -141,7 +163,7 @@ Double_t TKDInterpolatorBase::Eval(const Double_t *point, Double_t &result, Doub
     return 0.;
   }
   TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[nodeIndex];
-  if((fStatus&1) && node->Cov() && !force) return node->CookPDF(point, result, error);
+  if(node->Cov() && !force) return node->CookPDF(point, result, error);
 
   // Allocate memory
   if(!fBuffer) fBuffer = new Double_t[2*fLambda];
@@ -151,6 +173,7 @@ 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];
     }
+    Info("TKDInterpolatorBase::Eval()", Form("Build TKDTree(%d, %d, %d)", fNTNodes, fNSize, kNhelper));
     fKDhelper = new TKDTreeIF(fNTNodes, fNSize, kNhelper, fRefPoints);
     fKDhelper->Build();
     fKDhelper->MakeBoundariesExact();
@@ -174,6 +197,9 @@ Double_t TKDInterpolatorBase::Eval(const Double_t *point, Double_t &result, Doub
 
   Bool_t kDOWN = kFALSE;
   do{
+    if(npoints){
+      Info("TKDInterpolatorBase::Eval()", Form("Interpolation failed. Trying to increase the number of interpolation points from %d to %d.", npoints, npoints_new));
+    }
     if(npoints == npoints_new){
       Error("TKDInterpolatorBase::Eval()", Form("Interpolation failed and number of interpolation points (%d) Can not be increased further.", npoints));
       result = 0.;
@@ -196,7 +222,14 @@ Double_t TKDInterpolatorBase::Eval(const Double_t *point, Double_t &result, Doub
     for(int in=0; in<npoints; in++){
       tnode = (TKDNodeInfo*)(*fTNodes)[index[in]];
       //tnode->Print();
-      if(fStatus&1){ // INT
+      if(UseCOG()){ // COG
+        Float_t *p = &(tnode->Data()[0]);
+        ipar = 0;
+        for(int idim=0; idim<fNSize; idim++){
+          fBuffer[ipar++] = p[idim];
+          for(int jdim=idim; jdim<fNSize; jdim++) fBuffer[ipar++] = p[idim]*p[jdim];
+        }
+      } else { // INT
         Float_t *bounds = &(tnode->Data()[fNSize]);
         ipar = 0;
         for(int idim=0; idim<fNSize; idim++){
@@ -204,24 +237,18 @@ Double_t TKDInterpolatorBase::Eval(const Double_t *point, Double_t &result, Doub
           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<fNSize; jdim++) fBuffer[ipar++] = (bounds[2*idim] + bounds[2*idim+1]) * (bounds[2*jdim] + bounds[2*jdim+1]) * .25;
         }
-      } else { // COG
-        Float_t *p = &(tnode->Data()[0]);
-        ipar = 0;
-        for(int idim=0; idim<fNSize; idim++){
-          fBuffer[ipar++] = p[idim];
-          for(int jdim=idim; jdim<fNSize; jdim++) fBuffer[ipar++] = p[idim]*p[jdim];
-        }
       }
 
       // calculate tri-cubic weighting function
-      if(fStatus&4){
-        d = dist[in]/ dist[npoints];
+      if(UseWeights()){
+        d = dist[in]/dist[npoints];
         w0 = (1. - d*d*d); w = w0*w0*w0;
+        if(w<1.e-30) continue;
       } else w = 1.;
       
-//       printf("%2d x[", index[in]);
+//       printf("%2d d[%f] w[%f] x[", index[in], d, w);
 //       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("]\n");  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_new = npoints+ (kDOWN ? 0 : kdN);
@@ -237,7 +264,7 @@ Double_t TKDInterpolatorBase::Eval(const Double_t *point, Double_t &result, Doub
   Double_t chi2 = fFitter->GetChisquare()/(npoints - 4 - fLambda);
 
   // store results
-  if(fStatus&2 && fStatus&1) node->Store(par, cov);
+  if(HasStore()) node->Store(par, cov);
     
   // Build df/dpi|x values
   Double_t *fdfdp = &fBuffer[fLambda];
@@ -260,36 +287,35 @@ Double_t TKDInterpolatorBase::Eval(const Double_t *point, Double_t &result, Doub
 }
 
 //_________________________________________________________________
-void TKDInterpolatorBase::DrawBins(UInt_t ax1, UInt_t ax2, Float_t ax1min, Float_t ax1max, Float_t ax2min, Float_t ax2max)
+void TKDInterpolatorBase::DrawProjection(UInt_t ax1, UInt_t ax2)
 {
 // Draw nodes structure projected on plane "ax1:ax2". The parameter
 // "depth" specifies the bucket size per node. If depth == -1 draw only
 // terminal nodes and evaluation points (default -1 i.e. bucket size per node equal bucket size specified by the user)
 //
-// Observation:
-// This function creates the nodes (TBox) array for the specified depth
-// but don't delete it. Abusing this function may cause memory leaks !
-
-
   
+  Float_t ax1min, ax1max, ax2min, ax2max;
+  GetRange(ax1, ax1min, ax1max);
+  GetRange(ax2, ax2min, ax2max);
   TH2 *h2 = 0x0;
-  if(!(h2 = (TH2S*)gROOT->FindObject("hKDnodes"))){ 
+  if(!(h2 = (TH2S*)gROOT->FindObject("hKDnodes"))){
     h2 = new TH2S("hKDnodes", "", 100, ax1min, ax1max, 100, ax2min, ax2max);
-    h2->GetXaxis()->SetTitle(Form("x_{%d}", ax1));
-    h2->GetYaxis()->SetTitle(Form("x_{%d}", ax2));
   }
+  h2->GetXaxis()->SetRangeUser(ax1min, ax1max);
+  h2->GetXaxis()->SetTitle(Form("x_{%d}", ax1));
+  h2->GetYaxis()->SetRangeUser(ax2min, ax2max);
+  h2->GetYaxis()->SetTitle(Form("x_{%d}", ax2));
   h2->Draw();
-  
+
 
   if(!fTNodesDraw) fTNodesDraw = new TKDNodeInfo::TKDNodeDraw[fNTNodes]; 
   TKDNodeInfo::TKDNodeDraw *box = 0x0;
-  for(Int_t in=0; in<fNTNodes; in++){ 
-    TKDNodeInfo *node = (TKDNodeInfo*)((*fTNodes)[in]);
-
+  for(Int_t in=fNTNodes; in--; ){ 
     box = &(fTNodesDraw[in]);
-    box->SetNode(node, fNSize, ax1, ax2);
+    box->SetNode((TKDNodeInfo*)((*fTNodes)[in]), fNSize, ax1, ax2);
     box->Draw();
   }
+
   return;
 }
 
@@ -312,30 +338,3 @@ void TKDInterpolatorBase::SetAlpha(Float_t a)
   return;
 }
 
-//__________________________________________________________________
-void TKDInterpolatorBase::SetInterpolationMethod(Bool_t on)
-{
-// Set interpolation bit to "on".
-  
-  if(on) fStatus += fStatus&1 ? 0 : 1;
-  else fStatus += fStatus&1 ? -1 : 0;
-}
-
-
-//_________________________________________________________________
-void TKDInterpolatorBase::SetStore(Bool_t on)
-{
-// Set store bit to "on"
-  
-  if(on) fStatus += fStatus&2 ? 0 : 2;
-  else fStatus += fStatus&2 ? -2 : 0;
-}
-
-//_________________________________________________________________
-void TKDInterpolatorBase::SetWeights(Bool_t on)
-{
-// Set weights bit to "on"
-  
-  if(on) fStatus += fStatus&4 ? 0 : 4;
-  else fStatus += fStatus&4 ? -4 : 0;
-}
index 4b831bb..bb8c451 100644 (file)
@@ -33,7 +33,11 @@ public:
     kdN = 4       // increase in the number of PDF if fit failled
    ,kNhelper = 30 // bucket size in helper kdTree
   };
-
+  enum EKDInterpolatorBaseBits {
+    kCOG   = 0  // COG interpolation method
+   ,kSTORE = 1  // Store interpolation results
+   ,kWEIGHTS = 2 // use weights
+  };
   TKDInterpolatorBase(Int_t size = 0);
   virtual ~TKDInterpolatorBase();
 
@@ -41,18 +45,20 @@ public:
   virtual Int_t GetNodeIndex(const Float_t *p) = 0;
   Float_t    GetAlpha() const {return fAlpha;}
   Bool_t     GetCOGPoint(Int_t node, Float_t *&coord, Float_t &val, Float_t &error) const;
-  Bool_t     GetInterpolationMethod() const {return fStatus&1;}
   TKDNodeInfo* GetNodeInfo(Int_t inode) const;
   Int_t      GetNTNodes() const {return fNTNodes;}
-  void       GetStatus();
-  Bool_t     GetStore() const {return fStatus&2;}
-  Bool_t     GetWeights() const {return fStatus&4;}
+  Bool_t     GetRange(Int_t ax, Float_t &min, Float_t &max) const;
+  void       GetStatus(Option_t *opt="");
+
+  Bool_t     HasStore() const {return TESTBIT(fStatus, kSTORE);}
+  Bool_t     UseCOG() const {return TESTBIT(fStatus, kCOG);}
+  Bool_t     UseWeights() const {return TESTBIT(fStatus, kWEIGHTS);}
 
-  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       DrawProjection(UInt_t ax1 = 0, UInt_t ax2 = 1);
   void       SetAlpha(Float_t a);
-  void       SetInterpolationMethod(Bool_t on = kTRUE);
-  void       SetStore(Bool_t on = kTRUE);
-  void       SetWeights(Bool_t on = kTRUE);
+  void       SetCOG(Bool_t on = kTRUE) {on ? SETBIT(fStatus, kCOG) : CLRBIT(fStatus, kCOG);}
+  void       SetStore(Bool_t on = kTRUE) {on ? SETBIT(fStatus, kSTORE) : CLRBIT(fStatus, kSTORE);}
+  void       SetWeights(Bool_t on = kTRUE) {on ? SETBIT(fStatus, kWEIGHTS) : CLRBIT(fStatus, kWEIGHTS);}
 
 protected:
   virtual void      Build(Int_t nnodes);
index 9040a71..85b25eb 100644 (file)
@@ -1,3 +1,22 @@
+////////////////////////////////////////////////////////
+//
+// Bucket representation for TKDInterpolator(Base)
+//
+// The class store data and provides the interface to draw and print.
+// The bucket - generalized histogram bin in N dimensions is represented by unprocessed data like
+//   - experimental PDF value and statistical error 
+//   - COG position (n-tuplu)
+//   - boundaries
+// and interpolated data like
+//   - parameters of the local parabolic fit
+//   - their covariance matrix
+//  
+// For drawing 2D projections the helper class TKDNodeInfo::TKDNodeDraw is used.
+//
+// Author Alexandru Bercuci <A.Bercuci@gsi.de>
+//
+////////////////////////////////////////////////////////
+
 #include "TKDNodeInfo.h"
 
 #include "TRandom.h"
@@ -79,7 +98,7 @@ void TKDNodeInfo::Build(Int_t dim)
 //     Info("Build()", "...");
 
   if(!dim) return;
-  
+  fNDim = 3*dim;  
   Int_t lambda = Int_t(1 + dim + .5*dim*(dim+1));
   if(fData) delete [] fData;
   fData = new Float_t[fNDim];
@@ -91,6 +110,15 @@ void TKDNodeInfo::Build(Int_t dim)
 }
 
 //_________________________________________________________________
+void TKDNodeInfo::SetNode(Int_t ndim, Float_t *data, Float_t *pdf)
+{
+  Build(ndim);
+  memcpy(fData, data, fNDim*sizeof(Float_t));
+  fVal[0]=pdf[0]; fVal[1]=pdf[1];
+}
+
+
+//_________________________________________________________________
 void TKDNodeInfo::Print(const Option_t *) const
 {
   // Print the content of the node
@@ -99,9 +127,9 @@ void TKDNodeInfo::Print(const Option_t *) const
   for(int idim=0; idim<dim; idim++) printf("%f ", fData[idim]);
   printf("] f = [%f +- %f]\n", fVal[0], fVal[1]);
 
-  //   Float_t *bounds = &fData[dim];
+  Float_t *bounds = &fData[dim];
   printf("range[");
-  for(int idim=0; idim<dim; idim++) printf("(%f %f) ", fData[2*idim], fData[2*idim+1]);
+  for(int idim=0; idim<dim; idim++) printf("(%f %f) ", bounds[2*idim], bounds[2*idim+1]);
   printf("]\n");
   
   printf("Fit parameters : ");
@@ -160,7 +188,10 @@ Double_t TKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t
 
 
 //_________________________________________________________________
-TKDNodeInfo::TKDNodeDraw::TKDNodeDraw() : TBox(), fCOG()
+TKDNodeInfo::TKDNodeDraw::TKDNodeDraw() 
+  :TBox()
+  ,fCOG()
+  ,fNode(0x0)
 {
   SetFillStyle(3002);
   SetFillColor(50+Int_t(gRandom->Uniform()*50.));
@@ -181,17 +212,22 @@ void TKDNodeInfo::TKDNodeDraw::Draw(Option_t* option)
 //_________________________________________________________________
 void TKDNodeInfo::TKDNodeDraw::SetNode(TKDNodeInfo *node, UChar_t size, UChar_t ax1, UChar_t ax2)
 {
+  fNode=node;
   const Float_t kBorder = 0.;//1.E-4;
   Float_t *bounds = &(node->Data()[size]);
-  //printf("x1[%f] x2[%f] y1[%f] y2[%f]\n", bounds[2*ax1], bounds[2*ax1+1], bounds[2*ax2], bounds[2*ax2+1]);
   fX1=bounds[2*ax1]+kBorder;
   fX2=bounds[2*ax1+1]-kBorder;
   fY1=bounds[2*ax2]+kBorder;
   fY2=bounds[2*ax2+1]-kBorder;
   
   Float_t x(node->Data()[ax1]), y(node->Data()[ax2]);
-  //printf("x[%f] y[%f]\n", x, y);
   fCOG.SetX(x); fCOG.SetY(y);
 }
 
 
+//_________________________________________________________________
+void TKDNodeInfo::TKDNodeDraw::Print(const Option_t* option) const
+{
+  if(!fNode) return;
+  fNode->Print(option);
+}
index 2a5a988..198aaba 100644 (file)
@@ -1,6 +1,14 @@
 #ifndef ROOT_TKDNodeInfo
 #define ROOT_TKDNodeInfo
 
+////////////////////////////////////////////////////////
+//
+// Bucket representation for TKDInterpolator(Base)
+//
+// Author Alexandru Bercuci <A.Bercuci@gsi.de>
+//
+////////////////////////////////////////////////////////
+
 #ifndef ROOT_TObject
 #include "TObject.h"
 #endif
@@ -19,40 +27,47 @@ template <typename Value> class TMatrixT;
 typedef class TMatrixT<Double_t> TMatrixD;
 class TKDNodeInfo : public TObject
 {
-  //friend class TKDInterpolatorBase;
 public:
+  friend class TKDPDF;
+  friend class TKDInterpolatorBase;
   class TKDNodeDraw : public TBox {
   public:
     TKDNodeDraw();
     ~TKDNodeDraw() {;}
     void  Draw(Option_t* option = "");
+    void  Print(const Option_t* option = "") const; // *MENU*
     void  SetNode(TKDNodeInfo*, UChar_t s, UChar_t ax1=0, UChar_t ax2=1);
   private:
     TKDNodeDraw(const TKDNodeDraw & ref);
     TKDNodeDraw& operator=(const TKDNodeDraw & ref);
 
-    TMarker  fCOG;      // COG of the node
-
+    TMarker     fCOG;    // COG of the node
+    TKDNodeInfo *fNode;  //! node data
     ClassDef(TKDNodeDraw, 1)   // graphical representation of TKDNodeInfo
   };
 
   TKDNodeInfo(Int_t ndim = 0);
   TKDNodeInfo(const TKDNodeInfo & ref);
   TKDNodeInfo& operator=(const TKDNodeInfo & ref);
-  virtual  ~TKDNodeInfo();
-  Double_t  CookPDF(const Double_t *point, Double_t &result, Double_t &error);
-  inline Bool_t    Has(const Float_t *p) const;
-  void      Print(const Option_t * = "") const;
-  void      Store(const TVectorD &par, const TMatrixD &cov);
-
-  Int_t GetSize() const { return fNDim; }
-  Float_t *  Data() { return fData; } 
-  Float_t *  Val() { return fVal; } 
-  TMatrixD * Cov() { return fCov; }
-  TVectorD * Par() { return fPar; }
+  virtual ~TKDNodeInfo();
+  Double_t      CookPDF(const Double_t *point, Double_t &result, Double_t &error);
+  inline void   GetBoundary(Int_t ax, Float_t &min, Float_t &max) const;
+  inline void   GetCOG(Float_t* &p) const;
+  Int_t         GetDimension() const { return fNDim/3; }
+  void          GetPDF(Float_t &pdf, Float_t &error) const {pdf=fVal[0]; error=fVal[1];}
+  Int_t         GetSize() const { return fNDim; }
+  inline Bool_t Has(const Float_t *p) const;
+  void          Print(const Option_t * = "") const;
+  void          Store(const TVectorD &par, const TMatrixD &cov);
+
+  TMatrixD*     Cov() const  { return fCov; }
+  TVectorD*     Par() const  { return fPar; }
 
+  void          SetNode(Int_t ndim, Float_t *data, Float_t *pdf);
 protected:
-  void      Build(Int_t ndim);
+  Float_t*      Data() { return fData;}
+  Float_t*      Val()  { return &fVal[0]; }
+  void          Build(Int_t ndim);
 
 private:
   Int_t     fNDim;          // 3 times data dimension
@@ -65,15 +80,35 @@ private:
 };
 
 //_____________________________________________________________________
-Bool_t TKDNodeInfo::Has(const Float_t *p) const
+inline Bool_t TKDNodeInfo::Has(const Float_t *p) const
 {
-  Int_t n = 0;
   Int_t ndim = fNDim/3;
-  for(int id=0; id<ndim; id++) if(p[id]>=fData[ndim+2*id] && p[id]<fData[ndim+2*id+1]) n++;
+  
+  Float_t *it = &fData[ndim]; Int_t n(0);
+  for(int id=0; id<ndim; id++, it+=2)
+    if(p[id]>=it[0] && p[id]<it[1]) n++;
   if(n==ndim) return kTRUE;
   return kFALSE;
 }
 
+//_____________________________________________________________________
+inline void TKDNodeInfo::GetBoundary(Int_t ax, Float_t &min, Float_t &max) const
+{
+  Int_t ndim = fNDim/3;
+  if(ax<0 || ax>=ndim){
+    min=0.; max=0.;
+    return;
+  }
+  Float_t *it = &fData[ndim+(ax<<1)];
+  min = it[0]; max = it[1];
+}
+
+//_____________________________________________________________________
+inline void TKDNodeInfo::GetCOG(Float_t* &p) const
+{
+  Int_t ndim = fNDim/3;
+  memcpy(p, fData, ndim*sizeof(Float_t));
+}
 
 #endif