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