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