STAT/Macros/pdfIO.C

   1 TStopwatch timer;
   2 const Int_t nt = 100;
   3 //___________________________________________________________
   4 void pdfIO(const Int_t method = 0)
   5 {
   6 // Test interpolator IO response for several data dimensions.
   7 // method = 0 : use COG interpolator
   8 // method = 1 : use INT interpolator
   9
  10         Float_t tw, tr;
  11         TGraph *gw = new TGraph(5);
  12         gw->SetMarkerStyle(24);gw->SetMarkerColor(2);
  13         TGraph *gr = new TGraph(5);
  14         gr->SetMarkerStyle(25);gr->SetMarkerColor(4);
  15         for(int idim = 1; idim<6; idim++){
  16                 tw = interpolWrite(method, idim);
  17                 gw->SetPoint(idim-1, idim, tw);
  18                 tr = interpolRead(idim);
  19                 gr->SetPoint(idim-1, idim, tr);
  20         }
  21         gw->Draw("apl");
  22         gr->Draw("pl");
  23 }
  24
  25 //___________________________________________________________
  26 Float_t interpolWrite(const Int_t method, const Int_t ndim)
  27 {
  28         if(!gSystem->FindFile(".", "5D_Gauss.root")) build(5, 1000000);
  29         TFile::Open("5D_Gauss.root");
  30         printf("\nWriting %dD interpolator ...\n", ndim);
  31         TTree *db = (TTree*)gFile->Get("db");
  32         TString var = "x0"; for(int idim=1; idim<ndim; idim++) var += Form(":x%d", idim);
  33         printf("\tBuilding interpolator ...\n");
  34         TKDPDF pdf(db, var.Data(), "", 400, 100000);
  35         pdf.SetInterpolationMethod(method);
  36         pdf.SetStore();
  37
  38
  39         printf("\tSetting up the interpolator ...\n");
  40         Float_t *c, val, v_err;
  41         Double_t *p = new Double_t[ndim], res, r_err;
  42         timer.Start(kTRUE);
  43         for(int inode=0; inode<pdf.GetNTNodes(); inode++){
  44                 pdf.GetCOGPoint(inode, c, val, v_err);
  45                 for(int idim=0; idim<ndim; idim++) p[idim] = (Double_t)c[idim];
  46 //              printf("Evaluate for node [%d] ", inode);
  47 //              for(int idim=0; idim<ndim; idim++) printf("%5.3f ", p[idim]); printf("\n");
  48                 for(int it=0; it<nt; it++) pdf.Eval(p, res, r_err, kTRUE);
  49 //              printf("R = %6.4f +- %6.4f\n", res, r_err);
  50         }
  51         timer.Stop();
  52
  53         Float_t time = 1.E3*timer.CpuTime()/float(nt);
  54         printf("\tFinish interpolation in %6.2f [ms]\n", time);
  55
  56         printf("\tSaving interpolator ...\n");
  57         TFile *fi = TFile::Open(Form("%dD_Interpolator.root", ndim), "RECREATE");
  58         pdf.Write(Form("%dDgauss", ndim));
  59         fi->Close();
  60         delete fi;
  61
  62         delete [] p;
  63
  64         return time;
  65 }
  66
  67 //___________________________________________________________
  68 Float_t interpolRead(const Int_t ndim)
  69 {
  70         printf("Reading %dD interpolator ...\n", ndim);
  71         TFile::Open(Form("%dD_Interpolator.root", ndim));
  72         TKDPDF *pdf = (TKDPDF*)gFile->Get(Form("%dDgauss", ndim));
  73         gFile->Close();
  74
  75         printf("\tDoing interpolation ...\n");
  76         Float_t *c, v, ve;
  77         Double_t *p = new Double_t[ndim], r, re;
  78         timer.Start(kTRUE);
  79         for(int ip=0; ip<pdf->GetNTNodes(); ip++){
  80                 pdf->GetCOGPoint(ip, c, v, ve);
  81                 for(int idim=0; idim<ndim; idim++) p[idim] = (Double_t)c[idim];
  82                 for(int it=0; it<nt; it++) pdf->Eval(p, r, re);
  83         }
  84         timer.Stop();
  85         Float_t time = 1.E3*timer.CpuTime()/float(nt);
  86         printf("\tFinish interpolation in %6.2f [ms]\n", time);
  87         delete [] p;
  88         return time;
  89 }
  90
  91 //___________________________________________________________
  92 void build(const int ndim, const int npoints)
  93 {
  94         printf("Building DB ...\n");
  95         Float_t data[ndim];
  96         TFile::Open(Form("%dD_Gauss.root", ndim), "RECREATE");
  97         TTree *t = new TTree("db", Form("%dD data base for kD statistics", ndim));
  98         for(int idim=0; idim<ndim; idim++) t->Branch(Form("x%d", idim), &data[idim], Form("x%d/F", idim));
  99
 100         for (Int_t ip=0; ip<npoints; ip++){
 101                 for (Int_t id=0; id<ndim; id++) data[id]= gRandom->Gaus();
 102                 t->Fill();
 103         }
 104
 105         t->Write();
 106         gFile->Close();
 107         delete gFile;
 108 }
 109
 110