]> 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
Additional protection and warning message.
[u/mrichter/AliRoot.git] / STEER / AliTrackFitterRieman.cxx
CommitLineData
98937d93 1#include "TMatrixDSym.h"
2#include "TMatrixD.h"
3#include "AliTrackFitterRieman.h"
24d4520d 4#include "AliLog.h"
98937d93 5
6ClassImp(AliTrackFitterRieman)
7
8AliTrackFitterRieman::AliTrackFitterRieman():
9 AliTrackFitter()
10{
11 //
12 // default constructor
13 //
14 fAlpha = 0.;
15 for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
46ae650f 16 fSumYY = 0;
98937d93 17 for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
46ae650f 18 fSumZZ = 0;
19 fNUsed = 0;
98937d93 20 fConv = kFALSE;
21}
22
23
24AliTrackFitterRieman::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;
46ae650f 32 fSumYY = 0;
98937d93 33 for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
46ae650f 34 fSumZZ = 0;
35 fNUsed = 0;
98937d93 36 fConv = kFALSE;
37}
38
39AliTrackFitterRieman::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];
46ae650f 47 fSumYY = rieman.fSumYY;
98937d93 48 for (Int_t i=0;i<9;i++) fSumXZ[i] = rieman.fSumXZ[i];
46ae650f 49 fSumZZ = rieman.fSumZZ;
50 fNUsed = rieman.fNUsed;
98937d93 51 fConv = rieman.fConv;
52}
53
54//_____________________________________________________________________________
55AliTrackFitterRieman &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];
46ae650f 64 fSumYY = rieman.fSumYY;
98937d93 65 for (Int_t i=0;i<9;i++) fSumXZ[i] = rieman.fSumXZ[i];
46ae650f 66 fSumZZ = rieman.fSumZZ;
67 fNUsed = rieman.fNUsed;
98937d93 68 fConv = rieman.fConv;
69
70 return *this;
71}
72
73AliTrackFitterRieman::~AliTrackFitterRieman()
74{
75 // destructor
76 //
77}
78
79void 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;
46ae650f 86 fSumYY = 0;
98937d93 87 for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
46ae650f 88 fSumZZ = 0;
89 fNUsed = 0;
98937d93 90 fConv =kFALSE;
91}
92
cc345ce3 93Bool_t AliTrackFitterRieman::Fit(const TArrayI *volIds,const TArrayI *volIdsFit,
98937d93 94 AliAlignObj::ELayerID layerRangeMin,
95 AliAlignObj::ELayerID layerRangeMax)
96{
97 // Fit the track points. The method takes as an input
cc345ce3 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
98937d93 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
cc345ce3 104 // volume id's from volids list. The second array of points represents
105 // the track extrapolations corresponding to the space points
98937d93 106 // in the first array. The two arrays can be used to find
cc345ce3 107 // the residuals in the volids and consequently construct a
98937d93 108 // chi2 function to be minimized during the alignment
109 // procedures. For the moment the track extrapolation is taken
cc345ce3 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.
98937d93 113
46ae650f 114 Reset();
98937d93 115
116 Int_t npoints = fPoints->GetNPoints();
117 if (npoints < 3) return kFALSE;
118
cc345ce3 119 Bool_t isAlphaCalc = kFALSE;
46ae650f 120 AliTrackPoint p,plocal;
cc345ce3 121// fPoints->GetPoint(p,0);
122// fAlpha = TMath::ATan2(p.GetY(),p.GetX());
98937d93 123
124 Int_t npVolId = 0;
46ae650f 125 fNUsed = 0;
98937d93 126 Int_t *pindex = new Int_t[npoints];
cc345ce3 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];
98937d93 132 for (Int_t ipoint = 0; ipoint < npoints; ipoint++)
133 {
134 fPoints->GetPoint(p,ipoint);
135 UShort_t iVolId = p.GetVolumeID();
cc345ce3 136 if (FindVolId(volIds,iVolId)) {
98937d93 137 pindex[npVolId] = ipoint;
138 npVolId++;
139 }
cc345ce3 140 if (volIdsFit != 0x0) {
141 if (!FindVolId(volIdsFit,iVolId)) continue;
46ae650f 142 }
143 else {
144 if (iVolId < AliAlignObj::LayerToVolUID(layerRangeMin,0) ||
145 iVolId > AliAlignObj::LayerToVolUID(layerRangeMax,
146 AliAlignObj::LayerSize(layerRangeMax-
147 AliAlignObj::kFirstLayer))) continue;
148 }
cc345ce3 149 if (!isAlphaCalc) {
150 fAlpha = p.GetAngle();
151 isAlphaCalc = kTRUE;
152 }
46ae650f 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++;
98937d93 157 }
158
24d4520d 159 if (npVolId == 0 || fNUsed < 3) {
98937d93 160 delete [] pindex;
cc345ce3 161 delete [] fX;
162 delete [] fY;
163 delete [] fZ;
164 delete [] fSy;
165 delete [] fSz;
98937d93 166 return kFALSE;
167 }
168
169 Update();
170
cc345ce3 171 delete [] fX;
172 delete [] fY;
173 delete [] fZ;
174 delete [] fSy;
175 delete [] fSz;
176
98937d93 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
cc345ce3 188
189 if (fNUsed < fMinNPoints) {
190 delete [] pindex;
191 return kFALSE;
192 }
193
46ae650f 194 fPVolId = new AliTrackPointArray(npVolId);
195 fPTrack = new AliTrackPointArray(npVolId);
98937d93 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)) {
cc345ce3 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]);
46ae650f 205 fPVolId->AddPoint(ipoint,&p);
206 fPTrack->AddPoint(ipoint,&p2);
98937d93 207 }
208 }
209
210 delete [] pindex;
211
46ae650f 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
98937d93 250 return kTRUE;
251}
252
253void 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 //
46ae650f 265 // substitution t = 1/(x^2+y^2), u = 2*x*t, v = 2*y*t, D0 = R^2 - x0^2- y0^2
98937d93 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 //
cc345ce3 281 fX[fNUsed] = x; fY[fNUsed]=y; fZ[fNUsed]=z; fSy[fNUsed]=sy; fSz[fNUsed]=sz;
98937d93 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;
46ae650f 297 fSumYY += v*v*weight;
98937d93 298 //
299 // XZ part
300 //
cc345ce3 301 if (1) {
46ae650f 302 weight = 1./(sz*sz);
cc345ce3 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;
46ae650f 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 }
98937d93 318}
319
320void AliTrackFitterRieman::Update(){
321 //
322 // Rieman update
323 //
324 //
325 for (Int_t i=0;i<6;i++)fParams[i]=0;
46ae650f 326 fChi2 = 0;
327 fNdf = 0;
98937d93 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);
46ae650f 352 TMatrixD tmp = res*sums.T();
353 fChi2 += fSumYY - tmp(0,0);
354 fNdf += fNUsed - 3;
98937d93 355 conv++;
356 }
357 //
358 // XZ part
359 //
cc345ce3 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];
46ae650f 379 smatrixz.Invert();
cc345ce3 380 TMatrixD sumsxz(1,2);
46ae650f 381 if (smatrixz.IsValid()){
cc345ce3 382 sumsxz(0,0) = fSumXZ[3];
383 sumsxz(0,1) = fSumXZ[4];
46ae650f 384 TMatrixD res = sumsxz*smatrixz;
385 fParams[3] = res(0,0);
386 fParams[4] = res(0,1);
cc345ce3 387 fParams[5] = 0;
388 for (Int_t i=0;i<2;i++)
46ae650f 389 for (Int_t j=0;j<=i;j++){
cc345ce3 390 (*fCov)(i+3,j+3)=smatrixz(i,j);
46ae650f 391 }
392 TMatrixD tmp = res*sumsxz.T();
393 fChi2 += fSumZZ - tmp(0,0);
cc345ce3 394 fNdf += fNUsed - 2;
46ae650f 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++;
98937d93 423 }
98937d93 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
444Double_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
cc345ce3 455Double_t AliTrackFitterRieman::GetDYat(Double_t x) const {
98937d93 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
cc345ce3 468Double_t AliTrackFitterRieman::GetZat(Double_t x) const {
98937d93 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
cc345ce3 480Double_t AliTrackFitterRieman::GetDZat(Double_t x) const {
98937d93 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
490Double_t AliTrackFitterRieman::GetC(){
491 return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
492}
493
494Bool_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
521Bool_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;
cc345ce3 535 // Track parameters in the global coordinate system
98937d93 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
cc345ce3 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;
24d4520d 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 }
cc345ce3 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);
98937d93 593
594 return kTRUE;
595}