ITS: fix degenerate LSE for matrix solving

Comparing the output of dev and this PR, I saw plently of cases where
the system of equation was fully degenerate and produced to different
floating instructions and compiler optimizations slightly different
results. The solution is to discard the vertex cand. if the LSE becomes
degenerate as not to produce non-sense solutions.

Signed-off-by: Felix Schlepper <felix.schlepper@cern.ch>

Author: f3sch

Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/ClusterLinesGPU.h

@@ -14,150 +14,56 @@
 
 #include <array>
 #include <vector>
+#include <Math/SMatrix.h>
+#include <Math/SVector.h>
+#include "Framework/Logger.h"
 #include "ITStracking/Cluster.h"
 #include "ITStracking/Constants.h"
-#include "ITStracking/Definitions.h"
 #include "ITStracking/Tracklet.h"
-#include "GPUCommonRtypes.h"
-#include "GPUCommonMath.h"
-#include "GPUCommonDef.h"
-#include "GPUCommonLogger.h"
 
 namespace o2::its
 {
+
 struct Line final {
-  GPUhdDefault() Line() = default;
-  GPUhdDefault() Line(const Line&) = default;
-  Line(std::array<float, 3> firstPoint, std::array<float, 3> secondPoint);
-  GPUhd() Line(const Tracklet&, const Cluster*, const Cluster*);
+  using SVector3f = ROOT::Math::SVector<float, 3>;
+  using SMatrix3f = ROOT::Math::SMatrix<float, 3, 3, ROOT::Math::MatRepSym<float, 3>>;
+
+  Line() = default;
+  Line(const Tracklet&, const Cluster*, const Cluster*);
 
   static float getDistanceFromPoint(const Line& line, const std::array<float, 3>& point);
-  GPUhd() static float getDistanceFromPoint(const Line& line, const float point[3]);
-  static std::array<float, 6> getDCAComponents(const Line& line, const std::array<float, 3> point);
-  GPUhd() static void getDCAComponents(const Line& line, const float point[3], float destArray[6]);
-  GPUhd() static float getDCA(const Line&, const Line&, const float precision = constants::Tolerance);
-  static bool areParallel(const Line&, const Line&, const float precision = constants::Tolerance);
-  GPUhd() unsigned char isEmpty() const { return (originPoint[0] == 0.f && originPoint[1] == 0.f && originPoint[2] == 0.f) &&
-                                                 (cosinesDirector[0] == 0.f && cosinesDirector[1] == 0.f && cosinesDirector[2] == 0.f); }
-  bool operator==(const Line&) const;
-  bool operator!=(const Line&) const;
-  GPUh() void print() const
+  static SMatrix3f getDCAComponents(const Line& line, const std::array<float, 3>& point);
+  static float getDCA(const Line&, const Line&, const float precision = constants::Tolerance);
+  bool isEmpty() const { return ROOT::Math::Dot(originPoint, originPoint) == 0.f &&
+                                ROOT::Math::Dot(cosinesDirector, cosinesDirector) == 0.f; }
+  bool operator==(const Line&) const = default;
+  void print() const
   {
-    LOGP(info, "TRKLT: x={} y={} z={} dx={} dy={} dz={} ts:{}+/-{}", originPoint[0], originPoint[1], originPoint[2], cosinesDirector[0], cosinesDirector[1], cosinesDirector[2], mTime.getTimeStamp(), mTime.getTimeStampError());
+    LOGP(info, "TRKLT: x={} y={} z={} dx={} dy={} dz={} ts:{}+/-{}", originPoint(0), originPoint(1), originPoint(2), cosinesDirector(0), cosinesDirector(1), cosinesDirector(2), mTime.getTimeStamp(), mTime.getTimeStampError());
   }
 
-  float originPoint[3] = {0, 0, 0};
-  float cosinesDirector[3] = {0, 0, 0};
+  SVector3f originPoint;
+  SVector3f cosinesDirector;
   TimeEstBC mTime;
 
   ClassDefNV(Line, 1);
 };
 
-GPUhdi() Line::Line(const Tracklet& tracklet, const Cluster* innerClusters, const Cluster* outerClusters) : mTime(tracklet.mTime)
-{
-  originPoint[0] = innerClusters[tracklet.firstClusterIndex].xCoordinate;
-  originPoint[1] = innerClusters[tracklet.firstClusterIndex].yCoordinate;
-  originPoint[2] = innerClusters[tracklet.firstClusterIndex].zCoordinate;
-
-  cosinesDirector[0] = outerClusters[tracklet.secondClusterIndex].xCoordinate - innerClusters[tracklet.firstClusterIndex].xCoordinate;
-  cosinesDirector[1] = outerClusters[tracklet.secondClusterIndex].yCoordinate - innerClusters[tracklet.firstClusterIndex].yCoordinate;
-  cosinesDirector[2] = outerClusters[tracklet.secondClusterIndex].zCoordinate - innerClusters[tracklet.firstClusterIndex].zCoordinate;
-
-  float inverseNorm{1.f / o2::gpu::CAMath::Hypot(cosinesDirector[0], cosinesDirector[1], cosinesDirector[2])};
-  cosinesDirector[0] *= inverseNorm;
-  cosinesDirector[1] *= inverseNorm;
-  cosinesDirector[2] *= inverseNorm;
-}
-
-// static functions:
-inline float Line::getDistanceFromPoint(const Line& line, const std::array<float, 3>& point)
-{
-  float DCASquared{0};
-  float cdelta{0};
-  for (int i{0}; i < 3; ++i) {
-    cdelta -= line.cosinesDirector[i] * (line.originPoint[i] - point[i]);
-  }
-  for (int i{0}; i < 3; ++i) {
-    DCASquared += (line.originPoint[i] - point[i] + line.cosinesDirector[i] * cdelta) *
-                  (line.originPoint[i] - point[i] + line.cosinesDirector[i] * cdelta);
-  }
-  return o2::gpu::CAMath::Sqrt(DCASquared);
-}
-
-GPUhdi() float Line::getDistanceFromPoint(const Line& line, const float point[3])
-{
-  const float dx = point[0] - line.originPoint[0];
-  const float dy = point[1] - line.originPoint[1];
-  const float dz = point[2] - line.originPoint[2];
-  const float d = (dx * line.cosinesDirector[0]) + (dy * line.cosinesDirector[1]) + (dz * line.cosinesDirector[2]);
-
-  const float vx = dx - (d * line.cosinesDirector[0]);
-  const float vy = dy - (d * line.cosinesDirector[1]);
-  const float vz = dz - (d * line.cosinesDirector[2]);
-
-  return o2::gpu::CAMath::Hypot(vx, vy, vz);
-}
-
-GPUhdi() float Line::getDCA(const Line& firstLine, const Line& secondLine, const float precision)
-{
-  const float nx = (firstLine.cosinesDirector[1] * secondLine.cosinesDirector[2]) -
-                   (firstLine.cosinesDirector[2] * secondLine.cosinesDirector[1]);
-  const float ny = -(firstLine.cosinesDirector[0] * secondLine.cosinesDirector[2]) +
-                   (firstLine.cosinesDirector[2] * secondLine.cosinesDirector[0]);
-  const float nz = (firstLine.cosinesDirector[0] * secondLine.cosinesDirector[1]) -
-                   (firstLine.cosinesDirector[1] * secondLine.cosinesDirector[0]);
-  const float norm2 = (nx * nx) + (ny * ny) + (nz * nz);
-
-  if (norm2 <= precision * precision) {
-    return getDistanceFromPoint(firstLine, secondLine.originPoint);
-  }
-
-  const float dx = secondLine.originPoint[0] - firstLine.originPoint[0];
-  const float dy = secondLine.originPoint[1] - firstLine.originPoint[1];
-  const float dz = secondLine.originPoint[2] - firstLine.originPoint[2];
-  const float triple = (dx * nx) + (dy * ny) + (dz * nz);
-
-  return o2::gpu::CAMath::Abs(triple) / o2::gpu::CAMath::Sqrt(norm2);
-}
-
-GPUhdi() void Line::getDCAComponents(const Line& line, const float point[3], float destArray[6])
-{
-  float cdelta{0.};
-  for (int i{0}; i < 3; ++i) {
-    cdelta -= line.cosinesDirector[i] * (line.originPoint[i] - point[i]);
-  }
-
-  destArray[0] = line.originPoint[0] - point[0] + line.cosinesDirector[0] * cdelta;
-  destArray[3] = line.originPoint[1] - point[1] + line.cosinesDirector[1] * cdelta;
-  destArray[5] = line.originPoint[2] - point[2] + line.cosinesDirector[2] * cdelta;
-  destArray[1] = o2::gpu::CAMath::Hypot(destArray[0], destArray[3]);
-  destArray[2] = o2::gpu::CAMath::Hypot(destArray[0], destArray[5]);
-  destArray[4] = o2::gpu::CAMath::Hypot(destArray[3], destArray[5]);
-}
-
-inline bool Line::operator==(const Line& rhs) const
-{
-  bool val{false};
-  for (int i{0}; i < 3; ++i) {
-    val &= this->originPoint[i] == rhs.originPoint[i];
-  }
-  return val;
-}
-
-inline bool Line::operator!=(const Line& rhs) const
-{
-  return !(*this == rhs);
-}
-
 class ClusterLines final
 {
+  using SMatrix3 = ROOT::Math::SMatrix<double, 3, 3, ROOT::Math::MatRepSym<double, 3>>;
+  using SMatrix3f = ROOT::Math::SMatrix<float, 3, 3, ROOT::Math::MatRepSym<float, 3>>;
+  using SVector3 = ROOT::Math::SVector<double, 3>;
+
  public:
   ClusterLines() = default;
   ClusterLines(const int firstLabel, const Line& firstLine, const int secondLabel, const Line& secondLine);
   void add(const int lineLabel, const Line& line);
   void computeClusterCentroid();
+  void accumulate(const Line& line);
+  bool isValid() const noexcept { return mIsValid; }
   auto const& getVertex() const { return mVertex; }
-  std::array<float, 6> getRMS2() const { return mRMS2; }
+  const float* getRMS2() const { return mRMS2.Array(); }
   float getAvgDistance2() const { return mAvgDistance2; }
   auto getSize() const noexcept { return mLabels.size(); }
   auto& getLabels() noexcept { return mLabels; }
@@ -167,13 +73,14 @@ class ClusterLines final
   float getR() const noexcept { return std::sqrt(getR2()); }
 
  protected:
-  std::array<float, 6> mAMatrix = {}; // AX=B
-  std::array<float, 3> mBMatrix = {}; // AX=B
-  std::array<float, 3> mVertex = {};  // cluster centroid position
-  std::array<float, 6> mRMS2 = {};    // symmetric matrix: diagonal is RMS2
-  float mAvgDistance2 = 0.f;          // substitute for chi2
-  TimeEstBC mTime;                    // time stamp
-  std::vector<int> mLabels;           // contributing labels
+  SMatrix3 mAMatrix;                 // AX=B, symmetric normal matrix
+  SVector3 mBMatrix;                 // AX=B, right-hand side
+  std::array<float, 3> mVertex = {}; // cluster centroid position
+  SMatrix3f mRMS2;                   // symmetric matrix: diagonal is RMS2
+  float mAvgDistance2 = 0.f;         // substitute for chi2
+  bool mIsValid = false;             // true if linear system was solved successfully
+  TimeEstBC mTime;                   // time stamp
+  std::vector<int> mLabels;          // contributing labels
 };
 
 } // namespace o2::its

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

@@ -87,7 +87,6 @@ struct VertexingParameters {
   std::string asString() const;
 
   int nIterations = 1; // Number of vertexing passes to perform
-  bool allowSingleContribClusters = false;
   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};
   int ZBins{1};

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

@@ -22,10 +22,9 @@ namespace o2::its
 struct VertexerParamConfig : public o2::conf::ConfigurableParamHelper<VertexerParamConfig> {
   bool saveTimeBenchmarks = false; // dump metrics on file
 
-  int nIterations = 1;                     // Number of vertexing passes to perform.
-  int vertPerRofThreshold = 0;             // Maximum number of vertices per ROF to trigger second a iteration.
-  bool allowSingleContribClusters = false; // attempt to find vertices in case of a single tracklet found.
-  int deltaRof = 0;                        // Number of ROFs to be considered for the vertexing.
+  int nIterations = 1;         // Number of vertexing passes to perform.
+  int vertPerRofThreshold = 0; // Maximum number of vertices per ROF to trigger second a iteration.
+  int deltaRof = 0;            // Number of ROFs to be considered for the vertexing.
 
   // geometrical cuts for tracklet selection
   float zCut = 0.002f;