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
  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,
 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 ( !t.Filter( h.Y(), h.Z(), err2Y, err2Z ) ) continue;
 366
 367     first = 0;
 368
 369     hitsNew[nHitsNew++] = hits[jhit];
 370   }
 371
 372   if ( CAMath::Abs( t.QPt() ) < 1.e-8 ) t.SetQPt( 1.e-8 );
 373
 374   bool ok = 1;
 375
 376   const float *c = t.Cov();
 377   for ( int i = 0; i < 15; i++ ) ok = ok && finite( c[i] );
 378   for ( int i = 0; i < 5; i++ ) ok = ok && finite( t.Par()[i] );
 379   ok = ok && ( t.GetX() > 50 );
 380
 381   if ( c[0] <= 0 || c[2] <= 0 || c[5] <= 0 || c[9] <= 0 || c[14] <= 0 ) ok = 0;
 382   if ( c[0] > 5. || c[2] > 5. || c[5] > 2. || c[9] > 2 || c[14] > 2. ) ok = 0;
 383
 384   if ( CAMath::Abs( t.SinPhi() ) > .99 ) ok = 0;
 385   else if ( l.CosPhi() >= 0 ) t.SetSignCosPhi( 1 );
 386   else t.SetSignCosPhi( -1 );
 387
 388   if ( ok ) {
 389     T = t;
 390     Alpha = Alpha0;
 391     NTrackHits = nHitsNew;
 392     for ( int i = 0; i < NTrackHits; i++ ) {
 393       hits[dir ?( NTrackHits-1-i ) :i] = hitsNew[i];
 394     }
 395   }
 396   return ok;
 397 }
 398
 399
 400 float AliHLTTPCCAMerger::GetChi2( float x1, float y1, float a00, float a10, float a11,
 401                                   float x2, float y2, float b00, float b10, float b11  )
 402 {
 403   //* Calculate Chi2/ndf deviation
 404
 405   float d[2] = { x1 - x2, y1 - y2 };
 406
 407   float mSi[3] = { a00 + b00, a10 + b10, a11 + b11 };
 408
 409   float s = ( mSi[0] * mSi[2] - mSi[1] * mSi[1] );
 410
 411   if ( s < 1.E-10 ) return 10000.;
 412
 413   float mS[3] = { mSi[2], -mSi[1], mSi[0] };
 414
 415   return TMath::Abs( ( ( mS[0]*d[0] + mS[1]*d[1] )*d[0]
 416                        + ( mS[1]*d[0] + mS[2]*d[1] )*d[1] ) / s / 2 );
 417
 418 }
 419
 420
 421
 422 void AliHLTTPCCAMerger::MakeBorderTracks( int iSlice, int iBorder, AliHLTTPCCABorderTrack B[], int &nB )
 423 {
 424   //* prepare slice tracks for merging with next/previous/same sector
 425   //* each track transported to the border line,
 426   //* in some cases both inner and outer parameters of the track are transported
 427
 428   static int statAll = 0, statOK = 0;
 429   nB = 0;
 430   float dAlpha = fSliceParam.DAlpha() / 2;
 431   float x0 = 0;
 432
 433   if ( iBorder == 0 ) { // transport to the left age of the sector and rotate horisontally
 434     dAlpha = dAlpha - CAMath::Pi() / 2 ;
 435   } else if ( iBorder == 1 ) { //  transport to the right age of the sector and rotate horisontally
 436     dAlpha = -dAlpha - CAMath::Pi() / 2 ;
 437   } else if ( iBorder == 2 ) { // transport to the left age of the sector and rotate vertically
 438     dAlpha = dAlpha;
 439     x0 = fSliceParam.RowX( 63 );
 440   } else if ( iBorder == 3 ) { // transport to the right age of the sector and rotate vertically
 441     dAlpha = -dAlpha;
 442     x0 =  fSliceParam.RowX( 63 );
 443   } else if ( iBorder == 4 ) { // transport to the middle of the sector, w/o rotation
 444     dAlpha = 0;
 445     x0 = fSliceParam.RowX( 63 );
 446   }
 447
 448   for ( int itr = 0; itr < fSliceNTrackInfos[iSlice]; itr++ ) {
 449
 450     const AliHLTTPCCASliceTrackInfo &track = fTrackInfos[ fSliceTrackInfoStart[iSlice] + itr ];
 451
 452     AliHLTTPCCATrackParam t0 = track.InnerParam();
 453     AliHLTTPCCATrackParam t1 = track.OuterParam();
 454
 455     const float maxSin = CAMath::Sin( 60. / 180.*CAMath::Pi() );
 456
 457     bool ok0 = t0.Rotate( dAlpha, maxSin );
 458     bool ok1 = t1.Rotate( dAlpha, maxSin );
 459
 460     bool do0 = ok0;
 461     bool do1 = ok1 && ( !ok0 || t1.SignCosPhi() * t0.SignCosPhi() < 0 );
 462
 463     if ( ok0 && !do1 && ok1 && ( t1.X() < t0.X() ) ) {
 464       do0 = 0;
 465       do1 = 1;
 466     }
 467
 468     if ( do0 ) {
 469       AliHLTTPCCABorderTrack &b = B[nB];
 470       b.SetX( t0.GetX() );
 471       if ( t0.TransportToX( x0, fSliceParam.GetBz( t0 ), maxSin ) ) {
 472         b.SetOK( 1 );
 473         b.SetTrackID( itr );
 474         b.SetNClusters( track.NClusters() );
 475         b.SetIRow( fClusterInfos[ track.FirstClusterRef() + 0 ].IRow() );
 476         b.SetParam( t0 );
 477         nB++;
 478       }
 479     }
 480     if ( do1 ) {
 481       AliHLTTPCCABorderTrack &b = B[nB];
 482       b.SetX( t1.GetX() );
 483       if ( t1.TransportToX( x0, fSliceParam.GetBz( t1 ), maxSin ) ) {
 484         b.SetOK( 1 );
 485         b.SetTrackID( itr );
 486         b.SetNClusters( track.NClusters() );
 487         b.SetIRow( fClusterInfos[ track.FirstClusterRef() + track.NClusters()-1 ].IRow() );
 488         b.SetParam( t1 );
 489         nB++;
 490       }
 491     }
 492     if ( do0 || do1 ) statOK++;
 493     statAll++;
 494   }
 495 }
 496
 497
 498
 499 void AliHLTTPCCAMerger::SplitBorderTracks( int iSlice1, AliHLTTPCCABorderTrack B1[], int N1,
 500     int iSlice2, AliHLTTPCCABorderTrack B2[], int N2
 501                                          )
 502 {
 503   //* split two sets of tracks
 504
 505   float factor2ys = 1.;//1.5;//SG!!!
 506   float factor2zt = 1.;//1.5;//SG!!!
 507   float factor2k = 2.0;//2.2;
 508
 509   factor2k  = 3.5 * 3.5 * factor2k * factor2k;
 510   factor2ys = 3.5 * 3.5 * factor2ys * factor2ys;
 511   factor2zt = 3.5 * 3.5 * factor2zt * factor2zt;
 512
 513   int minNPartHits = 10;//SG!!!
 514   int minNTotalHits = 20;
 515
 516   //float maxDX = fSliceParam.RowX(40) -  fSliceParam.RowX(0);
 517
 518   for ( int i1 = 0; i1 < N1; i1++ ) {
 519     AliHLTTPCCABorderTrack &b1 = B1[i1];
 520     if ( !b1.OK() ) continue;
 521     if ( b1.NClusters() < minNPartHits ) continue;
 522     const AliHLTTPCCATrackParam &t1 = b1.Param();
 523     int iBest2 = -1;
 524     int lBest2 = 0;
 525     int start2 = ( iSlice1 != iSlice2 ) ? 0 : i1 + 1;
 526     for ( int i2 = start2; i2 < N2; i2++ ) {
 527       AliHLTTPCCABorderTrack &b2 = B2[i2];
 528       if ( !b2.OK() ) continue;
 529       if ( b2.NClusters() < minNPartHits ) continue;
 530       if ( b2.NClusters() < lBest2 ) continue;
 531       if ( b1.NClusters() + b2.NClusters() < minNTotalHits ) continue;
 532
 533       //if( TMath::Abs(b1.fX - b2.fX)>maxDX ) continue;
 534
 535       const AliHLTTPCCATrackParam &t2 = b2.Param();
 536
 537       float c = t2.SignCosPhi() * t1.SignCosPhi() >= 0 ? 1 : -1;
 538       float dk = t2.QPt() - c * t1.QPt();
 539       float s2k = t2.Err2QPt() + t1.Err2QPt();
 540
 541       if ( dk*dk > factor2k*s2k ) continue;
 542
 543
 544       float chi2ys = GetChi2( t1.Y(), c * t1.SinPhi(), t1.Cov()[0], c * t1.Cov()[3], t1.Cov()[5],
 545                               t2.Y(),  t2.SinPhi(), t2.Cov()[0],  t2.Cov()[3], t2.Cov()[5] );
 546
 547       if ( chi2ys > factor2ys ) continue;
 548
 549       float chi2zt = GetChi2( t1.Z(), c * t1.DzDs(), t1.Cov()[2], c * t1.Cov()[7], t1.Cov()[9],
 550                               t2.Z(),  t2.DzDs(), t2.Cov()[2],  t2.Cov()[7], t2.Cov()[9] );
 551
 552       if ( chi2zt > factor2zt ) continue;
 553
 554       lBest2 = b2.NClusters();
 555       iBest2 = b2.TrackID();
 556     }
 557
 558     if ( iBest2 < 0 ) continue;
 559
 560     AliHLTTPCCASliceTrackInfo &newTrack1 = fTrackInfos[fSliceTrackInfoStart[iSlice1] + b1.TrackID() ];
 561     AliHLTTPCCASliceTrackInfo &newTrack2 = fTrackInfos[fSliceTrackInfoStart[iSlice2] + iBest2 ];
 562
 563     int old1 = newTrack2.PrevNeighbour();
 564
 565     if ( old1 >= 0 ) {
 566       AliHLTTPCCASliceTrackInfo &oldTrack1 = fTrackInfos[fSliceTrackInfoStart[iSlice1] + old1];
 567       if ( oldTrack1.NClusters()  < newTrack1.NClusters() ) {
 568         newTrack2.SetPrevNeighbour( -1 );
 569         oldTrack1.SetNextNeighbour( -1 );
 570       } else continue;
 571     }
 572     int old2 = newTrack1.NextNeighbour();
 573     if ( old2 >= 0 ) {
 574       AliHLTTPCCASliceTrackInfo &oldTrack2 = fTrackInfos[fSliceTrackInfoStart[iSlice2] + old2];
 575       if ( oldTrack2.NClusters() < newTrack2.NClusters() ) {
 576         oldTrack2.SetPrevNeighbour( -1 );
 577       } else continue;
 578     }
 579     newTrack1.SetNextNeighbour( iBest2 );
 580     newTrack2.SetPrevNeighbour( b1.TrackID() );
 581   }
 582
 583 }
 584
 585
 586 void AliHLTTPCCAMerger::Merging()
 587 {
 588   //* track merging between slices
 589
 590   fOutput->SetNTracks( 0 );
 591   fOutput->SetNTrackClusters( 0 );
 592   fOutput->SetPointers();
 593
 594
 595   // for each slice set number of the next neighbouring slice
 596
 597   int nextSlice[100], prevSlice[100];
 598
 599   for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
 600     nextSlice[iSlice] = iSlice + 1;
 601     prevSlice[iSlice] = iSlice - 1;
 602   }
 603   int mid = fgkNSlices / 2 - 1 ;
 604   int last = fgkNSlices - 1 ;
 605   if ( mid < 0 ) mid = 0; // to avoid compiler warning
 606   if ( last < 0 ) last = 0; //
 607   nextSlice[ mid ] = 0;
 608   prevSlice[ 0 ] = mid;
 609   nextSlice[ last ] = fgkNSlices / 2;
 610   prevSlice[ fgkNSlices/2 ] = last;
 611
 612   int maxNSliceTracks = 0;
 613   for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
 614     if ( maxNSliceTracks < fSliceNTrackInfos[iSlice] ) maxNSliceTracks = fSliceNTrackInfos[iSlice];
 615   }
 616
 617   if ( 1 ) {// merging track segments withing one slice
 618
 619     AliHLTResizableArray<AliHLTTPCCABorderTrack> bord( maxNSliceTracks*10 );
 620
 621     AliHLTTPCCASliceTrackInfo *tmpT = new AliHLTTPCCASliceTrackInfo[maxNSliceTracks];
 622     AliHLTTPCCAClusterInfo *tmpH = new AliHLTTPCCAClusterInfo[fMaxClusterInfos];
 623
 624     for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
 625
 626       int nBord = 0;
 627       MakeBorderTracks( iSlice, 4, bord.Data(), nBord );
 628       SplitBorderTracks( iSlice, bord.Data(), nBord, iSlice, bord.Data(), nBord );
 629
 630       int nTr = 0, nH = 0;
 631       int sliceFirstClusterRef = 0;
 632       for ( int itr = 0; itr < fSliceNTrackInfos[iSlice]; itr++ ) {
 633         AliHLTTPCCASliceTrackInfo &track = fTrackInfos[ fSliceTrackInfoStart[iSlice] + itr];
 634         if ( itr == 0 ) sliceFirstClusterRef = track.FirstClusterRef();
 635         track.SetPrevNeighbour( -1 );
 636         if ( track.NextNeighbour() == -2 ) {
 637           track.SetNextNeighbour( -1 );
 638           continue;
 639         }
 640         AliHLTTPCCASliceTrackInfo &trackNew = tmpT[nTr];
 641         trackNew = track;
 642         trackNew.SetFirstClusterRef( sliceFirstClusterRef + nH );
 643
 644         for ( int ih = 0; ih < track.NClusters(); ih++ ) tmpH[nH+ih] = fClusterInfos[track.FirstClusterRef()+ih];
 645         nTr++;
 646         nH += track.NClusters();
 647
 648         int jtr =  track.NextNeighbour();
 649
 650         if ( jtr < 0 ) continue;
 651         AliHLTTPCCASliceTrackInfo &neighTrack = fTrackInfos[ fSliceTrackInfoStart[iSlice] + jtr];
 652
 653         track.SetNextNeighbour( -1 );
 654         neighTrack.SetNextNeighbour( -2 );
 655
 656         for ( int ih = 0; ih < neighTrack.NClusters(); ih++ )
 657           tmpH[nH+ih] = fClusterInfos[neighTrack.FirstClusterRef()+ih];
 658
 659         trackNew.SetNClusters( trackNew.NClusters() + neighTrack.NClusters() );
 660         trackNew.SetNextNeighbour( -1 );
 661         nH += neighTrack.NClusters();
 662         if ( neighTrack.InnerParam().X() < track.InnerParam().X() ) trackNew.SetInnerParam( neighTrack.InnerParam() );
 663         if ( neighTrack.OuterParam().X() > track.OuterParam().X() ) trackNew.SetOuterParam( neighTrack.OuterParam() );
 664       }
 665
 666       fSliceNTrackInfos[iSlice] = nTr;
 667       for ( int itr = 0; itr < nTr; itr++ ) fTrackInfos[ fSliceTrackInfoStart[iSlice] + itr] = tmpT[itr];
 668       for ( int ih = 0; ih < nH; ih++ ) fClusterInfos[sliceFirstClusterRef + ih] = tmpH[ih];
 669
 670     }
 671     delete[] tmpT;
 672     delete[] tmpH;
 673   }
 674
 675
 676   //* merging tracks between slices
 677
 678
 679   // arrays for the rotated track parameters
 680
 681   AliHLTTPCCABorderTrack
 682   *bCurr0 = new AliHLTTPCCABorderTrack[maxNSliceTracks*10],
 683   *bNext0 = new AliHLTTPCCABorderTrack[maxNSliceTracks*10],
 684   *bCurr = new AliHLTTPCCABorderTrack[maxNSliceTracks*10],
 685   *bNext = new AliHLTTPCCABorderTrack[maxNSliceTracks*10];
 686
 687   for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
 688
 689     int jSlice = nextSlice[iSlice];
 690
 691     int nCurr0 = 0, nNext0 = 0;
 692     int nCurr = 0, nNext = 0;
 693
 694     MakeBorderTracks( iSlice, 0, bCurr, nCurr );
 695     MakeBorderTracks( jSlice, 1, bNext, nNext );
 696     MakeBorderTracks( iSlice, 2, bCurr0, nCurr0 );
 697     MakeBorderTracks( jSlice, 3, bNext0, nNext0 );
 698
 699     SplitBorderTracks( iSlice, bCurr0, nCurr0, jSlice, bNext0, nNext0 );
 700     SplitBorderTracks( iSlice, bCurr, nCurr, jSlice, bNext, nNext );
 701   }
 702
 703   if ( bCurr0 ) delete[] bCurr0;
 704   if ( bNext0 ) delete[] bNext0;
 705   if ( bCurr  ) delete[] bCurr;
 706   if ( bNext  ) delete[] bNext;
 707
 708
 709   //static int nRejected = 0;
 710
 711   int nOutTracks = 0;
 712   int nOutTrackClusters = 0;
 713
 714   AliHLTTPCCAMergedTrack *outTracks = new AliHLTTPCCAMergedTrack[fMaxTrackInfos];
 715   unsigned int   *outClusterId = new unsigned int [fMaxClusterInfos];
 716   UChar_t  *outClusterPackedAmp = new UChar_t [fMaxClusterInfos];
 717
 718   for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
 719
 720     for ( int itr = 0; itr < fSliceNTrackInfos[iSlice]; itr++ ) {
 721
 722       AliHLTTPCCASliceTrackInfo &track = fTrackInfos[fSliceTrackInfoStart[iSlice] + itr];
 723
 724       if ( track.Used() ) continue;
 725       if ( track.PrevNeighbour() >= 0 ) continue;
 726
 727       AliHLTTPCCATrackParam startPoint = track.InnerParam(), endPoint = track.OuterParam();
 728       float startAlpha = track.InnerAlpha(), endAlpha = track.OuterAlpha();
 729
 730       int hits[2000];
 731       int firstHit = 1000;
 732       int nHits = 0;
 733       int jSlice = iSlice;
 734       int jtr = itr;
 735
 736       {
 737         track.SetUsed( 1 );
 738         for ( int jhit = 0; jhit < track.NClusters(); jhit++ ) {
 739           int id = track.FirstClusterRef() + jhit;
 740           hits[firstHit+jhit] = id;
 741         }
 742         nHits = track.NClusters();
 743         jtr = track.NextNeighbour();
 744         jSlice = nextSlice[iSlice];
 745       }
 746
 747       while ( jtr >= 0 ) {
 748         AliHLTTPCCASliceTrackInfo &segment = fTrackInfos[fSliceTrackInfoStart[jSlice] + jtr];
 749         if ( segment.Used() ) break;
 750         segment.SetUsed( 1 );
 751         bool dir = 0;
 752         int startHit = firstHit + nHits;
 753         float d00 = startPoint.GetDistXZ2( segment.InnerParam() );
 754         float d01 = startPoint.GetDistXZ2( segment.OuterParam() );
 755         float d10 = endPoint.GetDistXZ2( segment.InnerParam() );
 756         float d11 = endPoint.GetDistXZ2( segment.OuterParam() );
 757         if ( d00 <= d01 && d00 <= d10 && d00 <= d11 ) {
 758           startPoint = segment.OuterParam();
 759           startAlpha = segment.OuterAlpha();
 760           dir = 1;
 761           firstHit -= segment.NClusters();
 762           startHit = firstHit;
 763         } else if ( d01 <= d10 && d01 <= d11 ) {
 764           startPoint = segment.InnerParam();
 765           startAlpha = segment.InnerAlpha();
 766           dir = 0;
 767           firstHit -= segment.NClusters();
 768           startHit = firstHit;
 769         } else if ( d10 <= d11 ) {
 770           endPoint = segment.OuterParam();
 771           endAlpha = segment.OuterAlpha();
 772           dir = 0;
 773         } else {
 774           endPoint = segment.InnerParam();
 775           endAlpha = segment.InnerAlpha();
 776           dir = 1;
 777         }
 778
 779         for ( int jhit = 0; jhit < segment.NClusters(); jhit++ ) {
 780           int id = segment.FirstClusterRef() + jhit;
 781           hits[startHit+( dir ?( segment.NClusters()-1-jhit ) :jhit )] = id;
 782         }
 783         nHits += segment.NClusters();
 784         jtr = segment.NextNeighbour();
 785         jSlice = nextSlice[jSlice];
 786       }
 787
 788       if ( endPoint.X() < startPoint.X() ) { // swap
 789         for ( int i = 0; i < nHits; i++ ) hits[i] = hits[firstHit+nHits-1-i];
 790         firstHit = 0;
 791       }
 792
 793       if ( nHits < 30 ) continue;    //SG!!!
 794
 795       // refit
 796
 797       // need best t0!!!SG
 798
 799       endPoint = startPoint;
 800       if ( !FitTrack( endPoint, endAlpha, startPoint, startAlpha, hits + firstHit, nHits, 0 ) ) continue;
 801       if ( !FitTrack( startPoint, startAlpha, endPoint, endAlpha, hits + firstHit, nHits, 1 ) ) continue;
 802
 803       if ( nHits < 30 ) continue;    //SG!!!
 804
 805       AliHLTTPCCATrackParam &p = startPoint;
 806
 807       {
 808         double xTPC = 83.65; //SG!!!
 809         double dAlpha = 0.00609235;
 810         AliHLTTPCCATrackParam::AliHLTTPCCATrackFitParam fitPar;
 811         p.CalculateFitParameters( fitPar );
 812
 813         if ( p.TransportToXWithMaterial( xTPC, fitPar, fSliceParam.GetBz( p ) ) ) {
 814           double y = p.GetY();
 815           double ymax = xTPC * CAMath::Tan( dAlpha / 2. );
 816           if ( y > ymax ) {
 817             if ( p.Rotate( dAlpha ) ) { startAlpha += dAlpha;  p.TransportToXWithMaterial( xTPC, fitPar, fSliceParam.GetBz( p ) ); }
 818           } else if ( y < -ymax ) {
 819             if ( p.Rotate( -dAlpha ) ) {  startAlpha -= dAlpha; p.TransportToXWithMaterial( xTPC, fitPar, fSliceParam.GetBz( p ) );}
 820           }
 821         }
 822       }
 823
 824       {
 825         bool ok = 1;
 826
 827         const float *c = p.Cov();
 828         for ( int i = 0; i < 15; i++ ) ok = ok && finite( c[i] );
 829         for ( int i = 0; i < 5; i++ ) ok = ok && finite( p.Par()[i] );
 830         ok = ok && ( p.GetX() > 50 );
 831
 832         if ( c[0] <= 0 || c[2] <= 0 || c[5] <= 0 || c[9] <= 0 || c[14] <= 0 ) ok = 0;
 833         if ( c[0] > 5. || c[2] > 5. || c[5] > 2. || c[9] > 2 || c[14] > 2 ) ok = 0;
 834         if ( !ok ) continue;
 835       }
 836
 837       AliHLTTPCCAMergedTrack &mergedTrack = outTracks[nOutTracks];
 838       mergedTrack.SetNClusters( nHits );
 839       mergedTrack.SetFirstClusterRef( nOutTrackClusters );
 840       mergedTrack.SetInnerParam( startPoint );
 841       mergedTrack.SetInnerAlpha( startAlpha );
 842       mergedTrack.SetOuterParam( endPoint );
 843       mergedTrack.SetOuterAlpha( endAlpha );
 844
 845       for ( int i = 0; i < nHits; i++ ) {
 846         AliHLTTPCCAClusterInfo &clu = fClusterInfos[hits[firstHit+i]];
 847         outClusterId[nOutTrackClusters+i] = clu.Id();
 848         outClusterPackedAmp[nOutTrackClusters+i] = clu.PackedAmp();
 849       }
 850
 851       nOutTracks++;
 852       nOutTrackClusters += nHits;
 853     }
 854   }
 855
 856   fOutput->SetNTracks( nOutTracks );
 857   fOutput->SetNTrackClusters( nOutTrackClusters );
 858   fOutput->SetPointers();
 859
 860   for ( int itr = 0; itr < nOutTracks; itr++ ) fOutput->SetTrack( itr, outTracks[itr] );
 861
 862   for ( int ic = 0; ic < nOutTrackClusters; ic++ ) {
 863     fOutput->SetClusterId( ic, outClusterId[ic] );
 864     fOutput->SetClusterPackedAmp( ic, outClusterPackedAmp[ic] );
 865   }
 866
 867   delete[] outTracks;
 868   delete[] outClusterId;
 869   delete[] outClusterPackedAmp;
 870 }