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
  25 //#define DRAW
  26
  27 #ifdef DRAW
  28 #include "AliHLTTPCCADisplay.h"
  29 #endif
  30
  31 GPUd() void AliHLTTPCCANeighboursFinder::Thread
  32 ( int /*nBlocks*/, int nThreads, int iBlock, int iThread, int iSync,
  33   AliHLTTPCCASharedMemory &s, AliHLTTPCCATracker &tracker )
  34 {
  35   //* find neighbours
  36
  37   if ( iSync == 0 ) {
  38 #ifdef HLTCA_GPUCODE
  39         for (int i = iThread;i < sizeof(AliHLTTPCCARow) / sizeof(int);i += nThreads)
  40         {
  41                 reinterpret_cast<int*>(&s.fRow)[i] = reinterpret_cast<int*>(&tracker.SliceDataRows()[iBlock])[i];
  42                 if (iBlock >= 2 && iBlock <= tracker.Param().NRows() - 3)
  43                 {
  44                         reinterpret_cast<int*>(&s.fRowUp)[i] = reinterpret_cast<int*>(&tracker.SliceDataRows()[iBlock + 2])[i];
  45                         reinterpret_cast<int*>(&s.fRowDown)[i] = reinterpret_cast<int*>(&tracker.SliceDataRows()[iBlock - 2])[i];
  46                 }
  47         }
  48         __syncthreads();
  49 #endif
  50     if ( iThread == 0 ) {
  51       s.fNRows = tracker.Param().NRows();
  52       s.fIRow = iBlock;
  53       if ( s.fIRow < s.fNRows ) {
  54 #ifdef HLTCA_GPUCODE
  55                 const AliHLTTPCCARow &row = s.fRow;
  56 #else
  57                 const AliHLTTPCCARow &row = tracker.Row( s.fIRow );
  58 #endif
  59         s.fNHits = row.NHits();
  60
  61         if ( ( s.fIRow >= 2 ) && ( s.fIRow <= s.fNRows - 3 ) ) {
  62           s.fIRowUp = s.fIRow + 2;
  63           s.fIRowDn = s.fIRow - 2;
  64
  65           // references to the rows above and below
  66
  67 #ifdef HLTCA_GPUCODE
  68           const AliHLTTPCCARow &rowUp = s.fRowUp;
  69           const AliHLTTPCCARow &rowDn = s.fRowDown;
  70 #else
  71           const AliHLTTPCCARow &rowUp = tracker.Row( s.fIRowUp );
  72           const AliHLTTPCCARow &rowDn = tracker.Row( s.fIRowDn );
  73 #endif
  74           // the axis perpendicular to the rows
  75           const float xDn = rowDn.X();
  76           const float x   = row.X();
  77           const float xUp = rowUp.X();
  78
  79           // number of hits in rows above and below
  80           s.fUpNHits = tracker.Row( s.fIRowUp ).NHits();
  81           s.fDnNHits = tracker.Row( s.fIRowDn ).NHits();
  82
  83           // distance of the rows (absolute and relative)
  84           s.fUpDx = xUp - x;
  85           s.fDnDx = xDn - x;
  86           s.fUpTx = xUp / x;
  87           s.fDnTx = xDn / x;
  88           // UpTx/DnTx is used to move the HitArea such that central events are preferred (i.e. vertices
  89           // coming from y = 0, z = 0).
  90
  91           //s.fGridUp = tracker.Row( s.fIRowUp ).Grid();
  92           //s.fGridDn = tracker.Row( s.fIRowDn ).Grid();
  93         }
  94       }
  95     }
  96   } else if ( iSync == 1 ) {
  97     if ( s.fIRow < s.fNRows ) {
  98       if ( ( s.fIRow <= 1 ) || ( s.fIRow >= s.fNRows - 2 ) ) {
  99 #ifdef HLTCA_GPUCODE
 100                 const AliHLTTPCCARow &row = s.fRow;
 101 #else
 102                 const AliHLTTPCCARow &row = tracker.Row( s.fIRow );
 103 #endif
 104         for ( int ih = iThread; ih < s.fNHits; ih += nThreads ) {
 105           tracker.SetHitLinkUpData( row, ih, -1 );
 106           tracker.SetHitLinkDownData( row, ih, -1 );
 107         }
 108       } else {
 109 /*#ifdef HLTCA_GPUCODE
 110           const AliHLTTPCCARow &rowUp = s.fRowUp;
 111           const AliHLTTPCCARow &rowDn = s.fRowDown;
 112 #else
 113           const AliHLTTPCCARow &rowUp = tracker.Row( s.fIRowUp );
 114           const AliHLTTPCCARow &rowDn = tracker.Row( s.fIRowDn );
 115 #endif
 116
 117         for ( unsigned int ih = iThread; ih < s.fGridUp.N() + s.fGridUp.Ny() + 2; ih += nThreads ) {
 118           s.fGridContentUp[ih] = tracker.FirstHitInBin( rowUp, ih );
 119         }
 120         for ( unsigned int ih = iThread; ih < s.fGridDn.N() + s.fGridDn.Ny() + 2; ih += nThreads ) {
 121           s.fGridContentDn[ih] = tracker.FirstHitInBin( rowDn, ih );
 122         }*/
 123       }
 124     }
 125   } else if ( iSync == 2 ) {
 126     if ( ( s.fIRow <= 1 ) || ( s.fIRow >= s.fNRows - 2 ) ) return;
 127
 128     float chi2Cut = 3.*3.*4 * ( s.fUpDx * s.fUpDx + s.fDnDx * s.fDnDx );
 129     const float kAreaSize = tracker.Param().NeighboursSearchArea();
 130     //float chi2Cut = 3.*3.*(s.fUpDx*s.fUpDx + s.fDnDx*s.fDnDx ); //SG
 131 #define kMaxN 20
 132
 133 #ifdef HLTCA_GPUCODE
 134                   const AliHLTTPCCARow &row = s.fRow;
 135           const AliHLTTPCCARow &rowUp = s.fRowUp;
 136           const AliHLTTPCCARow &rowDn = s.fRowDown;
 137 #else
 138                   const AliHLTTPCCARow &row = tracker.Row( s.fIRow );
 139                   const AliHLTTPCCARow &rowUp = tracker.Row( s.fIRowUp );
 140           const AliHLTTPCCARow &rowDn = tracker.Row( s.fIRowDn );
 141 #endif
 142     const float y0 = row.Grid().YMin();
 143     const float z0 = row.Grid().ZMin();
 144     const float stepY = row.HstepY();
 145     const float stepZ = row.HstepZ();
 146
 147         for ( int ih = iThread; ih < s.fNHits; ih += nThreads ) {
 148
 149       unsigned short *neighUp = s.fB[iThread];
 150       float2 *yzUp = s.fA[iThread];
 151 #if defined(HLTCA_GPUCODE) & kMaxN > ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP
 152           unsigned short neighUp2[kMaxN - ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP];
 153           float2 yzUp2[kMaxN - ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP];
 154 #endif
 155       //unsigned short neighUp[5];
 156       //float2 yzUp[5];
 157
 158       int linkUp = -1;
 159       int linkDn = -1;
 160
 161       if ( s.fDnNHits > 0 && s.fUpNHits > 0 ) {
 162
 163         int nNeighUp = 0;
 164
 165         // coordinates of the hit in the current row
 166 #if defined(HLTCA_GPU_TEXTURE_FETCHa)
 167                 ushort2 tmpval = tex1Dfetch(gAliTexRefu2, ((char*) tracker.Data().HitData() - tracker.pData()->GPUTextureBase()) / sizeof(ushort2) + row.HitNumberOffset() + ih);
 168         const float y = y0 + tmpval.x * stepY;
 169         const float z = z0 + tmpval.y * stepZ;
 170 #else
 171         const float y = y0 + tracker.HitDataY( row, ih ) * stepY;
 172         const float z = z0 + tracker.HitDataZ( row, ih ) * stepZ;
 173 #endif
 174
 175         AliHLTTPCCAHitArea areaDn, areaUp;
 176         // TODO: for NVIDIA GPUs it should use the GridContentUp/-Dn that got copied into shared mem
 177         areaUp.Init( rowUp, tracker.Data(), y*s.fUpTx, z*s.fUpTx, kAreaSize, kAreaSize );
 178         areaDn.Init( rowDn, tracker.Data(), y*s.fDnTx, z*s.fDnTx, kAreaSize, kAreaSize );
 179
 180         do {
 181           AliHLTTPCCAHit h;
 182           int i = areaUp.GetNext( tracker, rowUp, tracker.Data(), &h );
 183           if ( i < 0 ) break;
 184 #if defined(HLTCA_GPUCODE) & kMaxN > ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP
 185                   if (nNeighUp >= ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP)
 186                   {
 187                         neighUp2[nNeighUp - ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP] = ( unsigned short ) i;
 188                         yzUp2[nNeighUp - ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP] = CAMath::MakeFloat2( s.fDnDx * ( h.Y() - y ), s.fDnDx * ( h.Z() - z ) );
 189                   }
 190                   else
 191 #endif
 192                   {
 193                         neighUp[nNeighUp] = ( unsigned short ) i;
 194                         yzUp[nNeighUp] = CAMath::MakeFloat2( s.fDnDx * ( h.Y() - y ), s.fDnDx * ( h.Z() - z ) );
 195                   }
 196           if ( ++nNeighUp >= kMaxN ) break;
 197         } while ( 1 );
 198
 199         int nNeighDn = 0;
 200
 201         if ( nNeighUp > 0 ) {
 202
 203           int bestDn = -1, bestUp = -1;
 204           float bestD = 1.e10;
 205
 206           do {
 207             AliHLTTPCCAHit h;
 208             int i = areaDn.GetNext( tracker, rowDn, tracker.Data(), &h );
 209             if ( i < 0 ) break;
 210
 211             nNeighDn++;
 212             float2 yzdn = CAMath::MakeFloat2( s.fUpDx * ( h.Y() - y ), s.fUpDx * ( h.Z() - z ) );
 213
 214             for ( int iUp = 0; iUp < nNeighUp; iUp++ ) {
 215 #if defined(HLTCA_GPUCODE) & kMaxN > ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP
 216                           float2 yzup = iUp >= ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP ? yzUp2[iUp - ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP] : yzUp[iUp];
 217 #else
 218               float2 yzup = yzUp[iUp];
 219 #endif
 220
 221               float dy = yzdn.x - yzup.x;
 222               float dz = yzdn.y - yzup.y;
 223               float d = dy * dy + dz * dz;
 224               if ( d < bestD ) {
 225                 bestD = d;
 226                 bestDn = i;
 227                 bestUp = iUp;
 228               }
 229             }
 230           } while ( 1 );
 231
 232           if ( bestD <= chi2Cut ) {
 233 #if defined(HLTCA_GPUCODE) & kMaxN > ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP
 234                         linkUp = bestUp >= ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP ? neighUp2[bestUp - ALIHLTTPCCANEIGHBOURS_FINDER_MAX_NNEIGHUP] : neighUp[bestUp];
 235 #else
 236             linkUp = neighUp[bestUp];
 237 #endif
 238             linkDn = bestDn;
 239           }
 240         }
 241 #ifdef DRAW
 242         std::cout << "n NeighUp = " << nNeighUp << ", n NeighDn = " << nNeighDn << std::endl;
 243 #endif
 244
 245       }
 246
 247       tracker.SetHitLinkUpData( row, ih, linkUp );
 248       tracker.SetHitLinkDownData( row, ih, linkDn );
 249 #ifdef DRAW
 250       std::cout << "Links for row " << s.fIRow << ", hit " << ih << ": " << linkUp << " " << linkDn << std::endl;
 251       if ( s.fIRow == 22 && ih == 5 ) {
 252         AliHLTTPCCADisplay::Instance().DrawSliceLink( s.fIRow, ih, -1, -1, 1 );
 253         AliHLTTPCCADisplay::Instance().DrawSliceHit( s.fIRow, ih, kBlue, 1. );
 254         AliHLTTPCCADisplay::Instance().Ask();
 255       }
 256 #endif
 257     }
 258   }
 259 }
 260