Total amount of used memory has been optimised
[u/mrichter/AliRoot.git] / 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
29 #include "TStopwatch.h"
30 #include "AliHLTTPCCAHitArea.h"
31 #include "AliHLTTPCCANeighboursFinder.h"
32 #include "AliHLTTPCCANeighboursCleaner.h"
33 #include "AliHLTTPCCAStartHitsFinder.h"
34 #include "AliHLTTPCCATrackletConstructor.h"
35 #include "AliHLTTPCCATrackletSelector.h"
36 #include "AliHLTTPCCAProcess.h"
37 #include "AliHLTTPCCAUsedHitsInitialiser.h"
38
39 #include "AliHLTTPCCATrackParam.h"
40
41 #if !defined(HLTCA_GPUCODE)
42 #include <iostream>
43 #endif
44
45 //#define DRAW
46
47 #ifdef DRAW
48   #include "AliHLTTPCCADisplay.h"
49 #endif //DRAW
50
51 #ifdef HLTCA_INTERNAL_PERFORMANCE
52  #include "AliHLTTPCCAPerformance.h"
53 #endif
54
55
56 ClassImp(AliHLTTPCCATracker)
57
58 #if !defined(HLTCA_GPUCODE)  
59
60 AliHLTTPCCATracker::AliHLTTPCCATracker()
61   :
62   fParam(),
63   fNHitsTotal(0),  
64   fCommonMemory(0), 
65   fCommonMemorySize(0),
66   fHitMemory(0), 
67   fHitMemorySize(0),
68   fTrackMemory(0), 
69   fTrackMemorySize(0),
70   fInputEvent(0),     
71   fInputEventSize(0), 
72   fRowData(0),     
73   fRowDataSize(0), 
74   fHitInputIDs(0), 
75   fHitWeights(0),  
76   fNTracklets(0),
77   fTrackletStartHits(0),
78   fTracklets(0), 
79   fNTracks(0),   
80   fTracks(0), 
81   fNTrackHits(0),
82   fTrackHits(0),
83   fNOutTracks(0),
84   fOutTracks(0), 
85   fNOutTrackHits(0),
86   fOutTrackHits(0),
87   fTmpHitInputIDs(0)
88 {
89   // constructor
90 }
91
92 AliHLTTPCCATracker::AliHLTTPCCATracker( const AliHLTTPCCATracker& )
93   :
94   fParam(),
95   fNHitsTotal(0),
96   fCommonMemory(0), 
97   fCommonMemorySize(0),
98   fHitMemory(0), 
99   fHitMemorySize(0),
100   fTrackMemory(0), 
101   fTrackMemorySize(0),
102   fInputEvent(0),     
103   fInputEventSize(0), 
104   fRowData(0),     
105   fRowDataSize(0), 
106   fHitInputIDs(0), 
107   fHitWeights(0),  
108   fNTracklets(0),
109   fTrackletStartHits(0),
110   fTracklets(0), 
111   fNTracks(0),   
112   fTracks(0), 
113   fNTrackHits(0),
114   fTrackHits(0),
115   fNOutTracks(0),
116   fOutTracks(0), 
117   fNOutTrackHits(0),
118   fOutTrackHits(0),
119   fTmpHitInputIDs(0)
120 {
121   // dummy
122 }
123
124 AliHLTTPCCATracker &AliHLTTPCCATracker::operator=( const AliHLTTPCCATracker& )
125 {
126   // dummy
127   fCommonMemory = 0;
128   fHitMemory = 0;
129   fTrackMemory = 0;
130   return *this;
131 }
132
133 GPUd() AliHLTTPCCATracker::~AliHLTTPCCATracker()
134 {
135   // destructor
136   if( fCommonMemory ) delete[] fCommonMemory;
137   if( fHitMemory ) delete[] fHitMemory;
138   if( fTrackMemory ) delete[] fTrackMemory;
139   if( fTmpHitInputIDs ) delete[] fTmpHitInputIDs;
140 }
141 #endif
142
143
144
145 // ----------------------------------------------------------------------------------
146 GPUd() void AliHLTTPCCATracker::Initialize( const AliHLTTPCCAParam &param )
147 {
148   // initialisation
149   fParam = param;
150   fParam.Update();
151   for( Int_t irow=0; irow<fParam.NRows(); irow++ ){
152     fRows[irow].SetX( fParam.RowX(irow) );
153     fRows[irow].SetMaxY( CAMath::Tan( fParam.DAlpha()/2.)*fRows[irow].X() );
154   }
155   StartEvent();
156 }
157
158 GPUd() void AliHLTTPCCATracker::StartEvent()
159 {
160   // start new event and fresh the memory  
161
162   if( !fCommonMemory ){
163     SetPointersCommon(); // just to calculate the size
164     fCommonMemory = reinterpret_cast<Char_t*> ( new uint4 [ fCommonMemorySize/sizeof(uint4) + 100] );
165     SetPointersCommon();// set pointers
166   }
167   
168   if( fHitMemory ) delete[] fHitMemory;
169   fHitMemory = 0;
170   if( fTrackMemory ) delete[] fTrackMemory;
171   fTrackMemory = 0;
172
173   fNHitsTotal = 0;
174   *fNTracklets = 0;
175   *fNTracks = 0 ;
176   *fNTrackHits = 0;
177   *fNOutTracks = 0;
178   *fNOutTrackHits = 0;
179   if( fTmpHitInputIDs ) delete[] fTmpHitInputIDs;
180   fTmpHitInputIDs = 0;
181 }
182
183
184 /*
185 GPUhd() void  AliHLTTPCCATracker::SetPointers()
186 {
187   // set all pointers to the event memory
188
189   Int_t gridSizeTotal = 2*(2*fNHitsTotal + 10*Param().NRows());
190   gridSizeTotal *=100;//SG!!!
191   ULong_t mem = (ULong_t) fCommonMemory;  
192   UInt_t sI = sizeof(Int_t);
193   UInt_t sF = sizeof(Float_t);
194   UInt_t sS = sizeof(Short_t);
195   UInt_t s4 = sizeof(uint4);
196
197   fInputEvent = (Char_t*) mem;  
198   fInputEventSize = (1+fParam.NRows()*2 + 1)*sI + (fNHitsTotal*2)*sF;
199   mem+= fInputEventSize;
200
201   mem = ( mem/s4 + 1 )*s4;
202   fRowData = (uint4*) mem;
203   fRowDataSize = ( 2*fNHitsTotal*sS +  //  yz
204                    gridSizeTotal*sS + // grid
205                    2*fNHitsTotal*sS +  // link up,link down
206                    fParam.NRows()*s4   // row alignment
207                    );
208   mem += fRowDataSize;
209
210   mem = ( mem/sI + 1 )*sI;
211
212   fHitInputIDs = (Int_t*) mem;
213   mem+= fNHitsTotal*sI;
214
215   fHitWeights = (Int_t*) mem;
216   mem+=  fNHitsTotal*sI;
217
218   fNTracklets = (Int_t*) mem;
219   mem+= sI;
220
221   fTrackletStartHits = (Int_t*) mem;
222   mem+= fNHitsTotal*sI;
223
224   mem = ( mem/sizeof(AliHLTTPCCATracklet) + 1 )*sizeof(AliHLTTPCCATracklet);
225
226   fTracklets = (AliHLTTPCCATracklet *) mem;
227   mem+= fNHitsTotal*sizeof(AliHLTTPCCATracklet);
228   
229   mem = ( mem/sI + 1 )*sI;
230
231   fNTracks = (Int_t*) mem;
232   mem+= sI;
233
234   mem = ( mem/sizeof(AliHLTTPCCATrack) + 1 )*sizeof(AliHLTTPCCATrack);
235
236   fTracks = (AliHLTTPCCATrack*) mem;
237   mem+= fNHitsTotal*sizeof(AliHLTTPCCATrack);
238   
239   mem = ( mem/sI + 1 )*sI;
240   fNTrackHits = (Int_t*) mem;
241   mem+= sI;
242
243   fTrackHits = (Int_t*) mem;
244   mem+= 100*fNHitsTotal*sI;//SG!!!
245
246   fNOutTracks = (Int_t*) mem;
247   mem+= sI;
248
249   mem = ( mem/sizeof(AliHLTTPCCAOutTrack) + 1 )*sizeof(AliHLTTPCCAOutTrack);
250   
251   fOutTracks = (AliHLTTPCCAOutTrack*) mem;
252   mem+= fNHitsTotal*sizeof(AliHLTTPCCAOutTrack);
253
254   mem = ( mem/sI + 1 )*sI;
255   fNOutTrackHits = (Int_t*) mem;
256   mem+= sI;
257   fOutTrackHits = (Int_t*) mem;
258   mem+= 100*fNHitsTotal*sI; //SG!!!
259
260   fCommonMemorySize = mem - (ULong_t) fCommonMemory;
261 }
262 */
263
264
265 GPUhd() void  AliHLTTPCCATracker::SetPointersCommon()
266 {
267   // set all pointers to the event memory
268
269   ULong_t mem = (ULong_t) fCommonMemory;  
270   UInt_t sI = sizeof(Int_t);
271   
272   // set common memory
273
274   fNTracklets = (Int_t*) mem;
275   mem+= sI;
276   fNTracks = (Int_t*) mem;
277   mem+= sI;
278   fNTrackHits = (Int_t*) mem;
279   mem+= sI;  
280   fNOutTracks = (Int_t*) mem;
281   mem+= sI;
282   fNOutTrackHits = (Int_t*) mem;
283   mem+= sI;
284
285   // calculate the size
286
287   fCommonMemorySize = mem - (ULong_t) fCommonMemory;
288 }
289
290
291 GPUhd() void  AliHLTTPCCATracker::SetPointersHits( Int_t MaxNHits )
292 {
293   // set all pointers to the event memory
294
295   Int_t gridSizeTotal = 2*(2*MaxNHits + 10*Param().NRows());
296   //gridSizeTotal *=100;//SG!!!
297
298   ULong_t mem = (ULong_t) fHitMemory;  
299   UInt_t sI = sizeof(Int_t);
300   UInt_t sF = sizeof(Float_t);
301   UInt_t sS = sizeof(Short_t);
302   UInt_t s4 = sizeof(uint4);
303    
304   // set input event
305
306   mem = ( mem/s4 + 1 )*s4;
307   fInputEvent = (Char_t*) mem;  
308   fInputEventSize = (1+fParam.NRows()*2 + 1)*sI + (MaxNHits*2)*sF;
309   mem+= fInputEventSize;
310
311   // set cluster data for TPC rows
312
313   mem = ( mem/s4 + 1 )*s4;
314   fRowData = (uint4*) mem;
315   fRowDataSize = ( 2*MaxNHits*sS +  //  yz
316                    gridSizeTotal*sS + // grid
317                    2*MaxNHits*sS +  // link up,link down
318                    fParam.NRows()*s4   // row alignment
319                    );
320   mem += fRowDataSize;
321
322   // extra arrays for tpc clusters
323   
324   mem = ( mem/sI + 1 )*sI;
325
326   fHitInputIDs = (Int_t*) mem;
327   mem+= MaxNHits*sI;
328
329   fTrackletStartHits = (Int_t*) mem;
330   mem+= MaxNHits*sI;
331
332   fHitWeights = (Int_t*) mem;
333   mem+=  MaxNHits*sI;
334
335   // arrays for track hits
336
337   fTrackHits = (Int_t*) mem;
338   mem+= 10*MaxNHits*sI;//SG!!!
339   
340   fOutTrackHits = (Int_t*) mem;
341   mem+= 10*MaxNHits*sI; //SG!!!
342
343   // calculate the size
344
345   fHitMemorySize = mem - (ULong_t) fHitMemory;
346 }
347
348
349 GPUhd() void  AliHLTTPCCATracker::SetPointersTracks( Int_t MaxNTracks )
350 {
351   // set all pointers to the tracks memory
352
353   ULong_t mem = (ULong_t) fTrackMemory;  
354   
355   // memory for tracklets
356
357   mem = ( mem/sizeof(AliHLTTPCCATracklet) + 1 )*sizeof(AliHLTTPCCATracklet);
358   fTracklets = (AliHLTTPCCATracklet *) mem;
359   mem+= MaxNTracks*sizeof(AliHLTTPCCATracklet);  
360   
361   // memory for selected tracks
362   
363   mem = ( mem/sizeof(AliHLTTPCCATrack) + 1 )*sizeof(AliHLTTPCCATrack);
364   fTracks = (AliHLTTPCCATrack*) mem;
365   mem+= MaxNTracks*sizeof(AliHLTTPCCATrack);  
366   
367   // memory for output tracks
368
369   mem = ( mem/sizeof(AliHLTTPCCAOutTrack) + 1 )*sizeof(AliHLTTPCCAOutTrack);
370   
371   fOutTracks = (AliHLTTPCCAOutTrack*) mem;
372   mem+= fNHitsTotal*sizeof(AliHLTTPCCAOutTrack);  
373
374   // calculate the size
375
376   fTrackMemorySize = mem - (ULong_t) fTrackMemory;
377 }
378
379
380
381 GPUd() void AliHLTTPCCATracker::ReadEvent( const Int_t *RowFirstHit, const Int_t *RowNHits, const Float_t *Y, const Float_t *Z, Int_t NHits )
382 {
383   //* Read event
384
385   StartEvent();
386   
387   fNHitsTotal = NHits;
388   
389   {
390     SetPointersHits(NHits); // to calculate the size
391     fHitMemory = reinterpret_cast<Char_t*> ( new uint4 [ fHitMemorySize/sizeof(uint4) + 100] );   
392     SetPointersHits(NHits); // set pointers for hits
393     *fNTracklets = 0;
394     *fNTracks = 0 ;
395     *fNOutTracks = 0;
396     *fNOutTrackHits = 0;
397   }
398
399   reinterpret_cast<Int_t*>( fInputEvent )[0] = fParam.NRows();
400   reinterpret_cast<Int_t*>( fInputEvent )[1+fParam.NRows()*2] = NHits;
401   Int_t *rowHeaders = reinterpret_cast<Int_t*>( fInputEvent ) +1;
402   Float_t *hitsYZ = reinterpret_cast<Float_t*>( fInputEvent ) + 1+fParam.NRows()*2+1;
403   for( Int_t iRow=0; iRow<fParam.NRows(); iRow++ ){
404     rowHeaders[iRow*2  ] = RowFirstHit[iRow];
405     rowHeaders[iRow*2+1] = RowNHits[iRow];
406   }
407   for( Int_t iHit=0; iHit<NHits; iHit++ ){
408     hitsYZ[iHit*2  ] = Y[iHit];
409     hitsYZ[iHit*2+1] = Z[iHit];
410   }
411
412   //SG cell finder - test code
413
414   if( fTmpHitInputIDs ) delete[] fTmpHitInputIDs;
415   fTmpHitInputIDs = new Int_t [NHits];
416   const Float_t areaY = .5;
417   const Float_t areaZ = .5;
418   Int_t newRowNHitsTotal = 0;
419   Bool_t *usedHits = new Bool_t [NHits];
420   for( Int_t iHit=0; iHit<NHits; iHit++ ) usedHits[iHit] = 0;
421   for( Int_t iRow=0; iRow<fParam.NRows(); iRow++ ){
422     rowHeaders[iRow*2  ] = newRowNHitsTotal; // new first hit
423     rowHeaders[iRow*2+1] = 0; // new N hits
424     Int_t newRowNHits = 0;
425     Int_t oldRowFirstHit = RowFirstHit[iRow];
426     Int_t oldRowLastHit = oldRowFirstHit + RowNHits[iRow];
427     for( Int_t iHit=oldRowFirstHit; iHit<oldRowLastHit; iHit++ ){
428       if( usedHits[iHit] ) continue;
429       Float_t y0 = Y[iHit];
430       Float_t z0 = Z[iHit];
431       Float_t cy = y0;
432       Float_t cz = z0;
433       Int_t nclu = 1;
434       usedHits[iHit]=1;
435       if(0) for( Int_t jHit=iHit+1; jHit<oldRowLastHit; jHit++ ){//SG!!!
436         //if( usedHits[jHit] ) continue;
437         Float_t dy = Y[jHit] - y0;
438         Float_t dz = Z[jHit] - z0;
439         if( CAMath::Abs(dy)<areaY && CAMath::Abs(dz)<areaZ ){
440           cy+=Y[jHit];
441           cz+=Z[jHit];
442           nclu++;
443           usedHits[jHit]=1;
444         }
445       }
446       Int_t id = newRowNHitsTotal+newRowNHits;
447       hitsYZ[id*2  ] = cy/nclu;
448       hitsYZ[id*2+1] = cz/nclu;
449       fTmpHitInputIDs[id] = iHit;
450       newRowNHits++;
451     }
452     rowHeaders[iRow*2+1] = newRowNHits;
453     newRowNHitsTotal+=newRowNHits;
454   }
455   fNHitsTotal = newRowNHitsTotal;
456   reinterpret_cast<Int_t*>( fInputEvent )[1+fParam.NRows()*2] = newRowNHitsTotal;
457
458   delete[] usedHits;  
459   SetupRowData();
460 }
461
462
463 GPUd() void AliHLTTPCCATracker::SetupRowData()
464 {
465   //* Convert input hits, create grids, etc.
466
467   fNHitsTotal = reinterpret_cast<Int_t*>( fInputEvent )[1+fParam.NRows()*2];
468   Int_t *rowHeaders = reinterpret_cast<Int_t*>( fInputEvent ) +1;
469   Float_t *hitsYZ = reinterpret_cast<Float_t*>( fInputEvent ) + 1+fParam.NRows()*2+1;
470
471   for( Int_t iRow=0; iRow<fParam.NRows(); iRow++ ){
472     
473     AliHLTTPCCARow &row = fRows[iRow];
474     row.SetFirstHit( rowHeaders[iRow*2] );
475     row.SetNHits( rowHeaders[iRow*2+1] );
476     Float_t yMin=1.e3, yMax=-1.e3, zMin=1.e3, zMax=-1.e3;
477     Int_t nGrid =  row.NHits();
478     for( Int_t i=0; i<row.NHits(); i++ ){
479       Int_t j = row.FirstHit()+i;
480       Float_t y = hitsYZ[j*2];
481       Float_t z = hitsYZ[j*2+1];
482       if( yMax < y ) yMax = y;
483       if( yMin > y ) yMin = y;
484       if( zMax < z ) zMax = z;
485       if( zMin > z ) zMin = z;
486     }
487     if( nGrid <= 0 ){
488       yMin = yMax = zMin = zMax = 0;
489       nGrid = 1;
490     }
491
492     AliHLTTPCCAGrid grid;
493     grid.Create( yMin, yMax, zMin, zMax, nGrid );
494     
495     Float_t sy = ( CAMath::Abs( grid.StepYInv() ) >1.e-4 ) ?1./grid.StepYInv() :1;
496     Float_t sz = ( CAMath::Abs( grid.StepZInv() ) >1.e-4 ) ?1./grid.StepZInv() :1;
497     
498     //cout<<"grid n = "<<row.Grid().N()<<" "<<sy<<" "<<sz<<" "<<yMin<<" "<<yMax<<" "<<zMin<<" "<<zMax<<endl;
499     
500     Bool_t recreate=0;
501     if( sy < 2. ) { recreate = 1; sy = 2; }
502     if( sz < 2. ) { recreate = 1; sz = 2; }
503     //recreate = 1;//SG!!!
504     //sy=2;
505     //sz=2;
506     if( recreate ) grid.Create( yMin, yMax, zMin, zMax, sy, sz );
507     row.SetGrid( grid );
508   }
509
510   AliHLTTPCCAHit ffHits[fNHitsTotal];
511
512   Int_t rowDataOffset = 0;
513
514   for( Int_t iRow=0; iRow<fParam.NRows(); iRow++ ){
515
516     AliHLTTPCCARow &row = fRows[iRow];
517     const AliHLTTPCCAGrid &grid = row.Grid();
518
519     Int_t c[grid.N()+3+10];
520     Int_t bins[row.NHits()];
521     Int_t filled[ row.Grid().N() +3+10 ];
522
523     for( UInt_t bin=0; bin<row.Grid().N()+3; bin++ ) filled[bin] = 0;  
524
525     for( Int_t i=0; i<row.NHits(); i++ ){
526       Int_t j = row.FirstHit()+i;
527       Int_t bin = row.Grid().GetBin( hitsYZ[2*j], hitsYZ[2*j+1] );
528       bins[i] = bin;
529       filled[bin]++;
530     }
531
532     {
533       Int_t n=0;
534       for( UInt_t bin=0; bin<row.Grid().N()+3; bin++ ){
535         c[bin] = n;
536         n+=filled[bin];
537       }
538     }
539
540     for( Int_t i=0; i<row.NHits(); i++ ){ 
541       Int_t bin = bins[i];
542       Int_t ind = c[bin] + filled[bin]-1;
543
544      AliHLTTPCCAHit &h = ffHits[row.FirstHit()+ind];
545       fHitInputIDs[row.FirstHit()+ind] = fTmpHitInputIDs[row.FirstHit()+i];
546       h.SetY( hitsYZ[2*(row.FirstHit()+i)] );
547       h.SetZ( hitsYZ[2*(row.FirstHit()+i)+1] );
548       filled[bin]--;
549     }
550
551     {
552       Float_t y0 = row.Grid().YMin();
553       Float_t stepY = (row.Grid().YMax() - y0)*(1./65535.);
554       Float_t z0 = row.Grid().ZMin();
555       Float_t stepZ = (row.Grid().ZMax() - z0)*(1./65535.);
556       if( stepY<1.e-4 ) stepY = 1.e-4;
557       if( stepZ<1.e-4 ) stepZ = 1.e-4;
558       Float_t stepYi = 1./stepY;
559       Float_t stepZi = 1./stepZ;
560       
561       row.SetHy0( y0 );
562       row.SetHz0( z0 );
563       row.SetHstepY( stepY );
564       row.SetHstepZ( stepZ );
565       row.SetHstepYi( stepYi );
566       row.SetHstepZi( stepZi );
567       
568       row.SetFullOffset( rowDataOffset );
569       ushort2 *p= (ushort2*)( fRowData + row.FullOffset() );
570       for( Int_t ih=0; ih<row.NHits(); ih++ ){
571         Int_t ihTot = row.FirstHit()+ih;
572         AliHLTTPCCAHit &hh = ffHits[ihTot];
573         Float_t xx = ((hh.Y() - y0)*stepYi);
574         Float_t yy = ((hh.Z() - z0)*stepZi);
575         if( xx<0 || yy<0 || xx>=65536 || yy>= 65536 ){
576           std::cout<<"!!!! hit packing error!!! "<<xx<<" "<<yy<<" "<<std::endl;
577         }
578         p[ih].x = (UShort_t) xx;
579         p[ih].y = (UShort_t) yy;
580       }
581       Int_t size = row.NHits()*sizeof(ushort2);
582
583       row.SetFullGridOffset( row.NHits()*2 );   
584       UShort_t *p1 = ((UShort_t *)p) + row.FullGridOffset();
585
586       Int_t n = grid.N();
587       for( Int_t i=0; i<n; i++ ){
588         p1[i] = c[i];
589       }     
590       UShort_t a = c[n];
591       Int_t nn = n+grid.Ny()+3;
592       for( Int_t i=n; i<nn; i++ ) p1[i] = a;
593
594       size+= (nn)*sizeof(UShort_t);
595       row.SetFullLinkOffset( row.NHits()*2 + nn );
596       size+= row.NHits()*2*sizeof(Short_t);
597
598       if( size%16 ) size = size/sizeof(uint4)+1;
599       else size = size/sizeof(uint4);
600       row.SetFullSize( size );
601       //cout<<iRow<<", "<<row.fNHits<<"= "<<size*16<<"b: "<<row.fFullOffset<<" "<<row.fFullSize<<" "<<row.fFullGridOffset<<" "<<row.fFullLinkOffset<<std::endl;
602
603       rowDataOffset+=size;
604     }  
605   } 
606 }
607
608
609 GPUh() void AliHLTTPCCATracker::Reconstruct()
610 {
611   //* reconstruction of event
612   //std::cout<<"Reconstruct slice "<<fParam.ISlice()<<", nHits="<<fNHitsTotal<<std::endl;
613
614   fTimers[0] = 0; // find neighbours
615   fTimers[1] = 0; // construct tracklets
616   fTimers[2] = 0; // fit tracklets
617   fTimers[3] = 0; // prolongation of tracklets
618   fTimers[4] = 0; // selection
619   fTimers[5] = 0; // write output
620   fTimers[6] = 0;
621   fTimers[7] = 0;
622   
623   //if( fParam.ISlice()<1 ) return; //SG!!!
624
625   TStopwatch timer0;
626
627   //SetupRowData();
628   if( fNHitsTotal < 1 ) return;
629
630 #ifdef DRAW
631
632   AliHLTTPCCADisplay::Instance().ClearView();  
633   AliHLTTPCCADisplay::Instance().SetSliceView();
634   AliHLTTPCCADisplay::Instance().SetCurrentSlice( this );
635   AliHLTTPCCADisplay::Instance().DrawSlice( this, 1 );  
636   if( fNHitsTotal>0 ){
637     AliHLTTPCCADisplay::Instance().DrawSliceHits( kRed, 1.);  
638     AliHLTTPCCADisplay::Instance().Ask();
639   }
640 #endif
641
642   *fNTracks = 0;
643   *fNTracklets = 0;
644
645 #if !defined(HLTCA_GPUCODE)  
646  
647   AliHLTTPCCAProcess<AliHLTTPCCANeighboursFinder>( Param().NRows(), 1, *this );
648
649 #ifdef HLTCA_INTERNAL_PERFORMANCE 
650   //if( Param().ISlice()<=2 ) 
651   //AliHLTTPCCAPerformance::Instance().LinkPerformance( Param().ISlice() );
652 #endif
653
654
655 #ifdef DRAW
656   if( fNHitsTotal>0 ){
657     AliHLTTPCCADisplay::Instance().DrawSliceLinks( -1, -1, 1);  
658     AliHLTTPCCADisplay::Instance().Ask();
659   }
660 #endif
661
662
663   AliHLTTPCCAProcess<AliHLTTPCCANeighboursCleaner>( Param().NRows()-2, 1, *this );
664   AliHLTTPCCAProcess<AliHLTTPCCAStartHitsFinder>( Param().NRows()-4, 1, *this );
665
666   Int_t nStartHits = *fNTracklets;
667   
668   Int_t nThreads = 128;
669   Int_t nBlocks = fNHitsTotal/nThreads + 1;
670   if( nBlocks<12 ){
671     nBlocks = 12; 
672     nThreads = fNHitsTotal/12+1;
673     if( nThreads%32 ) nThreads = (nThreads/32+1)*32;
674   }
675       
676   nThreads = fNHitsTotal;
677   nBlocks = 1;
678
679   AliHLTTPCCAProcess<AliHLTTPCCAUsedHitsInitialiser>(nBlocks, nThreads,*this);
680
681
682   {
683     SetPointersTracks(nStartHits); // to calculate the size
684     fTrackMemory = reinterpret_cast<Char_t*> ( new uint4 [ fTrackMemorySize/sizeof(uint4) + 100] );   
685     SetPointersTracks(nStartHits); // set pointers for hits
686   }
687
688   Int_t nMemThreads = AliHLTTPCCATrackletConstructor::NMemThreads();
689   nThreads = 256;//96;
690   nBlocks = nStartHits/nThreads + 1;
691   if( nBlocks<30 ){
692     nBlocks = 30;
693     nThreads = (nStartHits)/30+1;
694     if( nThreads%32 ) nThreads = (nThreads/32+1)*32;
695   }
696
697   nThreads = nStartHits;
698   nBlocks = 1;
699
700   AliHLTTPCCAProcess1<AliHLTTPCCATrackletConstructor>(nBlocks, nMemThreads+nThreads,*this);
701
702   //std::cout<<"Slice "<<Param().ISlice()<<": NHits="<<fNHitsTotal<<", NTracklets="<<*NTracklets()<<std::endl;
703
704   { 
705     nThreads = 128;
706     nBlocks = nStartHits/nThreads + 1;
707     if( nBlocks<12 ){
708       nBlocks = 12;  
709       nThreads = nStartHits/12+1;
710       nThreads = (nThreads/32+1)*32;
711     }
712         
713     *fNTrackHits = 0;
714
715     nThreads = nStartHits;
716     nBlocks = 1;
717
718
719     AliHLTTPCCAProcess<AliHLTTPCCATrackletSelector>(nBlocks, nThreads,*this);
720
721     //std::cout<<"Slice "<<Param().ISlice()<<": N start hits/tracklets/tracks = "<<nStartHits<<" "<<nStartHits<<" "<<*fNTracks<<std::endl;
722    }
723
724   //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;
725
726   
727   WriteOutput();      
728   
729
730 #endif
731
732 #ifdef DRAW
733  {
734    AliHLTTPCCADisplay &disp = AliHLTTPCCADisplay::Instance();
735    AliHLTTPCCATracker &slice = *this;
736    std::cout<<"N out tracks = "<<*slice.NOutTracks()<<std::endl;
737    //disp.Ask();
738    AliHLTTPCCADisplay::Instance().SetCurrentSlice( this );
739    AliHLTTPCCADisplay::Instance().DrawSlice( this, 1 );
740    disp.DrawSliceHits(-1,.5);
741    for( Int_t itr=0; itr<*slice.NOutTracks(); itr++ ){
742      std::cout<<"track N "<<itr<<", nhits="<<slice.OutTracks()[itr].NHits()<<std::endl;
743      disp.DrawSliceOutTrack( itr, kBlue );
744      //disp.Ask();
745      //int id = slice.OutTracks()[itr].OrigTrackID();
746      //AliHLTTPCCATrack &tr = Tracks()[id];
747      //for( Int_t ih=0; ih<tr.NHits(); ih++ ){       
748      //Int_t ic = (fTrackHits[tr.FirstHitID()+ih]);
749      //std::cout<<ih<<" "<<ID2IRow(ic)<<" "<<ID2IHit(ic)<<std::endl;
750      //}
751      //disp.DrawSliceTrack( id, kBlue );
752      //disp.Ask();
753    }
754    disp.Ask();
755  }
756 #endif
757
758   timer0.Stop();  
759   fTimers[0] = timer0.CpuTime()/100.;
760
761  }
762
763
764
765
766 GPUh() void AliHLTTPCCATracker::WriteOutput()
767 {
768   // write output
769
770   TStopwatch timer;
771  
772   *fNOutTrackHits = 0;
773   *fNOutTracks = 0;
774
775   //cout<<"output: nTracks = "<<*fNTracks<<", nHitsTotal="<<fNHitsTotal<<std::endl;  
776
777   for( Int_t iTr=0; iTr<*fNTracks; iTr++){
778
779     AliHLTTPCCATrack &iTrack = fTracks[iTr];
780
781     //cout<<"iTr = "<<iTr<<", nHits="<<iTrack.NHits()<<std::endl;
782
783     //if( !iTrack.Alive() ) continue;
784     if( iTrack.NHits()<3 ) continue;          
785     AliHLTTPCCAOutTrack &out = fOutTracks[*fNOutTracks];
786     out.SetFirstHitRef( *fNOutTrackHits );
787     out.SetNHits( 0 );
788     out.SetOrigTrackID( iTr );    
789     out.SetStartPoint( iTrack.Param() );
790     out.SetEndPoint( iTrack.Param() );    
791
792     Int_t iID = iTrack.FirstHitID();
793     Int_t nOutTrackHitsOld = *fNOutTrackHits;
794  
795     for( Int_t ith=0; ith<iTrack.NHits(); ith++ ){
796       Int_t ic = (fTrackHits[iID+ith]);
797       const AliHLTTPCCARow &row = ID2Row(ic);
798       Int_t ih = ID2IHit(ic);
799       fOutTrackHits[*fNOutTrackHits] = fHitInputIDs[row.FirstHit()+ih];      
800       (*fNOutTrackHits)++;
801       //cout<<"write i,row,hit,id="<<ith<<", "<<ID2IRow(ic)<<", "<<ih<<", "<<fHitInputIDs[row.FirstHit()+ih]<<std::endl;     
802       if( *fNOutTrackHits>=10*fNHitsTotal ){
803         std::cout<<"fNOutTrackHits>fNHitsTotal"<<std::endl;
804         //exit(0);
805         return;//SG!!!
806       }
807       out.SetNHits( out.NHits() + 1 );      
808     }    
809     if( out.NHits() >= 2 ){
810       (*fNOutTracks)++;
811     }else {
812       (*fNOutTrackHits) = nOutTrackHitsOld;
813     }
814   }
815   timer.Stop();
816   fTimers[5]+=timer.CpuTime();
817 }
818
819 GPUh() void AliHLTTPCCATracker::FitTrackFull( AliHLTTPCCATrack &/**/, Float_t * /**/ ) const 
820 {  
821   // fit track with material
822 #ifdef XXX    
823   //* Fit the track   
824   FitTrack( iTrack, tt0 );
825   if( iTrack.NHits()<=3 ) return;
826     
827   AliHLTTPCCATrackParam &t = iTrack.Param();
828   AliHLTTPCCATrackParam t0 = t;
829
830   t.Chi2() = 0;
831   t.NDF() = -5; 
832   Bool_t first = 1;
833
834   Int_t iID = iTrack.FirstHitID();
835   for( Int_t ih=0; ih<iTrack.NHits(); ih++, iID++ ){
836     Int_t *ic = &(fTrackHits[iID]);
837     Int_t iRow = ID2IRow(*ic);
838     AliHLTTPCCARow &row = fRows[iRow];      
839     if( !t0.TransportToX( row.X() ) ) continue;      
840     Float_t dy, dz;
841     AliHLTTPCCAHit &h = ID2Hit(*ic);
842
843     // check for wrong hits     
844     if(0){
845       dy = t0.GetY() - h.Y();
846       dz = t0.GetZ() - h.Z();
847       
848       //if( dy*dy > 3.5*3.5*(/*t0.GetErr2Y() + */h.ErrY()*h.ErrY() ) ) continue;//SG!!!
849       //if( dz*dz > 3.5*3.5*(/*t0.GetErr2Z() + */h.ErrZ()*h.ErrZ() ) ) continue;      
850     }
851
852     if( !t.TransportToX( row.X() ) ) continue;  
853
854     //* Update the track
855             
856     if( first ){
857       t.Cov()[ 0] = .5*.5;
858       t.Cov()[ 1] = 0;
859       t.Cov()[ 2] = .5*.5;
860       t.Cov()[ 3] = 0;
861       t.Cov()[ 4] = 0;
862       t.Cov()[ 5] = .2*.2;
863       t.Cov()[ 6] = 0;
864       t.Cov()[ 7] = 0;
865       t.Cov()[ 8] = 0;
866       t.Cov()[ 9] = .2*.2;
867       t.Cov()[10] = 0;
868       t.Cov()[11] = 0;
869       t.Cov()[12] = 0;
870       t.Cov()[13] = 0;
871       t.Cov()[14] = .2*.2;
872       t.Chi2() = 0;
873       t.NDF() = -5;
874     }
875     Float_t err2Y, err2Z;
876     GetErrors2( iRow, t, err2Y, err2Z );
877
878     if( !t.Filter2( h.Y(), h.Z(), err2Y, err2Z ) ) continue;
879
880     first = 0;      
881   }  
882   /*
883   Float_t cosPhi = iTrack.Param().GetCosPhi();
884   p0.Param().TransportToX(ID2Row( iTrack.PointID()[0] ).X());
885   p2.Param().TransportToX(ID2Row( iTrack.PointID()[1] ).X());  
886   if( p0.Param().GetCosPhi()*cosPhi<0 ){ // change direction
887   Float_t *par = p0.Param().Par();
888   Float_t *cov = p0.Param().Cov();
889   par[2] = -par[2]; // sin phi
890   par[3] = -par[3]; // DzDs
891     par[4] = -par[4]; // kappa
892     cov[3] = -cov[3];
893     cov[4] = -cov[4];
894     cov[6] = -cov[6];
895     cov[7] = -cov[7];
896     cov[10] = -cov[10];
897     cov[11] = -cov[11];
898     p0.Param().CosPhi() = -p0.Param().GetCosPhi();
899   }
900   */
901 #endif
902 }
903
904 GPUh() void AliHLTTPCCATracker::FitTrack( AliHLTTPCCATrack &/*track*/, Float_t */*t0[]*/ ) const 
905 {      
906   //* Fit the track   
907 #ifdef XXX
908   AliHLTTPCCAEndPoint &p2 = ID2Point(track.PointID()[1]);
909   AliHLTTPCCAHit &c0 = ID2Hit(fTrackHits[p0.TrackHitID()].HitID());     
910   AliHLTTPCCAHit &c1 = ID2Hit(fTrackHits[track.HitID()[1]].HitID());    
911   AliHLTTPCCAHit &c2 = ID2Hit(fTrackHits[p2.TrackHitID()].HitID());     
912   AliHLTTPCCARow &row0 = ID2Row(fTrackHits[p0.TrackHitID()].HitID());
913   AliHLTTPCCARow &row1 = ID2Row(fTrackHits[track.HitID()[1]].HitID());
914   AliHLTTPCCARow &row2 = ID2Row(fTrackHits[p2.TrackHitID()].HitID());
915   Float_t sp0[5] = {row0.X(), c0.Y(), c0.Z(), c0.ErrY(), c0.ErrZ() };
916   Float_t sp1[5] = {row1.X(), c1.Y(), c1.Z(), c1.ErrY(), c1.ErrZ() };
917   Float_t sp2[5] = {row2.X(), c2.Y(), c2.Z(), c2.ErrY(), c2.ErrZ() };
918   //cout<<"Fit track, points ="<<sp0[0]<<" "<<sp0[1]<<" / "<<sp1[0]<<" "<<sp1[1]<<" / "<<sp2[0]<<" "<<sp2[1]<<std::endl;
919   if( track.NHits()>=3 ){    
920     p0.Param().ConstructXYZ3(sp0,sp1,sp2,p0.Param().CosPhi(), t0);
921     p2.Param().ConstructXYZ3(sp2,sp1,sp0,p2.Param().CosPhi(), t0);
922     //p2.Param() = p0.Param();
923     //p2.Param().TransportToX(row2.X());
924     //p2.Param().Par()[1] = -p2.Param().Par()[1];
925     //p2.Param().Par()[4] = -p2.Param().Par()[4];
926   } else {
927     p0.Param().X() = row0.X();
928     p0.Param().Y() = c0.Y();
929     p0.Param().Z() = c0.Z();
930     p0.Param().Err2Y() = c0.ErrY()*c0.ErrY();
931     p0.Param().Err2Z() = c0.ErrZ()*c0.ErrZ();
932     p2.Param().X() = row2.X();
933     p2.Param().Y() = c2.Y();
934     p2.Param().Z() = c2.Z();
935     p2.Param().Err2Y() = c2.ErrY()*c2.ErrY();
936     p2.Param().Err2Z() = c2.ErrZ()*c2.ErrZ();
937   }
938 #endif
939 }
940
941
942
943 GPUd() void AliHLTTPCCATracker::GetErrors2( Int_t iRow, Float_t z, Float_t sinPhi, Float_t cosPhi, Float_t DzDs, Float_t &Err2Y, Float_t &Err2Z ) const
944 {
945   //
946   // Use calibrated cluster error from OCDB
947   //
948
949   z = CAMath::Abs((250.-0.275)-CAMath::Abs(z));
950   Int_t    type = (iRow<63) ? 0: ( (iRow>126) ? 1:2 );
951   Float_t cosPhiInv = CAMath::Abs(cosPhi)>1.e-2 ?1./cosPhi :0;
952   Float_t angleY = sinPhi*cosPhiInv ;
953   Float_t angleZ = DzDs*cosPhiInv ; // SG was bug??? 
954
955   Err2Y = fParam.GetClusterError2(0,type, z,angleY);  
956   Err2Z = fParam.GetClusterError2(1,type, z,angleZ);
957 }
958
959 GPUd() void AliHLTTPCCATracker::GetErrors2( Int_t iRow, const AliHLTTPCCATrackParam &t, Float_t &Err2Y, Float_t &Err2Z ) const
960 {
961   //
962   // Use calibrated cluster error from OCDB
963   //
964
965   Float_t z = CAMath::Abs((250.-0.275)-CAMath::Abs(t.GetZ()));
966   Int_t    type = (iRow<63) ? 0: ( (iRow>126) ? 1:2 );
967   Float_t cosPhiInv = CAMath::Abs(t.GetCosPhi())>1.e-2 ?1./t.GetCosPhi() :0;
968   Float_t angleY = t.GetSinPhi()*cosPhiInv ;
969   Float_t angleZ = t.GetDzDs()*cosPhiInv ; // SG was bug??? 
970
971   Err2Y = fParam.GetClusterError2(0,type, z,angleY);  
972   Err2Z = fParam.GetClusterError2(1,type, z,angleZ);
973 }
974
975
976 #if !defined(HLTCA_GPUCODE)  
977
978 GPUh() void AliHLTTPCCATracker::WriteEvent( std::ostream &out ) 
979 {
980   // write event to the file
981   for( Int_t iRow=0; iRow<fParam.NRows(); iRow++ ){
982     out<<fRows[iRow].FirstHit()<<" "<<fRows[iRow].NHits()<<std::endl;
983   } 
984   out<<fNHitsTotal<<std::endl;
985
986   Float_t y[fNHitsTotal], z[fNHitsTotal];
987  
988   for( Int_t iRow=0; iRow<fParam.NRows(); iRow++){
989     AliHLTTPCCARow &row = fRows[iRow];
990     Float_t y0 = row.Grid().YMin();
991     Float_t z0 = row.Grid().ZMin();
992     Float_t stepY = row.HstepY();
993     Float_t stepZ = row.HstepZ();
994     const uint4* tmpint4 = RowData() + row.FullOffset();
995     const ushort2 *hits = reinterpret_cast<const ushort2*>(tmpint4);
996     for( Int_t ih=0; ih<fRows[iRow].NHits(); ih++ ){
997       Int_t ihTot = row.FirstHit() + ih;
998       Int_t id = fHitInputIDs[ihTot];
999       ushort2 hh = hits[ih];
1000       y[id] = y0 + hh.x*stepY;
1001       z[id] = z0 + hh.y*stepZ;
1002     }
1003   }
1004   for( Int_t ih=0; ih<fNHitsTotal; ih++ ){
1005     out<<y[ih]<<" "<<z[ih]<<std::endl;
1006   }
1007 }
1008
1009 GPUh() void AliHLTTPCCATracker::ReadEvent( std::istream &in ) 
1010 {
1011   //* Read event from file 
1012
1013   Int_t rowFirstHit[Param().NRows()], rowNHits[Param().NRows()];  
1014   for( Int_t iRow=0; iRow<Param().NRows(); iRow++ ){
1015     in>>rowFirstHit[iRow]>>rowNHits[iRow];
1016   }
1017   Int_t nHits;
1018   in >> nHits;
1019   Float_t y[nHits], z[nHits];
1020   for( Int_t ih=0; ih<nHits; ih++ ){
1021     in>>y[ih]>>z[ih];
1022   }
1023   ReadEvent( rowFirstHit, rowNHits, y, z, nHits );
1024
1025
1026 GPUh() void AliHLTTPCCATracker::WriteTracks( std::ostream &out ) 
1027 {
1028   //* Write tracks to file 
1029  
1030   out<<fTimers[0]<<std::endl;
1031   out<<*fNOutTrackHits<<std::endl;
1032   for( Int_t ih=0; ih<*fNOutTrackHits; ih++ ){
1033     out<< fOutTrackHits[ih]<<" ";
1034   }
1035   out<<std::endl;
1036   
1037   out<<*fNOutTracks<<std::endl;
1038
1039   for( Int_t itr=0; itr<*fNOutTracks; itr++ ){
1040     AliHLTTPCCAOutTrack &t = fOutTracks[itr];    
1041     AliHLTTPCCATrackParam p1 = t.StartPoint();  
1042     AliHLTTPCCATrackParam p2 = t.EndPoint();    
1043     out<< t.NHits()<<" ";
1044     out<< t.FirstHitRef()<<" ";
1045     out<< t.OrigTrackID()<<" ";
1046     out<<std::endl;
1047     out<< p1.X()<<" ";
1048     out<< p1.CosPhi()<<" ";
1049     out<< p1.Chi2()<<" ";
1050     out<< p1.NDF()<<std::endl;
1051     for( Int_t i=0; i<5; i++ ) out<<p1.Par()[i]<<" ";
1052     out<<std::endl;
1053     for( Int_t i=0; i<15; i++ ) out<<p1.Cov()[i]<<" ";
1054     out<<std::endl;
1055     out<< p2.X()<<" ";
1056     out<< p2.CosPhi()<<" ";
1057     out<< p2.Chi2()<<" ";
1058     out<< p2.NDF()<<std::endl;
1059     for( Int_t i=0; i<5; i++ ) out<<p2.Par()[i]<<" ";
1060     out<<std::endl;
1061     for( Int_t i=0; i<15; i++ ) out<<p2.Cov()[i]<<" ";
1062     out<<std::endl;
1063   }
1064 }
1065
1066 GPUh() void AliHLTTPCCATracker::ReadTracks( std::istream &in )
1067 {
1068   //* Read tracks  from file 
1069   in>>fTimers[0];
1070   in>>*fNOutTrackHits;  
1071
1072   for( Int_t ih=0; ih<*fNOutTrackHits; ih++ ){
1073     in>>fOutTrackHits[ih];
1074   }
1075   in>>*fNOutTracks;
1076
1077   for( Int_t itr=0; itr<*fNOutTracks; itr++ ){
1078     AliHLTTPCCAOutTrack &t = fOutTracks[itr];    
1079     AliHLTTPCCATrackParam p1, p2;
1080     Int_t i;
1081     Float_t f;
1082     in>> i; t.SetNHits( i );
1083     in>> i; t.SetFirstHitRef( i );
1084     in>> i; t.SetOrigTrackID( i );
1085     in>> f; p1.SetX( f );
1086     in>> f; p1.SetCosPhi( f );
1087     in>> f; p1.SetChi2( f );
1088     in>> i; p1.SetNDF( i );
1089     for( Int_t j=0; j<5; j++ ){ in>>f; p1.SetPar(j,f); }
1090     for( Int_t j=0; j<15; j++ ){ in>>f; p1.SetCov(j,f); }
1091     in>> f; p2.SetX( f );
1092     in>> f; p2.SetCosPhi( f );
1093     in>> f; p2.SetChi2( f );
1094     in>> i; p2.SetNDF( i );
1095     for( Int_t j=0; j<5; j++ ){ in>>f; p2.SetPar(j,f); }
1096     for( Int_t j=0; j<15; j++ ){ in>>f; p2.SetCov(j,f); }
1097     t.SetStartPoint( p1 );
1098     t.SetEndPoint( p2 );
1099   }
1100 }
1101 #endif