STAT/Macros/pdfCompare.C

   1 void pdfCompare(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 int bucket = nstat/100; //400;
   8         const int nevals = 100;
   9         Float_t *data[1];
  10         Float_t data1[nstat];
  11         for(int istat=0; istat<nstat; istat++) data1[istat] = gRandom->Gaus();
  12         data[0] = &data1[0];
  13
  14         TLegend *leg = new TLegend(0.7157258,0.8280255,0.9637097,0.9596603);
  15         leg->SetBorderSize(1);
  16         TF1 *f=new TF1("f1", "gaus(0);x;f(x)", -5., 5.);
  17         f->SetLineWidth(1);
  18         f->SetLineColor(2);
  19         f->SetParameter(0, 1.);
  20         f->SetParameter(1, 0.);
  21         f->SetParameter(2, 1.);
  22         f->SetParameter(0, 1./f->Integral(-10., 10.));
  23         f->Draw();      leg->AddEntry(f, "model", "l");
  24
  25
  26         Float_t cog, val, val_err;
  27         Double_t x, chi2, CHI2, result, err;
  28         Int_t npoints;
  29         TKDPDF pdf(nstat, 1, bucket, data);
  30         pdf.SetInterpolationMethod(method);
  31         pdf.SetStore();
  32         pdf.SetAlpha(1.);
  33         Float_t *bounds;// = in.GetBoundary(0);
  34         Double_t dx = 10./nevals;
  35
  36         TGraph *gRE = new TGraph(npoints);
  37         gRE->SetMarkerStyle(4);
  38         gRE->SetMarkerColor(4);
  39         gRE->SetLineColor(4);
  40
  41         // do a preevaluation for INT method on COG points
  42         Float_t *pcog;
  43         for(int icog=0; icog<pdf.GetNTNodes(); icog++){
  44                 pdf.GetCOGPoint(icog, pcog, val, val_err);
  45                 x = pcog[0];
  46                 chi2 = pdf.Eval(&x, result, err);
  47         }
  48         //pdf.GetStatus();
  49
  50         TGraphErrors *gINT = new TGraphErrors(npoints);
  51         gINT->SetMarkerStyle(7);
  52         x = 5.-dx/2.;
  53         CHI2 = 0.; npoints = 0;
  54         for(int ip=0; ip<nevals; ip++){
  55                 chi2 = pdf.Eval(&x, result, err);
  56                 gINT->SetPoint(ip, x, result);
  57                 gINT->SetPointError(ip, 0., err);
  58                 gRE->SetPoint(ip, x, .01*err/TMath::Abs(result));
  59                 if(TMath::Abs(x)<2.){
  60                         CHI2 += (result - f->Eval(x))*(result - f->Eval(x))/err/err;
  61                         npoints++;
  62                 }
  63                 x -= dx;
  64         }
  65         gINT->Draw("p"); leg->AddEntry(gINT, Form("interpolation [%s]", method ? "INT" : "COG"), "pl");
  66
  67         gRE->Draw("pl"); leg->AddEntry(gRE, "1% relative errors", "pl");
  68         CHI2 /= npoints;
  69         printf("chi2(%d) = %f\n", npoints, CHI2);
  70
  71         leg->Draw();
  72 }
  73
  74
  75