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