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>
48 #include "AliRieman.h"
54 AliRieman::AliRieman() :
73 // default constructor
75 for (Int_t i=0;i<6;i++) fParams[i] = 0;
76 for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
77 for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
78 for (Int_t i=0;i<6;i++) {
85 AliRieman::AliRieman(Int_t capacity) :
89 fX(new Double_t[fCapacity]),
90 fY(new Double_t[fCapacity]),
91 fZ(new Double_t[fCapacity]),
92 fSy(new Double_t[fCapacity]),
93 fSz(new Double_t[fCapacity]),
94 fCovar(new TMatrixDSym(6)),
95 fCovarPolY(new TMatrixDSym(3)),
96 fCovarPolZ(new TMatrixDSym(2)),
104 // default constructor
106 for (Int_t i=0;i<6;i++) fParams[i] = 0;
107 for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
108 for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
109 for (Int_t i=0;i<6;i++) {
115 AliRieman::AliRieman(const AliRieman &rieman):
117 fCapacity(rieman.fN),
119 fX(new Double_t[fN]),
120 fY(new Double_t[fN]),
121 fZ(new Double_t[fN]),
122 fSy(new Double_t[fN]),
123 fSz(new Double_t[fN]),
124 fCovar(new TMatrixDSym(*(rieman.fCovar))),
125 fCovarPolY(new TMatrixDSym(*(rieman.fCovarPolY))),
126 fCovarPolZ(new TMatrixDSym(*(rieman.fCovarPolZ))),
127 fSumZZ(rieman.fSumZZ),
129 fChi2Y(rieman.fChi2Y),
130 fChi2Z(rieman.fChi2Z),
137 for (Int_t i=0;i<6;i++) fParams[i] = rieman.fParams[i];
138 for (Int_t i=0;i<9;i++) fSumXY[i] = rieman.fSumXY[i];
139 for (Int_t i=0;i<9;i++) fSumXZ[i] = rieman.fSumXZ[i];
140 for (Int_t i=0;i<6;i++) {
141 fSumPolY[i]=rieman.fSumPolY[i];
142 fSumPolZ[i]=rieman.fSumPolZ[i];
144 for (Int_t i=0;i<rieman.fN;i++){
145 fX[i] = rieman.fX[i];
146 fY[i] = rieman.fY[i];
147 fZ[i] = rieman.fZ[i];
148 fSy[i] = rieman.fSy[i];
149 fSz[i] = rieman.fSz[i];
153 AliRieman::~AliRieman()
168 void AliRieman::Reset()
171 // Reset all the data members
174 for (Int_t i=0;i<6;i++) fParams[i] = 0;
175 for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
176 for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
177 for (Int_t i=0;i<6;i++) {
186 void AliRieman::AddPoint(Double_t x, Double_t y, Double_t z, Double_t sy, Double_t sz)
191 //------------------------------------------------------
194 // (x-x0)^2+(y-y0)^2-R^2=0 ===>
196 // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==>
198 // substitution t = 1/(x^2+y^2), u = 2*x*t, v = 2*y*t, D0 = R^2 - x0^2- y0^2
200 // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation
202 // next substition a = 1/y0 b = -x0/y0 c = -D0/y0
204 // final linear equation : a + u*b +t*c - v =0;
208 // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min;
210 // where sigmai is the error of maesurement (a + ui*b +ti*c - vi)
212 // neglecting error of xi, and supposing xi>>yi sigmai ~ sigmaVi ~ 2*sigmay*t
214 if (fN==fCapacity-1) return; // out of capacity
215 fX[fN] = x; fY[fN]=y; fZ[fN]=z; fSy[fN]=sy; fSz[fN]=sz;
219 Double_t t = x*x+y*y;
224 Double_t error = 2.*sy*t;
226 Double_t weight = 1./error;
228 fSumXY[1] +=u*weight; fSumXY[2]+=v*weight; fSumXY[3]+=t*weight;
229 fSumXY[4] +=u*u*weight; fSumXY[5]+=t*t*weight;
230 fSumXY[6] +=u*v*weight; fSumXY[7]+=u*t*weight; fSumXY[8]+=v*t*weight;
236 Double_t r = TMath::Sqrt(x*x+y*y);
237 fSumXZ[1] +=weight*r; fSumXZ[2] +=weight*r*r; fSumXZ[3] +=weight*z; fSumXZ[4] += weight*r*z;
238 // Now the auxulary sums
239 fSumXZ[5] +=weight*r*r*r/24; fSumXZ[6] +=weight*r*r*r*r/12; fSumXZ[7] +=weight*r*r*r*z/24;
240 fSumXZ[8] +=weight*r*r*r*r*r*r/(24*24);
241 fSumZZ += z*z*weight;
243 // sum accumulation for rough error estimates of the track extrapolation error
245 Double_t maty = 1./(sy*sy);
246 Double_t matz = 1./(sz*sz);
247 for (Int_t i=0;i<5; i++){
257 void AliRieman::UpdatePol(){
263 for (Int_t i=0;i<6;i++)fParams[i]=0;
268 TMatrixDSym smatrix(3);
271 // smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2;
272 // smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut;
273 // sums(0,0) = sv; sums(0,1)=suv; sums(0,2)=svt;
275 smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5];
276 smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7];
277 smatrix(1,0) = fSumXY[1]; smatrix(2,0)=fSumXY[3]; smatrix(2,1)=fSumXY[7];
278 sums(0,0) = fSumXY[2]; sums(0,1) =fSumXY[6]; sums(0,2) =fSumXY[8];
280 if (smatrix.IsValid()){
281 for (Int_t i=0;i<3;i++)
282 for (Int_t j=0;j<=i;j++){
283 (*fCovar)(i,j)=smatrix(i,j);
285 TMatrixD res = sums*smatrix;
286 fParams[0] = res(0,0);
287 fParams[1] = res(0,1);
288 fParams[2] = res(0,2);
294 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
296 // fSumXZ[1] += fSumXZ[5]*rm1*rm1;
297 // fSumXZ[2] += fSumXZ[6]*rm1*rm1 + fSumXZ[8]*rm1*rm1*rm1*rm1;
298 // fSumXZ[4] += fSumXZ[7]*rm1*rm1;
299 Double_t sum1 = fSumXZ[1] + fSumXZ[5]*rm1*rm1;
300 Double_t sum2 = fSumXZ[2] + fSumXZ[6]*rm1*rm1 + fSumXZ[8]*rm1*rm1*rm1*rm1;
301 Double_t sum4 = fSumXZ[4] + fSumXZ[7]*rm1*rm1;
303 TMatrixDSym smatrixz(2);
304 // smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = fSumXZ[1]; smatrixz(1,1) = fSumXZ[2];
305 //smatrixz(1,0) = fSumXZ[1];
306 smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = sum1; smatrixz(1,1) = sum2;
307 smatrixz(1,0) = sum1;
309 TMatrixD sumsxz(1,2);
310 if (smatrixz.IsValid()){
311 sumsxz(0,0) = fSumXZ[3];
312 // sumsxz(0,1) = fSumXZ[4];
314 TMatrixD res = sumsxz*smatrixz;
315 fParams[3] = res(0,0);
316 fParams[4] = res(0,1);
318 for (Int_t i=0;i<2;i++)
319 for (Int_t j=0;j<=i;j++){
320 (*fCovar)(i+3,j+3)=smatrixz(i,j);
324 UpdateCovariancePol();
325 // (x-x0)^2+(y-y0)^2-R^2=0 ===>
327 // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==>
328 // substitution t = 1/(x^2+y^2), u = 2*x*t, y = 2*y*t, D0 = R^2 - x0^2- y0^2
330 // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation
332 // next substition a = 1/y0 b = -x0/y0 c = -D0/y0
333 // final linear equation : a + u*b +t*c - v =0;
336 // fParam[1] = -x0/y0
337 // fParam[2] = - (R^2 - x0^2 - y0^2)/y0
338 if (conv>1) fConv =kTRUE;
343 void AliRieman::Update(){
349 for (Int_t i=0;i<6;i++)fParams[i]=0;
354 TMatrixDSym smatrix(3);
357 // smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2;
358 // smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut;
359 // sums(0,0) = sv; sums(0,1)=suv; sums(0,2)=svt;
361 smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5];
362 smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7];
364 smatrix(1,0) = fSumXY[1];
365 smatrix(2,0) = fSumXY[3];
366 smatrix(2,1) = fSumXY[7];
368 sums(0,0) = fSumXY[2]; sums(0,1) =fSumXY[6]; sums(0,2) =fSumXY[8];
369 //TDecompChol decomp(smatrix);
370 //decomp.SetTol(1.0e-14);
374 if (smatrix.IsValid()){
375 for (Int_t i=0;i<3;i++)
376 for (Int_t j=0;j<3;j++){
377 (*fCovar)(i,j)=smatrix(i,j);
379 TMatrixD res = sums*smatrix;
380 fParams[0] = res(0,0);
381 fParams[1] = res(0,1);
382 fParams[2] = res(0,2);
386 fConv = kFALSE; //non converged
389 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.<0){
390 fConv = kFALSE; //non converged
396 Double_t x0 = -fParams[1]/fParams[0];
397 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.);
400 for (Int_t i=0;i<9;i++) sumXZ[i]=0;
401 for (Int_t i=0;i<fN;i++){
402 Double_t phi = TMath::ASin((fX[i]-x0)*rm1);
403 Double_t phi0 = TMath::ASin((0.-x0)*rm1);
404 Double_t weight = 1/fSz[i];
406 Double_t dphi = (phi-phi0)/rm1;
408 sumXZ[1] +=weight*dphi;
409 sumXZ[2] +=weight*dphi*dphi;
410 sumXZ[3] +=weight*fZ[i];
411 sumXZ[4] +=weight*dphi*fZ[i];
415 TMatrixDSym smatrixz(2);
417 smatrixz(0,0) = sumXZ[0]; smatrixz(0,1) = sumXZ[1]; smatrixz(1,1) = sumXZ[2];
418 smatrixz(1,0) = sumXZ[1];
420 if (smatrixz.IsValid()){
421 sumsz(0,0) = sumXZ[3];
422 sumsz(0,1) = sumXZ[4];
423 TMatrixD res = sumsz*smatrixz;
424 fParams[3] = res(0,0);
425 fParams[4] = res(0,1);
426 for (Int_t i=0;i<2;i++)
427 for (Int_t j=0;j<2;j++){
428 (*fCovar)(i+3,j+3)=smatrixz(i,j);
432 UpdateCovariancePol();
433 // (x-x0)^2+(y-y0)^2-R^2=0 ===>
435 // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==>
436 // substitution t = 1/(x^2+y^2), u = 2*x*t, v = 2*y*t, D0 = R^2 - x0^2- y0^2
438 // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation
440 // next substition a = 1/y0 b = -x0/y0 c = -D0/y0
441 // final linear equation : a + u*b +t*c - v =0;
444 // fParam[1] = -x0/y0
445 // fParam[2] = - (R^2 - x0^2 - y0^2)/y0
446 if (conv>1) fConv =kTRUE;
449 fChi2Y = CalcChi2Y();
450 fChi2Z = CalcChi2Z();
451 fChi2 = fChi2Y +fChi2Z;
454 void AliRieman::UpdateCovariancePol(){
456 // covariance matrices for rough extrapolation error estimate
457 // covariance matrix - get error estimates at point x = 0
458 // ! NOTE - error estimates is very rough and it is valid only if dl<<R
459 // when dl is the distance between first and last point
462 (*fCovarPolY)(0,0) = fSumPolY[0]; (*fCovarPolY)(0,1) = fSumPolY[1]; (*fCovarPolY)(0,2) = fSumPolY[2];
463 (*fCovarPolY)(1,0) = fSumPolY[1]; (*fCovarPolY)(1,1) = fSumPolY[2]; (*fCovarPolY)(1,2) = fSumPolY[3];
464 (*fCovarPolY)(2,0) = fSumPolY[2]; (*fCovarPolY)(2,1) = fSumPolY[3]; (*fCovarPolY)(2,2) = fSumPolY[4];
465 fCovarPolY->Invert();
468 (*fCovarPolZ)(0,0) = fSumPolZ[0]; (*fCovarPolZ)(0,1) = fSumPolZ[1];
469 (*fCovarPolZ)(1,0) = fSumPolZ[1]; (*fCovarPolZ)(1,1) = fSumPolZ[2];
470 fCovarPolZ->Invert();
474 Double_t AliRieman::GetErrY(Double_t x) const{
476 // P0' = P0 + P1 * x + P2 * x^2
487 TMatrixD covarX(trans,TMatrixD::kMult,*fCovarPolY);
492 Double_t AliRieman::GetErrZ(Double_t x) const{
494 // assumption error of curvature determination neggligible
503 TMatrixD covarX(trans,TMatrixD::kMult,*fCovarPolZ);
508 Bool_t AliRieman::GetExternalParameters(Double_t xref, Double_t *params, Double_t * covar){
510 // Get external parameters
511 // + approximative covariance
512 // 0 1 2 3 4 // y - local y
513 // 5 6 7 8 // z - local z
514 // 9 10 11 // snp - local sin fi
515 // 12 13 // tgl - deep angle
516 // 14 // C - curvature
517 Double_t sign = (GetC()>0) ? 1.:-1.;
519 params[0] = GetYat(xref);
520 params[1] = GetZat(xref);
521 params[2] = TMath::Sin(TMath::ATan(GetDYat(xref)));
522 params[3] = sign*fParams[4];
525 // covariance matrix in y
527 TMatrixD transY(3,3);
530 transY(0,2) = xref*xref;
534 TMatrixD covarY(transY,TMatrixD::kMult,*fCovarPolY);
537 TMatrixD transZ(2,2);
541 TMatrixD covarZ(transZ,TMatrixD::kMult,*fCovarPolZ);
544 // C ~ 2*P2 - in rotated frame
546 covar[0] = covarY(0,0); covar[1] = 0; covar[2] = covarY(0,1); covar[3] = 0; covar[4] = TMath::Sqrt(2.)*covarY(0,2);
547 covar[5] = covarZ(0,0); covar[6] = 0; covar[7] = sign*covarZ(0,1); covar[8] = 0;
548 covar[9] = covarY(1,1); covar[10]= 0; covar[11]= TMath::Sqrt(2.)*covarY(1,2);
549 covar[12] = covarZ(1,1); covar[13]= 0;
550 covar[14] = 2.*covarY(2,2);
558 Double_t AliRieman::GetYat(Double_t x) const {
560 // get y at given x position
562 if (!fConv) return 0.;
563 Double_t res2 = (x*fParams[0]+fParams[1]);
565 res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
566 if (res2<0) return 0;
567 res2 = TMath::Sqrt(res2);
568 res2 = (1-res2)/fParams[0];
572 Double_t AliRieman::GetDYat(Double_t x) const {
574 // get dy/dx at given x position
576 if (!fConv) return 0.;
577 Double_t x0 = -fParams[1]/fParams[0];
578 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<0) return 0;
579 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
580 if ( 1./(rm1*rm1)-(x-x0)*(x-x0)<=0) return 0;
581 Double_t res = (x-x0)/TMath::Sqrt(1./(rm1*rm1)-(x-x0)*(x-x0));
582 if (fParams[0]<0) res*=-1.;
588 Double_t AliRieman::GetZat(Double_t x) const {
590 // get z at given x position
592 if (!fConv) return 0.;
593 Double_t x0 = -fParams[1]/fParams[0];
594 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
595 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
596 Double_t phi = TMath::ASin((x-x0)*rm1);
597 Double_t phi0 = TMath::ASin((0.-x0)*rm1);
598 Double_t dphi = (phi-phi0);
599 Double_t res = fParams[3]+fParams[4]*dphi/rm1;
603 Double_t AliRieman::GetDZat(Double_t x) const {
605 // get dz/dx at given x postion
607 if (!fConv) return 0.;
608 Double_t x0 = -fParams[1]/fParams[0];
609 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
610 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
611 Double_t res = fParams[4]/TMath::Cos(TMath::ASin((x-x0)*rm1));
616 //_____________________________________________________________________________
617 Bool_t AliRieman::GetXYZat(Double_t r, Double_t alpha, Float_t *xyz) const
620 // Returns position given radius
622 if (!fConv) return kFALSE;
623 Double_t res2 = (r*fParams[0]+fParams[1]);
625 res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
626 if (res2<0) return kFALSE;
627 res2 = TMath::Sqrt(res2);
628 res2 = (1-res2)/fParams[0];
630 if (!fConv) return kFALSE;
631 Double_t x0 = -fParams[1]/fParams[0];
632 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
633 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
634 Double_t phi = TMath::ASin((r-x0)*rm1);
635 Double_t phi0 = TMath::ASin((0.-x0)*rm1);
636 Double_t dphi = (phi-phi0);
637 Double_t res = fParams[3]+fParams[4]*dphi/rm1;
639 Double_t sin = TMath::Sin(alpha);
640 Double_t cos = TMath::Cos(alpha);
641 xyz[0] = r*cos - res2*sin;
642 xyz[1] = res2*cos + r*sin;
649 Double_t AliRieman::GetC() const{
653 return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
656 Double_t AliRieman::CalcChi2Y() const{
658 // calculate chi2 for Y
660 Double_t sumchi2 = 0;
661 for (Int_t i=0;i<fN;i++){
662 Double_t chi2 = (fY[i] - GetYat(fX[i]))/fSy[i];
669 Double_t AliRieman::CalcChi2Z() const{
671 // calculate chi2 for Z
673 Double_t sumchi2 = 0;
674 for (Int_t i=0;i<fN;i++){
675 Double_t chi2 = (fZ[i] - GetZat(fX[i]))/fSy[i];
681 Double_t AliRieman::CalcChi2() const{
683 // sum chi2 in both coord - supposing Y and Z independent
685 return CalcChi2Y()+CalcChi2Z();
688 AliRieman * AliRieman::MakeResiduals() const{
690 // create residual structure - ONLY for debug purposes
692 AliRieman *rieman = new AliRieman(fN);
693 for (Int_t i=0;i<fN;i++){
694 rieman->AddPoint(fX[i],fY[i]-GetYat(fX[i]),fZ[i]-GetZat(fX[i]), fSy[i],fSz[i]);