AliRoot Core  v5-06-30 (35d6c57)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliTPCkalmanFit.cxx
Go to the documentation of this file.
1 
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 
66 
67 AliTPCkalmanFit* AliTPCkalmanFit::fgInstance = 0;
68 
70  TNamed(),
71  fCalibration(0),
72  fCalibParam(0),
73  fCalibCovar(0),
74  fLinearParam(0),
75  fLinearCovar(0),
76  fLastTimeStamp(-1),
77  fCurrentAlpha(0),
78  fCA(0),
79  fSA(0)
80 {
82 
83  for (Int_t ihis=0; ihis<12; ihis++){
84  fLinearTrackDelta[ihis]=0;
85  fLinearTrackPull[ihis]=0;
86  }
87 }
88 
91 
92  Int_t ncalibs = fCalibration->GetEntries();
93  for (Int_t icalib=0;icalib<ncalibs; icalib++){
94  AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
95  transform->Init();
96  }
97 }
98 
101 
102  Update(kalman);
103  for (Int_t i=0;i<12;i++){
104  if (fLinearTrackDelta[i] && kalman->fLinearTrackDelta[i]){
105  fLinearTrackDelta[i]->Add(kalman->fLinearTrackDelta[i]);
106  }
107  if (fLinearTrackPull[i] && kalman->fLinearTrackPull[i]){
108  fLinearTrackPull[i]->Add(kalman->fLinearTrackPull[i]);
109  }
110  }
111 
112 }
113 
114 
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->GetParam();
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->GetSigma()*transform->GetSigma();
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 
196 void AliTPCkalmanFit::SetStatus(const char * mask, Bool_t setOn, Bool_t isOr){
200 
201  Int_t ncalibs = fCalibration->GetEntries();
202  if (mask==0) {
203  for (Int_t i=0; i<ncalibs;i++){
205  transform->SetActive(setOn);
206  }
207  }
208  else{
209  for (Int_t i=0; i<ncalibs;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 
223 
224  Int_t ncalibs = fCalibration->GetEntries();
225  TMatrixD vecXk=*fCalibParam;
226  TMatrixD covXk=*fCalibCovar;
227  TMatrixD &vecZk = *(kalman->fCalibParam);
228  TMatrixD &measR = *(kalman->fCalibCovar);
229 
230  TMatrixD matHk(ncalibs,ncalibs);
231  TMatrixD vecYk(ncalibs,1);
232  TMatrixD matHkT(ncalibs,ncalibs);
233  TMatrixD matSk(ncalibs,ncalibs);
234  TMatrixD matKk(ncalibs,ncalibs);
235  TMatrixD covXk2(ncalibs,ncalibs);
236  TMatrixD covXk3(ncalibs,ncalibs);
237  //
238  for (Int_t i=0;i<ncalibs;i++){
239  for (Int_t j=0;j<ncalibs;j++) matHk(i,j)=0;
240  matHk(i,i)=1;
241  }
242  vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
243  matHkT=matHk.T(); matHk.T();
244  matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
245  matSk.Invert();
246  matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
247  vecXk += matKk*vecYk; // updated vector
248  covXk2= (matHk-(matKk*matHk));
249  covXk3 = covXk2*covXk;
250  covXk = covXk3;
251  (*fCalibParam) = vecXk;
252  (*fCalibCovar) = covXk;
253 }
254 
255 
256 
257 
258 
259 void AliTPCkalmanFit::FitTrackLinear(AliTrackPointArray &points, TTreeSRedirector *debug, Float_t scalingRMSY, Float_t scalingRMSZ){
261 
262  if (!fCalibParam) {
263  AliError("Kalman Fit not initialized");
264  return;
265  }
266  //
267  // 1. Make initial track hypothesis
268  //
269  TLinearFitter lfitY(2,"pol1");
270  TLinearFitter lfitZ(2,"pol1");
271  TVectorD vecZ(2);
272  TVectorD vecY(2);
273  //
274  lfitY.StoreData(kTRUE);
275  lfitZ.StoreData(kTRUE);
276  lfitY.ClearPoints();
277  lfitZ.ClearPoints();
278  Int_t npoints = points.GetNPoints();
279  if (npoints<2) return;
280  const Double_t kFac=npoints*npoints*100;
281  const Double_t kOff0=0.01*0.01;
282  const Double_t kOff1=kOff0/(250.*250.);
283  //
284  // 0. Make seed
285  //
286  // AliTrackPointArray pointsc(points);
287  //ApplyCalibration(&pointsc, -1.);
288  //
289  // 1.a Choosing reference rotation alpha
290  //
291  fCurrentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]);
292  fCA = TMath::Cos(fCurrentAlpha);
293  fSA = TMath::Sin(fCurrentAlpha);
294  //
295  // 1.b Fit the track in the rotated frame - MakeSeed
296  //
297  for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
298  Double_t rx = fCA*points.GetX()[ipoint]+fSA*points.GetY()[ipoint];
299  Double_t ry = -fSA*points.GetX()[ipoint]+fCA*points.GetY()[ipoint];
300  Double_t rz = points.GetZ()[ipoint];
301  lfitY.AddPoint(&rx,ry,1);
302  lfitZ.AddPoint(&rx,rz,1);
303  }
304  lfitY.Eval();
305  lfitZ.Eval();
306  lfitY.GetParameters(vecY);
307  lfitZ.GetParameters(vecZ);
308  //
309  // 2. Set initial parameters and the covariance matrix
310  //
311  Int_t ncalibs = fCalibration->GetEntries();
312  if (!fLinearParam) {
313  fLinearParam = new TMatrixD(ncalibs+4,1);
314  fLinearCovar = new TMatrixD(ncalibs+4,ncalibs+4);
315  }
316  //
317  for (Int_t i=0; i<ncalibs+4;i++){
318  (*fLinearParam)(i,0)=0;
319  if (i<ncalibs) (*fLinearParam)(i,0) = (*fCalibParam)(i,0);
320  for (Int_t j=0; j<ncalibs+4;j++){
321  (*fLinearCovar)(i,j)=0;
322  if (i<ncalibs&&j<ncalibs) (*fLinearCovar)(i,j) = (*fCalibCovar)(i,j);
323  }
324  }
325  Double_t chi2Y = lfitY.GetChisquare()/lfitY.GetNpoints();
326  Double_t chi2Z = lfitZ.GetChisquare()/lfitZ.GetNpoints();
327  (*fLinearParam)(ncalibs+0,0) = lfitY.GetParameter(0);
328  (*fLinearCovar)(ncalibs+0,ncalibs+0)= lfitY.GetCovarianceMatrixElement(0,0)*chi2Y*kFac+kOff0;
329  (*fLinearParam)(ncalibs+1,0) = lfitY.GetParameter(1);
330  (*fLinearCovar)(ncalibs+1,ncalibs+1)= lfitY.GetCovarianceMatrixElement(1,1)*chi2Y*kFac+kOff1;
331  //
332  (*fLinearParam)(ncalibs+2,0) = lfitZ.GetParameter(0);
333  (*fLinearCovar)(ncalibs+2,ncalibs+2)= lfitZ.GetCovarianceMatrixElement(0,0)*chi2Z*kFac+kOff0;
334  (*fLinearParam)(ncalibs+3,0) = lfitZ.GetParameter(1);
335  (*fLinearCovar)(ncalibs+3,ncalibs+3)= lfitZ.GetCovarianceMatrixElement(1,1)*chi2Z*kFac+kOff1;
336  //
337  // Fit thetrack together with correction
338  //
339  AliTrackPoint point;
340  for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
341  //
342  if (!points.GetPoint(point,ipoint)) continue;
343  Double_t erry2 = chi2Y;
344  Double_t errz2 = chi2Z;
345  //set error - it is hack
346  Float_t *cov = (Float_t*) point.GetCov();
347  cov[1] = (erry2+kOff0)*scalingRMSY;
348  cov[2] = (errz2+kOff0)*scalingRMSZ;
349  UpdateLinear(point,debug);
350  if (!points.GetPoint(point,npoints-1-ipoint)) continue;
351  //set error - it is hack
352  cov = (Float_t*) point.GetCov();
353  cov[1] = (erry2+kOff0)*scalingRMSY;
354  cov[2] = (errz2+kOff0)*scalingRMSZ;
355  UpdateLinear(point,debug);
356  }
357  //
358  // save current param and covariance
359  for (Int_t i=0; i<ncalibs;i++){
361  transform->SetParam( (*fLinearParam)(i,0));
362  (*fCalibParam)(i,0) = (*fLinearParam)(i,0);
363  for (Int_t j=0; j<ncalibs;j++){
364  (*fCalibCovar)(i,j) = (*fLinearCovar)(i,j);
365  }
366  }
367  if (debug){ // dump debug info
368  (*debug)<<"fitLinear"<<
369  "vecY.="<<&vecY<<
370  "vecZ.="<<&vecZ<<
371  "chi2Y="<<chi2Y<<
372  "chi2Z="<<chi2Z<<
373  "lP.="<<fLinearParam<<
374  "cP.="<<fCalibParam<<
375  "lC.="<<fLinearCovar<<
376  "cC.="<<fCalibCovar<<
377  "\n";
378  }
379 }
380 
381 void AliTPCkalmanFit::DumpTrackLinear(AliTrackPointArray &points, TTreeSRedirector *debug){
386 
387  if (!fCalibParam) {
388  AliError("Kalman Fit not initialized");
389  return;
390  }
391  //
392  //
393  //
394  TLinearFitter *fitters[16];
395  TVectorD *params[16];
396  TVectorD *errs[16];
397  TVectorD chi2N(16);
398  Int_t npoints = points.GetNPoints();
399  AliTrackPointArray pointsTrans(points);
400  ApplyCalibration(&pointsTrans,-1.);
401  for (Int_t ifit=0; ifit<8;ifit++){
402  fitters[ifit] = new TLinearFitter(2,"pol1");
403  params[ifit] = new TVectorD(2);
404  fitters[ifit+8]= new TLinearFitter(3,"pol2");
405  params[ifit+8] = new TVectorD(3);
406  errs[ifit] = new TVectorD(2);
407  errs[ifit+8] = new TVectorD(3);
408  }
409  //
410  // calculate mean x point and corrdinate frame
411  //
412  fCurrentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]);
413  fCA = TMath::Cos(fCurrentAlpha);
414  fSA = TMath::Sin(fCurrentAlpha);
415  Double_t xmean=0, sum=0;
416  for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
417  Double_t rx = fCA*points.GetX()[ipoint]+fSA*points.GetY()[ipoint];
418  xmean+=rx;
419  sum++;
420  }
421  xmean/=sum;
422  //
423  for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
424  Double_t rx = fCA*points.GetX()[ipoint]+fSA*points.GetY()[ipoint];
425  Double_t ry = -fSA*points.GetX()[ipoint]+fCA*points.GetY()[ipoint];
426  Double_t rz = points.GetZ()[ipoint];
427  //
428  Double_t rxT = fCA*pointsTrans.GetX()[ipoint]+fSA*pointsTrans.GetY()[ipoint];
429  Double_t ryT = -fSA*pointsTrans.GetX()[ipoint]+fCA*pointsTrans.GetY()[ipoint];
430  Double_t rzT = pointsTrans.GetZ()[ipoint];
431  if (rx>xmean){
432  fitters[0]->AddPoint(&rxT, ryT,1);
433  fitters[2]->AddPoint(&rxT, rzT,1);
434  fitters[4]->AddPoint(&rx, ry,1);
435  fitters[6]->AddPoint(&rx, rz,1);
436  fitters[8]->AddPoint(&rxT, ryT,1);
437  fitters[10]->AddPoint(&rxT, rzT,1);
438  fitters[12]->AddPoint(&rx, ry,1);
439  fitters[14]->AddPoint(&rx, rz,1);
440  }else{
441  fitters[1]->AddPoint(&rxT, ryT,1);
442  fitters[3]->AddPoint(&rxT, rzT,1);
443  fitters[5]->AddPoint(&rx, ry,1);
444  fitters[7]->AddPoint(&rx, rz,1);
445  fitters[9]->AddPoint(&rxT, ryT,1);
446  fitters[11]->AddPoint(&rxT, rzT,1);
447  fitters[13]->AddPoint(&rx, ry,1);
448  fitters[15]->AddPoint(&rx, rz,1);
449  }
450  }
451  for (Int_t ifit=0;ifit<16;ifit++){
452  fitters[ifit]->Eval();
453  fitters[ifit]->GetParameters(*(params[ifit]));
454  fitters[ifit]->GetErrors(*(errs[ifit]));
455  chi2N[ifit] = TMath::Sqrt(fitters[ifit]->GetChisquare()/(fitters[ifit]->GetNpoints()-1));
456  (*(errs[ifit]))[0]*=chi2N[ifit];
457  (*(errs[ifit]))[1]*=chi2N[ifit];
458  if (ifit>=8) (*(errs[ifit]))[2]*=chi2N[ifit]; // second derivative
459  }
460  if (debug){
461  (*debug)<<"dumpLinear"<<
462  "alpha="<<fCurrentAlpha<<
463  "xmean="<<xmean<<
464  "y0T.="<<params[0]<<
465  "y1T.="<<params[1]<<
466  "z0T.="<<params[2]<<
467  "z1T.="<<params[3]<<
468  "y0O.="<<params[4]<<
469  "y1O.="<<params[5]<<
470  "z0O.="<<params[6]<<
471  "z1O.="<<params[7]<<
472  "y0T2.="<<params[8]<<
473  "y1T2.="<<params[9]<<
474  "z0T2.="<<params[10]<<
475  "z1T2.="<<params[11]<<
476  "y0O2.="<<params[12]<<
477  "y1O2.="<<params[13]<<
478  "z0O2.="<<params[14]<<
479  "z1O2.="<<params[15]<<
480  "y0TErr.="<<errs[0]<<
481  "y1TErr.="<<errs[1]<<
482  "z0TErr.="<<errs[2]<<
483  "z1TErr.="<<errs[3]<<
484  "y0OErr.="<<errs[4]<<
485  "y1OErr.="<<errs[5]<<
486  "z0OErr.="<<errs[6]<<
487  "z1OErr.="<<errs[7]<<
488  "y0T2Err.="<<errs[8]<<
489  "y1T2Err.="<<errs[9]<<
490  "z0T2Err.="<<errs[10]<<
491  "z1T2Err.="<<errs[11]<<
492  "y0O2Err.="<<errs[12]<<
493  "y1O2Err.="<<errs[13]<<
494  "z0O2Err.="<<errs[14]<<
495  "z1O2Err.="<<errs[15]<<
496  "chi2.="<<&chi2N<<
497  "\n";
498  }
499  //
500  // delta alpha y0 z0 ky kz
501  Double_t x[6]={0,0,0,0,0,0};
502  x[1]=fCurrentAlpha;
503  x[2]=(*params[0])[0];
504  x[3]=(*params[2])[0];
505  x[4]=(*params[0])[1];
506  x[5]=(*params[2])[1];
507  //
508  //Delta y
509  //
510  x[0]= (*params[1])[0]-(*params[0])[0];
511  fLinearTrackDelta[0]->Fill(x);
512  x[0]/=TMath::Sqrt((*errs[1])[0]*(*errs[1])[0]+(*errs[0])[0]*(*errs[0])[0]);
513  fLinearTrackPull[0]->Fill(x);
514  //
515  //Delta ky
516  //
517  x[0]= 100*((*params[1])[1]-(*params[0])[1]);
518  fLinearTrackDelta[1]->Fill(x);
519  x[0]/=100*TMath::Sqrt((*errs[1])[1]*(*errs[1])[1]+(*errs[0])[1]*(*errs[0])[1]);
520  fLinearTrackPull[1]->Fill(x);
521  //
522  // ky2_0
523  //
524  x[0]= 100*100*((*params[8])[2]);
525  fLinearTrackDelta[2]->Fill(x);
526  x[0]/=100*100*TMath::Sqrt((*errs[8])[2]*(*errs[8])[2]);
527  fLinearTrackPull[2]->Fill(x);
528  //
529  // ky2_1
530  //
531  x[0]= 100*100*((*params[9])[2]);
532  fLinearTrackDelta[3]->Fill(x);
533  x[0]/=100*100*TMath::Sqrt((*errs[9])[2]*(*errs[9])[2]);
534  fLinearTrackPull[3]->Fill(x);
535  //
536  //
537  //Delta z
538  //
539  x[0]= (*params[3])[0]-(*params[2])[0];
540  fLinearTrackDelta[4]->Fill(x);
541  x[0]/=TMath::Sqrt((*errs[3])[0]*(*errs[3])[0]+(*errs[2])[0]*(*errs[2])[0]);
542  fLinearTrackPull[4]->Fill(x);
543  //
544  //Delta kz
545  //
546  x[0]= 100*((*params[3])[1]-(*params[2])[1]);
547  fLinearTrackDelta[5]->Fill(x);
548  x[0]/=100*TMath::Sqrt((*errs[3])[1]*(*errs[3])[1]+(*errs[2])[1]*(*errs[2])[1]);
549  fLinearTrackPull[5]->Fill(x);
550  //
551  // kz2_0
552  //
553  x[0]= 100*100*((*params[10])[2]);
554  fLinearTrackDelta[6]->Fill(x);
555  x[0]/=100*100*TMath::Sqrt((*errs[10])[2]*(*errs[10])[2]);
556  fLinearTrackPull[6]->Fill(x);
557  //
558  // kz2_1
559  //
560  x[0]= 100*100*((*params[11])[2]);
561  fLinearTrackDelta[7]->Fill(x);
562  x[0]/=100*100*TMath::Sqrt((*errs[11])[2]*(*errs[11])[2]);
563  fLinearTrackPull[7]->Fill(x);
564 
565  //
566  // rms of track
567  //
568  x[0]= chi2N[0];
569  fLinearTrackDelta[8]->Fill(x);
570  x[0]= chi2N[1];
571  fLinearTrackDelta[9]->Fill(x);
572  x[0]= chi2N[2];
573  fLinearTrackDelta[10]->Fill(x);
574  x[0]= chi2N[3];
575  fLinearTrackDelta[11]->Fill(x);
576  //
577 
578 
579  for (Int_t ifit=0; ifit<8;ifit++){
580  delete fitters[ifit];
581  delete params[ifit];
582  delete fitters[ifit+8];
583  delete params[ifit+8];
584  delete errs[ifit];
585  delete errs[ifit+8];
586  }
587 }
588 
589 
590 void AliTPCkalmanFit::Propagate(TTreeSRedirector * /*debug*/){
592 
593 }
594 
595 void AliTPCkalmanFit::AddCovariance(const char * varName, Double_t sigma){
597 
598  if (!fCalibCovar) return;
599  if (!fCalibration) return;
600  if (!fCalibration->FindObject(varName)) return;
601  Int_t ncalibs = fCalibration->GetEntries();
602  TString strVar(varName);
603  for (Int_t icalib=0;icalib<ncalibs; icalib++){
604  AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
605  if (strVar.CompareTo(transform->GetName())==0){
606  (*fCalibCovar)(icalib,icalib)+=sigma*sigma;
607  }
608  }
609 }
610 
611 
615 
616  if (!fCalibCovar) return;
617  Int_t ncalibs = fCalibration->GetEntries();
618  Double_t deltaT = 0;
619  if (fLastTimeStamp>0) deltaT = (fLastTimeStamp-time)/(60*60.); // delta T in hours
620  fLastTimeStamp = time;
621  for (Int_t icalib=0;icalib<ncalibs; icalib++){
622  AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
623  if ((*fCalibCovar)(icalib,icalib)<transform->GetSigmaMax()*transform->GetSigmaMax())
624  (*fCalibCovar)(icalib,icalib)+= transform->GetSigma2Time()*TMath::Abs(deltaT);
625  }
626 }
627 
628 
629 void AliTPCkalmanFit::UpdateLinear(AliTrackPoint &point, TTreeSRedirector *debug){
631 
632  Int_t ncalibs = fCalibration->GetEntries();
633  Int_t kNmeas = 2;
634  Int_t nelem = ncalibs+4;
635  TMatrixD &vecXk=*fLinearParam;
636  TMatrixD &covXk=*fLinearCovar;
637  //
638  TMatrixD matHk(kNmeas,nelem);
639  TMatrixD vecYk(kNmeas,1);
640  TMatrixD vecZk(kNmeas,1);
641  TMatrixD measR(kNmeas,kNmeas);
642  TMatrixD matHkT(nelem,kNmeas);
643  TMatrixD matSk(kNmeas,kNmeas);
644  TMatrixD matKk(nelem,kNmeas);
645  TMatrixD covXk2(nelem,nelem);
646  TMatrixD covXk3(nelem,nelem);
647  TMatrixD mat1(nelem,nelem);
648  //
649  // reset matHk
650  for (Int_t iel=0;iel<nelem;iel++){
651  for (Int_t ip=0;ip<kNmeas;ip++) {
652  matHk(ip,iel)=0;
653  }
654  }
655  for (Int_t iel0=0;iel0<nelem;iel0++)
656  for (Int_t iel1=0;iel1<nelem;iel1++){
657  if (iel0!=iel1) mat1(iel0,iel1)=0;
658  if (iel0==iel1) mat1(iel0,iel1)=1;
659  }
660  //
661  // fill matrix Hk
662  //
663  Int_t volId = point.GetVolumeID();
664  Double_t gxyz[3]={point.GetX(), point.GetY(),point.GetZ()};
665  Double_t rxyz[3]={fCA*gxyz[0]+fSA*gxyz[1],-fSA*gxyz[0]+fCA*gxyz[1] ,point.GetZ()};
666  //
667  Double_t dxdydz[3]={0,0,0};
668  Double_t rdxdydz[3]={0,0,0};
669 
670  for (Int_t icalib=0;icalib<ncalibs;icalib++){
671  AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
672  dxdydz[0] = transform->GetDeltaXYZ(0,volId, 1, gxyz[0], gxyz[1], gxyz[2]);
673  dxdydz[1] = transform->GetDeltaXYZ(1,volId, 1, gxyz[0], gxyz[1], gxyz[2]);
674  dxdydz[2] = transform->GetDeltaXYZ(2,volId, 1, gxyz[0], gxyz[1], gxyz[2]);
675  rdxdydz[0] = fCA*dxdydz[0]+fSA*dxdydz[1];
676  rdxdydz[1] = -fSA*dxdydz[0]+fCA*dxdydz[1];
677  rdxdydz[2] = dxdydz[2];
678  //
679  matHk(0,icalib)= rdxdydz[1]-rdxdydz[0]*(*fLinearParam)(ncalibs+1,0); // shift y + shift x * angleY
680  matHk(1,icalib)= rdxdydz[2]-rdxdydz[0]*(*fLinearParam)(ncalibs+3,0); // shift z + shift x * angleZ
681  }
682  matHk(0,ncalibs+0)=1;
683  matHk(0,ncalibs+1)=rxyz[0];
684  matHk(1,ncalibs+2)=1;
685  matHk(1,ncalibs+3)=rxyz[0];
686  //
687  //
688  //
689  vecZk(0,0) = rxyz[1];
690  vecZk(1,0) = rxyz[2];
691  measR(0,0) = point.GetCov()[1]; measR(0,1)=0;
692  measR(1,1) = point.GetCov()[2]; measR(1,0)=0;
693  //
694  vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
695  matHkT=matHk.T(); matHk.T();
696  matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
697  matSk.Invert();
698  matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
699  //
700  covXk2= (mat1-(matKk*matHk));
701  covXk3 = covXk2*covXk;
702  //
703  CheckCovariance(covXk3,100);
704  vecXk += matKk*vecYk; // updated vector
705  covXk = covXk3;
706  if (debug){ // dump debug info
707  (*debug)<<"updateLinear"<<
708  //
709  "point.="<<&point<<
710  "vecYk.="<<&vecYk<<
711  "vecZk.="<<&vecZk<<
712  "measR.="<<&measR<<
713  "matHk.="<<&matHk<<
714  "matHkT.="<<&matHkT<<
715  "matSk.="<<&matSk<<
716  "matKk.="<<&matKk<<
717  "covXk2.="<<&covXk2<<
718  "covXk.="<<&covXk<<
719  "vecXk.="<<&vecXk<<
720  "\n";
721  }
722 }
723 
724 
725 AliTrackPointArray * AliTPCkalmanFit::SortPoints(AliTrackPointArray &points){
729 
730  Int_t npoints = points.GetNPoints();
731  if (npoints<1) return 0;
732  Double_t currentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]);
733  Double_t ca = TMath::Cos(currentAlpha);
734  Double_t sa = TMath::Sin(currentAlpha);
735  //
736  // 1. sort the points
737  //
738  Double_t *rxvector = new Double_t[npoints];
739  Int_t *indexes = new Int_t[npoints];
740  for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
741  rxvector[ipoint]=ca*points.GetX()[ipoint]+sa*points.GetY()[ipoint];
742  }
743  TMath::Sort(npoints, rxvector,indexes,kFALSE);
744  AliTrackPoint point;
745  AliTrackPointArray *pointsSorted= new AliTrackPointArray(npoints);
746  for (Int_t ipoint=0; ipoint<npoints; ipoint++){
747  if (!points.GetPoint(point,indexes[ipoint])) continue;
748  pointsSorted->AddPoint(ipoint,&point);
749  }
750  delete [] rxvector;
751  delete [] indexes;
752  return pointsSorted;
753 }
754 
755 
756 
757 AliTrackPointArray * AliTPCkalmanFit::MakePointArrayLinear(Double_t alpha, Double_t y0, Double_t z0, Double_t ky, Double_t kz, Double_t err){
759 
760  AliTrackPointArray array(500);
761  Float_t cov[10];
762  Int_t npoints=0;
763  for (Int_t i=0;i<6;i++) cov[i]=0.001;
764  for (Int_t i=0;i<500;i++){
765  AliTrackPoint point(0, 0, 0, cov, 0,0,0);
766  array.AddPoint(npoints, &point);
767  npoints++;
768  }
769  npoints=0;
770  for (Float_t ir = -245; ir<245; ir++){
771  //
772  //
773  if (TMath::Abs(ir)<80) continue;
774  Double_t ca = TMath::Cos(alpha);
775  Double_t sa = TMath::Sin(alpha);
776  Double_t lx = ir;
777  Double_t ly = y0+lx*ky+gRandom->Gaus(0,err);
778  Double_t lz = z0+lx*kz+gRandom->Gaus(0,err);
779  Double_t gx = lx*ca-ly*sa;
780  Double_t gy = lx*sa+ly*ca;
781  Double_t gz = lz;
782  Double_t galpha = TMath::ATan2(gy,gx);
783  Int_t isec = TMath::Nint((9*galpha/TMath::Pi()+0.5));
784  if (isec<0) isec+=18;
785  if (gz<0) isec+=18;
786  if (ir>130) isec+=36;
787  //
788  AliTrackPoint point(gx, gy, gz, cov, isec,0,0);
789  array.AddPoint(npoints, &point);
790  npoints++;
791  }
792  AliTrackPointArray *parray = new AliTrackPointArray(npoints);
793  AliTrackPoint point;
794  for (Int_t i=0;i<npoints;i++){
795  array.GetPoint(point,i);
796  parray->AddPoint(i,&point);
797  }
798 
799  return parray;
800 }
801 
804 
805  if (!fCalibration) fCalibration = new TObjArray(2000);
806  fCalibration->AddLast(calib);
807 }
808 
809 Int_t AliTPCkalmanFit::GetTransformationIndex(const char * trName){
811 
812  if (!fCalibration) return -1;
813  if (!fCalibration->FindObject(trName)) return -1;
814  //
815  Int_t ncalibs = fCalibration->GetEntries();
816  TString strVar(trName);
817  for (Int_t icalib=0;icalib<ncalibs; icalib++){
818  AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
819  if (strVar.CompareTo(transform->GetName())==0){
820  return icalib;
821  }
822  }
823  return -1;
824 }
825 
826 
827 
828 void AliTPCkalmanFit::ApplyCalibration(AliTrackPointArray *array, Double_t csign){
830 
831  if (!fCalibration) return;
832  Int_t ncalibs = fCalibration->GetEntries();
833  if (ncalibs==0) return;
834  Int_t npoints = array->GetNPoints();
835  for (Int_t ipoint=0; ipoint<npoints; ipoint++){
836  Int_t volId = array->GetVolumeID()[ipoint];
837  Double_t xyz[3]={array->GetX()[ipoint], array->GetY()[ipoint],array->GetZ()[ipoint]};
838  Double_t dxdydz[3]={0,0,0};
839  for (Int_t icalib=0; icalib<ncalibs; icalib++){
840  AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
841  dxdydz[0] += transform->GetDeltaXYZ(0,volId, transform->GetParam(), xyz[0], xyz[1], xyz[2]);
842  dxdydz[1] += transform->GetDeltaXYZ(1,volId, transform->GetParam(), xyz[0], xyz[1], xyz[2]);
843  dxdydz[2] += transform->GetDeltaXYZ(2,volId, transform->GetParam(), xyz[0], xyz[1], xyz[2]);
844  }
845  ((Float_t*)array->GetX())[ipoint]+=csign*dxdydz[0];
846  ((Float_t*)array->GetY())[ipoint]+=csign*dxdydz[1];
847  ((Float_t*)array->GetZ())[ipoint]+=csign*dxdydz[2];
848  }
849 }
850 
851 Bool_t AliTPCkalmanFit::DumpCorelation(Double_t threshold, const char *mask0, const char *mask1){
853 
854  TMatrixD &mat = *fCalibCovar;
855  Int_t nrow= mat.GetNrows();
856  for (Int_t irow=0; irow<nrow; irow++){
857  AliTPCTransformation * trans0 = GetTransformation(irow);
858  TString strName0(trans0->GetName());
859  if (mask0){
860  if (!strName0.Contains(mask0)) continue;
861  }
862  for (Int_t icol=irow+1; icol<nrow; icol++){
863  AliTPCTransformation * trans1 = GetTransformation(icol);
864  TString strName1(trans1->GetName());
865  if (mask1){
866  if (!strName1.Contains(mask1)) continue;
867  }
868  //
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  continue;
873  }
874  Double_t corr0 = mat(irow,icol)/diag;
875  if (TMath::Abs(corr0)>threshold){
876  printf("%d\t%d\t%s\t%s\t%f\t%f\t%f\n", irow,icol, trans0->GetName(), trans1->GetName(),
877  TMath::Sqrt(mat(irow,irow)), TMath::Sqrt(mat(icol,icol)), corr0);
878  }
879  }
880  }
881  return (nrow>0);
882 }
883 
884 Bool_t AliTPCkalmanFit::DumpCalib(const char *mask, Float_t correlationCut){
886 
887  TMatrixD &mat = *fCalibCovar;
888  Int_t nrow= mat.GetNrows();
889  TString strMask(mask);
890  TVectorD vecCorrSum(nrow);
891  for (Int_t irow=0; irow<nrow; irow++){
892  vecCorrSum[irow]=0;
893  for (Int_t icol=0; icol<nrow; icol++){
894  if (icol==irow) continue;
895  Double_t diag = TMath::Sqrt(TMath::Abs(mat(irow,irow)*mat(icol,icol)));
896  Double_t corr0 = mat(irow,icol)/diag;
897  vecCorrSum[irow]+=TMath::Abs(corr0);
898  }
899  vecCorrSum[irow]*=0.5;
900  }
901 
902  for (Int_t irow=0; irow<nrow; irow++){
903  AliTPCTransformation * trans0 = GetTransformation(irow);
904  TString strName(trans0->GetName());
905  if (mask){
906  if (!strName.Contains(mask)) continue;
907  }
908  if (vecCorrSum[irow]<correlationCut) continue;
909  printf("%d\t%s\t%f\t%f\t%f\n",
910  irow,
911  trans0->GetName(),
912  (*fCalibParam)(irow,0),
913  TMath::Sqrt(mat(irow,irow)), vecCorrSum[irow]);
914  }
915  return (nrow>0);
916 }
917 
918 
919 Bool_t AliTPCkalmanFit::CheckCovariance(TMatrixD &mat, Float_t /*maxEl*/){
922 
923  Bool_t isOK=kTRUE;
924  Int_t nrow= mat.GetNrows();
925  for (Int_t irow=0; irow<nrow; irow++){
926  if (mat(irow,irow)<=0){
927  printf("Negative covariance\t%d\t%f\n",irow,mat(irow,irow));
928  isOK=kFALSE;
929  }
930 // if (mat(irow,irow)>maxEl){
931 // printf("Too big covariance\t%d\t%f\n",irow,mat(irow,irow));
932 // isOK=kFALSE;
933 // }
934 
935  for (Int_t icol=0; icol<nrow; icol++){
936  // if (mat(irow,irow)
937  Double_t diag = TMath::Sqrt(TMath::Abs(mat(irow,irow)*mat(icol,icol)));
938  if (diag<=0){
939  printf("Negative covariance\t%d\t%d\t%f\n",irow,icol, mat(irow,icol));
940  isOK=kFALSE;
941  continue;
942  }
943  Double_t cov0 = mat(irow,icol)/diag;
944  Double_t cov1 = mat(icol,irow)/diag;
945  if (TMath::Abs(cov0)>1 || TMath::Abs(cov1)>1 ){
946  printf("Covariance Problem %d\t%d\t%f\t%f\n",irow,icol, cov0, cov1);
947  isOK=kFALSE;
948  }
949  if (TMath::Abs(cov0-cov1)>0.0000001){
950  printf("Asymetry problem %d\t%d\t%f\t%f\n",irow,icol, cov0, cov1);
951  isOK=kFALSE;
952  }
953  //
954  // symetrize the covariance matrix
955  Double_t mean = (mat(irow,icol)+ mat(icol,irow))*0.5;
956  mat(irow,icol)=mean;
957  mat(icol,irow)=mean;
958  }
959  }
960  return isOK;
961 }
962 
963 
966 
967  //
968  // 1. Setup transformation
969  //
970 
971 
972  TVectorD fpar(10);
973  AliTPCTransformation * transformation=0;
974  AliTPCkalmanFit * kalmanFit0 = new AliTPCkalmanFit;
975  AliTPCkalmanFit * kalmanFit2 = new AliTPCkalmanFit;
976  //
977  // Radial scaling
978  //
979  for (Int_t iside=0; iside<=1; iside++)
980  for (Int_t ipar0=0; ipar0<3; ipar0++)
981  for (Int_t ipar1=0; ipar1<3; ipar1++){
982  fpar[0]=ipar0;
983  fpar[1]=ipar1;
984  if (ipar0+ipar1==0) continue;
985  Double_t param = (gRandom->Rndm()-0.5)*0.5;
986  char tname[100];
987  snprintf(tname,100,"tscalingR%d%dSide%d",ipar0,ipar1,iside);
988  transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingRPol",0,0,1);
989  transformation->SetParams(0,5*0.25,0,&fpar);
990  kalmanFit0->AddCalibration(transformation);
991  transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingRPol",0,0,1);
992  transformation->SetParams(param,0.25,0,&fpar);
993  kalmanFit2->AddCalibration(transformation);
994  }
995 
996 
997  //
998  // 2. Init transformation
999  //
1000  kalmanFit2->Init();
1001  kalmanFit0->Init();
1002 
1003  //
1004  // 3. Run kalman filter
1005  //
1006  Int_t ncalibs = kalmanFit0->fCalibration->GetEntries();
1007  TVectorD err(ncalibs);
1008  TTreeSRedirector *pcstream = new TTreeSRedirector("kalmanfitTest.root");
1009  for (Int_t i=0;i<ntracks;i++){
1010  if (i%100==0) printf("%d\n",i);
1011  Double_t alpha = gRandom->Rndm()*TMath::TwoPi();
1012  Double_t y0 = (gRandom->Rndm()-0.5)*180;
1013  Double_t z0 = (gRandom->Rndm()-0.5)*250*2;
1014  Double_t ky = (gRandom->Rndm()-0.5)*1;
1015  Double_t kz = (gRandom->Rndm()-0.5)*1;
1016  //generate random TPC track
1017  AliTrackPointArray * array = AliTPCkalmanFit::MakePointArrayLinear(alpha,y0,z0, ky, kz,0.04);
1018  AliTrackPointArray * arrayB = new AliTrackPointArray(*array);
1019  kalmanFit2->ApplyCalibration(array,1.); // misalign ideal track
1020  for (Int_t icalib=0; icalib<ncalibs; icalib++){
1021  err[icalib] = TMath::Sqrt((*kalmanFit0->fCalibCovar)(icalib,icalib));
1022  }
1023  //
1024  (*pcstream)<<"dump0"<<
1025  "alpha="<<alpha<<
1026  "y0="<<y0<<
1027  "z0="<<z0<<
1028  "ky="<<ky<<
1029  "lz="<<kz<<
1030  "p.="<<array<<
1031  "pB.="<<arrayB<<
1032  "cparam.="<<kalmanFit0->fCalibParam<<
1033  "ccovar.="<<kalmanFit0->fCalibCovar<<
1034  "err.="<<&err<<
1035  "gparam.="<<kalmanFit2->fCalibParam<<
1036  "gcovar.="<<kalmanFit2->fCalibCovar<<
1037  "\n";
1038  if (i%20==0) {
1039  kalmanFit0->FitTrackLinear(*array,pcstream); // fit track - dump intermediate results
1040  }else{
1041  kalmanFit0->FitTrackLinear(*array,0); // fit track + calibration
1042  }
1043  kalmanFit0->DumpTrackLinear(*array,pcstream); // dump track residuals to the tree + fill histograms
1044  }
1045  pcstream->GetFile()->cd();
1046  kalmanFit0->Write("kalmanFit");
1047  delete pcstream;
1048  return kalmanFit0;
1049 }
1050 
1051 
1052 // Double_t AliTPCkalmanFit::GetTPCDeltaXYZ(Int_t coord, Int_t volID, Double_t x, Double_t y, Double_t z){
1053 // //
1054 // // function for visualization purposes
1055 // //
1056 // if (!fCalibration) return 0;
1057 // Int_t ncalibs = fCalibration->GetEntries();
1058 // if (ncalibs==0) return 0;
1059 // Double_t dxdydz[3]={0,0,0};
1060 // //
1061 // if (volID<0){
1062 // Double_t alpha = TMath::ATan2(y,x);
1063 // Double_t r = TMath::Sqrt(y*y+x*x);
1064 // volID = TMath::Nint(9*alpha/TMath::Pi()-0.5);
1065 // if (volID<0) volID+=18;
1066 // if (z<0) volID+=18;
1067 // if (r>120) volID+=36;
1068 // }
1069 // for (Int_t icalib=0; icalib<ncalibs; icalib++){
1070 // AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
1071 // Double_t param = (*fCalibParam)(icalib,0);
1072 // dxdydz[coord] += transform->GetDeltaXYZ(coord,volID, param, x, y,z);
1073 // }
1074 // return dxdydz[coord];
1075 // }
1076 
1077 // Double_t AliTPCkalmanFit::SGetTPCDeltaXYZ(Int_t coord, Int_t volID, Double_t x, Double_t y, Double_t z){
1078 // //
1079 // //
1080 // //
1081 // if (AliTPCkalmanFit::fgInstance==0) return 0;
1082 // return AliTPCkalmanFit::fgInstance->GetTPCDeltaXYZ(coord, volID,x,y,z);
1083 // }
1084 
1085 
1086 
1087 
1088 
1089 Double_t AliTPCkalmanFit::GetTPCDeltaXYZ(Int_t coord, Int_t volID, Int_t icoordsys, Double_t x, Double_t y, Double_t z){
1103 
1104  if (!fCalibration) return 0;
1105  Int_t ncalibs = fCalibration->GetEntries();
1106  if (ncalibs==0) return 0;
1107  Double_t xyz[3]={0,0,0};
1108  Double_t dxdydz[6]={0,0,0,0,0,0};
1109  Double_t alpha=0;
1110  Double_t r=0;
1111  if(icoordsys==0){alpha=TMath::ATan2(y,x); r =TMath::Sqrt(y*y+x*x);}
1112  if(icoordsys==1){alpha=y; r=x;}
1113  Double_t ca = TMath::Cos(alpha);
1114  Double_t sa = TMath::Sin(alpha);
1115  if(icoordsys==0){xyz[0]=x; xyz[1]=y; xyz[2]=z;}
1116  if(icoordsys==1){xyz[0]=x*ca; xyz[1]=x*sa; xyz[2]=z;}
1117  //
1118  if (volID<0){
1119  // Double_t alpha = TMath::ATan2(y,x);
1120  //Double_t r = TMath::Sqrt(y*y+x*x);
1121  volID = TMath::Nint(9*alpha/TMath::Pi()-0.5);
1122  if (volID<0) volID+=18;
1123  if (z<0) volID+=18;
1124  if (r>120) volID+=36;
1125  }
1126  for (Int_t icalib=0; icalib<ncalibs; icalib++){
1127  AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
1128  Double_t param = (*fCalibParam)(icalib,0);
1129  for (Int_t icoord=0;icoord<6;icoord++){
1130  dxdydz[icoord] += transform->GetDeltaXYZ(icoord,volID, param, xyz[0],xyz[1],xyz[2]);
1131  }
1132  }
1133 
1134  return dxdydz[coord];
1135 }
1136 
1137 
1138 Double_t AliTPCkalmanFit::SGetTPCDeltaXYZ(Int_t coord, Int_t volID, Int_t icoordsys, Double_t x, Double_t y, Double_t z){
1140 
1141  if (AliTPCkalmanFit::fgInstance==0) return 0;
1142  return AliTPCkalmanFit::fgInstance->GetTPCDeltaXYZ(coord, volID, icoordsys,x,y,z);
1143 }
1144 
1145 
1146 
1147 
1148 
1149 Double_t AliTPCkalmanFit::GetTPCtransXYZ(Int_t coord, Int_t volID, Int_t calibID, Int_t icoordsys, Double_t x, Double_t y, Double_t z){
1150 
1151  Int_t ncalibs = fCalibration->GetEntries();
1152  if (calibID>=ncalibs) return 0;
1153  //Int_t volID=-1;
1154  //Double_t xyz[3]={x,y,z};
1155  Double_t r=0;
1156  Double_t alpha=0;
1157  if(icoordsys==0){r=TMath::Sqrt(x*x+y*y); alpha = TMath::ATan2(y,x);}
1158  if(icoordsys==1){r=x; alpha = y;}
1159  Double_t ca = TMath::Cos(alpha);
1160  Double_t sa = TMath::Sin(alpha);
1161  Double_t xyz[3]={0,0,0};
1162  if(icoordsys==0){xyz[0]=x;xyz[1]=y;xyz[2]=z;}
1163  if(icoordsys==1){xyz[0]=x*ca; xyz[1]=x*sa; xyz[2]=z;}
1164  //xyz[3]=param; xyz[4]=volID;
1165 
1166  if (volID<0){
1167  //Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
1168  //Double_t r = TMath::Sqrt(xyz[1]*xyz[1]+xyz[0]*xyz[0]);
1169  volID = TMath::Nint(9*alpha/TMath::Pi()-0.5);
1170  if (volID<0) volID+=18;
1171  if (xyz[2]<0) volID+=18;
1172  if (r>120) volID+=36;
1173  }
1174  AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(calibID);
1175  //transform->SetInstance(transform);
1176  Double_t param = (*fCalibParam)(calibID,0);
1177  Double_t delta = (Double_t)transform->GetDeltaXYZ(coord,volID, param, xyz[0],xyz[1],xyz[2]);
1178 
1179  return delta;
1180 }
1181 
1182 Double_t AliTPCkalmanFit::SGetTPCtransXYZ(Int_t coord, Int_t volID, Int_t calibID, Int_t icoordsys, Double_t x, Double_t y, Double_t z){
1184 
1185  if (AliTPCkalmanFit::fgInstance==0) return 0;
1186  return AliTPCkalmanFit::fgInstance->GetTPCtransXYZ(coord, volID, calibID,icoordsys,x,y,z);
1187 }
1188 
1189 
1190 void AliTPCkalmanFit::MakeTreeTrans(TTreeSRedirector *debug, const char *treeName){
1192 
1193  if (!fCalibParam) {
1194  AliError("Kalman Fit not initialized");
1195  return;
1196  }
1197  //
1198  //
1199  //
1200  const Int_t ncalibs = fCalibration->GetEntries();
1201  TMatrixD dxdydz(ncalibs,5);
1202  Double_t * adx = new Double_t[ncalibs];
1203  Double_t * ady = new Double_t[ncalibs];
1204  Double_t * adz = new Double_t[ncalibs];
1205  Double_t * adr = new Double_t[ncalibs];
1206  Double_t * adrphi = new Double_t[ncalibs];
1207 
1208  Double_t x[3];
1209  for (x[0]=-250.;x[0]<=250.;x[0]+=10.){
1210  for (x[1]=-250.;x[1]<=250.;x[1]+=10.){
1211  for (x[2]=-250.;x[2]<=250.;x[2]+=20.) {
1212  Double_t r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
1213  if (r<20) continue;
1214  if (r>260) continue;
1215  //Double_t z = x[2];
1216  Double_t phi=TMath::ATan2(x[1],x[0]);
1217  Double_t ca=TMath::Cos(phi);
1218  Double_t sa=TMath::Sin(phi);
1219  Double_t dx=0;
1220  Double_t dy=0;
1221  Double_t dz=0;
1222  Double_t dr=0;
1223  Double_t rdphi=0;
1224 
1225  Int_t volID= TMath::Nint(9*phi/TMath::Pi()-0.5);
1226  if (volID<0) volID+=18;
1227  if (x[2]<0) volID+=18; //C-side
1228  if (r>120) volID+=36; //outer
1229  Double_t volalpha=(volID+0.5)*TMath::Pi()/9;
1230  Double_t cva=TMath::Cos(volalpha);
1231  Double_t sva=TMath::Sin(volalpha);
1232 
1233  Double_t lx=x[0]*cva+x[1]*sva;
1234  Double_t ly=-x[0]*sva+x[1]*cva;
1235 
1236 
1237  for(Int_t icalib=0;icalib<ncalibs;icalib++){
1238  for(Int_t icoord=0;icoord<5;icoord++){
1239  dxdydz(icalib,icoord)= GetTPCtransXYZ(icoord, -1, icalib, 0, x[0], x[1], x[2]);
1240  }
1241  }
1242  dx=GetTPCDeltaXYZ(0, -1, 0, x[0], x[1], x[2]);
1243  dy=GetTPCDeltaXYZ(1, -1, 0, x[0], x[1], x[2]);
1244  dz=GetTPCDeltaXYZ(2, -1, 0, x[0], x[1], x[2]);
1245  dr=GetTPCDeltaXYZ(3, -1, 0, x[0], x[1], x[2]);
1246  rdphi=GetTPCDeltaXYZ(4, -1, 0, x[0], x[1], x[2]);
1247 
1248 
1249  if(debug){
1250  TTreeStream &cstream=
1251  (*debug)<<treeName<<
1252  "x="<<x[0]<<
1253  "y="<<x[1]<<
1254  "z="<<x[2]<<
1255  "r="<<r<<
1256  "ca="<<ca<<
1257  "sa="<<sa<<
1258  "lx="<<lx<<
1259  "ly="<<ly<<
1260  "sector="<<volID<<
1261  "phi="<<phi<<
1262  "dx="<<dx<<
1263  "dy="<<dy<<
1264  "dz="<<dz<<
1265  "dr="<<dr<<
1266  "rdphi="<<rdphi<<
1267  "dxdydz.="<<&dxdydz;
1268  for (Int_t icalib=0;icalib<ncalibs;icalib++){
1269  AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
1270  char tname[1000];
1271  //
1272  snprintf(tname,1000,"dx%s=",transform->GetName());
1273  adx[icalib] =dxdydz(icalib,0);
1274  cstream<<tname<<adx[icalib];
1275  snprintf(tname,1000,"dy%s=",transform->GetName());
1276  ady[icalib] =dxdydz(icalib,1);
1277  cstream<<tname<<ady[icalib];
1278  snprintf(tname,1000,"dz%s=",transform->GetName());
1279  adz[icalib] =dxdydz(icalib,2);
1280  cstream<<tname<<adz[icalib];
1281  //
1282  snprintf(tname,1000,"dr%s=",transform->GetName());
1283  adr[icalib] =dxdydz(icalib,3);
1284  cstream<<tname<<adr[icalib];
1285  snprintf(tname,1000,"rdphi%s=",transform->GetName());
1286  adrphi[icalib] =dxdydz(icalib,4);
1287  cstream<<tname<<adrphi[icalib];
1288  }
1289  cstream<<"\n";
1290  }
1291  }
1292  }
1293  Printf("x0=%f finished",x[0]);
1294  }
1295  delete [] adx;// = new Double_t[ncalibs];
1296  delete [] ady;// = new Double_t[ncalibs];
1297  delete [] adz;// = new Double_t[ncalibs];
1298  delete [] adr;// = new Double_t[ncalibs];
1299  delete [] adrphi;// = new Double_t[ncalibs];
1300 
1301 }
THnSparse * fLinearTrackDelta[12]
linear tracks matching residuals - delta
AliTPCkalmanFit * Test(Int_t ntracks)
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
void FitTrackLinear(AliTrackPointArray &points, TTreeSRedirector *debug=0, Float_t scalingRMSY=1., Float_t scalingRMSZ=1.)
Double_t GetTPCDeltaXYZ(Int_t coord, Int_t volID, Int_t icoordsys, Double_t x, Double_t y, Double_t z)
void Propagate(TTreeSRedirector *debug=0)
#define TObjArray
Double_t GetSigma() const
Double_t fCA
! cosine of current angle
static AliTrackPointArray * MakePointArrayLinear(Double_t alpha, Double_t y0, Double_t z0, Double_t ky, Double_t kz, Double_t err=0.02)
TMatrixD * fCalibCovar
calibration parameters
TMatrixD * fLinearCovar
linear covariance
TMatrixD mat
Definition: AnalyzeLaser.C:9
Int_t fLastTimeStamp
last time stamp - used for propagation of parameters
TTreeSRedirector * pcstream
static Double_t SGetTPCDeltaXYZ(Int_t coord, Int_t volID, Int_t icoordsys, Double_t x, Double_t y, Double_t z)
void AddCovariance(const char *varName, Double_t sigma)
npoints
Definition: driftITSTPC.C:85
void PropagateTime(Int_t time)
TObjArray * array
Definition: AnalyzeLaser.C:12
ClassImp(TPCGenInfo)
Definition: AliTPCCmpNG.C:254
void ApplyCalibration(AliTrackPointArray *array, Double_t csign)
void UpdateLinear(AliTrackPoint &point, TTreeSRedirector *debug=0)
void sum()
Bool_t DumpCorelation(Double_t threshold, const char *mask0=0, const char *mask1=0)
Int_t GetTransformationIndex(const char *trName)
static Double_t SGetTPCtransXYZ(Int_t coord, Int_t volID, Int_t calibID, Int_t icoordsys, Double_t x, Double_t y, Double_t z)
THnSparse * fLinearTrackPull[12]
linear tracks matching residuals - pull
Double_t GetTPCtransXYZ(Int_t coord, Int_t volID, Int_t calibID, Int_t icoordsys, Double_t x, Double_t y, Double_t z)
void SetParams(Double_t param, Double_t sigma, Double_t sigma2Time, const TVectorD *const fixedParams)
void SetActive(Bool_t flag)
void SetParam(Double_t param)
Should represent general non linear transformation.
static AliTPCkalmanFit * fgInstance
! Instance of this class (singleton implementation)
Bool_t DumpCalib(const char *mask=0, Float_t correlationCut=-1)
Double_t GetSigmaMax() const
void DumpTrackLinear(AliTrackPointArray &points, TTreeSRedirector *debug)
Double_t GetParam() const
static AliTrackPointArray * SortPoints(AliTrackPointArray &points)
void AddCalibration(AliTPCTransformation *calib)
Double_t GetSigma2Time() const
virtual Double_t GetDeltaXYZ(Int_t coord, Int_t volID, Double_t param, Double_t x, Double_t y, Double_t z)
TMatrixD * fCalibParam
calibration parameters
AliTPCTransformation * GetTransformation(Int_t i)
Double_t fSA
! sinus of current angle
void Update(const AliTPCkalmanFit *kalman)
static TBits * BitsSide(Bool_t aside)
TObjArray * fCalibration
array of calibrations
TMatrixD * fLinearParam
linear parameters
Int_t debug
Double_t fCurrentAlpha
! current rotation frame
void Add(const AliTPCkalmanFit *kalman)
void SetStatus(const char *name, Bool_t setOn, Bool_t isOr=kTRUE)
void MakeTreeTrans(TTreeSRedirector *debug, const char *treeName="all")
Bool_t CheckCovariance(TMatrixD &covar, Float_t maxEl)