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