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. *
6 // Primary Authors: Sergey Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de> *
7 // for The ALICE HLT Project. *
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. *
17 //***************************************************************************
19 #include "AliHLTTPCGMSliceTrack.h"
20 #include "AliHLTTPCCAMath.h"
21 #include "AliHLTTPCGMBorderTrack.h"
22 #include "AliHLTTPCGMTrackLinearisation.h"
23 #include "Riostream.h"
24 #include "AliHLTTPCCAParam.h"
29 bool AliHLTTPCGMSliceTrack::FilterErrors( AliHLTTPCCAParam ¶m, float maxSinPhi )
31 float lastX = fOrigTrack->Cluster(fOrigTrack->NClusters()-1 ).GetX();
34 const float *cy = param.GetParamS0Par(0,0);
35 const float *cz = param.GetParamS0Par(1,0);
37 float bz = -param.ConstBz();
38 const float kZLength = 250.f - 0.275f;
40 float trDzDs2 = fDzDs*fDzDs;
42 float dx = .33*(lastX - fX);
45 float kdx205 = 2.f+kdx*kdx*0.5f;
48 float secPhi2 = fSecPhi*fSecPhi;
49 float zz = fabs( kZLength - fabs(fZ) );
51 float angleY2 = secPhi2 - 1.f;
52 float angleZ2 = trDzDs2 * secPhi2 ;
54 float cy0 = cy[0] + cy[1]*zz + cy[3]*zz2;
55 float cy1 = cy[2] + cy[5]*zz;
57 float cz0 = cz[0] + cz[1]*zz + cz[3]*zz2;
58 float cz1 = cz[2] + cz[5]*zz;
61 fC0 = fabs( cy0 + angleY2 * ( cy1 + angleY2*cy2 ) );
62 fC2 = fabs( cz0 + angleZ2 * ( cz1 + angleZ2*cz2 ) );
73 for( int iStep=0; iStep<N; iStep++ ){
82 if( fabs( ey1 ) > maxSinPhi ) return 0;
85 float ex1 = sqrt(1.f - ey1*ey1);
88 float dxcci = dx / cc;
90 float dy = dxcci * ss;
91 float norm2 = 1.f + ey*ey1 + ex*ex1;
92 float dl = dxcci * sqrt( norm2 + norm2 );
96 float dSin = 0.5f*k*dl;
98 const float k2 = 1.f/6.f;
99 const float k4 = 3.f/40.f;
100 dS = dl + dl*a*(k2 + a*(k4 ));//+ k6*a) );
103 float dz = dS * fDzDs;
107 float secPhi2 = ex1i*ex1i;
108 float zz = fabs( kZLength - fabs(z) );
110 float angleY2 = secPhi2 - 1.f;
111 float angleZ2 = trDzDs2 * secPhi2 ;
113 float cy0 = cy[0] + cy[1]*zz + cy[3]*zz2;
114 float cy1 = cy[2] + cy[5]*zz;
116 float cz0 = cz[0] + cz[1]*zz + cz[3]*zz2;
117 float cz1 = cz[2] + cz[5]*zz;
120 err2Y = fabs( cy0 + angleY2 * ( cy1 + angleY2*cy2 ) );
121 err2Z = fabs( cz0 + angleZ2 * ( cz1 + angleZ2*cz2 ) );
124 float hh = kdx205 * dxcci*ex1i;
125 float h2 = hh * fSecPhi;
134 float h4 = bz*dxcci*hh;
144 float c20ph4c42 = c20 + h4*c42;
145 float h2c22 = h2*c22;
146 float h4c44 = h4*c44;
147 float n7 = c31 + dS*c33;
148 float n10 = c40 + h2*c42 + h4c44;
149 float n12 = c42 + dxBz*c44;
152 fC0+= h2*h2c22 + h4*h4c44 + 2.f*( h2*c20ph4c42 + h4*c40 );
154 fC3 = c20ph4c42 + h2c22 + dxBz*n10;
157 fC5 = c22 + dxBz*( c42 + n12 );
163 } // end transport block
175 float mS0 = 1.f/(err2Y + c00);
176 float mS2 = 1.f/(err2Z + c11);
180 float k00, k11, k20, k31, k40;
201 //* Check that the track parameters and covariance matrix are reasonable
205 const float *c = &fX;
206 for ( int i = 0; i < 17; i++ ) ok = ok && finite( c[i] );
208 if ( fC0 <= 0.f || fC2 <= 0.f || fC5 <= 0.f || fC9 <= 0.f || fC14 <= 0.f
209 || fC0 > 5.f || fC2 > 5.f || fC5 > 2.f || fC9 > 2.f ) ok = 0;
213 && ( fC3*fC3<=fC5*fC0 )
214 && ( fC7*fC7<=fC9*fC2 )
215 && ( fC10*fC10<=fC14*fC0 )
216 && ( fC12*fC12<=fC14*fC5 );
224 bool AliHLTTPCGMSliceTrack::TransportToX( float x, float Bz, AliHLTTPCGMBorderTrack &b, float maxSinPhi ) const
231 float ey1 = k*dx + ey;
233 if( fabs( ey1 ) > maxSinPhi ) return 0;
235 float ex1 = sqrt( 1.f - ey1 * ey1 );
236 float dxBz = dx * Bz;
240 float dxcci = dx / cc;
241 float norm2 = 1.f + ey*ey1 + ex*ex1;
243 float dy = dxcci * ss;
247 float dl = dxcci * sqrt( norm2 + norm2 );
248 float dSin = 0.5f*k*dl;
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 dS = dl + dl*a*(k2 + a*(k4 ));//+ k6*a) );
256 float ex1i = 1.f/ex1;
257 float dz = dS * fDzDs;
259 float hh = dxcci*ex1i*norm2;
260 float h2 = hh *fSecPhi;
261 float h4 = Bz*dxcci*hh;
271 float c20ph4c42 = c20 + h4*c42;
272 float h2c22 = h2*c22;
273 float h4c44 = h4*c44;
274 float n7 = c31 + dS*c33;
276 b.SetPar(0, fY + dy );
277 b.SetPar(1, fZ + dz );
282 b.SetCov(0, fC0 + h2*h2c22 + h4*h4c44 + 2.f*( h2*c20ph4c42 + h4*c40 ));
283 b.SetCov(1, fC2+ dS*(c31 + n7) );
284 b.SetCov(2, c22 + dxBz*( c42 + c42 + dxBz*c44 ));
287 b.SetCovD(0, c20ph4c42 + h2c22 + dxBz*(c40 + h2*c42 + h4c44) );
294 bool AliHLTTPCGMSliceTrack::TransportToXAlpha( float newX, float sinAlpha, float cosAlpha, float Bz, AliHLTTPCGMBorderTrack &b, float maxSinPhi ) const
310 float sinPhi = fSinPhi;
311 float cosPhi = fCosPhi;
312 float secPhi = fSecPhi;
316 // Rotate the coordinate system in XY on the angle alpha
318 float sP = sinPhi, cP = cosPhi;
319 cosPhi = cP * cosAlpha + sP * sinAlpha;
320 sinPhi = -cP * sinAlpha + sP * cosAlpha;
322 if ( CAMath::Abs( sinPhi ) > .999 || CAMath::Abs( cP ) < 1.e-2 ) return 0;
325 float j0 = cP *secPhi;
326 float j2 = cosPhi / cP;
327 x = fX*cosAlpha + fY*sinAlpha ;
328 y = -fX*sinAlpha + fY*cosAlpha ;
335 if( cosPhi < 0.f ){ // rotate to 180'
352 float ey1 = k*dx + ey;
354 if( fabs( ey1 ) > maxSinPhi ) return 0;
356 float ex1 = sqrt( 1.f - ey1 * ey1 );
358 float dxBz = dx * Bz;
362 float dxcci = dx / cc;
363 float norm2 = 1.f + ey*ey1 + ex*ex1;
365 float dy = dxcci * ss;
369 float dl = dxcci * sqrt( norm2 + norm2 );
370 float dSin = 0.5f*k*dl;
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 dS = dl + dl*a*(k2 + a*(k4 ));//+ k6*a) );
378 float ex1i = 1.f/ex1;
379 float dz = dS * dzds;
381 float hh = dxcci*ex1i*norm2;
382 float h2 = hh * secPhi;
383 float h4 = Bz*dxcci*hh;
385 float c20ph4c42 = c20 + h4*c42;
386 float h2c22 = h2*c22;
387 float h4c44 = h4*c44;
388 float n7 = c31 + dS*c33;
390 b.SetPar(0, y + dy );
391 b.SetPar(1, z + dz );
396 b.SetCov(0, c00 + h2*h2c22 + h4*h4c44 + 2.f*( h2*c20ph4c42 + h4*c40 ));
397 b.SetCov(1, c11 + dS*(c31 + n7) );
398 b.SetCov(2, c22 + dxBz*( c42 + c42 + dxBz*c44 ));
401 b.SetCovD(0, c20ph4c42 + h2c22 + dxBz*(c40 + h2*c42 + h4c44) );