Adding interface for delta R, RPhi and Phi
[u/mrichter/AliRoot.git] / TPC / AliTPCkalmanFit.cxx
CommitLineData
0bfb86a0 1/*
2 marian.ivanov@cern.ch
3
4 AliTPCkalmanFit:
5
6 Kalman filter(s) for fitting of the tracks together with calibration/transformation
7 parameters.
8
9 Correction/Transformation are currently described by set (TObjArray) of primitive
10 correction/transformatio - AliTPCTransformation. Currently we assume that transformation
11 comute (in first order). AliTPCTransformation describe general non linear transformation.
12
13 Current calibration parameters and covariance stored (fCalibParam, fCalibCovar).
14
15 Currenly only linear track model implemented.
16 Fits to be implemented:
17 0. Plane fitting (Laser CE)
18 1. Primary vertex fitting.
19 2. Propagation in magnetic field and fit of planes
20
21
22
23 How to use it - see AliTPCkalmanFit::Test function
24
25 Simple test (see AliTPCkalmanFit::Test)
26 AliTPCTransformation::BuildBasicFormulas();
27 AliTPCkalmanFit *kalmanFit0 = AliTPCkalmanFit::Test(2000);
28 TFile f("kalmanfitTest.root");
29
30 Transformation visualization:
31 Transforamtion can be visualized using TFn (TF1,TF2 ...) and using tree->Draw()
32 e.g:
33 kalmanFit0->SetInstance(kalmanFit0); //
34 kalmanFit0->InitTransformation(); //
2497be21 35 //
0e9efcbe 36 TF2 fxRZdz("fxRZdz","sign(y)*1000*(AliTPCkalmanFit::SGetTPCDeltaXYZ(0,-1,x,0,y)-AliTPCkalmanFit::SGetTPCDeltaXYZ(0,-1,x,0,y-1))",85,245,-250,250);
0bfb86a0 37 fxRZdz->Draw("");
38
2497be21 39 TF2 fxRZ("fxRZ","sign(y)*10*(AliTPCkalmanFit::SGetTPCDeltaXYZ(0,-1,x,0,y))",85,245,-250,250);
40 fxRZ->Draw("");
41
0bfb86a0 42
43
44*/
45
46#include "TRandom.h"
47#include "TMath.h"
48#include "TBits.h"
49#include "TFormula.h"
50#include "TF1.h"
51#include "TLinearFitter.h"
52#include "TFile.h"
53#include "THnSparse.h"
54#include "TAxis.h"
55
56
57#include "TTreeStream.h"
58#include "AliTrackPointArray.h"
59#include "AliLog.h"
60#include "AliTPCTransformation.h"
61#include "AliTPCkalmanFit.h"
62
63ClassImp(AliTPCkalmanFit)
64
65AliTPCkalmanFit* AliTPCkalmanFit::fgInstance = 0;
66
67AliTPCkalmanFit::AliTPCkalmanFit():
68 TNamed(),
69 fCalibration(0),
70 fCalibParam(0),
71 fCalibCovar(0),
72 fLinearParam(0),
73 fLinearCovar(0),
74 fLastTimeStamp(-1),
75 fCurrentAlpha(0), //! current rotation frame
76 fCA(0), //! cosine of current angle
77 fSA(0) //! sinus of current angle
78{
79 //
80 // Default constructor
81 //
c64c82db 82 for (Int_t ihis=0; ihis<12; ihis++){
0bfb86a0 83 fLinearTrackDelta[ihis]=0;
84 fLinearTrackPull[ihis]=0;
85 }
86}
87
88void AliTPCkalmanFit::InitTransformation(){
89 //
90 // Initialize pointers to the transforamtion functions
91 //
92 //
93 Int_t ncalibs = fCalibration->GetEntries();
94 for (Int_t icalib=0;icalib<ncalibs; icalib++){
95 AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
96 transform->Init();
97 }
98}
99
0e9efcbe 100void AliTPCkalmanFit::Add(const AliTPCkalmanFit * kalman){
101 //
102 //
103 //
104 Update(kalman);
c64c82db 105 for (Int_t i=0;i<12;i++){
0e9efcbe 106 if (fLinearTrackDelta[i] && kalman->fLinearTrackDelta[i]){
107 fLinearTrackDelta[i]->Add(kalman->fLinearTrackDelta[i]);
108 }
109 if (fLinearTrackPull[i] && kalman->fLinearTrackPull[i]){
110 fLinearTrackPull[i]->Add(kalman->fLinearTrackPull[i]);
111 }
112 }
113
114}
115
116
0bfb86a0 117void AliTPCkalmanFit::Init(){
118 //
119 // Initialize parameter vector and covariance matrix
120 // To be called after initialization of all of the transformations
121 //
122 //
123 Int_t ncalibs = fCalibration->GetEntries();
124 fCalibParam = new TMatrixD(ncalibs,1);
125 fCalibCovar = new TMatrixD(ncalibs,ncalibs);
126 for (Int_t icalib=0;icalib<ncalibs; icalib++){
127 AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
128 (*fCalibParam)(icalib,0) = transform->fParam;
129 for (Int_t jcalib=0;jcalib<ncalibs; jcalib++){
130 if (icalib!=jcalib) (*fCalibCovar)(icalib,jcalib)= 0;
131 if (icalib==jcalib) (*fCalibCovar)(icalib,jcalib) = transform->fSigma*transform->fSigma;
132 }
133 }
134 //
135 // Build QA histograms
136 //
137 Double_t mpi = TMath::Pi();
138 //
139 // delta alpha y0 z0 ky kz
140 Int_t binsQA[6] = {300, 36*4, 30, 25, 20, 20};
141 Double_t xminQA[6] = {-0.5, -mpi, -120, -250, -1, -1.};
142 Double_t xmaxQA[6] = { 0.5, mpi, 120, 250, 1, 1.};
143 TString axisName[6]={"#Delta",
144 "#alpha",
145 "y_{0}",
146 "z_{0}",
147 "tan(#phi)",
148 "tan(theta)"};
149 //
c64c82db 150 TString deltaName[12]={"#Delta_{y}(cm)",
151 "100*#Delta_{#phi}(cm)",
152 "100^{2}dy0^{2}/dx0^{2}(cm)",
153 "100^{2}dy1^{2}/dx1^{2}(cm)",
154 "#Delta_{z}(cm)",
155 "100*#Delta_{#theta}(cm)",
156 "100^{2}*dz0^{2}/dx0^{2}(cm)",
157 "100^{2}*dz1^{2}/dx1^{2}(cm)",
158 "RMSy_{0} (cm)",
159 "RMSy_{1} (cm)",
160 "RMSz_{0} (cm)",
161 "RMSz_{1} (cm)"};
162
163 TString pullName[12]={"#Delta_{y}(unit)",
164 "100*#Delta_{#phi}(unit)",
165 "100^{2}dy0^{2}/dx0^{2}(unit)",
166 "100^{2}dy1^{2}/dx1^{2}(unit)",
167 "#Delta_{z}(unit)",
168 "100*#Delta_{#theta}(unit)",
169 "100^{2}*dz0^{2}/dx0^{2}(unit)",
170 "100^{2}*dz1^{2}/dx1^{2}(unit)"
171 "RMSy_{0} (cm)",
172 "RMSy_{1} (cm)",
173 "RMSz_{0} (cm)",
174 "RMSz_{1} (cm)"};
175 //
176 //
177 //
178 for (Int_t ihis=0; ihis<12; ihis++){
0bfb86a0 179 fLinearTrackDelta[ihis]=0;
180 fLinearTrackPull[ihis]=0;
181 xminQA[0]=-0.5; xmaxQA[0] = 0.5;
182 fLinearTrackDelta[ihis] = new THnSparseS(deltaName[ihis],deltaName[ihis], 6, binsQA,xminQA, xmaxQA);
183 xminQA[0]=-10; xmaxQA[0] = 10;
184 fLinearTrackPull[ihis] = new THnSparseS(pullName[ihis],pullName[ihis], 6, binsQA,xminQA, xmaxQA);
185 for (Int_t iaxis=1; iaxis<6; iaxis++){
186 fLinearTrackDelta[ihis]->GetAxis(iaxis)->SetName(axisName[iaxis]);
187 fLinearTrackDelta[ihis]->GetAxis(iaxis)->SetTitle(axisName[iaxis]);
188 fLinearTrackPull[ihis]->GetAxis(iaxis)->SetName(axisName[iaxis]);
189 fLinearTrackPull[ihis]->GetAxis(iaxis)->SetTitle(axisName[iaxis]);
190 }
191 fLinearTrackDelta[ihis]->GetAxis(0)->SetName(deltaName[ihis]);
192 fLinearTrackDelta[ihis]->GetAxis(0)->SetTitle(deltaName[ihis]);
193 fLinearTrackPull[ihis]->GetAxis(0)->SetName(deltaName[ihis]);
194 fLinearTrackPull[ihis]->GetAxis(0)->SetTitle(deltaName[ihis]);
195 }
196
197
198}
199
0e9efcbe 200void AliTPCkalmanFit::SetStatus(const char * mask, Bool_t setOn, Bool_t isOr){
201 //
202 // 0. To activate all transforamtion call SetStatus(0,kTRUE)
203 // 1. To disable everything SetStatus(0,kFALSE)
204 // 2. To activate/desactivate SetStatus("xxx",kTRUE/kFALSE,kFALSE)
205 Int_t ncalibs = fCalibration->GetEntries();
206 if (mask==0) {
207 for (Int_t i=0; i<ncalibs;i++){
208 AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(i);
209 transform->SetActive(setOn);
210 }
211 }
212 else{
213 for (Int_t i=0; i<ncalibs;i++){
214 AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(i);
215 TString strName(transform->GetName());
216 if (strName.Contains(mask)){
217 if (!isOr) transform->SetActive( transform->IsActive() && setOn);
218 if (isOr) transform->SetActive( transform->IsActive() || setOn);
219 }
220 }
221 }
222}
223
224
225void AliTPCkalmanFit::Update(const AliTPCkalmanFit * kalman){
226 //
227 // Update Kalman filter
228 //
229 Int_t ncalibs = fCalibration->GetEntries();
230 TMatrixD vecXk=*fCalibParam; // X vector
231 TMatrixD covXk=*fCalibCovar; // X covariance
232 TMatrixD &vecZk = *(kalman->fCalibParam);
233 TMatrixD &measR = *(kalman->fCalibCovar);
234
235 TMatrixD matHk(ncalibs,ncalibs); // vector to mesurement
236 TMatrixD vecYk(ncalibs,1); // Innovation or measurement residual
237 TMatrixD matHkT(ncalibs,ncalibs); // helper matrix Hk transpose
238 TMatrixD matSk(ncalibs,ncalibs); // Innovation (or residual) covariance
239 TMatrixD matKk(ncalibs,ncalibs); // Optimal Kalman gain
240 TMatrixD covXk2(ncalibs,ncalibs); // helper matrix
241 TMatrixD covXk3(ncalibs,ncalibs); // helper matrix
242 //
243 for (Int_t i=0;i<ncalibs;i++){
244 for (Int_t j=0;j<ncalibs;j++) matHk(i,j)=0;
245 matHk(i,i)=1;
246 }
247 vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
248 matHkT=matHk.T(); matHk.T();
249 matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
250 matSk.Invert();
251 matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
252 vecXk += matKk*vecYk; // updated vector
253 covXk2= (matHk-(matKk*matHk));
254 covXk3 = covXk2*covXk;
255 covXk = covXk3;
256 (*fCalibParam) = vecXk;
257 (*fCalibCovar) = covXk;
258}
259
260
261
262
0bfb86a0 263
2497be21 264void AliTPCkalmanFit::FitTrackLinear(AliTrackPointArray &points, TTreeSRedirector *debug, Float_t scalingRMSY, Float_t scalingRMSZ){
0bfb86a0 265 //
266 //
267
268 if (!fCalibParam) {
269 AliError("Kalman Fit not initialized");
270 return;
271 }
272 //
273 // 1. Make initial track hypothesis
274 //
275 TLinearFitter lfitY(2,"pol1");
276 TLinearFitter lfitZ(2,"pol1");
277 TVectorD vecZ(2);
278 TVectorD vecY(2);
279 //
280 lfitY.StoreData(kTRUE);
281 lfitZ.StoreData(kTRUE);
282 lfitY.ClearPoints();
283 lfitZ.ClearPoints();
284 Int_t npoints = points.GetNPoints();
285 if (npoints<2) return;
286 const Double_t kFac=npoints*npoints*100;
287 const Double_t kOff0=0.01*0.01;
288 const Double_t kOff1=kOff0/(250.*250.);
289 //
290 // 0. Make seed
291 //
292 // AliTrackPointArray pointsc(points);
293 //ApplyCalibration(&pointsc, -1.);
294 //
295 // 1.a Choosing reference rotation alpha
296 //
297 fCurrentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]);
298 fCA = TMath::Cos(fCurrentAlpha);
299 fSA = TMath::Sin(fCurrentAlpha);
300 //
301 // 1.b Fit the track in the rotated frame - MakeSeed
302 //
c64c82db 303 for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
0bfb86a0 304 Double_t rx = fCA*points.GetX()[ipoint]+fSA*points.GetY()[ipoint];
305 Double_t ry = -fSA*points.GetX()[ipoint]+fCA*points.GetY()[ipoint];
306 Double_t rz = points.GetZ()[ipoint];
307 lfitY.AddPoint(&rx,ry,1);
308 lfitZ.AddPoint(&rx,rz,1);
309 }
310 lfitY.Eval();
311 lfitZ.Eval();
312 lfitY.GetParameters(vecY);
313 lfitZ.GetParameters(vecZ);
314 //
315 // 2. Set initial parameters and the covariance matrix
316 //
317 Int_t ncalibs = fCalibration->GetEntries();
318 if (!fLinearParam) {
319 fLinearParam = new TMatrixD(ncalibs+4,1);
320 fLinearCovar = new TMatrixD(ncalibs+4,ncalibs+4);
321 }
322 //
323 for (Int_t i=0; i<ncalibs+4;i++){
324 (*fLinearParam)(i,0)=0;
325 if (i<ncalibs) (*fLinearParam)(i,0) = (*fCalibParam)(i,0);
326 for (Int_t j=0; j<ncalibs+4;j++){
327 (*fLinearCovar)(i,j)=0;
328 if (i<ncalibs&&j<ncalibs) (*fLinearCovar)(i,j) = (*fCalibCovar)(i,j);
329 }
330 }
331 Double_t chi2Y = lfitY.GetChisquare()/lfitY.GetNpoints();
332 Double_t chi2Z = lfitZ.GetChisquare()/lfitZ.GetNpoints();
333 (*fLinearParam)(ncalibs+0,0) = lfitY.GetParameter(0);
334 (*fLinearCovar)(ncalibs+0,ncalibs+0)= lfitY.GetCovarianceMatrixElement(0,0)*chi2Y*kFac+kOff0;
335 (*fLinearParam)(ncalibs+1,0) = lfitY.GetParameter(1);
336 (*fLinearCovar)(ncalibs+1,ncalibs+1)= lfitY.GetCovarianceMatrixElement(1,1)*chi2Y*kFac+kOff1;
337 //
338 (*fLinearParam)(ncalibs+2,0) = lfitZ.GetParameter(0);
339 (*fLinearCovar)(ncalibs+2,ncalibs+2)= lfitZ.GetCovarianceMatrixElement(0,0)*chi2Z*kFac+kOff0;
340 (*fLinearParam)(ncalibs+3,0) = lfitZ.GetParameter(1);
341 (*fLinearCovar)(ncalibs+3,ncalibs+3)= lfitZ.GetCovarianceMatrixElement(1,1)*chi2Z*kFac+kOff1;
342 //
343 // Fit thetrack together with correction
344 //
345 AliTrackPoint point;
c64c82db 346 for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
0bfb86a0 347 //
348 if (!points.GetPoint(point,ipoint)) continue;
349 Double_t erry2 = chi2Y;
350 Double_t errz2 = chi2Z;
351 //set error - it is hack
352 Float_t *cov = (Float_t*) point.GetCov();
2497be21 353 cov[1] = (erry2+kOff0)*scalingRMSY;
354 cov[2] = (errz2+kOff0)*scalingRMSZ;
0bfb86a0 355 UpdateLinear(point,debug);
356 if (!points.GetPoint(point,npoints-1-ipoint)) continue;
357 //set error - it is hack
358 cov = (Float_t*) point.GetCov();
2497be21 359 cov[1] = (erry2+kOff0)*scalingRMSY;
360 cov[2] = (errz2+kOff0)*scalingRMSZ;
0bfb86a0 361 UpdateLinear(point,debug);
362 }
363 //
364 // save current param and covariance
365 for (Int_t i=0; i<ncalibs;i++){
366 AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(i);
367 transform->fParam= (*fLinearParam)(i,0);
368 (*fCalibParam)(i,0) = (*fLinearParam)(i,0);
369 for (Int_t j=0; j<ncalibs;j++){
370 (*fCalibCovar)(i,j) = (*fLinearCovar)(i,j);
371 }
372 }
373 if (debug){ // dump debug info
374 (*debug)<<"fitLinear"<<
375 "vecY.="<<&vecY<<
376 "vecZ.="<<&vecZ<<
377 "chi2Y="<<chi2Y<<
378 "chi2Z="<<chi2Z<<
379 "lP.="<<fLinearParam<<
380 "cP.="<<fCalibParam<<
381 "lC.="<<fLinearCovar<<
382 "cC.="<<fCalibCovar<<
383 "\n";
384 }
385}
386
387void AliTPCkalmanFit::DumpTrackLinear(AliTrackPointArray &points, TTreeSRedirector *debug){
388 //
389 // Dump the track parameters before and after current calibration
390 //
391 // Track is divided to two parts -
392 // X mean is defined as middle point
393 //
394
395 if (!fCalibParam) {
396 AliError("Kalman Fit not initialized");
397 return;
398 }
399 //
400 //
401 //
402 TLinearFitter *fitters[16];
403 TVectorD *params[16];
404 TVectorD *errs[16];
405 TVectorD chi2N(16);
406 Int_t npoints = points.GetNPoints();
407 AliTrackPointArray pointsTrans(points);
408 ApplyCalibration(&pointsTrans,-1.);
409 for (Int_t ifit=0; ifit<8;ifit++){
410 fitters[ifit] = new TLinearFitter(2,"pol1");
411 params[ifit] = new TVectorD(2);
412 fitters[ifit+8]= new TLinearFitter(3,"pol2");
413 params[ifit+8] = new TVectorD(3);
414 errs[ifit] = new TVectorD(2);
415 errs[ifit+8] = new TVectorD(3);
416 }
417 //
418 // calculate mean x point and corrdinate frame
419 //
420 fCurrentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]);
421 fCA = TMath::Cos(fCurrentAlpha);
422 fSA = TMath::Sin(fCurrentAlpha);
423 Double_t xmean=0, sum=0;
424 for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
425 Double_t rx = fCA*points.GetX()[ipoint]+fSA*points.GetY()[ipoint];
426 xmean+=rx;
427 sum++;
428 }
429 xmean/=sum;
430 //
431 for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
432 Double_t rx = fCA*points.GetX()[ipoint]+fSA*points.GetY()[ipoint];
433 Double_t ry = -fSA*points.GetX()[ipoint]+fCA*points.GetY()[ipoint];
434 Double_t rz = points.GetZ()[ipoint];
435 //
436 Double_t rxT = fCA*pointsTrans.GetX()[ipoint]+fSA*pointsTrans.GetY()[ipoint];
437 Double_t ryT = -fSA*pointsTrans.GetX()[ipoint]+fCA*pointsTrans.GetY()[ipoint];
438 Double_t rzT = pointsTrans.GetZ()[ipoint];
439 if (rx>xmean){
440 fitters[0]->AddPoint(&rxT, ryT,1);
441 fitters[2]->AddPoint(&rxT, rzT,1);
442 fitters[4]->AddPoint(&rx, ry,1);
443 fitters[6]->AddPoint(&rx, rz,1);
444 fitters[8]->AddPoint(&rxT, ryT,1);
445 fitters[10]->AddPoint(&rxT, rzT,1);
446 fitters[12]->AddPoint(&rx, ry,1);
447 fitters[14]->AddPoint(&rx, rz,1);
448 }else{
449 fitters[1]->AddPoint(&rxT, ryT,1);
450 fitters[3]->AddPoint(&rxT, rzT,1);
451 fitters[5]->AddPoint(&rx, ry,1);
452 fitters[7]->AddPoint(&rx, rz,1);
453 fitters[9]->AddPoint(&rxT, ryT,1);
454 fitters[11]->AddPoint(&rxT, rzT,1);
455 fitters[13]->AddPoint(&rx, ry,1);
456 fitters[15]->AddPoint(&rx, rz,1);
457 }
458 }
459 for (Int_t ifit=0;ifit<16;ifit++){
460 fitters[ifit]->Eval();
461 fitters[ifit]->GetParameters(*(params[ifit]));
462 fitters[ifit]->GetErrors(*(errs[ifit]));
463 chi2N[ifit] = TMath::Sqrt(fitters[ifit]->GetChisquare()/(fitters[ifit]->GetNpoints()-1));
464 (*(errs[ifit]))[0]*=chi2N[ifit];
465 (*(errs[ifit]))[1]*=chi2N[ifit];
466 if (ifit>=8) (*(errs[ifit]))[2]*=chi2N[ifit]; // second derivative
467 }
468 if (debug){
469 (*debug)<<"dumpLinear"<<
470 "alpha="<<fCurrentAlpha<<
471 "xmean="<<xmean<<
472 "y0T.="<<params[0]<<
473 "y1T.="<<params[1]<<
474 "z0T.="<<params[2]<<
475 "z1T.="<<params[3]<<
476 "y0O.="<<params[4]<<
477 "y1O.="<<params[5]<<
478 "z0O.="<<params[6]<<
479 "z1O.="<<params[7]<<
480 "y0T2.="<<params[8]<<
481 "y1T2.="<<params[9]<<
482 "z0T2.="<<params[10]<<
483 "z1T2.="<<params[11]<<
484 "y0O2.="<<params[12]<<
485 "y1O2.="<<params[13]<<
486 "z0O2.="<<params[14]<<
487 "z1O2.="<<params[15]<<
488 "y0TErr.="<<errs[0]<<
489 "y1TErr.="<<errs[1]<<
490 "z0TErr.="<<errs[2]<<
491 "z1TErr.="<<errs[3]<<
492 "y0OErr.="<<errs[4]<<
493 "y1OErr.="<<errs[5]<<
494 "z0OErr.="<<errs[6]<<
495 "z1OErr.="<<errs[7]<<
496 "y0T2Err.="<<errs[8]<<
497 "y1T2Err.="<<errs[9]<<
498 "z0T2Err.="<<errs[10]<<
499 "z1T2Err.="<<errs[11]<<
500 "y0O2Err.="<<errs[12]<<
501 "y1O2Err.="<<errs[13]<<
502 "z0O2Err.="<<errs[14]<<
503 "z1O2Err.="<<errs[15]<<
504 "chi2.="<<&chi2N<<
505 "\n";
506 }
507 //
508 // delta alpha y0 z0 ky kz
509 Double_t x[6]={0,0,0,0,0,0};
510 x[1]=fCurrentAlpha;
511 x[2]=(*params[0])[0];
512 x[3]=(*params[2])[0];
513 x[4]=(*params[0])[1];
514 x[5]=(*params[2])[1];
515 //
516 //Delta y
517 //
518 x[0]= (*params[1])[0]-(*params[0])[0];
519 fLinearTrackDelta[0]->Fill(x);
520 x[0]/=TMath::Sqrt((*errs[1])[0]*(*errs[1])[0]+(*errs[0])[0]*(*errs[0])[0]);
521 fLinearTrackPull[0]->Fill(x);
522 //
523 //Delta ky
524 //
525 x[0]= 100*((*params[1])[1]-(*params[0])[1]);
526 fLinearTrackDelta[1]->Fill(x);
527 x[0]/=100*TMath::Sqrt((*errs[1])[1]*(*errs[1])[1]+(*errs[0])[1]*(*errs[0])[1]);
528 fLinearTrackPull[1]->Fill(x);
529 //
530 // ky2_0
531 //
532 x[0]= 100*100*((*params[8])[2]);
533 fLinearTrackDelta[2]->Fill(x);
534 x[0]/=100*100*TMath::Sqrt((*errs[8])[2]*(*errs[8])[2]);
535 fLinearTrackPull[2]->Fill(x);
536 //
537 // ky2_1
538 //
539 x[0]= 100*100*((*params[9])[2]);
540 fLinearTrackDelta[3]->Fill(x);
541 x[0]/=100*100*TMath::Sqrt((*errs[9])[2]*(*errs[9])[2]);
542 fLinearTrackPull[3]->Fill(x);
543 //
544 //
545 //Delta z
546 //
547 x[0]= (*params[3])[0]-(*params[2])[0];
548 fLinearTrackDelta[4]->Fill(x);
549 x[0]/=TMath::Sqrt((*errs[3])[0]*(*errs[3])[0]+(*errs[2])[0]*(*errs[2])[0]);
550 fLinearTrackPull[4]->Fill(x);
551 //
552 //Delta kz
553 //
554 x[0]= 100*((*params[3])[1]-(*params[2])[1]);
555 fLinearTrackDelta[5]->Fill(x);
556 x[0]/=100*TMath::Sqrt((*errs[3])[1]*(*errs[3])[1]+(*errs[2])[1]*(*errs[2])[1]);
557 fLinearTrackPull[5]->Fill(x);
558 //
559 // kz2_0
560 //
561 x[0]= 100*100*((*params[10])[2]);
562 fLinearTrackDelta[6]->Fill(x);
563 x[0]/=100*100*TMath::Sqrt((*errs[10])[2]*(*errs[10])[2]);
564 fLinearTrackPull[6]->Fill(x);
565 //
566 // kz2_1
567 //
568 x[0]= 100*100*((*params[11])[2]);
569 fLinearTrackDelta[7]->Fill(x);
570 x[0]/=100*100*TMath::Sqrt((*errs[11])[2]*(*errs[11])[2]);
571 fLinearTrackPull[7]->Fill(x);
572
c64c82db 573 //
574 // rms of track
575 //
576 x[0]= chi2N[0];
577 fLinearTrackDelta[8]->Fill(x);
578 x[0]= chi2N[1];
579 fLinearTrackDelta[9]->Fill(x);
580 x[0]= chi2N[2];
581 fLinearTrackDelta[10]->Fill(x);
582 x[0]= chi2N[3];
583 fLinearTrackDelta[11]->Fill(x);
584 //
0bfb86a0 585
586
587 for (Int_t ifit=0; ifit<8;ifit++){
588 delete fitters[ifit];
589 delete params[ifit];
590 delete fitters[ifit+8];
591 delete params[ifit+8];
592 delete errs[ifit];
593 delete errs[ifit+8];
594 }
595}
596
597
598void AliTPCkalmanFit::Propagate(TTreeSRedirector */*debug*/){
599 //
600 // Propagate the Kalman
601 //
602}
603
2497be21 604void AliTPCkalmanFit::AddCovariance(const char * varName, Double_t sigma){
605 //
606 //
607 //
608 if (fCalibCovar) return;
609 if (!fCalibration) return;
610 if (!fCalibration->FindObject(varName)) return;
611 Int_t ncalibs = fCalibration->GetEntries();
612 TString strVar(varName);
613 for (Int_t icalib=0;icalib<ncalibs; icalib++){
614 AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
615 if (strVar.CompareTo(transform->GetName())==0){
616 (*fCalibCovar)(icalib,icalib)+=sigma*sigma;
617 }
618 }
619}
620
0bfb86a0 621
622void AliTPCkalmanFit::PropagateTime(Int_t time){
623 //
624 // Propagate the calibration in time
625 // - Increase covariance matrix
626 //
627 if (fCalibCovar) return;
628 Int_t ncalibs = fCalibration->GetEntries();
629 Double_t deltaT = 0;
630 if (fLastTimeStamp>0) deltaT = (fLastTimeStamp-time)/(60*60.); // delta T in hours
631 fLastTimeStamp = time;
632 for (Int_t icalib=0;icalib<ncalibs; icalib++){
633 AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
634 (*fCalibCovar)(icalib,icalib)+= transform->fSigma2Time*deltaT;
635 }
636}
637
638
639void AliTPCkalmanFit::UpdateLinear(AliTrackPoint &point, TTreeSRedirector *debug){
640 //
641 // Update Kalman
642 //
643 //
644 Int_t ncalibs = fCalibration->GetEntries();
645 Int_t kNmeas = 2;
646 Int_t nelem = ncalibs+4;
647 TMatrixD &vecXk=*fLinearParam; // X vector
648 TMatrixD &covXk=*fLinearCovar; // X covariance
649 //
650 TMatrixD matHk(kNmeas,nelem); // vector to mesurement
651 TMatrixD vecYk(kNmeas,1); // Innovation or measurement residual
652 TMatrixD vecZk(kNmeas,1); // Innovation or measurement residual
653 TMatrixD measR(kNmeas,kNmeas);
654 TMatrixD matHkT(nelem,kNmeas); // helper matrix Hk transpose
655 TMatrixD matSk(kNmeas,kNmeas); // Innovation (or residual) covariance
656 TMatrixD matKk(nelem,kNmeas); // Optimal Kalman gain
657 TMatrixD covXk2(nelem,nelem); // helper matrix
658 TMatrixD covXk3(nelem,nelem); // helper matrix
659 TMatrixD mat1(nelem,nelem);
660 //
661 // reset matHk
662 for (Int_t iel=0;iel<nelem;iel++){
663 for (Int_t ip=0;ip<kNmeas;ip++) {
664 matHk(ip,iel)=0;
665 }
666 }
667 for (Int_t iel0=0;iel0<nelem;iel0++)
668 for (Int_t iel1=0;iel1<nelem;iel1++){
669 if (iel0!=iel1) mat1(iel0,iel1)=0;
670 if (iel0==iel1) mat1(iel0,iel1)=1;
671 }
672 //
673 // fill matrix Hk
674 //
675 Int_t volId = point.GetVolumeID();
676 Double_t gxyz[3]={point.GetX(), point.GetY(),point.GetZ()};
677 Double_t rxyz[3]={fCA*gxyz[0]+fSA*gxyz[1],-fSA*gxyz[0]+fCA*gxyz[1] ,point.GetZ()};
678 //
679 Double_t dxdydz[3]={0,0,0};
680 Double_t rdxdydz[3]={0,0,0};
681
682 for (Int_t icalib=0;icalib<ncalibs;icalib++){
683 AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
684 dxdydz[0] = transform->GetDeltaXYZ(0,volId, 1, gxyz[0], gxyz[1], gxyz[2]);
685 dxdydz[1] = transform->GetDeltaXYZ(1,volId, 1, gxyz[0], gxyz[1], gxyz[2]);
686 dxdydz[2] = transform->GetDeltaXYZ(2,volId, 1, gxyz[0], gxyz[1], gxyz[2]);
687 rdxdydz[0] = fCA*dxdydz[0]+fSA*dxdydz[1];
688 rdxdydz[1] = -fSA*dxdydz[0]+fCA*dxdydz[1];
689 rdxdydz[2] = dxdydz[2];
690 //
691 matHk(0,icalib)= rdxdydz[1]-rdxdydz[0]*(*fLinearParam)(ncalibs+1,0); // shift y + shift x * angleY
692 matHk(1,icalib)= rdxdydz[2]-rdxdydz[0]*(*fLinearParam)(ncalibs+3,0); // shift z + shift x * angleZ
693 }
694 matHk(0,ncalibs+0)=1;
695 matHk(0,ncalibs+1)=rxyz[0];
696 matHk(1,ncalibs+2)=1;
697 matHk(1,ncalibs+3)=rxyz[0];
698 //
699 //
700 //
701 vecZk(0,0) = rxyz[1];
702 vecZk(1,0) = rxyz[2];
703 measR(0,0) = point.GetCov()[1]; measR(0,1)=0;
704 measR(1,1) = point.GetCov()[2]; measR(1,0)=0;
c64c82db 705 //
0bfb86a0 706 vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
707 matHkT=matHk.T(); matHk.T();
708 matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
709 matSk.Invert();
710 matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
711 //
712 covXk2= (mat1-(matKk*matHk));
713 covXk3 = covXk2*covXk;
714 //
715 CheckCovariance(covXk3,100);
716 vecXk += matKk*vecYk; // updated vector
717 covXk = covXk3;
718 if (debug){ // dump debug info
719 (*debug)<<"updateLinear"<<
720 //
721 "point.="<<&point<<
722 "vecYk.="<<&vecYk<<
723 "vecZk.="<<&vecZk<<
724 "measR.="<<&measR<<
725 "matHk.="<<&matHk<<
726 "matHkT.="<<&matHkT<<
727 "matSk.="<<&matSk<<
728 "matKk.="<<&matKk<<
729 "covXk2.="<<&covXk2<<
730 "covXk.="<<&covXk<<
731 "vecXk.="<<&vecXk<<
732 "\n";
733 }
734}
735
736
ddf81107 737AliTrackPointArray * AliTPCkalmanFit::SortPoints(AliTrackPointArray &points){
738 //
739 //Creates the array - points sorted according radius - neccessay for kalman fit
740 //
741 //
742 // 0. choose the frame - rotation angle
743 //
744 Int_t npoints = points.GetNPoints();
745 if (npoints<1) return 0;
746 Double_t currentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]);
747 Double_t ca = TMath::Cos(currentAlpha);
748 Double_t sa = TMath::Sin(currentAlpha);
749 //
750 // 1. sort the points
751 //
752 Double_t *rxvector = new Double_t[npoints];
753 Int_t *indexes = new Int_t[npoints];
754 for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
755 rxvector[ipoint]=ca*points.GetX()[ipoint]+sa*points.GetY()[ipoint];
756 }
757 TMath::Sort(npoints, rxvector,indexes,kFALSE);
758 AliTrackPoint point;
759 AliTrackPointArray *pointsSorted= new AliTrackPointArray(npoints);
760 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
761 if (!points.GetPoint(point,indexes[ipoint])) continue;
762 pointsSorted->AddPoint(ipoint,&point);
763 }
764 delete [] rxvector;
765 delete [] indexes;
766 return pointsSorted;
767}
768
769
770
0bfb86a0 771AliTrackPointArray * AliTPCkalmanFit::MakePointArrayLinear(Double_t alpha, Double_t y0, Double_t z0, Double_t ky, Double_t kz, Double_t err){
772 //
773 //
774 //
775 AliTrackPointArray array(500);
776 Float_t cov[10]; // dummy covariance
777 Int_t npoints=0;
778 for (Int_t i=0;i<6;i++) cov[i]=0.001;
779 for (Int_t i=0;i<500;i++){
780 AliTrackPoint point(0, 0, 0, cov, -1,0,0);
781 array.AddPoint(npoints, &point);
782 npoints++;
783 }
784 npoints=0;
785 for (Float_t ir = -245; ir<245; ir++){
786 //
787 //
788 if (TMath::Abs(ir)<80) continue;
789 Double_t ca = TMath::Cos(alpha);
790 Double_t sa = TMath::Sin(alpha);
791 Double_t lx = ir;
792 Double_t ly = y0+lx*ky+gRandom->Gaus(0,err);
793 Double_t lz = z0+lx*kz+gRandom->Gaus(0,err);
794 Double_t gx = lx*ca-ly*sa;
795 Double_t gy = lx*sa+ly*ca;
796 Double_t gz = lz;
797 Double_t galpha = TMath::ATan2(gy,gx);
798 Int_t isec = TMath::Nint((9*galpha/TMath::Pi()+0.5));
799 if (isec<0) isec+=18;
800 if (gz<0) isec+=18;
801 if (ir>130) isec+=36;
802 //
803 AliTrackPoint point(gx, gy, gz, cov, isec,0,0);
804 array.AddPoint(npoints, &point);
805 npoints++;
806 }
807 AliTrackPointArray *parray = new AliTrackPointArray(npoints);
808 AliTrackPoint point;
809 for (Int_t i=0;i<npoints;i++){
810 array.GetPoint(point,i);
811 parray->AddPoint(i,&point);
812 }
813
814 return parray;
815}
816
817void AliTPCkalmanFit::AddCalibration(AliTPCTransformation * calib){
818 //
819 // Add calibration
820 //
821 if (!fCalibration) fCalibration = new TObjArray(2000);
822 fCalibration->AddLast(calib);
823}
824
825void AliTPCkalmanFit::ApplyCalibration(AliTrackPointArray *array, Double_t csign){
826 //
827 //
828 //
829 if (!fCalibration) return;
830 Int_t ncalibs = fCalibration->GetEntries();
831 if (ncalibs==0) return;
832 Int_t npoints = array->GetNPoints();
833 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
834 Int_t volId = array->GetVolumeID()[ipoint];
835 Double_t xyz[3]={array->GetX()[ipoint], array->GetY()[ipoint],array->GetZ()[ipoint]};
836 Double_t dxdydz[3]={0,0,0};
837 for (Int_t icalib=0; icalib<ncalibs; icalib++){
838 AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
839 dxdydz[0] += transform->GetDeltaXYZ(0,volId, transform->fParam, xyz[0], xyz[1], xyz[2]);
840 dxdydz[1] += transform->GetDeltaXYZ(1,volId, transform->fParam, xyz[0], xyz[1], xyz[2]);
841 dxdydz[2] += transform->GetDeltaXYZ(2,volId, transform->fParam, xyz[0], xyz[1], xyz[2]);
842 }
843 ((Float_t*)array->GetX())[ipoint]+=csign*dxdydz[0];
844 ((Float_t*)array->GetY())[ipoint]+=csign*dxdydz[1];
845 ((Float_t*)array->GetZ())[ipoint]+=csign*dxdydz[2];
846 }
847}
848
c64c82db 849Bool_t AliTPCkalmanFit::DumpCorelation(Double_t threshold, const char *mask0, const char *mask1){
0bfb86a0 850 //
851 //
852 //
853 TMatrixD &mat = *fCalibCovar;
854 Int_t nrow= mat.GetNrows();
855 for (Int_t irow=0; irow<nrow; irow++){
c64c82db 856 AliTPCTransformation * trans0 = GetTransformation(irow);
857 TString strName0(trans0->GetName());
858 if (mask0){
859 if (!strName0.Contains(mask0)) continue;
860 }
0bfb86a0 861 for (Int_t icol=irow+1; icol<nrow; icol++){
c64c82db 862 AliTPCTransformation * trans1 = GetTransformation(icol);
863 TString strName1(trans1->GetName());
864 if (mask1){
865 if (!strName1.Contains(mask1)) continue;
866 }
867 //
0bfb86a0 868 Double_t diag = TMath::Sqrt(TMath::Abs(mat(irow,irow)*mat(icol,icol)));
869 if (diag<=0){
870 printf("Negative covariance\t%d\t%d\t%f\n",irow,icol, mat(irow,icol));
871 continue;
872 }
873 Double_t corr0 = mat(irow,icol)/diag;
874 if (TMath::Abs(corr0)>threshold){
0bfb86a0 875 printf("%d\t%d\t%s\t%s\t%f\t%f\t%f\n", irow,icol, trans0->GetName(), trans1->GetName(),
876 TMath::Sqrt(mat(irow,irow)), TMath::Sqrt(mat(icol,icol)), corr0);
877 }
878 }
879 }
2497be21 880 return (nrow>0);
0bfb86a0 881}
882
2497be21 883Bool_t AliTPCkalmanFit::DumpCalib(const char *mask, Float_t correlationCut){
0bfb86a0 884 //
885 // Print calibration entries - name, value, error
886 //
887 TMatrixD &mat = *fCalibCovar;
888 Int_t nrow= mat.GetNrows();
0e9efcbe 889 TString strMask(mask);
2497be21 890 TVectorD vecCorrSum(nrow);
891 for (Int_t irow=0; irow<nrow; irow++){
892 vecCorrSum[irow]=0;
893 for (Int_t icol=0; icol<nrow; icol++){
894 Double_t diag = TMath::Sqrt(TMath::Abs(mat(irow,irow)*mat(icol,icol)));
895 Double_t corr0 = mat(irow,icol)/diag;
896 vecCorrSum[irow]+=TMath::Abs(corr0);
897 }
898 vecCorrSum[irow]-=1.;
899 vecCorrSum[irow]*=0.5;
900 }
901
0bfb86a0 902 for (Int_t irow=0; irow<nrow; irow++){
903 AliTPCTransformation * trans0 = GetTransformation(irow);
0e9efcbe 904 TString strName(trans0->GetName());
905 if (mask){
906 if (!strName.Contains(mask)) continue;
907 }
2497be21 908 if (vecCorrSum[irow]<correlationCut) continue;
909 printf("%d\t%s\t%f\t%f\t%f\n",
0bfb86a0 910 irow,
911 trans0->GetName(),
912 (*fCalibParam)(irow,0),
2497be21 913 TMath::Sqrt(mat(irow,irow)), vecCorrSum[irow]);
0bfb86a0 914 }
915 return (nrow>0);
916}
917
918
919Bool_t AliTPCkalmanFit::CheckCovariance(TMatrixD &mat, Float_t /*maxEl*/){
920 //
921 // Check consistency of covariance matrix
922 // + symetrize coavariance matrix
923 Bool_t isOK=kTRUE;
924 Int_t nrow= mat.GetNrows();
925 for (Int_t irow=0; irow<nrow; irow++){
926 if (mat(irow,irow)<=0){
927 printf("Negative covariance\t%d\t%f\n",irow,mat(irow,irow));
928 isOK=kFALSE;
929 }
930// if (mat(irow,irow)>maxEl){
931// printf("Too big covariance\t%d\t%f\n",irow,mat(irow,irow));
932// isOK=kFALSE;
933// }
934
935 for (Int_t icol=0; icol<nrow; icol++){
936 // if (mat(irow,irow)
937 Double_t diag = TMath::Sqrt(TMath::Abs(mat(irow,irow)*mat(icol,icol)));
938 if (diag<=0){
939 printf("Negative covariance\t%d\t%d\t%f\n",irow,icol, mat(irow,icol));
940 isOK=kFALSE;
941 continue;
942 }
943 Double_t cov0 = mat(irow,icol)/diag;
944 Double_t cov1 = mat(icol,irow)/diag;
945 if (TMath::Abs(cov0)>1 || TMath::Abs(cov1)>1 ){
946 printf("Covariance Problem %d\t%d\t%f\t%f\n",irow,icol, cov0, cov1);
947 isOK=kFALSE;
948 }
949 if (TMath::Abs(cov0-cov1)>0.0000001){
950 printf("Asymetry problem %d\t%d\t%f\t%f\n",irow,icol, cov0, cov1);
951 isOK=kFALSE;
952 }
953 //
954 // symetrize the covariance matrix
955 Double_t mean = (mat(irow,icol)+ mat(icol,irow))*0.5;
956 mat(irow,icol)=mean;
957 mat(icol,irow)=mean;
958 }
959 }
960 return isOK;
961}
962
963
964AliTPCkalmanFit * AliTPCkalmanFit::Test(Int_t ntracks){
965 //
966 // This is test example
967 //
968
969 //
970 // 1. Setup transformation
971 //
972
973
974 TVectorD fpar(10);
975 AliTPCTransformation * transformation=0;
976 AliTPCkalmanFit * kalmanFit0 = new AliTPCkalmanFit;
977 AliTPCkalmanFit * kalmanFit2 = new AliTPCkalmanFit;
978 //
979 // Radial scaling
980 //
981 for (Int_t iside=0; iside<=1; iside++)
982 for (Int_t ipar0=0; ipar0<3; ipar0++)
983 for (Int_t ipar1=0; ipar1<3; ipar1++){
984 fpar[0]=ipar0;
985 fpar[1]=ipar1;
986 if (ipar0+ipar1==0) continue;
987 Double_t param = (gRandom->Rndm()-0.5)*0.5; // generate random parameters
988 char tname[100];
989 sprintf(tname,"tscalingR%d%dSide%d",ipar0,ipar1,iside);
990 transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingRPol",0,0,1);
991 transformation->SetParams(0,5*0.25,0,&fpar);
992 kalmanFit0->AddCalibration(transformation);
993 transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingRPol",0,0,1);
994 transformation->SetParams(param,0.25,0,&fpar);
995 kalmanFit2->AddCalibration(transformation);
996 }
997
998
999 //
1000 // 2. Init transformation
1001 //
1002 kalmanFit2->Init();
1003 kalmanFit0->Init();
1004
1005 //
1006 // 3. Run kalman filter
1007 //
1008 Int_t ncalibs = kalmanFit0->fCalibration->GetEntries();
1009 TVectorD err(ncalibs);
1010 TTreeSRedirector *pcstream = new TTreeSRedirector("kalmanfitTest.root");
1011 for (Int_t i=0;i<ntracks;i++){
1012 if (i%100==0) printf("%d\n",i);
1013 Double_t alpha = gRandom->Rndm()*TMath::TwoPi();
1014 Double_t y0 = (gRandom->Rndm()-0.5)*180;
1015 Double_t z0 = (gRandom->Rndm()-0.5)*250*2;
1016 Double_t ky = (gRandom->Rndm()-0.5)*1;
1017 Double_t kz = (gRandom->Rndm()-0.5)*1;
1018 //generate random TPC track
1019 AliTrackPointArray * array = AliTPCkalmanFit::MakePointArrayLinear(alpha,y0,z0, ky, kz,0.04);
1020 AliTrackPointArray * arrayB = new AliTrackPointArray(*array); // backup
1021 kalmanFit2->ApplyCalibration(array,1.); // misalign ideal track
1022 for (Int_t icalib=0; icalib<ncalibs; icalib++){
1023 err[icalib] = TMath::Sqrt((*kalmanFit0->fCalibCovar)(icalib,icalib));
1024 }
1025 //
1026 (*pcstream)<<"dump0"<<
1027 "alpha="<<alpha<<
1028 "y0="<<y0<<
1029 "z0="<<z0<<
1030 "ky="<<ky<<
1031 "lz="<<kz<<
1032 "p.="<<array<<
1033 "pB.="<<arrayB<<
1034 "cparam.="<<kalmanFit0->fCalibParam<<
1035 "ccovar.="<<kalmanFit0->fCalibCovar<<
1036 "err.="<<&err<<
1037 "gparam.="<<kalmanFit2->fCalibParam<<
1038 "gcovar.="<<kalmanFit2->fCalibCovar<<
1039 "\n";
1040 if (i%20==0) {
c64c82db 1041 kalmanFit0->FitTrackLinear(*array,pcstream); // fit track - dump intermediate results
0bfb86a0 1042 }else{
c64c82db 1043 kalmanFit0->FitTrackLinear(*array,0); // fit track + calibration
0bfb86a0 1044 }
1045 kalmanFit0->DumpTrackLinear(*array,pcstream); // dump track residuals to the tree + fill histograms
1046 }
1047 pcstream->GetFile()->cd();
1048 kalmanFit0->Write("kalmanFit");
1049 delete pcstream;
1050 return kalmanFit0;
1051}
1052
1053
1054Double_t AliTPCkalmanFit::GetTPCDeltaXYZ(Int_t coord, Int_t volID, Double_t x, Double_t y, Double_t z){
1055 //
1056 // function for visualization purposes
1057 //
1058 if (!fCalibration) return 0;
1059 Int_t ncalibs = fCalibration->GetEntries();
1060 if (ncalibs==0) return 0;
1061 Double_t dxdydz[3]={0,0,0};
1062 //
1063 if (volID<0){
1064 Double_t alpha = TMath::ATan2(y,x);
1065 Double_t r = TMath::Sqrt(y*y+x*x);
1066 volID = TMath::Nint(9*alpha/TMath::Pi()-0.5);
1067 if (volID<0) volID+=18;
1068 if (z<0) volID+=18;
1069 if (r>120) volID+=36;
1070 }
1071 for (Int_t icalib=0; icalib<ncalibs; icalib++){
1072 AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
1073 Double_t param = (*fCalibParam)(icalib,0);
1074 dxdydz[coord] += transform->GetDeltaXYZ(coord,volID, param, x, y,z);
1075 }
1076 return dxdydz[coord];
1077}
1078
1079Double_t AliTPCkalmanFit::SGetTPCDeltaXYZ(Int_t coord, Int_t volID, Double_t x, Double_t y, Double_t z){
1080 //
1081 //
1082 //
1083 if (AliTPCkalmanFit::fgInstance==0) return 0;
1084 return AliTPCkalmanFit::fgInstance->GetTPCDeltaXYZ(coord, volID,x,y,z);
1085}
1086