Improve CA selection algorithm

Improve CA selection algorithm

Author: mpuccio

Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackerTraitsGPU.cu

@@ -240,7 +240,7 @@ GPUd() void computeLayerCells(DeviceStoreNV& devStore, const int layerIndex,
 
                 cellsVector.emplace(currentIndex, currentTracklet.firstClusterIndex,
                                     nextTracklet.firstClusterIndex, nextTracklet.secondClusterIndex, currentTrackletIndex,
-                                    iNextLayerTracklet, normalizedPlaneVector, 1.0f / cellTrajectoryRadius);
+                                    iNextLayerTracklet, averageTanLambda);
                 ++trackletCellsNum;
               }
             }

Detectors/ITSMFT/ITS/tracking/include/ITStracking/Cell.h

@@ -31,16 +31,15 @@ namespace its
 class Cell final
 {
  public:
-  GPUd() Cell(const int, const int, const int, const int, const int, const float3&, const float);
+  GPUd() Cell(const int, const int, const int, const int, const int, const float);
 
   GPUhdni() int getFirstClusterIndex() const { return mFirstClusterIndex; };
   GPUhdni() int getSecondClusterIndex() const { return mSecondClusterIndex; };
   GPUhdni() int getThirdClusterIndex() const { return mThirdClusterIndex; };
-  GPUhd() int getFirstTrackletIndex() const { return mFirstTrackletIndex; };
-  int getSecondTrackletIndex() const { return mSecondTrackletIndex; };
+  GPUhdni() int getFirstTrackletIndex() const { return mFirstTrackletIndex; };
+  GPUhdni() int getSecondTrackletIndex() const { return mSecondTrackletIndex; };
+  GPUhdni() float getTanLambda() const { return mTanLambda; };
   int getLevel() const { return mLevel; };
-  float getCurvature() const { return mCurvature; };
-  const float3& getNormalVectorCoordinates() const { return mNormalVectorCoordinates; };
   void setLevel(const int level) { mLevel = level; };
 
  private:
@@ -49,21 +48,18 @@ class Cell final
   const int mThirdClusterIndex;
   const int mFirstTrackletIndex;
   const int mSecondTrackletIndex;
-  const float3 mNormalVectorCoordinates;
-  const float mCurvature;
+  const float mTanLambda;
   int mLevel;
 };
 
 GPUdi() Cell::Cell(const int firstClusterIndex, const int secondClusterIndex, const int thirdClusterIndex,
-                   const int firstTrackletIndex, const int secondTrackletIndex,
-                   const float3& normalVectorCoordinates, const float curvature)
+                   const int firstTrackletIndex, const int secondTrackletIndex, const float tanL)
   : mFirstClusterIndex{firstClusterIndex},
     mSecondClusterIndex{secondClusterIndex},
     mThirdClusterIndex{thirdClusterIndex},
     mFirstTrackletIndex(firstTrackletIndex),
     mSecondTrackletIndex(secondTrackletIndex),
-    mNormalVectorCoordinates(normalVectorCoordinates),
-    mCurvature{curvature},
+    mTanLambda{tanL},
     mLevel{1}
 {
   // Nothing to do

Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h

@@ -51,17 +51,7 @@ struct TrackingParameters {
   void CopyCuts(TrackingParameters& other, float scale = 1.)
   {
     TrackletMaxDeltaPhi = other.TrackletMaxDeltaPhi * scale;
-    for (unsigned int ii{0}; ii < TrackletMaxDeltaZ.size(); ++ii) {
-      TrackletMaxDeltaZ[ii] = other.TrackletMaxDeltaZ[ii] * scale;
-    }
-    CellMaxDeltaTanLambda = other.CellMaxDeltaTanLambda * scale;
-    for (unsigned int ii{0}; ii < CellMaxDCA.size(); ++ii) {
-      CellMaxDCA[ii] = other.CellMaxDCA[ii] * scale;
-    }
-    for (unsigned int ii{0}; ii < NeighbourMaxDeltaCurvature.size(); ++ii) {
-      NeighbourMaxDeltaCurvature[ii] = other.NeighbourMaxDeltaCurvature[ii] * scale;
-      NeighbourMaxDeltaN[ii] = other.NeighbourMaxDeltaN[ii] * scale;
-    }
+    NSigmaCut = other.NSigmaCut * scale;
   }
 
   int CellMinimumLevel();
@@ -72,23 +62,21 @@ struct TrackingParameters {
   int DeltaROF = 0;
   std::vector<float> LayerZ = {16.333f + 1, 16.333f + 1, 16.333f + 1, 42.140f + 1, 42.140f + 1, 73.745f + 1, 73.745f + 1};
   std::vector<float> LayerRadii = {2.33959f, 3.14076f, 3.91924f, 19.6213f, 24.5597f, 34.388f, 39.3329f};
+  std::vector<float> LayerxX0 = {5.e-3f, 5.e-3f, 5.e-3f, 1.e-2f, 1.e-2f, 1.e-2f, 1.e-2f};
+  std::vector<float> LayerResolution = {5.e-4f, 5.e-4f, 5.e-4f, 5.e-4f, 5.e-4f, 5.e-4f, 5.e-4f};
+  std::vector<float> LayerMisalignment = {0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f};
   int ZBins{256};
   int PhiBins{128};
 
   /// General parameters
   int ClusterSharing = 0;
   int MinTrackLength = 7;
+  float NSigmaCut = 5;
+  float PVres = 1.e-2f; /// FIXME: this has to be taken directly from the reconstructed vertices
   /// Trackleting cuts
   float TrackletMaxDeltaPhi = 0.3f;
-  std::vector<float> TrackletMaxDeltaZ = {0.1f, 0.1f, 0.3f, 0.3f, 0.3f, 0.3f};
   /// Cell finding cuts
-  float CellMaxDeltaTanLambda = 0.025f;
-  std::vector<float> CellMaxDCA = {0.05f, 0.04f, 0.05f, 0.2f, 0.4f};
-  float CellMaxDeltaPhi = 0.14f;
-  std::vector<float> CellMaxDeltaZ = {0.2f, 0.4f, 0.5f, 0.6f, 3.0f};
-  /// Neighbour finding cuts
-  std::vector<float> NeighbourMaxDeltaCurvature = {0.008f, 0.0025f, 0.003f, 0.0035f};
-  std::vector<float> NeighbourMaxDeltaN = {0.002f, 0.0090f, 0.002f, 0.005f};
+  float CellDeltaTanLambdaSigma = 0.007f;
   /// Fitter parameters
   bool UseMatBudLUT = false;
   std::array<float, 2> FitIterationMaxChi2 = {100, 50};

Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h

@@ -78,6 +78,8 @@ class TimeFrame final
 
   float getMinR(int layer) const { return mMinR[layer]; }
   float getMaxR(int layer) const { return mMaxR[layer]; }
+  float getMSangle(int layer) const { return mMSangles[layer]; }
+  float getPositionResolution(int layer) const { return mPositionResolution[layer]; }
 
   gsl::span<Cluster> getClustersOnLayer(int rofId, int layerId);
   gsl::span<const Cluster> getClustersOnLayer(int rofId, int layerId) const;
@@ -142,6 +144,8 @@ class TimeFrame final
   float mBeamPos[2] = {0.f, 0.f};
   std::vector<float> mMinR;
   std::vector<float> mMaxR;
+  std::vector<float> mMSangles;
+  std::vector<float> mPositionResolution;
   std::vector<bool> mMultiplicityCutMask;
   std::vector<int> mROframesPV = {0};
   std::vector<std::vector<int>> mROframesClusters;

Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h

@@ -83,7 +83,7 @@ class Tracker
 
  private:
   track::TrackParCov buildTrackSeed(const Cluster& cluster1, const Cluster& cluster2, const Cluster& cluster3,
-                                    const TrackingFrameInfo& tf3);
+                                    const TrackingFrameInfo& tf3, float resolution);
   template <typename... T>
   void initialiseTimeFrame(T&&... args);
   void computeTracklets();

Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h

@@ -32,6 +32,8 @@
 #include "ITStracking/TimeFrame.h"
 #include "ITStracking/Road.h"
 
+// #define OPTIMISATION_OUTPUT
+
 namespace o2
 {
 namespace gpu

Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h

@@ -40,6 +40,13 @@ struct TrackerParamConfig : public o2::conf::ConfigurableParamHelper<TrackerPara
 
   // Use TGeo for mat. budget
   bool useMatCorrTGeo = false;
+  float misalignment[7] = {-1.f, -1.f, -1.f, -1.f, -1.f, -1.f, -1.f};
+  float nSigmaCut = -1.f;
+  float deltaTanLres = -1.f;
+  float phiCut = -1.f;
+  float pvRes = -1.f;
+  int LUTbinsPhi = -1;
+  int LUTbinsZ = -1;
 
   O2ParamDef(TrackerParamConfig, "ITSCATrackerParam");
 };

Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracklet.h

@@ -29,6 +29,7 @@ namespace its
 struct Tracklet final {
   Tracklet();
   GPUdi() Tracklet(const int, const int, const Cluster&, const Cluster&, int rof0, int rof1);
+  GPUdi() Tracklet(const int, const int, float tanL, float phi, int rof0, int rof1);
 #ifdef _ALLOW_DEBUG_TREES_ITS_
   unsigned char isEmpty() const;
   void dump();
@@ -60,6 +61,16 @@ GPUdi() Tracklet::Tracklet(const int firstClusterOrderingIndex, const int second
   // Nothing to do
 }
 
+GPUdi() Tracklet::Tracklet(const int idx0, const int idx1, float tanL, float phi, int rof0, int rof1)
+  : firstClusterIndex{idx0},
+    secondClusterIndex{idx1},
+    tanLambda{tanL},
+    phi{phi},
+    rof{static_cast<unsigned short>(rof0), static_cast<unsigned short>(rof1)}
+{
+  // Nothing to do
+}
+
 #ifdef _ALLOW_DEBUG_TREES_ITS_
 inline unsigned char Tracklet::isEmpty() const
 {

Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx

@@ -20,6 +20,7 @@
 #include "DataFormatsITSMFT/TopologyDictionary.h"
 #include "ITSBase/GeometryTGeo.h"
 #include "ITSMFTBase/SegmentationAlpide.h"
+#include "ITStracking/TrackingConfigParam.h"
 
 #include <iostream>
 
@@ -31,6 +32,13 @@ struct ClusterHelper {
   int bin;
   int ind;
 };
+
+float MSangle(float mass, float p, float xX0)
+{
+  float beta = p / std::hypot(mass, p);
+  return 0.0136f * std::sqrt(xX0) * (1.f + 0.038f * std::log(xX0)) / (beta * p);
+}
+
 } // namespace
 
 namespace o2
@@ -180,6 +188,8 @@ void TimeFrame::initialise(const int iteration, const MemoryParameters& memParam
     mTrackletsLookupTable.resize(trkParam.CellsPerRoad());
     mIndexTables.clear();
     mIndexTableUtils.setTrackingParameters(trkParam);
+    mMSangles.resize(trkParam.NLayers);
+    mPositionResolution.resize(trkParam.NLayers);
 
     for (unsigned int iLayer{0}; iLayer < mClusters.size(); ++iLayer) {
       if (mClusters[iLayer].size()) {
@@ -189,6 +199,8 @@ void TimeFrame::initialise(const int iteration, const MemoryParameters& memParam
       mClusters[iLayer].resize(mUnsortedClusters[iLayer].size());
       mUsedClusters[iLayer].clear();
       mUsedClusters[iLayer].resize(mUnsortedClusters[iLayer].size(), false);
+      mMSangles[iLayer] = MSangle(0.14f, 0.3f, trkParam.LayerxX0[iLayer]);
+      mPositionResolution[iLayer] = std::hypot(trkParam.LayerMisalignment[iLayer], trkParam.LayerResolution[iLayer]);
     }
 
     mIndexTables.resize(mNrof);

Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx

@@ -97,6 +97,9 @@ void Tracker::computeCells()
 
 void Tracker::findCellsNeighbours(int& iteration)
 {
+#ifdef OPTIMISATION_OUTPUT
+  std::ofstream off(fmt::format("cellneighs{}.txt", iteration));
+#endif
   for (int iLayer{0}; iLayer < mTrkParams[iteration].CellsPerRoad() - 1; ++iLayer) {
 
     if (mTimeFrame->getCells()[iLayer + 1].empty() ||
@@ -114,35 +117,26 @@ void Tracker::findCellsNeighbours(int& iteration)
       const int nextLayerTrackletIndex{currentCell.getSecondTrackletIndex()};
       const int nextLayerFirstCellIndex{mTimeFrame->getCellsLookupTable()[iLayer][nextLayerTrackletIndex]};
       const int nextLayerLastCellIndex{mTimeFrame->getCellsLookupTable()[iLayer][nextLayerTrackletIndex + 1]};
-      for (int iNextLayerCell{nextLayerFirstCellIndex}; iNextLayerCell < nextLayerLastCellIndex; ++iNextLayerCell) {
+      for (int iNextCell{nextLayerFirstCellIndex}; iNextCell < nextLayerLastCellIndex; ++iNextCell) {
 
-        Cell& nextCell{mTimeFrame->getCells()[iLayer + 1][iNextLayerCell]};
+        Cell& nextCell{mTimeFrame->getCells()[iLayer + 1][iNextCell]};
         if (nextCell.getFirstTrackletIndex() != nextLayerTrackletIndex) {
           break;
         }
 
-        const float3 currentCellNormalVector{currentCell.getNormalVectorCoordinates()};
-        const float3 nextCellNormalVector{nextCell.getNormalVectorCoordinates()};
-        const float3 normalVectorsDeltaVector{currentCellNormalVector.x - nextCellNormalVector.x,
-                                              currentCellNormalVector.y - nextCellNormalVector.y,
-                                              currentCellNormalVector.z - nextCellNormalVector.z};
-
-        const float deltaNormalVectorsModulus{(normalVectorsDeltaVector.x * normalVectorsDeltaVector.x) +
-                                              (normalVectorsDeltaVector.y * normalVectorsDeltaVector.y) +
-                                              (normalVectorsDeltaVector.z * normalVectorsDeltaVector.z)};
-        const float deltaCurvature{std::abs(currentCell.getCurvature() - nextCell.getCurvature())};
-
-        if (deltaNormalVectorsModulus < mTrkParams[iteration].NeighbourMaxDeltaN[iLayer] &&
-            deltaCurvature < mTrkParams[iteration].NeighbourMaxDeltaCurvature[iLayer]) {
-
-          mTimeFrame->getCellsNeighbours()[iLayer][iNextLayerCell].push_back(iCell);
+#ifdef OPTIMISATION_OUTPUT
+        bool good{mTimeFrame->getCellsLabel(iLayer)[iCell] == mTimeFrame->getCellsLabel(iLayer + 1)[iNextCell]};
+        float signedDelta{currentCell.getTanLambda() - nextCell.getTanLambda()};
+        off << fmt::format("{}\t{:d}\t{}\t{}", iLayer, good, signedDelta, signedDelta / mTrkParams[iteration].CellDeltaTanLambdaSigma) << std::endl;
+#endif
+        mTimeFrame->getCellsNeighbours()[iLayer][iNextCell].push_back(iCell);
 
-          const int currentCellLevel{currentCell.getLevel()};
+        const int currentCellLevel{currentCell.getLevel()};
 
-          if (currentCellLevel >= nextCell.getLevel()) {
-            nextCell.setLevel(currentCellLevel + 1);
-          }
+        if (currentCellLevel >= nextCell.getLevel()) {
+          nextCell.setLevel(currentCellLevel + 1);
         }
+        // }
       }
     }
   }
@@ -283,7 +277,7 @@ void Tracker::findTracks()
     const auto& cluster3_tf = mTimeFrame->getTrackingFrameInfoOnLayer(lastCellLevel).at(clusters[lastCellLevel]);
 
     /// FIXME!
-    TrackITSExt temporaryTrack{buildTrackSeed(cluster1_glo, cluster2_glo, cluster3_glo, cluster3_tf)};
+    TrackITSExt temporaryTrack{buildTrackSeed(cluster1_glo, cluster2_glo, cluster3_glo, cluster3_tf, mTimeFrame->getPositionResolution(lastCellLevel))};
     for (size_t iC = 0; iC < clusters.size(); ++iC) {
       temporaryTrack.setExternalClusterIndex(iC, clusters[iC], clusters[iC] != constants::its::UnusedIndex);
     }
@@ -376,12 +370,15 @@ bool Tracker::fitTrack(TrackITSExt& track, int start, int end, int step, const f
       return false;
     }
 
-    auto predChi2{track.getPredictedChi2(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)};
+    GPUArray<float, 3> cov{trackingHit.covarianceTrackingFrame};
+    cov[0] = std::hypot(cov[0], mTrkParams[0].LayerMisalignment[iLayer]);
+    cov[2] = std::hypot(cov[2], mTrkParams[0].LayerMisalignment[iLayer]);
+    auto predChi2{track.getPredictedChi2(trackingHit.positionTrackingFrame, cov)};
     if (nCl >= 3 && predChi2 > chi2cut * (nCl * 2 - 5)) {
       return false;
     }
     track.setChi2(track.getChi2() + predChi2);
-    if (!track.o2::track::TrackParCov::update(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)) {
+    if (!track.o2::track::TrackParCov::update(trackingHit.positionTrackingFrame, cov)) {
       return false;
     }
     nCl++;
@@ -608,7 +605,7 @@ void Tracker::rectifyClusterIndices()
 /// whereas the others are referred to the global frame. This function is almost a clone of CookSeed, adapted to return
 /// a TrackParCov
 track::TrackParCov Tracker::buildTrackSeed(const Cluster& cluster1, const Cluster& cluster2,
-                                           const Cluster& cluster3, const TrackingFrameInfo& tf3)
+                                           const Cluster& cluster3, const TrackingFrameInfo& tf3, float resolution)
 {
   const float ca = std::cos(tf3.alphaTrackingFrame), sa = std::sin(tf3.alphaTrackingFrame);
   const float x1 = cluster1.xCoordinate * ca + cluster1.yCoordinate * sa;
@@ -628,10 +625,10 @@ track::TrackParCov Tracker::buildTrackSeed(const Cluster& cluster1, const Cluste
 
   const float fy = 1. / (cluster2.radius - cluster3.radius);
   const float& tz = fy;
-  const float cy = (math_utils::computeCurvature(x1, y1, x2, y2 + constants::its::Resolution, x3, y3) - crv) /
-                   (constants::its::Resolution * getBz() * o2::constants::math::B2C) *
+  const float cy = (math_utils::computeCurvature(x1, y1, x2, y2 + resolution, x3, y3) - crv) /
+                   (resolution * getBz() * o2::constants::math::B2C) *
                    20.f; // FIXME: MS contribution to the cov[14] (*20 added)
-  constexpr float s2 = constants::its::Resolution * constants::its::Resolution;
+  const float s2 = resolution;
 
   return track::TrackParCov(tf3.xTrackingFrame, tf3.alphaTrackingFrame,
                             {y3, z3, crv * (x3 - x0), 0.5f * (tgl12 + tgl23),
@@ -647,6 +644,19 @@ void Tracker::getGlobalConfiguration()
   if (tc.useMatCorrTGeo) {
     setCorrType(o2::base::PropagatorImpl<float>::MatCorrType::USEMatCorrTGeo);
   }
+  for (auto& params : mTrkParams) {
+    if (params.NLayers == 7) {
+      for (int i{0}; i < 7; ++i) {
+        params.LayerMisalignment[i] = tc.misalignment[i] > 0 ? tc.misalignment[i] : params.LayerMisalignment[i];
+      }
+    }
+    params.PhiBins = tc.LUTbinsPhi > 0 ? tc.LUTbinsPhi : params.PhiBins;
+    params.ZBins = tc.LUTbinsZ > 0 ? tc.LUTbinsZ : params.ZBins;
+    params.PVres = tc.pvRes > 0 ? tc.pvRes : params.PVres;
+    params.NSigmaCut *= tc.nSigmaCut > 0 ? tc.nSigmaCut : 1.f;
+    params.CellDeltaTanLambdaSigma *= tc.deltaTanLres > 0 ? tc.deltaTanLres : 1.f;
+    params.TrackletMaxDeltaPhi *= tc.phiCut > 0 ? tc.phiCut : 1.f;
+  }
 }
 
 void Tracker::adoptTimeFrame(TimeFrame& tf)