HLT/TPCLib/merger-ca/AliHLTTPCGMMerger.cxx

   1 // $Id: AliHLTTPCGMMerger.cxx 30732 2009-01-22 23:02:02Z sgorbuno $
   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 <stdio.h>
  21 #include "AliHLTTPCCASliceOutTrack.h"
  22 #include "AliHLTTPCCATracker.h"
  23 #include "AliHLTTPCCATrackParam.h"
  24 #include "AliHLTTPCGMCluster.h"
  25
  26 #include "AliHLTTPCGMMerger.h"
  27
  28 #include "AliHLTTPCCAMath.h"
  29 #include "TStopwatch.h"
  30
  31 #include "AliHLTTPCCATrackParam.h"
  32 #include "AliHLTTPCCASliceOutput.h"
  33 #include "AliHLTTPCGMMergedTrack.h"
  34 #include "AliHLTTPCCADataCompressor.h"
  35 #include "AliHLTTPCCAParam.h"
  36 #include "AliHLTTPCCATrackLinearisation.h"
  37 #include "AliHLTTPCCADataCompressor.h"
  38
  39 #include "AliHLTTPCGMTrackParam.h"
  40 #include "AliHLTTPCGMTrackLinearisation.h"
  41 #include "AliHLTTPCGMSliceTrack.h"
  42 #include "AliHLTTPCGMBorderTrack.h"
  43 #include <cmath>
  44
  45 #include <algorithm>
  46
  47 #include "AliHLTTPCCAGPUConfig.h"
  48 #include "MemoryAssignmentHelpers.h"
  49
  50 #define GLOBAL_TRACKS_SPECIAL_TREATMENT
  51
  52 AliHLTTPCGMMerger::AliHLTTPCGMMerger()
  53   :
  54   fSliceParam(),
  55   fNOutputTracks( 0 ),
  56   fNOutputTrackClusters( 0 ),
  57   fOutputTracks( 0 ),
  58   fOutputClusterIds(0),
  59   fSliceTrackInfos( 0 ),
  60   fMaxSliceTracks(0),
  61   fClusterX(0),
  62   fClusterY(0),
  63   fClusterZ(0),
  64   fClusterRowType(0),
  65   fClusterAngle(0),
  66   fBorderMemory(0),
  67   fBorderRangeMemory(0),
  68   fGPUTracker(NULL),
  69   fDebugLevel(0),
  70   fNClusters(0)
  71 {
  72   //* constructor
  73
  74   for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
  75     fNextSliceInd[iSlice] = iSlice + 1;
  76     fPrevSliceInd[iSlice] = iSlice - 1;
  77   }
  78   int mid = fgkNSlices / 2 - 1 ;
  79   int last = fgkNSlices - 1 ;
  80   fNextSliceInd[ mid ] = 0;
  81   fPrevSliceInd[ 0 ] = mid;  fNextSliceInd[ last ] = fgkNSlices / 2;
  82   fPrevSliceInd[ fgkNSlices/2 ] = last;
  83   {
  84     const double kCLight = 0.000299792458;
  85     double constBz = fSliceParam.BzkG() * kCLight;
  86
  87     fPolinomialFieldBz[0] = constBz * (  0.999286   );
  88     fPolinomialFieldBz[1] = constBz * ( -4.54386e-7 );
  89     fPolinomialFieldBz[2] = constBz * (  2.32950e-5 );
  90     fPolinomialFieldBz[3] = constBz * ( -2.99912e-7 );
  91     fPolinomialFieldBz[4] = constBz * ( -2.03442e-8 );
  92     fPolinomialFieldBz[5] = constBz * (  9.71402e-8 );
  93   }
  94
  95   Clear();
  96 }
  97
  98
  99 AliHLTTPCGMMerger::AliHLTTPCGMMerger(const AliHLTTPCGMMerger&)
 100   :
 101   fSliceParam(),
 102   fNOutputTracks( 0 ),
 103   fNOutputTrackClusters( 0 ),
 104   fOutputTracks( 0 ),
 105   fOutputClusterIds(0),
 106   fSliceTrackInfos( 0 ),
 107   fMaxSliceTracks(0),
 108   fClusterX(0),
 109   fClusterY(0),
 110   fClusterZ(0),
 111   fClusterRowType(0),
 112   fClusterAngle(0),
 113   fBorderMemory(0),
 114   fBorderRangeMemory(0),
 115   fGPUTracker(NULL),
 116   fDebugLevel(0),
 117   fNClusters(0)
 118 {
 119   //* dummy
 120   for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
 121     fNextSliceInd[iSlice] = 0;
 122     fPrevSliceInd[iSlice] = 0;
 123   }
 124   {
 125     const double kCLight = 0.000299792458;
 126     double constBz = fSliceParam.BzkG() * kCLight;
 127
 128     fPolinomialFieldBz[0] = constBz * (  0.999286   );
 129     fPolinomialFieldBz[1] = constBz * ( -4.54386e-7 );
 130     fPolinomialFieldBz[2] = constBz * (  2.32950e-5 );
 131     fPolinomialFieldBz[3] = constBz * ( -2.99912e-7 );
 132     fPolinomialFieldBz[4] = constBz * ( -2.03442e-8 );
 133     fPolinomialFieldBz[5] = constBz * (  9.71402e-8 );
 134   }
 135   Clear();
 136 }
 137
 138 const AliHLTTPCGMMerger &AliHLTTPCGMMerger::operator=(const AliHLTTPCGMMerger&) const
 139 {
 140   //* dummy
 141   return *this;
 142 }
 143
 144
 145 AliHLTTPCGMMerger::~AliHLTTPCGMMerger()
 146 {
 147   //* destructor
 148   ClearMemory();
 149 }
 150
 151 void AliHLTTPCGMMerger::Clear()
 152 {
 153   for ( int i = 0; i < fgkNSlices; ++i ) {
 154     fkSlices[i] = 0;
 155     fSliceNTrackInfos[ i ] = 0;
 156     fSliceTrackInfoStart[ i ] = 0;
 157   }
 158   ClearMemory();
 159 }
 160
 161
 162 void AliHLTTPCGMMerger::ClearMemory()
 163 {
 164   if (fOutputClusterIds) delete[] fOutputClusterIds;
 165   if (fSliceTrackInfos) delete[] fSliceTrackInfos;
 166   if (!fGPUTracker)
 167   {
 168           if (fOutputTracks) delete[] fOutputTracks;
 169           if (fClusterX) delete[] fClusterX;
 170           if (fClusterY) delete[] fClusterY;
 171           if (fClusterZ) delete[] fClusterZ;
 172           if (fClusterRowType) delete[] fClusterRowType;
 173           if (fClusterAngle) delete[] fClusterAngle;
 174   }
 175   if (fBorderMemory) delete[] fBorderMemory;
 176   if (fBorderRangeMemory) delete[] fBorderRangeMemory;
 177
 178   fNOutputTracks = 0;
 179   fOutputTracks = 0;
 180   fOutputClusterIds = 0;
 181   fSliceTrackInfos = 0;
 182   fMaxSliceTracks = 0;
 183   fClusterX = 0;
 184   fClusterY = 0;
 185   fClusterZ = 0;
 186   fClusterRowType = 0;
 187   fClusterAngle = 0;
 188   fBorderMemory = 0;
 189   fBorderRangeMemory = 0;
 190 }
 191
 192
 193 void AliHLTTPCGMMerger::SetSliceData( int index, const AliHLTTPCCASliceOutput *sliceData )
 194 {
 195   fkSlices[index] = sliceData;
 196 }
 197
 198
 199 bool AliHLTTPCGMMerger::Reconstruct()
 200 {
 201   //* main merging routine
 202
 203   {
 204     const double kCLight = 0.000299792458;
 205     double constBz = fSliceParam.BzkG() * kCLight;
 206
 207     fPolinomialFieldBz[0] = constBz * (  0.999286   );
 208     fPolinomialFieldBz[1] = constBz * ( -4.54386e-7 );
 209     fPolinomialFieldBz[2] = constBz * (  2.32950e-5 );
 210     fPolinomialFieldBz[3] = constBz * ( -2.99912e-7 );
 211     fPolinomialFieldBz[4] = constBz * ( -2.03442e-8 );
 212     fPolinomialFieldBz[5] = constBz * (  9.71402e-8 );
 213   }
 214
 215   int nIter = 1;
 216   TStopwatch timer;
 217 #ifdef HLTCA_STANDALONE
 218   unsigned long long int a, b, c, d, e, f, g;
 219   AliHLTTPCCATracker::StandaloneQueryFreq(&g);
 220 #endif
 221   //cout<<"Merger..."<<endl;
 222   for( int iter=0; iter<nIter; iter++ ){
 223     if( !AllocateMemory() ) return 0;
 224 #ifdef HLTCA_STANDALONE
 225         AliHLTTPCCATracker::StandaloneQueryTime(&a);
 226 #endif
 227     UnpackSlices();
 228 #ifdef HLTCA_STANDALONE
 229         AliHLTTPCCATracker::StandaloneQueryTime(&b);
 230 #endif
 231     MergeWithingSlices();
 232 #ifdef HLTCA_STANDALONE
 233         AliHLTTPCCATracker::StandaloneQueryTime(&c);
 234 #endif
 235     MergeSlices();
 236 #ifdef HLTCA_STANDALONE
 237         AliHLTTPCCATracker::StandaloneQueryTime(&d);
 238 #endif
 239     CollectMergedTracks();
 240 #ifdef HLTCA_STANDALONE
 241         AliHLTTPCCATracker::StandaloneQueryTime(&e);
 242 #endif
 243     Refit();
 244 #ifdef HLTCA_STANDALONE
 245         AliHLTTPCCATracker::StandaloneQueryTime(&f);
 246         if (fDebugLevel > 0)
 247         {
 248                 printf("Merge Time:\tUnpack Slices:\t%lld us\n", (b - a) * 1000000 / g);
 249                 printf("\t\tMerge Within:\t%lld us\n", (c - b) * 1000000 / g);
 250                 printf("\t\tMerge Slices:\t%lld us\n", (d - c) * 1000000 / g);
 251                 printf("\t\tCollect:\t%lld us\n", (e - d) * 1000000 / g);
 252                 printf("\t\tRefit:\t\t%lld us\n", (f - e) * 1000000 / g);
 253         }
 254         int newTracks = 0;
 255         for (int i = 0;i < fNOutputTracks;i++) if (fOutputTracks[i].OK()) newTracks++;
 256         printf("Output Tracks: %d\n", newTracks);
 257 #endif
 258   }
 259   timer.Stop();
 260   //cout<<"\nMerger time = "<<timer.CpuTime()*1.e3/nIter<<" ms\n"<<endl;
 261
 262   return 1;
 263 }
 264
 265
 266
 267 bool AliHLTTPCGMMerger::AllocateMemory()
 268 {
 269   //* memory allocation
 270
 271   ClearMemory();
 272
 273   int nTracks = 0;
 274   fNClusters = 0;
 275   fMaxSliceTracks  = 0;
 276
 277   for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
 278     if ( !fkSlices[iSlice] ) continue;
 279     nTracks += fkSlices[iSlice]->NTracks();
 280     fNClusters += fkSlices[iSlice]->NTrackClusters();
 281     if( fMaxSliceTracks < fkSlices[iSlice]->NTracks() ) fMaxSliceTracks = fkSlices[iSlice]->NTracks();
 282   }
 283
 284   //cout<<"\nMerger: input "<<nTracks<<" tracks, "<<nClusters<<" clusters"<<endl;
 285
 286   fOutputClusterIds = new UInt_t[fNClusters];
 287   fSliceTrackInfos = new AliHLTTPCGMSliceTrack[nTracks];
 288   if (fGPUTracker)
 289   {
 290         char* basemem = fGPUTracker->MergerBaseMemory();
 291         AssignMemory(fClusterX, basemem, fNClusters);
 292         AssignMemory(fClusterY, basemem, fNClusters);
 293         AssignMemory(fClusterZ, basemem, fNClusters);
 294         AssignMemory(fClusterAngle, basemem, fNClusters);
 295         AssignMemory(fClusterRowType, basemem, fNClusters);
 296         AssignMemory(fOutputTracks, basemem, nTracks);
 297   }
 298   else
 299   {
 300           fOutputTracks = new AliHLTTPCGMMergedTrack[nTracks];
 301           fClusterX = new float[fNClusters];
 302           fClusterY = new float[fNClusters];
 303           fClusterZ = new float[fNClusters];
 304           fClusterRowType = new UInt_t[fNClusters];
 305           fClusterAngle = new float[fNClusters];
 306   }
 307   fBorderMemory = new AliHLTTPCGMBorderTrack[fMaxSliceTracks*2];
 308   fBorderRangeMemory = new AliHLTTPCGMBorderTrack::Range[fMaxSliceTracks*2];
 309
 310   return ( ( fOutputTracks!=NULL )
 311            && ( fOutputClusterIds!=NULL )
 312            && ( fSliceTrackInfos!=NULL )
 313            && ( fClusterX!=NULL )
 314            && ( fClusterY!=NULL )
 315            && ( fClusterZ!=NULL )
 316            && ( fClusterRowType!=NULL )
 317            && ( fClusterAngle!=NULL )
 318            && ( fBorderMemory!=NULL )
 319            && ( fBorderRangeMemory!=NULL )
 320            );
 321 }
 322
 323
 324
 325 void AliHLTTPCGMMerger::UnpackSlices()
 326 {
 327   //* unpack the cluster information from the slice tracks and initialize track info array
 328
 329   int nTracksCurrent = 0;
 330
 331   const AliHLTTPCCASliceOutTrack* firstGlobalTracks[fgkNSlices];
 332
 333   for( int i=0; i<fgkNSlices; i++) firstGlobalTracks[i] = 0;
 334
 335 #ifdef GLOBAL_TRACKS_SPECIAL_TREATMENT
 336   const int kMaxTrackIdInSlice = AliHLTTPCCASliceOutTrack::MaxTrackId();
 337   int* TrackIds = (int*) malloc(kMaxTrackIdInSlice * fgkNSlices * sizeof(int));
 338   for (int i = 0;i < kMaxTrackIdInSlice * fgkNSlices;i++) TrackIds[i] = -1;
 339 #endif
 340
 341   for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
 342
 343     fSliceTrackInfoStart[ iSlice ] = nTracksCurrent;
 344
 345     fSliceNTrackInfos[ iSlice ] = 0;
 346
 347     if ( !fkSlices[iSlice] ) continue;
 348
 349     float alpha = fSliceParam.Alpha( iSlice );
 350
 351     const AliHLTTPCCASliceOutput &slice = *( fkSlices[iSlice] );
 352     const AliHLTTPCCASliceOutTrack *sliceTr = slice.GetFirstTrack();
 353
 354 #ifdef GLOBAL_TRACKS_SPECIAL_TREATMENT
 355     for ( int itr = 0; itr < slice.NLocalTracks(); itr++, sliceTr = sliceTr->GetNextTrack() ) {
 356 #else
 357     for ( int itr = 0; itr < slice.NTracks(); itr++, sliceTr = sliceTr->GetNextTrack() ) {
 358 #endif
 359       AliHLTTPCGMSliceTrack &track = fSliceTrackInfos[nTracksCurrent];
 360       track.Set( sliceTr, alpha );
 361       if( !track.FilterErrors( fSliceParam, .999 ) ) continue;
 362       track.SetPrevNeighbour( -1 );
 363       track.SetNextNeighbour( -1 );
 364       track.SetSliceNeighbour( -1 );
 365       track.SetUsed( 0 );
 366 #ifdef GLOBAL_TRACKS_SPECIAL_TREATMENT
 367           track.SetGlobalTrackId(0, -1);
 368           track.SetGlobalTrackId(1, -1);
 369           TrackIds[iSlice * kMaxTrackIdInSlice + sliceTr->LocalTrackId()] = nTracksCurrent;
 370 #endif
 371           nTracksCurrent++;
 372       fSliceNTrackInfos[ iSlice ]++;
 373     }
 374     firstGlobalTracks[iSlice] = sliceTr;
 375
 376     //std::cout<<"Unpack slice "<<iSlice<<": ntracks "<<slice.NTracks()<<"/"<<fSliceNTrackInfos[iSlice]<<std::endl;
 377     }
 378 #ifdef GLOBAL_TRACKS_SPECIAL_TREATMENT
 379   for (int iSlice = 0;iSlice < fgkNSlices;iSlice++)
 380   {
 381           fSliceTrackGlobalInfoStart[iSlice] = nTracksCurrent;
 382           fSliceNGlobalTrackInfos[iSlice] = 0;
 383
 384           if ( !fkSlices[iSlice] ) continue;
 385       float alpha = fSliceParam.Alpha( iSlice );
 386
 387           const AliHLTTPCCASliceOutput &slice = *( fkSlices[iSlice] );
 388           const AliHLTTPCCASliceOutTrack *sliceTr = firstGlobalTracks[iSlice];
 389           for (int itr = slice.NLocalTracks();itr < slice.NTracks();itr++, sliceTr = sliceTr->GetNextTrack())
 390           {
 391                   if (TrackIds[sliceTr->LocalTrackId()] == -1) continue;
 392                   AliHLTTPCGMSliceTrack &track = fSliceTrackInfos[nTracksCurrent];
 393                   track.Set( sliceTr, alpha );
 394                   //if( !track.FilterErrors( fSliceParam, .999 ) ) continue;
 395                   track.SetPrevNeighbour( -1 );
 396                   track.SetNextNeighbour( -1 );
 397                   track.SetSliceNeighbour( -1 );
 398                   track.SetUsed( 0 );
 399                   track.SetLocalTrackId(TrackIds[sliceTr->LocalTrackId()]);
 400                   nTracksCurrent++;
 401                   fSliceNGlobalTrackInfos[ iSlice ]++;
 402           }
 403   }
 404
 405   free(TrackIds);
 406 #endif
 407 }
 408
 409
 410
 411
 412
 413 void AliHLTTPCGMMerger::MakeBorderTracks( int iSlice, int iBorder, AliHLTTPCGMBorderTrack B[], int &nB )
 414 {
 415   //* prepare slice tracks for merging with next/previous/same sector
 416   //* each track transported to the border line
 417
 418   float fieldBz = fSliceParam.ConstBz();
 419
 420   nB = 0;
 421
 422   float dAlpha = fSliceParam.DAlpha() / 2;
 423   float x0 = 0;
 424
 425   if ( iBorder == 0 ) { // transport to the left age of the sector and rotate horisontally
 426     dAlpha = dAlpha - CAMath::Pi() / 2 ;
 427   } else if ( iBorder == 1 ) { //  transport to the right age of the sector and rotate horisontally
 428     dAlpha = -dAlpha - CAMath::Pi() / 2 ;
 429   } else if ( iBorder == 2 ) { // transport to the left age of the sector and rotate vertically
 430     dAlpha = dAlpha;
 431     x0 = fSliceParam.RowX( 63 );
 432   } else if ( iBorder == 3 ) { // transport to the right age of the sector and rotate vertically
 433     dAlpha = -dAlpha;
 434     x0 =  fSliceParam.RowX( 63 );
 435   } else if ( iBorder == 4 ) { // transport to the middle of the sector, w/o rotation
 436     dAlpha = 0;
 437     x0 = fSliceParam.RowX( 63 );
 438   }
 439
 440   const float maxSin = CAMath::Sin( 60. / 180.*CAMath::Pi() );
 441   float cosAlpha = AliHLTTPCCAMath::Cos( dAlpha );
 442   float sinAlpha = AliHLTTPCCAMath::Sin( dAlpha );
 443
 444   for ( int itr = 0; itr < fSliceNTrackInfos[iSlice]; itr++ ) {
 445
 446     AliHLTTPCGMSliceTrack &track = fSliceTrackInfos[ fSliceTrackInfoStart[iSlice] + itr ];
 447
 448     if( track.Used() ) continue;
 449     AliHLTTPCGMBorderTrack &b = B[nB];
 450
 451     if(  track.TransportToXAlpha( x0, sinAlpha, cosAlpha, fieldBz, b, maxSin)){
 452       b.SetTrackID( itr );
 453       b.SetNClusters( track.NClusters() );
 454           for (int i = 0;i < 4;i++) if (fabs(b.Cov()[i]) >= 5.0) b.SetCov(i, 5.0);
 455           if (fabs(b.Cov()[4]) >= 0.5) b.SetCov(4, 0.5);
 456       nB++;
 457     }
 458   }
 459 }
 460
 461
 462 void AliHLTTPCGMMerger::MergeBorderTracks ( int iSlice1, AliHLTTPCGMBorderTrack B1[],  int N1,
 463                                             int iSlice2, AliHLTTPCGMBorderTrack B2[],  int N2 )
 464 {
 465   //* merge two sets of tracks
 466
 467   //std::cout<<" Merge slices "<<iSlice1<<"+"<<iSlice2<<": tracks "<<N1<<"+"<<N2<<std::endl;
 468   int statAll=0, statMerged=0;
 469   float factor2ys = 1.5;//1.5;//SG!!!
 470   float factor2zt = 1.5;//1.5;//SG!!!
 471   float factor2k = 2.0;//2.2;
 472
 473   factor2k  = 3.5 * 3.5 * factor2k * factor2k;
 474   factor2ys = 3.5 * 3.5 * factor2ys * factor2ys;
 475   factor2zt = 3.5 * 3.5 * factor2zt * factor2zt;
 476
 477   int minNPartHits = 10;//SG!!!
 478   int minNTotalHits = 20;
 479
 480   AliHLTTPCGMBorderTrack::Range *range1 = fBorderRangeMemory;
 481   AliHLTTPCGMBorderTrack::Range *range2 = fBorderRangeMemory + N1;
 482
 483   bool sameSlice = (iSlice1 == iSlice2);
 484   {
 485     for ( int itr = 0; itr < N1; itr++ ){
 486       AliHLTTPCGMBorderTrack &b = B1[itr];
 487       //   if( iSlice1==7 && iSlice2==8 ){
 488       //cout<<b.TrackID()<<": "<<b.Cov()[0]<<" "<<b.Cov()[1]<<endl;
 489       //}
 490       float d = 3.5*sqrt(b.Cov()[1]);
 491       range1[itr].fId = itr;
 492       range1[itr].fMin = b.Par()[1] - d;
 493       range1[itr].fMax = b.Par()[1] + d;
 494     }
 495     std::sort(range1,range1+N1,AliHLTTPCGMBorderTrack::Range::CompMin);
 496     if( sameSlice ){
 497       for(int i=0; i<N1; i++) range2[i]= range1[i];
 498       std::sort(range2,range2+N1,AliHLTTPCGMBorderTrack::Range::CompMax);
 499       N2 = N1;
 500       B2 = B1;
 501     }else{
 502       for ( int itr = 0; itr < N2; itr++ ){
 503         AliHLTTPCGMBorderTrack &b = B2[itr];
 504         float d = 3.5*sqrt(b.Cov()[1]);
 505         range2[itr].fId = itr;
 506         range2[itr].fMin = b.Par()[1] - d;
 507         range2[itr].fMax = b.Par()[1] + d;
 508       }
 509       std::sort(range2,range2+N2,AliHLTTPCGMBorderTrack::Range::CompMax);
 510     }
 511   }
 512
 513   int i2 = 0;
 514   for ( int i1 = 0; i1 < N1; i1++ ) {
 515
 516     AliHLTTPCGMBorderTrack::Range r1 = range1[i1];
 517     while( i2<N2 && range2[i2].fMax< r1.fMin ) i2++;
 518
 519     AliHLTTPCGMBorderTrack &b1 = B1[r1.fId];
 520     if ( b1.NClusters() < minNPartHits ) continue;
 521     int iBest2 = -1;
 522     int lBest2 = 0;
 523     statAll++;
 524     for( int k2 = i2; k2<N2; k2++){
 525
 526       AliHLTTPCGMBorderTrack::Range r2 = range2[k2];
 527       if( r2.fMin > r1.fMax ) break;
 528       if( sameSlice && (r1.fId >= r2.fId) ) continue;
 529       // do check
 530       AliHLTTPCGMBorderTrack &b2 = B2[r2.fId];
 531       if ( b2.NClusters() < lBest2 ) continue;
 532
 533       if( !b1.CheckChi2Y(b2, factor2ys ) ) continue;
 534       //if( !b1.CheckChi2Z(b2, factor2zt ) ) continue;
 535       if( !b1.CheckChi2QPt(b2, factor2k ) ) continue;
 536       if( !b1.CheckChi2YS(b2, factor2ys ) ) continue;
 537       if( !b1.CheckChi2ZT(b2, factor2zt ) ) continue;
 538       if ( b2.NClusters() < minNPartHits ) continue;
 539       if ( b1.NClusters() + b2.NClusters() < minNTotalHits ) continue;
 540
 541       lBest2 = b2.NClusters();
 542       iBest2 = b2.TrackID();
 543     }
 544
 545     if ( iBest2 < 0 ) continue;
 546     statMerged++;
 547     AliHLTTPCGMSliceTrack &newTrack1 = fSliceTrackInfos[fSliceTrackInfoStart[iSlice1] + b1.TrackID() ];
 548     AliHLTTPCGMSliceTrack &newTrack2 = fSliceTrackInfos[fSliceTrackInfoStart[iSlice2] + iBest2 ];
 549
 550     int old1 = newTrack2.PrevNeighbour();
 551
 552     if ( old1 >= 0 ) {
 553       AliHLTTPCGMSliceTrack &oldTrack1 = fSliceTrackInfos[fSliceTrackInfoStart[iSlice1] + old1];
 554       if ( oldTrack1.NClusters()  < newTrack1.NClusters() ) {
 555         newTrack2.SetPrevNeighbour( -1 );
 556         oldTrack1.SetNextNeighbour( -1 );
 557       } else continue;
 558     }
 559     int old2 = newTrack1.NextNeighbour();
 560     if ( old2 >= 0 ) {
 561       AliHLTTPCGMSliceTrack &oldTrack2 = fSliceTrackInfos[fSliceTrackInfoStart[iSlice2] + old2];
 562       if ( oldTrack2.NClusters() < newTrack2.NClusters() ) {
 563         oldTrack2.SetPrevNeighbour( -1 );
 564       } else continue;
 565     }
 566     newTrack1.SetNextNeighbour( iBest2 );
 567     newTrack2.SetPrevNeighbour( b1.TrackID() );
 568   }
 569   //cout<<"slices "<<iSlice1<<","<<iSlice2<<": all "<<statAll<<" merged "<<statMerged<<endl;
 570 }
 571
 572
 573 void AliHLTTPCGMMerger::MergeWithingSlices()
 574 {
 575   //* merge track segments withing one slice
 576
 577   float x0 = fSliceParam.RowX( 63 );
 578   const float maxSin = CAMath::Sin( 60. / 180.*CAMath::Pi() );
 579
 580   for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
 581
 582     int nBord = 0;
 583     for ( int itr = 0; itr < fSliceNTrackInfos[iSlice]; itr++ ) {
 584       AliHLTTPCGMSliceTrack &track = fSliceTrackInfos[ fSliceTrackInfoStart[iSlice] + itr ];
 585       //track.SetPrevNeighbour( -1 );
 586       //track.SetNextNeighbour( -1 );
 587       //track.SetSliceNeighbour( -1 );
 588       //track.SetUsed(0);
 589
 590       AliHLTTPCGMBorderTrack &b = fBorderMemory[nBord];
 591       if( track.TransportToX( x0, fSliceParam.ConstBz(), b, maxSin) ){
 592         b.SetTrackID( itr );
 593         b.SetNClusters( track.NClusters() );
 594         nBord++;
 595       }
 596     }
 597
 598     MergeBorderTracks( iSlice, fBorderMemory, nBord, iSlice, fBorderMemory, nBord );
 599
 600     for ( int itr = 0; itr < fSliceNTrackInfos[iSlice]; itr++ ) {
 601       AliHLTTPCGMSliceTrack &track = fSliceTrackInfos[ fSliceTrackInfoStart[iSlice] + itr];
 602       if( track.PrevNeighbour()>=0 || track.Used() ) continue;
 603       int jtr = track.NextNeighbour();
 604       track.SetSliceNeighbour( jtr );
 605       track.SetNextNeighbour(-1);
 606       while( jtr>=0 ){
 607         AliHLTTPCGMSliceTrack &trackN = fSliceTrackInfos[ fSliceTrackInfoStart[iSlice] + jtr];
 608         if( trackN.NClusters()>track.NClusters() ) track.CopyParamFrom(trackN);
 609         trackN.SetUsed(2);
 610         jtr = trackN.NextNeighbour();
 611         trackN.SetSliceNeighbour( jtr );
 612         trackN.SetNextNeighbour(-1);
 613         trackN.SetPrevNeighbour(-1);
 614       }
 615     }
 616   }
 617 }
 618
 619
 620
 621
 622 void AliHLTTPCGMMerger::MergeSlices()
 623 {
 624   //* track merging between slices
 625
 626   //for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
 627   //for ( int itr = 0; itr < fSliceNTrackInfos[iSlice]; itr++ ) {
 628   //AliHLTTPCGMSliceTrack &track = fSliceTrackInfos[ fSliceTrackInfoStart[iSlice] + itr ];
 629   //track.SetPrevNeighbour( -1 );
 630   //track.SetNextNeighbour( -1 );
 631   //}
 632   //}
 633
 634   AliHLTTPCGMBorderTrack
 635     *bCurr = fBorderMemory,
 636     *bNext = fBorderMemory + fMaxSliceTracks;
 637
 638   for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
 639     int jSlice = fNextSliceInd[iSlice];
 640     int nCurr = 0, nNext = 0;
 641     MakeBorderTracks( iSlice, 2, bCurr, nCurr );
 642     MakeBorderTracks( jSlice, 3, bNext, nNext );
 643     MergeBorderTracks( iSlice, bCurr, nCurr, jSlice, bNext, nNext );
 644     MakeBorderTracks( iSlice, 0, bCurr, nCurr );
 645     MakeBorderTracks( jSlice, 1, bNext, nNext );
 646     MergeBorderTracks( iSlice, bCurr, nCurr, jSlice, bNext, nNext );
 647   }
 648 }
 649
 650
 651
 652
 653
 654 void AliHLTTPCGMMerger::CollectMergedTracks()
 655 {
 656   //*
 657
 658   //for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
 659   //for ( int itr = 0; itr < fSliceNTrackInfos[iSlice]; itr++ ) {
 660   //AliHLTTPCGMSliceTrack &track = fSliceTrackInfos[ fSliceTrackInfoStart[iSlice] + itr ];
 661   //if( track.Used()!=2 ) track.SetUsed(0);
 662   //}
 663   //}
 664
 665   //Resolve connections for global tracks first
 666 #ifdef GLOBAL_TRACKS_SPECIAL_TREATMENT
 667   for (int iSlice = 0;iSlice < fgkNSlices;iSlice++)
 668   {
 669     for (int itr = 0;itr < fSliceNGlobalTrackInfos[iSlice];itr++)
 670         {
 671                 AliHLTTPCGMSliceTrack &globalTrack = fSliceTrackInfos[fSliceTrackGlobalInfoStart[iSlice] + itr];
 672                 AliHLTTPCGMSliceTrack &localTrack = fSliceTrackInfos[globalTrack.LocalTrackId()];
 673                 localTrack.SetGlobalTrackId(localTrack.GlobalTrackId(0) != -1, fSliceTrackGlobalInfoStart[iSlice] + itr);
 674         }
 675   }
 676 #endif
 677
 678   //Now collect the merged tracks
 679   fNOutputTracks = 0;
 680   int nOutTrackClusters = 0;
 681   const int kMaxParts = 400;
 682   const int kMaxClusters = 1000;
 683
 684   const AliHLTTPCGMSliceTrack *trackParts[kMaxParts];
 685
 686   for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
 687
 688     for ( int itr = 0; itr < fSliceNTrackInfos[iSlice]; itr++ ) {
 689
 690       AliHLTTPCGMSliceTrack &track = fSliceTrackInfos[fSliceTrackInfoStart[iSlice] + itr];
 691
 692       if ( track.Used() ) continue;
 693       if ( track.PrevNeighbour() >= 0 ) continue;
 694       int nParts = 0;
 695       int jSlice = iSlice;
 696       AliHLTTPCGMSliceTrack *trbase = &track, *tr = &track;
 697       tr->SetUsed( 1 );
 698       do{
 699             if( nParts >= kMaxParts ) break;
 700             trackParts[nParts++] = tr;
 701 #ifdef GLOBAL_TRACKS_SPECIAL_TREATMENT
 702                 for (int i = 0;i < 2;i++) if (tr->GlobalTrackId(i) != -1) trackParts[nParts++] = &fSliceTrackInfos[tr->GlobalTrackId(i)];
 703 #endif
 704             int jtr = tr->SliceNeighbour();
 705             if( jtr >= 0 ) {
 706               tr = &(fSliceTrackInfos[fSliceTrackInfoStart[jSlice] + jtr]);
 707               tr->SetUsed( 2 );
 708               continue;
 709             }
 710         jtr = trbase->NextNeighbour();
 711             if( jtr>=0 ){
 712               jSlice = fNextSliceInd[jSlice];
 713               trbase = &(fSliceTrackInfos[fSliceTrackInfoStart[jSlice] + jtr]);
 714               tr = trbase;
 715               if( tr->Used() ) break;
 716               tr->SetUsed( 1 );
 717               continue;
 718             }
 719             break;
 720       }while(1);
 721
 722       // unpack and sort clusters
 723
 724           std::sort(trackParts, trackParts+nParts, CompareTrackParts );
 725
 726       AliHLTTPCCASliceOutCluster tmp[kMaxClusters];
 727       int currCluster = nOutTrackClusters;
 728       int nHits = 0;
 729       for( int ipart=0; ipart<nParts; ipart++ ){
 730         const AliHLTTPCGMSliceTrack *t = trackParts[ipart];
 731         int nTrackHits = t->NClusters();
 732         if( nHits + nTrackHits >= kMaxClusters ) break;
 733         const AliHLTTPCCASliceOutCluster *c= t->OrigTrack()->Clusters();
 734         AliHLTTPCCASliceOutCluster *c2 = tmp+nHits + nTrackHits-1;
 735         for( int i=0; i<nTrackHits; i++, c++, c2-- ) *c2 = *c;
 736         float alpha =  t->Alpha();
 737         for( int i=0; i<nTrackHits; i++) fClusterAngle[currCluster++] = alpha;
 738         nHits+=nTrackHits;
 739       }
 740
 741       if ( nHits < 30 ) continue;
 742
 743           float *clX = fClusterX + nOutTrackClusters;
 744           clX[0] = tmp[0].GetX();
 745           int ordered = 1;
 746       for( int i=1; i<nHits; i++ )
 747           {
 748                   clX[i] = tmp[i].GetX();
 749                   if (clX[i] > clX[i - 1]) ordered = 0;
 750           }
 751
 752           if (ordered == 0)
 753           {
 754                   int2 tmpFilter[kMaxClusters];
 755                   for( int i = 0;i < nHits;i++)
 756                   {
 757                           tmpFilter[i].x = i;
 758                           tmpFilter[i].y = tmp[i].GetId();
 759                   }
 760                   qsort(tmpFilter, nHits, sizeof(int2), CompareClusterIds);
 761                   bool tmpHitsOk[kMaxClusters];
 762                   tmpHitsOk[tmpFilter[0].x] = true;
 763                   for (int i = 1;i < nHits;i++)
 764                   {
 765                         tmpHitsOk[tmpFilter[i].x] = (tmpFilter[i].y != tmpFilter[i - 1].y);
 766                   }
 767                   int nFilteredHits = 0;
 768                   float *clA = fClusterAngle + nOutTrackClusters;
 769                   for (int i = 0;i < nHits;i++)
 770                   {
 771                         if (tmpHitsOk[i])
 772                         {
 773                                 if (nFilteredHits < i)
 774                                 {
 775                                         tmp[nFilteredHits] = tmp[i];
 776                                         clA[nFilteredHits] = clA[i];
 777                                 }
 778                                 nFilteredHits++;
 779                         }
 780                   }
 781
 782                   nHits = nFilteredHits;
 783                   for( int i=0; i<nHits; i++ ) clX[i] = tmp[i].GetX();
 784           }
 785
 786           UInt_t *clId = fOutputClusterIds + nOutTrackClusters;
 787       for( int i=0; i<nHits; i++ ) clId[i] = tmp[i].GetId();
 788
 789       UInt_t *clT  = fClusterRowType + nOutTrackClusters;
 790       for( int i=0; i<nHits; i++ ) clT[i] = tmp[i].GetRowType();
 791
 792       float *clY = fClusterY + nOutTrackClusters;
 793       for( int i=0; i<nHits; i++ ) clY[i] = tmp[i].GetY();
 794
 795       float *clZ = fClusterZ + nOutTrackClusters;
 796       for( int i=0; i<nHits; i++ ) clZ[i] = tmp[i].GetZ();
 797
 798       AliHLTTPCGMMergedTrack &mergedTrack = fOutputTracks[fNOutputTracks];
 799       mergedTrack.SetOK(1);
 800       mergedTrack.SetNClusters( nHits );
 801       mergedTrack.SetFirstClusterRef( nOutTrackClusters );
 802       AliHLTTPCGMTrackParam &p1 = mergedTrack.Param();
 803       const AliHLTTPCGMSliceTrack &p2 = *(trackParts[0]);
 804
 805       p1.X() = p2.X();
 806       p1.Y() = p2.Y();
 807       p1.Z() = p2.Z();
 808       p1.SinPhi()  = p2.SinPhi();
 809       p1.DzDs()  = p2.DzDs();
 810       p1.QPt()  = p2.QPt();
 811       mergedTrack.SetAlpha( p2.Alpha() );
 812
 813       fNOutputTracks++;
 814       nOutTrackClusters += nHits;
 815     }
 816   }
 817   fNOutputTrackClusters = nOutTrackClusters;
 818 }
 819
 820 void AliHLTTPCGMMerger::Refit()
 821 {
 822   //* final refit
 823 #ifdef HLTCA_GPU_MERGER
 824         if (fGPUTracker)
 825         {
 826                 fGPUTracker->RefitMergedTracks(this);
 827         }
 828         else
 829 #endif
 830         {
 831 #ifdef HLTCA_STANDALONE
 832 #pragma omp parallel for
 833 #endif
 834           for ( int itr = 0; itr < fNOutputTracks; itr++ ) {
 835
 836                 AliHLTTPCGMMergedTrack &track = fOutputTracks[itr];
 837                 if( !track.OK() ) continue;
 838
 839                 int nTrackHits = track.NClusters();
 840
 841                 AliHLTTPCGMTrackParam t = track.Param();
 842                 float Alpha = track.Alpha();
 843
 844                 t.Fit( fPolinomialFieldBz,
 845                    fClusterX+track.FirstClusterRef(),
 846                    fClusterY+track.FirstClusterRef(),
 847                    fClusterZ+track.FirstClusterRef(),
 848                    fClusterRowType+track.FirstClusterRef(),
 849                    fClusterAngle+track.FirstClusterRef(),
 850                    fSliceParam, nTrackHits, Alpha, 0 );
 851
 852                 if ( fabs( t.QPt() ) < 1.e-4 ) t.QPt() = 1.e-4 ;
 853
 854                 bool ok = nTrackHits >= 30 && t.CheckNumericalQuality() && fabs( t.SinPhi() ) <= .999;
 855
 856                 track.SetOK(ok);
 857                 if (!ok) continue;
 858
 859                 if( 1 ){//SG!!!
 860                   track.SetNClusters( nTrackHits );
 861                   track.Param() = t;
 862                   track.Alpha() = Alpha;
 863                 }
 864
 865                 {
 866                   int ind = track.FirstClusterRef();
 867                   float alpha = fClusterAngle[ind];
 868                   float x = fClusterX[ind];
 869                   float y = fClusterY[ind];
 870                   float z = fClusterZ[ind];
 871                   float sinA = AliHLTTPCCAMath::Sin( alpha - track.Alpha());
 872                   float cosA = AliHLTTPCCAMath::Cos( alpha - track.Alpha());
 873                   track.SetLastX( x*cosA - y*sinA );
 874                   track.SetLastY( x*sinA + y*cosA );
 875                   track.SetLastZ( z );
 876                 }
 877           }
 878         }
 879 }
 880