STEER/AliTrackFitterRieman.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
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //
  20 // Class to the track points on the Riemann sphere. Inputs are
  21 // the set of id's (volids) of the volumes in which residuals are
  22 // calculated to construct a chi2 function to be minimized during
  23 // the alignment procedures. For the moment the track extrapolation is
  24 // taken at the space-point reference plane. The reference plane is
  25 // found using the covariance matrix of the point
  26 // (assuming sigma(x)=0 at the reference coordinate system.
  27 //
  28 //////////////////////////////////////////////////////////////////////////////
  29
  30 #include "TMatrixDSym.h"
  31 #include "TMatrixD.h"
  32 #include "AliTrackFitterRieman.h"
  33 #include "AliLog.h"
  34
  35 ClassImp(AliTrackFitterRieman)
  36
  37 AliTrackFitterRieman::AliTrackFitterRieman():
  38   AliTrackFitter()
  39 {
  40   //
  41   // default constructor
  42   //
  43   fAlpha = 0.;
  44   for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
  45   fSumYY = 0;
  46   for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
  47   fSumZZ = 0;
  48   fNUsed = 0;
  49   fConv = kFALSE;
  50 }
  51
  52
  53 AliTrackFitterRieman::AliTrackFitterRieman(AliTrackPointArray *array, Bool_t owner):
  54   AliTrackFitter(array,owner)
  55 {
  56   //
  57   // Constructor
  58   //
  59   fAlpha = 0.;
  60   for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
  61   fSumYY = 0;
  62   for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
  63   fSumZZ = 0;
  64   fNUsed = 0;
  65   fConv = kFALSE;
  66 }
  67
  68 AliTrackFitterRieman::AliTrackFitterRieman(const AliTrackFitterRieman &rieman):
  69   AliTrackFitter(rieman)
  70 {
  71   //
  72   // copy constructor
  73   //
  74   fAlpha = rieman.fAlpha;
  75   for (Int_t i=0;i<9;i++) fSumXY[i]  = rieman.fSumXY[i];
  76   fSumYY = rieman.fSumYY;
  77   for (Int_t i=0;i<9;i++) fSumXZ[i]  = rieman.fSumXZ[i];
  78   fSumZZ = rieman.fSumZZ;
  79   fNUsed = rieman.fNUsed;
  80   fConv = rieman.fConv;
  81 }
  82
  83 //_____________________________________________________________________________
  84 AliTrackFitterRieman &AliTrackFitterRieman::operator =(const AliTrackFitterRieman& rieman)
  85 {
  86   //
  87   // Assignment operator
  88   //
  89   if(this==&rieman) return *this;
  90   ((AliTrackFitter *)this)->operator=(rieman);
  91
  92   fAlpha = rieman.fAlpha;
  93   for (Int_t i=0;i<9;i++) fSumXY[i]  = rieman.fSumXY[i];
  94   fSumYY = rieman.fSumYY;
  95   for (Int_t i=0;i<9;i++) fSumXZ[i]  = rieman.fSumXZ[i];
  96   fSumZZ = rieman.fSumZZ;
  97   fNUsed = rieman.fNUsed;
  98   fConv = rieman.fConv;
  99
 100   return *this;
 101 }
 102
 103 //_____________________________________________________________________________
 104 void AliTrackFitterRieman::Reset()
 105 {
 106   //
 107   // Reset the track parameters and
 108   // rieman sums
 109   //
 110   AliTrackFitter::Reset();
 111   fAlpha = 0.;
 112   for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
 113   fSumYY = 0;
 114   for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
 115   fSumZZ = 0;
 116   fNUsed = 0;
 117   fConv =kFALSE;
 118 }
 119
 120 Bool_t AliTrackFitterRieman::Fit(const TArrayI *volIds,const TArrayI *volIdsFit,
 121                                  AliAlignObj::ELayerID layerRangeMin,
 122                                  AliAlignObj::ELayerID layerRangeMax)
 123 {
 124   // Fit the track points. The method takes as an input
 125   // the set of id's (volids) of the volumes in which
 126   // one wants to calculate the residuals.
 127   // The following parameters are used to define the
 128   // range of volumes to be used in the fitting
 129   // As a result two AliTrackPointArray's obects are filled.
 130   // The first one contains the space points with
 131   // volume id's from volids list. The second array of points represents
 132   // the track extrapolations corresponding to the space points
 133   // in the first array. The two arrays can be used to find
 134   // the residuals in the volids and consequently construct a
 135   // chi2 function to be minimized during the alignment
 136   // procedures. For the moment the track extrapolation is taken
 137   // at the space-point reference plane. The reference plane is
 138   // found using the covariance matrix of the point
 139   // (assuming sigma(x)=0 at the reference coordinate system.
 140
 141   Reset();
 142
 143   Int_t npoints = fPoints->GetNPoints();
 144   if (npoints < 3) return kFALSE;
 145
 146   Bool_t isAlphaCalc = kFALSE;
 147   AliTrackPoint p,plocal;
 148 //   fPoints->GetPoint(p,0);
 149 //   fAlpha = TMath::ATan2(p.GetY(),p.GetX());
 150
 151   Int_t npVolId = 0;
 152   fNUsed = 0;
 153   Int_t *pindex = new Int_t[npoints];
 154   fX  = new Float_t[npoints];
 155   fY  = new Float_t[npoints];
 156   fZ  = new Float_t[npoints];
 157   fSy = new Float_t[npoints];
 158   fSz = new Float_t[npoints];
 159   for (Int_t ipoint = 0; ipoint < npoints; ipoint++)
 160     {
 161       fPoints->GetPoint(p,ipoint);
 162       UShort_t iVolId = p.GetVolumeID();
 163       if (FindVolId(volIds,iVolId)) {
 164         pindex[npVolId] = ipoint;
 165         npVolId++;
 166       }
 167       if (volIdsFit != 0x0) {
 168         if (!FindVolId(volIdsFit,iVolId)) continue;
 169       }
 170       else {
 171         if (iVolId < AliAlignObj::LayerToVolUID(layerRangeMin,0) ||
 172             iVolId > AliAlignObj::LayerToVolUID(layerRangeMax,
 173                                                 AliAlignObj::LayerSize(layerRangeMax))) continue;
 174       }
 175       if (!isAlphaCalc) {
 176         fAlpha = p.GetAngle();
 177         isAlphaCalc = kTRUE;
 178       }
 179       plocal = p.Rotate(fAlpha);
 180       AddPoint(plocal.GetX(),plocal.GetY(),plocal.GetZ(),
 181                TMath::Sqrt(plocal.GetCov()[3]),TMath::Sqrt(plocal.GetCov()[5]));
 182       fNUsed++;
 183     }
 184
 185   if (npVolId == 0 || fNUsed < 3) {
 186     delete [] pindex;
 187     delete [] fX;
 188     delete [] fY;
 189     delete [] fZ;
 190     delete [] fSy;
 191     delete [] fSz;
 192     return kFALSE;
 193   }
 194
 195   Update();
 196
 197   delete [] fX;
 198   delete [] fY;
 199   delete [] fZ;
 200   delete [] fSy;
 201   delete [] fSz;
 202
 203   if (!fConv) {
 204     delete [] pindex;
 205     return kFALSE;
 206   }
 207
 208   if ((fParams[0] == 0) ||
 209       ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0)) {
 210     delete [] pindex;
 211     return kFALSE;
 212   }
 213
 214
 215   if (fNUsed < fMinNPoints) {
 216     delete [] pindex;
 217     return kFALSE;
 218   }
 219
 220   fPVolId = new AliTrackPointArray(npVolId);
 221   fPTrack = new AliTrackPointArray(npVolId);
 222   AliTrackPoint p2;
 223   for (Int_t ipoint = 0; ipoint < npVolId; ipoint++)
 224     {
 225       Int_t index = pindex[ipoint];
 226       fPoints->GetPoint(p,index);
 227       if (GetPCA(p,p2)) {
 228         Float_t xyz[3],xyz2[3];
 229         p.GetXYZ(xyz); p2.GetXYZ(xyz2);
 230         //      printf("residuals %f %d %d %f %f %f %f %f %f\n",fChi2,fNUsed,fConv,xyz[0],xyz[1],xyz[2],xyz2[0]-xyz[0],xyz2[1]-xyz[1],xyz2[2]-xyz[2]);
 231         fPVolId->AddPoint(ipoint,&p);
 232         fPTrack->AddPoint(ipoint,&p2);
 233       }
 234     }
 235
 236   delete [] pindex;
 237
 238   // debug info
 239 //   Float_t chi2 = 0, chi22 = 0;
 240 //   for (Int_t ipoint = 0; ipoint < npoints; ipoint++)
 241 //     {
 242 //       fPoints->GetPoint(p,ipoint);
 243 //       UShort_t iVolId = p.GetVolumeID();
 244 //       if (volIdFit != 0) {
 245 //      if (iVolId != volIdFit) continue;
 246 //       }
 247 //       else {
 248 //      if (iVolId < AliAlignObj::LayerToVolUID(layerRangeMin,0) ||
 249 //          iVolId > AliAlignObj::LayerToVolUID(layerRangeMax,AliAlignObj::LayerSize(layerRangeMax))) continue;
 250 //       }
 251 //       plocal = p.Rotate(fAlpha);
 252 //       Float_t delta = (fParams[0]*(plocal.GetX()*plocal.GetX()+plocal.GetY()*plocal.GetY())+
 253 //                     2.*plocal.GetX()*fParams[1]+
 254 //                     fParams[2]-
 255 //                     2.*plocal.GetY())/
 256 //                    (2.*TMath::Sqrt(plocal.GetCov()[3]));
 257 // //       Float_t delta2 = (fParams[3]+
 258 // //                  plocal.GetX()*fParams[4]+
 259 // //                  plocal.GetX()*plocal.GetX()*fParams[5]-
 260 // //                  plocal.GetZ())/
 261 // //                 (TMath::Sqrt(plocal.GetCov()[5]));
 262 //       Double_t r = TMath::Sqrt(plocal.GetX()*plocal.GetX()+plocal.GetY()*plocal.GetY());
 263 //       Double_t Rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 264 //       Float_t delta2 = (fParams[3]+
 265 //                      r*fParams[4]+r*r*r*fParams[4]*Rm1*Rm1/24-
 266 //                      plocal.GetZ())/
 267 //                     (TMath::Sqrt(plocal.GetCov()[5]));
 268 //       chi2 += delta*delta;
 269 //       chi22 += delta2*delta2;
 270 //       //      printf("pulls %d %d %f %f\n",ipoint,iVolId,delta,delta2);
 271
 272 //     }
 273 //   printf("My chi2 = %f + %f = %f\n",chi2,chi22,chi2+chi22);
 274
 275   return kTRUE;
 276 }
 277
 278 void AliTrackFitterRieman::AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy, Float_t sz)
 279 {
 280   //
 281   //  Rieman update
 282   //
 283   //------------------------------------------------------
 284   // XY direction
 285   //
 286   //  (x-x0)^2+(y-y0)^2-R^2=0 ===>
 287   //
 288   //  (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0;  ==>
 289   //
 290   //   substitution t = 1/(x^2+y^2),   u = 2*x*t, v = 2*y*t,  D0 = R^2 - x0^2- y0^2
 291   //
 292   //  1 - u*x0 - v*y0 - t *D0 =0 ;  - linear equation
 293   //
 294   //  next substition   a = 1/y0    b = -x0/y0   c = -D0/y0
 295   //
 296   //  final linear equation :   a + u*b +t*c - v =0;
 297   //
 298   // Minimization :
 299   //
 300   // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min;
 301   //
 302   // where sigmai is the error of  maesurement  (a + ui*b +ti*c - vi)
 303   //
 304   // neglecting error of xi, and supposing  xi>>yi    sigmai ~ sigmaVi ~ 2*sigmay*t
 305   //
 306   fX[fNUsed] = x; fY[fNUsed]=y; fZ[fNUsed]=z; fSy[fNUsed]=sy; fSz[fNUsed]=sz;
 307   //
 308   // XY part
 309   //
 310   Double_t  t  =  x*x+y*y;
 311   if (t<2) return;
 312   t            = 1./t;
 313   Double_t  u  =  2.*x*t;
 314   Double_t  v  =  2.*y*t;
 315   Double_t  error = 2.*sy*t;
 316   error *=error;
 317   Double_t weight = 1./error;
 318   fSumXY[0] +=weight;
 319   fSumXY[1] +=u*weight;      fSumXY[2]+=v*weight;  fSumXY[3]+=t*weight;
 320   fSumXY[4] +=u*u*weight;    fSumXY[5]+=t*t*weight;
 321   fSumXY[6] +=u*v*weight;    fSumXY[7]+=u*t*weight; fSumXY[8]+=v*t*weight;
 322   fSumYY += v*v*weight;
 323   //
 324   // XZ part
 325   //
 326   if (1) {
 327     weight = 1./(sz*sz);
 328 //     fSumXZ[0] +=weight;
 329 //     fSumXZ[1] +=weight*x;   fSumXZ[2] +=weight*x*x; fSumXZ[3] +=weight*x*x*x; fSumXZ[4] += weight*x*x*x*x;
 330 //     fSumXZ[5] +=weight*z;   fSumXZ[6] +=weight*x*z; fSumXZ[7] +=weight*x*x*z;
 331     fSumZZ += z*z*weight;
 332   }
 333   else {
 334     weight = 1./(sz*sz);
 335     fSumXZ[0] +=weight;
 336     Double_t r = TMath::Sqrt(x*x+y*y);
 337     fSumXZ[1] +=weight*r;   fSumXZ[2] +=weight*r*r; fSumXZ[3] +=weight*z; fSumXZ[4] += weight*r*z;
 338     // Now the auxulary sums
 339     fSumXZ[5] +=weight*r*r*r/24; fSumXZ[6] +=weight*r*r*r*r/12; fSumXZ[7] +=weight*r*r*r*z/24;
 340     fSumXZ[8] +=weight*r*r*r*r*r*r/(24*24);
 341     fSumZZ += z*z*weight;
 342   }
 343 }
 344
 345 void AliTrackFitterRieman::Update(){
 346   //
 347   //  Rieman update
 348   //
 349   //
 350   for (Int_t i=0;i<6;i++)fParams[i]=0;
 351   fChi2 = 0;
 352   fNdf = 0;
 353   Int_t conv=0;
 354   //
 355   // XY part
 356   //
 357   TMatrixDSym     smatrix(3);
 358   TMatrixD        sums(1,3);
 359   //
 360   //   smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2;
 361   //   smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut;
 362   //   sums(0,0)    = sv; sums(0,1)=suv; sums(0,2)=svt;
 363
 364   smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5];
 365   smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7];
 366   sums(0,0)    = fSumXY[2]; sums(0,1)   =fSumXY[6]; sums(0,2)   =fSumXY[8];
 367   smatrix.Invert();
 368   if (smatrix.IsValid()){
 369     for (Int_t i=0;i<3;i++)
 370       for (Int_t j=0;j<=i;j++){
 371         (*fCov)(i,j)=smatrix(i,j);
 372       }
 373     TMatrixD  res = sums*smatrix;
 374     fParams[0] = res(0,0);
 375     fParams[1] = res(0,1);
 376     fParams[2] = res(0,2);
 377     TMatrixD  tmp = res*sums.T();
 378     fChi2 += fSumYY - tmp(0,0);
 379     fNdf  += fNUsed - 3;
 380     conv++;
 381   }
 382   //
 383   // XZ part
 384   //
 385   if (1) {
 386     Double_t x0 = -fParams[1]/fParams[0];
 387     Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 388
 389     for (Int_t i=0;i<fNUsed;i++){
 390       Double_t phi  = TMath::ASin((fX[i]-x0)*rm1);
 391       Double_t phi0 = TMath::ASin((0.-x0)*rm1);
 392       Double_t weight = 1/fSz[i];
 393       weight *=weight;
 394       Double_t dphi = (phi-phi0)/rm1;
 395       fSumXZ[0] +=weight;
 396       fSumXZ[1] +=weight*dphi;
 397       fSumXZ[2] +=weight*dphi*dphi;
 398       fSumXZ[3] +=weight*fZ[i];
 399       fSumXZ[4] +=weight*dphi*fZ[i];
 400     }
 401
 402     TMatrixDSym     smatrixz(2);
 403     smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = fSumXZ[1]; smatrixz(1,1) = fSumXZ[2];
 404     smatrixz.Invert();
 405     TMatrixD        sumsxz(1,2);
 406     if (smatrixz.IsValid()){
 407       sumsxz(0,0)    = fSumXZ[3];
 408       sumsxz(0,1)    = fSumXZ[4];
 409       TMatrixD res = sumsxz*smatrixz;
 410       fParams[3] = res(0,0);
 411       fParams[4] = res(0,1);
 412       fParams[5] = 0;
 413       for (Int_t i=0;i<2;i++)
 414         for (Int_t j=0;j<=i;j++){
 415           (*fCov)(i+3,j+3)=smatrixz(i,j);
 416         }
 417       TMatrixD  tmp = res*sumsxz.T();
 418       fChi2 += fSumZZ - tmp(0,0);
 419       fNdf  += fNUsed - 2;
 420       conv++;
 421     }
 422   }
 423   else {
 424     Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 425     fSumXZ[1] += fSumXZ[5]*rm1*rm1;
 426     fSumXZ[2] += fSumXZ[6]*rm1*rm1 + fSumXZ[8]*rm1*rm1*rm1*rm1;
 427     fSumXZ[4] += fSumXZ[7]*rm1*rm1;
 428
 429     TMatrixDSym     smatrixz(2);
 430     smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = fSumXZ[1]; smatrixz(1,1) = fSumXZ[2];
 431     smatrixz.Invert();
 432     TMatrixD        sumsxz(1,2);
 433     if (smatrixz.IsValid()){
 434       sumsxz(0,0)    = fSumXZ[3];
 435       sumsxz(0,1)    = fSumXZ[4];
 436       TMatrixD res = sumsxz*smatrixz;
 437       fParams[3] = res(0,0);
 438       fParams[4] = res(0,1);
 439       fParams[5] = 0;
 440       for (Int_t i=0;i<2;i++)
 441         for (Int_t j=0;j<=i;j++){
 442           (*fCov)(i+3,j+3)=smatrixz(i,j);
 443         }
 444       TMatrixD  tmp = res*sumsxz.T();
 445       fChi2 += fSumZZ - tmp(0,0);
 446       fNdf  += fNUsed - 2;
 447       conv++;
 448     }
 449   }
 450
 451   //  (x-x0)^2+(y-y0)^2-R^2=0 ===>
 452   //
 453   //  (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0;  ==>
 454   //   substitution t = 1/(x^2+y^2),   u = 2*x*t, y = 2*y*t,  D0 = R^2 - x0^2- y0^2
 455   //
 456   //  1 - u*x0 - v*y0 - t *D0 =0 ;  - linear equation
 457   //
 458   //  next substition   a = 1/y0    b = -x0/y0   c = -D0/y0
 459   //  final linear equation :   a + u*b +t*c - v =0;
 460   //
 461   //  fParam[0]  = 1/y0
 462   //  fParam[1]  = -x0/y0
 463   //  fParam[2]  = - (R^2 - x0^2 - y0^2)/y0
 464   if (conv>1) fConv =kTRUE;
 465   else
 466     fConv=kFALSE;
 467 }
 468
 469 //_____________________________________________________________________________
 470 Double_t AliTrackFitterRieman::GetYat(Double_t x) const
 471 {
 472   //
 473   // Returns value of Y at given X
 474   //
 475   if (!fConv) return 0.;
 476   Double_t res2 = (x*fParams[0]+fParams[1]);
 477   res2*=res2;
 478   res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
 479   if (res2<0) return 0;
 480   res2 = TMath::Sqrt(res2);
 481   res2 = (1-res2)/fParams[0];
 482   return res2;
 483 }
 484
 485 //_____________________________________________________________________________
 486 Double_t AliTrackFitterRieman::GetDYat(Double_t x) const
 487 {
 488   //
 489   // Returns value of dY/dX at given X
 490   //
 491   if (!fConv) return 0.;
 492   Double_t x0 = -fParams[1]/fParams[0];
 493   if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<0) return 0;
 494   Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 495   if ( 1./(rm1*rm1)-(x-x0)*(x-x0)<=0) return 0;
 496   Double_t res = (x-x0)/TMath::Sqrt(1./(rm1*rm1)-(x-x0)*(x-x0));
 497   if (fParams[0]<0) res*=-1.;
 498   return res;
 499 }
 500
 501 //_____________________________________________________________________________
 502 Double_t AliTrackFitterRieman::GetZat(Double_t x) const
 503 {
 504   //
 505   // Returns value of Z given X
 506   //
 507   if (!fConv) return 0.;
 508   Double_t x0 = -fParams[1]/fParams[0];
 509   if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
 510   Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 511   Double_t phi  = TMath::ASin((x-x0)*rm1);
 512   Double_t phi0 = TMath::ASin((0.-x0)*rm1);
 513   Double_t dphi = (phi-phi0);
 514   Double_t res = fParams[3]+fParams[4]*dphi/rm1;
 515   return res;
 516 }
 517
 518 //_____________________________________________________________________________
 519 Double_t AliTrackFitterRieman::GetDZat(Double_t x) const
 520 {
 521   //
 522   // Returns value of dZ/dX given X
 523   //
 524   if (!fConv) return 0.;
 525   Double_t x0 = -fParams[1]/fParams[0];
 526   if  (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
 527   Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 528   Double_t res = fParams[4]/TMath::Cos(TMath::ASin((x-x0)*rm1));
 529   return res;
 530 }
 531
 532
 533 //_____________________________________________________________________________
 534 Double_t AliTrackFitterRieman::GetC() const
 535 {
 536   //
 537   // Returns curvature
 538   //
 539   return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 540 }
 541
 542 //_____________________________________________________________________________
 543 Bool_t AliTrackFitterRieman::GetXYZat(Double_t r, Float_t *xyz) const
 544 {
 545   //
 546   // Returns position given radius
 547   //
 548   if (!fConv) return kFALSE;
 549   Double_t res2 = (r*fParams[0]+fParams[1]);
 550   res2*=res2;
 551   res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
 552   if (res2<0) return kFALSE;
 553   res2 = TMath::Sqrt(res2);
 554   res2 = (1-res2)/fParams[0];
 555
 556   if (!fConv) return kFALSE;
 557   Double_t x0 = -fParams[1]/fParams[0];
 558   if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
 559   Double_t rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 560   Double_t phi  = TMath::ASin((r-x0)*rm1);
 561   Double_t phi0 = TMath::ASin((0.-x0)*rm1);
 562   Double_t dphi = (phi-phi0);
 563   Double_t res = fParams[3]+fParams[4]*dphi/rm1;
 564
 565   Double_t sin = TMath::Sin(fAlpha);
 566   Double_t cos = TMath::Cos(fAlpha);
 567   xyz[0] = r*cos - res2*sin;
 568   xyz[1] = res2*cos + r*sin;
 569   xyz[2] = res;
 570
 571   return kTRUE;
 572 }
 573
 574 //_____________________________________________________________________________
 575 Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) const
 576 {
 577   //
 578   // Get the closest to a given spacepoint track trajectory point
 579   // Look for details in the description of the Fit() method
 580   //
 581   if (!fConv) return kFALSE;
 582
 583   // First X and Y coordinates
 584   Double_t sin = TMath::Sin(fAlpha);
 585   Double_t cos = TMath::Cos(fAlpha);
 586   //  fParam[0]  = 1/y0
 587   //  fParam[1]  = -x0/y0
 588   //  fParam[2]  = - (R^2 - x0^2 - y0^2)/y0
 589   if (fParams[0] == 0) return kFALSE;
 590   // Track parameters in the global coordinate system
 591   Double_t x0 = -fParams[1]/fParams[0]*cos -         1./fParams[0]*sin;
 592   Double_t y0 =          1./fParams[0]*cos - fParams[1]/fParams[0]*sin;
 593   if ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0) return kFALSE;
 594   Double_t r  = TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1)/
 595                 fParams[0];
 596
 597   // Define space-point refence plane
 598   Double_t alphap = p.GetAngle();
 599   Double_t sinp = TMath::Sin(alphap);
 600   Double_t cosp = TMath::Cos(alphap);
 601   Double_t x  = p.GetX()*cosp + p.GetY()*sinp;
 602   Double_t y  = p.GetY()*cosp - p.GetX()*sinp;
 603   Double_t x0p= x0*cosp + y0*sinp;
 604   Double_t y0p= y0*cosp - x0*sinp;
 605   if ((r*r - (x-x0p)*(x-x0p))<0) {
 606     AliWarning(Form("Track extrapolation failed ! (Track radius = %f, track circle x = %f, space-point x = %f, reference plane angle = %f\n",r,x0p,x,alphap));
 607     return kFALSE;
 608   }
 609   Double_t temp = TMath::Sqrt(r*r - (x-x0p)*(x-x0p));
 610   Double_t y1 = y0p + temp;
 611   Double_t y2 = y0p - temp;
 612   Double_t yprime = y1;
 613   if(TMath::Abs(y2-y) < TMath::Abs(y1-y)) yprime = y2;
 614
 615   // Back to the global coordinate system
 616   Double_t xsecond = x*cosp - yprime*sinp;
 617   Double_t ysecond = yprime*cosp + x*sinp;
 618
 619   // Now Z coordinate and track angles
 620   Double_t x2 = xsecond*cos + ysecond*sin;
 621   Double_t zsecond = GetZat(x2);
 622   Double_t dydx = GetDYat(x2);
 623   Double_t dzdx = GetDZat(x2);
 624
 625   // Fill the cov matrix of the track extrapolation point
 626   Double_t cov[6] = {0,0,0,0,0,0};
 627   Double_t sigmax = 100*100.;
 628   cov[0] = sigmax;           cov[1] = sigmax*dydx;      cov[2] = sigmax*dzdx;
 629   cov[3] = sigmax*dydx*dydx; cov[4] = sigmax*dydx*dzdx;
 630   cov[5] = sigmax*dzdx*dzdx;
 631
 632   Float_t  newcov[6];
 633   newcov[0] = cov[0]*cos*cos-
 634             2*cov[1]*sin*cos+
 635               cov[3]*sin*sin;
 636   newcov[1] = cov[1]*(cos*cos-sin*sin)-
 637              (cov[3]-cov[0])*sin*cos;
 638   newcov[2] = cov[2]*cos-
 639               cov[4]*sin;
 640   newcov[3] = cov[0]*sin*sin+
 641             2*cov[1]*sin*cos+
 642               cov[3]*cos*cos;
 643   newcov[4] = cov[4]*cos+
 644               cov[2]*sin;
 645   newcov[5] = cov[5];
 646
 647   p2.SetXYZ(xsecond,ysecond,zsecond,newcov);
 648
 649   return kTRUE;
 650 }