HLT/TPCLib/tracking-ca/AliHLTTPCCATrackParam.cxx

   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   float dx = GetX() - t.GetX();
  43   float dy = GetY() - t.GetY();
  44   float dz = GetZ() - t.GetZ();
  45   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   float dx = GetX() - t.GetX();
  53   float dz = GetZ() - t.GetZ();
  54   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   float k  = GetKappa( Bz );
  63   float ex = GetCosPhi();
  64   float ey = GetSinPhi();
  65   x -= GetX();
  66   y -= GetY();
  67   float dS = x * ex + y * ey;
  68   if ( CAMath::Abs( k ) > 1.e-4 ) dS = CAMath::ATan2( k * dS, 1 + k * ( x * ey - y * ex ) ) / k;
  69   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   float x0 = GetX();
  79   float y0 = GetY();
  80   float k  = GetKappa( Bz );
  81   float ex = GetCosPhi();
  82   float ey = GetSinPhi();
  83   float dx = x - x0;
  84   float dy = y - y0;
  85   float ax = dx * k + ey;
  86   float ay = dy * k - ex;
  87   float a = sqrt( ax * ax + ay * ay );
  88   xp = x0 + ( dx - ey * ( ( dx * dx + dy * dy ) * k - 2 * ( -dx * ey + dy * ex ) ) / ( a + 1 ) ) / a;
  89   yp = y0 + ( dy + ex * ( ( dx * dx + dy * dy ) * k - 2 * ( -dx * ey + dy * ex ) ) / ( a + 1 ) ) / a;
  90   float s = GetS( x, y, Bz );
  91   zp = GetZ() + GetDzDs() * s;
  92   if ( CAMath::Abs( k ) > 1.e-2 ) {
  93     float dZ = CAMath::Abs( GetDzDs() * CAMath::TwoPi() / k );
  94     if ( dZ > .1 ) {
  95       zp += CAMath::Nint( ( z - zp ) / dZ ) * dZ;
  96     }
  97   }
  98 }
  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   float ex = t0.CosPhi();
 115   float ey = t0.SinPhi();
 116   float k  =-t0.QPt() * Bz;
 117   float dx = x - X();
 118
 119   float ey1 = k * dx + ey;
 120   float ex1;
 121
 122   // check for intersection with X=x
 123
 124   if ( CAMath::Abs( ey1 ) > maxSinPhi ) return 0;
 125
 126   ex1 = CAMath::Sqrt( 1 - ey1 * ey1 );
 127   if ( ex < 0 ) ex1 = -ex1;
 128
 129   float dx2 = dx * dx;
 130   float ss = ey + ey1;
 131   float cc = ex + ex1;
 132
 133   if ( CAMath::Abs( cc ) < 1.e-4 || CAMath::Abs( ex ) < 1.e-4 || CAMath::Abs( ex1 ) < 1.e-4 ) return 0;
 134
 135   float tg = ss / cc; // tan((phi1+phi)/2)
 136
 137   float dy = dx * tg;
 138   float dl = dx * CAMath::Sqrt( 1 + tg * tg );
 139
 140   if ( cc < 0 ) dl = -dl;
 141   float dSin = dl * k / 2;
 142   if ( dSin > 1 ) dSin = 1;
 143   if ( dSin < -1 ) dSin = -1;
 144   float dS = ( CAMath::Abs( k ) > 1.e-4 )  ? ( 2 * CAMath::ASin( dSin ) / k ) : dl;
 145   float dz = dS * t0.DzDs();
 146
 147   if ( DL ) *DL = -dS * CAMath::Sqrt( 1 + t0.DzDs() * t0.DzDs() );
 148
 149   float cci = 1. / cc;
 150   float exi = 1. / ex;
 151   float ex1i = 1. / ex1;
 152
 153   float d[5] = { 0,
 154                  0,
 155                  GetPar(2) - t0.SinPhi(),
 156                  GetPar(3) - t0.DzDs(),
 157                  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   float h2 = dx * ( 1 + ey * ey1 + ex * ex1 ) * exi * ex1i * cci;
 167   float h4 = dx2 * ( cc + ss * ey1 * ex1i ) * cci * cci * (-Bz);
 168   float dxBz = dx * (-Bz);
 169
 170   t0.SetCosPhi( ex1 );
 171   t0.SetSinPhi( ey1 );
 172
 173   SetX(X() + dx);
 174   SetPar(0, Y() + dy     + h2 * d[2]           +   h4 * d[4]);
 175   SetPar(1, Z() + dz               + dS * d[3]);
 176   SetPar(2, t0.SinPhi() +     d[2]           + dxBz * d[4]);
 177
 178   float c00 = fC[0];
 179   float c10 = fC[1];
 180   float c11 = fC[2];
 181   float c20 = fC[3];
 182   float c21 = fC[4];
 183   float c22 = fC[5];
 184   float c30 = fC[6];
 185   float c31 = fC[7];
 186   float c32 = fC[8];
 187   float c33 = fC[9];
 188   float c40 = fC[10];
 189   float c41 = fC[11];
 190   float c42 = fC[12];
 191   float c43 = fC[13];
 192   float c44 = fC[14];
 193
 194   fC[0] = ( c00  + h2 * h2 * c22 + h4 * h4 * c44
 195             + 2 * ( h2 * c20 + h4 * c40 + h2 * h4 * c42 )  );
 196
 197   fC[1] = c10 + h2 * c21 + h4 * c41 + dS * ( c30 + h2 * c32 + h4 * c43 );
 198   fC[2] = c11 + 2 * dS * c31 + dS * dS * c33;
 199
 200   fC[3] = c20 + h2 * c22 + h4 * c42 + dxBz * ( c40 + h2 * c42 + h4 * c44 );
 201   fC[4] = c21 + dS * c32 + dxBz * ( c41 + dS * c43 );
 202   fC[5] = c22 + 2 * dxBz * c42 + dxBz * dxBz * c44;
 203
 204   fC[6] = c30 + h2 * c32 + h4 * c43;
 205   fC[7] = c31 + dS * c33;
 206   fC[8] = c32 + dxBz * c43;
 207   fC[9] = c33;
 208
 209   fC[10] = c40 + h2 * c42 + h4 * c44;
 210   fC[11] = c41 + dS * c43;
 211   fC[12] = c42 + dxBz * c44;
 212   fC[13] = c43;
 213   fC[14] = c44;
 214
 215   return 1;
 216 }
 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   float dx = x - X();
 231
 232   if ( CAMath::Abs( ex ) < 1.e-4 ) return 0;
 233   float exi = 1. / ex;
 234
 235   float dxBz = dx * (-Bz);
 236   float dS = dx * exi;
 237   float h2 = dS * exi * exi;
 238   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   float sinPhi = SinPhi() + dxBz * QPt();
 247   if ( maxSinPhi > 0 && CAMath::Abs( sinPhi ) > maxSinPhi ) return 0;
 248
 249   SetX(X() + dx);
 250   SetPar(0, GetPar(0) + dS * ey + h2 * ( SinPhi() - ey )  +   h4 * QPt());
 251   SetPar(1, GetPar(1) + dS * DzDs());
 252   SetPar(2, sinPhi);
 253
 254
 255   float c00 = fC[0];
 256   float c10 = fC[1];
 257   float c11 = fC[2];
 258   float c20 = fC[3];
 259   float c21 = fC[4];
 260   float c22 = fC[5];
 261   float c30 = fC[6];
 262   float c31 = fC[7];
 263   float c32 = fC[8];
 264   float c33 = fC[9];
 265   float c40 = fC[10];
 266   float c41 = fC[11];
 267   float c42 = fC[12];
 268   float c43 = fC[13];
 269   float c44 = fC[14];
 270
 271
 272   fC[0] = ( c00  + h2 * h2 * c22 + h4 * h4 * c44
 273             + 2 * ( h2 * c20 + h4 * c40 + h2 * h4 * c42 )  );
 274
 275   fC[1] = c10 + h2 * c21 + h4 * c41 + dS * ( c30 + h2 * c32 + h4 * c43 );
 276   fC[2] = c11 + 2 * dS * c31 + dS * dS * c33;
 277
 278   fC[3] = c20 + h2 * c22 + h4 * c42 + dxBz * ( c40 + h2 * c42 + h4 * c44 );
 279   fC[4] = c21 + dS * c32 + dxBz * ( c41 + dS * c43 );
 280   fC[5] = c22 + 2 * dxBz * c42 + dxBz * dxBz * c44;
 281
 282   fC[6] = c30 + h2 * c32 + h4 * c43;
 283   fC[7] = c31 + dS * c33;
 284   fC[8] = c32 + dxBz * c43;
 285   fC[9] = c33;
 286
 287   fC[10] = c40 + h2 * c42 + h4 * c44;
 288   fC[11] = c41 + dS * c43;
 289   fC[12] = c42 + dxBz * c44;
 290   fC[13] = c43;
 291   fC[14] = c44;
 292
 293   return 1;
 294 }
 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   AliHLTTPCCATrackLinearisation t0( *this );
 306
 307   return TransportToX( x, t0, Bz, maxSinPhi );
 308 }
 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   float dl;
 320
 321   if ( !TransportToX( x, t0, Bz,  maxSinPhi, &dl ) ) return 0;
 322
 323   CorrectForMeanMaterial( dl*kRhoOverRadLen, dl*kRho, par );
 324   return 1;
 325 }
 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   AliHLTTPCCATrackLinearisation t0( *this );
 333   return TransportToXWithMaterial( x, t0, par, Bz, maxSinPhi );
 334 }
 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   AliHLTTPCCATrackFitParam par;
 341   CalculateFitParameters( par );
 342   return TransportToXWithMaterial( x, par, Bz, maxSinPhi );
 343 }
 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   const float x0  = kp1 * 2.303;
 376   const float x1  = kp2 * 2.303;
 377   const float mI  = kp3;
 378   const float mZA = kp4;
 379   const float maxT = 2 * me * bg2;    // neglecting the electron mass
 380
 381   //*** Density effect
 382   float d2 = 0.;
 383   const float x = 0.5 * AliHLTTPCCAMath::Log( bg2 );
 384   const float lhwI = AliHLTTPCCAMath::Log( 28.816 * 1e-9 * AliHLTTPCCAMath::Sqrt( rho * mZA ) / mI );
 385   if ( x > x1 ) {
 386     d2 = lhwI + x - 0.5;
 387   } else if ( x > x0 ) {
 388     const float r = ( x1 - x ) / ( x1 - x0 );
 389     d2 = lhwI + x - 0.5 + ( 0.5 - lhwI - x0 ) * r * r * r;
 390   }
 391
 392   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   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   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   if ( beta2 >= 1 ) return 0;
 436
 437   if ( beta2 / ( 1 - beta2 ) > 3.5*3.5 )
 438     return 0.153e-3 / beta2*( log( 3.5*5940 ) + 0.5*log( beta2 / ( 1 - beta2 ) ) - beta2 );
 439   return 0.153e-3 / beta2*( log( 5940*beta2 / ( 1 - beta2 ) ) - beta2 );
 440 }
 441
 442
 443 MEM_CLASS_PRE() GPUdi() void MEM_LG(AliHLTTPCCATrackParam)::CalculateFitParameters( AliHLTTPCCATrackFitParam &par, float mass )
 444 {
 445   //*!
 446
 447   float qpt = GetPar(4);
 448   if( fC[14]>=1. ) qpt = 1./0.35;
 449
 450   float p2 = ( 1. + GetPar(3) * GetPar(3) );
 451   float k2 = qpt * qpt;
 452   float mass2 = mass * mass;
 453   float beta2 = p2 / ( p2 + mass2 * k2 );
 454
 455   float pp2 = ( k2 > 1.e-8 ) ? p2 / k2 : 10000; // impuls 2
 456
 457   //par.fBethe = BetheBlochGas( pp2/mass2);
 458   par.fBethe = ApproximateBetheBloch( pp2 / mass2 );
 459   par.fE = CAMath::Sqrt( pp2 + mass2 );
 460   par.fTheta2 = 14.1 * 14.1 / ( beta2 * pp2 * 1e6 );
 461   par.fEP2 = par.fE / pp2;
 462
 463   // Approximate energy loss fluctuation (M.Ivanov)
 464
 465   const float knst = 0.07; // To be tuned.
 466   par.fSigmadE2 = knst * par.fEP2 * qpt;
 467   par.fSigmadE2 = par.fSigmadE2 * par.fSigmadE2;
 468
 469   par.fK22 = ( 1. + GetPar(3) * GetPar(3) );
 470   par.fK33 = par.fK22 * par.fK22;
 471   par.fK43 = 0;
 472   par.fK44 = GetPar(3) * GetPar(3) * k2;
 473
 474 }
 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   float &fC22 = fC[5];
 486   float &fC33 = fC[9];
 487   float &fC40 = fC[10];
 488   float &fC41 = fC[11];
 489   float &fC42 = fC[12];
 490   float &fC43 = fC[13];
 491   float &fC44 = fC[14];
 492
 493   //Energy losses************************
 494
 495   float dE = par.fBethe * xTimesRho;
 496   if ( CAMath::Abs( dE ) > 0.3*par.fE ) return 0; //30% energy loss is too much!
 497   float corr = ( 1. - par.fEP2 * dE );
 498   if ( corr < 0.3 || corr > 1.3 ) return 0;
 499
 500   SetPar(4, GetPar(4) * corr);
 501   fC40 *= corr;
 502   fC41 *= corr;
 503   fC42 *= corr;
 504   fC43 *= corr;
 505   fC44 *= corr * corr;
 506   fC44 += par.fSigmadE2 * CAMath::Abs( dE );
 507
 508   //Multiple scattering******************
 509
 510   float theta2 = par.fTheta2 * CAMath::Abs( xOverX0 );
 511   fC22 += theta2 * par.fK22 * (1.-GetPar(2))*(1.+GetPar(2));
 512   fC33 += theta2 * par.fK33;
 513   fC43 += theta2 * par.fK43;
 514   fC44 += theta2 * par.fK44;
 515
 516   return 1;
 517 }
 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   float cA = CAMath::Cos( alpha );
 530   float sA = CAMath::Sin( alpha );
 531   float x = X(), y = Y(), sP = SinPhi(), cP = GetCosPhi();
 532   float cosPhi = cP * cA + sP * sA;
 533   float sinPhi = -cP * sA + sP * cA;
 534
 535   if ( CAMath::Abs( sinPhi ) > maxSinPhi || CAMath::Abs( cosPhi ) < 1.e-2 || CAMath::Abs( cP ) < 1.e-2  ) return 0;
 536
 537   float j0 = cP / cosPhi;
 538   float j2 = cosPhi / cP;
 539
 540   SetX( x*cA +  y*sA );
 541   SetY( -x*sA +  y*cA );
 542   SetSignCosPhi( cosPhi );
 543   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   fC[0] *= j0 * j0;
 554   fC[1] *= j0;
 555   fC[3] *= j0;
 556   fC[6] *= j0;
 557   fC[10] *= j0;
 558
 559   fC[3] *= j2;
 560   fC[4] *= j2;
 561   fC[5] *= j2 * j2;
 562   fC[8] *= j2;
 563   fC[12] *= j2;
 564
 565         if (cosPhi < 0)
 566         {
 567                 SetSinPhi(-SinPhi());
 568                 SetDzDs(-DzDs());
 569                 SetQPt(-QPt());
 570                 fC[3] = - fC[3];
 571                 fC[4] = - fC[4];
 572                 fC[6] = - fC[6];
 573                 fC[7] = - fC[7];
 574                 fC[10] = -fC[10];
 575                 fC[11] = -fC[11];
 576         }
 577
 578   //cout<<"      "<<fC[0]<<" "<<fC[1]<<" "<<fC[6]<<" "<<fC[10]<<" "<<fC[4]<<" "<<fC[5]<<" "<<fC[8]<<" "<<fC[12]<<endl;
 579   return 1;
 580 }
 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   float cA = CAMath::Cos( alpha );
 587   float sA = CAMath::Sin( alpha );
 588   float x0 = X(), y0 = Y(), sP = t0.SinPhi(), cP = t0.CosPhi();
 589   float cosPhi = cP * cA + sP * sA;
 590   float sinPhi = -cP * sA + sP * cA;
 591
 592   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   float j0 = cP / cosPhi;
 601   float j2 = cosPhi / cP;
 602   float d[2] = {Y() - y0, SinPhi() - sP};
 603
 604   SetX( x0*cA +  y0*sA );
 605   SetY( -x0*sA +  y0*cA + j0*d[0] );
 606   t0.SetCosPhi( cosPhi );
 607   t0.SetSinPhi( sinPhi );
 608
 609   SetSinPhi( sinPhi + j2*d[1] );
 610
 611   fC[0] *= j0 * j0;
 612   fC[1] *= j0;
 613   fC[3] *= j0;
 614   fC[6] *= j0;
 615   fC[10] *= j0;
 616
 617   fC[3] *= j2;
 618   fC[4] *= j2;
 619   fC[5] *= j2 * j2;
 620   fC[8] *= j2;
 621   fC[12] *= j2;
 622
 623   return 1;
 624 }
 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   c00 = fC[ 0],
 632         c11 = fC[ 2],
 633               c20 = fC[ 3],
 634                     c31 = fC[ 7],
 635                           c40 = fC[10];
 636
 637   err2Y += c00;
 638   err2Z += c11;
 639
 640   float
 641   z0 = y - GetPar(0),
 642        z1 = z - GetPar(1);
 643
 644   if ( err2Y < 1.e-8 || err2Z < 1.e-8 ) return 0;
 645
 646   float mS0 = 1. / err2Y;
 647   float mS2 = 1. / err2Z;
 648
 649   // K = CHtS
 650
 651   float k00, k11, k20, k31, k40;
 652
 653   k00 = c00 * mS0;
 654   k20 = c20 * mS0;
 655   k40 = c40 * mS0;
 656
 657   k11 = c11 * mS2;
 658   k31 = c31 * mS2;
 659
 660   float sinPhi = GetPar(2) + k20 * z0  ;
 661
 662   if ( maxSinPhi > 0 && CAMath::Abs( sinPhi ) >= maxSinPhi ) return 0;
 663
 664   fNDF  += 2;
 665   fChi2 += mS0 * z0 * z0 + mS2 * z1 * z1 ;
 666
 667   SetPar(0, GetPar(0) + k00 * z0);
 668   SetPar(1, GetPar(1) + k11 * z1);
 669   SetPar(2, sinPhi);
 670   SetPar(3, GetPar(3) + k31 * z1);
 671   SetPar(4, GetPar(4) + k40 * z0);
 672
 673   fC[ 0] -= k00 * c00 ;
 674   fC[ 3] -= k20 * c00 ;
 675   fC[ 5] -= k20 * c20 ;
 676   fC[10] -= k40 * c00 ;
 677   fC[12] -= k40 * c20 ;
 678   fC[14] -= k40 * c40 ;
 679
 680   fC[ 2] -= k11 * c11 ;
 681   fC[ 7] -= k31 * c11 ;
 682   fC[ 9] -= k31 * c31 ;
 683
 684   return 1;
 685 }
 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   bool ok = AliHLTTPCCAMath::Finite( GetX() ) && AliHLTTPCCAMath::Finite( fSignCosPhi ) && AliHLTTPCCAMath::Finite( fChi2 ) && AliHLTTPCCAMath::Finite( fNDF );
 692
 693   const float *c = Cov();
 694   for ( int i = 0; i < 15; i++ ) ok = ok && AliHLTTPCCAMath::Finite( c[i] );
 695   for ( int i = 0; i < 5; i++ ) ok = ok && AliHLTTPCCAMath::Finite( Par()[i] );
 696
 697   if ( c[0] <= 0 || c[2] <= 0 || c[5] <= 0 || c[9] <= 0 || c[14] <= 0 ) ok = 0;
 698   if ( c[0] > 5. || c[2] > 5. || c[5] > 2. || c[9] > 2
 699        //|| ( CAMath::Abs( QPt() ) > 1.e-2 && c[14] > 2. )
 700        ) ok = 0;
 701
 702   if ( CAMath::Abs( SinPhi() ) > .99 ) ok = 0;
 703   if ( CAMath::Abs( QPt() ) > 1. / 0.05 ) ok = 0;
 704   if( ok ){
 705     ok = ok
 706       && ( c[1]*c[1]<=c[2]*c[0] )
 707       && ( c[3]*c[3]<=c[5]*c[0] )
 708       && ( c[4]*c[4]<=c[5]*c[2] )
 709       && ( c[6]*c[6]<=c[9]*c[0] )
 710       && ( c[7]*c[7]<=c[9]*c[2] )
 711       && ( c[8]*c[8]<=c[9]*c[5] )
 712       && ( c[10]*c[10]<=c[14]*c[0] )
 713       && ( c[11]*c[11]<=c[14]*c[2] )
 714       && ( c[12]*c[12]<=c[14]*c[5] )
 715       && ( c[13]*c[13]<=c[14]*c[9] );
 716   }
 717   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   std::cout << "track: x=" << GetX() << " c=" << GetSignCosPhi() << ", P= " << GetY() << " " << GetZ() << " " << GetSinPhi() << " " << GetDzDs() << " " << GetQPt() << std::endl;
 731   std::cout << "errs2: " << GetErr2Y() << " " << GetErr2Z() << " " << GetErr2SinPhi() << " " << GetErr2DzDs() << " " << GetErr2QPt() << std::endl;
 732 #endif
 733 }
 734