HLT/TPCLib/tracking-ca/AliHLTTPCCANeighboursFinder.cxx

   1 // @(#) $Id: AliHLTTPCCANeighboursFinder1.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 "AliHLTTPCCANeighboursFinder.h"
  21 #include "AliHLTTPCCAMath.h"
  22 #include "AliHLTTPCCAHitArea.h"
  23 #include "AliHLTTPCCATracker.h"
  24 #include "AliHLTTPCCAHit.h"
  25
  26 //#define DRAW
  27
  28 #ifdef DRAW
  29 #include "AliHLTTPCCADisplay.h"
  30 #endif //DRAW
  31
  32 GPUd() void AliHLTTPCCANeighboursFinder::Thread
  33 ( int /*nBlocks*/, int nThreads, int iBlock, int iThread, int iSync,
  34   AliHLTTPCCASharedMemory &s, AliHLTTPCCATracker &tracker )
  35 {
  36   //* find neighbours
  37
  38   if ( iSync == 0 ) {
  39 #ifdef HLTCA_GPUCODE
  40         for (int i = iThread;i < sizeof(AliHLTTPCCARow) / sizeof(int);i += nThreads)
  41         {
  42                 reinterpret_cast<int*>(&s.fRow)[i] = reinterpret_cast<int*>(&tracker.SliceDataRows()[iBlock])[i];
  43                 if (iBlock >= 2 && iBlock <= tracker.Param().NRows() - 3)
  44                 {
  45                         reinterpret_cast<int*>(&s.fRowUp)[i] = reinterpret_cast<int*>(&tracker.SliceDataRows()[iBlock + 2])[i];
  46                         reinterpret_cast<int*>(&s.fRowDown)[i] = reinterpret_cast<int*>(&tracker.SliceDataRows()[iBlock - 2])[i];
  47                 }
  48         }
  49         __syncthreads();
  50 #endif
  51     if ( iThread == 0 ) {
  52       s.fNRows = tracker.Param().NRows();
  53       s.fIRow = iBlock;
  54       if ( s.fIRow < s.fNRows ) {
  55 #ifdef HLTCA_GPUCODE
  56                 const AliHLTTPCCARow &row = s.fRow;
  57 #else
  58                 const AliHLTTPCCARow &row = tracker.Row( s.fIRow );
  59 #endif
  60         s.fNHits = row.NHits();
  61
  62         if ( ( s.fIRow >= 2 ) && ( s.fIRow <= s.fNRows - 3 ) ) {
  63           s.fIRowUp = s.fIRow + 2;
  64           s.fIRowDn = s.fIRow - 2;
  65
  66           // references to the rows above and below
  67
  68 #ifdef HLTCA_GPUCODE
  69           const AliHLTTPCCARow &rowUp = s.fRowUp;
  70           const AliHLTTPCCARow &rowDn = s.fRowDown;
  71 #else
  72           const AliHLTTPCCARow &rowUp = tracker.Row( s.fIRowUp );
  73           const AliHLTTPCCARow &rowDn = tracker.Row( s.fIRowDn );
  74 #endif
  75           // the axis perpendicular to the rows
  76           const float xDn = rowDn.X();
  77           const float x   = row.X();
  78           const float xUp = rowUp.X();
  79
  80           // number of hits in rows above and below
  81           s.fUpNHits = tracker.Row( s.fIRowUp ).NHits();
  82           s.fDnNHits = tracker.Row( s.fIRowDn ).NHits();
  83
  84           // distance of the rows (absolute and relative)
  85           s.fUpDx = xUp - x;
  86           s.fDnDx = xDn - x;
  87           s.fUpTx = xUp / x;
  88           s.fDnTx = xDn / x;
  89           // UpTx/DnTx is used to move the HitArea such that central events are preferred (i.e. vertices
  90           // coming from y = 0, z = 0).
  91
  92           //s.fGridUp = tracker.Row( s.fIRowUp ).Grid();
  93           //s.fGridDn = tracker.Row( s.fIRowDn ).Grid();
  94         }
  95       }
  96     }
  97   } else if ( iSync == 1 ) {
  98     if ( s.fIRow < s.fNRows ) {
  99       if ( ( s.fIRow <= 1 ) || ( s.fIRow >= s.fNRows - 2 ) ) {
 100 #ifdef HLTCA_GPUCODE
 101                 const AliHLTTPCCARow &row = s.fRow;
 102 #else
 103                 const AliHLTTPCCARow &row = tracker.Row( s.fIRow );
 104 #endif
 105         for ( int ih = iThread; ih < s.fNHits; ih += nThreads ) {
 106           tracker.SetHitLinkUpData( row, ih, -1 );
 107           tracker.SetHitLinkDownData( row, ih, -1 );
 108         }
 109       } else {
 110 /*#ifdef HLTCA_GPUCODE
 111           const AliHLTTPCCARow &rowUp = s.fRowUp;
 112           const AliHLTTPCCARow &rowDn = s.fRowDown;
 113 #else
 114           const AliHLTTPCCARow &rowUp = tracker.Row( s.fIRowUp );
 115           const AliHLTTPCCARow &rowDn = tracker.Row( s.fIRowDn );
 116 #endif
 117
 118         for ( unsigned int ih = iThread; ih < s.fGridUp.N() + s.fGridUp.Ny() + 2; ih += nThreads ) {
 119           s.fGridContentUp[ih] = tracker.FirstHitInBin( rowUp, ih );
 120         }
 121         for ( unsigned int ih = iThread; ih < s.fGridDn.N() + s.fGridDn.Ny() + 2; ih += nThreads ) {
 122           s.fGridContentDn[ih] = tracker.FirstHitInBin( rowDn, ih );
 123         }*/
 124       }
 125     }
 126   } else if ( iSync == 2 ) {
 127     if ( ( s.fIRow <= 1 ) || ( s.fIRow >= s.fNRows - 2 ) ) return;
 128
 129     float chi2Cut = 3.*3.*4 * ( s.fUpDx * s.fUpDx + s.fDnDx * s.fDnDx );
 130     const float kAreaSize = tracker.Param().NeighboursSearchArea();
 131     //float chi2Cut = 3.*3.*(s.fUpDx*s.fUpDx + s.fDnDx*s.fDnDx ); //SG
 132 #define kMaxN 20
 133
 134 #ifdef HLTCA_GPUCODE
 135                   const AliHLTTPCCARow &row = s.fRow;
 136           const AliHLTTPCCARow &rowUp = s.fRowUp;
 137           const AliHLTTPCCARow &rowDn = s.fRowDown;
 138 #else
 139                   const AliHLTTPCCARow &row = tracker.Row( s.fIRow );
 140                   const AliHLTTPCCARow &rowUp = tracker.Row( s.fIRowUp );
 141           const AliHLTTPCCARow &rowDn = tracker.Row( s.fIRowDn );
 142 #endif
 143     const float y0 = row.Grid().YMin();
 144     const float z0 = row.Grid().ZMin();
 145     const float stepY = row.HstepY();
 146     const float stepZ = row.HstepZ();
 147
 148         for ( int ih = iThread; ih < s.fNHits; ih += nThreads ) {
 149
 150       int linkUp = -1;
 151       int linkDn = -1;
 152
 153       if ( s.fDnNHits > 0 && s.fUpNHits > 0 ) {
 154
 155
 156
 157         // coordinates of the hit in the current row
 158 #if defined(HLTCA_GPU_TEXTURE_FETCHa)
 159                 ushort2 tmpval = tex1Dfetch(gAliTexRefu2, ((char*) tracker.Data().HitData() - tracker.pData()->GPUTextureBase()) / sizeof(ushort2) + row.HitNumberOffset() + ih);
 160         const float y = y0 + tmpval.x * stepY;
 161         const float z = z0 + tmpval.y * stepZ;
 162 #else
 163         const float y = y0 + tracker.HitDataY( row, ih ) * stepY;
 164         const float z = z0 + tracker.HitDataZ( row, ih ) * stepZ;
 165 #endif
 166
 167 #ifdef FAST_NEIGHBOURS_FINDER
 168
 169 //#define NFDEBUG
 170
 171 #ifdef NFDEBUG
 172         printf("\nSearching Neighbours for: %f %f\n", y, z);
 173 #endif
 174
 175         const float y0Up = rowUp.Grid().YMin();
 176     const float z0Up = rowUp.Grid().ZMin();
 177     const float stepYUp = rowUp.HstepY();
 178     const float stepZUp = rowUp.HstepZ();
 179
 180     const float y0Dn = rowDn.Grid().YMin();
 181     const float z0Dn = rowDn.Grid().ZMin();
 182     const float stepYDn = rowDn.HstepY();
 183     const float stepZDn = rowDn.HstepZ();
 184
 185                 const float yMinUp = y*s.fUpTx - kAreaSize;
 186                 const float yMaxUp = y*s.fUpTx + kAreaSize;
 187                 const float yMinDn = y*s.fDnTx - kAreaSize;
 188                 const float yMaxDn = y*s.fDnTx + kAreaSize;
 189                 const float zMinUp = z*s.fUpTx - kAreaSize;
 190                 const float zMaxUp = z*s.fUpTx + kAreaSize;
 191                 const float zMinDn = z*s.fDnTx - kAreaSize;
 192                 const float zMaxDn = z*s.fDnTx + kAreaSize;
 193
 194                 int bYUpMin, bZUpMin, bYDnMin, bZDnMin;
 195                 rowUp.Grid().GetBin(yMinUp, zMinUp, &bYUpMin, &bZUpMin);
 196                 rowDn.Grid().GetBin(yMinDn, zMinDn, &bYDnMin, &bZDnMin);
 197
 198                 int bYUpMax, bZUpMax, bYDnMax, bZDnMax;
 199                 rowUp.Grid().GetBin(yMaxUp, zMaxUp, &bYUpMax, &bZUpMax);
 200                 rowDn.Grid().GetBin(yMaxDn, zMaxDn, &bYDnMax, &bZDnMax);
 201
 202                 int nYUp = rowUp.Grid().Ny();
 203                 int nYDn = rowDn.Grid().Ny();
 204
 205                 int ihUp = tracker.Data().FirstHitInBin(rowUp, bZUpMin * nYUp + bYUpMin);
 206                 int ihDn = bZDnMax * nYDn + bYDnMax >= rowDn.Grid().N() ? (rowDn.NHits() - 1) : (tracker.Data().FirstHitInBin(rowDn, bZDnMax * nYDn + bYDnMax + 1) - 1);
 207
 208                 int ihUpMax = tracker.Data().FirstHitInBin(rowUp, bZUpMin * nYUp + bYUpMax + 1) - 1;
 209                 int ihDnMin = tracker.Data().FirstHitInBin(rowDn, bZDnMax * nYDn + bYDnMin);
 210
 211                 float bestD = 1.e10;
 212                 int bestDn, bestUp;
 213
 214                 int lastUp = 0, lastDn = 0;
 215
 216                 while (true)
 217                 {
 218                         float yUp = y0Up + tracker.HitDataY(rowUp, ihUp) * stepYUp;
 219                         float zUp = z0Up + tracker.HitDataZ(rowUp, ihUp) * stepZUp;
 220
 221                         float yDn = y0Dn + tracker.HitDataY(rowDn, ihDn) * stepYDn;
 222                         float zDn = z0Dn + tracker.HitDataZ(rowDn, ihDn) * stepZDn;
 223
 224                         int ihUpNext, ihDnNext;
 225                         if (ihUp >= ihUpMax)
 226                         {
 227                                 if (bZUpMin < bZUpMax)
 228                                 {
 229                                         ihUpNext = tracker.Data().FirstHitInBin(rowUp, (bZUpMin + 1) * nYUp + bYUpMin);
 230                                 }
 231                                 else
 232                                 {
 233                                         lastUp = 1;
 234                                 }
 235                         }
 236                         else
 237                         {
 238                                 ihUpNext = ihUp + 1;
 239                         }
 240
 241                         if (ihDn <= ihDnMin)
 242                         {
 243                                 if (bZDnMax > bZDnMin)
 244                                 {
 245                                         ihDnNext = tracker.Data().FirstHitInBin(rowDn, (bZDnMax - 1) * nYDn + bYDnMax);
 246                                 }
 247                                 else
 248                                 {
 249                                         lastDn = 1;
 250                                 }
 251                         }
 252                         else
 253                         {
 254                                 ihDnNext = ihDn - 1;
 255                         }
 256
 257
 258
 259
 260                         float dUp, dDn;
 261                         if (!lastUp)
 262                         {
 263                                 const float yUpNext = y0Up + tracker.HitDataY(rowUp, ihUpNext) * stepYUp;
 264                                 const float zUpNext = z0Up + tracker.HitDataZ(rowUp, ihUpNext) * stepZUp;
 265                                 const float dYUp = s.fUpDx * (yUpNext - y) - s.fDnDx * (yDn - y);
 266                                 const float dZUp = s.fUpDx * (zUpNext - y) - s.fDnDx * (zDn - z);
 267 #ifdef NFDEBUG
 268                                 printf("Checking Up y: %f nexty: %f z: %f nextz: %f\n", yUp, yUpNext, zUp, zUpNext);
 269 #endif
 270                                 dUp = dYUp * dYUp + dZUp * dZUp;
 271                         }
 272
 273                         if (!lastDn)
 274                         {
 275                                 const float yDnNext = y0Dn + tracker.HitDataY(rowDn, ihDnNext) * stepYDn;
 276                                 const float zDnNext = z0Dn + tracker.HitDataZ(rowDn, ihDnNext) * stepZDn;
 277                                 const float dYDn = s.fDnDx * (yDnNext - y) - s.fUpDx * (yUp - y);
 278                                 const float dZDn = s.fDnDx * (zDnNext - y) - s.fUpDx * (zUp - z);
 279 #ifdef NFDEBUG
 280                                 printf("Checking Dn y: %f nexty: %f z: %f nextz: %f\n", yDn, yDnNext, zDn, zDnNext);
 281 #endif
 282                                 dDn = dYDn * dYDn + dZDn * dZDn;
 283                         }
 284
 285                         float d;
 286                         if (lastDn || (dUp < dDn && !lastUp))
 287                         {
 288                                 if (lastUp) break;
 289                                 d = dUp;
 290                                 if (ihUp >= ihUpMax)
 291                                 {
 292                                         bZUpMin++;
 293                                         ihUpMax = tracker.Data().FirstHitInBin(rowUp, bZUpMin * nYUp + bYUpMax + 1) - 1;
 294                                 }
 295                                 ihUp = ihUpNext;
 296                         }
 297                         else
 298                         {
 299                                 d = dDn;
 300                                 if (ihDn <= ihDnMin)
 301                                 {
 302                                         bZDnMax--;
 303                                         ihDnMin = tracker.Data().FirstHitInBin(rowDn, bZDnMax * nYDn + bYDnMin);
 304                                 }
 305                                 ihDn = ihDnNext;
 306                         }
 307                         if (d < bestD)
 308                         {
 309                                 bestD = d;
 310                                 bestUp = ihUp;
 311                                 bestDn = ihDn;
 312                         }
 313                 }
 314
 315                 if (bestD < chi2Cut)
 316                 {
 317                         linkUp = bestUp;
 318                         linkDn = bestDn;
 319                 }
 320
 321 #else
 322                 //Old Slow Neighbours finder
 323       unsigned short *neighUp = s.fB[iThread];
 324       float2 *yzUp = s.fA[iThread];
 325 #if defined(HLTCA_GPUCODE) & kMaxN > ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP
 326           unsigned short neighUp2[kMaxN - ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP];
 327           float2 yzUp2[kMaxN - ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP];
 328 #endif
 329       //unsigned short neighUp[5];
 330       //float2 yzUp[5];
 331
 332                 int nNeighUp = 0;
 333         AliHLTTPCCAHitArea areaDn, areaUp;
 334         areaUp.Init( rowUp, tracker.Data(), y*s.fUpTx, z*s.fUpTx, kAreaSize, kAreaSize );
 335         areaDn.Init( rowDn, tracker.Data(), y*s.fDnTx, z*s.fDnTx, kAreaSize, kAreaSize );
 336
 337         do {
 338           AliHLTTPCCAHit h;
 339           int i = areaUp.GetNext( tracker, rowUp, tracker.Data(), &h );
 340           if ( i < 0 ) break;
 341 #if defined(HLTCA_GPUCODE) & kMaxN > ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP
 342                   if (nNeighUp >= ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP)
 343                   {
 344                         neighUp2[nNeighUp - ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP] = ( unsigned short ) i;
 345                         yzUp2[nNeighUp - ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP] = CAMath::MakeFloat2( s.fDnDx * ( h.Y() - y ), s.fDnDx * ( h.Z() - z ) );
 346                   }
 347                   else
 348 #endif
 349                   {
 350                         neighUp[nNeighUp] = ( unsigned short ) i;
 351                         yzUp[nNeighUp] = CAMath::MakeFloat2( s.fDnDx * ( h.Y() - y ), s.fDnDx * ( h.Z() - z ) );
 352                   }
 353           if ( ++nNeighUp >= kMaxN ) break;
 354         } while ( 1 );
 355
 356         int nNeighDn = 0;
 357
 358         if ( nNeighUp > 0 ) {
 359
 360           int bestDn = -1, bestUp = -1;
 361           float bestD = 1.e10;
 362
 363           do {
 364             AliHLTTPCCAHit h;
 365             int i = areaDn.GetNext( tracker, rowDn, tracker.Data(), &h );
 366             if ( i < 0 ) break;
 367
 368             nNeighDn++;
 369             float2 yzdn = CAMath::MakeFloat2( s.fUpDx * ( h.Y() - y ), s.fUpDx * ( h.Z() - z ) );
 370
 371             for ( int iUp = 0; iUp < nNeighUp; iUp++ ) {
 372 #if defined(HLTCA_GPUCODE) & kMaxN > ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP
 373                           float2 yzup = iUp >= ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP ? yzUp2[iUp - ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP] : yzUp[iUp];
 374 #else
 375               float2 yzup = yzUp[iUp];
 376 #endif
 377
 378               float dy = yzdn.x - yzup.x;
 379               float dz = yzdn.y - yzup.y;
 380               float d = dy * dy + dz * dz;
 381               if ( d < bestD ) {
 382                 bestD = d;
 383                 bestDn = i;
 384                 bestUp = iUp;
 385               }
 386             }
 387           } while ( 1 );
 388
 389           if ( bestD <= chi2Cut ) {
 390 #if defined(HLTCA_GPUCODE) & kMaxN > ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP
 391                         linkUp = bestUp >= ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP ? neighUp2[bestUp - ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP] : neighUp[bestUp];
 392 #else
 393             linkUp = neighUp[bestUp];
 394 #endif
 395             linkDn = bestDn;
 396           }
 397         }
 398 #ifdef DRAW
 399         std::cout << "n NeighUp = " << nNeighUp << ", n NeighDn = " << nNeighDn << std::endl;
 400 #endif //DRAW
 401 #endif //FAST_NEIGHBOURS_FINDER
 402       }
 403
 404       tracker.SetHitLinkUpData( row, ih, linkUp );
 405       tracker.SetHitLinkDownData( row, ih, linkDn );
 406 #ifdef DRAW
 407       std::cout << "Links for row " << s.fIRow << ", hit " << ih << ": " << linkUp << " " << linkDn << std::endl;
 408       if ( s.fIRow == 22 && ih == 5 ) {
 409         AliHLTTPCCADisplay::Instance().DrawSliceLink( s.fIRow, ih, -1, -1, 1 );
 410         AliHLTTPCCADisplay::Instance().DrawSliceHit( s.fIRow, ih, kBlue, 1. );
 411         AliHLTTPCCADisplay::Instance().Ask();
 412       }
 413 #endif //DRAW
 414     }
 415   }
 416 }
 417