added offline wrapper for HLT TPC CA tracker (Sergey)
[u/mrichter/AliRoot.git] / HLT / TPCLib / tracking-ca / AliHLTTPCCATrackParam.cxx
CommitLineData
d54804bf 1// $Id$
2//***************************************************************************
3// This file is property of and copyright by the ALICE HLT Project *
4// ALICE Experiment at CERN, All rights reserved. *
5// *
6// Primary Authors: Sergey Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de> *
7// Ivan Kisel <kisel@kip.uni-heidelberg.de> *
8// for The ALICE HLT Project. *
9// *
10// Permission to use, copy, modify and distribute this software and its *
11// documentation strictly for non-commercial purposes is hereby granted *
12// without fee, provided that the above copyright notice appears in all *
13// copies and that both the copyright notice and this permission notice *
14// appear in the supporting documentation. The authors make no claims *
15// about the suitability of this software for any purpose. It is *
16// provided "as is" without express or implied warranty. *
17//***************************************************************************
18
19#include "AliHLTTPCCATrackParam.h"
20#include "TMath.h"
21#include "AliExternalTrackParam.h"
22
23//ClassImp(AliHLTTPCCATrackParam)
24
25//
26// Circle in XY:
27//
28// R = 1/TMath::Abs(Kappa);
29// Xc = X - sin(Phi)/Kappa;
30// Yc = Y + cos(Phi)/Kappa;
31//
32
33
34
35void AliHLTTPCCATrackParam::ConstructXY3( const Float_t x[3], const Float_t y[3],
36 const Float_t sigmaY2[3], Float_t CosPhi0 )
37{
38 //* Construct the track in XY plane by 3 points
39
40 Float_t x0 = x[0];
41 Float_t y0 = y[0];
42 Float_t x1 = x[1] - x0;
43 Float_t y1 = y[1] - y0;
44 Float_t x2 = x[2] - x0;
45 Float_t y2 = y[2] - y0;
46
47 Float_t a1 = x1*x1 + y1*y1;
48 Float_t a2 = x2*x2 + y2*y2;
49 Float_t a = 2*(x1*y2 - y1*x2);
50 Float_t lx = a1*y2 - a2*y1;
51 Float_t ly = -a1*x2 + a2*x1;
52 Float_t l = TMath::Sqrt(lx*lx + ly*ly);
53
54 Float_t li = 1./l;
55 Float_t li2 = li*li;
56 Float_t li3 = li2*li;
57 Float_t cosPhi = ly*li;
58
59 Float_t sinPhi = -lx*li;
60 Float_t kappa = a/l;
61
62 Float_t dlx = a2 - a1; // D lx / D y0
63 Float_t dly = -a; // D ly / D y0
64 Float_t dA = 2*(x2 - x1); // D a / D y0
65 Float_t dl = (lx*dlx + ly*dly)*li;
66
67 // D sinPhi,kappa / D y0
68
69 Float_t d0[2] = { -(dlx*ly-lx*dly)*ly*li3, (dA*l-a*dl)*li2 };
70
71 // D sinPhi,kappa / D y1
72
73 dlx = -a2 + 2*y1*y2;
74 dly = -2*x2*y1;
75 dA = -2*x2;
76 dl = (lx*dlx + ly*dly)*li;
77
78 Float_t d1[2] = { -(dlx*ly-lx*dly)*ly*li3, (dA*l-a*dl)*li2 };
79
80 // D sinPhi,kappa / D y2
81
82 dlx = a1 - 2*y1*y2;
83 dly = -2*x1*y2;
84 dA = 2*x1;
85 dl = (lx*dlx + ly*dly)*li;
86
87 Float_t d2[2] = { -(dlx*ly-lx*dly)*ly*li3, (dA*l-a*dl)*li2 };
88
89 if( CosPhi0*cosPhi <0 ){
90 cosPhi = -cosPhi;
91 sinPhi = -sinPhi;
92 kappa = -kappa;
93 d0[0] = -d0[0];
94 d0[1] = -d0[1];
95 d1[0] = -d1[0];
96 d1[1] = -d1[1];
97 d2[0] = -d2[0];
98 d2[1] = -d2[1];
99 }
100
101 X() = x0;
102 Y() = y0;
103 SinPhi() = sinPhi;
104 Kappa() = kappa;
105 CosPhi() = cosPhi;
106
107 Float_t s0 = sigmaY2[0];
108 Float_t s1 = sigmaY2[1];
109 Float_t s2 = sigmaY2[2];
110
111 fC[0] = s0;
112 fC[1] = 0;
113 fC[2] = 0;
114
115 fC[3] = d0[0]*s0;
116 fC[4] = 0;
117 fC[5] = d0[0]*d0[0]*s0 + d1[0]*d1[0]*s1 + d2[0]*d2[0]*s2;
118
119 fC[6] = 0;
120 fC[7] = 0;
121 fC[8] = 0;
122 fC[9] = 0;
123
124 fC[10] = d0[1]*s0;
125 fC[11] = 0;
126 fC[12] = d0[0]*d0[1]*s0 + d1[0]*d1[1]*s1 + d2[0]*d2[1]*s2;
127 fC[13] = 0;
128 fC[14] = d0[1]*d0[1]*s0 + d1[1]*d1[1]*s1 + d2[1]*d2[1]*s2;
129}
130
131
132Float_t AliHLTTPCCATrackParam::GetS( Float_t x, Float_t y ) const
133{
134 //* Get XY path length to the given point
135
136 Float_t k = GetKappa();
137 Float_t ex = GetCosPhi();
138 Float_t ey = GetSinPhi();
139 x-= GetX();
140 y-= GetY();
141 Float_t dS = x*ex + y*ey;
142 if( TMath::Abs(k)>1.e-4 ) dS = TMath::ATan2( k*dS, 1+k*(x*ey-y*ex) )/k;
143 return dS;
144}
145
146void AliHLTTPCCATrackParam::GetDCAPoint( Float_t x, Float_t y, Float_t z,
147 Float_t &xp, Float_t &yp, Float_t &zp ) const
148{
149 //* Get the track point closest to the (x,y,z)
150
151 Float_t x0 = GetX();
152 Float_t y0 = GetY();
153 Float_t k = GetKappa();
154 Float_t ex = GetCosPhi();
155 Float_t ey = GetSinPhi();
156 Float_t dx = x - x0;
157 Float_t dy = y - y0;
158 Float_t ax = dx*k+ey;
159 Float_t ay = dy*k-ex;
160 Float_t a = sqrt( ax*ax+ay*ay );
161 xp = x0 + (dx - ey*( (dx*dx+dy*dy)*k - 2*(-dx*ey+dy*ex) )/(a+1) )/a;
162 yp = y0 + (dy + ex*( (dx*dx+dy*dy)*k - 2*(-dx*ey+dy*ex) )/(a+1) )/a;
163 Float_t s = GetS(x,y);
164 zp = GetZ() + GetDzDs()*s;
165 if( TMath::Abs(k)>1.e-2 ){
166 Float_t dZ = TMath::Abs( GetDzDs()*TMath::TwoPi()/k );
167 if( dZ>.1 ){
168 zp+= TMath::Nint((z-zp)/dZ)*dZ;
169 }
170 }
171}
172
173void AliHLTTPCCATrackParam::ConstructXYZ3( const Float_t p0[5], const Float_t p1[5],
174 const Float_t p2[5],
175 Float_t CosPhi0, Float_t t0[] )
176{
177 //* Construct the track in XYZ by 3 points
178
179 Float_t px[3] = { p0[0], p1[0], p2[0] };
180 Float_t py[3] = { p0[1], p1[1], p2[1] };
181 Float_t pz[3] = { p0[2], p1[2], p2[2] };
182 Float_t ps2y[3] = { p0[3]*p0[3], p1[3]*p1[3], p2[3]*p2[3] };
183 Float_t ps2z[3] = { p0[4]*p0[4], p1[4]*p1[4], p2[4]*p2[4] };
184
185 Float_t kold = t0 ?t0[4] :0;
186 ConstructXY3( px, py, ps2y, CosPhi0 );
187
188 Float_t pS[3] = { GetS(px[0],py[0]), GetS(px[1],py[1]), GetS(px[2],py[2]) };
189 Float_t k = Kappa();
190 if( TMath::Abs(k)>1.e-2 ){
191 Float_t dS = TMath::Abs( TMath::TwoPi()/k );
192 pS[1]+= TMath::Nint( (pS[0]-pS[1])/dS )*dS; // not more than half turn
193 pS[2]+= TMath::Nint( (pS[1]-pS[2])/dS )*dS;
194 if( t0 ){
195 Float_t dZ = TMath::Abs(t0[3]*dS);
196 if( TMath::Abs(dZ)>1. ){
197 Float_t dsDz = 1./t0[3];
198 if( kold*k<0 ) dsDz = -dsDz;
199 Float_t s0 = (pz[0]-t0[1])*dsDz;
200 Float_t s1 = (pz[1]-t0[1])*dsDz;
201 Float_t s2 = (pz[2]-t0[1])*dsDz;
202 pS[0]+= TMath::Nint( (s0-pS[0])/dS )*dS ;
203 pS[1]+= TMath::Nint( (s1-pS[1])/dS )*dS ;
204 pS[2]+= TMath::Nint( (s2-pS[2])/dS )*dS ;
205 }
206 }
207 }
208
209 Float_t s = pS[0] + pS[1] + pS[2];
210 Float_t z = pz[0] + pz[1] + pz[2];
211 Float_t sz = pS[0]*pz[0] + pS[1]*pz[1] + pS[2]*pz[2];
212 Float_t ss = pS[0]*pS[0] + pS[1]*pS[1] + pS[2]*pS[2];
213
214 Float_t a = 3*ss-s*s;
215 Z() = (z*ss-sz*s)/a; // z0
216 DzDs() = (3*sz-z*s)/a; // t = dz/ds
217
218 Float_t dz0[3] = {ss - pS[0]*s,ss - pS[1]*s,ss - pS[2]*s };
219 Float_t dt [3] = {3*pS[0] - s, 3*pS[1] - s, 3*pS[2] - s };
220
221 fC[2] = (dz0[0]*dz0[0]*ps2z[0] + dz0[1]*dz0[1]*ps2z[1] + dz0[2]*dz0[2]*ps2z[2])/a/a;
222 fC[7]= (dz0[0]*dt [0]*ps2z[0] + dz0[1]*dt [1]*ps2z[1] + dz0[2]*dt [2]*ps2z[2])/a/a;
223 fC[9]= (dt [0]*dt [0]*ps2z[0] + dt [1]*dt [1]*ps2z[1] + dt [2]*dt [2]*ps2z[2])/a/a;
224}
225
226
227Bool_t AliHLTTPCCATrackParam::TransportToX( Float_t x )
228{
229 //* Transport the track parameters to X=x
230
231 Bool_t ret = 1;
232
233 Float_t x0 = X();
234 //Float_t y0 = Y();
235 Float_t k = Kappa();
236 Float_t ex = CosPhi();
237 Float_t ey = SinPhi();
238 Float_t dx = x - x0;
239
240 Float_t ey1 = k*dx + ey;
241 Float_t ex1;
242 if( TMath::Abs(ey1)>1 ){ // no intersection -> check the border
243 ey1 = ( ey1>0 ) ?1 :-1;
244 ex1 = 0;
245 dx = ( TMath::Abs(k)>1.e-4) ? ( (ey1-ey)/k ) :0;
246
247 Float_t ddx = TMath::Abs(x0+dx - x)*k*k;
248 Float_t hx[] = {0, -k, 1+ey };
249 Float_t sx2 = hx[1]*hx[1]*fC[ 3] + hx[2]*hx[2]*fC[ 5];
250 if( ddx*ddx>3.5*3.5*sx2 ) ret = 0; // x not withing the error
251 ret = 0; // any case
252 return ret;
253 }else{
254 ex1 = TMath::Sqrt(1 - ey1*ey1);
255 if( ex<0 ) ex1 = -ex1;
256 }
257
258 Float_t dx2 = dx*dx;
259 CosPhi() = ex1;
260 Float_t ss = ey+ey1;
261 Float_t cc = ex+ex1;
262 Float_t tg = 0;
263 if( TMath::Abs(cc)>1.e-4 ) tg = ss/cc; // tan((phi1+phi)/2)
264 else ret = 0;
265 Float_t dy = dx*tg;
266 Float_t dl = dx*TMath::Sqrt(1+tg*tg);
267
268 if( cc<0 ) dl = -dl;
269 Float_t dSin = dl*k/2;
270 if( dSin > 1 ) dSin = 1;
271 if( dSin <-1 ) dSin = -1;
272 Float_t dS = ( TMath::Abs(k)>1.e-4) ? (2*TMath::ASin(dSin)/k) :dl;
273 Float_t dz = dS*DzDs();
274
275 Float_t cci = 0, exi = 0, ex1i = 0;
276 if( TMath::Abs(cc)>1.e-4 ) cci = 1./cc;
277 else ret = 0;
278 if( TMath::Abs(ex)>1.e-4 ) exi = 1./ex;
279 else ret = 0;
280 if( TMath::Abs(ex1)>1.e-4 ) ex1i = 1./ex1;
281 else ret = 0;
282
283 if( !ret ) return ret;
284 X() += dx;
285 fP[0]+= dy;
286 fP[1]+= dz;
287 fP[2] = ey1;
288 fP[3] = fP[3];
289 fP[4] = fP[4];
290
291 Float_t h2 = dx*(1+ ex*ex1 + ey*ey1 )*cci*exi*ex1i;
292 Float_t h4 = dx2*(cc + ss*ey1*ex1i )*cci*cci;
293
294 Float_t c00 = fC[0];
295 Float_t c10 = fC[1];
296 Float_t c11 = fC[2];
297 Float_t c20 = fC[3];
298 Float_t c21 = fC[4];
299 Float_t c22 = fC[5];
300 Float_t c30 = fC[6];
301 Float_t c31 = fC[7];
302 Float_t c32 = fC[8];
303 Float_t c33 = fC[9];
304 Float_t c40 = fC[10];
305 Float_t c41 = fC[11];
306 Float_t c42 = fC[12];
307 Float_t c43 = fC[13];
308 Float_t c44 = fC[14];
309
310 //Float_t H0[5] = { 1,0, h2, 0, h4 };
311 //Float_t H1[5] = { 0, 1, 0, dS, 0 };
312 //Float_t H2[5] = { 0, 0, 1, 0, dx };
313 //Float_t H3[5] = { 0, 0, 0, 1, 0 };
314 //Float_t H4[5] = { 0, 0, 0, 0, 1 };
315
316
317 fC[0]=( c00 + h2*h2*c22 + h4*h4*c44
318 + 2*( h2*c20 + h4*c40 + h2*h4*c42 ) );
319
320 fC[1]= c10 + h2*c21 + h4*c41 + dS*(c30 + h2*c32 + h4*c43);
321 fC[2]= c11 + 2*dS*c31 + dS*dS*c33;
322
323 fC[3]= c20 + h2*c22 + h4*c42 + dx*( c40 + h2*c42 + h4*c44);
324 fC[4]= c21 + dS*c32 + dx*(c41 + dS*c43);
325 fC[5]= c22 +2*dx*c42 + dx2*c44;
326
327 fC[6]= c30 + h2*c32 + h4*c43;
328 fC[7]= c31 + dS*c33;
329 fC[8]= c32 + dx*c43;
330 fC[9]= c33;
331
332 fC[10]= c40 + h2*c42 + h4*c44;
333 fC[11]= c41 + dS*c43;
334 fC[12]= c42 + dx*c44;
335 fC[13]= c43;
336 fC[14]= c44;
337
338 return ret;
339}
340
341
342
343Bool_t AliHLTTPCCATrackParam::Rotate( Float_t alpha )
344{
345 //* Rotate the coordinate system in XY on the angle alpha
346
347 Bool_t ret = 1;
348
349 Float_t cA = TMath::Cos( alpha );
350 Float_t sA = TMath::Sin( alpha );
351 Float_t x = X(), y= Y(), sP= SinPhi(), cP= CosPhi();
352
353 X() = x*cA + y*sA;
354 Y() = -x*sA + y*cA;
355 CosPhi() = cP*cA + sP*sA;
356 SinPhi() = -cP*sA + sP*cA;
357
358 Float_t j0 = 0, j2 = 0;
359
360 if( TMath::Abs(CosPhi())>1.e-4 ) j0 = cP/CosPhi(); else ret = 0;
361 if( TMath::Abs(cP)>1.e-4 ) j2 = CosPhi()/cP; else ret = 0;
362
363 //Float_t J[5][5] = { { j0, 0, 0, 0, 0 }, // Y
364 // { 0, 1, 0, 0, 0 }, // Z
365 // { 0, 0, j2, 0, 0 }, // SinPhi
366 // { 0, 0, 0, 1, 0 }, // DzDs
367 // { 0, 0, 0, 0, 1 } }; // Kappa
368
369 fC[0]*= j0*j0;
370 fC[1]*= j0;
371 //fC[3]*= j0;
372 fC[6]*= j0;
373 fC[10]*= j0;
374
375 //fC[3]*= j2;
376 fC[4]*= j2;
377 fC[5]*= j2*j2;
378 fC[8]*= j2;
379 fC[12]*= j2;
380 return ret;
381}
382
383
384void AliHLTTPCCATrackParam::GetExtParam( AliExternalTrackParam &T, Double_t alpha, Double_t Bz ) const
385{
386 //* Convert to AliExternalTrackParam parameterisation,
387 //* the angle alpha is the global angle of the local X axis
388
389 Double_t par[5], cov[15];
390 for( Int_t i=0; i<5; i++ ) par[i] = fP[i];
391 for( Int_t i=0; i<15; i++ ) cov[i] = fC[i];
392
393 if(par[2]>.999 ) par[2]=.999;
394 if(par[2]<-.999 ) par[2]=-.999;
395
396 const Double_t kCLight = 0.000299792458;
397 Double_t c = 1./(Bz*kCLight);
398 { // kappa => 1/pt
399 par[4] *= c;
400 cov[10]*= c;
401 cov[11]*= c;
402 cov[12]*= c;
403 cov[13]*= c;
404 cov[14]*= c*c;
405 }
406 if( GetCosPhi()<0 ){ // change direction
407 par[2] = -par[2]; // sin phi
408 par[3] = -par[3]; // DzDs
409 par[4] = -par[4]; // kappa
410 cov[3] = -cov[3];
411 cov[4] = -cov[4];
412 cov[6] = -cov[6];
413 cov[7] = -cov[7];
414 cov[10] = -cov[10];
415 cov[11] = -cov[11];
416 }
417 T.Set(GetX(),alpha,par,cov);
418}
419
420void AliHLTTPCCATrackParam::SetExtParam( const AliExternalTrackParam &T, Double_t Bz )
421{
422 //* Convert from AliExternalTrackParam parameterisation
423
424 for( Int_t i=0; i<5; i++ ) fP[i] = T.GetParameter()[i];
425 for( Int_t i=0; i<15; i++ ) fC[i] = T.GetCovariance()[i];
426 X() = T.GetX();
427 if(SinPhi()>.999 ) SinPhi()=.999;
428 if(SinPhi()<-.999 ) SinPhi()=-.999;
429 CosPhi() = TMath::Sqrt(1.-SinPhi()*SinPhi());
430 const Double_t kCLight = 0.000299792458;
431 Double_t c = Bz*kCLight;
432 { // 1/pt -> kappa
433 fP[4] *= c;
434 fC[10]*= c;
435 fC[11]*= c;
436 fC[12]*= c;
437 fC[13]*= c;
438 fC[14]*= c*c;
439 }
440}
441
442
443void AliHLTTPCCATrackParam::Filter( Float_t y, Float_t z, Float_t erry, Float_t errz )
444{
445 //* Add the y,z measurement with the Kalman filter
446
447 Float_t
448 c00 = fC[ 0],
449 c10 = fC[ 1], c11 = fC[ 2],
450 c20 = fC[ 3], c21 = fC[ 4],
451 c30 = fC[ 6], c31 = fC[ 7],
452 c40 = fC[10], c41 = fC[11];
453
454 Float_t
455 z0 = y-fP[0],
456 z1 = z-fP[1];
457
458 Float_t v[3] = {erry*erry, 0, errz*errz};
459
460 Float_t mS[3] = { c00+v[0], c10+v[1], c11+v[2] };
461
462 Float_t mSi[3];
463 Float_t det = (mS[0]*mS[2] - mS[1]*mS[1]);
464
465 if( TMath::Abs(det)<1.e-8 ) return;
466 det = 1./det;
467 mSi[0] = mS[2]*det;
468 mSi[1] = -mS[1]*det;
469 mSi[2] = mS[0]*det;
470
471 fNDF += 2;
472 fChi2 += ( +(mSi[0]*z0 + mSi[1]*z1 )*z0
473 +(mSi[1]*z0 + mSi[2]*z1 )*z1 );
474
475 // K = CHtS
476
477 Float_t k00, k01 , k10, k11, k20, k21, k30, k31, k40, k41;
478
479 k00 = c00*mSi[0] + c10*mSi[1]; k01 = c00*mSi[1] + c10*mSi[2];
480 k10 = c10*mSi[0] + c11*mSi[1]; k11 = c10*mSi[1] + c11*mSi[2];
481 k20 = c20*mSi[0] + c21*mSi[1]; k21 = c20*mSi[1] + c21*mSi[2];
482 k30 = c30*mSi[0] + c31*mSi[1]; k31 = c30*mSi[1] + c31*mSi[2] ;
483 k40 = c40*mSi[0] + c41*mSi[1]; k41 = c40*mSi[1] + c41*mSi[2] ;
484
485 Float_t sinPhi = fP[2] + k20*z0 + k21*z1 ;
486 if( TMath::Abs(sinPhi)>=0.99 ) return;
487
488 fP[ 0]+= k00*z0 + k01*z1 ;
489 fP[ 1]+= k10*z0 + k11*z1 ;
490 fP[ 2]+= k20*z0 + k21*z1 ;
491 fP[ 3]+= k30*z0 + k31*z1 ;
492 fP[ 4]+= k40*z0 + k41*z1 ;
493
494
495 fC[ 0]-= k00*c00 + k01*c10 ;
496
497 fC[ 1]-= k10*c00 + k11*c10 ;
498 fC[ 2]-= k10*c10 + k11*c11 ;
499
500 fC[ 3]-= k20*c00 + k21*c10 ;
501 fC[ 4]-= k20*c10 + k21*c11 ;
502 fC[ 5]-= k20*c20 + k21*c21 ;
503
504 fC[ 6]-= k30*c00 + k31*c10 ;
505 fC[ 7]-= k30*c10 + k31*c11 ;
506 fC[ 8]-= k30*c20 + k31*c21 ;
507 fC[ 9]-= k30*c30 + k31*c31 ;
508
509 fC[10]-= k40*c00 + k41*c10 ;
510 fC[11]-= k40*c10 + k41*c11 ;
511 fC[12]-= k40*c20 + k41*c21 ;
512 fC[13]-= k40*c30 + k41*c31 ;
513 fC[14]-= k40*c40 + k41*c41 ;
514
515 if( CosPhi()>=0 ){
516 CosPhi() = TMath::Sqrt(1-SinPhi()*SinPhi());
517 }else{
518 CosPhi() = -TMath::Sqrt(1-SinPhi()*SinPhi());
519 }
520
521}