refactor: Use particle outcome for alive check

Fatras/include/ActsFatras/EventData/Particle.hpp

@@ -78,6 +78,16 @@ class Particle {
     return p;
   }
 
+  /// Kill the particle by setting the outcome to a non-alive status.
+  /// @param outcome The outcome status to set for the killed particle (must not be Alive)
+  /// @throws std::invalid_argument if the provided outcome is Alive
+  void killParticle(SimulationOutcome outcome) {
+    if (outcome == SimulationOutcome::Alive) {
+      throw std::invalid_argument("Cannot kill particle with outcome 'Alive'.");
+    }
+    m_outcome = outcome;
+  }
+
   /// Set the process type that generated this particle.
   /// @param proc Process type that generated this particle
   /// @return Reference to this particle for method chaining
@@ -165,22 +175,33 @@ class Particle {
   /// Set the absolute momentum.
   /// @param absMomentum Absolute momentum magnitude
   /// @return Reference to this particle for method chaining
+  /// @throws std::invalid_argument if absMomentum is negative
   Particle &setAbsoluteMomentum(double absMomentum) {
+    if (absMomentum < 0) {
+      throw std::invalid_argument("Absolute momentum cannot be negative.");
+    }
     m_absMomentum = absMomentum;
     return *this;
   }
 
-  /// Change the energy by the given amount.
-  ///
-  /// 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.
-  /// @param delta Energy change (negative for energy loss)
+  /// Reduce the energy by the given amount. If the energy loss exceeds the
+  /// current energy, the particle is killed with the specified outcome. If
+  /// stopping the particle was not expected (stoppedOutcome is Alive), an
+  /// exception is thrown.
+  /// @param delta The energy loss amount to subtract from the current energy
+  /// @param stoppedOutcome The outcome to set if the energy loss exceeds the current energy
   /// @return Reference to this particle for method chaining
-  Particle &correctEnergy(double delta) {
-    const auto newEnergy = std::hypot(m_mass, m_absMomentum) + delta;
+  /// @throws std::invalid_argument if the energy loss exceeds the current energy and stoppedOutcome is Alive
+  Particle &loseEnergy(double delta, SimulationOutcome stoppedOutcome =
+                                         SimulationOutcome::Alive) {
+    const double newEnergy = energy() - delta;
     if (newEnergy <= m_mass) {
-      m_absMomentum = 0.;
+      if (stoppedOutcome == SimulationOutcome::Alive) {
+        throw std::invalid_argument(
+            "Energy loss cannot exceed the current energy of the particle if "
+            "the particle is to remain alive.");
+      }
+      killParticle(stoppedOutcome);
     } else {
       m_absMomentum = Acts::fastCathetus(newEnergy, m_mass);
     }
@@ -270,7 +291,7 @@ class Particle {
 
   /// Check if the particle is alive, i.e. is not at rest.
   /// @return True if particle has non-zero momentum, false otherwise
-  bool isAlive() const { return 0. < m_absMomentum; }
+  bool isAlive() const { return m_outcome == SimulationOutcome::Alive; }
 
   /// Check if this is a secondary particle.
   /// @return True if particle is a secondary (has non-zero vertex secondary, generation, or sub-particle), false otherwise
@@ -388,17 +409,17 @@ class Particle {
   /// PDG particle number.
   Acts::PdgParticle m_pdg = Acts::PdgParticle::eInvalid;
   // Particle charge and mass.
-  double m_charge = 0.;
-  double m_mass = 0.;
+  double m_charge = 0;
+  double m_mass = 0;
   // kinematics, i.e. things that change over the particle lifetime.
   Acts::Vector3 m_direction = Acts::Vector3::UnitZ();
-  double m_absMomentum = 0.;
+  double m_absMomentum = 0;
   Acts::Vector4 m_position4 = Acts::Vector4::Zero();
   /// proper time in the particle rest frame
-  double m_properTime = 0.;
+  double m_properTime = 0;
   // accumulated material
-  double m_pathInX0 = 0.;
-  double m_pathInL0 = 0.;
+  double m_pathInX0 = 0;
+  double m_pathInL0 = 0;
   /// number of hits
   std::uint32_t m_numberOfHits = 0;
   /// reference surface

Fatras/include/ActsFatras/Physics/ElectroMagnetic/BetheBloch.hpp

@@ -59,7 +59,7 @@ struct BetheBloch {
     const double loss = lossDistribution(generator);
 
     // Apply the energy loss
-    particle.correctEnergy(-loss);
+    particle.loseEnergy(loss, SimulationOutcome::KilledInteraction);
 
     // Generates no new particles
     return {};

Fatras/include/ActsFatras/Physics/ElectroMagnetic/BetheHeitler.hpp

@@ -76,7 +76,8 @@ struct BetheHeitler {
                           photon.energy() * photon.direction());
 
     // apply the energy loss
-    particle.correctEnergy(-sampledEnergyLoss);
+    particle.loseEnergy(sampledEnergyLoss,
+                        SimulationOutcome::KilledInteraction);
 
     return {photon};
   }

Python/Examples/tests/root_file_hashes.txt

@@ -7,7 +7,7 @@ test_seeding__estimatedparams.root: 18375210d597e956ff4c37d7c424e409522008c5eedb
 test_seeding__performance_seeding.root: 992f9c611d30dde0d3f3ab676bab19ada61ab6a4442828e27b65ec5e5b7a2880
 test_seeding__particles.root: a0744d4373572ff91bf85d2f52ca18421c092325bb3ec1cf06f17211da5c649b
 test_seeding__particles_simulation.root: 7fed798e6219f6c33698885ce5b5cb736b5f760f393d03378efb89666f018611
-test_hashing_seeding__estimatedparams.root: 24ae339846e5bba75fadd68c4b05a20f4d06dcf5d15eb5e9219a3d28983eca18
+test_hashing_seeding__estimatedparams.root: ee85052f5b178f1ea0ac2c64abf93ad8660d79aca448678ddae00d6a68325ef9
 test_seeding_orthogonal__estimatedparams.root: 18375210d597e956ff4c37d7c424e409522008c5eedb94189fb51342020a6d2a
 test_seeding_orthogonal__performance_seeding.root: 60fbedcf5cb2b37cd8e526251940564432890d3a159d231ed819e915a904682c
 test_seeding_orthogonal__particles.root: a0744d4373572ff91bf85d2f52ca18421c092325bb3ec1cf06f17211da5c649b
@@ -49,14 +49,14 @@ test_ckf_tracks_example[generic-truth_estimated]__tracksummary_ckf.root: 30a84c6
 test_ckf_tracks_example[generic-truth_estimated]__performance_seeding.root: 1facb05c066221f6361b61f015cdf0918e94d9f3fce2269ec7b6a4dffeb2bc7e
 test_ckf_tracks_example[generic-truth_smeared]__trackstates_ckf.root: 4249f9a792262cd283ca33862ce6ebdc112af156443b81064cfed913a35ee7d7
 test_ckf_tracks_example[generic-truth_smeared]__tracksummary_ckf.root: 8339d60e8669c090ed554ce720bf10640ecd017f6e09df7c95c6cc0f4ca1dca4
-test_ckf_tracks_example[odd-full_seeding]__trackstates_ckf.root: 087e89cb8961ae3b778bc7aed630f6fff45da88712854e5871b4eb04b5bf1669
-test_ckf_tracks_example[odd-full_seeding]__tracksummary_ckf.root: 8cafcdb0da421e9947e1eecbf6231b2951490c635d7fd147f25a7f0229fe6835
+test_ckf_tracks_example[odd-full_seeding]__trackstates_ckf.root: 904ef7c85108abf42e70301eaa79b483f5b210e6966ac740b96252f19e67ff32
+test_ckf_tracks_example[odd-full_seeding]__tracksummary_ckf.root: 30114df40f3304e4bd6dd4335e28d9e09fa6168647614f8e9d547e34bda1dc12
 test_ckf_tracks_example[odd-full_seeding]__performance_seeding_trees.root: 43c58577aafe07645e5660c4f43904efadf91d8cda45c5c04c248bbe0f59814f
-test_ckf_tracks_example[odd-truth_estimated]__trackstates_ckf.root: 28648b9d4716af831ae802add2bcd9e49e96c28c4585dd1055db5355ae75ec03
-test_ckf_tracks_example[odd-truth_estimated]__tracksummary_ckf.root: fecfb72945bbb8a26c5ca11ed99fb5c9d788b2350a4c846c43dd274d112c02e4
+test_ckf_tracks_example[odd-truth_estimated]__trackstates_ckf.root: ae12cb24a0f06d0166d4d53202ef2288a55641d87361b6376f5ad6e9593a0c66
+test_ckf_tracks_example[odd-truth_estimated]__tracksummary_ckf.root: ec9132c5ad7cefa2335a8d90a983dcba561b98601a6778387bea1b6f8b5f517a
 test_ckf_tracks_example[odd-truth_estimated]__performance_seeding.root: 1a36b7017e59f1c08602ef3c2cb0483c51df248f112e3780c66594110719c575
-test_ckf_tracks_example[odd-truth_smeared]__trackstates_ckf.root: 6874579122581c6d3f96f13384cd037383f340097d1657d75b7d00196b41980d
-test_ckf_tracks_example[odd-truth_smeared]__tracksummary_ckf.root: 5ec3587cf438358c7dd86b9a926a122b3e188f9be42b0aed7c2a9c5baee395f9
+test_ckf_tracks_example[odd-truth_smeared]__trackstates_ckf.root: 9ad92736b71b86340285b341e3139a23a1ddc4e0e9ecce41835ef4df32fef125
+test_ckf_tracks_example[odd-truth_smeared]__tracksummary_ckf.root: 52a8907722beaa5b9cdd97364c1772a2fea426c56649decb02bbf840bde67d23
 test_vertex_fitting_reading[Truth-False-100]__performance_vertexing.root: 76ef6084d758dfdfc0151ddec2170e12d73394424e3dac4ffe46f0f339ec8293
 test_vertex_fitting_reading[Iterative-False-100]__performance_vertexing.root: 60372210c830a04f95ceb78c6c68a9b0de217746ff59e8e73053750c837b57eb
 test_vertex_fitting_reading[Iterative-True-100]__performance_vertexing.root: e34f217d524a5051dbb04a811d3407df3ebe2cc4bb7f54f6bda0847dbd7b52c3

Tests/UnitTests/Fatras/EventData/ParticleTests.cpp

@@ -10,11 +10,12 @@
 
 #include "Acts/Definitions/PdgParticle.hpp"
 #include "Acts/Definitions/Units.hpp"
+#include "Acts/Utilities/MathHelpers.hpp"
 #include "ActsFatras/EventData/Barcode.hpp"
 #include "ActsFatras/EventData/Particle.hpp"
+#include "ActsFatras/EventData/SimulationOutcome.hpp"
 #include "ActsTests/CommonHelpers/FloatComparisons.hpp"
 
-#include <cmath>
 #include <limits>
 
 using Acts::PdgParticle;
@@ -51,8 +52,6 @@ BOOST_AUTO_TEST_CASE(Construct) {
   CHECK_CLOSE_REL(particle.direction().norm(), 1, eps);
   BOOST_CHECK_EQUAL(particle.transverseMomentum(), 0.);
   BOOST_CHECK_EQUAL(particle.absoluteMomentum(), 0.);
-  // particle is created at rest and thus not alive
-  BOOST_CHECK(!particle.isAlive());
 }
 
 BOOST_AUTO_TEST_CASE(CorrectEnergy) {
@@ -66,48 +65,42 @@ BOOST_AUTO_TEST_CASE(CorrectEnergy) {
   BOOST_CHECK_EQUAL(particle.fourMomentum().x(), 2_GeV);
   BOOST_CHECK_EQUAL(particle.fourMomentum().y(), 0_GeV);
   BOOST_CHECK_EQUAL(particle.fourMomentum().z(), 0_GeV);
-  BOOST_CHECK_EQUAL(particle.fourMomentum().w(), std::hypot(1_GeV, 2_GeV));
+  BOOST_CHECK_EQUAL(particle.fourMomentum().w(), Acts::fastHypot(1_GeV, 2_GeV));
   BOOST_CHECK_EQUAL(particle.transverseMomentum(), 2_GeV);
   BOOST_CHECK_EQUAL(particle.absoluteMomentum(), 2_GeV);
-  BOOST_CHECK_EQUAL(particle.energy(), std::hypot(1_GeV, 2_GeV));
+  BOOST_CHECK_EQUAL(particle.energy(), Acts::fastHypot(1_GeV, 2_GeV));
   // particle direction must be normalized
   CHECK_CLOSE_REL(particle.direction().norm(), 1, eps);
 
   // lose some energy
-  particle.correctEnergy(-100_MeV);
+  particle.loseEnergy(100_MeV);
   BOOST_CHECK_LT(particle.transverseMomentum(), 2_GeV);
+  BOOST_CHECK_GT(particle.transverseMomentum(), 0_GeV);
   BOOST_CHECK_LT(particle.absoluteMomentum(), 2_GeV);
-  CHECK_CLOSE_REL(particle.energy(), std::hypot(1_GeV, 2_GeV) - 100_MeV, eps);
+  BOOST_CHECK_GT(particle.absoluteMomentum(), 0_GeV);
+  CHECK_CLOSE_REL(particle.energy(), Acts::fastHypot(1_GeV, 2_GeV) - 100_MeV,
+                  eps);
   CHECK_CLOSE_REL(particle.direction().norm(), 1, eps);
-  // particle is still alive
-  BOOST_CHECK(particle.isAlive());
 
   // lose some more energy
-  particle.correctEnergy(-200_MeV);
+  particle.loseEnergy(200_MeV);
   BOOST_CHECK_LT(particle.transverseMomentum(), 2_GeV);
+  BOOST_CHECK_GT(particle.transverseMomentum(), 0_GeV);
   BOOST_CHECK_LT(particle.absoluteMomentum(), 2_GeV);
-  CHECK_CLOSE_REL(particle.energy(), std::hypot(1_GeV, 2_GeV) - 300_MeV, eps);
+  BOOST_CHECK_GT(particle.absoluteMomentum(), 0_GeV);
+  CHECK_CLOSE_REL(particle.energy(), Acts::fastHypot(1_GeV, 2_GeV) - 300_MeV,
+                  eps);
   CHECK_CLOSE_REL(particle.direction().norm(), 1, eps);
-  // particle is still alive
-  BOOST_CHECK(particle.isAlive());
 
-  // lose a lot of energy
-  particle.correctEnergy(-3_GeV);
-  BOOST_CHECK_EQUAL(particle.transverseMomentum(), 0.);
-  BOOST_CHECK_EQUAL(particle.absoluteMomentum(), 0.);
-  BOOST_CHECK_EQUAL(particle.energy(), particle.mass());
-  CHECK_CLOSE_REL(particle.direction().norm(), 1, eps);
-  // particle is not alive anymore
-  BOOST_CHECK(!particle.isAlive());
+  // lose too much energy without a stopping strategy
+  BOOST_CHECK_THROW(particle.loseEnergy(10_GeV), std::invalid_argument);
 
-  // losing even more energy does nothing
-  particle.correctEnergy(-10_GeV);
-  BOOST_CHECK_EQUAL(particle.transverseMomentum(), 0.);
-  BOOST_CHECK_EQUAL(particle.absoluteMomentum(), 0.);
-  BOOST_CHECK_EQUAL(particle.energy(), particle.mass());
+  // lose too much energy with a stopping strategy
+  particle.loseEnergy(10_GeV, SimulationOutcome::KilledInteraction);
+  BOOST_CHECK_GT(particle.transverseMomentum(), 0.);
+  BOOST_CHECK_GT(particle.absoluteMomentum(), 0.);
   CHECK_CLOSE_REL(particle.direction().norm(), 1, eps);
-  // particle is still not alive
-  BOOST_CHECK(!particle.isAlive());
+  BOOST_CHECK_EQUAL(particle.outcome(), SimulationOutcome::KilledInteraction);
 }
 
 BOOST_AUTO_TEST_SUITE_END()

Tests/UnitTests/Fatras/Kernel/SimulationActorTests.cpp

@@ -77,7 +77,7 @@ struct MockInteractionList {
                      Particle &particle,
                      std::vector<Particle> &generated) const {
     generated.push_back(particle);
-    particle.correctEnergy(-energyLoss);
+    particle.loseEnergy(energyLoss);
     // break if particle is not alive anymore
     return !particle.isAlive();
   }

Tests/UnitTests/Fatras/Physics/BetheBlochTests.cpp

@@ -27,16 +27,22 @@ BOOST_AUTO_TEST_SUITE(PhysicsSuite)
 BOOST_DATA_TEST_CASE(FatrasBetheBloch, Dataset::parameters, pdg, phi, theta, p,
                      seed) {
   Generator gen(seed);
-  ActsFatras::Particle before = Dataset::makeParticle(pdg, phi, theta, p);
+  const ActsFatras::Particle before = Dataset::makeParticle(pdg, phi, theta, p);
   ActsFatras::Particle after = before;
 
   ActsFatras::BetheBloch process;
   const auto outgoing = process(gen, makeUnitSlab(), after);
-  // energy loss changes momentum and energy
-  BOOST_CHECK_LT(after.absoluteMomentum(), before.absoluteMomentum());
-  BOOST_CHECK_LT(after.energy(), before.energy());
+
   // energy loss creates no new particles
   BOOST_CHECK(outgoing.empty());
+  // energy loss changes momentum and energy
+  if (after.isAlive()) {
+    BOOST_CHECK_LT(after.absoluteMomentum(), before.absoluteMomentum());
+    BOOST_CHECK_LT(after.energy(), before.energy());
+  } else {
+    BOOST_CHECK_EQUAL(after.outcome(),
+                      ActsFatras::SimulationOutcome::KilledInteraction);
+  }
 }
 
 BOOST_AUTO_TEST_SUITE_END()