HLT/TPCLib/merger-ca/AliHLTTPCGMTrackParam.cxx

   1 // $Id: AliHLTTPCGMTrackParam.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 "AliHLTTPCGMTrackParam.h"
  20 #include "AliHLTTPCCAMath.h"
  21 #include "AliHLTTPCGMTrackLinearisation.h"
  22 #include "AliHLTTPCGMBorderTrack.h"
  23 #include "Riostream.h"
  24 #ifndef HLTCA_STANDALONE
  25 #include "AliExternalTrackParam.h"
  26 #endif
  27 #include "AliHLTTPCCAParam.h"
  28 #include <cmath>
  29
  30 GPUd() void AliHLTTPCGMTrackParam::Fit
  31 (
  32  float* PolinomialFieldBz,
  33  float x[], float y[], float z[], unsigned int rowType[], float alpha[], AliHLTTPCCAParam &param,
  34  int &N,
  35  float &Alpha,
  36  bool UseMeanPt,
  37  float maxSinPhi
  38  ){
  39
  40   const float kRho = 1.025e-3;//0.9e-3;
  41   const float kRadLen = 29.532;//28.94;
  42   const float kRhoOverRadLen = kRho / kRadLen;
  43
  44   AliHLTTPCGMTrackLinearisation t0(*this);
  45
  46   const float kZLength = 250.f - 0.275f;
  47   float trDzDs2 = t0.DzDs()*t0.DzDs();
  48
  49   AliHLTTPCGMTrackFitParam par;
  50   CalculateFitParameters( par, kRhoOverRadLen, kRho, UseMeanPt );
  51
  52   int maxN = N;
  53
  54   bool first = 1;
  55   N = 0;
  56
  57   for( int ihit=0; ihit<maxN; ihit++ ){
  58
  59     float sliceAlpha =  alpha[ihit];
  60
  61     if ( fabs( sliceAlpha - Alpha ) > 1.e-4 ) {
  62       if( !Rotate(  sliceAlpha - Alpha, t0, .999 ) ) break;
  63       Alpha = sliceAlpha;
  64     }
  65
  66     float dL=0;
  67     float bz =  GetBz(x[ihit], y[ihit],z[ihit], PolinomialFieldBz);
  68
  69     float err2Y, err2Z;
  70
  71     { // transport block
  72
  73       bz = -bz;
  74
  75       float ex = t0.CosPhi();
  76
  77       float ey = t0.SinPhi();
  78       float k  = t0.QPt()*bz;
  79       float dx = x[ihit] - X();
  80       float kdx = k*dx;
  81       float ey1 = kdx + ey;
  82
  83       if( fabs( ey1 ) >= maxSinPhi ) break;
  84
  85       float ss = ey + ey1;
  86       float ex1 = sqrt(1 - ey1*ey1);
  87
  88       float dxBz = dx * bz;
  89
  90       float cc = ex + ex1;
  91       float dxcci = dx * Reciprocal(cc);
  92       float kdx205 = kdx*kdx*0.5f;
  93
  94       float dy = dxcci * ss;
  95       float norm2 = float(1.f) + ey*ey1 + ex*ex1;
  96       float dl = dxcci * sqrt( norm2 + norm2 );
  97
  98       float dS;
  99       {
 100         float dSin = float(0.5f)*k*dl;
 101         float a = dSin*dSin;
 102         const float k2 = 1.f/6.f;
 103         const float k4 = 3.f/40.f;
 104         //const float k6 = 5.f/112.f;
 105         dS = dl + dl*a*(k2 + a*(k4 ));//+ k6*a) );
 106       }
 107
 108       float ex1i = Reciprocal(ex1);
 109       float dz = dS * t0.DzDs();
 110
 111       dL = -dS * t0.DlDs();
 112
 113       float hh = dxcci*ex1i*(2.f+kdx205);
 114       float h2 = hh * t0.SecPhi();
 115       float h4 = bz*dxcci*hh;
 116
 117       float d2 = fP[2] - t0.SinPhi();
 118       float d3 = fP[3] - t0.DzDs();
 119       float d4 = fP[4] - t0.QPt();
 120
 121
 122       fX+=dx;
 123       fP[0]+= dy     + h2 * d2           +   h4 * d4;
 124       fP[1]+= dz               + dS * d3;
 125       fP[2] = ey1 +     d2           + dxBz * d4;
 126
 127       t0.CosPhi() = ex1;
 128       t0.SecPhi() = ex1i;
 129       t0.SinPhi() = ey1;
 130
 131       {
 132         const float *cy = param.GetParamS0Par(0,rowType[ihit]);
 133         const float *cz = param.GetParamS0Par(1,rowType[ihit]);
 134
 135         float secPhi2 = ex1i*ex1i;
 136         float zz = fabs( kZLength - fabs(fP[2]) );
 137         float zz2 = zz*zz;
 138         float angleY2 = secPhi2 - 1.f;
 139         float angleZ2 = trDzDs2 * secPhi2 ;
 140
 141         float cy0 = cy[0] + cy[1]*zz + cy[3]*zz2;
 142         float cy1 = cy[2] + cy[5]*zz;
 143         float cy2 = cy[4];
 144         float cz0 = cz[0] + cz[1]*zz + cz[3]*zz2;
 145         float cz1 = cz[2] + cz[5]*zz;
 146         float cz2 = cz[4];
 147
 148         err2Y = fabs( cy0 + angleY2 * ( cy1 + angleY2*cy2 ) );
 149         err2Z = fabs( cz0 + angleZ2 * ( cz1 + angleZ2*cz2 ) );
 150      }
 151
 152
 153       if ( first ) {
 154         fP[0] = y[ihit];
 155         fP[1] = z[ihit];
 156         SetCov( 0, err2Y );
 157         SetCov( 1,  0 );
 158         SetCov( 2, err2Z);
 159         SetCov( 3,  0 );
 160         SetCov( 4,  0 );
 161         SetCov( 5,  1 );
 162         SetCov( 6,  0 );
 163         SetCov( 7,  0 );
 164         SetCov( 8,  0 );
 165         SetCov( 9,  1 );
 166         SetCov( 10,  0 );
 167         SetCov( 11,  0 );
 168         SetCov( 12,  0 );
 169         SetCov( 13,  0 );
 170         SetCov( 14,  10 );
 171         SetChi2( 0 );
 172         SetNDF( -3 );
 173         CalculateFitParameters( par, kRhoOverRadLen, kRho, UseMeanPt );
 174         first = 0;
 175         N+=1;
 176         continue;
 177       }
 178
 179       float c20 = fC[3];
 180       float c21 = fC[4];
 181       float c22 = fC[5];
 182       float c30 = fC[6];
 183       float c31 = fC[7];
 184       float c32 = fC[8];
 185       float c33 = fC[9];
 186       float c40 = fC[10];
 187       float c41 = fC[11];
 188       float c42 = fC[12];
 189       float c43 = fC[13];
 190       float c44 = fC[14];
 191
 192       float c20ph4c42 =  c20 + h4*c42;
 193       float h2c22 = h2*c22;
 194       float h4c44 = h4*c44;
 195
 196       float n6 = c30 + h2*c32 + h4*c43;
 197       float n7 = c31 + dS*c33;
 198       float n10 = c40 + h2*c42 + h4c44;
 199       float n11 = c41 + dS*c43;
 200       float n12 = c42 + dxBz*c44;
 201
 202       fC[8] = c32 + dxBz * c43;
 203
 204       fC[0]+= h2*h2c22 + h4*h4c44 + float(2.f)*( h2*c20ph4c42  + h4*c40 );
 205
 206       fC[1]+= h2*c21 + h4*c41 + dS*n6;
 207       fC[6] = n6;
 208
 209       fC[2]+= dS*(c31 + n7);
 210       fC[7] = n7;
 211
 212       fC[3] = c20ph4c42 + h2c22  + dxBz*n10;
 213       fC[10] = n10;
 214
 215       fC[4] = c21 + dS*c32 + dxBz*n11;
 216       fC[11] = n11;
 217
 218       fC[5] = c22 + dxBz*( c42 + n12 );
 219       fC[12] = n12;
 220
 221     } // end transport block
 222
 223
 224     float &fC22 = fC[5];
 225     float &fC33 = fC[9];
 226     float &fC40 = fC[10];
 227     float &fC41 = fC[11];
 228     float &fC42 = fC[12];
 229     float &fC43 = fC[13];
 230     float &fC44 = fC[14];
 231
 232     float
 233       c00 = fC[ 0],
 234       c11 = fC[ 2],
 235       c20 = fC[ 3],
 236       c31 = fC[ 7];
 237
 238
 239     // MS block
 240
 241     float dLmask = 0.f;
 242     bool maskMS = ( fabs( dL ) < par.fDLMax );
 243
 244
 245     // Filter block
 246
 247     float mS0 = Reciprocal(err2Y + c00);
 248
 249     // MS block
 250     Assign( dLmask, maskMS, dL );
 251
 252     // Filter block
 253
 254     float  z0 = y[ihit] - fP[0];
 255     float mS2 = Reciprocal(err2Z + c11);
 256
 257     if( fabs( fP[2] + z0*c20*mS0  ) > maxSinPhi ) break;
 258
 259     // MS block
 260
 261     float dLabs = fabs( dLmask);
 262     float corr = float(1.f) - par.fEP2* dLmask ;
 263
 264     fP[4]*= corr;
 265     fC40 *= corr;
 266     fC41 *= corr;
 267     fC42 *= corr;
 268     fC43 *= corr;
 269     fC44  = fC44*corr*corr + dLabs*par.fSigmadE2;
 270
 271     fC22 += dLabs * par.fK22 * (float(1.f)-fP[2]*fP[2]);
 272     fC33 += dLabs * par.fK33;
 273     fC43 += dLabs * par.fK43;
 274
 275
 276     // Filter block
 277
 278     float c40 = fC40;
 279
 280     // K = CHtS
 281
 282     float k00, k11, k20, k31, k40;
 283
 284     k00 = c00 * mS0;
 285     k20 = c20 * mS0;
 286     k40 = c40 * mS0;
 287     fChi2  += mS0*z0*z0;
 288     fP[0] += k00 * z0;
 289     fP[2] += k20 * z0;
 290     fP[4] += k40 * z0;
 291     fC[ 0] -= k00 * c00 ;
 292     fC[ 5] -= k20 * c20 ;
 293     fC[10] -= k00 * c40 ;
 294     fC[12] -= k40 * c20 ;
 295     fC[ 3] -= k20 * c00 ;
 296     fC[14] -= k40 * c40 ;
 297
 298     float  z1 = z[ihit] - fP[1];
 299
 300     k11 = c11 * mS2;
 301     k31 = c31 * mS2;
 302
 303     fChi2  +=  mS2*z1*z1 ;
 304     fNDF  += 2;
 305     N+=1;
 306
 307     fP[1] += k11 * z1;
 308     fP[3] += k31 * z1;
 309
 310     fC[ 7] -= k31 * c11;
 311     fC[ 2] -= k11 * c11;
 312     fC[ 9] -= k31 * c31;
 313   }
 314 }
 315
 316 GPUd() bool AliHLTTPCGMTrackParam::CheckNumericalQuality() const
 317 {
 318   //* Check that the track parameters and covariance matrix are reasonable
 319
 320   bool ok = AliHLTTPCCAMath::Finite(fX) && AliHLTTPCCAMath::Finite( fChi2 ) && AliHLTTPCCAMath::Finite( fNDF );
 321
 322   const float *c = fC;
 323   for ( int i = 0; i < 15; i++ ) ok = ok && AliHLTTPCCAMath::Finite( c[i] );
 324   for ( int i = 0; i < 5; i++ ) ok = ok && AliHLTTPCCAMath::Finite( fP[i] );
 325
 326   if ( c[0] <= 0 || c[2] <= 0 || c[5] <= 0 || c[9] <= 0 || c[14] <= 0 ) ok = 0;
 327   if ( c[0] > 5. || c[2] > 5. || c[5] > 2. || c[9] > 2.
 328        //|| ( CAMath::Abs( QPt() ) > 1.e-2 && c[14] > 2. )
 329        ) ok = 0;
 330
 331   if ( fabs( fP[2] ) > .999 ) ok = 0;
 332   if ( fabs( fP[4] ) > 1. / 0.05 ) ok = 0;
 333   if( ok ){
 334     ok = ok
 335       && ( c[1]*c[1]<=c[2]*c[0] )
 336       && ( c[3]*c[3]<=c[5]*c[0] )
 337       && ( c[4]*c[4]<=c[5]*c[2] )
 338       && ( c[6]*c[6]<=c[9]*c[0] )
 339       && ( c[7]*c[7]<=c[9]*c[2] )
 340       && ( c[8]*c[8]<=c[9]*c[5] )
 341       && ( c[10]*c[10]<=c[14]*c[0] )
 342       && ( c[11]*c[11]<=c[14]*c[2] )
 343       && ( c[12]*c[12]<=c[14]*c[5] )
 344       && ( c[13]*c[13]<=c[14]*c[9] );
 345   }
 346   return ok;
 347 }
 348
 349 //*
 350 //*  Multiple scattering and energy losses
 351 //*
 352
 353 GPUd() float AliHLTTPCGMTrackParam::ApproximateBetheBloch( float beta2 )
 354 {
 355   //------------------------------------------------------------------
 356   // This is an approximation of the Bethe-Bloch formula with
 357   // the density effect taken into account at beta*gamma > 3.5
 358   // (the approximation is reasonable only for solid materials)
 359   //------------------------------------------------------------------
 360
 361   const float log0 = log( float(5940.f));
 362   const float log1 = log( float(3.5f*5940.f) );
 363
 364   bool bad = (beta2 >= .999f)||( beta2 < 1.e-8f );
 365
 366   Assign( beta2, bad, 0.5f);
 367
 368   float a = beta2 / ( 1.f - beta2 );
 369   float b = 0.5*log(a);
 370   float d =  0.153e-3 / beta2;
 371   float c = b - beta2;
 372
 373   float ret = d*(log0 + b + c );
 374   float case1 = d*(log1 + c );
 375
 376   Assign( ret, ( a > 3.5*3.5  ), case1);
 377   Assign( ret,  bad, 0. );
 378
 379   return ret;
 380 }
 381
 382  GPUd() void AliHLTTPCGMTrackParam::CalculateFitParameters( AliHLTTPCGMTrackFitParam &par, float RhoOverRadLen,  float Rho, bool NoField, float mass )
 383 {
 384   //*!
 385
 386   float qpt = fP[4];
 387   if( NoField ) qpt = 1./0.35;
 388
 389   float p2 = ( 1. + fP[3] * fP[3] );
 390   float k2 = qpt * qpt;
 391   Assign( k2, (  k2 < 1.e-4f ), 1.e-4f );
 392
 393   float mass2 = mass * mass;
 394   float beta2 = p2 / ( p2 + mass2 * k2 );
 395
 396   float pp2 = p2 / k2; // impuls 2
 397
 398   //par.fBethe = BetheBlochGas( pp2/mass2);
 399   par.fBetheRho = ApproximateBetheBloch( pp2 / mass2 )*Rho;
 400   par.fE = sqrt( pp2 + mass2 );
 401   par.fTheta2 = ( 14.1*14.1/1.e6 ) / ( beta2 * pp2 )*RhoOverRadLen;
 402   par.fEP2 = par.fE / pp2;
 403
 404   // Approximate energy loss fluctuation (M.Ivanov)
 405
 406   const float knst = 0.07; // To be tuned.
 407   par.fSigmadE2 = knst * par.fEP2 * qpt;
 408   par.fSigmadE2 = par.fSigmadE2 * par.fSigmadE2;
 409
 410   float k22 = 1. + fP[3] * fP[3];
 411   par.fK22 = par.fTheta2*k22;
 412   par.fK33 = par.fK22 * k22;
 413   par.fK43 = 0.;
 414   par.fK44 =  par.fTheta2*fP[3] * fP[3] * k2;
 415
 416   float br=1.e-8f;
 417   Assign( br, ( par.fBetheRho>1.e-8f ), par.fBetheRho );
 418   par.fDLMax = 0.3*par.fE * Reciprocal( br );
 419
 420   par.fEP2*=par.fBetheRho;
 421   par.fSigmadE2 = par.fSigmadE2*par.fBetheRho+par.fK44;
 422 }
 423
 424
 425
 426
 427 //*
 428 //* Rotation
 429 //*
 430
 431
 432 GPUd() bool AliHLTTPCGMTrackParam::Rotate( float alpha, AliHLTTPCGMTrackLinearisation &t0, float maxSinPhi )
 433 {
 434   //* Rotate the coordinate system in XY on the angle alpha
 435
 436   float cA = CAMath::Cos( alpha );
 437   float sA = CAMath::Sin( alpha );
 438   float x0 = X(), y0 = Y(), sP = t0.SinPhi(), cP = t0.CosPhi();
 439   float cosPhi = cP * cA + sP * sA;
 440   float sinPhi = -cP * sA + sP * cA;
 441
 442   if ( CAMath::Abs( sinPhi ) > maxSinPhi || CAMath::Abs( cosPhi ) < 1.e-2 || CAMath::Abs( cP ) < 1.e-2  ) return 0;
 443
 444   //float J[5][5] = { { j0, 0, 0,  0,  0 }, // Y
 445   //                    {  0, 1, 0,  0,  0 }, // Z
 446   //                    {  0, 0, j2, 0,  0 }, // SinPhi
 447   //                  {  0, 0, 0,  1,  0 }, // DzDs
 448   //                  {  0, 0, 0,  0,  1 } }; // Kappa
 449
 450   float j0 = cP / cosPhi;
 451   float j2 = cosPhi / cP;
 452   float d[2] = {Y() - y0, SinPhi() - sP};
 453
 454   X() = ( x0*cA +  y0*sA );
 455   Y() = ( -x0*sA +  y0*cA + j0*d[0] );
 456   t0.CosPhi() = fabs( cosPhi );
 457   t0.SecPhi() = ( 1./t0.CosPhi() );
 458   t0.SinPhi() = ( sinPhi );
 459
 460   SinPhi() = ( sinPhi + j2*d[1] );
 461
 462   fC[0] *= j0 * j0;
 463   fC[1] *= j0;
 464   fC[3] *= j0;
 465   fC[6] *= j0;
 466   fC[10] *= j0;
 467
 468   fC[3] *= j2;
 469   fC[4] *= j2;
 470   fC[5] *= j2 * j2;
 471   fC[8] *= j2;
 472   fC[12] *= j2;
 473   if( cosPhi <0 ){ // change direction
 474     t0.SinPhi() = -sinPhi;
 475     t0.DzDs() = -t0.DzDs();
 476     t0.DlDs() = -t0.DlDs();
 477     t0.QPt() = -t0.QPt();
 478     SinPhi() = -SinPhi();
 479     DzDs() = -DzDs();
 480     QPt() = -QPt();
 481     fC[3] = -fC[3];
 482     fC[4] = -fC[4];
 483     fC[6] = -fC[6];
 484     fC[7] = -fC[7];
 485     fC[10] = -fC[10];
 486     fC[11] = -fC[11];
 487   }
 488
 489   return 1;
 490 }
 491
 492 #if !defined(HLTCA_STANDALONE) & !defined(HLTCA_GPUCODE)
 493 bool AliHLTTPCGMTrackParam::GetExtParam( AliExternalTrackParam &T, double alpha ) const
 494 {
 495   //* Convert from AliHLTTPCGMTrackParam to AliExternalTrackParam parameterisation,
 496   //* the angle alpha is the global angle of the local X axis
 497
 498   bool ok = CheckNumericalQuality();
 499
 500   double par[5], cov[15];
 501   for ( int i = 0; i < 5; i++ ) par[i] = fP[i];
 502   for ( int i = 0; i < 15; i++ ) cov[i] = fC[i];
 503
 504   if ( par[2] > .99 ) par[2] = .99;
 505   if ( par[2] < -.99 ) par[2] = -.99;
 506
 507   if ( fabs( par[4] ) < 1.e-5 ) par[4] = 1.e-5; // some other software will crash if q/Pt==0
 508   if ( fabs( par[4] ) > 1./0.08 ) ok = 0; // some other software will crash if q/Pt is too big
 509
 510   T.Set( (double) fX, alpha, par, cov );
 511   return ok;
 512 }
 513
 514
 515
 516 void AliHLTTPCGMTrackParam::SetExtParam( const AliExternalTrackParam &T )
 517 {
 518   //* Convert from AliExternalTrackParam parameterisation
 519
 520   for ( int i = 0; i < 5; i++ ) fP[i] = T.GetParameter()[i];
 521   for ( int i = 0; i < 15; i++ ) fC[i] = T.GetCovariance()[i];
 522   fX = T.GetX();
 523   if ( fP[2] > .999 ) fP[2] = .999;
 524   if ( fP[2] < -.999 ) fP[2] = -.999;
 525 }
 526 #endif
 527
 528 #ifdef HLTCA_GPUCODE
 529
 530 #include "AliHLTTPCGMMergedTrack.h"
 531
 532 GPUg() void RefitTracks(AliHLTTPCGMMergedTrack* tracks, int nTracks, float* PolinomialFieldBz, float* x, float* y, float* z, unsigned int* rowType, float* alpha, AliHLTTPCCAParam* param)
 533 {
 534         for (int i = blockDim.x * blockIdx.x + threadIdx.x;i < nTracks;i += blockDim.x * gridDim.x)
 535         {
 536                 //This is in fact a copy of ReFit() in AliHLTTPCGMMerger.cxx
 537                 AliHLTTPCGMMergedTrack& track = tracks[i];
 538                 float Alpha = track.Alpha();
 539                 int N = track.NClusters();
 540                 AliHLTTPCGMTrackParam t = track.Param();
 541
 542                 t.Fit(PolinomialFieldBz,
 543                         x+track.FirstClusterRef(),
 544                         y+track.FirstClusterRef(),
 545                         z+track.FirstClusterRef(),
 546                         rowType+track.FirstClusterRef(),
 547                         alpha+track.FirstClusterRef(),
 548                         *param,
 549                         N,
 550                         Alpha,
 551                         0
 552                         );
 553
 554                 if ( fabs( t.QPt() ) < 1.e-4 ) t.QPt() = 1.e-4 ;
 555
 556                 bool ok = N >= 30 && t.CheckNumericalQuality() && fabs( t.SinPhi() ) <= .999;
 557                 track.SetOK(ok);
 558                 if( !ok ) continue;
 559
 560                 if( 1 ){//SG!!!
 561                   track.SetNClusters( N );
 562                   track.Param() = t;
 563                   track.Alpha() = Alpha;
 564                 }
 565
 566                 {
 567                   int ind = track.FirstClusterRef();
 568                   float alphaalpha = alpha[ind];
 569                   float xx = x[ind];
 570                   float yy = y[ind];
 571                   float zz = z[ind];
 572                   float sinA = AliHLTTPCCAMath::Sin( alphaalpha - track.Alpha());
 573                   float cosA = AliHLTTPCCAMath::Cos( alphaalpha - track.Alpha());
 574                   track.SetLastX( xx*cosA - yy*sinA );
 575                   track.SetLastY( xx*sinA + yy*cosA );
 576                   track.SetLastZ( zz );
 577                 }
 578         }
 579 }
 580
 581 #endif