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