Coding conventions
[u/mrichter/AliRoot.git] / STAT / TKDNodeInfo.cxx
CommitLineData
a3408ed3 1#include "TKDNodeInfo.h"
2
3#include "TVectorD.h"
4#include "TMatrixD.h"
5#include "TMath.h"
6
7ClassImp(TKDNodeInfo)
8
9
10//_________________________________________________________________
1a439134 11TKDNodeInfo::TKDNodeInfo(Int_t dim):
a3408ed3 12 TObject()
13 ,fNDim(3*dim)
14 ,fData(0x0)
15 ,fCov(0x0)
16 ,fPar(0x0)
17{
1a439134 18 // Default constructor
a3408ed3 19 fVal[0] = 0.; fVal[1] = 0.;
20 Build(dim);
21}
22
23//_________________________________________________________________
24TKDNodeInfo::TKDNodeInfo(const TKDNodeInfo &ref):
25 TObject((TObject&) ref)
26 ,fNDim(fNDim)
27 ,fData(0x0)
28 ,fCov(0x0)
29 ,fPar(0x0)
30{
1a439134 31 // Copy constructor
a3408ed3 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//_________________________________________________________________
43TKDNodeInfo::~TKDNodeInfo()
44{
1a439134 45 // Destructor
a3408ed3 46 if(fData) delete [] fData;
47 if(fCov){
48 delete fPar;
49 delete fCov;
50 }
51}
52
53//_________________________________________________________________
54TKDNodeInfo& 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//_________________________________________________________________
1a439134 73void TKDNodeInfo::Build(Int_t dim)
a3408ed3 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//_________________________________________________________________
1a439134 92void TKDNodeInfo::Print(const Option_t *) const
a3408ed3 93{
1a439134 94 // Print the content of the node
a3408ed3 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
1a439134 100 // Float_t *bounds = &fData[dim];
a3408ed3 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
1a439134 111 // Int_t lambda = Int_t(1 + dim + .5*dim*(dim+1));
a3408ed3 112 for(int ip=0; ip<3; ip++) printf("p%d[%f] ", ip, (*fPar)(ip));
113 printf("\n");
114}
115
116//_________________________________________________________________
117void TKDNodeInfo::Store(const TVectorD &par, const TMatrixD &cov)
118{
1a439134 119 // Store the parameters and the covariance in the node
a3408ed3 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//_________________________________________________________________
129Double_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