new macros for testing the STAT package
[u/mrichter/AliRoot.git] / 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