ci: Add reco, reco after vertex fit, and truth track momenta to physm…

…on (#2221) Comparing the reconstructed momenta to their truth values in the physics monitoring. Both the momenta before and after the vertex fit are considered.
acts-project · Aug 14, 2023 · d0759c0 · d0759c0
1 parent e9977db
commit d0759c0
Show file tree

Hide file tree

Showing 15 changed files with 471 additions and 127 deletions.
diff --git a/CI/physmon/reference/performance_amvf_orthogonal_hist.root b/CI/physmon/reference/performance_amvf_orthogonal_hist.root
diff --git a/CI/physmon/reference/performance_amvf_seeded_hist.root b/CI/physmon/reference/performance_amvf_seeded_hist.root
diff --git a/CI/physmon/reference/performance_amvf_truth_estimated_hist.root b/CI/physmon/reference/performance_amvf_truth_estimated_hist.root
diff --git a/CI/physmon/reference/performance_amvf_truth_smeared_hist.root b/CI/physmon/reference/performance_amvf_truth_smeared_hist.root
diff --git a/CI/physmon/reference/performance_amvf_ttbar_hist.root b/CI/physmon/reference/performance_amvf_ttbar_hist.root
diff --git a/CI/physmon/reference/performance_ivf_orthogonal_hist.root b/CI/physmon/reference/performance_ivf_orthogonal_hist.root
diff --git a/CI/physmon/reference/performance_ivf_seeded_hist.root b/CI/physmon/reference/performance_ivf_seeded_hist.root
diff --git a/CI/physmon/reference/performance_ivf_truth_estimated_hist.root b/CI/physmon/reference/performance_ivf_truth_estimated_hist.root
diff --git a/CI/physmon/reference/performance_ivf_truth_smeared_hist.root b/CI/physmon/reference/performance_ivf_truth_smeared_hist.root
diff --git a/CI/physmon/vertexing_config.yml b/CI/physmon/vertexing_config.yml
@@ -4,11 +4,22 @@ histograms:
     min: -0.1
     max: 0.1
 
+  "trk_res.*":
+    nbins: 100
+    min: -0.1
+    max: 0.1
+
   "pull.*":
     nbins: 50
     min: -6
     max: 6
 
+
+  "trk_pull.*":
+    nbins: 50
+    min: -10
+    max: 10
+
   "cov.*":
     nbins: 100
     min: -0.0005
@@ -34,6 +45,26 @@ histograms:
     min: -2000
     max: 2000
 
+  "trk_truthPhi|trk_recoPhi|trk_recoPhiFitted":
+    nbins: 100
+    min: -3.2
+    max: 3.2
+
+  "trk_truthTheta|trk_recoTheta|trk_recoThetaFitted":
+    nbins: 100
+    min: -0.1
+    max: 3.1
+
+  "trk_truthQOverP|trk_recoQOverP|trk_recoQOverPFitted":
+    nbins: 100
+    min: -1
+    max: 1
+
+  "trk_momOverlap|trk_momOverlapFitted":
+    nbins: 100
+    min: 0.99
+    max: 1
+
 extra_histograms:
   - expression: df["nRecoVtx"] / df["nTrueVtx"]
     name: "recoOverTrue"

diff --git a/Core/include/Acts/Utilities/UnitVectors.hpp b/Core/include/Acts/Utilities/UnitVectors.hpp
@@ -55,6 +55,22 @@ inline Eigen::Matrix<T, 3, 1> makeDirectionUnitFromPhiTheta(T phi, T theta) {
   };
 }
 
+/// Construct a phi and theta angle from a direction vector.
+///
+/// @param unitDir 3D vector indicating a direction
+///
+template <typename T>
+inline Eigen::Matrix<T, 2, 1> makePhiThetaFromDirectionUnit(
+    Eigen::Matrix<T, 3, 1> unitDir) {
+  unitDir.normalize();
+  T phi = std::atan2(unitDir[1], unitDir[0]);
+  T theta = std::acos(unitDir[2]);
+  return {
+      phi,
+      theta,
+  };
+}
+
 /// Construct the first curvilinear unit vector `U` for the given direction.
 ///
 /// @param direction is the input direction vector

diff --git a/Examples/Io/Performance/ActsExamples/Io/Performance/VertexPerformanceWriter.cpp b/Examples/Io/Performance/ActsExamples/Io/Performance/VertexPerformanceWriter.cpp
diff --git a/Examples/Io/Performance/ActsExamples/Io/Performance/VertexPerformanceWriter.hpp b/Examples/Io/Performance/ActsExamples/Io/Performance/VertexPerformanceWriter.hpp
@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2021 CERN for the benefit of the Acts project
+// Copyright (C) 2019-2023 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
@@ -10,6 +10,7 @@
 
 #include "Acts/Definitions/TrackParametrization.hpp"
 #include "Acts/EventData/TrackParameters.hpp"
+#include "Acts/MagneticField/MagneticFieldProvider.hpp"
 #include "Acts/Utilities/Logger.hpp"
 #include "Acts/Vertexing/Vertex.hpp"
 #include "ActsExamples/EventData/Index.hpp"
@@ -60,6 +61,8 @@ class VertexPerformanceWriter final
     std::string inputMeasurementParticlesMap;
     /// Input vertex collection.
     std::string inputVertices;
+    /// Magnetic field
+    std::shared_ptr<Acts::MagneticFieldProvider> bField;
     /// Output filename.
     std::string filePath = "vertexingperformance.root";
     /// Name of the output tree.
@@ -104,42 +107,89 @@ class VertexPerformanceWriter final
   TFile* m_outputFile{nullptr};  ///< The output file
   TTree* m_outputTree{nullptr};  ///< The output tree
 
-  std::vector<float> m_truthX;
-  std::vector<float> m_truthY;
-  std::vector<float> m_truthZ;
-  std::vector<float> m_truthT;
-
-  std::vector<float> m_recoX;
-  std::vector<float> m_recoY;
-  std::vector<float> m_recoZ;
-  std::vector<float> m_recoT;
-
-  /// Difference in x position between reco and true vtx
-  std::vector<float> m_resX;
-  std::vector<float> m_resY;
-  std::vector<float> m_resZ;
-  std::vector<float> m_resT;
-
-  std::vector<float> m_pullX;
-  std::vector<float> m_pullY;
-  std::vector<float> m_pullZ;
-  std::vector<float> m_pullT;
-
-  std::vector<float> m_covXX;
-  std::vector<float> m_covYY;
-  std::vector<float> m_covZZ;
-  std::vector<float> m_covTT;
-  std::vector<float> m_covXY;
-  std::vector<float> m_covXZ;
-  std::vector<float> m_covXT;
-  std::vector<float> m_covYZ;
-  std::vector<float> m_covYT;
-  std::vector<float> m_covZT;
-
+  // True 4D vertex position
+  std::vector<double> m_truthX;
+  std::vector<double> m_truthY;
+  std::vector<double> m_truthZ;
+  std::vector<double> m_truthT;
+
+  // Reconstructed 4D vertex position
+  std::vector<double> m_recoX;
+  std::vector<double> m_recoY;
+  std::vector<double> m_recoZ;
+  std::vector<double> m_recoT;
+
+  /// Difference of reconstructed and true vertex 4D position
+  std::vector<double> m_resX;
+  std::vector<double> m_resY;
+  std::vector<double> m_resZ;
+  std::vector<double> m_resT;
+
+  // pull(X) = (X_reco - X_true)/Var(X_reco)^(1/2)
+  std::vector<double> m_pullX;
+  std::vector<double> m_pullY;
+  std::vector<double> m_pullZ;
+  std::vector<double> m_pullT;
+
+  // Vertex covariance
+  std::vector<double> m_covXX;
+  std::vector<double> m_covYY;
+  std::vector<double> m_covZZ;
+  std::vector<double> m_covTT;
+  std::vector<double> m_covXY;
+  std::vector<double> m_covXZ;
+  std::vector<double> m_covXT;
+  std::vector<double> m_covYZ;
+  std::vector<double> m_covYT;
+  std::vector<double> m_covZT;
+
+  //--------------------------------------------------------------
+  // Track-related variables are contained in a vector of vectors: The inner
+  // vectors contain the values of all tracks corresponding to one vertex. The
+  // outer vector can then have the same length as the flat vectors of
+  // vertex-related variables (see above). E.g.,
+  // m_truthPhi = ((truthPhi of 1st trk belonging to vtx 1,
+  //                truthPhi of 2nd trk belonging to vtx 1, ...),
+  //               (truthPhi of 1st trk belonging to vtx 2,
+  //                truthPhi of 2nd trk belonging to vtx 2, ...),
+  //                ...)
+  //
+  // True track momenta at the vertex
+  std::vector<std::vector<double>> m_truthPhi;
+  std::vector<std::vector<double>> m_truthTheta;
+  std::vector<std::vector<double>> m_truthQOverP;
+
+  // Reconstructed track momenta at the vertex before and after the vertex fit
+  std::vector<std::vector<double>> m_recoPhi;
+  std::vector<std::vector<double>> m_recoPhiFitted;
+  std::vector<std::vector<double>> m_recoTheta;
+  std::vector<std::vector<double>> m_recoThetaFitted;
+  std::vector<std::vector<double>> m_recoQOverP;
+  std::vector<std::vector<double>> m_recoQOverPFitted;
+
+  // Difference between reconstructed momenta and true momenta
+  std::vector<std::vector<double>> m_resPhi;
+  std::vector<std::vector<double>> m_resPhiFitted;
+  std::vector<std::vector<double>> m_resTheta;
+  std::vector<std::vector<double>> m_resThetaFitted;
+  std::vector<std::vector<double>> m_resQOverP;
+  std::vector<std::vector<double>> m_resQOverPFitted;
+  std::vector<std::vector<double>> m_momOverlap;
+  std::vector<std::vector<double>> m_momOverlapFitted;
+
+  // Pulls
+  std::vector<std::vector<double>> m_pullPhi;
+  std::vector<std::vector<double>> m_pullPhiFitted;
+  std::vector<std::vector<double>> m_pullTheta;
+  std::vector<std::vector<double>> m_pullThetaFitted;
+  std::vector<std::vector<double>> m_pullQOverP;
+  std::vector<std::vector<double>> m_pullQOverPFitted;
+
+  // Number of tracks associated with truth/reconstructed vertex
   std::vector<int> m_nTracksOnTruthVertex;
   std::vector<int> m_nTracksOnRecoVertex;
 
-  std::vector<float> m_trackVtxMatchFraction;
+  std::vector<double> m_trackVtxMatchFraction;
 
   /// Number of reconstructed vertices
   int m_nRecoVtx = -1;

diff --git a/Examples/Python/python/acts/examples/reconstruction.py b/Examples/Python/python/acts/examples/reconstruction.py
@@ -1554,6 +1554,7 @@ def addVertexFitting(
                 inputTrackParameters=trackParameters,
                 inputAssociatedTruthParticles=associatedParticles,
                 inputVertices=outputVertices,
+                bField=field,
                 minTrackVtxMatchFraction=0.5 if associatedParticles else 0.0,
                 treeName="vertexing",
                 filePath=str(outputDirRoot / "performance_vertexing.root"),

diff --git a/Examples/Python/src/Output.cpp b/Examples/Python/src/Output.cpp
@@ -312,8 +312,8 @@ void addOutput(Context& ctx) {
       ActsExamples::VertexPerformanceWriter, mex, "VertexPerformanceWriter",
       inputAllTruthParticles, inputSelectedTruthParticles, inputTrackParameters,
       inputAssociatedTruthParticles, inputTrackParameters, inputTrajectories,
-      inputMeasurementParticlesMap, inputVertices, filePath, treeName, fileMode,
-      minTrackVtxMatchFraction, truthMatchProbMin, useTracks);
+      inputMeasurementParticlesMap, inputVertices, bField, filePath, treeName,
+      fileMode, minTrackVtxMatchFraction, truthMatchProbMin, useTracks);
 
   // CSV WRITERS
   ACTS_PYTHON_DECLARE_WRITER(ActsExamples::CsvParticleWriter, mex,