STEER/AliRieman.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   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  **************************************************************************/
  15
  16 /* $Id$ */
  17 // Class for global helix fit of a track
  18 // Author: M.Ivanov
  19 // The method uses the following idea:
  20 //------------------------------------------------------
  21 // XY direction
  22 //
  23 //  (x-x0)^2+(y-y0)^2-R^2=0 ===>
  24 //
  25 //  (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0;  ==>
  26 //
  27 //   substitution t = 1/(x^2+y^2),   u = 2*x*t, v = 2*y*t,  D0 = R^2 - x0^2- y0^2
  28 //
  29 //  1 - u*x0 - v*y0 - t *D0 =0 ;  - linear equation
  30 //
  31 //  next substition   a = 1/y0    b = -x0/y0   c = -D0/y0
  32 //
  33 //  final linear equation :   a + u*b +t*c - v =0;
  34 //
  35 // Minimization :
  36 //
  37 // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min;
  38 //
  39 // where sigmai is the error of  maesurement  (a + ui*b +ti*c - vi)
  40 //
  41 // neglecting error of xi, and supposing  xi>>yi    sigmai ~ sigmaVi ~ 2*sigmay*t
  42
  43
  44 #include "TMatrixDSym.h"
  45 //#include "TDecompChol.h"
  46 #include "TMatrixD.h"
  47 #include "AliRieman.h"
  48
  49 ClassImp(AliRieman)
  50
  51
  52
  53 AliRieman::AliRieman() :
  54   TObject(),
  55   fCapacity(0),
  56   fN(0),
  57   fX(0),
  58   fY(0),
  59   fZ(0),
  60   fSy(0),
  61   fSz(0),
  62   fCovar(0),
  63   fCovarPolY(0),
  64   fCovarPolZ(0),
  65   fSumZZ(0),
  66   fChi2(0),
  67   fChi2Y(0),
  68   fChi2Z(0),
  69   fConv(kFALSE)
  70 {
  71   //
  72   // default constructor
  73   //
  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++) {
  78     fSumPolY[i]=0;
  79     fSumPolZ[i]=0;
  80   }
  81 }
  82
  83
  84 AliRieman::AliRieman(Int_t capacity) :
  85   TObject(),
  86   fCapacity(capacity),
  87   fN(0),
  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)),
  96   fSumZZ(0),
  97   fChi2(0),
  98   fChi2Y(0),
  99   fChi2Z(0),
 100   fConv(kFALSE)
 101 {
 102   //
 103   // default constructor
 104   //
 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++) {
 109     fSumPolY[i]=0;
 110     fSumPolZ[i]=0;
 111   }
 112 }
 113
 114 AliRieman::AliRieman(const AliRieman &rieman):
 115   TObject(rieman),
 116   fCapacity(rieman.fN),
 117   fN(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),
 127   fChi2(rieman.fChi2),
 128   fChi2Y(rieman.fChi2Y),
 129   fChi2Z(rieman.fChi2Z),
 130   fConv(rieman.fConv)
 131
 132 {
 133   //
 134   // copy constructor
 135   //
 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];
 142   }
 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];
 149   }
 150 }
 151
 152 AliRieman::~AliRieman()
 153 {
 154   //
 155   // Destructor
 156   //
 157   delete[]fX;
 158   delete[]fY;
 159   delete[]fZ;
 160   delete[]fSy;
 161   delete[]fSz;
 162   delete fCovar;
 163   delete fCovarPolY;
 164   delete fCovarPolZ;
 165 }
 166
 167 void AliRieman::Reset()
 168 {
 169   //
 170   // Reset all the data members
 171   //
 172   fN=0;
 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++) {
 177     fSumPolY[i]=0;
 178     fSumPolZ[i]=0;
 179   }
 180   fSumZZ =0;
 181   fConv =kFALSE;
 182 }
 183
 184
 185 void AliRieman::AddPoint(Double_t x, Double_t y, Double_t z, Double_t sy, Double_t sz)
 186 {
 187   //
 188   //  Rieman update
 189   //
 190   //------------------------------------------------------
 191   // XY direction
 192   //
 193   //  (x-x0)^2+(y-y0)^2-R^2=0 ===>
 194   //
 195   //  (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0;  ==>
 196   //
 197   //   substitution t = 1/(x^2+y^2),   u = 2*x*t, v = 2*y*t,  D0 = R^2 - x0^2- y0^2
 198   //
 199   //  1 - u*x0 - v*y0 - t *D0 =0 ;  - linear equation
 200   //
 201   //  next substition   a = 1/y0    b = -x0/y0   c = -D0/y0
 202   //
 203   //  final linear equation :   a + u*b +t*c - v =0;
 204   //
 205   // Minimization :
 206   //
 207   // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min;
 208   //
 209   // where sigmai is the error of  maesurement  (a + ui*b +ti*c - vi)
 210   //
 211   // neglecting error of xi, and supposing  xi>>yi    sigmai ~ sigmaVi ~ 2*sigmay*t
 212   //
 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;
 215   //
 216   // XY part
 217   //
 218   Double_t  t  =  x*x+y*y;
 219   if (t<2) return;
 220   t            = 1./t;
 221   Double_t  u  =  2.*x*t;
 222   Double_t  v  =  2.*y*t;
 223   Double_t  error = 2.*sy*t;
 224   error *=error;
 225   Double_t weight = 1./error;
 226   fSumXY[0] +=weight;
 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;
 230   //
 231   // XZ part
 232   //
 233   weight = 1./(sz*sz);
 234   fSumXZ[0] +=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;
 241   //
 242   // sum accumulation for rough error estimates of the track extrapolation error
 243   //
 244   Double_t maty = 1./(sy*sy);
 245   Double_t matz = 1./(sz*sz);
 246   for (Int_t i=0;i<5; i++){
 247     fSumPolY[i] += maty;
 248     fSumPolZ[i] += matz;
 249     maty *= x;
 250     matz *= x;
 251   }
 252   fN++;
 253 }
 254
 255
 256 void AliRieman::UpdatePol(){
 257   //
 258   //  Rieman update
 259   //
 260   //
 261   if (fN==0) return;
 262   for (Int_t i=0;i<6;i++)fParams[i]=0;
 263   Int_t conv=0;
 264   //
 265   // XY part
 266   //
 267   TMatrixDSym     smatrix(3);
 268   TMatrixD        sums(1,3);
 269   //
 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;
 273
 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];
 278   smatrix.Invert();
 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);
 283       }
 284     TMatrixD  res = sums*smatrix;
 285     fParams[0] = res(0,0);
 286     fParams[1] = res(0,1);
 287     fParams[2] = res(0,2);
 288     conv++;
 289   }
 290   //
 291   // XZ part
 292   //
 293   Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 294
 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;
 301
 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;
 307   smatrixz.Invert();
 308   TMatrixD        sumsxz(1,2);
 309   if (smatrixz.IsValid()){
 310     sumsxz(0,0)    = fSumXZ[3];
 311     //    sumsxz(0,1)    = fSumXZ[4];
 312     sumsxz(0,1)    = sum4;
 313     TMatrixD res = sumsxz*smatrixz;
 314     fParams[3] = res(0,0);
 315     fParams[4] = res(0,1);
 316     fParams[5] = 0;
 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);
 320       }
 321     conv++;
 322   }
 323   UpdateCovariancePol();
 324   //  (x-x0)^2+(y-y0)^2-R^2=0 ===>
 325   //
 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
 328   //
 329   //  1 - u*x0 - v*y0 - t *D0 =0 ;  - linear equation
 330   //
 331   //  next substition   a = 1/y0    b = -x0/y0   c = -D0/y0
 332   //  final linear equation :   a + u*b +t*c - v =0;
 333   //
 334   //  fParam[0]  = 1/y0
 335   //  fParam[1]  = -x0/y0
 336   //  fParam[2]  = - (R^2 - x0^2 - y0^2)/y0
 337   if (conv>1) fConv =kTRUE;
 338   else
 339     fConv=kFALSE;
 340 }
 341
 342 void AliRieman::Update(){
 343   //
 344   //  Rieman update
 345   //
 346   //
 347   if (fN==0) return;
 348   for (Int_t i=0;i<6;i++)fParams[i]=0;
 349   Int_t conv=0;
 350   //
 351   // XY part
 352   //
 353   TMatrixDSym     smatrix(3);
 354   TMatrixD        sums(1,3);
 355   //
 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;
 359
 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];
 362   //
 363   smatrix(1,0) = fSumXY[1];
 364   smatrix(2,0) = fSumXY[3];
 365   smatrix(2,1) = fSumXY[7];
 366
 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);
 370   //smatrix =
 371   //decomp.Invert();
 372   smatrix.Invert();
 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);
 377       }
 378     TMatrixD  res = sums*smatrix;
 379     fParams[0] = res(0,0);
 380     fParams[1] = res(0,1);
 381     fParams[2] = res(0,2);
 382     conv++;
 383   }
 384   if (conv==0){
 385     fConv = kFALSE;   //non converged
 386     return;
 387   }
 388   if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.<0){
 389     fConv = kFALSE;   //non converged
 390     return;
 391   }
 392   //
 393   // XZ part
 394   //
 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.);
 397   Double_t sumXZ[9];
 398
 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];
 404     weight *=weight;
 405     Double_t dphi = (phi-phi0)/rm1;
 406     sumXZ[0] +=weight;
 407     sumXZ[1] +=weight*dphi;
 408     sumXZ[2] +=weight*dphi*dphi;
 409     sumXZ[3] +=weight*fZ[i];
 410     sumXZ[4] +=weight*dphi*fZ[i];
 411
 412   }
 413
 414   TMatrixDSym     smatrixz(2);
 415   TMatrixD        sumsz(1,2);
 416   smatrixz(0,0) = sumXZ[0]; smatrixz(0,1) = sumXZ[1]; smatrixz(1,1) = sumXZ[2];
 417   smatrixz(1,0) = sumXZ[1];
 418   smatrixz.Invert();
 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);
 428     }
 429     conv++;
 430   }
 431   UpdateCovariancePol();
 432   //  (x-x0)^2+(y-y0)^2-R^2=0 ===>
 433   //
 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
 436   //
 437   //  1 - u*x0 - v*y0 - t *D0 =0 ;  - linear equation
 438   //
 439   //  next substition   a = 1/y0    b = -x0/y0   c = -D0/y0
 440   //  final linear equation :   a + u*b +t*c - v =0;
 441   //
 442   //  fParam[0]  = 1/y0
 443   //  fParam[1]  = -x0/y0
 444   //  fParam[2]  = - (R^2 - x0^2 - y0^2)/y0
 445   if (conv>1) fConv =kTRUE;
 446   else
 447     fConv=kFALSE;
 448   fChi2Y = CalcChi2Y();
 449   fChi2Z = CalcChi2Z();
 450   fChi2  = fChi2Y +fChi2Z;
 451 }
 452
 453 void AliRieman::UpdateCovariancePol(){
 454   //
 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
 459   //
 460   //
 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();
 465   //
 466
 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();
 470   //
 471 }
 472
 473 Double_t  AliRieman::GetErrY(Double_t x) const{
 474   //
 475   //    P0'  = P0 + P1 * x +  P2 * x^2
 476   //    P1'  =      P1     +  P2 * x
 477   //    P2'  =             +  P2
 478   TMatrixD trans(3,3);
 479   trans(0,0) = 1;
 480   trans(0,1) = x;
 481   trans(0,2) = x*x;
 482   trans(1,1) = 1;
 483   trans(1,2) = x;
 484   trans(2,2) = 1;
 485   //
 486   TMatrixD covarX(trans,TMatrixD::kMult,*fCovarPolY);
 487   covarX*=trans.T();
 488   return covarX(0,0);
 489 }
 490
 491 Double_t  AliRieman::GetErrZ(Double_t x) const{
 492   //
 493   //    assumption error of curvature determination neggligible
 494   //
 495   //    P0'  = P0 + P1 * x
 496   //    P1'  =      P1
 497   TMatrixD trans(2,2);
 498   trans(0,0) = 1;
 499   trans(0,1) = x;
 500   trans(1,1) = 1;
 501   //
 502   TMatrixD covarX(trans,TMatrixD::kMult,*fCovarPolZ);
 503   covarX*=trans.T();
 504   return covarX(0,0);
 505 }
 506
 507 Bool_t AliRieman::GetExternalParameters(Double_t xref, Double_t *params, Double_t * covar){
 508   //
 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.;
 517
 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];
 522   params[4] = GetC();
 523   //
 524   // covariance matrix in y
 525   //
 526   TMatrixD transY(3,3);
 527   transY(0,0) = 1;
 528   transY(0,1) = xref;
 529   transY(0,2) = xref*xref;
 530   transY(1,1) = 1;
 531   transY(1,2) = xref;
 532   transY(2,2) = 1;
 533   TMatrixD covarY(transY,TMatrixD::kMult,*fCovarPolY);
 534   covarY*=transY.T();
 535   //
 536   TMatrixD transZ(2,2);
 537   transZ(0,0) = 1;
 538   transZ(0,1) = xref;
 539   transZ(1,1) = 1;
 540   TMatrixD covarZ(transZ,TMatrixD::kMult,*fCovarPolZ);
 541   covarZ*=transZ.T();
 542   //
 543   // C ~ 2*P2 - in rotated frame
 544   //
 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);
 550   //
 551   return fConv;
 552 }
 553
 554
 555
 556
 557 Double_t AliRieman::GetYat(Double_t x) const {
 558   //
 559   // get y at given x position
 560   //
 561   if (!fConv) return 0.;
 562   Double_t res2 = (x*fParams[0]+fParams[1]);
 563   res2*=res2;
 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];
 568   return res2;
 569 }
 570
 571 Double_t AliRieman::GetDYat(Double_t x) const {
 572   //
 573   // get dy/dx at given x position
 574   //
 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.;
 582   return res;
 583 }
 584
 585
 586
 587 Double_t AliRieman::GetZat(Double_t x) const {
 588   //
 589   // get z at given x position
 590   //
 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;
 599   return res;
 600 }
 601
 602 Double_t AliRieman::GetDZat(Double_t x) const {
 603   //
 604   // get dz/dx at given x postion
 605   //
 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));
 611   return res;
 612 }
 613
 614
 615 //_____________________________________________________________________________
 616 Bool_t AliRieman::GetXYZat(Double_t r, Double_t alpha, Float_t *xyz) const
 617 {
 618   //
 619   // Returns position given radius
 620   //
 621   if (!fConv) return kFALSE;
 622   Double_t res2 = (r*fParams[0]+fParams[1]);
 623   res2*=res2;
 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];
 628
 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;
 637
 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;
 642   xyz[2] = res;
 643
 644   return kTRUE;
 645 }
 646
 647
 648 Double_t AliRieman::GetC() const{
 649   //
 650   // get curvature
 651   //
 652   return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 653 }
 654
 655 Double_t AliRieman::CalcChi2Y() const{
 656   //
 657   // calculate chi2 for Y
 658   //
 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];
 662     sumchi2+=chi2*chi2;
 663   }
 664   return sumchi2;
 665 }
 666
 667
 668 Double_t AliRieman::CalcChi2Z() const{
 669   //
 670   // calculate chi2 for Z
 671   //
 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];
 675     sumchi2+=chi2*chi2;
 676   }
 677   return sumchi2;
 678 }
 679
 680 Double_t AliRieman::CalcChi2() const{
 681   //
 682   // sum chi2 in both coord - supposing Y and Z independent
 683   //
 684   return CalcChi2Y()+CalcChi2Z();
 685 }
 686
 687 AliRieman * AliRieman::MakeResiduals() const{
 688   //
 689   // create residual structure - ONLY for debug purposes
 690   //
 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]);
 694   }
 695   return rieman;
 696 }
 697
 698