From 05edc63d6e702afb93b7c2600e14667fc6b64227 Mon Sep 17 00:00:00 2001
From: Tobias Boeckh <tobias.boeckh@cern.ch>
Date: Sun, 13 Mar 2022 23:18:02 +0100
Subject: [PATCH] added example for ACTS propagation
---
.../MyExtrapolationExample/CMakeLists.txt | 12 ++
.../python/MyExtrapolationExampleConfig.py | 30 +++++
.../MyExtrapolationExample/python/__init__.py | 0
.../src/MyExtrapolationExample.cxx | 114 ++++++++++++++++++
.../src/MyExtrapolationExample.h | 42 +++++++
.../MyExtrapolationExample_entries.cxx | 3 +
.../test/MyExtrapolationExampleDbg.py | 30 +++++
7 files changed, 231 insertions(+)
create mode 100644 Tracker/TrackerRecAlgs/MyExtrapolationExample/CMakeLists.txt
create mode 100644 Tracker/TrackerRecAlgs/MyExtrapolationExample/python/MyExtrapolationExampleConfig.py
create mode 100644 Tracker/TrackerRecAlgs/MyExtrapolationExample/python/__init__.py
create mode 100644 Tracker/TrackerRecAlgs/MyExtrapolationExample/src/MyExtrapolationExample.cxx
create mode 100644 Tracker/TrackerRecAlgs/MyExtrapolationExample/src/MyExtrapolationExample.h
create mode 100644 Tracker/TrackerRecAlgs/MyExtrapolationExample/src/components/MyExtrapolationExample_entries.cxx
create mode 100644 Tracker/TrackerRecAlgs/MyExtrapolationExample/test/MyExtrapolationExampleDbg.py
diff --git a/Tracker/TrackerRecAlgs/MyExtrapolationExample/CMakeLists.txt b/Tracker/TrackerRecAlgs/MyExtrapolationExample/CMakeLists.txt
new file mode 100644
index 000000000..72f4a7ac2
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/MyExtrapolationExample/CMakeLists.txt
@@ -0,0 +1,12 @@
+atlas_subdir(MyExtrapolationExample)
+
+atlas_add_component(
+ MyExtrapolationExample
+ src/MyExtrapolationExample.h
+ src/MyExtrapolationExample.cxx
+ src/components/MyExtrapolationExample_entries.cxx
+ LINK_LIBRARIES AthenaBaseComps StoreGateLib GeneratorObjects FaserActsGeometryLib TrackerSimEvent TrackerIdentifier TrackerReadoutGeometry
+)
+
+atlas_install_python_modules(python/*.py)
+atlas_install_scripts(test/*.py)
diff --git a/Tracker/TrackerRecAlgs/MyExtrapolationExample/python/MyExtrapolationExampleConfig.py b/Tracker/TrackerRecAlgs/MyExtrapolationExample/python/MyExtrapolationExampleConfig.py
new file mode 100644
index 000000000..6c2738946
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/MyExtrapolationExample/python/MyExtrapolationExampleConfig.py
@@ -0,0 +1,30 @@
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+from FaserSCT_GeoModel.FaserSCT_GeoModelConfig import FaserSCT_GeometryCfg
+from MagFieldServices.MagFieldServicesConfig import MagneticFieldSvcCfg
+
+# def FaserActsTrackingGeometrySvcCfg(flags, **kwargs):
+# acc = ComponentAccumulator()
+# FaserActsTrackingGeometrySvc = CompFactory.FaserActsTrackingGeometrySvc
+# acc.addService(FaserActsTrackingGeometrySvc(name="FaserActsTrackingGeometrySvc", **kwargs))
+# return acc
+
+# def FaserActsAlignmentCondAlgCfg(flags, **kwargs):
+# acc = ComponentAccumulator()
+# acc.addCondAlgo(CompFactory.FaserActsAlignmentCondAlg(name="FaserActsAlignmentCondAlg", **kwargs))
+# return acc
+
+def MyExtrapolationExampleCfg(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")
+
+ my_extrapolation_example = CompFactory.Tracker.MyExtrapolationExample(**kwargs)
+ my_extrapolation_example.ExtrapolationTool = actsExtrapolationTool
+ acc.addEventAlgo(my_extrapolation_example)
+ return acc
diff --git a/Tracker/TrackerRecAlgs/MyExtrapolationExample/python/__init__.py b/Tracker/TrackerRecAlgs/MyExtrapolationExample/python/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/Tracker/TrackerRecAlgs/MyExtrapolationExample/src/MyExtrapolationExample.cxx b/Tracker/TrackerRecAlgs/MyExtrapolationExample/src/MyExtrapolationExample.cxx
new file mode 100644
index 000000000..aecebc372
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/MyExtrapolationExample/src/MyExtrapolationExample.cxx
@@ -0,0 +1,114 @@
+#include "MyExtrapolationExample.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 {
+
+MyExtrapolationExample::MyExtrapolationExample(const std::string &name, ISvcLocator *pSvcLocator)
+ : AthReentrantAlgorithm(name, pSvcLocator) {}
+
+StatusCode MyExtrapolationExample::initialize() {
+ ATH_CHECK(m_mcEventCollectionKey.initialize());
+ ATH_CHECK(m_faserSiHitKey.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 MyExtrapolationExample::execute(const EventContext &ctx) const {
+ const Acts::GeometryContext gctx =
+ m_extrapolationTool->trackingGeometryTool()->getNominalGeometryContext().context();
+
+ std::vector<double> z_positions {};
+ SG::ReadHandle<FaserSiHitCollection> siHitCollection {m_faserSiHitKey, ctx};
+ ATH_CHECK(siHitCollection.isValid());
+ for (const FaserSiHit& hit : *siHitCollection) {
+ if ((hit.getStation() == 1) && (hit.getPlane() == 0) && (hit.getSensor() == 0)) {
+ if (hit.particleLink()) {
+ if (std::abs(hit.particleLink()->pdg_id()) != 13)
+ continue;
+ Identifier id = m_idHelper->wafer_id(hit.getStation(), hit.getPlane(), hit.getRow(), hit.getModule(), hit.getSensor());
+ const TrackerDD::SiDetectorElement* element = m_detMgr->getDetectorElement(id);
+ const HepGeom::Point3D<double> globalStartPosition =
+ Amg::EigenTransformToCLHEP(element->transformHit()) * hit.localStartPosition();
+ z_positions.push_back(globalStartPosition.z());
+ ATH_MSG_DEBUG("SiHit: " << globalStartPosition.x() << ", " << globalStartPosition.y() << ", " << globalStartPosition.z());
+ }
+ }
+ }
+
+ double z_mean = 0;
+ for (double z : z_positions) {
+ z_mean += z;
+ }
+ z_mean /= z_positions.size();
+
+ 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;
+ }
+
+
+ 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 MyExtrapolationExample::finalize() {
+ return StatusCode::SUCCESS;
+}
+
+} // Tracker
diff --git a/Tracker/TrackerRecAlgs/MyExtrapolationExample/src/MyExtrapolationExample.h b/Tracker/TrackerRecAlgs/MyExtrapolationExample/src/MyExtrapolationExample.h
new file mode 100644
index 000000000..cd5a2ed07
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/MyExtrapolationExample/src/MyExtrapolationExample.h
@@ -0,0 +1,42 @@
+/*
+Copyright (C) 2022 CERN for the benefit of the FASER collaboration
+*/
+
+#ifndef MYEXTRAPOLATIONEXAMPLE_MYEXTRAPOLATIONEXAMPLE_H
+#define MYEXTRAPOLATIONEXAMPLE_MYEXTRAPOLATIONEXAMPLE_H
+
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "FaserActsGeometryInterfaces/IFaserActsExtrapolationTool.h"
+#include "GeneratorObjects/McEventCollection.h"
+#include "TrackerSimEvent/FaserSiHitCollection.h"
+
+class FaserSCT_ID;
+namespace TrackerDD {
+class SCT_DetectorManager;
+
+}
+namespace Tracker {
+
+class MyExtrapolationExample : public AthReentrantAlgorithm {
+ public:
+ MyExtrapolationExample(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 <FaserSiHitCollection> m_faserSiHitKey {
+ this, "FaserSiHitCollection", "SCT_Hits"};
+ ToolHandle<IFaserActsExtrapolationTool> m_extrapolationTool {
+ this, "ExtrapolationTool", "FaserActsExtrapolationTool"};
+};
+
+} // namespace Tracker
+
+#endif // MYEXTRAPOLATIONEXAMPLE_MYEXTRAPOLATIONEXAMPLE_H
diff --git a/Tracker/TrackerRecAlgs/MyExtrapolationExample/src/components/MyExtrapolationExample_entries.cxx b/Tracker/TrackerRecAlgs/MyExtrapolationExample/src/components/MyExtrapolationExample_entries.cxx
new file mode 100644
index 000000000..6601e0a3f
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/MyExtrapolationExample/src/components/MyExtrapolationExample_entries.cxx
@@ -0,0 +1,3 @@
+#include "../MyExtrapolationExample.h"
+
+DECLARE_COMPONENT(Tracker::MyExtrapolationExample)
\ No newline at end of file
diff --git a/Tracker/TrackerRecAlgs/MyExtrapolationExample/test/MyExtrapolationExampleDbg.py b/Tracker/TrackerRecAlgs/MyExtrapolationExample/test/MyExtrapolationExampleDbg.py
new file mode 100644
index 000000000..fe83560fa
--- /dev/null
+++ b/Tracker/TrackerRecAlgs/MyExtrapolationExample/test/MyExtrapolationExampleDbg.py
@@ -0,0 +1,30 @@
+#!/usr/bin/env python
+
+import sys
+from AthenaCommon.Logging import log
+from AthenaCommon.Constants import DEBUG, VERBOSE
+from AthenaCommon.Configurable import Configurable
+from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+from MyExtrapolationExample.MyExtrapolationExampleConfig import MyExtrapolationExampleCfg
+
+
+log.setLevel(DEBUG)
+Configurable.configurableRun3Behavior = True
+
+# Configure
+ConfigFlags.Input.Files = ['my.HITS.pool.root']
+ConfigFlags.Output.ESDFileName = "MyExtrapolationExample.ESD.root"
+ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01"
+ConfigFlags.GeoModel.Align.Dynamic = False
+ConfigFlags.Beam.NumberOfCollisions = 0.
+ConfigFlags.lock()
+
+acc = MainServicesCfg(ConfigFlags)
+acc.merge(PoolReadCfg(ConfigFlags))
+acc.merge(MyExtrapolationExampleCfg(ConfigFlags))
+acc.getEventAlgo("Tracker::MyExtrapolationExample").OutputLevel = VERBOSE
+
+sc = acc.run(maxEvents=10)
+sys.exit(not sc.isSuccess())
--
GitLab