diff --git a/Tracker/TrackerRecAlgs/PairVertex/CMakeLists.txt b/Tracker/TrackerRecAlgs/PairVertex/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..85053a1f49bd551ea136c5cd61f6cee40543ed26
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/PairVertex/CMakeLists.txt
@@ -0,0 +1,12 @@
+atlas_subdir(PairVertex)
+
+atlas_add_component(
+ PairVertex
+ src/PairVertexAlg.h
+ src/PairVertexAlg.cxx
+ src/component/PairVertex_entries.cxx
+ LINK_LIBRARIES AthenaBaseComps StoreGateLib GeneratorObjects FaserActsGeometryLib TrackerSimEvent TrackerIdentifier TrackerReadoutGeometry TrkTrack GeoPrimitives
+)
+
+atlas_install_python_modules(python/*.py)
+# atlas_install_scripts(test/*.py)
diff --git a/Tracker/TrackerRecAlgs/PairVertex/python/PairVertexAlgConfig.py b/Tracker/TrackerRecAlgs/PairVertex/python/PairVertexAlgConfig.py
new file mode 100644
index 0000000000000000000000000000000000000000..f98af1629a39bfc590ef3ab0929da4b849ad5ce6
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/PairVertex/python/PairVertexAlgConfig.py
@@ -0,0 +1,90 @@
+"""
+ Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
+"""
+
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+from FaserSCT_GeoModel.FaserSCT_GeoModelConfig import FaserSCT_GeometryCfg
+from MagFieldServices.MagFieldServicesConfig import MagneticFieldSvcCfg
+
+def PairVertexAlgCfg(flags, **kwargs):
+ acc = FaserSCT_GeometryCfg(flags)
+ acc.merge(MagneticFieldSvcCfg(flags))
+ # acc.merge(FaserActsTrackingGeometrySvcCfg(flags))
+ # acc.merge(FaserActsAlignmentCondAlgCfg(flags))
+
+ actsExtrapolationTool = CompFactory.FaserActsExtrapolationTool("FaserActsExtrapolationTool")
+ actsExtrapolationTool.MaxSteps = 1000
+ actsExtrapolationTool.TrackingGeometryTool = CompFactory.FaserActsTrackingGeometryTool("TrackingGeometryTool")
+
+ PairVertexAlg = CompFactory.Tracker.PairVertexAlg("PairVertexAlg",**kwargs)
+ PairVertexAlg.ExtrapolationTool = actsExtrapolationTool
+ acc.addEventAlgo(PairVertexAlg)
+ return acc
+
+if __name__ == "__main__":
+
+ import sys
+ from AthenaCommon.Logging import log, logging
+ from AthenaCommon.Constants import DEBUG, VERBOSE, INFO
+ from AthenaCommon.Configurable import Configurable
+ from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+ from AthenaConfiguration.TestDefaults import defaultTestFiles
+ from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+ from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+ from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+
+ # Set up logging and new style config
+ log.setLevel(DEBUG)
+ Configurable.configurableRun3Behavior = True
+
+ # Configure
+ ConfigFlags.Input.Files = [
+ '/eos/experiment/faser/sim/mdc/foresee/110004/rec/r0007/FaserMC-MDC_FS_Amm_316MeV_2Em6-110004-00000-s0003-r0007-xAOD.root'
+ ]
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02" # Always needed; must match FaserVersionS
+ ConfigFlags.IOVDb.DatabaseInstance = "OFLP200" # Use MC conditions for now
+ ConfigFlags.Input.ProjectName = "data21" # Needed to bypass autoconfig
+ ConfigFlags.Input.isMC = True # Needed to bypass autoconfig
+ ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # FASER geometry
+ ConfigFlags.Common.isOnline = False
+ ConfigFlags.GeoModel.Align.Dynamic = False
+ ConfigFlags.Beam.NumberOfCollisions = 0.
+
+ ConfigFlags.Detector.GeometryFaserSCT = True
+
+ ConfigFlags.lock()
+
+ # Core components
+ acc = MainServicesCfg(ConfigFlags)
+ acc.merge(PoolReadCfg(ConfigFlags))
+
+ # algorithm
+ acc.merge(PairVertexAlgCfg(ConfigFlags))
+
+ # # Hack to avoid problem with our use of MC databases when isMC = False
+ # replicaSvc = acc.getService("DBReplicaSvc")
+ # replicaSvc.COOLSQLiteVetoPattern = ""
+ # replicaSvc.UseCOOLSQLite = True
+ # replicaSvc.UseCOOLFrontier = False
+ # replicaSvc.UseGeomSQLite = True
+
+ # Timing
+ #acc.merge(MergeRecoTimingObjCfg(ConfigFlags))
+
+ # Dump config
+ # logging.getLogger('forcomps').setLevel(VERBOSE)
+ # acc.foreach_component("*").OutputLevel = VERBOSE
+ # acc.foreach_component("*ClassID*").OutputLevel = INFO
+ # acc.getCondAlgo("FaserSCT_AlignCondAlg").OutputLevel = VERBOSE
+ # acc.getCondAlgo("FaserSCT_DetectorElementCondAlg").OutputLevel = VERBOSE
+ # acc.getService("StoreGateSvc").Dump = True
+ # acc.getService("ConditionStore").Dump = True
+ # acc.printConfig(withDetails=True)
+ # ConfigFlags.dump()
+
+ # Execute and finish
+ sc = acc.run(maxEvents=-1)
+
+ # Success should be 0
+ sys.exit(not sc.isSuccess())
\ No newline at end of file
diff --git a/Tracker/TrackerRecAlgs/PairVertex/python/__init__.py b/Tracker/TrackerRecAlgs/PairVertex/python/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.cxx b/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..99058166c7252ddbff852b0a09ec8d49a15e54fb
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.cxx
@@ -0,0 +1,163 @@
+#include "PairVertexAlg.h"
+
+#include "FaserActsGeometryInterfaces/IFaserActsTrackingGeometryTool.h"
+#include "GeoPrimitives/CLHEPtoEigenConverter.h"
+#include "Acts/Surfaces/PerigeeSurface.hpp"
+#include "TrackerIdentifier/FaserSCT_ID.h"
+#include "TrackerReadoutGeometry/SCT_DetectorManager.h"
+#include "TrackerReadoutGeometry/SiDetectorElement.h"
+#include <cmath>
+
+
+namespace Tracker {
+
+PairVertexAlg::PairVertexAlg(const std::string &name, ISvcLocator *pSvcLocator)
+ : AthReentrantAlgorithm(name, pSvcLocator) {}
+
+StatusCode PairVertexAlg::initialize()
+{
+ ATH_CHECK(m_mcEventCollectionKey.initialize());
+ ATH_CHECK(m_trackCollectionKey.initialize());
+ ATH_CHECK(m_extrapolationTool.retrieve());
+ ATH_CHECK(detStore()->retrieve(m_idHelper, "FaserSCT_ID"));
+ ATH_CHECK(detStore()->retrieve(m_detMgr, "SCT"));
+ return StatusCode::SUCCESS;
+}
+
+StatusCode PairVertexAlg::execute(const EventContext &ctx) const
+{
+ m_numberOfEvents++;
+ const Acts::GeometryContext gctx =
+ m_extrapolationTool->trackingGeometryTool()->getNominalGeometryContext().context();
+
+ std::vector<double> z_positions {};
+
+ SG::ReadHandle<McEventCollection> mcEvents {m_mcEventCollectionKey, ctx};
+ ATH_CHECK(mcEvents.isValid());
+ if (mcEvents->size() != 1)
+ {
+ ATH_MSG_ERROR("There should be exactly one event in the McEventCollection.");
+ return StatusCode::FAILURE;
+ }
+
+ SG::ReadHandle<TrackCollection> tracks { m_trackCollectionKey, ctx };
+ ATH_CHECK(tracks.isValid());
+
+ const Trk::TrackParameters* positive {nullptr};
+ const Trk::TrackParameters* negative {nullptr};
+
+ for (auto trk : *tracks)
+ {
+ const Trk::TrackParameters* upstream {nullptr};
+ if (trk == nullptr)
+ {
+ ATH_MSG_WARNING("Null pointer from TrackContainer.");
+ m_invalidTrackContainerEntries++;
+ continue;
+ }
+ auto parameters = trk->trackParameters();
+ if (parameters == nullptr || parameters->empty())
+ {
+ m_numberOfNoParameterTracks++;
+ ATH_MSG_WARNING("Track without parameters found.");
+ return StatusCode::SUCCESS;
+ }
+ for (auto state : *parameters)
+ {
+ if (state == nullptr)
+ {
+ m_numberOfNullParameters++;
+ ATH_MSG_WARNING("Null track parameters returned.");
+ continue;
+ }
+ if (!state->hasSurface())
+ {
+ m_numberOfParametersWithoutSurface++;
+ ATH_MSG_WARNING("Track state without surface found.");
+ continue;
+ }
+ if (!upstream || upstream->associatedSurface().center().z() > state->associatedSurface().center().z())
+ upstream = state;
+ }
+ if (!upstream)
+ {
+ m_numberOfNoUpstreamTracks++;
+ ATH_MSG_WARNING("Unable to find track parameters on any surface");
+ continue;
+ }
+ auto momentum = upstream->momentum();
+ double charge = upstream->charge();
+ if (charge > 0 || momentum.mag() > positive->momentum().mag())
+ {
+ positive = upstream;
+ }
+ else if (charge < 0 || momentum.mag() > negative->momentum().mag())
+ {
+ negative = upstream;
+ }
+ }
+
+ if (positive != nullptr) m_numberOfPositiveTracks++;
+ if (negative != nullptr) m_numberOfNegativeTracks++;
+
+ if (positive != nullptr && negative != nullptr) m_numberOfOppositeSignPairs++;
+
+ // for (const HepMC::GenParticle* particle : mcEvents->front()->particle_range())
+ // {
+ // if ((std::abs(particle->pdg_id()) != 13)) continue;
+ // const HepMC::FourVector& vertex = particle->production_vertex()->position();
+ // if (vertex.z() > 0) continue;
+ // const HepMC::FourVector& momentum = particle->momentum();
+ // double phi = momentum.phi();
+ // double theta = momentum.theta();
+ // double charge = particle->pdg_id() > 0 ? -1 : 1;
+ // double abs_momentum = momentum.rho() * m_MeV2GeV;
+ // double qop = charge / abs_momentum;
+ // // The coordinate system of the Acts::PlaneSurface is defined as
+ // // T = Z = normal, U = X x T = -Y, V = T x U = x
+ // // Acts::BoundVector pars;
+ // // pars << -vertex.y(), vertex.x(), phi, theta, qop, vertex.t();
+ // Acts::BoundVector pars = Acts::BoundVector::Zero();
+ // pars[Acts::eBoundLoc0] = -vertex.y();
+ // pars[Acts::eBoundLoc1] = vertex.x();
+ // pars[Acts::eBoundPhi] = phi;
+ // pars[Acts::eBoundTheta] = theta;
+ // pars[Acts::eBoundQOverP] = qop;
+ // pars[Acts::eBoundTime] = vertex.t();
+
+ // auto startSurface = Acts::Surface::makeShared<Acts::PlaneSurface>(
+ // Acts::Vector3(0, 0, vertex.z()), Acts::Vector3(0, 0, 1));
+ // auto targetSurface = Acts::Surface::makeShared<Acts::PlaneSurface>(
+ // Acts::Vector3(0, 0, z_mean), Acts::Vector3(0, 0, 1));
+ // Acts::BoundTrackParameters startParameters(
+ // std::move(startSurface), pars, charge);
+ // ATH_MSG_DEBUG("vertex: " << vertex.x() << ", " << vertex.y() << ", " << vertex.z());
+ // ATH_MSG_DEBUG("vertex momentum: " << momentum.x() * m_MeV2GeV << ", " << momentum.y() * m_MeV2GeV << ", " << momentum.z() * m_MeV2GeV);
+ // std::unique_ptr<const Acts::BoundTrackParameters> targetParameters =
+ // m_extrapolationTool->propagate(ctx, startParameters, *targetSurface);
+ // if (targetParameters)
+ // {
+ // Acts::Vector3 targetPosition = targetParameters->position(gctx);
+ // Acts::Vector3 targetMomentum = targetParameters->momentum();
+ // ATH_MSG_DEBUG("vertex: " << vertex.x() << ", " << vertex.y() << ", " << vertex.z());
+ // ATH_MSG_DEBUG("origin: " << targetPosition.x() << ", " << targetPosition.y() << ", " << targetPosition.z());
+ // ATH_MSG_DEBUG("vertex momentum: " << momentum.x() * m_MeV2GeV << ", " << momentum.y() * m_MeV2GeV << ", " << momentum.z() * m_MeV2GeV);
+ // ATH_MSG_DEBUG("origin momentum: " << targetMomentum.x() << ", " << targetMomentum.y() << ", " << targetMomentum.z());
+ // }
+ // }
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode PairVertexAlg::finalize()
+{
+ ATH_MSG_INFO("Found " << m_numberOfOppositeSignPairs << " opposite sign pairs in " << m_numberOfEvents << " total events.");
+ ATH_MSG_INFO(m_numberOfPositiveTracks << " events had a positive track, and " << m_numberOfNegativeTracks << " had a negative track.");
+ ATH_MSG_INFO(m_numberOfNoUpstreamTracks << " tracks could not locate an upstream position.");
+ ATH_MSG_INFO(m_numberOfNoParameterTracks << " tracks in track pairs had no track parameters.");
+ ATH_MSG_INFO(m_numberOfNullParameters << " track parameters were null and " << m_numberOfParametersWithoutSurface << " had no associated surface.");
+ ATH_MSG_INFO(m_invalidTrackContainerEntries << " invalid TrackContainer entries found.");
+ return StatusCode::SUCCESS;
+}
+
+} // Tracker
diff --git a/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.h b/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..d15fddf6e288828d143218fa6d2f686a50076ba0
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.h
@@ -0,0 +1,51 @@
+/*
+Copyright (C) 2022 CERN for the benefit of the FASER collaboration
+*/
+
+#pragma once
+
+#include "GeoPrimitives/GeoPrimitives.h"
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "FaserActsGeometryInterfaces/IFaserActsExtrapolationTool.h"
+#include "GeneratorObjects/McEventCollection.h"
+#include "TrkTrack/TrackCollection.h"
+
+class FaserSCT_ID;
+namespace TrackerDD
+{
+ class SCT_DetectorManager;
+}
+
+namespace Tracker
+{
+ class PairVertexAlg : public AthReentrantAlgorithm
+ {
+ public:
+ PairVertexAlg(const std::string& name, ISvcLocator* pSvcLocator);
+
+ virtual StatusCode initialize() override;
+ virtual StatusCode execute(const EventContext& ctx) const override;
+ virtual StatusCode finalize() override;
+
+ private:
+ double m_MeV2GeV = 1e-3;
+ const FaserSCT_ID* m_idHelper {nullptr};
+ const TrackerDD::SCT_DetectorManager* m_detMgr {nullptr};
+ SG::ReadHandleKey<McEventCollection> m_mcEventCollectionKey { this, "McEventCollection", "TruthEvent" };
+
+ SG::ReadHandleKey<TrackCollection> m_trackCollectionKey { this, "TrackCollection", "CKFTrackCollection" };
+
+ ToolHandle<IFaserActsExtrapolationTool> m_extrapolationTool { this, "ExtrapolationTool", "FaserActsExtrapolationTool" };
+
+ mutable std::atomic<size_t> m_numberOfEvents{0};
+ mutable std::atomic<size_t> m_numberOfPositiveTracks{0};
+ mutable std::atomic<size_t> m_numberOfNegativeTracks{0};
+ mutable std::atomic<size_t> m_numberOfOppositeSignPairs{0};
+ mutable std::atomic<size_t> m_numberOfNoParameterTracks{0};
+ mutable std::atomic<size_t> m_invalidTrackContainerEntries{0};
+ mutable std::atomic<size_t> m_numberOfNullParameters{0};
+ mutable std::atomic<size_t> m_numberOfParametersWithoutSurface{0};
+ mutable std::atomic<size_t> m_numberOfNoUpstreamTracks{0};
+ };
+
+}
\ No newline at end of file
diff --git a/Tracker/TrackerRecAlgs/PairVertex/src/component/PairVertex_entries.cxx b/Tracker/TrackerRecAlgs/PairVertex/src/component/PairVertex_entries.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..db08d4e764775c0a28da9e53538e2e5daae46932
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/PairVertex/src/component/PairVertex_entries.cxx
@@ -0,0 +1,3 @@
+#include "../PairVertexAlg.h"
+
+DECLARE_COMPONENT(Tracker::PairVertexAlg)
\ No newline at end of file