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 #include "AliHLTTPCCATracker.h"
  20
  21 #include "AliHLTTPCCAHit.h"
  22 #include "AliHLTTPCCACell.h"
  23 #include "AliHLTTPCCAEndPoint.h"
  24 #include "AliHLTTPCCAOutTrack.h"
  25 #include "AliHLTTPCCAGrid.h"
  26 #include "AliHLTTPCCARow.h"
  27 #include "AliHLTTPCCATrack.h"
  28
  29 #include "TMath.h"
  30 #include "Riostream.h"
  31 //#include <algo.h>
  32 #include "TStopwatch.h"
  33
  34 //#define DRAW
  35
  36 #ifdef DRAW
  37   #include "AliHLTTPCCADisplay.h"
  38   #include "TApplication.h"
  39 #endif //DRAW
  40
  41 ClassImp(AliHLTTPCCATracker)
  42
  43
  44 AliHLTTPCCATracker::AliHLTTPCCATracker()
  45   :fParam(),fRows(0),fOutTrackHits(0),fNOutTrackHits(0),fOutTracks(0),fNOutTracks(0),fNHitsTotal(0),fTracks(0),fNTracks(0),fCellHitPointers(0),fCells(0),fEndPoints(0)
  46 {
  47   // constructor
  48   fRows = new AliHLTTPCCARow[fParam.NRows()];
  49   Initialize( fParam );
  50 }
  51
  52 AliHLTTPCCATracker::AliHLTTPCCATracker( const AliHLTTPCCATracker& )
  53   :fParam(),fRows(0),fOutTrackHits(0),fNOutTrackHits(0),fOutTracks(0),fNOutTracks(0),fNHitsTotal(0),fTracks(0),fNTracks(0),fCellHitPointers(0),fCells(0),fEndPoints(0)
  54 {
  55   // dummy
  56 }
  57
  58 AliHLTTPCCATracker &AliHLTTPCCATracker::operator=( const AliHLTTPCCATracker& )
  59 {
  60   // dummy
  61   fRows=0;
  62   fOutTrackHits=0;
  63   fOutTracks=0;
  64   fNOutTracks=0;
  65   return *this;
  66 }
  67
  68 AliHLTTPCCATracker::~AliHLTTPCCATracker()
  69 {
  70   // destructor
  71   StartEvent();
  72   delete[] fRows;
  73 }
  74
  75 // ----------------------------------------------------------------------------------
  76 void AliHLTTPCCATracker::Initialize( AliHLTTPCCAParam &param )
  77 {
  78   // initialosation
  79   StartEvent();
  80   delete[] fRows;
  81   fRows = 0;
  82   fParam = param;
  83   fParam.Update();
  84   fRows = new AliHLTTPCCARow[fParam.NRows()];
  85   Float_t xStep = 1;
  86   Float_t deltaY = TMath::Tan(fParam.CellConnectionAngleXY());
  87   Float_t deltaZ = TMath::Tan(fParam.CellConnectionAngleXZ());
  88   for( Int_t irow=0; irow<fParam.NRows(); irow++ ){
  89     fRows[irow].X() = fParam.RowX(irow);
  90     if( irow < fParam.NRows()-1 ) xStep = fParam.RowX(irow+1) - fParam.RowX(irow);
  91     fRows[irow].DeltaY() = xStep*deltaY;
  92     fRows[irow].DeltaZ() = xStep*deltaZ;
  93     fRows[irow].MaxY() = TMath::Tan( fParam.DAlpha()/2.)*fRows[irow].X();
  94   }
  95   StartEvent();
  96 }
  97
  98 void AliHLTTPCCATracker::StartEvent()
  99 {
 100   // start new event and fresh the memory
 101
 102   if( fTracks ) delete[] fTracks;
 103   if( fOutTrackHits ) delete[] fOutTrackHits;
 104   if( fOutTracks ) delete[] fOutTracks;
 105   if( fCellHitPointers ) delete[] fCellHitPointers;
 106   if( fCells ) delete[] fCells;
 107   if( fEndPoints ) delete[] fEndPoints;
 108   fTracks = 0;
 109   fOutTrackHits = 0;
 110   fOutTracks = 0;
 111   fCellHitPointers = 0;
 112   fCells = 0;
 113   fEndPoints = 0;
 114   fNTracks = 0;
 115   fNOutTrackHits = 0;
 116   fNOutTracks = 0;
 117   fNHitsTotal = 0;
 118   for( Int_t irow=0; irow<fParam.NRows(); irow++ ){
 119     fRows[irow].Clear();
 120   }
 121 }
 122
 123
 124 void AliHLTTPCCATracker::ReadHitRow( Int_t iRow, AliHLTTPCCAHit *Row, Int_t NHits )
 125 {
 126   // read row of hits
 127   AliHLTTPCCARow &row = fRows[iRow];
 128   row.Hits() = new AliHLTTPCCAHit[NHits];
 129   for( Int_t i=0; i<NHits; i++ ){
 130     row.Hits()[i]=Row[i];
 131     row.Hits()[i].ErrY()*= fParam.YErrorCorrection();
 132     row.Hits()[i].ErrZ()*= fParam.ZErrorCorrection();
 133   }
 134   row.NHits() = NHits;
 135   fNHitsTotal += NHits;
 136 }
 137
 138 void AliHLTTPCCATracker::Reconstruct()
 139 {
 140   //* reconstruction of event
 141
 142 #ifdef DRAW
 143   if( !gApplication ){
 144     TApplication *myapp = new TApplication("myapp",0,0);
 145   }
 146   //AliHLTTPCCADisplay::Instance().Init();
 147
 148   AliHLTTPCCADisplay::Instance().SetCurrentSlice( this );
 149   AliHLTTPCCADisplay::Instance().SetSliceView();
 150   AliHLTTPCCADisplay::Instance().DrawSlice( this );
 151   //for( Int_t iRow=0; iRow<fParam.NRows(); iRow++ )
 152   //for (Int_t i = 0; i<fRows[iRow].NHits(); i++)
 153   //AliHLTTPCCADisplay::Instance().DrawHit( iRow, i );
 154   //AliHLTTPCCADisplay::Instance().Ask();
 155 #endif
 156
 157   fTimers[0] = 0;
 158   fTimers[1] = 0;
 159   fTimers[2] = 0;
 160   fTimers[3] = 0;
 161   fTimers[4] = 0;
 162   fTimers[5] = 0;
 163   fTimers[6] = 0;
 164   fTimers[7] = 0;
 165   if( fNHitsTotal < 1 ) return;
 166
 167   //cout<<"Find Cells..."<<endl;
 168   FindCells();
 169   //cout<<"Merge Cells..."<<endl;
 170   MergeCells();
 171   //cout<<"Find Tracks..."<<endl;
 172   FindTracks();
 173   //cout<<"Find Tracks OK"<<endl;
 174  }
 175
 176
 177 void AliHLTTPCCATracker::FindCells()
 178 {
 179   //* cell finder - neighbouring hits are grouped to cells
 180
 181   TStopwatch timer;
 182
 183   fCellHitPointers = new Int_t [fNHitsTotal];
 184   fCells = new AliHLTTPCCACell[fNHitsTotal];
 185   fEndPoints = new AliHLTTPCCAEndPoint[fNHitsTotal];
 186
 187   struct THitCont{
 188     Float_t Ymin, Ymax, Zmin, Zmax;
 189     Int_t binYmin, binYmax, binZmin, binZmax;
 190     Bool_t used;
 191     AliHLTTPCCAHit *h;
 192     THitCont *next;
 193   };
 194   THitCont *hitCont = new THitCont[fNHitsTotal];
 195
 196   Int_t lastCellHitPointer = 0;
 197   Int_t lastCell = 0;
 198
 199   for( Int_t irow=0; irow<fParam.NRows(); irow++ ){
 200     AliHLTTPCCARow &row=fRows[irow];
 201     Int_t nHits = row.NHits();
 202     //cout<<"row = "<<irow<<", x="<<row.X()<<endl;
 203     if( nHits<1 ) continue;
 204     //cout<<nHits*sizeof(AliHLTTPCCAHit)/1024.<<endl;
 205
 206     Float_t deltaY = row.DeltaY();
 207     Float_t deltaZ = row.DeltaZ();
 208
 209     Float_t yMin = 1.e20, zMin = 1.e20, yMax = -1.e20, zMax = -1.e20;
 210     for (Int_t ih = 0; ih<nHits; ih++){
 211       AliHLTTPCCAHit &h  = row.Hits()[ih];
 212       if( yMin> h.Y() ) yMin = h.Y();
 213       if( yMax< h.Y() ) yMax = h.Y();
 214       if( zMin> h.Z() ) zMin = h.Z();
 215       if( zMax< h.Z() ) zMax = h.Z();
 216     }
 217     AliHLTTPCCAGrid grid;
 218     grid.Create( yMin, yMax, zMin, zMax, nHits );
 219     //cout<<"row "<<irow<<", delta = "<<delta<<" :\n"<<endl;
 220     for (Int_t ih = 0; ih<nHits; ih++){
 221       AliHLTTPCCAHit &h  = row.Hits()[ih];
 222       THitCont &cont = hitCont[ih];
 223       THitCont *&bin = * ((THitCont **) grid.GetNoCheck(h.Y(),h.Z()));
 224       cont.h = &h;
 225       cont.used = 0;
 226       Float_t y = h.Y();
 227       //cout<<"ih = "<<ih<<", y= "<<y<<endl;
 228       Float_t dY = 3.5*h.ErrY() + deltaY;
 229       cont.Ymin = y-dY;
 230       cont.Ymax = y+dY;
 231       Float_t z = h.Z();
 232       Float_t dZ = 3.5*h.ErrZ() + deltaZ;
 233       cont.Zmin = z-dZ;
 234       cont.Zmax = z+dZ;
 235       cont.binYmin = (Int_t ) ( (cont.Ymin-dY-grid.YMin())*grid.StepYInv() );
 236       cont.binYmax = (Int_t ) ( (cont.Ymax+dY-grid.YMin())*grid.StepYInv() );
 237       cont.binZmin = (Int_t ) ( (cont.Zmin-dZ-grid.ZMin())*grid.StepZInv() );
 238       cont.binZmax = (Int_t ) ( (cont.Zmax+dZ-grid.ZMin())*grid.StepZInv() );
 239       if( cont.binYmin<0 ) cont.binYmin = 0;
 240       if( cont.binYmin>=grid.Ny() ) cont.binYmin = grid.Ny()-1;
 241       if( cont.binYmax<0 ) cont.binYmax = 0;
 242       if( cont.binYmax>=grid.Ny() ) cont.binYmax = grid.Ny()-1;
 243       if( cont.binZmin<0 ) cont.binZmin = 0;
 244       if( cont.binZmin>=grid.Nz() ) cont.binZmin = grid.Nz()-1;
 245       if( cont.binZmax<0 ) cont.binZmax = 0;
 246       if( cont.binZmax>=grid.Nz() ) cont.binZmax = grid.Nz()-1;
 247       cont.next = bin;
 248       bin = &cont;
 249     }
 250
 251     row.CellHitPointers() = fCellHitPointers + lastCellHitPointer;
 252     row.Cells() = fCells + lastCell;
 253     Int_t nPointers = 0;
 254     Int_t nCells = 0;
 255
 256     //Int_t statMaxBins = 0;
 257     //Int_t statMaxHits = 0;
 258     for (Int_t ih = 0; ih<nHits; ih++){
 259       THitCont &cont = hitCont[ih];
 260       if( cont.used ) continue;
 261       // cell start
 262       AliHLTTPCCACell &cell = row.Cells()[nCells++];
 263       cell.FirstHitRef() = nPointers;
 264       cell.NHits() = 1;
 265       cell.Link() = -1;
 266       cell.Status() = 0;
 267       cell.TrackID() = -1;
 268       row.CellHitPointers()[nPointers++] = ih;
 269       cont.used = 1;
 270
 271 #ifdef DRAW
 272       //AliHLTTPCCADisplay::Instance().DrawHit( irow, ih, kRed );
 273 #endif
 274   // cell finder - neighbouring hits are grouped to cells
 275
 276
 277       Float_t ymin = cont.Ymin;
 278       Float_t ymax = cont.Ymax;
 279       Float_t zmin = cont.Zmin;
 280       Float_t zmax = cont.Zmax;
 281       Int_t binYmin = cont.binYmin;
 282       Int_t binYmax = cont.binYmax;
 283       Int_t binZmin = cont.binZmin;
 284       Int_t binZmax = cont.binZmax;
 285
 286       Bool_t repeat = 1;
 287       while( repeat ){
 288         repeat = 0;
 289         THitCont ** startY = (THitCont **) grid.Grid() + binZmin*grid.Ny();
 290         //Float_t Ymax1 = Ymax;
 291         //Float_t Ymin1 = Ymin;
 292         //Float_t Zmax1 = Zmax;
 293         //Float_t Zmin1 = Zmin;
 294         Int_t binYmax1 = binYmax;
 295         Int_t binYmin1 = binYmin;
 296         Int_t binZmax1 = binZmax;
 297         Int_t binZmin1 = binZmin;
 298 #ifdef DRAW
 299         //cell.Y() = .5*(Ymin+Ymax);
 300         //cell.Z() = .5*(Zmin+Zmax);
 301         //cell.ErrY() = .5*( Ymax - Ymin )/3.5;
 302         //cell.ErrZ() = .5*( Zmax - Zmin )/3.5;
 303         //cell.YMin() = Ymin;
 304         //cell.YMax() = Ymax;
 305         //AliHLTTPCCADisplay::Instance().DrawCell( irow, nCells-1, 1,kRed );
 306         //AliHLTTPCCADisplay::Instance().Ask();
 307 #endif
 308         for( Int_t iGridZ=binZmin1; iGridZ<=binZmax1; iGridZ++, startY += grid.Ny() ){
 309           for( Int_t iGridY=binYmin1; iGridY<=binYmax1; iGridY++ ){
 310             for( THitCont *bin = *(startY + iGridY); bin; bin=bin->next ){
 311               Int_t jh = bin->h-row.Hits();
 312               THitCont &cont1 = hitCont[jh];
 313               if( cont1.used ) continue;
 314               //cout<<"["<<Ymin<<","<<Ymax<<"]: ["<<cont1.Ymin<<","<<cont1.Ymax<<"]"<<endl;
 315               if( cont1.Ymax < ymin ) continue;
 316               if( cont1.Ymin > ymax ) continue;
 317               if( cont1.Zmax < zmin ) break;// in the grid cell hit Y is decreasing
 318               if( cont1.Zmin > zmax ) continue;
 319
 320               if( cont1.Ymin < ymin ){ ymin = cont1.Ymin; repeat = 1; }
 321               if( cont1.Ymax > ymax ){ ymax = cont1.Ymax; repeat = 1; }
 322               if( cont1.Zmin < zmin ){ zmin = cont1.Zmin; repeat = 1; }
 323               if( cont1.Zmax > zmax ){ zmax = cont1.Zmax; repeat = 1; }
 324               if( cont1.binYmin < binYmin ){ binYmin = cont1.binYmin; repeat = 1; }
 325               if( cont1.binYmax > binYmax ){ binYmax = cont1.binYmax; repeat = 1; }
 326               if( cont1.binZmin < binZmin ){ binZmin = cont1.binZmin; repeat = 1; }
 327               if( cont1.binZmax > binZmax ){ binZmax = cont1.binZmax; repeat = 1; }
 328
 329               row.CellHitPointers()[nPointers++] = jh;
 330               cell.NHits()++;
 331               cont1.used = 1;
 332 #ifdef DRAW
 333               //AliHLTTPCCADisplay::Instance().DrawHit( irow, jh, kRed );
 334               //AliHLTTPCCADisplay::Instance().Ask();
 335 #endif
 336             }
 337           }
 338         }
 339       }
 340
 341       cell.Y() = .5*(ymin+ymax);
 342       cell.Z() = .5*(zmin+zmax);
 343       cell.ErrY() = .5*( ymax - ymin - 2*deltaY)/3.5;
 344       cell.ErrZ() = .5*( zmax - zmin -2*deltaZ)/3.5;
 345       cell.ZMin() = zmin;
 346       cell.ZMax() = zmax;
 347 #ifdef DRAW
 348       //AliHLTTPCCADisplay::Instance().DrawCell( irow, nCells-1 );
 349       //AliHLTTPCCADisplay::Instance().Ask();
 350 #endif
 351     }
 352     //cout<<statMaxBins<<"/"<<grid.N()<<" "<<statMaxHits<<"/"<<nHits<<endl;
 353
 354
 355     row.NCells() = nCells;
 356     lastCellHitPointer += nPointers;
 357     lastCell += nCells;
 358   }
 359   delete[] hitCont;
 360   timer.Stop();
 361   fTimers[0] = timer.CpuTime();
 362 }
 363
 364
 365 void AliHLTTPCCATracker::MergeCells()
 366 {
 367   // First step:
 368   //  for each Cell find one neighbour in the next row (irow+1)
 369   //  when there are no neighbours, look to the rows (irow+2),(irow+3)
 370   //
 371   // Initial state:
 372   //  cell.Link  =-1
 373   //  cell.Link1  =-1
 374   //  cell.Track      =-1
 375   //
 376   // Intermediate state: same as final
 377   //
 378   // Final state:
 379   //  cell.Link = Neighbour ID, if there is a neighbour
 380   //             = -1, if no neighbours found
 381   //             = -2, if there was more than one neighbour
 382   //  cell.Link1 = ID of the cell which has this Cell as a forward neighbour
 383   //             = -1 there are no backward neighbours
 384   //             = -2 there are more than one neighbour
 385   //  cell.Track = -1
 386   //
 387
 388   TStopwatch timer;
 389
 390   Int_t nStartCells = 0;
 391   for( Int_t iRow1=0; iRow1<fParam.NRows(); iRow1++ ){
 392     AliHLTTPCCARow &row1 = fRows[iRow1];
 393
 394     Float_t deltaY = row1.DeltaY();
 395     Float_t deltaZ = row1.DeltaZ();
 396     Float_t xStep = 1;
 397     if( iRow1 < fParam.NRows()-1 ) xStep = fParam.RowX(iRow1+1) - fParam.RowX(iRow1);
 398     Float_t tx = xStep/row1.X();
 399
 400     Int_t lastRow2 = iRow1+3;
 401     if( lastRow2>=fParam.NRows() ) lastRow2 = fParam.NRows()-1;
 402
 403     for (Int_t i1 = 0; i1<row1.NCells(); i1++){
 404       AliHLTTPCCACell &c1  = row1.Cells()[i1];
 405       //cout<<"row, cell= "<<iRow1<<" "<<i1<<" "<<c1.Y()<<" "<<c1.ErrY()<<" "<<c1.Z()<<" "<<c1.ErrZ()<<endl;
 406       //Float_t sy1 = c1.ErrY()*c1.ErrY();
 407       Float_t yy = c1.Y() +tx*c1.Y();
 408       Float_t zz = c1.Z() +tx*c1.Z();
 409
 410       Float_t yMin = yy - 3.5*c1.ErrY() - deltaY;
 411       Float_t yMax = yy + 3.5*c1.ErrY() + deltaY;
 412       Float_t zMin = zz - 3.5*c1.ErrZ() - deltaZ;
 413       Float_t zMax = zz + 3.5*c1.ErrZ() + deltaZ;
 414       //Float_t sz1 = c1.ErrZ()*c1.ErrZ();
 415       if( c1.Status()<=0 ) nStartCells++;
 416
 417       // looking for neighbour for the Cell c1
 418       Bool_t found = 0;
 419       for( Int_t iRow2=iRow1+1; iRow2<=lastRow2&&(!found); iRow2++ ){
 420         AliHLTTPCCARow &row2 = fRows[iRow2];
 421         AliHLTTPCCACell *cc2 = lower_bound(row2.Cells(),row2.Cells()+row2.NCells(),zMin,AliHLTTPCCARow::CompareCellZMax);
 422         for (Int_t i2 = (cc2 - row2.Cells()); i2<row2.NCells(); i2++){
 423           //cout<<"   candidat = "<<iRow2<<" "<<i2<<endl;
 424
 425           AliHLTTPCCACell &c2  = row2.Cells()[i2];
 426           Float_t y2Min = c2.Y() - 3.5*c2.ErrY();
 427           Float_t y2Max = c2.Y() + 3.5*c2.ErrY();
 428           Float_t z2Min = c2.Z() - 3.5*c2.ErrZ();
 429           Float_t z2Max = c2.Z() + 3.5*c2.ErrZ();
 430
 431           if( y2Min > yMax ) continue;
 432           if( y2Max < yMin ) continue;
 433           if( z2Min > zMax ) break;
 434           if( z2Max < zMin ) continue;
 435
 436           // c1 & c2 are neighbours
 437
 438           found = 1;
 439
 440           if( c1.Link() ==-1 && c2.Status()==0 ){
 441             // one-to-one connection - OK
 442             c1.Link() = IRowICell2ID(iRow2,i2);
 443             c2.Status() = 1;
 444           }else{
 445             // multi-connection - break all links
 446             if( c1.Link()>=0 ) ID2Cell(c1.Link()).Status() = -1;
 447             c1.Link()  = -2;
 448             c2.Status() = -1;
 449           }
 450         }
 451       }//row2
 452     }
 453   }//row1
 454
 455   timer.Stop();
 456   fTimers[1] = timer.CpuTime();
 457
 458   // Second step: create tracks
 459   //  for each sequence of neighbouring Cells create Track object
 460   //
 461   // Final state:
 462   //  cell.Track     = TrackNumber for first and last track cell
 463   //                 = -1 for other cells
 464   //  cell.Link     = Neighbour ID, if there is a neighbour
 465   //                 = -1, if no neighbour (last cell on the track )
 466   //  cell.Link1     = backward neighbour ID, if there is a neighbour
 467   //                 = -1 for first and last track cells
 468
 469   TStopwatch timer2;
 470
 471   fTracks = new AliHLTTPCCATrack[nStartCells];
 472   fNTracks = 0;
 473
 474   for( Int_t iRow=0; iRow<fParam.NRows(); iRow++ ){
 475     AliHLTTPCCARow &row = fRows[iRow];
 476     for( Int_t iCell = 0; iCell<row.NCells(); iCell++){
 477       AliHLTTPCCACell &c  = row.Cells()[iCell];
 478       if( c.Status()>0 ) continue; // not a starting cell
 479
 480       Int_t firstID = IRowICell2ID( iRow, iCell );
 481       Int_t midID = firstID;
 482       Int_t lastID = firstID;
 483
 484       AliHLTTPCCATrack &track = fTracks[fNTracks];
 485       track.Alive() = 1;
 486       track.NCells() = 1;
 487       AliHLTTPCCACell *last = &c;
 488       while( last->Link() >=0 ){
 489         Int_t nextID = last->Link();
 490         AliHLTTPCCACell *next = & ID2Cell(nextID);
 491         if(next->Status()!=1 ){
 492           last->Link() = -1;
 493           break;
 494         }
 495         track.NCells()++;
 496         last = next;
 497         lastID = nextID;
 498       }
 499       Int_t nCells05 = (track.NCells()-1)/2;
 500       for( Int_t i=0; i<nCells05; i++ ) midID = ID2Cell(midID).Link();
 501       //cout<<fNTracks<<", NCells="<<track.NCells()<<" "<<nCells05<<"id="<<firstID<<" "<<midID<<" "<<lastID<<endl;
 502       c.TrackID() = fNTracks;
 503       last->TrackID() = fNTracks;
 504       track.FirstCellID() = firstID;
 505       track.CellID()[0] = firstID;
 506       track.CellID()[1] = midID;
 507       track.CellID()[2] = lastID;
 508       track.PointID()[0] = -1;
 509       track.PointID()[1] = -1;
 510       //cout<<"Track N "<<fNTracks<<", NCells="<<track.NCells()<<endl;
 511
 512       fNTracks++;
 513     }
 514   }
 515   if( fNTracks != nStartCells ){
 516     //cout<<"fNTracks="<<fNTracks<<", NStrartCells="<<nStartCells<<endl;
 517     //exit(0);
 518     return;
 519   }
 520
 521   // create endpoints
 522
 523   Int_t nEndPointsTotal = 0;
 524   for( Int_t iRow=0; iRow<fParam.NRows(); iRow++ ){
 525     AliHLTTPCCARow &row = fRows[iRow];
 526     row.EndPoints()= fEndPoints + nEndPointsTotal;
 527     row.NEndPoints()=0;
 528     for( Int_t iCell = 0; iCell<row.NCells(); iCell++){
 529       AliHLTTPCCACell &c  = row.Cells()[iCell];
 530       if( c.TrackID()< 0 ) continue; // not an endpoint
 531       AliHLTTPCCAEndPoint &p = row.EndPoints()[row.NEndPoints()];
 532       p.CellID() = IRowICell2ID(iRow,iCell);
 533       p.TrackID() = c.TrackID();
 534       p.Link() = -1;
 535       AliHLTTPCCATrack &track = fTracks[c.TrackID()];
 536       if( c.Link()>=0 ){
 537         track.PointID()[0] = IRowICell2ID(iRow,row.NEndPoints());
 538       }else{
 539         track.PointID()[1] = IRowICell2ID(iRow,row.NEndPoints());
 540         if( track.PointID()[0]<0 )track.PointID()[0] = track.PointID()[1];
 541       }
 542       row.NEndPoints()++;
 543     }
 544     nEndPointsTotal += row.NEndPoints();
 545   }
 546   timer2.Stop();
 547   fTimers[2] = timer2.CpuTime();
 548 }
 549
 550 void AliHLTTPCCATracker::FindTracks()
 551 {
 552   // the Cellular Automaton track finder
 553   TStopwatch timer3;
 554   //cout<<"combine & fit tracks"<<endl;
 555   for( Int_t itr=0; itr<fNTracks; itr++ ){
 556     AliHLTTPCCATrack &iTrack = fTracks[itr];
 557     //if( iTrack.NCells()<3 ) continue;
 558     //cout<<" fit track "<<itr<<", NCells="<<iTrack.NCells()<<endl;
 559     ID2Point(iTrack.PointID()[0]).Param().CosPhi() = -1;
 560     ID2Point(iTrack.PointID()[1]).Param().CosPhi() = 1;
 561     FitTrack( iTrack );
 562     //if( iTrack.Param().Chi2() > fParam.TrackChi2Cut()*iTrack.Param().NDF() ){
 563       //iTrack.Alive() = 0;
 564     //}
 565   }
 566   timer3.Stop();
 567   fTimers[3] = timer3.CpuTime();
 568
 569 #ifdef DRAW
 570   if( !gApplication ){
 571     TApplication *myapp = new TApplication("myapp",0,0);
 572   }
 573   //AliHLTTPCCADisplay::Instance().Init();
 574
 575   AliHLTTPCCADisplay::Instance().SetCurrentSlice( this );
 576   AliHLTTPCCADisplay::Instance().DrawSlice( this );
 577   cout<<"draw hits..."<<endl;
 578   for( Int_t iRow=0; iRow<fParam.NRows(); iRow++ )
 579     for (Int_t i = 0; i<fRows[iRow].NHits(); i++)
 580       AliHLTTPCCADisplay::Instance().DrawHit( iRow, i );
 581   AliHLTTPCCADisplay::Instance().Ask();
 582   AliHLTTPCCADisplay::Instance().ClearView();
 583   AliHLTTPCCADisplay::Instance().DrawSlice( this );
 584   cout<<"draw cells..."<<endl;
 585   for( Int_t iRow=0; iRow<fParam.NRows(); iRow++ )
 586     for (Int_t i = 0; i<fRows[iRow].NCells(); i++)
 587       AliHLTTPCCADisplay::Instance().DrawCell( iRow, i );
 588   AliHLTTPCCADisplay::Instance().Ask();
 589
 590   Int_t nConnectedCells = 0;
 591
 592   cout<<"draw merged cells..."<<endl;
 593
 594   for( Int_t iRow=0; iRow<fParam.NRows(); iRow++ )
 595     for (Int_t i = 0; i<fRows[iRow].NCells(); i++)
 596       {
 597         AliHLTTPCCACell &c = fRows[iRow].Cells()[i];
 598         Int_t id = c.Link();
 599         if( id<0 ) continue;
 600         AliHLTTPCCADisplay::Instance().ConnectCells( iRow,c,ID2IRow(id),ID2Cell(id) );
 601         nConnectedCells++;
 602       }
 603   if( nConnectedCells>0 ){
 604     AliHLTTPCCADisplay::Instance().Ask();
 605   }
 606
 607
 608   AliHLTTPCCADisplay::Instance().ClearView();
 609   AliHLTTPCCADisplay::Instance().DrawSlice( this );
 610   for( Int_t iRow=0; iRow<fParam.NRows(); iRow++ )
 611     for (Int_t i = 0; i<fRows[iRow].NCells(); i++)
 612       AliHLTTPCCADisplay::Instance().DrawCell( iRow, i );
 613
 614   cout<<"draw initial tracks"<<endl;
 615
 616   for( Int_t itr=0; itr<fNTracks; itr++ ){
 617     AliHLTTPCCATrack &iTrack = fTracks[itr];
 618     if( iTrack.NCells()<3 ) continue;
 619     if( !iTrack.Alive() ) continue;
 620     AliHLTTPCCADisplay::Instance().DrawTrack( iTrack );
 621   }
 622   //if( fNTracks>0 )
 623   AliHLTTPCCADisplay::Instance().Ask();
 624 #endif
 625
 626   TStopwatch timer4;
 627
 628   Bool_t doMerging=1;//SG!!!
 629
 630   Float_t factor2 = fParam.TrackConnectionFactor()*fParam.TrackConnectionFactor();
 631   Int_t *refEndPoints = new Int_t[fNHitsTotal];
 632   Int_t nRefEndPoints = 0;
 633   for( Int_t iRow=0; iRow<fParam.NRows(); iRow++ ){
 634     AliHLTTPCCARow &irow = fRows[iRow];
 635     for (Int_t iPoint = 0; iPoint<irow.NEndPoints(); iPoint++){
 636       refEndPoints[nRefEndPoints++] = IRowICell2ID(iRow,iPoint);
 637     }
 638   }
 639   bool first = 1;
 640   while( doMerging ){
 641
 642     //cout<<"do merging"<<endl;
 643
 644     doMerging = 0;
 645
 646     // find nejghbouring tracks
 647     TStopwatch timer5;
 648     //cout<<"nRefEndPoints = "<<nRefEndPoints<<endl;
 649
 650     for( Int_t iRef=0; iRef<nRefEndPoints; iRef++ ){
 651       Int_t iRow = ID2IRow(refEndPoints[iRef]);
 652       AliHLTTPCCARow &irow = fRows[iRow];
 653       Int_t iPoint = ID2ICell(refEndPoints[iRef]);
 654       AliHLTTPCCAEndPoint &ip  = irow.EndPoints()[iPoint];
 655       if( ip.TrackID()<0 ) continue; // not active endpoint
 656
 657       AliHLTTPCCATrack &iTrack = fTracks[ip.TrackID()];
 658       if( iTrack.NCells()<3 ) continue;
 659
 660       Int_t jRowMin = iRow - fParam.MaxTrackMatchDRow();
 661       Int_t jRowMax = iRow + fParam.MaxTrackMatchDRow();
 662       if( jRowMin<0 ) jRowMin = 0;
 663       if( jRowMax>=fParam.NRows() ) jRowMax = fParam.NRows()-1;
 664
 665       if( ip.Param().CosPhi()>=0 ){ //TMath::Abs([2])<TMath::Pi()/2. ){
 666         jRowMin = iRow;
 667       }else{
 668         jRowMax = iRow;
 669       }
 670
 671       Int_t bestNeighbourN = -1;
 672       Float_t bestDist2 = 1.e20;
 673       Int_t bestLink = -1;
 674       if( ip.Link()>=0 ){
 675         bestNeighbourN = fTracks[ID2Point(ip.Link()).TrackID()].NCells();
 676       }
 677
 678       Float_t y0 = ip.Param().GetY();
 679       Float_t z0 = ip.Param().GetZ();
 680
 681  #ifdef DRAW
 682       //AliHLTTPCCADisplay::Instance().DrawTrackletPoint(ip.Param(),kBlue);
 683       //AliHLTTPCCADisplay::Instance().Ask();
 684 #endif
 685
 686       for( Int_t jRow = jRowMin; jRow<=jRowMax; jRow++){
 687         AliHLTTPCCARow &jrow = fRows[jRow];
 688         Float_t dx2 = (irow.X()-jrow.X())*(irow.X()-jrow.X());
 689
 690         // extrapolate the track to row jRow
 691
 692         //for( int ii=0; ii<10; ii++){
 693         //AliHLTTPCCATrackParam iPar = ip.Param();
 694         //iPar.TransportToX( jrow.X());
 695         //}
 696         AliHLTTPCCATrackParam iPar = ip.Param();
 697         Bool_t ok = iPar.TransportToX( jrow.X());
 698         if( !ok ) continue;
 699 #ifdef DRAW
 700         //AliHLTTPCCADisplay::Instance().DrawTrackletPoint(iPar, kBlack);
 701         //AliHLTTPCCADisplay::Instance().Ask();
 702 #endif
 703         Float_t zMin = iPar.GetZ() - 3.5*TMath::Sqrt(iPar.GetErr2Z()) - .2*3.5;
 704         AliHLTTPCCAEndPoint *jjp = lower_bound(jrow.EndPoints(),jrow.EndPoints()+jrow.NEndPoints(),zMin,AliHLTTPCCARow::CompareEndPointZ);
 705
 706         for (Int_t jPoint = jjp-jrow.EndPoints(); jPoint<jrow.NEndPoints(); jPoint++){
 707
 708           AliHLTTPCCAEndPoint &jp  = jrow.EndPoints()[jPoint];
 709
 710           if( jp.TrackID()<0 ) continue; // endpoint not active
 711           if( jp.TrackID()==ip.TrackID() ) continue; // same track
 712
 713           AliHLTTPCCATrack &jTrack = fTracks[jp.TrackID()];
 714
 715           if( bestNeighbourN > jTrack.NCells() ){
 716             continue; // there is already better neighbour found
 717           }
 718           if( jp.Link()>=0 &&
 719               fTracks[ID2Point(jp.Link()).TrackID()].NCells()>=iTrack.NCells() ){
 720             continue; // jTrack is already linked to a better track or to iTrack
 721           }
 722
 723           Float_t dy2, dz2;
 724           AliHLTTPCCATrackParam &jPar = jp.Param();
 725
 726           // check direction
 727           {
 728             if( jPar.GetCosPhi()*iPar.GetCosPhi()>=0 ) continue;
 729           }
 730           // check for neighbouring
 731           {
 732             float d = jPar.GetY() - iPar.GetY();
 733             float s2 = jPar.GetErr2Y() + iPar.GetErr2Y();
 734             if( d*d>factor2*s2 ){
 735               continue;
 736             }
 737             //cout<<"\ndy="<<TMath::Sqrt(d*d)<<", err="<<TMath::Sqrt(factor2*s2)<<endl;
 738             d = jPar.GetZ() - iPar.GetZ();
 739             s2 = jPar.GetErr2Z() + iPar.GetErr2Z();
 740             if( d*d>factor2*s2 ){
 741               if( d>0 ) break;
 742               continue;
 743             }
 744             //cout<<"dz="<<TMath::Sqrt(d*d)<<", err="<<TMath::Sqrt(factor2*s2)<<endl;
 745             if( iTrack.NCells()>=3 && jTrack.NCells()>=3 ){
 746               d = jPar.GetSinPhi() + iPar.GetSinPhi(); //! phi sign is different
 747               s2 = jPar.GetErr2SinPhi() + iPar.GetErr2SinPhi();
 748               if( d*d>factor2*s2 ) continue;
 749               //cout<<"dphi="<<TMath::Sqrt(d*d)<<", err="<<TMath::Sqrt(factor2*s2)<<endl;
 750               d = jPar.GetKappa() + iPar.GetKappa(); // ! kappa sign iz different
 751               s2 = jPar.GetErr2Kappa() + iPar.GetErr2Kappa();
 752               if( d*d>factor2*s2 ) continue;
 753               //cout<<"dk="<<TMath::Sqrt(d*d)<<", err="<<TMath::Sqrt(factor2*s2)<<endl;
 754               d = jPar.GetDzDs() + iPar.GetDzDs(); // ! DzDs sign iz different
 755               s2 = jPar.GetErr2DzDs() + iPar.GetErr2DzDs();
 756               if( d*d>factor2*s2 ) continue;
 757               //cout<<"dlam="<<TMath::Sqrt(d*d)<<", err="<<TMath::Sqrt(factor2*s2)<<endl;
 758             }
 759           }
 760
 761           Float_t dy = jPar.GetY() - y0;
 762           dy2 = dy*dy;
 763           Float_t dz = jPar.GetZ() - z0;
 764           dz2 = dz*dz;
 765           Float_t dist2 = dx2+dy2+dz2;
 766
 767           if( ( bestNeighbourN == jTrack.NCells() ) && dist2>bestDist2 ) continue;
 768
 769           // tracks can be matched
 770
 771           bestLink = IRowICell2ID( jRow, jPoint );
 772           bestNeighbourN = jTrack.NCells();
 773           bestDist2 = dist2;
 774         }
 775       }
 776
 777       if( bestLink < 0 ) continue; // no neighbours found
 778
 779       AliHLTTPCCAEndPoint &jp = ID2Point(bestLink);
 780
 781       if( ip.Link()>=0 ){ // break existing link of iTrack
 782         ID2Point(ip.Link()).Link() = -1;
 783       }
 784       if( jp.Link()>=0 ){ // break existing link of jTrack
 785         ID2Point(jp.Link()).Link() = -1;
 786       }
 787       ip.Link() = bestLink;
 788       jp.Link()= IRowICell2ID( iRow, iPoint );
 789
 790       //cout<<"create link ("<<jp.Link()<<","<<ip.TrackID()<<")->("<<ip.Link()<<","<<jp.TrackID()<<")"<<endl;
 791
 792     }
 793
 794     timer5.Stop();
 795     if(first) fTimers[5] += timer5.CpuTime();
 796
 797     //cout<<"merge neighbours"<<endl;
 798     // merge neighbours
 799
 800     TStopwatch timer6;
 801
 802     Int_t nRefEndPointsNew = 0;
 803     for( Int_t iRef=0; iRef<nRefEndPoints; iRef++ ){
 804
 805       Int_t iRow = ID2IRow(refEndPoints[iRef]);
 806       Int_t iPoint = ID2ICell(refEndPoints[iRef]);
 807       AliHLTTPCCARow &irow = fRows[iRow];
 808       AliHLTTPCCAEndPoint &ip  = irow.EndPoints()[iPoint];
 809       if( ip.TrackID()<0 ) continue; // not active endpoint
 810       if( ip.Link()<0 ) continue; // no neighbours found
 811
 812       Int_t ipID = IRowICell2ID(iRow,iPoint);
 813       Int_t jpID = ip.Link();
 814       AliHLTTPCCAEndPoint &jp  = ID2Point(jpID);
 815
 816       if( jp.Link()!=ipID ){
 817         //cout<<"broken link: jp.Link()!=iID"<<endl;
 818         //exit(0);
 819         return;
 820       }
 821       if( jp.TrackID()<0 ){
 822         //cout<<"broken link: jp.TrackID()<=0"<<endl;
 823         //exit(0);
 824         return;
 825       }
 826       if( jp.TrackID()==ip.TrackID() ){
 827         //cout<<"broken link: jp.TrackID()==ip.TrackID()"<<endl;
 828         //exit(0);
 829         return;
 830       }
 831
 832       //cout<<"Merge neighbours ("<<ipID<<","<<ip.TrackID()<<")->("<<jpID<<","<<jp.TrackID()<<")"<<endl;
 833
 834       AliHLTTPCCATrack &iTrack = fTracks[ip.TrackID()];
 835       AliHLTTPCCATrack &jTrack = fTracks[jp.TrackID()];
 836
 837       // rotate cell link direction for jTrack if necessary
 838
 839       Int_t icID = ip.CellID();
 840       Int_t jcID = jp.CellID();
 841       AliHLTTPCCACell &ic  = ID2Cell(icID);
 842       AliHLTTPCCACell &jc  = ID2Cell(jcID);
 843
 844       if( ic.Link()<0 && jc.Link()<0 || ic.Link()>=0 && jc.Link()>=0 ){
 845
 846         Int_t currID =  jTrack.CellID()[0];
 847         jTrack.CellID()[0] = jTrack.CellID()[2];
 848         jTrack.CellID()[2] = currID;
 849
 850         Int_t pID =  jTrack.PointID()[0];
 851         jTrack.PointID()[0] = jTrack.PointID()[1];
 852         jTrack.PointID()[1] = pID;
 853
 854         currID = jTrack.FirstCellID();
 855         Int_t lastID = -1;
 856         while( currID>=0 ){
 857           AliHLTTPCCACell &c = ID2Cell( currID );
 858           Int_t nextID = c.Link();
 859           c.Link() = lastID;
 860           lastID = currID;
 861           currID = nextID;
 862         }
 863         jTrack.FirstCellID() = lastID;
 864       }
 865       //cout<<"track i "<<ip.TrackID()<<", points="<<ipID<<", "<<iTrack.PointID()[0]<<", "<<iTrack.PointID()[1]<<endl;
 866       //cout<<"track j "<<jp.TrackID()<<", points="<<jpID<<", "<<jTrack.PointID()[0]<<", "<<jTrack.PointID()[1]<<endl;
 867       Int_t itr = ip.TrackID();
 868       ip.TrackID() = -1;
 869       jp.TrackID() = -1;
 870       ip.Link()  = -1;
 871       jp.Link()  = -1;
 872       jTrack.Alive() = 0;
 873
 874       //cout<<"iTrack ID: "<<iTrack.CellID()[0]<<" "<<iTrack.CellID()[1]<<" "<<iTrack.CellID()[2]<<endl;
 875       //cout<<"jTrack ID: "<<jTrack.CellID()[0]<<" "<<jTrack.CellID()[1]<<" "<<jTrack.CellID()[2]<<endl;
 876       if( ic.Link()<0 ){ //match iTrack->jTrack
 877         ic.Link() = jcID;
 878         iTrack.PointID()[1] = jTrack.PointID()[1];
 879         ID2Point(iTrack.PointID()[1]).TrackID() = itr;
 880         if( jTrack.NCells()<3 ){
 881           refEndPoints[nRefEndPointsNew++] = iTrack.PointID()[1];
 882           doMerging = 1;
 883           ID2Point(iTrack.PointID()[1]).Param() = ip.Param();// just to set phi direction
 884         }
 885         if( iTrack.NCells()<3 ){
 886           refEndPoints[nRefEndPointsNew++] = iTrack.PointID()[0];
 887           doMerging = 1;
 888           ID2Point(iTrack.PointID()[0]).Param() = jp.Param();// just to set phi direction
 889         }
 890
 891         if( TMath::Abs(ID2Cell(jTrack.CellID()[2]).Z()-ID2Cell(iTrack.CellID()[0]).Z())>
 892             TMath::Abs(ID2Cell(iTrack.CellID()[2]).Z()-ID2Cell(iTrack.CellID()[0]).Z())  ){
 893           iTrack.CellID()[2] = jTrack.CellID()[2];
 894         }
 895       }else{ //match jTrack->iTrack
 896         jc.Link() = icID;
 897         iTrack.FirstCellID()=jTrack.FirstCellID();
 898         iTrack.PointID()[0] = jTrack.PointID()[0];
 899         ID2Point(iTrack.PointID()[0]).TrackID() = itr;
 900         if( jTrack.NCells()<3 ){
 901           refEndPoints[nRefEndPointsNew++] = iTrack.PointID()[0];
 902           doMerging = 1;
 903           ID2Point(iTrack.PointID()[0]).Param() = ip.Param(); // just to set phi direction
 904         }
 905         if( iTrack.NCells()<3 ){
 906           refEndPoints[nRefEndPointsNew++] = iTrack.PointID()[1];
 907           doMerging = 1;
 908           ID2Point(iTrack.PointID()[1]).Param() = jp.Param();// just to set phi direction
 909         }
 910         if( TMath::Abs(ID2Cell(jTrack.CellID()[0]).Z()-ID2Cell(iTrack.CellID()[2]).Z())>
 911             TMath::Abs(ID2Cell(iTrack.CellID()[0]).Z()-ID2Cell(iTrack.CellID()[2]).Z())  ){
 912           iTrack.CellID()[0] = jTrack.CellID()[0];
 913         }
 914       }
 915
 916       //cout<<"merged ID: "<<iTrack.CellID()[0]<<" "<<iTrack.CellID()[1]<<" "<<iTrack.CellID()[2]<<endl;
 917
 918       if( jTrack.NCells()>iTrack.NCells() ){
 919         iTrack.CellID()[1] = jTrack.CellID()[1];
 920       }
 921
 922       AliHLTTPCCAEndPoint &p0 = ID2Point(iTrack.PointID()[0]);
 923       AliHLTTPCCAEndPoint &p1 = ID2Point(iTrack.PointID()[1]);
 924
 925       if( p0.Link() == iTrack.PointID()[1] ){
 926         p0.Link() = -1;
 927         p1.Link() = -1;
 928       }
 929       //cout<<" NCells itr/jtr= "<<iTrack.NCells()<<" "<<jTrack.NCells()<<endl;
 930       //cout<<" fit merged track "<<itr<<", NCells="<<iTrack.NCells()<<endl;
 931       Float_t *t0 = ( jTrack.NCells()>iTrack.NCells() ) ?jp.Param().Par() :ip.Param().Par();
 932       iTrack.NCells()+=jTrack.NCells();
 933       FitTrack(iTrack,t0);
 934
 935 #ifdef DRAW
 936       cout<<"merged points..."<<ipID<<"/"<<jpID<<endl;
 937       //AliHLTTPCCADisplay::Instance().ConnectCells( iRow,ic,ID2IRow(jcID),jc,kRed );
 938       AliHLTTPCCADisplay::Instance().ConnectEndPoints( ipID,jpID,1.,2,kRed );
 939       AliHLTTPCCADisplay::Instance().DrawEndPoint( ipID,1.,2,kRed );
 940       AliHLTTPCCADisplay::Instance().DrawEndPoint( jpID,1.,2,kRed );
 941       AliHLTTPCCADisplay::Instance().Ask();
 942       cout<<"merged track"<<endl;
 943       AliHLTTPCCADisplay::Instance().DrawTrack(iTrack);
 944       AliHLTTPCCADisplay::Instance().Ask();
 945 #endif
 946       /*
 947       static int ntr=0;
 948       if( ntr++==1 ){
 949         doMerging = 0;
 950         break;
 951       }
 952       */
 953       //doMerging = 1;
 954     }
 955
 956     timer6.Stop();
 957     if(first)fTimers[6] += timer6.CpuTime();
 958
 959     nRefEndPoints = nRefEndPointsNew;
 960
 961     //cout<<"merging ok"<<endl;
 962     //first = 0;
 963   }// do merging
 964
 965   delete[] refEndPoints;
 966   timer4.Stop();
 967   fTimers[4] = timer4.CpuTime();
 968
 969 #ifdef DRAW
 970   //if( fNTracks>0 ) AliHLTTPCCADisplay::Instance().Ask();
 971 #endif
 972
 973
 974
 975
 976 #ifdef DRAWXX
 977   AliHLTTPCCADisplay::Instance().ClearView();
 978   AliHLTTPCCADisplay::Instance().DrawSlice( this );
 979   for( Int_t iRow=0; iRow<fParam.NRows(); iRow++ )
 980     for (Int_t i = 0; i<fRows[iRow].NCells(); i++)
 981       AliHLTTPCCADisplay::Instance().DrawCell( iRow, i );
 982
 983   cout<<"draw final tracks"<<endl;
 984
 985   for( Int_t itr=0; itr<fNTracks; itr++ ){
 986     AliHLTTPCCATrack &iTrack = fTracks[itr];
 987     if( iTrack.NCells()<3 ) continue;
 988     if( !iTrack.Alive() ) continue;
 989     AliHLTTPCCADisplay::Instance().DrawTrack( iTrack );
 990     cout<<"final track "<<itr<<", ncells="<<iTrack.NCells()<<endl;
 991     AliHLTTPCCADisplay::Instance().Ask();
 992   }
 993   AliHLTTPCCADisplay::Instance().Ask();
 994 #endif
 995
 996   // write output
 997   TStopwatch timer7;
 998
 999   //cout<<"write output"<<endl;
1000 #ifdef DRAW
1001   AliHLTTPCCADisplay::Instance().ClearView();
1002   AliHLTTPCCADisplay::Instance().DrawSlice( this );
1003   for( Int_t iRow=0; iRow<fParam.NRows(); iRow++ )
1004     for (Int_t i = 0; i<fRows[iRow].NCells(); i++)
1005       AliHLTTPCCADisplay::Instance().DrawCell( iRow, i );
1006
1007   cout<<"draw out tracks"<<endl;
1008 #endif
1009
1010   fOutTrackHits = new Int_t[fNHitsTotal];
1011   fOutTracks = new AliHLTTPCCAOutTrack[fNTracks];
1012   fNOutTrackHits = 0;
1013   fNOutTracks = 0;
1014
1015   for( Int_t iTr=0; iTr<fNTracks; iTr++){
1016     AliHLTTPCCATrack &iTrack = fTracks[iTr];
1017     if( !iTrack.Alive() ) continue;
1018     if( iTrack.NCells()<3 ) continue;
1019     AliHLTTPCCAOutTrack &out = fOutTracks[fNOutTracks];
1020     out.FirstHitRef() = fNOutTrackHits;
1021     out.NHits() = 0;
1022     out.OrigTrackID() = iTr;
1023     {
1024       AliHLTTPCCAEndPoint &p0 = ID2Point(iTrack.PointID()[0]);
1025       AliHLTTPCCAEndPoint &p2 = ID2Point(iTrack.PointID()[1]);
1026       out.StartPoint() = p0.Param();
1027       out.EndPoint() = p2.Param();
1028     }
1029     AliHLTTPCCATrackParam &t = out.StartPoint();//SG!!!
1030     AliHLTTPCCATrackParam t0 = t;
1031
1032     t.Chi2() = 0;
1033     t.NDF() = -5;
1034     Bool_t first = 1;
1035
1036     Int_t iID = iTrack.FirstCellID();
1037     Int_t fNOutTrackHitsOld = fNOutTrackHits;
1038     for( AliHLTTPCCACell *ic = &ID2Cell(iID); ic->Link()>=0; iID = ic->Link(), ic = &ID2Cell(iID) ){
1039       //cout<<"itr="<<iTr<<", cell ="<<ID2IRow(iID)<<" "<<ID2ICell(iID)<<endl;
1040       AliHLTTPCCARow &row = ID2Row(iID);
1041       if( !t0.TransportToX( row.X() ) ) continue;
1042       Int_t jHit = -1;
1043       Float_t dy, dz, d = 1.e10;
1044       for( Int_t iHit=0; iHit<ic->NHits(); iHit++ ){
1045         AliHLTTPCCAHit &h = row.GetCellHit(*ic,iHit);
1046
1047         // check for wrong hits
1048         {
1049           Float_t ddy = t0.GetY() - h.Y();
1050           Float_t ddz = t0.GetZ() - h.Z();
1051           Float_t dd = ddy*ddy+ddz*ddz;
1052           if( dd<d ){
1053             d = dd;
1054             dy = ddy;
1055             dz = ddz;
1056             jHit = iHit;
1057           }
1058         }
1059       }
1060       if( jHit<0 ) continue;
1061       AliHLTTPCCAHit &h = row.GetCellHit(*ic,jHit);
1062       //if( dy*dy > 3.5*3.5*(/*t0.GetErr2Y() + */h.ErrY()*h.ErrY() ) ) continue;//SG!!!
1063       //if( dz*dz > 3.5*3.5*(/*t0.GetErr2Z() + */h.ErrZ()*h.ErrZ() ) ) continue;
1064       //if( !t0.Filter2( h.Y(), h.Z(), h.ErrY()*h.ErrY(), h.ErrZ()*h.ErrZ() ) ) continue;
1065
1066
1067       if( !t.TransportToX( row.X() ) ) continue;
1068
1069       //* Update the track
1070
1071           if( first ){
1072             t.Cov()[ 0] = .5*.5;
1073             t.Cov()[ 1] = 0;
1074             t.Cov()[ 2] = .5*.5;
1075             t.Cov()[ 3] = 0;
1076             t.Cov()[ 4] = 0;
1077             t.Cov()[ 5] = .2*.2;
1078             t.Cov()[ 6] = 0;
1079             t.Cov()[ 7] = 0;
1080             t.Cov()[ 8] = 0;
1081             t.Cov()[ 9] = .2*.2;
1082             t.Cov()[10] = 0;
1083             t.Cov()[11] = 0;
1084             t.Cov()[12] = 0;
1085             t.Cov()[13] = 0;
1086             t.Cov()[14] = .2*.2;
1087             t.Chi2() = 0;
1088             t.NDF() = -5;
1089           }
1090
1091       if( t.Filter2( h.Y(), h.Z(), h.ErrY()*h.ErrY(), h.ErrZ()*h.ErrZ() ) ) first = 0;
1092       else continue;
1093
1094       fOutTrackHits[fNOutTrackHits] = h.ID();
1095       fNOutTrackHits++;
1096       if( fNOutTrackHits>fNHitsTotal ){
1097         cout<<"fNOutTrackHits>fNHitsTotal"<<endl;
1098         exit(0);//SG!!!
1099         return;
1100       }
1101       out.NHits()++;
1102     }
1103     //cout<<fNOutTracks<<": itr = "<<iTr<<", n outhits = "<<out.NHits()<<endl;
1104     if( out.NHits() > 3 ){
1105       fNOutTracks++;
1106 #ifdef DRAW
1107       AliHLTTPCCADisplay::Instance().DrawTrack( iTrack );
1108       cout<<"out track "<<(fNOutTracks-1)<<", orig = "<<iTr<<", nhits="<<out.NHits()<<endl;
1109       AliHLTTPCCADisplay::Instance().Ask();
1110 #endif
1111
1112     }else {
1113       fNOutTrackHits = fNOutTrackHitsOld;
1114     }
1115   }
1116   //cout<<"end writing"<<endl;
1117   timer7.Stop();
1118   fTimers[7] = timer7.CpuTime();
1119
1120 #ifdef DRAW
1121   AliHLTTPCCADisplay::Instance().Ask();
1122   //AliHLTTPCCADisplay::Instance().DrawMCTracks(fParam.fISec);
1123   //AliHLTTPCCADisplay::Instance().Update();
1124   //AliHLTTPCCADisplay::Instance().Ask();
1125 #endif
1126 }
1127
1128
1129 void AliHLTTPCCATracker::FitTrack( AliHLTTPCCATrack &track, Float_t t0[] ) const
1130 {
1131   //* Fit the track
1132
1133   AliHLTTPCCAEndPoint &p0 = ID2Point(track.PointID()[0]);
1134   AliHLTTPCCAEndPoint &p2 = ID2Point(track.PointID()[1]);
1135   AliHLTTPCCACell &c0 = ID2Cell(p0.CellID());
1136   AliHLTTPCCACell &c1 = ID2Cell(track.CellID()[1]);
1137   AliHLTTPCCACell &c2 = ID2Cell(p2.CellID());
1138   AliHLTTPCCARow &row0 = ID2Row(p0.CellID());
1139   AliHLTTPCCARow &row1 = ID2Row(track.CellID()[1]);
1140   AliHLTTPCCARow &row2 = ID2Row(p2.CellID());
1141
1142
1143   Float_t sp0[5] = {row0.X(), c0.Y(), c0.Z(), c0.ErrY(), c0.ErrZ() };
1144   Float_t sp1[5] = {row1.X(), c1.Y(), c1.Z(), c1.ErrY(), c1.ErrZ() };
1145   Float_t sp2[5] = {row2.X(), c2.Y(), c2.Z(), c2.ErrY(), c2.ErrZ() };
1146   if( track.NCells()>=3 ){
1147     p0.Param().ConstructXYZ3(sp0,sp1,sp2,p0.Param().CosPhi(), t0);
1148     p2.Param().ConstructXYZ3(sp2,sp1,sp0,p2.Param().CosPhi(), t0);
1149     //p2.Param() = p0.Param();
1150     //p2.Param().TransportToX(row2.X());
1151     //p2.Param().Par()[1] = -p2.Param().Par()[1];
1152     //p2.Param().Par()[4] = -p2.Param().Par()[4];
1153   } else {
1154     p0.Param().X() = row0.X();
1155     p0.Param().Y() = c0.Y();
1156     p0.Param().Z() = c0.Z();
1157     p0.Param().Err2Y() = c0.ErrY()*c0.ErrY();
1158     p0.Param().Err2Z() = c0.ErrZ()*c0.ErrZ();
1159     p2.Param().X() = row2.X();
1160     p2.Param().Y() = c2.Y();
1161     p2.Param().Z() = c2.Z();
1162     p2.Param().Err2Y() = c2.ErrY()*c2.ErrY();
1163     p2.Param().Err2Z() = c2.ErrZ()*c2.ErrZ();
1164   }
1165 }