CMake: Removing all lhapdf link dependencies
[u/mrichter/AliRoot.git] / STAT / TKDNodeInfo.cxx
1 ////////////////////////////////////////////////////////
2 //
3 // Bucket representation for TKDInterpolator(Base)
4 //
5 // The class store data and provides the interface to draw and print.
6 // The bucket - generalized histogram bin in N dimensions is represented by unprocessed data like
7 //   - experimental PDF value and statistical error 
8 //   - COG position (n-tuplu)
9 //   - boundaries
10 // and interpolated data like
11 //   - parameters of the local parabolic fit
12 //   - their covariance matrix
13 //  
14 // For drawing 2D projections the helper class TKDNodeInfo::TKDNodeDraw is used.
15 //
16 // Author Alexandru Bercuci <A.Bercuci@gsi.de>
17 //
18 ////////////////////////////////////////////////////////
19
20 #include "TKDNodeInfo.h"
21
22 #include "TVectorD.h"
23 #include "TMatrixD.h"
24 #include "TRandom.h"
25 #include "TMath.h"
26
27 ClassImp(TKDNodeInfo)
28 ClassImp(TKDNodeInfo::TKDNodeDraw)
29
30
31 //_________________________________________________________________
32 TKDNodeInfo::TKDNodeInfo(Int_t dim):
33   TObject()
34   ,fNDim(3*dim)
35   ,fData(NULL)
36   ,fNpar(0)
37   ,fNcov(0)
38   ,fPar(NULL)
39   ,fCov(NULL)
40 {
41   // Default constructor
42   fVal[0] = 0.; fVal[1] = 0.;
43   Build(dim);
44 }
45
46 //_________________________________________________________________
47 TKDNodeInfo::TKDNodeInfo(const TKDNodeInfo &ref):
48   TObject((TObject&) ref)
49   ,fNDim(ref.fNDim)
50   ,fData(NULL)
51   ,fNpar(0)
52   ,fNcov(0)
53   ,fPar(NULL)
54   ,fCov(NULL)
55 {
56   // Copy constructor
57   Build(fNDim/3);
58
59   fData = new Float_t[fNDim];
60   memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
61   fVal[0] = ref.fVal[0];
62   fVal[1] = ref.fVal[1];
63   if(ref.fNpar&&ref.fPar){ 
64     fNpar = ref.fNpar;
65     fPar=new Double_t[fNpar];
66     memcpy(fPar, ref.fPar, fNpar*sizeof(Double_t));
67   }
68   if(ref.fNcov && ref.fCov){ 
69     fNcov = ref.fNcov;
70     fCov=new Double_t[fNcov];
71     memcpy(fCov, ref.fCov, fNcov*sizeof(Double_t));
72   }
73 }
74
75
76 //_________________________________________________________________
77 TKDNodeInfo::~TKDNodeInfo()
78 {
79   // Destructor
80   if(fData) delete [] fData;
81   if(fPar) delete [] fPar;
82   if(fCov) delete [] fCov;
83 }
84
85 //_________________________________________________________________
86 TKDNodeInfo& TKDNodeInfo::operator=(const TKDNodeInfo & ref)
87 {
88 //      Info("operator==()", "...");
89   
90   if(this == &ref) return *this;
91   Int_t ndim = fNDim/3;
92   if(fNDim != ref.fNDim){
93     fNDim = ref.fNDim;
94     Build(ndim);
95   }
96   memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
97   fVal[0] = ref.fVal[0];
98   fVal[1] = ref.fVal[1];
99   if(ref.fNpar&&ref.fPar){ 
100     fNpar = ref.fNpar;
101     fPar=new Double_t[fNpar];
102     memcpy(fPar, ref.fPar, fNpar*sizeof(Double_t));
103   }
104   if(ref.fNcov && ref.fCov){ 
105     fNcov = ref.fNcov;
106     fCov=new Double_t[fNcov];
107     memcpy(fCov, ref.fCov, fNcov*sizeof(Double_t));
108   }
109   return *this;
110 }
111
112 //_________________________________________________________________
113 void TKDNodeInfo::Build(Int_t dim)
114 {
115 // Allocate/Reallocate space for this node.
116
117 //      Info("Build()", "...");
118
119   if(!dim) return;
120   fNDim = 3*dim;
121   if(fData) delete [] fData;
122   fData = new Float_t[fNDim];
123   return;
124 }
125
126 //_________________________________________________________________
127 void TKDNodeInfo::Bootstrap()
128 {
129   if(!fNpar || !fPar) return;
130
131   Int_t ndim = Int_t(.5*(TMath::Sqrt(1.+8.*fNpar)-1.))-1;
132   fNDim = 3*ndim;
133 }
134
135 //_________________________________________________________________
136 void TKDNodeInfo::SetNode(Int_t ndim, Float_t *data, Float_t *pdf)
137 {
138   Build(ndim);
139   memcpy(fData, data, fNDim*sizeof(Float_t));
140   fVal[0]=pdf[0]; fVal[1]=pdf[1];
141 }
142
143
144 //_________________________________________________________________
145 void TKDNodeInfo::Print(const Option_t *opt) const
146 {
147   // Print the content of the node
148   Int_t dim = Int_t(fNDim/3.);
149   printf("x[");
150   for(int idim=0; idim<dim; idim++) printf("%f ", fData?fData[idim]:0.);
151   printf("] f = [%f +- %f]\n", fVal[0], fVal[1]);
152
153   if(fData){
154     Float_t *bounds = &fData[dim];
155     printf("range[");
156     for(int idim=0; idim<dim; idim++) printf("(%f %f) ", bounds[2*idim], bounds[2*idim+1]);
157     printf("]\n");
158   }
159   if(strcmp(opt, "a")!=0) return;
160
161   if(fNpar){ 
162     printf("Fit parameters : \n");
163     for(int ip=0; ip<fNpar; ip++) printf("p%d[%f] ", ip, fPar[ip]);
164     printf("\n");
165   }
166   if(!fNcov) return;
167   for(int ip(0), n(0); ip<fNpar; ip++){
168     for(int jp(ip); jp<fNpar; jp++) printf("c(%d %d)[%f] ", ip, jp, fCov[n++]);
169     printf("\n");
170   }
171 }
172
173 //_________________________________________________________________
174 void TKDNodeInfo::Store(TVectorD const *par, TMatrixD const *cov)
175 {
176 // Store the parameters and the covariance in the node
177
178   if(!fPar){SetNpar(); fPar = new Double_t[fNpar];}
179   for(int ip=0; ip<fNpar; ip++) fPar[ip] = (*par)[ip];
180
181   if(!cov) return;
182   if(!fCov){SetNcov(); fCov = new Double_t[fNcov];}
183   for(int ip(0), np(0); ip<fNpar; ip++)
184     for(int jp=ip; jp<fNpar; jp++) fCov[np++] = (*cov)(ip,jp);
185 }
186
187 //_________________________________________________________________
188 Bool_t TKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t &error) const
189 {
190 // Recalculate the PDF for one node from the results of interpolation (parameters and covariance matrix)
191
192   Int_t ndim = Int_t(fNDim/3.);
193   if(ndim>10) return kFALSE; // support only up to 10 dimensions
194   //printf("ndim[%d] npar[%d] ncov[%d]\n", ndim, fNpar, fNcov);
195
196   Double_t fdfdp[66]; memset(fdfdp, 0, ndim*sizeof(Double_t));
197   Int_t ipar = 0;
198   fdfdp[ipar++] = 1.;
199   for(int idim=0; idim<ndim; idim++){
200     fdfdp[ipar++] = point[idim];
201     for(int jdim=idim; jdim<ndim; jdim++) fdfdp[ipar++] = point[idim]*point[jdim];
202   }
203
204   // calculate estimation
205   result =0.; error = 0.;
206   for(int i=0; i<fNpar; i++) result += fdfdp[i]*fPar[i];
207   if(!fNcov) return kTRUE;
208
209   for(int i(0), n(0); i<fNpar; i++){
210     error += fdfdp[i]*fdfdp[i]*fCov[n++];
211     for(int j(i+1); j<fNpar; j++) error += 2.*fdfdp[i]*fdfdp[j]*fCov[n++];
212   }     
213   error = TMath::Sqrt(error);
214   
215   //printf("TKDNodeInfo::CookPDF() : %6.3f +- %6.3f\n", result, error);
216
217   return kTRUE;
218 }
219
220
221
222 //_________________________________________________________________
223 TKDNodeInfo::TKDNodeDraw::TKDNodeDraw() 
224   :TBox()
225   ,fCOG()
226   ,fNode(NULL)
227 {
228   SetFillStyle(3002);
229   SetFillColor(50+Int_t(gRandom->Uniform()*50.));
230
231   fCOG.SetMarkerStyle(3);
232   fCOG.SetMarkerSize(.7);
233   fCOG.SetMarkerColor(2);
234 }
235
236
237 //_________________________________________________________________
238 void TKDNodeInfo::TKDNodeDraw::Draw(Option_t* option)
239 {
240   TBox::Draw(option);
241   fCOG.Draw("p");
242 }
243
244 //_________________________________________________________________
245 void TKDNodeInfo::TKDNodeDraw::SetNode(TKDNodeInfo *node, UChar_t size, UChar_t ax1, UChar_t ax2)
246 {
247   fNode=node;
248   const Float_t kBorder = 0.;//1.E-4;
249   Float_t *bounds = &(node->Data()[size]);
250   fX1=bounds[2*ax1]+kBorder;
251   fX2=bounds[2*ax1+1]-kBorder;
252   fY1=bounds[2*ax2]+kBorder;
253   fY2=bounds[2*ax2+1]-kBorder;
254   
255   Float_t x(node->Data()[ax1]), y(node->Data()[ax2]);
256   fCOG.SetX(x); fCOG.SetY(y);
257 }
258
259
260 //_________________________________________________________________
261 void TKDNodeInfo::TKDNodeDraw::Print(const Option_t* option) const
262 {
263   if(!fNode) return;
264   fNode->Print(option);
265 }