]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCkalmanFit.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
* fixed bug while getting leading jet candidate
[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 //
c64c82db 78 for (Int_t ihis=0; ihis<12; ihis++){
0bfb86a0 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);
c64c82db 101 for (Int_t i=0;i<12;i++){
0e9efcbe 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 //
c64c82db 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++){
0bfb86a0 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
0e9efcbe 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
0bfb86a0 259
c64c82db 260void AliTPCkalmanFit::FitTrackLinear(AliTrackPointArray &points, TTreeSRedirector *debug){
0bfb86a0 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 //
c64c82db 299 for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
0bfb86a0 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;
c64c82db 342 for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
0bfb86a0 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
c64c82db 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 //
0bfb86a0 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;
c64c82db 684 //
0bfb86a0 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::MakePointArrayLinear(Double_t alpha, Double_t y0, Double_t z0, Double_t ky, Double_t kz, Double_t err){
717 //
718 //
719 //
720 AliTrackPointArray array(500);
721 Float_t cov[10]; // dummy covariance
722 Int_t npoints=0;
723 for (Int_t i=0;i<6;i++) cov[i]=0.001;
724 for (Int_t i=0;i<500;i++){
725 AliTrackPoint point(0, 0, 0, cov, -1,0,0);
726 array.AddPoint(npoints, &point);
727 npoints++;
728 }
729 npoints=0;
730 for (Float_t ir = -245; ir<245; ir++){
731 //
732 //
733 if (TMath::Abs(ir)<80) continue;
734 Double_t ca = TMath::Cos(alpha);
735 Double_t sa = TMath::Sin(alpha);
736 Double_t lx = ir;
737 Double_t ly = y0+lx*ky+gRandom->Gaus(0,err);
738 Double_t lz = z0+lx*kz+gRandom->Gaus(0,err);
739 Double_t gx = lx*ca-ly*sa;
740 Double_t gy = lx*sa+ly*ca;
741 Double_t gz = lz;
742 Double_t galpha = TMath::ATan2(gy,gx);
743 Int_t isec = TMath::Nint((9*galpha/TMath::Pi()+0.5));
744 if (isec<0) isec+=18;
745 if (gz<0) isec+=18;
746 if (ir>130) isec+=36;
747 //
748 AliTrackPoint point(gx, gy, gz, cov, isec,0,0);
749 array.AddPoint(npoints, &point);
750 npoints++;
751 }
752 AliTrackPointArray *parray = new AliTrackPointArray(npoints);
753 AliTrackPoint point;
754 for (Int_t i=0;i<npoints;i++){
755 array.GetPoint(point,i);
756 parray->AddPoint(i,&point);
757 }
758
759 return parray;
760}
761
762void AliTPCkalmanFit::AddCalibration(AliTPCTransformation * calib){
763 //
764 // Add calibration
765 //
766 if (!fCalibration) fCalibration = new TObjArray(2000);
767 fCalibration->AddLast(calib);
768}
769
770void AliTPCkalmanFit::ApplyCalibration(AliTrackPointArray *array, Double_t csign){
771 //
772 //
773 //
774 if (!fCalibration) return;
775 Int_t ncalibs = fCalibration->GetEntries();
776 if (ncalibs==0) return;
777 Int_t npoints = array->GetNPoints();
778 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
779 Int_t volId = array->GetVolumeID()[ipoint];
780 Double_t xyz[3]={array->GetX()[ipoint], array->GetY()[ipoint],array->GetZ()[ipoint]};
781 Double_t dxdydz[3]={0,0,0};
782 for (Int_t icalib=0; icalib<ncalibs; icalib++){
783 AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
784 dxdydz[0] += transform->GetDeltaXYZ(0,volId, transform->fParam, xyz[0], xyz[1], xyz[2]);
785 dxdydz[1] += transform->GetDeltaXYZ(1,volId, transform->fParam, xyz[0], xyz[1], xyz[2]);
786 dxdydz[2] += transform->GetDeltaXYZ(2,volId, transform->fParam, xyz[0], xyz[1], xyz[2]);
787 }
788 ((Float_t*)array->GetX())[ipoint]+=csign*dxdydz[0];
789 ((Float_t*)array->GetY())[ipoint]+=csign*dxdydz[1];
790 ((Float_t*)array->GetZ())[ipoint]+=csign*dxdydz[2];
791 }
792}
793
c64c82db 794Bool_t AliTPCkalmanFit::DumpCorelation(Double_t threshold, const char *mask0, const char *mask1){
0bfb86a0 795 //
796 //
797 //
798 TMatrixD &mat = *fCalibCovar;
799 Int_t nrow= mat.GetNrows();
800 for (Int_t irow=0; irow<nrow; irow++){
c64c82db 801 AliTPCTransformation * trans0 = GetTransformation(irow);
802 TString strName0(trans0->GetName());
803 if (mask0){
804 if (!strName0.Contains(mask0)) continue;
805 }
0bfb86a0 806 for (Int_t icol=irow+1; icol<nrow; icol++){
c64c82db 807 AliTPCTransformation * trans1 = GetTransformation(icol);
808 TString strName1(trans1->GetName());
809 if (mask1){
810 if (!strName1.Contains(mask1)) continue;
811 }
812 //
0bfb86a0 813 Double_t diag = TMath::Sqrt(TMath::Abs(mat(irow,irow)*mat(icol,icol)));
814 if (diag<=0){
815 printf("Negative covariance\t%d\t%d\t%f\n",irow,icol, mat(irow,icol));
816 continue;
817 }
818 Double_t corr0 = mat(irow,icol)/diag;
819 if (TMath::Abs(corr0)>threshold){
0bfb86a0 820 printf("%d\t%d\t%s\t%s\t%f\t%f\t%f\n", irow,icol, trans0->GetName(), trans1->GetName(),
821 TMath::Sqrt(mat(irow,irow)), TMath::Sqrt(mat(icol,icol)), corr0);
822 }
823 }
824 }
825 return (nrow>0);
826}
827
0e9efcbe 828Bool_t AliTPCkalmanFit::DumpCalib(const char *mask){
0bfb86a0 829 //
830 // Print calibration entries - name, value, error
831 //
832 TMatrixD &mat = *fCalibCovar;
833 Int_t nrow= mat.GetNrows();
0e9efcbe 834 TString strMask(mask);
0bfb86a0 835 for (Int_t irow=0; irow<nrow; irow++){
836 AliTPCTransformation * trans0 = GetTransformation(irow);
0e9efcbe 837 TString strName(trans0->GetName());
838 if (mask){
839 if (!strName.Contains(mask)) continue;
840 }
0bfb86a0 841 printf("%d\t%s\t%f\t%f\n",
842 irow,
843 trans0->GetName(),
844 (*fCalibParam)(irow,0),
845 TMath::Sqrt(mat(irow,irow)));
846 }
847 return (nrow>0);
848}
849
850
851Bool_t AliTPCkalmanFit::CheckCovariance(TMatrixD &mat, Float_t /*maxEl*/){
852 //
853 // Check consistency of covariance matrix
854 // + symetrize coavariance matrix
855 Bool_t isOK=kTRUE;
856 Int_t nrow= mat.GetNrows();
857 for (Int_t irow=0; irow<nrow; irow++){
858 if (mat(irow,irow)<=0){
859 printf("Negative covariance\t%d\t%f\n",irow,mat(irow,irow));
860 isOK=kFALSE;
861 }
862// if (mat(irow,irow)>maxEl){
863// printf("Too big covariance\t%d\t%f\n",irow,mat(irow,irow));
864// isOK=kFALSE;
865// }
866
867 for (Int_t icol=0; icol<nrow; icol++){
868 // if (mat(irow,irow)
869 Double_t diag = TMath::Sqrt(TMath::Abs(mat(irow,irow)*mat(icol,icol)));
870 if (diag<=0){
871 printf("Negative covariance\t%d\t%d\t%f\n",irow,icol, mat(irow,icol));
872 isOK=kFALSE;
873 continue;
874 }
875 Double_t cov0 = mat(irow,icol)/diag;
876 Double_t cov1 = mat(icol,irow)/diag;
877 if (TMath::Abs(cov0)>1 || TMath::Abs(cov1)>1 ){
878 printf("Covariance Problem %d\t%d\t%f\t%f\n",irow,icol, cov0, cov1);
879 isOK=kFALSE;
880 }
881 if (TMath::Abs(cov0-cov1)>0.0000001){
882 printf("Asymetry problem %d\t%d\t%f\t%f\n",irow,icol, cov0, cov1);
883 isOK=kFALSE;
884 }
885 //
886 // symetrize the covariance matrix
887 Double_t mean = (mat(irow,icol)+ mat(icol,irow))*0.5;
888 mat(irow,icol)=mean;
889 mat(icol,irow)=mean;
890 }
891 }
892 return isOK;
893}
894
895
896AliTPCkalmanFit * AliTPCkalmanFit::Test(Int_t ntracks){
897 //
898 // This is test example
899 //
900
901 //
902 // 1. Setup transformation
903 //
904
905
906 TVectorD fpar(10);
907 AliTPCTransformation * transformation=0;
908 AliTPCkalmanFit * kalmanFit0 = new AliTPCkalmanFit;
909 AliTPCkalmanFit * kalmanFit2 = new AliTPCkalmanFit;
910 //
911 // Radial scaling
912 //
913 for (Int_t iside=0; iside<=1; iside++)
914 for (Int_t ipar0=0; ipar0<3; ipar0++)
915 for (Int_t ipar1=0; ipar1<3; ipar1++){
916 fpar[0]=ipar0;
917 fpar[1]=ipar1;
918 if (ipar0+ipar1==0) continue;
919 Double_t param = (gRandom->Rndm()-0.5)*0.5; // generate random parameters
920 char tname[100];
921 sprintf(tname,"tscalingR%d%dSide%d",ipar0,ipar1,iside);
922 transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingRPol",0,0,1);
923 transformation->SetParams(0,5*0.25,0,&fpar);
924 kalmanFit0->AddCalibration(transformation);
925 transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingRPol",0,0,1);
926 transformation->SetParams(param,0.25,0,&fpar);
927 kalmanFit2->AddCalibration(transformation);
928 }
929
930
931 //
932 // 2. Init transformation
933 //
934 kalmanFit2->Init();
935 kalmanFit0->Init();
936
937 //
938 // 3. Run kalman filter
939 //
940 Int_t ncalibs = kalmanFit0->fCalibration->GetEntries();
941 TVectorD err(ncalibs);
942 TTreeSRedirector *pcstream = new TTreeSRedirector("kalmanfitTest.root");
943 for (Int_t i=0;i<ntracks;i++){
944 if (i%100==0) printf("%d\n",i);
945 Double_t alpha = gRandom->Rndm()*TMath::TwoPi();
946 Double_t y0 = (gRandom->Rndm()-0.5)*180;
947 Double_t z0 = (gRandom->Rndm()-0.5)*250*2;
948 Double_t ky = (gRandom->Rndm()-0.5)*1;
949 Double_t kz = (gRandom->Rndm()-0.5)*1;
950 //generate random TPC track
951 AliTrackPointArray * array = AliTPCkalmanFit::MakePointArrayLinear(alpha,y0,z0, ky, kz,0.04);
952 AliTrackPointArray * arrayB = new AliTrackPointArray(*array); // backup
953 kalmanFit2->ApplyCalibration(array,1.); // misalign ideal track
954 for (Int_t icalib=0; icalib<ncalibs; icalib++){
955 err[icalib] = TMath::Sqrt((*kalmanFit0->fCalibCovar)(icalib,icalib));
956 }
957 //
958 (*pcstream)<<"dump0"<<
959 "alpha="<<alpha<<
960 "y0="<<y0<<
961 "z0="<<z0<<
962 "ky="<<ky<<
963 "lz="<<kz<<
964 "p.="<<array<<
965 "pB.="<<arrayB<<
966 "cparam.="<<kalmanFit0->fCalibParam<<
967 "ccovar.="<<kalmanFit0->fCalibCovar<<
968 "err.="<<&err<<
969 "gparam.="<<kalmanFit2->fCalibParam<<
970 "gcovar.="<<kalmanFit2->fCalibCovar<<
971 "\n";
972 if (i%20==0) {
c64c82db 973 kalmanFit0->FitTrackLinear(*array,pcstream); // fit track - dump intermediate results
0bfb86a0 974 }else{
c64c82db 975 kalmanFit0->FitTrackLinear(*array,0); // fit track + calibration
0bfb86a0 976 }
977 kalmanFit0->DumpTrackLinear(*array,pcstream); // dump track residuals to the tree + fill histograms
978 }
979 pcstream->GetFile()->cd();
980 kalmanFit0->Write("kalmanFit");
981 delete pcstream;
982 return kalmanFit0;
983}
984
985
986Double_t AliTPCkalmanFit::GetTPCDeltaXYZ(Int_t coord, Int_t volID, Double_t x, Double_t y, Double_t z){
987 //
988 // function for visualization purposes
989 //
990 if (!fCalibration) return 0;
991 Int_t ncalibs = fCalibration->GetEntries();
992 if (ncalibs==0) return 0;
993 Double_t dxdydz[3]={0,0,0};
994 //
995 if (volID<0){
996 Double_t alpha = TMath::ATan2(y,x);
997 Double_t r = TMath::Sqrt(y*y+x*x);
998 volID = TMath::Nint(9*alpha/TMath::Pi()-0.5);
999 if (volID<0) volID+=18;
1000 if (z<0) volID+=18;
1001 if (r>120) volID+=36;
1002 }
1003 for (Int_t icalib=0; icalib<ncalibs; icalib++){
1004 AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
1005 Double_t param = (*fCalibParam)(icalib,0);
1006 dxdydz[coord] += transform->GetDeltaXYZ(coord,volID, param, x, y,z);
1007 }
1008 return dxdydz[coord];
1009}
1010
1011Double_t AliTPCkalmanFit::SGetTPCDeltaXYZ(Int_t coord, Int_t volID, Double_t x, Double_t y, Double_t z){
1012 //
1013 //
1014 //
1015 if (AliTPCkalmanFit::fgInstance==0) return 0;
1016 return AliTPCkalmanFit::fgInstance->GetTPCDeltaXYZ(coord, volID,x,y,z);
1017}
1018