diff --git a/Tracking/Acts/ActsGeometry/ActsGeometry/ActsExtrapolationTool.h b/Tracking/Acts/ActsGeometry/ActsGeometry/ActsExtrapolationTool.h
index e794e2c9fc55c4ed57ed5d5ea1c232cbafab7134..1b1b5e3a82b276b96bb451b7458bf3372cb8489a 100644
--- a/Tracking/Acts/ActsGeometry/ActsGeometry/ActsExtrapolationTool.h
+++ b/Tracking/Acts/ActsGeometry/ActsGeometry/ActsExtrapolationTool.h
@@ -68,8 +68,32 @@ private:
public:
virtual
std::vector<Acts::detail::Step>
+ propagationSteps(const EventContext& ctx,
+ const Acts::BoundParameters& startParameters,
+ Acts::NavigationDirection navDir = Acts::forward,
+ double pathLimit = std::numeric_limits<double>::max()) const override;
+
+ virtual
+ std::unique_ptr<const Acts::CurvilinearParameters>
+ propagate(const EventContext& ctx,
+ const Acts::BoundParameters& startParameters,
+ Acts::NavigationDirection navDir = Acts::forward,
+ double pathLimit = std::numeric_limits<double>::max()) const override;
+
+ virtual
+ std::vector<Acts::detail::Step>
+ propagationSteps(const EventContext& ctx,
+ const Acts::BoundParameters& startParameters,
+ const Acts::Surface& target,
+ Acts::NavigationDirection navDir = Acts::forward,
+ double pathLimit = std::numeric_limits<double>::max()) const override;
+
+ virtual
+ std::unique_ptr<const Acts::BoundParameters>
propagate(const EventContext& ctx,
const Acts::BoundParameters& startParameters,
+ const Acts::Surface& target,
+ Acts::NavigationDirection navDir = Acts::forward,
double pathLimit = std::numeric_limits<double>::max()) const override;
virtual
diff --git a/Tracking/Acts/ActsGeometry/CMakeLists.txt b/Tracking/Acts/ActsGeometry/CMakeLists.txt
index 44bf468ff03d9a5826419aa8d7acbe1cc793efc5..9bd53b33422c69c7f0fd27bc6a88ab12978c06e4 100755
--- a/Tracking/Acts/ActsGeometry/CMakeLists.txt
+++ b/Tracking/Acts/ActsGeometry/CMakeLists.txt
@@ -9,9 +9,9 @@ atlas_depends_on_subdirs( PUBLIC
DetectorDescription/Identifier
InnerDetector/InDetDetDescr/InDetIdentifier
InnerDetector/InDetDetDescr/InDetReadoutGeometry
- InnerDetector/InDetDetDescr/PixelReadoutGeometry
- InnerDetector/InDetDetDescr/SCT_ReadoutGeometry
- InnerDetector/InDetDetDescr/TRT_ReadoutGeometry
+ InnerDetector/InDetDetDescr/PixelReadoutGeometry
+ InnerDetector/InDetDetDescr/SCT_ReadoutGeometry
+ InnerDetector/InDetDetDescr/TRT_ReadoutGeometry
Control/AthenaBaseComps
AthenaKernel
DetectorDescription/GeoModel/GeoModelUtilities
@@ -46,9 +46,9 @@ atlas_add_library( ActsGeometryLib
ActsInteropLib
ActsGeometryInterfacesLib
ActsCore
- PixelReadoutGeometry
- SCT_ReadoutGeometry
- TRT_ReadoutGeometry)
+ PixelReadoutGeometry
+ SCT_ReadoutGeometry
+ TRT_ReadoutGeometry)
atlas_add_component( ActsGeometry
src/ActsExtrapolationAlg.cxx
diff --git a/Tracking/Acts/ActsGeometry/python/ActsGeometryConfig.py b/Tracking/Acts/ActsGeometry/python/ActsGeometryConfig.py
index d838f1599ed485bed7a878e335ac6a3631bddabc..2a07d14ba361401689ca8c9d504f86a2ed974174 100644
--- a/Tracking/Acts/ActsGeometry/python/ActsGeometryConfig.py
+++ b/Tracking/Acts/ActsGeometry/python/ActsGeometryConfig.py
@@ -1,9 +1,61 @@
-from AthenaCommon.Logging import logging
-logging.getLogger().info("Importing %s", __name__)
+# Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
-from AthenaCommon.AthenaCommonFlags import athenaCommonFlags
-from AthenaCommon.GlobalFlags import GlobalFlags
-from AthenaCommon import CfgMgr
-def ActsTrackingGeometrySvc(name="ActsTrackingGeometrySvc", **kwargs):
- return CfgMgr.ActsTrackingGeometrySvc(name, **kwargs)
+def ActsTrackingGeometrySvcCfg(configFlags, name = "ActsTrackingGeometrySvc" ) :
+ result = ComponentAccumulator()
+
+ Acts_ActsTrackingGeometrySvc = CompFactory.ActsTrackingGeometrySvc
+ subDetectors = []
+ if configFlags.Detector.GeometryPixel:
+ subDetectors += ["Pixel"]
+ if configFlags.Detector.GeometrySCT:
+ subDetectors += ["SCT"]
+ if configFlags.Detector.GeometryTRT:
+ subDetectors += ["TRT"]
+ if configFlags.Detector.GeometryCalo:
+ subDetectors += ["Calo"]
+
+ actsTrackingGeometrySvc = Acts_ActsTrackingGeometrySvc(name, BuildSubDetectors = subDetectors)
+ result.addService(actsTrackingGeometrySvc)
+ return result
+
+def ActsTrackingGeometryToolCfg(configFlags, name = "ActsTrackingGeometryTool" ) :
+ result = ComponentAccumulator()
+
+ acc = ActsTrackingGeometrySvcCfg(configFlags)
+ result.merge(acc)
+
+ Acts_ActsTrackingGeometryTool = CompFactory.ActsTrackingGeometryTool
+ actsTrackingGeometryTool = Acts_ActsTrackingGeometryTool(name)
+ result.addPublicTool(actsTrackingGeometryTool)
+
+ return result, actsTrackingGeometryTool
+
+def NominalAlignmentCondAlgCfg(configFlags, name = "NominalAlignmentCondAlg", **kwargs) :
+ result = ComponentAccumulator()
+
+ acc = ActsTrackingGeometrySvcCfg(configFlags)
+ result.merge(acc)
+
+ Acts_NominalAlignmentCondAlg = CompFactory.NominalAlignmentCondAlg
+ nominalAlignmentCondAlg = Acts_NominalAlignmentCondAlg(name, **kwargs)
+ result.addCondAlgo(nominalAlignmentCondAlg)
+
+ return result
+
+from MagFieldServices.MagFieldServicesConfig import MagneticFieldSvcCfg
+def ActsExtrapolationToolCfg(configFlags, name = "ActsExtrapolationTool" ) :
+ result=ComponentAccumulator()
+
+ acc = MagneticFieldSvcCfg(configFlags)
+ result.merge(acc)
+
+ acc, actsTrackingGeometryTool = ActsTrackingGeometryToolCfg(configFlags)
+ result.merge(acc)
+
+ Acts_ActsExtrapolationTool = CompFactory.ActsExtrapolationTool
+ actsExtrapolationTool = Acts_ActsExtrapolationTool(name)
+ result.addPublicTool(actsExtrapolationTool, primary=True)
+ return result
diff --git a/Tracking/Acts/ActsGeometry/src/ActsExtrapolationAlg.cxx b/Tracking/Acts/ActsGeometry/src/ActsExtrapolationAlg.cxx
index c5804e38955135adbb3d617e931789c17490eb70..beb72522c2ef2eab451bf0d393948deca1267ee2 100755
--- a/Tracking/Acts/ActsGeometry/src/ActsExtrapolationAlg.cxx
+++ b/Tracking/Acts/ActsGeometry/src/ActsExtrapolationAlg.cxx
@@ -114,7 +114,7 @@ StatusCode ActsExtrapolationAlg::execute(const EventContext &ctx) const {
auto anygctx = gctx.any();
Acts::BoundParameters startParameters(
anygctx, std::move(cov), std::move(pars), std::move(surface));
- steps = m_extrapolationTool->propagate(ctx, startParameters);
+ steps = m_extrapolationTool->propagationSteps(ctx, startParameters);
m_propStepWriterSvc->write(steps);
}
diff --git a/Tracking/Acts/ActsGeometry/src/ActsExtrapolationTool.cxx b/Tracking/Acts/ActsGeometry/src/ActsExtrapolationTool.cxx
index 6f518b3f150622fc00a14fd5d11bcb70a741bac0..328085f7bd849af49f67b8db3faac313522e473c 100644
--- a/Tracking/Acts/ActsGeometry/src/ActsExtrapolationTool.cxx
+++ b/Tracking/Acts/ActsGeometry/src/ActsExtrapolationTool.cxx
@@ -109,9 +109,10 @@ ActsExtrapolationTool::initialize()
std::vector<Acts::detail::Step>
-ActsExtrapolationTool::propagate(const EventContext& ctx,
- const Acts::BoundParameters& startParameters,
- double pathLimit /*= std::numeric_limits<double>::max()*/) const
+ActsExtrapolationTool::propagationSteps(const EventContext& ctx,
+ const Acts::BoundParameters& startParameters,
+ Acts::NavigationDirection navDir /*= Acts::forward*/,
+ double pathLimit /*= std::numeric_limits<double>::max()*/) const
{
using namespace Acts::UnitLiterals;
ATH_MSG_VERBOSE(name() << "::" << __FUNCTION__ << " begin");
@@ -137,6 +138,7 @@ ActsExtrapolationTool::propagate(const EventContext& ctx,
= (Acts::VectorHelpers::perp(startParameters.momentum())
< m_ptLoopers * 1_MeV);
options.maxStepSize = m_maxStepSize * 1_m;
+ options.direction = navDir;
std::vector<Acts::detail::Step> steps;
DebugOutput::result_type debugOutput;
@@ -170,3 +172,164 @@ ActsExtrapolationTool::propagate(const EventContext& ctx,
return steps;
}
+
+
+std::unique_ptr<const Acts::CurvilinearParameters>
+ActsExtrapolationTool::propagate(const EventContext& ctx,
+ const Acts::BoundParameters& startParameters,
+ Acts::NavigationDirection navDir /*= Acts::forward*/,
+ double pathLimit /*= std::numeric_limits<double>::max()*/) const
+{
+ using namespace Acts::UnitLiterals;
+ ATH_MSG_VERBOSE(name() << "::" << __FUNCTION__ << " begin");
+
+ Acts::MagneticFieldContext mctx;
+ const ActsGeometryContext& gctx
+ = m_trackingGeometryTool->getGeometryContext(ctx);
+
+ auto anygctx = gctx.any();
+
+ // Action list and abort list
+ using ActionList = Acts::ActionList<DebugOutput>;
+ using AbortConditions = Acts::AbortList<EndOfWorld>;
+ using Options = Acts::PropagatorOptions<ActionList, AbortConditions>;
+
+ Options options(anygctx, mctx);
+ options.pathLimit = pathLimit;
+ bool debug = msg().level() == MSG::VERBOSE;
+ options.debug = debug;
+
+ options.loopProtection
+ = (Acts::VectorHelpers::perp(startParameters.momentum())
+ < m_ptLoopers * 1_MeV);
+ options.maxStepSize = m_maxStepSize * 1_m;
+ options.direction = navDir;
+
+ std::vector<Acts::detail::Step> steps;
+ DebugOutput::result_type debugOutput;
+
+ auto parameters = boost::apply_visitor([&](const auto& propagator) -> std::unique_ptr<const Acts::CurvilinearParameters> {
+ auto result = propagator.propagate(startParameters, options);
+ if (!result.ok()) {
+ ATH_MSG_ERROR("Got error during propagation:" << result.error()
+ << ". Returning empty parameters.");
+ return nullptr;
+ }
+ return std::move(result.value().endParameters);
+ }, *m_varProp);
+
+ return parameters;
+}
+
+std::vector<Acts::detail::Step>
+ActsExtrapolationTool::propagationSteps(const EventContext& ctx,
+ const Acts::BoundParameters& startParameters,
+ const Acts::Surface& target,
+ Acts::NavigationDirection navDir /*= Acts::forward*/,
+ double pathLimit /*= std::numeric_limits<double>::max()*/) const
+{
+ using namespace Acts::UnitLiterals;
+ ATH_MSG_VERBOSE(name() << "::" << __FUNCTION__ << " begin");
+
+ Acts::MagneticFieldContext mctx;
+ const ActsGeometryContext& gctx
+ = m_trackingGeometryTool->getGeometryContext(ctx);
+
+ auto anygctx = gctx.any();
+
+ // Action list and abort list
+ using ActionList = Acts::ActionList<SteppingLogger, DebugOutput>;
+ using AbortConditions = Acts::AbortList<EndOfWorld>;
+ using Options = Acts::PropagatorOptions<ActionList, AbortConditions>;
+
+ Options options(anygctx, mctx);
+ options.pathLimit = pathLimit;
+ bool debug = msg().level() == MSG::VERBOSE;
+ options.debug = debug;
+
+ options.loopProtection
+ = (Acts::VectorHelpers::perp(startParameters.momentum())
+ < m_ptLoopers * 1_MeV);
+ options.maxStepSize = m_maxStepSize * 1_m;
+ options.direction = navDir;
+
+ std::vector<Acts::detail::Step> steps;
+ DebugOutput::result_type debugOutput;
+
+ auto res = boost::apply_visitor([&](const auto& propagator) -> ResultType {
+ auto result = propagator.propagate(startParameters, target, options);
+ if (!result.ok()) {
+ return result.error();
+ }
+ auto& propRes = *result;
+
+ auto steppingResults = propRes.template get<SteppingLogger::result_type>();
+ auto debugOutput = propRes.template get<DebugOutput::result_type>();
+ // try to force return value optimization, not sure this is necessary
+ return std::make_pair(std::move(steppingResults.steps), std::move(debugOutput));
+ }, *m_varProp);
+
+ if (!res.ok()) {
+ ATH_MSG_ERROR("Got error during propagation:" << res.error()
+ << ". Returning empty step vector.");
+ return {};
+ }
+ std::tie(steps, debugOutput) = std::move(*res);
+
+ if(debug) {
+ ATH_MSG_VERBOSE(debugOutput.debugString);
+ }
+
+ ATH_MSG_VERBOSE("Collected " << steps.size() << " steps");
+ ATH_MSG_VERBOSE(name() << "::" << __FUNCTION__ << " end");
+
+ return steps;
+}
+
+std::unique_ptr<const Acts::BoundParameters>
+ActsExtrapolationTool::propagate(const EventContext& ctx,
+ const Acts::BoundParameters& startParameters,
+ const Acts::Surface& target,
+ Acts::NavigationDirection navDir /*= Acts::forward*/,
+ double pathLimit /*= std::numeric_limits<double>::max()*/) const
+{
+ using namespace Acts::UnitLiterals;
+ ATH_MSG_VERBOSE(name() << "::" << __FUNCTION__ << " begin");
+
+ Acts::MagneticFieldContext mctx;
+ const ActsGeometryContext& gctx
+ = m_trackingGeometryTool->getGeometryContext(ctx);
+
+ auto anygctx = gctx.any();
+
+ // Action list and abort list
+ using ActionList = Acts::ActionList<DebugOutput>;
+ using AbortConditions = Acts::AbortList<EndOfWorld>;
+ using Options = Acts::PropagatorOptions<ActionList, AbortConditions>;
+
+ Options options(anygctx, mctx);
+ options.pathLimit = pathLimit;
+ bool debug = msg().level() == MSG::VERBOSE;
+ options.debug = debug;
+
+ options.loopProtection
+ = (Acts::VectorHelpers::perp(startParameters.momentum())
+ < m_ptLoopers * 1_MeV);
+ options.maxStepSize = m_maxStepSize * 1_m;
+ options.direction = navDir;
+
+ std::vector<Acts::detail::Step> steps;
+ DebugOutput::result_type debugOutput;
+
+ auto parameters = boost::apply_visitor([&](const auto& propagator) -> std::unique_ptr<const Acts::BoundParameters> {
+ auto result = propagator.propagate(startParameters, target, options);
+ if (!result.ok()) {
+ ATH_MSG_ERROR("Got error during propagation: " << result.error()
+ << ". Returning empty parameters.");
+ return nullptr;
+ }
+ return std::move(result.value().endParameters);
+ }, *m_varProp);
+
+ return parameters;
+}
diff --git a/Tracking/Acts/ActsGeometryInterfaces/ActsGeometryInterfaces/IActsExtrapolationTool.h b/Tracking/Acts/ActsGeometryInterfaces/ActsGeometryInterfaces/IActsExtrapolationTool.h
index 7b3c790f162444e9c1670a70c428f6aeff837a67..03ecfe769983bd51d4b3327db61b98c2b92661a2 100644
--- a/Tracking/Acts/ActsGeometryInterfaces/ActsGeometryInterfaces/IActsExtrapolationTool.h
+++ b/Tracking/Acts/ActsGeometryInterfaces/ActsGeometryInterfaces/IActsExtrapolationTool.h
@@ -28,8 +28,32 @@ class IActsExtrapolationTool : virtual public IAlgTool {
virtual
std::vector<Acts::detail::Step>
+ propagationSteps(const EventContext& ctx,
+ const Acts::BoundParameters& startParameters,
+ Acts::NavigationDirection navDir = Acts::forward,
+ double pathLimit = std::numeric_limits<double>::max()) const = 0;
+
+ virtual
+ std::unique_ptr<const Acts::CurvilinearParameters>
+ propagate(const EventContext& ctx,
+ const Acts::BoundParameters& startParameters,
+ Acts::NavigationDirection navDir = Acts::forward,
+ double pathLimit = std::numeric_limits<double>::max()) const = 0;
+
+ virtual
+ std::vector<Acts::detail::Step>
+ propagationSteps(const EventContext& ctx,
+ const Acts::BoundParameters& startParameters,
+ const Acts::Surface& target,
+ Acts::NavigationDirection navDir = Acts::forward,
+ double pathLimit = std::numeric_limits<double>::max()) const = 0;
+
+ virtual
+ std::unique_ptr<const Acts::BoundParameters>
propagate(const EventContext& ctx,
const Acts::BoundParameters& startParameters,
+ const Acts::Surface& target,
+ Acts::NavigationDirection navDir = Acts::forward,
double pathLimit = std::numeric_limits<double>::max()) const = 0;
virtual
diff --git a/Tracking/TrkExtrapolation/TrkExAlgs/CMakeLists.txt b/Tracking/TrkExtrapolation/TrkExAlgs/CMakeLists.txt
index 7b2e37026a48651565105d5234225a783e7482aa..d93a99657fb0dd0309d67258a9a549e052c5eb2d 100644
--- a/Tracking/TrkExtrapolation/TrkExAlgs/CMakeLists.txt
+++ b/Tracking/TrkExtrapolation/TrkExAlgs/CMakeLists.txt
@@ -22,18 +22,25 @@ atlas_depends_on_subdirs( PUBLIC
Tracking/TrkEvent/TrkEventPrimitives
Tracking/TrkEvent/TrkTrack
Tracking/TrkExtrapolation/TrkExInterfaces
- Tracking/TrkExtrapolation/TrkExUtils )
+ Tracking/TrkExtrapolation/TrkExUtils
+ ActsGeometryInterfaces
+ ActsGeometry
+ ActsInterop)
# External dependencies:
find_package( Eigen )
find_package( ROOT COMPONENTS Core Tree MathCore Hist RIO pthread )
+find_package( Acts COMPONENTS Core )
+
# Component(s) in the package:
atlas_add_component( TrkExAlgs
src/*.cxx
src/components/*.cxx
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${EIGEN_INCLUDE_DIRS}
- LINK_LIBRARIES ${ROOT_LIBRARIES} ${EIGEN_LIBRARIES} AthenaBaseComps AthContainers GeoPrimitives EventPrimitives GaudiKernel MagFieldInterfaces TrkSurfaces TrkParameters StoreGateLib SGtests TrkGeometry TrkVolumes TrkEventPrimitives TrkTrack TrkExInterfaces TrkExUtils )
+ LINK_LIBRARIES ${ROOT_LIBRARIES} ${EIGEN_LIBRARIES} AthenaBaseComps AthContainers GeoPrimitives EventPrimitives GaudiKernel MagFieldInterfaces TrkSurfaces TrkParameters StoreGateLib SGtests TrkGeometry TrkVolumes TrkEventPrimitives TrkTrack TrkExInterfaces TrkExUtils ActsGeometryInterfacesLib ActsInteropLib ActsGeometryLib ActsCore)
+
+atlas_install_python_modules( python/*.py )
# Install files from the package:
atlas_install_headers( TrkExAlgs )
diff --git a/Tracking/TrkExtrapolation/TrkExAlgs/TrkExAlgs/ExtrapolatorComparisonTest.h b/Tracking/TrkExtrapolation/TrkExAlgs/TrkExAlgs/ExtrapolatorComparisonTest.h
new file mode 100644
index 0000000000000000000000000000000000000000..aa969d304ab7560608f5adf3ba6dad5e0698ada4
--- /dev/null
+++ b/Tracking/TrkExtrapolation/TrkExAlgs/TrkExAlgs/ExtrapolatorComparisonTest.h
@@ -0,0 +1,130 @@
+/*
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// ExtrapolatorComparisonTest.h, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+#ifndef TRKEXALGS_TRKEXTRAPOLATORCOMPARISONTEST_H
+#define TRKEXALGS_TRKEXTRAPOLATORCOMPARISONTEST_H
+
+// Gaudi includes
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#include "GaudiKernel/MsgStream.h"
+
+// ATHENA
+#include "GaudiKernel/ServiceHandle.h"
+#include "GaudiKernel/Property.h" /*no forward decl: typedef*/
+#include "GaudiKernel/ISvcLocator.h"
+#include "GaudiKernel/EventContext.h"
+#include "GaudiKernel/RndmGenerators.h"
+#include "TrkExInterfaces/IExtrapolator.h"
+// ACTS
+#include "ActsGeometryInterfaces/IActsExtrapolationTool.h"
+#include "ActsGeometry/ActsGeometryContext.h"
+#include "TrkExAlgs/PropResultRootWriterSvc.h"
+// STL
+#include <memory>
+#include <vector>
+#include <fstream>
+#include <mutex>
+
+namespace Acts {
+ class TrackingGeometry;
+
+ namespace detail {
+ class Step;
+ }
+}
+
+struct perigeeParameters {
+ double m_d0 ;
+ double m_z0 ;
+ double m_phi ;
+ double m_eta ;
+ double m_pt ;
+ double m_charge;
+ perigeeParameters(double d0, double z0, double phi, double eta, double pt, double charge):
+ m_d0(d0), m_z0(z0), m_phi(phi), m_eta(eta), m_pt(pt), m_charge(charge) {}
+};
+
+namespace Trk
+{
+ /** @class ExtrapolatorComparisonTest
+
+ The ExtrapolatorComparisonTest Algorithm runs a number of n test extrapolations
+ from randomly distributed Track Parameters to reference surfcas within
+
+ - a) the Inner Detector if DetFlags.ID_On()
+ - b) the Calorimeter if DetFlags.Calo_On()
+ - c) the Muon System if DetFlags.Muon_On()
+
+ It uses and compares the output obtaining using
+ the ATLAS Extrapolator and the ACTS Extrapolation
+
+ */
+
+ class ExtrapolatorComparisonTest : public AthReentrantAlgorithm
+ {
+ public:
+
+ /** Standard Athena-Algorithm Constructor */
+ ExtrapolatorComparisonTest(const std::string& name, ISvcLocator* pSvcLocator);
+ /** Default Destructor */
+ ~ExtrapolatorComparisonTest();
+
+ /** standard Athena-Algorithm method */
+ StatusCode initialize() override;
+ /** standard Athena-Algorithm method */
+ StatusCode execute(const EventContext& ctx) const override;
+ /** standard Athena-Algorithm method */
+ StatusCode finalize() override;
+
+ private:
+
+ void generatePerigee(std::vector<perigeeParameters>& parameters);
+
+ /** The ACTS ExtrapolationTool to be retrieved */
+ ToolHandle<IActsExtrapolationTool> m_extrapolationTool{this, "ExtrapolationTool", "ActsExtrapolationTool"};
+
+ /** The ATLAS Extrapolator to be retrieved */
+ ToolHandle<Trk::IExtrapolator> m_atlasExtrapolator {this, "Extrapolator", "Trk::Extrapolator/AtlasExtrapolator"};
+
+
+ /** Random Number setup */
+ Rndm::Numbers* m_gaussDist;
+ Rndm::Numbers* m_flatDist;
+
+ double m_sigmaD0; //!< Sigma of distribution for D0
+ double m_sigmaZ0; //!< Sigma of distribution for Z0
+ double m_minPhi; //!< Minimal phi value
+ double m_maxPhi; //!< Maximal phi value
+ double m_minEta; //!< Minimal eta value
+ double m_maxEta; //!< Maximal eta value
+ double m_minPt; //!< Minimal p value
+ double m_maxPt; //!< Maximal p value
+
+ int m_particleType; //!< the particle typre for the extrap.
+
+ /** member variables for algorithm properties: */
+ unsigned int m_referenceSurfaces;
+
+ std::vector<double> m_referenceSurfaceRadius;
+ std::vector<double> m_referenceSurfaceHalflength;
+
+ std::vector<double> m_referenceSurfaceNegativeBoundary;
+ std::vector<double> m_referenceSurfacePositiveBoundary;
+
+ std::vector< std::vector<const Trk::Surface*> > m_atlasReferenceSurfaceTriples;
+ std::vector< std::vector< std::shared_ptr<const Acts::Surface> > > m_actsReferenceSurfaceTriples;
+
+ int m_eventsPerExecute;
+
+ ServiceHandle<PropResultRootWriterSvc> m_atlasPropResultWriterSvc;
+ ServiceHandle<PropResultRootWriterSvc> m_actsPropResultWriterSvc;
+
+ };
+} // end of namespace
+
+#endif
diff --git a/Tracking/TrkExtrapolation/TrkExAlgs/TrkExAlgs/PropResultRootWriterSvc.h b/Tracking/TrkExtrapolation/TrkExAlgs/TrkExAlgs/PropResultRootWriterSvc.h
new file mode 100644
index 0000000000000000000000000000000000000000..edbc78a276b0c752198a7ab0148f3a813f046828
--- /dev/null
+++ b/Tracking/TrkExtrapolation/TrkExAlgs/TrkExAlgs/PropResultRootWriterSvc.h
@@ -0,0 +1,211 @@
+
+/*
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef TRKEXALGS_PROPRESULTROOTWRITERSVC_H
+#define TRKEXALGS_PROPRESULTROOTWRITERSVC_H
+
+#include "GaudiKernel/IInterface.h"
+#include "GaudiKernel/Property.h" /*no forward decl: typedef*/
+#include "GaudiKernel/ITHistSvc.h"
+#include "AthenaBaseComps/AthService.h"
+
+#include "TTree.h"
+
+#include <mutex>
+
+class TFile;
+
+namespace Trk {
+
+class PropResultRootWriterSvc : public AthService {
+ public:
+ /// constructor
+ PropResultRootWriterSvc(const std::string& name, ISvcLocator* svcloc);
+ /// destructor
+ virtual ~PropResultRootWriterSvc();
+
+ /// Retrieve interface ID
+ static const InterfaceID& interfaceID() {
+ /// Declaration of the interface ID ( interface id, major version, minor version)
+ static const InterfaceID IID_PropResultRootWriterSvc("PropResultRootWriterSvc", 1, 0);
+ return IID_PropResultRootWriterSvc;
+ }
+
+ virtual StatusCode queryInterface(const InterfaceID& riid, void** ppvInterface) {
+ ATH_MSG_DEBUG("in queryInterface()");
+ if (PropResultRootWriterSvc::interfaceID().versionMatch(riid)) {
+ ATH_MSG_DEBUG("matched PropResultRootWriterSvc");
+ *ppvInterface=(PropResultRootWriterSvc*)this;
+ } else {
+ ATH_MSG_DEBUG("Default to Service interface");
+ return AthService::queryInterface(riid,ppvInterface);
+ }
+ return StatusCode::SUCCESS;
+ }
+
+ virtual StatusCode initialize() override;
+ virtual StatusCode finalize() override;
+
+ template <typename T>
+ void write(const T* initialPerigee,
+ const T* fwdParameters=nullptr, double fwdtime=std::numeric_limits<float>::quiet_NaN(),
+ const T* bkwParameters=nullptr, double bkwtime=std::numeric_limits<float>::quiet_NaN() );
+
+private:
+ std::mutex m_writeMutex;
+
+ // jobOptions properties
+ ServiceHandle<ITHistSvc> m_thistSvc;
+ TTree* m_tree;
+ Gaudi::Property<std::string> m_ntupleDirName{this, "DirName", "/ExtrapolationStudies/", ""};
+ Gaudi::Property<std::string> m_treeName{this, "TreeName", "ATLAS", ""};
+
+ int m_eventNum;
+
+ float m_start_d0 ;
+ float m_start_z0 ;
+ float m_start_phi ;
+ float m_start_theta ;
+ float m_start_qop ;
+
+ float m_fwd_success ;
+ float m_fwd_time ;
+ float m_fwd_final_l0 ;
+ float m_fwd_final_l1 ;
+ float m_fwd_final_x ;
+ float m_fwd_final_y ;
+ float m_fwd_final_z ;
+ float m_fwd_final_phi ;
+ float m_fwd_final_theta ;
+ float m_fwd_final_qop ;
+ float m_fwd_final_sigma_l0 ;
+ float m_fwd_final_sigma_l1 ;
+ float m_fwd_final_sigma_phi ;
+ float m_fwd_final_sigma_theta ;
+ float m_fwd_final_sigma_qop ;
+
+ float m_bkw_success ;
+ float m_bkw_time ;
+ float m_bkw_final_d0 ;
+ float m_bkw_final_z0 ;
+ float m_bkw_final_phi ;
+ float m_bkw_final_theta ;
+ float m_bkw_final_qop ;
+ float m_bkw_final_sigma_d0 ;
+ float m_bkw_final_sigma_z0 ;
+ float m_bkw_final_sigma_phi ;
+ float m_bkw_final_sigma_theta ;
+ float m_bkw_final_sigma_qop ;
+
+ };
+
+ template <typename T>
+ void PropResultRootWriterSvc::write(const T* initialPerigee,
+ const T* fwdParameters, double fwdtime,
+ const T* bkwParameters, double bkwtime ) {
+
+ auto ctx = Gaudi::Hive::currentContext();
+
+ m_eventNum = ctx.eventID().event_number();
+
+ if (initialPerigee) {
+ m_start_d0 = initialPerigee->parameters()[0];
+ m_start_z0 = initialPerigee->parameters()[1];
+ m_start_phi = initialPerigee->parameters()[2];
+ m_start_theta= initialPerigee->parameters()[3];
+ m_start_qop = initialPerigee->parameters()[4];
+ } else {
+ m_start_d0 = std::numeric_limits<float>::quiet_NaN();
+ m_start_z0 = std::numeric_limits<float>::quiet_NaN();
+ m_start_phi = std::numeric_limits<float>::quiet_NaN();
+ m_start_theta= std::numeric_limits<float>::quiet_NaN();
+ m_start_qop = std::numeric_limits<float>::quiet_NaN();
+ }
+
+ if (fwdParameters) {
+ m_fwd_success = 1;
+ m_fwd_time = fwdtime;
+ m_fwd_final_l0 = fwdParameters->parameters()[0];
+ m_fwd_final_l1 = fwdParameters->parameters()[1];
+ m_fwd_final_phi = fwdParameters->parameters()[2];
+ m_fwd_final_theta = fwdParameters->parameters()[3];
+ m_fwd_final_qop = fwdParameters->parameters()[4];
+ m_fwd_final_x = fwdParameters->position()[0];
+ m_fwd_final_y = fwdParameters->position()[1];
+ m_fwd_final_z = fwdParameters->position()[2];
+ if (fwdParameters->covariance()) {
+ m_fwd_final_sigma_l0 = (*fwdParameters->covariance())(0,0);
+ m_fwd_final_sigma_l1 = (*fwdParameters->covariance())(1,1);
+ m_fwd_final_sigma_phi = (*fwdParameters->covariance())(2,2);
+ m_fwd_final_sigma_theta = (*fwdParameters->covariance())(3,3);
+ m_fwd_final_sigma_qop = (*fwdParameters->covariance())(4,4);
+ } else {
+ m_fwd_final_sigma_l0 = std::numeric_limits<float>::quiet_NaN();
+ m_fwd_final_sigma_l1 = std::numeric_limits<float>::quiet_NaN();
+ m_fwd_final_sigma_phi = std::numeric_limits<float>::quiet_NaN();
+ m_fwd_final_sigma_theta = std::numeric_limits<float>::quiet_NaN();
+ m_fwd_final_sigma_qop = std::numeric_limits<float>::quiet_NaN();
+ }
+ } else {
+ m_fwd_success = std::numeric_limits<float>::quiet_NaN();
+ m_fwd_time = std::numeric_limits<float>::quiet_NaN();
+ m_fwd_final_l0 = std::numeric_limits<float>::quiet_NaN();
+ m_fwd_final_l1 = std::numeric_limits<float>::quiet_NaN();
+ m_fwd_final_phi = std::numeric_limits<float>::quiet_NaN();
+ m_fwd_final_theta = std::numeric_limits<float>::quiet_NaN();
+ m_fwd_final_qop = std::numeric_limits<float>::quiet_NaN();
+ m_fwd_final_x = std::numeric_limits<float>::quiet_NaN();
+ m_fwd_final_y = std::numeric_limits<float>::quiet_NaN();
+ m_fwd_final_z = std::numeric_limits<float>::quiet_NaN();
+ m_fwd_final_sigma_l0 = std::numeric_limits<float>::quiet_NaN();
+ m_fwd_final_sigma_l1 = std::numeric_limits<float>::quiet_NaN();
+ m_fwd_final_sigma_phi = std::numeric_limits<float>::quiet_NaN();
+ m_fwd_final_sigma_theta = std::numeric_limits<float>::quiet_NaN();
+ m_fwd_final_sigma_qop = std::numeric_limits<float>::quiet_NaN();
+ }
+
+ if (bkwParameters) {
+ m_bkw_success = 1;
+ m_bkw_time = bkwtime;
+ m_bkw_final_d0 = bkwParameters->parameters()[0];
+ m_bkw_final_z0 = bkwParameters->parameters()[1];
+ m_bkw_final_phi = bkwParameters->parameters()[2];
+ m_bkw_final_theta = bkwParameters->parameters()[3];
+ m_bkw_final_qop = bkwParameters->parameters()[4];
+ if (bkwParameters->covariance()) {
+ m_bkw_final_sigma_d0 = (*bkwParameters->covariance())(0,0);
+ m_bkw_final_sigma_z0 = (*bkwParameters->covariance())(1,1);
+ m_bkw_final_sigma_phi = (*bkwParameters->covariance())(2,2);
+ m_bkw_final_sigma_theta = (*bkwParameters->covariance())(3,3);
+ m_bkw_final_sigma_qop = (*bkwParameters->covariance())(4,4);
+ } else {
+ m_bkw_final_sigma_d0 = std::numeric_limits<float>::quiet_NaN();
+ m_bkw_final_sigma_z0 = std::numeric_limits<float>::quiet_NaN();
+ m_bkw_final_sigma_phi = std::numeric_limits<float>::quiet_NaN();
+ m_bkw_final_sigma_theta = std::numeric_limits<float>::quiet_NaN();
+ m_bkw_final_sigma_qop = std::numeric_limits<float>::quiet_NaN();
+ }
+ } else {
+ m_bkw_success = std::numeric_limits<float>::quiet_NaN();
+ m_bkw_time = std::numeric_limits<float>::quiet_NaN();
+ m_bkw_final_d0 = std::numeric_limits<float>::quiet_NaN();
+ m_bkw_final_z0 = std::numeric_limits<float>::quiet_NaN();
+ m_bkw_final_phi = std::numeric_limits<float>::quiet_NaN();
+ m_bkw_final_theta = std::numeric_limits<float>::quiet_NaN();
+ m_bkw_final_qop = std::numeric_limits<float>::quiet_NaN();
+ m_bkw_final_sigma_d0 = std::numeric_limits<float>::quiet_NaN();
+ m_bkw_final_sigma_z0 = std::numeric_limits<float>::quiet_NaN();
+ m_bkw_final_sigma_phi = std::numeric_limits<float>::quiet_NaN();
+ m_bkw_final_sigma_theta = std::numeric_limits<float>::quiet_NaN();
+ m_bkw_final_sigma_qop = std::numeric_limits<float>::quiet_NaN();
+ }
+ m_tree->Fill();
+ }
+
+
+
+}
+
+#endif
diff --git a/Tracking/TrkExtrapolation/TrkExAlgs/python/TrkExAlgsConfig.py b/Tracking/TrkExtrapolation/TrkExAlgs/python/TrkExAlgsConfig.py
new file mode 100644
index 0000000000000000000000000000000000000000..390f038c8eb32d9d404c90dd3136befdaa6012c1
--- /dev/null
+++ b/Tracking/TrkExtrapolation/TrkExAlgs/python/TrkExAlgsConfig.py
@@ -0,0 +1,123 @@
+# Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+
+from ActsGeometry.ActsGeometryConfig import ActsExtrapolationToolCfg
+from ActsGeometry.ActsGeometryConfig import NominalAlignmentCondAlgCfg
+from TrkConfig.AtlasExtrapolatorConfig import AtlasExtrapolatorCfg
+
+def PropResultRootWriterSvcCfg(configFlags, name="PropResultRootWriterSvc", **kwargs) :
+ result = ComponentAccumulator()
+
+ Trk__PropResultRootWriterSvc = CompFactory.Trk.PropResultRootWriterSvc
+ propResultRootWriterSvc = Trk__PropResultRootWriterSvc(name, **kwargs)
+ result.addService(propResultRootWriterSvc)
+
+ return result, propResultRootWriterSvc
+
+
+def ExtrapolatorComparisonTestCfg(configFlags, name = "ExtrapolatorComparisonTest", **kwargs ) :
+ result=ComponentAccumulator()
+
+ if "Extrapolator" not in kwargs:
+ kwargs["Extrapolator"] = result.popToolsAndMerge(AtlasExtrapolatorCfg(configFlags, 'AtlasExtrapolator'))
+
+ if "ExtrapolationTool" not in kwargs:
+ actsextrapolation = ActsExtrapolationToolCfg(configFlags)
+ kwargs["ExtrapolationTool"] = actsextrapolation.getPrimary()
+ result.merge(actsextrapolation)
+
+ Trk__ExtrapolatorComparisonTest = CompFactory.Trk.ExtrapolatorComparisonTest
+ extrapolationTest = Trk__ExtrapolatorComparisonTest(name, **kwargs)
+ result.addEventAlgo(extrapolationTest)
+
+ return result
+
+if __name__=="__main__":
+ from AthenaCommon.Configurable import Configurable
+ from AthenaCommon.Logging import log
+ from AthenaCommon.Constants import VERBOSE
+ from AthenaConfiguration.AllConfigFlags import ConfigFlags
+ from AthenaConfiguration.MainServicesConfig import MainServicesThreadedCfg
+
+ #log.setLevel(VERBOSE)
+
+ Configurable.configurableRun3Behavior = True
+
+ ## Just enable ID for the moment.
+ ConfigFlags.Input.isMC = True
+ ConfigFlags.Beam.Type = ''
+ ConfigFlags.GeoModel.AtlasVersion = "ATLAS-R2-2016-01-00-01"
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-SIM-00-00-00"
+ ConfigFlags.Detector.SimulateBpipe = True
+ ConfigFlags.Detector.SimulateID = True
+ ConfigFlags.Detector.GeometryBpipe = True
+ ConfigFlags.Detector.GeometryID = True
+ ConfigFlags.Detector.GeometryPixel = True
+ ConfigFlags.Detector.GeometrySCT = True
+ ConfigFlags.Detector.GeometryCalo = False
+ ConfigFlags.Detector.GeometryMuon = False
+ #ConfigFlags.Detector.GeometryTRT = True
+ ConfigFlags.TrackingGeometry.MaterialSource = "Input"
+
+ # This should run serially for the moment.
+ ConfigFlags.Concurrency.NumThreads = 1
+ ConfigFlags.Concurrency.NumConcurrentEvents = 1
+
+ log.debug('Lock config flags now.')
+ ConfigFlags.lock()
+
+ ConfigFlags.dump()
+
+ cfg=MainServicesThreadedCfg(ConfigFlags)
+
+ from AtlasGeoModel.InDetGMConfig import InDetGeometryCfg
+ cfg.merge(InDetGeometryCfg(ConfigFlags))
+
+ from BeamPipeGeoModel.BeamPipeGMConfig import BeamPipeGeometryCfg
+ cfg.merge(BeamPipeGeometryCfg(ConfigFlags))
+
+ nominalAlignmentCondAlg = NominalAlignmentCondAlgCfg(ConfigFlags, OutputLevel=VERBOSE)
+ cfg.merge(nominalAlignmentCondAlg)
+
+ histSvc = CompFactory.THistSvc(Output = ["ExtrapolationStudies DATAFILE='ExtrapolationStudies.root' OPT='RECREATE'"])
+ histSvc.OutputLevel=VERBOSE
+ cfg.addService( histSvc )
+
+ accAtlasPropResult, atlasPropResultRootWriterSvc = PropResultRootWriterSvcCfg(ConfigFlags,
+ name="ATLASPropResultRootWriterSvc",
+ TreeName="ATLAS")
+ cfg.merge(accAtlasPropResult)
+
+ accActsPropResult, actsPropResultRootWriterSvc = PropResultRootWriterSvcCfg(ConfigFlags,
+ name="ACTSPropResultRootWriterSvc",
+ TreeName="ACTS")
+ cfg.merge(accActsPropResult)
+
+ topoAcc=ExtrapolatorComparisonTestCfg(ConfigFlags,
+ EventsPerExecute = 1000,
+ StartPerigeeMinPt = 10000,
+ StartPerigeeMaxPt = 10000,
+ ReferenceSurfaceRadius = [80],
+ ReferenceSurfaceHalfZ = [500],
+ ATLASPropResultRootWriter = atlasPropResultRootWriterSvc,
+ ACTSPropResultRootWriter = actsPropResultRootWriterSvc)
+ cfg.merge(topoAcc)
+
+ cfg.printConfig()
+
+ #cfg.getPublicTool("InDetLayerArrayCreator").OutputLevel=VERBOSE
+ #cfg.getPublicTool("InDetTrackingVolumeArrayCreator").OutputLevel=VERBOSE
+ #cfg.getPublicTool("InDetCylinderVolumeCreator").OutputLevel=VERBOSE
+ cfg.getPublicTool("InDetCylinderVolumeCreator").LayerArrayCreator=cfg.getPublicTool("InDetLayerArrayCreator")
+ #cfg.getPublicTool("InDetTrackingGeometryBuilder").OutputLevel=VERBOSE
+ #cfg.getService("AtlasTrackingGeometrySvc").OutputLevel=VERBOSE
+
+ ## Building TRT straws
+ #cfg.getPublicTool("InDetTrackingGeometryBuilder").LayerBinningType = [2, 2, 2]
+ #cfg.getPublicTool("TRT_LayerBuilder").ModelLayersOnly = False
+
+ cfg.run(10)
+ f=open("ExtrapolatorComparisonTestConfig.pkl","w")
+ cfg.store(f)
+ f.close()
diff --git a/Tracking/TrkExtrapolation/TrkExAlgs/src/ExtrapolatorComparisonTest.cxx b/Tracking/TrkExtrapolation/TrkExAlgs/src/ExtrapolatorComparisonTest.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..692a3da6dbe6ba8ccdfcdd2649bcb0e0d301bed6
--- /dev/null
+++ b/Tracking/TrkExtrapolation/TrkExAlgs/src/ExtrapolatorComparisonTest.cxx
@@ -0,0 +1,358 @@
+/*
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// ExtrapolatorComparisonTest.cxx, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+#include "TrkExAlgs/ExtrapolatorComparisonTest.h"
+
+// Tracking
+#include "TrkSurfaces/PerigeeSurface.h"
+#include "TrkSurfaces/CylinderSurface.h"
+#include "TrkSurfaces/DiscSurface.h"
+#include "TrkEventPrimitives/ParticleHypothesis.h"
+#include "TrkParametersBase/ParametersT.h"
+#include "EventPrimitives/EventPrimitives.h"
+#include "MagFieldInterfaces/IMagFieldSvc.h"
+
+#include "GaudiKernel/SystemOfUnits.h"
+#include "GaudiKernel/ISvcLocator.h"
+#include "AthenaKernel/RNGWrapper.h"
+
+// ACTS
+#include "Acts/EventData/TrackParameters.hpp"
+#include "Acts/Surfaces/Surface.hpp"
+#include "Acts/Surfaces/PerigeeSurface.hpp"
+#include "Acts/Surfaces/DiscSurface.hpp"
+#include "Acts/Surfaces/CylinderSurface.hpp"
+#include "Acts/Utilities/Helpers.hpp"
+#include "Acts/Utilities/Logger.hpp"
+#include "ActsInterop/Logger.h"
+#include "Acts/Propagator/detail/SteppingLogger.hpp"
+#include "ActsGeometry/ActsGeometryContext.h"
+#include "ActsGeometryInterfaces/IActsTrackingGeometryTool.h"
+#include "Acts/Geometry/GeometryContext.hpp"
+#include "Acts/Utilities/Definitions.hpp"
+#include "Acts/Utilities/ParameterDefinitions.hpp"
+#include "Acts/Propagator/detail/SteppingLogger.hpp"
+
+// OTHER
+#include "CLHEP/Random/RandomEngine.h"
+
+// STL
+#include <string>
+#include <fstream>
+#include <chrono>
+#include <cmath>
+
+using xclock = std::chrono::steady_clock;
+
+//================ Constructor =================================================
+
+Trk::ExtrapolatorComparisonTest::ExtrapolatorComparisonTest(const std::string& name, ISvcLocator* pSvcLocator)
+ :
+ AthReentrantAlgorithm(name,pSvcLocator),
+ m_gaussDist(0),
+ m_flatDist(0),
+ m_sigmaD0(17.*Gaudi::Units::micrometer),
+ m_sigmaZ0(50.*Gaudi::Units::millimeter),
+ m_minPhi(-M_PI),
+ m_maxPhi(M_PI),
+ m_minEta(-3.),
+ m_maxEta(3.),
+ m_minPt(0.5*Gaudi::Units::GeV),
+ m_maxPt(50000.*Gaudi::Units::GeV),
+ m_particleType(2),
+ m_referenceSurfaces(0),
+ m_eventsPerExecute(1),
+ m_atlasPropResultWriterSvc("ATLASPropResultRootWriterSvc", name),
+ m_actsPropResultWriterSvc("ACTSPropResultRootWriterSvc", name)
+{
+ // algorithm steering
+ declareProperty("StartPerigeeSigmaD0" , m_sigmaD0 );
+ declareProperty("StartPerigeeSigmaZ0" , m_sigmaZ0 );
+ declareProperty("StartPerigeeMinPhi" , m_minPhi );
+ declareProperty("StartPerigeeMaxPhi" , m_maxPhi );
+ declareProperty("StartPerigeeMinEta" , m_minEta );
+ declareProperty("StartPerigeeMaxEta" , m_maxEta );
+ declareProperty("StartPerigeeMinPt" , m_minPt );
+ declareProperty("StartPerigeeMaxPt" , m_maxPt );
+ declareProperty("ParticleType" , m_particleType );
+ declareProperty("ReferenceSurfaceRadius" , m_referenceSurfaceRadius );
+ declareProperty("ReferenceSurfaceHalfZ" , m_referenceSurfaceHalflength);
+ declareProperty("EventsPerExecute" , m_eventsPerExecute );
+ declareProperty("ATLASPropResultRootWriter" , m_atlasPropResultWriterSvc );
+ declareProperty("ACTSPropResultRootWriter" , m_actsPropResultWriterSvc );
+}
+
+//================ Destructor =================================================
+
+Trk::ExtrapolatorComparisonTest::~ExtrapolatorComparisonTest()
+{
+ delete m_gaussDist;
+ delete m_flatDist;
+
+ // cleanup of the Trk::Surfaces
+ std::vector< std::vector< const Trk::Surface* > >::const_iterator atlasSurfaceTripleIter = m_atlasReferenceSurfaceTriples.begin();
+ std::vector< std::vector< const Trk::Surface* > >::const_iterator atlasSurfaceTripleIterEnd = m_atlasReferenceSurfaceTriples.end();
+ for ( ; atlasSurfaceTripleIter != atlasSurfaceTripleIterEnd; ++atlasSurfaceTripleIter) {
+ std::vector<const Trk::Surface*>::const_iterator surfIter = (*atlasSurfaceTripleIter).begin();
+ std::vector<const Trk::Surface*>::const_iterator surfIterEnd = (*atlasSurfaceTripleIter).end();
+ for ( ; surfIter != surfIterEnd; delete (*surfIter), ++surfIter);
+ }
+}
+
+
+//================ Initialisation =================================================
+
+StatusCode Trk::ExtrapolatorComparisonTest::initialize()
+{
+ // Code entered here will be executed once at program start.
+
+ ATH_MSG_INFO(" initialize()");
+
+ ATH_CHECK( m_extrapolationTool.retrieve() );
+ ATH_CHECK( m_atlasExtrapolator.retrieve() );
+
+ // Create the destination surfaces for extrapolation
+ // --> you need the Trk::Surfaces and the Acts::Surfaces
+ if (m_referenceSurfaceRadius.size() == m_referenceSurfaceHalflength.size()) {
+ // assign the size
+ m_referenceSurfaces = m_referenceSurfaceRadius.size();
+ // loop over it and create the
+ for (unsigned int surface = 0; surface < m_referenceSurfaces; surface++) {
+
+ double radius = m_referenceSurfaceRadius.at(surface);
+ double halfZ = m_referenceSurfaceHalflength.at(surface);
+
+ // create the Surface triplet
+ std::vector< const Trk::Surface*> trkSurfaceTriplet;
+ trkSurfaceTriplet.push_back(new Trk::DiscSurface (new Amg::Transform3D(Amg::Translation3D(0.,0., halfZ)), 0.,radius));
+ trkSurfaceTriplet.push_back(new Trk::CylinderSurface(new Amg::Transform3D(Amg::Translation3D(0.,0., 0.)),radius, halfZ));
+ trkSurfaceTriplet.push_back(new Trk::DiscSurface (new Amg::Transform3D(Amg::Translation3D(0.,0.,-halfZ)), 0.,radius));
+ ATH_MSG_INFO("Creating Trk::Surface at: R " << radius << " Z " << halfZ);
+ m_atlasReferenceSurfaceTriples.push_back(trkSurfaceTriplet);
+
+ // create the Surface triplet
+ std::vector<std::shared_ptr<const Acts::Surface>> actsSurfaceTriplet;
+
+ const Acts::Transform3D * posTransf = new Acts::Transform3D(Acts::Transform3D::Identity()*Acts::Translation3D(Acts::Vector3D(0.,0., halfZ)));
+ const Acts::Transform3D * cTransf = new Acts::Transform3D(Acts::Transform3D::Identity()*Acts::Translation3D(Acts::Vector3D(0.,0., 0.)));
+ const Acts::Transform3D * negTransf = new Acts::Transform3D(Acts::Transform3D::Identity()*Acts::Translation3D(Acts::Vector3D(0.,0.,-halfZ)));
+
+ auto posSurface = Acts::Surface::makeShared<Acts::DiscSurface> (std::shared_ptr<const Acts::Transform3D>(posTransf), 0.,radius);
+ auto cSurface = Acts::Surface::makeShared<Acts::CylinderSurface>(std::shared_ptr<const Acts::Transform3D>( cTransf),radius, halfZ);
+ auto negSurface = Acts::Surface::makeShared<Acts::DiscSurface> (std::shared_ptr<const Acts::Transform3D>(negTransf), 0.,radius);
+
+ actsSurfaceTriplet.push_back(posSurface);
+ actsSurfaceTriplet.push_back(cSurface );
+ actsSurfaceTriplet.push_back(negSurface);
+ ATH_MSG_INFO("Creating Acts::Surface at: R " << radius << " Z " << halfZ);
+ m_actsReferenceSurfaceTriples.push_back(actsSurfaceTriplet);
+
+ m_referenceSurfaceNegativeBoundary.push_back(atan2(radius,-halfZ));
+ m_referenceSurfacePositiveBoundary.push_back(atan2(radius, halfZ));
+ }
+ } else {
+ ATH_MSG_WARNING("Not compatible size of ReferenceSurfaceRadius and ReferenceSurfaceHalfZ!! Returning FAILURE!");
+ return StatusCode::FAILURE;
+ }
+
+ m_gaussDist = new Rndm::Numbers(randSvc(), Rndm::Gauss(0.,1.));
+ m_flatDist = new Rndm::Numbers(randSvc(), Rndm::Flat(0.,1.));
+
+ ATH_CHECK( m_atlasPropResultWriterSvc.retrieve() );
+ ATH_CHECK( m_actsPropResultWriterSvc.retrieve() );
+
+ ATH_MSG_INFO("initialize() successful!!");
+
+ return StatusCode::SUCCESS;
+}
+
+//================ Finalisation =================================================
+
+StatusCode Trk::ExtrapolatorComparisonTest::finalize()
+{
+ // Code entered here will be executed once at the end of the program run.
+ return StatusCode::SUCCESS;
+}
+
+//================ Execution ====================================================
+
+StatusCode Trk::ExtrapolatorComparisonTest::execute(const EventContext& ctx) const {
+
+ float milliseconds_to_seconds = 1000.;
+
+ // generate perigees with random number generator
+
+ std::vector<perigeeParameters> parameters = {};
+ for (int ext = 0; ext<m_eventsPerExecute; ext++) {
+ // generate with random number generator
+ double d0 = m_gaussDist->shoot() * m_sigmaD0;
+ double z0 = m_gaussDist->shoot() * m_sigmaZ0;
+ double phi = m_minPhi + (m_maxPhi-m_minPhi)* m_flatDist->shoot();
+ double eta = m_minEta + m_flatDist->shoot()*(m_maxEta-m_minEta);
+ double pt = m_minPt + m_flatDist->shoot()*(m_maxPt-m_minPt);
+ double charge = (m_flatDist->shoot() > 0.5 ) ? -1. : 1.;
+ parameters.push_back(perigeeParameters(d0, z0, phi, eta, pt, charge));
+ }
+
+ int n_extraps = 0;
+ auto start = xclock::now();
+ for (auto& perigee : parameters) {
+
+ Acts::Vector3D momentum(perigee.m_pt * std::cos(perigee.m_phi), perigee.m_pt * std::sin(perigee.m_phi), perigee.m_pt * std::sinh(perigee.m_eta));
+ double theta = Acts::VectorHelpers::theta(momentum);
+ double qOverP = perigee.m_charge / momentum.norm();
+
+ const Trk::PerigeeSurface atlPerigeeSurface;
+ const Trk::Perigee * atlPerigee = new const Trk::Perigee(perigee.m_d0, perigee.m_z0, perigee.m_phi, theta, qOverP, atlPerigeeSurface);
+
+ for (unsigned int surface = 0; surface < m_atlasReferenceSurfaceTriples.size(); surface++) {
+ n_extraps++;
+
+ double negRef = m_referenceSurfaceNegativeBoundary.at(surface);
+ double posRef = m_referenceSurfacePositiveBoundary.at(surface);
+
+ // decide which reference surface to take
+ int refSurface = theta < posRef ? 2 : 1;
+ refSurface = theta > negRef ? 0 : 1;
+
+ const Trk::Surface* destinationSurface = m_atlasReferenceSurfaceTriples.at(surface).at(refSurface);
+
+ ATH_MSG_VERBOSE("Starting extrapolation " << n_extraps << " from : " << *atlPerigee << " to : " << *destinationSurface);
+
+ auto start_fwd = xclock::now();
+ const Trk::TrackParameters* destParameters = m_atlasExtrapolator->extrapolate(*atlPerigee,
+ *destinationSurface,
+ Trk::alongMomentum,
+ true,
+ (Trk::ParticleHypothesis)m_particleType);
+ auto end_fwd = xclock::now();
+ float ms_fwd = std::chrono::duration_cast<std::chrono::milliseconds>(end_fwd-start_fwd).count();
+
+ if (destParameters) {
+ ATH_MSG_VERBOSE(" ATLAS Extrapolator succeded!! --> Forward" );
+ ATH_MSG_VERBOSE(" [ intersection ] with surface at (x,y,z) = " << destParameters->position().x() << ", " << destParameters->position().y() << ", " << destParameters->position().z() );
+ ATH_MSG_VERBOSE(" [ intersection ] parameters: " << destParameters->parameters() );
+ ATH_MSG_VERBOSE(" [ intersection ] cov matrix: " << destParameters->covariance() );
+
+ // now try backward extrapolation
+ auto start_bkw = xclock::now();
+ const Trk::TrackParameters* finalperigee = m_atlasExtrapolator->extrapolate(*destParameters,
+ atlPerigee->associatedSurface(),
+ Trk::oppositeMomentum,
+ true,
+ (Trk::ParticleHypothesis)m_particleType);
+ auto end_bkw = xclock::now();
+ float ms_bkw = std::chrono::duration_cast<std::chrono::milliseconds>(end_bkw-start_bkw).count();
+
+ if (finalperigee) {
+ ATH_MSG_VERBOSE(" ATLAS Extrapolator succeded!! --> Backward" );
+ ATH_MSG_VERBOSE(" [extrapolation to perigee] input: " << atlPerigee->parameters() );
+ ATH_MSG_VERBOSE(" [extrapolation to perigee] output: " << finalperigee->parameters() );
+ ATH_MSG_VERBOSE(" [extrapolation to perigee] cov matrix: " << finalperigee->covariance() );
+
+ } else if (!finalperigee) {
+ ATH_MSG_DEBUG(" ATLAS Extrapolation to perigee failed for input parameters: " << destParameters->parameters());
+ }
+
+ m_atlasPropResultWriterSvc->write<Trk::TrackParameters>(atlPerigee, destParameters, ms_fwd, finalperigee, ms_bkw);
+ delete finalperigee;
+ } else if (!destParameters) {
+ ATH_MSG_DEBUG(" ATLAS Extrapolation not successful! " );
+ m_atlasPropResultWriterSvc->write<Trk::TrackParameters>(atlPerigee);
+ }
+ delete destParameters;
+ }
+ delete atlPerigee;
+ }
+ auto end = xclock::now();
+ auto secs = std::chrono::duration_cast<std::chrono::milliseconds>(end-start).count() / milliseconds_to_seconds;
+ double secs_per_ex = secs / n_extraps;
+
+ ATH_MSG_INFO("ATLAS : Time for " << n_extraps << " iterations: " << secs << "s (" << secs_per_ex << "s per extrapolation)");
+
+ n_extraps = 0;
+ start = xclock::now();
+ for (auto& perigee : parameters) {
+ Acts::Vector3D momentum(perigee.m_pt * std::cos(perigee.m_phi), perigee.m_pt * std::sin(perigee.m_phi), perigee.m_pt * std::sinh(perigee.m_eta));
+ double theta = Acts::VectorHelpers::theta(momentum);
+ double qOverP = perigee.m_charge / momentum.norm();
+
+ std::shared_ptr<Acts::PerigeeSurface> actsPerigeeSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3D(0, 0, 0));
+ double t = 0.;
+ Acts::BoundVector pars;
+ pars << perigee.m_d0, perigee.m_z0, perigee.m_phi, theta, qOverP, t;
+ std::optional<Acts::BoundSymMatrix> cov = std::nullopt;
+
+ // Perigee, no alignment -> default geo context
+ ActsGeometryContext gctx = m_extrapolationTool->trackingGeometryTool()->getNominalGeometryContext();
+ auto anygctx = gctx.any();
+ const Acts::BoundParameters* startParameters = new const Acts::BoundParameters(anygctx, std::move(cov), std::move(pars), std::move(actsPerigeeSurface));
+
+ for (unsigned int surface = 0; surface < m_actsReferenceSurfaceTriples.size(); surface++) {
+ n_extraps++;
+
+ double negRef = m_referenceSurfaceNegativeBoundary.at(surface);
+ double posRef = m_referenceSurfacePositiveBoundary.at(surface);
+
+ // decide which reference surface to take
+ int refSurface = theta < posRef ? 2 : 1;
+ refSurface = theta > negRef ? 0 : 1;
+
+ auto destinationSurface = m_actsReferenceSurfaceTriples.at(surface).at(refSurface);
+
+ ATH_MSG_VERBOSE("Starting extrapolation " << n_extraps << " from : " << pars << " to : " << destinationSurface);
+
+ auto start_fwd = xclock::now();
+ auto destParameters = m_extrapolationTool->propagate(ctx, *startParameters, *destinationSurface, Acts::forward);
+ auto end_fwd = xclock::now();
+ float ms_fwd = std::chrono::duration_cast<std::chrono::milliseconds>(end_fwd-start_fwd).count();
+
+ if (destParameters) {
+ ATH_MSG_VERBOSE(" ACTS Extrapolator succeded!! --> Forward" );
+ ATH_MSG_VERBOSE(" [ intersection ] with surface at (x,y,z) = " << destParameters->position().x() << ", " << destParameters->position().y() << ", " << destParameters->position().z() );
+ ATH_MSG_VERBOSE(" [ intersection ] parameters: " << destParameters->parameters() );
+ ATH_MSG_VERBOSE(" [ intersection ] cov matrix: " << *destParameters->covariance() );
+
+ // now try backward extrapolation
+ auto start_bkw = xclock::now();
+ auto finalperigee = m_extrapolationTool->propagate(ctx, *destParameters, startParameters->referenceSurface(), Acts::backward);
+ auto end_bkw = xclock::now();
+ float ms_bkw = std::chrono::duration_cast<std::chrono::milliseconds>(end_bkw-start_bkw).count();
+
+ if (finalperigee) {
+ ATH_MSG_VERBOSE(" ACTS Extrapolator succeded!! --> Backward" );
+ ATH_MSG_VERBOSE(" [extrapolation to perigee] input: " << startParameters->parameters() );
+ ATH_MSG_VERBOSE(" [extrapolation to perigee] output: " << finalperigee->parameters() );
+ ATH_MSG_VERBOSE(" [extrapolation to perigee] cov matrix: " << *finalperigee->covariance() );
+
+ } else if (!finalperigee) {
+ ATH_MSG_DEBUG(" ACTS Extrapolation to perigee failed for input parameters: " << destParameters->parameters());
+ }
+
+ m_actsPropResultWriterSvc->write<Acts::BoundParameters>(startParameters, destParameters.release(), ms_fwd, finalperigee.release(), ms_bkw);
+ } else if (!destParameters) {
+ ATH_MSG_DEBUG(" ACTS Extrapolation not successful! " );
+ m_actsPropResultWriterSvc->write<Acts::BoundParameters>(startParameters);
+ }
+ }
+ delete startParameters;
+ }
+
+ end = xclock::now();
+ secs = std::chrono::duration_cast<std::chrono::milliseconds>(end-start).count() / milliseconds_to_seconds;
+ secs_per_ex = secs / n_extraps;
+
+ ATH_MSG_INFO("ACTS : Time for " << n_extraps << " iterations: " << secs << "s (" << secs_per_ex << "s per extrapolation)");
+
+ return StatusCode::SUCCESS;
+}
+
+//============================================================================================
+
+
+
diff --git a/Tracking/TrkExtrapolation/TrkExAlgs/src/PropResultRootWriterSvc.cxx b/Tracking/TrkExtrapolation/TrkExAlgs/src/PropResultRootWriterSvc.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..a9905d4a8f90cac4332ec40dbd4d8bdad4f67b3b
--- /dev/null
+++ b/Tracking/TrkExtrapolation/TrkExAlgs/src/PropResultRootWriterSvc.cxx
@@ -0,0 +1,84 @@
+
+/*
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "TrkExAlgs/PropResultRootWriterSvc.h"
+
+#include <vector>
+#include <deque>
+#include <mutex>
+#include <thread>
+
+#include "TFile.h"
+
+Trk::PropResultRootWriterSvc::PropResultRootWriterSvc( const std::string& name, ISvcLocator* svc )
+: AthService(name, svc),
+ m_thistSvc("THistSvc", name),
+ m_tree(0) {
+}
+
+Trk::PropResultRootWriterSvc::~PropResultRootWriterSvc() {}
+
+StatusCode Trk::PropResultRootWriterSvc::initialize()
+{
+ ATH_CHECK( AthService::initialize() );
+
+ CHECK(m_thistSvc.retrieve());
+
+ const std::string treeName = m_treeName;
+ m_tree = new TTree(treeName.data(), "MyStudies");
+
+ std::string fullNtupleName = m_ntupleDirName + m_treeName;
+ CHECK(m_thistSvc->regTree(fullNtupleName, m_tree));
+
+ if(m_tree == nullptr) {
+ ATH_MSG_ERROR("Unable to create TTree");
+ return StatusCode::FAILURE;
+ }
+
+ m_tree->Branch("event_nr" , &m_eventNum );
+
+ m_tree->Branch("start_d0" , &m_start_d0 );
+ m_tree->Branch("start_z0" , &m_start_z0 );
+ m_tree->Branch("start_phi" , &m_start_phi );
+ m_tree->Branch("start_theta", &m_start_theta);
+ m_tree->Branch("start_qop" , &m_start_qop );
+
+ m_tree->Branch("fwd_success" , &m_fwd_success );
+ m_tree->Branch("fwd_time" , &m_fwd_time );
+ m_tree->Branch("fwd_final_l0" , &m_fwd_final_l0 );
+ m_tree->Branch("fwd_final_l1" , &m_fwd_final_l1 );
+ m_tree->Branch("fwd_final_x" , &m_fwd_final_x );
+ m_tree->Branch("fwd_final_y" , &m_fwd_final_y );
+ m_tree->Branch("fwd_final_z" , &m_fwd_final_z );
+ m_tree->Branch("fwd_final_phi" , &m_fwd_final_phi );
+ m_tree->Branch("fwd_final_theta", &m_fwd_final_theta);
+ m_tree->Branch("fwd_final_qop" , &m_fwd_final_qop );
+ m_tree->Branch("fwd_final_sigma_l0" , &m_fwd_final_sigma_l0 );
+ m_tree->Branch("fwd_final_sigma_l1" , &m_fwd_final_sigma_l1 );
+ m_tree->Branch("fwd_final_sigma_phi" , &m_fwd_final_sigma_phi );
+ m_tree->Branch("fwd_final_sigma_theta", &m_fwd_final_sigma_theta);
+ m_tree->Branch("fwd_final_sigma_qop" , &m_fwd_final_sigma_qop );
+
+ m_tree->Branch("bkw_success" , &m_bkw_success );
+ m_tree->Branch("bkw_time" , &m_bkw_time );
+ m_tree->Branch("bkw_final_d0" , &m_bkw_final_d0 );
+ m_tree->Branch("bkw_final_z0" , &m_bkw_final_z0 );
+ m_tree->Branch("bkw_final_phi" , &m_bkw_final_phi );
+ m_tree->Branch("bkw_final_theta", &m_bkw_final_theta);
+ m_tree->Branch("bkw_final_qop" , &m_bkw_final_qop );
+ m_tree->Branch("bkw_final_sigma_d0" , &m_bkw_final_sigma_d0 );
+ m_tree->Branch("bkw_final_sigma_z0" , &m_bkw_final_sigma_z0 );
+ m_tree->Branch("bkw_final_sigma_phi" , &m_bkw_final_sigma_phi );
+ m_tree->Branch("bkw_final_sigma_theta", &m_bkw_final_sigma_theta);
+ m_tree->Branch("bkw_final_sigma_qop" , &m_bkw_final_sigma_qop );
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode Trk::PropResultRootWriterSvc::finalize()
+{
+ return StatusCode::SUCCESS;
+}
+
diff --git a/Tracking/TrkExtrapolation/TrkExAlgs/src/components/TrkExAlgs_entries.cxx b/Tracking/TrkExtrapolation/TrkExAlgs/src/components/TrkExAlgs_entries.cxx
index ee87e754693843c217a3ff82e1c8ece6df873da0..b87faee0acdd8b826f816beacf05717907e13487 100644
--- a/Tracking/TrkExtrapolation/TrkExAlgs/src/components/TrkExAlgs_entries.cxx
+++ b/Tracking/TrkExtrapolation/TrkExAlgs/src/components/TrkExAlgs_entries.cxx
@@ -1,16 +1,20 @@
#include "TrkExAlgs/ExtrapolatorTest.h"
+#include "TrkExAlgs/ExtrapolatorComparisonTest.h"
#include "TrkExAlgs/CombinedExtrapolatorTest.h"
#include "TrkExAlgs/CETmaterial.h"
#include "TrkExAlgs/ExtrapolationValidation.h"
#include "TrkExAlgs/EnergyLossExtrapolationValidation.h"
#include "TrkExAlgs/RiddersAlgorithm.h"
+#include "TrkExAlgs/PropResultRootWriterSvc.h"
+
using namespace Trk;
DECLARE_COMPONENT( RiddersAlgorithm )
DECLARE_COMPONENT( ExtrapolatorTest )
+DECLARE_COMPONENT( ExtrapolatorComparisonTest )
DECLARE_COMPONENT( ExtrapolationValidation )
DECLARE_COMPONENT( EnergyLossExtrapolationValidation )
DECLARE_COMPONENT( CombinedExtrapolatorTest )
DECLARE_COMPONENT( CETmaterial )
-
+DECLARE_COMPONENT( PropResultRootWriterSvc )