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