]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HLT/TPCLib/tracking-ca/AliHLTTPCCATrackParam.cxx
unnecessary include removed
[u/mrichter/AliRoot.git] / HLT / TPCLib / tracking-ca / AliHLTTPCCATrackParam.cxx
CommitLineData
d54804bf 1// $Id$
ce565086 2// **************************************************************************
fbb9b71b 3// This file is property of and copyright by the ALICE HLT Project *
d54804bf 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. *
ce565086 17// *
d54804bf 18//***************************************************************************
19
ce565086 20
d54804bf 21#include "AliHLTTPCCATrackParam.h"
00d07bcd 22#include "AliHLTTPCCAMath.h"
15d2e9cf 23#include "AliHLTTPCCATrackLinearisation.h"
d54804bf 24
25//
26// Circle in XY:
15d2e9cf 27//
fbb9b71b 28// kCLight = 0.000299792458;
a8714ffa 29// Kappa = -Bz*kCLight*QPt;
d54804bf 30// R = 1/TMath::Abs(Kappa);
31// Xc = X - sin(Phi)/Kappa;
32// Yc = Y + cos(Phi)/Kappa;
33//
34
fbb9b71b 35GPUd() float AliHLTTPCCATrackParam::GetDist2( const AliHLTTPCCATrackParam &t ) const
693d2443 36{
fbb9b71b 37 // get squared distance between tracks
693d2443 38
fbb9b71b 39 float dx = GetX() - t.GetX();
40 float dy = GetY() - t.GetY();
41 float dz = GetZ() - t.GetZ();
693d2443 42 return dx*dx + dy*dy + dz*dz;
43}
44
fbb9b71b 45GPUd() float AliHLTTPCCATrackParam::GetDistXZ2( const AliHLTTPCCATrackParam &t ) const
693d2443 46{
fbb9b71b 47 // get squared distance between tracks in X&Z
693d2443 48
fbb9b71b 49 float dx = GetX() - t.GetX();
50 float dz = GetZ() - t.GetZ();
693d2443 51 return dx*dx + dz*dz;
52}
d54804bf 53
54
fbb9b71b 55GPUd() float AliHLTTPCCATrackParam::GetS( float x, float y, float Bz ) const
d54804bf 56{
57 //* Get XY path length to the given point
58
fbb9b71b 59 float k = GetKappa( Bz );
60 float ex = GetCosPhi();
61 float ey = GetSinPhi();
62 x -= GetX();
63 y -= GetY();
64 float dS = x * ex + y * ey;
65 if ( CAMath::Abs( k ) > 1.e-4 ) dS = CAMath::ATan2( k * dS, 1 + k * ( x * ey - y * ex ) ) / k;
d54804bf 66 return dS;
67}
68
fbb9b71b 69GPUd() void AliHLTTPCCATrackParam::GetDCAPoint( float x, float y, float z,
70 float &xp, float &yp, float &zp,
71 float Bz ) const
d54804bf 72{
73 //* Get the track point closest to the (x,y,z)
74
fbb9b71b 75 float x0 = GetX();
76 float y0 = GetY();
77 float k = GetKappa( Bz );
78 float ex = GetCosPhi();
79 float ey = GetSinPhi();
80 float dx = x - x0;
81 float dy = y - y0;
82 float ax = dx * k + ey;
83 float ay = dy * k - ex;
84 float a = sqrt( ax * ax + ay * ay );
85 xp = x0 + ( dx - ey * ( ( dx * dx + dy * dy ) * k - 2 * ( -dx * ey + dy * ex ) ) / ( a + 1 ) ) / a;
86 yp = y0 + ( dy + ex * ( ( dx * dx + dy * dy ) * k - 2 * ( -dx * ey + dy * ex ) ) / ( a + 1 ) ) / a;
87 float s = GetS( x, y, Bz );
88 zp = GetZ() + GetDzDs() * s;
89 if ( CAMath::Abs( k ) > 1.e-2 ) {
90 float dZ = CAMath::Abs( GetDzDs() * CAMath::TwoPi() / k );
91 if ( dZ > .1 ) {
92 zp += CAMath::Nint( ( z - zp ) / dZ ) * dZ;
d54804bf 93 }
94 }
95}
96
d54804bf 97
15d2e9cf 98//*
99//* Transport routines
100//*
d54804bf 101
d54804bf 102
fbb9b71b 103GPUd() bool AliHLTTPCCATrackParam::TransportToX( float x, AliHLTTPCCATrackLinearisation &t0, float Bz, float maxSinPhi, float *DL )
d54804bf 104{
15d2e9cf 105 //* Transport the track parameters to X=x, using linearization at t0, and the field value Bz
106 //* maxSinPhi is the max. allowed value for |t0.SinPhi()|
fbb9b71b 107 //* linearisation of trajectory t0 is also transported to X=x,
15d2e9cf 108 //* returns 1 if OK
109 //*
d54804bf 110
fbb9b71b 111 float ex = t0.CosPhi();
112 float ey = t0.SinPhi();
a8714ffa 113 float k =-t0.QPt() * Bz;
fbb9b71b 114 float dx = x - X();
115
116 float ey1 = k * dx + ey;
117 float ex1;
118
119 // check for intersection with X=x
120
121 if ( CAMath::Abs( ey1 ) > maxSinPhi ) return 0;
122
123 ex1 = CAMath::Sqrt( 1 - ey1 * ey1 );
124 if ( ex < 0 ) ex1 = -ex1;
125
126 float dx2 = dx * dx;
127 float ss = ey + ey1;
128 float cc = ex + ex1;
129
130 if ( CAMath::Abs( cc ) < 1.e-4 || CAMath::Abs( ex ) < 1.e-4 || CAMath::Abs( ex1 ) < 1.e-4 ) return 0;
131
132 float tg = ss / cc; // tan((phi1+phi)/2)
133
134 float dy = dx * tg;
135 float dl = dx * CAMath::Sqrt( 1 + tg * tg );
136
137 if ( cc < 0 ) dl = -dl;
138 float dSin = dl * k / 2;
139 if ( dSin > 1 ) dSin = 1;
140 if ( dSin < -1 ) dSin = -1;
141 float dS = ( CAMath::Abs( k ) > 1.e-4 ) ? ( 2 * CAMath::ASin( dSin ) / k ) : dl;
142 float dz = dS * t0.DzDs();
143
144 if ( DL ) *DL = -dS * CAMath::Sqrt( 1 + t0.DzDs() * t0.DzDs() );
145
146 float cci = 1. / cc;
147 float exi = 1. / ex;
148 float ex1i = 1. / ex1;
149
150 float d[5] = { 0,
4acc2401 151 0,
73a33d2e 152 GetPar(2) - t0.SinPhi(),
153 GetPar(3) - t0.DzDs(),
154 GetPar(4) - t0.QPt()
4acc2401 155 };
fbb9b71b 156
157 //float H0[5] = { 1,0, h2, 0, h4 };
158 //float H1[5] = { 0, 1, 0, dS, 0 };
159 //float H2[5] = { 0, 0, 1, 0, dxBz };
160 //float H3[5] = { 0, 0, 0, 1, 0 };
161 //float H4[5] = { 0, 0, 0, 0, 1 };
162
163 float h2 = dx * ( 1 + ey * ey1 + ex * ex1 ) * exi * ex1i * cci;
a8714ffa 164 float h4 = dx2 * ( cc + ss * ey1 * ex1i ) * cci * cci * (-Bz);
165 float dxBz = dx * (-Bz);
15d2e9cf 166
693d2443 167 t0.SetCosPhi( ex1 );
fbb9b71b 168 t0.SetSinPhi( ey1 );
693d2443 169
73a33d2e 170 SetX(X() + dx);
171 SetPar(0, Y() + dy + h2 * d[2] + h4 * d[4]);
172 SetPar(1, Z() + dz + dS * d[3]);
173 SetPar(2, t0.SinPhi() + d[2] + dxBz * d[4]);
fbb9b71b 174
175 float c00 = fC[0];
176 float c10 = fC[1];
177 float c11 = fC[2];
178 float c20 = fC[3];
179 float c21 = fC[4];
180 float c22 = fC[5];
181 float c30 = fC[6];
182 float c31 = fC[7];
183 float c32 = fC[8];
184 float c33 = fC[9];
185 float c40 = fC[10];
186 float c41 = fC[11];
187 float c42 = fC[12];
188 float c43 = fC[13];
189 float c44 = fC[14];
190
191 fC[0] = ( c00 + h2 * h2 * c22 + h4 * h4 * c44
192 + 2 * ( h2 * c20 + h4 * c40 + h2 * h4 * c42 ) );
193
194 fC[1] = c10 + h2 * c21 + h4 * c41 + dS * ( c30 + h2 * c32 + h4 * c43 );
195 fC[2] = c11 + 2 * dS * c31 + dS * dS * c33;
196
197 fC[3] = c20 + h2 * c22 + h4 * c42 + dxBz * ( c40 + h2 * c42 + h4 * c44 );
198 fC[4] = c21 + dS * c32 + dxBz * ( c41 + dS * c43 );
199 fC[5] = c22 + 2 * dxBz * c42 + dxBz * dxBz * c44;
200
201 fC[6] = c30 + h2 * c32 + h4 * c43;
202 fC[7] = c31 + dS * c33;
203 fC[8] = c32 + dxBz * c43;
204 fC[9] = c33;
205
206 fC[10] = c40 + h2 * c42 + h4 * c44;
207 fC[11] = c41 + dS * c43;
208 fC[12] = c42 + dxBz * c44;
209 fC[13] = c43;
210 fC[14] = c44;
693d2443 211
4687b8fc 212 return 1;
213}
214
693d2443 215
fbb9b71b 216GPUd() bool AliHLTTPCCATrackParam::TransportToX( float x, float sinPhi0, float cosPhi0, float Bz, float maxSinPhi )
15d2e9cf 217{
fbb9b71b 218 //* Transport the track parameters to X=x, using linearization at phi0 with 0 curvature,
15d2e9cf 219 //* and the field value Bz
220 //* maxSinPhi is the max. allowed value for |t0.SinPhi()|
fbb9b71b 221 //* linearisation of trajectory t0 is also transported to X=x,
15d2e9cf 222 //* returns 1 if OK
223 //*
224
fbb9b71b 225 float ex = cosPhi0;
226 float ey = sinPhi0;
227 float dx = x - X();
228
229 if ( CAMath::Abs( ex ) < 1.e-4 ) return 0;
230 float exi = 1. / ex;
693d2443 231
a8714ffa 232 float dxBz = dx * (-Bz);
fbb9b71b 233 float dS = dx * exi;
234 float h2 = dS * exi * exi;
235 float h4 = .5 * h2 * dxBz;
693d2443 236
fbb9b71b 237 //float H0[5] = { 1,0, h2, 0, h4 };
238 //float H1[5] = { 0, 1, 0, dS, 0 };
239 //float H2[5] = { 0, 0, 1, 0, dxBz };
240 //float H3[5] = { 0, 0, 0, 1, 0 };
241 //float H4[5] = { 0, 0, 0, 0, 1 };
693d2443 242
fbb9b71b 243 float sinPhi = SinPhi() + dxBz * QPt();
244 if ( maxSinPhi > 0 && CAMath::Abs( sinPhi ) > maxSinPhi ) return 0;
4687b8fc 245
73a33d2e 246 SetX(X() + dx);
247 SetPar(0, GetPar(0) + dS * ey + h2 * ( SinPhi() - ey ) + h4 * QPt());
248 SetPar(1, GetPar(1) + dS * DzDs());
249 SetPar(2, sinPhi);
eb30eb49 250
fbb9b71b 251
252 float c00 = fC[0];
253 float c10 = fC[1];
254 float c11 = fC[2];
255 float c20 = fC[3];
256 float c21 = fC[4];
257 float c22 = fC[5];
258 float c30 = fC[6];
259 float c31 = fC[7];
260 float c32 = fC[8];
261 float c33 = fC[9];
262 float c40 = fC[10];
263 float c41 = fC[11];
264 float c42 = fC[12];
265 float c43 = fC[13];
266 float c44 = fC[14];
267
268
269 fC[0] = ( c00 + h2 * h2 * c22 + h4 * h4 * c44
270 + 2 * ( h2 * c20 + h4 * c40 + h2 * h4 * c42 ) );
271
272 fC[1] = c10 + h2 * c21 + h4 * c41 + dS * ( c30 + h2 * c32 + h4 * c43 );
273 fC[2] = c11 + 2 * dS * c31 + dS * dS * c33;
274
275 fC[3] = c20 + h2 * c22 + h4 * c42 + dxBz * ( c40 + h2 * c42 + h4 * c44 );
276 fC[4] = c21 + dS * c32 + dxBz * ( c41 + dS * c43 );
277 fC[5] = c22 + 2 * dxBz * c42 + dxBz * dxBz * c44;
278
279 fC[6] = c30 + h2 * c32 + h4 * c43;
280 fC[7] = c31 + dS * c33;
281 fC[8] = c32 + dxBz * c43;
282 fC[9] = c33;
283
284 fC[10] = c40 + h2 * c42 + h4 * c44;
285 fC[11] = c41 + dS * c43;
286 fC[12] = c42 + dxBz * c44;
287 fC[13] = c43;
288 fC[14] = c44;
eb30eb49 289
15d2e9cf 290 return 1;
eb30eb49 291}
292
293
294
693d2443 295
693d2443 296
693d2443 297
fbb9b71b 298GPUd() bool AliHLTTPCCATrackParam::TransportToX( float x, float Bz, float maxSinPhi )
15d2e9cf 299{
fbb9b71b 300 //* Transport the track parameters to X=x
693d2443 301
fbb9b71b 302 AliHLTTPCCATrackLinearisation t0( *this );
693d2443 303
15d2e9cf 304 return TransportToX( x, t0, Bz, maxSinPhi );
305}
693d2443 306
693d2443 307
693d2443 308
fbb9b71b 309GPUd() bool AliHLTTPCCATrackParam::TransportToXWithMaterial( float x, AliHLTTPCCATrackLinearisation &t0, AliHLTTPCCATrackFitParam &par, float Bz, float maxSinPhi )
15d2e9cf 310{
311 //* Transport the track parameters to X=x taking into account material budget
693d2443 312
fbb9b71b 313 const float kRho = 1.025e-3;//0.9e-3;
314 const float kRadLen = 29.532;//28.94;
315 const float kRhoOverRadLen = kRho / kRadLen;
316 float dl;
693d2443 317
fbb9b71b 318 if ( !TransportToX( x, t0, Bz, maxSinPhi, &dl ) ) return 0;
693d2443 319
fbb9b71b 320 CorrectForMeanMaterial( dl*kRhoOverRadLen, dl*kRho, par );
15d2e9cf 321 return 1;
322}
693d2443 323
693d2443 324
fbb9b71b 325GPUd() bool AliHLTTPCCATrackParam::TransportToXWithMaterial( float x, AliHLTTPCCATrackFitParam &par, float Bz, float maxSinPhi )
15d2e9cf 326{
327 //* Transport the track parameters to X=x taking into account material budget
693d2443 328
fbb9b71b 329 AliHLTTPCCATrackLinearisation t0( *this );
15d2e9cf 330 return TransportToXWithMaterial( x, t0, par, Bz, maxSinPhi );
331}
693d2443 332
fbb9b71b 333GPUd() bool AliHLTTPCCATrackParam::TransportToXWithMaterial( float x, float Bz, float maxSinPhi )
15d2e9cf 334{
335 //* Transport the track parameters to X=x taking into account material budget
336
337 AliHLTTPCCATrackFitParam par;
338 CalculateFitParameters( par );
fbb9b71b 339 return TransportToXWithMaterial( x, par, Bz, maxSinPhi );
15d2e9cf 340}
693d2443 341
693d2443 342
15d2e9cf 343//*
344//* Multiple scattering and energy losses
345//*
693d2443 346
15d2e9cf 347
fbb9b71b 348float AliHLTTPCCATrackParam::BetheBlochGeant( float bg2,
349 float kp0,
350 float kp1,
351 float kp2,
352 float kp3,
353 float kp4 )
15d2e9cf 354{
355 //
356 // This is the parameterization of the Bethe-Bloch formula inspired by Geant.
357 //
358 // bg2 - (beta*gamma)^2
359 // kp0 - density [g/cm^3]
360 // kp1 - density effect first junction point
361 // kp2 - density effect second junction point
362 // kp3 - mean excitation energy [GeV]
363 // kp4 - mean Z/A
364 //
fbb9b71b 365 // The default values for the kp* parameters are for silicon.
15d2e9cf 366 // The returned value is in [GeV/(g/cm^2)].
fbb9b71b 367 //
368
369 const float mK = 0.307075e-3; // [GeV*cm^2/g]
370 const float me = 0.511e-3; // [GeV/c^2]
371 const float rho = kp0;
372 const float x0 = kp1 * 2.303;
373 const float x1 = kp2 * 2.303;
374 const float mI = kp3;
375 const float mZA = kp4;
376 const float maxT = 2 * me * bg2; // neglecting the electron mass
377
15d2e9cf 378 //*** Density effect
fbb9b71b 379 float d2 = 0.;
7be9b0d7 380 const float x = 0.5 * AliHLTTPCCAMath::Log( bg2 );
381 const float lhwI = AliHLTTPCCAMath::Log( 28.816 * 1e-9 * AliHLTTPCCAMath::Sqrt( rho * mZA ) / mI );
fbb9b71b 382 if ( x > x1 ) {
15d2e9cf 383 d2 = lhwI + x - 0.5;
fbb9b71b 384 } else if ( x > x0 ) {
385 const float r = ( x1 - x ) / ( x1 - x0 );
386 d2 = lhwI + x - 0.5 + ( 0.5 - lhwI - x0 ) * r * r * r;
693d2443 387 }
388
7be9b0d7 389 return mK*mZA*( 1 + bg2 ) / bg2*( 0.5*AliHLTTPCCAMath::Log( 2*me*bg2*maxT / ( mI*mI ) ) - bg2 / ( 1 + bg2 ) - d2 );
15d2e9cf 390}
391
fbb9b71b 392float AliHLTTPCCATrackParam::BetheBlochSolid( float bg )
15d2e9cf 393{
394 //------------------------------------------------------------------
fbb9b71b 395 // This is an approximation of the Bethe-Bloch formula,
396 // reasonable for solid materials.
15d2e9cf 397 // All the parameters are, in fact, for Si.
398 // The returned value is in [GeV]
399 //------------------------------------------------------------------
400
fbb9b71b 401 return BetheBlochGeant( bg );
15d2e9cf 402}
403
fbb9b71b 404float AliHLTTPCCATrackParam::BetheBlochGas( float bg )
15d2e9cf 405{
406 //------------------------------------------------------------------
fbb9b71b 407 // This is an approximation of the Bethe-Bloch formula,
15d2e9cf 408 // reasonable for gas materials.
409 // All the parameters are, in fact, for Ne.
410 // The returned value is in [GeV]
411 //------------------------------------------------------------------
412
fbb9b71b 413 const float rho = 0.9e-3;
414 const float x0 = 2.;
415 const float x1 = 4.;
416 const float mI = 140.e-9;
417 const float mZA = 0.49555;
15d2e9cf 418
fbb9b71b 419 return BetheBlochGeant( bg, rho, x0, x1, mI, mZA );
693d2443 420}
421
422
15d2e9cf 423
424
fbb9b71b 425GPUd() float AliHLTTPCCATrackParam::ApproximateBetheBloch( float beta2 )
eb30eb49 426{
427 //------------------------------------------------------------------
fbb9b71b 428 // This is an approximation of the Bethe-Bloch formula with
eb30eb49 429 // the density effect taken into account at beta*gamma > 3.5
fbb9b71b 430 // (the approximation is reasonable only for solid materials)
eb30eb49 431 //------------------------------------------------------------------
fbb9b71b 432 if ( beta2 >= 1 ) return 0;
eb30eb49 433
fbb9b71b 434 if ( beta2 / ( 1 - beta2 ) > 3.5*3.5 )
435 return 0.153e-3 / beta2*( log( 3.5*5940 ) + 0.5*log( beta2 / ( 1 - beta2 ) ) - beta2 );
436 return 0.153e-3 / beta2*( log( 5940*beta2 / ( 1 - beta2 ) ) - beta2 );
eb30eb49 437}
438
439
fbb9b71b 440GPUd() void AliHLTTPCCATrackParam::CalculateFitParameters( AliHLTTPCCATrackFitParam &par, float mass )
eb30eb49 441{
442 //*!
443
73a33d2e 444 float qpt = GetPar(4);
a371e863 445 if( fC[14]>=1. ) qpt = 1./0.35;
446
73a33d2e 447 float p2 = ( 1. + GetPar(3) * GetPar(3) );
a371e863 448 float k2 = qpt * qpt;
fbb9b71b 449 float mass2 = mass * mass;
450 float beta2 = p2 / ( p2 + mass2 * k2 );
15d2e9cf 451
fbb9b71b 452 float pp2 = ( k2 > 1.e-8 ) ? p2 / k2 : 10000; // impuls 2
eb30eb49 453
15d2e9cf 454 //par.fBethe = BetheBlochGas( pp2/mass2);
fbb9b71b 455 par.fBethe = ApproximateBetheBloch( pp2 / mass2 );
456 par.fE = CAMath::Sqrt( pp2 + mass2 );
457 par.fTheta2 = 14.1 * 14.1 / ( beta2 * pp2 * 1e6 );
458 par.fEP2 = par.fE / pp2;
eb30eb49 459
460 // Approximate energy loss fluctuation (M.Ivanov)
fbb9b71b 461
462 const float knst = 0.07; // To be tuned.
a371e863 463 par.fSigmadE2 = knst * par.fEP2 * qpt;
eb30eb49 464 par.fSigmadE2 = par.fSigmadE2 * par.fSigmadE2;
fbb9b71b 465
73a33d2e 466 par.fK22 = ( 1. + GetPar(3) * GetPar(3) );
fbb9b71b 467 par.fK33 = par.fK22 * par.fK22;
b4b9adac 468 par.fK43 = 0;
73a33d2e 469 par.fK44 = GetPar(3) * GetPar(3) * k2;
15d2e9cf 470
eb30eb49 471}
472
473
fbb9b71b 474GPUd() bool AliHLTTPCCATrackParam::CorrectForMeanMaterial( float xOverX0, float xTimesRho, const AliHLTTPCCATrackFitParam &par )
eb30eb49 475{
476 //------------------------------------------------------------------
477 // This function corrects the track parameters for the crossed material.
478 // "xOverX0" - X/X0, the thickness in units of the radiation length.
fbb9b71b 479 // "xTimesRho" - is the product length*density (g/cm^2).
eb30eb49 480 //------------------------------------------------------------------
481
fbb9b71b 482 float &fC22 = fC[5];
483 float &fC33 = fC[9];
484 float &fC40 = fC[10];
485 float &fC41 = fC[11];
486 float &fC42 = fC[12];
487 float &fC43 = fC[13];
488 float &fC44 = fC[14];
eb30eb49 489
490 //Energy losses************************
fbb9b71b 491
492 float dE = par.fBethe * xTimesRho;
493 if ( CAMath::Abs( dE ) > 0.3*par.fE ) return 0; //30% energy loss is too much!
494 float corr = ( 1. - par.fEP2 * dE );
495 if ( corr < 0.3 || corr > 1.3 ) return 0;
496
73a33d2e 497 SetPar(4, GetPar(4) * corr);
fbb9b71b 498 fC40 *= corr;
499 fC41 *= corr;
500 fC42 *= corr;
501 fC43 *= corr;
502 fC44 *= corr * corr;
503 fC44 += par.fSigmadE2 * CAMath::Abs( dE );
eb30eb49 504
505 //Multiple scattering******************
fbb9b71b 506
507 float theta2 = par.fTheta2 * CAMath::Abs( xOverX0 );
bfd20868 508 fC22 += theta2 * par.fK22 * (1.-GetPar(2))*(1.+GetPar(2));
fbb9b71b 509 fC33 += theta2 * par.fK33;
510 fC43 += theta2 * par.fK43;
511 fC44 += theta2 * par.fK44;
15d2e9cf 512
eb30eb49 513 return 1;
514}
515
516
15d2e9cf 517//*
518//* Rotation
519//*
520
eb30eb49 521
fbb9b71b 522GPUd() bool AliHLTTPCCATrackParam::Rotate( float alpha, float maxSinPhi )
eb30eb49 523{
524 //* Rotate the coordinate system in XY on the angle alpha
fbb9b71b 525
526 float cA = CAMath::Cos( alpha );
527 float sA = CAMath::Sin( alpha );
528 float x = X(), y = Y(), sP = SinPhi(), cP = GetCosPhi();
529 float cosPhi = cP * cA + sP * sA;
530 float sinPhi = -cP * sA + sP * cA;
531
532 if ( CAMath::Abs( sinPhi ) > maxSinPhi || CAMath::Abs( cosPhi ) < 1.e-2 || CAMath::Abs( cP ) < 1.e-2 ) return 0;
533
534 float j0 = cP / cosPhi;
535 float j2 = cosPhi / cP;
536
693d2443 537 SetX( x*cA + y*sA );
fbb9b71b 538 SetY( -x*sA + y*cA );
15d2e9cf 539 SetSignCosPhi( cosPhi );
693d2443 540 SetSinPhi( sinPhi );
d54804bf 541
d54804bf 542
fbb9b71b 543 //float J[5][5] = { { j0, 0, 0, 0, 0 }, // Y
d54804bf 544 // { 0, 1, 0, 0, 0 }, // Z
545 // { 0, 0, j2, 0, 0 }, // SinPhi
fbb9b71b 546 // { 0, 0, 0, 1, 0 }, // DzDs
547 // { 0, 0, 0, 0, 1 } }; // Kappa
eb30eb49 548 //cout<<"alpha="<<alpha<<" "<<x<<" "<<y<<" "<<sP<<" "<<cP<<" "<<j0<<" "<<j2<<endl;
549 //cout<<" "<<fC[0]<<" "<<fC[1]<<" "<<fC[6]<<" "<<fC[10]<<" "<<fC[4]<<" "<<fC[5]<<" "<<fC[8]<<" "<<fC[12]<<endl;
fbb9b71b 550 fC[0] *= j0 * j0;
551 fC[1] *= j0;
552 fC[3] *= j0;
553 fC[6] *= j0;
554 fC[10] *= j0;
555
556 fC[3] *= j2;
557 fC[4] *= j2;
558 fC[5] *= j2 * j2;
559 fC[8] *= j2;
560 fC[12] *= j2;
eb30eb49 561 //cout<<" "<<fC[0]<<" "<<fC[1]<<" "<<fC[6]<<" "<<fC[10]<<" "<<fC[4]<<" "<<fC[5]<<" "<<fC[8]<<" "<<fC[12]<<endl;
562 return 1;
d54804bf 563}
564
fbb9b71b 565GPUd() bool AliHLTTPCCATrackParam::Rotate( float alpha, AliHLTTPCCATrackLinearisation &t0, float maxSinPhi )
693d2443 566{
567 //* Rotate the coordinate system in XY on the angle alpha
fbb9b71b 568
569 float cA = CAMath::Cos( alpha );
570 float sA = CAMath::Sin( alpha );
571 float x0 = X(), y0 = Y(), sP = t0.SinPhi(), cP = t0.CosPhi();
572 float cosPhi = cP * cA + sP * sA;
573 float sinPhi = -cP * sA + sP * cA;
574
575 if ( CAMath::Abs( sinPhi ) > maxSinPhi || CAMath::Abs( cosPhi ) < 1.e-2 || CAMath::Abs( cP ) < 1.e-2 ) return 0;
576
577 //float J[5][5] = { { j0, 0, 0, 0, 0 }, // Y
15d2e9cf 578 // { 0, 1, 0, 0, 0 }, // Z
579 // { 0, 0, j2, 0, 0 }, // SinPhi
fbb9b71b 580 // { 0, 0, 0, 1, 0 }, // DzDs
581 // { 0, 0, 0, 0, 1 } }; // Kappa
15d2e9cf 582
fbb9b71b 583 float j0 = cP / cosPhi;
584 float j2 = cosPhi / cP;
585 float d[2] = {Y() - y0, SinPhi() - sP};
693d2443 586
15d2e9cf 587 SetX( x0*cA + y0*sA );
fbb9b71b 588 SetY( -x0*sA + y0*cA + j0*d[0] );
693d2443 589 t0.SetCosPhi( cosPhi );
590 t0.SetSinPhi( sinPhi );
591
693d2443 592 SetSinPhi( sinPhi + j2*d[1] );
593
fbb9b71b 594 fC[0] *= j0 * j0;
595 fC[1] *= j0;
596 fC[3] *= j0;
597 fC[6] *= j0;
598 fC[10] *= j0;
693d2443 599
fbb9b71b 600 fC[3] *= j2;
601 fC[4] *= j2;
602 fC[5] *= j2 * j2;
603 fC[8] *= j2;
604 fC[12] *= j2;
693d2443 605
606 return 1;
607}
608
fbb9b71b 609GPUd() bool AliHLTTPCCATrackParam::Filter( float y, float z, float err2Y, float err2Z, float maxSinPhi )
4687b8fc 610{
fbb9b71b 611 //* Add the y,z measurement with the Kalman filter
612
613 float
614 c00 = fC[ 0],
615 c11 = fC[ 2],
616 c20 = fC[ 3],
617 c31 = fC[ 7],
618 c40 = fC[10];
619
620 err2Y += c00;
621 err2Z += c11;
622
623 float
73a33d2e 624 z0 = y - GetPar(0),
625 z1 = z - GetPar(1);
fbb9b71b 626
627 if ( err2Y < 1.e-8 || err2Z < 1.e-8 ) return 0;
628
629 float mS0 = 1. / err2Y;
630 float mS2 = 1. / err2Z;
631
4687b8fc 632 // K = CHtS
fbb9b71b 633
634 float k00, k11, k20, k31, k40;
635
636 k00 = c00 * mS0;
637 k20 = c20 * mS0;
638 k40 = c40 * mS0;
639
640 k11 = c11 * mS2;
641 k31 = c31 * mS2;
642
73a33d2e 643 float sinPhi = GetPar(2) + k20 * z0 ;
fbb9b71b 644
645 if ( maxSinPhi > 0 && CAMath::Abs( sinPhi ) >= maxSinPhi ) return 0;
646
693d2443 647 fNDF += 2;
fbb9b71b 648 fChi2 += mS0 * z0 * z0 + mS2 * z1 * z1 ;
693d2443 649
73a33d2e 650 SetPar(0, GetPar(0) + k00 * z0);
651 SetPar(1, GetPar(1) + k11 * z1);
652 SetPar(2, sinPhi);
653 SetPar(3, GetPar(3) + k31 * z1);
654 SetPar(4, GetPar(4) + k40 * z0);
fbb9b71b 655
656 fC[ 0] -= k00 * c00 ;
657 fC[ 3] -= k20 * c00 ;
658 fC[ 5] -= k20 * c20 ;
659 fC[10] -= k40 * c00 ;
660 fC[12] -= k40 * c20 ;
661 fC[14] -= k40 * c40 ;
662
663 fC[ 2] -= k11 * c11 ;
664 fC[ 7] -= k31 * c11 ;
665 fC[ 9] -= k31 * c31 ;
666
693d2443 667 return 1;
668}
669
4acc2401 670GPUd() bool AliHLTTPCCATrackParam::CheckNumericalQuality() const
671{
672 //* Check that the track parameters and covariance matrix are reasonable
673
73a33d2e 674 bool ok = AliHLTTPCCAMath::Finite( GetX() ) && AliHLTTPCCAMath::Finite( fSignCosPhi ) && AliHLTTPCCAMath::Finite( fChi2 ) && AliHLTTPCCAMath::Finite( fNDF );
4acc2401 675
676 const float *c = Cov();
7be9b0d7 677 for ( int i = 0; i < 15; i++ ) ok = ok && AliHLTTPCCAMath::Finite( c[i] );
678 for ( int i = 0; i < 5; i++ ) ok = ok && AliHLTTPCCAMath::Finite( Par()[i] );
693d2443 679
4acc2401 680 if ( c[0] <= 0 || c[2] <= 0 || c[5] <= 0 || c[9] <= 0 || c[14] <= 0 ) ok = 0;
8f4bbe22 681 if ( c[0] > 5. || c[2] > 5. || c[5] > 2. || c[9] > 2
682 //|| ( CAMath::Abs( QPt() ) > 1.e-2 && c[14] > 2. )
683 ) ok = 0;
4acc2401 684
685 if ( CAMath::Abs( SinPhi() ) > .99 ) ok = 0;
686 if ( CAMath::Abs( QPt() ) > 1. / 0.05 ) ok = 0;
21ea647f 687 if( ok ){
688 ok = ok
689 && ( c[1]*c[1]<=c[2]*c[0] )
690 && ( c[3]*c[3]<=c[5]*c[0] )
691 && ( c[4]*c[4]<=c[5]*c[2] )
692 && ( c[6]*c[6]<=c[9]*c[0] )
693 && ( c[7]*c[7]<=c[9]*c[2] )
694 && ( c[8]*c[8]<=c[9]*c[5] )
695 && ( c[10]*c[10]<=c[14]*c[0] )
696 && ( c[11]*c[11]<=c[14]*c[2] )
697 && ( c[12]*c[12]<=c[14]*c[5] )
698 && ( c[13]*c[13]<=c[14]*c[9] );
699 }
4acc2401 700 return ok;
701}
eb30eb49 702
eb30eb49 703
00d07bcd 704#if !defined(HLTCA_GPUCODE)
4687b8fc 705#include <iostream>
00d07bcd 706#endif
707
708GPUd() void AliHLTTPCCATrackParam::Print() const
eb30eb49 709{
00d07bcd 710 //* print parameters
fbb9b71b 711
00d07bcd 712#if !defined(HLTCA_GPUCODE)
fbb9b71b 713 std::cout << "track: x=" << GetX() << " c=" << GetSignCosPhi() << ", P= " << GetY() << " " << GetZ() << " " << GetSinPhi() << " " << GetDzDs() << " " << GetQPt() << std::endl;
714 std::cout << "errs2: " << GetErr2Y() << " " << GetErr2Z() << " " << GetErr2SinPhi() << " " << GetErr2DzDs() << " " << GetErr2QPt() << std::endl;
00d07bcd 715#endif
d54804bf 716}
15d2e9cf 717