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(){
55 // default constructor
57 for (Int_t i=0;i<6;i++) fParams[i] = 0;
58 for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
59 for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
60 for (Int_t i=0;i<6;i++) {
82 AliRieman::AliRieman(Int_t capacity){
84 // default constructor
86 for (Int_t i=0;i<6;i++) fParams[i] = 0;
87 for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
88 for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
89 for (Int_t i=0;i<6;i++) {
96 fX = new Double_t[fCapacity];
97 fY = new Double_t[fCapacity];
98 fZ = new Double_t[fCapacity];
99 fSy = new Double_t[fCapacity];
100 fSz = new Double_t[fCapacity];
101 fCovar = new TMatrixDSym(6);
102 fCovarPolY = new TMatrixDSym(3);
103 fCovarPolZ = new TMatrixDSym(2);
110 AliRieman::AliRieman(const AliRieman &rieman):TObject(rieman){
114 for (Int_t i=0;i<6;i++) fParams[i] = rieman.fParams[i];
115 for (Int_t i=0;i<9;i++) fSumXY[i] = rieman.fSumXY[i];
116 for (Int_t i=0;i<9;i++) fSumXZ[i] = rieman.fSumXZ[i];
117 for (Int_t i=0;i<6;i++) {
118 fSumPolY[i]=rieman.fSumPolY[i];
119 fSumPolZ[i]=rieman.fSumPolZ[i];
121 fSumZZ = rieman.fSumZZ;
122 fCapacity = rieman.fN;
124 fCovar = new TMatrixDSym(*(rieman.fCovar));
125 fCovarPolY = new TMatrixDSym(*(rieman.fCovarPolY));
126 fCovarPolZ = new TMatrixDSym(*(rieman.fCovarPolZ));
127 fConv = rieman.fConv;
128 fX = new Double_t[fN];
129 fY = new Double_t[fN];
130 fZ = new Double_t[fN];
131 fSy = new Double_t[fN];
132 fSz = new Double_t[fN];
133 for (Int_t i=0;i<rieman.fN;i++){
134 fX[i] = rieman.fX[i];
135 fY[i] = rieman.fY[i];
136 fZ[i] = rieman.fZ[i];
137 fSy[i] = rieman.fSy[i];
138 fSz[i] = rieman.fSz[i];
142 AliRieman::~AliRieman()
157 void AliRieman::Reset()
160 // Reset all the data members
163 for (Int_t i=0;i<6;i++) fParams[i] = 0;
164 for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
165 for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
166 for (Int_t i=0;i<6;i++) {
175 void AliRieman::AddPoint(Double_t x, Double_t y, Double_t z, Double_t sy, Double_t sz)
180 //------------------------------------------------------
183 // (x-x0)^2+(y-y0)^2-R^2=0 ===>
185 // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==>
187 // substitution t = 1/(x^2+y^2), u = 2*x*t, v = 2*y*t, D0 = R^2 - x0^2- y0^2
189 // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation
191 // next substition a = 1/y0 b = -x0/y0 c = -D0/y0
193 // final linear equation : a + u*b +t*c - v =0;
197 // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min;
199 // where sigmai is the error of maesurement (a + ui*b +ti*c - vi)
201 // neglecting error of xi, and supposing xi>>yi sigmai ~ sigmaVi ~ 2*sigmay*t
203 if (fN==fCapacity-1) return; // out of capacity
204 fX[fN] = x; fY[fN]=y; fZ[fN]=z; fSy[fN]=sy; fSz[fN]=sz;
208 Double_t t = x*x+y*y;
213 Double_t error = 2.*sy*t;
215 Double_t weight = 1./error;
217 fSumXY[1] +=u*weight; fSumXY[2]+=v*weight; fSumXY[3]+=t*weight;
218 fSumXY[4] +=u*u*weight; fSumXY[5]+=t*t*weight;
219 fSumXY[6] +=u*v*weight; fSumXY[7]+=u*t*weight; fSumXY[8]+=v*t*weight;
225 Double_t r = TMath::Sqrt(x*x+y*y);
226 fSumXZ[1] +=weight*r; fSumXZ[2] +=weight*r*r; fSumXZ[3] +=weight*z; fSumXZ[4] += weight*r*z;
227 // Now the auxulary sums
228 fSumXZ[5] +=weight*r*r*r/24; fSumXZ[6] +=weight*r*r*r*r/12; fSumXZ[7] +=weight*r*r*r*z/24;
229 fSumXZ[8] +=weight*r*r*r*r*r*r/(24*24);
230 fSumZZ += z*z*weight;
232 // sum accumulation for rough error estimates of the track extrapolation error
234 Double_t maty = 1./(sy*sy);
235 Double_t matz = 1./(sz*sz);
236 for (Int_t i=0;i<5; i++){
246 void AliRieman::UpdatePol(){
252 for (Int_t i=0;i<6;i++)fParams[i]=0;
257 TMatrixDSym smatrix(3);
260 // smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2;
261 // smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut;
262 // sums(0,0) = sv; sums(0,1)=suv; sums(0,2)=svt;
264 smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5];
265 smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7];
266 smatrix(1,0) = fSumXY[1]; smatrix(2,0)=fSumXY[3]; smatrix(2,1)=fSumXY[7];
267 sums(0,0) = fSumXY[2]; sums(0,1) =fSumXY[6]; sums(0,2) =fSumXY[8];
269 if (smatrix.IsValid()){
270 for (Int_t i=0;i<3;i++)
271 for (Int_t j=0;j<=i;j++){
272 (*fCovar)(i,j)=smatrix(i,j);
274 TMatrixD res = sums*smatrix;
275 fParams[0] = res(0,0);
276 fParams[1] = res(0,1);
277 fParams[2] = res(0,2);
283 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
285 // fSumXZ[1] += fSumXZ[5]*rm1*rm1;
286 // fSumXZ[2] += fSumXZ[6]*rm1*rm1 + fSumXZ[8]*rm1*rm1*rm1*rm1;
287 // fSumXZ[4] += fSumXZ[7]*rm1*rm1;
288 Double_t sum1 = fSumXZ[1] + fSumXZ[5]*rm1*rm1;
289 Double_t sum2 = fSumXZ[2] + fSumXZ[6]*rm1*rm1 + fSumXZ[8]*rm1*rm1*rm1*rm1;
290 Double_t sum4 = fSumXZ[4] + fSumXZ[7]*rm1*rm1;
292 TMatrixDSym smatrixz(2);
293 // smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = fSumXZ[1]; smatrixz(1,1) = fSumXZ[2];
294 //smatrixz(1,0) = fSumXZ[1];
295 smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = sum1; smatrixz(1,1) = sum2;
296 smatrixz(1,0) = sum1;
298 TMatrixD sumsxz(1,2);
299 if (smatrixz.IsValid()){
300 sumsxz(0,0) = fSumXZ[3];
301 // sumsxz(0,1) = fSumXZ[4];
303 TMatrixD res = sumsxz*smatrixz;
304 fParams[3] = res(0,0);
305 fParams[4] = res(0,1);
307 for (Int_t i=0;i<2;i++)
308 for (Int_t j=0;j<=i;j++){
309 (*fCovar)(i+3,j+3)=smatrixz(i,j);
313 UpdateCovariancePol();
314 // (x-x0)^2+(y-y0)^2-R^2=0 ===>
316 // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==>
317 // substitution t = 1/(x^2+y^2), u = 2*x*t, y = 2*y*t, D0 = R^2 - x0^2- y0^2
319 // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation
321 // next substition a = 1/y0 b = -x0/y0 c = -D0/y0
322 // final linear equation : a + u*b +t*c - v =0;
325 // fParam[1] = -x0/y0
326 // fParam[2] = - (R^2 - x0^2 - y0^2)/y0
327 if (conv>1) fConv =kTRUE;
332 void AliRieman::Update(){
338 for (Int_t i=0;i<6;i++)fParams[i]=0;
343 TMatrixDSym smatrix(3);
346 // smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2;
347 // smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut;
348 // sums(0,0) = sv; sums(0,1)=suv; sums(0,2)=svt;
350 smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5];
351 smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7];
353 smatrix(1,0) = fSumXY[1];
354 smatrix(2,0) = fSumXY[3];
355 smatrix(2,1) = fSumXY[7];
357 sums(0,0) = fSumXY[2]; sums(0,1) =fSumXY[6]; sums(0,2) =fSumXY[8];
358 //TDecompChol decomp(smatrix);
359 //decomp.SetTol(1.0e-14);
363 if (smatrix.IsValid()){
364 for (Int_t i=0;i<3;i++)
365 for (Int_t j=0;j<3;j++){
366 (*fCovar)(i,j)=smatrix(i,j);
368 TMatrixD res = sums*smatrix;
369 fParams[0] = res(0,0);
370 fParams[1] = res(0,1);
371 fParams[2] = res(0,2);
375 fConv = kFALSE; //non converged
378 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.<0){
379 fConv = kFALSE; //non converged
385 Double_t x0 = -fParams[1]/fParams[0];
386 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.);
389 for (Int_t i=0;i<9;i++) sumXZ[i]=0;
390 for (Int_t i=0;i<fN;i++){
391 Double_t phi = TMath::ASin((fX[i]-x0)*rm1);
392 Double_t phi0 = TMath::ASin((0.-x0)*rm1);
393 Double_t weight = 1/fSz[i];
395 Double_t dphi = (phi-phi0)/rm1;
397 sumXZ[1] +=weight*dphi;
398 sumXZ[2] +=weight*dphi*dphi;
399 sumXZ[3] +=weight*fZ[i];
400 sumXZ[4] +=weight*dphi*fZ[i];
404 TMatrixDSym smatrixz(2);
406 smatrixz(0,0) = sumXZ[0]; smatrixz(0,1) = sumXZ[1]; smatrixz(1,1) = sumXZ[2];
407 smatrixz(1,0) = sumXZ[1];
409 if (smatrixz.IsValid()){
410 sumsz(0,0) = sumXZ[3];
411 sumsz(0,1) = sumXZ[4];
412 TMatrixD res = sumsz*smatrixz;
413 fParams[3] = res(0,0);
414 fParams[4] = res(0,1);
415 for (Int_t i=0;i<2;i++)
416 for (Int_t j=0;j<2;j++){
417 (*fCovar)(i+3,j+3)=smatrixz(i,j);
421 UpdateCovariancePol();
422 // (x-x0)^2+(y-y0)^2-R^2=0 ===>
424 // (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0; ==>
425 // substitution t = 1/(x^2+y^2), u = 2*x*t, v = 2*y*t, D0 = R^2 - x0^2- y0^2
427 // 1 - u*x0 - v*y0 - t *D0 =0 ; - linear equation
429 // next substition a = 1/y0 b = -x0/y0 c = -D0/y0
430 // final linear equation : a + u*b +t*c - v =0;
433 // fParam[1] = -x0/y0
434 // fParam[2] = - (R^2 - x0^2 - y0^2)/y0
435 if (conv>1) fConv =kTRUE;
438 fChi2Y = CalcChi2Y();
439 fChi2Z = CalcChi2Z();
440 fChi2 = fChi2Y +fChi2Z;
443 void AliRieman::UpdateCovariancePol(){
445 // covariance matrices for rough extrapolation error estimate
446 // covariance matrix - get error estimates at point x = 0
447 // ! NOTE - error estimates is very rough and it is valid only if dl<<R
448 // when dl is the distance between first and last point
451 (*fCovarPolY)(0,0) = fSumPolY[0]; (*fCovarPolY)(0,1) = fSumPolY[1]; (*fCovarPolY)(0,2) = fSumPolY[2];
452 (*fCovarPolY)(1,0) = fSumPolY[1]; (*fCovarPolY)(1,1) = fSumPolY[2]; (*fCovarPolY)(1,2) = fSumPolY[3];
453 (*fCovarPolY)(2,0) = fSumPolY[2]; (*fCovarPolY)(2,1) = fSumPolY[3]; (*fCovarPolY)(2,2) = fSumPolY[4];
454 fCovarPolY->Invert();
457 (*fCovarPolZ)(0,0) = fSumPolZ[0]; (*fCovarPolZ)(0,1) = fSumPolZ[1];
458 (*fCovarPolZ)(1,0) = fSumPolZ[1]; (*fCovarPolZ)(1,1) = fSumPolZ[2];
459 fCovarPolZ->Invert();
463 Double_t AliRieman::GetErrY(Double_t x) const{
465 // P0' = P0 + P1 * x + P2 * x^2
476 TMatrixD covarX(trans,TMatrixD::kMult,*fCovarPolY);
481 Double_t AliRieman::GetErrZ(Double_t x) const{
483 // assumption error of curvature determination neggligible
492 TMatrixD covarX(trans,TMatrixD::kMult,*fCovarPolZ);
497 Bool_t AliRieman::GetExternalParameters(Double_t xref, Double_t *params, Double_t * covar){
499 // Get external parameters
500 // + approximative covariance
501 // 0 1 2 3 4 // y - local y
502 // 5 6 7 8 // z - local z
503 // 9 10 11 // snp - local sin fi
504 // 12 13 // tgl - deep angle
505 // 14 // C - curvature
506 Double_t sign = (GetC()>0) ? 1.:-1.;
508 params[0] = GetYat(xref);
509 params[1] = GetZat(xref);
510 params[2] = TMath::Sin(TMath::ATan(GetDYat(xref)));
511 params[3] = sign*fParams[4];
514 // covariance matrix in y
516 TMatrixD transY(3,3);
519 transY(0,2) = xref*xref;
523 TMatrixD covarY(transY,TMatrixD::kMult,*fCovarPolY);
526 TMatrixD transZ(2,2);
530 TMatrixD covarZ(transZ,TMatrixD::kMult,*fCovarPolZ);
533 // C ~ 2*P2 - in rotated frame
535 covar[0] = covarY(0,0); covar[1] = 0; covar[2] = covarY(0,1); covar[3] = 0; covar[4] = TMath::Sqrt(2.)*covarY(0,2);
536 covar[5] = covarZ(0,0); covar[6] = 0; covar[7] = sign*covarZ(0,1); covar[8] = 0;
537 covar[9] = covarY(1,1); covar[10]= 0; covar[11]= TMath::Sqrt(2.)*covarY(1,2);
538 covar[12] = covarZ(1,1); covar[13]= 0;
539 covar[14] = 2.*covarY(2,2);
547 Double_t AliRieman::GetYat(Double_t x) const {
549 // get y at given x position
551 if (!fConv) return 0.;
552 Double_t res2 = (x*fParams[0]+fParams[1]);
554 res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
555 if (res2<0) return 0;
556 res2 = TMath::Sqrt(res2);
557 res2 = (1-res2)/fParams[0];
561 Double_t AliRieman::GetDYat(Double_t x) const {
563 // get dy/dx at given x position
565 if (!fConv) return 0.;
566 Double_t x0 = -fParams[1]/fParams[0];
567 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<0) return 0;
568 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
569 if ( 1./(rm1*rm1)-(x-x0)*(x-x0)<=0) return 0;
570 Double_t res = (x-x0)/TMath::Sqrt(1./(rm1*rm1)-(x-x0)*(x-x0));
571 if (fParams[0]<0) res*=-1.;
577 Double_t AliRieman::GetZat(Double_t x) const {
579 // get z at given x position
581 if (!fConv) return 0.;
582 Double_t x0 = -fParams[1]/fParams[0];
583 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
584 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
585 Double_t phi = TMath::ASin((x-x0)*rm1);
586 Double_t phi0 = TMath::ASin((0.-x0)*rm1);
587 Double_t dphi = (phi-phi0);
588 Double_t res = fParams[3]+fParams[4]*dphi/rm1;
592 Double_t AliRieman::GetDZat(Double_t x) const {
594 // get dz/dx at given x postion
596 if (!fConv) return 0.;
597 Double_t x0 = -fParams[1]/fParams[0];
598 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
599 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
600 Double_t res = fParams[4]/TMath::Cos(TMath::ASin((x-x0)*rm1));
605 //_____________________________________________________________________________
606 Bool_t AliRieman::GetXYZat(Double_t r, Double_t alpha, Float_t *xyz) const
609 // Returns position given radius
611 if (!fConv) return kFALSE;
612 Double_t res2 = (r*fParams[0]+fParams[1]);
614 res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
615 if (res2<0) return kFALSE;
616 res2 = TMath::Sqrt(res2);
617 res2 = (1-res2)/fParams[0];
619 if (!fConv) return kFALSE;
620 Double_t x0 = -fParams[1]/fParams[0];
621 if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
622 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
623 Double_t phi = TMath::ASin((r-x0)*rm1);
624 Double_t phi0 = TMath::ASin((0.-x0)*rm1);
625 Double_t dphi = (phi-phi0);
626 Double_t res = fParams[3]+fParams[4]*dphi/rm1;
628 Double_t sin = TMath::Sin(alpha);
629 Double_t cos = TMath::Cos(alpha);
630 xyz[0] = r*cos - res2*sin;
631 xyz[1] = res2*cos + r*sin;
638 Double_t AliRieman::GetC() const{
642 return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
645 Double_t AliRieman::CalcChi2Y() const{
647 // calculate chi2 for Y
649 Double_t sumchi2 = 0;
650 for (Int_t i=0;i<fN;i++){
651 Double_t chi2 = (fY[i] - GetYat(fX[i]))/fSy[i];
658 Double_t AliRieman::CalcChi2Z() const{
660 // calculate chi2 for Z
662 Double_t sumchi2 = 0;
663 for (Int_t i=0;i<fN;i++){
664 Double_t chi2 = (fZ[i] - GetZat(fX[i]))/fSy[i];
670 Double_t AliRieman::CalcChi2() const{
672 // sum chi2 in both coord - supposing Y and Z independent
674 return CalcChi2Y()+CalcChi2Z();
677 AliRieman * AliRieman::MakeResiduals() const{
679 // create residual structure - ONLY for debug purposes
681 AliRieman *rieman = new AliRieman(fN);
682 for (Int_t i=0;i<fN;i++){
683 rieman->AddPoint(fX[i],fY[i]-GetYat(fX[i]),fZ[i]-GetZat(fX[i]), fSy[i],fSz[i]);