output of the slice tracker optimised
[u/mrichter/AliRoot.git] / HLT / TPCLib / tracking-ca / AliHLTTPCCAMerger.cxx
1 // $Id: AliHLTTPCCAMerger.cxx 30732 2009-01-22 23:02:02Z sgorbuno $
2 // **************************************************************************
3 // This file is property of and copyright by the ALICE HLT Project          *
4 // ALICE Experiment at CERN, All rights reserved.                           *
5 //                                                                          *
6 // Primary Authors: Sergey Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de> *
7 //                  Ivan Kisel <kisel@kip.uni-heidelberg.de>                *
8 //                  for The ALICE HLT Project.                              *
9 //                                                                          *
10 // Permission to use, copy, modify and distribute this software and its     *
11 // documentation strictly for non-commercial purposes is hereby granted     *
12 // without fee, provided that the above copyright notice appears in all     *
13 // copies and that both the copyright notice and this permission notice     *
14 // appear in the supporting documentation. The authors make no claims       *
15 // about the suitability of this software for any purpose. It is            *
16 // provided "as is" without express or implied warranty.                    *
17 //                                                                          *
18 //***************************************************************************
19
20 #include <stdio.h>
21 #include "AliHLTTPCCASliceOutTrack.h"
22 #include "AliHLTTPCCATracker.h"
23 #include "AliHLTTPCCATrackParam.h"
24
25 #include "AliHLTTPCCAMerger.h"
26
27 #include "AliHLTTPCCAMath.h"
28 #include "TStopwatch.h"
29
30 #include "AliHLTTPCCATrackParam.h"
31 #include "AliHLTTPCCASliceOutput.h"
32 #include "AliHLTTPCCAMergedTrack.h"
33 #include "AliHLTTPCCAMergerOutput.h"
34 #include "AliHLTTPCCADataCompressor.h"
35 #include "AliHLTTPCCAParam.h"
36 #include "AliHLTTPCCATrackLinearisation.h"
37 #include "AliHLTTPCCADataCompressor.h"
38
39 class AliHLTTPCCAMerger::AliHLTTPCCASliceTrackInfo
40 {
41
42   public:
43
44     const AliHLTTPCCATrackParam &InnerParam() const { return fInnerParam;      }
45     const AliHLTTPCCATrackParam &OuterParam() const { return fOuterParam;      }
46     float InnerAlpha()                      const { return fInnerAlpha;      }
47     float OuterAlpha()                      const { return fOuterAlpha;      }
48     int   NClusters()                       const { return fNClusters;       }
49     int   FirstClusterRef()                 const { return fFirstClusterRef; }
50     int   PrevNeighbour()                   const { return fPrevNeighbour;   }
51     int   NextNeighbour()                   const { return fNextNeighbour;   }
52     bool  Used()                            const { return fUsed;            }
53
54     void SetInnerParam( const AliHLTTPCCATrackParam &v ) { fInnerParam = v;      }
55     void SetOuterParam( const AliHLTTPCCATrackParam &v ) { fOuterParam = v;      }
56     void SetInnerAlpha( float v )                      { fInnerAlpha = v;      }
57     void SetOuterAlpha( float v )                      { fOuterAlpha = v;      }
58     void SetNClusters ( int v )                        { fNClusters = v;       }
59     void SetFirstClusterRef( int v )                   { fFirstClusterRef = v; }
60     void SetPrevNeighbour( int v )                     { fPrevNeighbour = v;   }
61     void SetNextNeighbour( int v )                     { fNextNeighbour = v;   }
62     void SetUsed( bool v )                             { fUsed = v;            }
63
64   private:
65
66     AliHLTTPCCATrackParam fInnerParam; // inner parameters
67     AliHLTTPCCATrackParam fOuterParam; // outer parameters
68     float fInnerAlpha;               // alpha angle for inner parameters
69     float fOuterAlpha;               // alpha angle for outer parameters
70     int fNClusters;                  // N clusters
71     int fFirstClusterRef;  // index of the first track cluster in the global cluster array
72     int fPrevNeighbour;    // neighbour in the previous slise
73     int fNextNeighbour;    // neighbour in the next slise
74     bool fUsed;            // is the slice track already merged
75
76 };
77
78
79 class AliHLTTPCCAMerger::AliHLTTPCCABorderTrack
80 {
81
82   public:
83
84     const AliHLTTPCCATrackParam &Param() const { return fParam;     }
85     int   TrackID()                    const { return fTrackID;   }
86     int   NClusters()                  const { return fNClusters; }
87     int   IRow()                       const { return fIRow;      }
88     float X()                          const { return fX;         }
89     bool  OK()                         const { return fOK;        }
90
91     void SetParam     ( const AliHLTTPCCATrackParam &v ) { fParam     = v; }
92     void SetTrackID   ( int v )                        { fTrackID   = v; }
93     void SetNClusters ( int v )                        { fNClusters = v; }
94     void SetIRow      ( int v )                        { fIRow      = v; }
95     void SetX         ( float v )                      { fX         = v; }
96     void SetOK        ( bool v )                       { fOK        = v; }
97
98   private:
99
100     AliHLTTPCCATrackParam fParam;  // track parameters at the border
101     int   fTrackID;              // track index
102     int   fNClusters;            // n clusters
103     int   fIRow;                 // row number of the closest cluster
104     float fX;                    // X coordinate of the closest cluster
105     bool  fOK;                   // is the track rotated and extrapolated correctly
106
107 };
108
109
110
111 AliHLTTPCCAMerger::AliHLTTPCCAMerger()
112     :
113     fSliceParam(),
114     fOutput( 0 ),
115     fTrackInfos( 0 ),
116     fMaxTrackInfos( 0 ),
117     fClusterInfos( 0 ),
118     fMaxClusterInfos( 0 )
119 {
120   //* constructor
121   Clear();
122 }
123
124 /*
125 AliHLTTPCCAMerger::AliHLTTPCCAMerger(const AliHLTTPCCAMerger&)
126   :
127   fSliceParam(),
128   fkSlices(0),
129   fOutput(0),
130   fTrackInfos(0),
131   fMaxTrackInfos(0),
132   fClusterInfos(0),
133   fMaxClusterInfos(0)
134 {
135 }
136
137 const AliHLTTPCCAMerger &AliHLTTPCCAMerger::operator=(const AliHLTTPCCAMerger&) const
138 {
139   return *this;
140 }
141 */
142
143 AliHLTTPCCAMerger::~AliHLTTPCCAMerger()
144 {
145   //* destructor
146   if ( fTrackInfos ) delete[] fTrackInfos;
147   if ( fClusterInfos ) delete[] fClusterInfos;
148   if ( fOutput ) delete[] ( ( float2* )( fOutput ) );
149 }
150
151 void AliHLTTPCCAMerger::Clear()
152 {
153   for ( int i = 0; i < fgkNSlices; ++i ) {
154     fkSlices[i] = 0;
155     fSliceNTrackInfos[ i ] = 0;
156     fSliceTrackInfoStart[ i ] = 0;
157   }
158   if ( fOutput ) delete[] ( ( float2* )( fOutput ) );
159   if ( fTrackInfos ) delete[] fTrackInfos;
160   if ( fClusterInfos ) delete[] fClusterInfos;
161   fOutput = 0;
162   fTrackInfos = 0;
163   fClusterInfos = 0;
164   fMaxTrackInfos = 0;
165   fMaxClusterInfos = 0;
166 }
167
168
169 void AliHLTTPCCAMerger::SetSliceData( int index, const AliHLTTPCCASliceOutput *SliceData )
170 {
171   fkSlices[index] = SliceData;
172 }
173
174 void AliHLTTPCCAMerger::Reconstruct()
175 {
176   //* main merging routine
177
178   UnpackSlices();
179   Merging();
180 }
181
182 void AliHLTTPCCAMerger::UnpackSlices()
183 {
184   //* unpack the cluster information from the slice tracks and initialize track info array
185
186   // get N tracks and N clusters in event
187
188   int nTracksTotal = 0;
189   int nTrackClustersTotal = 0;
190   for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
191     if ( !fkSlices[iSlice] ) continue;
192     nTracksTotal += fkSlices[iSlice]->NTracks();
193     nTrackClustersTotal += fkSlices[iSlice]->NTrackClusters();
194   }
195
196   // book/clean memory if necessary
197   {
198     if ( nTracksTotal > fMaxTrackInfos || ( fMaxTrackInfos > 100 && nTracksTotal < 0.5*fMaxTrackInfos ) ) {
199       if ( fTrackInfos ) delete[] fTrackInfos;
200       fMaxTrackInfos = ( int ) ( nTracksTotal * 1.2 );
201       fTrackInfos = new AliHLTTPCCASliceTrackInfo[fMaxTrackInfos];
202     }
203
204     if ( nTrackClustersTotal > fMaxClusterInfos || ( fMaxClusterInfos > 1000 && nTrackClustersTotal < 0.5*fMaxClusterInfos ) ) {
205       if ( fClusterInfos ) delete[] fClusterInfos;
206       fMaxClusterInfos = ( int ) ( nTrackClustersTotal * 1.2 );
207       fClusterInfos = new AliHLTTPCCAClusterInfo [fMaxClusterInfos];
208     }
209
210     if ( fOutput ) delete[] ( ( float2* )( fOutput ) );
211     int size = fOutput->EstimateSize( nTracksTotal, nTrackClustersTotal );
212     fOutput = ( AliHLTTPCCAMergerOutput* )( new float2[size/sizeof( float2 )+1] );
213   }
214
215   // unpack track and cluster information
216
217   int nTracksCurrent = 0;
218   int nClustersCurrent = 0;
219
220   for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
221
222     fSliceTrackInfoStart[ iSlice ] = nTracksCurrent;
223     fSliceNTrackInfos[ iSlice ] = 0;
224
225     if ( !fkSlices[iSlice] ) continue;
226
227     const AliHLTTPCCASliceOutput &slice = *( fkSlices[iSlice] );
228
229     const AliHLTTPCCASliceOutTrack *sliceTr = slice.GetFirstTrack();
230     
231     for ( int itr = 0; itr < slice.NTracks(); itr++ ) {
232
233       AliHLTTPCCATrackParam t0;
234       t0.InitParam();
235       t0.SetParam(sliceTr->Param());
236       int nCluNew = 0;
237       int nCluOld = sliceTr->NClusters();
238       for ( int iTrClu = 0; iTrClu < nCluOld; iTrClu++ ) {
239         
240         // unpack cluster information
241
242         AliHLTTPCCAClusterInfo &clu = fClusterInfos[nClustersCurrent + nCluNew];
243         int id, row;
244         float x,y,z;
245         sliceTr->Cluster( iTrClu ).Get(id,row,x,y,z);
246         
247         clu.SetISlice( iSlice );
248         clu.SetIRow( row );
249         clu.SetId( id );
250         clu.SetPackedAmp( 0 );
251         clu.SetX( x );
252         clu.SetY( y );
253         clu.SetZ( z );
254
255         if ( !t0.TransportToX( clu.X(), fSliceParam.GetBz( t0 ), .999 ) ) continue;
256
257         float err2Y, err2Z;
258         fSliceParam.GetClusterErrors2( clu.IRow(), clu.Z(), t0.SinPhi(), t0.GetCosPhi(), t0.DzDs(), err2Y, err2Z );
259
260         clu.SetErr2Y( err2Y );
261         clu.SetErr2Z( err2Z );
262         nCluNew++ ;
263       }
264
265       // refit the track
266
267       int hits[1000];
268       int nHits = nCluNew;
269       for ( int i = 0; i < nHits; i++ ) hits[i] = nClustersCurrent + i;
270
271       AliHLTTPCCATrackParam startPoint;
272       startPoint.InitParam();
273       startPoint.SetParam(sliceTr->Param());
274       AliHLTTPCCATrackParam endPoint = startPoint;
275       float startAlpha = fSliceParam.Alpha( iSlice );
276       float endAlpha = startAlpha;
277   
278       sliceTr = sliceTr->GetNextTrack();
279       
280       if ( nCluNew < .8*nCluOld ) continue;
281
282       if ( !FitTrack( endPoint, endAlpha, startPoint, startAlpha, hits, nHits, 0 ) ) continue;
283
284       startPoint = endPoint;
285       startAlpha = endAlpha;
286       if ( !FitTrack( startPoint, startAlpha, endPoint, endAlpha, hits, nHits, 1 ) ) continue;
287
288       if ( nHits < .8*nCluOld ) continue;
289
290       // store the track
291
292       AliHLTTPCCASliceTrackInfo &track = fTrackInfos[nTracksCurrent];
293
294       track.SetInnerParam( startPoint );
295       track.SetInnerAlpha( startAlpha );
296       track.SetOuterParam( endPoint );
297       track.SetOuterAlpha( endAlpha );
298       track.SetFirstClusterRef( nClustersCurrent );
299       track.SetNClusters( nHits );
300       track.SetPrevNeighbour( -1 );
301       track.SetNextNeighbour( -1 );
302       track.SetUsed( 0 );
303
304       for ( int i = 0; i < nHits; i++ )
305         fClusterInfos[nClustersCurrent + i] = fClusterInfos[hits[i]];
306       nTracksCurrent++;
307       fSliceNTrackInfos[ iSlice ]++;
308       nClustersCurrent += nHits;
309     }
310     //std::cout<<"Unpack slice "<<iSlice<<": ntracks "<<slice.NTracks()<<"/"<<fSliceNTrackInfos[iSlice]<<std::endl;
311   }
312 }
313
314
315
316 bool AliHLTTPCCAMerger::FitTrack( AliHLTTPCCATrackParam &T, float &Alpha,
317                                   AliHLTTPCCATrackParam t0, float Alpha0,
318                                   int hits[], int &NTrackHits, bool dir, bool final,
319                                   AliHLTTPCCAClusterInfo *infoArray )
320 {
321   // Fit the track
322
323   AliHLTTPCCATrackParam::AliHLTTPCCATrackFitParam fitPar;
324   AliHLTTPCCATrackParam t = t0;
325   AliHLTTPCCATrackLinearisation l( t0 );
326
327   bool first = 1;
328   bool doErrors = 1;
329   if ( !infoArray ) {
330     infoArray = fClusterInfos;
331     doErrors = 0;
332   }
333
334   t.CalculateFitParameters( fitPar );
335
336   int hitsNew[1000];
337   int nHitsNew = 0;
338
339   for ( int ihit = 0; ihit < NTrackHits; ihit++ ) {
340
341     int jhit = dir ? ( NTrackHits - 1 - ihit ) : ihit;
342     AliHLTTPCCAClusterInfo &h = infoArray[hits[jhit]];
343
344     int iSlice = h.ISlice();
345
346     float sliceAlpha =  fSliceParam.Alpha( iSlice );
347
348     if ( CAMath::Abs( sliceAlpha - Alpha0 ) > 1.e-4 ) {
349       if ( ! t.Rotate(  sliceAlpha - Alpha0, l, .999 ) ) continue;
350       Alpha0 = sliceAlpha;
351     }
352
353     //float x = fSliceParam.RowX( h.IRow() );
354     float x = h.X();
355
356     if ( !t.TransportToXWithMaterial( x, l, fitPar, fSliceParam.GetBz( t ) ) ) continue;
357
358     if ( first ) {
359       t.SetCov( 0, 10 );
360       t.SetCov( 1,  0 );
361       t.SetCov( 2, 10 );
362       t.SetCov( 3,  0 );
363       t.SetCov( 4,  0 );
364       t.SetCov( 5,  1 );
365       t.SetCov( 6,  0 );
366       t.SetCov( 7,  0 );
367       t.SetCov( 8,  0 );
368       t.SetCov( 9,  1 );
369       t.SetCov( 10,  0 );
370       t.SetCov( 11,  0 );
371       t.SetCov( 12,  0 );
372       t.SetCov( 13,  0 );
373       t.SetCov( 14,  10 );
374       t.SetChi2( 0 );
375       t.SetNDF( -5 );
376       t.CalculateFitParameters( fitPar );
377     }
378
379     float err2Y = h.Err2Y();
380     float err2Z = h.Err2Z();
381     if ( doErrors ) fSliceParam.GetClusterErrors2( h.IRow(), h.Z(), l.SinPhi(), l.CosPhi(), l.DzDs(), err2Y, err2Z );
382     if( !final ){
383       err2Y*= fSliceParam.ClusterError2CorrectionY();
384       err2Z*= fSliceParam.ClusterError2CorrectionZ();
385     }
386
387     if ( !t.Filter( h.Y(), h.Z(), err2Y, err2Z ) ) continue;
388
389     first = 0;
390
391     hitsNew[nHitsNew++] = hits[jhit];
392   }
393
394   if ( CAMath::Abs( t.QPt() ) < 1.e-4 ) t.SetQPt( 1.e-4 );
395
396   bool ok = t.CheckNumericalQuality();
397
398   if ( CAMath::Abs( t.SinPhi() ) > .99 ) ok = 0;
399   else if ( l.CosPhi() >= 0 ) t.SetSignCosPhi( 1 );
400   else t.SetSignCosPhi( -1 );
401
402   if ( ok ) {
403     T = t;
404     Alpha = Alpha0;
405     NTrackHits = nHitsNew;
406     for ( int i = 0; i < NTrackHits; i++ ) {
407       hits[dir ?( NTrackHits-1-i ) :i] = hitsNew[i];
408     }
409   }
410   return ok;
411 }
412
413
414 float AliHLTTPCCAMerger::GetChi2( float x1, float y1, float a00, float a10, float a11,
415                                   float x2, float y2, float b00, float b10, float b11  )
416 {
417   //* Calculate Chi2/ndf deviation
418
419   float d[2] = { x1 - x2, y1 - y2 };
420
421   float mSi[3] = { a00 + b00, a10 + b10, a11 + b11 };
422
423   float s = ( mSi[0] * mSi[2] - mSi[1] * mSi[1] );
424
425   if ( s < 1.E-10 ) return 10000.;
426
427   float mS[3] = { mSi[2], -mSi[1], mSi[0] };
428
429   return AliHLTTPCCAMath::Abs( ( ( mS[0]*d[0] + mS[1]*d[1] )*d[0]
430                        + ( mS[1]*d[0] + mS[2]*d[1] )*d[1] ) / s / 2 );
431
432 }
433
434
435
436 void AliHLTTPCCAMerger::MakeBorderTracks( int iSlice, int iBorder, AliHLTTPCCABorderTrack B[], int &nB )
437 {
438   //* prepare slice tracks for merging with next/previous/same sector
439   //* each track transported to the border line,
440   //* in some cases both inner and outer parameters of the track are transported
441
442   static int statAll = 0, statOK = 0;
443   nB = 0;
444   float dAlpha = fSliceParam.DAlpha() / 2;
445   float x0 = 0;
446
447   if ( iBorder == 0 ) { // transport to the left age of the sector and rotate horisontally
448     dAlpha = dAlpha - CAMath::Pi() / 2 ;
449   } else if ( iBorder == 1 ) { //  transport to the right age of the sector and rotate horisontally
450     dAlpha = -dAlpha - CAMath::Pi() / 2 ;
451   } else if ( iBorder == 2 ) { // transport to the left age of the sector and rotate vertically
452     dAlpha = dAlpha;
453     x0 = fSliceParam.RowX( 63 );
454   } else if ( iBorder == 3 ) { // transport to the right age of the sector and rotate vertically
455     dAlpha = -dAlpha;
456     x0 =  fSliceParam.RowX( 63 );
457   } else if ( iBorder == 4 ) { // transport to the middle of the sector, w/o rotation
458     dAlpha = 0;
459     x0 = fSliceParam.RowX( 63 );
460   }
461
462   for ( int itr = 0; itr < fSliceNTrackInfos[iSlice]; itr++ ) {
463
464     const AliHLTTPCCASliceTrackInfo &track = fTrackInfos[ fSliceTrackInfoStart[iSlice] + itr ];
465
466     AliHLTTPCCATrackParam t0 = track.InnerParam();
467     AliHLTTPCCATrackParam t1 = track.OuterParam();
468
469     const float maxSin = CAMath::Sin( 60. / 180.*CAMath::Pi() );
470
471     bool ok0 = t0.Rotate( dAlpha, maxSin );
472     bool ok1 = t1.Rotate( dAlpha, maxSin );
473
474     bool do0 = ok0;
475     bool do1 = ok1 && ( !ok0 || t1.SignCosPhi() * t0.SignCosPhi() < 0 );
476
477     if ( ok0 && !do1 && ok1 && ( t1.X() < t0.X() ) ) {
478       do0 = 0;
479       do1 = 1;
480     }
481
482     if ( do0 ) {
483       AliHLTTPCCABorderTrack &b = B[nB];
484       b.SetX( t0.GetX() );
485       
486       if ( t0.TransportToX( x0, fSliceParam.GetBz( t0 ), maxSin ) ) {
487         b.SetOK( 1 );
488         b.SetTrackID( itr );
489         b.SetNClusters( track.NClusters() );
490         b.SetIRow( fClusterInfos[ track.FirstClusterRef() + 0 ].IRow() );
491         b.SetParam( t0 );
492         nB++;
493       }
494     }
495     if ( do1 ) {
496       AliHLTTPCCABorderTrack &b = B[nB];
497       b.SetX( t1.GetX() );
498       
499       if ( t1.TransportToX( x0, fSliceParam.GetBz( t1 ), maxSin ) ) {
500         b.SetOK( 1 );
501         b.SetTrackID( itr );
502         b.SetNClusters( track.NClusters() );
503         b.SetIRow( fClusterInfos[ track.FirstClusterRef() + track.NClusters()-1 ].IRow() );
504         b.SetParam( t1 );
505         nB++;
506       }
507     }
508     if ( do0 || do1 ) statOK++;
509     statAll++;
510   }
511 }
512
513
514
515 void AliHLTTPCCAMerger::MergeBorderTracks( int iSlice1, AliHLTTPCCABorderTrack B1[], int N1,
516     int iSlice2, AliHLTTPCCABorderTrack B2[], int N2
517                                          )
518 {
519   //* merge two sets of tracks
520
521   //std::cout<<" Merge slices "<<iSlice1<<"+"<<iSlice2<<": tracks "<<N1<<"+"<<N2<<std::endl;
522   
523   float factor2ys = 1.;//1.5;//SG!!!
524   float factor2zt = 1.;//1.5;//SG!!!
525   float factor2k = 2.0;//2.2;
526
527   factor2k  = 3.5 * 3.5 * factor2k * factor2k;
528   factor2ys = 3.5 * 3.5 * factor2ys * factor2ys;
529   factor2zt = 3.5 * 3.5 * factor2zt * factor2zt;
530
531   int minNPartHits = 10;//SG!!!
532   int minNTotalHits = 20;
533
534   //float maxDX = fSliceParam.RowX(40) -  fSliceParam.RowX(0);
535
536   for ( int i1 = 0; i1 < N1; i1++ ) {
537     AliHLTTPCCABorderTrack &b1 = B1[i1];
538     if ( !b1.OK() ) continue;
539     if ( b1.NClusters() < minNPartHits ) continue;
540     const AliHLTTPCCATrackParam &t1 = b1.Param();
541     int iBest2 = -1;
542     int lBest2 = 0;
543     int start2 = ( iSlice1 != iSlice2 ) ? 0 : i1 + 1;
544     for ( int i2 = start2; i2 < N2; i2++ ) {
545       AliHLTTPCCABorderTrack &b2 = B2[i2];
546       if ( !b2.OK() ) continue;
547       if ( b2.NClusters() < minNPartHits ) continue;
548       if ( b2.NClusters() < lBest2 ) continue;
549       if ( b1.NClusters() + b2.NClusters() < minNTotalHits ) continue;
550
551       //if( TMath::Abs(b1.fX - b2.fX)>maxDX ) continue;
552
553       const AliHLTTPCCATrackParam &t2 = b2.Param();
554
555       float c = t2.SignCosPhi() * t1.SignCosPhi() >= 0 ? 1 : -1;
556       float dk = t2.QPt() - c * t1.QPt();
557       float s2k = t2.Err2QPt() + t1.Err2QPt();
558       //std::cout<<" check 1.. "<<dk/sqrt(factor2k)<<std::endl;
559       if ( dk*dk > factor2k*s2k ) continue;
560
561
562       float chi2ys = GetChi2( t1.Y(), c * t1.SinPhi(), t1.Cov()[0], c * t1.Cov()[3], t1.Cov()[5],
563                               t2.Y(),  t2.SinPhi(), t2.Cov()[0],  t2.Cov()[3], t2.Cov()[5] );
564
565       //std::cout<<" check 2.. "<<sqrt(chi2ys/factor2ys)<<std::endl;
566       if ( chi2ys > factor2ys ) continue;
567       
568
569       float chi2zt = GetChi2( t1.Z(), c * t1.DzDs(), t1.Cov()[2], c * t1.Cov()[7], t1.Cov()[9],
570                               t2.Z(),  t2.DzDs(), t2.Cov()[2],  t2.Cov()[7], t2.Cov()[9] );
571
572       //std::cout<<" check 3.. "<<sqrt(chi2zt/factor2zt)<<std::endl;
573       if ( chi2zt > factor2zt ) continue;
574
575       lBest2 = b2.NClusters();
576       iBest2 = b2.TrackID();
577     }
578
579     if ( iBest2 < 0 ) continue;
580
581     //std::cout<<"Neighbour found for "<<i1<<": "<<iBest2<<std::endl;
582     AliHLTTPCCASliceTrackInfo &newTrack1 = fTrackInfos[fSliceTrackInfoStart[iSlice1] + b1.TrackID() ];
583     AliHLTTPCCASliceTrackInfo &newTrack2 = fTrackInfos[fSliceTrackInfoStart[iSlice2] + iBest2 ];
584
585     int old1 = newTrack2.PrevNeighbour();
586
587     if ( old1 >= 0 ) {
588       AliHLTTPCCASliceTrackInfo &oldTrack1 = fTrackInfos[fSliceTrackInfoStart[iSlice1] + old1];
589       if ( oldTrack1.NClusters()  < newTrack1.NClusters() ) {
590         newTrack2.SetPrevNeighbour( -1 );
591         oldTrack1.SetNextNeighbour( -1 );
592       } else continue;
593     }
594     int old2 = newTrack1.NextNeighbour();
595     if ( old2 >= 0 ) {
596       AliHLTTPCCASliceTrackInfo &oldTrack2 = fTrackInfos[fSliceTrackInfoStart[iSlice2] + old2];
597       if ( oldTrack2.NClusters() < newTrack2.NClusters() ) {
598         oldTrack2.SetPrevNeighbour( -1 );
599       } else continue;
600     }
601     newTrack1.SetNextNeighbour( iBest2 );
602     newTrack2.SetPrevNeighbour( b1.TrackID() );
603     //std::cout<<"Neighbourhood is set"<<std::endl;
604   }
605
606 }
607
608
609 void AliHLTTPCCAMerger::Merging()
610 {
611   //* track merging between slices
612
613   fOutput->SetNTracks( 0 );
614   fOutput->SetNTrackClusters( 0 );
615   fOutput->SetPointers();
616
617
618   // for each slice set number of the next neighbouring slice
619
620   int nextSlice[100], prevSlice[100];
621
622   for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
623     nextSlice[iSlice] = iSlice + 1;
624     prevSlice[iSlice] = iSlice - 1;
625   }
626   int mid = fgkNSlices / 2 - 1 ;
627   int last = fgkNSlices - 1 ;
628   if ( mid < 0 ) mid = 0; // to avoid compiler warning
629   if ( last < 0 ) last = 0; //
630   nextSlice[ mid ] = 0;
631   prevSlice[ 0 ] = mid;
632   nextSlice[ last ] = fgkNSlices / 2;
633   prevSlice[ fgkNSlices/2 ] = last;
634
635   int maxNSliceTracks = 0;
636   for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
637     if ( maxNSliceTracks < fSliceNTrackInfos[iSlice] ) maxNSliceTracks = fSliceNTrackInfos[iSlice];
638   }
639
640   if ( 1 ) {// merging track segments withing one slice
641
642     AliHLTResizableArray<AliHLTTPCCABorderTrack> bord( maxNSliceTracks*2 );
643
644     AliHLTTPCCASliceTrackInfo *tmpT = new AliHLTTPCCASliceTrackInfo[maxNSliceTracks];
645     AliHLTTPCCAClusterInfo *tmpH = new AliHLTTPCCAClusterInfo[fMaxClusterInfos];
646
647     for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
648
649       int nBord = 0;
650       MakeBorderTracks( iSlice, 4, bord.Data(), nBord );
651       MergeBorderTracks( iSlice, bord.Data(), nBord, iSlice, bord.Data(), nBord );
652
653       int nTr = 0, nH = 0;
654       int sliceFirstClusterRef = 0;
655       for ( int itr = 0; itr < fSliceNTrackInfos[iSlice]; itr++ ) {
656         AliHLTTPCCASliceTrackInfo &track = fTrackInfos[ fSliceTrackInfoStart[iSlice] + itr];
657         if ( itr == 0 ) sliceFirstClusterRef = track.FirstClusterRef();
658         track.SetPrevNeighbour( -1 );
659         if ( track.NextNeighbour() == -2 ) {
660           track.SetNextNeighbour( -1 );
661           continue;
662         }
663         AliHLTTPCCASliceTrackInfo &trackNew = tmpT[nTr];
664         trackNew = track;
665         trackNew.SetFirstClusterRef( sliceFirstClusterRef + nH );
666
667         for ( int ih = 0; ih < track.NClusters(); ih++ ) tmpH[nH+ih] = fClusterInfos[track.FirstClusterRef()+ih];
668         nTr++;
669         nH += track.NClusters();
670
671         int jtr =  track.NextNeighbour();
672
673         if ( jtr < 0 ) continue;
674         AliHLTTPCCASliceTrackInfo &neighTrack = fTrackInfos[ fSliceTrackInfoStart[iSlice] + jtr];
675
676         track.SetNextNeighbour( -1 );
677         neighTrack.SetNextNeighbour( -2 );
678
679         for ( int ih = 0; ih < neighTrack.NClusters(); ih++ )
680           tmpH[nH+ih] = fClusterInfos[neighTrack.FirstClusterRef()+ih];
681
682         trackNew.SetNClusters( trackNew.NClusters() + neighTrack.NClusters() );
683         trackNew.SetNextNeighbour( -1 );
684         nH += neighTrack.NClusters();
685         if ( neighTrack.InnerParam().X() < track.InnerParam().X() ) trackNew.SetInnerParam( neighTrack.InnerParam() );
686         if ( neighTrack.OuterParam().X() > track.OuterParam().X() ) trackNew.SetOuterParam( neighTrack.OuterParam() );
687       }
688
689       fSliceNTrackInfos[iSlice] = nTr;
690       for ( int itr = 0; itr < nTr; itr++ ) fTrackInfos[ fSliceTrackInfoStart[iSlice] + itr] = tmpT[itr];
691       for ( int ih = 0; ih < nH; ih++ ) fClusterInfos[sliceFirstClusterRef + ih] = tmpH[ih];
692
693     }
694     delete[] tmpT;
695     delete[] tmpH;
696   }
697
698
699   //* merging tracks between slices
700
701
702   // arrays for the rotated track parameters
703
704   AliHLTTPCCABorderTrack
705   *bCurr0 = new AliHLTTPCCABorderTrack[maxNSliceTracks*2],
706   *bNext0 = new AliHLTTPCCABorderTrack[maxNSliceTracks*2],
707   *bCurr = new AliHLTTPCCABorderTrack[maxNSliceTracks*2],
708   *bNext = new AliHLTTPCCABorderTrack[maxNSliceTracks*2];
709
710   for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
711
712     int jSlice = nextSlice[iSlice];
713
714     int nCurr0 = 0, nNext0 = 0;
715     int nCurr = 0, nNext = 0;
716
717     MakeBorderTracks( iSlice, 0, bCurr, nCurr );
718     MakeBorderTracks( jSlice, 1, bNext, nNext );
719     MakeBorderTracks( iSlice, 2, bCurr0, nCurr0 );
720     MakeBorderTracks( jSlice, 3, bNext0, nNext0 );
721
722     MergeBorderTracks( iSlice, bCurr0, nCurr0, jSlice, bNext0, nNext0 );
723     MergeBorderTracks( iSlice, bCurr, nCurr, jSlice, bNext, nNext );
724   }
725
726   if ( bCurr0 ) delete[] bCurr0;
727   if ( bNext0 ) delete[] bNext0;
728   if ( bCurr  ) delete[] bCurr;
729   if ( bNext  ) delete[] bNext;
730
731
732   //static int nRejected = 0;
733
734   int nOutTracks = 0;
735   int nOutTrackClusters = 0;
736
737   AliHLTTPCCAMergedTrack *outTracks = new AliHLTTPCCAMergedTrack[fMaxTrackInfos];
738   unsigned int   *outClusterId = new unsigned int [fMaxClusterInfos];
739   UChar_t  *outClusterPackedAmp = new UChar_t [fMaxClusterInfos];
740
741   for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
742
743     for ( int itr = 0; itr < fSliceNTrackInfos[iSlice]; itr++ ) {
744
745       AliHLTTPCCASliceTrackInfo &track = fTrackInfos[fSliceTrackInfoStart[iSlice] + itr];
746
747       if ( track.Used() ) continue;
748       if ( track.PrevNeighbour() >= 0 ) continue;
749       //std::cout<<"Merged track candidate, nhits "<<track.NClusters()<<std::endl;
750       AliHLTTPCCATrackParam startPoint = track.InnerParam(), endPoint = track.OuterParam();
751       float startAlpha = track.InnerAlpha(), endAlpha = track.OuterAlpha();
752
753       int hits[2000];
754       int firstHit = 1000;
755       int nHits = 0;
756       int jSlice = iSlice;
757       int jtr = itr;
758
759       {
760         track.SetUsed( 1 );
761         for ( int jhit = 0; jhit < track.NClusters(); jhit++ ) {
762           int id = track.FirstClusterRef() + jhit;
763           hits[firstHit+jhit] = id;
764         }
765         nHits = track.NClusters();
766         jtr = track.NextNeighbour();
767         jSlice = nextSlice[iSlice];
768       }
769
770       while ( jtr >= 0 ) {
771         AliHLTTPCCASliceTrackInfo &segment = fTrackInfos[fSliceTrackInfoStart[jSlice] + jtr];
772         if ( segment.Used() ) break;
773         segment.SetUsed( 1 );
774         bool dir = 0;
775         int startHit = firstHit + nHits;
776         float d00 = startPoint.GetDistXZ2( segment.InnerParam() );
777         float d01 = startPoint.GetDistXZ2( segment.OuterParam() );
778         float d10 = endPoint.GetDistXZ2( segment.InnerParam() );
779         float d11 = endPoint.GetDistXZ2( segment.OuterParam() );
780         if ( d00 <= d01 && d00 <= d10 && d00 <= d11 ) {
781           startPoint = segment.OuterParam();
782           startAlpha = segment.OuterAlpha();
783           dir = 1;
784           firstHit -= segment.NClusters();
785           startHit = firstHit;
786         } else if ( d01 <= d10 && d01 <= d11 ) {
787           startPoint = segment.InnerParam();
788           startAlpha = segment.InnerAlpha();
789           dir = 0;
790           firstHit -= segment.NClusters();
791           startHit = firstHit;
792         } else if ( d10 <= d11 ) {
793           endPoint = segment.OuterParam();
794           endAlpha = segment.OuterAlpha();
795           dir = 0;
796         } else {
797           endPoint = segment.InnerParam();
798           endAlpha = segment.InnerAlpha();
799           dir = 1;
800         }
801
802         for ( int jhit = 0; jhit < segment.NClusters(); jhit++ ) {
803           int id = segment.FirstClusterRef() + jhit;
804           hits[startHit+( dir ?( segment.NClusters()-1-jhit ) :jhit )] = id;
805         }
806         nHits += segment.NClusters();
807         jtr = segment.NextNeighbour();
808         jSlice = nextSlice[jSlice];
809       }
810
811       if ( endPoint.X() < startPoint.X() ) { // swap
812         for ( int i = 0; i < nHits; i++ ) hits[i] = hits[firstHit+nHits-1-i];
813         firstHit = 0;
814       }
815
816       if ( nHits < 30 ) continue;    //SG!!!
817
818       // refit
819
820       // need best t0!!!SG
821
822       endPoint = startPoint;
823
824       if ( !FitTrack( endPoint, endAlpha, startPoint, startAlpha, hits + firstHit, nHits, 0,1 ) ) continue;
825       if ( !FitTrack( startPoint, startAlpha, endPoint, endAlpha, hits + firstHit, nHits, 1,1 ) ) continue;
826       if ( nHits < 30 ) continue;    //SG!!!
827
828       AliHLTTPCCATrackParam p = startPoint;
829       
830       if(0){
831         double xTPC = 83.65; //SG!!!
832         double dAlpha = 0.349066;
833         double ymax = 2.* xTPC * CAMath::Tan( dAlpha / 2. );
834         
835         double dRot = 0;
836         if ( p.TransportToXWithMaterial( xTPC, fSliceParam.GetBz( p ) ) ) {
837           double y = p.GetY();
838           if ( y > ymax ) {         
839             if ( p.Rotate( dAlpha ) ){
840               dRot = dAlpha;
841               p.TransportToXWithMaterial( xTPC, fSliceParam.GetBz( p ) );
842             }
843           } else if( y< -ymax ){
844             if ( p.Rotate( -dAlpha ) ){
845               dRot = -dAlpha;
846               p.TransportToXWithMaterial( xTPC, fSliceParam.GetBz( p ) );
847             }
848           }
849         }
850         
851         if ( -ymax <= p.GetY() && p.GetY() <= ymax && p.CheckNumericalQuality() ){
852           startPoint = p;
853           startAlpha+=dRot;
854         }     
855       }
856
857       if ( !startPoint.CheckNumericalQuality() ) continue;
858
859       AliHLTTPCCAMergedTrack &mergedTrack = outTracks[nOutTracks];
860       mergedTrack.SetNClusters( nHits );
861       mergedTrack.SetFirstClusterRef( nOutTrackClusters );
862       mergedTrack.SetInnerParam( startPoint );
863       mergedTrack.SetInnerAlpha( startAlpha );
864       mergedTrack.SetOuterParam( endPoint );
865       mergedTrack.SetOuterAlpha( endAlpha );
866
867       for ( int i = 0; i < nHits; i++ ) {
868         AliHLTTPCCAClusterInfo &clu = fClusterInfos[hits[firstHit+i]];
869         outClusterId[nOutTrackClusters+i] = clu.Id();
870         outClusterPackedAmp[nOutTrackClusters+i] = clu.PackedAmp();
871       }
872
873       nOutTracks++;
874       nOutTrackClusters += nHits;
875     }
876   }
877
878   fOutput->SetNTracks( nOutTracks );
879   #ifdef HLTCA_STANDALONE
880   printf("Tracks Output: %d\n", nOutTracks);
881   #endif
882   fOutput->SetNTrackClusters( nOutTrackClusters );
883   fOutput->SetPointers();
884
885   for ( int itr = 0; itr < nOutTracks; itr++ ) fOutput->SetTrack( itr, outTracks[itr] );
886
887   for ( int ic = 0; ic < nOutTrackClusters; ic++ ) {
888     fOutput->SetClusterId( ic, outClusterId[ic] );
889     fOutput->SetClusterPackedAmp( ic, outClusterPackedAmp[ic] );
890   }
891
892   delete[] outTracks;
893   delete[] outClusterId;
894   delete[] outClusterPackedAmp;
895 }