STEER/AliTrackFitterRieman.cxx

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