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