AliRoot Core  3dc7879 (3dc7879)
interpolTest.C
Go to the documentation of this file.
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 }
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
TFile * Open(const char *filename, Long64_t &nevents)
void interpolTest(const Int_t ndim=1)
Definition: interpolTest.C:24
Definition: TKDPDF.h:20
npoints
Definition: driftITSTPC.C:85
Double_t interpolation(const int nstat, const int ndim, const int first)
Definition: interpolTest.C:106
Double_t chi2
Definition: AnalyzeLaser.C:7
TF1 * f
Definition: interpolTest.C:21
void test()
Definition: interpolTest.C:100
void build(const int ndim, const int npoints)
Definition: interpolTest.C:225