code of fast version of GlobalMerger added (it is not yet called by the HLT component)
authorsgorbuno <sgorbuno@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 24 Nov 2010 15:54:05 +0000 (15:54 +0000)
committersgorbuno <sgorbuno@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 24 Nov 2010 15:54:05 +0000 (15:54 +0000)
22 files changed:
HLT/TPCLib/merger-ca/AliHLTTPCGMBorderTrack.h [new file with mode: 0644]
HLT/TPCLib/merger-ca/AliHLTTPCGMCluster.h [new file with mode: 0644]
HLT/TPCLib/merger-ca/AliHLTTPCGMMergedTrack.h [new file with mode: 0644]
HLT/TPCLib/merger-ca/AliHLTTPCGMMerger.cxx [new file with mode: 0644]
HLT/TPCLib/merger-ca/AliHLTTPCGMMerger.h [new file with mode: 0644]
HLT/TPCLib/merger-ca/AliHLTTPCGMSliceTrack.cxx [new file with mode: 0644]
HLT/TPCLib/merger-ca/AliHLTTPCGMSliceTrack.h [new file with mode: 0644]
HLT/TPCLib/merger-ca/AliHLTTPCGMTrackLinearisation.h [new file with mode: 0644]
HLT/TPCLib/merger-ca/AliHLTTPCGMTrackParam.cxx [new file with mode: 0644]
HLT/TPCLib/merger-ca/AliHLTTPCGMTrackParam.h [new file with mode: 0644]
HLT/TPCLib/tracking-ca/AliHLTTPCCAMerger.cxx
HLT/TPCLib/tracking-ca/AliHLTTPCCAMerger.h
HLT/TPCLib/tracking-ca/AliHLTTPCCAParam.cxx
HLT/TPCLib/tracking-ca/AliHLTTPCCAParam.h
HLT/TPCLib/tracking-ca/AliHLTTPCCAPerformance.cxx
HLT/TPCLib/tracking-ca/AliHLTTPCCASliceOutCluster.h
HLT/TPCLib/tracking-ca/AliHLTTPCCASliceOutTrack.h
HLT/TPCLib/tracking-ca/AliHLTTPCCAStandaloneFramework.cxx
HLT/TPCLib/tracking-ca/AliHLTTPCCAStandaloneFramework.h
HLT/TPCLib/tracking-ca/AliHLTTPCCATrackerOutputConverter.cxx
HLT/TPCLib/tracking-ca/AliTPCtrackerCA.cxx
HLT/libAliHLTTPC.pkg

diff --git a/HLT/TPCLib/merger-ca/AliHLTTPCGMBorderTrack.h b/HLT/TPCLib/merger-ca/AliHLTTPCGMBorderTrack.h
new file mode 100644 (file)
index 0000000..8f2325a
--- /dev/null
@@ -0,0 +1,95 @@
+//-*- Mode: C++ -*-
+// ************************************************************************
+// This file is property of and copyright by the ALICE HLT Project        *
+// ALICE Experiment at CERN, All rights reserved.                         *
+// See cxx source for full Copyright notice                               *
+//                                                                        *
+//*************************************************************************
+
+
+#ifndef ALIHLTTPCGMBORDERTRACK_H
+#define ALIHLTTPCGMBORDERTRACK_H
+
+#include "TMath.h"
+
+/**
+ * @class AliHLTTPCGMBorderTrack
+ *
+ * The class describes TPC slice tracks at sector borders. 
+ * Used in AliHLTTPCGMMerger
+ *
+ */
+class AliHLTTPCGMBorderTrack
+{
+
+ public:
+
+  struct Range{
+    int fId;
+    float fMin, fMax;
+    static bool CompMin(const Range &a, const Range &b)  { return a.fMin<b.fMin; }
+    static bool CompMax(const Range &a, const Range &b)  { return a.fMax<b.fMax; }
+  };
+
+
+  int   TrackID()                    const { return fTrackID;   }
+  int   NClusters()                  const { return fNClusters; }  
+  const float *Par() const { return fP; }
+  const float *Cov() const { return fC; }
+  const float *CovD() const { return fD; }
+
+  void SetTrackID   ( int v )                        { fTrackID   = v; }
+  void SetNClusters ( int v )                        { fNClusters = v; }
+  void SetPar( int i, float x ) { fP[i] = x; }
+  void SetCov( int i, float x ) { fC[i] = x; }
+  void SetCovD( int i, float x ) { fD[i] = x; }
+  static float CheckChi2( float x1, float y1, float cx1, float cxy1, float cy1,
+                         float x2, float y2, float cx2, float cxy2, float cy2, float chi2cut  )
+  {
+    //* Calculate Chi2/ndf deviation
+    float dx = x1 - x2;
+    float dy = y1 - y2;
+    float cx = cx1 + cx2;
+    float cxy = cxy1 + cxy2;
+    float cy = cy1 + cy2;
+    float det = cx*cy - cxy*cxy ;
+    return ( ( cy*dx - (cxy+cxy)*dy )*dx + cx*dy*dy < (det+det)*chi2cut );
+  }
+
+  bool CheckChi2Y( const AliHLTTPCGMBorderTrack &t, float chi2cut ) const {
+    float d = fP[0]-t.fP[0];
+    return ( d*d < chi2cut*(fC[0] + t.fC[0]) );
+  }
+
+  bool CheckChi2Z( const AliHLTTPCGMBorderTrack &t, float chi2cut ) const {
+    float d = fP[1]-t.fP[1];
+    return ( d*d < chi2cut *(fC[1] + t.fC[1]) );
+  }
+
+  bool CheckChi2QPt( const AliHLTTPCGMBorderTrack &t, float chi2cut  ) const {
+    float d = fP[4]-t.fP[4];
+    return ( d*d < chi2cut*(fC[4] + t.fC[4]) );
+  }
+  
+  bool CheckChi2YS( const AliHLTTPCGMBorderTrack &t, float chi2cut ) const {
+    return CheckChi2(   fP[0],   fP[2],   fC[0],   fD[0], fC[2],
+                     t.fP[0], t.fP[2], t.fC[0], t.fD[0], t.fC[2], chi2cut );
+  }
+  bool CheckChi2ZT( const AliHLTTPCGMBorderTrack &t, float chi2cut ) const {
+    return  CheckChi2(   fP[1],   fP[3],   fC[1],   fD[1], fC[3],
+                      t.fP[1], t.fP[3], t.fC[1], t.fD[1], t.fC[3], chi2cut );
+    
+  }
+
+ private:
+
+  int   fTrackID;              // track index
+  int   fNClusters;            // n clusters
+  float fP[5];
+  float fC[5];
+  float fD[2];
+};
+
+#endif
diff --git a/HLT/TPCLib/merger-ca/AliHLTTPCGMCluster.h b/HLT/TPCLib/merger-ca/AliHLTTPCGMCluster.h
new file mode 100644 (file)
index 0000000..1daa4d3
--- /dev/null
@@ -0,0 +1,57 @@
+//-*- Mode: C++ -*-
+// $Id: AliHLTTPCGMCluster.h 39008 2010-02-18 17:33:32Z sgorbuno $
+// ************************************************************************
+// This file is property of and copyright by the ALICE HLT Project        *
+// ALICE Experiment at CERN, All rights reserved.                         *
+// See cxx source for full Copyright notice                               *
+//                                                                        *
+//*************************************************************************
+
+
+#ifndef ALIHLTTPCGMCLUSTER_H
+#define ALIHLTTPCGMCLUSTER_H
+#include "TMath.h"
+/**
+ * @class AliHLTTPCGMCluster
+ *
+ * The class describes TPC clusters used in AliHLTTPCGMMerger
+ */
+class AliHLTTPCGMCluster
+{
+  
+ public:
+
+  unsigned char  ISlice()  const { return fISlice;    }
+  unsigned char  IRow()    const { return fIRow;    }
+  int  Id()      const { return fId;      }
+  UChar_t PackedAmp() const { return fPackedAmp; }
+  float X()         const { return fX;         }
+  float Y()         const { return fY;         }
+  float Z()         const { return fZ;         }
+  float Err2Y()     const { return fErr2Y;     }
+  float Err2Z()     const { return fErr2Z;     }
+  
+  void SetISlice    ( unsigned char v  ) { fISlice    = v; }
+  void SetIRow    ( unsigned char v  ) { fIRow    = v; }
+  void SetId      (  int v  ) { fId      = v; }
+  void SetPackedAmp ( UChar_t v ) { fPackedAmp = v; }
+  void SetX         ( float v ) { fX         = v; }
+  void SetY         ( float v ) { fY         = v; }
+  void SetZ         ( float v ) { fZ         = v; }
+  void SetErr2Y     ( float v ) { fErr2Y     = v; }
+  void SetErr2Z     ( float v ) { fErr2Z     = v; }
+
+  //private:
+
+  unsigned char fISlice;   // slice number
+  unsigned char fIRow;     // row number
+  int fId;                 // cluster Id
+  UChar_t fPackedAmp;      // packed amplitude
+  float fX;                // x position (slice coord.system)
+  float fY;                // y position (slice coord.system)
+  float fZ;                // z position (slice coord.system)
+  float fErr2Y;            // Squared measurement error of y position
+  float fErr2Z;            // Squared measurement error of z position
+};
+
+#endif
diff --git a/HLT/TPCLib/merger-ca/AliHLTTPCGMMergedTrack.h b/HLT/TPCLib/merger-ca/AliHLTTPCGMMergedTrack.h
new file mode 100644 (file)
index 0000000..229f8e1
--- /dev/null
@@ -0,0 +1,57 @@
+//-*- Mode: C++ -*-
+// ************************************************************************
+// This file is property of and copyright by the ALICE HLT Project        *
+// ALICE Experiment at CERN, All rights reserved.                         *
+// See cxx source for full Copyright notice                               *
+//                                                                        *
+//*************************************************************************
+
+
+#ifndef ALIHLTTPCGMMERGEDTRACK_H
+#define ALIHLTTPCGMMERGEDTRACK_H
+
+#include "AliHLTTPCGMTrackParam.h"
+
+/**
+ * @class AliHLTTPCGMMergedTrack
+ * 
+ * The class is used to store merged tracks in AliHLTTPCGMMerger
+ */
+class AliHLTTPCGMMergedTrack
+{
+ public:
+
+  int NClusters()                      const { return fNClusters;       }
+  int FirstClusterRef()                const { return fFirstClusterRef; }
+  const AliHLTTPCGMTrackParam &GetParam() const { return fParam;      }
+  float GetAlpha()                        const { return fAlpha;      }
+  AliHLTTPCGMTrackParam &Param() { return fParam;      }
+  float &Alpha()                 { return fAlpha;      }
+  float LastX()                        const { return fLastX; }
+  float LastY()                        const { return fLastY; }
+  float LastZ()                        const { return fLastZ; }
+  bool OK() const{ return fOK; }
+
+  void SetNClusters      ( int v )                { fNClusters = v;       }
+  void SetFirstClusterRef( int v )                { fFirstClusterRef = v; }
+  void SetParam( const AliHLTTPCGMTrackParam &v ) { fParam = v;      }     
+  void SetAlpha( float v )                        { fAlpha = v;      }  
+  void SetLastX( float v )                        { fLastX = v; }
+  void SetLastY( float v )                        { fLastY = v; }
+  void SetLastZ( float v )                        { fLastZ = v; }
+  void SetOK( bool v ) {fOK = v;}
+ private:
+
+  AliHLTTPCGMTrackParam fParam; //* fitted track parameters 
+
+  float fAlpha;                 //* alpha angle 
+  float fLastX; //* outer X
+  float fLastY; //* outer Y
+  float fLastZ; //* outer Z
+  int fFirstClusterRef;         //* index of the first track cluster in corresponding cluster arrays
+  int fNClusters;               //* number of track clusters
+  bool fOK;//
+};
+
+
+#endif 
diff --git a/HLT/TPCLib/merger-ca/AliHLTTPCGMMerger.cxx b/HLT/TPCLib/merger-ca/AliHLTTPCGMMerger.cxx
new file mode 100644 (file)
index 0000000..72934bb
--- /dev/null
@@ -0,0 +1,673 @@
+// $Id: AliHLTTPCGMMerger.cxx 30732 2009-01-22 23:02:02Z sgorbuno $
+// **************************************************************************
+// This file is property of and copyright by the ALICE HLT Project          *
+// ALICE Experiment at CERN, All rights reserved.                           *
+//                                                                          *
+// Primary Authors: Sergey Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de> *
+//                  Ivan Kisel <kisel@kip.uni-heidelberg.de>                *
+//                  for The ALICE HLT Project.                              *
+//                                                                          *
+// Permission to use, copy, modify and distribute this software and its     *
+// documentation strictly for non-commercial purposes is hereby granted     *
+// without fee, provided that the above copyright notice appears in all     *
+// copies and that both the copyright notice and this permission notice     *
+// appear in the supporting documentation. The authors make no claims       *
+// about the suitability of this software for any purpose. It is            *
+// provided "as is" without express or implied warranty.                    *
+//                                                                          *
+//***************************************************************************
+
+#include <stdio.h>
+#include "AliHLTTPCCASliceOutTrack.h"
+#include "AliHLTTPCCATracker.h"
+#include "AliHLTTPCCATrackParam.h"
+#include "AliHLTTPCGMCluster.h"
+
+#include "AliHLTTPCGMMerger.h"
+
+#include "AliHLTTPCCAMath.h"
+#include "TStopwatch.h"
+
+#include "AliHLTTPCCATrackParam.h"
+#include "AliHLTTPCCASliceOutput.h"
+#include "AliHLTTPCGMMergedTrack.h"
+#include "AliHLTTPCCADataCompressor.h"
+#include "AliHLTTPCCAParam.h"
+#include "AliHLTTPCCATrackLinearisation.h"
+#include "AliHLTTPCCADataCompressor.h"
+
+#include "AliHLTTPCGMTrackParam.h"
+#include "AliHLTTPCGMTrackLinearisation.h"
+#include "AliHLTTPCGMSliceTrack.h"
+#include "AliHLTTPCGMBorderTrack.h"
+
+
+
+AliHLTTPCGMMerger::AliHLTTPCGMMerger()
+  :
+  fSliceParam(),
+  fNOutputTracks( 0 ),
+  fOutputTracks( 0 ),
+  fOutputClusterIds(0),
+  fSliceTrackInfos( 0 ),  
+  fMaxSliceTracks(0),
+  fClusterX(0),
+  fClusterY(0),
+  fClusterZ(0),
+  fClusterRowType(0),
+  fClusterAngle(0),
+  fBorderMemory(0),
+  fBorderRangeMemory(0)
+{
+  //* constructor
+  
+  for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
+    fNextSliceInd[iSlice] = iSlice + 1;
+    fPrevSliceInd[iSlice] = iSlice - 1;
+  }
+  int mid = fgkNSlices / 2 - 1 ;
+  int last = fgkNSlices - 1 ;
+  if ( mid < 0 ) mid = 0; // to avoid compiler warning
+  if ( last < 0 ) last = 0; //
+  fNextSliceInd[ mid ] = 0;
+  fPrevSliceInd[ 0 ] = mid;  fNextSliceInd[ last ] = fgkNSlices / 2;
+  fPrevSliceInd[ fgkNSlices/2 ] = last;
+  
+  Clear();
+}
+
+
+AliHLTTPCGMMerger::AliHLTTPCGMMerger(const AliHLTTPCGMMerger&)
+  :
+  fSliceParam(),
+  fNOutputTracks( 0 ),
+  fOutputTracks( 0 ),
+  fOutputClusterIds(0),
+  fSliceTrackInfos( 0 ),  
+  fMaxSliceTracks(0),
+  fClusterX(0),
+  fClusterY(0),
+  fClusterZ(0),
+  fClusterRowType(0),
+  fClusterAngle(0),
+  fBorderMemory(0),
+  fBorderRangeMemory(0)
+{
+  //* dummy
+}
+
+const AliHLTTPCGMMerger &AliHLTTPCGMMerger::operator=(const AliHLTTPCGMMerger&) const
+{
+  //* dummy
+  return *this;
+}
+
+
+AliHLTTPCGMMerger::~AliHLTTPCGMMerger()
+{
+  //* destructor
+  ClearMemory();
+}
+
+void AliHLTTPCGMMerger::Clear()
+{
+  for ( int i = 0; i < fgkNSlices; ++i ) {
+    fkSlices[i] = 0;
+    fSliceNTrackInfos[ i ] = 0;
+    fSliceTrackInfoStart[ i ] = 0;
+  }
+  ClearMemory();
+}
+
+void AliHLTTPCGMMerger::ClearMemory()
+{
+  delete[] fOutputTracks;
+  delete[] fOutputClusterIds;
+  delete[] fSliceTrackInfos;  
+  delete[] fClusterX;
+  delete[] fClusterY;
+  delete[] fClusterZ;
+  delete[] fClusterRowType;
+  delete[] fClusterAngle;
+  delete[] fBorderMemory;
+  delete[] fBorderRangeMemory;
+
+  fNOutputTracks = 0;
+  fOutputTracks = 0;
+  fOutputClusterIds = 0;
+  fSliceTrackInfos = 0;
+  fMaxSliceTracks = 0;
+  fClusterX = 0;
+  fClusterY = 0;
+  fClusterZ = 0;
+  fClusterRowType = 0;
+  fClusterAngle = 0;
+  fBorderMemory = 0;  
+  fBorderRangeMemory = 0;
+}
+
+
+void AliHLTTPCGMMerger::SetSliceData( int index, const AliHLTTPCCASliceOutput *sliceData )
+{
+  fkSlices[index] = sliceData;
+}
+
+
+bool AliHLTTPCGMMerger::Reconstruct()
+{
+  //* main merging routine
+
+  {
+    const double kCLight = 0.000299792458;
+    double constBz = fSliceParam.BzkG() * kCLight;
+
+    AliHLTTPCGMTrackParam::fPolinomialFieldBz[0] = constBz * (  0.999286   );
+    AliHLTTPCGMTrackParam::fPolinomialFieldBz[1] = constBz * ( -4.54386e-7 );
+    AliHLTTPCGMTrackParam::fPolinomialFieldBz[2] = constBz * (  2.32950e-5 );
+    AliHLTTPCGMTrackParam::fPolinomialFieldBz[3] = constBz * ( -2.99912e-7 );
+    AliHLTTPCGMTrackParam::fPolinomialFieldBz[4] = constBz * ( -2.03442e-8 );
+    AliHLTTPCGMTrackParam::fPolinomialFieldBz[5] = constBz * (  9.71402e-8 );    
+  }
+  
+  int nIter = 100;
+  TStopwatch timer;
+  cout<<"Merger..."<<endl;
+  for( int iter=0; iter<nIter; iter++ ){
+    if( !AllocateMemory() ) return 0;
+    UnpackSlices();
+    MergeWithingSlices();
+    MergeSlices();
+    CollectMergedTracks();
+    Refit();
+  }  
+  timer.Stop();  
+  cout<<"\nMerger time = "<<timer.CpuTime()*1.e3/nIter<<" ms\n"<<endl;
+
+  return 1;
+}
+
+
+
+bool AliHLTTPCGMMerger::AllocateMemory()
+{
+  //* memory allocation
+  
+  ClearMemory();
+
+  int nTracks = 0;
+  int nClusters = 0;
+  fMaxSliceTracks  = 0;
+  
+  for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
+    if ( !fkSlices[iSlice] ) continue;
+    nTracks += fkSlices[iSlice]->NTracks();
+    nClusters += fkSlices[iSlice]->NTrackClusters();
+    if( fMaxSliceTracks < fkSlices[iSlice]->NTracks() ) fMaxSliceTracks = fkSlices[iSlice]->NTracks();
+  }
+
+  //cout<<"\nMerger: input "<<nTracks<<" tracks, "<<nClusters<<" clusters"<<endl;
+  
+  fOutputTracks = new AliHLTTPCGMMergedTrack[nTracks];
+  fOutputClusterIds = new UInt_t[nClusters];
+  fSliceTrackInfos = new AliHLTTPCGMSliceTrack[nTracks];  
+  fClusterX = new float[nClusters];
+  fClusterY = new float[nClusters];
+  fClusterZ = new float[nClusters];
+  fClusterRowType = new UInt_t[nClusters];
+  fClusterAngle = new float[nClusters];        
+  fBorderMemory = new AliHLTTPCGMBorderTrack[fMaxSliceTracks*2];
+  fBorderRangeMemory = new AliHLTTPCGMBorderTrack::Range[fMaxSliceTracks*2];  
+
+  return ( ( fOutputTracks!=NULL )
+          && ( fOutputClusterIds!=NULL )
+          && ( fSliceTrackInfos!=NULL )
+          && ( fClusterX!=NULL )
+          && ( fClusterY!=NULL )
+          && ( fClusterZ!=NULL )
+          && ( fClusterRowType!=NULL )
+          && ( fClusterAngle!=NULL )
+          && ( fBorderMemory!=NULL )
+          && ( fBorderRangeMemory!=NULL )
+          );
+}
+
+
+
+void AliHLTTPCGMMerger::UnpackSlices()
+{
+  //* unpack the cluster information from the slice tracks and initialize track info array
+  
+  int nTracksCurrent = 0;
+  for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
+
+    fSliceTrackInfoStart[ iSlice ] = nTracksCurrent;
+
+    fSliceNTrackInfos[ iSlice ] = 0;
+
+    if ( !fkSlices[iSlice] ) continue;
+
+    float alpha = fSliceParam.Alpha( iSlice );
+
+    const AliHLTTPCCASliceOutput &slice = *( fkSlices[iSlice] );
+    const AliHLTTPCCASliceOutTrack *sliceTr = slice.GetFirstTrack();    
+
+    for ( int itr = 0; itr < slice.NTracks(); itr++, sliceTr = sliceTr->GetNextTrack() ) {
+      AliHLTTPCGMSliceTrack &track = fSliceTrackInfos[nTracksCurrent];
+      track.Set( sliceTr, alpha );
+      if( !track.FilterErrors( fSliceParam, .999 ) ) continue;
+      track.SetPrevNeighbour( -1 );
+      track.SetNextNeighbour( -1 );
+      track.SetSliceNeighbour( -1 );
+      track.SetUsed( 0 );           
+      nTracksCurrent++;
+      fSliceNTrackInfos[ iSlice ]++;
+    }
+    
+    //std::cout<<"Unpack slice "<<iSlice<<": ntracks "<<slice.NTracks()<<"/"<<fSliceNTrackInfos[iSlice]<<std::endl;
+  } 
+}
+
+
+
+
+
+void AliHLTTPCGMMerger::MakeBorderTracks( int iSlice, int iBorder, AliHLTTPCGMBorderTrack B[], int &nB )
+{
+  //* prepare slice tracks for merging with next/previous/same sector
+  //* each track transported to the border line
+
+  float fieldBz = fSliceParam.ConstBz();
+  
+  nB = 0;
+  
+  float dAlpha = fSliceParam.DAlpha() / 2;
+  float x0 = 0;
+
+  if ( iBorder == 0 ) { // transport to the left age of the sector and rotate horisontally
+    dAlpha = dAlpha - CAMath::Pi() / 2 ;
+  } else if ( iBorder == 1 ) { //  transport to the right age of the sector and rotate horisontally
+    dAlpha = -dAlpha - CAMath::Pi() / 2 ;
+  } else if ( iBorder == 2 ) { // transport to the left age of the sector and rotate vertically
+    dAlpha = dAlpha;
+    x0 = fSliceParam.RowX( 63 );
+  } else if ( iBorder == 3 ) { // transport to the right age of the sector and rotate vertically
+    dAlpha = -dAlpha;
+    x0 =  fSliceParam.RowX( 63 );
+  } else if ( iBorder == 4 ) { // transport to the middle of the sector, w/o rotation
+    dAlpha = 0;
+    x0 = fSliceParam.RowX( 63 );
+  }
+
+  const float maxSin = CAMath::Sin( 60. / 180.*CAMath::Pi() );
+  float cosAlpha = TMath::Cos( dAlpha );
+  float sinAlpha = TMath::Sin( dAlpha );
+
+  for ( int itr = 0; itr < fSliceNTrackInfos[iSlice]; itr++ ) {
+
+    AliHLTTPCGMSliceTrack &track = fSliceTrackInfos[ fSliceTrackInfoStart[iSlice] + itr ];
+
+    if( track.Used() ) continue;
+    AliHLTTPCGMBorderTrack &b = B[nB];
+    
+    if(  track.TransportToXAlpha( x0, sinAlpha, cosAlpha, fieldBz, b, maxSin)){
+      b.SetTrackID( itr );
+      b.SetNClusters( track.NClusters() );
+      nB++; 
+    }
+  }
+}
+
+
+void AliHLTTPCGMMerger::MergeBorderTracks ( int iSlice1, AliHLTTPCGMBorderTrack B1[],  int N1,
+                                           int iSlice2, AliHLTTPCGMBorderTrack B2[],  int N2 )
+{
+  //* merge two sets of tracks
+
+  //std::cout<<" Merge slices "<<iSlice1<<"+"<<iSlice2<<": tracks "<<N1<<"+"<<N2<<std::endl;
+  int statAll=0, statMerged=0;
+  float factor2ys = 1.5;//1.5;//SG!!!
+  float factor2zt = 1.5;//1.5;//SG!!!
+  float factor2k = 2.0;//2.2;
+
+  factor2k  = 3.5 * 3.5 * factor2k * factor2k;
+  factor2ys = 3.5 * 3.5 * factor2ys * factor2ys;
+  factor2zt = 3.5 * 3.5 * factor2zt * factor2zt;
+  int minNPartHits = 10;//SG!!!
+  int minNTotalHits = 20;
+
+  AliHLTTPCGMBorderTrack::Range *range1 = fBorderRangeMemory;
+  AliHLTTPCGMBorderTrack::Range *range2 = fBorderRangeMemory + N1;
+
+  bool sameSlice = (iSlice1 == iSlice2);
+  {
+    for ( int itr = 0; itr < N1; itr++ ){
+      AliHLTTPCGMBorderTrack &b = B1[itr];
+      //   if( iSlice1==7 && iSlice2==8 ){
+      //cout<<b.TrackID()<<": "<<b.Cov()[0]<<" "<<b.Cov()[1]<<endl;
+      //}
+      float d = 3.5*sqrt(b.Cov()[1]);
+      range1[itr].fId = itr;
+      range1[itr].fMin = b.Par()[1] - d;
+      range1[itr].fMax = b.Par()[1] + d;
+    }
+    std::sort(range1,range1+N1,AliHLTTPCGMBorderTrack::Range::CompMin);
+    if( sameSlice ){
+      for(int i=0; i<N1; i++) range2[i]= range1[i];
+      std::sort(range2,range2+N1,AliHLTTPCGMBorderTrack::Range::CompMax);
+      N2 = N1;
+      B2 = B1;
+    }else{
+      for ( int itr = 0; itr < N2; itr++ ){
+       AliHLTTPCGMBorderTrack &b = B2[itr];
+       float d = 3.5*sqrt(b.Cov()[1]);
+       range2[itr].fId = itr;
+       range2[itr].fMin = b.Par()[1] - d;
+       range2[itr].fMax = b.Par()[1] + d;
+      }        
+      std::sort(range2,range2+N2,AliHLTTPCGMBorderTrack::Range::CompMax);
+    }
+  }
+
+  int i2 = 0;
+  for ( int i1 = 0; i1 < N1; i1++ ) {
+
+    AliHLTTPCGMBorderTrack::Range r1 = range1[i1];
+    while( i2<N2 && range2[i2].fMax< r1.fMin ) i2++;
+    AliHLTTPCGMBorderTrack &b1 = B1[r1.fId];   
+    if ( b1.NClusters() < minNPartHits ) continue;
+    int iBest2 = -1;
+    int lBest2 = 0;
+    statAll++;
+    for( int k2 = i2; k2<N2; k2++){
+      
+      AliHLTTPCGMBorderTrack::Range r2 = range2[k2];
+      if( r2.fMin > r1.fMax ) break;
+      if( sameSlice && (r1.fId >= r2.fId) ) continue;
+      // do check
+      AliHLTTPCGMBorderTrack &b2 = B2[r2.fId];
+      if ( b2.NClusters() < lBest2 ) continue;
+      
+      if( !b1.CheckChi2Y(b2, factor2ys ) ) continue;
+      //if( !b1.CheckChi2Z(b2, factor2zt ) ) continue;
+      if( !b1.CheckChi2QPt(b2, factor2k ) ) continue;
+      if( !b1.CheckChi2YS(b2, factor2ys ) ) continue;
+      if( !b1.CheckChi2ZT(b2, factor2zt ) ) continue;
+      if ( b2.NClusters() < minNPartHits ) continue;
+      if ( b1.NClusters() + b2.NClusters() < minNTotalHits ) continue;      
+
+      lBest2 = b2.NClusters();
+      iBest2 = b2.TrackID();
+    }
+
+    if ( iBest2 < 0 ) continue;
+    statMerged++;
+    AliHLTTPCGMSliceTrack &newTrack1 = fSliceTrackInfos[fSliceTrackInfoStart[iSlice1] + b1.TrackID() ];
+    AliHLTTPCGMSliceTrack &newTrack2 = fSliceTrackInfos[fSliceTrackInfoStart[iSlice2] + iBest2 ];
+
+    int old1 = newTrack2.PrevNeighbour();
+
+    if ( old1 >= 0 ) {
+      AliHLTTPCGMSliceTrack &oldTrack1 = fSliceTrackInfos[fSliceTrackInfoStart[iSlice1] + old1];
+      if ( oldTrack1.NClusters()  < newTrack1.NClusters() ) {
+        newTrack2.SetPrevNeighbour( -1 );
+        oldTrack1.SetNextNeighbour( -1 );
+      } else continue;
+    }
+    int old2 = newTrack1.NextNeighbour();
+    if ( old2 >= 0 ) {
+      AliHLTTPCGMSliceTrack &oldTrack2 = fSliceTrackInfos[fSliceTrackInfoStart[iSlice2] + old2];
+      if ( oldTrack2.NClusters() < newTrack2.NClusters() ) {
+        oldTrack2.SetPrevNeighbour( -1 );
+      } else continue;
+    }
+    newTrack1.SetNextNeighbour( iBest2 );
+    newTrack2.SetPrevNeighbour( b1.TrackID() );
+  }
+  //cout<<"slices "<<iSlice1<<","<<iSlice2<<": all "<<statAll<<" merged "<<statMerged<<endl;
+}
+
+
+void AliHLTTPCGMMerger::MergeWithingSlices()
+{
+  //* merge track segments withing one slice
+
+  float x0 = fSliceParam.RowX( 63 );  
+  const float maxSin = CAMath::Sin( 60. / 180.*CAMath::Pi() );
+
+  for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
+
+    int nBord = 0;
+    for ( int itr = 0; itr < fSliceNTrackInfos[iSlice]; itr++ ) {
+      AliHLTTPCGMSliceTrack &track = fSliceTrackInfos[ fSliceTrackInfoStart[iSlice] + itr ];
+      //track.SetPrevNeighbour( -1 );      
+      //track.SetNextNeighbour( -1 );
+      //track.SetSliceNeighbour( -1 );
+      //track.SetUsed(0);
+      
+      AliHLTTPCGMBorderTrack &b = fBorderMemory[nBord];
+      if( track.TransportToX( x0, fSliceParam.ConstBz(), b, maxSin) ){
+       b.SetTrackID( itr );
+       b.SetNClusters( track.NClusters() );
+       nBord++;
+      }
+    }  
+
+    MergeBorderTracks( iSlice, fBorderMemory, nBord, iSlice, fBorderMemory, nBord );
+    
+    for ( int itr = 0; itr < fSliceNTrackInfos[iSlice]; itr++ ) {
+      AliHLTTPCGMSliceTrack &track = fSliceTrackInfos[ fSliceTrackInfoStart[iSlice] + itr];
+      if( track.PrevNeighbour()>=0 || track.Used() ) continue;
+      int jtr = track.NextNeighbour();
+      track.SetSliceNeighbour( jtr );
+      track.SetNextNeighbour(-1);
+      while( jtr>=0 ){
+       AliHLTTPCGMSliceTrack &trackN = fSliceTrackInfos[ fSliceTrackInfoStart[iSlice] + jtr];
+       if( trackN.NClusters()>track.NClusters() ) track.CopyParamFrom(trackN);
+       trackN.SetUsed(2);
+       jtr = trackN.NextNeighbour();
+       trackN.SetSliceNeighbour( jtr );
+       trackN.SetNextNeighbour(-1);
+       trackN.SetPrevNeighbour(-1);
+      }
+    }
+  }
+}
+
+
+
+
+void AliHLTTPCGMMerger::MergeSlices()
+{
+  //* track merging between slices
+
+  //for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
+  //for ( int itr = 0; itr < fSliceNTrackInfos[iSlice]; itr++ ) {
+  //AliHLTTPCGMSliceTrack &track = fSliceTrackInfos[ fSliceTrackInfoStart[iSlice] + itr ];
+  //track.SetPrevNeighbour( -1 );
+  //track.SetNextNeighbour( -1 );
+  //}
+  //}
+
+  AliHLTTPCGMBorderTrack 
+    *bCurr = fBorderMemory,
+    *bNext = fBorderMemory + fMaxSliceTracks;
+  
+  for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {    
+    int jSlice = fNextSliceInd[iSlice];    
+    int nCurr = 0, nNext = 0;
+    MakeBorderTracks( iSlice, 2, bCurr, nCurr );
+    MakeBorderTracks( jSlice, 3, bNext, nNext );
+    MergeBorderTracks( iSlice, bCurr, nCurr, jSlice, bNext, nNext );
+    MakeBorderTracks( iSlice, 0, bCurr, nCurr );
+    MakeBorderTracks( jSlice, 1, bNext, nNext );
+    MergeBorderTracks( iSlice, bCurr, nCurr, jSlice, bNext, nNext );
+  }
+}
+
+
+
+
+
+void AliHLTTPCGMMerger::CollectMergedTracks()
+{
+  //* 
+
+  //for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
+  //for ( int itr = 0; itr < fSliceNTrackInfos[iSlice]; itr++ ) {
+  //AliHLTTPCGMSliceTrack &track = fSliceTrackInfos[ fSliceTrackInfoStart[iSlice] + itr ];
+  //if( track.Used()!=2 ) track.SetUsed(0);
+  //}
+  //}
+
+  fNOutputTracks = 0;
+  int nOutTrackClusters = 0;
+
+  const AliHLTTPCGMSliceTrack *trackParts[100];
+
+  for ( int iSlice = 0; iSlice < fgkNSlices; iSlice++ ) {
+
+    for ( int itr = 0; itr < fSliceNTrackInfos[iSlice]; itr++ ) {
+
+      AliHLTTPCGMSliceTrack &track = fSliceTrackInfos[fSliceTrackInfoStart[iSlice] + itr];
+
+      if ( track.Used() ) continue;
+      if ( track.PrevNeighbour() >= 0 ) continue;
+      int nParts = 0;
+      int jSlice = iSlice;
+      AliHLTTPCGMSliceTrack *trbase = &track, *tr = &track;
+      tr->SetUsed( 1 );
+      do{
+       trackParts[nParts++] = tr;
+       int jtr = tr->SliceNeighbour();
+       if( jtr >= 0 ) {
+         tr = &(fSliceTrackInfos[fSliceTrackInfoStart[jSlice] + jtr]);
+         tr->SetUsed( 2 );
+         continue;
+       }
+        jtr = trbase->NextNeighbour();
+       if( jtr>=0 ){
+         jSlice = fNextSliceInd[jSlice];
+         trbase = &(fSliceTrackInfos[fSliceTrackInfoStart[jSlice] + jtr]);
+         tr = trbase;
+         if( tr->Used() ) break;
+         tr->SetUsed( 1 );
+         continue;       
+       }
+       break;
+      }while(1);
+
+      // unpack and sort clusters
+      
+      sort(trackParts, trackParts+nParts, CompareTrackParts );
+
+      AliHLTTPCCASliceOutCluster tmp[400];
+      int currCluster = nOutTrackClusters;
+      int nHits = 0;
+      for( int ipart=0; ipart<nParts; ipart++ ){
+       const AliHLTTPCGMSliceTrack *t = trackParts[ipart];
+       int nTrackHits = t->NClusters();
+       const AliHLTTPCCASliceOutCluster *c= t->OrigTrack()->Clusters();
+       AliHLTTPCCASliceOutCluster *c2 = tmp+nHits + nTrackHits-1;
+       for( int i=0; i<nTrackHits; i++, c++, c2-- ) *c2 = *c;  
+       float alpha =  t->Alpha();
+       for( int i=0; i<nTrackHits; i++) fClusterAngle[currCluster++] = alpha;
+       nHits+=nTrackHits;
+      }
+
+      UInt_t *clId = fOutputClusterIds + nOutTrackClusters;      
+      for( int i=0; i<nHits; i++ ) clId[i] = tmp[i].GetId();      
+      
+      UInt_t *clT  = fClusterRowType + nOutTrackClusters;
+      for( int i=0; i<nHits; i++ ) clT[i] = tmp[i].GetRowType();  
+      
+      float *clX = fClusterX + nOutTrackClusters;
+      for( int i=0; i<nHits; i++ ) clX[i] = tmp[i].GetX();      
+
+      float *clY = fClusterY + nOutTrackClusters;
+      for( int i=0; i<nHits; i++ ) clY[i] = tmp[i].GetY();      
+
+      float *clZ = fClusterZ + nOutTrackClusters;
+      for( int i=0; i<nHits; i++ ) clZ[i] = tmp[i].GetZ();      
+
+      if ( nHits < 30 ) continue;   
+
+      AliHLTTPCGMMergedTrack &mergedTrack = fOutputTracks[fNOutputTracks];
+      mergedTrack.SetOK(1);
+      mergedTrack.SetNClusters( nHits );
+      mergedTrack.SetFirstClusterRef( nOutTrackClusters );
+      AliHLTTPCGMTrackParam &p1 = mergedTrack.Param();
+      const AliHLTTPCGMSliceTrack &p2 = *(trackParts[0]);
+
+      p1.X() = p2.X();
+      p1.Y() = p2.Y();
+      p1.Z() = p2.Z();
+      p1.SinPhi()  = p2.SinPhi();
+      p1.DzDs()  = p2.DzDs();
+      p1.QPt()  = p2.QPt();
+      mergedTrack.SetAlpha( p2.Alpha() );
+
+      fNOutputTracks++;
+      nOutTrackClusters += nHits;
+    }
+  }
+}
+
+
+
+void AliHLTTPCGMMerger::Refit()
+{
+  //* final refit
+
+  for ( int itr = 0; itr < fNOutputTracks; itr++ ) {
+    AliHLTTPCGMMergedTrack &track = fOutputTracks[itr];
+    if( !track.OK() ) continue;    
+    track.SetOK(0);
+
+    int nTrackHits = track.NClusters();
+       
+    AliHLTTPCGMTrackParam t = track.Param();
+    float Alpha = track.Alpha();  
+    
+    t.Fit( fClusterX+track.FirstClusterRef(),
+          fClusterY+track.FirstClusterRef(),
+          fClusterZ+track.FirstClusterRef(),
+          fClusterRowType+track.FirstClusterRef(),
+          fClusterAngle+track.FirstClusterRef(),
+          fSliceParam, nTrackHits, Alpha, 0 );      
+    
+    if ( nTrackHits < 30 ) continue;    //SG!!!
+
+    if ( fabs( t.QPt() ) < 1.e-4 ) t.QPt() = 1.e-4 ;
+    
+    bool ok = t.CheckNumericalQuality();
+    
+    if ( fabs( t.SinPhi() ) > .999 ) ok = 0;
+    
+    if( !ok ) continue;
+    
+    if( 1 ){//SG!!!
+      track.SetNClusters( nTrackHits );
+      track.Param() = t;
+      track.Alpha() = Alpha;
+    }
+        
+    track.SetOK(1);
+
+    {
+      int ind = track.FirstClusterRef();
+      float alpha = fClusterAngle[ind];
+      float x = fClusterX[ind];
+      float y = fClusterY[ind];
+      float z = fClusterZ[ind];
+      float sinA = TMath::Sin( alpha - track.Alpha());
+      float cosA = TMath::Cos( alpha - track.Alpha());
+      track.SetLastX( x*cosA - y*sinA );
+      track.SetLastY( x*sinA + y*cosA );
+      track.SetLastZ( z );
+    }
+  }
+}
+
diff --git a/HLT/TPCLib/merger-ca/AliHLTTPCGMMerger.h b/HLT/TPCLib/merger-ca/AliHLTTPCGMMerger.h
new file mode 100644 (file)
index 0000000..0128b08
--- /dev/null
@@ -0,0 +1,99 @@
+//-*- Mode: C++ -*-
+// ************************************************************************
+// This file is property of and copyright by the ALICE HLT Project        *
+// ALICE Experiment at CERN, All rights reserved.                         *
+// See cxx source for full Copyright notice                               *
+//                                                                        *
+//*************************************************************************
+
+#ifndef ALIHLTTPCGMMERGER_H
+#define ALIHLTTPCGMMERGER_H
+
+#include "AliHLTTPCCADef.h"
+#include "AliHLTTPCCAParam.h"
+#include "AliHLTTPCGMBorderTrack.h"
+#include "AliHLTTPCGMSliceTrack.h"
+
+#if !defined(HLTCA_GPUCODE)
+#include <iostream>
+#endif //HLTCA_GPUCODE
+
+class AliHLTTPCCASliceTrack;
+class AliHLTTPCCASliceOutput;
+class AliHLTTPCGMCluster;
+class AliHLTTPCGMTrackParam;
+class AliHLTTPCGMMergedTrack;
+
+/**
+ * @class AliHLTTPCGMMerger
+ *
+ */
+class AliHLTTPCGMMerger
+{
+  
+public:
+  
+  AliHLTTPCGMMerger();
+  ~AliHLTTPCGMMerger();
+  
+  void SetSliceParam( const AliHLTTPCCAParam &v ) { fSliceParam = v; }
+  
+  void Clear();
+  void SetSliceData( int index, const AliHLTTPCCASliceOutput *SliceData );
+  bool Reconstruct();
+  
+  Int_t NOutputTracks() const { return fNOutputTracks; }
+  const AliHLTTPCGMMergedTrack * OutputTracks() const { return fOutputTracks; }
+  const UInt_t * OutputClusterIds() const { return fOutputClusterIds; }
+   
+  
+  const AliHLTTPCCAParam &SliceParam() const { return fSliceParam; }
+  
+private:
+  
+  AliHLTTPCGMMerger( const AliHLTTPCGMMerger& );
+
+  const AliHLTTPCGMMerger &operator=( const AliHLTTPCGMMerger& ) const;
+  
+  void MakeBorderTracks( int iSlice, int iBorder, AliHLTTPCGMBorderTrack B[], int &nB );
+
+  void MergeBorderTracks( int iSlice1, AliHLTTPCGMBorderTrack B1[],  int N1,
+                         int iSlice2, AliHLTTPCGMBorderTrack B2[],  int N2 );
+  
+  static bool CompareTrackParts( const AliHLTTPCGMSliceTrack *t1, const AliHLTTPCGMSliceTrack *t2 ){
+    return (t1->X() > t2->X() );
+  }
+
+  void ClearMemory();
+  bool AllocateMemory();
+  void UnpackSlices();
+  void MergeWithingSlices();
+  void MergeSlices();
+  void CollectMergedTracks();
+  void Refit();
+  
+  static const int fgkNSlices = 36;       //* N slices
+  int fNextSliceInd[fgkNSlices];
+  int fPrevSliceInd[fgkNSlices];
+  
+  AliHLTTPCCAParam fSliceParam;           //* slice parameters (geometry, calibr, etc.)
+  const AliHLTTPCCASliceOutput *fkSlices[fgkNSlices]; //* array of input slice tracks
+
+  Int_t fNOutputTracks;
+  AliHLTTPCGMMergedTrack *fOutputTracks;       //* array of output merged tracks
+  UInt_t * fOutputClusterIds;
+  
+  AliHLTTPCGMSliceTrack *fSliceTrackInfos; //* additional information for slice tracks
+  int fSliceTrackInfoStart[fgkNSlices];   //* slice starting index in fTrackInfos array;
+  int fSliceNTrackInfos[fgkNSlices];      //* N of slice track infos in fTrackInfos array;
+  int fMaxSliceTracks;      // max N tracks in one slice
+  float *fClusterX;         // cluster X
+  float *fClusterY;         // cluster Y
+  float *fClusterZ;         // cluster Z
+  UInt_t *fClusterRowType;  // cluster row type
+  float *fClusterAngle;     // angle    
+  AliHLTTPCGMBorderTrack *fBorderMemory; // memory for border tracks
+  AliHLTTPCGMBorderTrack::Range *fBorderRangeMemory; // memory for border tracks
+};
+
+#endif //ALIHLTTPCCAMERGER_H
diff --git a/HLT/TPCLib/merger-ca/AliHLTTPCGMSliceTrack.cxx b/HLT/TPCLib/merger-ca/AliHLTTPCGMSliceTrack.cxx
new file mode 100644 (file)
index 0000000..6f36508
--- /dev/null
@@ -0,0 +1,404 @@
+// $Id: AliHLTTPCGMSliceTrack.cxx 41769 2010-06-16 13:58:00Z sgorbuno $
+// **************************************************************************
+// This file is property of and copyright by the ALICE HLT Project          *
+// ALICE Experiment at CERN, All rights reserved.                           *
+//                                                                          *
+// Primary Authors: Sergey Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de> *
+//                  for The ALICE HLT Project.                              *
+//                                                                          *
+// Permission to use, copy, modify and distribute this software and its     *
+// documentation strictly for non-commercial purposes is hereby granted     *
+// without fee, provided that the above copyright notice appears in all     *
+// copies and that both the copyright notice and this permission notice     *
+// appear in the supporting documentation. The authors make no claims       *
+// about the suitability of this software for any purpose. It is            *
+// provided "as is" without express or implied warranty.                    *
+//                                                                          *
+//***************************************************************************
+
+#include "AliHLTTPCGMSliceTrack.h"
+#include "AliHLTTPCCAMath.h"
+#include "AliHLTTPCGMBorderTrack.h"
+#include "AliHLTTPCGMTrackLinearisation.h"
+#include "Riostream.h"
+#include "AliHLTTPCCAParam.h"
+
+
+
+bool AliHLTTPCGMSliceTrack::FilterErrors( AliHLTTPCCAParam &param, float maxSinPhi )
+{
+  float lastX = fOrigTrack->Cluster(fOrigTrack->NClusters()-1 ).GetX();
+
+  const int N = 3;
+  const float *cy = param.GetParamS0Par(0,0);
+  const float *cz = param.GetParamS0Par(1,0);
+
+  float bz = -param.ConstBz();
+  const float kZLength = 250.f - 0.275f;
+
+  float trDzDs2 = fDzDs*fDzDs;
+  float k  = fQPt*bz;               
+  float dx = .33*(lastX - fX);  
+  float kdx = k*dx;
+  float dxBz = dx * bz;
+  float kdx205 = 2.f+kdx*kdx*0.5f;
+
+  {
+    float secPhi2 = fSecPhi*fSecPhi;
+    float zz = fabs( kZLength - fabs(fZ) );    
+    float zz2 = zz*zz;
+    float angleY2 = secPhi2 - 1.f; 
+    float angleZ2 = trDzDs2 * secPhi2 ;
+    
+    float cy0 = cy[0] + cy[1]*zz + cy[3]*zz2;
+    float cy1 = cy[2] + cy[5]*zz;
+    float cy2 = cy[4];
+    float cz0 = cz[0] + cz[1]*zz + cz[3]*zz2;
+    float cz1 = cz[2] + cz[5]*zz;
+    float cz2 = cz[4];
+    
+    fC0 = fabs( cy0 + angleY2 * ( cy1 + angleY2*cy2 ) );
+    fC2 = fabs( cz0 + angleZ2 * ( cz1 + angleZ2*cz2 ) );    
+
+    fC3 = 0;
+    fC5 = 1;
+    fC7 = 0;
+    fC9 = 1;
+    fC10 = 0;
+    fC12 = 0;
+    fC14 = 10;
+  }
+
+  for( int iStep=0; iStep<N; iStep++ ){    
+
+    float err2Y, err2Z;
+
+    { // transport block
+      float ex = fCosPhi; 
+      float ey = fSinPhi;
+      float ey1 = kdx + ey;
+      if( fabs( ey1 ) > maxSinPhi ) return 0;
+
+      float ss = ey + ey1;      
+      float ex1 = sqrt(1.f - ey1*ey1);
+          
+      float cc = ex + ex1;  
+      float dxcci = dx / cc;
+      
+      float dy = dxcci * ss;
+      float norm2 = 1.f + ey*ey1 + ex*ex1;
+      float dl = dxcci * sqrt( norm2 + norm2 );
+     
+      float dS;   
+      {
+       float dSin = 0.5f*k*dl;
+       float a = dSin*dSin;
+       const float k2 = 1.f/6.f;
+       const float k4 = 3.f/40.f;
+       dS = dl + dl*a*(k2 + a*(k4 ));//+ k6*a) );
+      }
+      float dz = dS * fDzDs;      
+      float ex1i =1.f/ex1;
+      float z = fZ+ dz;
+      {        
+       float secPhi2 = ex1i*ex1i;
+       float zz = fabs( kZLength - fabs(z) );  
+       float zz2 = zz*zz;
+       float angleY2 = secPhi2 - 1.f; 
+       float angleZ2 = trDzDs2 * secPhi2 ;
+
+       float cy0 = cy[0] + cy[1]*zz + cy[3]*zz2;
+       float cy1 = cy[2] + cy[5]*zz;
+       float cy2 = cy[4];
+       float cz0 = cz[0] + cz[1]*zz + cz[3]*zz2;
+       float cz1 = cz[2] + cz[5]*zz;
+       float cz2 = cz[4];
+       
+       err2Y = fabs( cy0 + angleY2 * ( cy1 + angleY2*cy2 ) );
+       err2Z = fabs( cz0 + angleZ2 * ( cz1 + angleZ2*cz2 ) );
+      }
+
+      float hh = kdx205 * dxcci*ex1i; 
+      float h2 = hh * fSecPhi;
+
+      fX+=dx;      
+      fY+= dy;
+      fZ+= dz;
+      fSinPhi = ey1;
+      fCosPhi = ex1;
+      fSecPhi = ex1i;
+    
+      float h4 = bz*dxcci*hh;
+      
+      float c20 = fC3;
+      float c22 = fC5;
+      float c31 = fC7;
+      float c33 = fC9;
+      float c40 = fC10;
+      float c42 = fC12;
+      float c44 = fC14;
+      
+      float c20ph4c42 =  c20 + h4*c42;
+      float h2c22 = h2*c22;
+      float h4c44 = h4*c44;
+      float n7 = c31 + dS*c33;
+      float n10 = c40 + h2*c42 + h4c44;
+      float n12 = c42 + dxBz*c44;
+      
+      
+      fC0+= h2*h2c22 + h4*h4c44 + 2.f*( h2*c20ph4c42  + h4*c40 );
+      
+      fC3 = c20ph4c42 + h2c22  + dxBz*n10;
+      fC10 = n10;
+      
+      fC5 = c22 + dxBz*( c42 + n12 );
+      fC12 = n12;
+      
+      fC2+= dS*(c31 + n7);
+      fC7 = n7; 
+      
+   } // end transport block 
+
+
+    // Filter block
+
+    float 
+      c00 = fC0,
+      c11 = fC2,
+      c20 = fC3,
+      c31 = fC7,
+      c40 = fC10;
+                    
+    float mS0 = 1.f/(err2Y + c00);    
+    float mS2 = 1.f/(err2Z + c11);
+            
+    // K = CHtS
+    
+    float k00, k11, k20, k31, k40;
+    
+    k00 = c00 * mS0;
+    k20 = c20 * mS0;
+    k40 = c40 * mS0;
+    
+    fC0 -= k00 * c00 ;
+    fC5 -= k20 * c20 ;
+    fC10 -= k00 * c40 ;
+    fC12 -= k40 * c20 ;
+    fC3 -= k20 * c00 ;
+    fC14 -= k40 * c40 ;
+        
+    k11 = c11 * mS2;
+    k31 = c31 * mS2;
+    
+    fC7 -= k31 * c11;
+    fC2 -= k11 * c11;
+    fC9 -= k31 * c31;   
+  }
+
+  //* Check that the track parameters and covariance matrix are reasonable
+
+  bool ok = 1;
+  
+  const float *c = &fX;
+  for ( int i = 0; i < 17; i++ ) ok = ok && finite( c[i] );
+
+  if ( fC0 <= 0.f || fC2 <= 0.f || fC5 <= 0.f || fC9 <= 0.f || fC14 <= 0.f 
+       || fC0 > 5.f || fC2 > 5.f || fC5 > 2.f || fC9 > 2.f             ) ok = 0;
+
+  if( ok ){
+    ok = ok 
+      && ( fC3*fC3<=fC5*fC0 )
+      && ( fC7*fC7<=fC9*fC2 )
+      && ( fC10*fC10<=fC14*fC0 )
+      && ( fC12*fC12<=fC14*fC5 );
+  }
+  return ok;
+}
+
+
+
+bool AliHLTTPCGMSliceTrack::TransportToX( float x, float Bz, AliHLTTPCGMBorderTrack &b, float maxSinPhi ) const 
+{
+  Bz = -Bz;
+  float ex = fCosPhi;
+  float ey = fSinPhi;
+  float k  = fQPt*Bz;
+  float dx = x - fX;
+  float ey1 = k*dx + ey;
+  
+  if( fabs( ey1 ) > maxSinPhi ) return 0;
+
+  float ex1 = sqrt( 1.f - ey1 * ey1 );
+  float dxBz = dx * Bz;
+    
+  float ss = ey + ey1;
+  float cc = ex + ex1;  
+  float dxcci = dx / cc;
+  float norm2 = 1.f + ey*ey1 + ex*ex1;
+
+  float dy = dxcci * ss;
+
+  float dS;    
+  {
+    float dl = dxcci * sqrt( norm2 + norm2 );
+    float dSin = 0.5f*k*dl;
+    float a = dSin*dSin;
+    const float k2 = 1.f/6.f;
+    const float k4 = 3.f/40.f;
+    //const float k6 = 5.f/112.f;
+    dS = dl + dl*a*(k2 + a*(k4 ));//+ k6*a) );
+  }
+  
+  float ex1i = 1.f/ex1;
+  float dz = dS * fDzDs;
+
+  float hh = dxcci*ex1i*norm2; 
+  float h2 = hh *fSecPhi;
+  float h4 = Bz*dxcci*hh;
+  
+  float c20 = fC3;
+  float c22 = fC5;  
+  float c31 = fC7;  
+  float c33 = fC9;
+  float c40 = fC10;  
+  float c42 = fC12;
+  float c44 = fC14;
+
+  float c20ph4c42 =  c20 + h4*c42;
+  float h2c22 = h2*c22;
+  float h4c44 = h4*c44;
+  float n7 = c31 + dS*c33;
+  
+  b.SetPar(0, fY + dy );
+  b.SetPar(1, fZ + dz );
+  b.SetPar(2, ey1 );
+  b.SetPar(3, fDzDs);
+  b.SetPar(4, fQPt);
+
+  b.SetCov(0, fC0 + h2*h2c22 + h4*h4c44 + 2.f*( h2*c20ph4c42  + h4*c40 ));
+  b.SetCov(1, fC2+ dS*(c31 + n7) );
+  b.SetCov(2, c22 + dxBz*( c42 + c42 + dxBz*c44 ));
+  b.SetCov(3, c33);
+  b.SetCov(4, c44);
+  b.SetCovD(0, c20ph4c42 + h2c22  + dxBz*(c40 + h2*c42 + h4c44) );
+  b.SetCovD(1, n7 );   
+  return 1;
+}
+
+
+
+bool AliHLTTPCGMSliceTrack::TransportToXAlpha( float newX, float sinAlpha, float cosAlpha, float Bz, AliHLTTPCGMBorderTrack &b, float maxSinPhi ) const 
+{
+  //* 
+
+  float c00 = fC0;
+  float c11 = fC2;
+  float c20 = fC3;
+  float c22 = fC5;  
+  float c31 = fC7;  
+  float c33 = fC9;
+  float c40 = fC10;  
+  float c42 = fC12;
+  float c44 = fC14;
+
+  float x,y;
+  float z = fZ;
+  float sinPhi = fSinPhi;
+  float cosPhi = fCosPhi;
+  float secPhi = fSecPhi;
+  float dzds = fDzDs;
+  float qpt = fQPt;
+
+  // Rotate the coordinate system in XY on the angle alpha
+  {
+    float sP = sinPhi, cP = cosPhi;
+    cosPhi =  cP * cosAlpha + sP * sinAlpha;
+    sinPhi = -cP * sinAlpha + sP * cosAlpha;
+    
+    if ( CAMath::Abs( sinPhi ) > .999 || CAMath::Abs( cP ) < 1.e-2  ) return 0;
+    
+    secPhi = 1./cosPhi;
+    float j0 = cP *secPhi;
+    float j2 = cosPhi / cP;    
+    x =   fX*cosAlpha +  fY*sinAlpha ;
+    y =  -fX*sinAlpha +  fY*cosAlpha ;    
+    
+    c00 *= j0 * j0;
+    c40 *= j0;
+    
+    c22 *= j2 * j2;    
+    c42 *= j2;    
+    if( cosPhi < 0.f ){ // rotate to 180'
+      cosPhi = -cosPhi;
+      secPhi = -secPhi;
+      sinPhi = -sinPhi;
+      dzds = -dzds;
+      qpt = -qpt;      
+      c20 = -c20;
+      c31 = -c31;
+      c40 = -c40;
+   }
+  }
+
+  Bz = -Bz;
+  float ex = cosPhi;
+  float ey = sinPhi;
+  float k  = qpt*Bz;
+  float dx = newX - x;
+  float ey1 = k*dx + ey;
+  
+  if( fabs( ey1 ) > maxSinPhi ) return 0;
+
+  float ex1 = sqrt( 1.f - ey1 * ey1 );
+
+  float dxBz = dx * Bz;
+    
+  float ss = ey + ey1;
+  float cc = ex + ex1;  
+  float dxcci = dx / cc;
+  float norm2 = 1.f + ey*ey1 + ex*ex1;
+
+  float dy = dxcci * ss;
+
+  float dS;    
+  {
+    float dl = dxcci * sqrt( norm2 + norm2 );
+    float dSin = 0.5f*k*dl;
+    float a = dSin*dSin;
+    const float k2 = 1.f/6.f;
+    const float k4 = 3.f/40.f;
+    //const float k6 = 5.f/112.f;
+    dS = dl + dl*a*(k2 + a*(k4 ));//+ k6*a) );
+  }
+  
+  float ex1i = 1.f/ex1;
+  float dz = dS * dzds;
+
+  float hh = dxcci*ex1i*norm2; 
+  float h2 = hh * secPhi;
+  float h4 = Bz*dxcci*hh;  
+
+  float c20ph4c42 =  c20 + h4*c42;
+  float h2c22 = h2*c22;
+  float h4c44 = h4*c44;
+  float n7 = c31 + dS*c33;
+  
+  b.SetPar(0, y + dy );
+  b.SetPar(1, z + dz );
+  b.SetPar(2, ey1 );
+  b.SetPar(3, dzds);
+  b.SetPar(4, qpt);
+  
+  b.SetCov(0, c00 + h2*h2c22 + h4*h4c44 + 2.f*( h2*c20ph4c42  + h4*c40 ));
+  b.SetCov(1, c11 + dS*(c31 + n7) );
+  b.SetCov(2, c22 + dxBz*( c42 + c42 + dxBz*c44 ));
+  b.SetCov(3, c33);
+  b.SetCov(4, c44);
+  b.SetCovD(0, c20ph4c42 + h2c22  + dxBz*(c40 + h2*c42 + h4c44) );
+  b.SetCovD(1, n7 ); 
+
+  return 1;
+}
diff --git a/HLT/TPCLib/merger-ca/AliHLTTPCGMSliceTrack.h b/HLT/TPCLib/merger-ca/AliHLTTPCGMSliceTrack.h
new file mode 100644 (file)
index 0000000..9fab349
--- /dev/null
@@ -0,0 +1,91 @@
+//-*- Mode: C++ -*-
+// ************************************************************************
+// This file is property of and copyright by the ALICE HLT Project        *
+// ALICE Experiment at CERN, All rights reserved.                         *
+// See cxx source for full Copyright notice                               *
+//                                                                        *
+//*************************************************************************
+
+
+#ifndef ALIHLTTPCGMSLICETRACK_H
+#define ALIHLTTPCGMSLICETRACK_H
+
+#include "TMath.h"
+#include "AliHLTTPCGMTrackParam.h"
+#include "AliHLTTPCCASliceOutTrack.h"
+
+/**
+ * @class AliHLTTPCGMSliceTrack
+ *
+ * The class describes TPC slice tracks used in AliHLTTPCGMMerger
+ */
+class AliHLTTPCGMSliceTrack
+{
+  
+ public:
+  
+  float Alpha()                      const { return fAlpha;      }
+  int   NClusters()                  const { return fNClusters;       }
+  int   PrevNeighbour()              const { return fPrevNeighbour;   }
+  int   NextNeighbour()              const { return fNextNeighbour;   }
+  int   SliceNeighbour()             const { return fSliceNeighbour; }
+  int   Used()                       const { return fUsed;            }
+  const AliHLTTPCCASliceOutTrack* OrigTrack() const { return fOrigTrack; }
+  float X()                      const { return fX;      }
+  float Y()                      const { return fY;      }
+  float Z()                      const { return fZ;      }
+  float SinPhi()                      const { return fSinPhi;      }
+  float CosPhi()                      const { return fCosPhi;      }
+  float SecPhi()                      const { return fSecPhi;      }
+  float DzDs()                      const { return fDzDs;      }
+  float QPt()                      const { return fQPt;      }
+
+  void Set( const AliHLTTPCCASliceOutTrack *sliceTr, float alpha ){
+    const AliHLTTPCCABaseTrackParam &t = sliceTr->Param();
+    fOrigTrack = sliceTr;
+    fX = t.GetX();
+    fY = t.GetY();
+    fZ = t.GetZ();
+    fDzDs = t.GetDzDs();
+    fSinPhi = t.GetSinPhi();
+    fQPt = t.GetQPt();
+    fCosPhi = sqrt(1.f - fSinPhi*fSinPhi);
+    fSecPhi = 1.f / fCosPhi;
+    fAlpha = alpha;
+    fNClusters = sliceTr->NClusters();    
+  }
+  
+  void SetNClusters ( int v )                        { fNClusters = v;       }
+  void SetPrevNeighbour( int v )                     { fPrevNeighbour = v;   }
+  void SetNextNeighbour( int v )                     { fNextNeighbour = v;   }
+  void SetUsed( int v )                             { fUsed = v;            }
+  void SetSliceNeighbour( int v )                    { fSliceNeighbour = v;            }
+
+
+  void CopyParamFrom( const AliHLTTPCGMSliceTrack &t ){
+    fX=t.fX; fY=t.fY; fZ=t.fZ;
+    fSinPhi=t.fSinPhi; fDzDs=t.fDzDs; fQPt=t.fQPt; fCosPhi=t.fCosPhi, fSecPhi=t.fSecPhi;
+    fAlpha = t.fAlpha;
+  }
+
+  bool FilterErrors( AliHLTTPCCAParam &param, float maxSinPhi =.999 );
+
+  bool TransportToX( float x, float Bz, AliHLTTPCGMBorderTrack &b, float maxSinPhi ) const ;
+
+  bool TransportToXAlpha( float x, float sinAlpha, float cosAlpha, float Bz, AliHLTTPCGMBorderTrack &b, float maxSinPhi ) const ;
+
+ private:
+
+  const AliHLTTPCCASliceOutTrack *fOrigTrack; // pointer to original slice track
+  float fX, fY, fZ, fSinPhi, fDzDs, fQPt, fCosPhi, fSecPhi; // parameters
+  float fC0, fC2, fC3, fC5, fC7, fC9, fC10, fC12, fC14; // covariances
+  float fAlpha;           // alpha angle 
+  int fNClusters;         // N clusters
+  int fPrevNeighbour;     // neighbour in the previous slise
+  int fNextNeighbour;     // neighbour in the next slise
+  int fSliceNeighbour;    // next neighbour withing the same slice;
+  int fUsed;              // is the slice track already merged
+  
+};
+
+#endif
diff --git a/HLT/TPCLib/merger-ca/AliHLTTPCGMTrackLinearisation.h b/HLT/TPCLib/merger-ca/AliHLTTPCGMTrackLinearisation.h
new file mode 100644 (file)
index 0000000..ed15975
--- /dev/null
@@ -0,0 +1,86 @@
+// ************************************************************************
+// This file is property of and copyright by the ALICE HLT Project        *
+// ALICE Experiment at CERN, All rights reserved.                         *
+// See cxx source for full Copyright notice                               *
+//                                                                        *
+//*************************************************************************
+
+
+#ifndef ALIHLTTPCGMTRACKLINEARISATION_H
+#define ALIHLTTPCGMTRACKLINEARISATION_H
+
+#include "AliHLTTPCGMTrackParam.h"
+
+
+/**
+ * @class AliHLTTPCGMTrackLinearisation
+ *
+ * AliHLTTPCGMTrackLinearisation class describes the parameters which are used
+ * to linearise the transport equations for the track trajectory.
+ *
+ * The class is used during track (re)fit, when the AliHLTTPCTrackParam track is only
+ * partially fitted, and there is some apriory knowledge about trajectory.
+ * This apriory knowledge is used to linearise the transport equations.
+ *
+ * In case the track is fully fitted, the best linearisation point is
+ * the track trajectory itself (AliHLTTPCGMTrackLinearisation = AliHLTTPCGMTrackParam ).
+ *
+ */
+class AliHLTTPCGMTrackLinearisation
+{
+
+  public:
+
+    AliHLTTPCGMTrackLinearisation()
+      : fSinPhi( 0. ), fCosPhi( 1. ), fSecPhi( 1. ), fDzDs( 0. ), fDlDs( 0.), fQPt( 0. ) {}
+
+    AliHLTTPCGMTrackLinearisation( float SinPhi1, float CosPhi1, float SecPhi1, float DzDs1, float DlDs1, float QPt1 )
+      : fSinPhi( SinPhi1 ), fCosPhi( CosPhi1 ), fSecPhi( SecPhi1 ), fDzDs( DzDs1 ), fDlDs( DlDs1), fQPt( QPt1 ) {}
+
+    AliHLTTPCGMTrackLinearisation( const AliHLTTPCGMTrackParam &t );
+
+    void Set( float SinPhi1, float CosPhi1, float SecPhi1, float DzDs1, float DlDs1, float QPt1 );
+
+    
+    float& SinPhi() { return fSinPhi; }
+    float& CosPhi() { return fCosPhi; }
+    float& SecPhi() { return fSecPhi; }
+    float& DzDs()   { return fDzDs; }
+    float& DlDs()   { return fDlDs; }
+    float& QPt()    { return fQPt; }
+
+
+ private:
+    float fSinPhi; // SinPhi
+    float fCosPhi; // CosPhi
+    float fSecPhi; // 1/cos(phi)
+    float fDzDs;   // DzDs
+    float fDlDs;   // DlDs
+    float fQPt;    // QPt
+};
+
+
+ inline AliHLTTPCGMTrackLinearisation::AliHLTTPCGMTrackLinearisation( const AliHLTTPCGMTrackParam &t )
+    : fSinPhi( t.GetSinPhi() ), fCosPhi( 0. ), fSecPhi( 0. ), fDzDs( t.GetDzDs() ), fDlDs( 0. ), fQPt( t.GetQPt() )
+{
+  fSinPhi = min( fSinPhi,  .999f );
+  fSinPhi = max( fSinPhi, -.999f );
+  fCosPhi = sqrt( 1. - fSinPhi * fSinPhi );
+  fSecPhi = 1./fCosPhi; //reciprocal(fCosPhi);
+  fDlDs = sqrt(1.+fDzDs*fDzDs);
+}
+
+
+ inline void AliHLTTPCGMTrackLinearisation::Set( float SinPhi1, float CosPhi1, float SecPhi1,
+    float DzDs1, float DlDs1, float QPt1 )
+{
+  fSinPhi = SinPhi1 ;
+  fCosPhi = CosPhi1 ;
+  fSecPhi = SecPhi1 ;
+  fDzDs = DzDs1;
+  fDlDs = DlDs1;
+  fQPt = QPt1;
+}
+
+
+#endif //ALIHLTTPCGMTRACKLINEARISATION_H
diff --git a/HLT/TPCLib/merger-ca/AliHLTTPCGMTrackParam.cxx b/HLT/TPCLib/merger-ca/AliHLTTPCGMTrackParam.cxx
new file mode 100644 (file)
index 0000000..435f773
--- /dev/null
@@ -0,0 +1,539 @@
+// $Id: AliHLTTPCGMTrackParam.cxx 41769 2010-06-16 13:58:00Z sgorbuno $
+// **************************************************************************
+// This file is property of and copyright by the ALICE HLT Project          *
+// ALICE Experiment at CERN, All rights reserved.                           *
+//                                                                          *
+// Primary Authors: Sergey Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de> *
+//                  for The ALICE HLT Project.                              *
+//                                                                          *
+// Permission to use, copy, modify and distribute this software and its     *
+// documentation strictly for non-commercial purposes is hereby granted     *
+// without fee, provided that the above copyright notice appears in all     *
+// copies and that both the copyright notice and this permission notice     *
+// appear in the supporting documentation. The authors make no claims       *
+// about the suitability of this software for any purpose. It is            *
+// provided "as is" without express or implied warranty.                    *
+//                                                                          *
+//***************************************************************************
+
+#include "AliHLTTPCGMTrackParam.h"
+#include "AliHLTTPCCAMath.h"
+#include "AliHLTTPCGMTrackLinearisation.h"
+#include "AliHLTTPCGMBorderTrack.h"
+#include "Riostream.h"
+#include "AliExternalTrackParam.h"
+#include "AliHLTTPCCAParam.h"
+
+
+float AliHLTTPCGMTrackParam::fPolinomialFieldBz[6];   // field coefficients
+
+
+void AliHLTTPCGMTrackParam::Fit
+(
+ float x[], float y[], float z[], unsigned int rowType[], float alpha[], AliHLTTPCCAParam &param,
+ int &N,
+ float &Alpha, 
+ bool UseMeanPt,
+ float maxSinPhi
+ ){
+  
+  const float kRho = 1.025e-3;//0.9e-3;
+  const float kRadLen = 29.532;//28.94;
+  const float kRhoOverRadLen = kRho / kRadLen;
+  
+  AliHLTTPCGMTrackLinearisation t0(*this);
+  const float kZLength = 250.f - 0.275f;
+  float trDzDs2 = t0.DzDs()*t0.DzDs();
+  AliHLTTPCGMTrackFitParam par;
+  CalculateFitParameters( par, kRhoOverRadLen, kRho, UseMeanPt );
+
+  int maxN = N;
+
+  bool first = 1;
+  N = 0;
+
+  for( int ihit=0; ihit<maxN; ihit++ ){
+    
+    float sliceAlpha =  alpha[ihit];
+    
+    if ( fabs( sliceAlpha - Alpha ) > 1.e-4 ) {
+      if( !Rotate(  sliceAlpha - Alpha, t0, .999 ) ) break;
+      Alpha = sliceAlpha;
+    }
+
+    float dL=0;    
+    float bz =  GetBz(x[ihit], y[ihit],z[ihit]);
+        
+    float err2Y, err2Z;
+
+    { // transport block
+      
+      bz = -bz;
+
+      float ex = t0.CosPhi();
+      
+      float ey = t0.SinPhi();
+      float k  = t0.QPt()*bz;
+      float dx = x[ihit] - X();
+      float kdx = k*dx;
+      float ey1 = kdx + ey;
+      
+      if( fabs( ey1 ) >= maxSinPhi ) break;
+
+      float ss = ey + ey1;   
+      float ex1 = sqrt(1 - ey1*ey1);
+      
+      float dxBz = dx * bz;
+    
+      float cc = ex + ex1;  
+      float dxcci = dx * Reciprocal(cc);
+      float kdx205 = kdx*kdx*0.5f;
+      
+      float dy = dxcci * ss;      
+      float norm2 = float(1.f) + ey*ey1 + ex*ex1;
+      float dl = dxcci * sqrt( norm2 + norm2 );
+
+      float dS;    
+      { 
+       float dSin = float(0.5f)*k*dl;
+       float a = dSin*dSin;
+       const float k2 = 1.f/6.f;
+       const float k4 = 3.f/40.f;
+       //const float k6 = 5.f/112.f;
+       dS = dl + dl*a*(k2 + a*(k4 ));//+ k6*a) );
+      }
+      
+      float ex1i = Reciprocal(ex1);
+      float dz = dS * t0.DzDs();  
+      
+      dL = -dS * t0.DlDs();
+      
+      float hh = dxcci*ex1i*(2.f+kdx205); 
+      float h2 = hh * t0.SecPhi();
+      float h4 = bz*dxcci*hh;
+
+      float d2 = fP[2] - t0.SinPhi();
+      float d3 = fP[3] - t0.DzDs();
+      float d4 = fP[4] - t0.QPt();
+      
+      
+      fX+=dx;
+      fP[0]+= dy     + h2 * d2           +   h4 * d4;
+      fP[1]+= dz               + dS * d3;
+      fP[2] = ey1 +     d2           + dxBz * d4;    
+      
+      t0.CosPhi() = ex1;
+      t0.SecPhi() = ex1i;
+      t0.SinPhi() = ey1;      
+
+      {
+       const float *cy = param.GetParamS0Par(0,rowType[ihit]);
+       const float *cz = param.GetParamS0Par(1,rowType[ihit]);
+
+       float secPhi2 = ex1i*ex1i;
+       float zz = fabs( kZLength - fabs(fP[2]) );      
+       float zz2 = zz*zz;
+       float angleY2 = secPhi2 - 1.f; 
+       float angleZ2 = trDzDs2 * secPhi2 ;
+
+       float cy0 = cy[0] + cy[1]*zz + cy[3]*zz2;
+       float cy1 = cy[2] + cy[5]*zz;
+       float cy2 = cy[4];
+       float cz0 = cz[0] + cz[1]*zz + cz[3]*zz2;
+       float cz1 = cz[2] + cz[5]*zz;
+       float cz2 = cz[4];
+       
+       err2Y = fabs( cy0 + angleY2 * ( cy1 + angleY2*cy2 ) );
+       err2Z = fabs( cz0 + angleZ2 * ( cz1 + angleZ2*cz2 ) );      
+     }
+
+
+      if ( first ) {
+       fP[0] = y[ihit];
+       fP[1] = z[ihit];
+       SetCov( 0, err2Y );
+       SetCov( 1,  0 );
+       SetCov( 2, err2Z);
+       SetCov( 3,  0 );
+       SetCov( 4,  0 );
+       SetCov( 5,  1 );
+       SetCov( 6,  0 );
+       SetCov( 7,  0 );
+       SetCov( 8,  0 );
+       SetCov( 9,  1 );
+       SetCov( 10,  0 );
+       SetCov( 11,  0 );
+       SetCov( 12,  0 );
+       SetCov( 13,  0 );
+       SetCov( 14,  10 );
+       SetChi2( 0 );
+       SetNDF( -3 );
+       CalculateFitParameters( par, kRhoOverRadLen, kRho, UseMeanPt );
+       first = 0;
+       N+=1;
+       continue;
+      }
+
+      float c20 = fC[3];
+      float c21 = fC[4];
+      float c22 = fC[5];
+      float c30 = fC[6];
+      float c31 = fC[7];
+      float c32 = fC[8];
+      float c33 = fC[9];
+      float c40 = fC[10];
+      float c41 = fC[11];
+      float c42 = fC[12];
+      float c43 = fC[13];
+      float c44 = fC[14];
+      
+      float c20ph4c42 =  c20 + h4*c42;
+      float h2c22 = h2*c22;
+      float h4c44 = h4*c44;
+      
+      float n6 = c30 + h2*c32 + h4*c43;
+      float n7 = c31 + dS*c33;
+      float n10 = c40 + h2*c42 + h4c44;
+      float n11 = c41 + dS*c43;
+      float n12 = c42 + dxBz*c44;
+      
+      fC[8] = c32 + dxBz * c43;
+      
+      fC[0]+= h2*h2c22 + h4*h4c44 + float(2.f)*( h2*c20ph4c42  + h4*c40 );
+      
+      fC[1]+= h2*c21 + h4*c41 + dS*n6;
+      fC[6] = n6;
+      
+      fC[2]+= dS*(c31 + n7);
+      fC[7] = n7; 
+      
+      fC[3] = c20ph4c42 + h2c22  + dxBz*n10;
+      fC[10] = n10;
+      
+      fC[4] = c21 + dS*c32 + dxBz*n11;
+      fC[11] = n11;
+      
+      fC[5] = c22 + dxBz*( c42 + n12 );
+      fC[12] = n12;
+      
+    } // end transport block 
+
+    float &fC22 = fC[5];
+    float &fC33 = fC[9];
+    float &fC40 = fC[10];
+    float &fC41 = fC[11];
+    float &fC42 = fC[12];
+    float &fC43 = fC[13];
+    float &fC44 = fC[14];
+    
+    float 
+      c00 = fC[ 0],
+      c11 = fC[ 2],
+      c20 = fC[ 3],
+      c31 = fC[ 7];
+    
+    
+    // MS block  
+    
+    float dLmask = 0.f;
+    bool maskMS = ( fabs( dL ) < par.fDLMax );
+
+    
+    // Filter block
+    
+    float mS0 = Reciprocal(err2Y + c00);
+    
+    // MS block
+    Assign( dLmask, maskMS, dL );
+    
+    // Filter block
+    
+    float  z0 = y[ihit] - fP[0];
+    float mS2 = Reciprocal(err2Z + c11);
+    
+    if( fabs( fP[2] + z0*c20*mS0  ) > maxSinPhi ) break;
+    
+    // MS block
+    
+    float dLabs = fabs( dLmask); 
+    float corr = float(1.f) - par.fEP2* dLmask ;
+    
+    fP[4]*= corr;
+    fC40 *= corr;
+    fC41 *= corr;
+    fC42 *= corr;
+    fC43 *= corr;
+    fC44  = fC44*corr*corr + dLabs*par.fSigmadE2;
+    
+    fC22 += dLabs * par.fK22 * (float(1.f)-fP[2]*fP[2]);
+    fC33 += dLabs * par.fK33;
+    fC43 += dLabs * par.fK43;
+        
+
+    // Filter block
+  
+    float c40 = fC40;
+    
+    // K = CHtS
+    
+    float k00, k11, k20, k31, k40;
+    
+    k00 = c00 * mS0;
+    k20 = c20 * mS0;
+    k40 = c40 * mS0;
+    fChi2  += mS0*z0*z0;
+    fP[0] += k00 * z0;
+    fP[2] += k20 * z0;
+    fP[4] += k40 * z0;
+    fC[ 0] -= k00 * c00 ;
+    fC[ 5] -= k20 * c20 ;
+    fC[10] -= k00 * c40 ;
+    fC[12] -= k40 * c20 ;
+    fC[ 3] -= k20 * c00 ;
+    fC[14] -= k40 * c40 ;
+  
+    float  z1 = z[ihit] - fP[1];
+    
+    k11 = c11 * mS2;
+    k31 = c31 * mS2;
+    
+    fChi2  +=  mS2*z1*z1 ;
+    fNDF  += 2;
+    N+=1;
+    
+    fP[1] += k11 * z1;
+    fP[3] += k31 * z1;
+    
+    fC[ 7] -= k31 * c11;
+    fC[ 2] -= k11 * c11;
+    fC[ 9] -= k31 * c31;    
+  } 
+}
+
+
+
+
+
+bool AliHLTTPCGMTrackParam::CheckNumericalQuality() const
+{
+  //* Check that the track parameters and covariance matrix are reasonable
+
+  bool ok = finite(fX) && finite( fChi2 ) && finite( fNDF );
+
+  const float *c = fC;
+  for ( int i = 0; i < 15; i++ ) ok = ok && finite( c[i] );
+  for ( int i = 0; i < 5; i++ ) ok = ok && finite( fP[i] );
+
+  if ( c[0] <= 0 || c[2] <= 0 || c[5] <= 0 || c[9] <= 0 || c[14] <= 0 ) ok = 0;
+  if ( c[0] > 5. || c[2] > 5. || c[5] > 2. || c[9] > 2. 
+       //|| ( CAMath::Abs( QPt() ) > 1.e-2 && c[14] > 2. ) 
+       ) ok = 0;
+
+  if ( fabs( fP[2] ) > .999 ) ok = 0;
+  if ( fabs( fP[4] ) > 1. / 0.05 ) ok = 0;
+  if( ok ){
+    ok = ok 
+      && ( c[1]*c[1]<=c[2]*c[0] )
+      && ( c[3]*c[3]<=c[5]*c[0] )
+      && ( c[4]*c[4]<=c[5]*c[2] )
+      && ( c[6]*c[6]<=c[9]*c[0] )
+      && ( c[7]*c[7]<=c[9]*c[2] )
+      && ( c[8]*c[8]<=c[9]*c[5] )
+      && ( c[10]*c[10]<=c[14]*c[0] )
+      && ( c[11]*c[11]<=c[14]*c[2] )
+      && ( c[12]*c[12]<=c[14]*c[5] )
+      && ( c[13]*c[13]<=c[14]*c[9] );      
+  }
+  return ok;
+}
+
+
+
+
+//*
+//*  Multiple scattering and energy losses
+//*
+
+float AliHLTTPCGMTrackParam::ApproximateBetheBloch( float beta2 )
+{
+  //------------------------------------------------------------------
+  // This is an approximation of the Bethe-Bloch formula with
+  // the density effect taken into account at beta*gamma > 3.5
+  // (the approximation is reasonable only for solid materials)
+  //------------------------------------------------------------------
+
+  const float log0 = log( float(5940.f));
+  const float log1 = log( float(3.5f*5940.f) );
+
+  bool bad = (beta2 >= .999f)||( beta2 < 1.e-8f );
+
+  Assign( beta2, bad, 0.5f);
+
+  float a = beta2 / ( 1.f - beta2 ); 
+  float b = 0.5*log(a);
+  float d =  0.153e-3 / beta2;
+  float c = b - beta2;
+
+  float ret = d*(log0 + b + c );
+  float case1 = d*(log1 + c );
+  
+  Assign( ret, ( a > 3.5*3.5  ), case1);
+  Assign( ret,  bad, 0. ); 
+
+  return ret;
+}
+
+
+ void AliHLTTPCGMTrackParam::CalculateFitParameters( AliHLTTPCGMTrackFitParam &par, float RhoOverRadLen,  float Rho, bool NoField, float mass )
+{
+  //*!
+
+  float qpt = fP[4];
+  if( NoField ) qpt = 1./0.35;
+
+  float p2 = ( 1. + fP[3] * fP[3] );
+  float k2 = qpt * qpt;
+  Assign( k2, (  k2 < 1.e-4f ), 1.e-4f );
+
+  float mass2 = mass * mass;
+  float beta2 = p2 / ( p2 + mass2 * k2 );
+  
+  float pp2 = p2 / k2; // impuls 2
+
+  //par.fBethe = BetheBlochGas( pp2/mass2);
+  par.fBetheRho = ApproximateBetheBloch( pp2 / mass2 )*Rho;
+  par.fE = sqrt( pp2 + mass2 );
+  par.fTheta2 = ( 14.1*14.1/1.e6 ) / ( beta2 * pp2 )*RhoOverRadLen;
+  par.fEP2 = par.fE / pp2;
+
+  // Approximate energy loss fluctuation (M.Ivanov)
+
+  const float knst = 0.07; // To be tuned.
+  par.fSigmadE2 = knst * par.fEP2 * qpt;
+  par.fSigmadE2 = par.fSigmadE2 * par.fSigmadE2;
+  
+  float k22 = 1. + fP[3] * fP[3];
+  par.fK22 = par.fTheta2*k22;
+  par.fK33 = par.fK22 * k22;
+  par.fK43 = 0.;
+  par.fK44 =  par.fTheta2*fP[3] * fP[3] * k2;
+  
+  float br=1.e-8f;
+  Assign( br, ( par.fBetheRho>1.e-8f ), par.fBetheRho );
+  par.fDLMax = 0.3*par.fE * Reciprocal( br );
+
+  par.fEP2*=par.fBetheRho;
+  par.fSigmadE2 = par.fSigmadE2*par.fBetheRho+par.fK44;  
+}
+
+
+
+
+//*
+//* Rotation
+//*
+
+
+bool AliHLTTPCGMTrackParam::Rotate( float alpha, AliHLTTPCGMTrackLinearisation &t0, float maxSinPhi )
+{
+  //* Rotate the coordinate system in XY on the angle alpha
+
+  float cA = CAMath::Cos( alpha );
+  float sA = CAMath::Sin( alpha );
+  float x0 = X(), y0 = Y(), sP = t0.SinPhi(), cP = t0.CosPhi();
+  float cosPhi = cP * cA + sP * sA;
+  float sinPhi = -cP * sA + sP * cA;
+
+  if ( CAMath::Abs( sinPhi ) > maxSinPhi || CAMath::Abs( cosPhi ) < 1.e-2 || CAMath::Abs( cP ) < 1.e-2  ) return 0;
+
+  //float J[5][5] = { { j0, 0, 0,  0,  0 }, // Y
+  //                    {  0, 1, 0,  0,  0 }, // Z
+  //                    {  0, 0, j2, 0,  0 }, // SinPhi
+  //                  {  0, 0, 0,  1,  0 }, // DzDs
+  //                  {  0, 0, 0,  0,  1 } }; // Kappa
+
+  float j0 = cP / cosPhi;
+  float j2 = cosPhi / cP;
+  float d[2] = {Y() - y0, SinPhi() - sP};
+
+  X() = ( x0*cA +  y0*sA );
+  Y() = ( -x0*sA +  y0*cA + j0*d[0] );
+  t0.CosPhi() = fabs( cosPhi );
+  t0.SecPhi() = ( 1./t0.CosPhi() );
+  t0.SinPhi() = ( sinPhi );
+
+  SinPhi() = ( sinPhi + j2*d[1] );
+
+  fC[0] *= j0 * j0;
+  fC[1] *= j0;
+  fC[3] *= j0;
+  fC[6] *= j0;
+  fC[10] *= j0;
+
+  fC[3] *= j2;
+  fC[4] *= j2;
+  fC[5] *= j2 * j2;
+  fC[8] *= j2;
+  fC[12] *= j2;
+  if( cosPhi <0 ){ // change direction
+    t0.SinPhi() = -sinPhi;
+    t0.DzDs() = -t0.DzDs();
+    t0.DlDs() = -t0.DlDs();
+    t0.QPt() = -t0.QPt();
+    SinPhi() = -SinPhi();
+    DzDs() = -DzDs();
+    QPt() = -QPt();
+    fC[3] = -fC[3];
+    fC[4] = -fC[4];
+    fC[6] = -fC[6];
+    fC[7] = -fC[7];
+    fC[10] = -fC[10];
+    fC[11] = -fC[11];
+  }
+
+  return 1;
+}
+
+
+
+
+bool AliHLTTPCGMTrackParam::GetExtParam( AliExternalTrackParam &T, double alpha ) const
+{
+  //* Convert from AliHLTTPCGMTrackParam to AliExternalTrackParam parameterisation,
+  //* the angle alpha is the global angle of the local X axis
+
+  bool ok = CheckNumericalQuality();
+
+  double par[5], cov[15];
+  for ( int i = 0; i < 5; i++ ) par[i] = fP[i];
+  for ( int i = 0; i < 15; i++ ) cov[i] = fC[i];
+
+  if ( par[2] > .99 ) par[2] = .99;
+  if ( par[2] < -.99 ) par[2] = -.99;
+
+  if ( fabs( par[4] ) < 1.e-5 ) par[4] = 1.e-5; // some other software will crash if q/Pt==0
+  if ( fabs( par[4] ) > 1./0.08 ) ok = 0; // some other software will crash if q/Pt is too big
+
+  T.Set( (double) fX, alpha, par, cov );
+  return ok;
+}
+
+
+void AliHLTTPCGMTrackParam::SetExtParam( const AliExternalTrackParam &T )
+{
+  //* Convert from AliExternalTrackParam parameterisation
+
+  for ( int i = 0; i < 5; i++ ) fP[i] = T.GetParameter()[i];
+  for ( int i = 0; i < 15; i++ ) fC[i] = T.GetCovariance()[i];
+  fX = T.GetX();
+  if ( fP[2] > .999 ) fP[2] = .999;
+  if ( fP[2] < -.999 ) fP[2] = -.999;
+}
+
+
+
+
diff --git a/HLT/TPCLib/merger-ca/AliHLTTPCGMTrackParam.h b/HLT/TPCLib/merger-ca/AliHLTTPCGMTrackParam.h
new file mode 100644 (file)
index 0000000..a4e964e
--- /dev/null
@@ -0,0 +1,143 @@
+//-*- Mode: C++ -*-
+// $Id: AliHLTTPCGMTrackParam.h 39008 2010-02-18 17:33:32Z sgorbuno $
+// ************************************************************************
+// This file is property of and copyright by the ALICE HLT Project        *
+// ALICE Experiment at CERN, All rights reserved.                         *
+// See cxx source for full Copyright notice                               *
+//                                                                        *
+//*************************************************************************
+
+
+#ifndef ALIHLTTPCGMTRACKPARAM_H
+#define ALIHLTTPCGMTRACKPARAM_H
+
+#include "TMath.h"
+
+class AliHLTTPCGMTrackLinearisation;
+class AliHLTTPCGMBorderTrack;
+class AliExternalTrackParam;
+class AliHLTTPCCAParam;
+
+/**
+ * @class AliHLTTPCGMTrackParam
+ *
+ * AliHLTTPCGMTrackParam class describes the track parametrisation
+ * which is used by the AliHLTTPCGMTracker slice tracker.
+ *
+ */
+class AliHLTTPCGMTrackParam
+{
+public:
+
+  struct AliHLTTPCGMTrackFitParam {
+    //float fBethe, fE, fTheta2, fEP2, fSigmadE2, fK22, fK33, fK43, fK44;// parameters
+    float fDLMax, fBetheRho, fE, fTheta2, fEP2, fSigmadE2, fK22, fK33, fK43, fK44;// parameters
+  };
+    
+  float& X()      { return fX;    }
+  float& Y()      { return fP[0]; }
+  float& Z()      { return fP[1]; }
+  float& SinPhi() { return fP[2]; }
+  float& DzDs()   { return fP[3]; }
+  float& QPt()    { return fP[4]; }
+  
+  float GetX()      const { return fX; }
+  float GetY()      const { return fP[0]; }
+  float GetZ()      const { return fP[1]; }
+  float GetSinPhi() const { return fP[2]; }
+  float GetDzDs()   const { return fP[3]; }
+  float GetQPt()    const { return fP[4]; }
+
+  float GetKappa( float Bz ) const { return -fP[4]*Bz; }
+
+  void SetX( float v ){ fX = v; }
+
+  float *Par() { return fP; }
+  const float *GetPar() const { return fP; }
+  float GetPar( int i) const { return(fP[i]); }
+  void SetPar( int i, float v ) { fP[i] = v; }
+
+  float& Chi2()  { return fChi2; }
+  int&   NDF()  { return fNDF; }
+  
+  float Err2Y()      const { return fC[0]; }
+  float Err2Z()      const { return fC[2]; }
+  float Err2SinPhi() const { return fC[5]; }
+  float Err2DzDs()   const { return fC[9]; }
+  float Err2QPt()    const { return fC[14]; }
+  
+  float GetChi2()   const { return fChi2; }
+  int   GetNDF()    const { return fNDF; }
+
+  float GetCosPhi() const { return sqrt( float(1.) - GetSinPhi()*GetSinPhi() ); }
+
+  float GetErr2Y()      const { return fC[0]; }
+  float GetErr2Z()      const { return fC[2]; }
+  float GetErr2SinPhi() const { return fC[5]; }
+  float GetErr2DzDs()   const { return fC[9]; }
+  float GetErr2QPt()    const { return fC[14]; }
+
+  float *Cov() { return fC; }
+
+  const float *GetCov() const { return fC; }
+  float GetCov(int i) const {return fC[i]; }
+
+
+  void SetCov( int i, float v ) { fC[i] = v; }
+  void SetChi2( float v )  {  fChi2 = v; }
+  void SetNDF( int v )   { fNDF = v; }
+  
+
+  static float ApproximateBetheBloch( float beta2 );
+
+  void CalculateFitParameters( AliHLTTPCGMTrackFitParam &par,float RhoOverRadLen,  float Rho,  bool NoField=0, float mass = 0.13957 );
+    
+
+  bool CheckNumericalQuality() const ;
+
+  void Fit
+  (
+   float x[], float y[], float z[], unsigned int rowType[], float alpha[], AliHLTTPCCAParam &param,
+   int &N, float &Alpha, 
+   bool UseMeanPt = 0,
+   float maxSinPhi = .999
+   );
+  
+  bool Rotate( float alpha, AliHLTTPCGMTrackLinearisation &t0, float maxSinPhi = .999 );
+  
+
+  static float fPolinomialFieldBz[6];   // field coefficients
+  static float GetBz( float x, float y, float z );
+  float GetBz() const{ return GetBz( fX, fP[0], fP[1] );}
+
+  static float Reciprocal( float x ){ return 1./x; }
+  static void Assign( float &x, bool mask, float v ){
+    if( mask ) x = v;
+  }
+
+  static void Assign( int &x, bool mask, int v ){
+    if( mask ) x = v;
+  }
+
+  bool GetExtParam( AliExternalTrackParam &T, double alpha ) const;
+  void SetExtParam( const AliExternalTrackParam &T );
+
+  private:
+  
+    float fX;      // x position
+    float fP[5];   // 'active' track parameters: Y, Z, SinPhi, DzDs, q/Pt
+    float fC[15];  // the covariance matrix for Y,Z,SinPhi,..
+    float fChi2;   // the chi^2 value
+    int   fNDF;    // the Number of Degrees of Freedom
+};
+
+
+inline  float AliHLTTPCGMTrackParam::GetBz( float x, float y, float z ) 
+{
+  float r2 = x * x + y * y;
+  float r  = sqrt( r2 );
+  const float *c = fPolinomialFieldBz;
+  return ( c[0] + c[1]*z  + c[2]*r  + c[3]*z*z + c[4]*z*r + c[5]*r2 );
+}
+
+#endif //ALIHLTTPCCATRACKPARAM_H
index 2873cb5..8c43b0e 100644 (file)
@@ -84,14 +84,12 @@ class AliHLTTPCCAMerger::AliHLTTPCCABorderTrack
     const AliHLTTPCCATrackParam &Param() const { return fParam;     }
     int   TrackID()                    const { return fTrackID;   }
     int   NClusters()                  const { return fNClusters; }
-    int   IRow()                       const { return fIRow;      }
     float X()                          const { return fX;         }
     bool  OK()                         const { return fOK;        }
 
     void SetParam     ( const AliHLTTPCCATrackParam &v ) { fParam     = v; }
     void SetTrackID   ( int v )                        { fTrackID   = v; }
     void SetNClusters ( int v )                        { fNClusters = v; }
-    void SetIRow      ( int v )                        { fIRow      = v; }
     void SetX         ( float v )                      { fX         = v; }
     void SetOK        ( bool v )                       { fOK        = v; }
 
@@ -100,7 +98,6 @@ class AliHLTTPCCAMerger::AliHLTTPCCABorderTrack
     AliHLTTPCCATrackParam fParam;  // track parameters at the border
     int   fTrackID;              // track index
     int   fNClusters;            // n clusters
-    int   fIRow;                 // row number of the closest cluster
     float fX;                    // X coordinate of the closest cluster
     bool  fOK;                   // is the track rotated and extrapolated correctly
 
@@ -240,22 +237,20 @@ void AliHLTTPCCAMerger::UnpackSlices()
         // unpack cluster information
 
         AliHLTTPCCAClusterInfo &clu = fClusterInfos[nClustersCurrent + nCluNew];
-       UInt_t id, row;
-       float x,y,z;
-       sliceTr->Cluster( iTrClu ).Get(iSlice,id,row,x,y,z);
+       const AliHLTTPCCASliceOutCluster &c = sliceTr->Cluster( iTrClu );
        
         clu.SetISlice( iSlice );
-        clu.SetIRow( row );
-        clu.SetId( id );
+        clu.SetRowType( c.GetRowType() );
+        clu.SetId( c.GetId() );
         clu.SetPackedAmp( 0 );
-        clu.SetX( x );
-        clu.SetY( y );
-        clu.SetZ( z );
+        clu.SetX( c.GetX() );
+        clu.SetY( c.GetY() );
+        clu.SetZ( c.GetZ() );
 
         if ( !t0.TransportToX( clu.X(), fSliceParam.GetBz( t0 ), .999 ) ) continue;
 
         float err2Y, err2Z;
-        fSliceParam.GetClusterErrors2( clu.IRow(), clu.Z(), t0.SinPhi(), t0.GetCosPhi(), t0.DzDs(), err2Y, err2Z );
+        fSliceParam.GetClusterErrors2v1( clu.RowType(), clu.Z(), t0.SinPhi(), t0.GetCosPhi(), t0.DzDs(), err2Y, err2Z );
 
         clu.SetErr2Y( err2Y );
         clu.SetErr2Z( err2Z );
@@ -378,7 +373,7 @@ bool AliHLTTPCCAMerger::FitTrack( AliHLTTPCCATrackParam &T, float &Alpha,
 
     float err2Y = h.Err2Y();
     float err2Z = h.Err2Z();
-    if ( doErrors ) fSliceParam.GetClusterErrors2( h.IRow(), h.Z(), l.SinPhi(), l.CosPhi(), l.DzDs(), err2Y, err2Z );
+    if ( doErrors ) fSliceParam.GetClusterErrors2v1( h.RowType(), h.Z(), l.SinPhi(), l.CosPhi(), l.DzDs(), err2Y, err2Z );
     if( !final ){
       err2Y*= fSliceParam.ClusterError2CorrectionY();
       err2Z*= fSliceParam.ClusterError2CorrectionZ();
@@ -486,8 +481,7 @@ void AliHLTTPCCAMerger::MakeBorderTracks( int iSlice, int iBorder, AliHLTTPCCABo
       if ( t0.TransportToX( x0, fSliceParam.GetBz( t0 ), maxSin ) ) {
         b.SetOK( 1 );
         b.SetTrackID( itr );
-        b.SetNClusters( track.NClusters() );
-        b.SetIRow( fClusterInfos[ track.FirstClusterRef() + 0 ].IRow() );
+        b.SetNClusters( track.NClusters() );        
         b.SetParam( t0 );
         nB++;
       }
@@ -500,7 +494,6 @@ void AliHLTTPCCAMerger::MakeBorderTracks( int iSlice, int iBorder, AliHLTTPCCABo
         b.SetOK( 1 );
         b.SetTrackID( itr );
         b.SetNClusters( track.NClusters() );
-        b.SetIRow( fClusterInfos[ track.FirstClusterRef() + track.NClusters()-1 ].IRow() );
         b.SetParam( t1 );
         nB++;
       }
index 95c2a45..bf13992 100644 (file)
@@ -37,7 +37,7 @@ class AliHLTTPCCAMerger
       public:
 
         unsigned char  ISlice()    const { return fISlice;    }
-        unsigned char  IRow()    const { return fIRow;    }
+        unsigned char  RowType()    const { return fRowType;    }
         int  Id()      const { return fId;      }
         UChar_t PackedAmp() const { return fPackedAmp; }
         float X()         const { return fX;         }
@@ -47,7 +47,7 @@ class AliHLTTPCCAMerger
         float Err2Z()     const { return fErr2Z;     }
 
         void SetISlice    ( unsigned char v  ) { fISlice    = v; }
-        void SetIRow    ( unsigned char v  ) { fIRow    = v; }
+        void SetRowType    ( unsigned char v  ) { fRowType    = v; }
         void SetId      (  int v  ) { fId      = v; }
         void SetPackedAmp ( UChar_t v ) { fPackedAmp = v; }
         void SetX         ( float v ) { fX         = v; }
@@ -59,7 +59,7 @@ class AliHLTTPCCAMerger
       private:
 
         unsigned char fISlice;            // slice number
-        unsigned char fIRow;            // row number
+        unsigned char fRowType;            // row type
         int fId;                 // cluster hlt number
         UChar_t fPackedAmp; // packed cluster amplitude
         float fX;                // x position (slice coord.system)
index 2ee748c..ad7e594 100644 (file)
@@ -180,6 +180,21 @@ GPUdi() void AliHLTTPCCAParam::GetClusterErrors2( int iRow, float z, float sinPh
   Err2Z = GetClusterError2( 1, type, z, angleZ );
 }
 
+GPUdi() void AliHLTTPCCAParam::GetClusterErrors2v1( int rowType, float z, float sinPhi, float cosPhi, float DzDs, float &Err2Y, float &Err2Z ) const
+{
+  //
+  // Use calibrated cluster error from OCDB
+  //
+
+  z = CAMath::Abs( ( 250. - 0.275 ) - CAMath::Abs( z ) );
+  float cosPhiInv = CAMath::Abs( cosPhi ) > 1.e-2 ? 1. / cosPhi : 0;
+  float angleY = sinPhi * cosPhiInv ; // dy/dx
+  float angleZ = DzDs * cosPhiInv ; // dz/dx
+
+  Err2Y = GetClusterError2( 0, rowType, z, angleY );
+  Err2Z = GetClusterError2( 1, rowType, z, angleZ );
+}
+
 #ifndef HLTCA_GPUCODE
 GPUh() void AliHLTTPCCAParam::WriteSettings( std::ostream &out ) const
 {
index bc980fc..921f573 100644 (file)
@@ -112,7 +112,8 @@ class AliHLTTPCCAParam
   GPUd() void SetMinTrackPt( float v ){ fMaxTrackQPt = 1./CAMath::Abs(v); }
 
     GPUd() float GetClusterError2( int yz, int type, float z, float angle ) const;
-    GPUd() void GetClusterErrors2( int iRow, float z, float sinPhi, float cosPhi, float DzDs, float &Err2Y, float &Err2Z ) const;
+    GPUd() void GetClusterErrors2( int row, float z, float sinPhi, float cosPhi, float DzDs, float &Err2Y, float &Err2Z ) const;
+    GPUd() void GetClusterErrors2v1( int rowType, float z, float sinPhi, float cosPhi, float DzDs, float &Err2Y, float &Err2Z ) const;
 
     void WriteSettings( std::ostream &out ) const;
     void ReadSettings( std::istream &in );
@@ -120,7 +121,9 @@ class AliHLTTPCCAParam
     GPUd() void SetParamS0Par( int i, int j, int k, float val ) {
       fParamS0Par[i][j][k] = val;
     }
-
+  
+    GPUd() const float *GetParamS0Par(int i, int j) const { return fParamS0Par[i][j]; }
     GPUd() float GetBzkG() const { return fBzkG;}
     GPUd() float GetConstBz() const { return fConstBz;}
     GPUd() float GetBz( float x, float y, float z ) const;
index 1d9b853..91f42c3 100644 (file)
@@ -1176,10 +1176,10 @@ void AliHLTTPCCAPerformance::SliceTrackCandPerformance( int /*iSlice*/, bool /*P
 
 
 
-void AliHLTTPCCAPerformance::SlicePerformance( int iSlice, bool PrintFlag )
+void AliHLTTPCCAPerformance::SlicePerformance( int /*iSlice*/, bool /*PrintFlag*/ )
 {
   //* calculate slice tracker performance
-
+#ifdef XXX
   AliHLTTPCCAStandaloneFramework &hlt = AliHLTTPCCAStandaloneFramework::Instance();
 
   int nRecTot = 0, nGhost = 0, nRecOut = 0;
@@ -1326,6 +1326,7 @@ void AliHLTTPCCAPerformance::SlicePerformance( int iSlice, bool PrintFlag )
     if ( nRecTot > 0 ) cout << nGhost / dRecTot; else cout << "_";
     cout << endl;
   }
+#endif
 }
 
 
@@ -1333,7 +1334,7 @@ void AliHLTTPCCAPerformance::SlicePerformance( int iSlice, bool PrintFlag )
 void AliHLTTPCCAPerformance::MergerPerformance()
 {
   // performance calculation for merged tracks
-
+#ifdef XXX
   int nRecTot = 0, nGhost = 0, nRecOut = 0;
   int nMCAll = 0, nRecAll = 0, nClonesAll = 0;
   int nMCRef = 0, nRecRef = 0, nClonesRef = 0;
@@ -1518,6 +1519,7 @@ void AliHLTTPCCAPerformance::MergerPerformance()
   fStatGBNMCRef  += nMCRef;
   fStatGBNRecRef  += nRecRef;
   fStatGBNClonesRef  += nClonesRef;
+#endif
 }
 
 
index ea78256..c380cf6 100644 (file)
@@ -23,23 +23,32 @@ class AliHLTTPCCASliceOutCluster
   public:
 
   GPUh() void Set( UInt_t id, UInt_t row, float x, float y, float z ){
-    UInt_t patch = (id>>1)&(0x7<<21); 
-    UInt_t cluster = id&0x1fffff;
-    fId = (row<<24)+ patch + cluster;
-    fXYZp = AliHLTTPCCADataCompressor::PackXYZ( row, x, y, z );
+    UInt_t rowtype;
+    //if( row<64 ) rowtype = 0;
+    //else if( row<128 ) rowtype = (UInt_t(2)<<30);
+    //else rowtype = (1<<30);
+    //fId = id|rowtype;
+    if( row<64 ) rowtype = 0;
+    else if( row<128 ) rowtype = 2;
+    else rowtype = 1;
+    fRowType = rowtype;
+    fId = id;
+    fX = x; fY = y; fZ = z;
   }
 
-  GPUh() void Get( int iSlice, UInt_t &Id, UInt_t &row, float &x, float &y, float &z ) const{
-    row = fId>>24;
-    UInt_t patch = (fId<<1)&(0x7<<22);
-    UInt_t cluster = fId&0x1fffff;
-    Id = (iSlice<<25) + patch + cluster;  
-    AliHLTTPCCADataCompressor::UnpackXYZ( row, fXYZp, x, y, z  );
-  }  
-  
+  GPUh() float GetX() const {return fX;}
+  GPUh() float GetY() const {return fY;}
+  GPUh() float GetZ() const {return fZ;}
+  GPUh() UInt_t GetId() const {return fId; } //fId & 0x3FFFFFFF;}
+  GPUh() UInt_t GetRowType() const {return fRowType; }//fId>>30;}
+
   private:
-    UInt_t fId; // Id ( row, patch, cluster )
-    AliHLTTPCCACompressedCluster fXYZp;// packed coordinates
+
+  UInt_t  fId; // Id ( slice, patch, cluster )    
+  UInt_t  fRowType; // row type
+  Float_t fX;// coordinates
+  Float_t fY;// coordinates
+  Float_t fZ;// coordinates
 };
 
 #endif 
index cc7b5c6..23e93ff 100644 (file)
@@ -33,6 +33,7 @@ class AliHLTTPCCASliceOutTrack
     GPUhd() int NClusters()                    const { return fNClusters;       }
     GPUhd() const AliHLTTPCCABaseTrackParam &Param() const { return fParam;           }
     GPUhd() const AliHLTTPCCASliceOutCluster &Cluster( int i ) const { return fClusters[i];           }
+    GPUhd() const AliHLTTPCCASliceOutCluster* Clusters() const { return fClusters;           }
 
     GPUhd() void SetNClusters( int v )                   { fNClusters = v;       }
     GPUhd() void SetParam( const AliHLTTPCCABaseTrackParam &v ) { fParam = v;           }
index 16e1b7e..f31ee67 100644 (file)
@@ -19,8 +19,6 @@
 
 #include "AliHLTTPCCAStandaloneFramework.h"
 #include "AliHLTTPCCATrackParam.h"
-#include "AliHLTTPCCAMerger.h"
-#include "AliHLTTPCCAMergerOutput.h"
 #include "AliHLTTPCCADataCompressor.h"
 #include "AliHLTTPCCAMath.h"
 #include "AliHLTTPCCAClusterData.h"
index b096256..e5022bc 100644 (file)
@@ -10,7 +10,7 @@
 #define ALIHLTTPCCASTANDALONEFRAMEWORK_H
 
 #include "AliHLTTPCCADef.h"
-#include "AliHLTTPCCAMerger.h"
+#include "AliHLTTPCGMMerger.h"
 #include "AliHLTTPCCAClusterData.h"
 #include "AliHLTTPCCATrackerFramework.h"
 #include <iostream>
@@ -41,7 +41,7 @@ class AliHLTTPCCAStandaloneFramework
        const AliHLTTPCCAParam &Param ( int iSlice ) const { return(fTracker.Param(iSlice)); }
        const AliHLTTPCCARow &Row ( int iSlice, int iRow ) const { return(fTracker.Row(iSlice, iRow)); }
     const AliHLTTPCCASliceOutput &Output( int iSlice ) const { return *fSliceOutput[iSlice]; }
-    AliHLTTPCCAMerger  &Merger()  { return fMerger; }
+    AliHLTTPCGMMerger  &Merger()  { return fMerger; }
     AliHLTTPCCAClusterData &ClusterData( int iSlice ) { return fClusterData[iSlice]; }
 
     /**
@@ -105,7 +105,7 @@ class AliHLTTPCCAStandaloneFramework
     AliHLTTPCCAStandaloneFramework( const AliHLTTPCCAStandaloneFramework& );
     const AliHLTTPCCAStandaloneFramework &operator=( const AliHLTTPCCAStandaloneFramework& ) const;
 
-    AliHLTTPCCAMerger fMerger;  //* global merger
+    AliHLTTPCGMMerger fMerger;  //* global merger
     AliHLTTPCCAClusterData fClusterData[fgkNSlices];
        AliHLTTPCCASliceOutput* fSliceOutput[fgkNSlices];
        AliHLTTPCCASliceOutput::outputControlStruct fOutputControl;
index 728d261..cfabf24 100644 (file)
@@ -327,14 +327,11 @@ int AliHLTTPCCATrackerOutputConverter::DoEvent( const AliHLTComponentEventData &
       currOutTrack->fFlags = 0;
       currOutTrack->fNPoints = nClu;    
       for( int i = 0; i< nClu; i++ ) { 
-       UInt_t id, row;
-       float x,y,z;
-       sliceTr->Cluster( i ).Get(slice,id,row,x,y,z);      
-       currOutTrack->fPointIDs[i] = id;
+       currOutTrack->fPointIDs[i] = sliceTr->Cluster( i ).GetId();
        if( i == nClu-1 ){
-         currOutTrack->fLastX = x;
-         currOutTrack->fLastY = y;
-         currOutTrack->fLastZ = z;
+         currOutTrack->fLastX = sliceTr->Cluster( i ).GetX();
+         currOutTrack->fLastY = sliceTr->Cluster( i ).GetY();
+         currOutTrack->fLastZ = sliceTr->Cluster( i ).GetZ();
        }
       }
       currOutTrack = ( AliHLTExternalTrackParam* )( (( Byte_t * )currOutTrack) + dSize );
index 0366660..5c5ad41 100644 (file)
@@ -32,9 +32,9 @@
 
 #include "AliHLTTPCCAPerformance.h"
 #include "AliHLTTPCCAParam.h"
-#include "AliHLTTPCCATrackConvertor.h"
 #include "AliHLTTPCCATracker.h"
 #include "AliHLTTPCCAStandaloneFramework.h"
+#include "AliHLTTPCGMMergedTrack.h"
 
 #include "TMath.h"
 #include "AliTPCLoader.h"
@@ -49,7 +49,6 @@
 #include "AliTrackReference.h"
 #include "TStopwatch.h"
 #include "AliTPCReconstructor.h"
-#include "AliHLTTPCCAMergerOutput.h"
 
 //#include <fstream.h>
 
@@ -134,7 +133,7 @@ AliTPCtrackerCA::AliTPCtrackerCA( const AliTPCParam *par ):
     param.SetMaxTrackMatchDRow( 5 );
     param.SetTrackConnectionFactor( 3.5 );
     param.SetMinNTrackClusters( 30 );
-    param.SetMinTrackPt(0.3);
+    param.SetMinTrackPt(0.2);
     AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
     for ( int iRow = 0; iRow < nRows; iRow++ ) {
       int    type = ( iRow < 63 ) ? 0 : ( ( iRow > 126 ) ? 1 : 2 );
@@ -378,18 +377,16 @@ int AliTPCtrackerCA::Clusters2Tracks( AliESDEvent *event )
 
   hlt.ProcessEvent();
 
-  if ( event ) {
+  if ( event ) {  
 
-    const AliHLTTPCCAMergerOutput &hltOut = *( hlt.Merger().Output() );
-
-    for ( int itr = 0; itr < hltOut.NTracks(); itr++ ) {
+    for ( int itr = 0; itr < hlt.Merger().NOutputTracks(); itr++ ) {
 
       AliTPCseed tTPC;
 
-      const AliHLTTPCCAMergedTrack &tCA = hltOut.Track( itr );
+      const AliHLTTPCGMMergedTrack &tCA = hlt.Merger().OutputTracks()[itr];
 
-      AliHLTTPCCATrackParam par = tCA.InnerParam();
-      AliHLTTPCCATrackConvertor::GetExtParam( par, tTPC, tCA.InnerAlpha() );
+      AliHLTTPCGMTrackParam par = tCA.GetParam();      
+      par.GetExtParam( tTPC, tCA.GetAlpha() );
       tTPC.SetMass( 0.13957 );
       if ( TMath::Abs( tTPC.GetSigned1Pt() ) > 1. / 0.02 ) continue;
       int nhits = tCA.NClusters();
@@ -399,17 +396,17 @@ int AliTPCtrackerCA::Clusters2Tracks( AliESDEvent *event )
         nhits = 160;
       }
       tTPC.SetNumberOfClusters( nhits );
-      float alpha = tCA.InnerAlpha();
-      AliHLTTPCCATrackParam t0 = par;
+      //float alpha = tCA.GetAlpha();
+      AliHLTTPCGMTrackParam t0 = par;
       for ( int ih = 0; ih < nhits; ih++ ) {
-        int index = hltOut.ClusterId( tCA.FirstClusterRef() + firstHit + ih );
-       index = index&0xffffff;
+        int index = hlt.Merger().OutputClusterIds()[ tCA.FirstClusterRef() + firstHit + ih ];
         tTPC.SetClusterIndex( ih, index );
         AliTPCclusterMI *c = &( fClusters[index] );
-        int iSlice = fClusterSliceRow[index] >> 8;
+        //int iSlice = fClusterSliceRow[index] >> 8;
         int row = fClusterSliceRow[index] & 0xff;
 
         tTPC.SetClusterPointer( row, c );
+       /*
         AliTPCTrackerPoint &point = *( tTPC.GetTrackPoint( row ) );
         {
           //AliHLTTPCCATracker &slice = hlt.SliceTracker( iSlice );
@@ -427,11 +424,12 @@ int AliTPCtrackerCA::Clusters2Tracks( AliESDEvent *event )
           point.SetAngleY( TMath::Abs( t0.GetSinPhi() / t0.GetCosPhi() ) );
           point.SetAngleZ( TMath::Abs( t0.GetDzDs() ) );
         }
+       */
       }
-      tTPC.CookdEdx( 0.02, 0.6 );
+      //tTPC.CookdEdx( 0.02, 0.6 );
 
       CookLabel( &tTPC, 0.1 );
-
+      
       if ( 1 ) { // correction like in off-line --- Adding systematic error
 
         const double *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
@@ -523,10 +521,10 @@ int AliTPCtrackerCA::Clusters2Tracks( AliESDEvent *event )
 }
 
 
-int AliTPCtrackerCA::RefitInward ( AliESDEvent *event )
+int AliTPCtrackerCA::RefitInward ( AliESDEvent * /*event*/ )
 {
   //* forward propagation of ESD tracks
-
+#ifdef XXX
   float xTPC = fkParam->GetInnerRadiusLow();
   float dAlpha = fkParam->GetInnerAngle() / 180.*TMath::Pi();
   float yMax = xTPC * TMath::Tan( dAlpha / 2. );
@@ -551,7 +549,9 @@ int AliTPCtrackerCA::RefitInward ( AliESDEvent *event )
       hits1[i] = i;
       int index = hits[i];
       infos[i].SetISlice( fClusterSliceRow[index] >> 8 );
-      infos[i].SetIRow( fClusterSliceRow[index] & 0xff );
+      int row = fClusterSliceRow[index] & 0xff;
+      int type = ( row < 63 ) ? 0 : ( ( row > 126 ) ? 1 : 2 );
+      infos[i].SetRowType( type );
       infos[i].SetId( index );
       infos[i].SetX( fClusters[index].GetX() );
       infos[i].SetY( fClusters[index].GetY() );
@@ -581,13 +581,14 @@ int AliTPCtrackerCA::RefitInward ( AliESDEvent *event )
       }
     }
   }
+#endif
   return 0;
 }
 
-int AliTPCtrackerCA::PropagateBack( AliESDEvent *event )
+int AliTPCtrackerCA::PropagateBack( AliESDEvent * /*event*/ )
 {
-
   //* backward propagation of ESD tracks
+#ifdef XXX
 
   AliHLTTPCCAStandaloneFramework &hlt = AliHLTTPCCAStandaloneFramework::Instance();
 
@@ -611,7 +612,9 @@ int AliTPCtrackerCA::PropagateBack( AliESDEvent *event )
       hits1[i] = i;
       int index = hits[i];
       infos[i].SetISlice( fClusterSliceRow[index] >> 8 );
-      infos[i].SetIRow( fClusterSliceRow[index] & 0xff );
+      int row = fClusterSliceRow[index] & 0xff;
+      int type = ( row < 63 ) ? 0 : ( ( row > 126 ) ? 1 : 2 );
+      infos[i].SetRowType( type );
       infos[i].SetId( index );
       infos[i].SetX( fClusters[index].GetX() );
       infos[i].SetY( fClusters[index].GetY() );
@@ -627,5 +630,6 @@ int AliTPCtrackerCA::PropagateBack( AliESDEvent *event )
       }
     }
   }
+#endif
   return 0;
 }
index 3eb6a60..fc32694 100644 (file)
@@ -85,6 +85,9 @@ CLASS_HDRS:=          AliHLTTPCTransform.h \
                 tracking-ca/AliHLTTPCCAInputDataCompressorComponent.h \
                 tracking-ca/AliHLTTPCCAGPUTracker.h \
                 tracking-ca/AliHLTTPCCATrackerOutputConverter.h \
+                merger-ca/AliHLTTPCGMTrackParam.h \
+                merger-ca/AliHLTTPCGMSliceTrack.h \
+                merger-ca/AliHLTTPCGMMerger.h \
                comp/AliHLTTPCCompDataCompressorHelper.h \
                comp/AliHLTTPCCompDumpComponent.h \
                comp/AliHLTTPCCompModelAnalysis.h \
@@ -148,6 +151,7 @@ EINCLUDE    := HLT/TPCLib \
                   HLT/TPCLib/tracking \
                   HLT/TPCLib/comp \
                   HLT/TPCLib/tracking-ca \
+                  HLT/TPCLib/merger-ca \
                   HLT/TPCLib/offline \
                   HLT/BASE \
                   HLT/BASE/util \