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