diff --git a/Core/include/Acts/Geometry/GeometryID.hpp b/Core/include/Acts/Geometry/GeometryID.hpp
index 62ff52245b42e51dc2b0ff28c61bcaf23897585e..878293081c77988f782e94d4b49d727a40df426b 100644
--- a/Core/include/Acts/Geometry/GeometryID.hpp
+++ b/Core/include/Acts/Geometry/GeometryID.hpp
@@ -9,6 +9,7 @@
#pragma once
#include <cstdint>
+#include <functional>
#include <iosfwd>
namespace Acts {
@@ -113,3 +114,13 @@ class GeometryID {
std::ostream& operator<<(std::ostream& os, GeometryID id);
} // namespace Acts
+
+// specialize std::hash so GeometryId can be used e.g. in an unordered_map
+namespace std {
+template <>
+struct hash<Acts::GeometryID> {
+ auto operator()(Acts::GeometryID gid) const noexcept {
+ return std::hash<Acts::GeometryID::Value>()(gid.value());
+ }
+};
+} // namespace std
diff --git a/Core/include/Acts/Utilities/MultiIndex.hpp b/Core/include/Acts/Utilities/MultiIndex.hpp
index 6520123d9feb47b755db2764aa07339e07581d0e..fbcca4e51b955b049a93f204ac936cc96a18ff90 100644
--- a/Core/include/Acts/Utilities/MultiIndex.hpp
+++ b/Core/include/Acts/Utilities/MultiIndex.hpp
@@ -148,15 +148,13 @@ class MultiIndex {
} // namespace Acts
+// specialize std::hash so MultiIndex can be used e.g. in an unordered_map
namespace std {
-
-// specialize std::hash so the MultiIndex can be used e.g. in an unordered_map
template <typename Storage, std::size_t... BitsPerLevel>
struct hash<Acts::MultiIndex<Storage, BitsPerLevel...>> {
- constexpr auto operator()(
- Acts::MultiIndex<Storage, BitsPerLevel...> idx) const noexcept {
+ auto operator()(Acts::MultiIndex<Storage, BitsPerLevel...> idx) const
+ noexcept {
return std::hash<Storage>()(idx.value());
}
};
-
} // namespace std
diff --git a/Core/include/Acts/Utilities/PdgParticle.hpp b/Core/include/Acts/Utilities/PdgParticle.hpp
index 3ad4722ca26f91a1e38f8316be049eec5587058a..56a60c7f06a2511a2b1b76d45be57ec61b4676dd 100644
--- a/Core/include/Acts/Utilities/PdgParticle.hpp
+++ b/Core/include/Acts/Utilities/PdgParticle.hpp
@@ -32,4 +32,10 @@ enum PdgParticle : int32_t {
eAntiProton = -eProton,
};
+/// Convert an anti-particle to its particle and leave particles as-is.
+constexpr inline PdgParticle makeAbsolutePdgParticle(PdgParticle pdg) {
+ const auto value = static_cast<int32_t>(pdg);
+ return static_cast<PdgParticle>((0 <= value) ? value : -value);
+}
+
} // namespace Acts
diff --git a/Core/include/Acts/Utilities/CurvilinearUnitVectors.hpp b/Core/include/Acts/Utilities/UnitVectors.hpp
similarity index 66%
rename from Core/include/Acts/Utilities/CurvilinearUnitVectors.hpp
rename to Core/include/Acts/Utilities/UnitVectors.hpp
index fe5338226306b29ca75e9da3d348148c47b8ab4a..bc502bc18783853e8b90b677bb549757f663b609 100644
--- a/Core/include/Acts/Utilities/CurvilinearUnitVectors.hpp
+++ b/Core/include/Acts/Utilities/UnitVectors.hpp
@@ -8,6 +8,7 @@
#pragma once
+#include <cmath>
#include <limits>
#include <utility>
@@ -15,6 +16,46 @@
namespace Acts {
+/// Construct a normalized direction vector from phi angle and pseudorapidity.
+///
+/// @param phi is the direction angle in the x-y plane.
+/// @param eta is the pseudorapidity towards the z-axis.
+///
+/// @note The input arguments intentionally use the same template type so that
+/// a compile error occurs if inconsistent input types are used. Avoids
+/// unexpected implicit type conversions and forces the user to
+/// explicitely cast missmatched input types.
+template <typename T>
+inline ActsVector<T, 3> makeDirectionUnitFromPhiEta(T phi, T eta) {
+ const auto coshEta = std::cosh(eta);
+ const auto tanhEta = std::tanh(eta);
+ return {
+ std::cos(phi) / coshEta,
+ std::sin(phi) / coshEta,
+ std::tanh(eta),
+ };
+}
+
+/// Construct a normalized direction vector from phi and theta angle.
+///
+/// @param phi is the direction angle in radian in the x-y plane.
+/// @param theta is the polar angle in radian towards the z-axis.
+///
+/// @note The input arguments intentionally use the same template type so that
+/// a compile error occurs if inconsistent input types are used. Avoids
+/// unexpected implicit type conversions and forces the user to
+/// explicitely cast missmatched input types.
+template <typename T>
+inline ActsVector<T, 3> makeDirectionUnitFromPhiTheta(T phi, T theta) {
+ const auto cosTheta = std::cos(theta);
+ const auto sinTheta = std::sin(theta);
+ return {
+ std::cos(phi) * sinTheta,
+ std::sin(phi) * sinTheta,
+ cosTheta,
+ };
+}
+
/// Construct the first curvilinear unit vector `U` for the given direction.
///
/// @param direction is the input direction vector
diff --git a/Fatras/CMakeLists.txt b/Fatras/CMakeLists.txt
index 8939d92b10fb2af7859c1d52af8ae682086ba92b..1390603518bc373a8467524aa28526936f842411 100644
--- a/Fatras/CMakeLists.txt
+++ b/Fatras/CMakeLists.txt
@@ -3,7 +3,8 @@ add_library(
src/LandauDistribution.cpp
src/Particle.cpp
src/ParticleData.cpp
- src/ProcessType.cpp)
+ src/ProcessType.cpp
+ src/StandardPhysicsLists.cpp)
target_compile_features(
ActsFatras
PUBLIC cxx_std_17)
diff --git a/Fatras/include/ActsFatras/EventData/Barcode.hpp b/Fatras/include/ActsFatras/EventData/Barcode.hpp
index f19e11be3e52901fd499d545537e5c0b0790590a..5a523b72650df7fcfdc16d5dc67b92b5f7e314f6 100644
--- a/Fatras/include/ActsFatras/EventData/Barcode.hpp
+++ b/Fatras/include/ActsFatras/EventData/Barcode.hpp
@@ -40,7 +40,7 @@ namespace ActsFatras {
/// a consequence is that only non-zero generation can have non-zero
/// sub-particle numbers. A non-zero generation indicates that the particle
/// is a descendant of the original particle, e.g. from interactions or decay,
-/// while the sub-particle number identifies the the descendat particle.
+/// while the sub-particle number identifies the descendant particle.
///
/// With this encoding, non-primary particles and their primary parent can
/// be easily identified at the expense of not storing the exact decay history.
@@ -92,7 +92,12 @@ namespace ActsFatras {
/// generation to contain unique values. However, this can only be done when all
/// particles are known.
class Barcode : public Acts::MultiIndex<uint64_t, 12, 12, 16, 8, 16> {
+ using Base = Acts::MultiIndex<uint64_t, 12, 12, 16, 8, 16>;
+
public:
+ using Base::Base;
+ using Base::Value;
+
/// Return the primary vertex identifier.
constexpr Value vertexPrimary() const { return level(0); }
/// Return the secondary vertex identifier.
@@ -124,3 +129,13 @@ class Barcode : public Acts::MultiIndex<uint64_t, 12, 12, 16, 8, 16> {
};
} // namespace ActsFatras
+
+// specialize std::hash so Barcode can be used e.g. in an unordered_map
+namespace std {
+template <>
+struct hash<ActsFatras::Barcode> {
+ auto operator()(ActsFatras::Barcode barcode) const noexcept {
+ return std::hash<ActsFatras::Barcode::Value>()(barcode.value());
+ }
+};
+} // namespace std
diff --git a/Fatras/include/ActsFatras/EventData/Hit.hpp b/Fatras/include/ActsFatras/EventData/Hit.hpp
index a4cad6e2bd1b26a72b86f7eec01fc8838bad97fc..95591338e8ac9bd7cb9501bb900c5975659b78b0 100644
--- a/Fatras/include/ActsFatras/EventData/Hit.hpp
+++ b/Fatras/include/ActsFatras/EventData/Hit.hpp
@@ -44,12 +44,12 @@ class Hit {
/// users responsibility to ensure that the position correspond to a
/// position on the given surface. The four-vector component order must be
/// [x,y,z,t] and [px,py,pz,E].
- template <typename Position4, typename Momentum4>
- Hit(Acts::GeometryID geoId, Barcode particleId,
+ template <typename Position4, typename Momentum40, typename Momentum41>
+ Hit(Acts::GeometryID geometryId, Barcode particleId,
const Eigen::MatrixBase<Position4>& pos4,
- const Eigen::MatrixBase<Momentum4>& before4,
- const Eigen::MatrixBase<Momentum4>& after4, int32_t index_ = -1)
- : m_geoId(geoId),
+ const Eigen::MatrixBase<Momentum40>& before4,
+ const Eigen::MatrixBase<Momentum41>& after4, int32_t index_ = -1)
+ : m_geometryId(geometryId),
m_particleId(particleId),
m_index(index_),
m_pos4(pos4),
@@ -61,13 +61,13 @@ class Hit {
Hit& operator=(Hit&&) = default;
/// Geometry identifier of the hit surface.
- Acts::GeometryID geoId() const { return m_geoId; }
+ constexpr Acts::GeometryID geometryId() const { return m_geometryId; }
/// Particle identifier of the particle that generated the hit.
- Barcode particleId() const { return m_particleId; }
+ constexpr Barcode particleId() const { return m_particleId; }
/// Hit index along the particle trajectory.
///
/// @retval negative if the hit index is undefined.
- int32_t index() const { return m_index; }
+ constexpr int32_t index() const { return m_index; }
/// Space-time position four-vector.
///
@@ -86,12 +86,14 @@ class Hit {
///
/// The component order is [px,py,pz,E].
const Vector4& momentum4After() const { return m_after4; }
- /// Particle direction before the hit.
- Vector3 directionBefore() const { return m_before4.head<3>().normalized(); }
- /// Particle direction after the hit.
- Vector3 directionAfter() const { return m_after4.head<3>().normalized(); }
- /// Average particle direction while passing through the surface.
- Vector3 direction() const {
+ /// Normalized particle direction vector before the hit.
+ Vector3 unitDirectionBefore() const {
+ return m_before4.head<3>().normalized();
+ }
+ /// Normalized particle direction vector the hit.
+ Vector3 unitDirectionAfter() const { return m_after4.head<3>().normalized(); }
+ /// Average normalized particle direction vector through the surface.
+ Vector3 unitDirection() const {
auto dir0 = m_before4 / (2 * m_before4.head<3>().norm());
auto dir1 = m_after4 / (2 * m_after4.head<3>().norm());
return (dir0 + dir1).head<3>().normalized();
@@ -100,11 +102,11 @@ class Hit {
///
/// @retval positive if the particle lost energy when it passed the surface
/// @retval negative if magic was involved
- Scalar depositedEnergy() const { return m_before4.w() - m_after4.w(); }
+ Scalar depositedEnergy() const { return m_before4[3] - m_after4[3]; }
private:
/// Identifier of the surface.
- Acts::GeometryID m_geoId;
+ Acts::GeometryID m_geometryId;
/// Identifier of the generating particle.
Barcode m_particleId;
/// Index of the hit along the particle trajectory.
diff --git a/Fatras/include/ActsFatras/EventData/Particle.hpp b/Fatras/include/ActsFatras/EventData/Particle.hpp
index 498d0e96022a141e21a0973eaec464f1f08e7933..4225500a0535cf771f5450cff8873f2ca53690bd 100644
--- a/Fatras/include/ActsFatras/EventData/Particle.hpp
+++ b/Fatras/include/ActsFatras/EventData/Particle.hpp
@@ -9,6 +9,7 @@
#pragma once
#include <cmath>
+#include <iosfwd>
#include <limits>
#include "Acts/Utilities/Definitions.hpp"
@@ -29,27 +30,39 @@ class Particle {
Particle() = default;
/// Construct a particle at rest with explicit mass and charge.
///
- /// @param id Particle identifier within an event
- /// @param pdg PDG particle number
+ /// @param particleId Particle identifier within an event
+ /// @param pdg PDG id
/// @param charge Particle charge in native units
- /// @param mass Particle mass in native units
+ /// @param mass Particle mass in native units
///
/// @warning It is the users responsibility that charge and mass match
/// the PDG particle number.
- Particle(Barcode id, Acts::PdgParticle pdg, Scalar charge, Scalar mass)
- : m_id(id), m_pdg(pdg), m_charge(charge), m_mass(mass) {}
+ Particle(Barcode particleId, Acts::PdgParticle pdg, Scalar charge,
+ Scalar mass)
+ : m_particleId(particleId), m_pdg(pdg), m_charge(charge), m_mass(mass) {}
/// Construct a particle at rest from a PDG particle number.
///
- /// @param id Particle identifier within an event
+ /// @param particleId Particle identifier within an event
/// @param pdg PDG particle number
///
/// Charge and mass are retrieved from the particle data table.
- Particle(Barcode id, Acts::PdgParticle pdg);
+ Particle(Barcode particleId, Acts::PdgParticle pdg);
Particle(const Particle &) = default;
Particle(Particle &&) = default;
Particle &operator=(const Particle &) = default;
Particle &operator=(Particle &&) = default;
+ /// Construct a new particle with a new identifier but same kinematics.
+ ///
+ /// @note This is intentionally not a regular setter. The particle id
+ /// is used to identify the whole particle. Setting it on an existing
+ /// particle is usually a mistake.
+ Particle withParticleId(Barcode particleId) const {
+ Particle p = *this;
+ p.m_particleId = particleId;
+ return p;
+ }
+
/// Set the process type that generated this particle.
Particle &setProcess(ProcessType proc) { return m_process = proc, *this; }
/// Set the space-time position four-vector.
@@ -75,21 +88,21 @@ class Particle {
}
/// Set the direction three-vector
Particle &setDirection(const Vector3 &direction) {
- m_direction = direction;
- m_direction.normalize();
+ m_unitDirection = direction;
+ m_unitDirection.normalize();
return *this;
}
/// Set the direction three-vector from scalar components.
Particle &setDirection(Scalar dx, Scalar dy, Scalar dz) {
- m_direction[0] = dx;
- m_direction[1] = dy;
- m_direction[2] = dz;
- m_direction.normalize();
+ m_unitDirection[0] = dx;
+ m_unitDirection[1] = dy;
+ m_unitDirection[2] = dz;
+ m_unitDirection.normalize();
return *this;
}
/// Set the absolute momentum.
- Particle &setMomentum(Scalar momentum) {
- m_momentum = momentum;
+ Particle &setAbsMomentum(Scalar absMomentum) {
+ m_absMomentum = absMomentum;
return *this;
}
/// Change the energy by the given amount.
@@ -98,30 +111,30 @@ class Particle {
/// would result in an unphysical value, the particle is put to rest, i.e.
/// its absolute momentum is set to zero.
Particle &correctEnergy(Scalar delta) {
- const auto newEnergy = std::hypot(m_mass, m_momentum) + delta;
+ const auto newEnergy = std::hypot(m_mass, m_absMomentum) + delta;
if (newEnergy <= m_mass) {
- m_momentum = Scalar(0);
+ m_absMomentum = Scalar(0);
} else {
- m_momentum = std::sqrt(newEnergy * newEnergy - m_mass * m_mass);
+ m_absMomentum = std::sqrt(newEnergy * newEnergy - m_mass * m_mass);
}
return *this;
}
/// Particle identifier within an event.
- Barcode id() const { return m_id; }
+ constexpr Barcode particleId() const { return m_particleId; }
/// Which type of process generated this particle.
- ProcessType process() const { return m_process; }
+ constexpr ProcessType process() const { return m_process; }
/// PDG particle number that identifies the type.
- Acts::PdgParticle pdg() const { return m_pdg; }
+ constexpr Acts::PdgParticle pdg() const { return m_pdg; }
/// Particle charge.
- Scalar charge() const { return m_charge; }
+ constexpr Scalar charge() const { return m_charge; }
/// Particle mass.
- Scalar mass() const { return m_mass; }
+ constexpr Scalar mass() const { return m_mass; }
/// Space-time position four-vector.
///
/// The component order is [x,y,z,t].
- const Vector4 &position4() const { return m_position4; }
+ constexpr const Vector4 &position4() const { return m_position4; }
/// Three-position, i.e. spatial coordinates without the time.
auto position() const { return m_position4.head<3>(); }
/// Time coordinate.
@@ -132,36 +145,33 @@ class Particle {
Vector4 momentum4() const {
Vector4 mom4;
// stored direction is always normalized
- mom4[0] = m_momentum * m_direction[0];
- mom4[1] = m_momentum * m_direction[1];
- mom4[2] = m_momentum * m_direction[2];
+ mom4[0] = m_absMomentum * m_unitDirection[0];
+ mom4[1] = m_absMomentum * m_unitDirection[1];
+ mom4[2] = m_absMomentum * m_unitDirection[2];
mom4[3] = energy();
return mom4;
}
- /// Three-direction, i.e. the normalized momentum three-vector.
- const Vector3 &direction() const { return m_direction; }
+ /// Unit three-direction, i.e. the normalized momentum three-vector.
+ const Vector3 &unitDirection() const { return m_unitDirection; }
+ /// Absolute momentum in the x-y plane.
+ Scalar transverseMomentum() const {
+ return m_absMomentum * m_unitDirection.head<2>().norm();
+ }
/// Absolute momentum.
- Scalar momentum() const { return m_momentum; }
+ constexpr Scalar absMomentum() const { return m_absMomentum; }
/// Total energy, i.e. norm of the four-momentum.
- Scalar energy() const { return std::hypot(m_mass, m_momentum); }
-
- /// Charge over absolute momentum.
- Scalar chargeOverMomentum() const { return m_charge / m_momentum; }
- /// Relativistic velocity.
- Scalar beta() const { return m_momentum / energy(); }
- /// Relativistic gamma factor.
- Scalar gamma() const { return std::hypot(1, m_momentum / m_mass); }
+ Scalar energy() const { return std::hypot(m_mass, m_absMomentum); }
/// Check if the particle is alive, i.e. is not at rest.
- operator bool() const { return Scalar(0) < m_momentum; }
+ constexpr operator bool() const { return Scalar(0) < m_absMomentum; }
/// Check if the particle is dead, i.e is at rest.
- bool operator!() const { return m_momentum <= Scalar(0); }
+ constexpr bool operator!() const { return m_absMomentum <= Scalar(0); }
/// Set the material that the particle has passed.
///
/// @param pathX0 passed material measured in radiation lengths
/// @param pathL0 passed thickness measured in interaction lengths
- Particle &setMaterialPassed(Scalar pathX0, Scalar pathL0) {
+ constexpr Particle &setMaterialPassed(Scalar pathX0, Scalar pathL0) {
m_pathX0 = pathX0;
m_pathL0 = pathL0;
return *this;
@@ -170,24 +180,24 @@ class Particle {
///
/// @param limitX0 maximum radiation lengths the particle can pass
/// @param limitL0 maximum interaction lengths the particle can pass
- Particle &setMaterialLimits(Scalar limitX0, Scalar limitL0) {
+ constexpr Particle &setMaterialLimits(Scalar limitX0, Scalar limitL0) {
m_limitX0 = limitX0;
m_limitL0 = limitL0;
return *this;
}
/// The passed material measured in radiation lengths.
- Scalar pathInX0() const { return m_pathX0; }
+ constexpr Scalar pathInX0() const { return m_pathX0; }
/// The passed material measured in interaction lengths.
- Scalar pathInL0() const { return m_pathL0; }
+ constexpr Scalar pathInL0() const { return m_pathL0; }
/// The maximum radation length the particle is allowed to pass.
- Scalar pathLimitX0() const { return m_limitX0; }
+ constexpr Scalar pathLimitX0() const { return m_limitX0; }
/// The maximum interaction length the particle is allowed to pass.
- Scalar pathLimitL0() const { return m_limitL0; }
+ constexpr Scalar pathLimitL0() const { return m_limitL0; }
private:
// identity, i.e. things that do not change over the particle lifetime.
/// Particle identifier within the event.
- Barcode m_id;
+ Barcode m_particleId;
/// Process type specifier.
ProcessType m_process = ProcessType::eUndefined;
/// PDG particle number.
@@ -196,8 +206,8 @@ class Particle {
Scalar m_charge = Scalar(0);
Scalar m_mass = Scalar(0);
// kinematics, i.e. things that change over the particle lifetime.
- Vector3 m_direction = Vector3::UnitZ();
- Scalar m_momentum = Scalar(0);
+ Vector3 m_unitDirection = Vector3::UnitZ();
+ Scalar m_absMomentum = Scalar(0);
Vector4 m_position4 = Vector4::Zero();
// simulation-specific X0/L0 information and limits
// these values are here to simplify the simulation of (nuclear) interactions.
@@ -213,4 +223,6 @@ class Particle {
Scalar m_limitL0 = std::numeric_limits<Scalar>::max();
};
+std::ostream &operator<<(std::ostream &os, const Particle &particle);
+
} // namespace ActsFatras
diff --git a/Fatras/include/ActsFatras/Kernel/Interactor.hpp b/Fatras/include/ActsFatras/Kernel/Interactor.hpp
index 5b77374ff5188e00f7d9cc920bcb546bc8f39a93..2d53d75e2aa62f4fa85c2aa4d4d5ef51602caaa8 100644
--- a/Fatras/include/ActsFatras/Kernel/Interactor.hpp
+++ b/Fatras/include/ActsFatras/Kernel/Interactor.hpp
@@ -100,7 +100,7 @@ struct Interactor {
.setPosition4(stepper.position(state.stepping),
stepper.time(state.stepping))
.setDirection(stepper.direction(state.stepping))
- .setMomentum(stepper.momentum(state.stepping));
+ .setAbsMomentum(stepper.momentum(state.stepping));
// we want to keep the particle state before and after the interaction.
// since the particle is modified in-place we need a copy.
auto after = before;
@@ -110,7 +110,7 @@ struct Interactor {
Acts::Vector2D local;
// TODO what to do in case of invalid return value?
surface.globalToLocal(state.geoContext, before.position(),
- before.direction(), local);
+ before.unitDirection(), local);
Acts::MaterialProperties slab =
surface.surfaceMaterial()->materialProperties(local);
@@ -120,7 +120,8 @@ struct Interactor {
auto normal = surface.normal(state.geoContext, local);
// dot-product(unit normal, direction) = cos(incidence angle)
// particle direction is normalized, not sure about surface normal
- auto cosIncidenceInv = normal.norm() / normal.dot(before.direction());
+ auto cosIncidenceInv =
+ normal.norm() / normal.dot(before.unitDirection());
slab.scaleThickness(cosIncidenceInv);
// TODO how to communicate the break condition to the propagator?
physics(*generator, slab, after, result.generatedParticles);
@@ -143,15 +144,15 @@ struct Interactor {
result.particle = after;
if (selectHitSurface(surface)) {
result.hits.emplace_back(
- surface.geoID(), before.id(),
+ surface.geoID(), before.particleId(),
// the interaction could potentially modify the particle position
Hit::Scalar(0.5) * (before.position4() + after.position4()),
before.momentum4(), after.momentum4(), result.hits.size());
}
// continue the propagation with the modified parameters
- stepper.update(state.stepping, after.position(), after.direction(),
- after.momentum(), after.time());
+ stepper.update(state.stepping, after.position(), after.unitDirection(),
+ after.absMomentum(), after.time());
}
/// Pure observer interface. Does not apply to the Fatras simulator.
diff --git a/Fatras/include/ActsFatras/Kernel/Simulator.hpp b/Fatras/include/ActsFatras/Kernel/Simulator.hpp
index 5a99b04c26868158cf737097e95aac24ff7481aa..c21849164488fd3b31c6ea7005c5c78ced6129c1 100644
--- a/Fatras/include/ActsFatras/Kernel/Simulator.hpp
+++ b/Fatras/include/ActsFatras/Kernel/Simulator.hpp
@@ -8,173 +8,319 @@
#pragma once
-#include <optional>
+#include <algorithm>
+#include <cassert>
+#include <iterator>
+#include <memory>
-#include "Acts/EventData/NeutralParameters.hpp"
-#include "Acts/EventData/TrackParameters.hpp"
+#include "Acts/EventData/ChargePolicy.hpp"
+#include "Acts/EventData/SingleCurvilinearTrackParameters.hpp"
+#include "Acts/Geometry/GeometryContext.hpp"
+#include "Acts/MagneticField/MagneticFieldContext.hpp"
#include "Acts/Propagator/AbortList.hpp"
#include "Acts/Propagator/ActionList.hpp"
#include "Acts/Propagator/Propagator.hpp"
#include "Acts/Propagator/detail/DebugOutputActor.hpp"
#include "Acts/Propagator/detail/StandardAborters.hpp"
#include "Acts/Utilities/Logger.hpp"
+#include "Acts/Utilities/Result.hpp"
+#include "ActsFatras/EventData/Hit.hpp"
+#include "ActsFatras/EventData/Particle.hpp"
+#include "ActsFatras/Kernel/Interactor.hpp"
namespace ActsFatras {
-struct VoidDetector {};
-
-/// @brief Fatras simulator
-///
-/// This is called from a Fatras steering algorithm
-/// Per call, the generator is provided which
-/// is then used to create a Acts propagator plugin.
+/// Single particle simulator with a fixed propagator and physics list.
///
-/// @tparam charged_propagator_t Type of the propagator for charged particles
-/// @tparam charged_selector_t Type of the slector (list) for charged particles
-/// @tparam charged_interactor_t Type of the dresser for chargerd particles
-///
-/// @tparam neutral_propagator_t Type of the propagator for neutral particles
-/// @tparam neutral_selector_t Type of the slector (list) for neutral particles
-/// @tparam neutral_interactor_t Type of the dresser for neutral particles
-template <typename charged_propagator_t, typename charged_selector_t,
- typename charged_interactor_t, typename neutral_propagator_t,
- typename neutral_selector_t, typename neutral_interactor_t>
-struct Simulator {
- Simulator(charged_propagator_t chpropagator, neutral_propagator_t npropagator)
- : chargedPropagator(std::move(chpropagator)),
- neutralPropagator(std::move(npropagator)),
- mlogger(Acts::getDefaultLogger("Simulator", Acts::Logging::INFO)) {}
+/// @tparam propagator_t is the type of the underlying propagator
+/// @tparam physics_list_t is the type of the simulated physics list
+/// @tparam hit_surface_selector_t is the type that selects hit surfaces
+template <typename propagator_t, typename physics_list_t,
+ typename hit_surface_selector_t>
+struct ParticleSimulator {
+ /// How and within which geometry to propagate the particle.
+ propagator_t propagator;
+ /// What should be simulated. Will be copied to the per-call interactor.
+ physics_list_t physics;
+ /// Where hits are registiered. Will be copied to the per-call interactor.
+ hit_surface_selector_t selectHitSurface;
+ /// Local logger for debug output.
+ std::shared_ptr<const Acts::Logger> localLogger = nullptr;
- using PhysicsList_t = typename charged_interactor_t::PhysicsList_t;
- charged_propagator_t chargedPropagator;
- charged_selector_t chargedSelector;
- PhysicsList_t physicsList;
+ /// Construct the simulator with the underlying propagator.
+ ParticleSimulator(propagator_t &&propagator_, Acts::Logging::Level lvl)
+ : propagator(propagator_),
+ localLogger(Acts::getDefaultLogger("Simulator", lvl)) {}
+
+ /// Provide access to the local logger instance, e.g. for logging macros.
+ const Acts::Logger &logger() const { return *localLogger; }
+
+ /// Simulate a single particle without secondaries.
+ ///
+ /// @param geoCtx is the geometry context to access surface geometries
+ /// @param magCtx is the magnetic field context to access field values
+ /// @param generator is the random number generator
+ /// @param particle is the initial particle state
+ /// @returns the result of the corresponding Interactor propagator action.
+ ///
+ /// @tparam generator_t is the type of the random number generator
+ template <typename generator_t>
+ Acts::Result<typename Interactor<generator_t, physics_list_t,
+ hit_surface_selector_t>::result_type>
+ simulate(const Acts::GeometryContext &geoCtx,
+ const Acts::MagneticFieldContext &magCtx, generator_t &generator,
+ const Particle &particle) const {
+ assert(localLogger and "Missing local logger");
- neutral_propagator_t neutralPropagator;
- neutral_selector_t neutralSelector;
+ // propagator-related additional types
+ using Interact =
+ Interactor<generator_t, physics_list_t, hit_surface_selector_t>;
+ using InteractResult = typename Interact::result_type;
+ using Actions = Acts::ActionList<Interact, Acts::detail::DebugOutputActor>;
+ using Abort = Acts::AbortList<Acts::detail::EndOfWorldReached>;
+ using PropagatorOptions = Acts::PropagatorOptions<Actions, Abort>;
- VoidDetector detector;
+ // Construct per-call options.
+ PropagatorOptions options(geoCtx, magCtx);
+ options.absPdgCode = particle.pdg();
+ options.mass = particle.mass();
+ options.debug = localLogger->doPrint(Acts::Logging::Level::DEBUG);
+ // setup the interactor as part of the propagator options
+ auto &interactor = options.actionList.template get<Interact>();
+ interactor.generator = &generator;
+ interactor.physics = physics;
+ interactor.selectHitSurface = selectHitSurface;
+ interactor.particle = particle;
- std::shared_ptr<const Acts::Logger> mlogger = nullptr;
+ // run with a start parameter type depending on the particle charge.
+ // TODO make track parameters consistently constructible regardless
+ // of the charge policy and template the class on the parameter.
+ if (particle.charge() != 0) {
+ Acts::SingleCurvilinearTrackParameters<Acts::ChargedPolicy> start(
+ std::nullopt, particle.position(),
+ particle.absMomentum() * particle.unitDirection(), particle.charge(),
+ particle.time());
+ auto result = propagator.propagate(start, options);
+ if (result.ok()) {
+ return result.value().template get<InteractResult>();
+ } else {
+ return result.error();
+ }
+ } else {
+ Acts::SingleCurvilinearTrackParameters<Acts::NeutralPolicy> start(
+ std::nullopt, particle.position(),
+ particle.absMomentum() * particle.unitDirection(), particle.time());
+ auto result = propagator.propagate(start, options);
+ if (result.ok()) {
+ return result.value().template get<InteractResult>();
+ } else {
+ return result.error();
+ }
+ }
+ }
+};
- bool debug = false;
+/// Fatras multi-particle simulator.
+///
+/// @tparam charged_selector_t Callable selector type for charged particles
+/// @tparam charged_simulator_t Single particle simulator for charged particles
+/// @tparam neutral_selector_t Callable selector type for neutral particles
+/// @tparam neutral_simulator_t Single particle simulator for neutral particles
+template <typename charged_selector_t, typename charged_simulator_t,
+ typename neutral_selector_t, typename neutral_simulator_t>
+struct Simulator {
+ charged_selector_t selectCharged;
+ neutral_selector_t selectNeutral;
+ charged_simulator_t charged;
+ neutral_simulator_t neutral;
- /// Private access to the logging instance
- const Acts::Logger &logger() const { return *mlogger; }
+ /// Construct from the single charged/neutral particle simulators.
+ Simulator(charged_simulator_t &&charged_, neutral_simulator_t &&neutral_)
+ : charged(std::move(charged_)), neutral(std::move(neutral_)) {}
- /// @brief call operator to the simulator
+ /// Simulate multiple particles and generated secondaries.
///
- /// @tparam context_t Type pf the context object
- /// @tparam generator_t Type of the generator object
- /// @tparam event_collection_t Type of the event collection
- /// @tparam hit_collection_t Type of the hit collection, needs insert()
+ /// @param geoCtx is the geometry context to access surface geometries
+ /// @param magCtx is the magnetic field context to access field values
+ /// @param inputParticles contains all particles that should be simulated
+ /// @param simulatedParticlesInitial contains initial particle states
+ /// @param simulatedParticlesFinal contains final particle states
+ /// @param hits contains all generated hits
///
- /// @param fatrasContext is the event-bound context
- /// @param fatrasGenerator is the event-bound random generator
- /// @param fatrasEvent is the truth event collection
- /// @param fatrasHits is the hit collection
- template <typename context_t, typename generator_t,
- typename event_collection_t, typename hit_collection_t>
- void operator()(context_t &fatrasContext, generator_t &fatrasGenerator,
- event_collection_t &fatrasEvent,
- hit_collection_t &fatrasHits) const {
- // if screen output is required
- typedef Acts::detail::DebugOutputActor DebugOutput;
-
- // Action list, abort list and options
- typedef Acts::ActionList<charged_interactor_t, DebugOutput>
- ChargedActionList;
- typedef Acts::AbortList<Acts::detail::EndOfWorldReached> ChargedAbortList;
- typedef Acts::PropagatorOptions<ChargedActionList, ChargedAbortList>
- ChargedOptions;
-
- // Action list, abort list and
- typedef Acts::ActionList<neutral_interactor_t, DebugOutput>
- NeutralActionList;
- typedef Acts::AbortList<Acts::detail::EndOfWorldReached> NeutralAbortList;
- typedef Acts::PropagatorOptions<NeutralActionList, NeutralAbortList>
- NeutralOptions;
-
- // loop over the input events
- // -> new secondaries will just be attached to that
- for (auto &vertex : fatrasEvent) {
- // take care here, the simulation can change the
- // particle collection
- for (std::size_t i = 0; i < vertex.outgoing.size(); i++) {
- // create a local copy since the collection can reallocate and
- // invalidate any reference.
- auto particle = vertex.outgoing[i];
- // charged particle detected and selected
- if (chargedSelector(detector, particle)) {
- // Need to construct them per call to set the particle
- // Options and configuration
- ChargedOptions chargedOptions(fatrasContext.geoContext,
- fatrasContext.magFieldContext);
- chargedOptions.debug = debug;
- // Get the charged interactor
- auto &chargedInteractor =
- chargedOptions.actionList.template get<charged_interactor_t>();
- // Result type typedef
- typedef typename charged_interactor_t::result_type ChargedResult;
- // Set the generator to guarantee event consistent entires
- chargedInteractor.generator = &fatrasGenerator;
- // Put all the additional information into the interactor
- chargedInteractor.initialParticle = particle;
- // Set the physics list
- chargedInteractor.physicsList = physicsList;
- // Create the kinematic start parameters
- Acts::CurvilinearParameters start(std::nullopt, particle.position(),
- particle.momentum(), particle.q(),
- particle.time());
- // Run the simulation
- const auto &result =
- chargedPropagator.propagate(start, chargedOptions).value();
- const auto &fatrasResult = result.template get<ChargedResult>();
- // a) Handle the hits
- // hits go to the hit collection, particle go to the particle
- // collection
- for (const auto &fHit : fatrasResult.simulatedHits) {
- fatrasHits.insert(fHit);
- }
- // b) deal with the particles
- const auto &simparticles = fatrasResult.outgoing;
- vertex.outgoing_insert(simparticles);
- // c) screen output if requested
- if (debug) {
- auto &fatrasDebug = result.template get<DebugOutput::result_type>();
- ACTS_INFO(fatrasDebug.debugString);
+ /// This takes all input particles and simulates those passing the selection
+ /// using the appropriate simulator. All selected particle states including
+ /// additional ones generated from interactions are stored in separate output
+ /// containers; both the initial state at the production vertex and the final
+ /// state after propagation are stored. Hits generated from selected input and
+ /// generated particles are stored in the hit container.
+ ///
+ /// @warning Particle-hit association is based on particle ids generated
+ /// during the simulation. This requires that all input particles
+ /// **must** have generation and sub-particle number set to zero.
+ ///
+ /// @tparam generator_t is the type of the random number generator
+ /// @tparam input_particles_t is a Container for particles
+ /// @tparam output_particles_t is a SequenceContainer for particles
+ /// @tparam hits_t is a SequenceContainer for hits
+ template <typename generator_t, typename input_particles_t,
+ typename output_particles_t, typename hits_t>
+ Acts::Result<void> simulate(const Acts::GeometryContext &geoCtx,
+ const Acts::MagneticFieldContext &magCtx,
+ generator_t &generator,
+ const input_particles_t &inputParticles,
+ output_particles_t &simulatedParticlesInitial,
+ output_particles_t &simulatedParticlesFinal,
+ hits_t &hits) const {
+ assert(
+ (simulatedParticlesInitial.size() == simulatedParticlesFinal.size()) and
+ "Inconsistent initial sizes of the simulated particle containers");
+
+ for (const Particle &inputParticle : inputParticles) {
+ if (not selectParticle(inputParticle)) {
+ continue;
+ }
+ // required to allow correct particle id numbering for secondaries later
+ if ((inputParticle.particleId().generation() != 0u) or
+ (inputParticle.particleId().subParticle() != 0u)) {
+ // TODO add meaningfull error code
+ return std::error_code(-1, std::generic_category());
+ }
+
+ // Do a *depth-first* simulation of the particle and its secondaries,
+ // i.e. we simulate all secondaries, tertiaries, ... before simulating
+ // the next primary particle. Use the end of the output container as
+ // a queue to store particles that should be simulated.
+ //
+ // WARNING the initial particle output container will be modified during
+ // iteration as new secondaries are added directly to it. to avoid
+ // issues, access must always occur via indices.
+ auto iinitial = simulatedParticlesInitial.size();
+ simulatedParticlesInitial.push_back(inputParticle);
+ for (; iinitial < simulatedParticlesInitial.size(); ++iinitial) {
+ // only simulatable particles are pushed to the container
+ // no additional check is necessary here.
+ const auto &initialParticle = simulatedParticlesInitial[iinitial];
+
+ if (selectCharged(initialParticle)) {
+ auto result =
+ charged.simulate(geoCtx, magCtx, generator, initialParticle);
+ if (not result.ok()) {
+ // do not keep unsimulated/ failed particles in the output
+ simulatedParticlesInitial.resize(iinitial);
+ return result.error();
}
- } else if (neutralSelector(detector, particle)) {
- // Options and configuration
- NeutralOptions neutralOptions(fatrasContext.geoContext,
- fatrasContext.magFieldContext);
- neutralOptions.debug = debug;
- // Get the charged interactor
- auto &neutralInteractor =
- neutralOptions.actionList.template get<neutral_interactor_t>();
- // Result type typedef
- typedef typename neutral_interactor_t::result_type NeutralResult;
- // Set the generator to guarantee event consistent entires
- neutralInteractor.generator = &fatrasGenerator;
- // Put all the additional information into the interactor
- neutralInteractor.initialParticle = particle;
- // Create the kinematic start parameters
- Acts::NeutralCurvilinearParameters start(
- std::nullopt, particle.position(), particle.momentum(), 0.);
- const auto &result =
- neutralPropagator.propagate(start, neutralOptions).value();
- auto &fatrasResult = result.template get<NeutralResult>();
- // a) deal with the particles
- const auto &simparticles = fatrasResult.outgoing;
- vertex.outgoing_insert(simparticles);
- // b) screen output if requested
- if (debug) {
- auto &fatrasDebug = result.template get<DebugOutput::result_type>();
- ACTS_INFO(fatrasDebug.debugString);
+ copyOutputs(result.value(), simulatedParticlesInitial,
+ simulatedParticlesFinal, hits);
+ } else if (selectNeutral(initialParticle)) {
+ auto result =
+ neutral.simulate(geoCtx, magCtx, generator, initialParticle);
+ if (not result.ok()) {
+ // do not keep unsimulated/ failed particles in the output
+ simulatedParticlesInitial.resize(iinitial);
+ return result.error();
}
- } // neutral processing
- } // loop over particles
- } // loop over events
+ copyOutputs(result.value(), simulatedParticlesInitial,
+ simulatedParticlesFinal, hits);
+ }
+
+ // since physics processes are independent, there can be particle id
+ // collisions within the generated secondaries. they can be resolved by
+ // renumbering within each sub-particle generation. this must happen
+ // before the particle is simulated since the particle id is used to
+ // associate generated hits back to the particle.
+ renumberTailParticleIds(simulatedParticlesInitial, iinitial);
+ }
+ }
+
+ return Acts::Result<void>::success();
+ }
+
+ private:
+ // Select if the particle should be simulated at all.
+ //
+ // This also enforces mutual-exclusivity of the two charge selections. If both
+ // charge selections evaluate true, they are probably not setup correctly and
+ // not simulating them at all is a reasonable fall-back.
+ bool selectParticle(const Particle &particle) const {
+ const bool isValidCharged = selectCharged(particle);
+ const bool isValidNeutral = selectNeutral(particle);
+ assert(not(isValidCharged and isValidNeutral) and
+ "Charge selectors are not mutually exclusive");
+ return isValidCharged xor isValidNeutral;
+ }
+
+ /// Copy Interactor results to output containers.
+ ///
+ /// @tparam interactor_result_t is an Interactor result struct
+ /// @tparam particles_t is a SequenceContainer for particles
+ /// @tparam hits_t is a SequenceContainer for hits
+ template <typename interactor_result_t, typename particles_t, typename hits_t>
+ void copyOutputs(const interactor_result_t &result,
+ particles_t &particlesInitial, particles_t &particlesFinal,
+ hits_t &hits) const {
+ // initial particle state was already pushed to the container before
+ // store final particle state at the end of the simulation
+ particlesFinal.push_back(result.particle);
+ // move generated secondaries that should be simulated to the output
+ std::copy_if(
+ result.generatedParticles.begin(), result.generatedParticles.end(),
+ std::back_inserter(particlesInitial),
+ [this](const Particle &particle) { return selectParticle(particle); });
+ std::copy(result.hits.begin(), result.hits.end(), std::back_inserter(hits));
+ }
+
+ /// Renumber particle ids in the tail of the container.
+ ///
+ /// Ensures particle ids are unique by modifying the sub-particle number
+ /// within each generation.
+ ///
+ /// @param particles particle container in which particles are renumbered
+ /// @param lastValid index of the last particle with a valid particle id
+ ///
+ /// @tparam particles_t is a SequenceContainer for particles
+ ///
+ /// @note This function assumes that all particles in the tail have the same
+ /// vertex numbers (primary/secondary) and particle number and are
+ /// ordered according to their generation number.
+ ///
+ template <typename particles_t>
+ static void renumberTailParticleIds(particles_t &particles,
+ std::size_t lastValid) {
+ // iterate over adjacent pairs; potentially modify the second element.
+ // assume e.g. a primary particle 2 with generation=subparticle=0 that
+ // generates two secondaries during simulation. we have the following
+ // ids (decoded as vertex|particle|generation|subparticle)
+ //
+ // v|2|0|0, v|2|1|0, v|2|1|0
+ //
+ // where v represents the vertex numbers. this will be renumbered to
+ //
+ // v|2|0|0, v|2|1|0, v|2|1|1
+ //
+ // if each secondary then generates two tertiaries we could have e.g.
+ //
+ // v|2|0|0, v|2|1|0, v|2|1|1, v|2|2|0, v|2|2|1, v|2|2|0, v|2|2|1
+ //
+ // which would then be renumbered to
+ //
+ // v|2|0|0, v|2|1|0, v|2|1|1, v|2|2|0, v|2|2|1, v|2|2|2, v|2|2|3
+ //
+ for (auto j = lastValid; (j + 1u) < particles.size(); ++j) {
+ const auto prevId = particles[j].particleId();
+ auto currId = particles[j + 1u].particleId();
+ // NOTE primary/secondary vertex and particle are assumed to be equal
+ // only renumber within one generation
+ if (prevId.generation() != currId.generation()) {
+ continue;
+ }
+ // ensure the sub-particle is strictly monotonic within a generation
+ if (prevId.subParticle() < currId.subParticle()) {
+ continue;
+ }
+ // sub-particle numbering must be non-zero
+ currId.setSubParticle(prevId.subParticle() + 1u);
+ particles[j + 1u] = particles[j + 1u].withParticleId(currId);
+ }
}
};
diff --git a/Fatras/include/ActsFatras/Physics/EnergyLoss/BetheBloch.hpp b/Fatras/include/ActsFatras/Physics/EnergyLoss/BetheBloch.hpp
index 8424547bdc115dbbb48870f600d59750f8622a32..07966800826c11ad8692b4c6d1b286197e41d721 100644
--- a/Fatras/include/ActsFatras/Physics/EnergyLoss/BetheBloch.hpp
+++ b/Fatras/include/ActsFatras/Physics/EnergyLoss/BetheBloch.hpp
@@ -43,7 +43,7 @@ struct BetheBloch {
// compute energy loss distribution parameters
const auto pdg = particle.pdg();
const auto m = particle.mass();
- const auto qOverP = particle.chargeOverMomentum();
+ const auto qOverP = particle.charge() / particle.absMomentum();
const auto q = particle.charge();
// most probable value
const auto energyLoss =
diff --git a/Fatras/include/ActsFatras/Physics/Scattering/GaussianMixture.hpp b/Fatras/include/ActsFatras/Physics/Scattering/GaussianMixture.hpp
index ceed5c179c5840f57104e2dbc76c8c4a303390da..a8078fbeeeb169a4c412bbe6cd836e229c2c3a7d 100644
--- a/Fatras/include/ActsFatras/Physics/Scattering/GaussianMixture.hpp
+++ b/Fatras/include/ActsFatras/Physics/Scattering/GaussianMixture.hpp
@@ -11,10 +11,10 @@
#include <random>
#include "Acts/Material/Interactions.hpp"
-#include "Acts/Material/MaterialProperties.hpp"
-#include "ActsFatras/EventData/Particle.hpp"
+#include "ActsFatras/Physics/Scattering/detail/Scattering.hpp"
namespace ActsFatras {
+namespace detail {
/// Generate scattering angles using a Gaussian mixture model.
struct GaussianMixture {
@@ -44,8 +44,8 @@ struct GaussianMixture {
Particle &particle) const {
/// Calculate the highland formula first
double sigma = Acts::computeMultipleScatteringTheta0(
- slab, particle.pdg(), particle.mass(), particle.chargeOverMomentum(),
- particle.charge());
+ slab, particle.pdg(), particle.mass(),
+ particle.charge() / particle.absMomentum(), particle.charge());
double sigma2 = sigma * sigma;
// Gauss distribution, will be sampled with generator
@@ -55,8 +55,11 @@ struct GaussianMixture {
// Now correct for the tail fraction
// d_0'
- double beta2 = particle.beta() * particle.beta();
- double dprime = slab.thickness() / (slab.material().X0() * beta2);
+ // beta² = (p/E)² = p²/(p² + m²) = 1/(1 + (m/p)²)
+ // 1/beta² = 1 + (m/p)²
+ double mOverP = particle.mass() / particle.absMomentum();
+ double beta2inv = 1 + mOverP * mOverP;
+ double dprime = slab.thicknessInX0() * beta2inv;
double log_dprime = std::log(dprime);
// d_0''
double log_dprimeprime =
@@ -76,7 +79,8 @@ struct GaussianMixture {
// G4 optimised / native double Gaussian model
if (optGaussianMixtureG4) {
- sigma2 = 225. * dprime / (particle.momentum() * particle.momentum());
+ sigma2 =
+ 225. * dprime / (particle.absMomentum() * particle.absMomentum());
}
// throw the random number core/tail
if (uniformDist(generator) < epsilon) {
@@ -87,4 +91,8 @@ struct GaussianMixture {
}
};
+} // namespace detail
+
+using GaussianMixtureScattering = detail::Scattering<detail::GaussianMixture>;
+
} // namespace ActsFatras
diff --git a/Fatras/include/ActsFatras/Physics/Scattering/GeneralMixture.hpp b/Fatras/include/ActsFatras/Physics/Scattering/GeneralMixture.hpp
index 022b897a463554bcd536c1fe8463680a5ed859dd..a3df2f254c4931cfde30384cb2e3f760e9e8b11d 100644
--- a/Fatras/include/ActsFatras/Physics/Scattering/GeneralMixture.hpp
+++ b/Fatras/include/ActsFatras/Physics/Scattering/GeneralMixture.hpp
@@ -11,11 +11,11 @@
#include <random>
#include "Acts/Material/Interactions.hpp"
-#include "Acts/Material/MaterialProperties.hpp"
#include "Acts/Utilities/PdgParticle.hpp"
-#include "ActsFatras/EventData/Particle.hpp"
+#include "ActsFatras/Physics/Scattering/detail/Scattering.hpp"
namespace ActsFatras {
+namespace detail {
/// Generate scattering angles using a general mixture model.
///
@@ -53,17 +53,22 @@ struct GeneralMixture {
//----------------------------------------------------------------------------
std::array<double, 4> scattering_params;
// Decide which mixture is best
- double beta2 = (particle.beta() * particle.beta());
+ // beta² = (p/E)² = p²/(p² + m²) = 1/(1 + (m/p)²)
+ // 1/beta² = 1 + (m/p)²
+ // beta = 1/sqrt(1 + (m/p)²)
+ double mOverP = particle.mass() / particle.absMomentum();
+ double beta2Inv = 1 + mOverP * mOverP;
+ double beta = 1 / std::sqrt(beta2Inv);
double tInX0 = slab.thicknessInX0();
- double tob2 = slab.thicknessInX0() / beta2;
+ double tob2 = tInX0 * beta2Inv;
if (tob2 > 0.6 / std::pow(slab.material().Z(), 0.6)) {
// Gaussian mixture or pure Gaussian
if (tob2 > 10) {
- scattering_params = getGaussian(particle.beta(), particle.momentum(),
- tInX0, genMixtureScalor);
+ scattering_params = getGaussian(beta, particle.absMomentum(), tInX0,
+ genMixtureScalor);
} else {
scattering_params =
- getGaussmix(particle.beta(), particle.momentum(), tInX0,
+ getGaussmix(beta, particle.absMomentum(), tInX0,
slab.material().Z(), genMixtureScalor);
}
// Simulate
@@ -71,7 +76,7 @@ struct GeneralMixture {
} else {
// Semigaussian mixture - get parameters
auto scattering_params_sg =
- getSemigauss(particle.beta(), particle.momentum(), tInX0,
+ getSemigauss(beta, particle.absMomentum(), tInX0,
slab.material().Z(), genMixtureScalor);
// Simulate
theta = semigauss(generator, scattering_params_sg);
@@ -80,8 +85,8 @@ struct GeneralMixture {
// for electrons we fall back to the Highland (extension)
// return projection factor times sigma times gauss random
const auto theta0 = Acts::computeMultipleScatteringTheta0(
- slab, particle.pdg(), particle.mass(), particle.chargeOverMomentum(),
- particle.charge());
+ slab, particle.pdg(), particle.mass(),
+ particle.charge() / particle.absMomentum(), particle.charge());
theta = std::normal_distribution<double>(0.0, theta0)(generator);
}
// scale from planar to 3d angle
@@ -195,4 +200,8 @@ struct GeneralMixture {
}
};
+} // namespace detail
+
+using GeneralMixtureScattering = detail::Scattering<detail::GeneralMixture>;
+
} // namespace ActsFatras
diff --git a/Fatras/include/ActsFatras/Physics/Scattering/Highland.hpp b/Fatras/include/ActsFatras/Physics/Scattering/Highland.hpp
index 211f2045c8ad3c680fb197e514cada4af3cf4aa2..b593327466327e8a6e23da4d4e5cb7cd5f6ecc46 100644
--- a/Fatras/include/ActsFatras/Physics/Scattering/Highland.hpp
+++ b/Fatras/include/ActsFatras/Physics/Scattering/Highland.hpp
@@ -11,10 +11,10 @@
#include <random>
#include "Acts/Material/Interactions.hpp"
-#include "Acts/Material/MaterialProperties.hpp"
-#include "ActsFatras/EventData/Particle.hpp"
+#include "ActsFatras/Physics/Scattering/detail/Scattering.hpp"
namespace ActsFatras {
+namespace detail {
/// Generate scattering angles using the Highland/PDG parametrization.
///
@@ -35,11 +35,15 @@ struct Highland {
Particle &particle) const {
// compute the planar scattering angle
const auto theta0 = Acts::computeMultipleScatteringTheta0(
- slab, particle.pdg(), particle.mass(), particle.chargeOverMomentum(),
- particle.charge());
+ slab, particle.pdg(), particle.mass(),
+ particle.charge() / particle.absMomentum(), particle.charge());
// draw from the normal distribution representing the 3d angle distribution
return std::normal_distribution<double>(0.0, M_SQRT2 * theta0)(generator);
}
};
+} // namespace detail
+
+using HighlandScattering = detail::Scattering<detail::Highland>;
+
} // namespace ActsFatras
diff --git a/Fatras/include/ActsFatras/Physics/Scattering/Scattering.hpp b/Fatras/include/ActsFatras/Physics/Scattering/detail/Scattering.hpp
similarity index 95%
rename from Fatras/include/ActsFatras/Physics/Scattering/Scattering.hpp
rename to Fatras/include/ActsFatras/Physics/Scattering/detail/Scattering.hpp
index b5fc3bd4fb19a512844c2b5600f3e339beea3005..04d1e51cb69c2a1c2eaf70ceabfeee36d70f4806 100644
--- a/Fatras/include/ActsFatras/Physics/Scattering/Scattering.hpp
+++ b/Fatras/include/ActsFatras/Physics/Scattering/detail/Scattering.hpp
@@ -12,11 +12,12 @@
#include <random>
#include "Acts/Material/MaterialProperties.hpp"
-#include "Acts/Utilities/CurvilinearUnitVectors.hpp"
#include "Acts/Utilities/Definitions.hpp"
+#include "Acts/Utilities/UnitVectors.hpp"
#include "ActsFatras/EventData/Particle.hpp"
namespace ActsFatras {
+namespace detail {
/// Simulate (multiple) scattering using a configurable scattering model.
///
@@ -54,7 +55,7 @@ struct Scattering {
// draw the scattering angle
const auto theta = angle(generator, slab, particle);
- Acts::Vector3D direction = particle.direction();
+ Acts::Vector3D direction = particle.unitDirection();
// construct the combined rotation to the scattered direction
Acts::RotationMatrix3D rotation(
// rotation of the scattering deflector axis relative to the reference
@@ -69,4 +70,5 @@ struct Scattering {
}
};
+} // namespace detail
} // namespace ActsFatras
diff --git a/Fatras/include/ActsFatras/Physics/StandardPhysicsLists.hpp b/Fatras/include/ActsFatras/Physics/StandardPhysicsLists.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ce8867bd7fb0b0f18d5fad66481a05bcfab5951a
--- /dev/null
+++ b/Fatras/include/ActsFatras/Physics/StandardPhysicsLists.hpp
@@ -0,0 +1,66 @@
+// 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 "ActsFatras/Kernel/PhysicsList.hpp"
+#include "ActsFatras/Kernel/Process.hpp"
+#include "ActsFatras/Physics/EnergyLoss/BetheBloch.hpp"
+#include "ActsFatras/Physics/EnergyLoss/BetheHeitler.hpp"
+#include "ActsFatras/Physics/Scattering/Highland.hpp"
+#include "ActsFatras/Selectors/KinematicCasts.hpp"
+#include "ActsFatras/Selectors/PdgSelectors.hpp"
+#include "ActsFatras/Selectors/SelectorHelpers.hpp"
+
+namespace ActsFatras {
+namespace detail {
+/// Select electrons and positrons only.
+using SelectElectronLike = AbsPdgSelector<Acts::PdgParticle::eElectron>;
+/// Select particles above a minimum absolute momentum.
+using SelectPMin = Min<Casts::P>;
+
+/// Highland multiple scattering that applies everywhere.
+using StandardScattering =
+ Process<HighlandScattering, EveryInput, EveryParticle, EveryParticle>;
+/// Ionisation/excitation energy loss with a lower p cut on output particles.
+///
+/// Bethe-Bloch generates no particles and the child selector has no effect.
+using StandardBetheBloch =
+ Process<BetheBloch, EveryInput, SelectPMin, EveryParticle>;
+/// Electron Bremsstrahlung energy loss with a lower p cut on output particles.
+///
+/// Only applies to electrons and positions. Bethe-Heitler generates no
+/// particles and the child selector has no effect.
+using StandardBetheHeitler =
+ Process<BetheHeitler, AsInputSelector<SelectElectronLike>, SelectPMin,
+ EveryParticle>;
+} // namespace detail
+
+/// Electro-magnetic interactions for charged particles.
+///
+/// Scattering must come first so it is computed with the unmodified initial
+/// energy before energy loss is applied.
+///
+/// @warning The list has no cuts on input particle charge or kinematics, i.e.
+/// it relies on the simulator to preselect relevant charged particles
+/// before application.
+/// @todo Bethe-Bloch does not describe electrons; add correct ionisation loss
+/// descriptions for electrons.
+/// @todo Bethe-Heitler is applied after energy loss and thus sees the wrong
+/// input energy.
+using ChargedElectroMagneticPhysicsList =
+ PhysicsList<detail::StandardScattering, detail::StandardBetheBloch,
+ detail::StandardBetheHeitler>;
+
+/// Construct the standard electro-magnetic physics list for charged particles.
+///
+/// @param minimumAbsMomentum lower p cut on output particles
+ChargedElectroMagneticPhysicsList makeChargedElectroMagneticPhysicsList(
+ double minimumAbsMomentum);
+
+} // namespace ActsFatras
diff --git a/Fatras/include/ActsFatras/Selectors/KinematicCasts.hpp b/Fatras/include/ActsFatras/Selectors/KinematicCasts.hpp
index 6229f607600821b6b1e2ae36528b44ab217877ca..67b3dd7051b387180f2de44c28fbd0bd76beda17 100644
--- a/Fatras/include/ActsFatras/Selectors/KinematicCasts.hpp
+++ b/Fatras/include/ActsFatras/Selectors/KinematicCasts.hpp
@@ -40,7 +40,7 @@ struct AbsVz {
struct Eta {
double operator()(const Particle& particle) const {
// particle direction is always normalized, i.e. dz = pz / p
- return std::atanh(particle.direction().z());
+ return std::atanh(particle.unitDirection().z());
}
};
@@ -48,7 +48,7 @@ struct Eta {
struct AbsEta {
double operator()(const Particle& particle) const {
// particle direction is always normalized, i.e. dz = pz / p
- return std::atanh(std::abs(particle.direction().z()));
+ return std::atanh(std::abs(particle.unitDirection().z()));
}
};
@@ -56,15 +56,15 @@ struct AbsEta {
struct Pt {
double operator()(const Particle& particle) const {
// particle direction is always normalized, i.e. dt²+dz²=1 w/ dt²=dx²+dy²
- return particle.momentum() *
- std::hypot(particle.direction().x(), particle.direction().y());
+ return particle.absMomentum() * std::hypot(particle.unitDirection().x(),
+ particle.unitDirection().y());
}
};
/// Retrieve the absolute momentum.
struct P {
double operator()(const Particle& particle) const {
- return particle.momentum();
+ return particle.absMomentum();
}
};
diff --git a/Fatras/include/ActsFatras/Selectors/PdgSelectors.hpp b/Fatras/include/ActsFatras/Selectors/PdgSelectors.hpp
index 2c120c2d3054a91a41406bb481a526a83a4d4943..dbe045be358830864cbb3b0351512d23e172343f 100644
--- a/Fatras/include/ActsFatras/Selectors/PdgSelectors.hpp
+++ b/Fatras/include/ActsFatras/Selectors/PdgSelectors.hpp
@@ -8,8 +8,7 @@
#pragma once
-#include <cstdlib>
-
+#include "Acts/Utilities/PdgParticle.hpp"
#include "ActsFatras/EventData/Particle.hpp"
namespace ActsFatras {
@@ -17,36 +16,38 @@ namespace ActsFatras {
/// Select particles of one specific type.
///
/// Particle and Antiparticle are treated as two separate types.
-template <int Pdg>
+template <Acts::PdgParticle Pdg>
struct PdgSelector {
bool operator()(const Particle &particle) const {
- return (static_cast<int>(particle.pdg()) == Pdg);
+ return (particle.pdg() == Pdg);
}
};
/// Select particles and antiparticles of one specific type.
-template <int Pdg>
+template <Acts::PdgParticle Pdg>
struct AbsPdgSelector {
bool operator()(const Particle &particle) const {
- return (std::abs(static_cast<int>(particle.pdg())) == std::abs(Pdg));
+ return (makeAbsolutePdgParticle(particle.pdg()) ==
+ makeAbsolutePdgParticle(Pdg));
}
};
/// Select all particles except one specific type.
///
/// Particle and Antiparticle are treated as two separate types.
-template <int Pdg>
+template <Acts::PdgParticle Pdg>
struct PdgExcluder {
bool operator()(const Particle &particle) const {
- return (static_cast<int>(particle.pdg()) != Pdg);
+ return (particle.pdg() != Pdg);
}
};
/// Select all particles except for (anti-)particles of one specific type.
-template <int Pdg>
+template <Acts::PdgParticle Pdg>
struct AbsPdgExcluder {
bool operator()(const Particle &particle) const {
- return (std::abs(static_cast<int>(particle.pdg())) != std::abs(Pdg));
+ return (makeAbsolutePdgParticle(particle.pdg()) !=
+ makeAbsolutePdgParticle(Pdg));
}
};
diff --git a/Fatras/src/Particle.cpp b/Fatras/src/Particle.cpp
index 58f91adf5a5b912691d388572927b58db01eada0..34c3283cc8109a4ffc7deb9008b993d6e00d10c6 100644
--- a/Fatras/src/Particle.cpp
+++ b/Fatras/src/Particle.cpp
@@ -10,5 +10,15 @@
#include "ActsFatras/Utilities/ParticleData.hpp"
-ActsFatras::Particle::Particle(Barcode id, Acts::PdgParticle pdg)
- : Particle(id, pdg, findCharge(pdg), findMass(pdg)) {}
+ActsFatras::Particle::Particle(Barcode particleId, Acts::PdgParticle pdg)
+ : Particle(particleId, pdg, findCharge(pdg), findMass(pdg)) {}
+
+std::ostream& ActsFatras::operator<<(std::ostream& os,
+ const ActsFatras::Particle& particle) {
+ os << "particle_id=" << particle.particleId();
+ os << " process=" << particle.process();
+ os << " pdg=" << particle.pdg();
+ os << " q=" << particle.charge();
+ os << " m=" << particle.mass();
+ return os;
+}
diff --git a/Fatras/src/StandardPhysicsLists.cpp b/Fatras/src/StandardPhysicsLists.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..044f1008bf157cf15282df975919fa21310ac9b4
--- /dev/null
+++ b/Fatras/src/StandardPhysicsLists.cpp
@@ -0,0 +1,19 @@
+// 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 "ActsFatras/Physics/StandardPhysicsLists.hpp"
+
+ActsFatras::ChargedElectroMagneticPhysicsList
+ActsFatras::makeChargedElectroMagneticPhysicsList(double minimumAbsMomentum) {
+ ChargedElectroMagneticPhysicsList pl;
+ pl.get<detail::StandardBetheBloch>().selectOutputParticle.valMin =
+ minimumAbsMomentum;
+ pl.get<detail::StandardBetheHeitler>().selectOutputParticle.valMin =
+ minimumAbsMomentum;
+ return pl;
+}
diff --git a/Tests/IntegrationTests/CMakeLists.txt b/Tests/IntegrationTests/CMakeLists.txt
index 8e8b0f55a42b001c8aa543f45d15ddae27bcf035..3198633b9a9a1258e08603ba0ca12e1ca81113b4 100644
--- a/Tests/IntegrationTests/CMakeLists.txt
+++ b/Tests/IntegrationTests/CMakeLists.txt
@@ -34,4 +34,5 @@ add_integrationtest(BVHNavigation BVHNavigationTest.cpp)
add_integrationtest(InterpolatedSolenoidBField InterpolatedSolenoidBFieldTest.cpp)
add_integrationtest(Propagation PropagationTests.cpp)
+add_subdirectory_if(Fatras ACTS_BUILD_FATRAS)
add_subdirectory_if(Legacy ACTS_BUILD_LEGACY)
diff --git a/Tests/IntegrationTests/Fatras/CMakeLists.txt b/Tests/IntegrationTests/Fatras/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..584bb46836e2364ee5a339e12a389d5c53bf8559
--- /dev/null
+++ b/Tests/IntegrationTests/Fatras/CMakeLists.txt
@@ -0,0 +1,3 @@
+set(integrationtest_extra_libraries ActsFatras)
+
+add_integrationtest(FatrasSimulation FatrasSimulationTests.cpp)
diff --git a/Tests/IntegrationTests/Fatras/FatrasSimulationTests.cpp b/Tests/IntegrationTests/Fatras/FatrasSimulationTests.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0b7c41a299b442855c4c5a2f538a3bf6f773b0c1
--- /dev/null
+++ b/Tests/IntegrationTests/Fatras/FatrasSimulationTests.cpp
@@ -0,0 +1,238 @@
+// 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 <algorithm>
+#include <boost/test/data/test_case.hpp>
+#include <boost/test/unit_test.hpp>
+#include <random>
+
+#include "Acts/EventData/NeutralParameters.hpp"
+#include "Acts/EventData/TrackParameters.hpp"
+#include "Acts/MagneticField/ConstantBField.hpp"
+#include "Acts/Propagator/EigenStepper.hpp"
+#include "Acts/Propagator/Navigator.hpp"
+#include "Acts/Propagator/StraightLineStepper.hpp"
+#include "Acts/Tests/CommonHelpers/CylindricalTrackingGeometry.hpp"
+#include "Acts/Utilities/UnitVectors.hpp"
+#include "ActsFatras/Kernel/PhysicsList.hpp"
+#include "ActsFatras/Kernel/Simulator.hpp"
+#include "ActsFatras/Physics/StandardPhysicsLists.hpp"
+#include "ActsFatras/Selectors/ChargeSelectors.hpp"
+#include "ActsFatras/Utilities/ParticleData.hpp"
+
+using namespace Acts::UnitLiterals;
+
+namespace {
+
+/// Mock-up process that splits a particle into two above a momentum threshold.
+struct SplitEnergyLoss {
+ double splitMomentumMin = 5_GeV;
+
+ template <typename generator_t>
+ bool operator()(generator_t&, const Acts::MaterialProperties&,
+ ActsFatras::Particle& particle,
+ std::vector<ActsFatras::Particle>& generated) const {
+ const auto p = particle.absMomentum();
+ if (splitMomentumMin < p) {
+ particle.setAbsMomentum(0.5 * p);
+ const auto pid = particle.particleId().makeDescendant();
+ generated.push_back(particle.withParticleId(pid).setAbsMomentum(0.5 * p));
+ }
+ // never break
+ return false;
+ }
+};
+
+// setup propagator-related types
+// use the default navigation
+using Navigator = Acts::Navigator;
+using MagneticField = Acts::ConstantBField;
+// propagate charged particles numerically in a constant magnetic field
+using ChargedStepper = Acts::EigenStepper<MagneticField>;
+using ChargedPropagator = Acts::Propagator<ChargedStepper, Navigator>;
+// propagate neutral particles with just straight lines
+using NeutralStepper = Acts::StraightLineStepper;
+using NeutralPropagator = Acts::Propagator<NeutralStepper, Navigator>;
+
+// setup simulator-related types
+// the random number generator type
+using Generator = std::ranlux48;
+// all charged particles w/ a mock-up physics list and hits everywhere
+using ChargedSelector = ActsFatras::ChargedSelector;
+using ChargedPhysicsList =
+ ActsFatras::PhysicsList<ActsFatras::detail::StandardScattering,
+ SplitEnergyLoss>;
+using ChargedSimulator =
+ ActsFatras::ParticleSimulator<ChargedPropagator, ChargedPhysicsList,
+ ActsFatras::EverySurface>;
+// all neutral particles w/o physics and no hits
+using NeutralSelector = ActsFatras::NeutralSelector;
+using NeutralSimulator =
+ ActsFatras::ParticleSimulator<NeutralPropagator, ActsFatras::PhysicsList<>,
+ ActsFatras::NoSurface>;
+// full simulator type for charged and neutrals
+using Simulator = ActsFatras::Simulator<ChargedSelector, ChargedSimulator,
+ NeutralSelector, NeutralSimulator>;
+
+// parameters for data-driven test cases
+
+const auto rangePdg = boost::unit_test::data::make({
+ Acts::PdgParticle::eElectron,
+ Acts::PdgParticle::eMuon,
+ Acts::PdgParticle::ePionPlus,
+ Acts::PdgParticle::ePionZero,
+});
+const auto rangePhi = boost::unit_test::data::make({
+ -135_degree,
+ -45_degree,
+ 45_degree,
+ 135_degree,
+});
+const auto rangeEta = boost::unit_test::data::make({
+ -1.0,
+ 0.0,
+ 1.0,
+ 3.0,
+});
+const auto rangeP = boost::unit_test::data::make({
+ 1_GeV,
+ 10_GeV,
+});
+const auto rangeNumParticles = boost::unit_test::data::make({
+ 1,
+ 2,
+});
+
+const auto dataset =
+ rangePdg * rangePhi * rangeEta * rangeP * rangeNumParticles;
+
+// helper functions for tests
+template <typename Container>
+void sortByParticleId(Container& container) {
+ std::sort(container.begin(), container.end(),
+ [](const auto& lhs, const auto& rhs) {
+ return lhs.particleId() < rhs.particleId();
+ });
+}
+template <typename Container>
+bool areParticleIdsUnique(const Container& sortedByParticleId) {
+ // assumes the container is sorted by particle id
+ auto ret =
+ std::adjacent_find(sortedByParticleId.begin(), sortedByParticleId.end(),
+ [](const auto& lhs, const auto& rhs) {
+ return lhs.particleId() == rhs.particleId();
+ });
+ return ret == sortedByParticleId.end();
+}
+template <typename Container, typename Value>
+bool containsParticleId(const Container& sortedByParticleId,
+ const Value& value) {
+ return std::binary_search(sortedByParticleId.begin(),
+ sortedByParticleId.end(), value,
+ [](const auto& lhs, const auto& rhs) {
+ return lhs.particleId() < rhs.particleId();
+ });
+}
+
+} // namespace
+
+BOOST_DATA_TEST_CASE(FatrasSimulation, dataset, pdg, phi, eta, p,
+ numParticles) {
+ using namespace Acts::UnitLiterals;
+
+ Acts::GeometryContext geoCtx;
+ Acts::MagneticFieldContext magCtx;
+ Acts::Logging::Level logLevel = Acts::Logging::Level::DEBUG;
+
+ // construct the example detector
+ Acts::Test::CylindricalTrackingGeometry geoBuilder(geoCtx);
+ auto trackingGeometry = geoBuilder();
+
+ // construct the propagators
+ Navigator navigator(trackingGeometry);
+ ChargedStepper chargedStepper(Acts::ConstantBField(0, 0, 1_T));
+ ChargedPropagator chargedPropagator(std::move(chargedStepper), navigator);
+ NeutralPropagator neutralPropagator(NeutralStepper(), navigator);
+
+ // construct the simulator
+ ChargedSimulator simulatorCharged(std::move(chargedPropagator), logLevel);
+ NeutralSimulator simulatorNeutral(std::move(neutralPropagator), logLevel);
+ Simulator simulator(std::move(simulatorCharged), std::move(simulatorNeutral));
+
+ // prepare simulation call parameters
+ // random number generator
+ Generator generator;
+ // input/ output particle and hits containers
+ std::vector<ActsFatras::Particle> input;
+ std::vector<ActsFatras::Particle> simulatedInitial;
+ std::vector<ActsFatras::Particle> simulatedFinal;
+ std::vector<ActsFatras::Hit> hits;
+
+ // create input particles. particle number should ne non-zero.
+ for (auto i = numParticles; 0 < i; --i) {
+ const auto pid = ActsFatras::Barcode().setVertexPrimary(42).setParticle(i);
+ const auto particle =
+ ActsFatras::Particle(pid, pdg)
+ .setDirection(Acts::makeDirectionUnitFromPhiEta(phi, eta))
+ .setAbsMomentum(p);
+ input.push_back(std::move(particle));
+ }
+ BOOST_TEST_INFO(input.front());
+ BOOST_TEST(input.size() == numParticles);
+
+ // run the simulation
+ auto result = simulator.simulate(geoCtx, magCtx, generator, input,
+ simulatedInitial, simulatedFinal, hits);
+
+ // should always succeed
+ BOOST_TEST(result.ok());
+
+ // ensure simulated particle containers have consistent content
+ BOOST_TEST(simulatedInitial.size() == simulatedFinal.size());
+ for (std::size_t i = 0; i < simulatedInitial.size(); ++i) {
+ const auto& initialParticle = simulatedInitial[i];
+ const auto& finalParticle = simulatedFinal[i];
+ // particle identify should not change during simulation
+ BOOST_TEST(initialParticle.particleId() == finalParticle.particleId());
+ BOOST_TEST(initialParticle.process() == finalParticle.process());
+ BOOST_TEST(initialParticle.pdg() == finalParticle.pdg());
+ BOOST_TEST(initialParticle.charge() == finalParticle.charge());
+ BOOST_TEST(initialParticle.mass() == finalParticle.mass());
+ }
+
+ // we have no particle cuts and should not loose any particles.
+ // might end up with more due to secondaries
+ BOOST_TEST(input.size() <= simulatedInitial.size());
+ BOOST_TEST(input.size() <= simulatedFinal.size());
+ // there should be some hits if we started with a charged particle
+ if (ActsFatras::findCharge(pdg) != 0) {
+ BOOST_TEST(0u < hits.size());
+ }
+
+ // sort all outputs by particle id to simply further tests
+ sortByParticleId(input);
+ sortByParticleId(simulatedInitial);
+ sortByParticleId(simulatedFinal);
+ sortByParticleId(hits);
+
+ // check that all particle ids are unique
+ BOOST_TEST(areParticleIdsUnique(input));
+ BOOST_TEST(areParticleIdsUnique(simulatedInitial));
+ BOOST_TEST(areParticleIdsUnique(simulatedFinal));
+ // hits must necessarily contain particle id duplicates
+ // check that every input particles is simulated
+ for (const auto& particle : input) {
+ BOOST_TEST(containsParticleId(simulatedInitial, particle));
+ BOOST_TEST(containsParticleId(simulatedFinal, particle));
+ }
+ // check that all hits can be associated to a particle
+ for (const auto& hit : hits) {
+ BOOST_TEST(containsParticleId(simulatedInitial, hit));
+ BOOST_TEST(containsParticleId(simulatedFinal, hit));
+ }
+}
diff --git a/Tests/UnitTests/Core/Utilities/CMakeLists.txt b/Tests/UnitTests/Core/Utilities/CMakeLists.txt
index 4ac191617214774cba59999c7b5eb343c542a923..5462c9b65ee126e349d0e4fda653fa90253a77c3 100644
--- a/Tests/UnitTests/Core/Utilities/CMakeLists.txt
+++ b/Tests/UnitTests/Core/Utilities/CMakeLists.txt
@@ -4,7 +4,6 @@ add_unittest(BFieldMapUtilsTests BFieldMapUtilsTests.cpp)
add_unittest(BinningDataTests BinningDataTests.cpp)
add_unittest(BinUtilityTests BinUtilityTests.cpp)
add_unittest(BoundingBoxTest BoundingBoxTest.cpp)
-add_unittest(CurvilinearUnitVectors CurvilinearUnitVectorsTests.cpp)
add_unittest(ExtendableTests ExtendableTests.cpp)
add_unittest(FiniteStateMachineTests FiniteStateMachineTests.cpp)
add_unittest(FrustumTest FrustumTest.cpp)
@@ -21,4 +20,5 @@ add_unittest(RealQuadraticEquationTests RealQuadraticEquationTests.cpp)
add_unittest(ResultTests ResultTests.cpp)
add_unittest(TypeTraitsTest TypeTraitsTest.cpp)
add_unittest(UnitConversionTests UnitConversionTests.cpp)
+add_unittest(UnitVectors UnitVectorsTests.cpp)
add_unittest(VisualizationTests VisualizationTests.cpp)
diff --git a/Tests/UnitTests/Core/Utilities/CurvilinearUnitVectorsTests.cpp b/Tests/UnitTests/Core/Utilities/CurvilinearUnitVectorsTests.cpp
deleted file mode 100644
index 55b4b9d0413cdfb7df4728acfedec17522ba96c2..0000000000000000000000000000000000000000
--- a/Tests/UnitTests/Core/Utilities/CurvilinearUnitVectorsTests.cpp
+++ /dev/null
@@ -1,71 +0,0 @@
-// 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 <boost/test/unit_test.hpp>
-
-#include <limits>
-
-#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
-#include "Acts/Utilities/CurvilinearUnitVectors.hpp"
-
-using Acts::Vector3D;
-
-namespace {
-constexpr auto eps = std::numeric_limits<double>::epsilon();
-
-template <typename Direction, typename RefUnitU, typename RefUnitV>
-void test(const Eigen::MatrixBase<Direction>& direction,
- const Eigen::MatrixBase<RefUnitU>& refU,
- const Eigen::MatrixBase<RefUnitV>& refV) {
- const auto u = Acts::makeCurvilinearUnitU(direction);
- const auto uv = Acts::makeCurvilinearUnitVectors(direction);
- // verify normalization
- CHECK_CLOSE_ABS(u.norm(), 1, eps);
- CHECK_CLOSE_ABS(uv.first.norm(), 1, eps);
- CHECK_CLOSE_ABS(uv.second.norm(), 1, eps);
- // verify orthonormality
- CHECK_SMALL(u.dot(direction), eps);
- CHECK_SMALL(uv.first.dot(direction), eps);
- CHECK_SMALL(uv.second.dot(direction), eps);
- CHECK_SMALL(uv.first.dot(uv.second), eps);
- // verify u is in the x-y plane
- CHECK_SMALL(u[2], eps);
- CHECK_SMALL(uv.first[2], eps);
- // verify references. do not use element-wise comparison to avoid issues with
- // small, epsilon-like differences.
- CHECK_CLOSE_ABS(u.dot(refU), 1, eps);
- CHECK_CLOSE_ABS(uv.first.dot(refU), 1, eps);
- CHECK_CLOSE_ABS(uv.second.dot(refV), 1, eps);
-}
-} // namespace
-
-BOOST_AUTO_TEST_SUITE(CurvilinearUnitVectors)
-
-BOOST_AUTO_TEST_CASE(DirectionTransverse) {
- // curvilinear system w/ direction in the transverse plane
- test(Vector3D(1, 1, 0), Vector3D(-1, 1, 0).normalized(),
- Vector3D(0, 0, 1).normalized());
-}
-
-BOOST_AUTO_TEST_CASE(DirectionPositiveZ) {
- // curvilinear system w/ direction along z
- test(Vector3D(0, 0, 1), Vector3D(1, 0, 0), Vector3D(0, 1, 0));
-}
-
-BOOST_AUTO_TEST_CASE(DirectionNegativeZ) {
- // curvilinear system w/ direction along z
- test(Vector3D(0, 0, -1), Vector3D(1, 0, 0), Vector3D(0, -1, 0));
-}
-
-BOOST_AUTO_TEST_CASE(DirectionCloseToZ) {
- // curvilinear system w/ direction close to z
- test(Vector3D(0, 32 * eps, 1 - 32 * eps), Vector3D(-1, 0, 0),
- Vector3D(0, -1, 32 * eps));
-}
-
-BOOST_AUTO_TEST_SUITE_END()
diff --git a/Tests/UnitTests/Core/Utilities/UnitVectorsTests.cpp b/Tests/UnitTests/Core/Utilities/UnitVectorsTests.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..86beb0610f09bf537005a47277c3054216de17bc
--- /dev/null
+++ b/Tests/UnitTests/Core/Utilities/UnitVectorsTests.cpp
@@ -0,0 +1,198 @@
+// 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 <boost/test/unit_test.hpp>
+
+#include <limits>
+
+#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
+#include "Acts/Utilities/UnitVectors.hpp"
+
+using Acts::Vector3D;
+
+namespace {
+constexpr auto eps = std::numeric_limits<double>::epsilon();
+}
+
+BOOST_AUTO_TEST_SUITE(UnitVectors)
+
+BOOST_AUTO_TEST_CASE(DirectionPhiEta) {
+ using Acts::makeDirectionUnitFromPhiEta;
+
+ // along positive x
+ const auto xPos1 = makeDirectionUnitFromPhiEta(0.0, 0.0);
+ CHECK_CLOSE_REL(xPos1.norm(), 1, eps);
+ CHECK_CLOSE_REL(xPos1.dot(Vector3D(1, 0, 0)), 1, eps);
+ const auto xPos2 = makeDirectionUnitFromPhiEta(2 * M_PI, 0.0);
+ CHECK_CLOSE_REL(xPos2.norm(), 1, eps);
+ CHECK_CLOSE_REL(xPos2.dot(Vector3D(1, 0, 0)), 1, eps);
+ // along negative x
+ const auto xNeg1 = makeDirectionUnitFromPhiEta(M_PI, 0.0);
+ CHECK_CLOSE_REL(xNeg1.norm(), 1, eps);
+ CHECK_CLOSE_REL(xNeg1.dot(Vector3D(-1, 0, 0)), 1, eps);
+ const auto xNeg2 = makeDirectionUnitFromPhiEta(-M_PI, 0.0);
+ CHECK_CLOSE_REL(xNeg2.norm(), 1, eps);
+ CHECK_CLOSE_REL(xNeg2.dot(Vector3D(-1, 0, 0)), 1, eps);
+ // along positive y
+ const auto yPos1 = makeDirectionUnitFromPhiEta(M_PI_2, 0.0);
+ CHECK_CLOSE_REL(yPos1.norm(), 1, eps);
+ CHECK_CLOSE_REL(yPos1.dot(Vector3D(0, 1, 0)), 1, eps);
+ const auto yPos2 = makeDirectionUnitFromPhiEta(-3 * M_PI_2, 0.0);
+ CHECK_CLOSE_REL(yPos2.norm(), 1, eps);
+ CHECK_CLOSE_REL(yPos2.dot(Vector3D(0, 1, 0)), 1, eps);
+ // along negative y
+ const auto yNeg1 = makeDirectionUnitFromPhiEta(-M_PI_2, 0.0);
+ CHECK_CLOSE_REL(yNeg1.norm(), 1, eps);
+ CHECK_CLOSE_REL(yNeg1.dot(Vector3D(0, -1, 0)), 1, eps);
+ const auto yNeg2 = makeDirectionUnitFromPhiEta(3 * M_PI_2, 0.0);
+ CHECK_CLOSE_REL(yNeg2.norm(), 1, eps);
+ CHECK_CLOSE_REL(yNeg2.dot(Vector3D(0, -1, 0)), 1, eps);
+
+ const auto inf = std::numeric_limits<double>::infinity();
+ // along positive z
+ const auto zPos1 = makeDirectionUnitFromPhiEta(0.0, inf);
+ CHECK_CLOSE_REL(zPos1.norm(), 1, eps);
+ CHECK_CLOSE_REL(zPos1.dot(Vector3D(0, 0, 1)), 1, eps);
+ const auto zPos2 = makeDirectionUnitFromPhiEta(M_PI_2, inf);
+ CHECK_CLOSE_REL(zPos2.norm(), 1, eps);
+ CHECK_CLOSE_REL(zPos2.dot(Vector3D(0, 0, 1)), 1, eps);
+ // along negative z
+ const auto zNeg1 = makeDirectionUnitFromPhiEta(0.0, -inf);
+ CHECK_CLOSE_REL(zNeg1.norm(), 1, eps);
+ CHECK_CLOSE_REL(zNeg1.dot(Vector3D(0, 0, -1)), 1, eps);
+ const auto zNeg2 = makeDirectionUnitFromPhiEta(M_PI_2, -inf);
+ CHECK_CLOSE_REL(zNeg2.norm(), 1, eps);
+ CHECK_CLOSE_REL(zNeg2.dot(Vector3D(0, 0, -1)), 1, eps);
+
+ // mixed direction
+ const auto mixed1 = makeDirectionUnitFromPhiEta(M_PI_4, 1.0);
+ CHECK_CLOSE_REL(mixed1.norm(), 1, eps);
+ CHECK_CLOSE_REL(
+ mixed1.dot(Vector3D(1, 1, M_SQRT2 * std::sinh(1.0)).normalized()), 1,
+ eps);
+ const auto mixed2 = makeDirectionUnitFromPhiEta(M_PI_4, -1.0);
+ CHECK_CLOSE_REL(mixed2.norm(), 1, eps);
+ CHECK_CLOSE_REL(
+ mixed2.dot(Vector3D(1, 1, M_SQRT2 * std::sinh(-1.0)).normalized()), 1,
+ eps);
+ const auto mixed3 = makeDirectionUnitFromPhiEta(-M_PI_4, -1.0);
+ CHECK_CLOSE_REL(mixed3.norm(), 1, eps);
+ CHECK_CLOSE_REL(
+ mixed3.dot(Vector3D(1, -1, M_SQRT2 * std::sinh(-1.0)).normalized()), 1,
+ eps);
+}
+
+BOOST_AUTO_TEST_CASE(DirectionPhiTheta) {
+ using Acts::makeDirectionUnitFromPhiTheta;
+
+ // along positive x
+ const auto xPos1 = makeDirectionUnitFromPhiTheta(0.0, M_PI_2);
+ CHECK_CLOSE_REL(xPos1.norm(), 1, eps);
+ CHECK_CLOSE_REL(xPos1.dot(Vector3D(1, 0, 0)), 1, eps);
+ const auto xPos2 = makeDirectionUnitFromPhiTheta(2 * M_PI, M_PI_2);
+ CHECK_CLOSE_REL(xPos2.norm(), 1, eps);
+ CHECK_CLOSE_REL(xPos2.dot(Vector3D(1, 0, 0)), 1, eps);
+ // along negative x
+ const auto xNeg1 = makeDirectionUnitFromPhiTheta(M_PI, M_PI_2);
+ CHECK_CLOSE_REL(xNeg1.norm(), 1, eps);
+ CHECK_CLOSE_REL(xNeg1.dot(Vector3D(-1, 0, 0)), 1, eps);
+ const auto xNeg2 = makeDirectionUnitFromPhiTheta(-M_PI, M_PI_2);
+ CHECK_CLOSE_REL(xNeg2.norm(), 1, eps);
+ CHECK_CLOSE_REL(xNeg2.dot(Vector3D(-1, 0, 0)), 1, eps);
+ // along positive y
+ const auto yPos1 = makeDirectionUnitFromPhiTheta(M_PI_2, M_PI_2);
+ CHECK_CLOSE_REL(yPos1.norm(), 1, eps);
+ CHECK_CLOSE_REL(yPos1.dot(Vector3D(0, 1, 0)), 1, eps);
+ const auto yPos2 = makeDirectionUnitFromPhiTheta(-3 * M_PI_2, M_PI_2);
+ CHECK_CLOSE_REL(yPos2.norm(), 1, eps);
+ CHECK_CLOSE_REL(yPos2.dot(Vector3D(0, 1, 0)), 1, eps);
+ // along negative y
+ const auto yNeg1 = makeDirectionUnitFromPhiTheta(-M_PI_2, M_PI_2);
+ CHECK_CLOSE_REL(yNeg1.norm(), 1, eps);
+ CHECK_CLOSE_REL(yNeg1.dot(Vector3D(0, -1, 0)), 1, eps);
+ const auto yNeg2 = makeDirectionUnitFromPhiTheta(3 * M_PI_2, M_PI_2);
+ CHECK_CLOSE_REL(yNeg2.norm(), 1, eps);
+ CHECK_CLOSE_REL(yNeg2.dot(Vector3D(0, -1, 0)), 1, eps);
+
+ // along positive z
+ const auto zPos1 = makeDirectionUnitFromPhiTheta(0.0, 0.0);
+ CHECK_CLOSE_REL(zPos1.norm(), 1, eps);
+ CHECK_CLOSE_REL(zPos1.dot(Vector3D(0, 0, 1)), 1, eps);
+ const auto zPos2 = makeDirectionUnitFromPhiTheta(M_PI_2, 0.0);
+ CHECK_CLOSE_REL(zPos2.norm(), 1, eps);
+ CHECK_CLOSE_REL(zPos2.dot(Vector3D(0, 0, 1)), 1, eps);
+ // along negative z
+ const auto zNeg1 = makeDirectionUnitFromPhiTheta(0.0, M_PI);
+ CHECK_CLOSE_REL(zNeg1.norm(), 1, eps);
+ CHECK_CLOSE_REL(zNeg1.dot(Vector3D(0, 0, -1)), 1, eps);
+ const auto zNeg2 = makeDirectionUnitFromPhiTheta(M_PI_2, M_PI);
+ CHECK_CLOSE_REL(zNeg2.norm(), 1, eps);
+ CHECK_CLOSE_REL(zNeg2.dot(Vector3D(0, 0, -1)), 1, eps);
+
+ // mixed direction
+ const auto mixed1 = makeDirectionUnitFromPhiTheta(M_PI_4, M_PI_4);
+ CHECK_CLOSE_REL(mixed1.norm(), 1, eps);
+ CHECK_CLOSE_REL(mixed1.dot(Vector3D(1, 1, M_SQRT2).normalized()), 1, eps);
+ const auto mixed2 = makeDirectionUnitFromPhiTheta(M_PI_4, 3 * M_PI_4);
+ CHECK_CLOSE_REL(mixed2.norm(), 1, eps);
+ CHECK_CLOSE_REL(mixed2.dot(Vector3D(1, 1, -M_SQRT2).normalized()), 1, eps);
+ const auto mixed3 = makeDirectionUnitFromPhiTheta(-M_PI_4, 3 * M_PI_4);
+ CHECK_CLOSE_REL(mixed3.norm(), 1, eps);
+ CHECK_CLOSE_REL(mixed3.dot(Vector3D(1, -1, -M_SQRT2).normalized()), 1, eps);
+}
+
+namespace {
+template <typename Direction, typename RefUnitU, typename RefUnitV>
+void testCurvilinear(const Eigen::MatrixBase<Direction>& direction,
+ const Eigen::MatrixBase<RefUnitU>& refU,
+ const Eigen::MatrixBase<RefUnitV>& refV) {
+ const auto u = Acts::makeCurvilinearUnitU(direction);
+ const auto uv = Acts::makeCurvilinearUnitVectors(direction);
+ // verify normalization
+ CHECK_CLOSE_ABS(u.norm(), 1, eps);
+ CHECK_CLOSE_ABS(uv.first.norm(), 1, eps);
+ CHECK_CLOSE_ABS(uv.second.norm(), 1, eps);
+ // verify orthonormality
+ CHECK_SMALL(u.dot(direction), eps);
+ CHECK_SMALL(uv.first.dot(direction), eps);
+ CHECK_SMALL(uv.second.dot(direction), eps);
+ CHECK_SMALL(uv.first.dot(uv.second), eps);
+ // verify u is in the x-y plane
+ CHECK_SMALL(u[2], eps);
+ CHECK_SMALL(uv.first[2], eps);
+ // verify references. do not use element-wise comparison to avoid issues with
+ // small, epsilon-like differences.
+ CHECK_CLOSE_ABS(u.dot(refU), 1, eps);
+ CHECK_CLOSE_ABS(uv.first.dot(refU), 1, eps);
+ CHECK_CLOSE_ABS(uv.second.dot(refV), 1, eps);
+}
+} // namespace
+
+BOOST_AUTO_TEST_CASE(CurvilinearTransverse) {
+ // curvilinear system w/ direction in the transverse plane
+ testCurvilinear(Vector3D(1, 1, 0), Vector3D(-1, 1, 0).normalized(),
+ Vector3D(0, 0, 1).normalized());
+}
+
+BOOST_AUTO_TEST_CASE(CurvilinearPositiveZ) {
+ // curvilinear system w/ direction along z
+ testCurvilinear(Vector3D(0, 0, 1), Vector3D(1, 0, 0), Vector3D(0, 1, 0));
+}
+
+BOOST_AUTO_TEST_CASE(CurvilinearNegativeZ) {
+ // curvilinear system w/ direction along z
+ testCurvilinear(Vector3D(0, 0, -1), Vector3D(1, 0, 0), Vector3D(0, -1, 0));
+}
+
+BOOST_AUTO_TEST_CASE(CurvilinearCloseToZ) {
+ // curvilinear system w/ direction close to z
+ testCurvilinear(Vector3D(0, 32 * eps, 1 - 32 * eps), Vector3D(-1, 0, 0),
+ Vector3D(0, -1, 32 * eps));
+}
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/Tests/UnitTests/Fatras/EventData/HitTests.cpp b/Tests/UnitTests/Fatras/EventData/HitTests.cpp
index 988a5e63ddc3b2fd7c573d75f7d1fcccaa9d1a1e..f9cee9a6f918362ff631e96d5ce06cb882c0136d 100644
--- a/Tests/UnitTests/Fatras/EventData/HitTests.cpp
+++ b/Tests/UnitTests/Fatras/EventData/HitTests.cpp
@@ -30,7 +30,7 @@ BOOST_AUTO_TEST_CASE(WithoutInteraction) {
auto m4 = Hit::Vector4(1, 1, 1, 4);
auto h = Hit(gid, pid, p4, m4, m4, 12u);
- BOOST_TEST(h.geoId() == gid);
+ BOOST_TEST(h.geometryId() == gid);
BOOST_TEST(h.particleId() == pid);
BOOST_TEST(h.index() == 12u);
CHECK_CLOSE_REL(h.position4(), p4, eps);
@@ -38,9 +38,11 @@ BOOST_AUTO_TEST_CASE(WithoutInteraction) {
CHECK_CLOSE_REL(h.time(), 4, eps);
CHECK_CLOSE_REL(h.momentum4Before(), m4, eps);
CHECK_CLOSE_REL(h.momentum4After(), m4, eps);
- CHECK_CLOSE_REL(h.directionBefore(), Hit::Vector3(1, 1, 1).normalized(), eps);
- CHECK_CLOSE_REL(h.directionAfter(), Hit::Vector3(1, 1, 1).normalized(), eps);
- CHECK_CLOSE_REL(h.direction(), Hit::Vector3(1, 1, 1).normalized(), eps);
+ CHECK_CLOSE_REL(h.unitDirectionBefore(), Hit::Vector3(1, 1, 1).normalized(),
+ eps);
+ CHECK_CLOSE_REL(h.unitDirectionAfter(), Hit::Vector3(1, 1, 1).normalized(),
+ eps);
+ CHECK_CLOSE_REL(h.unitDirection(), Hit::Vector3(1, 1, 1).normalized(), eps);
CHECK_SMALL(h.depositedEnergy(), eps);
}
@@ -52,7 +54,7 @@ BOOST_AUTO_TEST_CASE(WithEnergyLoss) {
auto m41 = Hit::Vector4(1.5, 0, 0, 1.5);
auto h = Hit(gid, pid, p4, m40, m41, 13u);
- BOOST_TEST(h.geoId() == gid);
+ BOOST_TEST(h.geometryId() == gid);
BOOST_TEST(h.particleId() == pid);
BOOST_TEST(h.index() == 13u);
CHECK_CLOSE_REL(h.position4(), p4, eps);
@@ -60,9 +62,10 @@ BOOST_AUTO_TEST_CASE(WithEnergyLoss) {
CHECK_CLOSE_REL(h.time(), 4, eps);
CHECK_CLOSE_OR_SMALL(h.momentum4Before(), m40, eps, eps);
CHECK_CLOSE_OR_SMALL(h.momentum4After(), m41, eps, eps);
- CHECK_CLOSE_OR_SMALL(h.directionBefore(), Hit::Vector3(1, 0, 0), eps, eps);
- CHECK_CLOSE_OR_SMALL(h.directionAfter(), Hit::Vector3(1, 0, 0), eps, eps);
- CHECK_CLOSE_OR_SMALL(h.direction(), Hit::Vector3(1, 0, 0), eps, eps);
+ CHECK_CLOSE_OR_SMALL(h.unitDirectionBefore(), Hit::Vector3(1, 0, 0), eps,
+ eps);
+ CHECK_CLOSE_OR_SMALL(h.unitDirectionAfter(), Hit::Vector3(1, 0, 0), eps, eps);
+ CHECK_CLOSE_OR_SMALL(h.unitDirection(), Hit::Vector3(1, 0, 0), eps, eps);
CHECK_CLOSE_REL(h.depositedEnergy(), 0.5, eps);
}
@@ -74,7 +77,7 @@ BOOST_AUTO_TEST_CASE(WithScattering) {
auto m41 = Hit::Vector4(0, -2, 2, 5);
auto h = Hit(gid, pid, p4, m40, m41, 42u);
- BOOST_TEST(h.geoId() == gid);
+ BOOST_TEST(h.geometryId() == gid);
BOOST_TEST(h.particleId() == pid);
BOOST_TEST(h.index() == 42u);
CHECK_CLOSE_REL(h.position4(), p4, eps);
@@ -82,11 +85,11 @@ BOOST_AUTO_TEST_CASE(WithScattering) {
CHECK_CLOSE_REL(h.time(), 4, eps);
CHECK_CLOSE_OR_SMALL(h.momentum4Before(), m40, eps, eps);
CHECK_CLOSE_OR_SMALL(h.momentum4After(), m41, eps, eps);
- CHECK_CLOSE_OR_SMALL(h.directionBefore(), Hit::Vector3(1, 0, 1).normalized(),
- eps, eps);
- CHECK_CLOSE_OR_SMALL(h.directionAfter(), Hit::Vector3(0, -1, 1).normalized(),
- eps, eps);
- CHECK_CLOSE_REL(h.direction(), Hit::Vector3(1, -1, 2).normalized(), eps);
+ CHECK_CLOSE_OR_SMALL(h.unitDirectionBefore(),
+ Hit::Vector3(1, 0, 1).normalized(), eps, eps);
+ CHECK_CLOSE_OR_SMALL(h.unitDirectionAfter(),
+ Hit::Vector3(0, -1, 1).normalized(), eps, eps);
+ CHECK_CLOSE_REL(h.unitDirection(), Hit::Vector3(1, -1, 2).normalized(), eps);
CHECK_SMALL(h.depositedEnergy(), eps);
}
@@ -98,7 +101,7 @@ BOOST_AUTO_TEST_CASE(WithEverything) {
auto m41 = Hit::Vector4(2, 1, 2, 4);
auto h = Hit(gid, pid, p4, m40, m41, 1u);
- BOOST_TEST(h.geoId() == gid);
+ BOOST_TEST(h.geometryId() == gid);
BOOST_TEST(h.particleId() == pid);
BOOST_TEST(h.index() == 1u);
CHECK_CLOSE_REL(h.position4(), p4, eps);
@@ -106,10 +109,12 @@ BOOST_AUTO_TEST_CASE(WithEverything) {
CHECK_CLOSE_REL(h.time(), 4, eps);
CHECK_CLOSE_OR_SMALL(h.momentum4Before(), m40, eps, eps);
CHECK_CLOSE_OR_SMALL(h.momentum4After(), m41, eps, eps);
- CHECK_CLOSE_REL(h.directionBefore(), Hit::Vector3(3, 2, 2).normalized(), eps);
- CHECK_CLOSE_REL(h.directionAfter(), Hit::Vector3(2, 1, 2).normalized(), eps);
+ CHECK_CLOSE_REL(h.unitDirectionBefore(), Hit::Vector3(3, 2, 2).normalized(),
+ eps);
+ CHECK_CLOSE_REL(h.unitDirectionAfter(), Hit::Vector3(2, 1, 2).normalized(),
+ eps);
CHECK_CLOSE_REL(
- h.direction(),
+ h.unitDirection(),
Hit::Vector3(0.7023994590205035, 0.41229136135810396, 0.5802161953247991),
eps);
CHECK_CLOSE_REL(h.depositedEnergy(), 1, eps);
diff --git a/Tests/UnitTests/Fatras/EventData/ParticleTests.cpp b/Tests/UnitTests/Fatras/EventData/ParticleTests.cpp
index 28e03564e3a9f1019c166b588d7447f9ed852c95..4ba1b078c260256a1ffb5f1a41b69e7bb23a2cf0 100644
--- a/Tests/UnitTests/Fatras/EventData/ParticleTests.cpp
+++ b/Tests/UnitTests/Fatras/EventData/ParticleTests.cpp
@@ -8,6 +8,9 @@
#include <boost/test/unit_test.hpp>
+#include <limits>
+
+#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
#include "Acts/Utilities/Units.hpp"
#include "ActsFatras/EventData/Particle.hpp"
@@ -16,13 +19,17 @@ using ActsFatras::Barcode;
using ActsFatras::Particle;
using namespace Acts::UnitLiterals;
+namespace {
+constexpr auto eps = std::numeric_limits<Particle::Scalar>::epsilon();
+}
+
BOOST_AUTO_TEST_SUITE(FatrasParticle)
BOOST_AUTO_TEST_CASE(Construct) {
- const auto id = Barcode().setVertexPrimary(1).setParticle(42);
- const auto particle = Particle(id, PdgParticle::eProton, 1_GeV, 1_e);
+ const auto pid = Barcode().setVertexPrimary(1).setParticle(42);
+ const auto particle = Particle(pid, PdgParticle::eProton, 1_GeV, 1_e);
- BOOST_TEST(particle.id() == id);
+ BOOST_TEST(particle.particleId() == pid);
BOOST_TEST(particle.pdg() == PdgParticle::eProton);
// particle is at rest at the origin
BOOST_TEST(particle.position4() == Particle::Vector4::Zero());
@@ -32,17 +39,19 @@ BOOST_AUTO_TEST_CASE(Construct) {
BOOST_TEST(particle.position4().y() == particle.position().y());
BOOST_TEST(particle.position4().z() == particle.position().z());
BOOST_TEST(particle.position4().w() == particle.time());
- // particle direction is undefined
- BOOST_TEST(particle.momentum() == Particle::Scalar(0));
+ // particle direction is undefined, but must be normalized
+ CHECK_CLOSE_REL(particle.unitDirection().norm(), 1, eps);
+ BOOST_TEST(particle.transverseMomentum() == Particle::Scalar(0));
+ BOOST_TEST(particle.absMomentum() == Particle::Scalar(0));
// particle is created at rest and thus not alive
BOOST_TEST(not particle);
}
BOOST_AUTO_TEST_CASE(CorrectEnergy) {
- const auto id = Barcode().setVertexPrimary(1).setParticle(42);
- auto particle = Particle(id, PdgParticle::eProton, 1_GeV, 1_e)
+ const auto pid = Barcode().setVertexPrimary(1).setParticle(42);
+ auto particle = Particle(pid, PdgParticle::eProton, 1_GeV, 1_e)
.setDirection(Particle::Vector3::UnitX())
- .setMomentum(2_GeV);
+ .setAbsMomentum(2_GeV);
BOOST_TEST(particle.mass() == 1_GeV);
// check that the particle has some input energy
@@ -50,23 +59,32 @@ BOOST_AUTO_TEST_CASE(CorrectEnergy) {
BOOST_TEST(particle.momentum4().y() == 0_GeV);
BOOST_TEST(particle.momentum4().z() == 0_GeV);
BOOST_TEST(particle.momentum4().w() == std::hypot(1_GeV, 2_GeV));
- BOOST_TEST(particle.momentum() == 2_GeV);
+ BOOST_TEST(particle.transverseMomentum() == 2_GeV);
+ BOOST_TEST(particle.absMomentum() == 2_GeV);
BOOST_TEST(particle.energy() == std::hypot(1_GeV, 2_GeV));
+ // particle direction must be normalized
+ CHECK_CLOSE_REL(particle.unitDirection().norm(), 1, eps);
// loose some energy
particle.correctEnergy(-100_MeV);
- BOOST_TEST(particle.momentum() < 2_GeV);
+ BOOST_TEST(particle.transverseMomentum() < 2_GeV);
+ BOOST_TEST(particle.absMomentum() < 2_GeV);
BOOST_TEST(particle.energy() ==
Particle::Scalar(std::hypot(1_GeV, 2_GeV) - 100_MeV));
+ CHECK_CLOSE_REL(particle.unitDirection().norm(), 1, eps);
// particle is still alive
BOOST_TEST(particle);
// loose a lot of energy
particle.correctEnergy(-3_GeV);
- BOOST_TEST(particle.momentum() == Particle::Scalar(0));
+ BOOST_TEST(particle.transverseMomentum() == Particle::Scalar(0));
+ BOOST_TEST(particle.absMomentum() == Particle::Scalar(0));
BOOST_TEST(particle.energy() == particle.mass());
+ CHECK_CLOSE_REL(particle.unitDirection().norm(), 1, eps);
// lossing even more energy does nothing
particle.correctEnergy(-10_GeV);
- BOOST_TEST(particle.momentum() == Particle::Scalar(0));
+ BOOST_TEST(particle.transverseMomentum() == Particle::Scalar(0));
+ BOOST_TEST(particle.absMomentum() == Particle::Scalar(0));
BOOST_TEST(particle.energy() == particle.mass());
+ CHECK_CLOSE_REL(particle.unitDirection().norm(), 1, eps);
// particle is not alive anymore
BOOST_TEST(not particle);
}
diff --git a/Tests/UnitTests/Fatras/Kernel/InteractorTests.cpp b/Tests/UnitTests/Fatras/Kernel/InteractorTests.cpp
index cdc1ab664e08679bda965996fe7c69fea02a596b..dd84e572f806237dd1f5582c0ad912e37dbb583c 100644
--- a/Tests/UnitTests/Fatras/Kernel/InteractorTests.cpp
+++ b/Tests/UnitTests/Fatras/Kernel/InteractorTests.cpp
@@ -94,15 +94,15 @@ struct Fixture {
Acts::PdgParticle::eProton, 0, 1)
.setPosition4(1, 2, 3, 4)
.setDirection(1, 0, 0)
- .setMomentum(100);
+ .setAbsMomentum(100);
interactor.generator = &generator;
interactor.physics.energyLoss = energyLoss;
interactor.particle = particle;
state.navigation.currentSurface = surface.get();
state.stepping.position = particle.position();
state.stepping.time = particle.time();
- state.stepping.direction = particle.direction();
- state.stepping.momentum = particle.momentum();
+ state.stepping.direction = particle.unitDirection();
+ state.stepping.momentum = particle.absMomentum();
}
};
@@ -142,7 +142,8 @@ BOOST_AUTO_TEST_CASE(HitsOnEmptySurface) {
BOOST_TEST(f.result.hits[1].index() == 1u);
// particle identity should be the same as the initial input
- BOOST_TEST(f.result.particle.id() == f.interactor.particle.id());
+ BOOST_TEST(f.result.particle.particleId() ==
+ f.interactor.particle.particleId());
BOOST_TEST(f.result.particle.process() == f.interactor.particle.process());
BOOST_TEST(f.result.particle.pdg() == f.interactor.particle.pdg());
BOOST_TEST(f.result.particle.charge() == f.interactor.particle.charge());
@@ -169,7 +170,8 @@ BOOST_AUTO_TEST_CASE(HitsOnMaterialSurface) {
BOOST_TEST(f.result.hits[1].index() == 1u);
// particle identity should be the same as the initial input
- BOOST_TEST(f.result.particle.id() == f.interactor.particle.id());
+ BOOST_TEST(f.result.particle.particleId() ==
+ f.interactor.particle.particleId());
BOOST_TEST(f.result.particle.process() == f.interactor.particle.process());
BOOST_TEST(f.result.particle.pdg() == f.interactor.particle.pdg());
BOOST_TEST(f.result.particle.charge() == f.interactor.particle.charge());
@@ -193,7 +195,8 @@ BOOST_AUTO_TEST_CASE(NoHitsEmptySurface) {
BOOST_TEST(f.result.hits.size() == 0u);
// particle identity should be the same as the initial input
- BOOST_TEST(f.result.particle.id() == f.interactor.particle.id());
+ BOOST_TEST(f.result.particle.particleId() ==
+ f.interactor.particle.particleId());
BOOST_TEST(f.result.particle.process() == f.interactor.particle.process());
BOOST_TEST(f.result.particle.pdg() == f.interactor.particle.pdg());
BOOST_TEST(f.result.particle.charge() == f.interactor.particle.charge());
@@ -217,7 +220,8 @@ BOOST_AUTO_TEST_CASE(NoHitsMaterialSurface) {
BOOST_TEST(f.result.hits.size() == 0u);
// particle identity should be the same as the initial input
- BOOST_TEST(f.result.particle.id() == f.interactor.particle.id());
+ BOOST_TEST(f.result.particle.particleId() ==
+ f.interactor.particle.particleId());
BOOST_TEST(f.result.particle.process() == f.interactor.particle.process());
BOOST_TEST(f.result.particle.pdg() == f.interactor.particle.pdg());
BOOST_TEST(f.result.particle.charge() == f.interactor.particle.charge());
diff --git a/Tests/UnitTests/Fatras/Kernel/ProcessTests.cpp b/Tests/UnitTests/Fatras/Kernel/ProcessTests.cpp
index 00d281ae683fbb609b6d30127e3c44493caf9c1a..7b4fe2a8c7e56ae061e6689214d0bbb63311be8f 100644
--- a/Tests/UnitTests/Fatras/Kernel/ProcessTests.cpp
+++ b/Tests/UnitTests/Fatras/Kernel/ProcessTests.cpp
@@ -30,10 +30,10 @@ struct MakeChildren {
ActsFatras::Particle &) const {
// create daughter particles
return {
- Particle().setMomentum(1_GeV),
- Particle().setMomentum(2_GeV),
- Particle().setMomentum(3_GeV),
- Particle().setMomentum(4_GeV),
+ Particle().setAbsMomentum(1_GeV),
+ Particle().setAbsMomentum(2_GeV),
+ Particle().setAbsMomentum(3_GeV),
+ Particle().setAbsMomentum(4_GeV),
};
}
};
@@ -43,14 +43,14 @@ struct HighP {
double minP = 10_GeV;
bool operator()(const ActsFatras::Particle &particle) const {
- return (minP <= particle.momentum());
+ return (minP <= particle.absMomentum());
}
};
struct Fixture {
std::default_random_engine generator;
Acts::MaterialProperties slab{Acts::Test::makeBeryllium(), 1_mm};
- Particle parent = Particle().setMomentum(10_GeV);
+ Particle parent = Particle().setAbsMomentum(10_GeV);
std::vector<Particle> children;
};
} // namespace
@@ -72,15 +72,15 @@ BOOST_AUTO_TEST_CASE(WithInputSelector) {
process.selectInput.minP = 10_GeV;
// above threshold should not abort
- f.parent.setMomentum(20_GeV);
+ f.parent.setAbsMomentum(20_GeV);
BOOST_TEST(not process(f.generator, f.slab, f.parent, f.children));
BOOST_TEST(f.children.size() == 4u);
// on threshold should still not abort
- f.parent.setMomentum(10_GeV);
+ f.parent.setAbsMomentum(10_GeV);
BOOST_TEST(not process(f.generator, f.slab, f.parent, f.children));
BOOST_TEST(f.children.size() == 8u);
// below threshold should abort and not run the process at all
- f.parent.setMomentum(2_GeV);
+ f.parent.setAbsMomentum(2_GeV);
BOOST_TEST(not process(f.generator, f.slab, f.parent, f.children));
// process did not run -> no new children
BOOST_TEST(f.children.size() == 8u);
@@ -92,15 +92,15 @@ BOOST_AUTO_TEST_CASE(WithOutputSelector) {
process.selectOutputParticle.minP = 10_GeV;
// above threshold should not abort
- f.parent.setMomentum(20_GeV);
+ f.parent.setAbsMomentum(20_GeV);
BOOST_TEST(not process(f.generator, f.slab, f.parent, f.children));
BOOST_TEST(f.children.size() == 4u);
// on threshold should still not abort
- f.parent.setMomentum(10_GeV);
+ f.parent.setAbsMomentum(10_GeV);
BOOST_TEST(not process(f.generator, f.slab, f.parent, f.children));
BOOST_TEST(f.children.size() == 8u);
// below threshold should abort but only after running the process
- f.parent.setMomentum(2_GeV);
+ f.parent.setAbsMomentum(2_GeV);
BOOST_TEST(process(f.generator, f.slab, f.parent, f.children));
// process did still run -> new children
BOOST_TEST(f.children.size() == 12u);
diff --git a/Tests/UnitTests/Fatras/Physics/Dataset.hpp b/Tests/UnitTests/Fatras/Physics/Dataset.hpp
index 63f4a072a5de01006d390aa0e56fe21c590b438b..eee0c82117753522cfbd3d9ade9a182fb58102bd 100644
--- a/Tests/UnitTests/Fatras/Physics/Dataset.hpp
+++ b/Tests/UnitTests/Fatras/Physics/Dataset.hpp
@@ -50,7 +50,7 @@ inline ActsFatras::Particle makeParticle(Acts::PdgParticle pdg, double phi,
.setPosition4(0, 0, 0, 0)
.setDirection(std::cos(lambda) * std::cos(phi),
std::cos(lambda) * std::sin(phi), std::sin(lambda))
- .setMomentum(p);
+ .setAbsMomentum(p);
}
} // namespace Dataset
diff --git a/Tests/UnitTests/Fatras/Physics/EnergyLossTests.cpp b/Tests/UnitTests/Fatras/Physics/EnergyLossTests.cpp
index ea73418011077b7797a718445f7578caa306db53..8d77fcac83b59b135a41e77da4afd513f471ad22 100644
--- a/Tests/UnitTests/Fatras/Physics/EnergyLossTests.cpp
+++ b/Tests/UnitTests/Fatras/Physics/EnergyLossTests.cpp
@@ -31,7 +31,7 @@ BOOST_DATA_TEST_CASE(BetheBloch, Dataset::parameters, pdg, phi, lambda, p,
ActsFatras::BetheBloch process;
const auto outgoing = process(gen, Acts::Test::makeUnitSlab(), after);
// energy loss changes momentum and energy
- BOOST_TEST(after.momentum() < before.momentum());
+ BOOST_TEST(after.absMomentum() < before.absMomentum());
BOOST_TEST(after.energy() < before.energy());
// energy loss creates no new particles
BOOST_TEST(outgoing.empty());
@@ -46,7 +46,7 @@ BOOST_DATA_TEST_CASE(BetheHeitler, Dataset::parameters, pdg, phi, lambda, p,
ActsFatras::BetheHeitler process;
const auto outgoing = process(gen, Acts::Test::makeUnitSlab(), after);
// energy loss changes momentum and energy
- BOOST_TEST(after.momentum() < before.momentum());
+ BOOST_TEST(after.absMomentum() < before.absMomentum());
BOOST_TEST(after.energy() < before.energy());
// energy loss creates no new particles
BOOST_TEST(outgoing.empty());
diff --git a/Tests/UnitTests/Fatras/Physics/ScatteringTests.cpp b/Tests/UnitTests/Fatras/Physics/ScatteringTests.cpp
index 159849f1a95ca3d8c060598bc3e20e19a4bc15a5..7a2cf875428c849ede728110c4622fa45f757aea 100644
--- a/Tests/UnitTests/Fatras/Physics/ScatteringTests.cpp
+++ b/Tests/UnitTests/Fatras/Physics/ScatteringTests.cpp
@@ -19,18 +19,11 @@
#include "ActsFatras/Physics/Scattering/GaussianMixture.hpp"
#include "ActsFatras/Physics/Scattering/GeneralMixture.hpp"
#include "ActsFatras/Physics/Scattering/Highland.hpp"
-#include "ActsFatras/Physics/Scattering/Scattering.hpp"
#include "Dataset.hpp"
namespace {
constexpr auto eps = std::numeric_limits<double>::epsilon();
-using GeneralMixtureScattering =
- ActsFatras::Scattering<ActsFatras::GeneralMixture>;
-using GaussianMixtureScattering =
- ActsFatras::Scattering<ActsFatras::GaussianMixture>;
-using HighlandScattering = ActsFatras::Scattering<ActsFatras::Highland>;
-
// Common test method that will be instantiated for each scattering model.
template <typename Scattering>
void test(const Scattering& scattering, uint32_t seed,
@@ -40,10 +33,10 @@ void test(const Scattering& scattering, uint32_t seed,
const auto outgoing = scattering(gen, Acts::Test::makePercentSlab(), after);
// scattering leaves absolute energy/momentum unchanged
- CHECK_CLOSE_REL(after.momentum(), before.momentum(), eps);
+ CHECK_CLOSE_REL(after.absMomentum(), before.absMomentum(), eps);
CHECK_CLOSE_REL(after.energy(), before.energy(), eps);
// scattering has changed the direction
- BOOST_TEST(before.direction().dot(after.direction()) < 1);
+ BOOST_TEST(before.unitDirection().dot(after.unitDirection()) < 1);
// scattering creates no new particles
BOOST_TEST(outgoing.empty());
}
@@ -53,18 +46,19 @@ BOOST_AUTO_TEST_SUITE(FatrasScattering)
BOOST_DATA_TEST_CASE(GeneralMixture, Dataset::parameters, pdg, phi, lambda, p,
seed) {
- test(GeneralMixtureScattering(), seed,
+ test(ActsFatras::GeneralMixtureScattering(), seed,
Dataset::makeParticle(pdg, phi, lambda, p));
}
BOOST_DATA_TEST_CASE(GaussianMixture, Dataset::parameters, pdg, phi, lambda, p,
seed) {
- test(GaussianMixtureScattering(), seed,
+ test(ActsFatras::GaussianMixtureScattering(), seed,
Dataset::makeParticle(pdg, phi, lambda, p));
}
BOOST_DATA_TEST_CASE(Highland, Dataset::parameters, pdg, phi, lambda, p, seed) {
- test(HighlandScattering(), seed, Dataset::makeParticle(pdg, phi, lambda, p));
+ test(ActsFatras::HighlandScattering(), seed,
+ Dataset::makeParticle(pdg, phi, lambda, p));
}
BOOST_AUTO_TEST_SUITE_END()
diff --git a/Tests/UnitTests/Fatras/Selectors/Dataset.hpp b/Tests/UnitTests/Fatras/Selectors/Dataset.hpp
index ddcca97a985ef094bea7fdec2f6b5405c1425eba..f725656d44bdfb47a2b78d3ffee9b10a756d76d6 100644
--- a/Tests/UnitTests/Fatras/Selectors/Dataset.hpp
+++ b/Tests/UnitTests/Fatras/Selectors/Dataset.hpp
@@ -21,7 +21,7 @@ ActsFatras::Particle makeParticle(Acts::PdgParticle pdg, double z, double eta) {
return ActsFatras::Particle(id, pdg)
.setPosition4(0.0, 0.0, z, 0.0)
.setDirection(1.0 / std::cosh(eta), 0.0, std::tanh(eta))
- .setMomentum(1.5_GeV);
+ .setAbsMomentum(1.5_GeV);
}
const auto centralElectron =