StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
AliHLTTPCCATrackletConstructor.cxx
1 // @(#) $Id: AliHLTTPCCATrackletConstructor.cxx,v 1.2 2016/07/15 14:43:33 fisyak Exp $
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 // Developed by: Igor Kulakov <I.Kulakov@gsi.de> *
11 // Maksym Zyzak <M.Zyzak@gsi.de> *
12 // *
13 // Permission to use, copy, modify and distribute this software and its *
14 // documentation strictly for non-commercial purposes is hereby granted *
15 // without fee, provided that the above copyright notice appears in all *
16 // copies and that both the copyright notice and this permission notice *
17 // appear in the supporting documentation. The authors make no claims *
18 // about the suitability of this software for any purpose. It is *
19 // provided "as is" without express or implied warranty. *
20 // *
21 //***************************************************************************
22 
23 
24 #include "AliHLTTPCCATracker.h"
25 #include "AliHLTTPCCATrackParamVector.h"
26 #include "AliHLTTPCCAParameters.h"
27 #include "AliHLTTPCCAGrid.h"
28 #include "AliHLTTPCCAMath.h"
29 #include "AliHLTTPCCADef.h"
30 #include "AliHLTTPCCATracklet.h"
31 #include "AliHLTTPCCATrackletConstructor.h"
32 #include "AliHLTArray.h"
33 #include "debug.h"
34 #include <iomanip>
35 
36 #ifdef MAIN_DRAW
37 #include "AliHLTTPCCADisplay.h"
38 #endif // MAIN_DRAW
39 
40 #include <limits>
41 #include <Vc/limits>
42 
43 using std::endl;
44 
45 /*
46  * Tracklets start as WaitingForFitLinkedHits. Once the start row is reached they go to
47  * FitLinkedHits. When the seed is fitted they go to ExtrapolateUp. At some row upwards
48  * extrapolation will stop and the tracklet goes into WaitingForExtrapolateDown. Once all tracklets
49  * in the vector are done with ExtrapolateUp down extrapolation starts from the highest row where
50  * fitting was done. Once a tracklet reaches its fEndRow it goes to the ExtrapolateDown state until
51  * extrapolation is done, when it will go into DoneStage.
52  *
53  * Tracklets in the vector whose trackIndex is < nTracks are in the NullStage.
54  */
55 enum Stage_t {
56  WaitingForFitLinkedHits = 0,
57  FitLinkedHits,
58  ExtrapolateUp,
59  WaitingForExtrapolateDown,
60  ExtrapolateDown,
61  DoneStage,
62  NullStage
63 };
64 
65 struct AliHLTTPCCATrackletConstructor::TrackMemory {
66  short_v fStartRow; // min row index with only the fitted part
67  ushort_v fEndRow; // max row index with only the fitted part
68  short_v fFirstRow; // min row index with the extrapolated parts
69  ushort_v fLastRow; // max row index with the extrapolated parts
70  short_v fCurrentHitIndex; // indef of the current hit (for chain-fit stage it is next candidate to filtrate, for extra stage it is last filtered hit)
71  short_v fStage; // reco stage
72  ushort_v fNHits; // n track hits
73  short_v fRemainingGap; // n missed hits during search
74  sfloat_v fLastY; // Y of the last fitted cluster
75  sfloat_v fLastZ; // Z of the last fitted cluster
76  TrackParamVector fParam;
77  ushort_m fIsFragile; // is extrapolation and finding of additional hits needed
78  // ushort_m fIsOnChain;
79 
80  inline TrackMemory()
81  : fStartRow( Vc::Zero ),
82  fCurrentHitIndex( -1 ),
83  fStage( Vc::Zero ),
84  fNHits( 2 ), // the first two hits are certain
85  fRemainingGap( AliHLTTPCCAParameters::MaximumExtrapolationRowGap ),
86  fLastY( Vc::Zero ),
87  fLastZ( Vc::Zero ),
88  fIsFragile( Vc::Zero )
89  // fIsOnChain( Vc::Zero )
90  {}
91 
92  short_m IsInvalid() {
93  return fNHits < 3 ||
94  static_cast<short_m>( CAMath::Abs( fParam.SinPhi() ) > .999f ||
95  fParam.Chi2() >= sfloat_v(15.f)*static_cast<sfloat_v>(fParam.NDF()) );
96  }
97 
98  std::ostream &operator<<( std::ostream &out ) {
99 #define SEP endl
100 //#define SEP " "
101  out << "SR " << fStartRow << SEP;
102  out << "ER " << fEndRow << SEP;
103  out << "FR " << fFirstRow << SEP;
104  out << "LR " << fLastRow << SEP;
105  out << "CI " << fCurrentHitIndex << SEP;
106  out << "St " << fStage << SEP;
107  out << "NH " << fNHits << SEP;
108  out << "RG " << fRemainingGap << SEP;
109  out << "LY " << fLastY << SEP;
110  out << "LZ " << fLastZ << SEP;
111  out << "IF " << fIsFragile << SEP;
112  out << fParam;
113  return out;
114 #undef SEP
115  };
116 };
117 
118 class InitTracklets
119 {
120  public:
121  InitTracklets( AliHLTTPCCATrackletConstructor::TrackMemory &_r, SliceData &data,
122  const Tracker &tracker,
123  TrackletVector &trackletVector, const ushort_v &_trackIndex )
124  : r( _r ), fData( data ), fTracker( tracker ),
125  fTrackletVector( trackletVector ), trackIndex( _trackIndex ) {}
126 
127  void operator()( int rowIndex );
128 
129  private:
130  AliHLTTPCCATrackletConstructor::TrackMemory &r;
131  SliceData &fData;
132  const Tracker &fTracker;
133  TrackletVector &fTrackletVector;
134  const ushort_v &trackIndex;
135 };
136 
140 void AliHLTTPCCATrackletConstructor::FitTracklet( TrackMemory &r, const int rowIndex,
141  const ushort_v trackIndex, TrackletVector &trackletVector )
142 {
143  const AliHLTTPCCARow &row = fData.Row( rowIndex );
144 
145  const short_m active = r.fStage == FitLinkedHits; // active chains (ones, which contains no fitted hits).
146 
147  assert( AliHLTTPCCAParameters::RowStep == 1 );
148  // if (rowStep == 2){
149  // active = active && ( r.fStartRow & short_v( Vc::One ) ) == ( rowIndex & 1 );
150  // }
151 
152  debugF() << "FitTracklet(" << rowIndex << ") stage: " << r.fStage << " r.fStartRow: " << r.fStartRow << " active: " << active << endl;
153 
154  assert( active == ( active && r.fStage == FitLinkedHits ) );
155  assert( active == ( active && r.fCurrentHitIndex >= 0 ) );
156  assert( active == ( active && rowIndex > r.fStartRow ) );
157  //assert( r.fActive == ( ( rowIndex - r.fStartRow ) % 2 == 0 ) );
158  ASSERT( active == ( active && r.fCurrentHitIndex < row.NHits() ),
159  active << r.fCurrentHitIndex << row.NHits() );
160 
161  const ushort_v oldHitIndex = static_cast<ushort_v>( r.fCurrentHitIndex );
162  // fData.SetHitAsUsedInTrack( fData.Row( rowIndex - 1 ), static_cast<ushort_v>( r.fCurrentHitIndex ), active ); // TODO mark first hit
163 
164 
165  // const short_v aaaaaItIsNotUsedButDoHelpsALotToSpeedUpAaaaaa = short_v(fData.HitDataIsUsed( row ), Vc::IndexesFromZero, short_m(Vc::Zero)); // TODO understand
166  const short_v isUsed(fData.HitDataIsUsed( row ), static_cast<ushort_v>(oldHitIndex), active);
167  const short_m fitMask = active && ( (isUsed != short_v(2)) && (isUsed != short_v(3)) ); // mask to add a new hit. // don't take hits from other tracks. In order to reduce calculations.
168 
169  const sfloat_v x = fData.RowX( rowIndex ); // convert to sfloat_v once now
170  const sfloat_v y( fData.HitDataY( row ), oldHitIndex, fitMask);
171  const sfloat_v z( fData.HitDataZ( row ), oldHitIndex, fitMask);
172 
173  debugF() << "x, y, z: " << x << y << z << endl;
174 
175  sfloat_m activeF( static_cast<sfloat_m>( fitMask ) ); // create float_v mask
176  // correct SinPhi if it is necessary
177  // calculate displacement to previous hit
178  const sfloat_v dx = x - r.fParam.X();
179  const sfloat_v dy = y - r.fLastY;
180  const sfloat_v dz = z - r.fLastZ;
181  debugF() << "dx, dy, dz: " << dx << dy << dz << endl;
182  r.fLastY( activeF ) = y;
183  r.fLastZ( activeF ) = z;
184 
185  sfloat_v err2Y, err2Z;
186  sfloat_v sinPhi = r.fParam.GetSinPhi();
187 
188 
189  const sfloat_m fragile =
190  static_cast<sfloat_m>( r.fNHits < ushort_v(AliHLTTPCCAParameters::MinimumHitsForFragileTracklet) ) || CAMath::Abs( r.fParam.SinPhi() ) >= .99f;
191  sinPhi( fragile ) = dy * CAMath::RSqrt( dx * dx + dy * dy );
192 
193 
194  assert( ( x == 0 && activeF ).isEmpty() );
195 
196  activeF = r.fParam.TransportToX( x, sinPhi, fTracker.Param().cBz(), -1.f, activeF );
197 
198  if ( !activeF.isEmpty() )
199  {
200  fTracker.GetErrors2( rowIndex, r.fParam, &err2Y, &err2Z );
201  const short_m hitAdded = static_cast<short_m>( r.fParam.Filter( activeF, y, z, err2Y, err2Z, .99f ) );
202  trackletVector.SetRowHits( rowIndex, trackIndex, static_cast<ushort_v>(r.fCurrentHitIndex), hitAdded );
203  ++r.fNHits( static_cast<ushort_m>(hitAdded) );
204  r.fEndRow( hitAdded ) = rowIndex;
205 
206  ASSERT( ( (row.NHits() > static_cast<ushort_v>(oldHitIndex) ) && hitAdded ) == hitAdded,
207  row.NHits() << static_cast<ushort_v>(oldHitIndex) << hitAdded );
208  fData.SetHitAsUsedInTrackFit( row, static_cast<ushort_v>( oldHitIndex ), hitAdded );
209  }
210 
211  r.fCurrentHitIndex.gather( fData.HitLinkUpData( row ), static_cast<ushort_v>( oldHitIndex ), active ); // set to next linked hit
212 
213  const short_m fittingDone = r.fCurrentHitIndex < 0 && active;
214  debugF() << "fittingDone = " << fittingDone << endl;
215 // if ( ISUNLIKELY( !fittingDone.isEmpty() ) ) {
216  ++r.fStage( fittingDone ); // goes to ExtrapolateUp if fitting is done (no other hit linked)
217  r.fLastRow( fittingDone ) = rowIndex;
218  debugF() << "hits: " << r.fNHits << ", sinphi: " << r.fParam.SinPhi() << endl;
219 
220 
221  const short_m invalidTracklet = r.IsInvalid() && fittingDone;
222 
223  r.fNHits.setZero( invalidTracklet );
224  r.fStage( invalidTracklet || (fittingDone && !r.fIsFragile) ) = DoneStage;
225 
226  debugF() << "r.fStage: " << r.fStage << endl;
227 // }
228 } // FitTracklet
229 
230 void AliHLTTPCCATrackletConstructor::FindNextHit( TrackMemory &r, const AliHLTTPCCARow &row,
231  sfloat_v::Memory &dy_best, sfloat_v::Memory &dz_best, short_m &active)
232 {
233  const sfloat_v fY = r.fParam.Y();
234  const sfloat_v fZ = r.fParam.Z();
235  const ushort_v fIndYmin = row.Grid().GetBinBounded( fY - AliHLTTPCCAParameters::MinCellSize*.5f, fZ - AliHLTTPCCAParameters::MinCellSize*.5f );
236 
237  short_v::Memory fHitIndexes;
238  for(int iv=0; iv<short_v::Size; iv++)
239  fHitIndexes[iv] = -1; // not member
240 
241  foreach_bit( int trackletIndex, active ) {
242  float minRadius2 = std::numeric_limits<float>::max();
243 
244  const float y = fY[trackletIndex];
245  const float z = fZ[trackletIndex];
246  const int indYmin = fIndYmin[trackletIndex];
247 
248  { // look at iY = 0,1; iZ = 0
249  const int end = fData.FirstHitInBin( row, indYmin + 2 );
250  for ( int hitIndex = fData.FirstHitInBin( row, indYmin ); hitIndex < end; hitIndex += float_v::Size ) {
251  float_v yz;
252  yz.load( fData.HitDataY( row ) + hitIndex, Vc::Unaligned);
253  const float_v dy = yz - y;
254  yz.load( fData.HitDataZ( row ) + hitIndex, Vc::Unaligned);
255 // const float_v dz = float_v::loadUnaligned( fData.HitDataZ( row ) + hitIndex ) - z;
256  const float_v dz = yz - z;
257  float_v radius2 = std::numeric_limits<float_v>::max();
258  radius2( uint_v( Vc::IndexesFromZero ) + hitIndex < end ) = dy * dy + dz * dz; // XXX Manhattan distance
259 // std::cout << y << " dy = " << dy << std::endl <<
260 // z << " dz = " << dz << std::endl <<
261 // " Manhattan distance = " << radius2 << std::endl;
262  const float min = radius2.min();
263  if ( min < minRadius2 ) {
264  minRadius2 = min;
265  int i = ( radius2 == min ).firstOne();
266  assert( minRadius2 == radius2[i] );
267  dy_best[trackletIndex] = dy[i];
268  dz_best[trackletIndex] = dz[i];
269  fHitIndexes[trackletIndex] = hitIndex + i;
270  }
271  }
272  }
273  { // look at iY = 0,1; iZ = 1
274  const int nY = row.Grid().Ny();
275  const int end = fData.FirstHitInBin( row, indYmin + nY + 2 );
276  for ( int hitIndex = fData.FirstHitInBin( row, indYmin + nY ); hitIndex < end; hitIndex += float_v::Size ) {
277  float_v yz;
278  yz.load( fData.HitDataY( row ) + hitIndex, Vc::Unaligned );
279  const float_v dy = yz - y;
280  yz.load( fData.HitDataZ( row ) + hitIndex, Vc::Unaligned );
281  const float_v dz = yz - z;
282 // const float_v dy = float_v::loadUnaligned( fData.HitDataY( row ) + hitIndex ) - y;
283 // const float_v dz = float_v::loadUnaligned( fData.HitDataZ( row ) + hitIndex ) - z;
284  float_v radius2 = std::numeric_limits<float_v>::max();
285  radius2( uint_v( Vc::IndexesFromZero ) + hitIndex < end ) = dy * dy + dz * dz; // XXX Manhattan distance
286  const float min = radius2.min();
287  if ( min < minRadius2 ) {
288  minRadius2 = min;
289  int i = ( radius2 == min ).firstOne();
290  assert( minRadius2 == radius2[i] );
291  dy_best[trackletIndex] = dy[i];
292  dz_best[trackletIndex] = dz[i];
293  fHitIndexes[trackletIndex] = hitIndex + i;
294  }
295  }
296  }
297  }
298 
299  r.fCurrentHitIndex( active ) = static_cast<short_v>( fHitIndexes );
300 // r.fCurrentHitIndex( static_cast<short_m>( row.NHits() <= 0 ) ) = -1; // sometimes row.NHits() == 0 but firstHitInBin != 0 CHECKME. But same checked in Extend\ExtrapolateTracklet
301  active &= (r.fCurrentHitIndex != short_v(-1));
302 } // FindNextHit
303 
307 short_m AliHLTTPCCATrackletConstructor::ExtrapolateTracklet( TrackMemory &r, const int rowIndex, const ushort_v trackIndex,
308  TrackletVector &trackletVector, const bool dir, const short_m &mask )
309 {
310  //assert( active == ( active && ( r.fStage == ExtrapolateUp || r.fStage == ExtrapolateDown ) ) );
311  // reconstruction of tracklets, tracklets update step
312 
313  const AliHLTTPCCARow &row = fData.Row( rowIndex );
314 
315  short_m active = mask;
316  --r.fRemainingGap( active );
317  assert( r.fRemainingGap >= 0 || !active );
318  debugF() << "ExtrapolateTracklet(" << rowIndex << ") fRemainingGap: " << r.fRemainingGap << endl;
319 
320 
321  const sfloat_v x = fData.RowX( rowIndex );
322 
323 // // correct SinPhi if it is necessary
324 // // calculate displacement to previous hit
325 // const sfloat_v dx = x - r.fParam.X();
326 // const sfloat_v dy = y - r.fLastY;
327 // const sfloat_v dz = z - r.fLastZ;
328 // debugF() << "dx, dy, dz: " << dx << dy << dz << endl;
329 // r.fLastY( activeF ) = y;
330 // r.fLastZ( activeF ) = z;
331 //
332 // sfloat_v err2Y, err2Z;
333 // sfloat_v sinPhi = r.fParam.GetSinPhi();
334 //
335 // const sfloat_v ri = sfloat_v( Vc::One ) / CAMath::Sqrt( dx * dx + dy * dy ); // RSqrt
336 // const sfloat_m fragile =
337 // static_cast<sfloat_m>( r.fNHits < AliHLTTPCCAParameters::MinimumHitsForFragileTracklet ) || CAMath::Abs( r.fParam.SinPhi() ) >= .99f;
338 // sinPhi( fragile ) = dy * ri;
339 
340 
341  sfloat_m activeF( active ); // get float mask
342  assert( ( x == 0 && activeF ).isEmpty() );
343 
344 // sfloat_m transport = r.fParam.TransportToXWithMaterial( x, fTracker.Param().Bz(), .99f );
345 // sfloat_m transport = r.fParam.TransportToX( x, sinPhi, fTracker.Param().cBz(), .99f, activeF );
346  sfloat_m transport = r.fParam.TransportToX( x, r.fParam.SinPhi(), fTracker.Param().cBz(), .99f, activeF );
347  activeF &= transport;
348 
349  active = static_cast<short_m>( activeF );
350 
351  if ( row.NHits() < 1 || active.isEmpty() ) {
352  ++r.fStage( r.fRemainingGap == 0 && mask ); // go to WaitingForExtrapolateDown or DoneStage if the gap got too big
353 #ifndef NODEBUG
354  debugF() << "r.fStage: " << r.fStage << " after ExtrapolateTracklet found empty row and fRemainingGap == " << r.fRemainingGap << endl;
355 #endif
356  // no hits in this row, skip (don't do this before TransportToX)
357  ASSERT( r.fRemainingGap >= 0,
358  r.fRemainingGap );
359  return short_m( Vc::Zero );
360  }
361 
362  //std::cout << std::setw( 4 ) << rowIndex << ": " << active << endl;
363 
364  // find next hit
365 
366 
367  short_m linkMask = active; linkMask &= short_m(Vc::Zero);
368  // if ( dir == 1 ) {
369  // if ( rowIndex-1 >= 0 ) {
370  // linkMask &= r.fLastRow == rowIndex-1; // check only hits on previous row. Others can't have link on current row.
371  // const AliHLTTPCCARow &prevRow = fData.Row( rowIndex-1 );
372  // r.fCurrentHitIndex.gather( fData.HitLinkUpData( prevRow ), static_cast<ushort_v>( r.fCurrentHitIndex ), linkMask ); // set to next linked hit
373  // }
374  // else linkMask &= short_m(Vc::Zero);
375  // }
376  // else {
377  // if( rowIndex+1 < fTracker.Param().NRows() ) {
378  // linkMask &= r.fFirstRow == rowIndex+1;
379  // const AliHLTTPCCARow &prevRow = fData.Row( rowIndex+1 );
380  // r.fCurrentHitIndex.gather( fData.HitLinkDownData( prevRow ), static_cast<ushort_v>( r.fCurrentHitIndex ), linkMask );
381  // }
382  // else linkMask &= short_m(Vc::Zero);
383  // }
384  // linkMask &= (r.fCurrentHitIndex != -1);
385 
386 
387  // sfloat_v y( fData.HitDataY( row ), static_cast<ushort_v>( r.fCurrentHitIndex ), linkMask );
388  // sfloat_v z( fData.HitDataZ( row ), static_cast<ushort_v>( r.fCurrentHitIndex ), linkMask );
389 
390  // foreach_bit( int trackletIndex, linkMask ) { // TODO optimize
391  // const float yT = r.fParam.Y()[trackletIndex];
392  // const float zT = r.fParam.Z()[trackletIndex];
393  // dy_tmp[trackletIndex] = y[trackletIndex] - yT;
394  // dz_tmp[trackletIndex] = z[trackletIndex] - zT;
395  // }
396 
397  // try to find hit by trying all closed hits
398  sfloat_v::Memory dy_tmp; dy_tmp = sfloat_v(Vc::Zero); // initialize to avoid arithmetic exception.
399  sfloat_v::Memory dz_tmp; dz_tmp = sfloat_v(Vc::Zero); // initialize to avoid arithmetic exception.
400  short_m findMask = active && !linkMask;
401 
402  FindNextHit(r, row, dy_tmp, dz_tmp, findMask);
403 
404  const sfloat_v dy( dy_tmp );
405  const sfloat_v dz( dz_tmp );
406  active = linkMask || findMask;
407 
408  // XXX use this branch? with 8 or 16 tracklets extrapolated at the same time this might be an
409  // unlikely return? Also see above.
410 //X if ( active.isEmpty() ) {
411 //X return;
412 //X }
413 
414  // const short_v isUsed(fData.HitDataIsUsed( row ), static_cast<ushort_v>(r.fCurrentHitIndex), active);
415  // active &= (isUsed != short_v( 2 )) && (isUsed != short_v( 1 )); // don't take hits from other chains. In order to reduce calculations.
416  // active &= (isUsed != short_v( 2 ));
417 
418  // check if found hit is acceptable
419  sfloat_v err2Y, err2Z;
420  fTracker.GetErrors2( rowIndex, r.fParam, &err2Y, &err2Z );
421 
422  const sfloat_v kFactor = AliHLTTPCCAParameters::HitPickUpFactor * AliHLTTPCCAParameters::HitPickUpFactor * 3.5f * 3.5f;
423  const sfloat_v two( 2.f );
424 
425  const sfloat_v sy2 = CAMath::Min( two, kFactor * ( r.fParam.GetErr2Y() + err2Y ) );
426  const sfloat_v sz2 = CAMath::Min( two, kFactor * ( r.fParam.GetErr2Z() + err2Z ) );
427 
428  activeF = static_cast<sfloat_m>( active );
429  debugF() << "activeF: " << activeF;
430  activeF &= dy * dy <= sy2 && dz * dz <= sz2;
431  debugF() << " -> " << activeF;
432 
433 
434  sfloat_m filtred = r.fParam.FilterDelta( activeF, dy, dz, err2Y, err2Z, .99f );
435  activeF &= filtred;
436  // r.fCurrentHitIndex( static_cast<short_m>(!activeF) && active ) = -1; // delete indexes of noaccepted hits
437  debugF() << " -> " << activeF;
438  active = static_cast<short_m>( activeF );
439  debugF() << " -> " << active << endl;
440 
441  debugF() << "SetRowHits( " << rowIndex << ", " << r.fCurrentHitIndex << ", " << active << ")" << endl;
442 
443  ASSERT( ( (row.NHits() > r.fCurrentHitIndex) && active) == active,
444  row.NHits() << r.fCurrentHitIndex << active );
445  fData.SetHitAsUsedInTrackExtend( row, static_cast<ushort_v>(r.fCurrentHitIndex), active );
446 
447  trackletVector.SetRowHits( rowIndex, trackIndex, static_cast<ushort_v>(r.fCurrentHitIndex), active );
448  ++r.fNHits( static_cast<ushort_m>(active) );
449  r.fRemainingGap( active ) = AliHLTTPCCAParameters::MaximumExtrapolationRowGap;
450  ++r.fStage( r.fRemainingGap == 0 && mask ); // go to WaitingForExtrapolateDown or DoneStage if the gap got too big
451 // r.fStage( r.fNHits > 10 && mask ) = DoneStage; // don't need long tracklets. merger will do work.
452 
453  r.fLastRow( active && short_m( dir) ) = rowIndex;
454  r.fFirstRow( active && short_m(!dir) ) = rowIndex;
455  return active;
456 } // ExtrapolateTracklet
457 
458 short_m AliHLTTPCCATrackletConstructor::ExtendTracklet( TrackMemory &r, const int rowIndex, const ushort_v trackIndex,
459  TrackletVector &trackletVector, const bool dir, const short_m &mask )
460 {
461 
462  //assert( activeExtraMask == ( activeExtraMask && ( r.fStage == ExtrapolateUp || r.fStage == ExtrapolateDown ) ) );
463  // reconstruction of tracklets, tracklets update step
464  assert( AliHLTTPCCAParameters::RowStep == 1 );
465 
466 
467  short_m activeExtraMask = mask;
468  const short_m activeFitMask = (r.fStage == FitLinkedHits); // active chains (ones, which contains no fitted hits).
469 
470  const AliHLTTPCCARow &row = fData.Row( rowIndex );
471 
472  assert( activeFitMask == ( activeFitMask && r.fStage == FitLinkedHits ) );
473  assert( activeFitMask == ( activeFitMask && r.fCurrentHitIndex >= 0 ) );
474  assert( activeFitMask == ( activeFitMask && rowIndex > r.fStartRow ) );
475  ASSERT( activeFitMask == ( activeFitMask && r.fCurrentHitIndex < row.NHits() ),
476  activeFitMask << r.fCurrentHitIndex << row.NHits() );
477  assert( (activeFitMask && activeExtraMask).isEmpty() );
478 
479 
480  const ushort_v oldHitIndex = static_cast<ushort_v>( r.fCurrentHitIndex );
481 
482 
483 
484  --r.fRemainingGap( activeExtraMask );
485  debugF() << "ExtrapolateTracklet(" << rowIndex << ") fRemainingGap: " << r.fRemainingGap << endl;
486 
487 
488  // -- TRANSPORT --
489 
490  const sfloat_v x = fData.RowX( rowIndex );
491 
492  // const short_v aaaaaItIsNotUsedButDoesHelpALotToSpeedUpAaaaaa = short_v(fData.HitDataIsUsed( row ), Vc::IndexesFromZero, short_m(Vc::Zero)); // todo understand. Here it takes a lot of time...
493  const short_v isUsed(fData.HitDataIsUsed( row ), static_cast<ushort_v>(oldHitIndex), activeFitMask); // Here it takes a lot of time...
494  const short_m fitMask = activeFitMask && (isUsed != short_v(2)) && (isUsed != short_v(3)); // mask to add a new hit. // don't take hits from other tracks. In order to reduce calculations.
495 
496  sfloat_m activeFitMaskF( static_cast<sfloat_m>( fitMask ) ); // create float_v mask
497  sfloat_m activeExtraMaskF( static_cast<sfloat_m>(activeExtraMask) );
498 
499  assert( row.NHits() > oldHitIndex || !fitMask );
500  sfloat_v y( fData.HitDataY( row ), oldHitIndex, activeFitMaskF);
501  sfloat_v z( fData.HitDataZ( row ), oldHitIndex, activeFitMaskF);
502 
503 
504  const sfloat_v dx = x - r.fParam.X();
505 
506  sfloat_v dy = y - r.fLastY;
507  sfloat_v dz = z - r.fLastZ;
508  r.fLastY( activeFitMaskF ) = y;
509  r.fLastZ( activeFitMaskF ) = z;
510 
511  sfloat_v sinPhi = r.fParam.GetSinPhi();
512 
513  const sfloat_m fragile = activeFitMaskF &&
514  ( static_cast<sfloat_m>( r.fNHits < ushort_v(AliHLTTPCCAParameters::MinimumHitsForFragileTracklet) ) || CAMath::Abs( r.fParam.SinPhi() ) >= .99f );
515  sinPhi( fragile ) = dy * CAMath::RSqrt( dx * dx + dy * dy );
516 
517  sfloat_v maxSinPhi = .99f;
518  maxSinPhi( activeFitMaskF ) = -1.f;
519 
520 
521 // sfloat_m transport = r.fParam.TransportToXWithMaterial( x, fTracker.Param().Bz(), .99f );
522 // sfloat_m transport = r.fParam.TransportToX( x, sinPhi, fTracker.Param().cBz(), .99f, activeExtraMaskF );
523 
524  assert( ( (x == 0) && (activeExtraMaskF || activeFitMaskF) ).isEmpty() );
525 
526  sfloat_m transport = r.fParam.TransportToX( x, sinPhi, fTracker.Param().cBz(), maxSinPhi, activeExtraMaskF || activeFitMaskF );
527  activeExtraMaskF &= transport;
528  activeFitMaskF &= transport;
529 
530 
531  if ( row.NHits() < 1 || ( activeExtraMaskF || activeFitMaskF ).isEmpty() ) {
532  // no hits in this row, skip (don't do this before TransportToX)
533  // end with extrapolation
534  ++r.fStage( r.fRemainingGap == 0 && mask ); // go to WaitingForExtrapolateDown or DoneStage if the gap got too big
535 #ifndef NODEBUG
536  debugF() << "r.fStage: " << r.fStage << " after ExtrapolateTracklet found empty row and fRemainingGap == " << r.fRemainingGap << endl;
537 #endif
538 
539  // end with chain fit
540  r.fCurrentHitIndex.gather( fData.HitLinkUpData( row ), static_cast<ushort_v>( oldHitIndex ), activeFitMask ); // prepare new hit for fit // TODO 2 dir??
541 
542  const short_m fittingDone = r.fCurrentHitIndex < 0 && activeFitMask;
543  ++r.fStage( fittingDone ); // goes to ExtrapolateUp if fitting is done (no other hit linked)
544  r.fLastRow( fittingDone ) = rowIndex;
545 
546  const short_m invalidTracklet = r.IsInvalid() && fittingDone;
547  r.fNHits.setZero( invalidTracklet );
548  r.fStage( invalidTracklet || (fittingDone && !r.fIsFragile) ) = DoneStage;
549 
550  return short_m( Vc::Zero );
551  }
552 
553  // -- FIND A NEXT HIT --
554 
555 
556  activeExtraMask = static_cast<short_m>( activeExtraMaskF );
557 
558  // try find hit by using links
559  short_m linkMask = short_m(Vc::Zero);//activeExtraMask && r.fIsOnChain; - // use links during extrapolation is turned off
560  // if ( !linkMask.isEmpty() ) {
561  // if ( dir == 1 ) {
562  // if ( rowIndex-1 >= 0 ) {
563  // linkMask &= r.fLastRow == rowIndex-1; // check only hits on previous row. Others can't have link on current row.
564  // const AliHLTTPCCARow &prevRow = fData.Row( rowIndex-1 );
565  // r.fCurrentHitIndex.gather( fData.HitLinkUpData( prevRow ), static_cast<ushort_v>( r.fCurrentHitIndex ), linkMask ); // set to next linked hit
566  // }
567  // else linkMask &= short_m(Vc::Zero);
568  // }
569  // else {
570  // if( rowIndex+1 < fTracker.Param().NRows() ) {
571  // linkMask &= r.fFirstRow == rowIndex+1;
572  // const AliHLTTPCCARow &prevRow = fData.Row( rowIndex+1 );
573  // r.fCurrentHitIndex.gather( fData.HitLinkDownData( prevRow ), static_cast<ushort_v>( r.fCurrentHitIndex ), linkMask );
574  // }
575  // else linkMask &= short_m(Vc::Zero);
576  // }
577  // r.fIsOnChain &= (r.fCurrentHitIndex != -1); // chain is ended
578  // linkMask &= (r.fCurrentHitIndex != -1);
579 
580  // assert( row.NHits() > static_cast<ushort_v>( r.fCurrentHitIndex ) || !linkMask );
581  // y.gather( fData.HitDataY( row ), static_cast<ushort_v>( r.fCurrentHitIndex ), static_cast<sfloat_m>(linkMask) );
582  // z.gather( fData.HitDataZ( row ), static_cast<ushort_v>( r.fCurrentHitIndex ), static_cast<sfloat_m>(linkMask) );
583  // }
584 
585  dy = y - r.fParam.Y();
586  dz = z - r.fParam.Z();
587 
588  // try to find hit by trying all closed hits
589  sfloat_v::Memory dy_tmp, dz_tmp;
590  short_m findMask = activeExtraMask && !linkMask; // try to find only if there is no chain
591  if ( !findMask.isEmpty() ) {
592  // const short_v prev(r.fCurrentHitIndex); // TODO oldHI
593 
594  FindNextHit(r, row, dy_tmp, dz_tmp, findMask);
595 
596  // if (dir == 1) {
597  // const short_v prev2( fData.HitLinkDownData( row ), static_cast<ushort_v>( r.fCurrentHitIndex ), findMask );
598  // r.fIsOnChain = prev == prev2 && findMask;
599  // }
600  // else {
601  // const short_v prev2( fData.HitLinkUpData( row ), static_cast<ushort_v>( r.fCurrentHitIndex ), findMask );
602  // r.fIsOnChain = prev == prev2 && findMask;
603  // }
604  }
605  activeExtraMask = linkMask || findMask;
606  activeExtraMaskF = static_cast<sfloat_m>( activeExtraMask );
607 
608  dy( static_cast<sfloat_m>( findMask ) ) = sfloat_v( dy_tmp );
609  dz( static_cast<sfloat_m>( findMask ) ) = sfloat_v( dz_tmp );
610 
611  // XXX use this branch? with 8 or 16 tracklets extrapolated at the same time this might be an
612  // unlikely return? Also see above.
613 //X if ( activeExtraMask.isEmpty() ) {
614 //X return;
615 //X }
616 
617 // const short_v isUsed2(fData.HitDataIsUsed( row ), static_cast<ushort_v>(r.fCurrentHitIndex), activeExtraMask);
618 // activeExtraMask &= (isUsed2 != short_v( 2 )) && (isUsed2 != short_v( 1 )); // don't take hits from other chains. In order to reduce calculations.
619 // activeExtraMask &= (isUsed2 != short_v( 2 ));
620 
621  // check if found hit is acceptable
622  sfloat_v err2Y, err2Z;
623  fTracker.GetErrors2( rowIndex, r.fParam, &err2Y, &err2Z );
624 
625  const sfloat_v kFactor = AliHLTTPCCAParameters::HitPickUpFactor * AliHLTTPCCAParameters::HitPickUpFactor * 3.5f * 3.5f;
626  const sfloat_v two( 2.f );
627 
628  const sfloat_v sy2 = CAMath::Min( two, kFactor * ( r.fParam.GetErr2Y() + err2Y ) );
629  const sfloat_v sz2 = CAMath::Min( two, kFactor * ( r.fParam.GetErr2Z() + err2Z ) );
630 
631 
632  activeExtraMaskF &= dy * dy <= sy2 && dz * dz <= sz2;
633 
634  activeFitMaskF &= (dy * dy <= sy2 && dz * dz <= sz2) || (r.fNHits < 3); // very strong cut for ghosts.
635 
636 
637  // -- FILTER THE NEXT HIT --
638 
639  const sfloat_m filtred = r.fParam.FilterDelta( activeExtraMaskF || activeFitMaskF, dy, dz, err2Y, err2Z, .99f );
640  activeExtraMaskF &= filtred;
641  activeExtraMask = static_cast<short_m>( activeExtraMaskF );
642  activeFitMaskF &= filtred;
643  const short_m hitAdded = static_cast<short_m>( activeFitMaskF );
644 
645 
646  // -- SAVE THE NEXT HIT --
647 
648 
649  trackletVector.SetRowHits( rowIndex, trackIndex, static_cast<ushort_v>(r.fCurrentHitIndex), activeExtraMask || hitAdded );
650 
651  r.fCurrentHitIndex.gather( fData.HitLinkUpData( row ), static_cast<ushort_v>( oldHitIndex ), activeFitMask ); // prepare new hit for fit
652 
653 
654  ASSERT( ( (row.NHits() > r.fCurrentHitIndex) && activeExtraMask) == activeExtraMask,
655  row.NHits() << r.fCurrentHitIndex << activeExtraMask );
656  fData.SetHitAsUsedInTrackExtend( row, static_cast<ushort_v>( r.fCurrentHitIndex ), activeExtraMask ); // TODO 1 function
657  fData.SetHitAsUsedInTrackFit( row, static_cast<ushort_v>( oldHitIndex ), hitAdded );
658 
659  ++r.fNHits( static_cast<ushort_m>( activeExtraMask || hitAdded ) );
660  r.fRemainingGap( activeExtraMask ) = AliHLTTPCCAParameters::MaximumExtrapolationRowGap;
661 
662  // -- MOVE TO NEXT STAGE --
663 
664  ++r.fStage( (r.fRemainingGap == 0) && mask ); // go to WaitingForExtrapolateDown or DoneStage if the gap got too big TODO not active
665  ASSERT( r.fRemainingGap >= 0,
666  r.fRemainingGap << mask << activeExtraMask );
667 // r.fStage( r.fNHits > 10 && mask ) = DoneStage; // don't need long tracklets. merdger will do work.
668 
669  const short_m fittingDone = r.fCurrentHitIndex < 0 && activeFitMask;
670  ++r.fStage( fittingDone ); // goes to ExtrapolateUp if fitting is done (no other hit linked)
671  r.fLastRow( fittingDone ) = rowIndex;
672 
673  const short_m invalidTracklet = r.IsInvalid() && fittingDone;
674  r.fNHits.setZero( invalidTracklet );
675  r.fStage( invalidTracklet || (fittingDone && !r.fIsFragile) ) = DoneStage; // TODO unmark used hits
676 
677  r.fEndRow( hitAdded ) = rowIndex;
678 
679  r.fLastRow( activeExtraMask && short_m( dir) ) = rowIndex;
680  r.fFirstRow( activeExtraMask && short_m(!dir) ) = rowIndex;
681  return activeExtraMask;
682 } // ExtendTracklet
683 
684 
685 void AliHLTTPCCATrackletConstructor::run()
686 {
687 
688  assert( *fTracker.NTracklets() < 32768 );
689  const short nTracks = *fTracker.NTracklets();
690 
691  for ( int trackIteration = 0; trackIteration * ushort_v::Size < nTracks; ++trackIteration ) {
692  const ushort_v trackIndex( ushort_v( Vc::IndexesFromZero ) + trackIteration * ushort_v::Size );
693  const ushort_m active = trackIndex < nTracks;
694  TrackMemory r;
695 
696  //std::cout << "----------------------------------------------------" << endl;
697  r.fStage( !active ) = NullStage;
698 
699  // if rowStep = 2 TrackletStartHits need to be sorted such that all even start rows come first. The odd start rows - last.
700  r.fStartRow.gather( fTracker.TrackletStartHits(), reinterpret_cast<short AliHLTTPCCAStartHitId::*>( &AliHLTTPCCAStartHitId::fRow ), trackIndex, active );
701  r.fCurrentHitIndex.gather( fTracker.TrackletStartHits(), reinterpret_cast<short int AliHLTTPCCAStartHitId::*>( &AliHLTTPCCAStartHitId::fHit ), trackIndex, active );
702  r.fEndRow = static_cast<ushort_v>( r.fStartRow );
703  r.fLastRow = static_cast<ushort_v>( r.fStartRow );
704  r.fStartRow( !active ) = std::numeric_limits<short_v>::max();
705  r.fFirstRow = r.fStartRow;
706  ushort_v length( fTracker.TrackletStartHits(), reinterpret_cast<unsigned short AliHLTTPCCAStartHitId::*>( &AliHLTTPCCAStartHitId::fLength ), trackIndex, active );
707  const ushort_v MaxNHitsForFragileTracklet(6);
708  r.fIsFragile = (length < MaxNHitsForFragileTracklet); // TODO parameter // TODO 5 is optimum for globalEff+timeSliceTracker. But one should take into account merger
709 
710  const sfloat_v zero( Vc::Zero );
711  const sfloat_v one( Vc::One );
712  r.fParam.SetSinPhi( zero );
713  r.fParam.SetDzDs( zero );
714  r.fParam.SetQPt( zero );
715  r.fParam.SetSignCosPhi( one );
716  r.fParam.SetChi2( zero );
717  r.fParam.SetNDF( -3 );
718  r.fParam.SetCov( 0, one );
719  r.fParam.SetCov( 1, zero );
720  r.fParam.SetCov( 2, one );
721  r.fParam.SetCov( 3, zero );
722  r.fParam.SetCov( 4, zero );
723  r.fParam.SetCov( 5, one );
724  r.fParam.SetCov( 6, zero );
725  r.fParam.SetCov( 7, zero );
726  r.fParam.SetCov( 8, zero );
727  r.fParam.SetCov( 9, one );
728  r.fParam.SetCov( 10, zero );
729  r.fParam.SetCov( 11, zero );
730  r.fParam.SetCov( 12, zero );
731  r.fParam.SetCov( 13, zero );
732  r.fParam.SetCov( 14, 10.f );
733 
734  // const short_m evenRows = ( r.fStartRow & short_v( Vc::One ) ) > short_v( Vc::Zero );
735  //
736  // // fStartRow should be sorted, with the one exception of odd/even boundary
737  // if ( r.fStartRow != r.fStartRow.sorted() ) {
738  // debugF() << r.fStartRow << r.fStartRow.sorted() << evenRows << evenRows.isMix() << endl;
739  // assert( r.fStartRow == r.fStartRow.sorted() || evenRows.isMix() );
740  // }
741  {
742  InitTracklets init( r, fData, fTracker, fTrackletVectors[trackIteration], trackIndex );
743  r.fStartRow.callWithValuesSorted( init );
744  }
745 
746 //#define USE_COUNTERS
747 #ifdef USE_COUNTERS
748  static unsigned int counters[6] = { 0, 0, 0, 0, 0, 0 };
749  counters[0]++;
750 #endif // USE_COUNTERS
751 
752  const int rowStep = AliHLTTPCCAParameters::RowStep;
753  { // fit and extrapolate upwards
754  int rowIndex = r.fStartRow.min() + rowStep*2;
755  const short_v activationRow = r.fStartRow + rowStep*2;
756  debugF() << "============================================= Start Fitting Upwards =============================================" << endl;
757  // TODO think about getting parameters of all tracks (fragile and !fragile) in the same point (begin or end of track). Needed by Merger
758 // #define DISABLE_HIT_SEARCH
759 #ifdef DISABLE_HIT_SEARCH
760  while ( rowIndex < fTracker.Param().NRows() && !( r.fStage <= FitLinkedHits ).isEmpty() ) {
761 #ifdef USE_COUNTERS
762  counters[1]++;
763 #endif // USE_COUNTERS
764 
765  ++r.fStage( rowIndex == activationRow ); // goes to FitLinkedHits on activation row
766  FitTracklet( r, rowIndex, trackIndex, fTrackletVectors[trackIteration] );
767  ++rowIndex;
768  }
769 
770 
771  const short_m ready = r.fStage <= DoneStage;
772  const sfloat_v x( fData.RowX(), sfloat_v::IndexType( r.fEndRow ) );
773 // const sfloat_v x( fData.RowX(), sfloat_v::IndexType( r.fStartRow ) );
774  debugF() << x << sfloat_v::IndexType( r.fEndRow ) << sfloat_v( fData.RowX(), sfloat_v::IndexType( Vc::IndexesFromZero ) ) << endl;
775  assert( ( x == 0 && static_cast<sfloat_m>( ready ) ).isEmpty() );
776  //assert ( ( (r.fParam.X() == x) && static_cast<sfloat_m>( ready )) == static_cast<sfloat_m>( ready ));
777 
778  const short_m transported = static_cast<short_m>( r.fParam.TransportToX(
779  x, fTracker.Param().cBz(), .999f, static_cast<sfloat_m>( ready ) ) );
780 
781  debugF() << "============================================= Stop Fitting Upwards ==============================================" << endl;
782 #ifdef MAIN_DRAW
783  if ( AliHLTTPCCADisplay::Instance().DrawType() == 1 ) {
784  foreach_bit( int ii, r.fStage < DoneStage ) {
785  TrackParam t( r.fParam, ii );
786  AliHLTTPCCADisplay::Instance().ClearView();
787  AliHLTTPCCADisplay::Instance().DrawSlice( &fTracker, 0 );
788  AliHLTTPCCADisplay::Instance().DrawSliceHits();
789  AliHLTTPCCADisplay::Instance().DrawTrackParam( t );
790  AliHLTTPCCADisplay::Instance().Ask();
791  }
792  }
793 #endif
794 #else // DISABLE_HIT_SEARCH
795  // fit all chains
796  while ( rowIndex < fTracker.Param().NRows() && ( r.fStage == ExtrapolateUp ).isEmpty() ) {
797 #ifdef USE_COUNTERS
798  counters[1]++;
799 #endif // USE_COUNTERS
800 
801  r.fStage( rowIndex == activationRow ) = FitLinkedHits; // goes to FitLinkedHits on activation row
802 
803  ExtendTracklet( r, rowIndex, trackIndex, fTrackletVectors[trackIteration], 1, short_m(Vc::Zero) );
804  // FitTracklet( r, rowIndex, trackIndex, fTrackletVectors[trackIteration] );
805  ++rowIndex;
806  }
807 
808  debugF() << "========================================== Start Extrapolating Upwards ==========================================" << endl;
809 
810  // some chains are fitted, some - are extrapolated
811  assert( r.fRemainingGap >= 0 );
812  //while ( rowIndex < fTracker.Param().NRows() && !( r.fStage <= FitLinkedHits ).isEmpty() ) {
813  while ( rowIndex < fTracker.Param().NRows() && !( r.fStage <= ExtrapolateUp ).isEmpty() ) {
814 #ifdef USE_COUNTERS
815  counters[2]++;
816 #endif // USE_COUNTERS
817 
818  r.fStage( rowIndex == activationRow ) = FitLinkedHits; // goes to FitLinkedHits on activation row
819  const short_m toExtrapolate = (r.fStage == ExtrapolateUp);
820  ExtendTracklet( r, rowIndex, trackIndex, fTrackletVectors[trackIteration], 1, toExtrapolate );
821  ++rowIndex;
822  }
823 
824  debugF() << "============================================= Stop Fitting Upwards ==============================================" << endl;
825 
826  // end extrapolate chains
827  short_m mask;
828  assert( r.fRemainingGap >= 0 );
829  while ( rowIndex < fTracker.Param().NRows() && !( mask = r.fStage == ExtrapolateUp ).isEmpty() ) {
830 #ifdef USE_COUNTERS
831  counters[3]++;
832 #endif // USE_COUNTERS
833 
834  assert ( r.fStage != FitLinkedHits );
835  ExtrapolateTracklet( r, rowIndex, trackIndex, fTrackletVectors[trackIteration], 1, mask );
836  ++rowIndex;
837  }
838  debugF() << "=========================================== Stop Extrapolating Upwards ==========================================" << endl;
839  }
840 
841  // r.fStage (r.fNHits < MaxNHitsForFragileTracklet) = DoneStage; // TODO investigate
842  { // extrapolate downwards
843  r.fRemainingGap = AliHLTTPCCAParameters::MaximumExtrapolationRowGap; // allow full gaps again
844  ++r.fStage( r.fStage == ExtrapolateUp ); // FitLinkedHits/ExtrapolateUp went so high that no gap put the tracklet into WaitingForExtrapolateDown
845  const short_m ready = r.fStage == WaitingForExtrapolateDown;
846  debugF() << "ready to extrapolate downwards: " << ready << endl;
847  ++r.fStage( ready ); // the wait is over
848 
849  if(1){ // set track parameters to x of end row of the fitting stage
850  //const sfloat_v x( fData.RowX(), sfloat_v::IndexType( r.fEndRow ) );
851  const sfloat_v x( fData.RowX(), sfloat_v::IndexType( r.fStartRow ), static_cast<sfloat_m>(ready));
852  debugF() << x << sfloat_v::IndexType( r.fEndRow ) << sfloat_v( fData.RowX(), sfloat_v::IndexType( Vc::IndexesFromZero ) ) << endl;
853  assert( ( x == 0 && static_cast<sfloat_m>( ready ) ).isEmpty() );
854  const short_m transported = static_cast<short_m>( r.fParam.TransportToX(
855  x, fTracker.Param().cBz(), .999f, static_cast<sfloat_m>( ready ) ) );
856  ++r.fStage( !transported ); // all those where transportation failed go to DoneStage
857  }
858 #ifdef MAIN_DRAW
859  if ( AliHLTTPCCADisplay::Instance().DrawType() == 1 ) {
860  foreach_bit( int ii, r.fStage < DoneStage ) {
861  TrackParam t( r.fParam, ii );
862  AliHLTTPCCADisplay::Instance().ClearView();
863  AliHLTTPCCADisplay::Instance().DrawSlice( &fTracker, 0 );
864  AliHLTTPCCADisplay::Instance().DrawSliceHits();
865  AliHLTTPCCADisplay::Instance().DrawTrackParam( t, 2 );
866  AliHLTTPCCADisplay::Instance().Ask();
867  }
868  }
869 #endif
870  debugF() << "========================================= Start Extrapolating Downwards =========================================" << endl;
871  // int rowIndex = r.fEndRow.max() - 1;
872  short_v tmpRow(r.fStartRow); // take only one hit from fitted chain
873  tmpRow(!active) = fTracker.Param().NRows()-1;
874  int rowIndex = tmpRow.max();
875  // int rowIndex = r.fLastRow.max() - 1;
876  // extrapolate along fitted chains
877  short_m mask;
878  const int minMaskedExtrRow = r.fStartRow.min();
879  while ( rowIndex >= minMaskedExtrRow && !( mask = r.fStage == ExtrapolateDown ).isEmpty() ) {
880 #ifdef USE_COUNTERS
881  counters[4]++;
882 #endif // USE_COUNTERS
883  // mask &= rowIndex < r.fEndRow;
884  mask &= rowIndex < r.fStartRow + 1;
885  mask &= ( ( rowIndex - r.fStartRow ) & short_v( std::numeric_limits<short_v::EntryType>::min() + 1 ) ) != short_v( Vc::Zero ); // CHECKME why do we need this?
886  ExtrapolateTracklet( r, rowIndex, trackIndex, fTrackletVectors[trackIteration], 0, mask );
887  --rowIndex;
888  }
889  // extrapolote downwards and find addition hits
890  while ( rowIndex >= 0 && !( mask = r.fStage == ExtrapolateDown ).isEmpty() ) {
891 #ifdef USE_COUNTERS
892  counters[5]++;
893 #endif // USE_COUNTERS
894  ExtrapolateTracklet( r, rowIndex, trackIndex, fTrackletVectors[trackIteration], 0, mask );
895  --rowIndex;
896  }
897  r.fFirstRow = CAMath::Min( r.fFirstRow, r.fStartRow );
898 #ifdef MAIN_DRAW
899  if ( AliHLTTPCCADisplay::Instance().DrawType() == 1 ) {
900  foreach_bit( int ii, r.fStage < NullStage ) {
901  TrackParam t( r.fParam, ii );
902  AliHLTTPCCADisplay::Instance().ClearView();
903  AliHLTTPCCADisplay::Instance().DrawSlice( &fTracker, 0 );
904  AliHLTTPCCADisplay::Instance().DrawSliceHits();
905  AliHLTTPCCADisplay::Instance().DrawTrackParam( t, 4 );
906  AliHLTTPCCADisplay::Instance().Ask();
907  }
908  }
909 #endif
910 
912  //std::cout << "fitted " << r.fParam.QPt() <<std::endl;
913  //r.fParam.SetQPt( qpt_approx);
915 
916 #endif // DISABLE_HIT_SEARCH
917  }
918 
919 #ifdef USE_COUNTERS
920  std::cout << "Counters= " << counters[0] << " " << counters[1] << " " << counters[2] << " " << counters[3] << " " << counters[4] << " " << counters[5] << std::endl;
921 #endif // USE_COUNTERS
922 
923  sfloat_m trackletOkF( r.fNHits >= ushort_v(AliHLTTPCCAParameters::MinimumHitsForTracklet) );
924  for ( int i = 0; i < 15; ++i ) {
925  trackletOkF &= CAMath::Finite( r.fParam.Cov()[i] );
926  }
927  for ( int i = 0; i < 5; ++i ) {
928  trackletOkF &= CAMath::Finite( r.fParam.Par()[i] );
929  }
930 
931  // 80 is row 0; if X is that small this track is garbage.
932  // XXX does this happen at all? If yes, under what conditions?
933  assert( ( r.fParam.X() > 50.f && trackletOkF ) == trackletOkF );
934  trackletOkF &= r.fParam.X() > 50.f // TODO: read from file!!!
935 
936 
937  // there must be errors and they must be positive or this track is garbage
938  && r.fParam.Err2QPt() > Vc::Zero;
939  trackletOkF &= r.fParam.Err2Y() > Vc::Zero && r.fParam.Err2Z() > Vc::Zero;
940  trackletOkF &= r.fParam.Err2SinPhi() > Vc::Zero && r.fParam.Err2DzDs() > Vc::Zero;
941 // trackletOkF &= ( r.fParam.Chi2()/static_cast<sfloat_v>(r.fParam.NDF()) < 25.f ); // TODO
942  debugF() << r.fParam << "-> trackletOk: " << trackletOkF << endl;
943 
944  const short_m trackletOk( trackletOkF );
945  r.fNHits.setZero( !trackletOk );
946 
948  //
950 
951  TrackletVector &tracklet = fTrackletVectors[trackIteration];
952  tracklet.SetNHits( r.fNHits );
953 
954  if ( !( r.fNHits > 0 ).isEmpty() ) {
955 #ifdef MAIN_DRAW
956  if ( AliHLTTPCCADisplay::Instance().DrawType() == 1 ) {
957  foreach_bit( int ii, r.fStage < DoneStage ) {
958  TrackParam t( r.fParam, ii );
959  AliHLTTPCCADisplay::Instance().ClearView();
960  AliHLTTPCCADisplay::Instance().DrawSlice( &fTracker, 0 );
961  AliHLTTPCCADisplay::Instance().DrawSliceHits();
962  AliHLTTPCCADisplay::Instance().DrawTrackParam( t, 2 );
963  AliHLTTPCCADisplay::Instance().Ask();
964  }
965  }
966 #endif
967  // start and end rows of the tracklet
968  tracklet.SetFirstRow( static_cast<ushort_v>( r.fFirstRow ) );
969  tracklet.SetLastRow( r.fLastRow );
970 
972  const sfloat_m MinQPt = CAMath::Abs(r.fParam.QPt()) < AliHLTTPCCAParameters::MinimumQPt;
973  sfloat_v NewQPt = r.fParam.QPt();
974  NewQPt(MinQPt) = AliHLTTPCCAParameters::MinimumQPt;
975  r.fParam.SetQPt(NewQPt);
976  // r.fParam.SetQPt( CAMath::Max( AliHLTTPCCAParameters::MinimumQPt, r.fParam.QPt() ) );
978  tracklet.SetParam( r.fParam );
979 
980  debugTS() << r.fFirstRow << r.fLastRow << endl;
981  debugTS() << "set hit weigths from row " << r.fFirstRow.min() << " until row " << r.fLastRow.max() << endl;
982  const ushort_v &weight = SliceData::CalculateHitWeight( r.fNHits, trackIndex );
983  // for all rows where we have a hit let the fTracker know what weight our hits have
984  for ( int rowIndex = r.fFirstRow.min(); rowIndex <= r.fLastRow.max(); ++rowIndex ) {
985  const ushort_v &hitIndex = tracklet.HitIndexAtRow( rowIndex );
986  debugTS() << "MaximizeHitWeight at " << rowIndex << hitIndex << " with " << weight << endl;
987  fData.MaximizeHitWeight( fData.Row( rowIndex ), hitIndex, weight );
988  }
989  }
990  }
991 }
992 
993 void InitTracklets::operator()( int rowIndex )
994 {
995  if ( ISUNLIKELY( rowIndex >= fTracker.Param().NRows() ) ) {
996  return;
997  }
998  debugF() << "InitTracklets(" << rowIndex << ")" << endl;
999  //std::cout<< "InitTracklets(" << rowIndex << ")" << std::endl;
1000  const int rowStep = AliHLTTPCCAParameters::RowStep;
1001  assert( rowIndex < fTracker.Param().NRows() - rowStep );
1002  // the first hit is a special case (easy)
1003  const ushort_m &mask = rowIndex == r.fStartRow;
1004  const sfloat_m maskF( mask );
1005  {
1006  const AliHLTTPCCARow &row = fData.Row( rowIndex );
1007  const sfloat_v x = fData.RowX( rowIndex );
1008  const ushort_v &hitIndex = static_cast<ushort_v>( r.fCurrentHitIndex );
1009  const sfloat_v y( fData.HitDataY( row ), hitIndex, mask );
1010  const sfloat_v z( fData.HitDataZ( row ), hitIndex, mask );
1011  r.fParam.SetX( x, maskF );
1012  r.fParam.SetY( y, maskF );
1013  r.fParam.SetZ( z, maskF );
1014  fTrackletVector.SetRowHits( rowIndex, trackIndex, hitIndex, mask );
1015 
1016  // mark first hit as used
1017  const short_v isUsed(fData.HitDataIsUsed( row ), static_cast<ushort_v>( r.fCurrentHitIndex ), mask);
1018  fData.SetHitAsUsedInTrackFit( row, static_cast<ushort_v>( r.fCurrentHitIndex ), static_cast<short_m>(mask) && ( isUsed == short_v(1) ) );
1019 
1020  r.fCurrentHitIndex.gather( fData.HitLinkUpData( row ), hitIndex, mask ); // set to next linked hit
1021  // the first hit in the Tracklet is guaranteed to have a link up, since StartHitsFinder
1022  // ensures it
1023  assert( ( r.fCurrentHitIndex >= 0 && mask ) == mask );
1024  }
1025 
1026  rowIndex += rowStep;
1027  {
1028  const AliHLTTPCCARow &row = fData.Row( rowIndex );
1029  const sfloat_v x = fData.RowX( rowIndex );
1030  const ushort_v &hitIndex = static_cast<ushort_v>( r.fCurrentHitIndex );
1031  const sfloat_v y( fData.HitDataY( row ), hitIndex, mask);
1032  const sfloat_v z( fData.HitDataZ( row ), hitIndex, mask);
1033  const sfloat_v &dx = x - r.fParam.X();
1034  const sfloat_v &dy = y - r.fParam.Y();
1035  const sfloat_v &dz = z - r.fParam.Z();
1036  const sfloat_v &ri = sfloat_v( Vc::One ) * CAMath::RSqrt( dx * dx + dy * dy );
1037  const sfloat_v &sinPhi = dy * ri;
1038  r.fParam.SetSinPhi( sinPhi, maskF );
1039  r.fParam.SetDzDs ( dz * ri, maskF );
1040  sfloat_v err2Y, err2Z;
1042 // fTracker.GetErrors2( rowIndex, r.fParam, &err2Y, &err2Z );
1043  fTracker.GetErrors2( rowIndex, r.fParam, &err2Y, &err2Z );
1045  r.fParam.SetCov( 0, err2Y, maskF );
1046  r.fParam.SetCov( 2, err2Z, maskF );
1047 
1048  const sfloat_m transported = r.fParam.TransportToX( x, sinPhi, fTracker.Param().cBz(), -1.f, maskF );
1049  // assert( transported == maskF );
1051 // fTracker.GetErrors2( rowIndex, r.fParam.GetZ(), sinPhi, r.fParam.GetDzDs(), &err2Y, &err2Z );
1052  fTracker.GetErrors2( rowIndex, r.fParam, &err2Y, &err2Z );
1054  const short_m hitAdded( r.fParam.Filter( maskF, y, z, err2Y, err2Z, .99f ) );
1055  // assert( hitAdded == mask );
1056  UNUSED_PARAM2( transported, hitAdded );
1057  fTrackletVector.SetRowHits( rowIndex, trackIndex, hitIndex, mask );
1058 
1059  r.fLastY( maskF ) = y;
1060  r.fLastZ( maskF ) = z;
1061 
1062  // mark 2-nd hit as used
1063  const short_v isUsed(fData.HitDataIsUsed( row ), static_cast<ushort_v>( r.fCurrentHitIndex ), mask);
1064  fData.SetHitAsUsedInTrackFit( row, static_cast<ushort_v>( r.fCurrentHitIndex ), static_cast<short_m>(mask) && ( isUsed == short_v(1) ) );
1065 
1066  r.fCurrentHitIndex.gather( fData.HitLinkUpData( row ), hitIndex, mask ); // set to next linked hit
1067 
1068  // the second hit in the Tracklet is also guaranteed to have a link up, since StartHitsFinder
1069  // ensures it
1070  assert( ( r.fCurrentHitIndex >= 0 && mask ) == mask );
1071  }
1072 }
1073 
AliHLTTPCCASliceData::FirstHitInBin
ushort_v FirstHitInBin(const AliHLTTPCCARow &row, ushort_v binIndexes) const
Definition: AliHLTTPCCASliceDataVector.h:371
AliHLTTPCCASliceData::Row
const AliHLTTPCCARow & Row(int rowIndex) const
Definition: AliHLTTPCCASliceDataVector.h:393
AliHLTTPCCAGrid::GetBinBounded
ushort_v GetBinBounded(const sfloat_v &Y, const sfloat_v &Z) const
Definition: AliHLTTPCCAGrid.h:115
AliHLTTPCCASliceData::CalculateHitWeight
static ushort_v CalculateHitWeight(ushort_v numberOfHits, ushort_v unique)
Definition: AliHLTTPCCASliceDataVector.h:398
AliHLTTPCCASliceData::MaximizeHitWeight
void MaximizeHitWeight(const AliHLTTPCCARow &row, const ushort_v &hitIndex, const ushort_v &weight)
Definition: AliHLTTPCCASliceDataVector.h:414
Generated by   doxygen 1.8.5