]>
Commit | Line | Data |
---|---|---|
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 | ||
18 | void pdfTest(const Int_t ndim = 1); | |
19 | Double_t interpolation(const int nstat, const int ndim, const int first); | |
20 | void build(const int ndim, const int npoints); | |
21 | Float_t Mem(); | |
22 | ||
23 | TF1 *f = 0x0; | |
24 | ||
25 | //__________________________________________________________ | |
26 | void 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 | //__________________________________________________________ | |
104 | Double_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 | //___________________________________________________________ | |
192 | void 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 | //______________________________________________________________ | |
211 | Float_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 |