Updated macros for PHOS alignment calculation
[u/mrichter/AliRoot.git] / STAT / TKDNodeInfo.cxx
CommitLineData
b273c7cb 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
a3408ed3 20#include "TKDNodeInfo.h"
21
22#include "TVectorD.h"
23#include "TMatrixD.h"
afb669c1 24#include "TRandom.h"
a3408ed3 25#include "TMath.h"
26
27ClassImp(TKDNodeInfo)
b6f6b31a 28ClassImp(TKDNodeInfo::TKDNodeDraw)
a3408ed3 29
30
31//_________________________________________________________________
1a439134 32TKDNodeInfo::TKDNodeInfo(Int_t dim):
b6f6b31a 33 TObject()
34 ,fNDim(3*dim)
afb669c1 35 ,fData(NULL)
36 ,fNpar(0)
37 ,fNcov(0)
38 ,fPar(NULL)
39 ,fCov(NULL)
a3408ed3 40{
1a439134 41 // Default constructor
b6f6b31a 42 fVal[0] = 0.; fVal[1] = 0.;
43 Build(dim);
a3408ed3 44}
45
46//_________________________________________________________________
47TKDNodeInfo::TKDNodeInfo(const TKDNodeInfo &ref):
b6f6b31a 48 TObject((TObject&) ref)
20932c8c 49 ,fNDim(ref.fNDim)
afb669c1 50 ,fData(NULL)
51 ,fNpar(0)
52 ,fNcov(0)
53 ,fPar(NULL)
54 ,fCov(NULL)
a3408ed3 55{
1a439134 56 // Copy constructor
b6f6b31a 57 Build(fNDim/3);
a3408ed3 58
20932c8c 59 fData = new Float_t[fNDim];
b6f6b31a 60 memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
61 fVal[0] = ref.fVal[0];
62 fVal[1] = ref.fVal[1];
afb669c1 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 }
a3408ed3 73}
74
75
76//_________________________________________________________________
77TKDNodeInfo::~TKDNodeInfo()
78{
1a439134 79 // Destructor
b6f6b31a 80 if(fData) delete [] fData;
afb669c1 81 if(fPar) delete [] fPar;
82 if(fCov) delete [] fCov;
a3408ed3 83}
84
85//_________________________________________________________________
86TKDNodeInfo& TKDNodeInfo::operator=(const TKDNodeInfo & ref)
87{
88// Info("operator==()", "...");
b6f6b31a 89
067a4e53 90 if(this == &ref) return *this;
b6f6b31a 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];
afb669c1 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 }
b6f6b31a 109 return *this;
a3408ed3 110}
111
112//_________________________________________________________________
1a439134 113void TKDNodeInfo::Build(Int_t dim)
a3408ed3 114{
115// Allocate/Reallocate space for this node.
116
117// Info("Build()", "...");
118
b6f6b31a 119 if(!dim) return;
afb669c1 120 fNDim = 3*dim;
b6f6b31a 121 if(fData) delete [] fData;
122 fData = new Float_t[fNDim];
b6f6b31a 123 return;
a3408ed3 124}
125
126//_________________________________________________________________
afb669c1 127void 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//_________________________________________________________________
b273c7cb 136void 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//_________________________________________________________________
afb669c1 145void TKDNodeInfo::Print(const Option_t *opt) const
a3408ed3 146{
1a439134 147 // Print the content of the node
afb669c1 148 Int_t dim = Int_t(fNDim/3.);
b6f6b31a 149 printf("x[");
afb669c1 150 for(int idim=0; idim<dim; idim++) printf("%f ", fData?fData[idim]:0.);
b6f6b31a 151 printf("] f = [%f +- %f]\n", fVal[0], fVal[1]);
152
afb669c1 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");
b6f6b31a 170 }
a3408ed3 171}
172
173//_________________________________________________________________
afb669c1 174void TKDNodeInfo::Store(TVectorD const *par, TMatrixD const *cov)
a3408ed3 175{
afb669c1 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);
a3408ed3 185}
186
187//_________________________________________________________________
afb669c1 188Bool_t TKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t &error) const
a3408ed3 189{
190// Recalculate the PDF for one node from the results of interpolation (parameters and covariance matrix)
191
afb669c1 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);
b6f6b31a 195
afb669c1 196 Double_t fdfdp[66]; memset(fdfdp, 0, ndim*sizeof(Double_t));
b6f6b31a 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.;
afb669c1 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++];
b6f6b31a 212 }
213 error = TMath::Sqrt(error);
214
215 //printf("TKDNodeInfo::CookPDF() : %6.3f +- %6.3f\n", result, error);
216
afb669c1 217 return kTRUE;
a3408ed3 218}
219
b6f6b31a 220
221
222//_________________________________________________________________
b273c7cb 223TKDNodeInfo::TKDNodeDraw::TKDNodeDraw()
224 :TBox()
225 ,fCOG()
afb669c1 226 ,fNode(NULL)
b6f6b31a 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//_________________________________________________________________
238void TKDNodeInfo::TKDNodeDraw::Draw(Option_t* option)
239{
240 TBox::Draw(option);
241 fCOG.Draw("p");
242}
243
244//_________________________________________________________________
245void TKDNodeInfo::TKDNodeDraw::SetNode(TKDNodeInfo *node, UChar_t size, UChar_t ax1, UChar_t ax2)
246{
b273c7cb 247 fNode=node;
b6f6b31a 248 const Float_t kBorder = 0.;//1.E-4;
249 Float_t *bounds = &(node->Data()[size]);
b6f6b31a 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]);
b6f6b31a 256 fCOG.SetX(x); fCOG.SetY(y);
257}
258
259
b273c7cb 260//_________________________________________________________________
261void TKDNodeInfo::TKDNodeDraw::Print(const Option_t* option) const
262{
263 if(!fNode) return;
264 fNode->Print(option);
265}