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