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))) continue;
 147       }
 148       if (!isAlphaCalc) {
 149         fAlpha = p.GetAngle();
 150         isAlphaCalc = kTRUE;
 151       }
 152       plocal = p.Rotate(fAlpha);
 153       AddPoint(plocal.GetX(),plocal.GetY(),plocal.GetZ(),
 154                TMath::Sqrt(plocal.GetCov()[3]),TMath::Sqrt(plocal.GetCov()[5]));
 155       fNUsed++;
 156     }
 157
 158   if (npVolId == 0 || fNUsed < 3) {
 159     delete [] pindex;
 160     delete [] fX;
 161     delete [] fY;
 162     delete [] fZ;
 163     delete [] fSy;
 164     delete [] fSz;
 165     return kFALSE;
 166   }
 167
 168   Update();
 169
 170   delete [] fX;
 171   delete [] fY;
 172   delete [] fZ;
 173   delete [] fSy;
 174   delete [] fSz;
 175
 176   if (!fConv) {
 177     delete [] pindex;
 178     return kFALSE;
 179   }
 180
 181   if ((fParams[0] == 0) ||
 182       ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0)) {
 183     delete [] pindex;
 184     return kFALSE;
 185   }
 186
 187
 188   if (fNUsed < fMinNPoints) {
 189     delete [] pindex;
 190     return kFALSE;
 191   }
 192
 193   fPVolId = new AliTrackPointArray(npVolId);
 194   fPTrack = new AliTrackPointArray(npVolId);
 195   AliTrackPoint p2;
 196   for (Int_t ipoint = 0; ipoint < npVolId; ipoint++)
 197     {
 198       Int_t index = pindex[ipoint];
 199       fPoints->GetPoint(p,index);
 200       if (GetPCA(p,p2)) {
 201         Float_t xyz[3],xyz2[3];
 202         p.GetXYZ(xyz); p2.GetXYZ(xyz2);
 203         //      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]);
 204         fPVolId->AddPoint(ipoint,&p);
 205         fPTrack->AddPoint(ipoint,&p2);
 206       }
 207     }
 208
 209   delete [] pindex;
 210
 211   // debug info
 212 //   Float_t chi2 = 0, chi22 = 0;
 213 //   for (Int_t ipoint = 0; ipoint < npoints; ipoint++)
 214 //     {
 215 //       fPoints->GetPoint(p,ipoint);
 216 //       UShort_t iVolId = p.GetVolumeID();
 217 //       if (volIdFit != 0) {
 218 //      if (iVolId != volIdFit) continue;
 219 //       }
 220 //       else {
 221 //      if (iVolId < AliAlignObj::LayerToVolUID(layerRangeMin,0) ||
 222 //          iVolId > AliAlignObj::LayerToVolUID(layerRangeMax,AliAlignObj::LayerSize(layerRangeMax))) continue;
 223 //       }
 224 //       plocal = p.Rotate(fAlpha);
 225 //       Float_t delta = (fParams[0]*(plocal.GetX()*plocal.GetX()+plocal.GetY()*plocal.GetY())+
 226 //                     2.*plocal.GetX()*fParams[1]+
 227 //                     fParams[2]-
 228 //                     2.*plocal.GetY())/
 229 //                    (2.*TMath::Sqrt(plocal.GetCov()[3]));
 230 // //       Float_t delta2 = (fParams[3]+
 231 // //                  plocal.GetX()*fParams[4]+
 232 // //                  plocal.GetX()*plocal.GetX()*fParams[5]-
 233 // //                  plocal.GetZ())/
 234 // //                 (TMath::Sqrt(plocal.GetCov()[5]));
 235 //       Double_t r = TMath::Sqrt(plocal.GetX()*plocal.GetX()+plocal.GetY()*plocal.GetY());
 236 //       Double_t Rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 237 //       Float_t delta2 = (fParams[3]+
 238 //                      r*fParams[4]+r*r*r*fParams[4]*Rm1*Rm1/24-
 239 //                      plocal.GetZ())/
 240 //                     (TMath::Sqrt(plocal.GetCov()[5]));
 241 //       chi2 += delta*delta;
 242 //       chi22 += delta2*delta2;
 243 //       //      printf("pulls %d %d %f %f\n",ipoint,iVolId,delta,delta2);
 244
 245 //     }
 246 //   printf("My chi2 = %f + %f = %f\n",chi2,chi22,chi2+chi22);
 247
 248   return kTRUE;
 249 }
 250
 251 void AliTrackFitterRieman::AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy, Float_t sz)
 252 {
 253   //
 254   //  Rieman update
 255   //
 256   //------------------------------------------------------
 257   // XY direction
 258   //
 259   //  (x-x0)^2+(y-y0)^2-R^2=0 ===>
 260   //
 261   //  (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0;  ==>
 262   //
 263   //   substitution t = 1/(x^2+y^2),   u = 2*x*t, v = 2*y*t,  D0 = R^2 - x0^2- y0^2
 264   //
 265   //  1 - u*x0 - v*y0 - t *D0 =0 ;  - linear equation
 266   //
 267   //  next substition   a = 1/y0    b = -x0/y0   c = -D0/y0
 268   //
 269   //  final linear equation :   a + u*b +t*c - v =0;
 270   //
 271   // Minimization :
 272   //
 273   // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min;
 274   //
 275   // where sigmai is the error of  maesurement  (a + ui*b +ti*c - vi)
 276   //
 277   // neglecting error of xi, and supposing  xi>>yi    sigmai ~ sigmaVi ~ 2*sigmay*t
 278   //
 279   fX[fNUsed] = x; fY[fNUsed]=y; fZ[fNUsed]=z; fSy[fNUsed]=sy; fSz[fNUsed]=sz;
 280   //
 281   // XY part
 282   //
 283   Double_t  t  =  x*x+y*y;
 284   if (t<2) return;
 285   t            = 1./t;
 286   Double_t  u  =  2.*x*t;
 287   Double_t  v  =  2.*y*t;
 288   Double_t  error = 2.*sy*t;
 289   error *=error;
 290   Double_t weight = 1./error;
 291   fSumXY[0] +=weight;
 292   fSumXY[1] +=u*weight;      fSumXY[2]+=v*weight;  fSumXY[3]+=t*weight;
 293   fSumXY[4] +=u*u*weight;    fSumXY[5]+=t*t*weight;
 294   fSumXY[6] +=u*v*weight;    fSumXY[7]+=u*t*weight; fSumXY[8]+=v*t*weight;
 295   fSumYY += v*v*weight;
 296   //
 297   // XZ part
 298   //
 299   if (1) {
 300     weight = 1./(sz*sz);
 301 //     fSumXZ[0] +=weight;
 302 //     fSumXZ[1] +=weight*x;   fSumXZ[2] +=weight*x*x; fSumXZ[3] +=weight*x*x*x; fSumXZ[4] += weight*x*x*x*x;
 303 //     fSumXZ[5] +=weight*z;   fSumXZ[6] +=weight*x*z; fSumXZ[7] +=weight*x*x*z;
 304     fSumZZ += z*z*weight;
 305   }
 306   else {
 307     weight = 1./(sz*sz);
 308     fSumXZ[0] +=weight;
 309     Double_t r = TMath::Sqrt(x*x+y*y);
 310     fSumXZ[1] +=weight*r;   fSumXZ[2] +=weight*r*r; fSumXZ[3] +=weight*z; fSumXZ[4] += weight*r*z;
 311     // Now the auxulary sums
 312     fSumXZ[5] +=weight*r*r*r/24; fSumXZ[6] +=weight*r*r*r*r/12; fSumXZ[7] +=weight*r*r*r*z/24;
 313     fSumXZ[8] +=weight*r*r*r*r*r*r/(24*24);
 314     fSumZZ += z*z*weight;
 315   }
 316 }
 317
 318 void AliTrackFitterRieman::Update(){
 319   //
 320   //  Rieman update
 321   //
 322   //
 323   for (Int_t i=0;i<6;i++)fParams[i]=0;
 324   fChi2 = 0;
 325   fNdf = 0;
 326   Int_t conv=0;
 327   //
 328   // XY part
 329   //
 330   TMatrixDSym     smatrix(3);
 331   TMatrixD        sums(1,3);
 332   //
 333   //   smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2;
 334   //   smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut;
 335   //   sums(0,0)    = sv; sums(0,1)=suv; sums(0,2)=svt;
 336
 337   smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5];
 338   smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7];
 339   sums(0,0)    = fSumXY[2]; sums(0,1)   =fSumXY[6]; sums(0,2)   =fSumXY[8];
 340   smatrix.Invert();
 341   if (smatrix.IsValid()){
 342     for (Int_t i=0;i<3;i++)
 343       for (Int_t j=0;j<=i;j++){
 344         (*fCov)(i,j)=smatrix(i,j);
 345       }
 346     TMatrixD  res = sums*smatrix;
 347     fParams[0] = res(0,0);
 348     fParams[1] = res(0,1);
 349     fParams[2] = res(0,2);
 350     TMatrixD  tmp = res*sums.T();
 351     fChi2 += fSumYY - tmp(0,0);
 352     fNdf  += fNUsed - 3;
 353     conv++;
 354   }
 355   //
 356   // XZ part
 357   //
 358   if (1) {
 359     Double_t x0 = -fParams[1]/fParams[0];
 360     Double_t Rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 361
 362     for (Int_t i=0;i<fNUsed;i++){
 363       Double_t phi  = TMath::ASin((fX[i]-x0)*Rm1);
 364       Double_t phi0 = TMath::ASin((0.-x0)*Rm1);
 365       Double_t weight = 1/fSz[i];
 366       weight *=weight;
 367       Double_t dphi = (phi-phi0)/Rm1;
 368       fSumXZ[0] +=weight;
 369       fSumXZ[1] +=weight*dphi;
 370       fSumXZ[2] +=weight*dphi*dphi;
 371       fSumXZ[3] +=weight*fZ[i];
 372       fSumXZ[4] +=weight*dphi*fZ[i];
 373     }
 374
 375     TMatrixDSym     smatrixz(2);
 376     smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = fSumXZ[1]; smatrixz(1,1) = fSumXZ[2];
 377     smatrixz.Invert();
 378     TMatrixD        sumsxz(1,2);
 379     if (smatrixz.IsValid()){
 380       sumsxz(0,0)    = fSumXZ[3];
 381       sumsxz(0,1)    = fSumXZ[4];
 382       TMatrixD res = sumsxz*smatrixz;
 383       fParams[3] = res(0,0);
 384       fParams[4] = res(0,1);
 385       fParams[5] = 0;
 386       for (Int_t i=0;i<2;i++)
 387         for (Int_t j=0;j<=i;j++){
 388           (*fCov)(i+3,j+3)=smatrixz(i,j);
 389         }
 390       TMatrixD  tmp = res*sumsxz.T();
 391       fChi2 += fSumZZ - tmp(0,0);
 392       fNdf  += fNUsed - 2;
 393       conv++;
 394     }
 395   }
 396   else {
 397     Double_t Rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 398     fSumXZ[1] += fSumXZ[5]*Rm1*Rm1;
 399     fSumXZ[2] += fSumXZ[6]*Rm1*Rm1 + fSumXZ[8]*Rm1*Rm1*Rm1*Rm1;
 400     fSumXZ[4] += fSumXZ[7]*Rm1*Rm1;
 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
 424   //  (x-x0)^2+(y-y0)^2-R^2=0 ===>
 425   //
 426   //  (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0;  ==>
 427   //   substitution t = 1/(x^2+y^2),   u = 2*x*t, y = 2*y*t,  D0 = R^2 - x0^2- y0^2
 428   //
 429   //  1 - u*x0 - v*y0 - t *D0 =0 ;  - linear equation
 430   //
 431   //  next substition   a = 1/y0    b = -x0/y0   c = -D0/y0
 432   //  final linear equation :   a + u*b +t*c - v =0;
 433   //
 434   //  fParam[0]  = 1/y0
 435   //  fParam[1]  = -x0/y0
 436   //  fParam[2]  = - (R^2 - x0^2 - y0^2)/y0
 437   if (conv>1) fConv =kTRUE;
 438   else
 439     fConv=kFALSE;
 440 }
 441
 442 Double_t AliTrackFitterRieman::GetYat(Double_t x){
 443   if (!fConv) return 0.;
 444   Double_t res2 = (x*fParams[0]+fParams[1]);
 445   res2*=res2;
 446   res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
 447   if (res2<0) return 0;
 448   res2 = TMath::Sqrt(res2);
 449   res2 = (1-res2)/fParams[0];
 450   return res2;
 451 }
 452
 453 Double_t AliTrackFitterRieman::GetDYat(Double_t x) const {
 454   if (!fConv) return 0.;
 455   Double_t x0 = -fParams[1]/fParams[0];
 456   if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<0) return 0;
 457   Double_t Rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 458   if ( 1./(Rm1*Rm1)-(x-x0)*(x-x0)<=0) return 0;
 459   Double_t res = (x-x0)/TMath::Sqrt(1./(Rm1*Rm1)-(x-x0)*(x-x0));
 460   if (fParams[0]<0) res*=-1.;
 461   return res;
 462 }
 463
 464
 465
 466 Double_t AliTrackFitterRieman::GetZat(Double_t x) const {
 467   if (!fConv) return 0.;
 468   Double_t x0 = -fParams[1]/fParams[0];
 469   if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
 470   Double_t Rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 471   Double_t phi  = TMath::ASin((x-x0)*Rm1);
 472   Double_t phi0 = TMath::ASin((0.-x0)*Rm1);
 473   Double_t dphi = (phi-phi0);
 474   Double_t res = fParams[3]+fParams[4]*dphi/Rm1;
 475   return res;
 476 }
 477
 478 Double_t AliTrackFitterRieman::GetDZat(Double_t x) const {
 479   if (!fConv) return 0.;
 480   Double_t x0 = -fParams[1]/fParams[0];
 481   if  (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
 482   Double_t Rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 483   Double_t res = fParams[4]/TMath::Cos(TMath::ASin((x-x0)*Rm1));
 484   return res;
 485 }
 486
 487
 488 Double_t AliTrackFitterRieman::GetC(){
 489   return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 490 }
 491
 492 Bool_t AliTrackFitterRieman::GetXYZat(Double_t r, Float_t *xyz){
 493   if (!fConv) return kFALSE;
 494   Double_t res2 = (r*fParams[0]+fParams[1]);
 495   res2*=res2;
 496   res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
 497   if (res2<0) return kFALSE;
 498   res2 = TMath::Sqrt(res2);
 499   res2 = (1-res2)/fParams[0];
 500
 501   if (!fConv) return kFALSE;
 502   Double_t x0 = -fParams[1]/fParams[0];
 503   if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
 504   Double_t Rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
 505   Double_t phi  = TMath::ASin((r-x0)*Rm1);
 506   Double_t phi0 = TMath::ASin((0.-x0)*Rm1);
 507   Double_t dphi = (phi-phi0);
 508   Double_t res = fParams[3]+fParams[4]*dphi/Rm1;
 509
 510   Double_t sin = TMath::Sin(fAlpha);
 511   Double_t cos = TMath::Cos(fAlpha);
 512   xyz[0] = r*cos - res2*sin;
 513   xyz[1] = res2*cos + r*sin;
 514   xyz[2] = res;
 515
 516   return kTRUE;
 517 }
 518
 519 Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) const
 520 {
 521   // Get the closest to a given spacepoint track trajectory point
 522   // Look for details in the description of the Fit() method
 523
 524   if (!fConv) return kFALSE;
 525
 526   // First X and Y coordinates
 527   Double_t sin = TMath::Sin(fAlpha);
 528   Double_t cos = TMath::Cos(fAlpha);
 529   //  fParam[0]  = 1/y0
 530   //  fParam[1]  = -x0/y0
 531   //  fParam[2]  = - (R^2 - x0^2 - y0^2)/y0
 532   if (fParams[0] == 0) return kFALSE;
 533   // Track parameters in the global coordinate system
 534   Double_t x0 = -fParams[1]/fParams[0]*cos -         1./fParams[0]*sin;
 535   Double_t y0 =          1./fParams[0]*cos - fParams[1]/fParams[0]*sin;
 536   if ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0) return kFALSE;
 537   Double_t R  = TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1)/
 538                 fParams[0];
 539
 540   // Define space-point refence plane
 541   Double_t alphap = p.GetAngle();
 542   Double_t sinp = TMath::Sin(alphap);
 543   Double_t cosp = TMath::Cos(alphap);
 544   Double_t x  = p.GetX()*cosp + p.GetY()*sinp;
 545   Double_t y  = p.GetY()*cosp - p.GetX()*sinp;
 546   Double_t x0p= x0*cosp + y0*sinp;
 547   Double_t y0p= y0*cosp - x0*sinp;
 548   if ((R*R - (x-x0p)*(x-x0p))<0) {
 549     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));
 550     return kFALSE;
 551   }
 552   Double_t temp = TMath::Sqrt(R*R - (x-x0p)*(x-x0p));
 553   Double_t y1 = y0p + temp;
 554   Double_t y2 = y0p - temp;
 555   Double_t yprime = y1;
 556   if(TMath::Abs(y2-y) < TMath::Abs(y1-y)) yprime = y2;
 557
 558   // Back to the global coordinate system
 559   Double_t xsecond = x*cosp - yprime*sinp;
 560   Double_t ysecond = yprime*cosp + x*sinp;
 561
 562   // Now Z coordinate and track angles
 563   Double_t x2 = xsecond*cos + ysecond*sin;
 564   Double_t zsecond = GetZat(x2);
 565   Double_t dydx = GetDYat(x2);
 566   Double_t dzdx = GetDZat(x2);
 567
 568   // Fill the cov matrix of the track extrapolation point
 569   Double_t cov[6] = {0,0,0,0,0,0};
 570   Double_t sigmax = 100*100.;
 571   cov[0] = sigmax;           cov[1] = sigmax*dydx;      cov[2] = sigmax*dzdx;
 572   cov[3] = sigmax*dydx*dydx; cov[4] = sigmax*dydx*dzdx;
 573   cov[5] = sigmax*dzdx*dzdx;
 574
 575   Float_t  newcov[6];
 576   newcov[0] = cov[0]*cos*cos-
 577             2*cov[1]*sin*cos+
 578               cov[3]*sin*sin;
 579   newcov[1] = cov[1]*(cos*cos-sin*sin)-
 580              (cov[3]-cov[0])*sin*cos;
 581   newcov[2] = cov[2]*cos-
 582               cov[4]*sin;
 583   newcov[3] = cov[0]*sin*sin+
 584             2*cov[1]*sin*cos+
 585               cov[3]*cos*cos;
 586   newcov[4] = cov[4]*cos+
 587               cov[2]*sin;
 588   newcov[5] = cov[5];
 589
 590   p2.SetXYZ(xsecond,ysecond,zsecond,newcov);
 591
 592   return kTRUE;
 593 }