]>
Commit | Line | Data |
---|---|---|
8a8f8274 | 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 |