HLT/TPCLib/tracking-ca/AliHLTTPCCATrackletConstructor.cxx

   1 // @(#) $Id: AliHLTTPCCATrackletConstructor.cxx 27042 2008-07-02 12:06:02Z richterm $
   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 "AliHLTTPCCATracker.h"
  22 #include "AliHLTTPCCATrackParam.h"
  23 #include "AliHLTTPCCATrackParam.h"
  24 #include "AliHLTTPCCAGrid.h"
  25 #include "AliHLTTPCCAHitArea.h"
  26 #include "AliHLTTPCCAMath.h"
  27 #include "AliHLTTPCCADef.h"
  28 #include "AliHLTTPCCATracklet.h"
  29 #include "AliHLTTPCCATrackletConstructor.h"
  30 //#include "AliHLTTPCCAPerformance.h"
  31 //#include "TH1D.h"
  32
  33 //#define DRAW
  34
  35 #ifdef DRAW
  36 #include "AliHLTTPCCADisplay.h"
  37 #endif
  38
  39
  40 GPUd() void AliHLTTPCCATrackletConstructor::Step0
  41 ( int nBlocks, int /*nThreads*/, int iBlock, int iThread,
  42   AliHLTTPCCASharedMemory &s, AliHLTTPCCAThreadMemory &r, AliHLTTPCCATracker &tracker, AliHLTTPCCATrackParam &/*tParam*/ )
  43 {
  44   // reconstruction of tracklets, step 0
  45
  46   r.fIsMemThread = ( iThread < NMemThreads() );
  47   if ( iThread == 0 ) {
  48     int nTracks = *tracker.NTracklets();
  49     int nTrPerBlock = nTracks / nBlocks + 1;
  50     s.fNRows = tracker.Param().NRows();
  51     s.fItr0 = nTrPerBlock * iBlock;
  52     s.fItr1 = s.fItr0 + nTrPerBlock;
  53     if ( s.fItr1 > nTracks ) s.fItr1 = nTracks;
  54     s.fUsedHits = tracker.HitWeights();
  55     s.fMinStartRow = 158;
  56     s.fMaxStartRow = 0;
  57   }
  58   if ( iThread < 32 ) {
  59     s.fMinStartRow32[iThread] = 158;
  60     s.fMaxStartRow32[iThread] = 0;
  61   }
  62 }
  63
  64
  65 GPUd() void AliHLTTPCCATrackletConstructor::Step1
  66 ( int /*nBlocks*/, int /*nThreads*/, int /*iBlock*/, int iThread,
  67   AliHLTTPCCASharedMemory &s, AliHLTTPCCAThreadMemory &r, AliHLTTPCCATracker &tracker, AliHLTTPCCATrackParam &tParam )
  68 {
  69   // reconstruction of tracklets, step 1
  70
  71   r.fItr = s.fItr0 + ( iThread - NMemThreads() );
  72   r.fGo = ( !r.fIsMemThread ) && ( r.fItr < s.fItr1 );
  73   r.fSave = r.fGo;
  74   r.fNHits = 0;
  75
  76   if ( !r.fGo ) return;
  77
  78   r.fStage = 0;
  79
  80   AliHLTTPCCATracklet &tracklet = tracker.Tracklets()[r.fItr];
  81
  82   unsigned int kThread = iThread % 32;//& 00000020;
  83   if ( SAVE() ) for ( int i = 0; i < 160; i++ ) tracklet.SetRowHit( i, -1 );
  84
  85   int id = tracker.TrackletStartHits()[r.fItr];
  86   r.fStartRow = AliHLTTPCCATracker::ID2IRow( id );
  87   r.fEndRow = r.fStartRow;
  88   r.fFirstRow = r.fStartRow;
  89   r.fLastRow = r.fFirstRow;
  90   r.fCurrIH =  AliHLTTPCCATracker::ID2IHit( id );
  91
  92   CAMath::AtomicMin( &s.fMinStartRow32[kThread], r.fStartRow );
  93   CAMath::AtomicMax( &s.fMaxStartRow32[kThread], r.fStartRow );
  94   tParam.SetSinPhi( 0 );
  95   tParam.SetDzDs( 0 );
  96   tParam.SetQPt( 0 );
  97   tParam.SetSignCosPhi( 1 );
  98   tParam.SetChi2( 0 );
  99   tParam.SetNDF( -3 );
 100   tParam.SetCov( 0, 1 );
 101   tParam.SetCov( 1, 0 );
 102   tParam.SetCov( 2, 1 );
 103   tParam.SetCov( 3, 0 );
 104   tParam.SetCov( 4, 0 );
 105   tParam.SetCov( 5, 1 );
 106   tParam.SetCov( 6, 0 );
 107   tParam.SetCov( 7, 0 );
 108   tParam.SetCov( 8, 0 );
 109   tParam.SetCov( 9, 1 );
 110   tParam.SetCov( 10, 0 );
 111   tParam.SetCov( 11, 0 );
 112   tParam.SetCov( 12, 0 );
 113   tParam.SetCov( 13, 0 );
 114   tParam.SetCov( 14, 10. );
 115
 116 }
 117
 118 GPUd() void AliHLTTPCCATrackletConstructor::Step2
 119 ( int /*nBlocks*/, int nThreads, int /*iBlock*/, int iThread,
 120   AliHLTTPCCASharedMemory &s, AliHLTTPCCAThreadMemory &/*r*/, AliHLTTPCCATracker &/*tracker*/, AliHLTTPCCATrackParam &/*tParam*/ )
 121 {
 122   // reconstruction of tracklets, step 2
 123
 124   if ( iThread == 0 ) {
 125     //CAMath::AtomicMinGPU(&s.fMinRow, s.fMinRow32[iThread]);
 126     int minStartRow = 158;
 127     int maxStartRow = 0;
 128     int n = ( nThreads > 32 ) ? 32 : nThreads;
 129     for ( int i = 0; i < n; i++ ) {
 130       if ( s.fMinStartRow32[i] < minStartRow ) minStartRow = s.fMinStartRow32[i];
 131       if ( s.fMaxStartRow32[i] > maxStartRow ) maxStartRow = s.fMaxStartRow32[i];
 132     }
 133     s.fMinStartRow = minStartRow;
 134     s.fMaxStartRow = maxStartRow;
 135   }
 136 }
 137
 138 GPUd() void AliHLTTPCCATrackletConstructor::ReadData
 139 ( int iThread, AliHLTTPCCASharedMemory &s, AliHLTTPCCAThreadMemory &r, AliHLTTPCCATracker &tracker, int iRow )
 140 {
 141   // reconstruction of tracklets, read data step
 142
 143   if ( r.fIsMemThread ) {
 144     const AliHLTTPCCARow &row = tracker.Row( iRow );
 145     bool jr = !r.fCurrentData;
 146     int n = row.FullSize();
 147     const uint4* gMem = tracker.RowData() + row.FullOffset();
 148     uint4 *sMem = s.fData[jr];
 149     for ( int i = iThread; i < n; i += NMemThreads() ) sMem[i] = gMem[i];
 150   }
 151 }
 152
 153
 154 GPUd() void AliHLTTPCCATrackletConstructor::StoreTracklet
 155 ( int /*nBlocks*/, int /*nThreads*/, int /*iBlock*/, int /*iThread*/,
 156   AliHLTTPCCASharedMemory &/*s*/, AliHLTTPCCAThreadMemory &r, AliHLTTPCCATracker &tracker, AliHLTTPCCATrackParam &tParam )
 157 {
 158   // reconstruction of tracklets, tracklet store step
 159
 160   if ( !r.fSave ) return;
 161
 162   //AliHLTTPCCAPerformance::Instance().HNHitsPerTrackCand()->Fill(r.fNHits);
 163
 164   do {
 165     {
 166       //std::cout<<"tracklet to store: "<<r.fItr<<", nhits = "<<r.fNHits<<std::endl;
 167     }
 168
 169     if ( r.fNHits < 5 ) {
 170       r.fNHits = 0;
 171       break;
 172     }
 173
 174     if ( 0 ) {
 175       if ( 1. / .5 < CAMath::Abs( tParam.QPt() ) ) { //SG!!!
 176         r.fNHits = 0;
 177         break;
 178       }
 179     }
 180
 181     {
 182       bool ok = 1;
 183       const float *c = tParam.Cov();
 184       for ( int i = 0; i < 15; i++ ) ok = ok && CAMath::Finite( c[i] );
 185       for ( int i = 0; i < 5; i++ ) ok = ok && CAMath::Finite( tParam.Par()[i] );
 186       ok = ok && ( tParam.X() > 50 );
 187
 188       if ( c[0] <= 0 || c[2] <= 0 || c[5] <= 0 || c[9] <= 0 || c[14] <= 0 ) ok = 0;
 189
 190       if ( !ok ) {
 191         r.fNHits = 0;
 192         break;
 193       }
 194     }
 195   } while ( 0 );
 196
 197   if ( !SAVE() ) return;
 198
 199   AliHLTTPCCATracklet &tracklet = tracker.Tracklets()[r.fItr];
 200
 201   tracklet.SetNHits( r.fNHits );
 202
 203   if ( r.fNHits > 0 ) {
 204 #ifdef DRAW
 205     if ( 0 ) {
 206       std::cout << "store tracklet " << r.fItr << ", nhits = " << r.fNHits << std::endl;
 207       if ( AliHLTTPCCADisplay::Instance().DrawTracklet( tParam, hitstore, kBlue, 1. ) ) {
 208         AliHLTTPCCADisplay::Instance().Ask();
 209       }
 210     }
 211 #endif
 212     if ( CAMath::Abs( tParam.Par()[4] ) < 1.e-8 ) tParam.SetPar( 4, 1.e-8 );
 213     tracklet.SetFirstRow( r.fFirstRow );
 214     tracklet.SetLastRow( r.fLastRow );
 215     tracklet.SetParam( tParam );
 216     int w = ( r.fNHits << 16 ) + r.fItr;
 217     for ( int iRow = 0; iRow < 160; iRow++ ) {
 218       int ih = tracklet.RowHit( iRow );
 219       if ( ih >= 0 ) {
 220         int ihTot = tracker.Row( iRow ).FirstHit() + ih;
 221         CAMath::AtomicMax( tracker.HitWeights() + ihTot, w );
 222       }
 223     }
 224   }
 225 }
 226
 227 GPUd() void AliHLTTPCCATrackletConstructor::UpdateTracklet
 228 ( int /*nBlocks*/, int /*nThreads*/, int /*iBlock*/, int /*iThread*/,
 229   AliHLTTPCCASharedMemory &s, AliHLTTPCCAThreadMemory &r, AliHLTTPCCATracker &tracker, AliHLTTPCCATrackParam &tParam, int iRow )
 230 {
 231   // reconstruction of tracklets, tracklets update step
 232
 233   //std::cout<<"Update tracklet: "<<r.fItr<<" "<<r.fGo<<" "<<r.fStage<<" "<<iRow<<std::endl;
 234   bool drawSearch = 0;//r.fItr==2;
 235   bool drawFit = 0;//r.fItr==2;
 236   bool drawFitted = drawFit ;//|| 1;//r.fItr==16;
 237
 238   if ( !r.fGo ) return;
 239
 240   const int kMaxRowGap = 4;
 241
 242   AliHLTTPCCATracklet &tracklet = tracker.Tracklets()[r.fItr];
 243
 244   const AliHLTTPCCARow &row = tracker.Row( iRow );
 245
 246   float y0 = row.Grid().YMin();
 247   float stepY = row.HstepY();
 248   float z0 = row.Grid().ZMin();
 249   float stepZ = row.HstepZ();
 250   float stepYi = row.HstepYi();
 251   float stepZi = row.HstepZi();
 252
 253   if ( r.fStage == 0 ) { // fitting part
 254     do {
 255
 256       if ( iRow < r.fStartRow || r.fCurrIH < 0  ) break;
 257
 258       if ( ( iRow - r.fStartRow ) % 2 != 0 ) break; // SG!!! - jump over the row
 259
 260       uint4 *tmpint4 = s.fData[r.fCurrentData];
 261       ushort2 hh = reinterpret_cast<ushort2*>( tmpint4 )[r.fCurrIH];
 262
 263       int oldIH = r.fCurrIH;
 264       r.fCurrIH = reinterpret_cast<short*>( tmpint4 )[row.FullLinkOffset() + r.fCurrIH];
 265
 266       float x = row.X();
 267       float y = y0 + hh.x * stepY;
 268       float z = z0 + hh.y * stepZ;
 269       if ( drawFit ) std::cout << " fit tracklet: new hit " << oldIH << ", xyz=" << x << " " << y << " " << z << std::endl;
 270
 271       if ( iRow == r.fStartRow ) {
 272         tParam.SetX( x );
 273         tParam.SetY( y );
 274         tParam.SetZ( z );
 275         r.fLastY = y;
 276         r.fLastZ = z;
 277         //#ifdef DRAW
 278         if ( drawFit ) std::cout << " fit tracklet " << r.fItr << ", row " << iRow << " first row" << std::endl;
 279         //#endif
 280       } else {
 281
 282         float err2Y, err2Z;
 283         float dx = x - tParam.X();
 284         float dy = y - r.fLastY;//tParam.Y();
 285         float dz = z - r.fLastZ;//tParam.Z();
 286         r.fLastY = y;
 287         r.fLastZ = z;
 288
 289         float ri = 1. / CAMath::Sqrt( dx * dx + dy * dy );
 290         if ( iRow == r.fStartRow + 1 ) {
 291           tParam.SetSinPhi( dy*ri );
 292           tParam.SetSignCosPhi( dx );
 293           tParam.SetDzDs( dz*ri );
 294           std::cout << "Init. errors... " << r.fItr << std::endl;
 295           tracker.GetErrors2( iRow, tParam, err2Y, err2Z );
 296           std::cout << "Init. errors = " << err2Y << " " << err2Z << std::endl;
 297           tParam.SetCov( 0, err2Y );
 298           tParam.SetCov( 2, err2Z );
 299         }
 300         if ( drawFit ) {
 301           //#ifdef DRAW
 302           std::cout << " fit tracklet " << r.fItr << ", row " << iRow << " transporting.." << std::endl;
 303           std::cout << " params before transport=" << std::endl;
 304           tParam.Print();
 305           //#endif
 306         }
 307         float sinPhi, cosPhi;
 308         if ( r.fNHits >= 10 && CAMath::Abs( tParam.SinPhi() ) < .99 ) {
 309           sinPhi = tParam.SinPhi();
 310           cosPhi = CAMath::Sqrt( 1 - sinPhi * sinPhi );
 311         } else {
 312           sinPhi = dy * ri;
 313           cosPhi = dx * ri;
 314         }
 315         //#ifdef DRAW
 316         if ( drawFit ) std::cout << "sinPhi0 = " << sinPhi << ", cosPhi0 = " << cosPhi << std::endl;
 317         //#endif
 318         if ( !tParam.TransportToX( x, sinPhi, cosPhi, tracker.Param().Bz(), -1 ) ) {
 319           //#ifdef DRAW
 320           if ( drawFit ) std::cout << " tracklet " << r.fItr << ", row " << iRow << ": can not transport!!" << std::endl;
 321 //#endif
 322           if ( SAVE() ) tracklet.SetRowHit( iRow, -1 );
 323           break;
 324         }
 325         //std::cout<<"mark1 "<<r.fItr<<std::endl;
 326         //tParam.Print();
 327         tracker.GetErrors2( iRow, tParam.GetZ(), sinPhi, cosPhi, tParam.GetDzDs(), err2Y, err2Z );
 328         //std::cout<<"mark2"<<std::endl;
 329
 330         if ( drawFit ) {
 331           //#ifdef DRAW
 332           std::cout << " params after transport=" << std::endl;
 333           tParam.Print();
 334           std::cout << "fit tracklet before filter: " << r.fItr << ", row " << iRow << " errs=" << err2Y << " " << err2Z << std::endl;
 335           //#endif
 336 #ifdef DRAW
 337           AliHLTTPCCADisplay::Instance().DrawTracklet( tParam, hitstore, kBlue, 2., 1 );
 338           AliHLTTPCCADisplay::Instance().Ask();
 339 #endif
 340         }
 341         if ( !tParam.Filter( y, z, err2Y, err2Z, .99 ) ) {
 342           //#ifdef DRAW
 343           if ( drawFit ) std::cout << " tracklet " << r.fItr << ", row " << iRow << ": can not filter!!" << std::endl;
 344           //#endif
 345           if ( SAVE() ) tracklet.SetRowHit( iRow, -1 );
 346           break;
 347         }
 348       }
 349       if ( SAVE() ) tracklet.SetRowHit( iRow, oldIH );
 350       if ( drawFit ) {
 351         //#ifdef DRAW
 352         std::cout << "fit tracklet after filter " << r.fItr << ", row " << iRow << std::endl;
 353         tParam.Print();
 354         //#endif
 355 #ifdef DRAW
 356         AliHLTTPCCADisplay::Instance().DrawTracklet( tParam, hitstore, kGreen, 2. );
 357         AliHLTTPCCADisplay::Instance().Ask();
 358 #endif
 359       }
 360       r.fNHits++;
 361       r.fLastRow = iRow;
 362       r.fEndRow = iRow;
 363       break;
 364     } while ( 0 );
 365
 366     if ( r.fCurrIH < 0 ) {
 367       //#ifdef DRAW
 368       if ( drawFitted ) std::cout << "fitted tracklet " << r.fItr << ", nhits=" << r.fNHits << std::endl;
 369       //#endif
 370       r.fStage = 1;
 371       //AliHLTTPCCAPerformance::Instance().HNHitsPerSeed()->Fill(r.fNHits);
 372       if ( r.fNHits < 3 ) { r.fNHits = 0; r.fGo = 0;}//SG!!!
 373       if ( CAMath::Abs( tParam.SinPhi() ) > .999 ) {
 374         //#ifdef DRAW
 375         if ( drawFitted ) std::cout << " fitted tracklet  error: sinPhi=" << tParam.SinPhi() << std::endl;
 376         //#endif
 377         r.fNHits = 0; r.fGo = 0;
 378       } else {
 379         //tParam.SetCosPhi( CAMath::Sqrt(1-tParam.SinPhi()*tParam.SinPhi()) );
 380       }
 381       if ( drawFitted ) {
 382         //#ifdef DRAW
 383         std::cout << "fitted tracklet " << r.fItr << " miss=" << r.fNMissed << " go=" << r.fGo << std::endl;
 384         tParam.Print();
 385 #ifdef DRAW
 386         AliHLTTPCCADisplay::Instance().DrawTracklet( tParam, hitstore, kBlue );
 387         AliHLTTPCCADisplay::Instance().Ask();
 388 #endif
 389       }
 390     }
 391   } else { // forward/backward searching part
 392     do {
 393       if ( drawSearch ) {
 394         //#ifdef DRAW
 395         std::cout << "search tracklet " << r.fItr << " row " << iRow << " miss=" << r.fNMissed << " go=" << r.fGo << " stage=" << r.fStage << std::endl;
 396         //#endif
 397       }
 398
 399       if ( r.fStage == 2 && ( ( iRow >= r.fEndRow ) ||
 400                               ( iRow >= r.fStartRow && ( iRow - r.fStartRow ) % 2 == 0 )
 401                             ) ) break;
 402       if ( r.fNMissed > kMaxRowGap  ) {
 403         break;
 404       }
 405
 406       r.fNMissed++;
 407
 408       float x = row.X();
 409       float err2Y, err2Z;
 410       if ( drawSearch ) {
 411         //#ifdef DRAW
 412         std::cout << "tracklet " << r.fItr << " before transport to row " << iRow << " : " << std::endl;
 413         tParam.Print();
 414         //#endif
 415       }
 416       if ( !tParam.TransportToX( x, tParam.SinPhi(), tParam.GetCosPhi(), tracker.Param().Bz(), .99 ) ) {
 417         //#ifdef DRAW
 418         if ( drawSearch ) std::cout << " tracklet " << r.fItr << ", row " << iRow << ": can not transport!!" << std::endl;
 419         //#endif
 420         break;
 421       }
 422       if ( drawSearch ) {
 423         //#ifdef DRAW
 424         std::cout << "tracklet " << r.fItr << " after transport to row " << iRow << " : " << std::endl;
 425         tParam.Print();
 426 #ifdef DRAW
 427         AliHLTTPCCADisplay::Instance().DrawTracklet( tParam, hitstore, kBlue, 2., 1 );
 428         AliHLTTPCCADisplay::Instance().Ask();
 429 #endif
 430       }
 431       uint4 *tmpint4 = s.fData[r.fCurrentData];
 432
 433       ushort2 *hits = reinterpret_cast<ushort2*>( tmpint4 );
 434
 435       float fY = tParam.GetY();
 436       float fZ = tParam.GetZ();
 437       int best = -1;
 438
 439       { // search for the closest hit
 440
 441         int ds;
 442         int fY0 = ( int ) ( ( fY - y0 ) * stepYi );
 443         int fZ0 = ( int ) ( ( fZ - z0 ) * stepZi );
 444         int ds0 = ( ( ( int )1 ) << 30 );
 445         ds = ds0;
 446
 447         unsigned int fIndYmin;
 448         unsigned int fHitYfst = 1, fHitYlst = 0, fHitYfst1 = 1, fHitYlst1 = 0;
 449
 450         fIndYmin = row.Grid().GetBin( ( float )( fY - 1. ), ( float )( fZ - 1. ) );
 451         if ( drawSearch ) {
 452 #ifdef DRAW
 453           std::cout << " tracklet " << r.fItr << ", row " << iRow << ": grid N=" << row.Grid().N() << std::endl;
 454           std::cout << " tracklet " << r.fItr << ", row " << iRow << ": minbin=" << fIndYmin << std::endl;
 455 #endif
 456         }
 457         {
 458           int nY = row.Grid().Ny();
 459
 460           unsigned short *sGridP = ( reinterpret_cast<unsigned short*>( tmpint4 ) ) + row.FullGridOffset();
 461           fHitYfst = sGridP[fIndYmin];
 462           fHitYlst = sGridP[fIndYmin+2];
 463           fHitYfst1 = sGridP[fIndYmin+nY];
 464           fHitYlst1 = sGridP[fIndYmin+nY+2];
 465           if ( drawSearch ) {
 466 #ifdef DRAW
 467             std::cout << " Grid, row " << iRow << ": nHits=" << row.NHits() << ", grid n=" << row.Grid().N() << ", c[n]=" << sGridP[row.Grid().N()] << std::endl;
 468             std::cout << "hit steps = " << stepY << " " << stepZ << std::endl;
 469             std::cout << " Grid bins:" << std::endl;
 470             for ( unsigned int i = 0; i < row.Grid().N(); i++ ) {
 471               std::cout << " bin " << i << ": ";
 472               for ( int j = sGridP[i]; j < sGridP[i+1]; j++ ) {
 473                 ushort2 hh = hits[j];
 474                 float y = y0 + hh.x * stepY;
 475                 float z = z0 + hh.y * stepZ;
 476                 std::cout << "[" << j << "|" << y << "," << z << "] ";
 477               }
 478               std::cout << std::endl;
 479             }
 480 #endif
 481           }
 482           if ( sGridP[row.Grid().N()] != row.NHits() ) {
 483 #ifdef DRAW
 484             std::cout << " grid, row " << iRow << ": nHits=" << row.NHits() << ", grid n=" << row.Grid().N() << ", c[n]=" << sGridP[row.Grid().N()] << std::endl;
 485             //exit(0);
 486 #endif
 487           }
 488         }
 489         if ( drawSearch ) {
 490           //#ifdef DRAW
 491           std::cout << " tracklet " << r.fItr << ", row " << iRow << ", yz= " << fY << "," << fZ << ": search hits=" << fHitYfst << " " << fHitYlst << " / " << fHitYfst1 << " " << fHitYlst1 << std::endl;
 492           std::cout << " hit search :" << std::endl;
 493           //#endif
 494         }
 495         for ( unsigned int fIh = fHitYfst; fIh < fHitYlst; fIh++ ) {
 496           ushort2 hh = hits[fIh];
 497           int ddy = ( int )( hh.x ) - fY0;
 498           int ddz = ( int )( hh.y ) - fZ0;
 499           int dds = CAMath::Abs( ddy ) + CAMath::Abs( ddz );
 500           if ( drawSearch ) {
 501             //#ifdef DRAW
 502             std::cout << fIh << ": hityz= " << hh.x << " " << hh.y << "(" << hh.x*stepY << " " << hh.y*stepZ << "), trackyz=" << fY0 << " " << fZ0 << "(" << fY0*stepY << " " << fZ0*stepZ << "), dy,dz,ds= " << ddy << " " << ddz << " " << dds << "(" << ddy*stepY << " " << ddz*stepZ << std::endl;
 503             //#endif
 504           }
 505           if ( dds < ds ) {
 506             ds = dds;
 507             best = fIh;
 508           }
 509         }
 510
 511         for ( unsigned int fIh = fHitYfst1; fIh < fHitYlst1; fIh++ ) {
 512           ushort2 hh = hits[fIh];
 513           int ddy = ( int )( hh.x ) - fY0;
 514           int ddz = ( int )( hh.y ) - fZ0;
 515           int dds = CAMath::Abs( ddy ) + CAMath::Abs( ddz );
 516           if ( drawSearch ) {
 517             //#ifdef DRAW
 518             std::cout << fIh << ": hityz= " << hh.x << " " << hh.y << "(" << hh.x*stepY << " " << hh.y*stepZ << "), trackyz=" << fY0 << " " << fZ0 << "(" << fY0*stepY << " " << fZ0*stepZ << "), dy,dz,ds= " << ddy << " " << ddz << " " << dds << "(" << ddy*stepY << " " << ddz*stepZ << std::endl;
 519             //#endif
 520           }
 521           if ( dds < ds ) {
 522             ds = dds;
 523             best = fIh;
 524           }
 525         }
 526       }// end of search for the closest hit
 527
 528       if ( best < 0 ) break;
 529       if ( drawSearch ) {
 530         //#ifdef DRAW
 531         std::cout << "hit search " << r.fItr << ", row " << iRow << " hit " << best << " found" << std::endl;
 532 #ifdef DRAW
 533         AliHLTTPCCADisplay::Instance().DrawSliceHit( iRow, best, kRed, 1. );
 534         AliHLTTPCCADisplay::Instance().Ask();
 535         AliHLTTPCCADisplay::Instance().DrawSliceHit( iRow, best, kWhite, 1 );
 536         AliHLTTPCCADisplay::Instance().DrawSliceHit( iRow, best );
 537 #endif
 538       }
 539
 540       ushort2 hh = hits[best];
 541
 542       //std::cout<<"mark 3, "<<r.fItr<<std::endl;
 543       //tParam.Print();
 544       tracker.GetErrors2( iRow, *( ( AliHLTTPCCATrackParam* )&tParam ), err2Y, err2Z );
 545       //std::cout<<"mark 4"<<std::endl;
 546
 547       float y = y0 + hh.x * stepY;
 548       float z = z0 + hh.y * stepZ;
 549       float dy = y - fY;
 550       float dz = z - fZ;
 551
 552       const float kFactor = tracker.Param().HitPickUpFactor() * tracker.Param().HitPickUpFactor() * 3.5 * 3.5;
 553       float sy2 = kFactor * ( tParam.GetErr2Y() +  err2Y );
 554       float sz2 = kFactor * ( tParam.GetErr2Z() +  err2Z );
 555       if ( sy2 > 2. ) sy2 = 2.;
 556       if ( sz2 > 2. ) sz2 = 2.;
 557
 558       if ( drawSearch ) {
 559         //#ifdef DRAW
 560         std::cout << "dy,sy= " << dy << " " << CAMath::Sqrt( sy2 ) << ", dz,sz= " << dz << " " << CAMath::Sqrt( sz2 ) << std::endl;
 561         std::cout << "dy,dz= " << dy << " " << dz << ", sy,sz= " << CAMath::Sqrt( sy2 ) << " " << CAMath::Sqrt( sz2 ) << ", sy,sz= " << CAMath::Sqrt( kFactor*( tParam.GetErr2Y() +  err2Y ) ) << " " << CAMath::Sqrt( kFactor*( tParam.GetErr2Z() +  err2Z ) ) << std::endl;
 562         //#endif
 563       }
 564       if ( CAMath::FMulRZ( dy, dy ) > sy2 || CAMath::FMulRZ( dz, dz ) > sz2  ) {
 565         if ( drawSearch ) {
 566           //#ifdef DRAW
 567
 568           std::cout << "found hit is out of the chi2 window\n " << std::endl;
 569           //#endif
 570         }
 571         break;
 572       }
 573 #ifdef DRAW
 574       //if( SAVE() ) hitstore[ iRow ] = best;
 575       //std::cout<<"hit search before filter: "<<r.fItr<<", row "<<iRow<<std::endl;
 576       //AliHLTTPCCADisplay::Instance().DrawTracklet(tParam, hitstore, kBlue);
 577       //AliHLTTPCCADisplay::Instance().Ask();
 578 #endif
 579       if ( !tParam.Filter( y, z, err2Y, err2Z, .99 ) ) {
 580         if ( drawSearch ) {
 581           //#ifdef DRAW
 582           std::cout << "tracklet " << r.fItr << " at row " << iRow << " : can not filter!!!! " << std::endl;
 583           //#endif
 584         }
 585         break;
 586       }
 587       if ( SAVE() ) tracklet.SetRowHit( iRow, best );
 588       if ( drawSearch ) {
 589         //#ifdef DRAW
 590         std::cout << "tracklet " << r.fItr << " after filter at row " << iRow << " : " << std::endl;
 591         tParam.Print();
 592 #ifdef DRAW
 593         AliHLTTPCCADisplay::Instance().DrawTracklet( tParam, hitstore, kRed );
 594         AliHLTTPCCADisplay::Instance().Ask();
 595 #endif
 596       }
 597       r.fNHits++;
 598       r.fNMissed = 0;
 599       if ( r.fStage == 1 ) r.fLastRow = iRow;
 600       else r.fFirstRow = iRow;
 601     } while ( 0 );
 602   }
 603 }
 604
 605
 606
 607 GPUd() void AliHLTTPCCATrackletConstructor::Thread
 608 ( int nBlocks, int nThreads, int iBlock, int iThread, int iSync,
 609   AliHLTTPCCASharedMemory &s, AliHLTTPCCAThreadMemory &r, AliHLTTPCCATracker &tracker, AliHLTTPCCATrackParam &tParam )
 610 {
 611
 612   // reconstruction of tracklets
 613   if ( iSync == 0 ) {
 614     Step0( nBlocks, nThreads, iBlock, iThread, s, r, tracker, tParam );
 615   } else if ( iSync == 1 ) {
 616     Step1( nBlocks, nThreads, iBlock, iThread, s, r, tracker, tParam );
 617   } else if ( iSync == 2 ) {
 618     Step2( nBlocks, nThreads, iBlock, iThread, s, r, tracker, tParam );
 619   }
 620
 621   else if ( iSync == 3 )
 622
 623   {
 624     r.fCurrentData = 1;
 625     ReadData( iThread, s, r, tracker, s.fMinStartRow );
 626     r.fCurrentData = 0;
 627     r.fNMissed = 0;
 628   } else if ( iSync == 3 + 159 + 1 ) {
 629     r.fCurrentData = 1;
 630     int nextRow = s.fMaxStartRow - 1;
 631     if ( nextRow < 0 ) nextRow = 0;
 632     ReadData( iThread, s, r, tracker, nextRow );
 633     r.fCurrentData = 0;
 634     r.fNMissed = 0;
 635     r.fStage = 2;
 636     if ( r.fGo ) {
 637       const AliHLTTPCCARow &row = tracker.Row( r.fEndRow );
 638       float x = row.X();
 639       if ( !tParam.TransportToX( x, tracker.Param().Bz(), .999 ) ) r.fGo = 0;
 640     }
 641   }
 642
 643   else if ( iSync <= 3 + 159 + 1 + 159 ) {
 644     int iRow, nextRow;
 645     if (  iSync <= 3 + 159 ) {
 646       iRow = iSync - 4;
 647       if ( iRow < s.fMinStartRow ) return;
 648       nextRow = iRow + 1;
 649       if ( nextRow > 158 ) nextRow = 158;
 650     } else {
 651       iRow = 158 - ( iSync - 4 - 159 - 1 );
 652       if ( iRow >= s.fMaxStartRow ) return;
 653       nextRow = iRow - 1;
 654       if ( nextRow < 0 ) nextRow = 0;
 655     }
 656
 657     if ( r.fIsMemThread ) {
 658       ReadData( iThread, s, r, tracker, nextRow );
 659     } else {
 660       UpdateTracklet( nBlocks, nThreads, iBlock, iThread,
 661                       s, r, tracker, tParam, iRow );
 662     }
 663     r.fCurrentData = !r.fCurrentData;
 664   }
 665
 666   else if ( iSync == 4 + 159*2 + 1 + 1 ) { //
 667     StoreTracklet( nBlocks, nThreads, iBlock, iThread,
 668                    s, r, tracker, tParam );
 669   }
 670 }
 671