diff --git a/.gitignore b/.gitignore
index 43c694310ebec2a7a87f52c645ea16d1feec7c28..ce66b6b9f0bb6191bfc40929b521f9e067e1c040 100644
--- a/.gitignore
+++ b/.gitignore
@@ -39,3 +39,6 @@ __pycache__/
# ROOT file
*.root
+*.sif
+**old
+install/*
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 966af4005239f58a34be7237bb79357c29ccb1f8..9c9515f7b228738620816e34dbd9dbe2bf269614 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -1,5 +1,5 @@
variables:
- EIC_TAG: "2.5-stable"
+ EIC_TAG: "2.6-stable"
default:
image: eicweb.phy.anl.gov:4567/containers/image_recipes/ubuntu_dind:latest
@@ -29,6 +29,10 @@ env:
fi
- echo "CI JUGGLER_TAG for this pipeline set to ${JUGGLER_TAG}"
- echo "JUGGLER_TAG=$JUGGLER_TAG" >> juggler.env
+ - |
+ if [ "${DETECTOR_TAG}" ]; then
+ echo "DETECTOR_TAG=$DETECTOR_TAG" >> juggler.env
+ fi
artifacts:
reports:
dotenv: juggler.env
diff --git a/JugBase/CMakeLists.txt b/JugBase/CMakeLists.txt
index 965766bc1f7462c2b82030c7bdbff11c7447dacd..f199ed02d8629cf60f649d903bfe710da52297f9 100644
--- a/JugBase/CMakeLists.txt
+++ b/JugBase/CMakeLists.txt
@@ -21,6 +21,10 @@ gaudi_depends_on_subdirs(GaudiAlg GaudiKernel)
gaudi_add_library(JugBase
src/*.cpp
+ src/Utilities/*.cpp
+ src/Plugins/BFieldOptions.cpp
+ #src/Plugins/BFieldScalor.cpp
+ src/Plugins/BFieldUtils.cpp
INCLUDE_DIRS EICD PODIO ROOT $ENV{HOME}/stow/podio/include
LINK_LIBRARIES GaudiAlgLib GaudiKernel ROOT DD4hep::DDG4IO
PUBLIC_HEADERS JugBase)
diff --git a/JugBase/JugBase/BField/BFieldOptions.hpp b/JugBase/JugBase/BField/BFieldOptions.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8e303da275ec51a97b5ba0648776263a3fe9a770
--- /dev/null
+++ b/JugBase/JugBase/BField/BFieldOptions.hpp
@@ -0,0 +1,67 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2017 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Utilities/detail/AxisFwd.hpp"
+#include "Acts/Utilities/detail/GridFwd.hpp"
+#include "JugBase/Utilities/OptionsFwd.hpp"
+
+#include <memory>
+#include <tuple>
+#include <variant>
+
+// Forward declarations
+namespace Acts {
+template <typename G>
+struct InterpolatedBFieldMapper;
+
+template <typename M>
+class InterpolatedBFieldMap;
+
+class ConstantBField;
+} // namespace Acts
+
+namespace Jug {
+namespace BField {
+class ScalableBField;
+}
+} // namespace Jug
+
+using InterpolatedMapper2D = Acts::InterpolatedBFieldMapper<
+ Acts::detail::Grid<Acts::Vector2, Acts::detail::EquidistantAxis,
+ Acts::detail::EquidistantAxis>>;
+
+using InterpolatedMapper3D = Acts::InterpolatedBFieldMapper<Acts::detail::Grid<
+ Acts::Vector3, Acts::detail::EquidistantAxis, Acts::detail::EquidistantAxis,
+ Acts::detail::EquidistantAxis>>;
+
+using InterpolatedBFieldMap2D =
+ Acts::InterpolatedBFieldMap<InterpolatedMapper2D>;
+using InterpolatedBFieldMap3D =
+ Acts::InterpolatedBFieldMap<InterpolatedMapper3D>;
+
+namespace Jug {
+
+namespace Options {
+
+using BFieldVariant =
+ std::variant<std::shared_ptr<InterpolatedBFieldMap2D>,
+ std::shared_ptr<InterpolatedBFieldMap3D>,
+ std::shared_ptr<Acts::ConstantBField>,
+ std::shared_ptr<Jug::BField::ScalableBField>>;
+
+// common bfield options, with a bf prefix
+void addBFieldOptions(boost::program_options::options_description& opt);
+
+// create the bfield maps
+BFieldVariant readBField(const boost::program_options::variables_map& vm);
+
+} // namespace Options
+} // namespace Jug
diff --git a/JugBase/JugBase/BField/BFieldScalor.hpp b/JugBase/JugBase/BField/BFieldScalor.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..90462d1990d87fb17d4bbb9abd4f120ca4e369b4
--- /dev/null
+++ b/JugBase/JugBase/BField/BFieldScalor.hpp
@@ -0,0 +1,66 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2019 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Geometry/GeometryContext.hpp"
+#include "Acts/Utilities/Logger.hpp"
+#include "Jug/Framework/AlgorithmContext.hpp"
+#include "Jug/Framework/IContextDecorator.hpp"
+
+#include <vector>
+
+namespace Jug {
+
+namespace BField {
+
+/// A mockup service that rotates a
+/// cylindrical geometry
+class BFieldScalor : public IContextDecorator {
+ public:
+ /// @brief nested configuration struct
+ struct Config {
+ /// Incremental scaling
+ double scalor = 1.2;
+ };
+
+ /// Constructor
+ ///
+ /// @param cfg Configuration struct
+ /// @param logger The logging framework
+ BFieldScalor(const Config& cfg,
+ std::unique_ptr<const Acts::Logger> logger =
+ Acts::getDefaultLogger("BFieldScalor", Acts::Logging::INFO));
+
+ /// Virtual destructor
+ virtual ~BFieldScalor() = default;
+
+ /// @brief decorates (adds, modifies) the AlgorithmContext
+ /// with a geometric rotation per event
+ ///
+ /// @note If decorators depend on each other, they have to be
+ /// added in order.
+ ///
+ /// @param context the bare (or at least non-const) Event context
+ ProcessCode decorate(AlgorithmContext& context) final override;
+
+ /// @brief decorator name() for screen output
+ const std::string& name() const final override { return m_name; }
+
+ private:
+ Config m_cfg; ///< the configuration class
+ std::unique_ptr<const Acts::Logger> m_logger; ///!< the logging instance
+ std::string m_name = "BFieldScalor";
+
+ /// Private access to the logging instance
+ const Acts::Logger& logger() const { return *m_logger; }
+};
+} // namespace BField
+
+} // namespace Jug
diff --git a/JugBase/JugBase/BField/BFieldUtils.hpp b/JugBase/JugBase/BField/BFieldUtils.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..48376e8df4384c76f2b41c192bcd69c850d87f3b
--- /dev/null
+++ b/JugBase/JugBase/BField/BFieldUtils.hpp
@@ -0,0 +1,225 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2017 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/Definitions/Units.hpp"
+#include "Acts/MagneticField/InterpolatedBFieldMap.hpp"
+#include "Acts/Utilities/Logger.hpp"
+#include "Acts/Utilities/detail/Axis.hpp"
+#include "Acts/Utilities/detail/Grid.hpp"
+
+namespace Jug {
+
+namespace BField {
+
+namespace txt {
+/// Method to setup the FieldMapper
+/// @param localToGlobalBin Function mapping the local bins of r,z to the
+/// global
+/// bin of the map magnetic field value e.g.: we have small grid with the
+/// values: r={2,3}, z ={4,5}, the corresponding indices are i(r) and j(z),
+/// the
+/// globalIndex is M and the field map is:
+///|| r | i || z | j || |B(r,z)| || M ||
+/// -----------------------------------
+///|| 2 | 0 || 4 | 0 || 2.323 || 0 ||
+///|| 2 | 0 || 5 | 1 || 2.334 || 1 ||
+///|| 3 | 1 || 4 | 0 || 2.325 || 2 ||
+///|| 3 | 1 || 5 | 1 || 2.331 || 3 ||
+///
+/// @code
+/// In this case the function would look like:
+/// [](std::array<size_t, 2> binsRZ, std::array<size_t, 2> nBinsRZ) {
+/// return (binsRZ.at(0) * nBinsRZ.at(1) + binsRZ.at(1));
+/// }
+/// @endcode
+/// @param[in] fieldMapFile Path to file containing field map in txt format
+/// @param[in] lengthUnit The unit of the grid points
+/// @param[in] BFieldUnit The unit of the magnetic field
+/// @param[in] nPoints Extimate of number of grid points in field map needed
+/// for allocation
+/// @note This information is only used as a hint for the required size of
+/// the internal vectors. A correct value is not needed, but will help
+/// to speed up the field map initialization process.
+/// @param[in] firstOctant Flag if set to true indicating that only the
+/// first
+/// quadrant of the grid points and the BField values has been given and
+/// that
+/// the BFieldMap should be created symmetrically for all quadrants.
+/// e.g. we have the grid values r={0,1} with BFieldValues={2,3} on the r
+/// axis.
+/// If the flag is set to true the r-axis grid values will be set to
+/// {-1,0,1}
+/// and the BFieldValues will be set to {3,2,3}.
+Acts::InterpolatedBFieldMapper<
+ Acts::detail::Grid<Acts::Vector2, Acts::detail::EquidistantAxis,
+ Acts::detail::EquidistantAxis>>
+fieldMapperRZ(std::function<size_t(std::array<size_t, 2> binsRZ,
+ std::array<size_t, 2> nBinsRZ)>
+ localToGlobalBin,
+ std::string fieldMapFile = "",
+ double lengthUnit = Acts::UnitConstants::mm,
+ double BFieldUnit = Acts::UnitConstants::T, size_t nPoints = 1000,
+ bool firstOctant = false);
+
+/// Method to setup the FieldMapper
+/// @param localToGlobalBin Function mapping the local bins of x,y,z to the
+/// global bin of the map magnetic field value e.g.: we have small grid with
+/// the
+/// values: x={2,3}, y={3,4}, z ={4,5}, the corresponding indices are i(x),
+/// j(y)
+/// and z(k), the globalIndex is M and the field map is:
+///|| x | i || y | j || z | k || |B(x,y,z)| || M ||
+/// --------------------------------------------
+///|| 2 | 0 || 3 | 0 || 4 | 0 || 2.323 || 0 ||
+///|| 2 | 0 || 3 | 0 || 5 | 1 || 2.334 || 1 ||
+///|| 2 | 0 || 4 | 1 || 4 | 0 || 2.325 || 2 ||
+///|| 2 | 0 || 4 | 1 || 5 | 1 || 2.331 || 3 ||
+///|| 3 | 1 || 3 | 0 || 4 | 0 || 2.323 || 4 ||
+///|| 3 | 1 || 3 | 0 || 5 | 1 || 2.334 || 5 ||
+///|| 3 | 1 || 4 | 1 || 4 | 0 || 2.325 || 6 ||
+///|| 3 | 1 || 4 | 1 || 5 | 1 || 2.331 || 7 ||
+///
+/// @code
+/// In this case the function would look like:
+/// [](std::array<size_t, 3> binsXYZ, std::array<size_t, 3> nBinsXYZ) {
+/// return (binsXYZ.at(0) * (nBinsXYZ.at(1) * nBinsXYZ.at(2))
+/// + binsXYZ.at(1) * nBinsXYZ.at(2)
+/// + binsXYZ.at(2));
+/// }
+/// @endcode
+/// @param[in] fieldMapFile Path to file containing field map in txt format
+/// @param[in] lengthUnit The unit of the grid points
+/// @param[in] BFieldUnit The unit of the magnetic field
+/// @param[in] nPoints Extimate of number of grid points in field map needed
+/// for allocation
+/// @note This information is only used as a hint for the required size of
+/// the internal vectors. A correct value is not needed, but will help
+/// to speed up the field map initialization process.
+/// @param[in] firstOctant Flag if set to true indicating that only the
+/// first
+/// octant of the grid points and the BField values has been given and that
+/// the BFieldMap should be created symmetrically for all quadrants.
+/// e.g. we have the grid values z={0,1} with BFieldValues={2,3} on the r
+/// axis.
+/// If the flag is set to true the z-axis grid values will be set to
+/// {-1,0,1}
+/// and the BFieldValues will be set to {3,2,3}.
+Acts::InterpolatedBFieldMapper<Acts::detail::Grid<
+ Acts::Vector3, Acts::detail::EquidistantAxis, Acts::detail::EquidistantAxis,
+ Acts::detail::EquidistantAxis>>
+fieldMapperXYZ(std::function<size_t(std::array<size_t, 3> binsXYZ,
+ std::array<size_t, 3> nBinsXYZ)>
+ localToGlobalBin,
+ std::string fieldMapFile = "",
+ double lengthUnit = Acts::UnitConstants::mm,
+ double BFieldUnit = Acts::UnitConstants::T,
+ size_t nPoints = 1000, bool firstOctant = false);
+} // namespace txt
+
+namespace root {
+/// Method to setup the FieldMapper
+/// @param localToGlobalBin Function mapping the local bins of r,z to the
+/// global
+/// bin of the map magnetic field value e.g.: we have small grid with the
+/// values: r={2,3}, z ={4,5}, the corresponding indices are i(r) and j(z),
+/// the
+/// globalIndex is M and the field map is:
+///|| r | i || z | j || |B(r,z)| || M ||
+/// -----------------------------------
+///|| 2 | 0 || 4 | 0 || 2.323 || 0 ||
+///|| 2 | 0 || 5 | 1 || 2.334 || 1 ||
+///|| 3 | 1 || 4 | 0 || 2.325 || 2 ||
+///|| 3 | 1 || 5 | 1 || 2.331 || 3 ||
+///
+/// @code
+/// In this case the function would look like:
+/// [](std::array<size_t, 2> binsRZ, std::array<size_t, 2> nBinsRZ) {
+/// return (binsRZ.at(0) * nBinsRZ.at(1) + binsRZ.at(1));
+/// }
+/// @endcode
+/// @param[in] fieldMapFile Path to file containing field map in txt format
+/// @param[in] treeName The name of the root tree
+/// @param[in] lengthUnit The unit of the grid points
+/// @param[in] BFieldUnit The unit of the magnetic field
+/// @param[in] firstQuadrant Flag if set to true indicating that only the
+/// first
+/// quadrant of the grid points and the BField values has been given and
+/// that
+/// the BFieldMap should be created symmetrically for all quadrants.
+/// e.g. we have the grid values r={0,1} with BFieldValues={2,3} on the r
+/// axis.
+/// If the flag is set to true the r-axis grid values will be set to
+/// {-1,0,1}
+/// and the BFieldValues will be set to {3,2,3}.
+Acts::InterpolatedBFieldMapper<
+ Acts::detail::Grid<Acts::Vector2, Acts::detail::EquidistantAxis,
+ Acts::detail::EquidistantAxis>>
+fieldMapperRZ(std::function<size_t(std::array<size_t, 2> binsRZ,
+ std::array<size_t, 2> nBinsRZ)>
+ localToGlobalBin,
+ std::string fieldMapFile = "", std::string treeName = "",
+ double lengthUnit = Acts::UnitConstants::mm,
+ double BFieldUnit = Acts::UnitConstants::T,
+ bool firstOctant = false);
+
+/// Method to setup the FieldMapper
+/// @param localToGlobalBin Function mapping the local bins of x,y,z to the
+/// global bin of the map magnetic field value e.g.: we have small grid with
+/// the
+/// values: x={2,3}, y={3,4}, z ={4,5}, the corresponding indices are i(x),
+/// j(y)
+/// and z(k), the globalIndex is M and the field map is:
+///|| x | i || y | j || z | k || |B(x,y,z)| || M ||
+/// --------------------------------------------
+///|| 2 | 0 || 3 | 0 || 4 | 0 || 2.323 || 0 ||
+///|| 2 | 0 || 3 | 0 || 5 | 1 || 2.334 || 1 ||
+///|| 2 | 0 || 4 | 1 || 4 | 0 || 2.325 || 2 ||
+///|| 2 | 0 || 4 | 1 || 5 | 1 || 2.331 || 3 ||
+///|| 3 | 1 || 3 | 0 || 4 | 0 || 2.323 || 4 ||
+///|| 3 | 1 || 3 | 0 || 5 | 1 || 2.334 || 5 ||
+///|| 3 | 1 || 4 | 1 || 4 | 0 || 2.325 || 6 ||
+///|| 3 | 1 || 4 | 1 || 5 | 1 || 2.331 || 7 ||
+///
+/// @code
+/// In this case the function would look like:
+/// [](std::array<size_t, 3> binsXYZ, std::array<size_t, 3> nBinsXYZ) {
+/// return (binsXYZ.at(0) * (nBinsXYZ.at(1) * nBinsXYZ.at(2))
+/// + binsXYZ.at(1) * nBinsXYZ.at(2)
+/// + binsXYZ.at(2));
+/// }
+/// @endcode
+/// @param[in] fieldMapFile Path to file containing field map in txt format
+/// @param[in] treeName The name of the root tree
+/// @param[in] lengthUnit The unit of the grid points
+/// @param[in] BFieldUnit The unit of the magnetic field
+/// @param[in] firstOctant Flag if set to true indicating that only the
+/// first
+/// octant of the grid points and the BField values has been given and that
+/// the BFieldMap should be created symmetrically for all quadrants.
+/// e.g. we have the grid values z={0,1} with BFieldValues={2,3} on the r
+/// axis.
+/// If the flag is set to true the z-axis grid values will be set to
+/// {-1,0,1}
+/// and the BFieldValues will be set to {3,2,3}.
+Acts::InterpolatedBFieldMapper<Acts::detail::Grid<
+ Acts::Vector3, Acts::detail::EquidistantAxis, Acts::detail::EquidistantAxis,
+ Acts::detail::EquidistantAxis>>
+fieldMapperXYZ(std::function<size_t(std::array<size_t, 3> binsXYZ,
+ std::array<size_t, 3> nBinsXYZ)>
+ localToGlobalBin,
+ std::string fieldMapFile = "", std::string treeName = "",
+ double lengthUnit = Acts::UnitConstants::mm,
+ double BFieldUnit = Acts::UnitConstants::T,
+ bool firstOctant = false);
+} // namespace root
+
+} // namespace BField
+
+} // namespace Jug
diff --git a/JugBase/JugBase/BField/ScalableBField.hpp b/JugBase/JugBase/BField/ScalableBField.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8307a5a4c9a3cd9d9c0f53a247db4631064a2898
--- /dev/null
+++ b/JugBase/JugBase/BField/ScalableBField.hpp
@@ -0,0 +1,137 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/Definitions/Algebra.hpp"
+#include "Acts/MagneticField/MagneticFieldContext.hpp"
+
+namespace Jug {
+
+namespace BField {
+
+/// The Context to be handed around
+struct ScalableBFieldContext {
+ double scalor = 1.;
+};
+
+/// @ingroup MagneticField
+///
+/// @brief returns a given constant field value at every point
+///
+/// This class is based on the constant magnetic field class
+/// but allows a event based context
+class ScalableBField final {
+ public:
+ struct Cache {
+ double scalor = 1.;
+
+ /// @brief constructor with context
+ Cache(const Acts::MagneticFieldContext& mcfg) {
+ scalor = std::any_cast<const ScalableBFieldContext>(mcfg).scalor;
+ }
+ };
+
+ /// @brief construct constant magnetic field from field vector
+ ///
+ /// @param [in] B magnetic field vector in global coordinate system
+ explicit ScalableBField(Acts::Vector3 B) : m_BField(std::move(B)) {}
+
+ /// @brief construct constant magnetic field from components
+ ///
+ /// @param [in] Bx magnetic field component in global x-direction
+ /// @param [in] By magnetic field component in global y-direction
+ /// @param [in] Bz magnetic field component in global z-direction
+ ScalableBField(double Bx = 0., double By = 0., double Bz = 0.)
+ : m_BField(Bx, By, Bz) {}
+
+ /// @brief retrieve magnetic field value
+ ///
+ /// @param [in] position global position
+ /// @return magnetic field vector
+ ///
+ /// @note The @p position is ignored and only kept as argument to provide
+ /// a consistent interface with other magnetic field services.
+ Acts::Vector3 getField(const Acts::Vector3& /*position*/) const {
+ return m_BField;
+ }
+
+ /// @brief retrieve magnetic field value
+ ///
+ /// @param [in] position global position
+ /// @param [in] cache Cache object (is ignored)
+ /// @return magnetic field vector
+ ///
+ /// @note The @p position is ignored and only kept as argument to provide
+ /// a consistent interface with other magnetic field services.
+ Acts::Vector3 getField(const Acts::Vector3& /*position*/,
+ Cache& cache) const {
+ return m_BField * cache.scalor;
+ }
+
+ /// @brief retrieve magnetic field value & its gradient
+ ///
+ /// @param [in] position global position
+ /// @param [out] derivative gradient of magnetic field vector as (3x3)
+ /// matrix
+ /// @return magnetic field vector
+ ///
+ /// @note The @p position is ignored and only kept as argument to provide
+ /// a consistent interface with other magnetic field services.
+ /// @note currently the derivative is not calculated
+ /// @todo return derivative
+ Acts::Vector3 getFieldGradient(const Acts::Vector3& /*position*/,
+ Acts::ActsMatrix<3, 3>& /*derivative*/) const {
+ return m_BField;
+ }
+
+ /// @brief retrieve magnetic field value & its gradient
+ ///
+ /// @param [in] position global position
+ /// @param [out] derivative gradient of magnetic field vector as (3x3)
+ /// matrix
+ /// @param [in] cache Cache object (is ignored)
+ /// @return magnetic field vector
+ ///
+ /// @note The @p position is ignored and only kept as argument to provide
+ /// a consistent interface with other magnetic field services.
+ /// @note currently the derivative is not calculated
+ /// @todo return derivative
+ Acts::Vector3 getFieldGradient(const Acts::Vector3& /*position*/,
+ Acts::ActsMatrix<3, 3>& /*derivative*/,
+ Cache& cache) const {
+ return m_BField * cache.scalor;
+ }
+
+ /// @brief check whether given 3D position is inside look-up domain
+ ///
+ /// @param [in] position global 3D position
+ /// @return @c true if position is inside the defined look-up grid,
+ /// otherwise @c false
+ /// @note The method will always return true for the constant B-Field
+ bool isInside(const Acts::Vector3& /*position*/) const { return true; }
+
+ /// @brief update magnetic field vector from components
+ ///
+ /// @param [in] Bx magnetic field component in global x-direction
+ /// @param [in] By magnetic field component in global y-direction
+ /// @param [in] Bz magnetic field component in global z-direction
+ void setField(double Bx, double By, double Bz) { m_BField << Bx, By, Bz; }
+
+ /// @brief update magnetic field vector
+ ///
+ /// @param [in] B magnetic field vector in global coordinate system
+ void setField(const Acts::Vector3& B) { m_BField = B; }
+
+ private:
+ /// magnetic field vector
+ Acts::Vector3 m_BField;
+};
+
+} // namespace BField
+} // namespace Jug
diff --git a/JugBase/JugBase/Utilities/GroupBy.hpp b/JugBase/JugBase/Utilities/GroupBy.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..399e18bf8a8d881420b33a0548df6794591fee9d
--- /dev/null
+++ b/JugBase/JugBase/Utilities/GroupBy.hpp
@@ -0,0 +1,141 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2020 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "JugBase/Utilities/Range.hpp"
+
+#include <algorithm>
+#include <iterator>
+#include <utility>
+
+namespace Jug {
+
+/// Proxy for iterating over groups of elements within a container.
+///
+/// @note Each group will contain at least one element.
+///
+/// Consecutive elements with the same key (as defined by the KeyGetter) are
+/// placed in one group. The proxy should always be used as part of a
+/// range-based for loop. In combination with structured bindings to reduce the
+/// boilerplate, the group iteration can be written as
+///
+/// for (auto&& [key, elements] : GroupBy<...>(...)) {
+/// // do something with just the key
+/// ...
+///
+/// // iterate over the group elements
+/// for (const auto& element : elements) {
+/// ...
+/// }
+/// }
+///
+template <typename Iterator, typename KeyGetter>
+class GroupBy {
+ public:
+ /// The key type that identifies elements within a group.
+ using Key = std::decay_t<decltype(KeyGetter()(*Iterator()))>;
+ /// A Group is an iterator range with the associated key.
+ using Group = std::pair<Key, Range<Iterator>>;
+ /// Iterator type representing the end of the groups.
+ ///
+ /// The end iterator will not be dereferenced in C++17 range-based loops. It
+ /// can thus be a simpler type without the overhead of the full group iterator
+ /// below.
+ using GroupEndIterator = Iterator;
+ /// Iterator type representing a group of elements.
+ class GroupIterator {
+ public:
+ using iterator_category = std::input_iterator_tag;
+ using value_type = Group;
+ using difference_type = std::ptrdiff_t;
+ using pointer = Group*;
+ using reference = Group&;
+
+ constexpr GroupIterator(const GroupBy& groupBy, Iterator groupBegin)
+ : m_groupBy(groupBy),
+ m_groupBegin(groupBegin),
+ m_groupEnd(groupBy.findEndOfGroup(groupBegin)) {}
+ /// Pre-increment operator to advance to the next group.
+ constexpr GroupIterator& operator++() {
+ // make the current end the new group beginning
+ std::swap(m_groupBegin, m_groupEnd);
+ // find the end of the next group starting from the new beginning
+ m_groupEnd = m_groupBy.findEndOfGroup(m_groupBegin);
+ return *this;
+ }
+ /// Post-increment operator to advance to the next group.
+ constexpr GroupIterator operator++(int) {
+ GroupIterator retval = *this;
+ ++(*this);
+ return retval;
+ }
+ /// Derefence operator that returns the pointed-to group of elements.
+ constexpr Group operator*() const {
+ const Key key = (m_groupBegin != m_groupEnd)
+ ? m_groupBy.m_keyGetter(*m_groupBegin)
+ : Key();
+ return {key, makeRange(m_groupBegin, m_groupEnd)};
+ }
+
+ private:
+ const GroupBy& m_groupBy;
+ Iterator m_groupBegin;
+ Iterator m_groupEnd;
+
+ friend constexpr bool operator==(const GroupIterator& lhs,
+ const GroupEndIterator& rhs) {
+ return lhs.m_groupBegin == rhs;
+ }
+ friend constexpr bool operator!=(const GroupIterator& lhs,
+ const GroupEndIterator& rhs) {
+ return not(lhs == rhs);
+ }
+ };
+
+ /// Construct the group-by proxy for an iterator range.
+ constexpr GroupBy(Iterator begin, Iterator end,
+ KeyGetter keyGetter = KeyGetter())
+ : m_begin(begin), m_end(end), m_keyGetter(std::move(keyGetter)) {}
+ constexpr GroupIterator begin() const {
+ return GroupIterator(*this, m_begin);
+ }
+ constexpr GroupEndIterator end() const { return m_end; }
+ constexpr bool empty() const { return m_begin == m_end; }
+
+ private:
+ Iterator m_begin;
+ Iterator m_end;
+ KeyGetter m_keyGetter;
+
+ /// Find the end of the group that starts at the given position.
+ ///
+ /// This uses a linear search from the start position and thus has linear
+ /// complexity in the group size. It does not assume any ordering of the
+ /// underlying container and is a cache-friendly access pattern.
+ constexpr Iterator findEndOfGroup(Iterator start) const {
+ // check for end so we can safely dereference the start iterator.
+ if (start == m_end) {
+ return start;
+ }
+ // search the first element that does not share a key with the start.
+ return std::find_if_not(std::next(start), m_end,
+ [this, start](const auto& x) {
+ return m_keyGetter(x) == m_keyGetter(*start);
+ });
+ }
+};
+
+/// Construct the group-by proxy for a container.
+template <typename Container, typename KeyGetter>
+auto makeGroupBy(const Container& container, KeyGetter keyGetter)
+ -> GroupBy<decltype(std::begin(container)), KeyGetter> {
+ return {std::begin(container), std::end(container), std::move(keyGetter)};
+}
+
+} // namespace Jug
diff --git a/JugBase/JugBase/Utilities/Helpers.hpp b/JugBase/JugBase/Utilities/Helpers.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..890bb25f4adae956c993d93b0bfe3640d20ea566
--- /dev/null
+++ b/JugBase/JugBase/Utilities/Helpers.hpp
@@ -0,0 +1,126 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2019 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include <string>
+
+#include "TEfficiency.h"
+#include "TFitResult.h"
+#include "TFitResultPtr.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TProfile.h"
+#include "TROOT.h"
+
+namespace Jug {
+
+namespace PlotHelpers {
+/// @brief Nested binning struct for booking plots
+struct Binning {
+ Binning(){};
+
+ Binning(std::string bTitle, int bins, float bMin, float bMax)
+ : title(bTitle), nBins(bins), min(bMin), max(bMax){};
+
+ std::string title; ///< title to be displayed
+ int nBins; ///< number of bins
+ float min; ///< minimum value
+ float max; ///< maximum value
+};
+
+/// @brief book a 1D histogram
+/// @param histName the name of histogram
+/// @param histTitle the title of histogram
+/// @param varBinning the binning info of variable
+/// @return histogram pointer
+TH1F* bookHisto(const char* histName, const char* histTitle,
+ const Binning& varBinning);
+
+/// @brief book a 2D histogram
+/// @param histName the name of histogram
+/// @param histTitle the title of histogram
+/// @param varXBinning the binning info of variable at x axis
+/// @param varYBinning the binning info of variable at y axis
+/// @return histogram pointer
+TH2F* bookHisto(const char* histName, const char* histTitle,
+ const Binning& varXBinning, const Binning& varYBinning);
+
+/// @brief fill a 1D histogram
+/// @param hist histogram to fill
+/// @param value value to fill
+/// @param weight weight to fill
+void fillHisto(TH1F* hist, float value, float weight = 1.0);
+
+/// @brief fill a 2D histogram
+/// @param hist histogram to fill
+/// @param xValue x value to fill
+/// @param yValue y value to fill
+/// @param weight weight to fill
+void fillHisto(TH2F* hist, float xValue, float yValue, float weight = 1.0);
+
+/// @brief extract details, i.e. mean and width of a 1D histogram and fill
+/// them into histograms
+/// @param inputHist histogram to investigate
+/// @param j which bin number of meanHist and widthHist to fill
+/// @param meanHist histogram to fill the mean value of inputHist
+/// @param widthHist histogram to fill the width value of inputHist
+///
+/// @todo write specialized helper class to extract details of hists
+void anaHisto(TH1D* inputHist, int j, TH1F* meanHist, TH1F* widthHist);
+
+/// @brief book a 1D efficiency plot
+/// @param effName the name of plot
+/// @param effTitle the title of plot
+/// @param varBinning the binning info of variable
+/// @return TEfficiency pointer
+TEfficiency* bookEff(const char* effName, const char* effTitle,
+ const Binning& varBinning);
+
+/// @brief book a 2D efficiency plot
+/// @param effName the name of plot
+/// @param effTitle the title of plot
+/// @param varXBinning the binning info of variable at x axis
+/// @param varYBinning the binning info of variable at y axis
+/// @return TEfficiency pointer
+TEfficiency* bookEff(const char* effName, const char* effTitle,
+ const Binning& varXBinning, const Binning& varYBinning);
+
+/// @brief fill a 1D efficiency plot
+/// @param efficiency plot to fill
+/// @param value value to fill
+/// @param status bool to denote passed or not
+void fillEff(TEfficiency* efficiency, float value, bool status);
+
+/// @brief fill a 2D efficiency plot
+/// @param efficiency plot to fill
+/// @param xValue x value to fill
+/// @param yValue y value to fill
+/// @param status bool to denote passed or not
+void fillEff(TEfficiency* efficiency, float xValue, float yValue, bool status);
+
+/// @brief book a TProfile plot
+/// @param profName the name of plot
+/// @param profTitle the title of plot
+/// @param varXBinning the binning info of variable at x axis
+/// @param varYBinning the binning info of variable at y axis
+/// @return TProfile pointer
+TProfile* bookProf(const char* profName, const char* profTitle,
+ const Binning& varXBinning, const Binning& varYBinning);
+
+/// @brief fill a TProfile plot
+/// @param profile plot to fill
+/// @param xValue xvalue to fill
+/// @param yValue yvalue to fill
+/// @param weight weight to fill
+void fillProf(TProfile* profile, float xValue, float yValue,
+ float weight = 1.0);
+
+} // namespace PlotHelpers
+
+} // namespace Jug
diff --git a/JugBase/JugBase/Utilities/Options.hpp b/JugBase/JugBase/Utilities/Options.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..c927f2d37f4a68e0bb2acb5d460d09f6bc8dd6b9
--- /dev/null
+++ b/JugBase/JugBase/Utilities/Options.hpp
@@ -0,0 +1,170 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2017-2020 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include <array>
+#include <iosfwd>
+#include <optional>
+#include <vector>
+
+namespace Jug {
+namespace Options {
+
+/// @defgroup option-types Additional types for program options
+///
+/// All types are intended as utility type for the user options and not as a
+/// variable type for the configuration structs. They should only be used where
+/// a single option can not be represented by an existing primitive types.
+///
+/// They also must be distinct types and can not just be typedefs; otherwise we
+/// can not define the required operator{<<,>>} overloads in this namespace.
+///
+/// @{
+
+/// Half open [lower,upper) interval type for a single user option.
+///
+/// A missing limit represents an unbounded upper or lower limit. With just
+/// one defined limit the interval is just a lower/upper bound; with both
+/// limits undefined, the interval is unbounded everywhere and thus contains
+/// all possible values.
+struct Interval {
+ std::optional<double> lower;
+ std::optional<double> upper;
+};
+
+/// A fixed number of real values as one user option.
+///
+/// @note Implemented as a subclass so it is distinct from `std::array`
+/// and we can provide overloads in the same namespace.
+template <size_t kSize>
+class Reals : public std::array<double, kSize> {};
+
+/// An arbitrary number of revaluesal as one user option.
+///
+/// @note Making this a `std::vector<double>` typedef or subclass confuses
+/// program options, since `std::vector<double>` is interpreted as a `double`
+/// option that can be provided multiple times.
+struct VariableReals {
+ std::vector<double> values;
+};
+
+/// A fixed number of integers as one user option.
+///
+/// @note Implemented as a subclass so it is distinct from `std::array`
+/// and we can provide overloads in the same namespace.
+template <size_t kSize>
+class Integers : public std::array<int, kSize> {};
+
+/// An arbitrary number of integers as one user option.
+///
+/// @note Making this a `std::vector<int>` typedef or subclass confuses
+/// program options, since `std::vector<int>` is interpreted as an `int`
+/// option that can be provided multiple times.
+struct VariableIntegers {
+ std::vector<int> values;
+};
+
+/// @}
+
+/// Extract an interval from an input of the form 'lower:upper'.
+///
+/// An input of the form `lower:` or `:upper` sets just one of the limits. Any
+/// other input leads to an unbounded interval.
+///
+/// @note The more common range notation uses `lower-upper` but the `-`
+/// separator complicates the parsing of negative values.
+std::istream& operator>>(std::istream& is, Interval& interval);
+
+/// Print an interval as `lower:upper`.
+std::ostream& operator<<(std::ostream& os, const Interval& interval);
+
+namespace detail {
+void parseDoublesFixed(std::istream& is, size_t size, double* values);
+void parseDoublesVariable(std::istream& is, std::vector<double>& values);
+void printDoubles(std::ostream& os, size_t size, const double* values);
+} // namespace detail
+
+/// Extract a fixed number of doubles from an input of the form 'x:y:z'.
+///
+/// @note If the values would be separated by whitespace, negative values
+/// and additional command line both start with `-` and would be
+/// undistinguishable.
+template <size_t kSize>
+inline std::istream& operator>>(std::istream& is, Reals<kSize>& values) {
+ detail::parseDoublesFixed(is, kSize, values.data());
+ return is;
+}
+
+/// Extract a variable number of doubles from an input of the form 'x:y:...'.
+///
+/// @note If the values would be separated by whitespace, negative values
+/// and additional command line both start with `-` and would be
+/// undistinguishable.
+inline std::istream& operator>>(std::istream& is, VariableReals& values) {
+ detail::parseDoublesVariable(is, values.values);
+ return is;
+}
+
+/// Print a fixed number of doubles as `x:y:z`.
+template <size_t kSize>
+inline std::ostream& operator<<(std::ostream& os, const Reals<kSize>& values) {
+ detail::printDoubles(os, kSize, values.data());
+ return os;
+}
+
+/// Print a variable number of doubles as `x:y:z:...`.
+inline std::ostream& operator<<(std::ostream& os, const VariableReals& values) {
+ detail::printDoubles(os, values.values.size(), values.values.data());
+ return os;
+}
+
+namespace detail {
+void parseIntegersFixed(std::istream& is, size_t size, int* values);
+void parseIntegersVariable(std::istream& is, std::vector<int>& values);
+void printIntegers(std::ostream& os, size_t size, const int* values);
+} // namespace detail
+
+/// Extract a fixed number of integers from an input of the form 'x:y:z'.
+///
+/// @note If the values would be separated by whitespace, negative values
+/// and additional command line both start with `-` and would be
+/// undistinguishable.
+template <size_t kSize>
+inline std::istream& operator>>(std::istream& is, Integers<kSize>& values) {
+ detail::parseIntegersFixed(is, kSize, values.data());
+ return is;
+}
+
+/// Extract a variable number of integers from an input of the form 'x:y:...'.
+///
+/// @note If the values would be separated by whitespace, negative values
+/// and additional command line both start with `-` and would be
+/// undistinguishable.
+inline std::istream& operator>>(std::istream& is, VariableIntegers& values) {
+ detail::parseIntegersVariable(is, values.values);
+ return is;
+}
+
+/// Print a fixed number of integers as `x:y:z`.
+template <size_t kSize>
+inline std::ostream& operator<<(std::ostream& os,
+ const Integers<kSize>& values) {
+ detail::printIntegers(os, kSize, values.data());
+ return os;
+}
+
+/// Print a variable number of integers as `x:y:z:...`.
+inline std::ostream& operator<<(std::ostream& os,
+ const VariableIntegers& values) {
+ detail::printIntegers(os, values.values.size(), values.values.data());
+ return os;
+}
+
+} // namespace Options
+} // namespace Jug
diff --git a/JugBase/JugBase/Utilities/OptionsFwd.hpp b/JugBase/JugBase/Utilities/OptionsFwd.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8f03f638979427a366e1e8d98fe32f94bac28cc9
--- /dev/null
+++ b/JugBase/JugBase/Utilities/OptionsFwd.hpp
@@ -0,0 +1,23 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2019 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+namespace boost {
+namespace program_options {
+class options_description;
+class variables_map;
+} // namespace program_options
+} // namespace boost
+
+namespace Jug {
+namespace Options {
+using Description = ::boost::program_options::options_description;
+using Variables = ::boost::program_options::variables_map;
+} // namespace Options
+} // namespace Jug
diff --git a/JugBase/JugBase/Utilities/Paths.hpp b/JugBase/JugBase/Utilities/Paths.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..555975100d371ddd6d3dfc13d78be206eeea21f0
--- /dev/null
+++ b/JugBase/JugBase/Utilities/Paths.hpp
@@ -0,0 +1,46 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2017 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include <string>
+#include <utility>
+#include <vector>
+
+namespace Jug {
+
+/// Ensure that the given directory exists and is writable.
+///
+/// @return Canonical path to the directory.
+///
+/// Will create missing directories and throw on any error.
+std::string ensureWritableDirectory(const std::string& dir);
+
+/// Join dir and name into one path with correct handling of empty dirs.
+std::string joinPaths(const std::string& dir, const std::string& name);
+
+/// Construct a file path of the form `[<dir>/]event<XXXXXXXXX>-<name>`.
+///
+/// @params dir output directory, current directory if empty
+/// @params name basic filename
+/// @params event event number
+std::string perEventFilepath(const std::string& dir, const std::string& name,
+ size_t event);
+
+/// Determine the range of available events in a directory of per-event files.
+///
+/// @params dir input directory, current directory if empty
+/// @params name base filename
+/// @return first and last+1 event number
+/// @returns {0, 0} when no matching files could be found
+///
+/// Event files must be named `[<dir>/]event<XXXXXXXXX>-<name>` to be considered
+std::pair<size_t, size_t> determineEventFilesRange(const std::string& dir,
+ const std::string& name);
+
+} // namespace Jug
diff --git a/JugBase/JugBase/Utilities/Range.hpp b/JugBase/JugBase/Utilities/Range.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0cd3588adeb18d2f1525bf05879c783426e7d756
--- /dev/null
+++ b/JugBase/JugBase/Utilities/Range.hpp
@@ -0,0 +1,56 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2019 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include <iterator>
+#include <utility>
+
+namespace Jug {
+
+/// A wrapper around a pair of iterators to simplify range-based loops.
+///
+/// Some standard library algorithms return pairs of iterators to identify
+/// a sub-range. This wrapper simplifies the iteration and should be used as
+/// follows:
+///
+/// for (auto x : makeRange(std::equal_range(...)) {
+/// ...
+/// }
+///
+template <typename Iterator>
+class Range {
+ public:
+ Range(Iterator b, Iterator e) : m_begin(b), m_end(e) {}
+ Range(Range&&) = default;
+ Range(const Range&) = default;
+ ~Range() = default;
+ Range& operator=(Range&&) = default;
+ Range& operator=(const Range&) = default;
+
+ Iterator begin() const { return m_begin; }
+ Iterator end() const { return m_end; }
+ bool empty() const { return m_begin == m_end; }
+ std::size_t size() const { return std::distance(m_begin, m_end); }
+
+ private:
+ Iterator m_begin;
+ Iterator m_end;
+};
+
+template <typename Iterator>
+Range<Iterator> makeRange(Iterator begin, Iterator end) {
+ return Range<Iterator>(begin, end);
+}
+
+template <typename Iterator>
+Range<Iterator> makeRange(std::pair<Iterator, Iterator> range) {
+ return Range<Iterator>(range.first, range.second);
+}
+
+} // namespace Jug
diff --git a/JugBase/src/Plugins/BFieldOptions.cpp b/JugBase/src/Plugins/BFieldOptions.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ca0831427b568434d8205885a4936b3cf504ab4e
--- /dev/null
+++ b/JugBase/src/Plugins/BFieldOptions.cpp
@@ -0,0 +1,228 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2019 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "JugBase/BField/BFieldOptions.hpp"
+
+#include "Acts/Definitions/Units.hpp"
+#include "Acts/MagneticField/ConstantBField.hpp"
+#include "Acts/MagneticField/InterpolatedBFieldMap.hpp"
+#include "Acts/Utilities/Logger.hpp"
+#include "JugBase/BField/BFieldUtils.hpp"
+#include "JugBase/BField/ScalableBField.hpp"
+#include "JugBase/Utilities/Options.hpp"
+
+#include <iostream>
+#include <tuple>
+#include <utility>
+
+//#include <boost/program_options.hpp>
+
+//namespace po = boost::program_options;
+
+using InterpolatedMapper2D = Acts::InterpolatedBFieldMapper<
+ Acts::detail::Grid<Acts::Vector2, Acts::detail::EquidistantAxis,
+ Acts::detail::EquidistantAxis>>;
+
+using InterpolatedMapper3D = Acts::InterpolatedBFieldMapper<Acts::detail::Grid<
+ Acts::Vector3, Acts::detail::EquidistantAxis, Acts::detail::EquidistantAxis,
+ Acts::detail::EquidistantAxis>>;
+
+using InterpolatedBFieldMap2D =
+ Acts::InterpolatedBFieldMap<InterpolatedMapper2D>;
+using InterpolatedBFieldMap3D =
+ Acts::InterpolatedBFieldMap<InterpolatedMapper3D>;
+
+namespace Jug {
+
+//namespace Options {
+//
+//// common bfield options, with a bf prefix
+//void addBFieldOptions(boost::program_options::options_description& opt) {
+// opt.add_options()("bf-map", po::value<std::string>()->default_value(""),
+// "Set this string to point to the bfield source file."
+// "That can either be a '.txt', a '.csv' or a '.root' file. "
+// "Omit for a constant magnetic field.")(
+// "bf-name", po::value<std::string>()->default_value("bField"),
+// "In case your field map file is given in root format, please specify "
+// "the "
+// "name of the TTree.")(
+// "bf-gridpoints", po::value<size_t>()->default_value(100000),
+// "Estimate of number of grid points, "
+// "needed for allocation, only for txt and csv files.")(
+// "bf-lscalor", po::value<double>()->default_value(1.),
+// "The default unit for the grid "
+// "points is mm. In case the grid points of your field map has another "
+// "unit, please set the scalor to mm.")(
+// "bf-bscalor", po::value<double>()->default_value(1.),
+// "The default unit for the magnetic field values is Tesla. In case the "
+// "grid points of your field map has another unit, please set the "
+// "scalor "
+// "to [T].")(
+// "bf-rz", po::value<bool>()->default_value(false),
+// "Please set this flag to true, if your grid points and your "
+// "magnetic field values are given in 'rz'. The default is 'xyz'.")(
+// "bf-foctant", po::value<bool>()->default_value(false),
+// "Please set this flag to true, if your field map is only given for the "
+// "first octant/quadrant and should be symmetrically created for all "
+// "other "
+// "octants/quadrants.")(
+// "bf-values", po::value<Reals<3>>()->default_value({{0., 0., 0.}}),
+// "In case no magnetic field map is handed over. A constant magnetic "
+// "field will be created automatically. The values can be set with this "
+// "options. Please hand over the coordinates in cartesian coordinates: "
+// "{Bx,By,Bz} in Tesla.")(
+// "bf-context-scalable", po::value<bool>()->default_value(false),
+// "This is for testing the event dependent magnetic field scaling.");
+//}
+
+// create the bfield maps
+//BFieldVariant readBField(const boost::program_options::variables_map& vm) {
+// std::string bfieldmap = "constfield";
+//
+// enum BFieldMapType { constant = 0, root = 1, text = 2 };
+//
+// std::shared_ptr<InterpolatedBFieldMap2D> map2D = nullptr;
+// std::shared_ptr<InterpolatedBFieldMap3D> map3D = nullptr;
+// std::shared_ptr<Acts::ConstantBField> mapConst = nullptr;
+// std::shared_ptr<Jug::BField::ScalableBField> mapScale = nullptr;
+//
+// int bfieldmaptype = constant;
+// if (vm.count("bf-map") && vm["bf-map"].template as<std::string>() != "") {
+// bfieldmap = vm["bf-map"].template as<std::string>();
+// std::cout << "- read in magnetic field map: "
+// << vm["bf-map"].template as<std::string>() << std::endl;
+// if (bfieldmap.find(".root") != std::string::npos) {
+// std::cout << "- root format for magnetic field detected" << std::endl;
+// bfieldmaptype = root;
+// } else if (bfieldmap.find(".txt") != std::string::npos ||
+// bfieldmap.find(".csv") != std::string::npos) {
+// std::cout << "- txt format for magnetic field detected" << std::endl;
+// bfieldmaptype = text;
+// } else {
+// std::cout << "- magnetic field format could not be detected";
+// std::cout << " use '.root', '.txt', or '.csv'." << std::endl;
+// throw std::runtime_error("Invalid BField options");
+// }
+// }
+// if (bfieldmaptype == text && vm.count("bf-gridpoints")) {
+// std::cout << "- number of points set to: "
+// << vm["bf-gridpoints"].template as<size_t>() << std::endl;
+// }
+// double lscalor = 1.;
+// if (bfieldmaptype != constant && vm.count("bf-lscalor")) {
+// lscalor = vm["bf-lscalor"].template as<double>();
+// std::cout << "- length scalor to mm set to: " << lscalor << std::endl;
+// }
+// double bscalor = 1.;
+// if (vm.count("bf-bscalor")) {
+// bscalor = vm["bf-bscalor"].template as<double>();
+// std::cout << "- BField (scalor to/in) Tesla set to: " << bscalor
+// << std::endl;
+// }
+// if (bfieldmaptype != constant && vm["bf-rz"].template as<bool>())
+// std::cout << "- BField map is given in 'rz' coordiantes." << std::endl;
+// else if (bfieldmaptype != constant)
+// std::cout << "- BField map is given in 'xyz' coordiantes." << std::endl;
+//
+// if (bfieldmaptype != constant && vm["bf-foctant"].template as<bool>()) {
+// std::cout
+// << "- Only the first octant/quadrant is given, bField map will be "
+// "symmetrically created for all other octants/quadrants"
+// << std::endl;
+// }
+//
+// // Declare the mapper
+// double lengthUnit = lscalor * Acts::UnitConstants::mm;
+// double BFieldUnit = bscalor * Acts::UnitConstants::T;
+//
+// // set the mapper - foort
+// if (bfieldmaptype == root) {
+// if (vm["bf-rz"].template as<bool>()) {
+// auto mapper2D = Jug::BField::root::fieldMapperRZ(
+// [](std::array<size_t, 2> binsRZ, std::array<size_t, 2> nBinsRZ) {
+// return (binsRZ.at(1) * nBinsRZ.at(0) + binsRZ.at(0));
+// },
+// vm["bf-map"].template as<std::string>(),
+// vm["bf-name"].template as<std::string>(), lengthUnit, BFieldUnit,
+// vm["bf-foctant"].template as<bool>());
+//
+// // create field mapping
+// InterpolatedBFieldMap2D::Config config2D(std::move(mapper2D));
+// config2D.scale = bscalor;
+// // create BField service
+// return std::make_shared<InterpolatedBFieldMap2D>(std::move(config2D));
+//
+// } else {
+// auto mapper3D = Jug::BField::root::fieldMapperXYZ(
+// [](std::array<size_t, 3> binsXYZ, std::array<size_t, 3> nBinsXYZ) {
+// return (binsXYZ.at(0) * (nBinsXYZ.at(1) * nBinsXYZ.at(2)) +
+// binsXYZ.at(1) * nBinsXYZ.at(2) + binsXYZ.at(2));
+// },
+// vm["bf-map"].template as<std::string>(),
+// vm["bf-name"].template as<std::string>(), lengthUnit, BFieldUnit,
+// vm["bf-foctant"].template as<bool>());
+//
+// // create field mapping
+// InterpolatedBFieldMap3D::Config config3D(std::move(mapper3D));
+// config3D.scale = bscalor;
+// // create BField service
+// return std::make_shared<InterpolatedBFieldMap3D>(std::move(config3D));
+// }
+// } else if (bfieldmaptype == text) {
+// if (vm["bf-rz"].template as<bool>()) {
+// auto mapper2D = Jug::BField::txt::fieldMapperRZ(
+// [](std::array<size_t, 2> binsRZ, std::array<size_t, 2> nBinsRZ) {
+// return (binsRZ.at(1) * nBinsRZ.at(0) + binsRZ.at(0));
+// },
+// vm["bf-map"].template as<std::string>(), lengthUnit, BFieldUnit,
+// vm["bf-gridpoints"].template as<size_t>(),
+// vm["bf-foctant"].template as<bool>());
+//
+// // create field mapping
+// InterpolatedBFieldMap2D::Config config2D(std::move(mapper2D));
+// config2D.scale = bscalor;
+// // create BField service
+// return std::make_shared<InterpolatedBFieldMap2D>(std::move(config2D));
+//
+// } else {
+// auto mapper3D = Jug::BField::txt::fieldMapperXYZ(
+// [](std::array<size_t, 3> binsXYZ, std::array<size_t, 3> nBinsXYZ) {
+// return (binsXYZ.at(0) * (nBinsXYZ.at(1) * nBinsXYZ.at(2)) +
+// binsXYZ.at(1) * nBinsXYZ.at(2) + binsXYZ.at(2));
+// },
+// vm["bf-map"].template as<std::string>(), lengthUnit, BFieldUnit,
+// vm["bf-gridpoints"].template as<size_t>(),
+// vm["bf-foctant"].template as<bool>());
+//
+// // create field mapping
+// InterpolatedBFieldMap3D::Config config3D(std::move(mapper3D));
+// config3D.scale = bscalor;
+// // create BField service
+// return std::make_shared<InterpolatedBFieldMap3D>(std::move(config3D));
+// }
+// } else { // constant
+// // No bfield map is handed over
+// // get the constant bField values
+// auto bFieldValues = vm["bf-values"].template as<Reals<3>>();
+// if (vm["bf-context-scalable"].template as<bool>()) {
+// // Create the scalable magnetic field
+// return std::make_shared<Jug::BField::ScalableBField>(
+// bFieldValues.at(0) * Acts::UnitConstants::T,
+// bFieldValues.at(1) * Acts::UnitConstants::T,
+// bFieldValues.at(2) * Acts::UnitConstants::T);
+// } else {
+// // Create the constant magnetic field
+// return std::make_shared<Acts::ConstantBField>(
+// bFieldValues.at(0) * Acts::UnitConstants::T,
+// bFieldValues.at(1) * Acts::UnitConstants::T,
+// bFieldValues.at(2) * Acts::UnitConstants::T);
+// }
+// }
+//}
+//} // namespace Options
+} // namespace Jug
diff --git a/JugBase/src/Plugins/BFieldScalor.cpp b/JugBase/src/Plugins/BFieldScalor.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c4f0e28a26407cc6193c759b7a2a80c8bc693a0c
--- /dev/null
+++ b/JugBase/src/Plugins/BFieldScalor.cpp
@@ -0,0 +1,27 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2019 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "JugBase/BField/BFieldScalor.hpp"
+
+#include "JugBase/BField/ScalableBField.hpp"
+
+#include <cmath>
+
+Jug::BField::BFieldScalor::BFieldScalor(
+ const Jug::BField::BFieldScalor::Config& cfg,
+ std::unique_ptr<const Acts::Logger> logger)
+ : m_cfg(cfg), m_logger(std::move(logger)) {}
+
+Jug::ProcessCode Jug::BField::BFieldScalor::decorate(
+ AlgorithmContext& context) {
+ ScalableBFieldContext bFieldContext{
+ std::pow(m_cfg.scalor, context.eventNumber)};
+ context.magFieldContext = std::make_any<ScalableBFieldContext>(bFieldContext);
+
+ return ProcessCode::SUCCESS;
+}
diff --git a/JugBase/src/Plugins/BFieldUtils.cpp b/JugBase/src/Plugins/BFieldUtils.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0e778609a9209b75fcaaa638172283334928cb44
--- /dev/null
+++ b/JugBase/src/Plugins/BFieldUtils.cpp
@@ -0,0 +1,201 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2017 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "JugBase/BField/BFieldUtils.hpp"
+
+#include "Acts/MagneticField/BFieldMapUtils.hpp"
+#include "Acts/Utilities/detail/Axis.hpp"
+#include "Acts/Utilities/detail/Grid.hpp"
+
+#include <fstream>
+
+#include "TFile.h"
+#include "TROOT.h"
+#include "TTree.h"
+
+Acts::InterpolatedBFieldMapper<
+ Acts::detail::Grid<Acts::Vector2, Acts::detail::EquidistantAxis,
+ Acts::detail::EquidistantAxis>>
+Jug::BField::txt::fieldMapperRZ(
+ std::function<size_t(std::array<size_t, 2> binsRZ,
+ std::array<size_t, 2> nBinsRZ)>
+ localToGlobalBin,
+ std::string fieldMapFile, double lengthUnit, double BFieldUnit,
+ size_t nPoints, bool firstQuadrant) {
+ /// [1] Read in field map file
+ // Grid position points in r and z
+ std::vector<double> rPos;
+ std::vector<double> zPos;
+ // components of magnetic field on grid points
+ std::vector<Acts::Vector2> bField;
+ // reserve estimated size
+ rPos.reserve(nPoints);
+ zPos.reserve(nPoints);
+ bField.reserve(nPoints);
+ // [1] Read in file and fill values
+ std::ifstream map_file(fieldMapFile.c_str(), std::ios::in);
+ std::string line;
+ double r = 0., z = 0.;
+ double br = 0., bz = 0.;
+ while (std::getline(map_file, line)) {
+ if (line.empty() || line[0] == '%' || line[0] == '#' ||
+ line.find_first_not_of(' ') == std::string::npos)
+ continue;
+
+ std::istringstream tmp(line);
+ tmp >> r >> z >> br >> bz;
+ rPos.push_back(r);
+ zPos.push_back(z);
+ bField.push_back(Acts::Vector2(br, bz));
+ }
+ map_file.close();
+ /// [2] use helper function in core
+ return Acts::fieldMapperRZ(localToGlobalBin, rPos, zPos, bField, lengthUnit,
+ BFieldUnit, firstQuadrant);
+}
+
+Acts::InterpolatedBFieldMapper<Acts::detail::Grid<
+ Acts::Vector3, Acts::detail::EquidistantAxis, Acts::detail::EquidistantAxis,
+ Acts::detail::EquidistantAxis>>
+Jug::BField::txt::fieldMapperXYZ(
+ std::function<size_t(std::array<size_t, 3> binsXYZ,
+ std::array<size_t, 3> nBinsXYZ)>
+ localToGlobalBin,
+ std::string fieldMapFile, double lengthUnit, double BFieldUnit,
+ size_t nPoints, bool firstOctant) {
+ /// [1] Read in field map file
+ // Grid position points in x, y and z
+ std::vector<double> xPos;
+ std::vector<double> yPos;
+ std::vector<double> zPos;
+ // components of magnetic field on grid points
+ std::vector<Acts::Vector3> bField;
+ // reserve estimated size
+ xPos.reserve(nPoints);
+ yPos.reserve(nPoints);
+ zPos.reserve(nPoints);
+ bField.reserve(nPoints);
+ // [1] Read in file and fill values
+ std::ifstream map_file(fieldMapFile.c_str(), std::ios::in);
+ std::string line;
+ double x = 0., y = 0., z = 0.;
+ double bx = 0., by = 0., bz = 0.;
+ while (std::getline(map_file, line)) {
+ if (line.empty() || line[0] == '%' || line[0] == '#' ||
+ line.find_first_not_of(' ') == std::string::npos)
+ continue;
+
+ std::istringstream tmp(line);
+ tmp >> x >> y >> z >> bx >> by >> bz;
+ xPos.push_back(x);
+ yPos.push_back(y);
+ zPos.push_back(z);
+ bField.push_back(Acts::Vector3(bx, by, bz));
+ }
+ map_file.close();
+
+ return Acts::fieldMapperXYZ(localToGlobalBin, xPos, yPos, zPos, bField,
+ lengthUnit, BFieldUnit, firstOctant);
+}
+
+Acts::InterpolatedBFieldMapper<
+ Acts::detail::Grid<Acts::Vector2, Acts::detail::EquidistantAxis,
+ Acts::detail::EquidistantAxis>>
+Jug::BField::root::fieldMapperRZ(
+ std::function<size_t(std::array<size_t, 2> binsRZ,
+ std::array<size_t, 2> nBinsRZ)>
+ localToGlobalBin,
+ std::string fieldMapFile, std::string treeName, double lengthUnit,
+ double BFieldUnit, bool firstQuadrant) {
+ /// [1] Read in field map file
+ // Grid position points in r and z
+ std::vector<double> rPos;
+ std::vector<double> zPos;
+ // components of magnetic field on grid points
+ std::vector<Acts::Vector2> bField;
+ // [1] Read in file and fill values
+ TFile* inputFile = TFile::Open(fieldMapFile.c_str());
+ TTree* tree = (TTree*)inputFile->Get(treeName.c_str());
+ Int_t entries = tree->GetEntries();
+
+ double r, z;
+ double Br, Bz;
+
+ tree->SetBranchAddress("r", &r);
+ tree->SetBranchAddress("z", &z);
+
+ tree->SetBranchAddress("Br", &Br);
+ tree->SetBranchAddress("Bz", &Bz);
+
+ // reserve size
+ rPos.reserve(entries);
+ zPos.reserve(entries);
+ bField.reserve(entries);
+
+ for (int i = 0; i < entries; i++) {
+ tree->GetEvent(i);
+ rPos.push_back(r);
+ zPos.push_back(z);
+ bField.push_back(Acts::Vector2(Br, Bz));
+ }
+ inputFile->Close();
+ /// [2] use helper function in core
+ return Acts::fieldMapperRZ(localToGlobalBin, rPos, zPos, bField, lengthUnit,
+ BFieldUnit, firstQuadrant);
+}
+
+Acts::InterpolatedBFieldMapper<Acts::detail::Grid<
+ Acts::Vector3, Acts::detail::EquidistantAxis, Acts::detail::EquidistantAxis,
+ Acts::detail::EquidistantAxis>>
+Jug::BField::root::fieldMapperXYZ(
+ std::function<size_t(std::array<size_t, 3> binsXYZ,
+ std::array<size_t, 3> nBinsXYZ)>
+ localToGlobalBin,
+ std::string fieldMapFile, std::string treeName, double lengthUnit,
+ double BFieldUnit, bool firstOctant) {
+ /// [1] Read in field map file
+ // Grid position points in x, y and z
+ std::vector<double> xPos;
+ std::vector<double> yPos;
+ std::vector<double> zPos;
+ // components of magnetic field on grid points
+ std::vector<Acts::Vector3> bField;
+ // [1] Read in file and fill values
+ TFile* inputFile = TFile::Open(fieldMapFile.c_str());
+ TTree* tree = (TTree*)inputFile->Get(treeName.c_str());
+ Int_t entries = tree->GetEntries();
+
+ double x, y, z;
+ double Bx, By, Bz;
+
+ tree->SetBranchAddress("x", &x);
+ tree->SetBranchAddress("y", &y);
+ tree->SetBranchAddress("z", &z);
+
+ tree->SetBranchAddress("Bx", &Bx);
+ tree->SetBranchAddress("By", &By);
+ tree->SetBranchAddress("Bz", &Bz);
+
+ // reserve size
+ xPos.reserve(entries);
+ yPos.reserve(entries);
+ zPos.reserve(entries);
+ bField.reserve(entries);
+
+ for (int i = 0; i < entries; i++) {
+ tree->GetEvent(i);
+ xPos.push_back(x);
+ yPos.push_back(y);
+ zPos.push_back(z);
+ bField.push_back(Acts::Vector3(Bx, By, Bz));
+ }
+ inputFile->Close();
+
+ return Acts::fieldMapperXYZ(localToGlobalBin, xPos, yPos, zPos, bField,
+ lengthUnit, BFieldUnit, firstOctant);
+}
diff --git a/JugBase/src/Utilities/Helpers.cpp b/JugBase/src/Utilities/Helpers.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a26d42067ed6123e270362b1b2d7f2395a2993b2
--- /dev/null
+++ b/JugBase/src/Utilities/Helpers.cpp
@@ -0,0 +1,103 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2019 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "JugBase/Utilities/Helpers.hpp"
+
+namespace Jug {
+
+namespace PlotHelpers {
+TH1F* bookHisto(const char* histName, const char* histTitle,
+ const Binning& varBinning) {
+ TH1F* hist = new TH1F(histName, histTitle, varBinning.nBins, varBinning.min,
+ varBinning.max);
+ hist->GetXaxis()->SetTitle(varBinning.title.c_str());
+ hist->GetYaxis()->SetTitle("Entries");
+ hist->Sumw2();
+ return hist;
+}
+
+TH2F* bookHisto(const char* histName, const char* histTitle,
+ const Binning& varXBinning, const Binning& varYBinning) {
+ TH2F* hist = new TH2F(histName, histTitle, varXBinning.nBins, varXBinning.min,
+ varXBinning.max, varYBinning.nBins, varYBinning.min,
+ varYBinning.max);
+ hist->GetXaxis()->SetTitle(varXBinning.title.c_str());
+ hist->GetYaxis()->SetTitle(varYBinning.title.c_str());
+ hist->Sumw2();
+ return hist;
+}
+
+void fillHisto(TH1F* hist, float value, float weight) {
+ assert(hist != nullptr);
+ hist->Fill(value, weight);
+}
+
+void fillHisto(TH2F* hist, float xValue, float yValue, float weight) {
+ assert(hist != nullptr);
+ hist->Fill(xValue, yValue, weight);
+}
+
+void anaHisto(TH1D* inputHist, int j, TH1F* meanHist, TH1F* widthHist) {
+ // evaluate mean and width via the Gauss fit
+ assert(inputHist != nullptr);
+ if (inputHist->GetEntries() > 0) {
+ TFitResultPtr r = inputHist->Fit("gaus", "QS0");
+ if (r.Get() and ((r->Status() % 1000) == 0)) {
+ // fill the mean and width into 'j'th bin of the meanHist and widthHist,
+ // respectively
+ meanHist->SetBinContent(j, r->Parameter(1));
+ meanHist->SetBinError(j, r->ParError(1));
+ widthHist->SetBinContent(j, r->Parameter(2));
+ widthHist->SetBinError(j, r->ParError(2));
+ }
+ }
+}
+
+TEfficiency* bookEff(const char* effName, const char* effTitle,
+ const Binning& varBinning) {
+ TEfficiency* efficiency = new TEfficiency(effName, effTitle, varBinning.nBins,
+ varBinning.min, varBinning.max);
+ return efficiency;
+}
+
+TEfficiency* bookEff(const char* effName, const char* effTitle,
+ const Binning& varXBinning, const Binning& varYBinning) {
+ TEfficiency* efficiency = new TEfficiency(
+ effName, effTitle, varXBinning.nBins, varXBinning.min, varXBinning.max,
+ varYBinning.nBins, varYBinning.min, varYBinning.max);
+ return efficiency;
+}
+
+void fillEff(TEfficiency* efficiency, float value, bool status) {
+ assert(efficiency != nullptr);
+ efficiency->Fill(status, value);
+}
+
+void fillEff(TEfficiency* efficiency, float xValue, float yValue, bool status) {
+ assert(efficiency != nullptr);
+ efficiency->Fill(status, xValue, yValue);
+}
+
+TProfile* bookProf(const char* profName, const char* profTitle,
+ const Binning& varXBinning, const Binning& varYBinning) {
+ TProfile* prof =
+ new TProfile(profName, profTitle, varXBinning.nBins, varXBinning.min,
+ varXBinning.max, varYBinning.min, varYBinning.max);
+ prof->GetXaxis()->SetTitle(varXBinning.title.c_str());
+ prof->GetYaxis()->SetTitle(varYBinning.title.c_str());
+ return prof;
+}
+
+void fillProf(TProfile* profile, float xValue, float yValue, float weight) {
+ assert(profile != nullptr);
+ profile->Fill(xValue, yValue, weight);
+}
+
+} // namespace PlotHelpers
+
+} // namespace Jug
diff --git a/JugBase/src/Utilities/Options.cpp b/JugBase/src/Utilities/Options.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..024fc9d3bcebc0735a6978d95932a1d8a4da58f8
--- /dev/null
+++ b/JugBase/src/Utilities/Options.cpp
@@ -0,0 +1,162 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2020 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "JugBase/Utilities/Options.hpp"
+
+#include <algorithm>
+#include <istream>
+#include <ostream>
+#include <stdexcept>
+
+namespace {
+constexpr char s_separator = ':';
+}
+
+// interval
+
+std::istream& Jug::Options::operator>>(
+ std::istream& is, Jug::Options::Interval& interval) {
+ std::string buf;
+ is >> buf;
+
+ // default to an unbounded interval
+ interval.lower.reset();
+ interval.upper.reset();
+
+ // find the limit separator
+ auto pos = buf.find_first_of(s_separator);
+ // no separator -> invalid input -> unbounded interval
+ if (pos == std::string::npos) {
+ return is;
+ }
+
+ // if it exists, parse limit before separator
+ if (0 < pos) {
+ auto lowerStr = buf.substr(0, pos);
+ interval.lower = std::stod(lowerStr);
+ }
+ // if it exists, parse limit after separator
+ if ((pos + 1) < buf.size()) {
+ auto upperStr = buf.substr(pos + 1);
+ interval.upper = std::stod(upperStr);
+ }
+
+ return is;
+}
+
+std::ostream& Jug::Options::operator<<(
+ std::ostream& os, const Jug::Options::Interval& interval) {
+ if (not interval.lower.has_value() and not interval.upper.has_value()) {
+ os << "unbounded";
+ } else {
+ if (interval.lower.has_value()) {
+ os << interval.lower.value();
+ }
+ os << s_separator;
+ if (interval.upper.has_value()) {
+ os << interval.upper.value();
+ }
+ }
+ return os;
+}
+
+// helper functions to parse and print multiple values
+
+namespace {
+
+template <typename value_t, typename converter_t>
+void parseVariable(std::istream& is, std::vector<value_t>& values,
+ converter_t&& convert) {
+ values.clear();
+
+ std::string buf;
+ is >> buf;
+ std::string bufValue;
+ std::string::size_type pos = 0;
+ std::string::size_type end = std::string::npos;
+ do {
+ end = buf.find_first_of(s_separator, pos);
+ if (end == std::string::npos) {
+ // last element; take the rest of the buffer
+ bufValue = buf.substr(pos);
+ } else {
+ bufValue = buf.substr(pos, end - pos);
+ pos = end + 1u;
+ }
+ values.push_back(convert(bufValue));
+ } while (end != std::string::npos);
+}
+
+template <typename value_t, typename converter_t>
+void parseFixed(std::istream& is, size_t size, value_t* values,
+ converter_t&& convert) {
+ // reserve space for the expected number of values
+ std::vector<value_t> tmp(size, 0);
+ parseVariable(is, tmp, std::forward<converter_t>(convert));
+ if (tmp.size() < size) {
+ throw std::invalid_argument(
+ "Not enough values for fixed-size user option, expected " +
+ std::to_string(size) + " received " + std::to_string(tmp.size()));
+ }
+ if (size < tmp.size()) {
+ throw std::invalid_argument(
+ "Too many values for fixed-size user option, expected " +
+ std::to_string(size) + " received " + std::to_string(tmp.size()));
+ }
+ std::copy(tmp.begin(), tmp.end(), values);
+}
+
+template <typename value_t>
+void print(std::ostream& os, size_t size, const value_t* values) {
+ for (size_t i = 0; i < size; ++i) {
+ if (0u < i) {
+ os << s_separator;
+ }
+ os << values[i];
+ }
+}
+
+} // namespace
+
+// fixed and variable number of generic values
+
+void Jug::Options::detail::parseDoublesFixed(std::istream& is,
+ size_t size,
+ double* values) {
+ parseFixed(is, size, values,
+ [](const std::string& s) { return std::stod(s); });
+}
+
+void Jug::Options::detail::parseDoublesVariable(
+ std::istream& is, std::vector<double>& values) {
+ parseVariable(is, values, [](const std::string& s) { return std::stod(s); });
+}
+
+void Jug::Options::detail::printDoubles(std::ostream& os, size_t size,
+ const double* values) {
+ print(os, size, values);
+}
+
+// fixed and variable number of integers
+
+void Jug::Options::detail::parseIntegersFixed(std::istream& is,
+ size_t size,
+ int* values) {
+ parseFixed(is, size, values,
+ [](const std::string& s) { return std::stoi(s); });
+}
+
+void Jug::Options::detail::parseIntegersVariable(
+ std::istream& is, std::vector<int>& values) {
+ parseVariable(is, values, [](const std::string& s) { return std::stoi(s); });
+}
+
+void Jug::Options::detail::printIntegers(std::ostream& os, size_t size,
+ const int* values) {
+ print(os, size, values);
+}
diff --git a/JugBase/src/Utilities/Paths.cpp b/JugBase/src/Utilities/Paths.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..cc4d3bff1d71a5788ca5dc32c6750300bfa2dc14
--- /dev/null
+++ b/JugBase/src/Utilities/Paths.cpp
@@ -0,0 +1,112 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2017 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "JugBase/Utilities/Paths.hpp"
+
+#include "Acts/Utilities/Logger.hpp"
+
+#include <charconv>
+#include <cstdio>
+#include <regex>
+#include <sstream>
+#include <stdexcept>
+
+#include <boost/filesystem.hpp>
+
+std::string Jug::ensureWritableDirectory(const std::string& dir) {
+ using boost::filesystem::current_path;
+ using boost::filesystem::path;
+
+ auto dir_path = dir.empty() ? current_path() : path(dir);
+ if (exists(dir_path) and not is_directory(dir_path)) {
+ throw std::runtime_error("'" + dir +
+ "' already exists but is not a directory");
+ }
+ create_directories(dir_path);
+ return canonical(dir_path).native();
+}
+
+std::string Jug::joinPaths(const std::string& dir,
+ const std::string& name) {
+ if (dir.empty()) {
+ return name;
+ } else {
+ return dir + '/' + name;
+ }
+}
+
+std::string Jug::perEventFilepath(const std::string& dir,
+ const std::string& name,
+ size_t event) {
+ char prefix[64];
+
+ snprintf(prefix, sizeof(prefix), "event%09zu-", event);
+
+ if (dir.empty()) {
+ return prefix + name;
+ } else {
+ return dir + '/' + prefix + name;
+ }
+}
+
+std::pair<size_t, size_t> Jug::determineEventFilesRange(
+ const std::string& dir, const std::string& name) {
+ using boost::filesystem::current_path;
+ using boost::filesystem::directory_iterator;
+ using boost::filesystem::path;
+
+ ACTS_LOCAL_LOGGER(
+ Acts::getDefaultLogger("EventFilesRange", Acts::Logging::VERBOSE));
+
+ // ensure directory path is valid
+ auto dir_path = dir.empty() ? current_path() : path(dir);
+ if (not exists(dir_path)) {
+ throw std::runtime_error("'" + dir_path.native() + "' does not exists");
+ }
+ if (not is_directory(dir_path)) {
+ throw std::runtime_error("'" + dir_path.native() + "' is not a directory");
+ }
+
+ // invalid default range that allows simple restriction later on
+ size_t eventMin = SIZE_MAX;
+ size_t eventMax = 0;
+
+ // filter matching event files from the directory listing
+ std::string filename;
+ std::regex re("^event([0-9]+)-" + name + "$");
+ std::cmatch match;
+
+ for (const auto& f : directory_iterator(dir_path)) {
+ if ((not exists(f.status())) or (not is_regular_file(f.status()))) {
+ continue;
+ }
+ // keep a copy so the match can refer to the underlying const char*
+ filename = f.path().filename().native();
+ if (std::regex_match(filename.c_str(), match, re)) {
+ ACTS_VERBOSE("Matching file " << filename);
+
+ // first sub_match is the whole string, second should be the event number
+ size_t event = 0;
+ auto ret = std::from_chars(match[1].first, match[1].second, event);
+ if (ret.ptr == match[1].first) {
+ throw std::runtime_error(
+ "Could not extract event number from filename");
+ }
+ // enlarge available event range
+ eventMin = std::min(eventMin, event);
+ eventMax = std::max(eventMax, event);
+ }
+ }
+ ACTS_VERBOSE("Detected event range [" << eventMin << "," << eventMax << "]");
+
+ // should only occur if no files matched and the initial values persisted.
+ if (eventMax < eventMin) {
+ return std::make_pair(0u, 0u);
+ }
+ return std::make_pair(eventMin, eventMax + 1);
+}
diff --git a/JugBase/src/components/GeoSvc.cpp b/JugBase/src/components/GeoSvc.cpp
index d2f5f85efc8bd60ee612b65d929a284c47203624..e85708feb013974ee90d42244547611766ced6c1 100644
--- a/JugBase/src/components/GeoSvc.cpp
+++ b/JugBase/src/components/GeoSvc.cpp
@@ -17,6 +17,44 @@
#include "Acts/Geometry/TrackingGeometry.hpp"
#include "Acts/Plugins/DD4hep/ConvertDD4hepDetector.hpp"
+#include "Acts/Surfaces/PlaneSurface.hpp"
+
+
+void draw_surfaces(std::shared_ptr<const Acts::TrackingGeometry> trk_geo, const std::string& fname)
+{
+ using namespace Acts;
+ Acts::GeometryContext tgContext = Acts::GeometryContext();
+ std::vector<const Surface*> surfaces;
+
+ trk_geo->visitSurfaces([&](const Acts::Surface* surface) {
+ // for now we just require a valid surface
+ if (not surface) {
+ return;
+ }
+ surfaces.push_back(surface);
+ });
+ std::ofstream os;
+ os.open(fname);
+ os << std::fixed << std::setprecision(4);
+ size_t nVtx = 0;
+ for (const auto& srfx : surfaces) {
+ const PlaneSurface* srf = dynamic_cast<const PlaneSurface*>(srfx);
+ const PlanarBounds* bounds = dynamic_cast<const PlanarBounds*>(&srf->bounds());
+ for (const auto& vtxloc : bounds->vertices()) {
+ Vector3 vtx = srf->transform(tgContext) * Vector3(vtxloc.x(), vtxloc.y(), 0);
+ os << "v " << vtx.x() << " " << vtx.y() << " " << vtx.z() << "\n";
+ }
+ // connect them
+ os << "f";
+ for (size_t i = 1; i <= bounds->vertices().size(); ++i) {
+ os << " " << nVtx + i;
+ }
+ os << "\n";
+ nVtx += bounds->vertices().size();
+ }
+ os.close();
+}
+
using namespace Gaudi;
DECLARE_COMPONENT(GeoSvc)
@@ -85,8 +123,9 @@ StatusCode GeoSvc::initialize() {
default:
geoMsgLevel = Acts::Logging::VERBOSE;
}
- m_trackingGeo = std::move(Acts::convertDD4hepDetector(m_dd4hepgeo->world(), geoMsgLevel, Acts::equidistant,
+ m_trackingGeo = std::move(Acts::convertDD4hepDetector(m_dd4hepgeo->world(), Acts::Logging::VERBOSE, Acts::equidistant,
Acts::equidistant, Acts::equidistant));
+ draw_surfaces(m_trackingGeo, "Derp.obj");
return StatusCode::SUCCESS;
}
diff --git a/JugBase/src/components/GeoSvc.h b/JugBase/src/components/GeoSvc.h
index 87fa21ce408bd11199fcfc935a08f44c1c10047f..2efeeec3c284fb31087fab16d1b2f0c42aa808c7 100644
--- a/JugBase/src/components/GeoSvc.h
+++ b/JugBase/src/components/GeoSvc.h
@@ -11,7 +11,7 @@
// Interface
#include "JugBase/IGeoSvc.h"
-#include "Acts/Utilities/Units.hpp"
+#include "Acts/Definitions/Units.hpp"
#include "DD4hep/DD4hepUnits.h"
//using namespace Acts::UnitLiterals;
@@ -27,17 +27,28 @@
#include "DDRec/SurfaceManager.h"
#include "DDRec/Surface.h"
+// Create a test context
+//#define CHECK_ROTATION_ANGLE(t, a, tolerance) \
+// { \
+// Vector3 v = (*t) * Vector3(1, 0, 0); \
+// CHECK_CLOSE_ABS(VectorHelpers::phi(v), (a), tolerance); \
+// }
+
+//using SrfVec = std::vector<std::shared_ptr<const Surface>>;
+void draw_surfaces(std::shared_ptr<const Acts::TrackingGeometry> trk_geo, const std::string& fname);
+
+
class GeoSvc : public extends<Service, IGeoSvc> {
public:
GeoSvc(const std::string& name, ISvcLocator* svc);
virtual ~GeoSvc();
-
+
virtual StatusCode initialize() final;
-
+
virtual StatusCode finalize() final;
-
+
/** Build the dd4hep geometry.
* This function generates the DD4hep geometry.
*/
@@ -64,8 +75,6 @@ public:
virtual double centralMagneticField() const {
return m_dd4hepgeo->field().magneticField({0,0,0}).z()*(Acts::UnitConstants::T/dd4hep::tesla);
}
-
-;
private:
@@ -76,9 +85,7 @@ private:
/// Pointer to the interface to the DD4hep geometry
dd4hep::Detector* m_dd4hepgeo;
-
std::shared_ptr<const dd4hep::rec::CellIDPositionConverter> m_cellid_converter = nullptr;//(*(m_geoSvc->detector()));
-
/// XML-files with the detector description
Gaudi::Property<std::vector<std::string>> m_xmlFileNames{this, "detectors", {}, "Detector descriptions XML-files"};
diff --git a/JugReco/CMakeLists.txt b/JugReco/CMakeLists.txt
index 5df1eed1d4cb4824d893837b837a8ee69c41f30a..9935c87515dc4313721cf2db06c9cc4058530c86 100644
--- a/JugReco/CMakeLists.txt
+++ b/JugReco/CMakeLists.txt
@@ -31,8 +31,10 @@ gaudi_add_module(JugRecoPlugins
src/components/TrackerHitReconstruction.cpp
src/components/TrackerSourceLinker.cpp
src/components/TrackerSourcesLinker.cpp
- src/components/TrackingHitsSourceLinker.cpp
+ #src/components/TrackingHitsSourceLinker.cpp
src/components/TrackFindingAlgorithm.cpp
+ #src/components/TrackFittingAlgorithm.cpp
+ #src/components/TrackFittingFunction.cpp
src/components/TestACTSLogger.cpp
src/components/TrackParamTruthInit.cpp
src/components/TrackParamClusterInit.cpp
diff --git a/JugReco/JugReco/BField.h b/JugReco/JugReco/BField.h
index b013d99194f9ed72ca052f3e273048c3938e4330..0d3bb15a1957f0062ff0289a21486b2c24ea8c5e 100644
--- a/JugReco/JugReco/BField.h
+++ b/JugReco/JugReco/BField.h
@@ -1,7 +1,7 @@
#ifndef Jug_BFIELD_HH
#define Jug_BFIELD_HH 1
-#include "Acts/Utilities//Definitions.hpp"
+#include "Acts//Definitions/Units.hpp"
#include "Acts/Utilities/detail/AxisFwd.hpp"
#include "Acts/Utilities/detail/GridFwd.hpp"
#include <memory>
@@ -9,7 +9,7 @@
#include<variant>
#include "Acts/MagneticField/MagneticFieldContext.hpp"
-#include "Acts/Utilities/Definitions.hpp"
+#include "Acts/Definitions/Common.hpp"
// Forward declarations
namespace Acts {
@@ -22,141 +22,138 @@ namespace Acts {
} // namespace Acts
namespace Jug {
- //namespace BField {
+ // namespace BField {
// class ScalableBField;
//}
-namespace BField {
-
-/// The Context to be handed around
-struct ScalableBFieldContext {
- double scalor = 1.;
-};
-
-/// @ingroup MagneticField
-///
-/// @brief returns a given constant field value at every point
-///
-/// This class is based on the constant magnetic field class
-/// but allows a event based context
-class ScalableBField final {
- public:
- struct Cache {
- double scalor = 1.;
-
- /// @brief constructor with context
- Cache(const Acts::MagneticFieldContext& mcfg) {
- scalor = std::any_cast<const ScalableBFieldContext>(mcfg).scalor;
- }
- };
-
- /// @brief construct constant magnetic field from field vector
- ///
- /// @param [in] B magnetic field vector in global coordinate system
- explicit ScalableBField(Acts::Vector3D B) : m_BField(std::move(B)) {}
-
- /// @brief construct constant magnetic field from components
- ///
- /// @param [in] Bx magnetic field component in global x-direction
- /// @param [in] By magnetic field component in global y-direction
- /// @param [in] Bz magnetic field component in global z-direction
- ScalableBField(double Bx = 0., double By = 0., double Bz = 0.)
- : m_BField(Bx, By, Bz) {}
-
- /// @brief retrieve magnetic field value
- ///
- /// @param [in] position global position
- /// @return magnetic field vector
- ///
- /// @note The @p position is ignored and only kept as argument to provide
- /// a consistent interface with other magnetic field services.
- Acts::Vector3D getField(const Acts::Vector3D& /*position*/) const {
- return m_BField;
- }
-
- /// @brief retrieve magnetic field value
- ///
- /// @param [in] position global position
- /// @param [in] cache Cache object (is ignored)
- /// @return magnetic field vector
- ///
- /// @note The @p position is ignored and only kept as argument to provide
- /// a consistent interface with other magnetic field services.
- Acts::Vector3D getField(const Acts::Vector3D& /*position*/,
- Cache& cache) const {
- return m_BField * cache.scalor;
- }
-
- /// @brief retrieve magnetic field value & its gradient
- ///
- /// @param [in] position global position
- /// @param [out] derivative gradient of magnetic field vector as (3x3)
- /// matrix
- /// @return magnetic field vector
- ///
- /// @note The @p position is ignored and only kept as argument to provide
- /// a consistent interface with other magnetic field services.
- /// @note currently the derivative is not calculated
- /// @todo return derivative
- Acts::Vector3D getFieldGradient(
- const Acts::Vector3D& /*position*/,
- Acts::ActsMatrixD<3, 3>& /*derivative*/) const {
- return m_BField;
- }
-
- /// @brief retrieve magnetic field value & its gradient
- ///
- /// @param [in] position global position
- /// @param [out] derivative gradient of magnetic field vector as (3x3)
- /// matrix
- /// @param [in] cache Cache object (is ignored)
- /// @return magnetic field vector
- ///
- /// @note The @p position is ignored and only kept as argument to provide
- /// a consistent interface with other magnetic field services.
- /// @note currently the derivative is not calculated
- /// @todo return derivative
- Acts::Vector3D getFieldGradient(const Acts::Vector3D& /*position*/,
- Acts::ActsMatrixD<3, 3>& /*derivative*/,
- Cache& cache) const {
- return m_BField * cache.scalor;
- }
-
- /// @brief check whether given 3D position is inside look-up domain
- ///
- /// @param [in] position global 3D position
- /// @return @c true if position is inside the defined look-up grid,
- /// otherwise @c false
- /// @note The method will always return true for the constant B-Field
- bool isInside(const Acts::Vector3D& /*position*/) const { return true; }
-
- /// @brief update magnetic field vector from components
- ///
- /// @param [in] Bx magnetic field component in global x-direction
- /// @param [in] By magnetic field component in global y-direction
- /// @param [in] Bz magnetic field component in global z-direction
- void setField(double Bx, double By, double Bz) { m_BField << Bx, By, Bz; }
-
- /// @brief update magnetic field vector
- ///
- /// @param [in] B magnetic field vector in global coordinate system
- void setField(const Acts::Vector3D& B) { m_BField = B; }
-
- private:
- /// magnetic field vector
- Acts::Vector3D m_BField;
-};
+ namespace BField {
+
+ /// The Context to be handed around
+ struct ScalableBFieldContext {
+ double scalor = 1.;
+ };
+
+ /// @ingroup MagneticField
+ ///
+ /// @brief returns a given constant field value at every point
+ ///
+ /// This class is based on the constant magnetic field class
+ /// but allows a event based context
+ class ScalableBField final {
+ public:
+ struct Cache {
+ double scalor = 1.;
+
+ /// @brief constructor with context
+ Cache(const Acts::MagneticFieldContext& mcfg)
+ {
+ scalor = std::any_cast<const ScalableBFieldContext>(mcfg).scalor;
+ }
+ };
+
+ /// @brief construct constant magnetic field from field vector
+ ///
+ /// @param [in] B magnetic field vector in global coordinate system
+ explicit ScalableBField(Acts::Vector3 B) : m_BField(std::move(B)) {}
+
+ /// @brief construct constant magnetic field from components
+ ///
+ /// @param [in] Bx magnetic field component in global x-direction
+ /// @param [in] By magnetic field component in global y-direction
+ /// @param [in] Bz magnetic field component in global z-direction
+ ScalableBField(double Bx = 0., double By = 0., double Bz = 0.) : m_BField(Bx, By, Bz) {}
+
+ /// @brief retrieve magnetic field value
+ ///
+ /// @param [in] position global position
+ /// @return magnetic field vector
+ ///
+ /// @note The @p position is ignored and only kept as argument to provide
+ /// a consistent interface with other magnetic field services.
+ Acts::Vector3 getField(const Acts::Vector3& /*position*/) const { return m_BField; }
+
+ /// @brief retrieve magnetic field value
+ ///
+ /// @param [in] position global position
+ /// @param [in] cache Cache object (is ignored)
+ /// @return magnetic field vector
+ ///
+ /// @note The @p position is ignored and only kept as argument to provide
+ /// a consistent interface with other magnetic field services.
+ Acts::Vector3 getField(const Acts::Vector3& /*position*/, Cache& cache) const
+ {
+ return m_BField * cache.scalor;
+ }
+
+ /// @brief retrieve magnetic field value & its gradient
+ ///
+ /// @param [in] position global position
+ /// @param [out] derivative gradient of magnetic field vector as (3x3)
+ /// matrix
+ /// @return magnetic field vector
+ ///
+ /// @note The @p position is ignored and only kept as argument to provide
+ /// a consistent interface with other magnetic field services.
+ /// @note currently the derivative is not calculated
+ /// @todo return derivative
+ Acts::Vector3 getFieldGradient(const Acts::Vector3& /*position*/, Acts::ActsMatrix<3, 3>& /*derivative*/) const
+ {
+ return m_BField;
+ }
+
+ /// @brief retrieve magnetic field value & its gradient
+ ///
+ /// @param [in] position global position
+ /// @param [out] derivative gradient of magnetic field vector as (3x3)
+ /// matrix
+ /// @param [in] cache Cache object (is ignored)
+ /// @return magnetic field vector
+ ///
+ /// @note The @p position is ignored and only kept as argument to provide
+ /// a consistent interface with other magnetic field services.
+ /// @note currently the derivative is not calculated
+ /// @todo return derivative
+ Acts::Vector3 getFieldGradient(const Acts::Vector3& /*position*/, Acts::ActsMatrix<3, 3>& /*derivative*/,
+ Cache& cache) const
+ {
+ return m_BField * cache.scalor;
+ }
+
+ /// @brief check whether given 3D position is inside look-up domain
+ ///
+ /// @param [in] position global 3D position
+ /// @return @c true if position is inside the defined look-up grid,
+ /// otherwise @c false
+ /// @note The method will always return true for the constant B-Field
+ bool isInside(const Acts::Vector3& /*position*/) const { return true; }
+
+ /// @brief update magnetic field vector from components
+ ///
+ /// @param [in] Bx magnetic field component in global x-direction
+ /// @param [in] By magnetic field component in global y-direction
+ /// @param [in] Bz magnetic field component in global z-direction
+ void setField(double Bx, double By, double Bz) { m_BField << Bx, By, Bz; }
+
+ /// @brief update magnetic field vector
+ ///
+ /// @param [in] B magnetic field vector in global coordinate system
+ void setField(const Acts::Vector3& B) { m_BField = B; }
+
+ private:
+ /// magnetic field vector
+ Acts::Vector3 m_BField;
+ };
} // namespace BField
} // namespace Jug
using InterpolatedMapper2D = Acts::InterpolatedBFieldMapper<
- Acts::detail::Grid<Acts::Vector2D, Acts::detail::EquidistantAxis, Acts::detail::EquidistantAxis>>;
+ Acts::detail::Grid<Acts::Vector2, Acts::detail::EquidistantAxis, Acts::detail::EquidistantAxis>>;
using InterpolatedMapper3D =
- Acts::InterpolatedBFieldMapper<Acts::detail::Grid<Acts::Vector3D, Acts::detail::EquidistantAxis,
+ Acts::InterpolatedBFieldMapper<Acts::detail::Grid<Acts::Vector3, Acts::detail::EquidistantAxis,
Acts::detail::EquidistantAxis, Acts::detail::EquidistantAxis>>;
-using InterpolatedBFieldMap2D = Acts::InterpolatedBFieldMap<InterpolatedMapper2D>;
-using InterpolatedBFieldMap3D = Acts::InterpolatedBFieldMap<InterpolatedMapper3D>;
+using InterpolatedBFieldMap2D = Acts::InterpolatedBFieldMap<InterpolatedMapper2D>;
+using InterpolatedBFieldMap3D = Acts::InterpolatedBFieldMap<InterpolatedMapper3D>;
using BFieldVariant = std::variant<std::shared_ptr<InterpolatedBFieldMap2D>, std::shared_ptr<InterpolatedBFieldMap3D>,
std::shared_ptr<Acts::ConstantBField>, std::shared_ptr<Jug::BField::ScalableBField>>;
diff --git a/JugReco/JugReco/Index.hpp b/JugReco/JugReco/Index.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..cb4eed81c098f7b79b12e9b6fb6f4a358e874204
--- /dev/null
+++ b/JugReco/JugReco/Index.hpp
@@ -0,0 +1,55 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2019-2020 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include <cstdint>
+
+#include <boost/container/flat_map.hpp>
+
+namespace Jug {
+
+/// Index type to reference elements in a container.
+///
+/// We do not expect to have more than 2^32 elements in any given container so a
+/// fixed sized integer type is sufficient.
+using Index = uint32_t;
+
+/// Store elements that are identified by an index, e.g. in another container.
+///
+/// Each index can have zero or more associated elements. A typical case could
+/// be to store all generating particles for a hit where the hit is identified
+/// by its index in the hit container.
+template <typename value_t>
+using IndexMultimap = boost::container::flat_multimap<Index, value_t>;
+
+/// Invert the multimap, i.e. from a -> {b...} to b -> {a...}.
+///
+/// @note This assumes that the value in the initial multimap is itself a
+/// sortable index-like object, as would be the case when mapping e.g.
+/// hit ids to particle ids/ barcodes.
+template <typename value_t>
+inline boost::container::flat_multimap<value_t, Index> invertIndexMultimap(
+ const IndexMultimap<value_t>& multimap) {
+ using InverseMultimap = boost::container::flat_multimap<value_t, Index>;
+
+ // switch key-value without enforcing the new ordering (linear copy)
+ typename InverseMultimap::sequence_type unordered;
+ unordered.reserve(multimap.size());
+ for (auto&& [index, value] : multimap) {
+ // value is now the key and the index is now the value
+ unordered.emplace_back(value, index);
+ }
+
+ // adopting the unordered sequence will reestablish the correct order
+ InverseMultimap inverse;
+ inverse.adopt_sequence(std::move(unordered));
+ return inverse;
+}
+
+} // namespace ActsExamples
diff --git a/JugReco/JugReco/IndexSourceLink.hpp b/JugReco/JugReco/IndexSourceLink.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..cfdb48b91ba5eb62c2da545c57938d1edc59b6f2
--- /dev/null
+++ b/JugReco/JugReco/IndexSourceLink.hpp
@@ -0,0 +1,62 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2020 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "JugReco/GeometryContainers.hpp"
+#include "JugReco/Index.hpp"
+
+#include <cassert>
+
+namespace Jug {
+
+ /// A source link that stores just an index.
+ ///
+ /// This is intentionally kept as barebones as possible. The source link
+ /// is just a reference and will be copied, moved around, etc. often.
+ /// Keeping it small and separate from the actual, potentially large,
+ /// measurement data should result in better overall performance.
+ /// Using an index instead of e.g. a pointer, means source link and
+ /// measurement are decoupled and the measurement represenation can be
+ /// easily changed without having to also change the source link.
+ class IndexSourceLink final {
+ public:
+ /// Construct from geometry identifier and index.
+ constexpr IndexSourceLink(Acts::GeometryIdentifier gid, Index idx) : m_geometryId(gid), m_index(idx) {}
+
+ // Construct an invalid source link. Must be default constructible to
+ /// satisfy SourceLinkConcept.
+ IndexSourceLink() = default;
+ IndexSourceLink(const IndexSourceLink&) = default;
+ IndexSourceLink(IndexSourceLink&&) = default;
+ IndexSourceLink& operator=(const IndexSourceLink&) = default;
+ IndexSourceLink& operator=(IndexSourceLink&&) = default;
+
+ /// Access the geometry identifier.
+ constexpr Acts::GeometryIdentifier geometryId() const { return m_geometryId; }
+ /// Access the index.
+ constexpr Index index() const { return m_index; }
+
+ public:
+ Acts::GeometryIdentifier m_geometryId;
+ Index m_index;
+
+ friend constexpr bool operator==(const IndexSourceLink& lhs, const IndexSourceLink& rhs)
+ {
+ return (lhs.m_geometryId == rhs.m_geometryId) and (lhs.m_index == rhs.m_index);
+ }
+ friend constexpr bool operator!=(const IndexSourceLink& lhs, const IndexSourceLink& rhs) { return not(lhs == rhs); }
+ };
+
+ /// Container of index source links.
+ ///
+ /// Since the source links provide a `.geometryId()` accessor, they can be
+ /// stored in an ordered geometry container.
+ using IndexSourceLinkContainer = GeometryIdMultiset<IndexSourceLink>;
+
+} // namespace Jug
diff --git a/JugReco/JugReco/Measurement.hpp b/JugReco/JugReco/Measurement.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ccd7725cbac029f5a1105751e3da6fefb1dbfce3
--- /dev/null
+++ b/JugReco/JugReco/Measurement.hpp
@@ -0,0 +1,52 @@
+// Copyright (C) 2020 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/EventData/Measurement.hpp"
+#include "JugReco/IndexSourceLink.hpp"
+
+#include <cassert>
+#include <vector>
+
+namespace Jug {
+
+ /// Variable measurement type that can contain all possible combinations.
+ using Measurement = ::Acts::BoundVariantMeasurement<IndexSourceLink>;
+ /// Container of measurements.
+ ///
+ /// In contrast to the source links, the measurements themself must not be
+ /// orderable. The source links stored in the measurements are treated
+ /// as opaque here and no ordering is enforced on the stored measurements.
+ using MeasurementContainer = std::vector<Measurement>;
+
+ /// Calibrator to convert an index source link to a measurement.
+ class MeasurementCalibrator {
+ public:
+ /// Construct an invalid calibrator. Required to allow copying.
+ MeasurementCalibrator() = default;
+ /// Construct using a user-provided container to chose measurements from.
+ MeasurementCalibrator(const MeasurementContainer& measurements) : m_measurements(&measurements) {}
+
+ /// Find the measurement corresponding to the source link.
+ ///
+ /// @tparam parameters_t Track parameters type
+ /// @param sourceLink Input source link
+ /// @param parameters Input track parameters (unused)
+ template <typename parameters_t>
+ const Measurement& operator()(const IndexSourceLink& sourceLink, const parameters_t& /* parameters */) const
+ {
+ assert(m_measurements and "Undefined measurement container in DigitizedCalibrator");
+ assert((sourceLink.index() < m_measurements->size()) and "Source link index is outside the container bounds");
+ return (*m_measurements)[sourceLink.m_index];
+ }
+
+ private:
+ // use pointer so the calibrator is copyable and default constructible.
+ const MeasurementContainer* m_measurements = nullptr;
+ };
+
+} // namespace Jug
diff --git a/JugReco/JugReco/SimMultiTrajectory.hpp b/JugReco/JugReco/SimMultiTrajectory.hpp
index 969a9e8f348767b080e9a88fb7e1583b58f65ec6..91b72ddcd4bc6cd425a1d2771fee33478cd2aae0 100644
--- a/JugReco/JugReco/SimMultiTrajectory.hpp
+++ b/JugReco/JugReco/SimMultiTrajectory.hpp
@@ -8,11 +8,11 @@
#pragma once
-#include "JugReco/SourceLinks.h"
//#include "ACTFW/Validation/ProtoTrackClassification.hpp"
#include "Acts/EventData/MultiTrajectory.hpp"
#include "Acts/EventData/TrackParameters.hpp"
+#include "JugReco/IndexSourceLink.hpp"
#include <unordered_map>
#include <utility>
@@ -48,7 +48,7 @@ struct SimMultiTrajectory {
/// @param tTips The entry indices for trajectories in multiTrajectory
/// @param parameters The fitted track parameters indexed by trajectory entry
/// index
- SimMultiTrajectory(const Acts::MultiTrajectory<SourceLink>& multiTraj,
+ SimMultiTrajectory(const Acts::MultiTrajectory<IndexSourceLink>& multiTraj,
const std::vector<size_t>& tTips,
const IndexedParams& parameters)
: m_multiTrajectory(multiTraj),
@@ -118,7 +118,7 @@ struct SimMultiTrajectory {
/// @return The multiTrajectory with trajectory entry indices
///
/// @note It could return an empty multiTrajectory
- std::pair<std::vector<size_t>, Acts::MultiTrajectory<SourceLink>>
+ std::pair<std::vector<size_t>, Acts::MultiTrajectory<IndexSourceLink>>
trajectory() const {
return std::make_pair(m_trackTips, m_multiTrajectory);
}
@@ -149,7 +149,7 @@ struct SimMultiTrajectory {
private:
// The multiTrajectory
- Acts::MultiTrajectory<SourceLink> m_multiTrajectory;
+ Acts::MultiTrajectory<IndexSourceLink> m_multiTrajectory;
// The entry indices of trajectories stored in multiTrajectory
std::vector<size_t> m_trackTips = {};
diff --git a/JugReco/JugReco/SourceLinks.h b/JugReco/JugReco/SourceLinks.h
index d53ab710a89a325dda099c43c48c47317c9be082..86442375c79f8b5ecc965ba3703720613957bcae 100644
--- a/JugReco/JugReco/SourceLinks.h
+++ b/JugReco/JugReco/SourceLinks.h
@@ -27,8 +27,9 @@ class SourceLink {
private:
Acts::BoundVector m_values;
Acts::BoundMatrix m_cov;
- size_t m_dim = 0u;
+ size_t m_dim = 2;
// store geo id copy to avoid indirection via truth hit
+ int32_t m_index;
Acts::GeometryIdentifier m_geometryId;
// need to store pointers to make the object copyable
const Acts::Surface* m_surface;
@@ -37,10 +38,11 @@ class SourceLink {
public:
SourceLink(const Acts::Surface& surface, //const ActsFatras::Hit& truthHit,
- size_t dim, Acts::BoundVector values, Acts::BoundMatrix cov)
+ size_t dim, int32_t index, Acts::BoundVector values, Acts::BoundMatrix cov)
: m_values(values),
m_cov(cov),
m_dim(dim),
+ m_index(index),
m_geometryId(0),//truthHit.geometryId()),
m_surface(&surface){}
//m_truthHit(&truthHit) {}
@@ -76,9 +78,14 @@ class SourceLink {
friend constexpr bool operator==(const SourceLink& lhs,
const SourceLink& rhs) {
- return lhs.geometryId() == rhs.geometryId();
+
+ return (lhs.geometryId() == rhs.geometryId()) && (lhs.m_index == rhs.m_index);
//lhs.m_truthHit == rhs.m_truthHit;
}
+ friend constexpr bool operator!=(const SourceLink& lhs,
+ const SourceLink& rhs) {
+ return not(lhs == rhs);
+ }
};
/// Store source links ordered by geometry identifier.
diff --git a/JugReco/JugReco/Track.hpp b/JugReco/JugReco/Track.hpp
index 595597599a1cc43b7b57123541ae3ba3ff0be521..57c3452cbea34ca76870dff618857a50cfc858b3 100644
--- a/JugReco/JugReco/Track.hpp
+++ b/JugReco/JugReco/Track.hpp
@@ -15,21 +15,19 @@
//#include "ACTFW/EventData/SimSourceLink.hpp"
#include "Acts/EventData/MultiTrajectory.hpp"
#include "Acts/EventData/TrackParameters.hpp"
-#include "JugReco/SourceLinks.h"
+#include "JugReco/IndexSourceLink.hpp"
#include <vector>
namespace Jug {
-
/// (Reconstructed) track parameters e.g. close to the vertex.
- using TrackParameters = Acts::CurvilinearTrackParameters;
-
+ using TrackParameters = ::Acts::BoundTrackParameters;
/// Container of reconstructed track states for multiple tracks.
using TrackParametersContainer = std::vector<TrackParameters>;
/// MultiTrajectory definition
- using Trajectory = Acts::MultiTrajectory<SourceLink>;
+ using Trajectory = Acts::MultiTrajectory<IndexSourceLink>;
/// Container for the truth fitting/finding track(s)
using TrajectoryContainer = std::vector<SimMultiTrajectory>;
diff --git a/JugReco/src/components/ParticlesFromTrackFit.cpp b/JugReco/src/components/ParticlesFromTrackFit.cpp
index 56b424be3ce4c26df72c0f8df8bd25066793b234..7b43f691de5388a07ee4cf46cb6ea11d2874356b 100644
--- a/JugReco/src/components/ParticlesFromTrackFit.cpp
+++ b/JugReco/src/components/ParticlesFromTrackFit.cpp
@@ -22,7 +22,7 @@
#include "eicd/ParticleCollection.h"
#include "eicd/TrackerHitCollection.h"
#include "eicd/TrackParametersCollection.h"
-#include "JugReco/SourceLinks.h"
+#include "JugReco/IndexSourceLink.hpp"
#include "JugReco/Track.hpp"
#include "Acts/Utilities/Helpers.hpp"
diff --git a/JugReco/src/components/TrackFindingAlgorithm.cpp b/JugReco/src/components/TrackFindingAlgorithm.cpp
index 52a6bacb3ef12b2c6efc4b2fd6151205c47cbaa6..c25e87666be7eb23693cfa21a1b373de581299a5 100644
--- a/JugReco/src/components/TrackFindingAlgorithm.cpp
+++ b/JugReco/src/components/TrackFindingAlgorithm.cpp
@@ -26,18 +26,20 @@
#include "Acts/Propagator/EigenStepper.hpp"
#include "Acts/Propagator/Navigator.hpp"
#include "Acts/Propagator/Propagator.hpp"
-#include "Acts/Utilities/Definitions.hpp"
+#include "Acts/Definitions/Common.hpp"
#include "Acts/Utilities/Helpers.hpp"
#include "Acts/Utilities/Logger.hpp"
-#include "Acts/Utilities/Units.hpp"
+#include "Acts/Definitions/Units.hpp"
#include "JugBase/DataHandle.h"
#include "JugBase/IGeoSvc.h"
#include "JugReco/GeometryContainers.hpp"
-#include "JugReco/SourceLinks.h"
+//#include "JugReco/SourceLinks.h"
+#include "JugReco/Measurement.hpp"
+#include "JugReco/Index.hpp"
+#include "JugReco/IndexSourceLink.hpp"
#include "JugReco/Track.hpp"
#include "JugReco/BField.h"
-#include "JugReco/SourceLinks.h"
#include "eicd/TrackerHitCollection.h"
@@ -47,19 +49,21 @@
#include <random>
#include <stdexcept>
-
-
namespace Jug::Reco {
using namespace Acts::UnitLiterals;
- TrackFindingAlgorithm::TrackFindingAlgorithm(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) {
+ TrackFindingAlgorithm::TrackFindingAlgorithm(const std::string& name, ISvcLocator* svcLoc)
+ : GaudiAlgorithm(name, svcLoc)
+ {
declareProperty("inputSourceLinks", m_inputSourceLinks, "");
+ declareProperty("inputMeasurements", m_inputMeasurements, "");
declareProperty("inputInitialTrackParameters", m_inputInitialTrackParameters, "");
declareProperty("outputTrajectories", m_outputTrajectories, "");
}
- StatusCode TrackFindingAlgorithm::initialize() {
+ StatusCode TrackFindingAlgorithm::initialize()
+ {
if (GaudiAlgorithm::initialize().isFailure())
return StatusCode::FAILURE;
m_geoSvc = service("GeoSvc");
@@ -68,24 +72,22 @@ namespace Jug::Reco {
<< "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
return StatusCode::FAILURE;
}
- m_BField = std::make_shared<Acts::ConstantBField>(Acts::Vector3D{0.0, 0.0, m_geoSvc->centralMagneticField()});
- m_fieldctx = BFieldVariant(m_BField);
-
+ m_BField = std::make_shared<Acts::ConstantBField>(Acts::Vector3{0.0, 0.0, m_geoSvc->centralMagneticField()});
+ m_fieldctx = BFieldVariant(m_BField);
// chi2 and #sourclinks per surface cutoffs
- m_sourcelinkSelectorCfg = { {Acts::GeometryIdentifier(), {15, 10}},
+ m_sourcelinkSelectorCfg = {
+ {Acts::GeometryIdentifier(), {15, 10}},
};
-
- findTracks = TrackFindingAlgorithm::makeTrackFinderFunction(m_geoSvc->trackingGeometry(), m_BField);
-
+ m_trackFinderFunc = TrackFindingAlgorithm::makeTrackFinderFunction(m_geoSvc->trackingGeometry(), m_BField);
return StatusCode::SUCCESS;
}
- StatusCode TrackFindingAlgorithm::execute() {
+ StatusCode TrackFindingAlgorithm::execute()
+ {
// Read input data
- const SourceLinkContainer* src_links = m_inputSourceLinks.get();
+ const IndexSourceLinkContainer* src_links = m_inputSourceLinks.get();
const TrackParametersContainer* init_trk_params = m_inputInitialTrackParameters.get();
- // const auto sourceLinks = ctx.eventStore.get<SourceLinkContainer>(m_cfg.inputSourceLinks);
- // const auto initialParameters = ctx.eventStore.get<TrackParametersContainer>(m_cfg.inputInitialTrackParameters);
+ const MeasurementContainer* measurements = m_inputMeasurements.get();
//// Prepare the output data with MultiTrajectory
// TrajectoryContainer trajectories;
@@ -93,7 +95,7 @@ namespace Jug::Reco {
trajectories->reserve(init_trk_params->size());
//// Construct a perigee surface as the target surface
- auto pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3D{0., 0., 0.});
+ auto pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3{0., 0., 0.});
ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("TrackFindingAlgorithm Logger", Acts::Logging::INFO));
@@ -104,13 +106,14 @@ namespace Jug::Reco {
Acts::PropagatorPlainOptions pOptions;
pOptions.maxSteps = 10000;
+
// Set the CombinatorialKalmanFilter options
- CKFOptions ckfOptions(m_geoctx, m_fieldctx, m_calibctx, m_sourcelinkSelectorCfg,
- Acts::LoggerWrapper{logger()}, pOptions,
- &(*pSurface));
+ TrackFindingAlgorithm::TrackFinderOptions options(
+ m_geoctx, m_fieldctx, m_calibctx, MeasurementCalibrator(*measurements),
+ Acts::CKFSourceLinkSelector(m_sourcelinkSelectorCfg), Acts::LoggerWrapper{logger()}, pOptions, &(*pSurface));
debug() << "Invoke track finding seeded by truth particle " << iseed << endmsg;
- auto result = findTracks(*src_links, initialParams, ckfOptions);
+ auto result = m_trackFinderFunc(*src_links, initialParams, options);
debug() << "finding done." << endmsg;
if (result.ok()) {
// Get the track finding output object
@@ -122,7 +125,7 @@ namespace Jug::Reco {
debug() << "Track finding failed for truth seed " << iseed << endmsg;
ACTS_WARNING("Track finding failed for truth seed " << iseed << " with error" << result.error());
// Track finding failed, but still create an empty SimMultiTrajectory
- //trajectories->push_back(SimMultiTrajectory());
+ // trajectories->push_back(SimMultiTrajectory());
}
}
@@ -141,8 +144,9 @@ namespace {
TrackFinderFunctionImpl(TrackFinder&& f) : trackFinder(std::move(f)) {}
Jug::Reco::TrackFindingAlgorithm::TrackFinderResult
- operator()(const Jug::SourceLinkContainer& sourceLinks, const Jug::TrackParameters& initialParameters,
- const Acts::CombinatorialKalmanFilterOptions<Acts::CKFSourceLinkSelector>& options) const {
+ operator()(const Jug::IndexSourceLinkContainer& sourceLinks, const Jug::TrackParameters& initialParameters,
+ const Jug::Reco::TrackFindingAlgorithm::TrackFinderOptions& options) const
+ {
return trackFinder.findTracks(sourceLinks, initialParameters, options);
};
};
@@ -169,7 +173,7 @@ namespace Jug::Reco {
using Navigator = Acts::Navigator;
using Propagator = Acts::Propagator<Stepper, Navigator>;
using SourceLinkSelector = Acts::CKFSourceLinkSelector;
- using CKF = Acts::CombinatorialKalmanFilter<Propagator, Updater, Smoother, SourceLinkSelector>;
+ using CKF = Acts::CombinatorialKalmanFilter<Propagator, Updater, Smoother>;
std::cout << " finding ...\n";
// construct all components for the track finder
@@ -179,14 +183,14 @@ namespace Jug::Reco {
navigator.resolvePassive = false;
navigator.resolveMaterial = true;
navigator.resolveSensitive = true;
- std::cout << " propagator\n" ;
+ std::cout << " propagator\n";
Propagator propagator(std::move(stepper), std::move(navigator));
CKF trackFinder(std::move(propagator));
+
// build the track finder functions. owns the track finder object.
return TrackFinderFunctionImpl<CKF>(std::move(trackFinder));
},
std::move(magneticField));
}
} // namespace Jug::Reco
-
diff --git a/JugReco/src/components/TrackFindingAlgorithm.h b/JugReco/src/components/TrackFindingAlgorithm.h
index 73d5054da23f1995dab5600dadc0ab1d001ef464..47cd1dc3dd5be7ee6180f736f5954dd9dc634eae 100644
--- a/JugReco/src/components/TrackFindingAlgorithm.h
+++ b/JugReco/src/components/TrackFindingAlgorithm.h
@@ -16,7 +16,7 @@
//#include "Acts/Geometry/TrackingGeometry.hpp"
//#include "Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp"
-//#include "Acts/Utilities/Definitions.hpp"
+//#include "Acts/Definitions/Common.hpp"
//#include "Acts/Utilities/Helpers.hpp"
//#include "Acts/Utilities/Logger.hpp"
@@ -24,9 +24,11 @@
#include <stdexcept>
#include <vector>
-#include "JugReco/SourceLinks.h"
+#include "JugReco/Index.hpp"
+#include "JugReco/IndexSourceLink.hpp"
#include "JugReco/Track.hpp"
#include "JugReco/BField.h"
+#include "JugReco/Measurement.hpp"
#include "eicd/TrackerHitCollection.h"
@@ -34,17 +36,12 @@
#include "Acts/Geometry/TrackingGeometry.hpp"
#include "Acts/TrackFinding/CKFSourceLinkSelector.hpp"
#include "Acts/TrackFinding/CombinatorialKalmanFilter.hpp"
-#include "Acts/Utilities/Definitions.hpp"
-
-//#include "Acts/Fitter/GainMatrixSmoother.hpp"
-//#include "Acts/Fitter/GainMatrixUpdater.hpp"
-//#include "Acts/MagneticField/ConstantBField.hpp"
-//#include "Acts/MagneticField/InterpolatedBFieldMap.hpp"
-//#include "Acts/MagneticField/SharedBField.hpp"
-//#include "Acts/Propagator/EigenStepper.hpp"
-//#include "Acts/Propagator/Navigator.hpp"
-//#include "Acts/Propagator/Propagator.hpp"
-//#include "Acts/Utilities/Units.hpp"
+#include "Acts/Definitions/Common.hpp"
+
+#include "Acts/Geometry/TrackingGeometry.hpp"
+#include "Acts/TrackFinding/CKFSourceLinkSelector.hpp"
+#include "Acts/TrackFinding/CombinatorialKalmanFilter.hpp"
+
#include <random>
#include <stdexcept>
@@ -53,17 +50,23 @@ namespace Jug::Reco {
class TrackFindingAlgorithm : public GaudiAlgorithm {
public:
- using TrackFinderResult = Acts::Result<Acts::CombinatorialKalmanFilterResult<SourceLink>>;
-
- /// Track finding function that takes input measurements, initial trackstate
- /// and track finder options and returns some track-finding-specific result.
- using CKFOptions = Acts::CombinatorialKalmanFilterOptions<Acts::CKFSourceLinkSelector>;
-
- using TrackFinderFunction =
- std::function<TrackFinderResult(const SourceLinkContainer&, const TrackParameters&, const CKFOptions&)>;
+ /// Track finder function that takes input measurements, initial trackstate
+ /// and track finder options and returns some track-finder-specific result.
+ using TrackFinderOptions = Acts::CombinatorialKalmanFilterOptions<MeasurementCalibrator, Acts::CKFSourceLinkSelector>;
+ using TrackFinderResult = Acts::Result<Acts::CombinatorialKalmanFilterResult<IndexSourceLink>>;
+ using TrackFinderFunction = std::function<TrackFinderResult(
+ const IndexSourceLinkContainer&, const TrackParameters&, const TrackFinderOptions&)>;
+
+ /// Create the track finder function implementation.
+ ///
+ /// The magnetic field is intentionally given by-value since the variant
+ /// contains shared_ptr anyways.
+ static TrackFinderFunction makeTrackFinderFunction(std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry,
+ BFieldVariant magneticField);
public:
- DataHandle<SourceLinkContainer> m_inputSourceLinks{"inputSourceLinks", Gaudi::DataHandle::Reader, this};
+ DataHandle<IndexSourceLinkContainer> m_inputSourceLinks{"inputSourceLinks", Gaudi::DataHandle::Reader, this};
+ DataHandle<MeasurementContainer> m_inputMeasurements{"inputMeasurements", Gaudi::DataHandle::Reader, this};
DataHandle<TrackParametersContainer> m_inputInitialTrackParameters{"inputInitialTrackParameters",
Gaudi::DataHandle::Reader, this};
DataHandle<TrajectoryContainer> m_outputTrajectories{"outputTrajectories", Gaudi::DataHandle::Writer, this};
@@ -79,15 +82,7 @@ namespace Jug::Reco {
TrackFindingAlgorithm(const std::string& name, ISvcLocator* svcLoc);
- /** Create the track finder function implementation.
- * The magnetic field is intentionally given by-value since the variant
- * contains shared_ptr anyways.
- */
- static TrackFinderFunction makeTrackFinderFunction(std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry,
- BFieldVariant magneticField);
- /// Type erased track finder function.
- TrackFinderFunction findTracks;
StatusCode initialize() override;
diff --git a/JugReco/src/components/TrackFittingAlgorithm.cpp b/JugReco/src/components/TrackFittingAlgorithm.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2a2d327a17687f00d5156f723a6b7e739068691f
--- /dev/null
+++ b/JugReco/src/components/TrackFittingAlgorithm.cpp
@@ -0,0 +1,199 @@
+//
+#include "TrackFittingAlgorithm.h"
+
+// Gaudi
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "GaudiKernel/ToolHandle.h"
+#include "GaudiAlg/Transformer.h"
+#include "GaudiAlg/GaudiTool.h"
+#include "GaudiKernel/RndmGenerators.h"
+#include "GaudiKernel/Property.h"
+
+#include "DDRec/CellIDPositionConverter.h"
+#include "DDRec/SurfaceManager.h"
+#include "DDRec/Surface.h"
+
+#include "Acts/Geometry/TrackingGeometry.hpp"
+#include "Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp"
+#include "Acts/Surfaces/PerigeeSurface.hpp"
+#include "Acts/MagneticField/ConstantBField.hpp"
+#include "Acts/MagneticField/InterpolatedBFieldMap.hpp"
+#include "Acts/MagneticField/SharedBField.hpp"
+#include "Acts/Propagator/EigenStepper.hpp"
+#include "Acts/Propagator/Navigator.hpp"
+#include "Acts/Propagator/Propagator.hpp"
+#include "Acts/Definitions/Common.hpp"
+#include "Acts/Utilities/Helpers.hpp"
+#include "Acts/Utilities/Logger.hpp"
+#include "Acts/Definitions/Units.hpp"
+
+#include "JugBase/DataHandle.h"
+#include "JugBase/IGeoSvc.h"
+#include "JugReco/GeometryContainers.hpp"
+#include "JugReco/SourceLinks.h"
+#include "JugReco/Track.hpp"
+#include "JugReco/BField.h"
+#include "JugReco/Measurement.hpp"
+#include "JugReco/SourceLinks.h"
+
+#include "eicd/TrackerHitCollection.h"
+
+#include <functional>
+#include <stdexcept>
+#include <vector>
+#include <random>
+#include <stdexcept>
+
+namespace Jug::Reco {
+
+ using namespace Acts::UnitLiterals;
+
+ TrackFittingAlgorithm::TrackFittingAlgorithm(const std::string& name, ISvcLocator* svcLoc)
+ : GaudiAlgorithm(name, svcLoc)
+ {
+ declareProperty("inputSourceLinks", m_inputSourceLinks, "");
+ declareProperty("initialTrackParameters", m_initialTrackParameters, "");
+ declareProperty("foundTracks", m_foundTracks, "");
+ declareProperty("outputTrajectories", m_outputTrajectories, "");
+ }
+
+ StatusCode TrackFittingAlgorithm::initialize()
+ {
+ if (GaudiAlgorithm::initialize().isFailure())
+ return StatusCode::FAILURE;
+ m_geoSvc = service("GeoSvc");
+ if (!m_geoSvc) {
+ error() << "Unable to locate Geometry Service. "
+ << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
+ return StatusCode::FAILURE;
+ }
+ m_BField = std::make_shared<Acts::ConstantBField>(Acts::Vector3{0.0, 0.0, m_geoSvc->centralMagneticField()});
+ m_fieldctx = BFieldVariant(m_BField);
+ // chi2 and #sourclinks per surface cutoffs
+ //m_sourcelinkSelectorCfg = {
+ // {Acts::GeometryIdentifier(), {15, 10}},
+ //};
+ m_trackFittingFunc = makeTrackFittingFunction(m_geoSvc->trackingGeometry(), m_BField);
+ return StatusCode::SUCCESS;
+ }
+
+ StatusCode TrackFittingAlgorithm::execute()
+ {
+ // Read input data
+ const SourceLinkContainer* src_links = m_inputSourceLinks.get();
+ const TrackParametersContainer* init_trk_params = m_initialTrackParameters.get();
+
+ // TrajectoryContainer trajectories;
+ auto trajectories = m_outputTrajectories.createAndPut();
+ trajectories->reserve(init_trk_params->size());
+
+ //// Construct a perigee surface as the target surface
+ auto pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3{0., 0., 0.});
+ ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("TrackFittingAlgorithm Logger", Acts::Logging::INFO));
+
+ // Perform the track finding for each starting parameter
+ // @TODO: use seeds from track seeding algorithm as starting parameter
+ for (std::size_t iseed = 0; iseed < init_trk_params->size(); ++iseed) {
+
+ const auto& initialParams = (*init_trk_params)[iseed];
+ Acts::PropagatorPlainOptions pOptions;
+ pOptions.maxSteps = 10000;
+
+ Acts::KalmanFitterOptions<Acts::VoidOutlierFinder> kfOptions(
+ m_geoctx, m_fieldctx, m_calibctx,
+ Acts::VoidOutlierFinder(), Acts::LoggerWrapper{logger()}, Acts::PropagatorPlainOptions(), &(*pSurface));
+
+ debug() << "Invoke track fitting ... " << iseed << endmsg;
+ auto result = m_trackFittingFunc(*src_links, initialParams, kfOptions);
+ debug() << "fitting done." << endmsg;
+ if (result.ok()) {
+ // Get the track finding output object
+ const auto& trackFindingOutput = result.value();
+ // Create a SimMultiTrajectory
+ trajectories->emplace_back(std::move(trackFindingOutput.fittedStates), std::move(trackFindingOutput.trackTips),
+ std::move(trackFindingOutput.fittedParameters));
+ } else {
+ debug() << "Track finding failed for truth seed " << iseed << endmsg;
+ ACTS_WARNING("Track finding failed for truth seed " << iseed << " with error" << result.error());
+ // Track finding failed, but still create an empty SimMultiTrajectory
+ // trajectories->push_back(SimMultiTrajectory());
+ }
+ }
+
+ // ctx.eventStore.add(m_cfg.outputTrajectories, std::move(trajectories));
+ return StatusCode::SUCCESS;
+
+
+ ///////////////////////////
+ // acts example
+
+ // Set the KalmanFitter options
+
+ // Perform the fit for each input track
+ //std::vector<IndexSourceLink> trackSourceLinks;
+ //for (std::size_t itrack = 0; itrack < protoTracks.size(); ++itrack) {
+ // // The list of hits and the initial start parameters
+ // const auto& protoTrack = protoTracks[itrack];
+ // const auto& initialParams = initialParameters[itrack];
+
+ // // We can have empty tracks which must give empty fit results so the number
+ // // of entries in input and output containers matches.
+ // if (protoTrack.empty()) {
+ // trajectories.push_back(Trajectories());
+ // ACTS_WARNING("Empty track " << itrack << " found.");
+ // continue;
+ // }
+
+ // // Clear & reserve the right size
+ // trackSourceLinks.clear();
+ // trackSourceLinks.reserve(protoTrack.size());
+
+ // // Fill the source links via their indices from the container
+ // for (auto hitIndex : protoTrack) {
+ // auto sourceLink = sourceLinks.nth(hitIndex);
+ // if (sourceLink == sourceLinks.end()) {
+ // ACTS_FATAL("Proto track " << itrack << " contains invalid hit index"
+ // << hitIndex);
+ // return ProcessCode::ABORT;
+ // }
+ // trackSourceLinks.push_back(*sourceLink);
+ // }
+
+ // ACTS_DEBUG("Invoke fitter");
+ // auto result = m_cfg.fit(trackSourceLinks, initialParams, kfOptions);
+ // if (result.ok()) {
+ // // Get the fit output object
+ // const auto& fitOutput = result.value();
+ // // The track entry indices container. One element here.
+ // std::vector<size_t> trackTips;
+ // trackTips.reserve(1);
+ // trackTips.emplace_back(fitOutput.trackTip);
+ // // The fitted parameters container. One element (at most) here.
+ // Trajectories::IndexedParameters indexedParams;
+ // if (fitOutput.fittedParameters) {
+ // const auto& params = fitOutput.fittedParameters.value();
+ // ACTS_VERBOSE("Fitted paramemeters for track " << itrack);
+ // ACTS_VERBOSE(" " << params.parameters().transpose());
+ // // Push the fitted parameters to the container
+ // indexedParams.emplace(fitOutput.trackTip, std::move(params));
+ // } else {
+ // ACTS_DEBUG("No fitted paramemeters for track " << itrack);
+ // }
+ // // store the result
+ // trajectories.emplace_back(std::move(fitOutput.fittedStates),
+ // std::move(trackTips), std::move(indexedParams));
+ // } else {
+ // ACTS_WARNING("Fit failed for track " << itrack << " with error"
+ // << result.error());
+ // // Fit failed. Add an empty result so the output container has
+ // // the same number of entries as the input.
+ // trajectories.push_back(Trajectories());
+ // }
+ //}
+
+
+ return StatusCode::SUCCESS;
+ }
+
+ DECLARE_COMPONENT(TrackFittingAlgorithm)
+} // namespace Jug::Reco
diff --git a/JugReco/src/components/TrackFittingAlgorithm.h b/JugReco/src/components/TrackFittingAlgorithm.h
new file mode 100644
index 0000000000000000000000000000000000000000..1386598b1b49d29744ce8f8a9a4f51a1fab4ab62
--- /dev/null
+++ b/JugReco/src/components/TrackFittingAlgorithm.h
@@ -0,0 +1,110 @@
+#ifndef JUGGLER_JUGRECO_TrackFittingAlgorithm_HH
+#define JUGGLER_JUGRECO_TrackFittingAlgorithm_HH
+
+#include "JugReco/GeometryContainers.hpp"
+
+// Gaudi
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "GaudiKernel/ToolHandle.h"
+//#include "GaudiAlg/Transformer.h"
+//#include "GaudiAlg/GaudiTool.h"
+//#include "GaudiKernel/RndmGenerators.h"
+#include "Gaudi/Property.h"
+
+#include "JugBase/DataHandle.h"
+#include "JugBase/IGeoSvc.h"
+
+//#include "Acts/Geometry/TrackingGeometry.hpp"
+//#include "Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp"
+//#include "Acts/Definitions/Common.hpp"
+//#include "Acts/Utilities/Helpers.hpp"
+//#include "Acts/Utilities/Logger.hpp"
+
+#include <functional>
+#include <stdexcept>
+#include <vector>
+
+#include "JugReco/SourceLinks.h"
+#include "JugReco/Track.hpp"
+#include "JugReco/BField.h"
+#include "JugReco/Measurement.hpp"
+
+#include "eicd/TrackerHitCollection.h"
+
+//#include "Acts/Surfaces/PerigeeSurface.hpp"
+
+#include "Acts/TrackFitting/KalmanFitter.hpp"
+#include "Acts/TrackFitting/GainMatrixSmoother.hpp"
+#include "Acts/TrackFitting/GainMatrixUpdater.hpp"
+#include "Acts/Geometry/TrackingGeometry.hpp"
+#include "Acts/Definitions/Common.hpp"
+
+//#include "Acts/Fitter/GainMatrixSmoother.hpp"
+//#include "Acts/Fitter/GainMatrixUpdater.hpp"
+//#include "Acts/MagneticField/ConstantBField.hpp"
+//#include "Acts/MagneticField/InterpolatedBFieldMap.hpp"
+//#include "Acts/MagneticField/SharedBField.hpp"
+//#include "Acts/Propagator/EigenStepper.hpp"
+//#include "Acts/Propagator/Navigator.hpp"
+//#include "Acts/Propagator/Propagator.hpp"
+//#include "Acts/Definitions/Units.hpp"
+
+#include <random>
+#include <stdexcept>
+
+namespace Jug::Reco {
+
+ class TrackFittingAlgorithm : public GaudiAlgorithm {
+ public:
+ //using TrackFinderResult = Acts::Result<Acts::CombinatorialKalmanFilterResult<SourceLink>>;
+ using FitterResult = Acts::Result<Acts::KalmanFitterResult<SourceLink>>;
+ /// Fit function that takes input measurements, initial trackstate and fitter
+
+ using FitterFunction = std::function<FitterResult(
+ const std::vector<SourceLink>&, const TrackParameters&,
+ const Acts::KalmanFitterOptions<Acts::VoidOutlierFinder>&)>;
+
+
+ /// Track fitter function that takes input measurements, initial trackstate
+ /// and fitter options and returns some track-fitter-specific result.
+ using TrackFitterOptions =
+ Acts::KalmanFitterOptions< Acts::VoidOutlierFinder>;
+ using TrackFitterResult =
+ Acts::Result<Acts::KalmanFitterResult<SourceLink>>;
+ using TrackFitterFunction = std::function<TrackFitterResult(
+ const std::vector<SourceLink>&, const TrackParameters&,
+ const TrackFitterOptions&)>;
+
+ public:
+ DataHandle<SourceLinkContainer> m_inputSourceLinks{"inputSourceLinks", Gaudi::DataHandle::Reader, this};
+ DataHandle<TrackParametersContainer> m_initialTrackParameters{"initialTrackParameters", Gaudi::DataHandle::Reader, this};
+ DataHandle<TrajectoryContainer> m_foundTracks{"foundTracks", Gaudi::DataHandle::Reader, this};
+ DataHandle<TrajectoryContainer> m_outputTrajectories{"outputTrajectories", Gaudi::DataHandle::Writer, this};
+
+ FitterFunction m_trackFittingFunc;
+ SmartIF<IGeoSvc> m_geoSvc;
+ std::shared_ptr<Acts::ConstantBField> m_BField = nullptr;
+ Acts::GeometryContext m_geoctx;
+ Acts::CalibrationContext m_calibctx;
+ Acts::MagneticFieldContext m_fieldctx;
+
+ //Acts::CKFSourceLinkSelector::Config m_sourcelinkSelectorCfg;
+
+ TrackFittingAlgorithm(const std::string& name, ISvcLocator* svcLoc);
+
+ /** Create the track finder function implementation.
+ * The magnetic field is intentionally given by-value since the variant
+ * contains shared_ptr anyways.
+ */
+ static FitterFunction makeTrackFittingFunction(std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry,
+ BFieldVariant magneticField);
+
+
+ StatusCode initialize() override;
+
+ StatusCode execute() override;
+ };
+
+} // namespace Jug::Reco
+
+#endif
diff --git a/JugReco/src/components/TrackFittingFunction.cpp b/JugReco/src/components/TrackFittingFunction.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9eafb46ffba71f01639e3a6368322266e9d02bc8
--- /dev/null
+++ b/JugReco/src/components/TrackFittingFunction.cpp
@@ -0,0 +1,77 @@
+#include "Acts/Geometry/GeometryIdentifier.hpp"
+#include "Acts/Geometry/TrackingGeometry.hpp"
+#include "Acts/MagneticField/ConstantBField.hpp"
+#include "Acts/MagneticField/InterpolatedBFieldMap.hpp"
+#include "Acts/MagneticField/SharedBField.hpp"
+#include "Acts/Propagator/EigenStepper.hpp"
+#include "Acts/Propagator/Navigator.hpp"
+#include "Acts/Propagator/Propagator.hpp"
+#include "Acts/Surfaces/Surface.hpp"
+#include "Acts/TrackFitting/GainMatrixSmoother.hpp"
+#include "Acts/TrackFitting/GainMatrixUpdater.hpp"
+#include "Acts/Utilities/Helpers.hpp"
+#include "Acts/Utilities/ParameterDefinitions.hpp"
+//#include "ActsExamples/Plugins/BField/ScalableBField.hpp"
+#include "TrackFittingAlgorithm.h"
+
+namespace {
+
+ template <typename track_fitter_t>
+ struct TrackFitterFunctionImpl {
+ track_fitter_t trackFitter;
+
+ TrackFitterFunctionImpl(track_fitter_t&& f) : trackFitter(std::move(f)) {}
+
+ Jug::Reco::TrackFittingAlgorithm::FitterResult operator()(
+ const std::vector<Jug::SourceLink>& sourceLinks,
+ const Jug::TrackParameters& initialParameters,
+ const Acts::KalmanFitterOptions<Acts::VoidOutlierFinder>& options) const {
+ return fitter.fit(sourceLinks, initialParameters, options);
+ }
+
+ //TrackFittingAlgorithm::TrackFitterResult
+ //operator()(const std::vector<SourceLink>& sourceLinks,
+ // const TrackParameters& initialParameters,
+ // const TrackFittingAlgorithm::TrackFitterOptions& options) const
+ //{
+ // return trackFitter.fit(sourceLinks, initialParameters, options);
+ //};
+ };
+
+} // namespace
+
+namespace Jug::Reco {
+ TrackFittingAlgorithm::TrackFitterFunction TrackFittingAlgorithm::makeTrackFittingFunction(
+ std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry, Options::BFieldVariant magneticField)
+ {
+ using Updater = Acts::GainMatrixUpdater;
+ using Smoother = Acts::GainMatrixSmoother;
+
+ // unpack the magnetic field variant and instantiate the corresponding fitter.
+ return std::visit(
+ [trackingGeometry](auto&& inputField) -> TrackFitterFunction {
+ // each entry in the variant is already a shared_ptr
+ // need ::element_type to get the real magnetic field type
+ using InputMagneticField = typename std::decay_t<decltype(inputField)>::element_type;
+ using MagneticField = Acts::SharedBField<InputMagneticField>;
+ using Stepper = Acts::EigenStepper<MagneticField>;
+ using Navigator = Acts::Navigator;
+ using Propagator = Acts::Propagator<Stepper, Navigator>;
+ using Fitter = Acts::KalmanFitter<Propagator, Updater, Smoother>;
+
+ // construct all components for the fitter
+ MagneticField field(std::move(inputField));
+ Stepper stepper(std::move(field));
+ Navigator navigator(trackingGeometry);
+ navigator.resolvePassive = false;
+ navigator.resolveMaterial = true;
+ navigator.resolveSensitive = true;
+ Propagator propagator(std::move(stepper), std::move(navigator));
+ Fitter trackFitter(std::move(propagator));
+
+ // build the fitter functions. owns the fitter object.
+ return TrackFitterFunctionImpl<Fitter>(std::move(trackFitter));
+ },
+ std::move(magneticField));
+ }
+} // namespace Jug::Reco
diff --git a/JugReco/src/components/TrackParamClusterInit.cpp b/JugReco/src/components/TrackParamClusterInit.cpp
index cd8010fdbae0db107de627e6fde213aeb6e66995..a9715615fabdb3002997c7c3c0313a21b9338c3b 100644
--- a/JugReco/src/components/TrackParamClusterInit.cpp
+++ b/JugReco/src/components/TrackParamClusterInit.cpp
@@ -10,12 +10,14 @@
#include "JugBase/DataHandle.h"
#include "JugBase/IGeoSvc.h"
#include "JugReco/Track.hpp"
-#include "Acts/Utilities/Units.hpp"
-#include "Acts/Utilities/Definitions.hpp"
+#include "Acts/Definitions/Units.hpp"
+#include "Acts/Definitions/Common.hpp"
#include "eicd/TrackerHitCollection.h"
#include "eicd/ClusterCollection.h"
+#include "Math/Vector3D.h"
+#include "Acts/Surfaces/PerigeeSurface.hpp"
///// (Reconstructed) track parameters e.g. close to the vertex.
//using TrackParameters = Acts::CurvilinearTrackParameters;
@@ -77,28 +79,52 @@ namespace Jug::Reco {
continue;
}
double len = std::hypot( c.x() , c.y() , c.z() );
+ ROOT::Math::XYZVector momentum(c.x() * p / len, c.y() * p / len, c.z() * p / len);
// build some track cov matrix
- Acts::BoundSymMatrix cov = Acts::BoundSymMatrix::Zero();
- cov(Acts::eBoundLoc0, Acts::eBoundLoc0) = 1.0 * mm*1.0 * mm;
- cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = 1.0 * mm*1.0 * mm;
- cov(Acts::eBoundPhi, Acts::eBoundPhi) = M_PI / 180.0;
- cov(Acts::eBoundTheta, Acts::eBoundTheta) = M_PI / 180.0;
- cov(Acts::eBoundQOverP, Acts::eBoundQOverP) = 1.0 / (p*p);
- cov(Acts::eBoundTime, Acts::eBoundTime) = Acts::UnitConstants::ns;
+ //Acts::BoundSymMatrix cov = Acts::BoundSymMatrix::Zero();
+ //cov(Acts::eBoundLoc0, Acts::eBoundLoc0) = 1.0 * mm*1.0 * mm;
+ //cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = 1.0 * mm*1.0 * mm;
+ //cov(Acts::eBoundPhi, Acts::eBoundPhi) = M_PI / 180.0;
+ //cov(Acts::eBoundTheta, Acts::eBoundTheta) = M_PI / 180.0;
+ //cov(Acts::eBoundQOverP, Acts::eBoundQOverP) = 1.0 / (p*p);
+ //cov(Acts::eBoundTime, Acts::eBoundTime) = Acts::UnitConstants::ns;
+
+ Acts::BoundVector params;
+ params(Acts::eBoundLoc0) = 0.0 * mm ;
+ params(Acts::eBoundLoc1) = 0.0 * mm ;
+ params(Acts::eBoundPhi) = momentum.Phi();
+ params(Acts::eBoundTheta) = momentum.Theta();
+ params(Acts::eBoundQOverP) = 1/p;
+ params(Acts::eBoundTime) = 0 * ns;
+
+ auto pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3{0,0,0});
debug() << "Invoke track finding seeded by truth particle with p = " << p/GeV << " GeV" << endmsg;
- // add all charges to the track candidate...
- init_trk_params->emplace_back(Acts::Vector4D(0 * mm, 0 * mm, 0 * mm, 0),
- Acts::Vector3D(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, -1,
- std::make_optional(cov));
- debug() << init_trk_params->back() << endmsg;
- init_trk_params->emplace_back(Acts::Vector4D(0 * mm, 0 * mm, 0 * mm, 0),
- Acts::Vector3D(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, 1,
- std::make_optional(cov));
- debug() << init_trk_params->back() << endmsg;
- //init_trk_params->emplace_back(Acts::Vector4D(0 * mm, 0 * mm, 0 * mm, 0),
- // Acts::Vector3D(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, 0,
+
+ // add both charges to the track candidate...
+ init_trk_params->push_back({pSurface, params, 1});
+
+ Acts::BoundVector params2;
+ params2(Acts::eBoundLoc0) = 0.0 * mm ;
+ params2(Acts::eBoundLoc1) = 0.0 * mm ;
+ params2(Acts::eBoundPhi) = momentum.Phi();
+ params2(Acts::eBoundTheta) = momentum.Theta();
+ params2(Acts::eBoundQOverP) = -1/p;
+ params2(Acts::eBoundTime) = 0 * ns;
+ init_trk_params->push_back({pSurface, params2, -1});
+
+ // acts v1.2.0:
+ //init_trk_params->emplace_back(Acts::Vector4(0 * mm, 0 * mm, 0 * mm, 0),
+ // Acts::Vector3(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, -1,
+ // std::make_optional(cov));
+ //debug() << init_trk_params->back() << endmsg;
+ //init_trk_params->emplace_back(Acts::Vector4(0 * mm, 0 * mm, 0 * mm, 0),
+ // Acts::Vector3(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, 1,
+ // std::make_optional(cov));
+ ////debug() << init_trk_params->back() << endmsg;
+ //init_trk_params->emplace_back(Acts::Vector4(0 * mm, 0 * mm, 0 * mm, 0),
+ // Acts::Vector3(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, 0,
// std::make_optional(cov));
}
return StatusCode::SUCCESS;
diff --git a/JugReco/src/components/TrackParamTruthInit.cpp b/JugReco/src/components/TrackParamTruthInit.cpp
index f12c86aeb5f462e6d4eacd164e7c21df5cd26d8c..34a8077c57454a72d7370056b3a77ec9aec86eef 100644
--- a/JugReco/src/components/TrackParamTruthInit.cpp
+++ b/JugReco/src/components/TrackParamTruthInit.cpp
@@ -10,11 +10,13 @@
#include "JugBase/DataHandle.h"
#include "JugBase/IGeoSvc.h"
#include "JugReco/Track.hpp"
-#include "Acts/Utilities/Units.hpp"
-#include "Acts/Utilities/Definitions.hpp"
+#include "Acts/Definitions/Units.hpp"
+#include "Acts/Definitions/Common.hpp"
#include "eicd/TrackerHitCollection.h"
#include "dd4pod/Geant4ParticleCollection.h"
+#include "Math/Vector3D.h"
+#include "Acts/Surfaces/PerigeeSurface.hpp"
///// (Reconstructed) track parameters e.g. close to the vertex.
@@ -70,11 +72,13 @@ namespace Jug::Reco {
if(part.genStatus() != 1 ) {
continue;
}
- using Acts::UnitConstants::MeV;
using Acts::UnitConstants::GeV;
+ using Acts::UnitConstants::MeV;
using Acts::UnitConstants::mm;
+ using Acts::UnitConstants::ns;
double p = std::hypot( part.psx() * GeV, part.psy() * GeV, part.psz() * GeV);
+ ROOT::Math::XYZVector momentum(part.psx() * GeV, part.psy() * GeV, part.psz() * GeV);
// build some track cov matrix
Acts::BoundSymMatrix cov = Acts::BoundSymMatrix::Zero();
@@ -83,15 +87,30 @@ namespace Jug::Reco {
cov(Acts::eBoundPhi, Acts::eBoundPhi) = M_PI / 180.0;
cov(Acts::eBoundTheta, Acts::eBoundTheta) = M_PI / 180.0;
cov(Acts::eBoundQOverP, Acts::eBoundQOverP) = 0.98 / (p*p);
- cov(Acts::eBoundTime, Acts::eBoundTime) = Acts::UnitConstants::ns;
-
- init_trk_params->emplace_back(Acts::Vector4D(part.vsx() * mm, part.vsy() * mm, part.vsz() * mm, part.time() * Acts::UnitConstants::ns),
- Acts::Vector3D(part.psx() * GeV, part.psy() * GeV, part.psz() * GeV),
- p+200*MeV,
- ((part.pdgID() > 0) ? -1 : 1),
- std::make_optional(std::move(cov))
- );
- //part .charge()
+ cov(Acts::eBoundTime, Acts::eBoundTime) = 1.0*ns;
+
+ Acts::BoundVector params;
+ params(Acts::eBoundLoc0) = 0.0 * mm ; // cylinder radius
+ params(Acts::eBoundLoc1) = 0.0 * mm ; // cylinder length
+ params(Acts::eBoundPhi) = momentum.Phi();
+ params(Acts::eBoundTheta) = momentum.Theta();
+ params(Acts::eBoundQOverP) = 1/p;
+ params(Acts::eBoundTime) = part.time() * ns;
+ /// \todo create or find better particle data interface.
+ // get the particle charge
+ int charge = ((part.pdgID() > 0) ? 1 : -1);
+ // electron is negative but positive pdg code
+ if( std::abs(part.pdgID()) == 11 ) {
+ charge *= -1;
+ }
+
+ //// Construct a perigee surface as the target surface
+ auto pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(
+ Acts::Vector3{part.vsx() * mm, part.vsy() * mm, part.vsz() * mm});
+
+ params(Acts::eBoundQOverP) = charge/p;
+ init_trk_params->push_back({pSurface, params, charge});
+ // std::make_optional(std::move(cov))
debug() << "Invoke track finding seeded by truth particle with p = " << p/GeV << " GeV" << endmsg;
//Acts::BoundMatrix cov = Acts::BoundMatrix::Zero();
diff --git a/JugReco/src/components/TrackParamVertexClusterInit.cpp b/JugReco/src/components/TrackParamVertexClusterInit.cpp
index 9b9fc91471f05dcf3f6eee3c815c625932cdb306..9855071d1594a963cb9f4f14099c2a6ff5429bff 100644
--- a/JugReco/src/components/TrackParamVertexClusterInit.cpp
+++ b/JugReco/src/components/TrackParamVertexClusterInit.cpp
@@ -12,11 +12,13 @@
#include "JugBase/DataHandle.h"
#include "JugBase/IGeoSvc.h"
#include "JugReco/Track.hpp"
-#include "Acts/Utilities/Units.hpp"
-#include "Acts/Utilities/Definitions.hpp"
+#include "Acts/Definitions/Units.hpp"
+#include "Acts/Definitions/Common.hpp"
#include "eicd/TrackerHitCollection.h"
#include "eicd/ClusterCollection.h"
+#include "Math/Vector3D.h"
+#include "Acts/Surfaces/PerigeeSurface.hpp"
using namespace Gaudi::Units;
@@ -79,6 +81,7 @@ namespace Jug::Reco {
continue;
}
+ auto pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3{0,0,0});
for (const auto& t : *vtx_hits) {
double len = std::hypot(t.x(), t.y(), t.z());
@@ -86,30 +89,54 @@ namespace Jug::Reco {
continue;
}
- // build some track cov matrix
- Acts::BoundSymMatrix cov = Acts::BoundSymMatrix::Zero();
- cov(Acts::eBoundLoc0, Acts::eBoundLoc0) = 1.0 * mm*1.0 * mm;
- cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = 1.0 * mm*1.0 * mm;
- cov(Acts::eBoundPhi, Acts::eBoundPhi) = M_PI / 180.0;
- cov(Acts::eBoundTheta, Acts::eBoundTheta) = M_PI / 180.0;
- cov(Acts::eBoundQOverP, Acts::eBoundQOverP) = 1.0 / (0.9*0.9*p_cluster*p_cluster);
- cov(Acts::eBoundTime, Acts::eBoundTime) = Acts::UnitConstants::ns;
-
- // add all charges to the track candidate...
- init_trk_params->emplace_back(Acts::Vector4D(0 * mm, 0 * mm, 0 * mm, 0),
- Acts::Vector3D(t.x() * p_cluster / len, t.y() * p_cluster / len, t.z() * p_cluster / len), p_cluster, -1,
- std::make_optional(cov));
- init_trk_params->emplace_back(Acts::Vector4D(0 * mm, 0 * mm, 0 * mm, 0),
- Acts::Vector3D(t.x() * p_cluster / len, t.y() * p_cluster / len, t.z() * p_cluster / len), p_cluster, 1,
- std::make_optional(cov));
+ //// build some track cov matrix
+ //Acts::BoundSymMatrix cov = Acts::BoundSymMatrix::Zero();
+ //cov(Acts::eBoundLoc0, Acts::eBoundLoc0) = 1.0 * mm*1.0 * mm;
+ //cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = 1.0 * mm*1.0 * mm;
+ //cov(Acts::eBoundPhi, Acts::eBoundPhi) = M_PI / 180.0;
+ //cov(Acts::eBoundTheta, Acts::eBoundTheta) = M_PI / 180.0;
+ //cov(Acts::eBoundQOverP, Acts::eBoundQOverP) = 1.0 / (0.9*0.9*p_cluster*p_cluster);
+ //cov(Acts::eBoundTime, Acts::eBoundTime) = Acts::UnitConstants::ns;
+
+ //// add all charges to the track candidate...
+ //init_trk_params->emplace_back(Acts::Vector4(0 * mm, 0 * mm, 0 * mm, 0),
+ // Acts::Vector3(t.x() * p_cluster / len, t.y() * p_cluster / len, t.z() * p_cluster / len), p_cluster, -1,
+ // std::make_optional(cov));
+ //init_trk_params->emplace_back(Acts::Vector4(0 * mm, 0 * mm, 0 * mm, 0),
+ // Acts::Vector3(t.x() * p_cluster / len, t.y() * p_cluster / len, t.z() * p_cluster / len), p_cluster, 1,
+ // std::make_optional(cov));
+
+ ROOT::Math::XYZVector momentum(t.x() * p_cluster / len, t.y() * p_cluster / len, t.z() * p_cluster / len);
+
+ Acts::BoundVector params;
+ params(Acts::eBoundLoc0) = 0.0 * mm;
+ params(Acts::eBoundLoc1) = 0.0 * mm;
+ params(Acts::eBoundPhi) = momentum.Phi();
+ params(Acts::eBoundTheta) = momentum.Theta();
+ params(Acts::eBoundQOverP) = 1 / p_cluster;
+ params(Acts::eBoundTime) = 0 * ns;
+
+ //debug() << "Invoke track finding seeded by truth particle with p = " << p / GeV << " GeV" << endmsg;
+
+ // add both charges to the track candidate...
+ init_trk_params->push_back({pSurface, params, 1});
+
+ Acts::BoundVector params2;
+ params2(Acts::eBoundLoc0) = 0.0 * mm;
+ params2(Acts::eBoundLoc1) = 0.0 * mm;
+ params2(Acts::eBoundPhi) = momentum.Phi();
+ params2(Acts::eBoundTheta) = momentum.Theta();
+ params2(Acts::eBoundQOverP) = -1 / p_cluster;
+ params2(Acts::eBoundTime) = 0 * ns;
+ init_trk_params->push_back({pSurface, params2, -1});
}
- //init_trk_params->emplace_back(Acts::Vector4D(0 * mm, 0 * mm, 0 * mm, 0),
- // Acts::Vector3D(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, 1,
+ // init_trk_params->emplace_back(Acts::Vector4(0 * mm, 0 * mm, 0 * mm, 0),
+ // Acts::Vector3(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, 1,
// std::make_optional(cov));
- //init_trk_params->emplace_back(Acts::Vector4D(0 * mm, 0 * mm, 0 * mm, 0),
- // Acts::Vector3D(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, 0,
+ // init_trk_params->emplace_back(Acts::Vector4(0 * mm, 0 * mm, 0 * mm, 0),
+ // Acts::Vector3(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, 0,
// std::make_optional(cov));
- debug() << "Invoke track finding seeded by truth particle with p = " << p_cluster/GeV << " GeV" << endmsg;
+ // debug() << "Invoke track finding seeded by truth particle with p = " << p_cluster/GeV << " GeV" << endmsg;
}
return StatusCode::SUCCESS;
}
diff --git a/JugReco/src/components/TrackerHitReconstruction.cpp b/JugReco/src/components/TrackerHitReconstruction.cpp
index 43e7190f117e72c7bac6b3ff5a21275fe5b38f51..d32c854fe8223e7d0efcf773af079b51cb35e9f5 100644
--- a/JugReco/src/components/TrackerHitReconstruction.cpp
+++ b/JugReco/src/components/TrackerHitReconstruction.cpp
@@ -2,11 +2,11 @@
// Gaudi
#include "GaudiAlg/GaudiAlgorithm.h"
-#include "GaudiKernel/ToolHandle.h"
+//#include "GaudiKernel/ToolHandle.h"
#include "GaudiAlg/Transformer.h"
#include "GaudiAlg/GaudiTool.h"
#include "GaudiKernel/RndmGenerators.h"
-#include "GaudiKernel/Property.h"
+#include "Gaudi/Property.h"
#include "DDRec/CellIDPositionConverter.h"
#include "DDRec/SurfaceManager.h"
@@ -20,7 +20,6 @@
//#include "GaudiExamples/MyTrack.h"
#include "eicd/RawTrackerHitCollection.h"
#include "eicd/TrackerHitCollection.h"
-#include "JugReco/SourceLinks.h"
#include "DD4hep/DD4hepUnits.h"
diff --git a/JugReco/src/components/TrackerSourceLinker.cpp b/JugReco/src/components/TrackerSourceLinker.cpp
index 0927d047782006423cde6308f84b25eca47c3a24..8e78cf5ac1a62ab4d60623f67ef68cc6a1c3537e 100644
--- a/JugReco/src/components/TrackerSourceLinker.cpp
+++ b/JugReco/src/components/TrackerSourceLinker.cpp
@@ -6,7 +6,7 @@
#include "GaudiAlg/Transformer.h"
#include "GaudiAlg/GaudiTool.h"
#include "GaudiKernel/RndmGenerators.h"
-#include "GaudiKernel/Property.h"
+#include "Gaudi/Property.h"
#include "JugBase/DataHandle.h"
#include "JugBase/IGeoSvc.h"
@@ -18,12 +18,15 @@
#include "DD4hep/DD4hepUnits.h"
-#include "Acts/Utilities/Units.hpp"
-#include "Acts/Utilities/Definitions.hpp"
+#include "Acts/Surfaces/Surface.hpp"
+#include "Acts/Definitions/Units.hpp"
+#include "Acts/Definitions/Common.hpp"
#include "Acts/Geometry/TrackingGeometry.hpp"
#include "Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp"
-#include "JugReco/SourceLinks.h"
+#include "JugReco/Index.hpp"
+#include "JugReco/IndexSourceLink.hpp"
+#include "JugReco/Measurement.hpp"
#include "eicd/TrackerHitCollection.h"
@@ -32,7 +35,8 @@ namespace Jug::Reco {
class TrackerSourceLinker : public GaudiAlgorithm {
public:
DataHandle<eic::TrackerHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
- DataHandle<SourceLinkContainer> m_outputSourceLinks{"outputSourceLinks", Gaudi::DataHandle::Writer, this};
+ DataHandle<IndexSourceLinkContainer> m_outputSourceLinks{"outputSourceLinks", Gaudi::DataHandle::Writer, this};
+ DataHandle<MeasurementContainer> m_outputMeasurements{"outputMeasurements", Gaudi::DataHandle::Writer, this};
/// Pointer to the geometry service
SmartIF<IGeoSvc> m_geoSvc;
@@ -44,6 +48,7 @@ namespace Jug::Reco {
: GaudiAlgorithm(name, svcLoc) {
declareProperty("inputHitCollection", m_inputHitCollection, "");
declareProperty("outputSourceLinks", m_outputSourceLinks, "");
+ declareProperty("outputMeasurements", m_outputMeasurements, "");
}
StatusCode initialize() override {
@@ -77,60 +82,60 @@ namespace Jug::Reco {
// input collection
const eic::TrackerHitCollection* hits = m_inputHitCollection.get();
// Create output collections
- auto source_links = m_outputSourceLinks.createAndPut();
- // setup local covariance
- // TODO add support for per volume/layer/module settings
+ auto sourceLinks = m_outputSourceLinks.createAndPut();
+ auto measurements = m_outputMeasurements.createAndPut();
+ // IndexMultimap<ActsFatras::Barcode> hitParticlesMap;
+ // IndexMultimap<Index> hitSimHitsMap;
+ sourceLinks->reserve(hits->size());
+ measurements->reserve(hits->size());
debug() << (*hits).size() << " hits " << endmsg;
+ int ihit = 0;
for(const auto& ahit : *hits) {
- Acts::BoundMatrix cov = Acts::BoundMatrix::Zero();
- cov(Acts::eBoundLoc0, Acts::eBoundLoc0) = ahit.covsym_xx()*Acts::UnitConstants::mm*ahit.covsym_xx()*Acts::UnitConstants::mm;
- cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = ahit.covsym_yy()*Acts::UnitConstants::mm*ahit.covsym_yy()*Acts::UnitConstants::mm;
+ Acts::SymMatrix2 cov = Acts::SymMatrix2::Zero();
+ cov(0,0) = ahit.covsym_xx()*Acts::UnitConstants::mm*ahit.covsym_xx()*Acts::UnitConstants::mm;
+ cov(1,1) = ahit.covsym_yy()*Acts::UnitConstants::mm*ahit.covsym_yy()*Acts::UnitConstants::mm;
- auto vol_ctx = m_geoSvc->cellIDPositionConverter()->findContext(ahit.cellID());
- auto vol_id = vol_ctx->identifier;
- //debug() << " hit : \n" << ahit << endmsg;
- //debug() << "cell ID : " << ahit.cellID() << endmsg;
- //debug() << " position : (" << ahit.position(0) << ", " << ahit.position(1) << ", "<< ahit.position(2) << ") " << endmsg;
- debug() << " vol_id : " << vol_id << endmsg;
- debug() << " placment pos : " << vol_ctx->volumePlacement().position() << endmsg;
-
- const auto is = m_surfaces.find(vol_id);
+ auto vol_ctx = m_geoSvc->cellIDPositionConverter()->findContext(ahit.cellID());
+ auto vol_id = vol_ctx->identifier;
+ const auto is = m_surfaces.find(vol_id);
if (is == m_surfaces.end()) {
debug() << " vol_id (" << vol_id << ") not found in m_surfaces." <<endmsg;
continue;
}
const Acts::Surface* surface = is->second;
-
debug() << " surface center : " << surface->center(Acts::GeometryContext()) << endmsg;
// transform global position into local coordinates
- Acts::Vector2D pos(0, 0);
+ Acts::Vector2 pos(0, 0);
// geometry context contains nothing here
pos = surface->globalToLocal(Acts::GeometryContext(), {ahit.x(), ahit.y(), ahit.z()}, {0, 0, 0}).value();//, pos);
- //// smear truth to create local measurement
- Acts::BoundVector loc = Acts::BoundVector::Zero();
+ Acts::Vector2 loc = Acts::Vector2::Zero();
loc[Acts::eBoundLoc0] = pos[0] ;//+ m_cfg.sigmaLoc0 * stdNormal(rng);
loc[Acts::eBoundLoc1] = pos[1] ;//+ m_cfg.sigmaLoc1 * stdNormal(rng);
+ //debug() << "loc : (" << loc[0] << ", " << loc[1] << ")" << endmsg;
- debug() << "loc : (" << loc[0] << ", " << loc[1] << ")" << endmsg;
+ //local position
+ //auto loc = {ahit.x(), ahit.y(), ahit.z()} - vol_ctx->volumePlacement().position()
+ //debug() << " hit : \n" << ahit << endmsg;
+ //debug() << " cell ID : " << ahit.cellID() << endmsg;
+ //debug() << " position : (" << ahit.position(0) << ", " << ahit.position(1) << ", "<< ahit.position(2) << ") " << endmsg;
+ //debug() << " vol_id : " << vol_id << endmsg;
+ //debug() << " placment pos : " << vol_ctx->volumePlacement().position() << endmsg;
- // create source link at the end of the container
- //auto it = source_links->emplace_hint(source_links->end(), *surface, hit, 2, loc, cov);
- auto it = source_links->emplace_hint(source_links->end(), *surface, 2, loc, cov);
- // ensure hits and links share the same order to prevent ugly surprises
- if (std::next(it) != source_links->end()) {
- error() << "The hit ordering broke. Run for your life." << endmsg;
- return StatusCode::FAILURE;
- }
+ // the measurement container is unordered and the index under which the
+ // measurement will be stored is known before adding it.
+ Index hitIdx = measurements->size();
+ IndexSourceLink sourceLink(vol_id, ihit);
+ auto meas = Acts::makeMeasurement(sourceLink, loc, cov, Acts::eBoundLoc0, Acts::eBoundLoc1);
+
+ // add to output containers. since the input is already geometry-order,
+ // new elements in geometry containers can just be appended at the end.
+ sourceLinks->emplace_hint(sourceLinks->end(), std::move(sourceLink));
+ measurements->emplace_back(std::move(meas));
- //std::array<double,3> posarr; pos.GetCoordinates(posarr);
- //std::array<double,3> dimarr; dim.GetCoordinates(posarr);
- //eic::TrackerHit hit;
- //eic::TrackerHit hit((long long)ahit.cellID0(), (long long)ahit.cellID(), (long long)ahit.time(),
- // (float)ahit.charge() / 10000.0, (float)0.0, {{pos.x(), pos.y(),pos.z()}},{{dim[0],dim[1],0.0}});
- //rec_hits->push_back(hit);
+ ihit++;
}
return StatusCode::SUCCESS;
}
@@ -140,82 +145,3 @@ namespace Jug::Reco {
} // namespace Jug::reco
-//HitSmearing::HitSmearing(const Config& cfg, Acts::Logging::Level lvl)
-// : BareAlgorithm("HitSmearing", lvl), m_cfg(cfg) {
-// if (m_cfg.inputSimulatedHits.empty()) {
-// throw std::invalid_argument("Missing input simulated hits collection");
-// }
-// if (m_cfg.outputSourceLinks.empty()) {
-// throw std::invalid_argument("Missing output source links collection");
-// }
-// if ((m_cfg.sigmaLoc0 < 0) or (m_cfg.sigmaLoc1 < 0)) {
-// throw std::invalid_argument("Invalid resolution setting");
-// }
-// if (not m_cfg.trackingGeometry) {
-// throw std::invalid_argument("Missing tracking geometry");
-// }
-// if (!m_cfg.randomNumbers) {
-// throw std::invalid_argument("Missing random numbers tool");
-// }
-// // fill the surface map to allow lookup by geometry id only
-// m_cfg.trackingGeometry->visitSurfaces([this](const Acts::Surface* surface) {
-// // for now we just require a valid surface
-// if (not surface) {
-// return;
-// }
-// this->m_surfaces.insert_or_assign(surface->associatedDetectorElement()->identifier(), surface);
-// });
-// }
-//}
-
-//ProcessCode HitSmearing::execute(const AlgorithmContext& ctx) const {
-// // setup input and output containers
-// const auto& hits =
-// ctx.eventStore.get<SimHitContainer>(m_cfg.inputSimulatedHits);
-// SimSourceLinkContainer sourceLinks;
-// sourceLinks.reserve(hits.size());
-//
-// // setup random number generator
-// auto rng = m_cfg.randomNumbers->spawnGenerator(ctx);
-// std::normal_distribution<double> stdNormal(0.0, 1.0);
-//
-// // setup local covariance
-// // TODO add support for per volume/layer/module settings
-// Acts::BoundMatrix cov = Acts::BoundMatrix::Zero();
-// cov(Acts::eLOC_0, Acts::eLOC_0) = m_cfg.sigmaLoc0 * m_cfg.sigmaLoc0;
-// cov(Acts::eLOC_1, Acts::eLOC_1) = m_cfg.sigmaLoc1 * m_cfg.sigmaLoc1;
-//
-// for (auto&& [moduleGeoId, moduleHits] : groupByModule(hits)) {
-// // check if we should create hits for this surface
-// const auto is = m_surfaces.find(moduleGeoId);
-// if (is == m_surfaces.end()) {
-// continue;
-// }
-//
-// // smear all truth hits for this module
-// const Acts::Surface* surface = is->second;
-// for (const auto& hit : moduleHits) {
-// // transform global position into local coordinates
-// Acts::Vector2D pos(0, 0);
-// surface->globalToLocal(ctx.geoContext, hit.position(),
-// hit.unitDirection(), pos);
-//
-// // smear truth to create local measurement
-// Acts::BoundVector loc = Acts::BoundVector::Zero();
-// loc[Acts::eLOC_0] = pos[0] + m_cfg.sigmaLoc0 * stdNormal(rng);
-// loc[Acts::eLOC_1] = pos[1] + m_cfg.sigmaLoc1 * stdNormal(rng);
-//
-// // create source link at the end of the container
-// auto it = sourceLinks.emplace_hint(sourceLinks.end(), *surface, hit, 2,
-// loc, cov);
-// // ensure hits and links share the same order to prevent ugly surprises
-// if (std::next(it) != sourceLinks.end()) {
-// ACTS_FATAL("The hit ordering broke. Run for your life.");
-// return ProcessCode::ABORT;
-// }
-// }
-// }
-//
-// ctx.eventStore.add(m_cfg.outputSourceLinks, std::move(sourceLinks));
-// return ProcessCode::SUCCESS;
-//}
diff --git a/JugReco/src/components/TrackerSourcesLinker.cpp b/JugReco/src/components/TrackerSourcesLinker.cpp
index d054487a2c25f00a336034bfc6d02375d8fe6353..c1ab51f2eb661f3335e70c66394222e4cbcfe08a 100644
--- a/JugReco/src/components/TrackerSourcesLinker.cpp
+++ b/JugReco/src/components/TrackerSourcesLinker.cpp
@@ -1,4 +1,5 @@
#include "JugReco/GeometryContainers.hpp"
+#include <ios>
// Gaudi
#include "GaudiAlg/GaudiAlgorithm.h"
@@ -15,15 +16,18 @@
#include "DDRec/SurfaceManager.h"
#include "DDRec/Surface.h"
#include "DD4hep/Volumes.h"
+#include "DD4hep/VolumeManager.h"
#include "DD4hep/DD4hepUnits.h"
-
-#include "Acts/Utilities/Units.hpp"
-#include "Acts/Utilities/Definitions.hpp"
+#include "Acts/Surfaces/Surface.hpp"
+#include "Acts/Definitions/Units.hpp"
+#include "Acts/Definitions/Common.hpp"
#include "Acts/Geometry/TrackingGeometry.hpp"
#include "Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp"
-#include "JugReco/SourceLinks.h"
+#include "JugReco/Index.hpp"
+#include "JugReco/IndexSourceLink.hpp"
+#include "JugReco/Measurement.hpp"
#include "eicd/TrackerHitCollection.h"
@@ -38,10 +42,12 @@ namespace Jug::Reco {
DataHandle<HitCol> m_ITrackerEndcapHits{"ITrackerEndcapHits", Gaudi::DataHandle::Reader, this};
DataHandle<HitCol> m_OTrackerBarrelHits{"OTrackerBarrelHits", Gaudi::DataHandle::Reader, this};
DataHandle<HitCol> m_OTrackerEndcapHits{"OTrackerEndcapHits", Gaudi::DataHandle::Reader, this};
+ DataHandle<HitCol> m_allTrackerHits{"allTrackerHits", Gaudi::DataHandle::Writer, this};
std::vector<DataHandle<HitCol>*> m_trackerHitCollections;
//Gaudi::Property<std::vector<std::string>> m_trackerHitCollections{this, "trackerHitCollections"};
- DataHandle<SourceLinkContainer> m_outputSourceLinks{"outputSourceLinks", Gaudi::DataHandle::Writer, this};
+ DataHandle<IndexSourceLinkContainer> m_outputSourceLinks{"outputSourceLinks", Gaudi::DataHandle::Writer, this};
+ DataHandle<MeasurementContainer> m_outputMeasurements{"outputMeasurements", Gaudi::DataHandle::Writer, this};
/// Pointer to the geometry service
SmartIF<IGeoSvc> m_geoSvc;
@@ -56,6 +62,8 @@ namespace Jug::Reco {
declareProperty("OTrackerBarrelHits", m_OTrackerBarrelHits, "");
declareProperty("OTrackerEndcapHits", m_OTrackerEndcapHits, "");
declareProperty("outputSourceLinks", m_outputSourceLinks, "");
+ declareProperty("outputMeasurements", m_outputMeasurements, "");
+ declareProperty("allTrackerHits", m_allTrackerHits, "");
}
@@ -64,7 +72,7 @@ namespace Jug::Reco {
return StatusCode::FAILURE;
m_trackerHitCollections.push_back(&m_ITrackerBarrelHits);
- m_trackerHitCollections.push_back(&m_ITrackerEndcapHits);
+ //m_trackerHitCollections.push_back(&m_ITrackerEndcapHits);
m_trackerHitCollections.push_back(&m_OTrackerBarrelHits);
m_trackerHitCollections.push_back(&m_OTrackerEndcapHits);
@@ -86,6 +94,9 @@ namespace Jug::Reco {
debug() << "invalid det_element!!! " << endmsg;
return;
}
+ if(det_element->identifier() == 542675589) {
+ debug() << "test volume FOUND!!!!!!!!!!! " << endmsg;
+ }
this->m_surfaces.insert_or_assign(det_element->identifier(), surface);
});
@@ -94,67 +105,79 @@ namespace Jug::Reco {
StatusCode execute() override {
// Create output collections
- auto source_links = m_outputSourceLinks.createAndPut();
- // setup local covariance
- // TODO add support for per volume/layer/module settings
+ auto sourceLinks = m_outputSourceLinks.createAndPut();
+ auto measurements = m_outputMeasurements.createAndPut();
+ auto allHits = m_allTrackerHits.createAndPut();
+ // IndexMultimap<ActsFatras::Barcode> hitParticlesMap;
+ // IndexMultimap<Index> hitSimHitsMap;
+ //int total_hits = 0;
+ //for(const auto col : m_trackerHitCollections) {
+ // total_hits += col->get();
+ //}
+ //sourceLinks->reserve(hits->size());
+ //measurements->reserve(hits->size());
debug() << " m_trackerHitCollections size : " << m_trackerHitCollections.size() << endmsg;
+ int ihit = 0;
for(const auto col : m_trackerHitCollections) {
// input collection
const eic::TrackerHitCollection* hits = col->get();
- //const eic::TrackerHitCollection* hits = get<eic::TrackerHitCollection>("/Event/"+ col);
+ // const eic::TrackerHitCollection* hits = get<eic::TrackerHitCollection>("/Event/"+ col);
debug() << (*hits).size() << " hits " << endmsg;
- for(const auto& ahit : *hits) {
+ for (const auto& ahit : *hits) {
+
+ //allHits->push_back(ahit.clone());
- Acts::BoundMatrix cov = Acts::BoundMatrix::Zero();
- cov(Acts::eBoundLoc0, Acts::eBoundLoc0) = ahit.covsym_xx()*Acts::UnitConstants::mm*ahit.covsym_xx()*Acts::UnitConstants::mm;
- cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = ahit.covsym_yy()*Acts::UnitConstants::mm*ahit.covsym_yy()*Acts::UnitConstants::mm;
+ Acts::SymMatrix2 cov = Acts::SymMatrix2::Zero();
+ cov(0, 0) = ahit.covsym_xx() * Acts::UnitConstants::mm * ahit.covsym_xx() * Acts::UnitConstants::mm;
+ cov(1, 1) = ahit.covsym_yy() * Acts::UnitConstants::mm * ahit.covsym_yy() * Acts::UnitConstants::mm;
- auto vol_ctx = m_geoSvc->cellIDPositionConverter()->findContext(ahit.cellID());
- auto vol_id = vol_ctx->identifier;
- //debug() << " hit : \n" << ahit << endmsg;
- //debug() << "cell ID : " << ahit.cellID() << endmsg;
- //debug() << " position : (" << ahit.position(0) << ", " << ahit.position(1) << ", "<< ahit.position(2) << ") " << endmsg;
- debug() << " vol_id : " << vol_id << endmsg;
- debug() << " placment pos : " << vol_ctx->volumePlacement().position() << endmsg;
+ auto vol_ctx = m_geoSvc->cellIDPositionConverter()->findContext(ahit.cellID());
+ auto vol_id = vol_ctx->identifier;
+ const auto is = m_surfaces.find(vol_id);
- const auto is = m_surfaces.find(vol_id);
+ debug() << " cell ID : 0x" << std::hex << ahit.cellID() << std::dec << endmsg;
+ debug() << " vol ID : 0x" << std::hex << vol_id << std::dec << endmsg;
if (is == m_surfaces.end()) {
- debug() << " vol_id (" << vol_id << ") not found in m_surfaces." <<endmsg;
+ auto detector = m_geoSvc->detector()->volumeManager().lookupDetector(vol_id);
+ debug() << " detector 0x" << std::hex << detector.id() <<std::dec << ", " << detector.path() << endmsg;
+ debug() << " vol_id " << vol_id << " (0x" << std::hex << vol_id <<std::dec<< ") not found in m_surfaces." << endmsg;
+ debug() << "hit " << ihit << endmsg;
continue;
}
const Acts::Surface* surface = is->second;
-
- debug() << " surface center : " << surface->center(Acts::GeometryContext()) << endmsg;
+ //debug() << " surface center : " << surface->center(Acts::GeometryContext()) << endmsg;
// transform global position into local coordinates
- Acts::Vector2D pos(0, 0);
+ Acts::Vector2 pos(0, 0);
// geometry context contains nothing here
- pos = surface->globalToLocal(Acts::GeometryContext(), {ahit.x(), ahit.y(), ahit.z()}, {0, 0, 0}).value();//, pos);
-
- //// smear truth to create local measurement
- Acts::BoundVector loc = Acts::BoundVector::Zero();
- loc[Acts::eBoundLoc0] = pos[0] ;//+ m_cfg.sigmaLoc0 * stdNormal(rng);
- loc[Acts::eBoundLoc0] = pos[1] ;//+ m_cfg.sigmaLoc1 * stdNormal(rng);
-
- debug() << "loc : (" << loc[0] << ", " << loc[1] << ")" << endmsg;
-
- // create source link at the end of the container
- //auto it = source_links->emplace_hint(source_links->end(), *surface, hit, 2, loc, cov);
- auto it = source_links->emplace_hint(source_links->end(), *surface, 2, loc, cov);
- // ensure hits and links share the same order to prevent ugly surprises
- if (std::next(it) != source_links->end()) {
- error() << "The hit ordering broke. Run for your life." << endmsg;
- return StatusCode::FAILURE;
- }
-
- //std::array<double,3> posarr; pos.GetCoordinates(posarr);
- //std::array<double,3> dimarr; dim.GetCoordinates(posarr);
- //eic::TrackerHit hit;
- //eic::TrackerHit hit((long long)ahit.cellID0(), (long long)ahit.cellID(), (long long)ahit.time(),
- // (float)ahit.charge() / 10000.0, (float)0.0, {{pos.x(), pos.y(),pos.z()}},{{dim[0],dim[1],0.0}});
- //rec_hits->push_back(hit);
+ pos = surface->globalToLocal(Acts::GeometryContext(), {ahit.x(), ahit.y(), ahit.z()}, {0, 0, 0})
+ .value(); //, pos);
+
+ Acts::Vector2 loc = Acts::Vector2::Zero();
+ loc[Acts::eBoundLoc0] = pos[0]; //+ m_cfg.sigmaLoc0 * stdNormal(rng);
+ loc[Acts::eBoundLoc1] = pos[1]; //+ m_cfg.sigmaLoc1 * stdNormal(rng);
+ // debug() << "loc : (" << loc[0] << ", " << loc[1] << ")" << endmsg;
+ // local position
+ // auto loc = {ahit.x(), ahit.y(), ahit.z()} - vol_ctx->volumePlacement().position()
+ //debug() << " hit : \n" << ahit << endmsg;
+ //debug() << " position : (" << ahit.x() << ", " << ahit.y() << ", "<< ahit.z() <<
+ //") " << endmsg; debug() << " vol_id : " << vol_id << endmsg; debug() << " placment pos : " <<
+ //vol_ctx->volumePlacement().position() << endmsg;
+
+ // the measurement container is unordered and the index under which the
+ // measurement will be stored is known before adding it.
+ Index hitIdx = measurements->size();
+ IndexSourceLink sourceLink(vol_id, ihit);
+ auto meas = Acts::makeMeasurement(sourceLink, loc, cov, Acts::eBoundLoc0, Acts::eBoundLoc1);
+
+ // add to output containers. since the input is already geometry-order,
+ // new elements in geometry containers can just be appended at the end.
+ sourceLinks->emplace_hint(sourceLinks->end(), std::move(sourceLink));
+ measurements->emplace_back(std::move(meas));
+
+ ihit++;
}
}
return StatusCode::SUCCESS;
diff --git a/JugReco/src/components/TrackingHitsSourceLinker.cpp b/JugReco/src/components/TrackingHitsSourceLinker.cpp
index e9c8d6b9c45b903245491c90a3f62ab241f0f2fc..ad7b147fc4153a0de1a3f1436ae6f614938f4d3a 100644
--- a/JugReco/src/components/TrackingHitsSourceLinker.cpp
+++ b/JugReco/src/components/TrackingHitsSourceLinker.cpp
@@ -20,8 +20,8 @@
#include "DD4hep/DD4hepUnits.h"
-#include "Acts/Utilities/Units.hpp"
-#include "Acts/Utilities/Definitions.hpp"
+#include "Acts/Definitions/Units.hpp"
+#include "Acts/Definitions/Common.hpp"
#include "Acts/Geometry/TrackingGeometry.hpp"
#include "Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp"
diff --git a/VERSION b/VERSION
index 26aaba0e86632e4d537006e45b0ec918d780b3b4..88c5fb891dcf1d1647d2b84bac0630cf9570d213 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-1.2.0
+1.4.0