1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 // Class for global helix fit of a track
19 // The method uses the following idea:
20 //------------------------------------------------------
23 // (x-x0)^2+(y-y0)^2-R^2=0 ===>
25 // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==>
27 // substitution t = 1/(x^2+y^2), u = 2*x*t, v = 2*y*t, D0 = R^2 - x0^2- y0^2
29 // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation
31 // next substition a = 1/y0 b = -x0/y0 c = -D0/y0
33 // final linear equation : a + u*b +t*c - v =0;
37 // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min;
39 // where sigmai is the error of maesurement (a + ui*b +ti*c - vi)
41 // neglecting error of xi, and supposing xi>>yi sigmai ~ sigmaVi ~ 2*sigmay*t
44 #include "TMatrixDSym.h"
45 //#include "TDecompChol.h"
47 #include "AliRieman.h"
53 AliRieman::AliRieman() :
72 // default constructor
74 for (Int_t i=0;i<6;i++) fParams[i] = 0;
75 for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
76 for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
77 for (Int_t i=0;i<6;i++) {
84 AliRieman::AliRieman(Int_t capacity) :
88 fX(new Double_t[fCapacity]),
89 fY(new Double_t[fCapacity]),
90 fZ(new Double_t[fCapacity]),
91 fSy(new Double_t[fCapacity]),
92 fSz(new Double_t[fCapacity]),
93 fCovar(new TMatrixDSym(6)),
94 fCovarPolY(new TMatrixDSym(3)),
95 fCovarPolZ(new TMatrixDSym(2)),
103 // default constructor
105 for (Int_t i=0;i<6;i++) fParams[i] = 0;
106 for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
107 for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
108 for (Int_t i=0;i<6;i++) {
114 AliRieman::AliRieman(const AliRieman &rieman):
116 fCapacity(rieman.fN),
118 fX(new Double_t[fN]),
119 fY(new Double_t[fN]),
120 fZ(new Double_t[fN]),
121 fSy(new Double_t[fN]),
122 fSz(new Double_t[fN]),
123 fCovar(new TMatrixDSym(*(rieman.fCovar))),
124 fCovarPolY(new TMatrixDSym(*(rieman.fCovarPolY))),
125 fCovarPolZ(new TMatrixDSym(*(rieman.fCovarPolZ))),
126 fSumZZ(rieman.fSumZZ),
128 fChi2Y(rieman.fChi2Y),
129 fChi2Z(rieman.fChi2Z),
136 for (Int_t i=0;i<6;i++) fParams[i] = rieman.fParams[i];
137 for (Int_t i=0;i<9;i++) fSumXY[i] = rieman.fSumXY[i];
138 for (Int_t i=0;i<9;i++) fSumXZ[i] = rieman.fSumXZ[i];
139 for (Int_t i=0;i<6;i++) {
140 fSumPolY[i]=rieman.fSumPolY[i];
141 fSumPolZ[i]=rieman.fSumPolZ[i];
143 for (Int_t i=0;i<rieman.fN;i++){
144 fX[i] = rieman.fX[i];
145 fY[i] = rieman.fY[i];
146 fZ[i] = rieman.fZ[i];
147 fSy[i] = rieman.fSy[i];
148 fSz[i] = rieman.fSz[i];
152 AliRieman::~AliRieman()
167 void AliRieman::Reset()
170 // Reset all the data members
173 for (Int_t i=0;i<6;i++) fParams[i] = 0;
174 for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
175 for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
176 for (Int_t i=0;i<6;i++) {
185 void AliRieman::AddPoint(Double_t x, Double_t y, Double_t z, Double_t sy, Double_t sz)
190 //------------------------------------------------------
193 // (x-x0)^2+(y-y0)^2-R^2=0 ===>
195 // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==>
197 // substitution t = 1/(x^2+y^2), u = 2*x*t, v = 2*y*t, D0 = R^2 - x0^2- y0^2
199 // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation
201 // next substition a = 1/y0 b = -x0/y0 c = -D0/y0
203 // final linear equation : a + u*b +t*c - v =0;
207 // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min;
209 // where sigmai is the error of maesurement (a + ui*b +ti*c - vi)
211 // neglecting error of xi, and supposing xi>>yi sigmai ~ sigmaVi ~ 2*sigmay*t
213 if (fN==fCapacity-1) return; // out of capacity
214 fX[fN] = x; fY[fN]=y; fZ[fN]=z; fSy[fN]=sy; fSz[fN]=sz;
218 Double_t t = x*x+y*y;
223 Double_t error = 2.*sy*t;
225 Double_t weight = 1./error;
227 fSumXY[1] +=u*weight; fSumXY[2]+=v*weight; fSumXY[3]+=t*weight;
228 fSumXY[4] +=u*u*weight; fSumXY[5]+=t*t*weight;
229 fSumXY[6] +=u*v*weight; fSumXY[7]+=u*t*weight; fSumXY[8]+=v*t*weight;
235 Double_t r = TMath::Sqrt(x*x+y*y);
236 fSumXZ[1] +=weight*r; fSumXZ[2] +=weight*r*r; fSumXZ[3] +=weight*z; fSumXZ[4] += weight*r*z;
237 // Now the auxulary sums
238 fSumXZ[5] +=weight*r*r*r/24; fSumXZ[6] +=weight*r*r*r*r/12; fSumXZ[7] +=weight*r*r*r*z/24;
239 fSumXZ[8] +=weight*r*r*r*r*r*r/(24*24);
240 fSumZZ += z*z*weight;
242 // sum accumulation for rough error estimates of the track extrapolation error
244 Double_t maty = 1./(sy*sy);
245 Double_t matz = 1./(sz*sz);
246 for (Int_t i=0;i<5; i++){
256 void AliRieman::UpdatePol(){
262 for (Int_t i=0;i<6;i++)fParams[i]=0;
267 TMatrixDSym smatrix(3);
270 // smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2;
271 // smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut;
272 // sums(0,0) = sv; sums(0,1)=suv; sums(0,2)=svt;
274 smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5];
275 smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7];
276 smatrix(1,0) = fSumXY[1]; smatrix(2,0)=fSumXY[3]; smatrix(2,1)=fSumXY[7];
277 sums(0,0) = fSumXY[2]; sums(0,1) =fSumXY[6]; sums(0,2) =fSumXY[8];
279 if (smatrix.IsValid()){
280 for (Int_t i=0;i<3;i++)
281 for (Int_t j=0;j<=i;j++){
282 (*fCovar)(i,j)=smatrix(i,j);
284 TMatrixD res = sums*smatrix;
285 fParams[0] = res(0,0);
286 fParams[1] = res(0,1);
287 fParams[2] = res(0,2);
293 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
295 // fSumXZ[1] += fSumXZ[5]*rm1*rm1;
296 // fSumXZ[2] += fSumXZ[6]*rm1*rm1 + fSumXZ[8]*rm1*rm1*rm1*rm1;
297 // fSumXZ[4] += fSumXZ[7]*rm1*rm1;
298 Double_t sum1 = fSumXZ[1] + fSumXZ[5]*rm1*rm1;
299 Double_t sum2 = fSumXZ[2] + fSumXZ[6]*rm1*rm1 + fSumXZ[8]*rm1*rm1*rm1*rm1;
300 Double_t sum4 = fSumXZ[4] + fSumXZ[7]*rm1*rm1;
302 TMatrixDSym smatrixz(2);
303 // smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = fSumXZ[1]; smatrixz(1,1) = fSumXZ[2];
304 //smatrixz(1,0) = fSumXZ[1];
305 smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = sum1; smatrixz(1,1) = sum2;
306 smatrixz(1,0) = sum1;
308 TMatrixD sumsxz(1,2);
309 if (smatrixz.IsValid()){
310 sumsxz(0,0) = fSumXZ[3];
311 // sumsxz(0,1) = fSumXZ[4];
313 TMatrixD res = sumsxz*smatrixz;
314 fParams[3] = res(0,0);
315 fParams[4] = res(0,1);
317 for (Int_t i=0;i<2;i++)
318 for (Int_t j=0;j<=i;j++){
319 (*fCovar)(i+3,j+3)=smatrixz(i,j);
323 UpdateCovariancePol();
324 // (x-x0)^2+(y-y0)^2-R^2=0 ===>
326 // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==>
327 // substitution t = 1/(x^2+y^2), u = 2*x*t, y = 2*y*t, D0 = R^2 - x0^2- y0^2
329 // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation
331 // next substition a = 1/y0 b = -x0/y0 c = -D0/y0
332 // final linear equation : a + u*b +t*c - v =0;
335 // fParam[1] = -x0/y0
336 // fParam[2] = - (R^2 - x0^2 - y0^2)/y0
337 if (conv>1) fConv =kTRUE;
342 void AliRieman::Update(){
348 for (Int_t i=0;i<6;i++)fParams[i]=0;
353 TMatrixDSym smatrix(3);
356 // smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2;
357 // smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut;
358 // sums(0,0) = sv; sums(0,1)=suv; sums(0,2)=svt;
360 smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5];
361 smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7];
363 smatrix(1,0) = fSumXY[1];
364 smatrix(2,0) = fSumXY[3];
365 smatrix(2,1) = fSumXY[7];
367 sums(0,0) = fSumXY[2]; sums(0,1) =fSumXY[6]; sums(0,2) =fSumXY[8];
368 //TDecompChol decomp(smatrix);
369 //decomp.SetTol(1.0e-14);
373 if (smatrix.IsValid()){
374 for (Int_t i=0;i<3;i++)
375 for (Int_t j=0;j<3;j++){
376 (*fCovar)(i,j)=smatrix(i,j);
378 TMatrixD res = sums*smatrix;
379 fParams[0] = res(0,0);
380 fParams[1] = res(0,1);
381 fParams[2] = res(0,2);
385 fConv = kFALSE; //non converged
388 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.<0){
389 fConv = kFALSE; //non converged
395 Double_t x0 = -fParams[1]/fParams[0];
396 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.);
399 for (Int_t i=0;i<9;i++) sumXZ[i]=0;
400 for (Int_t i=0;i<fN;i++){
401 Double_t phi = TMath::ASin((fX[i]-x0)*rm1);
402 Double_t phi0 = TMath::ASin((0.-x0)*rm1);
403 Double_t weight = 1/fSz[i];
405 Double_t dphi = (phi-phi0)/rm1;
407 sumXZ[1] +=weight*dphi;
408 sumXZ[2] +=weight*dphi*dphi;
409 sumXZ[3] +=weight*fZ[i];
410 sumXZ[4] +=weight*dphi*fZ[i];
414 TMatrixDSym smatrixz(2);
416 smatrixz(0,0) = sumXZ[0]; smatrixz(0,1) = sumXZ[1]; smatrixz(1,1) = sumXZ[2];
417 smatrixz(1,0) = sumXZ[1];
419 if (smatrixz.IsValid()){
420 sumsz(0,0) = sumXZ[3];
421 sumsz(0,1) = sumXZ[4];
422 TMatrixD res = sumsz*smatrixz;
423 fParams[3] = res(0,0);
424 fParams[4] = res(0,1);
425 for (Int_t i=0;i<2;i++)
426 for (Int_t j=0;j<2;j++){
427 (*fCovar)(i+3,j+3)=smatrixz(i,j);
431 UpdateCovariancePol();
432 // (x-x0)^2+(y-y0)^2-R^2=0 ===>
434 // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==>
435 // substitution t = 1/(x^2+y^2), u = 2*x*t, v = 2*y*t, D0 = R^2 - x0^2- y0^2
437 // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation
439 // next substition a = 1/y0 b = -x0/y0 c = -D0/y0
440 // final linear equation : a + u*b +t*c - v =0;
443 // fParam[1] = -x0/y0
444 // fParam[2] = - (R^2 - x0^2 - y0^2)/y0
445 if (conv>1) fConv =kTRUE;
448 fChi2Y = CalcChi2Y();
449 fChi2Z = CalcChi2Z();
450 fChi2 = fChi2Y +fChi2Z;
453 void AliRieman::UpdateCovariancePol(){
455 // covariance matrices for rough extrapolation error estimate
456 // covariance matrix - get error estimates at point x = 0
457 // ! NOTE - error estimates is very rough and it is valid only if dl<<R
458 // when dl is the distance between first and last point
461 (*fCovarPolY)(0,0) = fSumPolY[0]; (*fCovarPolY)(0,1) = fSumPolY[1]; (*fCovarPolY)(0,2) = fSumPolY[2];
462 (*fCovarPolY)(1,0) = fSumPolY[1]; (*fCovarPolY)(1,1) = fSumPolY[2]; (*fCovarPolY)(1,2) = fSumPolY[3];
463 (*fCovarPolY)(2,0) = fSumPolY[2]; (*fCovarPolY)(2,1) = fSumPolY[3]; (*fCovarPolY)(2,2) = fSumPolY[4];
464 fCovarPolY->Invert();
467 (*fCovarPolZ)(0,0) = fSumPolZ[0]; (*fCovarPolZ)(0,1) = fSumPolZ[1];
468 (*fCovarPolZ)(1,0) = fSumPolZ[1]; (*fCovarPolZ)(1,1) = fSumPolZ[2];
469 fCovarPolZ->Invert();
473 Double_t AliRieman::GetErrY(Double_t x) const{
475 // P0' = P0 + P1 * x + P2 * x^2
486 TMatrixD covarX(trans,TMatrixD::kMult,*fCovarPolY);
491 Double_t AliRieman::GetErrZ(Double_t x) const{
493 // assumption error of curvature determination neggligible
502 TMatrixD covarX(trans,TMatrixD::kMult,*fCovarPolZ);
507 Bool_t AliRieman::GetExternalParameters(Double_t xref, Double_t *params, Double_t * covar){
509 // Get external parameters
510 // + approximative covariance
511 // 0 1 2 3 4 // y - local y
512 // 5 6 7 8 // z - local z
513 // 9 10 11 // snp - local sin fi
514 // 12 13 // tgl - deep angle
515 // 14 // C - curvature
516 Double_t sign = (GetC()>0) ? 1.:-1.;
518 params[0] = GetYat(xref);
519 params[1] = GetZat(xref);
520 params[2] = TMath::Sin(TMath::ATan(GetDYat(xref)));
521 params[3] = sign*fParams[4];
524 // covariance matrix in y
526 TMatrixD transY(3,3);
529 transY(0,2) = xref*xref;
533 TMatrixD covarY(transY,TMatrixD::kMult,*fCovarPolY);
536 TMatrixD transZ(2,2);
540 TMatrixD covarZ(transZ,TMatrixD::kMult,*fCovarPolZ);
543 // C ~ 2*P2 - in rotated frame
545 covar[0] = covarY(0,0); covar[1] = 0; covar[2] = covarY(0,1); covar[3] = 0; covar[4] = TMath::Sqrt(2.)*covarY(0,2);
546 covar[5] = covarZ(0,0); covar[6] = 0; covar[7] = sign*covarZ(0,1); covar[8] = 0;
547 covar[9] = covarY(1,1); covar[10]= 0; covar[11]= TMath::Sqrt(2.)*covarY(1,2);
548 covar[12] = covarZ(1,1); covar[13]= 0;
549 covar[14] = 2.*covarY(2,2);
557 Double_t AliRieman::GetYat(Double_t x) const {
559 // get y at given x position
561 if (!fConv) return 0.;
562 Double_t res2 = (x*fParams[0]+fParams[1]);
564 res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
565 if (res2<0) return 0;
566 res2 = TMath::Sqrt(res2);
567 res2 = (1-res2)/fParams[0];
571 Double_t AliRieman::GetDYat(Double_t x) const {
573 // get dy/dx at given x position
575 if (!fConv) return 0.;
576 Double_t x0 = -fParams[1]/fParams[0];
577 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<0) return 0;
578 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
579 if ( 1./(rm1*rm1)-(x-x0)*(x-x0)<=0) return 0;
580 Double_t res = (x-x0)/TMath::Sqrt(1./(rm1*rm1)-(x-x0)*(x-x0));
581 if (fParams[0]<0) res*=-1.;
587 Double_t AliRieman::GetZat(Double_t x) const {
589 // get z at given x position
591 if (!fConv) return 0.;
592 Double_t x0 = -fParams[1]/fParams[0];
593 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
594 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
595 Double_t phi = TMath::ASin((x-x0)*rm1);
596 Double_t phi0 = TMath::ASin((0.-x0)*rm1);
597 Double_t dphi = (phi-phi0);
598 Double_t res = fParams[3]+fParams[4]*dphi/rm1;
602 Double_t AliRieman::GetDZat(Double_t x) const {
604 // get dz/dx at given x postion
606 if (!fConv) return 0.;
607 Double_t x0 = -fParams[1]/fParams[0];
608 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
609 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
610 Double_t res = fParams[4]/TMath::Cos(TMath::ASin((x-x0)*rm1));
615 //_____________________________________________________________________________
616 Bool_t AliRieman::GetXYZat(Double_t r, Double_t alpha, Float_t *xyz) const
619 // Returns position given radius
621 if (!fConv) return kFALSE;
622 Double_t res2 = (r*fParams[0]+fParams[1]);
624 res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
625 if (res2<0) return kFALSE;
626 res2 = TMath::Sqrt(res2);
627 res2 = (1-res2)/fParams[0];
629 if (!fConv) return kFALSE;
630 Double_t x0 = -fParams[1]/fParams[0];
631 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
632 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
633 Double_t phi = TMath::ASin((r-x0)*rm1);
634 Double_t phi0 = TMath::ASin((0.-x0)*rm1);
635 Double_t dphi = (phi-phi0);
636 Double_t res = fParams[3]+fParams[4]*dphi/rm1;
638 Double_t sin = TMath::Sin(alpha);
639 Double_t cos = TMath::Cos(alpha);
640 xyz[0] = r*cos - res2*sin;
641 xyz[1] = res2*cos + r*sin;
648 Double_t AliRieman::GetC() const{
652 return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
655 Double_t AliRieman::CalcChi2Y() const{
657 // calculate chi2 for Y
659 Double_t sumchi2 = 0;
660 for (Int_t i=0;i<fN;i++){
661 Double_t chi2 = (fY[i] - GetYat(fX[i]))/fSy[i];
668 Double_t AliRieman::CalcChi2Z() const{
670 // calculate chi2 for Z
672 Double_t sumchi2 = 0;
673 for (Int_t i=0;i<fN;i++){
674 Double_t chi2 = (fZ[i] - GetZat(fX[i]))/fSy[i];
680 Double_t AliRieman::CalcChi2() const{
682 // sum chi2 in both coord - supposing Y and Z independent
684 return CalcChi2Y()+CalcChi2Z();
687 AliRieman * AliRieman::MakeResiduals() const{
689 // create residual structure - ONLY for debug purposes
691 AliRieman *rieman = new AliRieman(fN);
692 for (Int_t i=0;i<fN;i++){
693 rieman->AddPoint(fX[i],fY[i]-GetYat(fX[i]),fZ[i]-GetZat(fX[i]), fSy[i],fSz[i]);