diff --git a/Core/include/Acts/Utilities/PdgParticle.hpp b/Core/include/Acts/Utilities/PdgParticle.hpp
index c8da3ec9f96fd8eb97554ba867ac78d86c97f6f7..3ad4722ca26f91a1e38f8316be049eec5587058a 100644
--- a/Core/include/Acts/Utilities/PdgParticle.hpp
+++ b/Core/include/Acts/Utilities/PdgParticle.hpp
@@ -8,10 +8,13 @@
#pragma once
+#include <cstdint>
+
namespace Acts {
/// Symbolic values for commonly used PDG particle numbers.
-enum PdgParticle : int {
+enum PdgParticle : int32_t {
+ eInvalid = 0,
eElectron = 11,
eAntiElectron = -eElectron,
ePositron = -eElectron,
diff --git a/Fatras/include/ActsFatras/EventData/Barcode.hpp b/Fatras/include/ActsFatras/EventData/Barcode.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..f19e11be3e52901fd499d545537e5c0b0790590a
--- /dev/null
+++ b/Fatras/include/ActsFatras/EventData/Barcode.hpp
@@ -0,0 +1,126 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2018-2020 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include <cstdint>
+
+#include "Acts/Utilities/MultiIndex.hpp"
+
+namespace ActsFatras {
+
+/// Particle identifier that encodes additional event information.
+///
+/// The barcode has to fulfill two separate requirements: be able to act as
+/// unique identifier for particles within an event and to encode details
+/// on the event structure for fast lookup. Since we only care about tracking
+/// here, we need to support two scenarios:
+///
+/// * Identify which primary/secondary vertex particles belong to. No
+/// information on intermediate/unstable/invisible particles needs to be
+/// retained. This information is already available in the underlying
+/// generator event and should not be duplicated.
+/// * If (visible) particles convert, decay, or interact with the detector,
+/// we need to be able to identify the initial (primary) particle. Typical
+/// examples are pion nuclear interactions or electron/gamma conversions.
+///
+/// The vertex information is encoded as two numbers that define the
+/// primary and secondary vertex. The primary vertex must be non-zero.
+/// Particles with a zero secondary vertex originate directly from the primary
+/// vertex.
+///
+/// Within one vertex (primary+secondary) each particle is identified by
+/// a particle, generation, and sub-particle number. Particles originating
+/// from the vertex must have zero generation and zero sub-particle number;
+/// 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.
+///
+/// With this encoding, non-primary particles and their primary parent can
+/// be easily identified at the expense of not storing the exact decay history.
+///
+/// A barcode with all elements set to zero (the default value) is an invalid
+/// value that can be used e.g. to mark missing or unknown particles.
+///
+/// ## Example
+///
+/// A particle generated in a primary interaction might have the barcode
+///
+/// 2|0|14|0|0 -> vertex=2 (primary), particle=14, generation=0, sub=0
+///
+/// A simulation module might generate an interaction and create two new
+/// particles. These are descendants of the initial particle and the simulation
+/// module can generate the new barcodes directly by increasing the
+/// generation number and chosing sub-particle identifiers:
+///
+/// 2|0|14|1|0 -> vertex=2 (primary), particle=14, generation=1, sub=0
+/// 2|0|14|1|1 -> vertex=2 (primary), particle=14, generation=1, sub=1
+///
+/// If these secondary particles generate further tertiary particles
+/// the barcode would be e.g.
+///
+/// 2|0|14|2|0 -> vertex=2 (primary), particle=14, generation=2, sub=0
+///
+/// ## Possible issues
+///
+/// The hierachical nature of the barcode allows barcode creation without
+/// a central service. Since the full history is not stored, generated barcodes
+/// for higher-generation particles can overlap when generated by independent
+/// interactions. Assuming an initial primary particle with barcode
+///
+/// 3|4|5|0|0 -> particle=5
+///
+/// a first interaction might create a secondary particle by increasing the
+/// generation number (without destroying the initial particle)
+///
+/// 3|4|5|1|0 -> particle=5, generation+=1, first sub-particle
+///
+/// The initial particle gets simulated further and at another step a second
+/// interaction also creates a new particle. Since it knows nothing about
+/// the previously created particle (no central service), it will generate
+///
+/// 3|4|5|1|0 -> particle=5, generation+=1, first sub-particle
+///
+/// which is identical to the previously create barcode. These cases can be
+/// easily solved by renumbering the sub-particle identifier within each
+/// 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> {
+ public:
+ /// Return the primary vertex identifier.
+ constexpr Value vertexPrimary() const { return level(0); }
+ /// Return the secondary vertex identifier.
+ constexpr Value vertexSecondary() const { return level(1); }
+ /// Return the particle identifier.
+ constexpr Value particle() const { return level(2); }
+ /// Return the generation identifier.
+ constexpr Value generation() const { return level(3); }
+ /// Return the sub-particle identifier.
+ constexpr Value subParticle() const { return level(4); }
+
+ /// Set the primary vertex identifier.
+ constexpr Barcode& setVertexPrimary(Value id) { return set(0, id), *this; }
+ /// Set the secondary vertex identifier.
+ constexpr Barcode& setVertexSecondary(Value id) { return set(1, id), *this; }
+ /// Set the parent particle identifier.
+ constexpr Barcode& setParticle(Value id) { return set(2, id), *this; }
+ /// Set the particle identifier.
+ constexpr Barcode& setGeneration(Value id) { return set(3, id), *this; }
+ /// Set the process identifier.
+ constexpr Barcode& setSubParticle(Value id) { return set(4, id), *this; }
+
+ /// Construct a new barcode representing a descendant particle.
+ ///
+ /// @param sub sub-particle index of the new barcode.
+ Barcode makeDescendant(Value sub = 0u) const {
+ return Barcode(*this).setGeneration(generation() + 1).setSubParticle(sub);
+ }
+};
+
+} // namespace ActsFatras
diff --git a/Fatras/include/ActsFatras/EventData/Particle.hpp b/Fatras/include/ActsFatras/EventData/Particle.hpp
index d5a08f2ef737e358e02e80c1e0b9c463dee99c87..0d207fec9ee404d36a45176a2f54d3bb96fc5810 100644
--- a/Fatras/include/ActsFatras/EventData/Particle.hpp
+++ b/Fatras/include/ActsFatras/EventData/Particle.hpp
@@ -9,206 +9,202 @@
#pragma once
#include <cmath>
+#include <cstdint>
+#include <limits>
-#include "Acts/Geometry/GeometryID.hpp"
-#include "Acts/Surfaces/Surface.hpp"
#include "Acts/Utilities/Definitions.hpp"
-#include "Acts/Utilities/Helpers.hpp"
-#include "Acts/Utilities/Units.hpp"
+#include "Acts/Utilities/PdgParticle.hpp"
+#include "ActsFatras/EventData/Barcode.hpp"
namespace ActsFatras {
-/// Typedef the pdg code
-typedef int pdg_type;
-
-/// Typedef the process code
-typedef unsigned int process_code;
-
-/// Typedef barcode
-typedef unsigned int barcode_type;
+/// Process type identifier.
+///
+/// Encodes the type of process that generated the particle.
+enum class ProcessType : uint32_t {
+ eUndefined = 0,
+};
/// Simulation particle information and kinematic state.
class Particle {
public:
- /// Default
- Particle() = default;
+ using Scalar = double;
+ using Vector3 = Acts::ActsVector<Scalar, 3>;
+ using Vector4 = Acts::ActsVector<Scalar, 4>;
- /// @brief Construct a particle consistently
+ /// Construct a default particle with invalid identity.
+ Particle() = default;
+ /// Construct a particle at rest with explicit mass and charge.
///
- /// @param pposition The particle position at construction
- /// @param pmomentum The particle momentum at construction
- /// @param pm The particle mass
- /// @param pq The partilce charge
- /// @param pbarcode The particle barcode
- Particle(const Acts::Vector3D &position, const Acts::Vector3D &momentum,
- double m, double q, pdg_type pdg = 0, barcode_type barcode = 0,
- double startTime = 0.)
- : m_position(position),
- m_momentum(momentum),
- m_m(m),
- m_q(q),
- m_p(momentum.norm()),
- m_pT(Acts::VectorHelpers::perp(momentum)),
- m_pdg(pdg),
- m_barcode(barcode),
- m_timeStamp(startTime) {
- m_E = std::sqrt(m_p * m_p + m_m * m_m);
- m_beta = (m_p / m_E);
- m_gamma = (m_E / m_m);
- }
-
- /// @brief Set the limits
+ /// @param id Encoded identifier within an event
+ /// @param pdg PDG particle number
+ /// @param mass Particle mass in native units
+ /// @param charge Particle charge in native units
///
- /// @param x0Limit the limit in X0 to be passed
- /// @param l0Limit the limit in L0 to be passed
- /// @param timeLimit the readout time limit to be passed
- void setLimits(double x0Limit, double l0Limit,
- double timeLimit = std::numeric_limits<double>::max()) {
- m_limitInX0 = x0Limit;
- m_limitInL0 = l0Limit;
- m_timeLimit = timeLimit;
+ /// @warning It is the users responsibility that charge and mass match
+ /// the PDG particle number.
+ Particle(Barcode id, Acts::PdgParticle pdg, Scalar mass, Scalar charge)
+ : m_id(id), m_pdg(pdg), m_charge(charge), m_mass(mass) {}
+ Particle(const Particle &) = default;
+ Particle(Particle &&) = default;
+ Particle &operator=(const Particle &) = default;
+ Particle &operator=(Particle &&) = default;
+
+ /// Set the process type that generated this particle.
+ Particle &process(ProcessType proc) { return m_process = proc, *this; }
+ /// Set the space-time position four-vector.
+ Particle &setPosition4(const Vector4 &pos4) {
+ m_position4 = pos4;
+ return *this;
}
-
- /// @brief Update the particle with applying energy loss
- ///
- /// @param deltaE is the energy loss to be applied
- void scatter(Acts::Vector3D nmomentum) {
- m_momentum = std::move(nmomentum);
- m_pT = Acts::VectorHelpers::perp(m_momentum);
+ /// Set the space-time position four-vector from three-position and time.
+ Particle &setPosition4(const Vector3 &position, Scalar time) {
+ m_position4.head<3>() = position;
+ m_position4[3] = time;
+ return *this;
}
-
- /// @brief Update the particle with applying energy loss
+ /// Set the space-time position four-vector from scalar components.
+ Particle &setPosition4(Scalar x, Scalar y, Scalar z, Scalar time) {
+ m_position4[0] = x;
+ m_position4[1] = y;
+ m_position4[2] = z;
+ m_position4[3] = time;
+ return *this;
+ }
+ /// Set the direction three-vector
+ Particle &setDirection(const Vector3 &direction) {
+ m_direction = direction;
+ m_direction.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();
+ return *this;
+ }
+ /// Set the absolute momentum.
+ Particle &setMomentum(Scalar momentum) {
+ m_momentum = momentum;
+ return *this;
+ }
+ /// Change the energy by the given amount.
///
- /// @param deltaE is the energy loss to be applied
- void energyLoss(double deltaE) {
- // particle falls to rest
- if (m_E - deltaE < m_m) {
- m_E = m_m;
- m_p = 0.;
- m_pT = 0.;
- m_beta = 0.;
- m_gamma = 1.;
- m_momentum = Acts::Vector3D(0., 0., 0.);
- m_alive = false;
+ /// Energy loss corresponds to a negative change. If the updated energy
+ /// 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;
+ if (newEnergy <= m_mass) {
+ m_momentum = Scalar(0);
+ } else {
+ m_momentum = std::sqrt(newEnergy * newEnergy - m_mass * m_mass);
}
- // updatet the parameters
- m_E -= deltaE;
- m_p = std::sqrt(m_E * m_E - m_m * m_m);
- m_momentum = m_p * m_momentum.normalized();
- m_pT = Acts::VectorHelpers::perp(m_momentum);
- m_beta = (m_p / m_E);
- m_gamma = (m_E / m_m);
+ return *this;
}
- /// @brief Update the particle with a new position and momentum,
- /// this corresponds to a step update
+ /// Encoded particle identifier within an event.
+ Barcode id() const { return m_id; }
+ /// Which type of process generated this particle.
+ ProcessType process() const { return m_process; }
+ /// PDG particle number that identifies the type.
+ Acts::PdgParticle pdg() const { return m_pdg; }
+ /// Particle charge.
+ Scalar charge() const { return m_charge; }
+ /// Particle mass.
+ Scalar mass() const { return m_mass; }
+
+ /// Space-time position four-vector.
+ 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.
+ Scalar time() const { return m_position4[3]; }
+ /// Energy-momentum four-vector.
+ 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[3] = energy();
+ return mom4;
+ }
+ /// Three-direction, i.e. the normalized momentum three-vector.
+ const Vector3 &direction() const { return m_direction; }
+ /// Absolute momentum.
+ Scalar momentum() const { return m_momentum; }
+ /// 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); }
+
+ /// Check if the particle is alive, i.e. is not at rest.
+ operator bool() const { return Scalar(0) < m_momentum; }
+ /// Check if the particle is dead, i.e is at rest.
+ bool operator!() const { return m_momentum <= Scalar(0); }
+
+ /// Register material that the particle has passed.
///
- /// @param position New position after update
- /// @param momentum New momentum after update
- /// @param deltaPathX0 passed since last step
- /// @param deltaPathL0 passed since last step
- /// @param deltaTime The time elapsed
+ /// @param thicknessX0 material thickness measured in radiation lengths
+ /// @param thicknessL0 material thickness measured in interaction lengths
+ Particle &addPassedMaterial(Scalar thicknessX0, Scalar thicknessL0) {
+ m_pathX0 += thicknessX0;
+ m_pathL0 += thicknessL0;
+ return *this;
+ }
+ /// Set the material limits.
///
- /// @return break condition
- bool update(const Acts::Vector3D &position, const Acts::Vector3D &momentum,
- double deltaPahtX0 = 0., double deltaPahtL0 = 0.,
- double deltaTime = 0.) {
- m_position = position;
- m_momentum = momentum;
- m_p = momentum.norm();
- if (m_p) {
- m_pT = Acts::VectorHelpers::perp(momentum);
- m_E = std::sqrt(m_p * m_p + m_m * m_m);
- m_timeStamp += deltaTime;
- m_beta = (m_p / m_E);
- m_gamma = (m_E / m_m);
-
- // set parameters and check limits
- m_pathInX0 += deltaPahtX0;
- m_pathInL0 += deltaPahtL0;
- m_timeStamp += deltaTime;
- if (m_pathInX0 >= m_limitInX0 || m_pathInL0 >= m_limitInL0 ||
- m_timeStamp > m_timeLimit) {
- m_alive = false;
- }
- }
- return !m_alive;
+ /// @param limitX0 maximum radiation lengths the particle can pass
+ /// @param limitL0 maximum interaction lengths the particle can pass
+ Particle &setMaterialLimits(Scalar limitX0, Scalar limitL0) {
+ m_limitX0 = limitX0;
+ m_limitL0 = limitL0;
+ return *this;
}
-
- /// @brief Access methods: position
- const Acts::Vector3D &position() const { return m_position; }
-
- /// @brief Access methods: momentum
- const Acts::Vector3D &momentum() const { return m_momentum; }
-
- /// @brief Access methods: p
- double p() const { return m_p; }
-
- /// @brief Access methods: pT
- double pT() const { return m_pT; }
-
- /// @brief Access methods: E
- double E() const { return m_E; }
-
- /// @brief Access methods: m
- double m() const { return m_m; }
-
- /// @brief Access methods: beta
- double beta() const { return m_beta; }
-
- /// @brief Access methods: gamma
- double gamma() const { return m_gamma; }
-
- /// @brief Access methods: charge
- double q() const { return m_q; }
-
- /// @brief Access methods: pdg code
- pdg_type pdg() const { return m_pdg; }
-
- /// @brief Access methods: barcode
- barcode_type barcode() const { return m_barcode; }
-
- /// @brief Access methods: path/X0
- double pathInX0() const { return m_pathInX0; }
-
- /// @brief Access methods: limit/X0
- double limitInX0() const { return m_limitInX0; }
-
- /// @brief Access methods: pdg code
- double pathInL0() const { return m_limitInX0; }
-
- /// @brief Access methods: barcode
- double limitInL0() const { return m_limitInL0; }
-
- /// @brief boolean operator indicating the particle to be alive
- operator bool() { return m_alive; }
+ /// The passed material measured in radiation lengths.
+ Scalar pathInX0() const { return m_pathX0; }
+ /// The passed material measured in interaction lengths.
+ Scalar pathInL0() const { return m_pathL0; }
+ /// The maximum radation length the particle is allowed to pass.
+ Scalar pathLimitX0() const { return m_limitX0; }
+ /// The maximum interaction length the particle is allowed to pass.
+ Scalar pathLimitL0() const { return m_limitL0; }
private:
- Acts::Vector3D m_position = Acts::Vector3D(0., 0., 0.); //!< kinematic info
- Acts::Vector3D m_momentum = Acts::Vector3D(0., 0., 0.); //!< kinematic info
-
- double m_m = 0.; //!< particle mass
- double m_E = 0.; //!< total energy
- double m_q = 0.; //!< the charge
- double m_beta = 0.; //!< relativistic beta factor
- double m_gamma = 1.; //!< relativistic gamma factor
- double m_p = 0.; //!< momentum magnitude
- double m_pT = 0.; //!< transverse momentum magnitude
- pdg_type m_pdg = 0; //!< pdg code of the particle
- barcode_type m_barcode = 0; //!< barcode of the particle
-
- double m_pathInX0 = 0.; //!< passed path in X0
- double m_limitInX0 =
- std::numeric_limits<double>::max(); //!< path limit in X0
-
- double m_pathInL0 = 0.; //!< passed path in L0
- double m_limitInL0 =
- std::numeric_limits<double>::max(); //!< path limit in X0
-
- double m_timeStamp = 0.; //!< passed time elapsed
- double m_timeLimit = std::numeric_limits<double>::max(); // time limit
-
- bool m_alive = true; //!< the particle is alive
+ // identity, i.e. things that do not change over the particle lifetime.
+ /// Particle identifier within the event.
+ Barcode m_id;
+ /// Process type specifier.
+ ProcessType m_process = ProcessType::eUndefined;
+ /// PDG particle number.
+ Acts::PdgParticle m_pdg = Acts::PdgParticle::eInvalid;
+ // Particle charge and mass.
+ 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);
+ Vector4 m_position4 = Vector4::Zero();
+ // simulation-specific X0/L0 information and limits
+ // these values are here to simplify the simulation of (nuclear) interactions.
+ // instead of checking at every surface whether an interaction should occur we
+ // can draw an overall limit once. the relevant interaction only needs to
+ // be executed once the limit is reached.
+ // this information is not really particle-specific and should probably be
+ // handled separately. for now, storing it directly here is the simplest
+ // solution.
+ Scalar m_pathX0 = Scalar(0);
+ Scalar m_pathL0 = Scalar(0);
+ Scalar m_limitX0 = std::numeric_limits<Scalar>::max();
+ Scalar m_limitL0 = std::numeric_limits<Scalar>::max();
};
} // namespace ActsFatras
diff --git a/Fatras/include/ActsFatras/Physics/EnergyLoss/BetheBloch.hpp b/Fatras/include/ActsFatras/Physics/EnergyLoss/BetheBloch.hpp
index e047119d09e1711728fc41c76b712c961499c4b3..f588440ea6a80027a618867ea53a1fddffcbeb1f 100644
--- a/Fatras/include/ActsFatras/Physics/EnergyLoss/BetheBloch.hpp
+++ b/Fatras/include/ActsFatras/Physics/EnergyLoss/BetheBloch.hpp
@@ -8,69 +8,68 @@
#pragma once
+#include <vector>
+
#include "Acts/Material/Interactions.hpp"
-#include "ActsFatras/Kernel/LandauDistribution.hpp"
+#include "Acts/Material/MaterialProperties.hpp"
+#include "ActsFatras/EventData/Particle.hpp"
+#include "ActsFatras/Utilities/LandauDistribution.hpp"
namespace ActsFatras {
-/// @brief The struct for the EnergyLoss physics list
-///
-/// This generates the energy loss according to the Bethe-Bloch
-/// description, applying a landau generated enery loss
+/// Simulate energy loss using the Bethe-Bloch/Landau description.
///
-/// It follows the interface of EnergyLoss samplers in Fatras
-/// that could return radiated photons for further processing,
-/// however, for the Bethe-Bloch application the return vector
-/// is always 0.
+/// Energy loss is computed using the most probable value and appropriate
+/// fluctuations from a Landau distribution. No secondaries are generated
+/// for the removed energy.
struct BetheBloch {
/// The flag to include BetheBloch process or not
bool betheBloch = true;
-
/// Scaling for most probable value
double scaleFactorMPV = 1.;
-
/// Scaling for Sigma
double scaleFactorSigma = 1.;
- /// @brief Call operator for the Bethe Bloch energy loss
- ///
- /// @tparam generator_t is a random number generator type
- /// @tparam detector_t is the detector information type
- /// @tparam particle_t is the particle information type
+ /// Simulate energy loss and update the particle parameters.
///
- /// @param[in] generator is the random number generator
- /// @param[in] detector the detector information
- /// @param[in] particle the particle which is being scattered
+ /// @param[in] generator is the random number generator
+ /// @param[in] slab defines the passed material
+ /// @param[in,out] particle is the particle being updated
+ /// @return Empty secondaries containers.
///
- /// @return empty vector for BetheBloch - no secondaries created
- template <typename generator_t, typename detector_t, typename particle_t>
- std::vector<particle_t> operator()(generator_t &generator,
- const detector_t &detector,
- particle_t &particle) const {
+ /// @tparam generator_t is a RandomNumberEngine
+ template <typename generator_t>
+ std::vector<Particle> operator()(generator_t &generator,
+ const Acts::MaterialProperties &slab,
+ Particle &particle) const {
// Do nothing if the flag is set to false
if (not betheBloch) {
return {};
}
- // Create a random landau distribution between in the intervall [0,1]
- LandauDistribution landauDist(0., 1.);
- double landau = landauDist(generator);
- double qop = particle.q() / particle.p();
-
- // @TODO Double investigate if we could do one call
- double energyLoss = Acts::computeEnergyLossLandau(
- detector, particle.pdg(), particle.m(), qop, particle.q());
- double energyLossSigma = Acts::computeEnergyLossLandauSigma(
- detector, particle.pdg(), particle.m(), qop, particle.q());
+ // compute energy loss distribution parameters
+ const auto pdg = particle.pdg();
+ const auto m = particle.mass();
+ const auto qOverP = particle.chargeOverMomentum();
+ const auto q = particle.charge();
+ // most probable value
+ const auto energyLoss =
+ Acts::computeEnergyLossLandau(slab, pdg, m, qOverP, q);
+ // Gaussian-equivalent sigma
+ const auto energyLossSigma =
+ Acts::computeEnergyLossLandauSigma(slab, pdg, m, qOverP, q);
// Simulate the energy loss
- double sampledEnergyLoss = scaleFactorMPV * std::fabs(energyLoss) +
- scaleFactorSigma * energyLossSigma * landau;
+ // TODO landau location and scale parameters are not identical to the most
+ // probable value and the Gaussian-equivalent sigma
+ LandauDistribution lossDistribution(scaleFactorMPV * energyLoss,
+ scaleFactorSigma * energyLossSigma);
+ const auto loss = lossDistribution(generator);
// Apply the energy loss
- particle.energyLoss(sampledEnergyLoss);
+ particle.correctEnergy(-loss);
- // return empty children
+ // Generates no new particles
return {};
}
};
diff --git a/Fatras/include/ActsFatras/Physics/EnergyLoss/BetheHeitler.hpp b/Fatras/include/ActsFatras/Physics/EnergyLoss/BetheHeitler.hpp
index 90e8e21147f770b7488f11cf06c49744575677e8..5890977d338426b6193360a9701f42978c0279bc 100644
--- a/Fatras/include/ActsFatras/Physics/EnergyLoss/BetheHeitler.hpp
+++ b/Fatras/include/ActsFatras/Physics/EnergyLoss/BetheHeitler.hpp
@@ -9,57 +9,54 @@
#pragma once
#include <random>
+#include <vector>
-namespace ActsFatras {
+#include "Acts/Material/MaterialProperties.hpp"
+#include "ActsFatras/EventData/Particle.hpp"
-const double log_2 = std::log(2.);
+namespace ActsFatras {
-/// The struct for the EnergyLoss physics list
+/// Simulate electron energy loss using the Bethe-Heitler description.
///
-/// Bethe-Heitler for electron brem description as described here:
+/// Bethe-Heitler for electron bremsstrahlung description as described here:
/// "A Gaussian-mixture approximation of the Bethe–Heitler model of electron
/// energy loss by bremsstrahlung" R. Frühwirth
-///
struct BetheHeitler {
/// The flag to include BetheHeitler process or not
bool betheHeitler = true;
-
/// A scaling factor to
double scaleFactor = 1.;
- /// @brief Call operator for the Bethe-Heitler energy loss
- ///
- /// @tparam generator_t is a random number generator type
- /// @tparam detector_t is the detector information type
- /// @tparam particle_t is the particle information type
+ /// Simulate energy loss and update the particle parameters.
///
- /// @param[in] generator is the random number generator
- /// @param[in] detector the detector information
- /// @param[in] particle the particle which is being scattered
+ /// @param[in] generator is the random number generator
+ /// @param[in] slab defines the passed material
+ /// @param[in,out] particle is the particle being updated
+ /// @return Empty secondaries containers.
///
- /// @return eventually produced photons
- template <typename generator_t, typename detector_t, typename particle_t>
- std::vector<particle_t> operator()(generator_t &generator,
- const detector_t &detector,
- particle_t &particle) const {
+ /// @tparam generator_t is a RandomNumberEngine
+ template <typename generator_t>
+ std::vector<Particle> operator()(generator_t &generator,
+ const Acts::MaterialProperties &slab,
+ Particle &particle) const {
// Do nothing if the flag is set to false
if (not betheHeitler) {
return {};
}
- double tInX0 = detector.thickness() / detector.material().X0();
-
// Take a random gamma-distributed value - depending on t/X0
- std::gamma_distribution<double> gDist(tInX0 / log_2, 1.);
+ std::gamma_distribution<double> gDist(slab.thicknessInX0() / std::log(2.0),
+ 1.0);
- double u = gDist(generator);
- double z = std::exp(-1. * u);
- double sampledEnergyLoss = std::abs(scaleFactor * particle.E() * (z - 1.));
+ const auto u = gDist(generator);
+ const auto z = std::exp(-u);
+ const auto sampledEnergyLoss =
+ std::abs(scaleFactor * particle.energy() * (z - 1.));
// apply the energy loss
- particle.energyLoss(sampledEnergyLoss);
+ particle.correctEnergy(-sampledEnergyLoss);
- // @TODO return photons, needs particle_creator_t
+ // TODO return the lost energy as a photon
return {};
}
};
diff --git a/Fatras/include/ActsFatras/Physics/Scattering/GaussianMixture.hpp b/Fatras/include/ActsFatras/Physics/Scattering/GaussianMixture.hpp
index adb985087f2735a62fee0a76812098f264fd6641..ceed5c179c5840f57104e2dbc76c8c4a303390da 100644
--- a/Fatras/include/ActsFatras/Physics/Scattering/GaussianMixture.hpp
+++ b/Fatras/include/ActsFatras/Physics/Scattering/GaussianMixture.hpp
@@ -11,15 +11,15 @@
#include <random>
#include "Acts/Material/Interactions.hpp"
+#include "Acts/Material/MaterialProperties.hpp"
+#include "ActsFatras/EventData/Particle.hpp"
namespace ActsFatras {
-/// @brief The struct to be provided to the Scatterer action
-/// This is the gaussian mixture
+/// Generate scattering angles using a Gaussian mixture model.
struct GaussianMixture {
/// Steering parameter
- bool log_include = true;
-
+ bool optGaussianMixtureG4 = false;
double gausMixSigma1_a0 = 8.471e-1;
double gausMixSigma1_a1 = 3.347e-2;
double gausMixSigma1_a2 = -1.843e-3;
@@ -30,43 +30,37 @@ struct GaussianMixture {
double gausMixEpsilon_b1 = 1.106e-1;
double gausMixEpsilon_b2 = -5.729e-3;
- bool optGaussianMixtureG4 = false;
-
- /// @brief Call operator to perform this scattering
+ /// Generate a single 3D scattering angle.
///
- /// @tparam generator_t is a random number generator type
- /// @tparam detector_t is the detector information type
- /// @tparam particle_t is the particle information type
+ /// @param[in] generator is the random number generator
+ /// @param[in] slab defines the passed material
+ /// @param[in,out] particle is the particle being scattered
+ /// @return a 3d scattering angle
///
- /// @param[in] generator is the random number generator
- /// @param[in] detector the detector information
- /// @param[in] particle the particle which is being scattered
- ///
- /// @return a scattering angle in 3D
- template <typename generator_t, typename detector_t, typename particle_t>
- double operator()(generator_t &generator, const detector_t &detector,
- particle_t &particle) const {
+ /// @tparam generator_t is a RandomNumberEngine
+ template <typename generator_t>
+ double operator()(generator_t &generator,
+ const Acts::MaterialProperties &slab,
+ Particle &particle) const {
/// Calculate the highland formula first
- double qop = particle.q() / particle.p();
double sigma = Acts::computeMultipleScatteringTheta0(
- detector, particle.pdg(), particle.m(), qop);
-
+ slab, particle.pdg(), particle.mass(), particle.chargeOverMomentum(),
+ particle.charge());
double sigma2 = sigma * sigma;
// Gauss distribution, will be sampled with generator
std::normal_distribution<double> gaussDist(0., 1.);
-
// Uniform distribution, will be sampled with generator
std::uniform_real_distribution<double> uniformDist(0., 1.);
// Now correct for the tail fraction
// d_0'
double beta2 = particle.beta() * particle.beta();
- double dprime = detector.thickness() / (detector.material().X0() * beta2);
+ double dprime = slab.thickness() / (slab.material().X0() * beta2);
double log_dprime = std::log(dprime);
// d_0''
double log_dprimeprime =
- std::log(std::pow(detector.material().Z(), 2.0 / 3.0) * dprime);
+ std::log(std::pow(slab.material().Z(), 2.0 / 3.0) * dprime);
// get epsilon
double epsilon =
@@ -82,7 +76,7 @@ struct GaussianMixture {
// G4 optimised / native double Gaussian model
if (optGaussianMixtureG4) {
- sigma2 = 225. * dprime / (particle.p() * particle.p());
+ sigma2 = 225. * dprime / (particle.momentum() * particle.momentum());
}
// throw the random number core/tail
if (uniformDist(generator) < epsilon) {
diff --git a/Fatras/include/ActsFatras/Physics/Scattering/GeneralMixture.hpp b/Fatras/include/ActsFatras/Physics/Scattering/GeneralMixture.hpp
index 1b7381c332adc5e9353c02e68ef95d51774bffae..022b897a463554bcd536c1fe8463680a5ed859dd 100644
--- a/Fatras/include/ActsFatras/Physics/Scattering/GeneralMixture.hpp
+++ b/Fatras/include/ActsFatras/Physics/Scattering/GeneralMixture.hpp
@@ -11,47 +11,40 @@
#include <random>
#include "Acts/Material/Interactions.hpp"
+#include "Acts/Material/MaterialProperties.hpp"
+#include "Acts/Utilities/PdgParticle.hpp"
+#include "ActsFatras/EventData/Particle.hpp"
namespace ActsFatras {
-/// This scatter emulated core and tail scattering
+/// Generate scattering angles using a general mixture model.
+///
+/// Emulates core and tail scattering as described in
+///
+/// General mixture model Fruehwirth, M. Liendl.
+/// Comp. Phys. Comm. 141 (2001) 230-246
///
-/// General mixture model Fruehwirth, M. Liendl.
-/// Comp. Phys. Comm. 141 (2001) 230-246
struct GeneralMixture {
/// Steering parameter
bool log_include = true;
-
- //- Scale the mixture level
+ /// Scale the mixture level
double genMixtureScalor = 1.;
- /// @brief Call operator to perform this scattering
+ /// Generate a single 3D scattering angle.
///
- /// @tparam generator_t is a random number generator type
- /// @tparam detector_t is the detector information type
- /// @tparam particle_t is the particle information type
+ /// @param[in] generator is the random number generator
+ /// @param[in] slab defines the passed material
+ /// @param[in,out] particle is the particle being scattered
+ /// @return a 3d scattering angle
///
- /// @param[in] generator is the random number generator
- /// @param[in] detector the detector information
- /// @param[in] particle the particle which is being scattered
- ///
- /// @return a scattering angle in 3D
- template <typename generator_t, typename detector_t, typename particle_t>
- double operator()(generator_t &generator, const detector_t &detector,
- particle_t &particle) const {
- // scale the path length to the radiation length
- // @todo path correction factor
- double tInX0 = detector.thickness() / detector.material().X0();
-
- // material properties
- double Z = detector.material().Z(); // charge layer material
-
- double theta(0.);
-
- if (std::abs(particle.pdg()) != 11) {
- /// Uniform distribution, will be sampled with generator
- std::uniform_real_distribution<double> uniformDist(0., 1.);
+ /// @tparam generator_t is a RandomNumberEngine
+ template <typename generator_t>
+ double operator()(generator_t &generator,
+ const Acts::MaterialProperties &slab,
+ Particle &particle) const {
+ double theta = 0.0;
+ if (std::abs(particle.pdg()) != Acts::PdgParticle::eElectron) {
//----------------------------------------------------------------------------
// see Mixture models of multiple scattering: computation and simulation.
// -
@@ -61,39 +54,37 @@ struct GeneralMixture {
std::array<double, 4> scattering_params;
// Decide which mixture is best
double beta2 = (particle.beta() * particle.beta());
- double tob2 = tInX0 / beta2;
- if (tob2 > 0.6 / std::pow(Z, 0.6)) {
+ double tInX0 = slab.thicknessInX0();
+ double tob2 = slab.thicknessInX0() / beta2;
+ 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.p(), tInX0,
- genMixtureScalor);
+ scattering_params = getGaussian(particle.beta(), particle.momentum(),
+ tInX0, genMixtureScalor);
} else {
scattering_params =
- getGaussmix(particle.beta(), particle.p(), tInX0,
- detector.material().Z(), genMixtureScalor);
+ getGaussmix(particle.beta(), particle.momentum(), tInX0,
+ slab.material().Z(), genMixtureScalor);
}
// Simulate
- theta = gaussmix(uniformDist, generator, scattering_params);
+ theta = gaussmix(generator, scattering_params);
} else {
// Semigaussian mixture - get parameters
auto scattering_params_sg =
- getSemigauss(particle.beta(), particle.p(), tInX0,
- detector.material().Z(), genMixtureScalor);
+ getSemigauss(particle.beta(), particle.momentum(), tInX0,
+ slab.material().Z(), genMixtureScalor);
// Simulate
- theta = semigauss(uniformDist, generator, scattering_params_sg);
+ theta = semigauss(generator, scattering_params_sg);
}
} else {
- /// Gauss distribution, will be sampled with generator
- std::normal_distribution<double> gaussDist(0., 1.);
-
// for electrons we fall back to the Highland (extension)
// return projection factor times sigma times gauss random
- double qop = particle.q() / particle.p();
- double theta0 = Acts::computeMultipleScatteringTheta0(
- detector, particle.pdg(), particle.m(), qop);
- theta = theta0 * gaussDist(generator);
+ const auto theta0 = Acts::computeMultipleScatteringTheta0(
+ slab, particle.pdg(), particle.mass(), particle.chargeOverMomentum(),
+ particle.charge());
+ theta = std::normal_distribution<double>(0.0, theta0)(generator);
}
- // return scaled by sqare root of two
+ // scale from planar to 3d angle
return M_SQRT2 * theta;
}
@@ -163,8 +154,9 @@ struct GeneralMixture {
///
/// @return a double value that represents the gaussian mixture
template <typename generator_t>
- double gaussmix(UniformDist &udist, generator_t &generator,
+ double gaussmix(generator_t &generator,
const std::array<double, 4> &scattering_params) const {
+ std::uniform_real_distribution<double> udist(0.0, 1.0);
double sigma_tot = scattering_params[0];
double var1 = scattering_params[1];
double var2 = scattering_params[2];
@@ -186,8 +178,9 @@ struct GeneralMixture {
///
/// @return a double value that represents the gaussian mixture
template <typename generator_t>
- double semigauss(UniformDist &udist, generator_t &generator,
+ double semigauss(generator_t &generator,
const std::array<double, 6> &scattering_params) const {
+ std::uniform_real_distribution<double> udist(0.0, 1.0);
double a = scattering_params[0];
double b = scattering_params[1];
double var1 = scattering_params[2];
diff --git a/Fatras/include/ActsFatras/Physics/Scattering/Highland.hpp b/Fatras/include/ActsFatras/Physics/Scattering/Highland.hpp
index 256c9c54e94ef55755aec477e73c518c803a51b6..211f2045c8ad3c680fb197e514cada4af3cf4aa2 100644
--- a/Fatras/include/ActsFatras/Physics/Scattering/Highland.hpp
+++ b/Fatras/include/ActsFatras/Physics/Scattering/Highland.hpp
@@ -11,36 +11,34 @@
#include <random>
#include "Acts/Material/Interactions.hpp"
+#include "Acts/Material/MaterialProperties.hpp"
+#include "ActsFatras/EventData/Particle.hpp"
namespace ActsFatras {
-/// @The struct for the Highland-based scattering
+/// Generate scattering angles using the Highland/PDG parametrization.
///
-/// This will scatter particles with a single gaussian distribution
-/// according to the highland formula.
+/// The angles are drawn from a single normal distribution with a width
+/// given by the Highland/PDG formula.
struct Highland {
- /// @brief Call operator to perform this scattering
+ /// Generate a single 3D scattering angle.
///
- /// @tparam generator_t is a random number generator type
- /// @tparam detector_t is the detector information type
- /// @tparam particle_t is the particle information type
+ /// @param[in] generator is the random number generator
+ /// @param[in] slab defines the passed material
+ /// @param[in,out] particle is the particle being scattered
+ /// @return a 3d scattering angle
///
- /// @param[in] generator is the random number generator
- /// @param[in] detector the detector information
- /// @param[in] particle the particle which is being scattered
- ///
- /// @return a scattering angle in 3D
- template <typename generator_t, typename detector_t, typename particle_t>
- double operator()(generator_t &generator, const detector_t &detector,
- particle_t &particle) const {
- // Gauss distribution, will be sampled sampled with generator
- std::normal_distribution<double> gaussDist(0., 1.);
-
- double qop = particle.q() / particle.p();
- double theta0 = Acts::computeMultipleScatteringTheta0(
- detector, particle.pdg(), particle.m(), qop, particle.q());
- // Return projection factor times sigma times grauss random
- return M_SQRT2 * theta0 * gaussDist(generator);
+ /// @tparam generator_t is a RandomNumberEngine
+ template <typename generator_t>
+ double operator()(generator_t &generator,
+ const Acts::MaterialProperties &slab,
+ Particle &particle) const {
+ // compute the planar scattering angle
+ const auto theta0 = Acts::computeMultipleScatteringTheta0(
+ slab, particle.pdg(), particle.mass(), particle.chargeOverMomentum(),
+ particle.charge());
+ // draw from the normal distribution representing the 3d angle distribution
+ return std::normal_distribution<double>(0.0, M_SQRT2 * theta0)(generator);
}
};
diff --git a/Fatras/include/ActsFatras/Physics/Scattering/Scattering.hpp b/Fatras/include/ActsFatras/Physics/Scattering/Scattering.hpp
index 3ede5644e18754efbe4900f1bd84b903418bd80f..96989f212f072006fdbfa2ce66603a55e9d9d827 100644
--- a/Fatras/include/ActsFatras/Physics/Scattering/Scattering.hpp
+++ b/Fatras/include/ActsFatras/Physics/Scattering/Scattering.hpp
@@ -10,53 +10,56 @@
#include <cmath>
#include <random>
+#include <vector>
+#include "Acts/Material/MaterialProperties.hpp"
#include "Acts/Utilities/Definitions.hpp"
#include "Acts/Utilities/Helpers.hpp"
+#include "ActsFatras/EventData/Particle.hpp"
namespace ActsFatras {
-/// @class Scatterering
+/// Simulate (multiple) scattering using a configurable scattering model.
///
-/// This is the (multiple) scattering plugin to the
-/// Physics list. It needs a scattering formula in order
-/// to provide the the scattering angle in 3D space.
-///
-/// There's two options to apply the scattering
-/// - a parametric action that relates phi and theta (default: off)
-/// - an actuall out of direction scattering applying two random numbers
-template <typename formula_t>
+/// @tparam scattering_model_t Model implementation to draw a scattering angle.
+template <typename scattering_model_t>
struct Scattering {
/// The flag to include scattering or not
bool scattering = true;
-
/// Include the log term
bool parametric = false;
- double projectionFactor = 1. / std::sqrt(2.);
-
/// The scattering formula
- formula_t angle;
-
- /// This is the scattering call operator
- template <typename generator_t, typename detector_t, typename particle_t>
- std::vector<particle_t> operator()(generator_t &gen, const detector_t &det,
- particle_t &in) const {
+ scattering_model_t angle;
+
+ /// Simulate scattering and update the particle parameters.
+ ///
+ /// @param[in] generator is the random number generator
+ /// @param[in] slab defines the passed material
+ /// @param[in,out] particle is the particle being updated
+ /// @return Empty secondaries containers.
+ ///
+ /// @tparam generator_t is a RandomNumberEngine
+ template <typename generator_t>
+ std::vector<Particle> operator()(generator_t &generator,
+ const Acts::MaterialProperties &slab,
+ Particle &particle) const {
// Do nothing if the flag is set to false
if (not scattering) {
return {};
}
// 3D scattering angle
- double angle3D = angle(gen, det, in);
+ double angle3D = angle(generator, slab, particle);
// parametric scattering
if (parametric) {
// the initial values
- double theta = Acts::VectorHelpers::theta(in.momentum());
- double phi = Acts::VectorHelpers::phi(in.momentum());
+ double theta = Acts::VectorHelpers::theta(particle.direction());
+ double phi = Acts::VectorHelpers::phi(particle.direction());
double sinTheta = (sin(theta) * sin(theta) > 10e-10) ? sin(theta) : 1.;
// sample them in an independent way
+ const double projectionFactor = 1. / std::sqrt(2.);
double deltaTheta = projectionFactor * angle3D;
double numDetlaPhi = 0.; //?? @THIS IS WRONG HERE !
double deltaPhi = projectionFactor * numDetlaPhi / sinTheta;
@@ -81,16 +84,17 @@ struct Scattering {
double ctheta = std::cos(theta);
// assign the new values
- in.scatter(in.p() * Acts::Vector3D(cphi * stheta, sphi * stheta, ctheta));
+ particle.setDirection(
+ Acts::Vector3D(cphi * stheta, sphi * stheta, ctheta));
} else {
/// uniform distribution
std::uniform_real_distribution<double> uniformDist(0., 1.);
// Create a random uniform distribution between in the intervall [0,1]
- double psi = 2. * M_PI * uniformDist(gen);
+ double psi = 2. * M_PI * uniformDist(generator);
// more complex but "more true"
- Acts::Vector3D pDirection(in.momentum().normalized());
+ Acts::Vector3D pDirection(particle.direction());
double x = -pDirection.y();
double y = pDirection.x();
double z = 0.;
@@ -107,7 +111,7 @@ struct Scattering {
rotation = Acts::AngleAxis3D(psi, pDirection) *
Acts::AngleAxis3D(angle3D, deflector);
// rotate and set a new direction to the cache
- in.scatter(in.p() * rotation * pDirection.normalized());
+ particle.setDirection(rotation * pDirection);
}
// scattering always returns an empty list
// - it is a non-distructive process
diff --git a/Fatras/include/ActsFatras/Selectors/ChargeSelectors.hpp b/Fatras/include/ActsFatras/Selectors/ChargeSelectors.hpp
index 254df7da56ae725ece0f334137e5181f0edf4a54..66cc9cd6f20d05bff7f68f2817f5294d4ce903ca 100644
--- a/Fatras/include/ActsFatras/Selectors/ChargeSelectors.hpp
+++ b/Fatras/include/ActsFatras/Selectors/ChargeSelectors.hpp
@@ -8,37 +8,39 @@
#pragma once
+#include "ActsFatras/EventData/Particle.hpp"
+
namespace ActsFatras {
-struct ChargedSelector {
- /// Return true for all particles with charge != 0.
- template <typename detector_t, typename particle_t>
- bool operator()(const detector_t &, const particle_t &particle) const {
- return (particle.q() * particle.q() > 0.);
+/// Select neutral particles.
+struct NeutralSelector {
+ template <typename detector_t>
+ bool operator()(const detector_t &, const Particle &particle) const {
+ return (particle.charge() == Particle::Scalar(0));
}
};
-struct NeutralSelector {
- /// Return true for all particles with charge == 0.
- template <typename detector_t, typename particle_t>
- bool operator()(const detector_t &, const particle_t &particle) const {
- return (particle.q() == 0.);
+/// Select all charged particles.
+struct ChargedSelector {
+ template <typename detector_t>
+ bool operator()(const detector_t &, const Particle &particle) const {
+ return (particle.charge() != Particle::Scalar(0));
}
};
+/// Select positively charged particles.
struct PositiveSelector {
- /// Return true for all particles with charge > 0.
- template <typename detector_t, typename particle_t>
- bool operator()(const detector_t &, const particle_t &particle) const {
- return (particle.q() > 0.);
+ template <typename detector_t>
+ bool operator()(const detector_t &, const Particle &particle) const {
+ return (Particle::Scalar(0) < particle.charge());
}
};
+/// Select negatively charged particles.
struct NegativeSelector {
- /// Return true for all particles with charge<> 0.
- template <typename detector_t, typename particle_t>
- bool operator()(const detector_t &, const particle_t &particle) const {
- return (particle.q() * particle.q() > 0.);
+ template <typename detector_t>
+ bool operator()(const detector_t &, const Particle &particle) const {
+ return (particle.charge() < Particle::Scalar(0));
}
};
diff --git a/Fatras/include/ActsFatras/Selectors/KinematicCasts.hpp b/Fatras/include/ActsFatras/Selectors/KinematicCasts.hpp
index 357a2d70ffb6f8e10d57528e2d62918dedcd3394..1c25b17ac36230fe6dc37aeab90cca1354556f9e 100644
--- a/Fatras/include/ActsFatras/Selectors/KinematicCasts.hpp
+++ b/Fatras/include/ActsFatras/Selectors/KinematicCasts.hpp
@@ -11,71 +11,65 @@
#include <cmath>
#include "Acts/Utilities/Helpers.hpp"
+#include "ActsFatras/EventData/Particle.hpp"
namespace ActsFatras {
namespace Casts {
-/// The Eta cast operator
-struct eta {
- template <typename particle_t>
- double operator()(const particle_t &particle) const {
- return Acts::VectorHelpers::eta(particle.momentum());
+/// Retrieve the transverse absolute distance of the vertex to the origin.
+struct Vrho {
+ double operator()(const Particle& particle) const {
+ return Acts::VectorHelpers::perp(particle.position());
}
};
-/// The Eta cast operator
-struct absEta {
- template <typename particle_t>
- double operator()(const particle_t &particle) const {
- return std::abs(Acts::VectorHelpers::eta(particle.momentum()));
+/// Retrieve the longitudinal distance of the vertex to the origin.
+struct Vz {
+ double operator()(const Particle& particle) const {
+ return particle.position().z();
}
};
-/// The Pt cast operator
-struct pT {
- template <typename particle_t>
- double operator()(const particle_t &particle) const {
- return Acts::VectorHelpers::perp(particle.momentum());
+/// Retrieve the longitudinal absolute distance of the vertex to the origin.
+struct AbsVz {
+ double operator()(const Particle& particle) const {
+ return std::abs(particle.position().z());
}
};
-/// The P cast operator
-struct p {
- template <typename particle_t>
- double operator()(const particle_t &particle) const {
- return particle.momentum().norm();
+/// Retrieve the direction pseudo-rapidity.
+struct Eta {
+ double operator()(const Particle& particle) const {
+ return Acts::VectorHelpers::eta(particle.direction());
}
};
-/// The E cast operator
-struct E {
- template <typename particle_t>
- double operator()(const particle_t &particle) const {
- return particle.E();
+/// Retrieve the direction absolute pseudo-rapidity.
+struct AbsEta {
+ double operator()(const Particle& particle) const {
+ return std::abs(Acts::VectorHelpers::eta(particle.direction()));
}
};
-/// The E cast operator
-struct vR {
- template <typename particle_t>
- double operator()(const particle_t &particle) const {
- return Acts::VectorHelpers::perp(particle.position());
+/// Retrieve the transverse momentum.
+struct Pt {
+ double operator()(const Particle& particle) const {
+ return particle.momentum() *
+ Acts::VectorHelpers::perp(particle.direction());
}
};
-/// The E cast operator
-struct vZ {
- template <typename particle_t>
- double operator()(const particle_t &particle) const {
- return particle.position().z();
+/// Retrieve the absolute momentum.
+struct P {
+ double operator()(const Particle& particle) const {
+ return particle.momentum();
}
};
-/// The E cast operator
-struct AbsVz {
- template <typename particle_t>
- double operator()(const particle_t &particle) const {
- return std::abs(particle.position().z());
+/// Retrieve the total energy.
+struct E {
+ double operator()(const Particle& particle) const {
+ return particle.energy();
}
};
diff --git a/Fatras/include/ActsFatras/Selectors/LimitSelectors.hpp b/Fatras/include/ActsFatras/Selectors/LimitSelectors.hpp
index 36d4eb23f1d9fc51a3f7e46b65994a799e3563cf..8ebad10eae8982db206e52a39690df73cb23dc38 100644
--- a/Fatras/include/ActsFatras/Selectors/LimitSelectors.hpp
+++ b/Fatras/include/ActsFatras/Selectors/LimitSelectors.hpp
@@ -8,31 +8,26 @@
#pragma once
+#include "Acts/Material/MaterialProperties.hpp"
+#include "ActsFatras/EventData/Particle.hpp"
+
namespace ActsFatras {
-struct X0Limit {
- /// Return true if the limit in X0 is reached
- /// @todo modify, return what's left from the detector: needs non-const
- /// detector
- template <typename detector_t, typename particle_t>
- bool operator()(const detector_t &detector,
- const particle_t &particle) const {
- return particle.pathInX0() +
- detector.thickness() / detector.material().X0() >=
- particle.limitInX0();
+/// Select particles whose X0 limit would be reached after material passage.
+struct PathLimitX0 {
+ bool operator()(const Acts::MaterialProperties &slab,
+ const Particle &particle) const {
+ return particle.pathLimitX0() <
+ (particle.pathInX0() + slab.thicknessInX0());
}
};
-struct L0Limit {
- /// Return true if the limit in X0 is reached
- /// @todo modify, return what's left from the detector: needs non-const
- /// detector
- template <typename detector_t, typename particle_t>
- bool operator()(const detector_t &detector,
- const particle_t &particle) const {
- return particle.pathInL0() +
- detector.thickness() / detector.material().L0() >=
- particle.limitInL0();
+/// Select particles whose L0 limit would be reached after material passage.
+struct PathLimitL0 {
+ bool operator()(const Acts::MaterialProperties &slab,
+ const Particle &particle) const {
+ return particle.pathLimitL0() <
+ (particle.pathInL0() + slab.thicknessInL0());
}
};
diff --git a/Fatras/include/ActsFatras/Selectors/PdgSelectors.hpp b/Fatras/include/ActsFatras/Selectors/PdgSelectors.hpp
index 99fcdfc8e8a075312cc8471c0d5f9f82a96286fb..9017340800601e46e9555c625d5ce6157406c0b4 100644
--- a/Fatras/include/ActsFatras/Selectors/PdgSelectors.hpp
+++ b/Fatras/include/ActsFatras/Selectors/PdgSelectors.hpp
@@ -8,57 +8,49 @@
#pragma once
-namespace ActsFatras {
+#include <cstdlib>
-template <int pdg_t>
-struct AbsPdgSelector {
- // absolute Pdg selection
- const int saPDG = pdg_t;
+#include "ActsFatras/EventData/Particle.hpp"
- /// Return true for all particles with | pdg | matching
- /// the selection criteria
- template <typename detector_t, typename particle_t>
- bool operator()(const detector_t &, const particle_t &particle) const {
- return (particle.pdg() * particle.pdg() == saPDG * saPDG);
- }
-};
+namespace ActsFatras {
-template <int pdg_t>
+/// Select particles of one specific type.
+///
+/// Particle and Antiparticle are treated as two separate types.
+template <int Pdg>
struct PdgSelector {
- // Pdg selection
- const int saPDG = pdg_t;
-
- /// Return true for all particles with pdg matching
- /// the selection criteria
- template <typename detector_t, typename particle_t>
- bool operator()(const detector_t &, const particle_t &particle) const {
- return (particle.pdg() == saPDG);
+ template <typename detector_t>
+ bool operator()(const detector_t &, const Particle &particle) const {
+ return (static_cast<int>(particle.pdg()) == Pdg);
}
};
-template <int pdg_t>
-struct AbsPdgExcluder {
- // absolute Pdg selection
- const int saPDG = pdg_t;
-
- /// Return true for all particles with | pdg | matching
- /// the selection criteria
- template <typename detector_t, typename particle_t>
- bool operator()(const detector_t &, const particle_t &particle) const {
- return !(particle.pdg() * particle.pdg() == saPDG * saPDG);
+/// Select particles and antiparticles of one specific type.
+template <int Pdg>
+struct AbsPdgSelector {
+ template <typename detector_t>
+ bool operator()(const detector_t &, const Particle &particle) const {
+ return (std::abs(static_cast<int>(particle.pdg())) == std::abs(Pdg));
}
};
-template <int pdg_t>
+/// Select all particles except one specific type.
+///
+/// Particle and Antiparticle are treated as two separate types.
+template <int Pdg>
struct PdgExcluder {
- // Pdg selection
- const int saPDG = pdg_t;
+ template <typename detector_t>
+ bool operator()(const detector_t &, const Particle &particle) const {
+ return (static_cast<int>(particle.pdg()) != Pdg);
+ }
+};
- /// Return true for all particles with pdg matching
- /// the selection criteria
- template <typename detector_t, typename particle_t>
- bool operator()(const detector_t &, const particle_t &particle) const {
- return !(particle.pdg() == saPDG);
+/// Select all particles except for (anti-)particles of one specific type.
+template <int Pdg>
+struct AbsPdgExcluder {
+ template <typename detector_t>
+ bool operator()(const detector_t &, const Particle &particle) const {
+ return (std::abs(static_cast<int>(particle.pdg())) != std::abs(Pdg));
}
};
diff --git a/Fatras/include/ActsFatras/Selectors/SelectorHelpers.hpp b/Fatras/include/ActsFatras/Selectors/SelectorHelpers.hpp
index 05f16fdb91db364bca04465982f9cbf2b0072626..c5d47b54719f4e4ad7805420531a7966d95cc846 100644
--- a/Fatras/include/ActsFatras/Selectors/SelectorHelpers.hpp
+++ b/Fatras/include/ActsFatras/Selectors/SelectorHelpers.hpp
@@ -10,49 +10,44 @@
#include <climits>
+#include "ActsFatras/EventData/Particle.hpp"
+
namespace ActsFatras {
-// static selectors
+/// Select all objects with an extracted value equal or larger than the cut.
template <typename cast_t>
struct Min {
- cast_t cast;
double valMin = 0.;
- /// Return true for all particles with transverse momentum
- /// bigger than the specified minimum value
- template <typename detector_t, typename particle_t>
- bool operator()(const detector_t &, const particle_t &particle) const {
- double val = cast(particle);
- return (val >= valMin);
+ template <typename detector_t, typename T>
+ bool operator()(const detector_t &, const T &thing) const {
+ return (valMin <= cast_t()(thing));
}
};
+/// Select all objects with an extracted value below the cut.
template <typename cast_t>
struct Max {
- cast_t cast;
double valMax = std::numeric_limits<double>::max();
- /// Return true for all particles with transverse momentum
- /// bigger than the specified minimum value
- template <typename detector_t, typename particle_t>
- bool operator()(const detector_t &, const particle_t &particle) const {
- double val = cast(particle);
- return (val <= valMax);
+ template <typename detector_t, typename T>
+ bool operator()(const detector_t &, const T &thing) const {
+ return (cast_t()(thing) < valMax);
}
};
+/// Select all objects with an extracted value within the range.
+///
+/// The range is defined as the left, half-open interval within the cuts.
template <typename cast_t>
struct Range {
- cast_t cast;
- double valMin = 0.;
+ double valMin = std::numeric_limits<double>::lowest();
double valMax = std::numeric_limits<double>::max();
- /// Return true for all particles with transverse momentum
- /// within the specified range
- template <typename detector_t, typename particle_t>
- bool operator()(const detector_t &, const particle_t &particle) const {
- double val = cast(particle);
- return (val >= valMin && val <= valMax);
+ template <typename detector_t, typename T>
+ bool operator()(const detector_t &, const T &thing) const {
+ const auto val = cast_t()(thing);
+ return ((valMin <= val) and (val < valMax));
}
};
diff --git a/Fatras/include/ActsFatras/Kernel/LandauDistribution.hpp b/Fatras/include/ActsFatras/Utilities/LandauDistribution.hpp
similarity index 80%
rename from Fatras/include/ActsFatras/Kernel/LandauDistribution.hpp
rename to Fatras/include/ActsFatras/Utilities/LandauDistribution.hpp
index 7f6becece237ded613145980da01870fde6681a4..54f1a0ef4b3be1102ce45f9a91717a7a44b3d9ab 100644
--- a/Fatras/include/ActsFatras/Kernel/LandauDistribution.hpp
+++ b/Fatras/include/ActsFatras/Utilities/LandauDistribution.hpp
@@ -22,13 +22,16 @@ class LandauDistribution {
/// Parameters must link back to the host distribution.
using distribution_type = LandauDistribution;
- /// Mean of the Landau distribution
- double mean = 0.;
- /// Scale of the Landau distribution
- double scale = 1.;
+ /// Location parameter.
+ ///
+ /// @warning This is neither the mean nor the most probable value.
+ double location = 0.0;
+ /// Scale parameter.
+ double scale = 1.0;
/// Construct from parameters.
- param_type(double mean_, double scale_) : mean(mean_), scale(scale_) {}
+ param_type(double location_, double scale_)
+ : location(location_), scale(scale_) {}
// Explicitlely defaulted construction and assignment
param_type() = default;
param_type(const param_type &) = default;
@@ -38,7 +41,7 @@ class LandauDistribution {
/// Parameters should be EqualityComparable
friend bool operator==(const param_type &lhs, const param_type &rhs) {
- return (lhs.mean == rhs.mean) and (lhs.scale == rhs.scale);
+ return (lhs.location == rhs.location) and (lhs.scale == rhs.scale);
}
friend bool operator!=(const param_type &lhs, const param_type &rhs) {
return not(lhs == rhs);
@@ -48,7 +51,7 @@ class LandauDistribution {
using result_type = double;
/// Construct directly from the distribution parameters.
- LandauDistribution(double mean, double scale) : m_cfg(mean, scale) {}
+ LandauDistribution(double location, double scale) : m_cfg(location, scale) {}
/// Construct from a parameter object.
LandauDistribution(const param_type &cfg) : m_cfg(cfg) {}
// Explicitlely defaulted construction and assignment
@@ -72,14 +75,14 @@ class LandauDistribution {
/// Generate a random number from the configured Landau distribution.
template <typename Generator>
- result_type operator()(Generator &engine) {
- return (*this)(engine, m_cfg);
+ result_type operator()(Generator &generator) {
+ return (*this)(generator, m_cfg);
}
/// Generate a random number from the given Landau distribution.
template <typename Generator>
- result_type operator()(Generator &engine, const param_type ¶ms) {
- double x = std::generate_canonical<double, 10>(engine);
- return params.mean + quantile(x, params.scale);
+ result_type operator()(Generator &generator, const param_type ¶ms) {
+ const auto z = std::uniform_real_distribution<double>()(generator);
+ return params.location + params.scale * quantile(z);
}
/// Provide standard comparison operators
@@ -95,7 +98,7 @@ class LandauDistribution {
private:
param_type m_cfg;
- static double quantile(double z, double xi);
+ static double quantile(double z);
};
} // namespace ActsFatras
diff --git a/Fatras/src/LandauDistribution.cpp b/Fatras/src/LandauDistribution.cpp
index 1c605f35d009e6e7b150b012e19538a17df0ab5f..c9166bc6fc3beb6f1c9605c7c034ed995455b10d 100644
--- a/Fatras/src/LandauDistribution.cpp
+++ b/Fatras/src/LandauDistribution.cpp
@@ -6,11 +6,10 @@
// 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/Kernel/LandauDistribution.hpp"
+#include "ActsFatras/Utilities/LandauDistribution.hpp"
-double ActsFatras::LandauDistribution::quantile(double z, double xi) {
+double ActsFatras::LandauDistribution::quantile(double z) {
// LANDAU quantile : algorithm from CERNLIB G110 ranlan
- // with scale parameter xi
// Converted by Rene Brun from CERNLIB routine ranlan(G110),
// Moved and adapted to QuantFuncMathCore by B. List 29.4.2010
@@ -180,8 +179,6 @@ double ActsFatras::LandauDistribution::quantile(double z, double xi) {
40.157721, 41.622399, 43.202525, 44.912465, 46.769077, 48.792279,
51.005773, 53.437996, 56.123356, 59.103894};
- if (xi <= 0)
- return 0;
if (z <= 0)
return -std::numeric_limits<double>::infinity();
if (z >= 1)
@@ -214,5 +211,5 @@ double ActsFatras::LandauDistribution::quantile(double z, double xi) {
((1 + 6.06511919E3 * u + 6.94021044E5 * v) * u);
}
}
- return xi * ranlan;
+ return ranlan;
}
diff --git a/Tests/UnitTests/Fatras/CMakeLists.txt b/Tests/UnitTests/Fatras/CMakeLists.txt
index f78d03b06ecb9402d278abd34673168b2e61e5be..d9bcfa9550ea8558f74ede0f1d1825e0c8e4ff85 100644
--- a/Tests/UnitTests/Fatras/CMakeLists.txt
+++ b/Tests/UnitTests/Fatras/CMakeLists.txt
@@ -1,3 +1,4 @@
+add_subdirectory(EventData)
add_subdirectory(Kernel)
add_subdirectory(Physics)
add_subdirectory(Selectors)
diff --git a/Tests/UnitTests/Fatras/EventData/BarcodeTests.cpp b/Tests/UnitTests/Fatras/EventData/BarcodeTests.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..254f5ace6fe134d64561c6381ed26e7210797b0c
--- /dev/null
+++ b/Tests/UnitTests/Fatras/EventData/BarcodeTests.cpp
@@ -0,0 +1,57 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2018-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 <random>
+
+#include "ActsFatras/EventData/Barcode.hpp"
+
+using ActsFatras::Barcode;
+
+BOOST_AUTO_TEST_SUITE(FatrasBarcode)
+
+BOOST_AUTO_TEST_CASE(MakeDescendant) {
+ // initial barcode with primary particle
+ auto p3 = Barcode().setVertexPrimary(1).setVertexSecondary(2).setParticle(3);
+ BOOST_TEST(p3.vertexPrimary() == 1u);
+ BOOST_TEST(p3.vertexSecondary() == 2u);
+ BOOST_TEST(p3.particle() == 3u);
+ BOOST_TEST(p3.generation() == 0u);
+ BOOST_TEST(p3.subParticle() == 0u);
+
+ // generate two first-generation descendants
+ auto p30 = p3.makeDescendant(0);
+ auto p31 = p3.makeDescendant(1);
+ BOOST_TEST(p30.vertexPrimary() == 1u);
+ BOOST_TEST(p30.vertexSecondary() == 2u);
+ BOOST_TEST(p30.particle() == 3u);
+ BOOST_TEST(p30.generation() == 1u);
+ BOOST_TEST(p30.subParticle() == 0u);
+ BOOST_TEST(p31.vertexPrimary() == 1u);
+ BOOST_TEST(p31.vertexSecondary() == 2u);
+ BOOST_TEST(p31.particle() == 3u);
+ BOOST_TEST(p31.generation() == 1u);
+ BOOST_TEST(p31.subParticle() == 1u);
+
+ // generate two (overlapping) second-generation descendants.
+ auto p300 = p30.makeDescendant(0);
+ auto p310 = p31.makeDescendant(0);
+ BOOST_TEST(p300.vertexPrimary() == 1u);
+ BOOST_TEST(p300.vertexSecondary() == 2u);
+ BOOST_TEST(p300.particle() == 3u);
+ BOOST_TEST(p300.generation() == 2u);
+ BOOST_TEST(p300.subParticle() == 0u);
+ BOOST_TEST(p310.vertexPrimary() == 1u);
+ BOOST_TEST(p310.vertexSecondary() == 2u);
+ BOOST_TEST(p310.particle() == 3u);
+ BOOST_TEST(p310.generation() == 2u);
+ BOOST_TEST(p310.subParticle() == 0u);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/Tests/UnitTests/Fatras/EventData/CMakeLists.txt b/Tests/UnitTests/Fatras/EventData/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f0f5147602f28c3dbef7b5d675e712e68ab00151
--- /dev/null
+++ b/Tests/UnitTests/Fatras/EventData/CMakeLists.txt
@@ -0,0 +1,4 @@
+set(unittest_extra_libraries ActsFatras)
+
+add_unittest(FatrasBarcode BarcodeTests.cpp)
+add_unittest(FatrasParticle ParticleTests.cpp)
diff --git a/Tests/UnitTests/Fatras/EventData/ParticleTests.cpp b/Tests/UnitTests/Fatras/EventData/ParticleTests.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c2be28f87499261488ce53562f0b9ef8561934cc
--- /dev/null
+++ b/Tests/UnitTests/Fatras/EventData/ParticleTests.cpp
@@ -0,0 +1,76 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2018-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 <random>
+
+#include "Acts/Utilities/Units.hpp"
+#include "ActsFatras/EventData/Particle.hpp"
+
+using Acts::PdgParticle;
+using ActsFatras::Barcode;
+using ActsFatras::Particle;
+using namespace Acts::UnitLiterals;
+
+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);
+
+ BOOST_TEST(particle.id() == id);
+ BOOST_TEST(particle.pdg() == PdgParticle::eProton);
+ // particle is at rest at the origin
+ BOOST_TEST(particle.position4() == Particle::Vector4::Zero());
+ BOOST_TEST(particle.position() == Particle::Vector3::Zero());
+ BOOST_TEST(particle.time() == Particle::Scalar(0));
+ BOOST_TEST(particle.position4().x() == particle.position().x());
+ 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 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)
+ .setDirection(Particle::Vector3::UnitX())
+ .setMomentum(2_GeV);
+
+ BOOST_TEST(particle.mass() == 1_GeV);
+ // check that the particle has some input energy
+ BOOST_TEST(particle.momentum4().x() == 2_GeV);
+ 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.energy() == std::hypot(1_GeV, 2_GeV));
+ // loose some energy
+ particle.correctEnergy(-100_MeV);
+ BOOST_TEST(particle.momentum() < 2_GeV);
+ BOOST_TEST(particle.energy() ==
+ Particle::Scalar(std::hypot(1_GeV, 2_GeV) - 100_MeV));
+ // 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.energy() == particle.mass());
+ // lossing even more energy does nothing
+ particle.correctEnergy(-10_GeV);
+ BOOST_TEST(particle.momentum() == Particle::Scalar(0));
+ BOOST_TEST(particle.energy() == particle.mass());
+ // particle is not alive anymore
+ BOOST_TEST(not particle);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/Tests/UnitTests/Fatras/Physics/Dataset.hpp b/Tests/UnitTests/Fatras/Physics/Dataset.hpp
index 9406efd00d7465cca9019dc8349cec36ea8c1686..7591d824c9b213bcc5fe0fd3efb8e2236975acea 100644
--- a/Tests/UnitTests/Fatras/Physics/Dataset.hpp
+++ b/Tests/UnitTests/Fatras/Physics/Dataset.hpp
@@ -57,12 +57,14 @@ const auto particleParameters = momentumPhi * momentumLambda * momentumAbs *
// utility function to build a particle from the dataset parameters
inline ActsFatras::Particle makeParticle(double phi, double lambda, double p,
- int pdg, double m, double q) {
- Acts::Vector3D pos(0, 0, 0);
- Acts::Vector3D mom(p * std::cos(lambda) * std::cos(phi),
- p * std::cos(lambda) * std::sin(phi),
- p * std::sin(lambda));
- return {pos, mom, m, q, pdg};
+ Acts::PdgParticle pdg, double m,
+ double q) {
+ const auto id = ActsFatras::Barcode().setVertexPrimary(1).setParticle(1);
+ return ActsFatras::Particle(id, pdg, m, q)
+ .setPosition4(0, 0, 0, 0)
+ .setDirection(std::cos(lambda) * std::cos(phi),
+ std::cos(lambda) * std::sin(phi), std::sin(lambda))
+ .setMomentum(p);
}
} // namespace
diff --git a/Tests/UnitTests/Fatras/Physics/EnergyLossTests.cpp b/Tests/UnitTests/Fatras/Physics/EnergyLossTests.cpp
index 18dfb89375d13d435d50857593c1c378b98cf349..06561bec6bd6a86bf3cacbbf6f94e4a4b011e0ef 100644
--- a/Tests/UnitTests/Fatras/Physics/EnergyLossTests.cpp
+++ b/Tests/UnitTests/Fatras/Physics/EnergyLossTests.cpp
@@ -30,10 +30,8 @@ BOOST_DATA_TEST_CASE(BetheBloch, Dataset::particleParameters, phi, lambda, p,
ActsFatras::BetheBloch process;
const auto outgoing = process(gen, Dataset::thickSlab, after);
// energy loss changes momentum and energy
- // TODO fix process computation
- // BOOST_TEST(after.pT() < before.pT());
- // BOOST_TEST(after.p() < before.p());
- // BOOST_TEST(after.E() < before.E());
+ BOOST_TEST(after.momentum() < before.momentum());
+ BOOST_TEST(after.energy() < before.energy());
// energy loss creates no new particles
BOOST_TEST(outgoing.empty());
}
@@ -48,10 +46,8 @@ BOOST_DATA_TEST_CASE(BetheHeitler, Dataset::particleParameters, phi, lambda, p,
ActsFatras::BetheHeitler process;
const auto outgoing = process(gen, Dataset::thickSlab, after);
// energy loss changes momentum and energy
- // TODO fix process computation
- // BOOST_TEST(after.pT() < before.pT());
- // BOOST_TEST(after.p() < before.p());
- // BOOST_TEST(after.E() < before.E());
+ BOOST_TEST(after.momentum() < before.momentum());
+ 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 8ff257b8fe9c246eba18fa2908f3881931891bcf..d8c8186dc63b0db6bbf35872f78825ccf7c45199 100644
--- a/Tests/UnitTests/Fatras/Physics/ScatteringTests.cpp
+++ b/Tests/UnitTests/Fatras/Physics/ScatteringTests.cpp
@@ -15,33 +15,55 @@
#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
#include "Acts/Tests/CommonHelpers/PredefinedMaterials.hpp"
#include "ActsFatras/EventData/Particle.hpp"
+#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"
-// TODO decrease this after removing computation round trip from particle
-static constexpr double eps = 1e-4;
+namespace {
-BOOST_AUTO_TEST_SUITE(FatrasScattering)
+constexpr double eps = 1e-10;
-/// Test the scattering implementation
-BOOST_DATA_TEST_CASE(HighlandScattering, Dataset::particleParameters, phi,
- lambda, p, pdg, m, q) {
- using HighlandScattering = ActsFatras::Scattering<ActsFatras::Highland>;
+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 run(const Scattering& scattering, const ActsFatras::Particle& before) {
std::default_random_engine gen;
- ActsFatras::Particle before =
- Dataset::makeParticle(phi, lambda, p, pdg, m, q);
ActsFatras::Particle after = before;
- HighlandScattering scattering;
const auto outgoing = scattering(gen, Dataset::thinSlab, after);
// scattering leaves absolute energy/momentum unchanged
- CHECK_CLOSE_REL(after.pT(), before.pT(), eps);
- CHECK_CLOSE_REL(after.p(), before.p(), eps);
- CHECK_CLOSE_REL(after.E(), before.E(), eps);
+ CHECK_CLOSE_REL(after.momentum(), before.momentum(), eps);
+ CHECK_CLOSE_REL(after.energy(), before.energy(), eps);
// scattering creates no new particles
BOOST_TEST(outgoing.empty());
}
+} // namespace
+
+BOOST_AUTO_TEST_SUITE(FatrasScattering)
+
+BOOST_DATA_TEST_CASE(GeneralMixture, Dataset::particleParameters, phi, lambda,
+ p, pdg, m, q) {
+ run(GeneralMixtureScattering(),
+ Dataset::makeParticle(phi, lambda, p, pdg, m, q));
+}
+
+BOOST_DATA_TEST_CASE(GaussianMixture, Dataset::particleParameters, phi, lambda,
+ p, pdg, m, q) {
+ run(GaussianMixtureScattering(),
+ Dataset::makeParticle(phi, lambda, p, pdg, m, q));
+}
+
+BOOST_DATA_TEST_CASE(Highland, Dataset::particleParameters, phi, lambda, p, pdg,
+ m, q) {
+ run(HighlandScattering(), Dataset::makeParticle(phi, lambda, p, pdg, m, q));
+}
+
BOOST_AUTO_TEST_SUITE_END()
diff --git a/Tests/UnitTests/Fatras/Selectors/CMakeLists.txt b/Tests/UnitTests/Fatras/Selectors/CMakeLists.txt
index 86a4a0f82cdcc77177c4f4d16770a59b49ab4100..50112c684dce66d6a490b3a9b0cf6f84245bbce3 100644
--- a/Tests/UnitTests/Fatras/Selectors/CMakeLists.txt
+++ b/Tests/UnitTests/Fatras/Selectors/CMakeLists.txt
@@ -1,5 +1,6 @@
set(unittest_extra_libraries ActsFatras)
+add_unittest(FatrasChargeSelectors ChargeSelectorsTests.cpp)
add_unittest(FatrasKinematicCasts KinematicCastsTests.cpp)
add_unittest(FatrasLimitSelectors LimitSelectorsTests.cpp)
add_unittest(FatrasPdgSelectors PdgSelectorsTests.cpp)
diff --git a/Tests/UnitTests/Fatras/Selectors/ChargeSelectorsTests.cpp b/Tests/UnitTests/Fatras/Selectors/ChargeSelectorsTests.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5695ebe886f5ebb9751bf1973b3497ab061e6877
--- /dev/null
+++ b/Tests/UnitTests/Fatras/Selectors/ChargeSelectorsTests.cpp
@@ -0,0 +1,47 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2018-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 "ActsFatras/Selectors/ChargeSelectors.hpp"
+#include "Dataset.hpp"
+
+using namespace ActsFatras;
+
+static const auto& detector = Dataset::thinSlab;
+
+BOOST_AUTO_TEST_SUITE(FatrasChargeSelectors)
+
+BOOST_AUTO_TEST_CASE(NegativeParticle) {
+ const auto& particle = Dataset::centralElectron;
+
+ BOOST_TEST(not NeutralSelector()(detector, particle));
+ BOOST_TEST(ChargedSelector()(detector, particle));
+ BOOST_TEST(not PositiveSelector()(detector, particle));
+ BOOST_TEST(NegativeSelector()(detector, particle));
+}
+
+BOOST_AUTO_TEST_CASE(NeutralParticle) {
+ const auto& particle = Dataset::centralNeutron;
+
+ BOOST_TEST(NeutralSelector()(detector, particle));
+ BOOST_TEST(not ChargedSelector()(detector, particle));
+ BOOST_TEST(not PositiveSelector()(detector, particle));
+ BOOST_TEST(not NegativeSelector()(detector, particle));
+}
+
+BOOST_AUTO_TEST_CASE(PositiveParticle) {
+ const auto& particle = Dataset::centralPositron;
+
+ BOOST_TEST(not NeutralSelector()(detector, particle));
+ BOOST_TEST(ChargedSelector()(detector, particle));
+ BOOST_TEST(PositiveSelector()(detector, particle));
+ BOOST_TEST(not NegativeSelector()(detector, particle));
+}
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/Tests/UnitTests/Fatras/Selectors/Dataset.hpp b/Tests/UnitTests/Fatras/Selectors/Dataset.hpp
index af9b15cce4ce8f260b6c95d9161f5e3ccc04a26c..1287661c020c2d19068c092fea97cb6eb33e87cb 100644
--- a/Tests/UnitTests/Fatras/Selectors/Dataset.hpp
+++ b/Tests/UnitTests/Fatras/Selectors/Dataset.hpp
@@ -27,25 +27,31 @@ constexpr auto massElectron = 0.51099891_MeV;
constexpr auto massMuon = 105.658367_MeV;
constexpr auto massPion = 134.9766_MeV;
-ActsFatras::Particle centralElectron({0_mm, 0_mm, 0_mm},
- {0_GeV, 1.5_GeV, 0_GeV}, massElectron,
- -1_e, Acts::PdgParticle::eElectron);
-ActsFatras::Particle centralPositron({0_mm, 0_mm, 0_mm},
- {0_GeV, 1.5_GeV, 0_GeV}, massElectron, 1_e,
- Acts::PdgParticle::ePositron);
-ActsFatras::Particle centralMuon({0_mm, 0_mm, 0_mm}, {0_GeV, 1.5_GeV, 0_GeV},
- massMuon, -1_e, Acts::PdgParticle::eMuon);
-ActsFatras::Particle centralAntiMuon({0_mm, 0_mm, 0_mm},
- {0_GeV, 1.5_GeV, 0_GeV}, massMuon, 1_e,
- Acts::PdgParticle::eAntiMuon);
-ActsFatras::Particle backwardPion({0_mm, 0_mm, -100_mm},
- {10_MeV, 10_MeV, -1.5_GeV}, massPion, -1_e,
- Acts::PdgParticle::ePionMinus);
-ActsFatras::Particle centralPion({0_mm, 0_mm, 0_mm}, {0_GeV, 1.5_GeV, 0_GeV},
- massPion, -1_e, Acts::PdgParticle::ePionMinus);
-ActsFatras::Particle forwardPion({0_mm, 0_mm, 100_mm},
- {10_MeV, 10_MeV, 1.5_GeV}, massPion, -1_e,
- Acts::PdgParticle::ePionMinus);
+ActsFatras::Particle makeParticle(double z, double eta, Acts::PdgParticle pdg,
+ double m, double q) {
+ const auto id = ActsFatras::Barcode().setVertexPrimary(1).setParticle(1);
+ return ActsFatras::Particle(id, pdg, m, q)
+ .setPosition4(0.0, 0.0, z, 0.0)
+ .setDirection(1.0 / std::cosh(eta), 0.0, std::tanh(eta))
+ .setMomentum(1.5_GeV);
+}
+
+const auto centralElectron =
+ makeParticle(0_mm, 0.0, Acts::PdgParticle::eElectron, massElectron, -1_e);
+const auto centralPositron =
+ makeParticle(0_mm, 0.0, Acts::PdgParticle::ePositron, massElectron, 1_e);
+const auto centralMuon =
+ makeParticle(0_mm, 0.0, Acts::PdgParticle::eMuon, massMuon, -1_e);
+const auto centralAntiMuon =
+ makeParticle(0_mm, 0.0, Acts::PdgParticle::eAntiMuon, massMuon, 1_e);
+const auto backwardPion =
+ makeParticle(-100_mm, -4.5, Acts::PdgParticle::ePionMinus, massPion, -1_e);
+const auto centralPion =
+ makeParticle(0_mm, 0.0, Acts::PdgParticle::ePionMinus, massPion, -1_e);
+const auto forwardPion =
+ makeParticle(100_mm, 4.5, Acts::PdgParticle::ePionMinus, massPion, -1_e);
+const auto centralNeutron =
+ makeParticle(1_mm, 0.1, Acts::PdgParticle::eNeutron, 1_GeV, 0_e);
} // namespace
} // namespace Dataset
diff --git a/Tests/UnitTests/Fatras/Selectors/KinematicCastsTests.cpp b/Tests/UnitTests/Fatras/Selectors/KinematicCastsTests.cpp
index eacafa5c2dca61a8599ae7a8372790b0c51a822e..c427323c670c1f49832bdf036cf879a0eccf4f82 100644
--- a/Tests/UnitTests/Fatras/Selectors/KinematicCastsTests.cpp
+++ b/Tests/UnitTests/Fatras/Selectors/KinematicCastsTests.cpp
@@ -19,31 +19,43 @@ static constexpr double eps = 1e-10;
BOOST_AUTO_TEST_SUITE(FatrasKinematicCasts)
+BOOST_AUTO_TEST_CASE(BackwardParticle) {
+ const auto& particle = Dataset::backwardPion;
+
+ CHECK_SMALL(Vrho()(particle), eps);
+ CHECK_CLOSE_REL(Vz()(particle), -100_mm, eps);
+ CHECK_CLOSE_REL(AbsVz()(particle), 100_mm, eps);
+ CHECK_CLOSE_REL(Eta()(particle), -4.5, eps);
+ CHECK_CLOSE_REL(AbsEta()(particle), 4.5, eps);
+ CHECK_CLOSE_REL(Pt()(particle), 1.5_GeV / std::cosh(4.5), eps);
+ CHECK_CLOSE_REL(P()(particle), 1.5_GeV, eps);
+ CHECK_CLOSE_REL(E()(particle), std::hypot(1.5_GeV, Dataset::massPion), eps);
+}
+
BOOST_AUTO_TEST_CASE(CentralParticle) {
const auto& particle = Dataset::centralPion;
- CHECK_SMALL(vR()(particle), eps);
- CHECK_SMALL(vZ()(particle), eps);
+ CHECK_SMALL(Vrho()(particle), eps);
+ CHECK_SMALL(Vz()(particle), eps);
CHECK_SMALL(AbsVz()(particle), eps);
- CHECK_SMALL(eta()(particle), eps);
- CHECK_SMALL(absEta()(particle), eps);
- CHECK_CLOSE_REL(pT()(particle), 1.5_GeV, eps);
- CHECK_CLOSE_REL(p()(particle), 1.5_GeV, eps);
- // TODO energy
+ CHECK_SMALL(Eta()(particle), eps);
+ CHECK_SMALL(AbsEta()(particle), eps);
+ CHECK_CLOSE_REL(Pt()(particle), 1.5_GeV, eps);
+ CHECK_CLOSE_REL(P()(particle), 1.5_GeV, eps);
+ CHECK_CLOSE_REL(E()(particle), std::hypot(1.5_GeV, Dataset::massPion), eps);
}
BOOST_AUTO_TEST_CASE(ForwardParticle) {
const auto& particle = Dataset::forwardPion;
- CHECK_SMALL(vR()(particle), eps);
- CHECK_CLOSE_REL(vZ()(particle), 100_mm, eps);
+ CHECK_SMALL(Vrho()(particle), eps);
+ CHECK_CLOSE_REL(Vz()(particle), 100_mm, eps);
CHECK_CLOSE_REL(AbsVz()(particle), 100_mm, eps);
- // TODO use value comparisons instead of relative checks
- BOOST_TEST(4.0 < eta()(particle));
- BOOST_TEST(4.0 < absEta()(particle));
- BOOST_TEST(10_MeV < pT()(particle));
- BOOST_TEST(1.5_GeV < p()(particle));
- // TODO energy
+ CHECK_CLOSE_REL(Eta()(particle), 4.5, eps);
+ CHECK_CLOSE_REL(AbsEta()(particle), 4.5, eps);
+ CHECK_CLOSE_REL(Pt()(particle), 1.5_GeV / std::cosh(4.5), eps);
+ CHECK_CLOSE_REL(P()(particle), 1.5_GeV, eps);
+ CHECK_CLOSE_REL(E()(particle), std::hypot(1.5_GeV, Dataset::massPion), eps);
}
BOOST_AUTO_TEST_SUITE_END()
diff --git a/Tests/UnitTests/Fatras/Selectors/LimitSelectorsTests.cpp b/Tests/UnitTests/Fatras/Selectors/LimitSelectorsTests.cpp
index 111148198e7be44efb9f11c09db01b287190e759..f032ee9bb5be29061c6651544e360c901333b513 100644
--- a/Tests/UnitTests/Fatras/Selectors/LimitSelectorsTests.cpp
+++ b/Tests/UnitTests/Fatras/Selectors/LimitSelectorsTests.cpp
@@ -16,21 +16,35 @@ using namespace Acts::UnitLiterals;
BOOST_AUTO_TEST_SUITE(LimitSelectors)
-BOOST_AUTO_TEST_CASE(X0L0Limits) {
- ActsFatras::X0Limit selectX0;
- ActsFatras::L0Limit selectL0;
-
+// Construct a particle that is close to its X0/L0 path limit.
+//
+// Passing a thin slab should still be Ok, but the thick slab should not.
+static ActsFatras::Particle makeParticleCloseToLimit() {
// create particle and move it close to the X0/L0 limit
auto particle = Dataset::centralPion;
- particle.setLimits(0.15, 0.45);
- particle.update(particle.position(), particle.momentum(), 0.10, 0.34);
+ particle.addPassedMaterial(0.125, 0.0125);
+ particle.setMaterialLimits(
+ 0.125 + 1.125 * Dataset::thinSlab.thicknessInX0(),
+ 0.0125 + 1.125 * Dataset::thinSlab.thicknessInL0());
+ return particle;
+}
+
+BOOST_AUTO_TEST_CASE(PathLimitX0) {
+ ActsFatras::PathLimitX0 select;
+ auto particle = makeParticleCloseToLimit();
+ // particle is still within limits for thin block
+ BOOST_TEST(not select(Dataset::thinSlab, particle));
+ // particle would pass limits for thick block
+ BOOST_TEST(select(Dataset::thickSlab, particle));
+}
+BOOST_AUTO_TEST_CASE(PathLimitL0) {
+ ActsFatras::PathLimitL0 select;
+ auto particle = makeParticleCloseToLimit();
// particle is still within limits for thin block
- BOOST_TEST(not selectX0(Dataset::thinSlab, particle));
- BOOST_TEST(not selectL0(Dataset::thinSlab, particle));
+ BOOST_TEST(not select(Dataset::thinSlab, particle));
// particle would pass limits for thick block
- BOOST_TEST(selectX0(Dataset::thickSlab, particle));
- BOOST_TEST(selectL0(Dataset::thickSlab, particle));
+ BOOST_TEST(select(Dataset::thickSlab, particle));
}
BOOST_AUTO_TEST_SUITE_END()
diff --git a/Tests/UnitTests/Fatras/Selectors/SelectorHelpersTests.cpp b/Tests/UnitTests/Fatras/Selectors/SelectorHelpersTests.cpp
index c43ea4505817486cdcb17ac6b8818ab4641e77b1..aa4a859d3e8e93edf1678b5b31f7f1e385bca5a1 100644
--- a/Tests/UnitTests/Fatras/Selectors/SelectorHelpersTests.cpp
+++ b/Tests/UnitTests/Fatras/Selectors/SelectorHelpersTests.cpp
@@ -23,17 +23,13 @@ BOOST_AUTO_TEST_SUITE(FatrasSelectorHelpers)
BOOST_AUTO_TEST_CASE(Min) {
// require a minimum eta value of 0.5
- ActsFatras::Min<ActsFatras::Casts::eta> minEta;
- minEta.valMin = 0.5;
-
+ ActsFatras::Min<ActsFatras::Casts::Eta> minEta{0.5};
BOOST_TEST(not minEta(detector, backward));
BOOST_TEST(not minEta(detector, central));
BOOST_TEST(minEta(detector, forward));
// require a mininum absolute eta value of 0.5
- ActsFatras::Min<ActsFatras::Casts::absEta> minAbsEta;
- minAbsEta.valMin = 0.5;
-
+ ActsFatras::Min<ActsFatras::Casts::AbsEta> minAbsEta{0.5};
BOOST_TEST(minAbsEta(detector, backward));
BOOST_TEST(not minAbsEta(detector, central));
BOOST_TEST(minAbsEta(detector, forward));
@@ -41,35 +37,25 @@ BOOST_AUTO_TEST_CASE(Min) {
BOOST_AUTO_TEST_CASE(Max) {
// require a maximum eta value of 0.5
- ActsFatras::Max<ActsFatras::Casts::eta> maxEta;
- maxEta.valMax = 0.5;
-
+ ActsFatras::Max<ActsFatras::Casts::Eta> maxEta{0.5};
BOOST_TEST(maxEta(detector, backward));
BOOST_TEST(maxEta(detector, central));
BOOST_TEST(not maxEta(detector, forward));
// require a maximum absolute eta value of 0.5
- ActsFatras::Max<ActsFatras::Casts::absEta> maxAbsEta;
- maxAbsEta.valMax = 0.5;
-
+ ActsFatras::Max<ActsFatras::Casts::AbsEta> maxAbsEta{0.5};
BOOST_TEST(not maxAbsEta(detector, backward));
BOOST_TEST(maxAbsEta(detector, central));
BOOST_TEST(not maxAbsEta(detector, forward));
}
BOOST_AUTO_TEST_CASE(Range) {
- ActsFatras::Range<ActsFatras::Casts::eta> rangeEta;
- rangeEta.valMin = -6.;
- rangeEta.valMax = -0.5;
-
+ ActsFatras::Range<ActsFatras::Casts::Eta> rangeEta{-6.0, -0.5};
BOOST_TEST(rangeEta(detector, backward));
BOOST_TEST(not rangeEta(detector, central));
BOOST_TEST(not rangeEta(detector, forward));
- ActsFatras::Range<ActsFatras::Casts::absEta> rangeAbsEta;
- rangeAbsEta.valMin = 0.5;
- rangeAbsEta.valMax = 6.0;
-
+ ActsFatras::Range<ActsFatras::Casts::AbsEta> rangeAbsEta{0.5, 6.0};
BOOST_TEST(rangeAbsEta(detector, backward));
BOOST_TEST(not rangeAbsEta(detector, central));
BOOST_TEST(rangeAbsEta(detector, forward));