Line data Source code
1 : // $Id: AliHLTTPCGMSliceTrack.cxx 41769 2010-06-16 13:58:00Z sgorbuno $
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 : // for The ALICE HLT Project. *
8 : // *
9 : // Permission to use, copy, modify and distribute this software and its *
10 : // documentation strictly for non-commercial purposes is hereby granted *
11 : // without fee, provided that the above copyright notice appears in all *
12 : // copies and that both the copyright notice and this permission notice *
13 : // appear in the supporting documentation. The authors make no claims *
14 : // about the suitability of this software for any purpose. It is *
15 : // provided "as is" without express or implied warranty. *
16 : // *
17 : //***************************************************************************
18 :
19 : #include "AliHLTTPCGMSliceTrack.h"
20 : #include "AliHLTTPCCAMath.h"
21 : #include "AliHLTTPCGMBorderTrack.h"
22 : #include "AliHLTTPCGMTrackLinearisation.h"
23 : #include "Riostream.h"
24 : #include "AliHLTTPCCAParam.h"
25 : #include <cmath>
26 :
27 :
28 :
29 : bool AliHLTTPCGMSliceTrack::FilterErrors( AliHLTTPCCAParam ¶m, float maxSinPhi )
30 : {
31 0 : float lastX = fOrigTrack->Cluster(fOrigTrack->NClusters()-1 ).GetX();
32 :
33 : const int N = 3;
34 0 : const float *cy = param.GetParamS0Par(0,0);
35 0 : const float *cz = param.GetParamS0Par(1,0);
36 :
37 0 : float bz = -param.ConstBz();
38 : const float kZLength = 250.f - 0.275f;
39 :
40 0 : float trDzDs2 = fDzDs*fDzDs;
41 0 : float k = fQPt*bz;
42 0 : float dx = .33*(lastX - fX);
43 0 : float kdx = k*dx;
44 0 : float dxBz = dx * bz;
45 0 : float kdx205 = 2.f+kdx*kdx*0.5f;
46 :
47 : {
48 0 : float secPhi2 = fSecPhi*fSecPhi;
49 0 : float zz = fabs( kZLength - fabs(fZ) );
50 0 : float zz2 = zz*zz;
51 0 : float angleY2 = secPhi2 - 1.f;
52 0 : float angleZ2 = trDzDs2 * secPhi2 ;
53 :
54 0 : float cy0 = cy[0] + cy[1]*zz + cy[3]*zz2;
55 0 : float cy1 = cy[2] + cy[5]*zz;
56 0 : float cy2 = cy[4];
57 0 : float cz0 = cz[0] + cz[1]*zz + cz[3]*zz2;
58 0 : float cz1 = cz[2] + cz[5]*zz;
59 0 : float cz2 = cz[4];
60 :
61 0 : fC0 = fabs( cy0 + angleY2 * ( cy1 + angleY2*cy2 ) );
62 0 : fC2 = fabs( cz0 + angleZ2 * ( cz1 + angleZ2*cz2 ) );
63 :
64 0 : fC3 = 0;
65 0 : fC5 = 1;
66 0 : fC7 = 0;
67 0 : fC9 = 1;
68 0 : fC10 = 0;
69 0 : fC12 = 0;
70 0 : fC14 = 10;
71 : }
72 :
73 0 : for( int iStep=0; iStep<N; iStep++ ){
74 :
75 : float err2Y, err2Z;
76 :
77 : { // transport block
78 :
79 0 : float ex = fCosPhi;
80 0 : float ey = fSinPhi;
81 0 : float ey1 = kdx + ey;
82 0 : if( fabs( ey1 ) > maxSinPhi ) return 0;
83 :
84 0 : float ss = ey + ey1;
85 0 : float ex1 = sqrt(1.f - ey1*ey1);
86 :
87 0 : float cc = ex + ex1;
88 0 : float dxcci = dx / cc;
89 :
90 0 : float dy = dxcci * ss;
91 0 : float norm2 = 1.f + ey*ey1 + ex*ex1;
92 0 : float dl = dxcci * sqrt( norm2 + norm2 );
93 :
94 : float dS;
95 : {
96 0 : float dSin = 0.5f*k*dl;
97 0 : float a = dSin*dSin;
98 : const float k2 = 1.f/6.f;
99 : const float k4 = 3.f/40.f;
100 0 : dS = dl + dl*a*(k2 + a*(k4 ));//+ k6*a) );
101 : }
102 :
103 0 : float dz = dS * fDzDs;
104 0 : float ex1i =1.f/ex1;
105 0 : float z = fZ+ dz;
106 : {
107 0 : float secPhi2 = ex1i*ex1i;
108 0 : float zz = fabs( kZLength - fabs(z) );
109 0 : float zz2 = zz*zz;
110 0 : float angleY2 = secPhi2 - 1.f;
111 0 : float angleZ2 = trDzDs2 * secPhi2 ;
112 :
113 0 : float cy0 = cy[0] + cy[1]*zz + cy[3]*zz2;
114 0 : float cy1 = cy[2] + cy[5]*zz;
115 0 : float cy2 = cy[4];
116 0 : float cz0 = cz[0] + cz[1]*zz + cz[3]*zz2;
117 0 : float cz1 = cz[2] + cz[5]*zz;
118 0 : float cz2 = cz[4];
119 :
120 0 : err2Y = fabs( cy0 + angleY2 * ( cy1 + angleY2*cy2 ) );
121 0 : err2Z = fabs( cz0 + angleZ2 * ( cz1 + angleZ2*cz2 ) );
122 : }
123 :
124 0 : float hh = kdx205 * dxcci*ex1i;
125 0 : float h2 = hh * fSecPhi;
126 :
127 0 : fX+=dx;
128 0 : fY+= dy;
129 0 : fZ+= dz;
130 0 : fSinPhi = ey1;
131 0 : fCosPhi = ex1;
132 0 : fSecPhi = ex1i;
133 :
134 0 : float h4 = bz*dxcci*hh;
135 :
136 0 : float c20 = fC3;
137 0 : float c22 = fC5;
138 0 : float c31 = fC7;
139 0 : float c33 = fC9;
140 0 : float c40 = fC10;
141 0 : float c42 = fC12;
142 0 : float c44 = fC14;
143 :
144 0 : float c20ph4c42 = c20 + h4*c42;
145 0 : float h2c22 = h2*c22;
146 0 : float h4c44 = h4*c44;
147 0 : float n7 = c31 + dS*c33;
148 0 : float n10 = c40 + h2*c42 + h4c44;
149 0 : float n12 = c42 + dxBz*c44;
150 :
151 :
152 0 : fC0+= h2*h2c22 + h4*h4c44 + 2.f*( h2*c20ph4c42 + h4*c40 );
153 :
154 0 : fC3 = c20ph4c42 + h2c22 + dxBz*n10;
155 0 : fC10 = n10;
156 :
157 0 : fC5 = c22 + dxBz*( c42 + n12 );
158 0 : fC12 = n12;
159 :
160 0 : fC2+= dS*(c31 + n7);
161 0 : fC7 = n7;
162 :
163 0 : } // end transport block
164 :
165 :
166 : // Filter block
167 :
168 : float
169 0 : c00 = fC0,
170 0 : c11 = fC2,
171 0 : c20 = fC3,
172 0 : c31 = fC7,
173 0 : c40 = fC10;
174 :
175 0 : float mS0 = 1.f/(err2Y + c00);
176 0 : float mS2 = 1.f/(err2Z + c11);
177 :
178 : // K = CHtS
179 :
180 : float k00, k11, k20, k31, k40;
181 :
182 0 : k00 = c00 * mS0;
183 0 : k20 = c20 * mS0;
184 0 : k40 = c40 * mS0;
185 :
186 0 : fC0 -= k00 * c00 ;
187 0 : fC5 -= k20 * c20 ;
188 0 : fC10 -= k00 * c40 ;
189 0 : fC12 -= k40 * c20 ;
190 0 : fC3 -= k20 * c00 ;
191 0 : fC14 -= k40 * c40 ;
192 :
193 0 : k11 = c11 * mS2;
194 0 : k31 = c31 * mS2;
195 :
196 0 : fC7 -= k31 * c11;
197 0 : fC2 -= k11 * c11;
198 0 : fC9 -= k31 * c31;
199 0 : }
200 :
201 : //* Check that the track parameters and covariance matrix are reasonable
202 :
203 : bool ok = 1;
204 :
205 : const float *c = &fX;
206 0 : for ( int i = 0; i < 17; i++ ) ok = ok && finite( c[i] );
207 :
208 0 : if ( fC0 <= 0.f || fC2 <= 0.f || fC5 <= 0.f || fC9 <= 0.f || fC14 <= 0.f
209 0 : || fC0 > 5.f || fC2 > 5.f || fC5 > 2.f || fC9 > 2.f ) ok = 0;
210 :
211 0 : if( ok ){
212 0 : ok = ok
213 0 : && ( fC3*fC3<=fC5*fC0 )
214 0 : && ( fC7*fC7<=fC9*fC2 )
215 0 : && ( fC10*fC10<=fC14*fC0 )
216 0 : && ( fC12*fC12<=fC14*fC5 );
217 0 : }
218 :
219 0 : return ok;
220 0 : }
221 :
222 :
223 :
224 : bool AliHLTTPCGMSliceTrack::TransportToX( float x, float Bz, AliHLTTPCGMBorderTrack &b, float maxSinPhi ) const
225 : {
226 0 : Bz = -Bz;
227 0 : float ex = fCosPhi;
228 0 : float ey = fSinPhi;
229 0 : float k = fQPt*Bz;
230 0 : float dx = x - fX;
231 0 : float ey1 = k*dx + ey;
232 :
233 0 : if( fabs( ey1 ) > maxSinPhi ) return 0;
234 :
235 0 : float ex1 = sqrt( 1.f - ey1 * ey1 );
236 0 : float dxBz = dx * Bz;
237 :
238 0 : float ss = ey + ey1;
239 0 : float cc = ex + ex1;
240 0 : float dxcci = dx / cc;
241 0 : float norm2 = 1.f + ey*ey1 + ex*ex1;
242 :
243 0 : float dy = dxcci * ss;
244 :
245 : float dS;
246 : {
247 0 : float dl = dxcci * sqrt( norm2 + norm2 );
248 0 : float dSin = 0.5f*k*dl;
249 0 : float a = dSin*dSin;
250 : const float k2 = 1.f/6.f;
251 : const float k4 = 3.f/40.f;
252 : //const float k6 = 5.f/112.f;
253 0 : dS = dl + dl*a*(k2 + a*(k4 ));//+ k6*a) );
254 : }
255 :
256 0 : float ex1i = 1.f/ex1;
257 0 : float dz = dS * fDzDs;
258 :
259 0 : float hh = dxcci*ex1i*norm2;
260 0 : float h2 = hh *fSecPhi;
261 0 : float h4 = Bz*dxcci*hh;
262 :
263 0 : float c20 = fC3;
264 0 : float c22 = fC5;
265 0 : float c31 = fC7;
266 0 : float c33 = fC9;
267 0 : float c40 = fC10;
268 0 : float c42 = fC12;
269 0 : float c44 = fC14;
270 :
271 0 : float c20ph4c42 = c20 + h4*c42;
272 0 : float h2c22 = h2*c22;
273 0 : float h4c44 = h4*c44;
274 0 : float n7 = c31 + dS*c33;
275 :
276 0 : b.SetPar(0, fY + dy );
277 0 : b.SetPar(1, fZ + dz );
278 0 : b.SetPar(2, ey1 );
279 0 : b.SetPar(3, fDzDs);
280 0 : b.SetPar(4, fQPt);
281 :
282 0 : b.SetCov(0, fC0 + h2*h2c22 + h4*h4c44 + 2.f*( h2*c20ph4c42 + h4*c40 ));
283 0 : b.SetCov(1, fC2+ dS*(c31 + n7) );
284 0 : b.SetCov(2, c22 + dxBz*( c42 + c42 + dxBz*c44 ));
285 0 : b.SetCov(3, c33);
286 0 : b.SetCov(4, c44);
287 0 : b.SetCovD(0, c20ph4c42 + h2c22 + dxBz*(c40 + h2*c42 + h4c44) );
288 0 : b.SetCovD(1, n7 );
289 : return 1;
290 0 : }
291 :
292 :
293 :
294 : bool AliHLTTPCGMSliceTrack::TransportToXAlpha( float newX, float sinAlpha, float cosAlpha, float Bz, AliHLTTPCGMBorderTrack &b, float maxSinPhi ) const
295 : {
296 : //*
297 :
298 0 : float c00 = fC0;
299 0 : float c11 = fC2;
300 0 : float c20 = fC3;
301 0 : float c22 = fC5;
302 0 : float c31 = fC7;
303 0 : float c33 = fC9;
304 0 : float c40 = fC10;
305 0 : float c42 = fC12;
306 0 : float c44 = fC14;
307 :
308 : float x,y;
309 0 : float z = fZ;
310 0 : float sinPhi = fSinPhi;
311 0 : float cosPhi = fCosPhi;
312 0 : float secPhi = fSecPhi;
313 0 : float dzds = fDzDs;
314 0 : float qpt = fQPt;
315 :
316 : // Rotate the coordinate system in XY on the angle alpha
317 : {
318 : float sP = sinPhi, cP = cosPhi;
319 0 : cosPhi = cP * cosAlpha + sP * sinAlpha;
320 0 : sinPhi = -cP * sinAlpha + sP * cosAlpha;
321 :
322 0 : if ( CAMath::Abs( sinPhi ) > .999 || CAMath::Abs( cP ) < 1.e-2 ) return 0;
323 :
324 0 : secPhi = 1./cosPhi;
325 0 : float j0 = cP *secPhi;
326 0 : float j2 = cosPhi / cP;
327 0 : x = fX*cosAlpha + fY*sinAlpha ;
328 0 : y = -fX*sinAlpha + fY*cosAlpha ;
329 :
330 0 : c00 *= j0 * j0;
331 0 : c40 *= j0;
332 :
333 0 : c22 *= j2 * j2;
334 0 : c42 *= j2;
335 0 : if( cosPhi < 0.f ){ // rotate to 180'
336 0 : cosPhi = -cosPhi;
337 0 : secPhi = -secPhi;
338 0 : sinPhi = -sinPhi;
339 0 : dzds = -dzds;
340 0 : qpt = -qpt;
341 0 : c20 = -c20;
342 0 : c31 = -c31;
343 0 : c40 = -c40;
344 0 : }
345 0 : }
346 :
347 0 : Bz = -Bz;
348 : float ex = cosPhi;
349 : float ey = sinPhi;
350 0 : float k = qpt*Bz;
351 0 : float dx = newX - x;
352 0 : float ey1 = k*dx + ey;
353 :
354 0 : if( fabs( ey1 ) > maxSinPhi ) return 0;
355 :
356 0 : float ex1 = sqrt( 1.f - ey1 * ey1 );
357 :
358 0 : float dxBz = dx * Bz;
359 :
360 0 : float ss = ey + ey1;
361 0 : float cc = ex + ex1;
362 0 : float dxcci = dx / cc;
363 0 : float norm2 = 1.f + ey*ey1 + ex*ex1;
364 :
365 0 : float dy = dxcci * ss;
366 :
367 : float dS;
368 : {
369 0 : float dl = dxcci * sqrt( norm2 + norm2 );
370 0 : float dSin = 0.5f*k*dl;
371 0 : float a = dSin*dSin;
372 : const float k2 = 1.f/6.f;
373 : const float k4 = 3.f/40.f;
374 : //const float k6 = 5.f/112.f;
375 0 : dS = dl + dl*a*(k2 + a*(k4 ));//+ k6*a) );
376 : }
377 :
378 0 : float ex1i = 1.f/ex1;
379 0 : float dz = dS * dzds;
380 :
381 0 : float hh = dxcci*ex1i*norm2;
382 0 : float h2 = hh * secPhi;
383 0 : float h4 = Bz*dxcci*hh;
384 :
385 0 : float c20ph4c42 = c20 + h4*c42;
386 0 : float h2c22 = h2*c22;
387 0 : float h4c44 = h4*c44;
388 0 : float n7 = c31 + dS*c33;
389 :
390 0 : b.SetPar(0, y + dy );
391 0 : b.SetPar(1, z + dz );
392 0 : b.SetPar(2, ey1 );
393 0 : b.SetPar(3, dzds);
394 0 : b.SetPar(4, qpt);
395 :
396 0 : b.SetCov(0, c00 + h2*h2c22 + h4*h4c44 + 2.f*( h2*c20ph4c42 + h4*c40 ));
397 0 : b.SetCov(1, c11 + dS*(c31 + n7) );
398 0 : b.SetCov(2, c22 + dxBz*( c42 + c42 + dxBz*c44 ));
399 0 : b.SetCov(3, c33);
400 0 : b.SetCov(4, c44);
401 0 : b.SetCovD(0, c20ph4c42 + h2c22 + dxBz*(c40 + h2*c42 + h4c44) );
402 0 : b.SetCovD(1, n7 );
403 :
404 : return 1;
405 0 : }
|
