new macros for testing the STAT package
authorabercuci <abercuci@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 25 Oct 2007 13:55:12 +0000 (13:55 +0000)
committerabercuci <abercuci@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 25 Oct 2007 13:55:12 +0000 (13:55 +0000)
STAT/Macros/interpolCompare.C [new file with mode: 0644]
STAT/Macros/interpolTest.C [new file with mode: 0644]
STAT/Macros/kDTreeTest.C [new file with mode: 0644]
STAT/Macros/kNNSpeed.C [new file with mode: 0644]
STAT/Macros/kNNTest.C [new file with mode: 0644]
STAT/Macros/pdfCompare.C [new file with mode: 0644]
STAT/Macros/pdfIO.C [new file with mode: 0644]
STAT/Macros/pdfTest.C [new file with mode: 0644]

diff --git a/STAT/Macros/interpolCompare.C b/STAT/Macros/interpolCompare.C
new file mode 100644 (file)
index 0000000..6a1ad2e
--- /dev/null
@@ -0,0 +1,127 @@
+void interpolCompare(const Int_t method = 0, const int nstat = 40000)
+{
+// Compare interpolation methods.
+// method = 0 : Select COG method
+// method = 1 : Select INT method
+
+       const Float_t alpha = 20.;
+       const int bucket = nstat/100; //400;
+       const int nevals = 100;
+       Float_t *data[1];
+       Float_t data1[nstat];
+       for(int istat=0; istat<nstat; istat++) data1[istat] = gRandom->Gaus();
+       data[0] = &data1[0];
+
+       TLegend *leg = new TLegend(0.7157258,0.8280255,0.9637097,0.9596603);
+       leg->SetBorderSize(1);
+       TF1 *f=new TF1("f1", "gaus(0);x;f(x)", -5., 5.);
+       f->SetLineWidth(1);
+       f->SetLineColor(2);
+       f->SetParameter(0, 1.);
+       f->SetParameter(1, 0.);
+       f->SetParameter(2, 1.);
+       f->SetParameter(0, 1./f->Integral(-10., 10.));
+       f->Draw();      leg->AddEntry(f, "model", "l");
+       
+       
+       Float_t cog, val, val_err;
+       Double_t x, chi2, chi2pdf, chi2int, result, err;
+       Int_t npPDF, npINT;
+       TKDPDF pdf(nstat, 1, bucket, data);
+       pdf.SetInterpolationMethod(method);
+       pdf.SetAlpha(alpha);
+       
+       TKDInterpolator in(1, pdf.GetNTNodes());
+       in.SetInterpolationMethod(method);
+       in.SetAlpha(alpha);
+       
+       // change to the reprezentation which is closer to the local
+       // interpolating polynom.
+       Double_t fx;
+       Int_t inode = 0;
+       TKDNodeInfo *node = 0x0;
+       while(node = pdf.GetNodeInfo(inode)){
+               if(node->fVal[0] > 0.){
+                       in.SetNode(inode, *node);
+                       node  = in.GetNodeInfo(inode);
+                       fx = node->fVal[0];
+                       node->fVal[0] = TMath::Log(fx);
+                       node->fVal[1] = node->fVal[1]/fx;
+               }
+               inode++;
+       }
+
+       // Do preprocessing for INT
+       Float_t *pcog;
+       if(method) for(int icog=0; icog<pdf.GetNTNodes(); icog++){
+               pdf.GetCOGPoint(icog, pcog, val, val_err);
+               x = pcog[0];
+               pdf.Eval(&x, result, err);
+               in.Eval(&x, result, err);
+       }
+
+
+       
+       // prepare drawing containers
+       TGraphErrors *gPDF = new TGraphErrors(nevals);
+       gPDF->SetMarkerStyle(7);
+       TGraph *gREP = new TGraph(nevals);
+       gREP->SetMarkerStyle(7);
+       gREP->SetMarkerColor(4);
+       gREP->SetLineColor(4);
+       
+       TGraphErrors *gINT = new TGraphErrors(nevals);
+       gINT->SetMarkerStyle(24);
+       gINT->SetMarkerSize(.6);
+       gINT->SetMarkerColor(3);
+       gINT->SetLineColor(3);
+       TGraph *gREI = new TGraph(nevals);
+       gREI->SetMarkerStyle(24);
+       gREI->SetMarkerColor(4);
+       gREI->SetLineColor(4);
+       
+       // Do the interpolation with PDF and Interpolator
+       Double_t dx = 10./nevals;
+       x = 5.-dx/2.;
+       chi2pdf = 0.; npPDF = 0;
+       chi2int = 0.; npINT = 0;
+       for(int ip=0; ip<nevals; ip++){
+               chi2 = pdf.Eval(&x, result, err);
+               gPDF->SetPoint(ip, x, result);
+               gPDF->SetPointError(ip, 0., err);
+               if(TMath::Abs(x)<3.){
+                       chi2pdf += (result - f->Eval(x))*(result - f->Eval(x))/f->Eval(x)/f->Eval(x);
+                       npPDF++;
+               }
+               gREP->SetPoint(ip, x, .1*err/TMath::Abs(result));
+               //printf("%2d x[%f] r[%f +- %f] chi2[%e]\n", ip, x, result, err, chi2);
+
+               chi2 = in.Eval(&x, result, err);
+               result = TMath::Exp(result);
+               err = err*result;
+               gINT->SetPoint(ip, x, result);
+               gINT->SetPointError(ip, 0., err);
+               if(TMath::Abs(x)<3.){
+                       chi2int += (result - f->Eval(x))*(result - f->Eval(x))/f->Eval(x)/f->Eval(x);
+                       npINT++;
+               }
+               gREI->SetPoint(ip, x, .1*err/TMath::Abs(result));
+               //printf("%2d x[%f] r[%f +- %f] chi2[%e]\n\n", ip, x, result, err, chi2);
+
+               x -= dx;
+       }
+       gINT->Draw("p"); leg->AddEntry(gINT, Form("Interpolator [%s]", method ? "INT" : "COG"), "pl");
+       gPDF->Draw("p"); leg->AddEntry(gPDF, Form("PDF [%s]", method ? "INT" : "COG"), "pl");   
+       gREI->Draw("pl"); leg->AddEntry(gREI, "1% #epsilon Interpolator", "pl");
+       gREP->Draw("pl"); leg->AddEntry(gREP, "1% #epsilon PDF", "pl");
+
+       chi2pdf /= npPDF;
+       printf("PDF : chi2(%d) = %f\n", npPDF, chi2pdf);
+       chi2int /= npINT;
+       printf("INT : chi2(%d) = %f\n", npINT, chi2int);
+
+       leg->Draw();
+}
+
+
+
diff --git a/STAT/Macros/interpolTest.C b/STAT/Macros/interpolTest.C
new file mode 100644 (file)
index 0000000..e7e1f9f
--- /dev/null
@@ -0,0 +1,241 @@
+// #include "TSystem.h"
+// #include "TFile.h"
+// #include "TTree.h"
+// #include "TProfile.h"
+// #include "TF1.h"
+// #include "TF2.h"
+// #include "TF3.h"
+// #include "TLegend.h"
+// #include "TRandom.h"
+// #include "TString.h"
+// #include "TGraph.h"
+// #include "TMarker.h"
+// #include "TStopwatch.h"
+// #include "../src/TKDPDF.h"
+// #include "../src/TKDInterpolator.h"
+
+void interpolTest(const Int_t ndim = 1);
+Double_t interpolation(const int nstat, const int ndim, const int first);
+void build(const int ndim, const int npoints);
+
+TF1 *f = 0x0;
+
+//__________________________________________________________
+void interpolTest(const Int_t ndim)
+{
+// Test interpolator on a variable dimension (ndim<=5) uncorrelated
+// Gauss model. The macro displays the chi2 between interpolation and
+// model for various statistics of experimental data.
+//
+// Check the worker function interpolation() to fine tune the
+// algorithm.
+// 
+// Attention:
+// The macro works in compiled mode !
+
+       // define model
+       switch(ndim){
+       case 1:
+               f=new TF1("f1", "gaus(0);x;f(x)", -5., 5.);
+               f->SetParameter(0, 1.);
+               f->SetParameter(1, 0.);
+               f->SetParameter(2, 1.);
+               f->SetParameter(0, 1./f->Integral(-5., 5.));
+               break;
+       case 2:
+               f=new TF2("f2", "[0]*exp(-0.5*((x-[1])/[2])**2)*exp(-0.5*((y-[3])/[4])**2)", -5., 5., -5., 5.);
+               f->SetParameter(0, 1.);
+               f->SetParameter(1, 0.);
+               f->SetParameter(2, 1.);
+               f->SetParameter(3, 0.);
+               f->SetParameter(4, 1.);
+               f->SetParameter(0, 1./f->Integral(-5., 5., -5., 5.));
+               break;
+       case 3:
+       case 4:
+       case 5:
+               f=new TF3("f3", "[0]*exp(-0.5*((x-[1])/[2])**2)*exp(-0.5*((y-[3])/[4])**2)*exp(-0.5*((z-[5])/[6])**2)", -5., 5., -5., 5., -5., 5.);
+               f->SetParameter(0, 1.);
+               f->SetParameter(1, 0.);
+               f->SetParameter(2, 1.);
+               f->SetParameter(3, 0.);
+               f->SetParameter(4, 1.);
+               f->SetParameter(5, 0.);
+               f->SetParameter(6, 1.);
+               f->SetParameter(0, 1./f->Integral(-5., 5., -5., 5., -5., 5.));
+               if(ndim>3) printf("Warning : chi2 results are unreliable !\n");
+               break;
+       default:
+               printf("%dD not supported in this macro.\n", ndim);
+               return;
+       }
+
+       Double_t chi2;
+       const Int_t nchi2 = 18;
+       const Int_t nfirst = 1;
+       //TProfile *hp = new TProfile("hp", "", 100, 1.E4, 1.E6);
+       //hp->SetMarkerStyle(24);
+       TGraph *gChi2 = new TGraph(nchi2);
+       gChi2->SetMarkerStyle(24);
+       Int_t ip = 0, nstat, first;
+       for(int ifirst=0; ifirst<nfirst; ifirst++){
+               first = 0;//Int_t(gRandom->Uniform(1.E4));
+               for(int i=2; i<10; i++){
+                       nstat = 10000*i;
+                       chi2 = interpolation(nstat, ndim, first);
+                       gChi2->SetPoint(ip++, TMath::Log10(float(nstat)), chi2);
+                       //hp->Fill(float(nstat), chi2);
+               }
+               for(int i=1; i<10; i++){
+                       nstat = 100000*i;
+                       chi2 = interpolation(nstat, ndim, first);
+                       gChi2->SetPoint(ip++, TMath::Log10(float(nstat)), chi2);
+                       //hp->Fill(float(nstat), chi2);
+               }
+               gChi2->Draw("apl");
+       }
+       //hp->Draw("pl");
+}
+
+void test()
+{
+       TClonesArray ar("TKDNodeInfo", 10);
+}
+
+//__________________________________________________________
+Double_t interpolation(const int nstat, const int ndim, const int first)
+{
+// Steer interpolation of a nD (n<=5) uncorrelated Gauss distribution.
+// The user is suppose to give the experimental statistics (nstat) the
+// dimension of the experimental point (ndim). The entry from which the
+// data base is being read is steered by "first".
+
+       const int bs = 400;
+       const int npoints = 100;
+       // switch on chi2 calculation between interpolation and model.
+       // Default calculates distance.
+       const Bool_t kCHI2 = kFALSE; 
+       const Bool_t kINT  = kFALSE; // switch on integral interpolator
+
+       // open data base
+       TString fn("5D_Gauss.root");
+       if(!gSystem->FindFile(".", fn)) build(5, 1000000);
+       TFile::Open("5D_Gauss.root");
+       TTree *t = (TTree*)gFile->Get("db");
+
+       // build interpolator
+       TString var = "x0"; for(int idim=1; idim<ndim; idim++) var += Form(":x%d", idim);
+       TKDPDF pdf(t, var.Data(), "", bs, nstat, first);
+       printf("\n\nFINISH BUILDING PDF\n\n");
+       
+       Double_t fx;
+       Int_t inode = 0;
+       TKDNodeInfo *node = 0x0;
+       TKDInterpolator in(ndim, pdf.GetNTNodes());
+/*     while(node = pdf.GetNodeInfo(inode)){
+               if(node->fVal[0] > 0.){
+                       fx = node->fVal[0];
+                       node->fVal[0] = TMath::Log(fx);
+                       node->fVal[1] = node->fVal[1]/fx;
+               }
+               in.SetNode(inode, *node);
+               inode++;
+       }*/
+       //pdf.SetInterpolationMethod(kINT);
+       pdf.SetAlpha(2.);
+       //in.SetStore();
+
+
+       Double_t x[2], r, e, chi2;
+       TH2 *h2 = new TH2F("h2", "", 50, -4., 4., 50, -4., 4.);
+       TAxis *ax = h2->GetXaxis(), *ay = h2->GetYaxis();
+       for(int ix=1; ix<=ax->GetNbins(); ix++){
+               x[0] = ax->GetBinCenter(ix);
+               for(int iy=1; iy<=ay->GetNbins(); iy++){
+                       x[1] = ay->GetBinCenter(iy);
+
+                       chi2 = pdf.Eval(x, r, e);
+                       printf("x[%2d] x[%2d] r[%f] e[%f] chi2[%f]\n", ix, iy, r, e, chi2);
+                       h2->SetBinContent(ix, iy, r/*TMath::Exp(r)*/);
+               }
+       }
+       h2->Draw("lego2");
+       return;
+       
+       // define testing grid
+       Double_t x[5], r, e;
+       Double_t dx[5] = {4., 4., 4., 4., 4.};
+       Double_t chi2, theor, result, err;
+       for(int idim=0; idim<ndim; idim++){
+               dx[idim] = (ndim<<2)/float(npoints);
+               x[idim]  = -2.+dx[idim]/2.;
+       }
+
+       
+       // setup the integral interpolator and do memory test
+       Float_t c[5], v, ve;
+       for(int ip=0; ip<in.GetNTNodes(); ip++){
+               in.GetCOGPoint(ip, c, v, ve);
+               for(int idim=0; idim<ndim; idim++) x[idim] = (Double_t)c[idim];
+               in.Eval(x, result, err);
+       }
+       printf("\nInterpolating %d data points in %dD ...\n", nstat, ndim);
+       printf("\tbucket size %d\n", bs);
+       printf("\tstarting interpolation ...\n");
+       return;
+
+       
+       // evaluate relative error
+       for(int idim=0; idim<ndim; idim++) x[idim] = 0.;
+       in.Eval(x, result, err);
+       Double_t threshold = err/result;
+       
+       // starting interpolation over the testing grid
+       chi2 = 0.; Int_t nsample = 0;
+       x[4]  = -2.+dx[4]/2.;
+       do{
+               x[3]  = -2.+dx[3]/2.;
+               do{
+                       x[2]  = -2.+dx[2]/2.;
+                       do{
+                               x[1]  = -2.+dx[1]/2.;
+                               do{
+                                       x[0]  = -2.+dx[0]/2.;
+                                       do{
+                                               in.Eval(x, result, err);
+                                               if(err/TMath::Abs(result) < 1.1*threshold){
+                                                       theor = f->Eval(x[0], x[1], x[2]);
+                                                       Double_t tmp = result - theor; tmp *= tmp;
+                                                       chi2 += tmp/(kCHI2 ? err*err : theor);
+                                                       nsample++;
+                                               }
+                                       } while((x[0] += dx[0]) < 2.);
+                               } while((x[1] += dx[1]) < 2.);
+                       } while((x[2] += dx[2]) < 2.);
+               } while((x[3] += dx[3]) < 2.);
+       } while((x[4] += dx[4]) < 2.);
+       chi2 /= float(nsample);
+       printf("\tinterpolation quality (chi2) %f\n", chi2);
+       
+       return kCHI2 ? chi2 : TMath::Sqrt(chi2);
+}
+
+
+//___________________________________________________________
+void build(const int ndim, const int npoints)
+{
+       printf("Building DB ...\n");
+       Float_t data[ndim];
+       TFile::Open(Form("%dD_Gauss.root", ndim), "RECREATE");
+       TTree *t = new TTree("db", Form("%dD data base for kD statistics", ndim));
+       for(int idim=0; idim<ndim; idim++) t->Branch(Form("x%d", idim), &data[idim], Form("x%d/F", idim));
+
+       for (Int_t ip=0; ip<npoints; ip++){
+               for (Int_t id=0; id<ndim; id++) data[id]= gRandom->Gaus();
+               t->Fill();
+       }
+
+       t->Write();
+       gFile->Close();
+       delete gFile;
+}
diff --git a/STAT/Macros/kDTreeTest.C b/STAT/Macros/kDTreeTest.C
new file mode 100644 (file)
index 0000000..80d8532
--- /dev/null
@@ -0,0 +1,229 @@
+#include <malloc.h>
+#include "TMatrixD.h"
+#include "TRandom.h"
+#include "TGraph.h"
+#include "TStopwatch.h"
+#include "../src/TKDTree.h"
+
+
+Float_t Mem()
+{
+  // size in KB
+  static struct mallinfo memdebug; 
+  memdebug = mallinfo();
+  return memdebug.uordblks/1024.;
+}
+
+
+void TestBuild(const Int_t npoints = 1000000, const Int_t bsize = 100);
+void TestSpeed(const Int_t npower2 = 20, const Int_t bsize = 10);
+
+void testindexes(Int_t nloop, Int_t npoints);
+void  testkdtreeIF(Int_t npoints=1000, Int_t bsize=9, Int_t nloop=1000, Int_t mode = 2);
+void TestSizeIF(Int_t npoints=1000, Int_t nrows = 150, Int_t nsec=18, Int_t mode = 1, Int_t bsize=9);
+
+//______________________________________________________________________
+void kDTreeTest()
+{
+       printf("\n\tTesting kDTree memory usage ...\n");
+       TestBuild();
+       
+       printf("\n\tTesting kDTree speed ...\n");
+       TestSpeed();
+}
+
+//______________________________________________________________________
+void TestBuild(const Int_t npoints, const Int_t bsize){  
+  Float_t *data0 =  new Float_t[npoints*2];
+  Float_t *data[2];
+  data[0] = &data0[0];
+  data[1] = &data0[npoints];
+  for (Int_t i=0;i<npoints;i++) {
+    data[1][i]= gRandom->Rndm();
+    data[0][i]= gRandom->Rndm();
+  }
+  Float_t before =Mem();
+  TKDTreeIF *kdtree = new TKDTreeIF(npoints, 2, bsize, data);
+       Float_t after = Mem();
+  printf("Memory usage %f KB\n",after-before);
+       delete kdtree;
+  Float_t end = Mem();
+  printf("Memory leak %f KB\n", end-before);
+       return; 
+}
+
+//______________________________________________________________________
+void TestSpeed(const Int_t npower2, const Int_t bsize)
+{
+       if(npower2 < 10){
+               printf("Please specify a power of 2 greater than 10\n");
+               return;
+       }
+       
+       Int_t npoints = Int_t(pow(2., npower2))*bsize;
+       Float_t *data0 =  new Float_t[npoints*2];
+  Float_t *data[2];
+  data[0] = &data0[0];
+  data[1] = &data0[npoints];
+  for (Int_t i=0;i<npoints;i++) {
+    data[1][i]= gRandom->Rndm();
+    data[0][i]= gRandom->Rndm();
+  }
+
+       TGraph *g = new TGraph(npower2-10);
+       g->SetMarkerStyle(7);
+       TStopwatch timer;
+       Int_t tpoints;
+  TKDTreeIF *kdtree = 0x0;
+       for(int i=10; i<npower2; i++){
+               tpoints = Int_t(pow(2., i))*bsize;
+               timer.Start(kTRUE);
+               kdtree = new TKDTreeIF(tpoints, 2, bsize, data);
+               timer.Stop();
+               g->SetPoint(i-10, i, timer.CpuTime());
+               printf("npoints [%d] nodes [%d] cpu time %f [s]\n", tpoints, kdtree->GetNNodes(), timer.CpuTime());
+               //timer.Print("u");
+               delete kdtree;
+       }
+       g->Draw("apl");
+       return;
+}
+
+//______________________________________________________________________
+void TestSizeIF(Int_t npoints, Int_t nrows, Int_t nsec, Int_t mode, Int_t bsize)
+{
+  //
+  // test size to build kdtree
+  //
+  Float_t before =Mem();
+  for (Int_t isec=0; isec<nsec;isec++)
+    for (Int_t irow=0;irow<nrows;irow++){
+      testkdtreeIF(npoints,0,mode,bsize);
+    }
+  Float_t after = Mem();
+  printf("Memory usage %f\n",after-before);
+}
+
+
+void testindexes(Int_t nloop, Int_t npoints){
+  //
+  // test indexing 
+  // 
+  TKDTree<Int_t, Float_t> *kdtree = new TKDTree<Int_t, Float_t>(0,0,0,0);
+  Int_t row =0;
+  Int_t collumn =0; 
+  TStopwatch timer;
+  timer.Start();
+  row = 10;
+  for (Int_t iloop=0;iloop<nloop;iloop++)
+  for (Int_t index=1;index<npoints;index++){
+    //row = TMath::Log2(index);
+    //row=0;
+   //  if (index< (16<<row)) row=0;
+//     for (; index>=(32<<row);row+=5);
+//     for (; index>=(2<<row);row++);
+//     collumn= index-(1<<row);
+      TKDTree<Int_t, Float_t>::GetCoord(index,row,collumn);
+    //
+    //Int_t index2=kdtree->GetIndex(row,collumn);
+    //printf("%d\t%d\t%d\t%d\n",index,index2,row,collumn);
+    if (kdtree->GetIndex(row,collumn)!=index || collumn<0) {
+      printf("Problem\n");
+    }
+  }
+  timer.Stop();
+  timer.Print();
+}
+
+
+void  testkdtreeIF(Int_t npoints, Int_t bsize, Int_t nloop, Int_t mode)
+{
+//
+// Test speed and functionality of 2D kdtree.
+// Input parametrs:
+// npoints - number of data points
+// bsize   - bucket size
+// nloop   - number of loops
+// mode    - tasks to be performed by the kdTree
+//         - 0  : time building the tree
+//
+
+  Float_t rangey  = 100;
+  Float_t rangez  = 100;
+  Float_t drangey = 0.1;
+  Float_t drangez = 0.1;
+
+  //
+  Float_t *data0 =  new Float_t[npoints*2];
+  Float_t *data[2];
+  data[0] = &data0[0];
+  data[1] = &data0[npoints];
+  Int_t i;   
+  for (i=0; i<npoints; i++){
+    data[0][i]          = gRandom->Uniform(-rangey, rangey);
+    data[1][i]          = gRandom->Uniform(-rangez, rangez);
+  }
+  TStopwatch timer;
+
+       // check time build
+       printf("building kdTree ...\n");
+       timer.Start(kTRUE);
+  TKDTreeIF *kdtree = new TKDTreeIF(npoints, 2, bsize, data);
+  timer.Stop();
+  timer.Print();
+       if(mode == 0) return;
+       
+
+
+       Float_t countern=0;
+       Float_t counteriter  = 0;
+       Float_t counterfound = 0;
+
+
+       if (mode ==2){
+               if (nloop) timer.Start(kTRUE);
+               Int_t res[npoints];
+               Int_t nfound = 0;
+               for (Int_t kloop = 0;kloop<nloop;kloop++){
+                       if (kloop==0){
+                               counteriter = 0;
+                               counterfound= 0;
+                               countern    = 0;
+                       }
+                       for (Int_t i=0;i<npoints;i++){
+                               Float_t point[2]={data[0][i],data[1][i]};
+                               Float_t delta[2]={drangey,drangez};
+                               Int_t iter  =0;
+                               nfound =0;
+                               Int_t bnode =0;
+                               //kdtree->FindBNode(point,delta, bnode);
+                               //continue;
+                               kdtree->FindInRangeA(point,delta,res,nfound,iter,bnode);
+                               if (kloop==0){
+                                       //Bool_t isOK = kTRUE;
+                                       Bool_t isOK = kFALSE;
+                                       for (Int_t ipoint=0;ipoint<nfound;ipoint++)
+                                               if (res[ipoint]==i) isOK =kTRUE;
+                                       counteriter+=iter;
+                                       counterfound+=nfound;
+                                       if (isOK) {
+                                               countern++;
+                                       }else{
+                                               printf("Bug\n");
+                                       }
+                               }
+                       }
+               }
+               
+               if (nloop){
+                       timer.Stop();
+                       timer.Print();
+               }
+       }
+       delete [] data0;
+
+       counteriter/=npoints;
+       counterfound/=npoints;
+       if (nloop) printf("Find nearest point:\t%f\t%f\t%f\n",countern, counteriter, counterfound);
+}
diff --git a/STAT/Macros/kNNSpeed.C b/STAT/Macros/kNNSpeed.C
new file mode 100644 (file)
index 0000000..cb8d133
--- /dev/null
@@ -0,0 +1,149 @@
+#include "TSystem.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TH2I.h"
+#include "TLegend.h"
+#include "TRandom.h"
+#include "TString.h"
+#include "TGraph.h"
+#include "TMarker.h"
+#include "TStopwatch.h"
+#include "../src/TKDInterpolator.h"
+#include "../src/TKDTree.h"
+
+const int ntimes = 500000; // times to repeate kNN search
+TStopwatch timer;
+
+Float_t kdTreeNN(Float_t **p, const int ndim = 5, const int np = 10000);
+Float_t standAloneNN(Float_t **p, const int ndim = 5, const int np = 10000);
+
+//______________________________________________________________
+TTree *db = 0x0;
+void kNNSpeed(const Int_t method = 0)
+{
+// Estimate kNN algorithm speed:
+// method = 0 : run kDTree kNN search
+// method = 1 : run kNN loop search
+
+       const Int_t ndim  = 5;
+       const Int_t nstat = 9;
+
+
+       // generate test sample
+       Float_t **p = new Float_t*[ntimes];
+       for(int ip =0 ; ip<ntimes; ip++){
+               p[ip] = new Float_t[ndim];
+               for(int idim=0; idim<ndim; idim++) p[ip][idim] = gRandom->Uniform(-.5, .5);
+       }
+
+       Float_t time;
+       Int_t stat;
+       TLegend *leg = new TLegend(.7, .7, .9, .9);
+       TH1 *h1 = new TH1I("h1", "", 100, 9.5, 20.5);
+       h1->SetMaximum(1.E2);
+       h1->SetMinimum(1.E-1);
+       h1->GetXaxis()->SetTitle("Statistics");
+       h1->GetYaxis()->SetTitle("CPU [#mus]");
+       h1->SetStats(0);
+       h1->Draw();
+       TGraph **g=new TGraph*[ndim];
+       for(int idim=2; idim<3/*ndim*/; idim++){
+               g[idim] = new TGraph(nstat);
+               g[idim]->SetMarkerStyle(20+idim);
+               g[idim]->SetMarkerColor(idim+1);
+               leg->AddEntry(g[idim], Form("%dD", idim+1), "pl");
+               stat = 16384;//1024;
+               for(int istat = 0; istat<nstat; istat++){
+                       time = kdTreeNN(p, idim+1, stat);
+                       //time = standAloneNN(p, idim+1, stat);
+                       g[idim]->SetPoint(istat, TMath::Log2(float(stat)), time);
+                       stat *= 2;
+               }
+               g[idim]->Draw("pl");
+       }
+       leg->Draw();
+       
+       for(int ip=0; ip<ntimes; ip++) delete [] p[ip];
+       delete [] p;
+}
+
+//______________________________________________________________
+Float_t kdTreeNN(Float_t **p, const int ndim, const int np)
+{
+       Float_t **x = new Float_t*[ndim];
+       for(int idim =0 ; idim<ndim; idim++){
+               x[idim] = new Float_t[np];
+               for(int ip=0; ip<np; ip++) x[idim][ip] = gRandom->Gaus();
+       }
+       TKDTreeIF nnFinder(np, ndim, 1, x);
+       nnFinder.MakeBoundaries();
+       
+       Int_t *index, fNN = 1;
+       Float_t *d;
+       timer.Start(kTRUE);
+       for(int itime=0; itime<ntimes; itime++){
+               nnFinder.FindNearestNeighbors(p[itime], fNN, index, d);
+       }
+       timer.Stop();
+       Float_t time = 1.E6*timer.CpuTime()/float(ntimes);
+
+       printf("np[%7d] nd[%d] time = %5.2f[mus]\n", np, ndim, time);
+       
+       // clean
+       for(int idim=0; idim<ndim; idim++) delete [] x[idim];
+       delete [] x;
+       
+
+       return time;
+}
+
+//______________________________________________________________
+Float_t standAloneNN(Float_t **p, const int ndim, const int np)
+{      
+       // STAND ALONE
+       const Int_t kNN = 1;
+       Float_t ftmp, gtmp, dist;
+       Int_t itmp, jtmp;
+       Int_t   fkNN[kNN];
+       Float_t fkNNdist[kNN];
+       for(int i=0; i<kNN; i++) fkNNdist[i] = 9999.;
+
+       Float_t **x = new Float_t*[ndim];
+       for(int idim =0 ; idim<ndim; idim++){
+               x[idim] = new Float_t[np];
+               for(int ip=0; ip<np; ip++) x[idim][ip] = gRandom->Gaus();
+       }
+       Int_t npoints = np;
+       
+       // calculate
+       Int_t ntimes = 100;
+       timer.Start(kTRUE);
+       for(int it=0; it<ntimes; it++){
+               for(int idx=0; idx<npoints; idx++){
+                       // calculate distance in the L1 metric
+                       dist = 0.;
+                       for(int idim=0; idim<ndim; idim++) dist += TMath::Abs(p[it*10][idim] - x[idim][idx]);
+                       if(dist >= fkNNdist[kNN-1]) continue;
+       
+                       //insert neighbor
+                       int iNN=0;
+                       while(dist >= fkNNdist[iNN]) iNN++;
+                       itmp = fkNN[iNN]; ftmp = fkNNdist[iNN];
+                       fkNN[iNN]     = idx;
+                       fkNNdist[iNN] = dist;
+                       for(int ireplace=iNN+1; ireplace<kNN; ireplace++){
+                               jtmp = fkNN[ireplace]; gtmp = fkNNdist[ireplace];
+                               fkNN[ireplace] = itmp; fkNNdist[ireplace] = ftmp;
+                               itmp = jtmp; ftmp = gtmp;
+                               if(ftmp == 9999.) break;
+                       }
+               }
+       }
+       timer.Stop();
+       
+       Float_t time = timer.CpuTime()/float(ntimes);
+       printf("np[%5d] nd[%d] time[%f] %6.2f\n", np, ndim, time, timer.CpuTime());
+
+       return time;
+}
+
diff --git a/STAT/Macros/kNNTest.C b/STAT/Macros/kNNTest.C
new file mode 100644 (file)
index 0000000..afcede3
--- /dev/null
@@ -0,0 +1,205 @@
+#include <malloc.h>
+
+#include "TSystem.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TH2I.h"
+#include "TLegend.h"
+#include "TRandom.h"
+#include "TString.h"
+#include "TGraph.h"
+#include "TMarker.h"
+#include "TStopwatch.h"
+#include "../src/TKDPDF.h"
+#include "../src/TKDTree.h"
+
+void kNNTest(const int np = 10000, const int ndim = 2);
+void kNNDraw(const Float_t *p, const int kNN=20);
+void build(const Int_t ndim = 2, const Int_t nstat = 1000000);
+Float_t Mem();
+
+
+//______________________________________________________________
+void kNNTest(const int np, const int ndim)
+{
+// Test speed and quality of nearest neighbors search.
+// The results are compared with a simple loop search through the data.
+// The macro should run in compiled mode.
+
+       const int ntimes = 1000; // times to repeate kNN search
+       const Int_t kNN = 1; // number of neighbors
+
+       Float_t **x = new Float_t*[ndim];
+       for(int idim =0 ; idim<ndim; idim++){
+               x[idim] = new Float_t[np];
+               for(int ip=0; ip<np; ip++) x[idim][ip] = gRandom->Gaus();
+       }
+       TKDTreeIF nnFinder(np, ndim, 1, x);
+       
+       // generate test sample
+       Float_t **p = new Float_t*[ntimes];
+       for(int ip =0 ; ip<ntimes; ip++){
+               p[ip] = new Float_t[ndim];
+               for(int idim=0; idim<ndim; idim++) p[ip][idim] = gRandom->Gaus();
+       }
+
+
+       Int_t *index;
+       Float_t *d;
+       TStopwatch timer;
+       timer.Start(kTRUE);
+       for(int itime=0; itime<ntimes; itime++) nnFinder.FindNearestNeighbors(p[itime], kNN, index, d);
+       timer.Stop();
+       printf("kDTree NN calculation.\n");
+       printf("cpu = %f [ms] real = %f [ms]\n", 1.E3*timer.CpuTime()/ntimes, 1.E3*timer.RealTime()/ntimes);
+       printf("Points indexes in increasing order of their distance to the target point.\n");
+       for(int i=0; i<kNN; i++){
+               printf("%5d [%7.4f] ", index[i], d[i]);
+               if(i%5==4) printf("\n");
+       }
+       printf("\n");
+       
+       // draw kNN
+       TLegend *leg = new TLegend(.7, .7, .9, .9);
+       leg->SetBorderSize(1);
+       leg->SetHeader("NN finders");
+       TH2 *h2 = new TH2I("h2", "", 100, -2., 2., 100, -2., 2.);
+       h2->Draw(); h2->SetStats(kFALSE);
+       TMarker *mp = new TMarker(p[ntimes-1][0], p[ntimes-1][1], 3);
+       mp->SetMarkerColor(4);
+       mp->Draw(); leg->AddEntry(mp, "Target", "p");
+       TGraph *gKD = new TGraph(kNN);
+       gKD->SetMarkerStyle(24);
+       gKD->SetMarkerColor(2);
+       gKD->SetMarkerSize(.5);
+       for(int ip=0; ip<kNN; ip++) gKD->SetPoint(ip, x[0][index[ip]], x[1][index[ip]]);
+       gKD->Draw("p"); leg->AddEntry(gKD, "kDTree", "p");
+       
+
+               
+       // STAND ALONE
+       Float_t ftmp, gtmp, dist;
+       Int_t itmp, jtmp;
+       Int_t   fkNN[kNN];
+       Float_t fkNNdist[kNN];
+       for(int i=0; i<kNN; i++) fkNNdist[i] = 9999.;
+
+       // calculate
+       timer.Start(kTRUE);
+       for(int idx=0; idx<np; idx++){
+               // calculate distance in the L1 metric
+               dist = 0.;
+               for(int idim=0; idim<ndim; idim++) dist += TMath::Abs(p[ntimes-1][idim] - x[idim][idx]);
+               if(dist >= fkNNdist[kNN-1]) continue;
+
+               //insert neighbor
+               int iNN=0;
+               while(dist >= fkNNdist[iNN]) iNN++;
+               itmp = fkNN[iNN]; ftmp = fkNNdist[iNN];
+               fkNN[iNN]     = idx;
+               fkNNdist[iNN] = dist;
+               for(int ireplace=iNN+1; ireplace<kNN; ireplace++){
+                       jtmp = fkNN[ireplace]; gtmp = fkNNdist[ireplace];
+                       fkNN[ireplace] = itmp; fkNNdist[ireplace] = ftmp;
+                       itmp = jtmp; ftmp = gtmp;
+                       if(ftmp == 9999.) break;
+               }
+       }
+       timer.Stop();
+       printf("\nStand Alone NN calculation.\n");
+       printf("cpu = %f [ms] real = %f [ms]\n", 1.E3*timer.CpuTime(), 1.E3*timer.RealTime());
+       printf("Points indexes in increasing order of their distance to the target point.\n");
+       for(int i=0; i<kNN; i++){
+               printf("%5d [%7.4f] ", fkNN[i], fkNNdist[i]);
+               if(i%5==4) printf("\n");
+       }
+       printf("\n");   
+
+       TGraph *gSA = new TGraph(kNN);
+       gSA->SetMarkerStyle(24);
+       for(int ip=0; ip<kNN; ip++) gSA->SetPoint(ip, x[0][fkNN[ip]], x[1][fkNN[ip]]);
+
+       gSA->Draw("p"); leg->AddEntry(gSA, "Stand Alone", "p");
+       leg->Draw();
+       
+       for(int ip=0; ip<ntimes; ip++) delete [] p[ip];
+       delete [] p;
+       for(int idim=0; idim<ndim; idim++) delete [] x[idim];
+       delete [] x;
+}
+
+
+//______________________________________________________________
+Float_t p[]={1.4, -.6}; //}{1.7, -.4};
+void kNNDraw(const Float_t *p, const int kNN)
+{
+// Test memory consumption and draw "kNN" nearest neighbours of point "p".
+// The distance (in the L1 metric) is encoded in the color code.
+
+       const Int_t npoints = 10000;
+       Float_t *data0 =  new Float_t[npoints*2];
+  Float_t *data[2];
+  data[0] = &data0[0];
+  data[1] = &data0[npoints];
+  for (Int_t i=0;i<npoints;i++) {
+    data[1][i]= gRandom->Gaus();
+    data[0][i]= gRandom->Gaus();
+  }
+
+       TKDPDF pdf(npoints, 2, 100, data);
+       pdf.DrawNode(pdf.GetNodeIndex(p));
+       
+       TMarker *mp = new TMarker(p[0], p[1], 3);
+       mp->SetMarkerColor(4);
+       mp->Draw();
+               
+       Int_t *index, color;
+       Float_t *d, d0, pknn[2];
+       Float_t start = Mem();
+       pdf.FindNearestNeighbors(p, kNN, index, d);
+       Float_t end = Mem();
+       printf("FindNearestNeighbors memory usage %fKB\n", end-start);
+       d0 = d[kNN-1];
+       TMarker *ma = new TMarker[kNN];
+       for(int ip=0; ip<kNN; ip++){
+               pdf.GetDataPoint(index[ip], pknn);
+               color = 101 - Int_t(50. * d[ip]/d0);
+               ma[ip].SetMarkerStyle(4);
+               ma[ip].SetMarkerColor(color);
+               ma[ip].DrawMarker(pknn[0], pknn[1]);
+       }
+}
+
+
+//______________________________________________________________
+Float_t Mem()
+{
+  // size in KB
+  static struct mallinfo memdebug; 
+  memdebug = mallinfo();
+  return memdebug.uordblks/1024.;
+}
+
+//______________________________________________________________
+void build(const Int_t ndim, const Int_t nstat)
+{
+// Build "nstat" data points in "ndim" dimensions taken from an
+// uncorrelated 2D Gauss distribution.
+       
+
+       printf("build data ... \n");
+       Double_t pntTrain[ndim];
+       TFile *f = TFile::Open(Form("%dD_Gauss.root", ndim), "RECREATE");
+       TTree *t = new TTree("db", "gauss database");
+       for(int idim=0; idim<ndim; idim++) t->Branch(Form("x%d", idim), &pntTrain[idim], Form("x%d/D", idim));
+       
+       for(int istat=0; istat<nstat; istat++){
+               for(int idim=0; idim<ndim; idim++) pntTrain[idim] = gRandom->Gaus();
+               t->Fill();
+       }
+       f->cd();
+       t->Write();
+       f->Close();
+       delete f;
+}
+
diff --git a/STAT/Macros/pdfCompare.C b/STAT/Macros/pdfCompare.C
new file mode 100644 (file)
index 0000000..75673b5
--- /dev/null
@@ -0,0 +1,75 @@
+void pdfCompare(const Int_t method = 0, const int nstat = 40000)
+{
+// Compare interpolation methods.
+// method = 0 : Select COG method
+// method = 1 : Select INT method
+
+       const int bucket = nstat/100; //400;
+       const int nevals = 100;
+       Float_t *data[1];
+       Float_t data1[nstat];
+       for(int istat=0; istat<nstat; istat++) data1[istat] = gRandom->Gaus();
+       data[0] = &data1[0];
+
+       TLegend *leg = new TLegend(0.7157258,0.8280255,0.9637097,0.9596603);
+       leg->SetBorderSize(1);
+       TF1 *f=new TF1("f1", "gaus(0);x;f(x)", -5., 5.);
+       f->SetLineWidth(1);
+       f->SetLineColor(2);
+       f->SetParameter(0, 1.);
+       f->SetParameter(1, 0.);
+       f->SetParameter(2, 1.);
+       f->SetParameter(0, 1./f->Integral(-10., 10.));
+       f->Draw();      leg->AddEntry(f, "model", "l");
+       
+       
+       Float_t cog, val, val_err;
+       Double_t x, chi2, CHI2, result, err;
+       Int_t npoints;
+       TKDPDF pdf(nstat, 1, bucket, data);
+       pdf.SetInterpolationMethod(method);
+       pdf.SetStore();
+       pdf.SetAlpha(1.);
+       Float_t *bounds;// = in.GetBoundary(0);
+       Double_t dx = 10./nevals;
+
+       TGraph *gRE = new TGraph(npoints);
+       gRE->SetMarkerStyle(4);
+       gRE->SetMarkerColor(4);
+       gRE->SetLineColor(4);
+       
+       // do a preevaluation for INT method on COG points
+       Float_t *pcog;
+       for(int icog=0; icog<pdf.GetNTNodes(); icog++){
+               pdf.GetCOGPoint(icog, pcog, val, val_err);
+               x = pcog[0];
+               chi2 = pdf.Eval(&x, result, err);
+       }
+       //pdf.GetStatus();
+
+       TGraphErrors *gINT = new TGraphErrors(npoints);
+       gINT->SetMarkerStyle(7);
+       x = 5.-dx/2.;
+       CHI2 = 0.; npoints = 0;
+       for(int ip=0; ip<nevals; ip++){
+               chi2 = pdf.Eval(&x, result, err);
+               gINT->SetPoint(ip, x, result);
+               gINT->SetPointError(ip, 0., err);
+               gRE->SetPoint(ip, x, .01*err/TMath::Abs(result));
+               if(TMath::Abs(x)<2.){
+                       CHI2 += (result - f->Eval(x))*(result - f->Eval(x))/err/err;
+                       npoints++;
+               }
+               x -= dx;
+       }
+       gINT->Draw("p"); leg->AddEntry(gINT, Form("interpolation [%s]", method ? "INT" : "COG"), "pl");
+       
+       gRE->Draw("pl"); leg->AddEntry(gRE, "1% relative errors", "pl");
+       CHI2 /= npoints;
+       printf("chi2(%d) = %f\n", npoints, CHI2);
+
+       leg->Draw();
+}
+
+
+
diff --git a/STAT/Macros/pdfIO.C b/STAT/Macros/pdfIO.C
new file mode 100644 (file)
index 0000000..9b80018
--- /dev/null
@@ -0,0 +1,110 @@
+TStopwatch timer;
+const Int_t nt = 100;
+//___________________________________________________________
+void pdfIO(const Int_t method = 0)
+{
+// Test interpolator IO response for several data dimensions.
+// method = 0 : use COG interpolator
+// method = 1 : use INT interpolator
+
+       Float_t tw, tr;
+       TGraph *gw = new TGraph(5);
+       gw->SetMarkerStyle(24);gw->SetMarkerColor(2);
+       TGraph *gr = new TGraph(5);
+       gr->SetMarkerStyle(25);gr->SetMarkerColor(4);
+       for(int idim = 1; idim<6; idim++){
+               tw = interpolWrite(method, idim);
+               gw->SetPoint(idim-1, idim, tw);
+               tr = interpolRead(idim);
+               gr->SetPoint(idim-1, idim, tr);
+       }
+       gw->Draw("apl");
+       gr->Draw("pl");
+}
+
+//___________________________________________________________
+Float_t interpolWrite(const Int_t method, const Int_t ndim)
+{
+       if(!gSystem->FindFile(".", "5D_Gauss.root")) build(5, 1000000);
+       TFile::Open("5D_Gauss.root");
+       printf("\nWriting %dD interpolator ...\n", ndim);
+       TTree *db = (TTree*)gFile->Get("db");
+       TString var = "x0"; for(int idim=1; idim<ndim; idim++) var += Form(":x%d", idim);
+       printf("\tBuilding interpolator ...\n");
+       TKDPDF pdf(db, var.Data(), "", 400, 100000);
+       pdf.SetInterpolationMethod(method);
+       pdf.SetStore();
+
+       
+       printf("\tSetting up the interpolator ...\n");
+       Float_t *c, val, v_err;
+       Double_t *p = new Double_t[ndim], res, r_err;
+       timer.Start(kTRUE);
+       for(int inode=0; inode<pdf.GetNTNodes(); inode++){
+               pdf.GetCOGPoint(inode, c, val, v_err);
+               for(int idim=0; idim<ndim; idim++) p[idim] = (Double_t)c[idim];
+//             printf("Evaluate for node [%d] ", inode);
+//             for(int idim=0; idim<ndim; idim++) printf("%5.3f ", p[idim]); printf("\n");
+               for(int it=0; it<nt; it++) pdf.Eval(p, res, r_err, kTRUE);
+//             printf("R = %6.4f +- %6.4f\n", res, r_err);
+       }
+       timer.Stop();
+       
+       Float_t time = 1.E3*timer.CpuTime()/float(nt);
+       printf("\tFinish interpolation in %6.2f [ms]\n", time);
+       
+       printf("\tSaving interpolator ...\n");
+       TFile *fi = TFile::Open(Form("%dD_Interpolator.root", ndim), "RECREATE");
+       pdf.Write(Form("%dDgauss", ndim));
+       fi->Close();
+       delete fi;
+       
+       delete [] p;
+
+       return time;
+}
+
+//___________________________________________________________
+Float_t interpolRead(const Int_t ndim)
+{
+       printf("Reading %dD interpolator ...\n", ndim);
+       TFile::Open(Form("%dD_Interpolator.root", ndim));
+       TKDPDF *pdf = (TKDPDF*)gFile->Get(Form("%dDgauss", ndim));
+       gFile->Close();
+       
+       printf("\tDoing interpolation ...\n");
+       Float_t *c, v, ve;
+       Double_t *p = new Double_t[ndim], r, re;
+       timer.Start(kTRUE);     
+       for(int ip=0; ip<pdf->GetNTNodes(); ip++){
+               pdf->GetCOGPoint(ip, c, v, ve);
+               for(int idim=0; idim<ndim; idim++) p[idim] = (Double_t)c[idim];
+               for(int it=0; it<nt; it++) pdf->Eval(p, r, re);
+       }
+       timer.Stop();
+       Float_t time = 1.E3*timer.CpuTime()/float(nt);
+       printf("\tFinish interpolation in %6.2f [ms]\n", time);
+       delete [] p;
+       return time;
+}
+
+//___________________________________________________________
+void build(const int ndim, const int npoints)
+{
+       printf("Building DB ...\n");
+       Float_t data[ndim];
+       TFile::Open(Form("%dD_Gauss.root", ndim), "RECREATE");
+       TTree *t = new TTree("db", Form("%dD data base for kD statistics", ndim));
+       for(int idim=0; idim<ndim; idim++) t->Branch(Form("x%d", idim), &data[idim], Form("x%d/F", idim));
+
+       for (Int_t ip=0; ip<npoints; ip++){
+               for (Int_t id=0; id<ndim; id++) data[id]= gRandom->Gaus();
+               t->Fill();
+       }
+
+       t->Write();
+       gFile->Close();
+       delete gFile;
+}
+
+
diff --git a/STAT/Macros/pdfTest.C b/STAT/Macros/pdfTest.C
new file mode 100644 (file)
index 0000000..318a7c1
--- /dev/null
@@ -0,0 +1,219 @@
+#include <malloc.h>
+
+#include "TSystem.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TProfile.h"
+#include "TF1.h"
+#include "TF2.h"
+#include "TF3.h"
+#include "TLegend.h"
+#include "TRandom.h"
+#include "TString.h"
+#include "TGraph.h"
+#include "TMarker.h"
+#include "TStopwatch.h"
+#include "../src/TKDPDF.h"
+
+void pdfTest(const Int_t ndim = 1);
+Double_t interpolation(const int nstat, const int ndim, const int first);
+void build(const int ndim, const int npoints);
+Float_t Mem();
+
+TF1 *f = 0x0;
+
+//__________________________________________________________
+void pdfTest(const Int_t ndim)
+{
+// Test interpolator on a variable dimension (ndim<=5) uncorrelated
+// Gauss model. The macro displays the chi2 between interpolation and
+// model for various statistics of experimental data.
+//
+// Check the worker function interpolation() to fine tune the
+// algorithm.
+// 
+// Attention:
+// The macro works in compiled mode !
+
+       // define model
+       switch(ndim){
+       case 1:
+               f=new TF1("f1", "gaus(0);x;f(x)", -5., 5.);
+               f->SetParameter(0, 1.);
+               f->SetParameter(1, 0.);
+               f->SetParameter(2, 1.);
+               f->SetParameter(0, 1./f->Integral(-5., 5.));
+               break;
+       case 2:
+               f=new TF2("f2", "[0]*exp(-0.5*((x-[1])/[2])**2)*exp(-0.5*((y-[3])/[4])**2)", -5., 5., -5., 5.);
+               f->SetParameter(0, 1.);
+               f->SetParameter(1, 0.);
+               f->SetParameter(2, 1.);
+               f->SetParameter(3, 0.);
+               f->SetParameter(4, 1.);
+               f->SetParameter(0, 1./f->Integral(-5., 5., -5., 5.));
+               break;
+       case 3:
+       case 4:
+       case 5:
+               f=new TF3("f3", "[0]*exp(-0.5*((x-[1])/[2])**2)*exp(-0.5*((y-[3])/[4])**2)*exp(-0.5*((z-[5])/[6])**2)", -5., 5., -5., 5., -5., 5.);
+               f->SetParameter(0, 1.);
+               f->SetParameter(1, 0.);
+               f->SetParameter(2, 1.);
+               f->SetParameter(3, 0.);
+               f->SetParameter(4, 1.);
+               f->SetParameter(5, 0.);
+               f->SetParameter(6, 1.);
+               f->SetParameter(0, 1./f->Integral(-5., 5., -5., 5., -5., 5.));
+               if(ndim>3) printf("Warning : chi2 results are unreliable !\n");
+               break;
+       default:
+               printf("%dD not supported in this macro.\n", ndim);
+               return;
+       }
+
+       Double_t chi2;
+       const Int_t nchi2 = 18;
+       const Int_t nfirst = 1;
+       //TProfile *hp = new TProfile("hp", "", 100, 1.E4, 1.E6);
+       //hp->SetMarkerStyle(24);
+       TGraph *gChi2 = new TGraph(nchi2);
+       gChi2->SetMarkerStyle(24);
+       Int_t ip = 0, nstat, first;
+       for(int ifirst=0; ifirst<nfirst; ifirst++){
+               first = 0;//Int_t(gRandom->Uniform(1.E4));
+               for(int i=2; i<10; i++){
+                       nstat = 10000*i;
+                       chi2 = interpolation(nstat, ndim, first);
+                       gChi2->SetPoint(ip++, TMath::Log10(float(nstat)), chi2);
+                       //hp->Fill(float(nstat), chi2);
+               }
+               for(int i=1; i<10; i++){
+                       nstat = 100000*i;
+                       chi2 = interpolation(nstat, ndim, first);
+                       gChi2->SetPoint(ip++, TMath::Log10(float(nstat)), chi2);
+                       //hp->Fill(float(nstat), chi2);
+               }
+               gChi2->Draw("apl");
+       }
+       //hp->Draw("pl");
+}
+
+
+//__________________________________________________________
+Double_t interpolation(const int nstat, const int ndim, const int first)
+{
+// Steer interpolation of a nD (n<=5) uncorrelated Gauss distribution.
+// The user is suppose to give the experimental statistics (nstat) the
+// dimension of the experimental point (ndim). The entry from which the
+// data base is being read is steered by "first".
+
+       const int bs = 400;
+       const int npoints = 100;
+       // switch on chi2 calculation between interpolation and model.
+       // Default calculates distance.
+       const Bool_t kCHI2 = kFALSE; 
+       const Bool_t kINT  = kTRUE; // switch on integral interpolator
+
+       // open data base
+       TString fn("5D_Gauss.root");
+       if(!gSystem->FindFile(".", fn)) build(5, 1000000);
+       TFile::Open("5D_Gauss.root");
+       TTree *t = (TTree*)gFile->Get("db");
+
+       // build interpolator
+       TString var = "x0"; for(int idim=1; idim<ndim; idim++) var += Form(":x%d", idim);
+       TKDPDF pdf(t, var.Data(), "", bs, nstat, first);
+       pdf.SetInterpolationMethod(kINT);
+       pdf.SetStore();
+       
+       // define testing grid
+       Double_t x[5];
+       Double_t dx[5] = {4., 4., 4., 4., 4.};
+       Double_t chi2, theor, result, err;
+       for(int idim=0; idim<ndim; idim++){
+               dx[idim] = (ndim<<2)/float(npoints);
+               x[idim]  = -2.+dx[idim]/2.;
+       }
+       
+       // setup the integral interpolator and do memory test
+       Float_t start, stop;
+       start = Mem();
+       Float_t *c, v, ve;
+       for(int ip=0; ip<pdf.GetNTNodes(); ip++){
+               pdf.GetCOGPoint(ip, c, v, ve);
+               for(int idim=0; idim<ndim; idim++) x[idim] = (Double_t)c[idim];
+               pdf.Eval(x, result, err);
+       }
+       stop  = Mem();
+       printf("\nInterpolating %d data points in %dD ...\n", nstat, ndim);
+       printf("\tbucket size %d\n", bs);
+       printf("\tmemory usage in Eval() %6.4fKB\n", stop-start);
+       printf("\tstarting interpolation ...\n");
+       //return stop-start;
+
+       // evaluate relative error
+       for(int idim=0; idim<ndim; idim++) x[idim] = 0.;
+       pdf.Eval(x, result, err);
+       Double_t threshold = err/result;
+       
+       // starting interpolation over the testing grid
+       chi2 = 0.; Int_t nsample = 0;
+       x[4]  = -2.+dx[4]/2.;
+       do{
+               x[3]  = -2.+dx[3]/2.;
+               do{
+                       x[2]  = -2.+dx[2]/2.;
+                       do{
+                               x[1]  = -2.+dx[1]/2.;
+                               do{
+                                       x[0]  = -2.+dx[0]/2.;
+                                       do{
+                                               pdf.Eval(x, result, err);
+                                               if(err/TMath::Abs(result) < 1.1*threshold){
+                                                       theor = f->Eval(x[0], x[1], x[2]);
+                                                       Double_t tmp = result - theor; tmp *= tmp;
+                                                       chi2 += tmp/(kCHI2 ? err*err : theor);
+                                                       nsample++;
+                                               }
+                                       } while((x[0] += dx[0]) < 2.);
+                               } while((x[1] += dx[1]) < 2.);
+                       } while((x[2] += dx[2]) < 2.);
+               } while((x[3] += dx[3]) < 2.);
+       } while((x[4] += dx[4]) < 2.);
+       chi2 /= float(nsample);
+       printf("\tinterpolation quality (chi2) %f\n", chi2);
+       
+       return kCHI2 ? chi2 : TMath::Sqrt(chi2);
+}
+
+
+//___________________________________________________________
+void build(const int ndim, const int npoints)
+{
+       printf("Building DB ...\n");
+       Float_t data[ndim];
+       TFile::Open(Form("%dD_Gauss.root", ndim), "RECREATE");
+       TTree *t = new TTree("db", Form("%dD data base for kD statistics", ndim));
+       for(int idim=0; idim<ndim; idim++) t->Branch(Form("x%d", idim), &data[idim], Form("x%d/F", idim));
+
+       for (Int_t ip=0; ip<npoints; ip++){
+               for (Int_t id=0; id<ndim; id++) data[id]= gRandom->Gaus();
+               t->Fill();
+       }
+
+       t->Write();
+       gFile->Close();
+       delete gFile;
+}
+
+//______________________________________________________________
+Float_t Mem()
+{
+  // size in KB
+  static struct mallinfo debug;
+  debug = mallinfo();
+       //printf("arena[%d] usmblks[%d] uordblks[%d] fordblks[%d]\n", debug.arena, debug.usmblks, debug.uordblks, debug.fordblks);
+  return debug.uordblks/1024.;
+}
+