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   //
  55   // default constructor
  56   //
  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++) {
  61     fSumPolY[i]=0;
  62     fSumPolZ[i]=0;
  63   }
  64   fSumZZ = 0;
  65   fCapacity = 0;
  66   fN =0;
  67   fX  = 0;
  68   fY  = 0;
  69   fZ  = 0;
  70   fSy = 0;
  71   fSz = 0;
  72   fChi2  = 0;
  73   fChi2Y = 0;
  74   fChi2Z = 0;
  75   fCovar = 0;
  76   fCovarPolY = 0;
  77   fCovarPolZ = 0;
  78   fConv = kFALSE;
  79 }
  80
  81
  82 AliRieman::AliRieman(Int_t capacity){
  83   //
  84   // default constructor
  85   //
  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++) {
  90     fSumPolY[i]=0;
  91     fSumPolZ[i]=0;
  92   }
  93   fSumZZ =0;
  94   fCapacity = capacity;
  95   fN =0;
  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);
 104   fChi2  = 0;
 105   fChi2Y = 0;
 106   fChi2Z = 0;
 107   fConv = kFALSE;
 108 }
 109
 110 AliRieman::AliRieman(const AliRieman &rieman):TObject(rieman){
 111   //
 112   // copy constructor
 113   //
 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];
 120   }
 121   fSumZZ    = rieman.fSumZZ;
 122   fCapacity = rieman.fN;
 123   fN =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];
 139   }
 140 }
 141
 142 AliRieman::~AliRieman()
 143 {
 144   //
 145   // Destructor
 146   //
 147   delete[]fX;
 148   delete[]fY;
 149   delete[]fZ;
 150   delete[]fSy;
 151   delete[]fSz;
 152   delete fCovar;
 153   delete fCovarPolY;
 154   delete fCovarPolZ;
 155 }
 156
 157 void AliRieman::Reset()
 158 {
 159   //
 160   // Reset all the data members
 161   //
 162   fN=0;
 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++) {
 167     fSumPolY[i]=0;
 168     fSumPolZ[i]=0;
 169   }
 170   fSumZZ =0;
 171   fConv =kFALSE;
 172 }
 173
 174
 175 void AliRieman::AddPoint(Double_t x, Double_t y, Double_t z, Double_t sy, Double_t sz)
 176 {
 177   //
 178   //  Rieman update
 179   //
 180   //------------------------------------------------------
 181   // XY direction
 182   //
 183   //  (x-x0)^2+(y-y0)^2-R^2=0 ===>
 184   //
 185   //  (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0;  ==>
 186   //
 187   //   substitution t = 1/(x^2+y^2),   u = 2*x*t, v = 2*y*t,  D0 = R^2 - x0^2- y0^2
 188   //
 189   //  1 - u*x0 - v*y0 - t *D0 =0 ;  - linear equation
 190   //
 191   //  next substition   a = 1/y0    b = -x0/y0   c = -D0/y0
 192   //
 193   //  final linear equation :   a + u*b +t*c - v =0;
 194   //
 195   // Minimization :
 196   //
 197   // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min;
 198   //
 199   // where sigmai is the error of  maesurement  (a + ui*b +ti*c - vi)
 200   //
 201   // neglecting error of xi, and supposing  xi>>yi    sigmai ~ sigmaVi ~ 2*sigmay*t
 202   //
 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;
 205   //
 206   // XY part
 207   //
 208   Double_t  t  =  x*x+y*y;
 209   if (t<2) return;
 210   t            = 1./t;
 211   Double_t  u  =  2.*x*t;
 212   Double_t  v  =  2.*y*t;
 213   Double_t  error = 2.*sy*t;
 214   error *=error;
 215   Double_t weight = 1./error;
 216   fSumXY[0] +=weight;
 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;
 220   //
 221   // XZ part
 222   //
 223   weight = 1./(sz*sz);
 224   fSumXZ[0] +=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;
 231   //
 232   // sum accumulation for rough error estimates of the track extrapolation error
 233   //
 234   Double_t maty = 1./(sy*sy);
 235   Double_t matz = 1./(sz*sz);
 236   for (Int_t i=0;i<5; i++){
 237     fSumPolY[i] += maty;
 238     fSumPolZ[i] += matz;
 239     maty *= x;
 240     matz *= x;
 241   }
 242   fN++;
 243 }
 244
 245
 246 void AliRieman::UpdatePol(){
 247   //
 248   //  Rieman update
 249   //
 250   //
 251   if (fN==0) return;
 252   for (Int_t i=0;i<6;i++)fParams[i]=0;
 253   Int_t conv=0;
 254   //
 255   // XY part
 256   //
 257   TMatrixDSym     smatrix(3);
 258   TMatrixD        sums(1,3);
 259   //
 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;
 263
 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];
 268   smatrix.Invert();
 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);
 273       }
 274     TMatrixD  res = sums*smatrix;
 275     fParams[0] = res(0,0);
 276     fParams[1] = res(0,1);
 277     fParams[2] = res(0,2);
 278     conv++;
 279   }
 280   //
 281   // XZ part
 282   //
 283   Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 284
 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;
 291
 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;
 297   smatrixz.Invert();
 298   TMatrixD        sumsxz(1,2);
 299   if (smatrixz.IsValid()){
 300     sumsxz(0,0)    = fSumXZ[3];
 301     //    sumsxz(0,1)    = fSumXZ[4];
 302     sumsxz(0,1)    = sum4;
 303     TMatrixD res = sumsxz*smatrixz;
 304     fParams[3] = res(0,0);
 305     fParams[4] = res(0,1);
 306     fParams[5] = 0;
 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);
 310       }
 311     conv++;
 312   }
 313   UpdateCovariancePol();
 314   //  (x-x0)^2+(y-y0)^2-R^2=0 ===>
 315   //
 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
 318   //
 319   //  1 - u*x0 - v*y0 - t *D0 =0 ;  - linear equation
 320   //
 321   //  next substition   a = 1/y0    b = -x0/y0   c = -D0/y0
 322   //  final linear equation :   a + u*b +t*c - v =0;
 323   //
 324   //  fParam[0]  = 1/y0
 325   //  fParam[1]  = -x0/y0
 326   //  fParam[2]  = - (R^2 - x0^2 - y0^2)/y0
 327   if (conv>1) fConv =kTRUE;
 328   else
 329     fConv=kFALSE;
 330 }
 331
 332 void AliRieman::Update(){
 333   //
 334   //  Rieman update
 335   //
 336   //
 337   if (fN==0) return;
 338   for (Int_t i=0;i<6;i++)fParams[i]=0;
 339   Int_t conv=0;
 340   //
 341   // XY part
 342   //
 343   TMatrixDSym     smatrix(3);
 344   TMatrixD        sums(1,3);
 345   //
 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;
 349
 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];
 352   //
 353   smatrix(1,0) = fSumXY[1];
 354   smatrix(2,0) = fSumXY[3];
 355   smatrix(2,1) = fSumXY[7];
 356
 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);
 360   //smatrix =
 361   //decomp.Invert();
 362   smatrix.Invert();
 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);
 367       }
 368     TMatrixD  res = sums*smatrix;
 369     fParams[0] = res(0,0);
 370     fParams[1] = res(0,1);
 371     fParams[2] = res(0,2);
 372     conv++;
 373   }
 374   if (conv==0){
 375     fConv = kFALSE;   //non converged
 376     return;
 377   }
 378   if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1.<0){
 379     fConv = kFALSE;   //non converged
 380     return;
 381   }
 382   //
 383   // XZ part
 384   //
 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.);
 387   Double_t sumXZ[9];
 388
 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];
 394     weight *=weight;
 395     Double_t dphi = (phi-phi0)/rm1;
 396     sumXZ[0] +=weight;
 397     sumXZ[1] +=weight*dphi;
 398     sumXZ[2] +=weight*dphi*dphi;
 399     sumXZ[3] +=weight*fZ[i];
 400     sumXZ[4] +=weight*dphi*fZ[i];
 401
 402   }
 403
 404   TMatrixDSym     smatrixz(2);
 405   TMatrixD        sumsz(1,2);
 406   smatrixz(0,0) = sumXZ[0]; smatrixz(0,1) = sumXZ[1]; smatrixz(1,1) = sumXZ[2];
 407   smatrixz(1,0) = sumXZ[1];
 408   smatrixz.Invert();
 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);
 418     }
 419     conv++;
 420   }
 421   UpdateCovariancePol();
 422   //  (x-x0)^2+(y-y0)^2-R^2=0 ===>
 423   //
 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
 426   //
 427   //  1 - u*x0 - v*y0 - t *D0 =0 ;  - linear equation
 428   //
 429   //  next substition   a = 1/y0    b = -x0/y0   c = -D0/y0
 430   //  final linear equation :   a + u*b +t*c - v =0;
 431   //
 432   //  fParam[0]  = 1/y0
 433   //  fParam[1]  = -x0/y0
 434   //  fParam[2]  = - (R^2 - x0^2 - y0^2)/y0
 435   if (conv>1) fConv =kTRUE;
 436   else
 437     fConv=kFALSE;
 438   fChi2Y = CalcChi2Y();
 439   fChi2Z = CalcChi2Z();
 440   fChi2  = fChi2Y +fChi2Z;
 441 }
 442
 443 void AliRieman::UpdateCovariancePol(){
 444   //
 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
 449   //
 450   //
 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();
 455   //
 456
 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();
 460   //
 461 }
 462
 463 Double_t  AliRieman::GetErrY(Double_t x) const{
 464   //
 465   //    P0'  = P0 + P1 * x +  P2 * x^2
 466   //    P1'  =      P1     +  P2 * x
 467   //    P2'  =             +  P2
 468   TMatrixD trans(3,3);
 469   trans(0,0) = 1;
 470   trans(0,1) = x;
 471   trans(0,2) = x*x;
 472   trans(1,1) = 1;
 473   trans(1,2) = x;
 474   trans(2,2) = 1;
 475   //
 476   TMatrixD covarX(trans,TMatrixD::kMult,*fCovarPolY);
 477   covarX*=trans.T();
 478   return covarX(0,0);
 479 }
 480
 481 Double_t  AliRieman::GetErrZ(Double_t x) const{
 482   //
 483   //    assumption error of curvature determination neggligible
 484   //
 485   //    P0'  = P0 + P1 * x
 486   //    P1'  =      P1
 487   TMatrixD trans(2,2);
 488   trans(0,0) = 1;
 489   trans(0,1) = x;
 490   trans(1,1) = 1;
 491   //
 492   TMatrixD covarX(trans,TMatrixD::kMult,*fCovarPolZ);
 493   covarX*=trans.T();
 494   return covarX(0,0);
 495 }
 496
 497 Bool_t AliRieman::GetExternalParameters(Double_t xref, Double_t *params, Double_t * covar){
 498   //
 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.;
 507
 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];
 512   params[4] = GetC();
 513   //
 514   // covariance matrix in y
 515   //
 516   TMatrixD transY(3,3);
 517   transY(0,0) = 1;
 518   transY(0,1) = xref;
 519   transY(0,2) = xref*xref;
 520   transY(1,1) = 1;
 521   transY(1,2) = xref;
 522   transY(2,2) = 1;
 523   TMatrixD covarY(transY,TMatrixD::kMult,*fCovarPolY);
 524   covarY*=transY.T();
 525   //
 526   TMatrixD transZ(2,2);
 527   transZ(0,0) = 1;
 528   transZ(0,1) = xref;
 529   transZ(1,1) = 1;
 530   TMatrixD covarZ(transZ,TMatrixD::kMult,*fCovarPolZ);
 531   covarZ*=transZ.T();
 532   //
 533   // C ~ 2*P2 - in rotated frame
 534   //
 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);
 540   //
 541   return fConv;
 542 }
 543
 544
 545
 546
 547 Double_t AliRieman::GetYat(Double_t x) const {
 548   //
 549   // get y at given x position
 550   //
 551   if (!fConv) return 0.;
 552   Double_t res2 = (x*fParams[0]+fParams[1]);
 553   res2*=res2;
 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];
 558   return res2;
 559 }
 560
 561 Double_t AliRieman::GetDYat(Double_t x) const {
 562   //
 563   // get dy/dx at given x position
 564   //
 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.;
 572   return res;
 573 }
 574
 575
 576
 577 Double_t AliRieman::GetZat(Double_t x) const {
 578   //
 579   // get z at given x position
 580   //
 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;
 589   return res;
 590 }
 591
 592 Double_t AliRieman::GetDZat(Double_t x) const {
 593   //
 594   // get dz/dx at given x postion
 595   //
 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));
 601   return res;
 602 }
 603
 604
 605 //_____________________________________________________________________________
 606 Bool_t AliRieman::GetXYZat(Double_t r, Double_t alpha, Float_t *xyz) const
 607 {
 608   //
 609   // Returns position given radius
 610   //
 611   if (!fConv) return kFALSE;
 612   Double_t res2 = (r*fParams[0]+fParams[1]);
 613   res2*=res2;
 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];
 618
 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;
 627
 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;
 632   xyz[2] = res;
 633
 634   return kTRUE;
 635 }
 636
 637
 638 Double_t AliRieman::GetC() const{
 639   //
 640   // get curvature
 641   //
 642   return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 643 }
 644
 645 Double_t AliRieman::CalcChi2Y() const{
 646   //
 647   // calculate chi2 for Y
 648   //
 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];
 652     sumchi2+=chi2*chi2;
 653   }
 654   return sumchi2;
 655 }
 656
 657
 658 Double_t AliRieman::CalcChi2Z() const{
 659   //
 660   // calculate chi2 for Z
 661   //
 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];
 665     sumchi2+=chi2*chi2;
 666   }
 667   return sumchi2;
 668 }
 669
 670 Double_t AliRieman::CalcChi2() const{
 671   //
 672   // sum chi2 in both coord - supposing Y and Z independent
 673   //
 674   return CalcChi2Y()+CalcChi2Z();
 675 }
 676
 677 AliRieman * AliRieman::MakeResiduals() const{
 678   //
 679   // create residual structure - ONLY for debug purposes
 680   //
 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]);
 684   }
 685   return rieman;
 686 }
 687
 688