Runtime configuration for the ITS CA Tracker

DataFormats/Detectors/ITSMFT/ITS/include/DataFormatsITS/TrackITS.h

@@ -55,7 +55,7 @@ class TrackITS : public o2::track::TrackParCov
   // Bool_t getPhiZat(Float_t r,Float_t &phi,Float_t &z) const;
 
   void setClusterIndex(Int_t layer, Int_t index);
-  void setExternalClusterIndex(Int_t layer, Int_t idx);
+  void setExternalClusterIndex(Int_t layer, Int_t idx, bool newCluster = false);
   void resetClusters();
 
   void setChi2(float chi2) { mChi2 = chi2; }

DataFormats/Detectors/ITSMFT/ITS/src/TrackITS.cxx

@@ -86,11 +86,13 @@ void TrackITS::setClusterIndex(Int_t l, Int_t i)
   mIndex[mNClusters++] = idx;
 }
 
-void TrackITS::setExternalClusterIndex(Int_t layer, Int_t idx)
+void TrackITS::setExternalClusterIndex(Int_t layer, Int_t idx, bool newCluster)
 {
   //--------------------------------------------------------------------
   // Set the cluster index within an external cluster array
   //--------------------------------------------------------------------
+  if (newCluster)
+    mNClusters++;
   mIndex[layer] = idx;
 }
 

DataFormats/Reconstruction/src/Track.cxx

@@ -141,7 +141,7 @@ bool TrackPar::rotateParam(float alpha)
   // RS: check if rotation does no invalidate track model (cos(local_phi)>=0, i.e. particle
   // direction in local frame is along the X axis
   if ((csp * ca + snp * sa) < 0) {
-    printf("Rotation failed: local cos(phi) would become %.2f\n", csp * ca + snp * sa);
+    //printf("Rotation failed: local cos(phi) would become %.2f\n", csp * ca + snp * sa);
     return false;
   }
   //
@@ -554,7 +554,7 @@ bool TrackParCov::rotate(float alpha)
   // RS: check if rotation does no invalidate track model (cos(local_phi)>=0, i.e. particle
   // direction in local frame is along the X axis
   if ((csp * ca + snp * sa) < 0) {
-    printf("Rotation failed: local cos(phi) would become %.2f\n", csp * ca + snp * sa);
+    //printf("Rotation failed: local cos(phi) would become %.2f\n", csp * ca + snp * sa);
     return false;
   }
   //

Detectors/ITSMFT/ITS/reconstruction/CMakeLists.txt

@@ -20,7 +20,6 @@ set(NO_DICT_SRCS # sources not for the dictionary
     src/CA/PrimaryVertexContext.cxx
     src/CA/Road.cxx
     src/CA/Tracker.cxx
-    src/CA/TrackingUtils.cxx
     src/CA/Tracklet.cxx
     src/CA/ClusterLines.cxx
     src/CA/Vertexer.cxx

Detectors/ITSMFT/ITS/reconstruction/include/ITSReconstruction/CA/Configuration.h

@@ -0,0 +1,112 @@
+// Copyright CERN and copyright holders of ALICE O2. This software is
+// distributed under the terms of the GNU General Public License v3 (GPL
+// Version 3), copied verbatim in the file "COPYING".
+//
+// See http://alice-o2.web.cern.ch/license for full licensing information.
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+///
+/// \file Configuration.h
+/// \brief
+///
+
+#ifndef TRACKINGITSU_INCLUDE_CONFIGURATION_H_
+#define TRACKINGITSU_INCLUDE_CONFIGURATION_H_
+
+#include <array>
+#include <climits>
+#include <vector>
+
+#include "json.h"
+
+#include "ITSReconstruction/CA/Constants.h"
+
+namespace o2
+{
+namespace ITS
+{
+namespace CA
+{
+
+template <typename Param>
+class Configuration : public Param
+{
+ public:
+  static Configuration<Param>& getInstance()
+  {
+    static Configuration<Param> instance;
+    return instance;
+  }
+
+ private:
+  Configuration() = default;
+  Configuration(const Configuration<Param>&) = delete;
+  const Configuration<Param>& operator=(const Configuration<Param>&) = delete;
+};
+
+struct TrackingParameters {
+  int CellMinimumLevel(int iteration);
+  int NumberOfIterations();
+
+  /// General parameters
+  int ClusterSharing = 0;
+  std::vector<int> MinTrackLength = { 7 };
+  /// Trackleting cuts
+  std::vector<float> TrackletMaxDeltaPhi = { 0.3f };
+  std::vector<std::array<float, Constants::ITS::TrackletsPerRoad>> TrackletMaxDeltaZ = { { 0.1f, 0.1f, 0.3f, 0.3f, 0.3f, 0.3f } };
+  /// Cell finding cuts
+  std::vector<float> CellMaxDeltaTanLambda = { 0.025f };
+  std::vector<std::array<float, Constants::ITS::CellsPerRoad>> CellMaxDCA = { { 0.05f, 0.04f, 0.05f, 0.2f, 0.4f } };
+  std::vector<float> CellMaxDeltaPhi = { 0.14f };
+  std::vector<std::array<float, Constants::ITS::CellsPerRoad>> CellMaxDeltaZ = { { 0.2f, 0.4f, 0.5f, 0.6f, 3.0f } };
+  /// Neighbour finding cuts
+  std::vector<std::array<float, Constants::ITS::CellsPerRoad - 1>> NeighbourMaxDeltaCurvature = { { 0.008f, 0.0025f, 0.003f, 0.0035f } };
+  std::vector<std::array<float, Constants::ITS::CellsPerRoad - 1>> NeighbourMaxDeltaN = { { 0.002f, 0.0090f, 0.002f, 0.005f } };
+};
+
+struct MemoryParameters {
+  /// Memory coefficients
+  int MemoryOffset = 256;
+  std::vector<std::array<float, Constants::ITS::CellsPerRoad>> CellsMemoryCoefficients = { { 2.3208e-08f, 2.104e-08f, 1.6432e-08f, 1.2412e-08f, 1.3543e-08f } };
+  std::vector<std::array<float, Constants::ITS::TrackletsPerRoad>> TrackletsMemoryCoefficients = { { 0.0016353f, 0.0013627f, 0.000984f, 0.00078135f, 0.00057934f, 0.00052217f } };
+};
+
+struct IndexTableParameters {
+  IndexTableParameters();
+  void ComputeInverseBinSizes();
+  int ZBins = 20;
+  int PhiBins = 20;
+  float InversePhiBinSize = 20 / Constants::Math::TwoPi;
+  std::array<float, Constants::ITS::LayersNumber> InverseZBinSize;
+};
+
+inline int TrackingParameters::NumberOfIterations()
+{
+  return MinTrackLength.size();
+}
+
+inline int TrackingParameters::CellMinimumLevel(int iteration)
+{
+  return MinTrackLength[iteration] - Constants::ITS::ClustersPerCell + 1;
+}
+
+inline IndexTableParameters::IndexTableParameters()
+{
+  ComputeInverseBinSizes();
+}
+
+inline void IndexTableParameters::ComputeInverseBinSizes()
+{
+  InversePhiBinSize = PhiBins / Constants::Math::TwoPi;
+  for (int iL = 0; iL < Constants::ITS::LayersNumber; ++iL) {
+    InverseZBinSize[iL] = 0.5f * ZBins / Constants::ITS::LayersZCoordinate()[iL];
+  }
+}
+
+} // namespace CA
+} // namespace ITS
+} // namespace o2
+
+#endif /* TRACKINGITSU_INCLUDE_CONFIGURATION_H_ */

Detectors/ITSMFT/ITS/reconstruction/include/ITSReconstruction/CA/Constants.h

@@ -16,6 +16,7 @@
 #define TRACKINGITSU_INCLUDE_CONSTANTS_H_
 
 #include <climits>
+#include <vector>
 
 #include "ITSReconstruction/CA/Definitions.h"
 
@@ -44,6 +45,7 @@ constexpr int LayersNumber{ 7 };
 constexpr int LayersNumberVertexer{ 3 };
 constexpr int TrackletsPerRoad{ LayersNumber - 1 };
 constexpr int CellsPerRoad{ LayersNumber - 2 };
+constexpr int ClustersPerCell{ 3 };
 constexpr int UnusedIndex{ -1 };
 constexpr float Resolution{ 0.0005f };
 
@@ -57,31 +59,6 @@ GPU_HOST_DEVICE constexpr GPUArray<float, LayersNumber> LayersRCoordinate()
 }
 } // namespace ITS
 
-namespace Thresholds
-{
-GPU_DEVICE constexpr GPUArray<float, ITS::TrackletsPerRoad> TrackletMaxDeltaZThreshold()
-{
-  return GPUArray<float, ITS::TrackletsPerRoad>{ { 0.1f, 0.1f, 0.3f, 0.3f, 0.3f, 0.3f } };
-}
-constexpr float CellMaxDeltaTanLambdaThreshold{ 0.025f };
-GPU_DEVICE constexpr GPUArray<float, ITS::CellsPerRoad> CellMaxDeltaZThreshold()
-{
-  return GPUArray<float, ITS::CellsPerRoad>{ { 0.2f, 0.4f, 0.5f, 0.6f, 3.0f } };
-}
-GPU_DEVICE constexpr GPUArray<float, ITS::CellsPerRoad> CellMaxDistanceOfClosestApproachThreshold()
-{
-  return GPUArray<float, ITS::CellsPerRoad>{ { 0.05f, 0.04f, 0.05f, 0.2f, 0.4f } };
-}
-constexpr float CellMaxDeltaPhiThreshold{ 0.14f };
-constexpr float ZCoordinateCut{ 0.5f };
-constexpr float PhiCoordinateCut{ 0.3f };
-constexpr GPUArray<float, ITS::CellsPerRoad - 1> NeighbourCellMaxNormalVectorsDelta{ { 0.002f, 0.009f, 0.002f,
-                                                                                       0.005f } };
-constexpr GPUArray<float, ITS::CellsPerRoad - 1> NeighbourCellMaxCurvaturesDelta{ { 0.008f, 0.0025f, 0.003f,
-                                                                                    0.0035f } };
-constexpr int CellsMinLevel{ 5 };
-} // namespace Thresholds
-
 namespace IndexTable
 {
 constexpr int ZBins{ 20 };
@@ -95,16 +72,6 @@ GPU_HOST_DEVICE constexpr GPUArray<float, ITS::LayersNumber> InverseZBinSize()
 }
 } // namespace IndexTable
 
-namespace Memory
-{
-constexpr GPUArray<float, ITS::TrackletsPerRoad> TrackletsMemoryCoefficients{
-  { 0.0016353f, 0.0013627f, 0.000984f, 0.00078135f, 0.00057934f, 0.00052217f }
-};
-constexpr GPUArray<float, ITS::CellsPerRoad> CellsMemoryCoefficients{ { 2.3208e-08f, 2.104e-08f, 1.6432e-08f,
-                                                                        1.2412e-08f, 1.3543e-08f } };
-constexpr int Offset = 256; /// Required for low multiplicity events
-} // namespace Memory
-
 namespace PDGCodes
 {
 constexpr int PionCode{ 211 };

Detectors/ITSMFT/ITS/reconstruction/include/ITSReconstruction/CA/Definitions.h

@@ -17,6 +17,17 @@
 
 #include <array>
 
+#ifdef CA_DEBUG
+#define CA_DEBUGGER(x) x
+#else
+#define CA_DEBUGGER(x) \
+  do {                 \
+  } while (0)
+#ifndef NDEBUG
+#define NDEBUG 1
+#endif
+#endif
+
 #if defined(TRACKINGITSU_CUDA_COMPILE)
 #define TRACKINGITSU_GPU_MODE true
 #else

Detectors/ITSMFT/ITS/reconstruction/include/ITSReconstruction/CA/IOUtils.h

@@ -20,7 +20,9 @@
 #include <unordered_map>
 #include <vector>
 
+#include "ITSReconstruction/CA/Configuration.h"
 #include "ITSReconstruction/CA/Event.h"
+#include "ITSReconstruction/CA/json.h"
 #include "ITSReconstruction/CA/Label.h"
 #include "ITSReconstruction/CA/Road.h"
 
@@ -45,8 +47,16 @@ namespace ITS
 namespace CA
 {
 
+void to_json(nlohmann::json& j, const TrackingParameters& par);
+void from_json(const nlohmann::json& j, TrackingParameters& par);
+void to_json(nlohmann::json& j, const MemoryParameters& par);
+void from_json(const nlohmann::json& j, MemoryParameters& par);
+void to_json(nlohmann::json& j, const IndexTableParameters& par);
+void from_json(const nlohmann::json& j, IndexTableParameters& par);
+
 namespace IOUtils
 {
+void loadConfigurations(const std::string&);
 std::vector<Event> loadEventData(const std::string&);
 void loadEventData(Event& events, const std::vector<ITSMFT::Cluster>* mClustersArray,
                    const dataformats::MCTruthContainer<MCCompLabel>* mClsLabels = nullptr);

Detectors/ITSMFT/ITS/reconstruction/include/ITSReconstruction/CA/MathUtils.h

@@ -33,6 +33,10 @@ float calculatePhiCoordinate(const float, const float);
 float calculateRCoordinate(const float, const float);
 GPU_HOST_DEVICE constexpr float getNormalizedPhiCoordinate(const float);
 GPU_HOST_DEVICE constexpr float3 crossProduct(const float3&, const float3&);
+float computeCurvature(float x1, float y1, float x2, float y2, float x3, float y3);
+float computeCurvatureCentreX(float x1, float y1, float x2, float y2, float x3, float y3);
+float computeTanDipAngle(float x1, float y1, float x2, float y2, float z1, float z2);
+
 } // namespace MathUtils
 
 inline float MathUtils::calculatePhiCoordinate(const float xCoordinate, const float yCoordinate)
@@ -59,6 +63,29 @@ GPU_HOST_DEVICE constexpr float3 MathUtils::crossProduct(const float3& firstVect
                  (firstVector.z * secondVector.x) - (firstVector.x * secondVector.z),
                  (firstVector.x * secondVector.y) - (firstVector.y * secondVector.x) };
 }
+
+inline float MathUtils::computeCurvature(float x1, float y1, float x2, float y2, float x3, float y3)
+{
+  const float d = (x2 - x1) * (y3 - y2) - (x3 - x2) * (y2 - y1);
+  const float a =
+    0.5f * ((y3 - y2) * (y2 * y2 - y1 * y1 + x2 * x2 - x1 * x1) - (y2 - y1) * (y3 * y3 - y2 * y2 + x3 * x3 - x2 * x2));
+  const float b =
+    0.5f * ((x2 - x1) * (y3 * y3 - y2 * y2 + x3 * x3 - x2 * x2) - (x3 - x2) * (y2 * y2 - y1 * y1 + x2 * x2 - x1 * x1));
+
+  return -1.f * d / std::sqrt((d * x1 - a) * (d * x1 - a) + (d * y1 - b) * (d * y1 - b));
+}
+
+inline float MathUtils::computeCurvatureCentreX(float x1, float y1, float x2, float y2, float x3, float y3)
+{
+  const float k1 = (y2 - y1) / (x2 - x1), k2 = (y3 - y2) / (x3 - x2);
+  return 0.5f * (k1 * k2 * (y1 - y3) + k2 * (x1 + x2) - k1 * (x2 + x3)) / (k2 - k1);
+}
+
+inline float MathUtils::computeTanDipAngle(float x1, float y1, float x2, float y2, float z1, float z2)
+{
+  return (z1 - z2) / std::sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
+}
+
 } // namespace CA
 } // namespace ITS
 } // namespace o2

Detectors/ITSMFT/ITS/reconstruction/include/ITSReconstruction/CA/PrimaryVertexContext.h

@@ -24,6 +24,7 @@
 #include "SimulationDataFormat/MCCompLabel.h"
 
 #include "ITSReconstruction/CA/Cell.h"
+#include "ITSReconstruction/CA/Configuration.h"
 #include "ITSReconstruction/CA/Constants.h"
 #include "ITSReconstruction/CA/Definitions.h"
 #include "ITSReconstruction/CA/Event.h"
@@ -50,7 +51,7 @@ class PrimaryVertexContext final
   PrimaryVertexContext(const PrimaryVertexContext&) = delete;
   PrimaryVertexContext& operator=(const PrimaryVertexContext&) = delete;
 
-  void initialise(const Event&, const int);
+  void initialise(const MemoryParameters& memParam, const Event& event, const int pvIndex, const int iteration);
   const float3& getPrimaryVertex() const;
   std::array<std::vector<Cluster>, Constants::ITS::LayersNumber>& getClusters();
   std::array<std::vector<Cell>, Constants::ITS::CellsPerRoad>& getCells();
@@ -144,7 +145,10 @@ inline bool PrimaryVertexContext::isClusterUsed(int layer, int clusterId) const
 inline void PrimaryVertexContext::markUsedCluster(int layer, int clusterId) { mUsedClusters[layer][clusterId] = true; }
 
 #if TRACKINGITSU_GPU_MODE
-inline GPU::PrimaryVertexContext& PrimaryVertexContext::getDeviceContext() { return *mGPUContextDevicePointer; }
+inline GPU::PrimaryVertexContext& PrimaryVertexContext::getDeviceContext()
+{
+  return *mGPUContextDevicePointer;
+}
 
 inline GPU::Array<GPU::Vector<Cluster>, Constants::ITS::LayersNumber>& PrimaryVertexContext::getDeviceClusters()
 {