rename Interpolator to PDF, add new class TKDInterpolator and
[u/mrichter/AliRoot.git] / STAT / TKDNodeInfo.cxx
1 #include "TKDNodeInfo.h"
2
3 #include "TVectorD.h"
4 #include "TMatrixD.h"
5 #include "TMath.h"
6
7 ClassImp(TKDNodeInfo)
8
9
10 //_________________________________________________________________
11 TKDNodeInfo::TKDNodeInfo(const Int_t dim):
12         TObject()
13         ,fNDim(3*dim)
14         ,fData(0x0)
15         ,fCov(0x0)
16         ,fPar(0x0)
17 {
18         fVal[0] = 0.; fVal[1] = 0.;
19         Build(dim);
20 }
21
22 //_________________________________________________________________
23 TKDNodeInfo::TKDNodeInfo(const TKDNodeInfo &ref):
24         TObject((TObject&) ref)
25         ,fNDim(fNDim)
26         ,fData(0x0)
27         ,fCov(0x0)
28         ,fPar(0x0)
29 {
30         Build(fNDim/3);
31
32         memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
33         fVal[0] = ref.fVal[0];
34         fVal[1] = ref.fVal[1];
35         if(ref.fCov) (*fCov) = (*ref.fCov);
36         if(ref.fPar) (*fPar) = (*ref.fPar);
37 }
38
39
40 //_________________________________________________________________
41 TKDNodeInfo::~TKDNodeInfo()
42 {
43         if(fData) delete [] fData;
44         if(fCov){
45                 delete fPar;
46                 delete fCov;
47         }
48 }
49
50 //_________________________________________________________________
51 TKDNodeInfo& TKDNodeInfo::operator=(const TKDNodeInfo & ref)
52 {
53 //      Info("operator==()", "...");
54         
55         Int_t ndim = fNDim/3;
56         if(fNDim != ref.fNDim){
57                 fNDim = ref.fNDim;
58                 Build(ndim);
59         }
60         memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
61         fVal[0] = ref.fVal[0];
62         fVal[1] = ref.fVal[1];
63         if(ref.fCov) (*fCov) = (*ref.fCov);
64         if(ref.fPar) (*fPar) = (*ref.fPar);
65         
66         return *this;
67 }
68
69 //_________________________________________________________________
70 void TKDNodeInfo::Build(const Int_t dim)
71 {
72 // Allocate/Reallocate space for this node.
73
74 //      Info("Build()", "...");
75
76         if(!dim) return;
77         
78         Int_t lambda = Int_t(1 + dim + .5*dim*(dim+1));
79         if(fData) delete [] fData;
80         fData = new Float_t[fNDim];
81         if(fCov){
82                 fCov->ResizeTo(lambda, lambda);
83                 fPar->ResizeTo(lambda);
84         }
85         return;
86 }
87
88 //_________________________________________________________________
89 void TKDNodeInfo::Print()
90 {
91         Int_t dim = fNDim/3;
92         printf("x[");
93         for(int idim=0; idim<dim; idim++) printf("%f ", fData[idim]);
94         printf("] f = [%f +- %f]\n", fVal[0], fVal[1]);
95
96         Float_t *bounds = &fData[dim];
97         printf("range[");
98         for(int idim=0; idim<dim; idim++) printf("(%f %f) ", fData[2*idim], fData[2*idim+1]);
99         printf("]\n");
100         
101         printf("Fit parameters : ");
102         if(!fCov){
103                 printf("Not defined.\n");
104                 return;
105         }
106         
107         Int_t lambda = Int_t(1 + dim + .5*dim*(dim+1));
108         for(int ip=0; ip<3; ip++) printf("p%d[%f] ", ip, (*fPar)(ip));
109         printf("\n");
110 }
111
112 //_________________________________________________________________
113 void TKDNodeInfo::Store(const TVectorD &par, const TMatrixD &cov)
114 {
115         if(!fCov){
116                 fCov = new TMatrixD(cov.GetNrows(), cov.GetNrows());
117                 fPar = new TVectorD(par.GetNrows());
118         }
119         (*fPar) = par;
120         (*fCov) = cov;
121 }
122
123 //_________________________________________________________________
124 Double_t TKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t &error)
125 {
126 // Recalculate the PDF for one node from the results of interpolation (parameters and covariance matrix)
127
128         Int_t ndim = fNDim/3;
129         if(ndim>10) return 0.; // support only up to 10 dimensions
130
131         Int_t lambda = 1 + ndim + (ndim*(ndim+1)>>1);
132         Double_t fdfdp[66];
133         Int_t ipar = 0;
134         fdfdp[ipar++] = 1.;
135         for(int idim=0; idim<ndim; idim++){
136                 fdfdp[ipar++] = point[idim];
137                 for(int jdim=idim; jdim<ndim; jdim++) fdfdp[ipar++] = point[idim]*point[jdim];
138         }
139
140         // calculate estimation
141         result =0.; error = 0.;
142         for(int i=0; i<lambda; i++){
143                 result += fdfdp[i]*(*fPar)(i);
144                 for(int j=0; j<lambda; j++) error += fdfdp[i]*fdfdp[j]*(*fCov)(i,j);
145         }       
146         error = TMath::Sqrt(error);
147         
148         //printf("TKDNodeInfo::CookPDF() : %6.3f +- %6.3f\n", result, error);
149
150         return 0.;
151 }
152