Merge branch 'master-pairvtx' into 'master'

Development of track pair analysis

See merge request faser/calypso!280

Tracker/TrackerRecAlgs/PairVertex/CMakeLists.txt

+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)

Tracker/TrackerRecAlgs/PairVertex/python/PairVertexAlgConfig.py

+"""
+    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

Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.cxx

+#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

Tracker/TrackerRecAlgs/PairVertex/src/PairVertexAlg.h

+/*
+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

Tracker/TrackerRecAlgs/PairVertex/src/component/PairVertex_entries.cxx

+#include "../PairVertexAlg.h"
+
+DECLARE_COMPONENT(Tracker::PairVertexAlg)
\ No newline at end of file