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