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