--- /dev/null
+#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.;
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff