Commit 4ff68c49 authored by Whitney Armstrong's avatar Whitney Armstrong
Browse files

Track finder is working

- Following the examples, the track finder is working with the
 hytop_test repo. There is still a lot of work to do but it is progress.
- Added source linker algorithm which uses just the reconstructed hits.
- Added track parameter initialization based on mc truth. This should be
replaced with proper seeding.
- Added TrackFindingAlgorithm following acts examples. This code is in a
single file and takes a long time to compile. This can be reduced with
some fwd declarations and headers.
parent da40d7d8
......@@ -28,10 +28,12 @@ target_link_libraries(JugBase
podio::podioRootIO
ActsCore ActsDD4hepPlugin
)
target_compile_options(JugBase PRIVATE -Wno-suggest-override)
gaudi_add_module(JugBasePlugins
src/components/*.cpp src/components/ReadTestConsumer.cxx
LINK_LIBRARIES GaudiAlgLib GaudiKernel JugBase ROOT NPDet::DD4podIO ActsCore ActsDD4hepPlugin)
target_compile_options(JugBasePlugins PRIVATE -Wno-suggest-override)
#gaudi_add_test(ProduceForReadTest
# WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
......
......@@ -29,6 +29,7 @@ gaudi_install_python_modules()
gaudi_add_module(JugDigiPlugins
src/components/*.cpp
LINK_LIBRARIES GaudiAlgLib GaudiKernel JugBase ROOT NPDet::DD4podIO EICD::eicd)
target_compile_options(JugDigiPlugins PRIVATE -Wno-suggest-override)
#gaudi_add_test(ProduceForReadTest
# WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
......
......@@ -30,7 +30,9 @@ gaudi_install_python_modules()
gaudi_add_module(JugRecoPlugins
src/components/TrackerHitReconstruction.cpp
src/components/TrackerSourceLinker.cpp
src/components/TrackFindingAlgorithm.cpp
src/components/TestACTSLogger.cpp
src/components/TrackParamTruthInit.cpp
LINK_LIBRARIES GaudiAlgLib GaudiKernel JugBase ROOT NPDet::DD4podIO EICD::eicd DDRec Acts)
target_compile_options(JugRecoPlugins PRIVATE -Wno-suggest-override)
......
#ifndef Jug_BFIELD_HH
#define Jug_BFIELD_HH 1
#include "Acts/Utilities//Definitions.hpp"
#include "Acts/Utilities/detail/AxisFwd.hpp"
#include "Acts/Utilities/detail/GridFwd.hpp"
#include <memory>
#include <tuple>
#include<variant>
#include "Acts/MagneticField/MagneticFieldContext.hpp"
#include "Acts/Utilities/Definitions.hpp"
// 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 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
} // namespace Jug
using InterpolatedMapper2D = Acts::InterpolatedBFieldMapper<
Acts::detail::Grid<Acts::Vector2D, Acts::detail::EquidistantAxis, Acts::detail::EquidistantAxis>>;
using InterpolatedMapper3D =
Acts::InterpolatedBFieldMapper<Acts::detail::Grid<Acts::Vector3D, Acts::detail::EquidistantAxis,
Acts::detail::EquidistantAxis, Acts::detail::EquidistantAxis>>;
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>>;
#endif
// 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 "JugReco/Utilities/GroupBy.hpp"
#include "JugReco/Utilities/Range.hpp"
#include "Acts/Geometry/GeometryID.hpp"
#include <algorithm>
#include <cstddef>
#include <utility>
#include <boost/container/flat_map.hpp>
#include <boost/container/flat_set.hpp>
namespace Jug {
namespace detail {
// extract the geometry identifier from a variety of types
struct GeometryIdGetter {
// explicit geometry identifier are just forwarded
constexpr Acts::GeometryID operator()(Acts::GeometryID geometryId) const {
return geometryId;
}
// encoded geometry ids are converted back to geometry identifiers.
constexpr Acts::GeometryID operator()(Acts::GeometryID::Value encoded) const {
return Acts::GeometryID(encoded);
}
// support elements in map-like structures.
template <typename T>
constexpr Acts::GeometryID operator()(
const std::pair<Acts::GeometryID, T>& mapItem) const {
return mapItem.first;
}
// support elements that implement `.geometryId()`.
template <typename T>
inline auto operator()(const T& thing) const
-> decltype(thing.geometryId(), Acts::GeometryID()) {
return thing.geometryId();
}
};
struct CompareGeometryId {
// indicate that comparisons between keys and full objects are allowed.
using is_transparent = void;
// compare two elements using the automatic key extraction.
template <typename Left, typename Right>
constexpr bool operator()(Left&& lhs, Right&& rhs) const {
return GeometryIdGetter()(lhs) < GeometryIdGetter()(rhs);
}
};
} // namespace detail
/// Store elements that know their detector geometry id, e.g. simulation hits.
///
/// @tparam T type to be stored, must be compatible with `CompareGeometryId`
///
/// The container stores an arbitrary number of elements for any geometry
/// id. Elements can be retrieved via the geometry id; elements can be selected
/// for a specific geometry id or for a larger range, e.g. a volume or a layer
/// within the geometry hierachy using the helper functions below. Elements can
/// also be accessed by index that uniquely identifies each element regardless
/// of geometry id.
template <typename T>
using GeometryIdMultiset =
boost::container::flat_multiset<T, detail::CompareGeometryId>;
/// Store elements indexed by an geometry id.
///
/// @tparam T type to be stored
///
/// The behaviour is the same as for the `GeometryIdMultiset` except that the
/// stored elements do not know their geometry id themself. When iterating
/// the iterator elements behave as for the `std::map`, i.e.
///
/// for (const auto& entry: elements) {
/// auto id = entry.first; // geometry id
/// const auto& el = entry.second; // stored element
/// }
///
template <typename T>
using GeometryIdMultimap = GeometryIdMultiset<std::pair<Acts::GeometryID, T>>;
/// Select all elements within the given volume.
template <typename T>
inline Range<typename GeometryIdMultiset<T>::const_iterator> selectVolume(
const GeometryIdMultiset<T>& container, Acts::GeometryID::Value volume) {
auto cmp = Acts::GeometryID().setVolume(volume);
auto beg = std::lower_bound(container.begin(), container.end(), cmp,
detail::CompareGeometryId{});
// WARNING overflows to volume==0 if the input volume is the last one
cmp = Acts::GeometryID().setVolume(volume + 1u);
// optimize search by using the lower bound as start point. also handles
// volume overflows since the geo id would be located before the start of
// the upper edge search window.
auto end =
std::lower_bound(beg, container.end(), cmp, detail::CompareGeometryId{});
return makeRange(beg, end);
}
template <typename T>
inline auto selectVolume(const GeometryIdMultiset<T>& container,
Acts::GeometryID id) {
return selectVolume(container, id.volume());
}
/// Select all elements within the given layer.
template <typename T>
inline Range<typename GeometryIdMultiset<T>::const_iterator> selectLayer(
const GeometryIdMultiset<T>& container, Acts::GeometryID::Value volume,
Acts::GeometryID::Value layer) {
auto cmp = Acts::GeometryID().setVolume(volume).setLayer(layer);
auto beg = std::lower_bound(container.begin(), container.end(), cmp,
detail::CompareGeometryId{});
// WARNING resets to layer==0 if the input layer is the last one
cmp = Acts::GeometryID().setVolume(volume).setLayer(layer + 1u);
// optimize search by using the lower bound as start point. also handles
// volume overflows since the geo id would be located before the start of
// the upper edge search window.
auto end =
std::lower_bound(beg, container.end(), cmp, detail::CompareGeometryId{});
return makeRange(beg, end);
}
template <typename T>
inline auto selectLayer(const GeometryIdMultiset<T>& container,
Acts::GeometryID id) {
return selectLayer(container, id.volume(), id.layer());
}
/// Select all elements for the given module / sensitive surface.
template <typename T>
inline Range<typename GeometryIdMultiset<T>::const_iterator> selectModule(
const GeometryIdMultiset<T>& container, Acts::GeometryID geoId) {
// module is the lowest level and defines a single geometry id value
return makeRange(container.equal_range(geoId));
}
template <typename T>
inline auto selectModule(const GeometryIdMultiset<T>& container,
Acts::GeometryID::Value volume,
Acts::GeometryID::Value layer,
Acts::GeometryID::Value module) {
return selectModule(
container,
Acts::GeometryID().setVolume(volume).setLayer(layer).setSensitive(
module));
}
/// Iterate over groups of elements belonging to each module/ sensitive surface.
template <typename T>
inline GroupBy<typename GeometryIdMultiset<T>::const_iterator,
detail::GeometryIdGetter>
groupByModule(const GeometryIdMultiset<T>& container) {
return makeGroupBy(container, detail::GeometryIdGetter());
}
} // namespace FW
// 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 <cstddef>
#include <utility>
#include <boost/container/flat_map.hpp>
namespace FW {
/// 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, typename Key = std::size_t>
using IndexMultimap = boost::container::flat_multimap<Key, Value>;
/// 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, typename Key>
inline IndexMultimap<Key, Value> invertIndexMultimap(
const IndexMultimap<Value, Key>& multimap) {
// switch key-value without enforcing the new ordering (linear copy)
typename IndexMultimap<Key, Value>::sequence_type unordered;
unordered.reserve(multimap.size());
for (const auto& keyValue : multimap) {
// value is now the key and the key is now the value
unordered.emplace_back(keyValue.second, keyValue.first);
}
// adopting the unordered sequence will reestablish the correct order
IndexMultimap<Key, Value> inverse;
inverse.adopt_sequence(std::move(unordered));
return inverse;
}
} // namespace FW
// 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 <cstddef>
#include <vector>
namespace FW {
/// A proto track is a collection of hits identified by their indices.
using ProtoTrack = std::vector<size_t>;
/// Container of proto tracks. Each proto track is identified by its index.
using ProtoTrackContainer = std::vector<ProtoTrack>;
} // namespace FW
// 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 "ACTFW/EventData/GeometryContainers.hpp"
//#include "ActsFatras/EventData/Hit.hpp"
//
//namespace FW {
//
//using SimHit = ::ActsFatras::Hit;
///// Store hits ordered by geometry identifier.
//using SimHitContainer = GeometryIdMultiset<::ActsFatras::Hit>;
//
//} // end of namespace FW
// 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/EventData/MeasurementHelpers.hpp"
//
//#include <cstddef>
//#include <vector>
//
//namespace FW {
//
///// A hit identifier with additional truth information.
/////
///// In an addition to a unique identifier, an list of additional indices is
///// stored. These can refer e.g. to particle indices or truth hit indices.
///// Using indices instead of pointers allows more flexibility, i.e. we are not
///// fixed to a specific object type when using e.g. pointers, and is more
///// robust since no requirements on stable memory locations of the pointed-to
///// objects are necessary (as would be the case for pointers).
//class SimIdentifier : public Acts::MinimalSourceLink {
// public:
// using Value = uint64_t;
// using Difference = int64_t;
//
// /// Constructor from encoded identifier value.
// ///
// /// @param value is the identifier value
// explicit SimIdentifier(Value value) : m_value(value) {}
// /// Constructor from encoded identifier value and truth information.
// ///
// /// @param value is the identifier value
// /// @param indices
// SimIdentifier(Value value, std::vector<std::size_t> indices)
// : m_value(value), m_indices(std::move(indices)) {}
//
// // Explicitely defaulted constructors and assignment operators
// SimIdentifier() = default;
// SimIdentifier(const SimIdentifier&) = default;
// SimIdentifier(SimIdentifier&&) = default;
// SimIdentifier& operator=(const SimIdentifier&) = default;
// SimIdentifier& operator=(SimIdentifier&&) = default;
//
// /// Assign from an identifier value.
// SimIdentifier& operator=(Value value);
// /// Cast to an identifier value.
// operator Value() const { return m_value; }
//
// /// Explicit access the underlying identifier value.
// Value value() const { return m_value; }
// /// Access all associated truth indices.
// const std::vector<std::size_t>& indices() const { return m_indices; }
//
// /// Attach a truth index to the identifier.
// void addIndex(std::size_t index) { m_indices.push_back(index); }
//
// private:
// /// The stored identifier value.
// Value m_value = 0u;
// /// Associated truth indices.
// std::vector<std::size_t> m_indices;
//
// friend constexpr bool operator<(const SimIdentifier& lhs,
// const SimIdentifier& rhs) {
// return lhs.m_value < rhs.m_value;
// }
// friend bool operator==(const SimIdentifier& lhs, const SimIdentifier& rhs) {
// return lhs.m_value == rhs.m_value;
// }
//};
//
//} // end of namespace FW
//
//using identifier_type = ::FW::SimIdentifier::Value;
//using identifier_diff = ::FW::SimIdentifier::Difference;
//using Identifier = ::FW::SimIdentifier;
// 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 "JugReco/SourceLinks.h"
//#include "ACTFW/Validation/ProtoTrackClassification.hpp"
#include "Acts/EventData/MultiTrajectory.hpp"
#include "Acts/EventData/TrackParameters.hpp"
#include <unordered_map>
#include <utility>
namespace Jug {
/// Associate a particle to its hit count within a proto track.
struct ParticleHitCount {
uint64_t particleId;
std::size_t hitCount;
};
using IndexedParams = std::unordered_map<size_t, Acts::BoundParameters>;
/// @brief Struct for truth track fitting/finding result with
/// Acts::KalmanFitter/Acts::CombinatorialKalmanFilter
///
/// It contains a MultiTrajectory with a vector of entry indices for individual
/// trajectories, and a map of fitted parameters indexed by the entry index.
/// In case of track fitting, there is at most one trajectory in the
/// MultiTrajectory; In case of track finding, there could be multiple
/// trajectories in the MultiTrajectory.
///
/// @note The MultiTrajectory is thought to be empty if there is no entry index
struct SimMultiTrajectory {
public:
/// @brief Default constructor
///
SimMultiTrajectory() = default;
/// @brief Constructor from multiTrajectory and fitted track parameters
///
/// @param multiTraj The multiTrajectory
/// @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,
const std::vector<size_t>& tTips,
const IndexedParams& parameters)
: m_multiTrajectory(multiTraj),
m_trackTips(tTips),
m_trackParameters(parameters) {}
/// @brief Copy constructor
///
/// @param rhs The source SimMultiTrajectory
SimMultiTrajectory(const SimMultiTrajectory& rhs)
: m_multiTrajectory(rhs.m_multiTrajectory),
m_trackTips(rhs.m_trackTips),
m_trackParameters(rhs.m_trackParameters) {}
/// @brief Copy move constructor
///
/// @param rhs The source SimMultiTrajectory
SimMultiTrajectory(SimMultiTrajectory&& rhs)
: m_multiTrajectory(std::move(rhs.m_multiTrajectory)),
m_trackTips(std::move(rhs.m_trackTips)),
m_trackParameters(std::move(rhs.m_trackParameters)) {}
/// @brief Default destructor
///
~SimMultiTrajectory() = default;
/// @brief assignment operator
///
/// @param rhs The source SimMultiTrajectory
SimMultiTrajectory& operator=(const SimMultiTrajectory& rhs) {