From Philippe & Laurent: new variant of MUON visualization.
[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
b6f6b31a 22#include "TRandom.h"
a3408ed3 23#include "TVectorD.h"
24#include "TMatrixD.h"
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)
35 ,fData(0x0)
36 ,fCov(0x0)
37 ,fPar(0x0)
a3408ed3 38{
1a439134 39 // Default constructor
b6f6b31a 40 fVal[0] = 0.; fVal[1] = 0.;
41 Build(dim);
a3408ed3 42}
43
44//_________________________________________________________________
45TKDNodeInfo::TKDNodeInfo(const TKDNodeInfo &ref):
b6f6b31a 46 TObject((TObject&) ref)
47 ,fNDim(fNDim)
48 ,fData(0x0)
49 ,fCov(0x0)
50 ,fPar(0x0)
a3408ed3 51{
1a439134 52 // Copy constructor
b6f6b31a 53 Build(fNDim/3);
a3408ed3 54
b6f6b31a 55 memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
56 fVal[0] = ref.fVal[0];
57 fVal[1] = ref.fVal[1];
58 if(ref.fCov) (*fCov) = (*ref.fCov);
59 if(ref.fPar) (*fPar) = (*ref.fPar);
a3408ed3 60}
61
62
63//_________________________________________________________________
64TKDNodeInfo::~TKDNodeInfo()
65{
1a439134 66 // Destructor
b6f6b31a 67 if(fData) delete [] fData;
68 if(fCov){
69 delete fPar;
70 delete fCov;
71 }
a3408ed3 72}
73
74//_________________________________________________________________
75TKDNodeInfo& TKDNodeInfo::operator=(const TKDNodeInfo & ref)
76{
77// Info("operator==()", "...");
b6f6b31a 78
79 Int_t ndim = fNDim/3;
80 if(fNDim != ref.fNDim){
81 fNDim = ref.fNDim;
82 Build(ndim);
83 }
84 memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
85 fVal[0] = ref.fVal[0];
86 fVal[1] = ref.fVal[1];
87 if(ref.fCov) (*fCov) = (*ref.fCov);
88 if(ref.fPar) (*fPar) = (*ref.fPar);
89
90 return *this;
a3408ed3 91}
92
93//_________________________________________________________________
1a439134 94void TKDNodeInfo::Build(Int_t dim)
a3408ed3 95{
96// Allocate/Reallocate space for this node.
97
98// Info("Build()", "...");
99
b6f6b31a 100 if(!dim) return;
b273c7cb 101 fNDim = 3*dim;
b6f6b31a 102 Int_t lambda = Int_t(1 + dim + .5*dim*(dim+1));
103 if(fData) delete [] fData;
104 fData = new Float_t[fNDim];
105 if(fCov){
106 fCov->ResizeTo(lambda, lambda);
107 fPar->ResizeTo(lambda);
108 }
109 return;
a3408ed3 110}
111
112//_________________________________________________________________
b273c7cb 113void TKDNodeInfo::SetNode(Int_t ndim, Float_t *data, Float_t *pdf)
114{
115 Build(ndim);
116 memcpy(fData, data, fNDim*sizeof(Float_t));
117 fVal[0]=pdf[0]; fVal[1]=pdf[1];
118}
119
120
121//_________________________________________________________________
1a439134 122void TKDNodeInfo::Print(const Option_t *) const
a3408ed3 123{
1a439134 124 // Print the content of the node
b6f6b31a 125 Int_t dim = fNDim/3;
126 printf("x[");
127 for(int idim=0; idim<dim; idim++) printf("%f ", fData[idim]);
128 printf("] f = [%f +- %f]\n", fVal[0], fVal[1]);
129
b273c7cb 130 Float_t *bounds = &fData[dim];
b6f6b31a 131 printf("range[");
b273c7cb 132 for(int idim=0; idim<dim; idim++) printf("(%f %f) ", bounds[2*idim], bounds[2*idim+1]);
b6f6b31a 133 printf("]\n");
134
135 printf("Fit parameters : ");
136 if(!fCov){
137 printf("Not defined.\n");
138 return;
139 }
140
141 // Int_t lambda = Int_t(1 + dim + .5*dim*(dim+1));
142 for(int ip=0; ip<3; ip++) printf("p%d[%f] ", ip, (*fPar)(ip));
143 printf("\n");
a3408ed3 144}
145
146//_________________________________________________________________
147void TKDNodeInfo::Store(const TVectorD &par, const TMatrixD &cov)
148{
1a439134 149 // Store the parameters and the covariance in the node
b6f6b31a 150 if(!fCov){
151 fCov = new TMatrixD(cov.GetNrows(), cov.GetNrows());
152 fPar = new TVectorD(par.GetNrows());
153 }
154 (*fPar) = par;
155 (*fCov) = cov;
a3408ed3 156}
157
158//_________________________________________________________________
159Double_t TKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t &error)
160{
161// Recalculate the PDF for one node from the results of interpolation (parameters and covariance matrix)
162
b6f6b31a 163 Int_t ndim = fNDim/3;
164 if(ndim>10) return 0.; // support only up to 10 dimensions
165
166 Int_t lambda = 1 + ndim + (ndim*(ndim+1)>>1);
167 Double_t fdfdp[66];
168 Int_t ipar = 0;
169 fdfdp[ipar++] = 1.;
170 for(int idim=0; idim<ndim; idim++){
171 fdfdp[ipar++] = point[idim];
172 for(int jdim=idim; jdim<ndim; jdim++) fdfdp[ipar++] = point[idim]*point[jdim];
173 }
174
175 // calculate estimation
176 result =0.; error = 0.;
177 for(int i=0; i<lambda; i++){
178 result += fdfdp[i]*(*fPar)(i);
179 for(int j=0; j<lambda; j++) error += fdfdp[i]*fdfdp[j]*(*fCov)(i,j);
180 }
181 error = TMath::Sqrt(error);
182
183 //printf("TKDNodeInfo::CookPDF() : %6.3f +- %6.3f\n", result, error);
184
185 return 0.;
a3408ed3 186}
187
b6f6b31a 188
189
190//_________________________________________________________________
b273c7cb 191TKDNodeInfo::TKDNodeDraw::TKDNodeDraw()
192 :TBox()
193 ,fCOG()
194 ,fNode(0x0)
b6f6b31a 195{
196 SetFillStyle(3002);
197 SetFillColor(50+Int_t(gRandom->Uniform()*50.));
198
199 fCOG.SetMarkerStyle(3);
200 fCOG.SetMarkerSize(.7);
201 fCOG.SetMarkerColor(2);
202}
203
204
205//_________________________________________________________________
206void TKDNodeInfo::TKDNodeDraw::Draw(Option_t* option)
207{
208 TBox::Draw(option);
209 fCOG.Draw("p");
210}
211
212//_________________________________________________________________
213void TKDNodeInfo::TKDNodeDraw::SetNode(TKDNodeInfo *node, UChar_t size, UChar_t ax1, UChar_t ax2)
214{
b273c7cb 215 fNode=node;
b6f6b31a 216 const Float_t kBorder = 0.;//1.E-4;
217 Float_t *bounds = &(node->Data()[size]);
b6f6b31a 218 fX1=bounds[2*ax1]+kBorder;
219 fX2=bounds[2*ax1+1]-kBorder;
220 fY1=bounds[2*ax2]+kBorder;
221 fY2=bounds[2*ax2+1]-kBorder;
222
223 Float_t x(node->Data()[ax1]), y(node->Data()[ax2]);
b6f6b31a 224 fCOG.SetX(x); fCOG.SetY(y);
225}
226
227
b273c7cb 228//_________________________________________________________________
229void TKDNodeInfo::TKDNodeDraw::Print(const Option_t* option) const
230{
231 if(!fNode) return;
232 fNode->Print(option);
233}