]>
Commit | Line | Data |
---|---|---|
db0e2c5f | 1 | #include "AliTRDTKDInterpolator.h" |
2 | ||
3 | #include "TClonesArray.h" | |
4 | #include "TLinearFitter.h" | |
5 | #include "TMath.h" | |
6 | #include "TRandom.h" | |
7 | ||
f2762b1c | 8 | #include "AliLog.h" |
9 | ||
db0e2c5f | 10 | #include "iostream" |
11 | using namespace std; | |
12 | ||
13 | ClassImp(AliTRDTKDInterpolator) | |
14 | ||
15 | //_________________________________________________________________ | |
16 | AliTRDTKDInterpolator::AliTRDTKDInterpolator() : | |
17 | TKDTreeIF(), | |
18 | fNDataNodes(0), | |
19 | fNodes(NULL), | |
20 | fLambda(0), | |
21 | fNPointsI(0), | |
22 | fUseHelperNodes(kFALSE), | |
23 | fUseWeights(kFALSE), | |
bd58d4b9 | 24 | fPDFMode(kInterpolation), |
25 | fStoreCov(kFALSE) | |
db0e2c5f | 26 | { |
27 | // default constructor | |
28 | } | |
29 | ||
30 | //_________________________________________________________________ | |
31 | AliTRDTKDInterpolator::AliTRDTKDInterpolator(Int_t npoints, Int_t ndim, UInt_t bsize, Float_t **data) : | |
32 | TKDTreeIF(npoints, ndim, bsize, data), | |
33 | fNDataNodes(0), | |
34 | fNodes(NULL), | |
35 | fLambda(1 + ndim + (ndim*(ndim+1)>>1)), | |
36 | fNPointsI(100), | |
37 | fUseHelperNodes(kFALSE), | |
38 | fUseWeights(kFALSE), | |
bd58d4b9 | 39 | fPDFMode(kInterpolation), |
40 | fStoreCov(kFALSE) | |
db0e2c5f | 41 | { |
42 | } | |
43 | ||
44 | //_________________________________________________________________ | |
45 | AliTRDTKDInterpolator::~AliTRDTKDInterpolator() | |
46 | { | |
47 | if(fNodes){ | |
48 | fNodes->Delete(); | |
49 | delete fNodes; | |
50 | fNodes=NULL; | |
51 | } | |
52 | } | |
53 | //_________________________________________________________________ | |
54 | AliTRDTKDInterpolator::AliTRDTKDInterpolator(const AliTRDTKDInterpolator &ref): | |
55 | TKDTreeIF(), | |
56 | fNDataNodes(ref.fNDataNodes), | |
57 | fNodes(ref.fNodes), | |
58 | fLambda(ref.fLambda), | |
59 | fNPointsI(ref.fNPointsI), | |
60 | fUseHelperNodes(ref.fUseHelperNodes), | |
61 | fUseWeights(ref.fUseWeights), | |
bd58d4b9 | 62 | fPDFMode(ref.fPDFMode), |
63 | fStoreCov(ref.fStoreCov) | |
db0e2c5f | 64 | { |
65 | // Copy constructor | |
f2762b1c | 66 | this->Print(""); |
db0e2c5f | 67 | } |
68 | ||
69 | //____________________________________________________________ | |
70 | ||
71 | AliTRDTKDInterpolator &AliTRDTKDInterpolator::operator=(const AliTRDTKDInterpolator &ref){ | |
72 | // | |
73 | // Assignment operator | |
74 | // | |
75 | if(this == &ref) return *this; | |
76 | ||
77 | // Make copy | |
78 | TObject::operator=(ref); | |
79 | ||
f2762b1c | 80 | this->Print(""); |
81 | ||
db0e2c5f | 82 | return *this; |
83 | } | |
84 | ||
85 | //_________________________________________________________________ | |
86 | Bool_t AliTRDTKDInterpolator::Build() | |
87 | { | |
88 | TKDTreeIF::Build(); | |
89 | if(!fBoundaries) MakeBoundaries(); | |
90 | ||
91 | // allocate interpolation nodes | |
f2762b1c | 92 | fNDataNodes = fNPoints/fBucketSize + ((fNPoints%fBucketSize)?1:0); |
db0e2c5f | 93 | |
94 | if(fNodes){ | |
95 | Warning("AliTRDTKDInterpolator::Build()", "Data already allocated."); | |
96 | fNodes->Delete(); | |
97 | } else { | |
98 | fNodes = new TClonesArray("AliTRDTKDInterpolator::AliTRDTKDNodeInfo", fNDataNodes); | |
99 | fNodes->SetOwner(); | |
100 | } | |
101 | for(int in=0; in<fNDataNodes; in++) new ((*fNodes)[in]) AliTRDTKDNodeInfo(fNDim); | |
102 | ||
103 | // Set Interpolator nodes | |
104 | ||
105 | for(int inode=0, tnode = fNNodes; inode<fNDataNodes-1; inode++, tnode++){ | |
106 | AliTRDTKDNodeInfo *node =GetNodeInfo(inode); | |
107 | memcpy(node->fBounds,GetBoundary(tnode),2*fNDim*sizeof(Float_t)); | |
108 | node->fVal[0] = Float_t(fBucketSize)/fNPoints; | |
109 | for(int idim=0; idim<fNDim; idim++) node->fVal[0] /= (node->fBounds[2*idim+1] - node->fBounds[2*idim]); | |
110 | node->fVal[1] = node->fVal[0]/TMath::Sqrt(float(fBucketSize)); | |
111 | ||
112 | Int_t *indexPoints = GetPointsIndexes(tnode); | |
113 | // loop points in this terminal node | |
114 | for(int idim=0; idim<fNDim; idim++){ | |
115 | node->fData[idim] = 0.; | |
116 | for(int ip = 0; ip<fBucketSize; ip++) node->fData[idim] += fData[idim][indexPoints[ip]]; | |
117 | node->fData[idim] /= fBucketSize; | |
118 | } | |
119 | } | |
120 | ||
121 | // Analyze last (incomplete) terminal node | |
122 | ||
123 | Int_t counts = fNPoints%fBucketSize; | |
124 | counts = counts ? counts : fBucketSize; | |
125 | Int_t inode = fNDataNodes - 1, tnode = inode + fNNodes; | |
126 | AliTRDTKDNodeInfo *ftnode = GetNodeInfo(inode); | |
127 | ftnode->fVal[0] = Float_t(counts)/fNPoints; | |
128 | memcpy(ftnode->fBounds,GetBoundary(tnode),2*fNDim*sizeof(Float_t)); | |
129 | ||
130 | for(int idim=0; idim<fNDim; idim++){ | |
131 | Float_t dx = ftnode->fBounds[2*idim+1]-ftnode->fBounds[2*idim]; | |
132 | if(dx < 1.e-30){ | |
133 | Warning("AliTRDTKDInterpolator::Build()", "Terminal bucket index[%d] too narrow on the %d dimension.", inode, idim); | |
134 | continue; | |
135 | } | |
136 | ftnode->fVal[0] /= (ftnode->fBounds[2*idim+1] - ftnode->fBounds[2*idim]); | |
137 | } | |
138 | ftnode->fVal[1] = ftnode->fVal[0]/TMath::Sqrt(float(counts)); | |
139 | ||
140 | // loop points in this terminal node | |
141 | Int_t *indexPoints = GetPointsIndexes(tnode); | |
142 | for(int idim=0; idim<fNDim; idim++){ | |
143 | ftnode->fData[idim] = 0.; | |
144 | for(int ip = 0; ip<counts; ip++) ftnode->fData[idim] += fData[idim][indexPoints[ip]]; | |
145 | ftnode->fData[idim] /= counts; | |
146 | } | |
147 | ||
148 | delete [] fBoundaries; | |
149 | fBoundaries = NULL; | |
150 | // Add Helper Nodes | |
151 | if(fUseHelperNodes){BuildBoundaryNodes();} | |
152 | ||
153 | if(fNPointsI>GetNTNodes()){fNPointsI=GetNTNodes();} | |
154 | ||
155 | BuildInterpolation(); | |
156 | ||
157 | return kTRUE; | |
158 | } | |
159 | ||
160 | ||
161 | //_________________________________________________________________ | |
162 | Bool_t AliTRDTKDInterpolator::Eval(const Double_t *point, Double_t &result, Double_t &error) | |
163 | { | |
f2762b1c | 164 | AliDebug(3,Form("Eval PDF Mode %d",fPDFMode)); |
165 | if((AliLog::GetDebugLevel("",IsA()->GetName()))>0){ | |
166 | printf("Point ["); | |
167 | for(int idim=0; idim<fNDim; idim++) printf("%f ", point[idim]); | |
168 | printf("] \n"); | |
169 | } | |
170 | ||
db0e2c5f | 171 | Float_t pointF[fNDim]; // local Float_t conversion for "point" |
172 | for(int idim=0; idim<fNDim; idim++) pointF[idim] = (Float_t)point[idim]; | |
173 | Int_t nodeIndex = GetNodeIndex(pointF); | |
174 | if(nodeIndex<0){ | |
f2762b1c | 175 | AliError("Can not retrieve node for data point"); |
db0e2c5f | 176 | result = 0.; |
177 | error = 1.E10; | |
178 | return kFALSE; | |
179 | } | |
180 | AliTRDTKDNodeInfo *node =GetNodeInfo(nodeIndex); | |
181 | ||
f2762b1c | 182 | if((AliLog::GetDebugLevel("",IsA()->GetName()))>0){ |
183 | printf("Node Info: \n"); | |
184 | node->Print("a"); | |
db0e2c5f | 185 | } |
f2762b1c | 186 | |
187 | return node->CookPDF(point, result, error,fPDFMode); | |
db0e2c5f | 188 | } |
189 | ||
190 | //__________________________________________________________________ | |
191 | void AliTRDTKDInterpolator::Print(const Option_t */*opt*/) const | |
192 | { | |
193 | for(Int_t i=GetNTNodes(); i--;){ | |
194 | printf("Node %d of %d: \n",i,GetNTNodes()); | |
195 | GetNodeInfo(i)->Print(); | |
196 | } | |
197 | ||
198 | } | |
199 | ||
200 | //__________________________________________________________________ | |
201 | Int_t AliTRDTKDInterpolator::GetNodeIndex(const Float_t *p) | |
202 | { | |
203 | Int_t inode=FindNode(p)-fNDataNodes+1; | |
f2762b1c | 204 | if(GetNodeInfo(inode)->Has(p)){ |
205 | AliDebug(2,Form("Find Node %d",inode)); | |
206 | return inode; | |
207 | } | |
db0e2c5f | 208 | |
209 | // Search extra nodes | |
210 | ||
211 | for(inode=fNDataNodes;inode<GetNTNodes();inode++){ | |
f2762b1c | 212 | if(GetNodeInfo(inode)->Has(p)){ |
213 | AliDebug(2,Form("Find Extra Node %d",inode)); | |
214 | return inode; | |
215 | } | |
db0e2c5f | 216 | } |
217 | ||
218 | // Search for nearest neighbor | |
219 | Float_t dist; | |
220 | Float_t closestdist=10000; | |
221 | inode=-1; | |
222 | for(Int_t ii=0;ii<GetNTNodes();ii++){ | |
223 | AliTRDTKDNodeInfo *node=GetNodeInfo(ii); | |
224 | dist=0; | |
225 | for(Int_t idim=0;idim<fNDim;idim++){ | |
226 | dist+=TMath::Power((node->fData[idim]-p[idim]),2); | |
227 | } | |
228 | dist=TMath::Sqrt(dist); | |
229 | if(dist<closestdist){closestdist=dist;inode=ii;} | |
230 | } | |
f2762b1c | 231 | AliDebug(2,Form("Find Nearest Neighbor Node %d",inode)); |
db0e2c5f | 232 | return inode; |
233 | } | |
234 | ||
235 | //_________________________________________________________________ | |
236 | AliTRDTKDInterpolator::AliTRDTKDNodeInfo* AliTRDTKDInterpolator::GetNodeInfo(Int_t inode) const | |
237 | { | |
238 | if(!fNodes || inode >= GetNTNodes()) return NULL; | |
239 | return (AliTRDTKDNodeInfo*)(*fNodes)[inode]; | |
240 | } | |
241 | ||
242 | //_________________________________________________________________ | |
243 | Int_t AliTRDTKDInterpolator::GetNTNodes() const | |
244 | { | |
245 | return fNodes?fNodes->GetEntriesFast():0; | |
246 | } | |
247 | ||
248 | //_________________________________________________________________ | |
249 | Bool_t AliTRDTKDInterpolator::GetRange(Int_t idim,Float_t range[2]) const | |
250 | { | |
251 | if(!fNodes) return kFALSE; | |
252 | if(idim<0 || idim>=fNDim){ | |
253 | range[0]=0.; range[1]=0.; | |
254 | return kFALSE; | |
255 | } | |
256 | range[0]=1.e10; range[1]=-1.e10; | |
257 | for(Int_t in=GetNTNodes(); in--; ){ | |
258 | AliTRDTKDNodeInfo *node = GetNodeInfo(in); | |
259 | ||
260 | if(node->fBounds[2*idim]<range[0]) range[0] = node->fBounds[2*idim]; | |
261 | if(node->fBounds[2*idim+1]>range[1]) range[1] = node->fBounds[2*idim+1]; | |
262 | } | |
263 | ||
264 | return kTRUE; | |
265 | } | |
266 | ||
267 | //_________________________________________________________________ | |
268 | TH2Poly *AliTRDTKDInterpolator::Projection(Int_t xdim,Int_t ydim) | |
269 | { | |
270 | Float_t rangex[2],rangey[2]; | |
271 | GetRange(xdim,rangex); | |
272 | GetRange(ydim,rangey); | |
273 | ||
274 | TH2Poly* h2 = new TH2Poly("hTKDnodes","hTKDnodes", rangex[0],rangex[1],rangey[0],rangey[1]); | |
275 | h2->GetXaxis()->SetTitle(Form("Q_{%d}", xdim)); | |
276 | h2->GetYaxis()->SetTitle(Form("Q_{%d}", ydim)); | |
277 | ||
278 | for(Int_t inode=0;inode<GetNTNodes();inode++){ | |
279 | ||
280 | AliTRDTKDNodeInfo* node=GetNodeInfo(inode); | |
281 | h2->AddBin(node->fBounds[2*xdim],node->fBounds[2*ydim],node->fBounds[2*xdim+1],node->fBounds[2*ydim+1]); | |
282 | h2->SetBinContent(inode+1, node->fVal[0]); | |
283 | } | |
284 | return h2; | |
285 | } | |
286 | ||
287 | //_________________________________________________________________ | |
288 | void AliTRDTKDInterpolator::BuildInterpolation() | |
289 | { | |
f2762b1c | 290 | AliInfo("Build Interpolation"); |
db0e2c5f | 291 | |
292 | // Calculate Interpolation | |
293 | ||
294 | Double_t buffer[fLambda]; | |
295 | ||
296 | Float_t **refPoints = new Float_t*[fNDim]; | |
297 | for(int id=0; id<fNDim; id++){ | |
298 | refPoints[id] = new Float_t[GetNTNodes()]; | |
299 | for(int in=0; in<GetNTNodes(); in++) refPoints[id][in] = GetNodeInfo(in)->fData[id]; | |
300 | } | |
301 | TKDTreeIF *KDhelper = new TKDTreeIF(GetNTNodes(), fNDim, 30, refPoints); | |
302 | KDhelper->Build(); | |
303 | KDhelper->MakeBoundariesExact(); | |
304 | ||
305 | Float_t dist[fNPointsI]; | |
306 | Int_t ind[fNPointsI]; | |
307 | ||
308 | TLinearFitter fitter(fLambda, Form("hyp%d", fLambda-1)); | |
309 | ||
310 | Int_t nodeIndex(0); Float_t param[6], *pp(NULL); | |
311 | nodeIndex=GetNTNodes(); pp=¶m[0]; | |
312 | while(nodeIndex--){ | |
313 | ||
314 | fitter.ClearPoints(); | |
315 | ||
316 | AliTRDTKDNodeInfo *node = GetNodeInfo(nodeIndex); | |
317 | // find nearest neighbors | |
318 | KDhelper->FindNearestNeighbors(node->fData,fNPointsI, &ind[0], &dist[0]); | |
319 | ||
320 | for(int in=0;in<fNPointsI;in++){ | |
321 | AliTRDTKDNodeInfo *nnode = GetNodeInfo(ind[in]); | |
322 | ||
323 | Float_t w=1; //weight | |
324 | // calculate tri-cubic weighting function | |
325 | if(fUseWeights){ | |
326 | Float_t d = dist[in]/dist[fNPointsI-1]; | |
327 | Float_t w0 = (1. - d*d*d); | |
328 | w = w0*w0*w0; | |
329 | if(w<1.e-30) continue; | |
330 | } | |
331 | Int_t ipar=0; | |
332 | for(int idim=0; idim<fNDim; idim++){ | |
333 | buffer[ipar++] = nnode->fData[idim]; | |
334 | for(int jdim=idim; jdim<fNDim; jdim++) buffer[ipar++] = nnode->fData[idim]*nnode->fData[jdim]; | |
335 | } | |
336 | fitter.AddPoint(buffer,nnode->fVal[0], nnode->fVal[1]/w); | |
337 | ||
338 | // Ensure Boundary Condition | |
339 | for(Int_t kdim=0;kdim<fNDim;kdim++){ | |
340 | if(node->fBounds[2*kdim]==0){ | |
341 | Float_t zdata[fNDim]; | |
342 | memcpy(&zdata[0],node->fData,fNDim*sizeof(Float_t)); | |
343 | zdata[kdim]=0; | |
344 | ipar=0; | |
345 | for(int idim=0; idim<fNDim; idim++){ | |
346 | buffer[ipar++] = zdata[idim]; | |
347 | for(int jdim=idim; jdim<fNDim; jdim++) buffer[ipar++] = zdata[idim]*zdata[jdim]; | |
348 | } | |
349 | fitter.AddPoint(buffer,0,1); | |
350 | } | |
351 | } | |
352 | } | |
353 | ||
f2762b1c | 354 | AliDebug(2,Form("Calculate Interpolation for Node %d",nodeIndex)); |
db0e2c5f | 355 | fitter.Eval(); |
356 | ||
357 | // retrive fitter results | |
358 | TMatrixD cov(fLambda, fLambda); | |
359 | TVectorD par(fLambda); | |
360 | fitter.GetCovarianceMatrix(cov); | |
361 | fitter.GetParameters(par); | |
362 | ||
363 | // store results | |
bd58d4b9 | 364 | node->Store(&par,&cov,fStoreCov); |
db0e2c5f | 365 | } |
366 | ||
367 | delete KDhelper; | |
368 | for(int id=0; id<fNDim; id++){ | |
369 | delete refPoints[id]; | |
370 | } | |
371 | delete[] refPoints; | |
372 | } | |
373 | ||
374 | //_________________________________________________________________ | |
375 | void AliTRDTKDInterpolator::BuildBoundaryNodes(){ | |
376 | ||
377 | Int_t nnew=0; | |
378 | ||
379 | Float_t treebounds[2*fNDim]; | |
380 | for(Int_t idim=0;idim<fNDim;idim++){ | |
381 | GetRange(idim,&treebounds[2*idim]); | |
382 | } | |
383 | ||
384 | for(int inode=0; inode<GetNTNodes(); inode++){ | |
385 | ||
386 | AliTRDTKDNodeInfo *node=GetNodeInfo(inode); | |
387 | ||
388 | for(Int_t vdim=0;vdim<fNDim;vdim++){ | |
389 | ||
390 | // Try expansion to lower and higher values | |
391 | for(Int_t iter=0;iter<2;iter++){ | |
392 | if(node->fBounds[2*vdim+iter]==treebounds[2*vdim+iter]){ | |
393 | ||
394 | // Add new Node | |
395 | new ((*fNodes)[GetNTNodes()]) AliTRDTKDNodeInfo(fNDim); | |
396 | ||
397 | AliTRDTKDNodeInfo *newnode = GetNodeInfo(GetNTNodes()-1); | |
398 | if(iter==0)newnode->fBounds[2*vdim+iter]=0; | |
399 | if(iter==1)newnode->fBounds[2*vdim+iter]=2*treebounds[2*vdim+iter]; | |
400 | newnode->fBounds[2*vdim+!iter]=node->fBounds[2*vdim+iter]; | |
401 | for(Int_t idim=0;idim<fNDim;idim++){ | |
402 | if(idim==vdim)continue; | |
403 | newnode->fBounds[2*idim]=node->fBounds[2*idim]; | |
404 | newnode->fBounds[2*idim+1]=node->fBounds[2*idim+1]; | |
405 | } | |
406 | newnode->fVal[0]=0; | |
407 | newnode->fVal[1]=Float_t(1)/fNPoints; | |
408 | for(int idim=0; idim<fNDim; idim++){ | |
409 | newnode->fVal[1] /= (newnode->fBounds[2*idim+1] - newnode->fBounds[2*idim]); | |
410 | newnode->fData[idim]=0.5*(newnode->fBounds[2*idim+1] + newnode->fBounds[2*idim]); | |
411 | } | |
412 | nnew++; | |
413 | } | |
414 | } | |
415 | } | |
416 | } | |
f2762b1c | 417 | AliInfo(Form("%d Boundary Nodes Added \n",nnew)); |
db0e2c5f | 418 | } |
419 | ||
420 | //_________________________________________________________________ | |
421 | AliTRDTKDInterpolator::AliTRDTKDNodeInfo::AliTRDTKDNodeInfo(Int_t ndim): | |
422 | TObject() | |
423 | ,fNDim(ndim) | |
424 | ,fNBounds(2*ndim) | |
425 | ,fNPar(1 + ndim + (ndim*(ndim+1)>>1)) | |
f2762b1c | 426 | ,fNCov(Int_t((fNPar+1)*fNPar/2)) |
db0e2c5f | 427 | ,fData(NULL) |
428 | ,fBounds(NULL) | |
429 | ,fPar(NULL) | |
430 | ,fCov(NULL) | |
431 | { | |
432 | // Default constructor | |
433 | fVal[0] = 0.; fVal[1] = 0.; | |
434 | fData=new Float_t[fNDim]; | |
435 | fBounds=new Float_t[fNBounds]; | |
436 | } | |
437 | ||
438 | //_________________________________________________________________ | |
439 | AliTRDTKDInterpolator::AliTRDTKDNodeInfo::AliTRDTKDNodeInfo(const AliTRDTKDInterpolator::AliTRDTKDNodeInfo &ref): | |
440 | TObject(ref) | |
441 | ,fNDim(ref.fNDim) | |
442 | ,fNBounds(ref.fNBounds) | |
443 | ,fNPar(ref.fNPar) | |
444 | ,fNCov(ref.fNCov) | |
445 | ,fData(NULL) | |
446 | ,fBounds(NULL) | |
447 | ,fPar(NULL) | |
448 | ,fCov(NULL) | |
449 | { | |
450 | // Copy constructor | |
451 | ||
452 | if(ref.fData){ | |
453 | fData = new Float_t[fNDim]; | |
454 | memcpy(fData, ref.fData, fNDim*sizeof(Float_t)); | |
455 | } | |
456 | if(ref.fBounds){ | |
457 | fBounds = new Float_t[2*fNDim]; | |
458 | memcpy(fBounds, ref.fBounds, 2*fNDim*sizeof(Float_t)); | |
459 | } | |
460 | ||
461 | fVal[0] = ref.fVal[0]; | |
462 | fVal[1] = ref.fVal[1]; | |
463 | ||
464 | if(ref.fPar){ | |
465 | fPar=new Double_t[fNPar]; | |
466 | memcpy(fPar, ref.fPar, fNPar*sizeof(Double_t)); | |
467 | } | |
468 | if(ref.fCov){ | |
469 | fCov=new Double_t[fNCov]; | |
470 | memcpy(fCov, ref.fCov, fNCov*sizeof(Double_t)); | |
471 | } | |
472 | } | |
473 | ||
474 | //_________________________________________________________________ | |
475 | AliTRDTKDInterpolator::AliTRDTKDNodeInfo::~AliTRDTKDNodeInfo() | |
476 | { | |
477 | // Destructor | |
478 | if(fData) delete [] fData; | |
479 | if(fBounds) delete [] fBounds; | |
480 | if(fPar) delete [] fPar; | |
481 | if(fCov) delete [] fCov; | |
482 | } | |
483 | ||
484 | //_________________________________________________________________ | |
485 | AliTRDTKDInterpolator::AliTRDTKDNodeInfo &AliTRDTKDInterpolator::AliTRDTKDNodeInfo::operator=(const AliTRDTKDInterpolator::AliTRDTKDNodeInfo &ref) | |
486 | { | |
487 | // | |
488 | // Assignment operator | |
489 | // | |
490 | ||
491 | if(&ref==this) return *this; | |
492 | memcpy(fData, ref.fData, fNDim*sizeof(Float_t)); | |
493 | memcpy(fBounds, ref.fBounds, 2*fNDim*sizeof(Float_t)); | |
494 | ||
495 | fVal[0] = ref.fVal[0]; | |
496 | fVal[1] = ref.fVal[1]; | |
497 | ||
498 | if(ref.fPar){ | |
499 | fPar=new Double_t[fNPar]; | |
500 | memcpy(fPar, ref.fPar, fNPar*sizeof(Double_t)); | |
501 | } | |
502 | if(ref.fCov){ | |
503 | fCov=new Double_t[fNCov]; | |
504 | memcpy(fCov, ref.fCov, fNCov*sizeof(Double_t)); | |
505 | } | |
506 | return *this; | |
507 | } | |
508 | ||
509 | ||
510 | //_________________________________________________________________ | |
511 | void AliTRDTKDInterpolator::AliTRDTKDNodeInfo::Print(const Option_t *opt) const | |
512 | { | |
513 | // Print the content of the node | |
514 | printf("x ["); | |
515 | for(int idim=0; idim<fNDim; idim++) printf("%f ", fData?fData[idim]:0.); | |
516 | printf("] f = [%e +- %e]\n", fVal[0], fVal[1]); | |
517 | ||
518 | if(fBounds){ | |
519 | printf("range ["); | |
520 | for(int idim=0; idim<fNDim; idim++) printf("(%f %f) ", fBounds[2*idim], fBounds[2*idim+1]); | |
521 | printf("]\n"); | |
522 | } | |
523 | if(strcmp(opt, "a")!=0) return; | |
524 | ||
525 | if(fPar){ | |
526 | printf("Fit parameters : \n"); | |
f2762b1c | 527 | for(int ip=0; ip<fNPar; ip++) printf("p%d[%e] ", ip, fPar[ip]); |
db0e2c5f | 528 | printf("\n"); |
529 | } | |
530 | if(!fCov) return; | |
531 | for(int ip(0), n(0); ip<fNPar; ip++){ | |
f2762b1c | 532 | for(int jp(ip); jp<fNPar; jp++) printf("c(%d %d)[%e] ", ip, jp, fCov[n++]); |
db0e2c5f | 533 | printf("\n"); |
534 | } | |
535 | } | |
536 | ||
537 | //_________________________________________________________________ | |
bd58d4b9 | 538 | void AliTRDTKDInterpolator::AliTRDTKDNodeInfo::Store(TVectorD const *par, TMatrixD const *cov,Bool_t storeCov) |
db0e2c5f | 539 | { |
f2762b1c | 540 | // Store the parameters and the covariance in the node |
db0e2c5f | 541 | |
f2762b1c | 542 | AliDebug(2,"Store Node Interpolation Parameters"); |
db0e2c5f | 543 | |
f2762b1c | 544 | if((AliLog::GetDebugLevel("",IsA()->GetName()))>=10){ |
545 | par->Print(""); | |
546 | cov->Print(""); | |
547 | } | |
548 | ||
549 | if(!fPar){fPar = new Double_t[fNPar];} | |
550 | for(int ip=0; ip<fNPar; ip++) fPar[ip] = (*par)[ip]; | |
bd58d4b9 | 551 | if(!cov||!storeCov) return; |
f2762b1c | 552 | if(!fCov){fCov = new Double_t[fNCov];} |
553 | for(int ip(0), np(0); ip<fNPar; ip++) | |
554 | for(int jp=ip; jp<fNPar; jp++) fCov[np++] = (*cov)(ip,jp); | |
555 | ||
556 | if((AliLog::GetDebugLevel("",IsA()->GetName()))>10){this->Print("a");} | |
db0e2c5f | 557 | } |
558 | ||
559 | //_________________________________________________________________ | |
f2762b1c | 560 | Bool_t AliTRDTKDInterpolator::AliTRDTKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t &error,TRDTKDMode mod) const |
db0e2c5f | 561 | { |
562 | // Recalculate the PDF for one node from the results of interpolation (parameters and covariance matrix) | |
563 | ||
564 | result =0.; error = 1.; | |
f2762b1c | 565 | |
566 | if(mod==kNodeVal){ | |
567 | error=fVal[1]; | |
568 | result=fVal[0]; | |
569 | return kTRUE; | |
570 | } | |
571 | ||
572 | if(!fPar){ | |
573 | AliDebug(1,"Requested Interpolation Parameters don't exist"); | |
574 | return kFALSE; | |
575 | } | |
db0e2c5f | 576 | |
577 | Double_t fdfdp[fNDim]; | |
578 | Int_t ipar = 0; | |
579 | fdfdp[ipar++] = 1.; | |
580 | for(int idim=0; idim<fNDim; idim++){ | |
581 | fdfdp[ipar++] = point[idim]; | |
582 | for(int jdim=idim; jdim<fNDim; jdim++) fdfdp[ipar++] = point[idim]*point[jdim]; | |
583 | } | |
584 | ||
585 | // calculate estimation | |
586 | for(int i=0; i<fNPar; i++) result += fdfdp[i]*fPar[i]; | |
587 | ||
f2762b1c | 588 | if(!fCov){ |
589 | AliDebug(3,"Interpolation Error cannot be estimated, Covariance Parameters don't exist"); | |
590 | return kTRUE; | |
591 | } | |
db0e2c5f | 592 | // calculate error |
593 | error=0; | |
594 | ||
595 | for(int i(0), n(0); i<fNPar; i++){ | |
596 | error += fdfdp[i]*fdfdp[i]*fCov[n++]; | |
597 | for(int j(i+1); j<fNPar; j++) error += 2.*fdfdp[i]*fdfdp[j]*fCov[n++]; | |
598 | } | |
599 | if(error>0)error = TMath::Sqrt(error); | |
600 | else{error=100;} | |
601 | ||
f2762b1c | 602 | if(mod==kMinError){ |
603 | if(error<fVal[1]){ | |
604 | return kTRUE; | |
605 | } | |
bd58d4b9 | 606 | if(error==1)AliWarning("Covariance not available, always choose bin content"); |
f2762b1c | 607 | error=fVal[1]; |
608 | result=fVal[0]; | |
609 | return kTRUE; | |
610 | } | |
611 | ||
db0e2c5f | 612 | // Boundary condition |
613 | if(result<0){ | |
f2762b1c | 614 | AliDebug(2,"Using Node Value to ensure Boundary Condition"); |
db0e2c5f | 615 | result=fVal[0]; |
616 | error=fVal[1]; | |
617 | } | |
618 | ||
f2762b1c | 619 | AliDebug(1,Form("Cook PDF Result: %e Error: %e",result,error)); |
620 | ||
db0e2c5f | 621 | return kTRUE; |
622 | } | |
623 | ||
624 | //_____________________________________________________________________ | |
625 | Bool_t AliTRDTKDInterpolator::AliTRDTKDNodeInfo::Has(const Float_t *p) const | |
626 | { | |
627 | Int_t n(0); | |
628 | for(int id=0; id<fNDim; id++) | |
629 | if(p[id]>=fBounds[2*id] && p[id]<fBounds[2*id+1]) n++; | |
630 | if(n==fNDim) return kTRUE; | |
631 | return kFALSE; | |
632 | } | |
633 |