1. Removed obsolete test functions
[u/mrichter/AliRoot.git] / STAT / Macros / interpolCompare.C
CommitLineData
8a8f8274 1void interpolCompare(const Int_t method = 0, const int nstat = 40000)
2{
3// Compare interpolation methods.
4// method = 0 : Select COG method
5// method = 1 : Select INT method
6
7 const Float_t alpha = 20.;
8 const int bucket = nstat/100; //400;
9 const int nevals = 100;
10 Float_t *data[1];
11 Float_t data1[nstat];
12 for(int istat=0; istat<nstat; istat++) data1[istat] = gRandom->Gaus();
13 data[0] = &data1[0];
14
15 TLegend *leg = new TLegend(0.7157258,0.8280255,0.9637097,0.9596603);
16 leg->SetBorderSize(1);
17 TF1 *f=new TF1("f1", "gaus(0);x;f(x)", -5., 5.);
18 f->SetLineWidth(1);
19 f->SetLineColor(2);
20 f->SetParameter(0, 1.);
21 f->SetParameter(1, 0.);
22 f->SetParameter(2, 1.);
23 f->SetParameter(0, 1./f->Integral(-10., 10.));
24 f->Draw(); leg->AddEntry(f, "model", "l");
25
26
27 Float_t cog, val, val_err;
28 Double_t x, chi2, chi2pdf, chi2int, result, err;
29 Int_t npPDF, npINT;
30 TKDPDF pdf(nstat, 1, bucket, data);
31 pdf.SetInterpolationMethod(method);
32 pdf.SetAlpha(alpha);
33
34 TKDInterpolator in(1, pdf.GetNTNodes());
35 in.SetInterpolationMethod(method);
36 in.SetAlpha(alpha);
37
38 // change to the reprezentation which is closer to the local
39 // interpolating polynom.
40 Double_t fx;
41 Int_t inode = 0;
42 TKDNodeInfo *node = 0x0;
43 while(node = pdf.GetNodeInfo(inode)){
44 if(node->fVal[0] > 0.){
45 in.SetNode(inode, *node);
46 node = in.GetNodeInfo(inode);
47 fx = node->fVal[0];
48 node->fVal[0] = TMath::Log(fx);
49 node->fVal[1] = node->fVal[1]/fx;
50 }
51 inode++;
52 }
53
54 // Do preprocessing for INT
55 Float_t *pcog;
56 if(method) for(int icog=0; icog<pdf.GetNTNodes(); icog++){
57 pdf.GetCOGPoint(icog, pcog, val, val_err);
58 x = pcog[0];
59 pdf.Eval(&x, result, err);
60 in.Eval(&x, result, err);
61 }
62
63
64
65 // prepare drawing containers
66 TGraphErrors *gPDF = new TGraphErrors(nevals);
67 gPDF->SetMarkerStyle(7);
68 TGraph *gREP = new TGraph(nevals);
69 gREP->SetMarkerStyle(7);
70 gREP->SetMarkerColor(4);
71 gREP->SetLineColor(4);
72
73 TGraphErrors *gINT = new TGraphErrors(nevals);
74 gINT->SetMarkerStyle(24);
75 gINT->SetMarkerSize(.6);
76 gINT->SetMarkerColor(3);
77 gINT->SetLineColor(3);
78 TGraph *gREI = new TGraph(nevals);
79 gREI->SetMarkerStyle(24);
80 gREI->SetMarkerColor(4);
81 gREI->SetLineColor(4);
82
83 // Do the interpolation with PDF and Interpolator
84 Double_t dx = 10./nevals;
85 x = 5.-dx/2.;
86 chi2pdf = 0.; npPDF = 0;
87 chi2int = 0.; npINT = 0;
88 for(int ip=0; ip<nevals; ip++){
89 chi2 = pdf.Eval(&x, result, err);
90 gPDF->SetPoint(ip, x, result);
91 gPDF->SetPointError(ip, 0., err);
92 if(TMath::Abs(x)<3.){
93 chi2pdf += (result - f->Eval(x))*(result - f->Eval(x))/f->Eval(x)/f->Eval(x);
94 npPDF++;
95 }
96 gREP->SetPoint(ip, x, .1*err/TMath::Abs(result));
97 //printf("%2d x[%f] r[%f +- %f] chi2[%e]\n", ip, x, result, err, chi2);
98
99 chi2 = in.Eval(&x, result, err);
100 result = TMath::Exp(result);
101 err = err*result;
102 gINT->SetPoint(ip, x, result);
103 gINT->SetPointError(ip, 0., err);
104 if(TMath::Abs(x)<3.){
105 chi2int += (result - f->Eval(x))*(result - f->Eval(x))/f->Eval(x)/f->Eval(x);
106 npINT++;
107 }
108 gREI->SetPoint(ip, x, .1*err/TMath::Abs(result));
109 //printf("%2d x[%f] r[%f +- %f] chi2[%e]\n\n", ip, x, result, err, chi2);
110
111 x -= dx;
112 }
113 gINT->Draw("p"); leg->AddEntry(gINT, Form("Interpolator [%s]", method ? "INT" : "COG"), "pl");
114 gPDF->Draw("p"); leg->AddEntry(gPDF, Form("PDF [%s]", method ? "INT" : "COG"), "pl");
115 gREI->Draw("pl"); leg->AddEntry(gREI, "1% #epsilon Interpolator", "pl");
116 gREP->Draw("pl"); leg->AddEntry(gREP, "1% #epsilon PDF", "pl");
117
118 chi2pdf /= npPDF;
119 printf("PDF : chi2(%d) = %f\n", npPDF, chi2pdf);
120 chi2int /= npINT;
121 printf("INT : chi2(%d) = %f\n", npINT, chi2int);
122
123 leg->Draw();
124}
125
126
127