HLT/TPCLib/tracking-ca/AliHLTTPCCASliceData.cxx

   1 // **************************************************************************
   2 // * This file is property of and copyright by the ALICE HLT Project        *
   3 // * All rights reserved.                                                   *
   4 // *                                                                        *
   5 // * Primary Authors:                                                       *
   6 // *     Copyright 2009       Matthias Kretz <kretz@kde.org>                *
   7 // *                                                                        *
   8 // * Permission to use, copy, modify and distribute this software and its   *
   9 // * documentation strictly for non-commercial purposes is hereby granted   *
  10 // * without fee, provided that the above copyright notice appears in all   *
  11 // * copies and that both the copyright notice and this permission notice   *
  12 // * appear in the supporting documentation. The authors make no claims     *
  13 // * about the suitability of this software for any purpose. It is          *
  14 // * provided "as is" without express or implied warranty.                  *
  15 // **************************************************************************
  16
  17 #include "AliHLTTPCCASliceData.h"
  18 #include "AliHLTTPCCAClusterData.h"
  19 #include "AliHLTTPCCAMath.h"
  20 #include "AliHLTArray.h"
  21 #include "AliHLTTPCCAHit.h"
  22 #include "AliHLTTPCCAParam.h"
  23 #include "MemoryAssignmentHelpers.h"
  24 #include "AliHLTTPCCAGPUConfig.h"
  25 #include "AliHLTTPCCAGPUTracker.h"
  26 #include <iostream>
  27
  28 // calculates an approximation for 1/sqrt(x)
  29 // Google for 0x5f3759df :)
  30 static inline float fastInvSqrt( float _x )
  31 {
  32   // the function calculates fast inverse sqrt
  33
  34   union { float f; int i; } x = { _x };
  35   const float xhalf = 0.5f * x.f;
  36   x.i = 0x5f3759df - ( x.i >> 1 );
  37   x.f = x.f * ( 1.5f - xhalf * x.f * x.f );
  38   return x.f;
  39 }
  40
  41 inline void AliHLTTPCCASliceData::CreateGrid( AliHLTTPCCARow *row, const AliHLTTPCCAClusterData &data, int ClusterDataHitNumberOffset )
  42 {
  43   // grid creation
  44
  45   if ( row->NHits() <= 0 ) { // no hits or invalid data
  46     // grid coordinates don't matter, since there are no hits
  47     row->fGrid.CreateEmpty();
  48     return;
  49   }
  50
  51   float yMin =  1.e3f;
  52   float yMax = -1.e3f;
  53   float zMin =  1.e3f;
  54   float zMax = -1.e3f;
  55   for ( int i = ClusterDataHitNumberOffset; i < ClusterDataHitNumberOffset + row->fNHits; ++i ) {
  56     const float y = data.Y( i );
  57     const float z = data.Z( i );
  58     if ( yMax < y ) yMax = y;
  59     if ( yMin > y ) yMin = y;
  60     if ( zMax < z ) zMax = z;
  61     if ( zMin > z ) zMin = z;
  62   }
  63
  64   const float norm = fastInvSqrt( row->fNHits );
  65   row->fGrid.Create( yMin, yMax, zMin, zMax,
  66                      CAMath::Max( ( yMax - yMin ) * norm, 2.f ),
  67                      CAMath::Max( ( zMax - zMin ) * norm, 2.f ) );
  68 }
  69
  70 inline void AliHLTTPCCASliceData::PackHitData( AliHLTTPCCARow* const row, const AliHLTArray<AliHLTTPCCAHit> &binSortedHits )
  71 {
  72   // hit data packing
  73
  74   static const float shortPackingConstant = 1.f / 65535.f;
  75   const float y0 = row->fGrid.YMin();
  76   const float z0 = row->fGrid.ZMin();
  77   const float stepY = ( row->fGrid.YMax() - y0 ) * shortPackingConstant;
  78   const float stepZ = ( row->fGrid.ZMax() - z0 ) * shortPackingConstant;
  79   const float stepYi = 1.f / stepY;
  80   const float stepZi = 1.f / stepZ;
  81
  82   row->fHy0 = y0;
  83   row->fHz0 = z0;
  84   row->fHstepY = stepY;
  85   row->fHstepZ = stepZ;
  86   row->fHstepYi = stepYi;
  87   row->fHstepZi = stepZi;
  88
  89   for ( int hitIndex = 0; hitIndex < row->fNHits; ++hitIndex ) {
  90     // bin sorted index!
  91     const int globalHitIndex = row->fHitNumberOffset + hitIndex;
  92     const AliHLTTPCCAHit &hh = binSortedHits[globalHitIndex];
  93     const float xx = ( ( hh.Y() - y0 ) * stepYi ) + .5 ;
  94     const float yy = ( ( hh.Z() - z0 ) * stepZi ) + .5 ;
  95     if ( xx < 0 || yy < 0 || xx >= 65536  || yy >= 65536 ) {
  96       std::cout << "!!!! hit packing error!!! " << xx << " " << yy << " " << std::endl;
  97     }
  98     // HitData is bin sorted
  99     fHitData[row->fHitNumberOffset + hitIndex].x = xx;
 100     fHitData[row->fHitNumberOffset + hitIndex].y = yy;
 101   }
 102 }
 103
 104 void AliHLTTPCCASliceData::Clear()
 105 {
 106   fNumberOfHits = 0;
 107 }
 108
 109 void AliHLTTPCCASliceData::InitializeRows( const AliHLTTPCCAParam &p )
 110 {
 111   // initialisation of rows
 112         if (!fRows) fRows = new AliHLTTPCCARow[HLTCA_ROW_COUNT + 1];
 113   for ( int i = 0; i < p.NRows(); ++i ) {
 114     fRows[i].fX = p.RowX( i );
 115     fRows[i].fMaxY = CAMath::Tan( p.DAlpha() / 2. ) * fRows[i].fX;
 116   }
 117 }
 118
 119 #ifndef HLTCA_GPUCODE
 120         AliHLTTPCCASliceData::~AliHLTTPCCASliceData()
 121         {
 122                 //Standard Destrcutor
 123                 if (fRows)
 124                 {
 125                         if (!fIsGpuSliceData) delete[] fRows;
 126                         fRows = NULL;
 127                 }
 128                 if (fMemory)
 129                 {
 130                         if (!fIsGpuSliceData) delete[] fMemory;
 131                         fMemory = NULL;
 132                 }
 133
 134         }
 135 #endif
 136
 137 GPUh() void AliHLTTPCCASliceData::SetGPUSliceDataMemory(void* const pSliceMemory, void* const pRowMemory)
 138 {
 139         //Set Pointer to slice data memory to external memory
 140         fMemory = (char*) pSliceMemory;
 141         fRows = (AliHLTTPCCARow*) pRowMemory;
 142 }
 143
 144 size_t AliHLTTPCCASliceData::SetPointers(const AliHLTTPCCAClusterData *data, bool allocate)
 145 {
 146         //Set slice data internal pointers
 147   int hitMemCount = 0;
 148   for ( int rowIndex = data->FirstRow(); rowIndex <= data->LastRow(); ++rowIndex )
 149   {
 150         hitMemCount += NextMultipleOf<sizeof(HLTCA_GPU_ROWALIGNMENT) / sizeof(ushort_v)>(data->NumberOfClusters( rowIndex ));
 151   }
 152         //Calculate Memory needed to store hits in rows
 153
 154   const int numberOfRows = data->LastRow() - data->FirstRow() + 1;
 155   const unsigned int kVectorAlignment = 256 /*sizeof( uint4 )*/ ;
 156   fNumberOfHitsPlusAlign = NextMultipleOf < ( kVectorAlignment > sizeof(HLTCA_GPU_ROWALIGNMENT) ? kVectorAlignment : sizeof(HLTCA_GPU_ROWALIGNMENT)) / sizeof( int ) > ( hitMemCount );
 157   fNumberOfHits = data->NumberOfClusters();
 158   const int firstHitInBinSize = (23 + sizeof(HLTCA_GPU_ROWALIGNMENT) / sizeof(int)) * numberOfRows + 4 * fNumberOfHits + 3;
 159   //FIXME: sizeof(HLTCA_GPU_ROWALIGNMENT) / sizeof(int) * numberOfRows is way to big and only to ensure to reserve enough memory for GPU Alignment.
 160   //Might be replaced by correct value
 161
 162   const int memorySize =
 163     // LinkData, HitData
 164     fNumberOfHitsPlusAlign * 4 * sizeof( short ) +
 165     // FirstHitInBin
 166     NextMultipleOf<kVectorAlignment>( ( firstHitInBinSize ) * sizeof( int ) ) +
 167     // HitWeights, ClusterDataIndex
 168     fNumberOfHitsPlusAlign * 2 * sizeof( int );
 169
 170   if ( 1||fMemorySize < memorySize ) { // release the memory on CPU
 171         fMemorySize = memorySize;
 172         if (allocate && !fIsGpuSliceData)
 173         {
 174                 if (fMemory)
 175                 {
 176                         delete[] fMemory;
 177                 }
 178           fMemory = new char[fMemorySize + 4];// kVectorAlignment];
 179         }
 180   }
 181
 182   char *mem = fMemory;
 183   AssignMemory( fLinkUpData,   mem, fNumberOfHitsPlusAlign );
 184   AssignMemory( fLinkDownData, mem, fNumberOfHitsPlusAlign );
 185   AssignMemory( fHitData,     mem, fNumberOfHitsPlusAlign );
 186   AssignMemory( fFirstHitInBin,  mem, firstHitInBinSize );
 187   fGpuMemorySize = mem - fMemory;
 188
 189   //Memory Allocated below will not be copied to GPU but instead be initialized on the gpu itself. Therefore it must not be copied to GPU!
 190   AssignMemory( fHitWeights,   mem, fNumberOfHitsPlusAlign );
 191   AssignMemory( fClusterDataIndex, mem, fNumberOfHitsPlusAlign );
 192   return(mem - fMemory);
 193 }
 194
 195 void AliHLTTPCCASliceData::InitFromClusterData( const AliHLTTPCCAClusterData &data )
 196 {
 197   // initialisation from cluster data
 198
 199   ////////////////////////////////////
 200   // 1. prepare arrays
 201   ////////////////////////////////////
 202
 203   const int numberOfRows = data.LastRow() - data.FirstRow() + 1;
 204   fNumberOfHits = data.NumberOfClusters();
 205
 206   /* TODO Vectorization
 207   for ( int rowIndex = data.FirstRow(); rowIndex <= data.LastRow(); ++rowIndex ) {
 208     int NumberOfClusters( int rowIndex ) const;
 209   }
 210   const int memorySize = fNumberOfHits * sizeof( short_v::Type )
 211   */
 212   SetPointers(&data, true);
 213
 214   ////////////////////////////////////
 215   // 2. fill HitData and FirstHitInBin
 216   ////////////////////////////////////
 217
 218   for ( int rowIndex = 0; rowIndex < data.FirstRow(); ++rowIndex ) {
 219     AliHLTTPCCARow &row = fRows[rowIndex];
 220     row.fGrid.CreateEmpty();
 221     row.fNHits = 0;
 222     row.fFullSize = 0;
 223     row.fHitNumberOffset = 0;
 224     row.fFirstHitInBinOffset = 0;
 225
 226     row.fHy0 = 0.f;
 227     row.fHz0 = 0.f;
 228     row.fHstepY = 1.f;
 229     row.fHstepZ = 1.f;
 230     row.fHstepYi = 1.f;
 231     row.fHstepZi = 1.f;
 232   }
 233   for ( int rowIndex = data.LastRow() + 1; rowIndex < HLTCA_ROW_COUNT + 1; ++rowIndex ) {
 234     AliHLTTPCCARow &row = fRows[rowIndex];
 235     row.fGrid.CreateEmpty();
 236     row.fNHits = 0;
 237     row.fFullSize = 0;
 238     row.fHitNumberOffset = 0;
 239     row.fFirstHitInBinOffset = 0;
 240
 241     row.fHy0 = 0.f;
 242     row.fHz0 = 0.f;
 243     row.fHstepY = 1.f;
 244     row.fHstepZ = 1.f;
 245     row.fHstepYi = 1.f;
 246     row.fHstepZi = 1.f;
 247   }
 248
 249
 250   AliHLTResizableArray<AliHLTTPCCAHit> binSortedHits( fNumberOfHits + sizeof(HLTCA_GPU_ROWALIGNMENT) / sizeof(ushort_v) * numberOfRows + 1 );
 251
 252   int gridContentOffset = 0;
 253   int hitOffset = 0;
 254
 255   int binCreationMemorySize = 103 * 2 + fNumberOfHits;
 256   AliHLTResizableArray<unsigned short> binCreationMemory( binCreationMemorySize );
 257
 258   fGPUSharedDataReq = 0;
 259
 260   for ( int rowIndex = data.FirstRow(); rowIndex <= data.LastRow(); ++rowIndex ) {
 261     AliHLTTPCCARow &row = fRows[rowIndex];
 262     row.fNHits = data.NumberOfClusters( rowIndex );
 263     assert( row.fNHits < ( 1 << sizeof( unsigned short ) * 8 ) );
 264         row.fHitNumberOffset = hitOffset;
 265         hitOffset += NextMultipleOf<sizeof(HLTCA_GPU_ROWALIGNMENT) / sizeof(ushort_v)>(data.NumberOfClusters( rowIndex ));
 266
 267     row.fFirstHitInBinOffset = gridContentOffset;
 268
 269     CreateGrid( &row, data, data.RowOffset( rowIndex ) );
 270     const AliHLTTPCCAGrid &grid = row.fGrid;
 271     const int numberOfBins = grid.N();
 272
 273     int binCreationMemorySizeNew;
 274     if ( ( binCreationMemorySizeNew = numberOfBins * 2 + 6 + row.fNHits + sizeof(HLTCA_GPU_ROWALIGNMENT) / sizeof(unsigned short) * numberOfRows + 1) > binCreationMemorySize ) {
 275       binCreationMemorySize = binCreationMemorySizeNew;
 276       binCreationMemory.Resize( binCreationMemorySize );
 277     }
 278
 279     AliHLTArray<unsigned short> c = binCreationMemory;           // number of hits in all previous bins
 280     AliHLTArray<unsigned short> bins = c + ( numberOfBins + 3 ); // cache for the bin index for every hit in this row
 281     AliHLTArray<unsigned short> filled = bins + row.fNHits;      // counts how many hits there are per bin
 282
 283     for ( unsigned int bin = 0; bin < row.fGrid.N() + 3; ++bin ) {
 284       filled[bin] = 0; // initialize filled[] to 0
 285     }
 286
 287     for ( int hitIndex = 0; hitIndex < row.fNHits; ++hitIndex ) {
 288       const int globalHitIndex = data.RowOffset( rowIndex ) + hitIndex;
 289       const unsigned short bin = row.fGrid.GetBin( data.Y( globalHitIndex ), data.Z( globalHitIndex ) );
 290
 291       bins[hitIndex] = bin;
 292       ++filled[bin];
 293     }
 294
 295     unsigned short n = 0;
 296     for ( int bin = 0; bin < numberOfBins + 3; ++bin ) {
 297       c[bin] = n;
 298       n += filled[bin];
 299     }
 300
 301     for ( int hitIndex = 0; hitIndex < row.fNHits; ++hitIndex ) {
 302       const unsigned short bin = bins[hitIndex];
 303       --filled[bin];
 304       const unsigned short ind = c[bin] + filled[bin]; // generate an index for this hit that is >= c[bin] and < c[bin + 1]
 305       const int globalBinsortedIndex = row.fHitNumberOffset + ind;
 306       const int globalHitIndex = data.RowOffset( rowIndex ) + hitIndex;
 307
 308       // allows to find the global hit index / coordinates from a global bin sorted hit index
 309       fClusterDataIndex[globalBinsortedIndex] = globalHitIndex;
 310       binSortedHits[globalBinsortedIndex].SetY( data.Y( globalHitIndex ) );
 311       binSortedHits[globalBinsortedIndex].SetZ( data.Z( globalHitIndex ) );
 312     }
 313
 314     PackHitData( &row, binSortedHits );
 315
 316     for ( int i = 0; i < numberOfBins; ++i ) {
 317       fFirstHitInBin[row.fFirstHitInBinOffset + i] = c[i]; // global bin-sorted hit index
 318     }
 319     const unsigned short a = c[numberOfBins];
 320     // grid.N is <= row.fNHits
 321     const int nn = numberOfBins + grid.Ny() + 3;
 322     for ( int i = numberOfBins; i < nn; ++i ) {
 323       assert( (signed) row.fFirstHitInBinOffset + i < 23 * numberOfRows + 4 * fNumberOfHits + 3 );
 324       fFirstHitInBin[row.fFirstHitInBinOffset + i] = a;
 325     }
 326
 327     row.fFullSize = nn;
 328     gridContentOffset += nn;
 329
 330         if (NextMultipleOf<sizeof(HLTCA_GPU_ROWALIGNMENT) / sizeof(ushort_v)>(row.fNHits) + nn > (unsigned) fGPUSharedDataReq)
 331                 fGPUSharedDataReq = NextMultipleOf<sizeof(HLTCA_GPU_ROWALIGNMENT) / sizeof(ushort_v)>(row.fNHits) + nn;
 332
 333         //Make pointer aligned
 334         gridContentOffset = NextMultipleOf<sizeof(HLTCA_GPU_ROWALIGNMENT) / sizeof(ushort_v)>(gridContentOffset);
 335   }
 336
 337 #if 0
 338   //SG cell finder - test code
 339
 340   if ( fTmpHitInputIDs ) delete[] fTmpHitInputIDs;
 341   fTmpHitInputIDs = new int [NHits];
 342   const float areaY = .5;
 343   const float areaZ = .5;
 344   int newRowNHitsTotal = 0;
 345   bool *usedHits = new bool [NHits];
 346   for ( int iHit = 0; iHit < NHits; iHit++ ) usedHits[iHit] = 0;
 347   for ( int iRow = 0; iRow < fParam.NRows(); iRow++ ) {
 348     rowHeaders[iRow*2  ] = newRowNHitsTotal; // new first hit
 349     rowHeaders[iRow*2+1] = 0; // new N hits
 350     int newRowNHits = 0;
 351     int oldRowFirstHit = RowFirstHit[iRow];
 352     int oldRowLastHit = oldRowFirstHit + RowNHits[iRow];
 353     for ( int iHit = oldRowFirstHit; iHit < oldRowLastHit; iHit++ ) {
 354       if ( usedHits[iHit] ) continue;
 355       float x0 = X[iHit];
 356       float y0 = Y[iHit];
 357       float z0 = Z[iHit];
 358       float cx = x0;
 359       float cy = y0;
 360       float cz = z0;
 361       int nclu = 1;
 362       usedHits[iHit] = 1;
 363       if ( 0 ) for ( int jHit = iHit + 1; jHit < oldRowLastHit; jHit++ ) {//SG!!!
 364           //if( usedHits[jHit] ) continue;
 365           float dy = Y[jHit] - y0;
 366           float dz = Z[jHit] - z0;
 367           if ( CAMath::Abs( dy ) < areaY && CAMath::Abs( dz ) < areaZ ) {
 368             cx += X[jHit];
 369             cy += Y[jHit];
 370             cz += Z[jHit];
 371             nclu++;
 372             usedHits[jHit] = 1;
 373           }
 374         }
 375       int id = newRowNHitsTotal + newRowNHits;
 376       hitsXYZ[id*3+0 ] = cx / nclu;
 377       hitsXYZ[id*3+1 ] = cy / nclu;
 378       hitsXYZ[id*3+2 ] = cz / nclu;
 379       fTmpHitInputIDs[id] = iHit;
 380       newRowNHits++;
 381     }
 382     rowHeaders[iRow*2+1] = newRowNHits;
 383     newRowNHitsTotal += newRowNHits;
 384   }
 385   NHitsTotal() = newRowNHitsTotal;
 386   reinterpret_cast<int*>( fInputEvent )[1+fParam.NRows()*2] = newRowNHitsTotal;
 387
 388   delete[] usedHits;
 389 #endif
 390 }
 391
 392 void AliHLTTPCCASliceData::ClearHitWeights()
 393 {
 394   // clear hit weights
 395
 396 #ifdef ENABLE_VECTORIZATION
 397   const int_v v0( Zero );
 398   const int *const end = fHitWeights + fNumberOfHits;
 399   for ( int *mem = fHitWeights; mem < end; mem += v0.Size ) {
 400     v0.store( mem );
 401   }
 402 #else
 403   for ( int i = 0; i < fNumberOfHitsPlusAlign; ++i ) {
 404     fHitWeights[i] = 0;
 405   }
 406 #endif
 407 }
 408
 409 void AliHLTTPCCASliceData::ClearLinks()
 410 {
 411   // link cleaning
 412
 413 #ifdef ENABLE_VECTORIZATION
 414   const short_v v0( -1 );
 415   const short *const end1 = fLinkUpData + fNumberOfHits;
 416   for ( short *mem = fLinkUpData; mem < end; mem += v0.Size ) {
 417     v0.store( mem );
 418   }
 419   const short *const end2 = fLinkDownData + fNumberOfHits;
 420   for ( short *mem = fLinkDownData; mem < end; mem += v0.Size ) {
 421     v0.store( mem );
 422   }
 423 #else
 424   for ( int i = 0; i < fNumberOfHits; ++i ) {
 425     fLinkUpData[i] = -1;
 426   }
 427   for ( int i = 0; i < fNumberOfHits; ++i ) {
 428     fLinkDownData[i] = -1;
 429   }
 430 #endif
 431 }
 432