6 Kalman filter(s) for fitting of the tracks together with calibration/transformation
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.
13 Current calibration parameters and covariance stored (fCalibParam, fCalibCovar).
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
23 How to use it - see AliTPCkalmanFit::Test function
25 Simple test (see AliTPCkalmanFit::Test)
26 AliTPCTransformation::BuildBasicFormulas();
27 AliTPCkalmanFit *kalmanFit0 = AliTPCkalmanFit::Test(2000);
28 TFile f("kalmanfitTest.root");
30 Transformation visualization:
31 Transforamtion can be visualized using TFn (TF1,TF2 ...) and using tree->Draw()
33 kalmanFit0->SetInstance(kalmanFit0); //
34 kalmanFit0->InitTransformation(); //
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);
39 TF2 fxRZ("fxRZ","sign(y)*10*(AliTPCkalmanFit::SGetTPCDeltaXYZ(0,-1,x,0,y))",85,245,-250,250);
51 #include "TLinearFitter.h"
53 #include "THnSparse.h"
57 #include "TTreeStream.h"
58 #include "AliTrackPointArray.h"
60 #include "AliTPCTransformation.h"
61 #include "AliTPCkalmanFit.h"
63 ClassImp(AliTPCkalmanFit)
65 AliTPCkalmanFit* AliTPCkalmanFit::fgInstance = 0;
67 AliTPCkalmanFit::AliTPCkalmanFit():
75 fCurrentAlpha(0), //! current rotation frame
76 fCA(0), //! cosine of current angle
77 fSA(0) //! sinus of current angle
80 // Default constructor
82 for (Int_t ihis=0; ihis<12; ihis++){
83 fLinearTrackDelta[ihis]=0;
84 fLinearTrackPull[ihis]=0;
88 void AliTPCkalmanFit::InitTransformation(){
90 // Initialize pointers to the transforamtion functions
93 Int_t ncalibs = fCalibration->GetEntries();
94 for (Int_t icalib=0;icalib<ncalibs; icalib++){
95 AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
100 void AliTPCkalmanFit::Add(const AliTPCkalmanFit * kalman){
105 for (Int_t i=0;i<12;i++){
106 if (fLinearTrackDelta[i] && kalman->fLinearTrackDelta[i]){
107 fLinearTrackDelta[i]->Add(kalman->fLinearTrackDelta[i]);
109 if (fLinearTrackPull[i] && kalman->fLinearTrackPull[i]){
110 fLinearTrackPull[i]->Add(kalman->fLinearTrackPull[i]);
117 void AliTPCkalmanFit::Init(){
119 // Initialize parameter vector and covariance matrix
120 // To be called after initialization of all of the transformations
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->GetParam();
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->GetSigma()*transform->GetSigma();
135 // Build QA histograms
137 Double_t mpi = TMath::Pi();
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",
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)",
155 "100*#Delta_{#theta}(cm)",
156 "100^{2}*dz0^{2}/dx0^{2}(cm)",
157 "100^{2}*dz1^{2}/dx1^{2}(cm)",
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)",
168 "100*#Delta_{#theta}(unit)",
169 "100^{2}*dz0^{2}/dx0^{2}(unit)",
170 "100^{2}*dz1^{2}/dx1^{2}(unit)"
178 for (Int_t ihis=0; ihis<12; ihis++){
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]);
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]);
200 void AliTPCkalmanFit::SetStatus(const char * mask, Bool_t setOn, Bool_t isOr){
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();
207 for (Int_t i=0; i<ncalibs;i++){
208 AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(i);
209 transform->SetActive(setOn);
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);
225 void AliTPCkalmanFit::Update(const AliTPCkalmanFit * kalman){
227 // Update Kalman filter
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);
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
243 for (Int_t i=0;i<ncalibs;i++){
244 for (Int_t j=0;j<ncalibs;j++) matHk(i,j)=0;
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
251 matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
252 vecXk += matKk*vecYk; // updated vector
253 covXk2= (matHk-(matKk*matHk));
254 covXk3 = covXk2*covXk;
256 (*fCalibParam) = vecXk;
257 (*fCalibCovar) = covXk;
264 void AliTPCkalmanFit::FitTrackLinear(AliTrackPointArray &points, TTreeSRedirector *debug, Float_t scalingRMSY, Float_t scalingRMSZ){
269 AliError("Kalman Fit not initialized");
273 // 1. Make initial track hypothesis
275 TLinearFitter lfitY(2,"pol1");
276 TLinearFitter lfitZ(2,"pol1");
280 lfitY.StoreData(kTRUE);
281 lfitZ.StoreData(kTRUE);
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.);
292 // AliTrackPointArray pointsc(points);
293 //ApplyCalibration(&pointsc, -1.);
295 // 1.a Choosing reference rotation alpha
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);
301 // 1.b Fit the track in the rotated frame - MakeSeed
303 for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
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);
312 lfitY.GetParameters(vecY);
313 lfitZ.GetParameters(vecZ);
315 // 2. Set initial parameters and the covariance matrix
317 Int_t ncalibs = fCalibration->GetEntries();
319 fLinearParam = new TMatrixD(ncalibs+4,1);
320 fLinearCovar = new TMatrixD(ncalibs+4,ncalibs+4);
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);
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;
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;
343 // Fit thetrack together with correction
346 for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
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();
353 cov[1] = (erry2+kOff0)*scalingRMSY;
354 cov[2] = (errz2+kOff0)*scalingRMSZ;
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();
359 cov[1] = (erry2+kOff0)*scalingRMSY;
360 cov[2] = (errz2+kOff0)*scalingRMSZ;
361 UpdateLinear(point,debug);
364 // save current param and covariance
365 for (Int_t i=0; i<ncalibs;i++){
366 AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(i);
367 transform->SetParam( (*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);
373 if (debug){ // dump debug info
374 (*debug)<<"fitLinear"<<
379 "lP.="<<fLinearParam<<
380 "cP.="<<fCalibParam<<
381 "lC.="<<fLinearCovar<<
382 "cC.="<<fCalibCovar<<
387 void AliTPCkalmanFit::DumpTrackLinear(AliTrackPointArray &points, TTreeSRedirector *debug){
389 // Dump the track parameters before and after current calibration
391 // Track is divided to two parts -
392 // X mean is defined as middle point
396 AliError("Kalman Fit not initialized");
402 TLinearFitter *fitters[16];
403 TVectorD *params[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);
418 // calculate mean x point and corrdinate frame
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];
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];
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];
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);
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);
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
469 (*debug)<<"dumpLinear"<<
470 "alpha="<<fCurrentAlpha<<
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]<<
508 // delta alpha y0 z0 ky kz
509 Double_t x[6]={0,0,0,0,0,0};
511 x[2]=(*params[0])[0];
512 x[3]=(*params[2])[0];
513 x[4]=(*params[0])[1];
514 x[5]=(*params[2])[1];
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);
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);
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);
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);
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);
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);
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);
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);
577 fLinearTrackDelta[8]->Fill(x);
579 fLinearTrackDelta[9]->Fill(x);
581 fLinearTrackDelta[10]->Fill(x);
583 fLinearTrackDelta[11]->Fill(x);
587 for (Int_t ifit=0; ifit<8;ifit++){
588 delete fitters[ifit];
590 delete fitters[ifit+8];
591 delete params[ifit+8];
598 void AliTPCkalmanFit::Propagate(TTreeSRedirector */*debug*/){
600 // Propagate the Kalman
604 void AliTPCkalmanFit::AddCovariance(const char * varName, Double_t sigma){
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;
622 void AliTPCkalmanFit::PropagateTime(Int_t time){
624 // Propagate the calibration in time
625 // - Increase covariance matrix
627 if (!fCalibCovar) return;
628 Int_t ncalibs = fCalibration->GetEntries();
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);
634 if ((*fCalibCovar)(icalib,icalib)<transform->GetSigmaMax()*transform->GetSigmaMax())
635 (*fCalibCovar)(icalib,icalib)+= transform->GetSigma2Time()*TMath::Abs(deltaT);
640 void AliTPCkalmanFit::UpdateLinear(AliTrackPoint &point, TTreeSRedirector *debug){
645 Int_t ncalibs = fCalibration->GetEntries();
647 Int_t nelem = ncalibs+4;
648 TMatrixD &vecXk=*fLinearParam; // X vector
649 TMatrixD &covXk=*fLinearCovar; // X covariance
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);
663 for (Int_t iel=0;iel<nelem;iel++){
664 for (Int_t ip=0;ip<kNmeas;ip++) {
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;
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()};
680 Double_t dxdydz[3]={0,0,0};
681 Double_t rdxdydz[3]={0,0,0};
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];
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
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];
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;
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
711 matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
713 covXk2= (mat1-(matKk*matHk));
714 covXk3 = covXk2*covXk;
716 CheckCovariance(covXk3,100);
717 vecXk += matKk*vecYk; // updated vector
719 if (debug){ // dump debug info
720 (*debug)<<"updateLinear"<<
727 "matHkT.="<<&matHkT<<
730 "covXk2.="<<&covXk2<<
738 AliTrackPointArray * AliTPCkalmanFit::SortPoints(AliTrackPointArray &points){
740 //Creates the array - points sorted according radius - neccessay for kalman fit
743 // 0. choose the frame - rotation angle
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);
751 // 1. sort the points
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];
758 TMath::Sort(npoints, rxvector,indexes,kFALSE);
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);
772 AliTrackPointArray * AliTPCkalmanFit::MakePointArrayLinear(Double_t alpha, Double_t y0, Double_t z0, Double_t ky, Double_t kz, Double_t err){
776 AliTrackPointArray array(500);
777 Float_t cov[10]; // dummy covariance
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);
786 for (Float_t ir = -245; ir<245; ir++){
789 if (TMath::Abs(ir)<80) continue;
790 Double_t ca = TMath::Cos(alpha);
791 Double_t sa = TMath::Sin(alpha);
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;
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;
802 if (ir>130) isec+=36;
804 AliTrackPoint point(gx, gy, gz, cov, isec,0,0);
805 array.AddPoint(npoints, &point);
808 AliTrackPointArray *parray = new AliTrackPointArray(npoints);
810 for (Int_t i=0;i<npoints;i++){
811 array.GetPoint(point,i);
812 parray->AddPoint(i,&point);
818 void AliTPCkalmanFit::AddCalibration(AliTPCTransformation * calib){
822 if (!fCalibration) fCalibration = new TObjArray(2000);
823 fCalibration->AddLast(calib);
826 Int_t AliTPCkalmanFit::GetTransformationIndex(const char * trName){
830 if (!fCalibration) return -1;
831 if (!fCalibration->FindObject(trName)) return -1;
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){
846 void AliTPCkalmanFit::ApplyCalibration(AliTrackPointArray *array, Double_t csign){
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->GetParam(), xyz[0], xyz[1], xyz[2]);
861 dxdydz[1] += transform->GetDeltaXYZ(1,volId, transform->GetParam(), xyz[0], xyz[1], xyz[2]);
862 dxdydz[2] += transform->GetDeltaXYZ(2,volId, transform->GetParam(), xyz[0], xyz[1], xyz[2]);
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];
870 Bool_t AliTPCkalmanFit::DumpCorelation(Double_t threshold, const char *mask0, const char *mask1){
874 TMatrixD &mat = *fCalibCovar;
875 Int_t nrow= mat.GetNrows();
876 for (Int_t irow=0; irow<nrow; irow++){
877 AliTPCTransformation * trans0 = GetTransformation(irow);
878 TString strName0(trans0->GetName());
880 if (!strName0.Contains(mask0)) continue;
882 for (Int_t icol=irow+1; icol<nrow; icol++){
883 AliTPCTransformation * trans1 = GetTransformation(icol);
884 TString strName1(trans1->GetName());
886 if (!strName1.Contains(mask1)) continue;
889 Double_t diag = TMath::Sqrt(TMath::Abs(mat(irow,irow)*mat(icol,icol)));
891 printf("Negative covariance\t%d\t%d\t%f\n",irow,icol, mat(irow,icol));
894 Double_t corr0 = mat(irow,icol)/diag;
895 if (TMath::Abs(corr0)>threshold){
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);
904 Bool_t AliTPCkalmanFit::DumpCalib(const char *mask, Float_t correlationCut){
906 // Print calibration entries - name, value, error
908 TMatrixD &mat = *fCalibCovar;
909 Int_t nrow= mat.GetNrows();
910 TString strMask(mask);
911 TVectorD vecCorrSum(nrow);
912 for (Int_t irow=0; irow<nrow; irow++){
914 for (Int_t icol=0; icol<nrow; icol++){
915 if (icol==irow) continue;
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);
920 vecCorrSum[irow]*=0.5;
923 for (Int_t irow=0; irow<nrow; irow++){
924 AliTPCTransformation * trans0 = GetTransformation(irow);
925 TString strName(trans0->GetName());
927 if (!strName.Contains(mask)) continue;
929 if (vecCorrSum[irow]<correlationCut) continue;
930 printf("%d\t%s\t%f\t%f\t%f\n",
933 (*fCalibParam)(irow,0),
934 TMath::Sqrt(mat(irow,irow)), vecCorrSum[irow]);
940 Bool_t AliTPCkalmanFit::CheckCovariance(TMatrixD &mat, Float_t /*maxEl*/){
942 // Check consistency of covariance matrix
943 // + symetrize coavariance matrix
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));
951 // if (mat(irow,irow)>maxEl){
952 // printf("Too big covariance\t%d\t%f\n",irow,mat(irow,irow));
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)));
960 printf("Negative covariance\t%d\t%d\t%f\n",irow,icol, mat(irow,icol));
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);
970 if (TMath::Abs(cov0-cov1)>0.0000001){
971 printf("Asymetry problem %d\t%d\t%f\t%f\n",irow,icol, cov0, cov1);
975 // symetrize the covariance matrix
976 Double_t mean = (mat(irow,icol)+ mat(icol,irow))*0.5;
985 AliTPCkalmanFit * AliTPCkalmanFit::Test(Int_t ntracks){
987 // This is test example
991 // 1. Setup transformation
996 AliTPCTransformation * transformation=0;
997 AliTPCkalmanFit * kalmanFit0 = new AliTPCkalmanFit;
998 AliTPCkalmanFit * kalmanFit2 = new AliTPCkalmanFit;
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++){
1007 if (ipar0+ipar1==0) continue;
1008 Double_t param = (gRandom->Rndm()-0.5)*0.5; // generate random parameters
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);
1021 // 2. Init transformation
1027 // 3. Run kalman filter
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));
1047 (*pcstream)<<"dump0"<<
1055 "cparam.="<<kalmanFit0->fCalibParam<<
1056 "ccovar.="<<kalmanFit0->fCalibCovar<<
1058 "gparam.="<<kalmanFit2->fCalibParam<<
1059 "gcovar.="<<kalmanFit2->fCalibCovar<<
1062 kalmanFit0->FitTrackLinear(*array,pcstream); // fit track - dump intermediate results
1064 kalmanFit0->FitTrackLinear(*array,0); // fit track + calibration
1066 kalmanFit0->DumpTrackLinear(*array,pcstream); // dump track residuals to the tree + fill histograms
1068 pcstream->GetFile()->cd();
1069 kalmanFit0->Write("kalmanFit");
1075 // Double_t AliTPCkalmanFit::GetTPCDeltaXYZ(Int_t coord, Int_t volID, Double_t x, Double_t y, Double_t z){
1077 // // function for visualization purposes
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};
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;
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);
1097 // return dxdydz[coord];
1100 // Double_t AliTPCkalmanFit::SGetTPCDeltaXYZ(Int_t coord, Int_t volID, Double_t x, Double_t y, Double_t z){
1104 // if (AliTPCkalmanFit::fgInstance==0) return 0;
1105 // return AliTPCkalmanFit::fgInstance->GetTPCDeltaXYZ(coord, volID,x,y,z);
1112 Double_t AliTPCkalmanFit::GetTPCDeltaXYZ(Int_t coord, Int_t volID, Int_t icoordsys, Double_t x, Double_t y, Double_t z){
1114 // function for visualization purposes
1116 // coord - coordinate for output
1124 //icoordsys - type of coordinate system for input
1128 if (!fCalibration) return 0;
1129 Int_t ncalibs = fCalibration->GetEntries();
1130 if (ncalibs==0) return 0;
1131 Double_t xyz[3]={0,0,0};
1132 Double_t dxdydz[6]={0,0,0,0,0,0};
1135 if(icoordsys==0){alpha=TMath::ATan2(y,x); r =TMath::Sqrt(y*y+x*x);}
1136 if(icoordsys==1){alpha=y; r=x;}
1137 Double_t ca = TMath::Cos(alpha);
1138 Double_t sa = TMath::Sin(alpha);
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;}
1143 // Double_t alpha = TMath::ATan2(y,x);
1144 //Double_t r = TMath::Sqrt(y*y+x*x);
1145 volID = TMath::Nint(9*alpha/TMath::Pi()-0.5);
1146 if (volID<0) volID+=18;
1148 if (r>120) volID+=36;
1150 for (Int_t icalib=0; icalib<ncalibs; icalib++){
1151 AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
1152 Double_t param = (*fCalibParam)(icalib,0);
1153 for (Int_t icoord=0;icoord<6;icoord++){
1154 dxdydz[icoord] += transform->GetDeltaXYZ(icoord,volID, param, xyz[0],xyz[1],xyz[2]);
1158 return dxdydz[coord];
1162 Double_t AliTPCkalmanFit::SGetTPCDeltaXYZ(Int_t coord, Int_t volID, Int_t icoordsys, Double_t x, Double_t y, Double_t z){
1166 if (AliTPCkalmanFit::fgInstance==0) return 0;
1167 return AliTPCkalmanFit::fgInstance->GetTPCDeltaXYZ(coord, volID, icoordsys,x,y,z);
1174 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){
1176 Int_t ncalibs = fCalibration->GetEntries();
1177 if (calibID>=ncalibs) return 0;
1179 //Double_t xyz[3]={x,y,z};
1182 if(icoordsys==0){r=TMath::Sqrt(x*x+y*y); alpha = TMath::ATan2(y,x);}
1183 if(icoordsys==1){r=x; alpha = y;}
1184 Double_t ca = TMath::Cos(alpha);
1185 Double_t sa = TMath::Sin(alpha);
1186 Double_t xyz[3]={0,0,0};
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;}
1189 //xyz[3]=param; xyz[4]=volID;
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;
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]);
1207 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){
1211 if (AliTPCkalmanFit::fgInstance==0) return 0;
1212 return AliTPCkalmanFit::fgInstance->GetTPCtransXYZ(coord, volID, calibID,icoordsys,x,y,z);
1216 void AliTPCkalmanFit::MakeTreeTrans(TTreeSRedirector *debug, const char *treeName){
1218 // Make the Tree before and after current calibration
1221 AliError("Kalman Fit not initialized");
1227 const Int_t ncalibs = fCalibration->GetEntries();
1228 TMatrixD dxdydz(ncalibs,5);
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];
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]);
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);
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);
1260 Double_t lx=x[0]*cva+x[1]*sva;
1261 Double_t ly=-x[0]*sva+x[1]*cva;
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]);
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]);
1277 TTreeStream &cstream=
1278 (*debug)<<treeName<<
1294 "dxdydz.="<<&dxdydz;
1295 for (Int_t icalib=0;icalib<ncalibs;icalib++){
1296 AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
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];
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];
1320 Printf("x0=%f finished",x[0]);