]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliTrackFitterRieman.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Return to the original methods from AliLoader (Marco)
[u/mrichter/AliRoot.git] / STEER / AliTrackFitterRieman.cxx
CommitLineData
6b6cba33 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
98937d93 30#include "TMatrixDSym.h"
31#include "TMatrixD.h"
32#include "AliTrackFitterRieman.h"
24d4520d 33#include "AliLog.h"
98937d93 34
35ClassImp(AliTrackFitterRieman)
36
37AliTrackFitterRieman::AliTrackFitterRieman():
38 AliTrackFitter()
39{
40 //
41 // default constructor
42 //
43 fAlpha = 0.;
44 for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
46ae650f 45 fSumYY = 0;
98937d93 46 for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
46ae650f 47 fSumZZ = 0;
48 fNUsed = 0;
98937d93 49 fConv = kFALSE;
50}
51
52
53AliTrackFitterRieman::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;
46ae650f 61 fSumYY = 0;
98937d93 62 for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
46ae650f 63 fSumZZ = 0;
64 fNUsed = 0;
98937d93 65 fConv = kFALSE;
66}
67
68AliTrackFitterRieman::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];
46ae650f 76 fSumYY = rieman.fSumYY;
98937d93 77 for (Int_t i=0;i<9;i++) fSumXZ[i] = rieman.fSumXZ[i];
46ae650f 78 fSumZZ = rieman.fSumZZ;
79 fNUsed = rieman.fNUsed;
98937d93 80 fConv = rieman.fConv;
81}
82
83//_____________________________________________________________________________
84AliTrackFitterRieman &AliTrackFitterRieman::operator =(const AliTrackFitterRieman& rieman)
85{
6b6cba33 86 //
87 // Assignment operator
98937d93 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];
46ae650f 94 fSumYY = rieman.fSumYY;
98937d93 95 for (Int_t i=0;i<9;i++) fSumXZ[i] = rieman.fSumXZ[i];
46ae650f 96 fSumZZ = rieman.fSumZZ;
97 fNUsed = rieman.fNUsed;
98937d93 98 fConv = rieman.fConv;
99
100 return *this;
101}
102
6b6cba33 103//_____________________________________________________________________________
98937d93 104void AliTrackFitterRieman::Reset()
105{
6b6cba33 106 //
98937d93 107 // Reset the track parameters and
108 // rieman sums
6b6cba33 109 //
98937d93 110 AliTrackFitter::Reset();
111 fAlpha = 0.;
112 for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
46ae650f 113 fSumYY = 0;
98937d93 114 for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
46ae650f 115 fSumZZ = 0;
116 fNUsed = 0;
98937d93 117 fConv =kFALSE;
118}
119
cc345ce3 120Bool_t AliTrackFitterRieman::Fit(const TArrayI *volIds,const TArrayI *volIdsFit,
98937d93 121 AliAlignObj::ELayerID layerRangeMin,
122 AliAlignObj::ELayerID layerRangeMax)
123{
124 // Fit the track points. The method takes as an input
cc345ce3 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
98937d93 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
cc345ce3 131 // volume id's from volids list. The second array of points represents
132 // the track extrapolations corresponding to the space points
98937d93 133 // in the first array. The two arrays can be used to find
cc345ce3 134 // the residuals in the volids and consequently construct a
98937d93 135 // chi2 function to be minimized during the alignment
136 // procedures. For the moment the track extrapolation is taken
cc345ce3 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.
98937d93 140
46ae650f 141 Reset();
98937d93 142
143 Int_t npoints = fPoints->GetNPoints();
144 if (npoints < 3) return kFALSE;
145
cc345ce3 146 Bool_t isAlphaCalc = kFALSE;
46ae650f 147 AliTrackPoint p,plocal;
cc345ce3 148// fPoints->GetPoint(p,0);
149// fAlpha = TMath::ATan2(p.GetY(),p.GetX());
98937d93 150
151 Int_t npVolId = 0;
46ae650f 152 fNUsed = 0;
98937d93 153 Int_t *pindex = new Int_t[npoints];
cc345ce3 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];
98937d93 159 for (Int_t ipoint = 0; ipoint < npoints; ipoint++)
160 {
161 fPoints->GetPoint(p,ipoint);
162 UShort_t iVolId = p.GetVolumeID();
cc345ce3 163 if (FindVolId(volIds,iVolId)) {
98937d93 164 pindex[npVolId] = ipoint;
165 npVolId++;
166 }
cc345ce3 167 if (volIdsFit != 0x0) {
168 if (!FindVolId(volIdsFit,iVolId)) continue;
46ae650f 169 }
170 else {
171 if (iVolId < AliAlignObj::LayerToVolUID(layerRangeMin,0) ||
172 iVolId > AliAlignObj::LayerToVolUID(layerRangeMax,
c041444f 173 AliAlignObj::LayerSize(layerRangeMax))) continue;
46ae650f 174 }
cc345ce3 175 if (!isAlphaCalc) {
176 fAlpha = p.GetAngle();
177 isAlphaCalc = kTRUE;
178 }
46ae650f 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++;
98937d93 183 }
184
24d4520d 185 if (npVolId == 0 || fNUsed < 3) {
98937d93 186 delete [] pindex;
cc345ce3 187 delete [] fX;
188 delete [] fY;
189 delete [] fZ;
190 delete [] fSy;
191 delete [] fSz;
98937d93 192 return kFALSE;
193 }
194
195 Update();
196
cc345ce3 197 delete [] fX;
198 delete [] fY;
199 delete [] fZ;
200 delete [] fSy;
201 delete [] fSz;
202
98937d93 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
cc345ce3 214
215 if (fNUsed < fMinNPoints) {
216 delete [] pindex;
217 return kFALSE;
218 }
219
46ae650f 220 fPVolId = new AliTrackPointArray(npVolId);
221 fPTrack = new AliTrackPointArray(npVolId);
98937d93 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)) {
cc345ce3 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]);
46ae650f 231 fPVolId->AddPoint(ipoint,&p);
232 fPTrack->AddPoint(ipoint,&p2);
98937d93 233 }
234 }
235
236 delete [] pindex;
237
46ae650f 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) ||
c041444f 249// iVolId > AliAlignObj::LayerToVolUID(layerRangeMax,AliAlignObj::LayerSize(layerRangeMax))) continue;
46ae650f 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
98937d93 275 return kTRUE;
276}
277
278void 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 //
46ae650f 290 // substitution t = 1/(x^2+y^2), u = 2*x*t, v = 2*y*t, D0 = R^2 - x0^2- y0^2
98937d93 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 //
cc345ce3 306 fX[fNUsed] = x; fY[fNUsed]=y; fZ[fNUsed]=z; fSy[fNUsed]=sy; fSz[fNUsed]=sz;
98937d93 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;
46ae650f 322 fSumYY += v*v*weight;
98937d93 323 //
324 // XZ part
325 //
cc345ce3 326 if (1) {
46ae650f 327 weight = 1./(sz*sz);
cc345ce3 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;
46ae650f 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 }
98937d93 343}
344
345void AliTrackFitterRieman::Update(){
346 //
347 // Rieman update
348 //
349 //
350 for (Int_t i=0;i<6;i++)fParams[i]=0;
46ae650f 351 fChi2 = 0;
352 fNdf = 0;
98937d93 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);
46ae650f 377 TMatrixD tmp = res*sums.T();
378 fChi2 += fSumYY - tmp(0,0);
379 fNdf += fNUsed - 3;
98937d93 380 conv++;
381 }
382 //
383 // XZ part
384 //
cc345ce3 385 if (1) {
386 Double_t x0 = -fParams[1]/fParams[0];
6b6cba33 387 Double_t rm1 = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
cc345ce3 388
389 for (Int_t i=0;i<fNUsed;i++){
6b6cba33 390 Double_t phi = TMath::ASin((fX[i]-x0)*rm1);
391 Double_t phi0 = TMath::ASin((0.-x0)*rm1);
cc345ce3 392 Double_t weight = 1/fSz[i];
393 weight *=weight;
6b6cba33 394 Double_t dphi = (phi-phi0)/rm1;
cc345ce3 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];
46ae650f 404 smatrixz.Invert();
cc345ce3 405 TMatrixD sumsxz(1,2);
46ae650f 406 if (smatrixz.IsValid()){
cc345ce3 407 sumsxz(0,0) = fSumXZ[3];
408 sumsxz(0,1) = fSumXZ[4];
46ae650f 409 TMatrixD res = sumsxz*smatrixz;
410 fParams[3] = res(0,0);
411 fParams[4] = res(0,1);
cc345ce3 412 fParams[5] = 0;
413 for (Int_t i=0;i<2;i++)
46ae650f 414 for (Int_t j=0;j<=i;j++){
cc345ce3 415 (*fCov)(i+3,j+3)=smatrixz(i,j);
46ae650f 416 }
417 TMatrixD tmp = res*sumsxz.T();
418 fChi2 += fSumZZ - tmp(0,0);
cc345ce3 419 fNdf += fNUsed - 2;
46ae650f 420 conv++;
421 }
422 }
423 else {
6b6cba33 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;
46ae650f 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++;
98937d93 448 }
98937d93 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
6b6cba33 469//_____________________________________________________________________________
470Double_t AliTrackFitterRieman::GetYat(Double_t x) const
471{
472 //
473 // Returns value of Y at given X
474 //
98937d93 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
6b6cba33 485//_____________________________________________________________________________
486Double_t AliTrackFitterRieman::GetDYat(Double_t x) const
487{
488 //
489 // Returns value of dY/dX at given X
490 //
98937d93 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;
6b6cba33 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));
98937d93 497 if (fParams[0]<0) res*=-1.;
498 return res;
499}
500
6b6cba33 501//_____________________________________________________________________________
502Double_t AliTrackFitterRieman::GetZat(Double_t x) const
503{
504 //
505 // Returns value of Z given X
506 //
98937d93 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;
6b6cba33 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);
98937d93 513 Double_t dphi = (phi-phi0);
6b6cba33 514 Double_t res = fParams[3]+fParams[4]*dphi/rm1;
98937d93 515 return res;
516}
517
6b6cba33 518//_____________________________________________________________________________
519Double_t AliTrackFitterRieman::GetDZat(Double_t x) const
520{
521 //
522 // Returns value of dZ/dX given X
523 //
98937d93 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;
6b6cba33 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));
98937d93 529 return res;
530}
531
532
6b6cba33 533//_____________________________________________________________________________
534Double_t AliTrackFitterRieman::GetC() const
535{
536 //
537 // Returns curvature
538 //
98937d93 539 return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
540}
541
6b6cba33 542//_____________________________________________________________________________
543Bool_t AliTrackFitterRieman::GetXYZat(Double_t r, Float_t *xyz) const
544{
545 //
546 // Returns position given radius
547 //
98937d93 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;
6b6cba33 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);
98937d93 562 Double_t dphi = (phi-phi0);
6b6cba33 563 Double_t res = fParams[3]+fParams[4]*dphi/rm1;
98937d93 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
6b6cba33 574//_____________________________________________________________________________
98937d93 575Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) const
576{
6b6cba33 577 //
98937d93 578 // Get the closest to a given spacepoint track trajectory point
579 // Look for details in the description of the Fit() method
6b6cba33 580 //
98937d93 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;
cc345ce3 590 // Track parameters in the global coordinate system
98937d93 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;
6b6cba33 594 Double_t r = TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1)/
98937d93 595 fParams[0];
596
cc345ce3 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;
6b6cba33 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));
24d4520d 607 return kFALSE;
608 }
6b6cba33 609 Double_t temp = TMath::Sqrt(r*r - (x-x0p)*(x-x0p));
cc345ce3 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);
98937d93 648
649 return kTRUE;
650}