Line data Source code
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 128 : ClassImp(TKDNodeInfo)
28 128 : ClassImp(TKDNodeInfo::TKDNodeDraw)
29 :
30 :
31 : //_________________________________________________________________
32 : TKDNodeInfo::TKDNodeInfo(Int_t dim):
33 22002 : TObject()
34 22002 : ,fNDim(3*dim)
35 22002 : ,fData(NULL)
36 22002 : ,fNpar(0)
37 22002 : ,fNcov(0)
38 22002 : ,fPar(NULL)
39 22002 : ,fCov(NULL)
40 110010 : {
41 : // Default constructor
42 22002 : fVal[0] = 0.; fVal[1] = 0.;
43 22002 : Build(dim);
44 44004 : }
45 :
46 : //_________________________________________________________________
47 : TKDNodeInfo::TKDNodeInfo(const TKDNodeInfo &ref):
48 0 : TObject((TObject&) ref)
49 0 : ,fNDim(ref.fNDim)
50 0 : ,fData(NULL)
51 0 : ,fNpar(0)
52 0 : ,fNcov(0)
53 0 : ,fPar(NULL)
54 0 : ,fCov(NULL)
55 0 : {
56 : // Copy constructor
57 0 : Build(fNDim/3);
58 :
59 0 : fData = new Float_t[fNDim];
60 0 : memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
61 0 : fVal[0] = ref.fVal[0];
62 0 : fVal[1] = ref.fVal[1];
63 0 : if(ref.fNpar&&ref.fPar){
64 0 : fNpar = ref.fNpar;
65 0 : fPar=new Double_t[fNpar];
66 0 : memcpy(fPar, ref.fPar, fNpar*sizeof(Double_t));
67 0 : }
68 0 : if(ref.fNcov && ref.fCov){
69 0 : fNcov = ref.fNcov;
70 0 : fCov=new Double_t[fNcov];
71 0 : memcpy(fCov, ref.fCov, fNcov*sizeof(Double_t));
72 0 : }
73 0 : }
74 :
75 :
76 : //_________________________________________________________________
77 : TKDNodeInfo::~TKDNodeInfo()
78 12 : {
79 : // Destructor
80 2 : if(fData) delete [] fData;
81 2 : if(fPar) delete [] fPar;
82 2 : if(fCov) delete [] fCov;
83 6 : }
84 :
85 : //_________________________________________________________________
86 : TKDNodeInfo& TKDNodeInfo::operator=(const TKDNodeInfo & ref)
87 : {
88 : // Info("operator==()", "...");
89 :
90 0 : if(this == &ref) return *this;
91 0 : Int_t ndim = fNDim/3;
92 0 : if(fNDim != ref.fNDim){
93 0 : fNDim = ref.fNDim;
94 0 : Build(ndim);
95 0 : }
96 0 : memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
97 0 : fVal[0] = ref.fVal[0];
98 0 : fVal[1] = ref.fVal[1];
99 0 : if(ref.fNpar&&ref.fPar){
100 0 : fNpar = ref.fNpar;
101 0 : fPar=new Double_t[fNpar];
102 0 : memcpy(fPar, ref.fPar, fNpar*sizeof(Double_t));
103 0 : }
104 0 : if(ref.fNcov && ref.fCov){
105 0 : fNcov = ref.fNcov;
106 0 : fCov=new Double_t[fNcov];
107 0 : memcpy(fCov, ref.fCov, fNcov*sizeof(Double_t));
108 0 : }
109 : return *this;
110 0 : }
111 :
112 : //_________________________________________________________________
113 : void TKDNodeInfo::Build(Int_t dim)
114 : {
115 : // Allocate/Reallocate space for this node.
116 :
117 : // Info("Build()", "...");
118 :
119 44004 : if(!dim) return;
120 0 : fNDim = 3*dim;
121 0 : if(fData) delete [] fData;
122 0 : fData = new Float_t[fNDim];
123 0 : return;
124 22002 : }
125 :
126 : //_________________________________________________________________
127 : void TKDNodeInfo::Bootstrap()
128 : {
129 0 : if(!fNpar || !fPar) return;
130 :
131 0 : Int_t ndim = Int_t(.5*(TMath::Sqrt(1.+8.*fNpar)-1.))-1;
132 0 : fNDim = 3*ndim;
133 0 : }
134 :
135 : //_________________________________________________________________
136 : void TKDNodeInfo::SetNode(Int_t ndim, Float_t *data, Float_t *pdf)
137 : {
138 0 : Build(ndim);
139 0 : memcpy(fData, data, fNDim*sizeof(Float_t));
140 0 : fVal[0]=pdf[0]; fVal[1]=pdf[1];
141 0 : }
142 :
143 :
144 : //_________________________________________________________________
145 : void TKDNodeInfo::Print(const Option_t *opt) const
146 : {
147 : // Print the content of the node
148 0 : Int_t dim = Int_t(fNDim/3.);
149 0 : printf("x[");
150 0 : for(int idim=0; idim<dim; idim++) printf("%f ", fData?fData[idim]:0.);
151 0 : printf("] f = [%f +- %f]\n", fVal[0], fVal[1]);
152 :
153 0 : if(fData){
154 0 : Float_t *bounds = &fData[dim];
155 0 : printf("range[");
156 0 : for(int idim=0; idim<dim; idim++) printf("(%f %f) ", bounds[2*idim], bounds[2*idim+1]);
157 0 : printf("]\n");
158 0 : }
159 0 : if(strcmp(opt, "a")!=0) return;
160 :
161 0 : if(fNpar){
162 0 : printf("Fit parameters : \n");
163 0 : for(int ip=0; ip<fNpar; ip++) printf("p%d[%f] ", ip, fPar[ip]);
164 0 : printf("\n");
165 0 : }
166 0 : if(!fNcov) return;
167 0 : for(int ip(0), n(0); ip<fNpar; ip++){
168 0 : for(int jp(ip); jp<fNpar; jp++) printf("c(%d %d)[%f] ", ip, jp, fCov[n++]);
169 0 : printf("\n");
170 : }
171 0 : }
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 0 : if(!fPar){SetNpar(); fPar = new Double_t[fNpar];}
179 0 : for(int ip=0; ip<fNpar; ip++) fPar[ip] = (*par)[ip];
180 :
181 0 : if(!cov) return;
182 0 : if(!fCov){SetNcov(); fCov = new Double_t[fNcov];}
183 0 : for(int ip(0), np(0); ip<fNpar; ip++)
184 0 : for(int jp=ip; jp<fNpar; jp++) fCov[np++] = (*cov)(ip,jp);
185 0 : }
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 0 : Int_t ndim = Int_t(fNDim/3.);
193 0 : 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 0 : Double_t fdfdp[66]; memset(fdfdp, 0, ndim*sizeof(Double_t));
197 : Int_t ipar = 0;
198 0 : fdfdp[ipar++] = 1.;
199 0 : for(int idim=0; idim<ndim; idim++){
200 0 : fdfdp[ipar++] = point[idim];
201 0 : for(int jdim=idim; jdim<ndim; jdim++) fdfdp[ipar++] = point[idim]*point[jdim];
202 : }
203 :
204 : // calculate estimation
205 0 : result =0.; error = 0.;
206 0 : for(int i=0; i<fNpar; i++) result += fdfdp[i]*fPar[i];
207 0 : if(!fNcov) return kTRUE;
208 :
209 0 : for(int i(0), n(0); i<fNpar; i++){
210 0 : error += fdfdp[i]*fdfdp[i]*fCov[n++];
211 0 : for(int j(i+1); j<fNpar; j++) error += 2.*fdfdp[i]*fdfdp[j]*fCov[n++];
212 : }
213 0 : error = TMath::Sqrt(error);
214 :
215 : //printf("TKDNodeInfo::CookPDF() : %6.3f +- %6.3f\n", result, error);
216 :
217 0 : return kTRUE;
218 0 : }
219 :
220 :
221 :
222 : //_________________________________________________________________
223 : TKDNodeInfo::TKDNodeDraw::TKDNodeDraw()
224 0 : :TBox()
225 0 : ,fCOG()
226 0 : ,fNode(NULL)
227 0 : {
228 0 : SetFillStyle(3002);
229 0 : SetFillColor(50+Int_t(gRandom->Uniform()*50.));
230 :
231 0 : fCOG.SetMarkerStyle(3);
232 0 : fCOG.SetMarkerSize(.7);
233 0 : fCOG.SetMarkerColor(2);
234 0 : }
235 :
236 :
237 : //_________________________________________________________________
238 : void TKDNodeInfo::TKDNodeDraw::Draw(Option_t* option)
239 : {
240 0 : TBox::Draw(option);
241 0 : fCOG.Draw("p");
242 0 : }
243 :
244 : //_________________________________________________________________
245 : void TKDNodeInfo::TKDNodeDraw::SetNode(TKDNodeInfo *node, UChar_t size, UChar_t ax1, UChar_t ax2)
246 : {
247 0 : fNode=node;
248 : const Float_t kBorder = 0.;//1.E-4;
249 0 : Float_t *bounds = &(node->Data()[size]);
250 0 : fX1=bounds[2*ax1]+kBorder;
251 0 : fX2=bounds[2*ax1+1]-kBorder;
252 0 : fY1=bounds[2*ax2]+kBorder;
253 0 : fY2=bounds[2*ax2+1]-kBorder;
254 :
255 0 : Float_t x(node->Data()[ax1]), y(node->Data()[ax2]);
256 0 : fCOG.SetX(x); fCOG.SetY(y);
257 0 : }
258 :
259 :
260 : //_________________________________________________________________
261 : void TKDNodeInfo::TKDNodeDraw::Print(const Option_t* option) const
262 : {
263 0 : if(!fNode) return;
264 0 : fNode->Print(option);
265 0 : }
|
