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