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 <iostream>
  25
  26 // calculates an approximation for 1/sqrt(x)
  27 // Google for 0x5f3759df :)
  28 static inline float fastInvSqrt( float _x )
  29 {
  30   union { float f; int i; } x = { _x };
  31   const float xhalf = 0.5f * x.f;
  32   x.i = 0x5f3759df - ( x.i >> 1 );
  33   x.f = x.f * ( 1.5f - xhalf * x.f * x.f );
  34   return x.f;
  35 }
  36
  37 inline void AliHLTTPCCASliceData::createGrid( AliHLTTPCCARow *row, const AliHLTTPCCAClusterData &data )
  38 {
  39   if ( row->NHits() <= 0 ) { // no hits or invalid data
  40     // grid coordinates don't matter, since there are no hits
  41     row->fGrid.CreateEmpty();
  42     return;
  43   }
  44
  45   float yMin =  1.e3f;
  46   float yMax = -1.e3f;
  47   float zMin =  1.e3f;
  48   float zMax = -1.e3f;
  49   for ( int i = row->fHitNumberOffset; i < row->fHitNumberOffset + row->fNHits; ++i ) {
  50     const float y = data.Y( i );
  51     const float z = data.Z( i );
  52     if ( yMax < y ) yMax = y;
  53     if ( yMin > y ) yMin = y;
  54     if ( zMax < z ) zMax = z;
  55     if ( zMin > z ) zMin = z;
  56   }
  57
  58   const float norm = fastInvSqrt( row->fNHits );
  59   row->fGrid.Create( yMin, yMax, zMin, zMax,
  60                      CAMath::Max( ( yMax - yMin ) * norm, 2.f ),
  61                      CAMath::Max( ( zMax - zMin ) * norm, 2.f ) );
  62 }
  63
  64 inline void AliHLTTPCCASliceData::PackHitData( AliHLTTPCCARow *row, const AliHLTArray<AliHLTTPCCAHit> &binSortedHits )
  65 {
  66   static const float shortPackingConstant = 1.f / 65535.f;
  67   const float y0 = row->fGrid.YMin();
  68   const float z0 = row->fGrid.ZMin();
  69   const float stepY = ( row->fGrid.YMax() - y0 ) * shortPackingConstant;
  70   const float stepZ = ( row->fGrid.ZMax() - z0 ) * shortPackingConstant;
  71   const float stepYi = 1.f / stepY;
  72   const float stepZi = 1.f / stepZ;
  73
  74   row->fHy0 = y0;
  75   row->fHz0 = z0;
  76   row->fHstepY = stepY;
  77   row->fHstepZ = stepZ;
  78   row->fHstepYi = stepYi;
  79   row->fHstepZi = stepZi;
  80
  81   for ( int hitIndex = 0; hitIndex < row->fNHits; ++hitIndex ) {
  82     // bin sorted index!
  83     const int globalHitIndex = row->fHitNumberOffset + hitIndex;
  84     const AliHLTTPCCAHit &hh = binSortedHits[globalHitIndex];
  85     const float xx = ( ( hh.Y() - y0 ) * stepYi ) + .5 ;
  86     const float yy = ( ( hh.Z() - z0 ) * stepZi ) + .5 ;
  87     if ( xx < 0 || yy < 0 || xx >= 65536  || yy >= 65536 ) {
  88       std::cout << "!!!! hit packing error!!! " << xx << " " << yy << " " << std::endl;
  89     }
  90     // HitData is bin sorted
  91     fHitDataY[row->fHitNumberOffset + hitIndex] = xx;
  92     fHitDataZ[row->fHitNumberOffset + hitIndex] = yy;
  93   }
  94 }
  95
  96 void AliHLTTPCCASliceData::Clear()
  97 {
  98   fNumberOfHits = 0;
  99 }
 100
 101 void AliHLTTPCCASliceData::InitializeRows( const AliHLTTPCCAParam &p )
 102 {
 103   for ( int i = 0; i < p.NRows(); ++i ) {
 104     fRows[i].fX = p.RowX( i );
 105     fRows[i].fMaxY = CAMath::Tan( p.DAlpha() / 2. ) * fRows[i].fX;
 106   }
 107 }
 108
 109 void AliHLTTPCCASliceData::InitFromClusterData( const AliHLTTPCCAClusterData &data )
 110 {
 111   ////////////////////////////////////
 112   // 1. prepare arrays
 113   ////////////////////////////////////
 114
 115   fNumberOfHits = data.NumberOfClusters();
 116
 117   /* TODO Vectorization
 118   for ( int rowIndex = data.FirstRow(); rowIndex <= data.LastRow(); ++rowIndex ) {
 119     int NumberOfClusters( int rowIndex ) const;
 120   }
 121   const int memorySize = fNumberOfHits * sizeof( short_v::Type )
 122   */
 123   const int numberOfRows = data.LastRow() - data.FirstRow();
 124   enum { VectorAlignment = sizeof( int ) };
 125   const int numberOfHitsPlusAlignment = NextMultipleOf < VectorAlignment / sizeof( int ) > ( fNumberOfHits );
 126   const int memorySize =
 127     // LinkData, HitData
 128     numberOfHitsPlusAlignment * 4 * sizeof( short ) +
 129     // FirstHitInBin
 130     NextMultipleOf<VectorAlignment>( ( 23 * numberOfRows + 4 * fNumberOfHits ) * sizeof( int ) ) +
 131     // HitWeights, ClusterDataIndex
 132     numberOfHitsPlusAlignment * 2 * sizeof( int );
 133
 134   if ( fMemorySize < memorySize ) {
 135     fMemorySize = memorySize;
 136     delete[] fMemory;
 137     fMemory = new char[fMemorySize + 4];// VectorAlignment];
 138   }
 139
 140   char *mem = fMemory;
 141   AssignMemory( fLinkUpData,   mem, numberOfHitsPlusAlignment );
 142   AssignMemory( fLinkDownData, mem, numberOfHitsPlusAlignment );
 143   AssignMemory( fHitDataY,     mem, numberOfHitsPlusAlignment );
 144   AssignMemory( fHitDataZ,     mem, numberOfHitsPlusAlignment );
 145   AssignMemory( fFirstHitInBin,  mem, 23 * numberOfRows + 4 * fNumberOfHits );
 146   AssignMemory( fHitWeights,   mem, numberOfHitsPlusAlignment );
 147   AssignMemory( fClusterDataIndex, mem, numberOfHitsPlusAlignment );
 148
 149   ////////////////////////////////////
 150   // 2. fill HitData and FirstHitInBin
 151   ////////////////////////////////////
 152
 153   for ( int rowIndex = 0; rowIndex < data.FirstRow(); ++rowIndex ) {
 154     AliHLTTPCCARow &row = fRows[rowIndex];
 155     row.fGrid.CreateEmpty();
 156     row.fNHits = 0;
 157     row.fFullSize = 0;
 158     row.fHitNumberOffset = 0;
 159     row.fFirstHitInBinOffset = 0;
 160
 161     row.fHy0 = 0.f;
 162     row.fHz0 = 0.f;
 163     row.fHstepY = 1.f;
 164     row.fHstepZ = 1.f;
 165     row.fHstepYi = 1.f;
 166     row.fHstepZi = 1.f;
 167   }
 168   for ( int rowIndex = data.LastRow() + 1; rowIndex < 160; ++rowIndex ) {
 169     AliHLTTPCCARow &row = fRows[rowIndex];
 170     row.fGrid.CreateEmpty();
 171     row.fNHits = 0;
 172     row.fFullSize = 0;
 173     row.fHitNumberOffset = 0;
 174     row.fFirstHitInBinOffset = 0;
 175
 176     row.fHy0 = 0.f;
 177     row.fHz0 = 0.f;
 178     row.fHstepY = 1.f;
 179     row.fHstepZ = 1.f;
 180     row.fHstepYi = 1.f;
 181     row.fHstepZi = 1.f;
 182   }
 183
 184   AliHLTResizableArray<AliHLTTPCCAHit> binSortedHits( fNumberOfHits );
 185
 186   int gridContentOffset = 0;
 187
 188   int binCreationMemorySize = 103 * 2 + fNumberOfHits;
 189   AliHLTResizableArray<unsigned short> binCreationMemory( binCreationMemorySize );
 190
 191   for ( int rowIndex = data.FirstRow(); rowIndex <= data.LastRow(); ++rowIndex ) {
 192     AliHLTTPCCARow &row = fRows[rowIndex];
 193
 194     row.fNHits = data.NumberOfClusters( rowIndex );
 195     assert( row.fNHits < ( 1 << sizeof( unsigned short ) * 8 ) );
 196     row.fHitNumberOffset = data.RowOffset( rowIndex );
 197     row.fFirstHitInBinOffset = gridContentOffset;
 198
 199     createGrid( &row, data );
 200     const AliHLTTPCCAGrid &grid = row.fGrid;
 201     const int numberOfBins = grid.N();
 202
 203     int binCreationMemorySizeNew;
 204     if ( ( binCreationMemorySizeNew = numberOfBins * 2 + 6 + row.fNHits ) > binCreationMemorySize ) {
 205       binCreationMemorySize = binCreationMemorySizeNew;
 206       binCreationMemory.Resize( binCreationMemorySize );
 207     }
 208
 209     AliHLTArray<unsigned short> c = binCreationMemory;           // number of hits in all previous bins
 210     AliHLTArray<unsigned short> bins = c + ( numberOfBins + 3 ); // cache for the bin index for every hit in this row
 211     AliHLTArray<unsigned short> filled = bins + row.fNHits;      // counts how many hits there are per bin
 212
 213     for ( unsigned int bin = 0; bin < row.fGrid.N() + 3; ++bin ) {
 214       filled[bin] = 0; // initialize filled[] to 0
 215     }
 216
 217     for ( int hitIndex = 0; hitIndex < row.fNHits; ++hitIndex ) {
 218       const int globalHitIndex = row.fHitNumberOffset + hitIndex;
 219       const unsigned short bin = row.fGrid.GetBin( data.Y( globalHitIndex ), data.Z( globalHitIndex ) );
 220       bins[hitIndex] = bin;
 221       ++filled[bin];
 222     }
 223
 224     unsigned short n = 0;
 225     for ( int bin = 0; bin < numberOfBins + 3; ++bin ) {
 226       c[bin] = n;
 227       n += filled[bin];
 228     }
 229
 230     for ( int hitIndex = 0; hitIndex < row.fNHits; ++hitIndex ) {
 231       const unsigned short bin = bins[hitIndex];
 232       --filled[bin];
 233       const unsigned short ind = c[bin] + filled[bin]; // generate an index for this hit that is >= c[bin] and < c[bin + 1]
 234       const int globalBinsortedIndex = row.fHitNumberOffset + ind;
 235       const int globalHitIndex = row.fHitNumberOffset + hitIndex;
 236
 237       // allows to find the global hit index / coordinates from a global bin sorted hit index
 238       fClusterDataIndex[globalBinsortedIndex] = globalHitIndex;
 239       binSortedHits[globalBinsortedIndex].SetY( data.Y( globalHitIndex ) );
 240       binSortedHits[globalBinsortedIndex].SetZ( data.Z( globalHitIndex ) );
 241     }
 242
 243     PackHitData( &row, binSortedHits );
 244
 245     for ( int i = 0; i < numberOfBins; ++i ) {
 246       fFirstHitInBin[row.fFirstHitInBinOffset + i] = c[i]; // global bin-sorted hit index
 247     }
 248     const unsigned short a = c[numberOfBins];
 249     // grid.N is <= row.fNHits
 250     const int nn = numberOfBins + grid.Ny() + 3;
 251     for ( int i = numberOfBins; i < nn; ++i ) {
 252       assert( row.fFirstHitInBinOffset + i < 23 * numberOfRows + 4 * fNumberOfHits );
 253       fFirstHitInBin[row.fFirstHitInBinOffset + i] = a;
 254     }
 255
 256     row.fFullSize = nn;
 257     gridContentOffset += nn;
 258   }
 259
 260 #if 0
 261   //SG cell finder - test code
 262
 263   if ( fTmpHitInputIDs ) delete[] fTmpHitInputIDs;
 264   fTmpHitInputIDs = new int [NHits];
 265   const float areaY = .5;
 266   const float areaZ = .5;
 267   int newRowNHitsTotal = 0;
 268   bool *usedHits = new bool [NHits];
 269   for ( int iHit = 0; iHit < NHits; iHit++ ) usedHits[iHit] = 0;
 270   for ( int iRow = 0; iRow < fParam.NRows(); iRow++ ) {
 271     rowHeaders[iRow*2  ] = newRowNHitsTotal; // new first hit
 272     rowHeaders[iRow*2+1] = 0; // new N hits
 273     int newRowNHits = 0;
 274     int oldRowFirstHit = RowFirstHit[iRow];
 275     int oldRowLastHit = oldRowFirstHit + RowNHits[iRow];
 276     for ( int iHit = oldRowFirstHit; iHit < oldRowLastHit; iHit++ ) {
 277       if ( usedHits[iHit] ) continue;
 278       float x0 = X[iHit];
 279       float y0 = Y[iHit];
 280       float z0 = Z[iHit];
 281       float cx = x0;
 282       float cy = y0;
 283       float cz = z0;
 284       int nclu = 1;
 285       usedHits[iHit] = 1;
 286       if ( 0 ) for ( int jHit = iHit + 1; jHit < oldRowLastHit; jHit++ ) {//SG!!!
 287           //if( usedHits[jHit] ) continue;
 288           float dy = Y[jHit] - y0;
 289           float dz = Z[jHit] - z0;
 290           if ( CAMath::Abs( dy ) < areaY && CAMath::Abs( dz ) < areaZ ) {
 291             cx += X[jHit];
 292             cy += Y[jHit];
 293             cz += Z[jHit];
 294             nclu++;
 295             usedHits[jHit] = 1;
 296           }
 297         }
 298       int id = newRowNHitsTotal + newRowNHits;
 299       hitsXYZ[id*3+0 ] = cx / nclu;
 300       hitsXYZ[id*3+1 ] = cy / nclu;
 301       hitsXYZ[id*3+2 ] = cz / nclu;
 302       fTmpHitInputIDs[id] = iHit;
 303       newRowNHits++;
 304     }
 305     rowHeaders[iRow*2+1] = newRowNHits;
 306     newRowNHitsTotal += newRowNHits;
 307   }
 308   NHitsTotal() = newRowNHitsTotal;
 309   reinterpret_cast<int*>( fInputEvent )[1+fParam.NRows()*2] = newRowNHitsTotal;
 310
 311   delete[] usedHits;
 312 #endif
 313 }
 314
 315 void AliHLTTPCCASliceData::ClearHitWeights()
 316 {
 317 #ifdef ENABLE_VECTORIZATION
 318   const int_v v0( Zero );
 319   const int *const end = fHitWeights + fNumberOfHits;
 320   for ( int *mem = fHitWeights; mem < end; mem += v0.Size ) {
 321     v0.store( mem );
 322   }
 323 #else
 324   for ( int i = 0; i < fNumberOfHits; ++i ) {
 325     fHitWeights[i] = 0;
 326   }
 327 #endif
 328 }
 329
 330 void AliHLTTPCCASliceData::ClearLinks()
 331 {
 332 #ifdef ENABLE_VECTORIZATION
 333   const short_v v0( -1 );
 334   const short *const end1 = fLinkUpData + fNumberOfHits;
 335   for ( short *mem = fLinkUpData; mem < end; mem += v0.Size ) {
 336     v0.store( mem );
 337   }
 338   const short *const end2 = fLinkDownData + fNumberOfHits;
 339   for ( short *mem = fLinkDownData; mem < end; mem += v0.Size ) {
 340     v0.store( mem );
 341   }
 342 #else
 343   for ( int i = 0; i < fNumberOfHits; ++i ) {
 344     fLinkUpData[i] = -1;
 345   }
 346   for ( int i = 0; i < fNumberOfHits; ++i ) {
 347     fLinkDownData[i] = -1;
 348   }
 349 #endif
 350 }
 351