ITS: cleanup lines and tracklet/clusterlines changes (#14511)

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

FullCi is green, merging.

Author: f3sch

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

@@ -15,25 +15,25 @@
 #include <array>
 #include <vector>
 #include "ITStracking/Cluster.h"
-#include "ITStracking/Definitions.h"
+#include "ITStracking/Constants.h"
 #include "ITStracking/Tracklet.h"
+#include "GPUCommonRtypes.h"
 #include "GPUCommonMath.h"
 
 namespace o2::its
 {
 struct Line final {
-  GPUhd() Line();
+  GPUhdDefault() Line() = default;
   GPUhd() Line(const Line&);
   Line(std::array<float, 3> firstPoint, std::array<float, 3> secondPoint);
-  GPUhd() Line(const float firstPoint[3], const float secondPoint[3]);
   GPUhd() 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 = 1e-14);
-  static bool areParallel(const Line&, const Line&, const float precision = 1e-14);
+  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); }
   GPUhdi() auto getDeltaROF() const { return rof[1] - rof[0]; }
@@ -42,56 +42,35 @@ struct Line final {
   bool operator!=(const Line&) const;
   short getMinROF() const { return rof[0] < rof[1] ? rof[0] : rof[1]; }
 
-  float originPoint[3], cosinesDirector[3];
-  float weightMatrix[6] = {1., 0., 0., 1., 0., 1.};
+  float originPoint[3] = {0};
+  float cosinesDirector[3] = {0};
+  // float weightMatrix[6] = {1., 0., 0., 1., 0., 1.};
   // weightMatrix is a symmetric matrix internally stored as
   //    0 --> row = 0, col = 0
   //    1 --> 0,1
   //    2 --> 0,2
   //    3 --> 1,1
   //    4 --> 1,2
   //    5 --> 2,2
-  short rof[2];
-};
+  short rof[2] = {-1, -1};
 
-GPUhdi() Line::Line() : weightMatrix{1., 0., 0., 1., 0., 1.}
-{
-  rof[0] = -1;
-  rof[1] = -1;
-}
+  ClassDefNV(Line, 1);
+};
 
 GPUhdi() Line::Line(const Line& other)
 {
   for (int i{0}; i < 3; ++i) {
     originPoint[i] = other.originPoint[i];
     cosinesDirector[i] = other.cosinesDirector[i];
   }
-  for (int i{0}; i < 6; ++i) {
-    weightMatrix[i] = other.weightMatrix[i];
-  }
+  // for (int i{0}; i < 6; ++i) {
+  //   weightMatrix[i] = other.weightMatrix[i];
+  // }
   for (int i{0}; i < 2; ++i) {
     rof[i] = other.rof[i];
   }
 }
 
-GPUhdi() Line::Line(const float firstPoint[3], const float secondPoint[3])
-{
-  for (int i{0}; i < 3; ++i) {
-    originPoint[i] = firstPoint[i];
-    cosinesDirector[i] = secondPoint[i] - firstPoint[i];
-  }
-
-  float inverseNorm{1.f / o2::gpu::CAMath::Sqrt(cosinesDirector[0] * cosinesDirector[0] + cosinesDirector[1] * cosinesDirector[1] +
-                                                cosinesDirector[2] * cosinesDirector[2])};
-
-  for (int index{0}; index < 3; ++index) {
-    cosinesDirector[index] *= inverseNorm;
-  }
-
-  rof[0] = -1;
-  rof[1] = -1;
-}
-
 GPUhdi() Line::Line(const Tracklet& tracklet, const Cluster* innerClusters, const Cluster* outerClusters)
 {
   originPoint[0] = innerClusters[tracklet.firstClusterIndex].xCoordinate;
@@ -102,12 +81,10 @@ GPUhdi() Line::Line(const Tracklet& tracklet, const Cluster* innerClusters, cons
   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::Sqrt(cosinesDirector[0] * cosinesDirector[0] + cosinesDirector[1] * cosinesDirector[1] +
-                                                cosinesDirector[2] * cosinesDirector[2])};
-
-  for (int index{0}; index < 3; ++index) {
-    cosinesDirector[index] *= inverseNorm;
-  }
+  float inverseNorm{1.f / o2::gpu::CAMath::Hypot(cosinesDirector[0], cosinesDirector[1], cosinesDirector[2])};
+  cosinesDirector[0] *= inverseNorm;
+  cosinesDirector[1] *= inverseNorm;
+  cosinesDirector[2] *= inverseNorm;
 
   rof[0] = tracklet.rof[0];
   rof[1] = tracklet.rof[1];
@@ -130,47 +107,38 @@ inline float Line::getDistanceFromPoint(const Line& line, const std::array<float
 
 GPUhdi() float Line::getDistanceFromPoint(const Line& line, const float point[3])
 {
-  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);
+  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)
 {
-  float normalVector[3];
-  normalVector[0] = firstLine.cosinesDirector[1] * secondLine.cosinesDirector[2] -
-                    firstLine.cosinesDirector[2] * secondLine.cosinesDirector[1];
-  normalVector[1] = -firstLine.cosinesDirector[0] * secondLine.cosinesDirector[2] +
-                    firstLine.cosinesDirector[2] * secondLine.cosinesDirector[0];
-  normalVector[2] = firstLine.cosinesDirector[0] * secondLine.cosinesDirector[1] -
-                    firstLine.cosinesDirector[1] * secondLine.cosinesDirector[0];
-
-  float norm{0.f}, distance{0.f};
-  for (int i{0}; i < 3; ++i) {
-    norm += normalVector[i] * normalVector[i];
-    distance += (secondLine.originPoint[i] - firstLine.originPoint[i]) * normalVector[i];
-  }
-  if (norm > precision) {
-    return o2::gpu::CAMath::Abs(distance / o2::gpu::CAMath::Sqrt(norm));
-  } else {
-#if defined(__CUDACC__) || defined(__HIPCC__)
-    float stdOriginPoint[3];
-    for (int i{0}; i < 3; ++i) {
-      stdOriginPoint[i] = secondLine.originPoint[1];
-    }
-#else
-    std::array<float, 3> stdOriginPoint = {};
-    std::copy_n(secondLine.originPoint, 3, stdOriginPoint.begin());
-#endif
-    return getDistanceFromPoint(firstLine, stdOriginPoint);
+  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])
@@ -199,11 +167,7 @@ inline bool Line::operator==(const Line& rhs) const
 
 inline bool Line::operator!=(const Line& rhs) const
 {
-  bool val;
-  for (int i{0}; i < 3; ++i) {
-    val &= this->originPoint[i] != rhs.originPoint[i];
-  }
-  return val;
+  return !(*this == rhs);
 }
 
 GPUhdi() void Line::print() const

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

@@ -20,53 +20,48 @@
 #include "GPUCommonRtypes.h"
 #include "GPUCommonMath.h"
 #include "GPUCommonDef.h"
+#include "GPUCommonLogger.h"
 
 #ifndef GPUCA_GPUCODE_DEVICE
+#ifndef GPU_NO_FMT
 #include <string>
+#include <fmt/format.h>
+#endif
 #endif
 
 namespace o2::its
 {
 
 struct Tracklet final {
-  GPUhdi() Tracklet();
+  GPUhdDefault() Tracklet() = default;
   GPUhdi() Tracklet(const int, const int, const Cluster&, const Cluster&, short rof0, short rof1);
   GPUhdi() Tracklet(const int, const int, float tanL, float phi, short rof0, short rof1);
-  GPUhdi() bool operator==(const Tracklet&) const;
-  GPUhdi() bool operator!=(const Tracklet&) const;
+  GPUhdDefault() bool operator==(const Tracklet&) const = default;
   GPUhdi() unsigned char isEmpty() const
   {
     return firstClusterIndex < 0 || secondClusterIndex < 0;
   }
   GPUhdi() auto getDeltaRof() const { return rof[1] - rof[0]; }
-  GPUhdi() void dump();
   GPUhdi() void dump() const;
-  GPUhdi() void dump(const int, const int);
   GPUhdi() void dump(const int, const int) const;
   GPUhdi() unsigned char operator<(const Tracklet&) const;
-#ifndef GPUCA_GPUCODE_DEVICE
+#if !defined(GPUCA_NO_FMT) && !defined(GPUCA_GPUCODE_DEVICE)
   std::string asString() const
   {
-    return "fClIdx: " + std::to_string(firstClusterIndex) + " sClIdx: " + std::to_string(secondClusterIndex) +
-           " rof1: " + std::to_string(rof[0]) + " rof2: " + std::to_string(rof[1]) + " delta: " + std::to_string(getDeltaRof());
+    return fmt::format("fClIdx:{} fROF:{} sClIdx:{} sROF:{} (DROF:{})", firstClusterIndex, rof[0], secondClusterIndex, rof[1], getDeltaRof());
   }
+  void print() const { LOG(info) << asString(); }
 #endif
 
-  int firstClusterIndex;
-  int secondClusterIndex;
-  float tanLambda;
-  float phi;
-  short rof[2];
+  int firstClusterIndex{-1};
+  int secondClusterIndex{-1};
+  float tanLambda{-999};
+  float phi{-999};
+  short rof[2] = {-1, -1};
 
   ClassDefNV(Tracklet, 1);
 };
 
-GPUhdi() Tracklet::Tracklet() : firstClusterIndex{-1}, secondClusterIndex{-1}, tanLambda{0.0f}, phi{0.0f}
-{
-  rof[0] = -1;
-  rof[1] = -1;
-}
-
 GPUhdi() Tracklet::Tracklet(const int firstClusterOrderingIndex, const int secondClusterOrderingIndex,
                             const Cluster& firstCluster, const Cluster& secondCluster, short rof0 = -1, short rof1 = -1)
   : firstClusterIndex{firstClusterOrderingIndex},
@@ -90,24 +85,6 @@ GPUhdi() Tracklet::Tracklet(const int idx0, const int idx1, float tanL, float ph
   // Nothing to do
 }
 
-GPUhdi() bool Tracklet::operator==(const Tracklet& rhs) const
-{
-  return this->firstClusterIndex == rhs.firstClusterIndex &&
-         this->secondClusterIndex == rhs.secondClusterIndex &&
-         this->tanLambda == rhs.tanLambda &&
-         this->phi == rhs.phi &&
-         this->rof[0] == rhs.rof[0] &&
-         this->rof[1] == rhs.rof[1];
-}
-
-GPUhdi() bool Tracklet::operator!=(const Tracklet& rhs) const
-{
-  return this->firstClusterIndex != rhs.firstClusterIndex ||
-         this->secondClusterIndex != rhs.secondClusterIndex ||
-         this->tanLambda != rhs.tanLambda ||
-         this->phi != rhs.phi;
-}
-
 GPUhdi() unsigned char Tracklet::operator<(const Tracklet& t) const
 {
   if (isEmpty()) {
@@ -116,21 +93,11 @@ GPUhdi() unsigned char Tracklet::operator<(const Tracklet& t) const
   return true;
 }
 
-GPUhdi() void Tracklet::dump(const int offsetFirst, const int offsetSecond)
-{
-  printf("fClIdx: %d sClIdx: %d  rof1: %hu rof2: %hu\n", firstClusterIndex + offsetFirst, secondClusterIndex + offsetSecond, rof[0], rof[1]);
-}
-
 GPUhdi() void Tracklet::dump(const int offsetFirst, const int offsetSecond) const
 {
   printf("fClIdx: %d sClIdx: %d  rof1: %hu rof2: %hu\n", firstClusterIndex + offsetFirst, secondClusterIndex + offsetSecond, rof[0], rof[1]);
 }
 
-GPUhdi() void Tracklet::dump()
-{
-  printf("fClIdx: %d sClIdx: %d  rof1: %hu rof2: %hu\n", firstClusterIndex, secondClusterIndex, rof[0], rof[1]);
-}
-
 GPUhdi() void Tracklet::dump() const
 {
   printf("fClIdx: %d sClIdx: %d  rof1: %hu rof2: %hu\n", firstClusterIndex, secondClusterIndex, rof[0], rof[1]);