HLT/TPCLib/tracking-ca/AliHLTTPCCATrackletConstructor.cxx

   1 // @(#) $Id: AliHLTTPCCATrackletConstructor.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 #include "AliHLTTPCCATracker.h"
  20 #include "AliHLTTPCCATrackParam.h"
  21 #include "AliHLTTPCCATrackParam1.h"
  22 #include "AliHLTTPCCAGrid.h"
  23 #include "AliHLTTPCCAHitArea.h"
  24 #include "AliHLTTPCCAMath.h"
  25 #include "AliHLTTPCCADef.h"
  26 #include "AliHLTTPCCATrackletConstructor.h"
  27
  28 //GPUd() void myprintf1(int i, int j){
  29   //printf("fwd: iS=%d, iRow=%d\n",i,j);
  30 //}
  31 //GPUd() void myprintf2(int i, int j){
  32   //printf("bck: iS=%d, iRow=%d\n",i,j);
  33 //}
  34
  35
  36 GPUd() void AliHLTTPCCATrackletConstructor::Step0
  37 ( Int_t nBlocks, Int_t /*nThreads*/, Int_t iBlock, Int_t iThread, Int_t /*iSync*/,
  38   AliHLTTPCCASharedMemory &s, AliHLTTPCCAThreadMemory &r, AliHLTTPCCATracker &tracker, AliHLTTPCCATrackParam1 &/*tParam*/ )
  39 {
  40   // reconstruction of tracklets, step 0
  41
  42   const Int_t kNMemThreads = 128;
  43
  44   r.fIsMemThread = ( iThread<kNMemThreads );
  45   if( iThread==0 ){
  46     int nTracks = tracker.StartHits()[0];
  47     if(iBlock==0) *tracker.Tracklets() = nTracks;
  48     int nTrPerBlock = nTracks/nBlocks+1;
  49     s.fNRows = tracker.Param().NRows();
  50     s.fItr0 = nTrPerBlock*iBlock;
  51     s.fItr1 = s.fItr0 + nTrPerBlock;
  52     if( s.fItr1> nTracks ) s.fItr1 = nTracks;
  53     s.fUsedHits = tracker.HitIsUsed();
  54     s.fMinStartRow = 158;
  55     s.fMaxStartRow = 0;
  56   }
  57   if( iThread<32 ){
  58     s.fMinStartRow32[iThread] = 158;
  59     s.fMaxStartRow32[iThread] = 0;
  60   }
  61 }
  62
  63
  64 GPUd() void AliHLTTPCCATrackletConstructor::Step1
  65 ( Int_t /*nBlocks*/, Int_t /*nThreads*/, Int_t /*iBlock*/, Int_t iThread, Int_t /*iSync*/,
  66   AliHLTTPCCASharedMemory &s, AliHLTTPCCAThreadMemory &r, AliHLTTPCCATracker &tracker, AliHLTTPCCATrackParam1 &tParam )
  67 {
  68   // reconstruction of tracklets, step 1
  69
  70   const Int_t kNMemThreads = 128;
  71
  72   r.fItr= s.fItr0 + ( iThread - kNMemThreads );
  73   r.fGo = (!r.fIsMemThread) && ( r.fItr<s.fItr1 );
  74   r.fSave = r.fGo;
  75   r.fNHits=0;
  76
  77   if( !r.fGo ) return;
  78
  79   r.fStage = 0;
  80
  81   r.fTrackStoreOffset = 1 + r.fItr*(5+ sizeof(AliHLTTPCCATrackParam)/4 + 160 );
  82   r.fHitStoreOffset = r.fTrackStoreOffset + 5+ sizeof(AliHLTTPCCATrackParam)/4 ;
  83
  84   int *hitstore = tracker.Tracklets() +r.fHitStoreOffset;
  85
  86   UInt_t kThread = iThread %32;//& 00000020;
  87   if( SAVE() ) for( int i=0; i<160; i++ ) hitstore[i] = -1;
  88
  89   int id = tracker.StartHits()[1 + r.fItr];
  90   r.fFirstRow = AliHLTTPCCATracker::ID2IRow(id);
  91   r.fCurrIH =  AliHLTTPCCATracker::ID2IHit(id);
  92   CAMath::atomicMin( &s.fMinStartRow32[kThread], r.fFirstRow);
  93   CAMath::atomicMax( &s.fMaxStartRow32[kThread], r.fFirstRow);
  94   tParam.SinPhi() = 0;
  95   tParam.DzDs() = 0;
  96   tParam.Kappa() = 0;
  97   tParam.CosPhi() = 1;
  98   tParam.Chi2() = 0;
  99   tParam.NDF() = -3;
 100
 101   tParam.Cov()[1] = 0;
 102   tParam.Cov()[4] = 0;
 103   tParam.Cov()[8] = 0;
 104   tParam.Cov()[13] = 0;
 105
 106   tParam.Cov()[3] = 0;
 107   tParam.Cov()[5] = 1.;
 108   tParam.Cov()[7] = 0;
 109   tParam.Cov()[9] = 1.;
 110   tParam.Cov()[10] = 0;
 111   tParam.Cov()[12] = 0;
 112   tParam.Cov()[14] = 1.;
 113   r.fLastRow = r.fFirstRow;
 114 }
 115
 116 GPUd() void AliHLTTPCCATrackletConstructor::Step2
 117 ( Int_t /*nBlocks*/, Int_t /*nThreads*/, Int_t /*iBlock*/, Int_t iThread, Int_t /*iSync*/,
 118   AliHLTTPCCASharedMemory &s, AliHLTTPCCAThreadMemory &/*r*/, AliHLTTPCCATracker &/*tracker*/, AliHLTTPCCATrackParam1 &/*tParam*/ )
 119 {
 120   // reconstruction of tracklets, step 2
 121   if( iThread==0 ){
 122     //CAMath::atomicMinGPU(&s.fMinRow, s.fMinRow32[iThread]);
 123     int minStartRow = 158;
 124     int maxStartRow = 0;
 125     for( int i=0; i<32; i++ ){
 126       if( s.fMinStartRow32[i]<minStartRow ) minStartRow = s.fMinStartRow32[i];
 127       if( s.fMaxStartRow32[i]>maxStartRow ) maxStartRow = s.fMaxStartRow32[i];
 128     }
 129     s.fMinStartRow = minStartRow;
 130     s.fMaxStartRow = maxStartRow;
 131   }
 132 }
 133
 134 GPUd() void AliHLTTPCCATrackletConstructor::ReadData
 135 ( Int_t iThread, AliHLTTPCCASharedMemory &s, AliHLTTPCCAThreadMemory &r, AliHLTTPCCATracker &tracker, Int_t iRow )
 136 {
 137
 138   // reconstruction of tracklets, read data step
 139   const Int_t kNMemThreads = 128;
 140   if( r.fIsMemThread ){
 141     AliHLTTPCCARow &row = tracker.Rows()[iRow];
 142     bool jr = !r.fCurrentData;
 143     Int_t n = row.FullSize();
 144     uint4* gMem = tracker.TexHitsFullData() + row.FullOffset();
 145     uint4 *sMem = s.fData[jr];
 146     for( int i=iThread; i<n; i+=kNMemThreads ) sMem[i] = gMem[i];
 147   }
 148 }
 149
 150
 151 GPUd() void AliHLTTPCCATrackletConstructor::UnpackGrid
 152 ( Int_t /*nBlocks*/, Int_t nThreads, Int_t /*iBlock*/, Int_t iThread, Int_t /*iSync*/,
 153   AliHLTTPCCASharedMemory &s, AliHLTTPCCAThreadMemory &r, AliHLTTPCCATracker &tracker, AliHLTTPCCATrackParam1 &/*tParam*/,Int_t iRow )
 154 {
 155   // reconstruction of tracklets, grid unpacking step
 156
 157   AliHLTTPCCARow &row = tracker.Rows()[iRow];
 158   int n = row.Grid().N()+1;
 159   int nY = row.Grid().Ny();
 160   uint4 *tmpint4 = s.fData[r.fCurrentData];
 161   UShort_t *sGridP = (reinterpret_cast<UShort_t*>(tmpint4)) + row.FullGridOffset();
 162
 163   UInt_t *sGrid = s.fGridContent1;
 164
 165   for( int i=iThread; i<n; i+=nThreads ){
 166     UInt_t s0 = sGridP[i];
 167     UInt_t e0 = sGridP[i+2];
 168     UInt_t s1 = sGridP[i+nY];
 169     UInt_t e1 = sGridP[i+nY+2];
 170     UInt_t nh0 = e0-s0;
 171     UInt_t nh1 = e1-s1;
 172     sGrid[i] = (nh1<<26)+(s1<<16)+( nh0<<10 ) + s0;
 173   }
 174 }
 175
 176
 177 GPUd() void AliHLTTPCCATrackletConstructor::StoreTracklet
 178 ( Int_t /*nBlocks*/, Int_t /*nThreads*/, Int_t /*iBlock*/, Int_t /*iThread*/, Int_t /*iSync*/,
 179   AliHLTTPCCASharedMemory &/*s*/, AliHLTTPCCAThreadMemory &r, AliHLTTPCCATracker &tracker, AliHLTTPCCATrackParam1 &tParam )
 180 {
 181   // reconstruction of tracklets, tracklet store step
 182
 183   if( !r.fSave ) return;
 184
 185   do{
 186     if( r.fNHits<10 ){
 187       r.fNHits = 0;
 188       break;
 189     }
 190
 191     {
 192       Bool_t ok=1;
 193       Float_t *c = tParam.Cov();
 194       for( int i=0; i<15; i++ ) ok = ok && CAMath::Finite(c[i]);
 195       for( int i=0; i<5; i++ ) ok = ok && CAMath::Finite(tParam.Par()[i]);
 196       ok = ok && (tParam.X()>50);
 197
 198       if( c[0]<=0 || c[2]<=0 || c[5]<=0 || c[9]<=0 || c[14]<=0 ) ok = 0;
 199
 200       if(!ok){
 201         r.fNHits = 0;
 202         break;
 203       }
 204     }
 205   }while(0);
 206
 207   if( !SAVE() ) return;
 208
 209   int *store = tracker.Tracklets() + r.fTrackStoreOffset;
 210   int *hitstore = tracker.Tracklets() +r.fHitStoreOffset;
 211   store[0] = r.fNHits;
 212
 213   if( r.fNHits>0 ){
 214     store[3] = r.fFirstRow;
 215     store[4] = r.fLastRow;
 216     if( CAMath::Abs(tParam.Par()[4])<1.e-8 ) tParam.Par()[4] = 1.e-8;
 217     *((AliHLTTPCCATrackParam1*)(store+5)) = tParam;
 218     int w = (r.fNHits<<16)+r.fItr;
 219     for( int iRow=0; iRow<160; iRow++ ){
 220       Int_t ih = hitstore[iRow];
 221       if( ih>=0 ){
 222         int ihTot = tracker.Rows()[iRow].FirstHit() + ih;
 223         CAMath::atomicMax( tracker.HitIsUsed() + ihTot, w );
 224       }
 225     }
 226   }
 227 }
 228
 229 GPUd() void AliHLTTPCCATrackletConstructor::UpdateTracklet
 230 ( Int_t /*nBlocks*/, Int_t /*nThreads*/, Int_t /*iBlock*/, Int_t /*iThread*/, Int_t /*iSync*/,
 231     AliHLTTPCCASharedMemory &s, AliHLTTPCCAThreadMemory &r, AliHLTTPCCATracker &tracker, AliHLTTPCCATrackParam1 &tParam,Int_t iRow )
 232 {
 233   // reconstruction of tracklets, tracklets update step
 234
 235   if( !r.fGo ) return;
 236
 237   const Int_t kMaxRowGap = 5;
 238
 239   int *hitstore = tracker.Tracklets() +r.fHitStoreOffset;
 240
 241   AliHLTTPCCARow &row = tracker.Rows()[iRow];
 242
 243   float y0 = row.Grid().YMin();
 244   float stepY = row.HstepY();
 245   float z0 = row.Grid().ZMin();
 246   float stepZ = row.HstepZ();
 247   float stepYi = row.HstepYi();
 248   float stepZi = row.HstepZi();
 249
 250   if( r.fStage == 0 ){ // fitting part
 251     do{
 252
 253       if( iRow<r.fFirstRow || r.fCurrIH<0  ) break;
 254
 255       uint4 *tmpint4 = s.fData[r.fCurrentData];
 256       ushort2 hh = reinterpret_cast<ushort2*>(tmpint4)[r.fCurrIH];
 257
 258       Int_t oldIH = r.fCurrIH;
 259       r.fCurrIH = reinterpret_cast<Short_t*>(tmpint4)[row.FullLinkOffset() + r.fCurrIH];
 260
 261       float x = row.X();
 262       float y = y0 + hh.x*stepY;
 263       float z = z0 + hh.y*stepZ;
 264
 265       if( iRow==r.fFirstRow ){
 266         tParam.X() = x;
 267         tParam.Y() = y;
 268         tParam.Z() = z;
 269         float err2Y, err2Z;
 270         tracker.GetErrors2( iRow, tParam, err2Y, err2Z );
 271         tParam.Cov()[0] = err2Y;
 272         tParam.Cov()[2] = err2Z;
 273       }else{
 274         if( !tParam.TransportToX0( x, .95 ) ){
 275           if( SAVE() ) hitstore[iRow] = -1;
 276           break;
 277         }
 278         float err2Y, err2Z;
 279         tracker.GetErrors2( iRow, *((AliHLTTPCCATrackParam*)&tParam), err2Y, err2Z );
 280         if( !tParam.Filter2( y, z, err2Y, err2Z, .95 ) ) {
 281           if( SAVE() ) hitstore[iRow] = -1;
 282           break;
 283         }
 284       }
 285       if( SAVE() ) hitstore[iRow] = oldIH;
 286       r.fNHits++;
 287       r.fLastRow = iRow;
 288       if( r.fCurrIH<0 ){
 289         r.fStage = 1;
 290         if( r.fNHits<3 ){ r.fNHits=0; r.fGo = 0;}
 291       }
 292       break;
 293     } while(0);
 294   }
 295   else // forward/backward searching part
 296     {
 297       do{
 298         if( r.fStage == 2 && iRow>=r.fFirstRow ) break;
 299         if( r.fNMissed>kMaxRowGap  ){
 300           r.fGo = 0;
 301           break;
 302         }
 303
 304         r.fNMissed++;
 305
 306         float x = row.X();
 307         float err2Y, err2Z;
 308         if( !tParam.TransportToX0( x, .95 ) ) break;
 309         uint4 *tmpint4 = s.fData[r.fCurrentData];
 310
 311         ushort2 *hits = reinterpret_cast<ushort2*>(tmpint4);
 312         UInt_t *gridContent1 = ((UInt_t*)(s.fGridContent1));
 313
 314         float fY = tParam.GetY();
 315         float fZ = tParam.GetZ();
 316         Int_t best = -1;
 317
 318         { // search for the closest hit
 319
 320           Int_t ds;
 321           Int_t fY0 = (Int_t) ((fY - y0)*stepYi);
 322           Int_t fZ0 = (Int_t) ((fZ - z0)*stepZi);
 323           Int_t ds0 = ( ((int)1)<<30);
 324           ds = ds0;
 325
 326           UInt_t fIndYmin;
 327           UInt_t fHitYfst=1, fHitYlst=0, fHitYfst1=1, fHitYlst1=0;
 328
 329           fIndYmin = row.Grid().GetBin( (float)(fY-1.), (float)(fZ-1.) );
 330           UInt_t c = gridContent1[fIndYmin];
 331           fHitYfst = c & 0x000003FF;
 332           fHitYlst = fHitYfst + ((c & 0x0000FC00)>>10);
 333           fHitYfst1 = ( c & 0x03FF0000 )>>16;
 334           fHitYlst1 = fHitYfst1 + ((c & 0xFC000000)>>26);
 335
 336           for( UInt_t fIh = fHitYfst; fIh<fHitYlst; fIh++ ){
 337             ushort2 hh = hits[fIh];
 338             Int_t ddy = (Int_t)(hh.x) - fY0;
 339             Int_t ddz = (Int_t)(hh.y) - fZ0;
 340             Int_t dds = CAMath::mul24(ddy,ddy) + CAMath::mul24(ddz,ddz);
 341             if( dds<ds ){
 342               ds = dds;
 343               best = fIh;
 344             }
 345           }
 346
 347           for( UInt_t fIh = fHitYfst1; fIh<fHitYlst1; fIh++ ){
 348             ushort2 hh = hits[fIh];
 349             Int_t ddy = (Int_t)(hh.x) - fY0;
 350             Int_t ddz = (Int_t)(hh.y) - fZ0;
 351             Int_t dds = CAMath::mul24(ddy,ddy) + CAMath::mul24(ddz,ddz);
 352             if( dds<ds ){
 353               ds = dds;
 354               best = fIh;
 355             }
 356           }
 357         }// end of search for the closest hit
 358
 359         if( best<0 ) break;
 360
 361         ushort2 hh = hits[best];
 362
 363         tracker.GetErrors2( iRow, *((AliHLTTPCCATrackParam*)&tParam), err2Y, err2Z );
 364
 365         float y = y0 + hh.x*stepY;
 366         float z = z0 + hh.y*stepZ;
 367         float dy = y - fY;
 368         float dz = z - fZ;
 369
 370         const Float_t kFactor = 3.5*3.5;
 371         Float_t sy2 = kFactor*( tParam.GetErr2Y() +  err2Y );
 372         Float_t sz2 = kFactor*( tParam.GetErr2Z() +  err2Z );
 373         if( sy2 > 1. ) sy2 = 1.;
 374         if( sz2 > 1. ) sz2 = 1.;
 375         if( iRow==63 || iRow==64 || iRow==65 ){
 376           if( sy2 < 4. ) sy2 = 4.;
 377           if( sz2 < 4. ) sz2 = 4.;
 378         }
 379
 380
 381         if( CAMath::fmul_rz(dy,dy)>sy2 || CAMath::fmul_rz(dz,dz)>sz2  ) break;
 382
 383         if( !tParam.Filter2( y, z, err2Y, err2Z, .95 ) ) break;
 384
 385         if( SAVE() ) hitstore[ iRow ] = best;
 386         r.fNHits++;
 387         r.fNMissed=0;
 388       }while(0);
 389     }
 390 }
 391
 392
 393
 394 GPUd() void AliHLTTPCCATrackletConstructor::Thread
 395 ( Int_t nBlocks, Int_t nThreads, Int_t iBlock, Int_t iThread, Int_t iSync,
 396   AliHLTTPCCASharedMemory &s, AliHLTTPCCAThreadMemory &r, AliHLTTPCCATracker &tracker, AliHLTTPCCATrackParam1 &tParam )
 397 {
 398
 399   // reconstruction of tracklets
 400   if( iSync==0 )
 401     {
 402       Step0( nBlocks, nThreads, iBlock, iThread, iSync, s, r, tracker, tParam );
 403     }
 404   else if( iSync==1 )
 405     {
 406       Step1( nBlocks, nThreads, iBlock, iThread, iSync, s, r, tracker, tParam );
 407     }
 408   else if( iSync==2 )
 409     {
 410       Step2( nBlocks, nThreads, iBlock, iThread, iSync, s, r, tracker, tParam );
 411     }
 412
 413   else if( iSync==3 )
 414
 415     {
 416       r.fCurrentData = 1;
 417       ReadData( iThread, s, r, tracker, s.fMinStartRow );
 418       r.fCurrentData = 0;
 419       r.fNMissed = 0;
 420     }
 421   else if( iSync==3+159*2+1 )//322
 422
 423     {
 424       r.fCurrentData = 1;
 425       Int_t nextRow = s.fMaxStartRow-1;
 426       if( nextRow<0 ) nextRow = 0;
 427       ReadData( iThread, s, r, tracker, nextRow );
 428       r.fCurrentData = 0;
 429       r.fNMissed = 0;
 430       r.fStage = 2;
 431     }
 432
 433   else if( iSync<=3+159*2+1+159*2 )
 434
 435     {
 436       int iRow, nextRow;
 437       if(  iSync<=3+159*2 ){
 438         iRow = (iSync -4)/2;
 439         //if( iBlock==0 && iThread==0 ) myprintf1(iSync,iRow);
 440         if( iRow < s.fMinStartRow ) return;
 441         nextRow = iRow+1;
 442         if( nextRow>158 ) nextRow = 158;
 443       }else{
 444         iRow = 158 - (iSync - 4-159*2-1)/2;
 445         //if( iBlock==0 && iThread==0 ) myprintf2(iSync,iRow);
 446         if( iRow >= s.fMaxStartRow ) return;
 447         nextRow = iRow-1;
 448         if( nextRow<0 ) nextRow = 0;
 449       }
 450
 451       if( iSync%2==0 ){
 452         UnpackGrid( nBlocks, nThreads, iBlock, iThread, iSync,
 453                     s, r, tracker, tParam, iRow );
 454       }else{
 455         if( r.fIsMemThread ){
 456           ReadData( iThread, s, r, tracker, nextRow );
 457         }else{
 458           UpdateTracklet( nBlocks, nThreads, iBlock, iThread, iSync,
 459                           s, r, tracker, tParam, iRow );
 460         }
 461         r.fCurrentData = !r.fCurrentData;
 462       }
 463     }
 464
 465   else if( iSync== 4+159*4 +1+1 ) // 642
 466
 467     {
 468       StoreTracklet( nBlocks, nThreads, iBlock, iThread, iSync, //SG!!!
 469                      s, r, tracker, tParam );
 470     }
 471 }
 472