new macros for testing the STAT package
[u/mrichter/AliRoot.git] / STAT / Macros / interpolCompare.C
1 void 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