diff --git a/PhysicsAnalysis/NtupleDumper/CMakeLists.txt b/PhysicsAnalysis/NtupleDumper/CMakeLists.txt
index 607c58d892d201731098f64829a82d556f1a302a..dd8faf032f48b2c3d78453b84494e1903d64367d 100644
--- a/PhysicsAnalysis/NtupleDumper/CMakeLists.txt
+++ b/PhysicsAnalysis/NtupleDumper/CMakeLists.txt
@@ -7,7 +7,7 @@ atlas_add_component(
src/NtupleDumperAlg.h
src/NtupleDumperAlg.cxx
src/component/NtupleDumper_entries.cxx
- LINK_LIBRARIES AthenaBaseComps StoreGateLib xAODFaserWaveform xAODFaserCalorimeter xAODFaserTrigger xAODFaserLHC ScintIdentifier FaserCaloIdentifier GeneratorObjects FaserActsGeometryLib TrackerSimEvent TrackerSimData TrackerIdentifier TrackerReadoutGeometry TrkTrack GeoPrimitives TrackerRIO_OnTrack TrackerSpacePoint FaserActsKalmanFilterLib
+ LINK_LIBRARIES AthenaBaseComps StoreGateLib xAODFaserWaveform xAODFaserCalorimeter xAODFaserTrigger xAODFaserLHC ScintIdentifier FaserCaloIdentifier GeneratorObjects FaserActsGeometryLib TrackerSimEvent TrackerSimData TrackerIdentifier TrackerReadoutGeometry TrkTrack GeoPrimitives TrackerRIO_OnTrack TrackerSpacePoint FaserActsKalmanFilterLib FaserActsmanVertexingLib
PRIVATE_LINK_LIBRARIES nlohmann_json::nlohmann_json
)
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
index 6e4c2aae531c6feee89df6102a90d118f6276dff..1568097245d10680bb9a9889b13ceb5b854f9e51 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
@@ -17,7 +17,9 @@
#include <cmath>
#include <TH1F.h>
#include <numeric>
+#include <limits>
+constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
using namespace Acts::UnitLiterals;
NtupleDumperAlg::NtupleDumperAlg(const std::string &name,
@@ -119,6 +121,7 @@ StatusCode NtupleDumperAlg::initialize()
ATH_CHECK(m_trackingGeometryTool.retrieve());
ATH_CHECK(m_trackTruthMatchingTool.retrieve());
ATH_CHECK(m_fiducialParticleTool.retrieve());
+ ATH_CHECK(m_vertexingTool.retrieve());
ATH_CHECK(m_spacePointContainerKey.initialize());
@@ -289,6 +292,8 @@ StatusCode NtupleDumperAlg::initialize()
//TRUTH
m_tree->Branch("t_pdg", &m_t_pdg);
m_tree->Branch("t_barcode", &m_t_barcode);
+ m_tree->Branch("t_pdg_parent", &m_t_pdg_parent);
+ m_tree->Branch("t_barcode_parent", &m_t_barcode_parent);
m_tree->Branch("t_truthHitRatio", &m_t_truthHitRatio);
m_tree->Branch("t_prodVtx_x", &m_t_prodVtx_x);
@@ -375,6 +380,19 @@ StatusCode NtupleDumperAlg::initialize()
m_tree->Branch("truthd1_y", &m_truthd1_y);
m_tree->Branch("truthd1_z", &m_truthd1_z);
+ m_tree->Branch("vertex_x", &m_vertex_x, "vertex_x/D");
+ m_tree->Branch("vertex_y", &m_vertex_y, "vertex_y/D");
+ m_tree->Branch("vertex_z", &m_vertex_z, "vertex_z/D");
+ m_tree->Branch("vertex_chi2", &m_vertex_chi2, "vertex_chi2/D");
+ m_tree->Branch("vertex_px0", &m_vertex_px0, "vertex_px0/D");
+ m_tree->Branch("vertex_py0", &m_vertex_py0, "vertex_py0/D");
+ m_tree->Branch("vertex_pz0", &m_vertex_pz0, "vertex_pz0/D");
+ m_tree->Branch("vertex_p0", &m_vertex_p0, "vertex_p0/D");
+ m_tree->Branch("vertex_px1", &m_vertex_px1, "vertex_px1/D");
+ m_tree->Branch("vertex_py1", &m_vertex_py1, "vertex_py1/D");
+ m_tree->Branch("vertex_pz1", &m_vertex_pz1, "vertex_pz1/D");
+ m_tree->Branch("vertex_p1", &m_vertex_p1, "vertex_p1/D");
+
ATH_CHECK(histSvc()->regTree("/HIST2/tree", m_tree));
// Register histograms
@@ -885,6 +903,8 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
}
ATH_CHECK(trackCollection.isValid());
+ // create list of good tracks in fiducial region (r < 95mm) and chi2/nDoF < 25, #Layers >= 7, #Hits >= 12
+ std::vector<const Trk::Track*> goodTracks {};
for (const Trk::Track* track : *trackCollection)
{
if (track == nullptr) continue;
@@ -951,6 +971,13 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_y_atMaxRadius.push_back(positionAtMaxRadius.y());
m_z_atMaxRadius.push_back(positionAtMaxRadius.z());
+ if ((m_Chi2.back() / m_DoF.back() < m_maxChi2NDOF) &&
+ (m_nLayers.back() >= m_minLayers) &&
+ (m_DoF.back() + 5 >= m_minHits) &&
+ (m_r_atMaxRadius.back() < 95)) {
+ goodTracks.push_back(track);
+ }
+
if (isMC) { // if simulation, store track truth info as well
auto [truthParticle, hitCount] = m_trackTruthMatchingTool->getTruthParticle(track);
if (truthParticle != nullptr) {
@@ -958,6 +985,13 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
truthParticleCount[truthParticle->barcode()] += 1;
m_t_pdg.push_back(truthParticle->pdgId());
m_t_barcode.push_back(truthParticle->barcode());
+ if (truthParticle->nParents() > 0){
+ m_t_pdg_parent.push_back(truthParticle->parent(0)->pdgId());
+ m_t_barcode_parent.push_back(truthParticle->parent(0)->barcode());
+ } else {
+ m_t_pdg_parent.push_back(0);
+ m_t_barcode_parent.push_back(-1);
+ }
// the track fit eats up 5 degrees of freedom, thus the number of hits on track is m_DoF + 5
m_t_truthHitRatio.push_back(hitCount / (m_DoF.back() + 5));
m_isFiducial.push_back(m_fiducialParticleTool->isFiducial(truthParticle->barcode()));
@@ -1152,6 +1186,41 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_longTracks++;
}
+ if (goodTracks.size() >= 2) {
+ // sort tracks my momentum and reconstruct vertex
+ std::sort(goodTracks.begin(), goodTracks.end(), [](const Trk::Track *lhs, const Trk::Track *rhs){
+ return lhs->trackParameters()->front()->momentum().z() < rhs->trackParameters()->front()->momentum().z();
+ });
+ std::vector<const Trk::Track*> lowMomentumTracks {goodTracks[0], goodTracks[1]};
+
+ std::optional<FaserTracking::POCA> vertex = m_vertexingTool->getVertex(lowMomentumTracks);
+ if (vertex) {
+ Eigen::Vector3d position = vertex->position;
+ m_vertex_x = position.x();
+ m_vertex_y = position.y();
+ m_vertex_z = position.z();
+ m_vertex_chi2 = vertex->chi2;
+
+ // get track parameters of two tracks with lowest momentum at vertex
+ auto vertexTrackParameters0 = m_vertexingTool->extrapolateTrack(goodTracks[0], m_vertex_z);
+ if (vertexTrackParameters0 != nullptr) {
+ // convert momentum from GeV to MeV
+ m_vertex_px0 = vertexTrackParameters0->momentum().x() * 1e3;
+ m_vertex_py0 = vertexTrackParameters0->momentum().y() * 1e3;
+ m_vertex_pz0 = vertexTrackParameters0->momentum().z() * 1e3;
+ m_vertex_p0 = std::sqrt(m_vertex_px0*m_vertex_px0 + m_vertex_py0*m_vertex_py0 + m_vertex_pz0*m_vertex_pz0);
+ }
+ auto vertexTrackParameters1 = m_vertexingTool->extrapolateTrack(goodTracks[1], m_vertex_z);
+ if (vertexTrackParameters1 != nullptr) {
+ // convert momentum from GeV to MeV
+ m_vertex_px1 = vertexTrackParameters1->momentum().x() * 1e3;
+ m_vertex_py1 = vertexTrackParameters1->momentum().y() * 1e3;
+ m_vertex_pz1 = vertexTrackParameters1->momentum().z() * 1e3;
+ m_vertex_p1 = std::sqrt(m_vertex_px1*m_vertex_px1 + m_vertex_py1*m_vertex_py1 + m_vertex_pz1*m_vertex_pz1);
+ }
+ }
+ }
+
if (!isMC) {
if (m_doTrackFilter) { // filter events: colliding bunches have at least one long track, non-colliding bunches have at least one long track or one calo module with raw_peak > 3 mV
if (m_longTracks == 0 && abs(m_distanceToCollidingBCID) <= 1) {
@@ -1400,6 +1469,8 @@ NtupleDumperAlg::clearTree() const
m_t_pdg.clear();
m_t_barcode.clear();
+ m_t_pdg_parent.clear();
+ m_t_barcode_parent.clear();
m_t_truthHitRatio.clear();
m_t_prodVtx_x.clear();
m_t_prodVtx_y.clear();
@@ -1466,11 +1537,25 @@ NtupleDumperAlg::clearTree() const
m_truthd1_y.clear();
m_truthd1_z.clear();
+ m_vertex_x = NaN;
+ m_vertex_y = NaN;
+ m_vertex_z = NaN;
+ m_vertex_chi2 = NaN;
+ m_vertex_px0 = NaN;
+ m_vertex_py0 = NaN;
+ m_vertex_pz0 = NaN;
+ m_vertex_p0 = NaN;
+ m_vertex_px1 = NaN;
+ m_vertex_py1 = NaN;
+ m_vertex_pz1 = NaN;
+ m_vertex_p1 = NaN;
}
void NtupleDumperAlg::clearTrackTruth() const {
m_t_pdg.push_back(0);
m_t_barcode.push_back(-1);
+ m_t_pdg_parent.push_back(0);
+ m_t_barcode_parent.push_back(-1);
m_t_truthHitRatio.push_back(0);
m_t_prodVtx_x.push_back(999999);
m_t_prodVtx_y.push_back(999999);
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
index 67f3bacf8c8f4a6ef134741637d093fdea60c24b..d94379ca4cdc62741cea0a0edd06089382718436 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
@@ -23,6 +23,7 @@
#include "TrackerSimEvent/FaserSiHitCollection.h"
#include "xAODEventInfo/EventInfo.h"
#include "StoreGate/ReadDecorHandle.h"
+#include "FaserActsVertexing/IVertexingTool.h"
#include <vector>
#include <nlohmann/json.hpp>
@@ -93,6 +94,7 @@ private:
ToolHandle<IFaserActsTrackingGeometryTool> m_trackingGeometryTool {this, "TrackingGeometryTool", "FaserActsTrackingGeometryTool"};
ToolHandle<ITrackTruthMatchingTool> m_trackTruthMatchingTool {this, "TrackTruthMatchingTool", "TrackTruthMatchingTool"};
ToolHandle<IFiducialParticleTool> m_fiducialParticleTool {this, "FiducialParticleTool", "FiducialParticleTool"};
+ ToolHandle<FaserTracking::IVertexingTool> m_vertexingTool { this, "VertexingTool", "FaserTracking::PointOfClosestApproachSearchTool"};
const TrackerDD::SCT_DetectorManager* m_detMgr {nullptr};
@@ -124,6 +126,11 @@ private:
BooleanProperty m_applyGoodRunsList {this, "ApplyGoodRunsList", false, "Only write out events passing GRL (data only)"};
StringProperty m_goodRunsList {this, "GoodRunsList", "", "GoodRunsList in json format"};
+ // track quality cuts
+ UnsignedIntegerProperty m_minLayers{this, "minLayers", 7, "Miminimum number of layers of a track."};
+ UnsignedIntegerProperty m_minHits{this, "minHits", 12, "Miminimum number of hits of a track."};
+ DoubleProperty m_maxChi2NDOF{this, "maxChi2NDOF", 25, "Maximum chi2 per degree of freedom."};
+
// json object to hold data read from GRL file (or empty if not)
nlohmann::json m_grl;
@@ -274,6 +281,8 @@ private:
mutable std::vector<int> m_t_pdg; // pdg code of the truth matched particle
mutable std::vector<int> m_t_barcode; // barcode of the truth matched particle
+ mutable std::vector<int> m_t_pdg_parent; // pdg code of the parent of the truth matched particle
+ mutable std::vector<int> m_t_barcode_parent; // barcode of the parent of the truth matched particle
mutable std::vector<double> m_t_truthHitRatio; // ratio of hits on track matched to the truth particle over all hits on track
mutable std::vector<double> m_t_prodVtx_x; // x component of the production vertex in mm
mutable std::vector<double> m_t_prodVtx_y; // y component of the production vertex in mm
@@ -355,6 +364,19 @@ private:
mutable int m_eventsPassed = 0;
mutable int m_eventsFailedGRL = 0;
+
+ mutable double m_vertex_x; // components of reconstructed vertex in mm
+ mutable double m_vertex_y;
+ mutable double m_vertex_z;
+ mutable double m_vertex_chi2; // chi2 value of vertex fit
+ mutable double m_vertex_px0; // components of lowest momentum track extrapolated to vertex in MeV
+ mutable double m_vertex_py0;
+ mutable double m_vertex_pz0;
+ mutable double m_vertex_p0;
+ mutable double m_vertex_px1; // components of second lowest momentum track extrapolated to vertex in MeV
+ mutable double m_vertex_py1;
+ mutable double m_vertex_pz1;
+ mutable double m_vertex_p1;
};
inline const ServiceHandle <ITHistSvc> &NtupleDumperAlg::histSvc() const {
diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsExtrapolationTool.cxx b/Tracking/Acts/FaserActsGeometry/src/FaserActsExtrapolationTool.cxx
index d31fe8d3c11e804e115aef2cc6dfc08a2354b9e7..14d2699822166993a61b0de47a74c1fc5b139f65 100644
--- a/Tracking/Acts/FaserActsGeometry/src/FaserActsExtrapolationTool.cxx
+++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsExtrapolationTool.cxx
@@ -71,7 +71,7 @@ FaserActsExtrapolationTool::initialize()
ATH_MSG_INFO("Initializing ACTS extrapolation");
- m_logger = Acts::getDefaultLogger("ExtrapolationTool", Acts::Logging::INFO);
+ m_logger = Acts::getDefaultLogger("ExtrapolationTool", Acts::Logging::ERROR);
// m_logger = makeActsAthenaLogger(this, "Prop", "FaserActsExtrapTool");
ATH_CHECK( m_trackingGeometryTool.retrieve() );
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/CMakeLists.txt b/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..61908a643b1d9f03eaa3b4ebf9c6eb2d34786a6a
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/CMakeLists.txt
@@ -0,0 +1,44 @@
+atlas_subdir(FaserActsVertexing)
+
+find_package( Eigen )
+find_package( Acts COMPONENTS Core )
+
+atlas_add_library( FaserActsmanVertexingLib
+ PUBLIC_HEADERS FaserActsVertexing
+ INCLUDE_DIRS ${EIGEN_INCLUDE_DIRS}
+ LINK_LIBRARIES ${EIGEN_LIBRARIES}
+ ActsCore
+ AthenaBaseComps
+ FaserActsGeometryLib
+ FaserActsKalmanFilterLib
+ MagFieldConditions
+ StoreGateLib
+ TrackerPrepRawData
+ TrackerRIO_OnTrack
+ TrackerSimEvent
+ TrackerSpacePoint
+ TrkTrack
+ xAODTruth
+)
+
+atlas_add_component(
+ FaserActsVertexing
+ FaserActsVertexing/IVertexingTool.h
+ src/PointOfClosestApproachSearchTool.h
+ src/PointOfClosestApproachSearchTool.cxx
+ src/components/FaserActsVertexing_entries.cxx
+ INCLUDE_DIRS ${EIGEN_INCLUDE_DIRS}
+ LINK_LIBRARIES ${EIGEN_LIBRARIES}
+ ActsCore
+ AthenaBaseComps
+ FaserActsGeometryLib
+ FaserActsKalmanFilterLib
+ MagFieldConditions
+ StoreGateLib
+ TrackerPrepRawData
+ TrackerRIO_OnTrack
+ TrackerSimEvent
+ TrackerSpacePoint
+ TrkTrack
+ xAODTruth
+)
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/FaserActsVertexing/IVertexingTool.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/FaserActsVertexing/IVertexingTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..3c6f93714b6cf07f2428a80c157318399d7b632e
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/FaserActsVertexing/IVertexingTool.h
@@ -0,0 +1,37 @@
+#ifndef FASERACTSVERTEXING_IVERTEXINGTOOL_H
+#define FASERACTSVERTEXING_IVERTEXINGTOOL_H
+
+#include "Acts/EventData/TrackParameters.hpp"
+#include "GaudiKernel/IAlgTool.h"
+#include "GaudiKernel/IInterface.h"
+#include "TrkTrack/Track.h"
+#include <vector>
+
+namespace FaserTracking {
+
+// Point of closest approach between n tracks. Contains the three-dimensional
+// reconstructed position, the chi2 value of the vertex fit and the tracks that
+// originate from this vertex.
+struct POCA {
+ POCA(Eigen::Vector3d vertex, double chi2, std::vector<const Trk::Track *> tracks)
+ : position(vertex), chi2(chi2), tracks(tracks) {}
+ Eigen::Vector3d position;
+ double chi2;
+ std::vector<const Trk::Track *> tracks;
+};
+
+class IVertexingTool : virtual public IAlgTool {
+public:
+ DeclareInterfaceID(IVertexingTool, 1, 0);
+
+ // Run vertex fit and return the fitted POCA.
+ virtual std::optional<POCA> getVertex(const std::vector<const Trk::Track *> &tracks) const = 0;
+
+ // Get bound track parameters at given z position.
+ virtual std::unique_ptr<const Acts::BoundTrackParameters>
+ extrapolateTrack(const Trk::Track *track, double targetPosition) const = 0;
+};
+
+} // namespace FaserTracking
+
+#endif // FASERACTSVERTEXING_IVERTEXINGTOOL_H
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/src/PointOfClosestApproachSearchTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/src/PointOfClosestApproachSearchTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..ea423d703372dbfcc954dc2089b8169282cdc634
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/src/PointOfClosestApproachSearchTool.cxx
@@ -0,0 +1,212 @@
+#include "PointOfClosestApproachSearchTool.h"
+#include "TrkParameters/TrackParameters.h"
+#include <climits>
+
+namespace FaserTracking {
+
+PointOfClosestApproachSearchTool::PointOfClosestApproachSearchTool(const std::string &type,
+ const std::string &name,
+ const IInterface *parent)
+ : base_class(type, name, parent) {}
+
+StatusCode PointOfClosestApproachSearchTool::initialize() {
+ ATH_CHECK(m_trackingGeometryTool.retrieve());
+ ATH_CHECK(m_extrapolationTool.retrieve());
+ return StatusCode::SUCCESS;
+}
+
+StatusCode PointOfClosestApproachSearchTool::finalize() { return StatusCode::SUCCESS; }
+
+std::optional<POCA>
+PointOfClosestApproachSearchTool::getVertex(const std::vector<const Trk::Track *> &tracks) const {
+
+ if (tracks.size() < 2) {
+ ATH_MSG_DEBUG("Found " << tracks.size() << " tracks, but require at least 2 tracks.");
+ return {};
+ } else if (tracks.size() > 3) {
+ ATH_MSG_DEBUG("Found " << tracks.size() << " tracks, but require at most 3 tracks.");
+ return {};
+ }
+
+ clear();
+
+ m_ctx = Gaudi::Hive::currentContext();
+ m_gctx = m_trackingGeometryTool->getNominalGeometryContext().context();
+
+ // convert to Acts::BoundTrackParameters
+ std::vector<Acts::BoundTrackParameters> trackParameters{};
+ for (const Trk::Track *track : tracks) {
+ trackParameters.push_back(getParametersFromTrack(track->trackParameters()->front()));
+ }
+
+ // extrapolate to nSteps z positions betweeen zMin and zMax and calculate center of tracks and
+ // chi2 for each position.
+ size_t nSteps = (m_zMax - m_zMin) / m_stepSize + 1;
+ m_allVertexCandidates.reserve(nSteps);
+ for (int z = m_zMin; z <= m_zMax; z += m_stepSize) {
+ auto vertexCandidate = getVertexCandidate(trackParameters, z);
+ if (vertexCandidate) {
+ m_allVertexCandidates.push_back(*vertexCandidate);
+ }
+ }
+
+ // find vertex candidates with minium (local) chi2 value
+ for (size_t i = 1; i < nSteps - 1; ++i) {
+ if ((m_allVertexCandidates[i - 1].chi2 > m_allVertexCandidates[i].chi2) and
+ (m_allVertexCandidates[i + 1].chi2 > m_allVertexCandidates[i].chi2)) {
+ m_goodVertexCandidates.push_back(m_allVertexCandidates[i]);
+ }
+ }
+
+ // Check if there are one or two local minima and calculate mean in the case of two minima.
+ if (m_goodVertexCandidates.size() == 0) {
+ ATH_MSG_DEBUG("Cannot find minimum. Extrapolation failed. This is likely an error.");
+ return {};
+ } else if (m_goodVertexCandidates.size() == 1) {
+ const VertexCandidate &vx = m_goodVertexCandidates.front();
+ return POCA(vx.position, vx.chi2, tracks);
+ } else if (m_goodVertexCandidates.size() == 2) {
+ const VertexCandidate &vx0 = m_goodVertexCandidates.at(0);
+ const VertexCandidate &vx1 = m_goodVertexCandidates.at(1);
+ return POCA(0.5 * (vx0.position + vx1.position), 0.5 * (vx0.chi2 + vx1.chi2), tracks);
+ } else {
+ ATH_MSG_DEBUG("Cannot find global minimum. Expect a maximum of two local chi2 minima but found "
+ << m_goodVertexCandidates.size() << " minima.");
+ return {};
+ }
+}
+
+// Calculate center of tracks
+Eigen::Vector2d PointOfClosestApproachSearchTool::getCenter(
+ const std::vector<PointOfClosestApproachSearchTool::TrackPosition> &trackPositions) const {
+ Eigen::Vector2d center{0, 0};
+ Eigen::Vector2d sumInverse{0, 0};
+ for (const TrackPosition &value : trackPositions) {
+ center(0) += value.inv(0, 0) * value.pos(0);
+ sumInverse(0) += value.inv(0, 0);
+ center(1) += value.inv(1, 1) * value.pos(1);
+ sumInverse(1) += value.inv(1, 1);
+ }
+ center(0) /= sumInverse(0);
+ center(1) /= sumInverse(1);
+ return center;
+}
+
+// Extrapolate track to given z positon.
+std::unique_ptr<const Acts::BoundTrackParameters>
+PointOfClosestApproachSearchTool::extrapolateTrack(const Trk::Track *track,
+ double targetPosition) const {
+ Acts::BoundTrackParameters startParameters =
+ getParametersFromTrack(track->trackParameters()->front());
+ auto targetSurface = Acts::Surface::makeShared<Acts::PlaneSurface>(
+ Acts::Vector3(0, 0, targetPosition), Acts::Vector3(0, 0, 1));
+ Acts::NavigationDirection navigationDirection =
+ targetPosition > startParameters.referenceSurface().center(m_gctx).z() ? Acts::forward
+ : Acts::backward;
+ return m_extrapolationTool->propagate(m_ctx, startParameters, *targetSurface,
+ navigationDirection);
+}
+
+std::optional<PointOfClosestApproachSearchTool::TrackPosition>
+PointOfClosestApproachSearchTool::extrapolateTrackParameters(
+ const Acts::BoundTrackParameters &startParameters, double z) const {
+ auto targetSurface =
+ Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, z), Acts::Vector3(0, 0, 1));
+ Acts::NavigationDirection navigationDirection =
+ z > startParameters.referenceSurface().center(m_gctx).z() ? Acts::forward : Acts::backward;
+ std::unique_ptr<const Acts::BoundTrackParameters> targetParameters =
+ m_extrapolationTool->propagate(m_ctx, startParameters, *targetSurface, navigationDirection);
+ if (targetParameters != nullptr && targetParameters->covariance().has_value()) {
+ // extrapolate covariance matrix?
+ // return TrackPosition{targetParameters->position(m_gctx),
+ // startParameters.covariance().value().topLeftCorner(2, 2)};
+ return TrackPosition(targetParameters->position(m_gctx),
+ targetParameters->covariance().value().topLeftCorner(2, 2));
+ } else {
+ ATH_MSG_DEBUG("Extrapolation failed.");
+ return std::nullopt;
+ }
+}
+
+std::optional<PointOfClosestApproachSearchTool::VertexCandidate>
+PointOfClosestApproachSearchTool::getVertexCandidate(
+ std::vector<Acts::BoundTrackParameters> &trackParameters, double z) const {
+ // extrapolate tracks to given z position
+ std::vector<TrackPosition> trackPositions{};
+ for (const Acts::BoundTrackParameters &startParameters : trackParameters) {
+ auto tp = extrapolateTrackParameters(startParameters, z);
+ if (tp) {
+ trackPositions.push_back(*tp);
+ }
+ }
+ if (trackPositions.size() < 2) {
+ ATH_MSG_DEBUG(
+ "Extrapolation failed for too many tracks. Require atleast 2 extrapolated positions.");
+ return std::nullopt;
+ }
+
+ // calculate center weighted by inverse covariance matrix
+ Eigen::Vector2d trackCenter = getCenter(trackPositions);
+
+ // either use weighted x and y position or only y position to calculate chi2 value
+ double vertexChi2 = chi2Y(trackPositions, trackCenter);
+
+ return PointOfClosestApproachSearchTool::VertexCandidate(z, vertexChi2, trackCenter);
+}
+
+double PointOfClosestApproachSearchTool::chi2(const std::vector<TrackPosition> &trackPositions,
+ const Eigen::Vector2d &trackCenter) const {
+ double chi2 = 0;
+ for (const TrackPosition &value : trackPositions) {
+ Eigen::Vector2d delta = trackCenter - value.pos;
+ chi2 += delta.dot(value.inv.diagonal().cwiseProduct(delta));
+ }
+ return chi2;
+}
+
+double PointOfClosestApproachSearchTool::chi2Y(const std::vector<TrackPosition> &trackPositions,
+ const Eigen::Vector2d &trackCenter) const {
+ double chi2 = 0;
+ for (const TrackPosition &value : trackPositions) {
+ Eigen::Vector2d delta = trackCenter - value.pos;
+ chi2 += delta.y() * value.inv(1, 1) * delta.y();
+ }
+ return chi2;
+}
+
+Acts::BoundTrackParameters PointOfClosestApproachSearchTool::getParametersFromTrack(
+ const Trk::TrackParameters *trackParameters) const {
+ using namespace Acts::UnitLiterals;
+
+ Acts::GeometryContext geometryContext =
+ m_trackingGeometryTool->getNominalGeometryContext().context();
+ Acts::AngleAxis3 rotZ(M_PI / 2, Acts::Vector3(0., 0., 1.));
+ Acts::Translation3 translation{0., 0., trackParameters->associatedSurface().center().z()};
+ auto transform = Acts::Transform3(translation * rotZ);
+ auto surface = Acts::Surface::makeShared<Acts::PlaneSurface>(transform);
+
+ auto bound = Acts::detail::transformFreeToBoundParameters(
+ trackParameters->position(), 0, trackParameters->momentum(),
+ trackParameters->charge() / trackParameters->momentum().mag() / 1_MeV, *surface,
+ geometryContext);
+ if (!bound.ok()) {
+ ATH_MSG_WARNING("FreeToBound parameter transformation failed");
+ }
+ // convert covariance matrix to GeV
+ Acts::BoundSymMatrix cov = Acts::BoundSymMatrix::Identity();
+ cov.topLeftCorner(5, 5) = *trackParameters->covariance();
+ for (int i = 0; i < cov.rows(); i++) {
+ cov(i, 4) = cov(i, 4) / 1_MeV;
+ }
+ for (int i = 0; i < cov.cols(); i++) {
+ cov(4, i) = cov(4, i) / 1_MeV;
+ }
+ return Acts::BoundTrackParameters{surface, bound.value(), trackParameters->charge(), cov};
+}
+
+void PointOfClosestApproachSearchTool::clear() const {
+ m_allVertexCandidates.clear();
+ m_goodVertexCandidates.clear();
+}
+
+} // namespace FaserTracking
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/src/PointOfClosestApproachSearchTool.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/src/PointOfClosestApproachSearchTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..9ed794948483ee01fcfb13203971ba7a897657da
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/src/PointOfClosestApproachSearchTool.h
@@ -0,0 +1,84 @@
+#ifndef FASERACTSVERTEXING_POINTOFCLOSESTAPPROACHSEARCHTOOL_H
+#define FASERACTSVERTEXING_POINTOFCLOSESTAPPROACHSEARCHTOOL_H
+
+#include "TrackerPrepRawData/FaserSCT_ClusterContainer.h"
+#include "Acts/EventData/TrackParameters.hpp"
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "FaserActsGeometryInterfaces/IFaserActsExtrapolationTool.h"
+#include "FaserActsGeometryInterfaces/IFaserActsTrackingGeometryTool.h"
+#include "FaserActsVertexing/IVertexingTool.h"
+#include "TrkTrack/Track.h"
+
+namespace FaserTracking {
+
+class PointOfClosestApproachSearchTool : public extends<AthAlgTool, IVertexingTool> {
+public:
+ PointOfClosestApproachSearchTool(const std::string &type, const std::string &name,
+ const IInterface *parent);
+ virtual ~PointOfClosestApproachSearchTool() = default;
+ virtual StatusCode initialize() override;
+ virtual StatusCode finalize() override;
+
+ std::optional<POCA> getVertex(const std::vector<const Trk::Track *> &tracks) const override;
+
+ std::unique_ptr<const Acts::BoundTrackParameters>
+ extrapolateTrack(const Trk::Track *track, double targetPosition) const override;
+
+private:
+ struct VertexCandidate {
+ VertexCandidate(double z, double chi2, const Eigen::Vector2d ¢er) : chi2(chi2) {
+ position = Eigen::Vector3d(center.x(), center.y(), z);
+ }
+ double chi2;
+ Eigen::Vector3d position;
+ };
+
+ struct TrackPosition {
+ TrackPosition(const Acts::Vector3 &position, const Acts::BoundSymMatrix &covariance) {
+ pos << position.x(), position.y();
+ cov << covariance(1, 1), 0, 0, covariance(0, 0);
+ inv = cov.inverse();
+ }
+ Eigen::Vector2d pos;
+ Eigen::Matrix2d cov;
+ Eigen::Matrix2d inv;
+ };
+
+ ToolHandle<IFaserActsTrackingGeometryTool> m_trackingGeometryTool{
+ this, "TrackingGeometryTool", "FaserActsTrackingGeometryTool"};
+ ToolHandle<IFaserActsExtrapolationTool> m_extrapolationTool{
+ this, "ExtrapolationTool", "FaserActsExtrapolationTool"};
+ DoubleProperty m_stepSize{this, "StepSize", 50};
+ DoubleProperty m_zMin{this, "ZMin", -2000};
+ DoubleProperty m_zMax{this, "ZMax", 500};
+ StringProperty m_debugLevel{this, "DebugLevel", "INFO"};
+
+ void clear() const;
+
+ std::optional<TrackPosition>
+ extrapolateTrackParameters(const Acts::BoundTrackParameters &startParameters, double z) const;
+
+ Eigen::Vector2d getCenter(const std::vector<TrackPosition> &trackPositions) const;
+
+ Acts::BoundTrackParameters
+ getParametersFromTrack(const Trk::TrackParameters *trackParameters) const;
+
+ double chi2(const std::vector<TrackPosition> &trackPositions,
+ const Eigen::Vector2d &trackCenter) const;
+
+ double chi2Y(const std::vector<TrackPosition> &trackPositions,
+ const Eigen::Vector2d &trackCenter) const;
+
+ std::optional<VertexCandidate>
+ getVertexCandidate(std::vector<Acts::BoundTrackParameters> &trackParameters, double z) const;
+
+ mutable EventContext m_ctx;
+ mutable Acts::GeometryContext m_gctx;
+
+ mutable std::vector<VertexCandidate> m_allVertexCandidates;
+ mutable std::vector<VertexCandidate> m_goodVertexCandidates;
+};
+
+} // namespace FaserTracking
+
+#endif /* FASERACTSVERTEXING_POINTOFCLOSESTAPPROACHSEARCHTOOL_H */
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/src/components/FaserActsVertexing_entries.cxx b/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/src/components/FaserActsVertexing_entries.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..4f27a5a3d95ba97b48d36ac0b43c634c1d99d9ae
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing/src/components/FaserActsVertexing_entries.cxx
@@ -0,0 +1,3 @@
+#include "../PointOfClosestApproachSearchTool.h"
+
+DECLARE_COMPONENT(FaserTracking::PointOfClosestApproachSearchTool)
\ No newline at end of file