HLT/TPCLib/tracking-ca/AliHLTTPCCATracker.cxx

   1 // @(#) $Id$
   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 "AliHLTTPCCATracker.h"
  21 #include "AliHLTTPCCAOutTrack.h"
  22 #include "AliHLTTPCCAGrid.h"
  23 #include "AliHLTTPCCARow.h"
  24 #include "AliHLTTPCCATrack.h"
  25 #include "AliHLTTPCCATracklet.h"
  26 #include "AliHLTTPCCAMath.h"
  27 #include "AliHLTTPCCAHit.h"
  28 #include "MemoryAssignmentHelpers.h"
  29
  30 #include "TStopwatch.h"
  31 #include "AliHLTTPCCAHitArea.h"
  32 #include "AliHLTTPCCANeighboursFinder.h"
  33 #include "AliHLTTPCCANeighboursCleaner.h"
  34 #include "AliHLTTPCCAStartHitsFinder.h"
  35 #include "AliHLTTPCCATrackletConstructor.h"
  36 #include "AliHLTTPCCATrackletSelector.h"
  37 #include "AliHLTTPCCAProcess.h"
  38 #include "AliHLTTPCCASliceTrack.h"
  39 #include "AliHLTTPCCASliceOutput.h"
  40 #include "AliHLTTPCCADataCompressor.h"
  41 #include "AliHLTTPCCAClusterData.h"
  42
  43 #include "AliHLTTPCCATrackParam.h"
  44
  45 #include "AliHLTTPCCAGPUConfig.h"
  46 #include "AliHLTTPCCAGPUTracker.h"
  47
  48 #if !defined(HLTCA_GPUCODE)
  49 #include <iostream>
  50 #include <string.h>
  51 #endif
  52
  53 //#define DRAW1
  54
  55 #ifdef DRAW1
  56 #include "AliHLTTPCCADisplay.h"
  57 #endif //DRAW
  58
  59 #ifdef HLTCA_INTERNAL_PERFORMANCE
  60 #include "AliHLTTPCCAPerformance.h"
  61 #endif
  62
  63 #ifdef HLTCA_STANDALONE
  64 #include "AliHLTTPCCAStandaloneFramework.h"
  65 #endif
  66
  67 ClassImp( AliHLTTPCCATracker )
  68
  69
  70 #if !defined(HLTCA_GPUCODE)
  71
  72 AliHLTTPCCATracker::~AliHLTTPCCATracker()
  73 {
  74   // destructor
  75         if (!fIsGPUTracker)
  76         {
  77                 if (fCommonMem) delete fCommonMem;
  78                 if (fHitMemory) delete[] fHitMemory;
  79                 if (fTrackletMemory) delete[] fTrackletMemory;
  80                 if (fTrackMemory) delete[] fTrackMemory;
  81                 fCommonMem = NULL;
  82                 fHitMemory = fTrackMemory = NULL;
  83         }
  84 }
  85
  86 // ----------------------------------------------------------------------------------
  87 void AliHLTTPCCATracker::Initialize( const AliHLTTPCCAParam &param )
  88 {
  89   // initialisation
  90   fParam = param;
  91   fParam.Update();
  92   fData.InitializeRows( fParam );
  93
  94   StartEvent();
  95 }
  96
  97 void AliHLTTPCCATracker::StartEvent()
  98 {
  99   // start new event and fresh the memory
 100
 101   SetupCommonMemory();
 102 }
 103
 104 void AliHLTTPCCATracker::SetGPUTracker()
 105 {
 106         fIsGPUTracker = true;
 107         fData.SetGpuSliceData();
 108 }
 109
 110 char* AliHLTTPCCATracker::SetGPUTrackerCommonMemory(char* pGPUMemory)
 111 {
 112         fCommonMem = (commonMemoryStruct*) pGPUMemory;
 113         return(pGPUMemory + sizeof(commonMemoryStruct));
 114 }
 115
 116
 117 char* AliHLTTPCCATracker::SetGPUTrackerHitsMemory(char* pGPUMemory, int MaxNHits)
 118 {
 119         fHitMemory = (char*) pGPUMemory;
 120         SetPointersHits(MaxNHits);
 121         pGPUMemory += fHitMemorySize;
 122         AssignMemory(fTrackletTmpStartHits, pGPUMemory, NHitsTotal());
 123         AssignMemory(fRowStartHitCountOffset, pGPUMemory, Param().NRows());
 124
 125         return(pGPUMemory);
 126 }
 127
 128 char* AliHLTTPCCATracker::SetGPUTrackerTrackletsMemory(char* pGPUMemory, int MaxNTracks)
 129 {
 130         fTrackletMemory = (char*) pGPUMemory;
 131         SetPointersTracklets(MaxNTracks);
 132         pGPUMemory += fTrackletMemorySize;
 133         AssignMemory(fGPUTrackletTemp, pGPUMemory, MaxNTracks);
 134         AssignMemory(fRowBlockTracklets, pGPUMemory, MaxNTracks * 2 * (Param().NRows() / HLTCA_GPU_SCHED_ROW_STEP + 1));
 135         AssignMemory(fRowBlockPos, pGPUMemory, 2 * (Param().NRows() / HLTCA_GPU_SCHED_ROW_STEP + 1));
 136         AssignMemory(fBlockStartingTracklet, pGPUMemory, HLTCA_GPU_BLOCK_COUNT);
 137
 138         return(pGPUMemory);
 139 }
 140
 141 char* AliHLTTPCCATracker::SetGPUTrackerTracksMemory(char* pGPUMemory, int MaxNTracks, int MaxNHits )
 142 {
 143         fTrackMemory = (char*) pGPUMemory;
 144         SetPointersTracks(MaxNTracks, MaxNHits);
 145         pGPUMemory += fTrackMemorySize;
 146
 147         return(pGPUMemory);
 148 }
 149
 150 void AliHLTTPCCATracker::DumpSliceData(std::ostream &out)
 151 {
 152         out << "Slice Data:" << endl;
 153         for (int i = 0;i < Param().NRows();i++)
 154         {
 155                 out << "Row: " << i << endl;
 156                 for (int j = 0;j < Row(i).NHits();j++)
 157                 {
 158                         if (j && j % 16 == 0) out << std::endl;
 159                         out << j << '-' << Data().HitDataY(Row(i), j) << '-' << Data().HitDataZ(Row(i), j) << ", ";
 160                 }
 161                 out << endl;
 162         }
 163 }
 164
 165 void AliHLTTPCCATracker::DumpLinks(std::ostream &out)
 166 {
 167         out << "Hit Links:" << endl;
 168         for (int i = 0;i < Param().NRows();i++)
 169         {
 170                 out << "Row: " << i << endl;
 171                 for (int j = 0;j < Row(i).NHits();j++)
 172                 {
 173                         if (j && j % 32 == 0) out << endl;
 174                         out << HitLinkUpData(Row(i), j) << "/" << HitLinkDownData(Row(i), j) << ", ";
 175                 }
 176                 out << endl;
 177         }
 178 }
 179
 180 void AliHLTTPCCATracker::DumpHitWeights(std::ostream &out)
 181 {
 182     out << "Hit Weights:" << endl;
 183     for (int i = 0;i < Param().NRows();i++)
 184     {
 185         out << "Row: " << i << ":" << endl;
 186         for (int j = 0;j < Row(i).NHits();j++)
 187         {
 188             if (j && j % 32 == 0) out << endl;
 189             out << HitWeight(Row(i), j) << ", ";
 190         }
 191         out << endl;
 192     }
 193 }
 194
 195 int AliHLTTPCCATracker::starthitSortComparison(const void*a, const void* b)
 196 {
 197         AliHLTTPCCAHitId* aa = (AliHLTTPCCAHitId*) a;
 198         AliHLTTPCCAHitId* bb = (AliHLTTPCCAHitId*) b;
 199
 200         if (aa->RowIndex() != bb->RowIndex()) return(aa->RowIndex() - bb->RowIndex());
 201         return(aa->HitIndex() - bb->HitIndex());
 202 }
 203
 204 void AliHLTTPCCATracker::DumpStartHits(std::ostream &out)
 205 {
 206         out << "Start Hits: (" << *NTracklets() << ")" << endl;
 207 #ifdef HLTCA_GPU_SORT_DUMPDATA
 208         qsort(TrackletStartHits(), *NTracklets(), sizeof(AliHLTTPCCAHitId), starthitSortComparison);
 209 #endif
 210         for (int i = 0;i < *NTracklets();i++)
 211         {
 212                 out << TrackletStartHit(i).RowIndex() << "-" << TrackletStartHit(i).HitIndex() << endl;
 213         }
 214         out << endl;
 215 }
 216
 217 void AliHLTTPCCATracker::DumpTrackHits(std::ostream &out)
 218 {
 219         out << "Tracks: (" << *NTracks() << ")" << endl;
 220 #ifdef HLTCA_GPU_SORT_DUMPDATA
 221         for (int k = 0;k < Param().NRows();k++)
 222         {
 223                 for (int l = 0;l < Row(k).NHits();l++)
 224                 {
 225 #endif
 226                         for (int j = 0;j < *NTracks();j++)
 227                         {
 228                                 if (Tracks()[j].NHits() == 0 || !Tracks()[j].Alive()) continue;
 229 #ifdef HLTCA_GPU_SORT_DUMPDATA
 230                                 if (TrackHits()[Tracks()[j].FirstHitID()].RowIndex() == k && TrackHits()[Tracks()[j].FirstHitID()].HitIndex() == l)
 231                                 {
 232 #endif
 233                                         for (int i = 0;i < Tracks()[j].NHits();i++)
 234                                         {
 235                                                 out << TrackHits()[Tracks()[j].FirstHitID() + i].RowIndex() << "-" << TrackHits()[Tracks()[j].FirstHitID() + i].HitIndex() << ", ";
 236                                         }
 237                                         out << "(Track: " << j << ")" << endl;
 238 #ifdef HLTCA_GPU_SORT_DUMPDATA
 239                                 }
 240                         }
 241 #endif
 242                 }
 243 #ifdef HLTCA_GPU_SORT_DUMPDATA
 244         }
 245 #endif
 246 }
 247
 248 void AliHLTTPCCATracker::DumpTrackletHits(std::ostream &out)
 249 {
 250         out << "Tracklets: (" << *NTracklets() << ")" << endl;
 251 #ifdef HLTCA_GPU_SORT_DUMPDATA
 252         AliHLTTPCCAHitId* tmpIds = new AliHLTTPCCAHitId[*NTracklets()];
 253         AliHLTTPCCATracklet* tmpTracklets = new AliHLTTPCCATracklet[*NTracklets()];
 254         memcpy(tmpIds, TrackletStartHits(), *NTracklets() * sizeof(AliHLTTPCCAHitId));
 255         memcpy(tmpTracklets, Tracklets(), *NTracklets() * sizeof(AliHLTTPCCATracklet));
 256 #ifdef EXTERN_ROW_HITS
 257         int* tmpHits = new int[*NTracklets() * Param().NRows()];
 258         memcpy(tmpHits, TrackletRowHits(), *NTracklets() * Param().NRows() * sizeof(int));
 259 #endif
 260         qsort(TrackletStartHits(), *NTracklets(), sizeof(AliHLTTPCCAHitId), starthitSortComparison);
 261         for (int i = 0;i < *NTracklets();i++)
 262         {
 263                 for (int j = 0;j < *NTracklets();j++)
 264                 {
 265                         if (tmpIds[i].RowIndex() == TrackletStartHit(j).RowIndex() && tmpIds[i].HitIndex() == TrackletStartHit(j).HitIndex())
 266                         {
 267                                 memcpy(&Tracklets()[j], &tmpTracklets[i], sizeof(AliHLTTPCCATracklet));
 268 #ifdef EXTERN_ROW_HITS
 269                                 if (tmpTracklets[i].NHits())
 270                                 {
 271                                         for (int k = tmpTracklets[i].FirstRow();k <= tmpTracklets[i].LastRow();k++)
 272                                         {
 273                                                 fTrackletRowHits[k * *NTracklets() + j] = tmpHits[k * *NTracklets() + i];
 274                                         }
 275                                 }
 276 #endif
 277                                 break;
 278                         }
 279                 }
 280         }
 281         delete[] tmpIds;
 282         delete[] tmpTracklets;
 283 #ifdef EXTERN_ROW_HITS
 284         delete[] tmpHits;
 285 #endif
 286 #endif
 287         for (int j = 0;j < *NTracklets();j++)
 288         {
 289                 out << "Tracklet " << j << " (Hits: " << setw(3) << Tracklets()[j].NHits() << ", Start: " << setw(3) << TrackletStartHit(j).RowIndex() << "-" << setw(3) << TrackletStartHit(j).HitIndex() << ") ";
 290                 if (Tracklets()[j].NHits() == 0);
 291                 else if (Tracklets()[j].LastRow() > Tracklets()[j].FirstRow() && (Tracklets()[j].FirstRow() >= Param().NRows() || Tracklets()[j].LastRow() >= Param().NRows()))
 292                 {
 293 #ifdef HLTCA_STANDALONE
 294                         printf("\nError: First %d Last %d Num %d", Tracklets()[j].FirstRow(), Tracklets()[j].LastRow(), Tracklets()[j].NHits());
 295 #endif
 296                 }
 297                 else if (Tracklets()[j].NHits() && Tracklets()[j].LastRow() > Tracklets()[j].FirstRow())
 298                 {
 299                         for (int i = Tracklets()[j].FirstRow();i <= Tracklets()[j].LastRow();i++)
 300                         {
 301                                 //if (Tracklets()[j].RowHit(i) != -1)
 302 #ifdef EXTERN_ROW_HITS
 303                                         out << i << "-" << fTrackletRowHits[i * fCommonMem->fNTracklets + j] << ", ";
 304 #else
 305                                         out << i << "-" << Tracklets()[j].RowHit(i) << ", ";
 306 #endif
 307                         }
 308                 }
 309                 out << endl;
 310         }
 311 }
 312
 313
 314 void AliHLTTPCCATracker::SetupCommonMemory()
 315 {
 316   // set up common memory
 317
 318   if (!fIsGPUTracker)
 319   {
 320     if ( !fCommonMem ) {
 321       // the 1600 extra bytes are not used unless fCommonMemorySize increases with a later event
 322       //fCommonMemory = reinterpret_cast<char*> ( new uint4 [ fCommonMemorySize/sizeof( uint4 ) + 100] );
 323           fCommonMem = new commonMemoryStruct;
 324     }
 325
 326         if (fHitMemory) delete[] fHitMemory;
 327         if (fTrackletMemory) delete[] fTrackletMemory;
 328         if (fTrackMemory) delete[] fTrackMemory;
 329   }
 330
 331   fHitMemory = fTrackletMemory = fTrackMemory = 0;
 332
 333   fData.Clear();
 334   fCommonMem->fNTracklets = 0;
 335   fCommonMem->fNTracks = 0 ;
 336   fCommonMem->fNTrackHits = 0;
 337 }
 338
 339 void AliHLTTPCCATracker::ReadEvent( AliHLTTPCCAClusterData *clusterData )
 340 {
 341   // read event
 342
 343   fClusterData = clusterData;
 344
 345   StartEvent();
 346
 347   //* Convert input hits, create grids, etc.
 348   fData.InitFromClusterData( *clusterData );
 349   {
 350     if (!fIsGPUTracker)
 351         {
 352                 SetPointersHits( fData.NumberOfHits() ); // to calculate the size
 353                 fHitMemory = reinterpret_cast<char*> ( new uint4 [ fHitMemorySize/sizeof( uint4 ) + 100] );
 354         }
 355     SetPointersHits( fData.NumberOfHits() ); // set pointers for hits
 356   }
 357 }
 358
 359 GPUhd() void  AliHLTTPCCATracker::SetPointersHits( int MaxNHits )
 360 {
 361   // set all pointers to the event memory
 362
 363   char *mem = fHitMemory;
 364
 365   // extra arrays for tpc clusters
 366
 367 #ifdef HLTCA_GPU_SORT_STARTHITS_2
 368   AssignMemory( fTrackletStartHits, mem, MaxNHits + 32);
 369 #else
 370   AssignMemory( fTrackletStartHits, mem, MaxNHits);
 371 #endif
 372
 373   // calculate the size
 374
 375   fHitMemorySize = mem - fHitMemory;
 376 }
 377
 378 GPUhd() void  AliHLTTPCCATracker::SetPointersTracklets( int MaxNTracklets )
 379 {
 380   // set all pointers to the tracklets memory
 381   char *mem = fTrackletMemory;
 382
 383   // memory for tracklets
 384
 385   AssignMemory( fTracklets, mem, MaxNTracklets );
 386 #ifdef EXTERN_ROW_HITS
 387   AssignMemory( fTrackletRowHits, mem, MaxNTracklets * Param().NRows());
 388 #endif
 389
 390   fTrackletMemorySize = mem - fTrackletMemory;
 391 }
 392
 393
 394 GPUhd() void  AliHLTTPCCATracker::SetPointersTracks( int MaxNTracks, int MaxNHits )
 395 {
 396   // set all pointers to the tracks memory
 397   char *mem = fTrackMemory;
 398
 399   // memory for selected tracks
 400
 401   AssignMemory( fTracks, mem, MaxNTracks );
 402   AssignMemory( fTrackHits, mem, 2 * MaxNHits );
 403
 404   // calculate the size
 405
 406   fTrackMemorySize = mem - fTrackMemory;
 407 }
 408
 409 GPUh() int AliHLTTPCCATracker::CheckEmptySlice()
 410 {
 411   if ( NHitsTotal() < 1 ) {
 412     {
 413           AliHLTTPCCASliceOutput::Allocate(*fOutput, 0, 0, fOutputControl);
 414           AliHLTTPCCASliceOutput* useOutput = *fOutput;
 415           if (fOutput == NULL) return(1);
 416       useOutput->SetNTracks( 0 );
 417       useOutput->SetNTrackClusters( 0 );
 418           useOutput->SetNOutTracks(0);
 419           useOutput->SetNOutTrackHits(0);
 420     }
 421
 422     return 1;
 423   }
 424   return 0;
 425 }
 426
 427 void AliHLTTPCCATracker::RunNeighboursFinder()
 428 {
 429         AliHLTTPCCAProcess<AliHLTTPCCANeighboursFinder>( Param().NRows(), 1, *this );
 430 }
 431
 432 void AliHLTTPCCATracker::RunNeighboursCleaner()
 433 {
 434         AliHLTTPCCAProcess<AliHLTTPCCANeighboursCleaner>( Param().NRows() - 2, 1, *this );
 435 }
 436
 437 void AliHLTTPCCATracker::RunStartHitsFinder()
 438 {
 439         AliHLTTPCCAProcess<AliHLTTPCCAStartHitsFinder>( Param().NRows() - 4, 1, *this );
 440 }
 441
 442 void AliHLTTPCCATracker::RunTrackletConstructor()
 443 {
 444     AliHLTTPCCATrackletConstructor::AliHLTTPCCATrackletConstructorNewCPU(*this);
 445 }
 446
 447 void AliHLTTPCCATracker::RunTrackletSelector()
 448 {
 449   AliHLTTPCCAProcess<AliHLTTPCCATrackletSelector>( 1, fCommonMem->fNTracklets, *this );
 450 }
 451
 452 #ifdef HLTCA_STANDALONE
 453 void AliHLTTPCCATracker::StandalonePerfTime(int i)
 454 {
 455   if (fGPUDebugLevel >= 1)
 456   {
 457           AliHLTTPCCAStandaloneFramework::StandaloneQueryTime(&fPerfTimers[i]);
 458   }
 459 }
 460 #else
 461 void AliHLTTPCCATracker::StandalonePerfTime(int /*i*/) {}
 462 #endif
 463
 464 GPUh() void AliHLTTPCCATracker::Reconstruct()
 465 {
 466   //* reconstruction of event
 467   //std::cout<<"Reconstruct slice "<<fParam.ISlice()<<", nHits="<<NHitsTotal()<<std::endl;
 468
 469   fTimers[0] = 0; // find neighbours
 470   fTimers[1] = 0; // construct tracklets
 471   fTimers[2] = 0; // fit tracklets
 472   fTimers[3] = 0; // prolongation of tracklets
 473   fTimers[4] = 0; // selection
 474   fTimers[5] = 0; // write output
 475   fTimers[6] = 0;
 476   fTimers[7] = 0;
 477
 478   //if( fParam.ISlice()<1 ) return; //SG!!!
 479
 480   TStopwatch timer0;
 481
 482   StandalonePerfTime(0);
 483
 484   if (CheckEmptySlice()) return;
 485
 486 #ifdef DRAW1
 487   //if( fParam.ISlice()==2 || fParam.ISlice()==3)
 488     {
 489   AliHLTTPCCADisplay::Instance().ClearView();
 490   AliHLTTPCCADisplay::Instance().SetSliceView();
 491   AliHLTTPCCADisplay::Instance().SetCurrentSlice( this );
 492   AliHLTTPCCADisplay::Instance().DrawSlice( this, 1 );
 493   if ( NHitsTotal() > 0 ) {
 494     AliHLTTPCCADisplay::Instance().DrawSliceHits( kRed, .5 );
 495     AliHLTTPCCADisplay::Instance().Ask();
 496   }
 497   }
 498 #endif
 499
 500         fCommonMem->fNTracklets = fCommonMem->fNTracks = fCommonMem->fNTrackHits = 0;
 501
 502 #if !defined(HLTCA_GPUCODE)
 503
 504   /*if (fGPUDebugLevel >= 6)
 505   {
 506           *fGPUDebugOut << endl << endl << "Slice: " << Param().ISlice() << endl;
 507           *fGPUDebugOut << "Slice Data:" << endl;
 508           DumpSliceData(*fGPUDebugOut);
 509   }*/
 510
 511   StandalonePerfTime(1);
 512
 513   RunNeighboursFinder();
 514
 515   StandalonePerfTime(2);
 516
 517   if (fGPUDebugLevel >= 6) DumpLinks(*fGPUDebugOut);
 518
 519 #ifdef HLTCA_INTERNAL_PERFORMANCE
 520   //if( Param().ISlice()<=2 )
 521   //AliHLTTPCCAPerformance::Instance().LinkPerformance( Param().ISlice() );
 522 #endif
 523
 524
 525 #ifdef DRAW1
 526   if ( NHitsTotal() > 0 ) {
 527     AliHLTTPCCADisplay::Instance().DrawSliceLinks( -1, -1, 1 );
 528     AliHLTTPCCADisplay::Instance().Ask();
 529   }
 530 #endif
 531
 532   RunNeighboursCleaner();
 533
 534   StandalonePerfTime(3);
 535
 536   if (fGPUDebugLevel >= 6) DumpLinks(*fGPUDebugOut);
 537
 538   RunStartHitsFinder();
 539
 540   StandalonePerfTime(4);
 541   StandalonePerfTime(5);
 542
 543   if (fGPUDebugLevel >= 6) DumpStartHits(*fGPUDebugOut);
 544
 545   fData.ClearHitWeights();
 546
 547   SetPointersTracklets( fCommonMem->fNTracklets * 2 ); // to calculate the size
 548   fTrackletMemory = reinterpret_cast<char*> ( new uint4 [ fTrackletMemorySize/sizeof( uint4 ) + 100] );
 549   SetPointersTracklets( fCommonMem->fNTracklets * 2 ); // set pointers for hits
 550
 551   SetPointersTracks( fCommonMem->fNTracklets * 2, NHitsTotal() ); // to calculate the size
 552   fTrackMemory = reinterpret_cast<char*> ( new uint4 [ fTrackMemorySize/sizeof( uint4 ) + 100] );
 553   SetPointersTracks( fCommonMem->fNTracklets * 2, NHitsTotal() ); // set pointers for hits
 554
 555   StandalonePerfTime(6);
 556   StandalonePerfTime(7);
 557
 558   RunTrackletConstructor();
 559
 560   StandalonePerfTime(8);
 561
 562   if (fGPUDebugLevel >= 6) DumpTrackletHits(*fGPUDebugOut);
 563   if (fGPUDebugLevel >= 6) DumpHitWeights(*fGPUDebugOut);
 564
 565   //std::cout<<"Slice "<<Param().ISlice()<<": NHits="<<NHitsTotal()<<", NTracklets="<<*NTracklets()<<std::endl;
 566
 567   RunTrackletSelector();
 568
 569   StandalonePerfTime(9);
 570
 571   //std::cout<<"Slice "<<Param().ISlice()<<": N start hits/tracklets/tracks = "<<nStartHits<<" "<<nStartHits<<" "<<*fNTracks<<std::endl;
 572
 573   if (fGPUDebugLevel >= 6) DumpTrackHits(*fGPUDebugOut);
 574
 575   //std::cout<<"Memory used for slice "<<fParam.ISlice()<<" : "<<fCommonMemorySize/1024./1024.<<" + "<<fHitMemorySize/1024./1024.<<" + "<<fTrackMemorySize/1024./1024.<<" = "<<( fCommonMemorySize+fHitMemorySize+fTrackMemorySize )/1024./1024.<<" Mb "<<std::endl;
 576
 577   WriteOutput();
 578
 579   StandalonePerfTime(10);
 580
 581 #endif
 582
 583 #ifdef DRAW1
 584   {
 585     AliHLTTPCCADisplay &disp = AliHLTTPCCADisplay::Instance();
 586     AliHLTTPCCATracker &slice = *this;
 587     std::cout << "N out tracks = " << slice.NOutTracks() << std::endl;
 588     AliHLTTPCCADisplay::Instance().SetSliceView();
 589     AliHLTTPCCADisplay::Instance().SetCurrentSlice( this );
 590     AliHLTTPCCADisplay::Instance().DrawSlice( this, 1 );
 591     disp.DrawSliceHits( kRed, .5 );
 592     disp.Ask();
 593     for ( int itr = 0; itr < slice.NOutTracks(); itr++ ) {
 594       std::cout << "track N " << itr << ", nhits=" << slice.OutTracks()[itr].NHits() << std::endl;
 595       disp.DrawSliceOutTrack( itr, kBlue );
 596       //disp.Ask();
 597       //int id = slice.OutTracks()[itr].OrigTrackID();
 598       //AliHLTTPCCATrack &tr = Tracks()[id];
 599       //for( int ih=0; ih<tr.NHits(); ih++ ){
 600       //int ic = (fTrackHits[tr.FirstHitID()+ih]);
 601       //std::cout<<ih<<" "<<ID2IRow(ic)<<" "<<ID2IHit(ic)<<std::endl;
 602       //}
 603       //disp.DrawSliceTrack( id, kBlue );
 604       //disp.Ask();
 605     }
 606     disp.Ask();
 607   }
 608 #endif
 609
 610   timer0.Stop();
 611   fTimers[0] = timer0.CpuTime() / 100.;
 612
 613 }
 614
 615 GPUh() void AliHLTTPCCATracker::WriteOutput()
 616 {
 617   // write output
 618
 619   TStopwatch timer;
 620
 621   //cout<<"output: nTracks = "<<*fNTracks<<", nHitsTotal="<<NHitsTotal()<<std::endl;
 622
 623   if (fOutputControl == NULL) fOutputControl = new AliHLTTPCCASliceOutput::outputControlStruct;
 624   AliHLTTPCCASliceOutput::Allocate(*fOutput, fCommonMem->fNTracks, fCommonMem->fNTrackHits, fOutputControl);
 625   AliHLTTPCCASliceOutput* useOutput = *fOutput;
 626   if (useOutput == NULL) return;
 627
 628   if (fOutputControl->fDefaultOutput)
 629   {
 630           useOutput->SetNTracks( fCommonMem->fNTracks );
 631           useOutput->SetNTrackClusters( fCommonMem->fNTrackHits );
 632
 633           int nStoredHits = 0;
 634
 635           for ( int iTr = 0; iTr < fCommonMem->fNTracks; iTr++ ) {
 636                 AliHLTTPCCATrack &iTrack = fTracks[iTr];
 637
 638                 AliHLTTPCCASliceTrack out;
 639                 out.SetFirstClusterRef( nStoredHits );
 640                 out.SetNClusters( iTrack.NHits() );
 641                 out.SetParam( iTrack.Param() );
 642
 643                 useOutput->SetTrack( iTr, out );
 644
 645                 int iID = iTrack.FirstHitID();
 646                 for ( int ith = 0; ith < iTrack.NHits(); ith++ ) {
 647                   const AliHLTTPCCAHitId &ic = fTrackHits[iID + ith];
 648                   int iRow = ic.RowIndex();
 649                   int ih = ic.HitIndex();
 650
 651                   const AliHLTTPCCARow &row = fData.Row( iRow );
 652
 653                   //float y0 = row.Grid().YMin();
 654                   //float z0 = row.Grid().ZMin();
 655                   //float stepY = row.HstepY();
 656                   //float stepZ = row.HstepZ();
 657                   //float x = row.X();
 658
 659                   //const uint4 *tmpint4 = RowData() + row.FullOffset();
 660                   //const ushort2 *hits = reinterpret_cast<const ushort2*>(tmpint4);
 661                   //ushort2 hh = hits[ih];
 662                   //float y = y0 + hh.x*stepY;
 663                   //float z = z0 + hh.y*stepZ;
 664
 665                   int clusterIndex = fData.ClusterDataIndex( row, ih );
 666                   int clusterRowIndex = clusterIndex - fClusterData->RowOffset( iRow );
 667
 668                   if ( clusterIndex < 0 || clusterIndex >= fClusterData->NumberOfClusters() ) {
 669                         //std::cout << inpIDtot << ", " << fClusterData->NumberOfClusters()
 670                         //<< "; " << inpID << ", " << fClusterData->NumberOfClusters( iRow ) << std::endl;
 671                         //abort();
 672                         continue;
 673                   }
 674                   if ( clusterRowIndex < 0 || clusterRowIndex >= fClusterData->NumberOfClusters( iRow ) ) {
 675                         //std::cout << inpIDtot << ", " << fClusterData->NumberOfClusters()
 676                         //<< "; " << inpID << ", " << fClusterData->NumberOfClusters( iRow ) << std::endl;
 677                         //abort();
 678                         continue;
 679                   }
 680
 681                   float origX = fClusterData->X( clusterIndex );
 682                   float origY = fClusterData->Y( clusterIndex );
 683                   float origZ = fClusterData->Z( clusterIndex );
 684
 685
 686                   int id = fClusterData->Id( clusterIndex );
 687
 688                   unsigned short hPackedYZ = 0;
 689                   UChar_t hPackedAmp = 0;
 690                   float2 hUnpackedYZ;
 691                   hUnpackedYZ.x = origY;
 692                   hUnpackedYZ.y = origZ;
 693                   float hUnpackedX = origX;
 694
 695                   useOutput->SetClusterId( nStoredHits, id  );
 696                   useOutput->SetClusterRow( nStoredHits, ( unsigned char ) iRow  );
 697                   useOutput->SetClusterPackedYZ( nStoredHits, hPackedYZ );
 698                   useOutput->SetClusterPackedAmp( nStoredHits, hPackedAmp );
 699                   useOutput->SetClusterUnpackedYZ( nStoredHits, hUnpackedYZ );
 700                   useOutput->SetClusterUnpackedX( nStoredHits, hUnpackedX );
 701                   nStoredHits++;
 702                 }
 703           }
 704   }
 705
 706
 707   // old stuff
 708   if (fOutputControl->fObsoleteOutput)
 709   {
 710           useOutput->SetNOutTrackHits(0);
 711           useOutput->SetNOutTracks(0);
 712
 713
 714           for ( int iTr = 0; iTr < fCommonMem->fNTracks; iTr++ ) {
 715
 716                 const AliHLTTPCCATrack &iTrack = fTracks[iTr];
 717
 718                 //std::cout<<"iTr = "<<iTr<<", nHits="<<iTrack.NHits()<<std::endl;
 719
 720                 //if( !iTrack.Alive() ) continue;
 721                 if ( iTrack.NHits() < 3 ) continue;
 722                 AliHLTTPCCAOutTrack &out = useOutput->OutTracks()[useOutput->NOutTracks()];
 723                 out.SetFirstHitRef( useOutput->NOutTrackHits() );
 724                 out.SetNHits( 0 );
 725                 out.SetOrigTrackID( iTr );
 726                 out.SetStartPoint( iTrack.Param() );
 727                 out.SetEndPoint( iTrack.Param() );
 728
 729                 int iID = iTrack.FirstHitID();
 730                 int nOutTrackHitsOld = useOutput->NOutTrackHits();
 731
 732                 for ( int ith = 0; ith < iTrack.NHits(); ith++ ) {
 733                   const AliHLTTPCCAHitId &ic = fTrackHits[iID + ith];
 734                   const AliHLTTPCCARow &row = Row( ic );
 735                   int ih = ic.HitIndex();
 736                   useOutput->SetOutTrackHit(useOutput->NOutTrackHits(), HitInputID( row, ih ));
 737                   useOutput->SetNOutTrackHits(useOutput->NOutTrackHits() + 1 );
 738                   //std::cout<<"write i,row,hit,id="<<ith<<", "<<ID2IRow(ic)<<", "<<ih<<", "<<HitInputID( row, ih )<<std::endl;
 739                   if ( useOutput->NOutTrackHits() >= 10*NHitsTotal() ) {
 740                         std::cout << "fNOutTrackHits>NHitsTotal()" << std::endl;
 741                         //exit(0);
 742                         return;//SG!!!
 743                   }
 744                   out.SetNHits( out.NHits() + 1 );
 745                 }
 746                 if ( out.NHits() >= 2 ) {
 747                   useOutput->SetNOutTracks(useOutput->NOutTracks() + 1);
 748                 } else {
 749                   useOutput->SetNOutTrackHits(nOutTrackHitsOld);
 750                 }
 751           }
 752   }
 753
 754   timer.Stop();
 755   fTimers[5] += timer.CpuTime();
 756 }
 757
 758 #endif
 759
 760 GPUh() void AliHLTTPCCATracker::FitTrackFull( const AliHLTTPCCATrack &/**/, float * /**/ ) const
 761 {
 762   // fit track with material
 763 #ifdef XXX
 764   //* Fit the track
 765   FitTrack( iTrack, tt0 );
 766   if ( iTrack.NHits() <= 3 ) return;
 767
 768   AliHLTTPCCATrackParam &t = iTrack.Param();
 769   AliHLTTPCCATrackParam t0 = t;
 770
 771   t.Chi2() = 0;
 772   t.NDF() = -5;
 773   bool first = 1;
 774
 775   int iID = iTrack.FirstHitID();
 776   for ( int ih = 0; ih < iTrack.NHits(); ih++, iID++ ) {
 777     const AliHLTTPCCAHitId &ic = fTrackHits[iID];
 778     int iRow = ic.rowIndex();
 779     const AliHLTTPCCARow &row = fData.Row( iRow );
 780     if ( !t0.TransportToX( row.X() ) ) continue;
 781     float dy, dz;
 782     const AliHLTTPCCAHit &h = ic.hitIndex();
 783
 784     // check for wrong hits
 785     if ( 0 ) {
 786       dy = t0.GetY() - h.Y();
 787       dz = t0.GetZ() - h.Z();
 788
 789       //if( dy*dy > 3.5*3.5*(/*t0.GetErr2Y() + */h.ErrY()*h.ErrY() ) ) continue;//SG!!!
 790       //if( dz*dz > 3.5*3.5*(/*t0.GetErr2Z() + */h.ErrZ()*h.ErrZ() ) ) continue;
 791     }
 792
 793     if ( !t.TransportToX( row.X() ) ) continue;
 794
 795     //* Update the track
 796
 797     if ( first ) {
 798       t.Cov()[ 0] = .5 * .5;
 799       t.Cov()[ 1] = 0;
 800       t.Cov()[ 2] = .5 * .5;
 801       t.Cov()[ 3] = 0;
 802       t.Cov()[ 4] = 0;
 803       t.Cov()[ 5] = .2 * .2;
 804       t.Cov()[ 6] = 0;
 805       t.Cov()[ 7] = 0;
 806       t.Cov()[ 8] = 0;
 807       t.Cov()[ 9] = .2 * .2;
 808       t.Cov()[10] = 0;
 809       t.Cov()[11] = 0;
 810       t.Cov()[12] = 0;
 811       t.Cov()[13] = 0;
 812       t.Cov()[14] = .2 * .2;
 813       t.Chi2() = 0;
 814       t.NDF() = -5;
 815     }
 816     float err2Y, err2Z;
 817     GetErrors2( iRow, t, err2Y, err2Z );
 818
 819     if ( !t.Filter2( h.Y(), h.Z(), err2Y, err2Z ) ) continue;
 820
 821     first = 0;
 822   }
 823   /*
 824   float cosPhi = iTrack.Param().GetCosPhi();
 825   p0.Param().TransportToX(ID2Row( iTrack.PointID()[0] ).X());
 826   p2.Param().TransportToX(ID2Row( iTrack.PointID()[1] ).X());
 827   if( p0.Param().GetCosPhi()*cosPhi<0 ){ // change direction
 828   float *par = p0.Param().Par();
 829   float *cov = p0.Param().Cov();
 830   par[2] = -par[2]; // sin phi
 831   par[3] = -par[3]; // DzDs
 832     par[4] = -par[4]; // kappa
 833     cov[3] = -cov[3];
 834     cov[4] = -cov[4];
 835     cov[6] = -cov[6];
 836     cov[7] = -cov[7];
 837     cov[10] = -cov[10];
 838     cov[11] = -cov[11];
 839     p0.Param().CosPhi() = -p0.Param().GetCosPhi();
 840   }
 841   */
 842 #endif
 843 }
 844
 845 GPUh() void AliHLTTPCCATracker::FitTrack( const AliHLTTPCCATrack &/*track*/, float * /*t0[]*/ ) const
 846 {
 847   //* Fit the track
 848 #ifdef XXX
 849   AliHLTTPCCAEndPoint &p2 = ID2Point( track.PointID()[1] );
 850   const AliHLTTPCCAHit &c0 = ID2Hit( fTrackHits[p0.TrackHitID()].HitID() );
 851   const AliHLTTPCCAHit &c1 = ID2Hit( fTrackHits[track.HitID()[1]].HitID() );
 852   const AliHLTTPCCAHit &c2 = ID2Hit( fTrackHits[p2.TrackHitID()].HitID() );
 853   const AliHLTTPCCARow &row0 = ID2Row( fTrackHits[p0.TrackHitID()].HitID() );
 854   const AliHLTTPCCARow &row1 = ID2Row( fTrackHits[track.HitID()[1]].HitID() );
 855   const AliHLTTPCCARow &row2 = ID2Row( fTrackHits[p2.TrackHitID()].HitID() );
 856   float sp0[5] = {row0.X(), c0.Y(), c0.Z(), c0.ErrY(), c0.ErrZ() };
 857   float sp1[5] = {row1.X(), c1.Y(), c1.Z(), c1.ErrY(), c1.ErrZ() };
 858   float sp2[5] = {row2.X(), c2.Y(), c2.Z(), c2.ErrY(), c2.ErrZ() };
 859   //std::cout<<"Fit track, points ="<<sp0[0]<<" "<<sp0[1]<<" / "<<sp1[0]<<" "<<sp1[1]<<" / "<<sp2[0]<<" "<<sp2[1]<<std::endl;
 860   if ( track.NHits() >= 3 ) {
 861     p0.Param().ConstructXYZ3( sp0, sp1, sp2, p0.Param().CosPhi(), t0 );
 862     p2.Param().ConstructXYZ3( sp2, sp1, sp0, p2.Param().CosPhi(), t0 );
 863     //p2.Param() = p0.Param();
 864     //p2.Param().TransportToX(row2.X());
 865     //p2.Param().Par()[1] = -p2.Param().Par()[1];
 866     //p2.Param().Par()[4] = -p2.Param().Par()[4];
 867   } else {
 868     p0.Param().X() = row0.X();
 869     p0.Param().Y() = c0.Y();
 870     p0.Param().Z() = c0.Z();
 871     p0.Param().Err2Y() = c0.ErrY() * c0.ErrY();
 872     p0.Param().Err2Z() = c0.ErrZ() * c0.ErrZ();
 873     p2.Param().X() = row2.X();
 874     p2.Param().Y() = c2.Y();
 875     p2.Param().Z() = c2.Z();
 876     p2.Param().Err2Y() = c2.ErrY() * c2.ErrY();
 877     p2.Param().Err2Z() = c2.ErrZ() * c2.ErrZ();
 878   }
 879 #endif
 880 }
 881
 882
 883 GPUd() void AliHLTTPCCATracker::GetErrors2( int iRow, float z, float sinPhi, float cosPhi, float DzDs, float &Err2Y, float &Err2Z ) const
 884 {
 885   //
 886   // Use calibrated cluster error from OCDB
 887   //
 888
 889   fParam.GetClusterErrors2( iRow, z, sinPhi, cosPhi, DzDs, Err2Y, Err2Z );
 890   Err2Y*=fParam.ClusterError2CorrectionY();
 891   Err2Z*=fParam.ClusterError2CorrectionZ();
 892 }
 893
 894 GPUd() void AliHLTTPCCATracker::GetErrors2( int iRow, const AliHLTTPCCATrackParam &t, float &Err2Y, float &Err2Z ) const
 895 {
 896   //
 897   // Use calibrated cluster error from OCDB
 898   //
 899
 900   fParam.GetClusterErrors2( iRow, t.GetZ(), t.SinPhi(), t.GetCosPhi(), t.DzDs(), Err2Y, Err2Z );
 901 }
 902
 903
 904 #if !defined(HLTCA_GPUCODE)
 905
 906 GPUh() void AliHLTTPCCATracker::WriteEvent( std::ostream &out )
 907 {
 908   // write event to the file
 909   for ( int iRow = 0; iRow < fParam.NRows(); iRow++ ) {
 910     out << fData.Row( iRow ).HitNumberOffset() << " " << fData.Row( iRow ).NHits() << std::endl;
 911   }
 912   out << NHitsTotal() << std::endl;
 913
 914   AliHLTResizableArray<float> y( NHitsTotal() ), z( NHitsTotal() );
 915
 916   for ( int iRow = 0; iRow < fParam.NRows(); iRow++ ) {
 917     const AliHLTTPCCARow &row = Row( iRow );
 918     float y0 = row.Grid().YMin();
 919     float z0 = row.Grid().ZMin();
 920     float stepY = row.HstepY();
 921     float stepZ = row.HstepZ();
 922     for ( int ih = 0; ih < fData.Row( iRow ).NHits(); ih++ ) {
 923       int id = HitInputID( row, ih );
 924       y[id] = y0 + HitDataY( row, ih ) * stepY;
 925       z[id] = z0 + HitDataZ( row, ih ) * stepZ;
 926     }
 927   }
 928   for ( int ih = 0; ih < NHitsTotal(); ih++ ) {
 929     out << y[ih] << " " << z[ih] << std::endl;
 930   }
 931 }
 932
 933 GPUh() void AliHLTTPCCATracker::WriteTracks( std::ostream &out )
 934 {
 935   //* Write tracks to file
 936   AliHLTTPCCASliceOutput* useOutput = *fOutput;
 937
 938   out << fTimers[0] << std::endl;
 939   out << useOutput->NOutTrackHits() << std::endl;
 940   for ( int ih = 0; ih < useOutput->NOutTrackHits(); ih++ ) {
 941     out << useOutput->OutTrackHit(ih) << " ";
 942   }
 943   out << std::endl;
 944
 945   out << useOutput->NOutTracks() << std::endl;
 946
 947   for ( int itr = 0; itr < useOutput->NOutTracks(); itr++ ) {
 948     const AliHLTTPCCAOutTrack &t = useOutput->OutTrack(itr);
 949     AliHLTTPCCATrackParam p1 = t.StartPoint();
 950     AliHLTTPCCATrackParam p2 = t.EndPoint();
 951     out << t.NHits() << " ";
 952     out << t.FirstHitRef() << " ";
 953     out << t.OrigTrackID() << " ";
 954     out << std::endl;
 955     out << p1.X() << " ";
 956     out << p1.SignCosPhi() << " ";
 957     out << p1.Chi2() << " ";
 958     out << p1.NDF() << std::endl;
 959     for ( int i = 0; i < 5; i++ ) out << p1.Par()[i] << " ";
 960     out << std::endl;
 961     for ( int i = 0; i < 15; i++ ) out << p1.Cov()[i] << " ";
 962     out << std::endl;
 963     out << p2.X() << " ";
 964     out << p2.SignCosPhi() << " ";
 965     out << p2.Chi2() << " ";
 966     out << p2.NDF() << std::endl;
 967     for ( int i = 0; i < 5; i++ ) out << p2.Par()[i] << " ";
 968     out << std::endl;
 969     for ( int i = 0; i < 15; i++ ) out << p2.Cov()[i] << " ";
 970     out << std::endl;
 971   }
 972 }
 973
 974 GPUh() void AliHLTTPCCATracker::ReadTracks( std::istream &in )
 975 {
 976   //* Read tracks  from file
 977   AliHLTTPCCASliceOutput::Allocate(*fOutput, 4096, 16384, fOutputControl);//Just some max values
 978   AliHLTTPCCASliceOutput* useOutput = *fOutput;
 979
 980   int tmpval;
 981   in >> fTimers[0];
 982   in >> tmpval;
 983   useOutput->SetNOutTrackHits(tmpval);
 984
 985   for ( int ih = 0; ih < useOutput->NOutTrackHits(); ih++ ) {
 986     in >> tmpval;
 987         useOutput->SetOutTrackHit(ih, tmpval);
 988   }
 989   in >> tmpval;
 990   useOutput->SetNOutTracks(tmpval);
 991
 992   for ( int itr = 0; itr < useOutput->NOutTracks(); itr++ ) {
 993     AliHLTTPCCAOutTrack &t = useOutput->OutTracks()[itr];
 994     AliHLTTPCCATrackParam p1, p2;
 995     int i;
 996     float f;
 997     in >> i; t.SetNHits( i );
 998     in >> i; t.SetFirstHitRef( i );
 999     in >> i; t.SetOrigTrackID( i );
1000     in >> f; p1.SetX( f );
1001     in >> f; p1.SetSignCosPhi( f );
1002     in >> f; p1.SetChi2( f );
1003     in >> i; p1.SetNDF( i );
1004     for ( int j = 0; j < 5; j++ ) { in >> f; p1.SetPar( j, f ); }
1005     for ( int j = 0; j < 15; j++ ) { in >> f; p1.SetCov( j, f ); }
1006     in >> f; p2.SetX( f );
1007     in >> f; p2.SetSignCosPhi( f );
1008     in >> f; p2.SetChi2( f );
1009     in >> i; p2.SetNDF( i );
1010     for ( int j = 0; j < 5; j++ ) { in >> f; p2.SetPar( j, f ); }
1011     for ( int j = 0; j < 15; j++ ) { in >> f; p2.SetCov( j, f ); }
1012     t.SetStartPoint( p1 );
1013     t.SetEndPoint( p2 );
1014   }
1015 }
1016 #endif