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 <TMath.h>
  46 #include <TMatrixD.h>
  47
  48 #include "AliRieman.h"
  49
  50 ClassImp(AliRieman)
  51
  52
  53
  54 AliRieman::AliRieman() :
  55   TObject(),
  56   fCapacity(0),
  57   fN(0),
  58   fX(0),
  59   fY(0),
  60   fZ(0),
  61   fSy(0),
  62   fSz(0),
  63   fCovar(0),
  64   fCovarPolY(0),
  65   fCovarPolZ(0),
  66   fSumZZ(0),
  67   fChi2(0),
  68   fChi2Y(0),
  69   fChi2Z(0),
  70   fConv(kFALSE)
  71 {
  72   //
  73   // default constructor
  74   //
  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++) {
  79     fSumPolY[i]=0;
  80     fSumPolZ[i]=0;
  81   }
  82 }
  83
  84
  85 AliRieman::AliRieman(Int_t capacity) :
  86   TObject(),
  87   fCapacity(capacity),
  88   fN(0),
  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)),
  97   fSumZZ(0),
  98   fChi2(0),
  99   fChi2Y(0),
 100   fChi2Z(0),
 101   fConv(kFALSE)
 102 {
 103   //
 104   // default constructor
 105   //
 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++) {
 110     fSumPolY[i]=0;
 111     fSumPolZ[i]=0;
 112   }
 113 }
 114
 115 AliRieman::AliRieman(const AliRieman &rieman):
 116   TObject(rieman),
 117   fCapacity(rieman.fN),
 118   fN(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),
 128   fChi2(rieman.fChi2),
 129   fChi2Y(rieman.fChi2Y),
 130   fChi2Z(rieman.fChi2Z),
 131   fConv(rieman.fConv)
 132
 133 {
 134   //
 135   // copy constructor
 136   //
 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];
 143   }
 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];
 150   }
 151 }
 152
 153 AliRieman::~AliRieman()
 154 {
 155   //
 156   // Destructor
 157   //
 158   delete[]fX;
 159   delete[]fY;
 160   delete[]fZ;
 161   delete[]fSy;
 162   delete[]fSz;
 163   delete fCovar;
 164   delete fCovarPolY;
 165   delete fCovarPolZ;
 166 }
 167
 168 void AliRieman::Reset()
 169 {
 170   //
 171   // Reset all the data members
 172   //
 173   fN=0;
 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++) {
 178     fSumPolY[i]=0;
 179     fSumPolZ[i]=0;
 180   }
 181   fSumZZ =0;
 182   fConv =kFALSE;
 183 }
 184
 185
 186 void AliRieman::AddPoint(Double_t x, Double_t y, Double_t z, Double_t sy, Double_t sz)
 187 {
 188   //
 189   //  Rieman update
 190   //
 191   //------------------------------------------------------
 192   // XY direction
 193   //
 194   //  (x-x0)^2+(y-y0)^2-R^2=0 ===>
 195   //
 196   //  (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0;  ==>
 197   //
 198   //   substitution t = 1/(x^2+y^2),   u = 2*x*t, v = 2*y*t,  D0 = R^2 - x0^2- y0^2
 199   //
 200   //  1 - u*x0 - v*y0 - t *D0 =0 ;  - linear equation
 201   //
 202   //  next substition   a = 1/y0    b = -x0/y0   c = -D0/y0
 203   //
 204   //  final linear equation :   a + u*b +t*c - v =0;
 205   //
 206   // Minimization :
 207   //
 208   // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min;
 209   //
 210   // where sigmai is the error of  maesurement  (a + ui*b +ti*c - vi)
 211   //
 212   // neglecting error of xi, and supposing  xi>>yi    sigmai ~ sigmaVi ~ 2*sigmay*t
 213   //
 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;
 216   //
 217   // XY part
 218   //
 219   Double_t  t  =  x*x+y*y;
 220   if (t<2) return;
 221   t            = 1./t;
 222   Double_t  u  =  2.*x*t;
 223   Double_t  v  =  2.*y*t;
 224   Double_t  error = 2.*sy*t;
 225   error *=error;
 226   Double_t weight = 1./error;
 227   fSumXY[0] +=weight;
 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;
 231   //
 232   // XZ part
 233   //
 234   weight = 1./(sz*sz);
 235   fSumXZ[0] +=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;
 242   //
 243   // sum accumulation for rough error estimates of the track extrapolation error
 244   //
 245   Double_t maty = 1./(sy*sy);
 246   Double_t matz = 1./(sz*sz);
 247   for (Int_t i=0;i<5; i++){
 248     fSumPolY[i] += maty;
 249     fSumPolZ[i] += matz;
 250     maty *= x;
 251     matz *= x;
 252   }
 253   fN++;
 254 }
 255
 256
 257 void AliRieman::UpdatePol(){
 258   //
 259   //  Rieman update
 260   //
 261   //
 262   if (fN==0) return;
 263   for (Int_t i=0;i<6;i++)fParams[i]=0;
 264   Int_t conv=0;
 265   //
 266   // XY part
 267   //
 268   TMatrixDSym     smatrix(3);
 269   TMatrixD        sums(1,3);
 270   //
 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;
 274
 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];
 279   smatrix.Invert();
 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);
 284       }
 285     TMatrixD  res = sums*smatrix;
 286     fParams[0] = res(0,0);
 287     fParams[1] = res(0,1);
 288     fParams[2] = res(0,2);
 289     conv++;
 290   }
 291   //
 292   // XZ part
 293   //
 294   Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 295
 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;
 302
 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;
 308   smatrixz.Invert();
 309   TMatrixD        sumsxz(1,2);
 310   if (smatrixz.IsValid()){
 311     sumsxz(0,0)    = fSumXZ[3];
 312     //    sumsxz(0,1)    = fSumXZ[4];
 313     sumsxz(0,1)    = sum4;
 314     TMatrixD res = sumsxz*smatrixz;
 315     fParams[3] = res(0,0);
 316     fParams[4] = res(0,1);
 317     fParams[5] = 0;
 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);
 321       }
 322     conv++;
 323   }
 324   UpdateCovariancePol();
 325   //  (x-x0)^2+(y-y0)^2-R^2=0 ===>
 326   //
 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
 329   //
 330   //  1 - u*x0 - v*y0 - t *D0 =0 ;  - linear equation
 331   //
 332   //  next substition   a = 1/y0    b = -x0/y0   c = -D0/y0
 333   //  final linear equation :   a + u*b +t*c - v =0;
 334   //
 335   //  fParam[0]  = 1/y0
 336   //  fParam[1]  = -x0/y0
 337   //  fParam[2]  = - (R^2 - x0^2 - y0^2)/y0
 338   if (conv>1) fConv =kTRUE;
 339   else
 340     fConv=kFALSE;
 341 }
 342
 343 void AliRieman::Update(){
 344   //
 345   //  Rieman update
 346   //
 347   //
 348   if (fN==0) return;
 349   for (Int_t i=0;i<6;i++)fParams[i]=0;
 350   Int_t conv=0;
 351   //
 352   // XY part
 353   //
 354   TMatrixDSym     smatrix(3);
 355   TMatrixD        sums(1,3);
 356   //
 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;
 360
 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];
 363   //
 364   smatrix(1,0) = fSumXY[1];
 365   smatrix(2,0) = fSumXY[3];
 366   smatrix(2,1) = fSumXY[7];
 367
 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);
 371   //smatrix =
 372   //decomp.Invert();
 373   smatrix.Invert();
 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);
 378       }
 379     TMatrixD  res = sums*smatrix;
 380     fParams[0] = res(0,0);
 381     fParams[1] = res(0,1);
 382     fParams[2] = res(0,2);
 383     conv++;
 384   }
 385   if (conv==0){
 386     fConv = kFALSE;   //non converged
 387     return;
 388   }
 389   if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.<0){
 390     fConv = kFALSE;   //non converged
 391     return;
 392   }
 393   //
 394   // XZ part
 395   //
 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.);
 398   Double_t sumXZ[9];
 399
 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];
 405     weight *=weight;
 406     Double_t dphi = (phi-phi0)/rm1;
 407     sumXZ[0] +=weight;
 408     sumXZ[1] +=weight*dphi;
 409     sumXZ[2] +=weight*dphi*dphi;
 410     sumXZ[3] +=weight*fZ[i];
 411     sumXZ[4] +=weight*dphi*fZ[i];
 412
 413   }
 414
 415   TMatrixDSym     smatrixz(2);
 416   TMatrixD        sumsz(1,2);
 417   smatrixz(0,0) = sumXZ[0]; smatrixz(0,1) = sumXZ[1]; smatrixz(1,1) = sumXZ[2];
 418   smatrixz(1,0) = sumXZ[1];
 419   smatrixz.Invert();
 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);
 429     }
 430     conv++;
 431   }
 432   UpdateCovariancePol();
 433   //  (x-x0)^2+(y-y0)^2-R^2=0 ===>
 434   //
 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
 437   //
 438   //  1 - u*x0 - v*y0 - t *D0 =0 ;  - linear equation
 439   //
 440   //  next substition   a = 1/y0    b = -x0/y0   c = -D0/y0
 441   //  final linear equation :   a + u*b +t*c - v =0;
 442   //
 443   //  fParam[0]  = 1/y0
 444   //  fParam[1]  = -x0/y0
 445   //  fParam[2]  = - (R^2 - x0^2 - y0^2)/y0
 446   if (conv>1) fConv =kTRUE;
 447   else
 448     fConv=kFALSE;
 449   fChi2Y = CalcChi2Y();
 450   fChi2Z = CalcChi2Z();
 451   fChi2  = fChi2Y +fChi2Z;
 452 }
 453
 454 void AliRieman::UpdateCovariancePol(){
 455   //
 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
 460   //
 461   //
 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();
 466   //
 467
 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();
 471   //
 472 }
 473
 474 Double_t  AliRieman::GetErrY(Double_t x) const{
 475   //
 476   //    P0'  = P0 + P1 * x +  P2 * x^2
 477   //    P1'  =      P1     +  P2 * x
 478   //    P2'  =             +  P2
 479   TMatrixD trans(3,3);
 480   trans(0,0) = 1;
 481   trans(0,1) = x;
 482   trans(0,2) = x*x;
 483   trans(1,1) = 1;
 484   trans(1,2) = x;
 485   trans(2,2) = 1;
 486   //
 487   TMatrixD covarX(trans,TMatrixD::kMult,*fCovarPolY);
 488   covarX*=trans.T();
 489   return covarX(0,0);
 490 }
 491
 492 Double_t  AliRieman::GetErrZ(Double_t x) const{
 493   //
 494   //    assumption error of curvature determination neggligible
 495   //
 496   //    P0'  = P0 + P1 * x
 497   //    P1'  =      P1
 498   TMatrixD trans(2,2);
 499   trans(0,0) = 1;
 500   trans(0,1) = x;
 501   trans(1,1) = 1;
 502   //
 503   TMatrixD covarX(trans,TMatrixD::kMult,*fCovarPolZ);
 504   covarX*=trans.T();
 505   return covarX(0,0);
 506 }
 507
 508 Bool_t AliRieman::GetExternalParameters(Double_t xref, Double_t *params, Double_t * covar){
 509   //
 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.;
 518
 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];
 523   params[4] = GetC();
 524   //
 525   // covariance matrix in y
 526   //
 527   TMatrixD transY(3,3);
 528   transY(0,0) = 1;
 529   transY(0,1) = xref;
 530   transY(0,2) = xref*xref;
 531   transY(1,1) = 1;
 532   transY(1,2) = xref;
 533   transY(2,2) = 1;
 534   TMatrixD covarY(transY,TMatrixD::kMult,*fCovarPolY);
 535   covarY*=transY.T();
 536   //
 537   TMatrixD transZ(2,2);
 538   transZ(0,0) = 1;
 539   transZ(0,1) = xref;
 540   transZ(1,1) = 1;
 541   TMatrixD covarZ(transZ,TMatrixD::kMult,*fCovarPolZ);
 542   covarZ*=transZ.T();
 543   //
 544   // C ~ 2*P2 - in rotated frame
 545   //
 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);
 551   //
 552   return fConv;
 553 }
 554
 555
 556
 557
 558 Double_t AliRieman::GetYat(Double_t x) const {
 559   //
 560   // get y at given x position
 561   //
 562   if (!fConv) return 0.;
 563   Double_t res2 = (x*fParams[0]+fParams[1]);
 564   res2*=res2;
 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];
 569   return res2;
 570 }
 571
 572 Double_t AliRieman::GetDYat(Double_t x) const {
 573   //
 574   // get dy/dx at given x position
 575   //
 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   Double_t arg = (1./rm1-(x-x0))*(1./rm1+(x-x0));
 581   if ( arg <= 0) return 0;
 582   Double_t res = (x-x0)/TMath::Sqrt(arg);
 583   if (fParams[0]<0) res*=-1.;
 584   return res;
 585 }
 586
 587
 588
 589 Double_t AliRieman::GetZat(Double_t x) const {
 590   //
 591   // get z at given x position
 592   //
 593   if (!fConv) return 0.;
 594   Double_t x0 = -fParams[1]/fParams[0];
 595   if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
 596   Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 597   Double_t phi  = TMath::ASin((x-x0)*rm1);
 598   Double_t phi0 = TMath::ASin((0.-x0)*rm1);
 599   Double_t dphi = (phi-phi0);
 600   Double_t res = fParams[3]+fParams[4]*dphi/rm1;
 601   return res;
 602 }
 603
 604 Double_t AliRieman::GetDZat(Double_t x) const {
 605   //
 606   // get dz/dx at given x postion
 607   //
 608   if (!fConv) return 0.;
 609   Double_t x0 = -fParams[1]/fParams[0];
 610   if  (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
 611   Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 612   Double_t res = fParams[4]/TMath::Cos(TMath::ASin((x-x0)*rm1));
 613   return res;
 614 }
 615
 616
 617 //_____________________________________________________________________________
 618 Bool_t AliRieman::GetXYZat(Double_t r, Double_t alpha, Float_t *xyz) const
 619 {
 620   //
 621   // Returns position given radius
 622   //
 623   if (!fConv) return kFALSE;
 624   Double_t res2 = (r*fParams[0]+fParams[1]);
 625   res2*=res2;
 626   res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
 627   if (res2<0) return kFALSE;
 628   res2 = TMath::Sqrt(res2);
 629   res2 = (1-res2)/fParams[0];
 630
 631   if (!fConv) return kFALSE;
 632   Double_t x0 = -fParams[1]/fParams[0];
 633   if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
 634   Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 635   Double_t phi  = TMath::ASin((r-x0)*rm1);
 636   Double_t phi0 = TMath::ASin((0.-x0)*rm1);
 637   Double_t dphi = (phi-phi0);
 638   Double_t res = fParams[3]+fParams[4]*dphi/rm1;
 639
 640   Double_t sin = TMath::Sin(alpha);
 641   Double_t cos = TMath::Cos(alpha);
 642   xyz[0] = r*cos - res2*sin;
 643   xyz[1] = res2*cos + r*sin;
 644   xyz[2] = res;
 645
 646   return kTRUE;
 647 }
 648
 649
 650 Double_t AliRieman::GetC() const{
 651   //
 652   // get curvature
 653   //
 654   return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 655 }
 656
 657 Double_t AliRieman::CalcChi2Y() const{
 658   //
 659   // calculate chi2 for Y
 660   //
 661   Double_t sumchi2 = 0;
 662   for (Int_t i=0;i<fN;i++){
 663     Double_t chi2 = (fY[i] - GetYat(fX[i]))/fSy[i];
 664     sumchi2+=chi2*chi2;
 665   }
 666   return sumchi2;
 667 }
 668
 669
 670 Double_t AliRieman::CalcChi2Z() const{
 671   //
 672   // calculate chi2 for Z
 673   //
 674   Double_t sumchi2 = 0;
 675   for (Int_t i=0;i<fN;i++){
 676     Double_t chi2 = (fZ[i] - GetZat(fX[i]))/fSy[i];
 677     sumchi2+=chi2*chi2;
 678   }
 679   return sumchi2;
 680 }
 681
 682 Double_t AliRieman::CalcChi2() const{
 683   //
 684   // sum chi2 in both coord - supposing Y and Z independent
 685   //
 686   return CalcChi2Y()+CalcChi2Z();
 687 }
 688
 689 AliRieman * AliRieman::MakeResiduals() const{
 690   //
 691   // create residual structure - ONLY for debug purposes
 692   //
 693   AliRieman *rieman = new AliRieman(fN);
 694   for (Int_t i=0;i<fN;i++){
 695     rieman->AddPoint(fX[i],fY[i]-GetYat(fX[i]),fZ[i]-GetZat(fX[i]), fSy[i],fSz[i]);
 696   }
 697   return rieman;
 698 }
 699
 700