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