]>
Commit | Line | Data |
---|---|---|
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 | ||
17 | void interpolTest(const Int_t ndim = 1); | |
18 | Double_t interpolation(const int nstat, const int ndim, const int first); | |
19 | void build(const int ndim, const int npoints); | |
20 | ||
21 | TF1 *f = 0x0; | |
22 | ||
23 | //__________________________________________________________ | |
24 | void 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 | ||
100 | void test() | |
101 | { | |
102 | TClonesArray ar("TKDNodeInfo", 10); | |
103 | } | |
104 | ||
105 | //__________________________________________________________ | |
106 | Double_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 | //___________________________________________________________ | |
225 | void 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 | } |