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