new macros for testing the STAT package
[u/mrichter/AliRoot.git] / STAT / Macros / pdfTest.C
CommitLineData
8a8f8274 1#include <malloc.h>
2
3#include "TSystem.h"
4#include "TFile.h"
5#include "TTree.h"
6#include "TProfile.h"
7#include "TF1.h"
8#include "TF2.h"
9#include "TF3.h"
10#include "TLegend.h"
11#include "TRandom.h"
12#include "TString.h"
13#include "TGraph.h"
14#include "TMarker.h"
15#include "TStopwatch.h"
16#include "../src/TKDPDF.h"
17
18void pdfTest(const Int_t ndim = 1);
19Double_t interpolation(const int nstat, const int ndim, const int first);
20void build(const int ndim, const int npoints);
21Float_t Mem();
22
23TF1 *f = 0x0;
24
25//__________________________________________________________
26void pdfTest(const Int_t ndim)
27{
28// Test interpolator on a variable dimension (ndim<=5) uncorrelated
29// Gauss model. The macro displays the chi2 between interpolation and
30// model for various statistics of experimental data.
31//
32// Check the worker function interpolation() to fine tune the
33// algorithm.
34//
35// Attention:
36// The macro works in compiled mode !
37
38 // define model
39 switch(ndim){
40 case 1:
41 f=new TF1("f1", "gaus(0);x;f(x)", -5., 5.);
42 f->SetParameter(0, 1.);
43 f->SetParameter(1, 0.);
44 f->SetParameter(2, 1.);
45 f->SetParameter(0, 1./f->Integral(-5., 5.));
46 break;
47 case 2:
48 f=new TF2("f2", "[0]*exp(-0.5*((x-[1])/[2])**2)*exp(-0.5*((y-[3])/[4])**2)", -5., 5., -5., 5.);
49 f->SetParameter(0, 1.);
50 f->SetParameter(1, 0.);
51 f->SetParameter(2, 1.);
52 f->SetParameter(3, 0.);
53 f->SetParameter(4, 1.);
54 f->SetParameter(0, 1./f->Integral(-5., 5., -5., 5.));
55 break;
56 case 3:
57 case 4:
58 case 5:
59 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.);
60 f->SetParameter(0, 1.);
61 f->SetParameter(1, 0.);
62 f->SetParameter(2, 1.);
63 f->SetParameter(3, 0.);
64 f->SetParameter(4, 1.);
65 f->SetParameter(5, 0.);
66 f->SetParameter(6, 1.);
67 f->SetParameter(0, 1./f->Integral(-5., 5., -5., 5., -5., 5.));
68 if(ndim>3) printf("Warning : chi2 results are unreliable !\n");
69 break;
70 default:
71 printf("%dD not supported in this macro.\n", ndim);
72 return;
73 }
74
75 Double_t chi2;
76 const Int_t nchi2 = 18;
77 const Int_t nfirst = 1;
78 //TProfile *hp = new TProfile("hp", "", 100, 1.E4, 1.E6);
79 //hp->SetMarkerStyle(24);
80 TGraph *gChi2 = new TGraph(nchi2);
81 gChi2->SetMarkerStyle(24);
82 Int_t ip = 0, nstat, first;
83 for(int ifirst=0; ifirst<nfirst; ifirst++){
84 first = 0;//Int_t(gRandom->Uniform(1.E4));
85 for(int i=2; i<10; i++){
86 nstat = 10000*i;
87 chi2 = interpolation(nstat, ndim, first);
88 gChi2->SetPoint(ip++, TMath::Log10(float(nstat)), chi2);
89 //hp->Fill(float(nstat), chi2);
90 }
91 for(int i=1; i<10; i++){
92 nstat = 100000*i;
93 chi2 = interpolation(nstat, ndim, first);
94 gChi2->SetPoint(ip++, TMath::Log10(float(nstat)), chi2);
95 //hp->Fill(float(nstat), chi2);
96 }
97 gChi2->Draw("apl");
98 }
99 //hp->Draw("pl");
100}
101
102
103//__________________________________________________________
104Double_t interpolation(const int nstat, const int ndim, const int first)
105{
106// Steer interpolation of a nD (n<=5) uncorrelated Gauss distribution.
107// The user is suppose to give the experimental statistics (nstat) the
108// dimension of the experimental point (ndim). The entry from which the
109// data base is being read is steered by "first".
110
111 const int bs = 400;
112 const int npoints = 100;
113 // switch on chi2 calculation between interpolation and model.
114 // Default calculates distance.
115 const Bool_t kCHI2 = kFALSE;
116 const Bool_t kINT = kTRUE; // switch on integral interpolator
117
118 // open data base
119 TString fn("5D_Gauss.root");
120 if(!gSystem->FindFile(".", fn)) build(5, 1000000);
121 TFile::Open("5D_Gauss.root");
122 TTree *t = (TTree*)gFile->Get("db");
123
124 // build interpolator
125 TString var = "x0"; for(int idim=1; idim<ndim; idim++) var += Form(":x%d", idim);
126 TKDPDF pdf(t, var.Data(), "", bs, nstat, first);
127 pdf.SetInterpolationMethod(kINT);
128 pdf.SetStore();
129
130 // define testing grid
131 Double_t x[5];
132 Double_t dx[5] = {4., 4., 4., 4., 4.};
133 Double_t chi2, theor, result, err;
134 for(int idim=0; idim<ndim; idim++){
135 dx[idim] = (ndim<<2)/float(npoints);
136 x[idim] = -2.+dx[idim]/2.;
137 }
138
139 // setup the integral interpolator and do memory test
140 Float_t start, stop;
141 start = Mem();
142 Float_t *c, v, ve;
143 for(int ip=0; ip<pdf.GetNTNodes(); ip++){
144 pdf.GetCOGPoint(ip, c, v, ve);
145 for(int idim=0; idim<ndim; idim++) x[idim] = (Double_t)c[idim];
146 pdf.Eval(x, result, err);
147 }
148 stop = Mem();
149 printf("\nInterpolating %d data points in %dD ...\n", nstat, ndim);
150 printf("\tbucket size %d\n", bs);
151 printf("\tmemory usage in Eval() %6.4fKB\n", stop-start);
152 printf("\tstarting interpolation ...\n");
153 //return stop-start;
154
155 // evaluate relative error
156 for(int idim=0; idim<ndim; idim++) x[idim] = 0.;
157 pdf.Eval(x, result, err);
158 Double_t threshold = err/result;
159
160 // starting interpolation over the testing grid
161 chi2 = 0.; Int_t nsample = 0;
162 x[4] = -2.+dx[4]/2.;
163 do{
164 x[3] = -2.+dx[3]/2.;
165 do{
166 x[2] = -2.+dx[2]/2.;
167 do{
168 x[1] = -2.+dx[1]/2.;
169 do{
170 x[0] = -2.+dx[0]/2.;
171 do{
172 pdf.Eval(x, result, err);
173 if(err/TMath::Abs(result) < 1.1*threshold){
174 theor = f->Eval(x[0], x[1], x[2]);
175 Double_t tmp = result - theor; tmp *= tmp;
176 chi2 += tmp/(kCHI2 ? err*err : theor);
177 nsample++;
178 }
179 } while((x[0] += dx[0]) < 2.);
180 } while((x[1] += dx[1]) < 2.);
181 } while((x[2] += dx[2]) < 2.);
182 } while((x[3] += dx[3]) < 2.);
183 } while((x[4] += dx[4]) < 2.);
184 chi2 /= float(nsample);
185 printf("\tinterpolation quality (chi2) %f\n", chi2);
186
187 return kCHI2 ? chi2 : TMath::Sqrt(chi2);
188}
189
190
191//___________________________________________________________
192void build(const int ndim, const int npoints)
193{
194 printf("Building DB ...\n");
195 Float_t data[ndim];
196 TFile::Open(Form("%dD_Gauss.root", ndim), "RECREATE");
197 TTree *t = new TTree("db", Form("%dD data base for kD statistics", ndim));
198 for(int idim=0; idim<ndim; idim++) t->Branch(Form("x%d", idim), &data[idim], Form("x%d/F", idim));
199
200 for (Int_t ip=0; ip<npoints; ip++){
201 for (Int_t id=0; id<ndim; id++) data[id]= gRandom->Gaus();
202 t->Fill();
203 }
204
205 t->Write();
206 gFile->Close();
207 delete gFile;
208}
209
210//______________________________________________________________
211Float_t Mem()
212{
213 // size in KB
214 static struct mallinfo debug;
215 debug = mallinfo();
216 //printf("arena[%d] usmblks[%d] uordblks[%d] fordblks[%d]\n", debug.arena, debug.usmblks, debug.uordblks, debug.fordblks);
217 return debug.uordblks/1024.;
218}
219