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