[u/mrichter/AliRoot.git] / HLT / TPCLib / merger-ca / AliHLTTPCGMSliceTrack.cxx

// $Id: AliHLTTPCGMSliceTrack.cxx 41769 2010-06-16 13:58:00Z sgorbuno $
// **************************************************************************
// This file is property of and copyright by the ALICE HLT Project          *
// ALICE Experiment at CERN, All rights reserved.                           *
//                                                                          *
// Primary Authors: Sergey Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de> *
//                  for The ALICE HLT Project.                              *
//                                                                          *
// Permission to use, copy, modify and distribute this software and its     *
// documentation strictly for non-commercial purposes is hereby granted     *
// without fee, provided that the above copyright notice appears in all     *
// copies and that both the copyright notice and this permission notice     *
// appear in the supporting documentation. The authors make no claims       *
// about the suitability of this software for any purpose. It is            *
// provided "as is" without express or implied warranty.                    *
//                                                                          *
//***************************************************************************

#include "AliHLTTPCGMSliceTrack.h"
#include "AliHLTTPCCAMath.h"
#include "AliHLTTPCGMBorderTrack.h"
#include "AliHLTTPCGMTrackLinearisation.h"
#include "Riostream.h"
#include "AliHLTTPCCAParam.h"
#include <cmath>


bool AliHLTTPCGMSliceTrack::FilterErrors( AliHLTTPCCAParam &param, float maxSinPhi )
{
  float lastX = fOrigTrack->Cluster(fOrigTrack->NClusters()-1 ).GetX();

  const int N = 3;
  const float *cy = param.GetParamS0Par(0,0);
  const float *cz = param.GetParamS0Par(1,0);

  float bz = -param.ConstBz();
  const float kZLength = 250.f - 0.275f;

  float trDzDs2 = fDzDs*fDzDs;
  float k  = fQPt*bz;               
  float dx = .33*(lastX - fX);  
  float kdx = k*dx;
  float dxBz = dx * bz;
  float kdx205 = 2.f+kdx*kdx*0.5f;

  {
    float secPhi2 = fSecPhi*fSecPhi;
    float zz = fabs( kZLength - fabs(fZ) );	
    float zz2 = zz*zz;
    float angleY2 = secPhi2 - 1.f; 
    float angleZ2 = trDzDs2 * secPhi2 ;
    
    float cy0 = cy[0] + cy[1]*zz + cy[3]*zz2;
    float cy1 = cy[2] + cy[5]*zz;
    float cy2 = cy[4];
    float cz0 = cz[0] + cz[1]*zz + cz[3]*zz2;
    float cz1 = cz[2] + cz[5]*zz;
    float cz2 = cz[4];
    
    fC0 = fabs( cy0 + angleY2 * ( cy1 + angleY2*cy2 ) );
    fC2 = fabs( cz0 + angleZ2 * ( cz1 + angleZ2*cz2 ) );    

    fC3 = 0;
    fC5 = 1;
    fC7 = 0;
    fC9 = 1;
    fC10 = 0;
    fC12 = 0;
    fC14 = 10;
  }

  for( int iStep=0; iStep<N; iStep++ ){    

    float err2Y, err2Z;

    { // transport block
 
      float ex = fCosPhi; 
      float ey = fSinPhi;
      float ey1 = kdx + ey;
      if( fabs( ey1 ) > maxSinPhi ) return 0;

      float ss = ey + ey1;      
      float ex1 = sqrt(1.f - ey1*ey1);
          
      float cc = ex + ex1;  
      float dxcci = dx / cc;
      
      float dy = dxcci * ss;
      float norm2 = 1.f + ey*ey1 + ex*ex1;
      float dl = dxcci * sqrt( norm2 + norm2 );
     
      float dS;   
      {
	float dSin = 0.5f*k*dl;
	float a = dSin*dSin;
	const float k2 = 1.f/6.f;
	const float k4 = 3.f/40.f;
	dS = dl + dl*a*(k2 + a*(k4 ));//+ k6*a) );
      }
 
      float dz = dS * fDzDs;      
      float ex1i =1.f/ex1;
      float z = fZ+ dz;
      {	
	float secPhi2 = ex1i*ex1i;
	float zz = fabs( kZLength - fabs(z) );	
	float zz2 = zz*zz;
	float angleY2 = secPhi2 - 1.f; 
	float angleZ2 = trDzDs2 * secPhi2 ;

	float cy0 = cy[0] + cy[1]*zz + cy[3]*zz2;
	float cy1 = cy[2] + cy[5]*zz;
	float cy2 = cy[4];
	float cz0 = cz[0] + cz[1]*zz + cz[3]*zz2;
	float cz1 = cz[2] + cz[5]*zz;
	float cz2 = cz[4];
	
	err2Y = fabs( cy0 + angleY2 * ( cy1 + angleY2*cy2 ) );
	err2Z = fabs( cz0 + angleZ2 * ( cz1 + angleZ2*cz2 ) );
      }

      float hh = kdx205 * dxcci*ex1i; 
      float h2 = hh * fSecPhi;

      fX+=dx;      
      fY+= dy;
      fZ+= dz;
      fSinPhi = ey1;
      fCosPhi = ex1;
      fSecPhi = ex1i;
    
      float h4 = bz*dxcci*hh;
      
      float c20 = fC3;
      float c22 = fC5;
      float c31 = fC7;
      float c33 = fC9;
      float c40 = fC10;
      float c42 = fC12;
      float c44 = fC14;
      
      float c20ph4c42 =  c20 + h4*c42;
      float h2c22 = h2*c22;
      float h4c44 = h4*c44;
      float n7 = c31 + dS*c33;
      float n10 = c40 + h2*c42 + h4c44;
      float n12 = c42 + dxBz*c44;
      
      
      fC0+= h2*h2c22 + h4*h4c44 + 2.f*( h2*c20ph4c42  + h4*c40 );
      
      fC3 = c20ph4c42 + h2c22  + dxBz*n10;
      fC10 = n10;
      
      fC5 = c22 + dxBz*( c42 + n12 );
      fC12 = n12;
      
      fC2+= dS*(c31 + n7);
      fC7 = n7; 
      
   } // end transport block 


    // Filter block

    float 
      c00 = fC0,
      c11 = fC2,
      c20 = fC3,
      c31 = fC7,
      c40 = fC10;
                    
    float mS0 = 1.f/(err2Y + c00);    
    float mS2 = 1.f/(err2Z + c11);
            
    // K = CHtS
    
    float k00, k11, k20, k31, k40;
    
    k00 = c00 * mS0;
    k20 = c20 * mS0;
    k40 = c40 * mS0;
    
    fC0 -= k00 * c00 ;
    fC5 -= k20 * c20 ;
    fC10 -= k00 * c40 ;
    fC12 -= k40 * c20 ;
    fC3 -= k20 * c00 ;
    fC14 -= k40 * c40 ;
        
    k11 = c11 * mS2;
    k31 = c31 * mS2;
    
    fC7 -= k31 * c11;
    fC2 -= k11 * c11;
    fC9 -= k31 * c31;   
  }

  //* Check that the track parameters and covariance matrix are reasonable

  bool ok = 1;
  
  const float *c = &fX;
  for ( int i = 0; i < 17; i++ ) ok = ok && finite( c[i] );

  if ( fC0 <= 0.f || fC2 <= 0.f || fC5 <= 0.f || fC9 <= 0.f || fC14 <= 0.f 
       || fC0 > 5.f || fC2 > 5.f || fC5 > 2.f || fC9 > 2.f             ) ok = 0;

  if( ok ){
    ok = ok 
      && ( fC3*fC3<=fC5*fC0 )
      && ( fC7*fC7<=fC9*fC2 )
      && ( fC10*fC10<=fC14*fC0 )
      && ( fC12*fC12<=fC14*fC5 );
  }
 
  return ok;
}


bool AliHLTTPCGMSliceTrack::TransportToX( float x, float Bz, AliHLTTPCGMBorderTrack &b, float maxSinPhi ) const 
{
  Bz = -Bz;
  float ex = fCosPhi;
  float ey = fSinPhi;
  float k  = fQPt*Bz;
  float dx = x - fX;
  float ey1 = k*dx + ey;
  
  if( fabs( ey1 ) > maxSinPhi ) return 0;

  float ex1 = sqrt( 1.f - ey1 * ey1 );
  float dxBz = dx * Bz;
    
  float ss = ey + ey1;
  float cc = ex + ex1;  
  float dxcci = dx / cc;
  float norm2 = 1.f + ey*ey1 + ex*ex1;

  float dy = dxcci * ss;

  float dS;    
  {
    float dl = dxcci * sqrt( norm2 + norm2 );
    float dSin = 0.5f*k*dl;
    float a = dSin*dSin;
    const float k2 = 1.f/6.f;
    const float k4 = 3.f/40.f;
    //const float k6 = 5.f/112.f;
    dS = dl + dl*a*(k2 + a*(k4 ));//+ k6*a) );
  }
  
  float ex1i = 1.f/ex1;
  float dz = dS * fDzDs;

  float hh = dxcci*ex1i*norm2; 
  float h2 = hh *fSecPhi;
  float h4 = Bz*dxcci*hh;
  
  float c20 = fC3;
  float c22 = fC5;  
  float c31 = fC7;  
  float c33 = fC9;
  float c40 = fC10;  
  float c42 = fC12;
  float c44 = fC14;

  float c20ph4c42 =  c20 + h4*c42;
  float h2c22 = h2*c22;
  float h4c44 = h4*c44;
  float n7 = c31 + dS*c33;
  
  b.SetPar(0, fY + dy );
  b.SetPar(1, fZ + dz );
  b.SetPar(2, ey1 );
  b.SetPar(3, fDzDs);
  b.SetPar(4, fQPt);

  b.SetCov(0, fC0 + h2*h2c22 + h4*h4c44 + 2.f*( h2*c20ph4c42  + h4*c40 ));
  b.SetCov(1, fC2+ dS*(c31 + n7) );
  b.SetCov(2, c22 + dxBz*( c42 + c42 + dxBz*c44 ));
  b.SetCov(3, c33);
  b.SetCov(4, c44);
  b.SetCovD(0, c20ph4c42 + h2c22  + dxBz*(c40 + h2*c42 + h4c44) );
  b.SetCovD(1, n7 );   
  return 1;
}


bool AliHLTTPCGMSliceTrack::TransportToXAlpha( float newX, float sinAlpha, float cosAlpha, float Bz, AliHLTTPCGMBorderTrack &b, float maxSinPhi ) const 
{
  //* 

  float c00 = fC0;
  float c11 = fC2;
  float c20 = fC3;
  float c22 = fC5;  
  float c31 = fC7;  
  float c33 = fC9;
  float c40 = fC10;  
  float c42 = fC12;
  float c44 = fC14;

  float x,y;
  float z = fZ;
  float sinPhi = fSinPhi;
  float cosPhi = fCosPhi;
  float secPhi = fSecPhi;
  float dzds = fDzDs;
  float qpt = fQPt;

  // Rotate the coordinate system in XY on the angle alpha
  {
    float sP = sinPhi, cP = cosPhi;
    cosPhi =  cP * cosAlpha + sP * sinAlpha;
    sinPhi = -cP * sinAlpha + sP * cosAlpha;
    
    if ( CAMath::Abs( sinPhi ) > .999 || CAMath::Abs( cP ) < 1.e-2  ) return 0;
    
    secPhi = 1./cosPhi;
    float j0 = cP *secPhi;
    float j2 = cosPhi / cP;    
    x =   fX*cosAlpha +  fY*sinAlpha ;
    y =  -fX*sinAlpha +  fY*cosAlpha ;    
    
    c00 *= j0 * j0;
    c40 *= j0;
    
    c22 *= j2 * j2;    
    c42 *= j2;    
    if( cosPhi < 0.f ){ // rotate to 180'
      cosPhi = -cosPhi;
      secPhi = -secPhi;
      sinPhi = -sinPhi;
      dzds = -dzds;
      qpt = -qpt;      
      c20 = -c20;
      c31 = -c31;
      c40 = -c40;
   }
  }

  Bz = -Bz;
  float ex = cosPhi;
  float ey = sinPhi;
  float k  = qpt*Bz;
  float dx = newX - x;
  float ey1 = k*dx + ey;
  
  if( fabs( ey1 ) > maxSinPhi ) return 0;

  float ex1 = sqrt( 1.f - ey1 * ey1 );

  float dxBz = dx * Bz;
    
  float ss = ey + ey1;
  float cc = ex + ex1;  
  float dxcci = dx / cc;
  float norm2 = 1.f + ey*ey1 + ex*ex1;

  float dy = dxcci * ss;

  float dS;    
  {
    float dl = dxcci * sqrt( norm2 + norm2 );
    float dSin = 0.5f*k*dl;
    float a = dSin*dSin;
    const float k2 = 1.f/6.f;
    const float k4 = 3.f/40.f;
    //const float k6 = 5.f/112.f;
    dS = dl + dl*a*(k2 + a*(k4 ));//+ k6*a) );
  }
  
  float ex1i = 1.f/ex1;
  float dz = dS * dzds;

  float hh = dxcci*ex1i*norm2; 
  float h2 = hh * secPhi;
  float h4 = Bz*dxcci*hh;  

  float c20ph4c42 =  c20 + h4*c42;
  float h2c22 = h2*c22;
  float h4c44 = h4*c44;
  float n7 = c31 + dS*c33;
  
  b.SetPar(0, y + dy );
  b.SetPar(1, z + dz );
  b.SetPar(2, ey1 );
  b.SetPar(3, dzds);
  b.SetPar(4, qpt);
  
  b.SetCov(0, c00 + h2*h2c22 + h4*h4c44 + 2.f*( h2*c20ph4c42  + h4*c40 ));
  b.SetCov(1, c11 + dS*(c31 + n7) );
  b.SetCov(2, c22 + dxBz*( c42 + c42 + dxBz*c44 ));
  b.SetCov(3, c33);
  b.SetCov(4, c44);
  b.SetCovD(0, c20ph4c42 + h2c22  + dxBz*(c40 + h2*c42 + h4c44) );
  b.SetCovD(1, n7 ); 

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