STAT/Macros/pdfTest.C

   1 #include <malloc.h>
   2
   3 #include "TSystem.h"
   4 #include "TFile.h"
   5 #include "TTree.h"
   6 #include "TProfile.h"
   7 #include "TF1.h"
   8 #include "TF2.h"
   9 #include "TF3.h"
  10 #include "TLegend.h"
  11 #include "TRandom.h"
  12 #include "TString.h"
  13 #include "TGraph.h"
  14 #include "TMarker.h"
  15 #include "TStopwatch.h"
  16 #include "../src/TKDPDF.h"
  17
  18 void pdfTest(const Int_t ndim = 1);
  19 Double_t interpolation(const int nstat, const int ndim, const int first);
  20 void build(const int ndim, const int npoints);
  21 Float_t Mem();
  22
  23 TF1 *f = 0x0;
  24
  25 //__________________________________________________________
  26 void pdfTest(const Int_t ndim)
  27 {
  28 // Test interpolator on a variable dimension (ndim<=5) uncorrelated
  29 // Gauss model. The macro displays the chi2 between interpolation and
  30 // model for various statistics of experimental data.
  31 //
  32 // Check the worker function interpolation() to fine tune the
  33 // algorithm.
  34 //
  35 // Attention:
  36 // The macro works in compiled mode !
  37
  38         // define model
  39         switch(ndim){
  40         case 1:
  41                 f=new TF1("f1", "gaus(0);x;f(x)", -5., 5.);
  42                 f->SetParameter(0, 1.);
  43                 f->SetParameter(1, 0.);
  44                 f->SetParameter(2, 1.);
  45                 f->SetParameter(0, 1./f->Integral(-5., 5.));
  46                 break;
  47         case 2:
  48                 f=new TF2("f2", "[0]*exp(-0.5*((x-[1])/[2])**2)*exp(-0.5*((y-[3])/[4])**2)", -5., 5., -5., 5.);
  49                 f->SetParameter(0, 1.);
  50                 f->SetParameter(1, 0.);
  51                 f->SetParameter(2, 1.);
  52                 f->SetParameter(3, 0.);
  53                 f->SetParameter(4, 1.);
  54                 f->SetParameter(0, 1./f->Integral(-5., 5., -5., 5.));
  55                 break;
  56         case 3:
  57         case 4:
  58         case 5:
  59                 f=new TF3("f3", "[0]*exp(-0.5*((x-[1])/[2])**2)*exp(-0.5*((y-[3])/[4])**2)*exp(-0.5*((z-[5])/[6])**2)", -5., 5., -5., 5., -5., 5.);
  60                 f->SetParameter(0, 1.);
  61                 f->SetParameter(1, 0.);
  62                 f->SetParameter(2, 1.);
  63                 f->SetParameter(3, 0.);
  64                 f->SetParameter(4, 1.);
  65                 f->SetParameter(5, 0.);
  66                 f->SetParameter(6, 1.);
  67                 f->SetParameter(0, 1./f->Integral(-5., 5., -5., 5., -5., 5.));
  68                 if(ndim>3) printf("Warning : chi2 results are unreliable !\n");
  69                 break;
  70         default:
  71                 printf("%dD not supported in this macro.\n", ndim);
  72                 return;
  73         }
  74
  75         Double_t chi2;
  76         const Int_t nchi2 = 18;
  77         const Int_t nfirst = 1;
  78         //TProfile *hp = new TProfile("hp", "", 100, 1.E4, 1.E6);
  79         //hp->SetMarkerStyle(24);
  80         TGraph *gChi2 = new TGraph(nchi2);
  81         gChi2->SetMarkerStyle(24);
  82         Int_t ip = 0, nstat, first;
  83         for(int ifirst=0; ifirst<nfirst; ifirst++){
  84                 first = 0;//Int_t(gRandom->Uniform(1.E4));
  85                 for(int i=2; i<10; i++){
  86                         nstat = 10000*i;
  87                         chi2 = interpolation(nstat, ndim, first);
  88                         gChi2->SetPoint(ip++, TMath::Log10(float(nstat)), chi2);
  89                         //hp->Fill(float(nstat), chi2);
  90                 }
  91                 for(int i=1; i<10; i++){
  92                         nstat = 100000*i;
  93                         chi2 = interpolation(nstat, ndim, first);
  94                         gChi2->SetPoint(ip++, TMath::Log10(float(nstat)), chi2);
  95                         //hp->Fill(float(nstat), chi2);
  96                 }
  97                 gChi2->Draw("apl");
  98         }
  99         //hp->Draw("pl");
 100 }
 101
 102
 103 //__________________________________________________________
 104 Double_t interpolation(const int nstat, const int ndim, const int first)
 105 {
 106 // Steer interpolation of a nD (n<=5) uncorrelated Gauss distribution.
 107 // The user is suppose to give the experimental statistics (nstat) the
 108 // dimension of the experimental point (ndim). The entry from which the
 109 // data base is being read is steered by "first".
 110
 111         const int bs = 400;
 112         const int npoints = 100;
 113         // switch on chi2 calculation between interpolation and model.
 114         // Default calculates distance.
 115         const Bool_t kCHI2 = kFALSE;
 116         const Bool_t kINT  = kTRUE; // switch on integral interpolator
 117
 118         // open data base
 119         TString fn("5D_Gauss.root");
 120         if(!gSystem->FindFile(".", fn)) build(5, 1000000);
 121         TFile::Open("5D_Gauss.root");
 122         TTree *t = (TTree*)gFile->Get("db");
 123
 124         // build interpolator
 125         TString var = "x0"; for(int idim=1; idim<ndim; idim++) var += Form(":x%d", idim);
 126         TKDPDF pdf(t, var.Data(), "", bs, nstat, first);
 127         pdf.SetInterpolationMethod(kINT);
 128         pdf.SetStore();
 129
 130         // define testing grid
 131         Double_t x[5];
 132         Double_t dx[5] = {4., 4., 4., 4., 4.};
 133         Double_t chi2, theor, result, err;
 134         for(int idim=0; idim<ndim; idim++){
 135                 dx[idim] = (ndim<<2)/float(npoints);
 136                 x[idim]  = -2.+dx[idim]/2.;
 137         }
 138
 139         // setup the integral interpolator and do memory test
 140         Float_t start, stop;
 141         start = Mem();
 142         Float_t *c, v, ve;
 143         for(int ip=0; ip<pdf.GetNTNodes(); ip++){
 144                 pdf.GetCOGPoint(ip, c, v, ve);
 145                 for(int idim=0; idim<ndim; idim++) x[idim] = (Double_t)c[idim];
 146                 pdf.Eval(x, result, err);
 147         }
 148         stop  = Mem();
 149         printf("\nInterpolating %d data points in %dD ...\n", nstat, ndim);
 150         printf("\tbucket size %d\n", bs);
 151         printf("\tmemory usage in Eval() %6.4fKB\n", stop-start);
 152         printf("\tstarting interpolation ...\n");
 153         //return stop-start;
 154
 155         // evaluate relative error
 156         for(int idim=0; idim<ndim; idim++) x[idim] = 0.;
 157         pdf.Eval(x, result, err);
 158         Double_t threshold = err/result;
 159
 160         // starting interpolation over the testing grid
 161         chi2 = 0.; Int_t nsample = 0;
 162         x[4]  = -2.+dx[4]/2.;
 163         do{
 164                 x[3]  = -2.+dx[3]/2.;
 165                 do{
 166                         x[2]  = -2.+dx[2]/2.;
 167                         do{
 168                                 x[1]  = -2.+dx[1]/2.;
 169                                 do{
 170                                         x[0]  = -2.+dx[0]/2.;
 171                                         do{
 172                                                 pdf.Eval(x, result, err);
 173                                                 if(err/TMath::Abs(result) < 1.1*threshold){
 174                                                         theor = f->Eval(x[0], x[1], x[2]);
 175                                                         Double_t tmp = result - theor; tmp *= tmp;
 176                                                         chi2 += tmp/(kCHI2 ? err*err : theor);
 177                                                         nsample++;
 178                                                 }
 179                                         } while((x[0] += dx[0]) < 2.);
 180                                 } while((x[1] += dx[1]) < 2.);
 181                         } while((x[2] += dx[2]) < 2.);
 182                 } while((x[3] += dx[3]) < 2.);
 183         } while((x[4] += dx[4]) < 2.);
 184         chi2 /= float(nsample);
 185         printf("\tinterpolation quality (chi2) %f\n", chi2);
 186
 187         return kCHI2 ? chi2 : TMath::Sqrt(chi2);
 188 }
 189
 190
 191 //___________________________________________________________
 192 void build(const int ndim, const int npoints)
 193 {
 194         printf("Building DB ...\n");
 195         Float_t data[ndim];
 196         TFile::Open(Form("%dD_Gauss.root", ndim), "RECREATE");
 197         TTree *t = new TTree("db", Form("%dD data base for kD statistics", ndim));
 198         for(int idim=0; idim<ndim; idim++) t->Branch(Form("x%d", idim), &data[idim], Form("x%d/F", idim));
 199
 200         for (Int_t ip=0; ip<npoints; ip++){
 201                 for (Int_t id=0; id<ndim; id++) data[id]= gRandom->Gaus();
 202                 t->Fill();
 203         }
 204
 205         t->Write();
 206         gFile->Close();
 207         delete gFile;
 208 }
 209
 210 //______________________________________________________________
 211 Float_t Mem()
 212 {
 213   // size in KB
 214   static struct mallinfo debug;
 215   debug = mallinfo();
 216         //printf("arena[%d] usmblks[%d] uordblks[%d] fordblks[%d]\n", debug.arena, debug.usmblks, debug.uordblks, debug.fordblks);
 217   return debug.uordblks/1024.;
 218 }
 219