feat: save vertex seed (acts-project#2885)

Saves 4D of the vertex seed (its x- and y- coordinate are not estimated in current seeding algorithms and should be set to 0).

Co-authored-by: Andreas Stefl <487211+andiwand@users.noreply.github.com>

CI/physmon/vertexing_config.yml

@@ -14,12 +14,12 @@ histograms:
     min: 0.3
     max: 1.01
 
-  "resX|resY|resZ":
+  "resX|resY|resZ|resSeedZ":
     nbins: 100
     min: -0.1
     max: 0.1
 
-  "resT":
+  "resT|resSeedT":
     nbins: 100
     min: -50
     max: 50
@@ -50,17 +50,17 @@ histograms:
     min: 0.0
     max: 100.0
 
-  "truthX|truthY|recoX|recoY":
+  "truthX|truthY|recoX|recoY|seedX|seedY":
     nbins: 100
     min: -0.06
     max: 0.06
 
-  "truthZ|recoZ":
+  "truthZ|recoZ|seedZ":
     nbins: 100
     min: -200
     max: 200
 
-  "truthT|recoT":
+  "truthT|recoT|seedT":
     nbins: 100
     min: -2000
     max: 2000

CI/physmon/vertexing_ttbar_config.yml

@@ -14,12 +14,12 @@ histograms:
     min: 0.0
     max: 1.01
 
-  "resX|resY|resZ":
+  "resX|resY|resZ|resSeedZ":
     nbins: 100
     min: -0.2
     max: 0.2
 
-  "resT":
+  "resT|resSeedT":
     nbins: 100
     min: -50
     max: 50
@@ -49,17 +49,17 @@ histograms:
     min: 0.0
     max: 100.0
 
-  "truthX|truthY|recoX|recoY":
+  "truthX|truthY|recoX|recoY|seedX|seedY":
     nbins: 100
     min: -0.1
     max: 0.1
 
-  "truthZ|recoZ":
+  "truthZ|recoZ|seedZ":
     nbins: 100
     min: -200
     max: 200
 
-  "truthT|recoT":
+  "truthT|recoT|seedT":
     nbins: 100
     min: -2000
     max: 2000

Core/include/Acts/Vertexing/Vertex.hpp

@@ -56,6 +56,10 @@ class Vertex {
   const Vector4& fullPosition() const;
   Vector4& fullPosition();
 
+  /// @return Returns 4D position of the vertex seed
+  const Vector4& fullSeedPosition() const;
+  Vector4& fullSeedPosition();
+
   /// @return Returns position covariance
   SquareMatrix3 covariance() const;
 
@@ -107,6 +111,7 @@ class Vertex {
 
  private:
   Vector4 m_position = Vector4::Zero();
+  Vector4 m_seedPosition = Vector4::Zero();
   SquareMatrix4 m_covariance = SquareMatrix4::Zero();
   std::vector<TrackAtVertex> m_tracksAtVertex;
   double m_chiSquared = 0.;  // chi2 of the fit

Core/src/Vertexing/Vertex.cpp

@@ -8,99 +8,109 @@
 
 #include "Acts/Vertexing/Vertex.hpp"
 
-Acts::Vertex::Vertex(const Vector3& position) {
-  m_position[ePos0] = position[ePos0];
-  m_position[ePos1] = position[ePos1];
-  m_position[ePos2] = position[ePos2];
+namespace Acts {
+
+Vertex::Vertex(const Vector3& position) {
+  m_position.head<3>() = position;
+  m_seedPosition.head<3>() = position;
 }
 
-Acts::Vertex::Vertex(const Vector4& position) : m_position(position) {}
+Vertex::Vertex(const Vector4& position)
+    : m_position(position), m_seedPosition(position) {}
 
-Acts::Vertex::Vertex(const Vector3& position, const SquareMatrix3& covariance,
-                     std::vector<TrackAtVertex> tracks)
+Vertex::Vertex(const Vector3& position, const SquareMatrix3& covariance,
+               std::vector<TrackAtVertex> tracks)
     : m_tracksAtVertex(std::move(tracks)) {
-  m_position[ePos0] = position[ePos0];
-  m_position[ePos1] = position[ePos1];
-  m_position[ePos2] = position[ePos2];
+  m_position.head<3>() = position;
+  m_seedPosition.head<3>() = position;
   m_covariance.block<3, 3>(ePos0, ePos0) = covariance;
 }
 
-Acts::Vertex::Vertex(const Vector4& position, const SquareMatrix4& covariance,
-                     std::vector<TrackAtVertex> tracks)
+Vertex::Vertex(const Vector4& position, const SquareMatrix4& covariance,
+               std::vector<TrackAtVertex> tracks)
     : m_position(position),
+      m_seedPosition(position),
       m_covariance(covariance),
       m_tracksAtVertex(std::move(tracks)) {}
 
-Acts::Vector3 Acts::Vertex::position() const {
+Vector3 Vertex::position() const {
   return VectorHelpers::position(m_position);
 }
 
-Acts::ActsScalar Acts::Vertex::time() const {
+ActsScalar Vertex::time() const {
   return m_position[eTime];
 }
 
-const Acts::Vector4& Acts::Vertex::fullPosition() const {
+const Vector4& Vertex::fullPosition() const {
   return m_position;
 }
 
-Acts::Vector4& Acts::Vertex::fullPosition() {
+Vector4& Vertex::fullPosition() {
   return m_position;
 }
 
-Acts::SquareMatrix3 Acts::Vertex::covariance() const {
+const Vector4& Vertex::fullSeedPosition() const {
+  return m_seedPosition;
+}
+
+Vector4& Vertex::fullSeedPosition() {
+  return m_seedPosition;
+}
+
+SquareMatrix3 Vertex::covariance() const {
   return m_covariance.block<3, 3>(ePos0, ePos0);
 }
 
-const Acts::SquareMatrix4& Acts::Vertex::fullCovariance() const {
+const SquareMatrix4& Vertex::fullCovariance() const {
   return m_covariance;
 }
 
-Acts::SquareMatrix4& Acts::Vertex::fullCovariance() {
+SquareMatrix4& Vertex::fullCovariance() {
   return m_covariance;
 }
 
-const std::vector<Acts::TrackAtVertex>& Acts::Vertex::tracks() const {
+const std::vector<TrackAtVertex>& Vertex::tracks() const {
   return m_tracksAtVertex;
 }
 
-std::pair<double, double> Acts::Vertex::fitQuality() const {
+std::pair<double, double> Vertex::fitQuality() const {
   return std::pair<double, double>(m_chiSquared, m_numberDoF);
 }
 
-void Acts::Vertex::setPosition(const Vector3& position, ActsScalar time) {
-  m_position[ePos0] = position[ePos0];
-  m_position[ePos1] = position[ePos1];
-  m_position[ePos2] = position[ePos2];
+void Vertex::setPosition(const Vector3& position, ActsScalar time) {
+  m_position.head<3>() = position;
   m_position[eTime] = time;
 }
 
-void Acts::Vertex::setFullPosition(const Vector4& fullPosition) {
+void Vertex::setFullPosition(const Vector4& fullPosition) {
   m_position = fullPosition;
 }
 
-void Acts::Vertex::setTime(ActsScalar time) {
+void Vertex::setTime(ActsScalar time) {
   m_position[eTime] = time;
 }
 
-void Acts::Vertex::setCovariance(const SquareMatrix3& covariance) {
+void Vertex::setCovariance(const SquareMatrix3& covariance) {
   m_covariance.setZero();
   m_covariance.block<3, 3>(ePos0, ePos0) = covariance;
 }
 
-void Acts::Vertex::setFullCovariance(const SquareMatrix4& covariance) {
+void Vertex::setFullCovariance(const SquareMatrix4& covariance) {
   m_covariance = covariance;
 }
 
-void Acts::Vertex::setTracksAtVertex(std::vector<TrackAtVertex> tracks) {
+void Vertex::setTracksAtVertex(std::vector<TrackAtVertex> tracks) {
   m_tracksAtVertex = std::move(tracks);
 }
 
-void Acts::Vertex::setFitQuality(double chiSquared, double numberDoF) {
+void Vertex::setFitQuality(double chiSquared, double numberDoF) {
   m_chiSquared = chiSquared;
   m_numberDoF = numberDoF;
 }
 
-void Acts::Vertex::setFitQuality(std::pair<double, double> fitQuality) {
+void Vertex::setFitQuality(std::pair<double, double> fitQuality) {
   m_chiSquared = fitQuality.first;
   m_numberDoF = fitQuality.second;
 }
+
+}  // namespace Acts

Examples/Io/Performance/ActsExamples/Io/Performance/VertexPerformanceWriter.cpp

@@ -108,11 +108,19 @@ ActsExamples::VertexPerformanceWriter::VertexPerformanceWriter(
     m_outputTree->Branch("recoZ", &m_recoZ);
     m_outputTree->Branch("recoT", &m_recoT);
 
+    m_outputTree->Branch("seedX", &m_seedX);
+    m_outputTree->Branch("seedY", &m_seedY);
+    m_outputTree->Branch("seedZ", &m_seedZ);
+    m_outputTree->Branch("seedT", &m_seedT);
+
     m_outputTree->Branch("resX", &m_resX);
     m_outputTree->Branch("resY", &m_resY);
     m_outputTree->Branch("resZ", &m_resZ);
     m_outputTree->Branch("resT", &m_resT);
 
+    m_outputTree->Branch("resSeedZ", &m_resSeedZ);
+    m_outputTree->Branch("resSeedT", &m_resSeedT);
+
     m_outputTree->Branch("pullX", &m_pullX);
     m_outputTree->Branch("pullY", &m_pullY);
     m_outputTree->Branch("pullZ", &m_pullZ);
@@ -554,12 +562,26 @@ ActsExamples::ProcessCode ActsExamples::VertexPerformanceWriter::writeT(
             m_recoZ.push_back(vtx.fullPosition()[Acts::FreeIndices::eFreePos2]);
             m_recoT.push_back(vtx.fullPosition()[Acts::FreeIndices::eFreeTime]);
 
+            m_seedX.push_back(
+                vtx.fullSeedPosition()[Acts::FreeIndices::eFreePos0]);
+            m_seedY.push_back(
+                vtx.fullSeedPosition()[Acts::FreeIndices::eFreePos1]);
+            m_seedZ.push_back(
+                vtx.fullSeedPosition()[Acts::FreeIndices::eFreePos2]);
+            m_seedT.push_back(
+                vtx.fullSeedPosition()[Acts::FreeIndices::eFreeTime]);
+
             const Acts::ActsVector<4> diffPos = vtx.fullPosition() - truePos;
             m_resX.push_back(diffPos[Acts::FreeIndices::eFreePos0]);
             m_resY.push_back(diffPos[Acts::FreeIndices::eFreePos1]);
             m_resZ.push_back(diffPos[Acts::FreeIndices::eFreePos2]);
             m_resT.push_back(diffPos[Acts::FreeIndices::eFreeTime]);
 
+            const Acts::ActsVector<4> diffSeedPos =
+                vtx.fullSeedPosition() - truePos;
+            m_resSeedZ.push_back(diffSeedPos[Acts::FreeIndices::eFreePos2]);
+            m_resSeedT.push_back(diffSeedPos[Acts::FreeIndices::eFreeTime]);
+
             Acts::ActsScalar varX = vtx.fullCovariance()(
                 Acts::FreeIndices::eFreePos0, Acts::FreeIndices::eFreePos0);
             Acts::ActsScalar varY = vtx.fullCovariance()(
@@ -742,10 +764,16 @@ ActsExamples::ProcessCode ActsExamples::VertexPerformanceWriter::writeT(
   m_recoY.clear();
   m_recoZ.clear();
   m_recoT.clear();
+  m_seedX.clear();
+  m_seedY.clear();
+  m_seedZ.clear();
+  m_seedT.clear();
   m_resX.clear();
   m_resY.clear();
   m_resZ.clear();
   m_resT.clear();
+  m_resSeedZ.clear();
+  m_resSeedT.clear();
   m_pullX.clear();
   m_pullY.clear();
   m_pullZ.clear();

Examples/Io/Performance/ActsExamples/Io/Performance/VertexPerformanceWriter.hpp

@@ -119,12 +119,23 @@ class VertexPerformanceWriter final
   std::vector<double> m_recoZ;
   std::vector<double> m_recoT;
 
+  // 4D position of the vertex seed. x and y coordinate are 0 in current
+  // implementations, we save them here as a check.
+  std::vector<double> m_seedX;
+  std::vector<double> m_seedY;
+  std::vector<double> m_seedZ;
+  std::vector<double> m_seedT;
+
   // 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;
 
+  // Difference between the seed and the true vertex z and t coordinate
+  std::vector<double> m_resSeedZ;
+  std::vector<double> m_resSeedT;
+
   // pull(X) = (X_reco - X_true)/Var(X_reco)^(1/2)
   std::vector<double> m_pullX;
   std::vector<double> m_pullY;