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