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