diff --git a/Fatras/include/ActsFatras/Kernel/Simulation.hpp b/Fatras/include/ActsFatras/Kernel/Simulation.hpp index cbc1115225a..2d8280737f0 100644 --- a/Fatras/include/ActsFatras/Kernel/Simulation.hpp +++ b/Fatras/include/ActsFatras/Kernel/Simulation.hpp @@ -245,6 +245,9 @@ struct Simulation { continue; } + assert(result->particle.particleId() == initialParticle.particleId() && + "Particle id must not change during simulation"); + copyOutputs(result.value(), simulatedParticlesInitial, simulatedParticlesFinal, hits); // since physics processes are independent, there can be particle id @@ -256,6 +259,10 @@ struct Simulation { } } + assert( + (simulatedParticlesInitial.size() == simulatedParticlesFinal.size()) && + "Inconsistent final sizes of the simulated particle containers"); + // the overall function call succeeded, i.e. no fatal errors occurred. // yet, there might have been some particle for which the propagation // failed. thus, the successful result contains a list of failed particles. @@ -284,12 +291,13 @@ struct Simulation { // initial particle state was already pushed to the container before // store final particle state at the end of the simulation particlesFinal.push_back(result.particle); + std::copy(result.hits.begin(), result.hits.end(), std::back_inserter(hits)); + // move generated secondaries that should be simulated to the output std::copy_if( result.generatedParticles.begin(), result.generatedParticles.end(), std::back_inserter(particlesInitial), [this](const Particle &particle) { return selectParticle(particle); }); - std::copy(result.hits.begin(), result.hits.end(), std::back_inserter(hits)); } /// Renumber particle ids in the tail of the container. diff --git a/Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp b/Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp index cc6d7200868..acbaaaa713d 100644 --- a/Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp +++ b/Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp @@ -69,7 +69,17 @@ struct SimulationActor { void act(propagator_state_t &state, stepper_t &stepper, navigator_t &navigator, result_type &result, const Acts::Logger &logger) const { - assert(generator && "The generator pointer must be valid"); + assert(generator != nullptr && "The generator pointer must be valid"); + + if (state.stage == Acts::PropagatorStage::prePropagation) { + // first step is special: there is no previous state and we need to arm + // the decay simulation for all future steps. + result.particle = + makeParticle(initialParticle, state, stepper, navigator); + result.properTimeLimit = + decay.generateProperTimeLimit(*generator, initialParticle); + return; + } // actors are called once more after the propagation terminated if (!result.isAlive) { @@ -82,28 +92,11 @@ struct SimulationActor { return; } - // check if we are still on the start surface and skip if so - if ((navigator.startSurface(state.navigation) != nullptr) && - (navigator.startSurface(state.navigation) == - navigator.currentSurface(state.navigation))) { - return; - } - // update the particle state first. this also computes the proper time which // needs the particle state from the previous step for reference. that means // this must happen for every step (not just on surface) and before // everything, e.g. any interactions that could modify the state. - if (std::isnan(result.properTimeLimit)) { - // first step is special: there is no previous state and we need to arm - // the decay simulation for all future steps. - result.particle = - makeParticle(initialParticle, state, stepper, navigator); - result.properTimeLimit = - decay.generateProperTimeLimit(*generator, initialParticle); - } else { - result.particle = - makeParticle(result.particle, state, stepper, navigator); - } + result.particle = makeParticle(result.particle, state, stepper, navigator); // decay check. needs to happen at every step, not just on surfaces. if (std::isfinite(result.properTimeLimit) && diff --git a/Tests/UnitTests/Fatras/Kernel/SimulationActorTests.cpp b/Tests/UnitTests/Fatras/Kernel/SimulationActorTests.cpp index e52d586c30b..3306c5702d6 100644 --- a/Tests/UnitTests/Fatras/Kernel/SimulationActorTests.cpp +++ b/Tests/UnitTests/Fatras/Kernel/SimulationActorTests.cpp @@ -244,7 +244,13 @@ BOOST_AUTO_TEST_CASE(HitsOnEmptySurface) { BOOST_CHECK_EQUAL(f.actor.initialParticle.absoluteMomentum(), f.p); BOOST_CHECK_EQUAL(f.actor.initialParticle.energy(), f.e); + // call.actor: pre propagation + f.state.stage = Acts::PropagatorStage::prePropagation; + f.actor.act(f.state, f.stepper, f.navigator, f.result, + Acts::getDummyLogger()); + // call.actor: surface selection -> one hit, no material -> no secondary + f.state.stage = Acts::PropagatorStage::postStep; f.actor.act(f.state, f.stepper, f.navigator, f.result, Acts::getDummyLogger()); BOOST_CHECK(f.result.isAlive); @@ -271,6 +277,7 @@ BOOST_AUTO_TEST_CASE(HitsOnEmptySurface) { BOOST_CHECK_EQUAL(f.state.stepping.p, f.result.particle.absoluteMomentum()); // call.actor again: one more hit, still no secondary + f.state.stage = Acts::PropagatorStage::postStep; f.actor.act(f.state, f.stepper, f.navigator, f.result, Acts::getDummyLogger()); BOOST_CHECK(f.result.isAlive); @@ -317,7 +324,13 @@ BOOST_AUTO_TEST_CASE(HitsOnMaterialSurface) { BOOST_CHECK_EQUAL(f.actor.initialParticle.absoluteMomentum(), f.p); BOOST_CHECK_EQUAL(f.actor.initialParticle.energy(), f.e); + // call.actor: pre propagation + f.state.stage = Acts::PropagatorStage::prePropagation; + f.actor.act(f.state, f.stepper, f.navigator, f.result, + Acts::getDummyLogger()); + // call.actor: surface selection -> one hit, material -> one secondary + f.state.stage = Acts::PropagatorStage::postStep; f.actor.act(f.state, f.stepper, f.navigator, f.result, Acts::getDummyLogger()); BOOST_CHECK(f.result.isAlive); @@ -346,6 +359,7 @@ BOOST_AUTO_TEST_CASE(HitsOnMaterialSurface) { tol); // call.actor again: one more hit, one more secondary + f.state.stage = Acts::PropagatorStage::postStep; f.actor.act(f.state, f.stepper, f.navigator, f.result, Acts::getDummyLogger()); BOOST_CHECK(f.result.isAlive); @@ -392,7 +406,13 @@ BOOST_AUTO_TEST_CASE(NoHitsEmptySurface) { BOOST_CHECK_EQUAL(f.actor.initialParticle.absoluteMomentum(), f.p); BOOST_CHECK_EQUAL(f.actor.initialParticle.energy(), f.e); + // call.actor: pre propagation + f.state.stage = Acts::PropagatorStage::prePropagation; + f.actor.act(f.state, f.stepper, f.navigator, f.result, + Acts::getDummyLogger()); + // call.actor: no surface sel. -> no hit, no material -> no secondary + f.state.stage = Acts::PropagatorStage::postStep; f.actor.act(f.state, f.stepper, f.navigator, f.result, Acts::getDummyLogger()); BOOST_CHECK(f.result.isAlive); @@ -419,6 +439,7 @@ BOOST_AUTO_TEST_CASE(NoHitsEmptySurface) { BOOST_CHECK_EQUAL(f.state.stepping.p, f.result.particle.absoluteMomentum()); // call.actor again: no hit, still no secondary + f.state.stage = Acts::PropagatorStage::postStep; f.actor.act(f.state, f.stepper, f.navigator, f.result, Acts::getDummyLogger()); BOOST_CHECK(f.result.isAlive); @@ -455,7 +476,13 @@ BOOST_AUTO_TEST_CASE(NoHitsEmptySurface) { BOOST_AUTO_TEST_CASE(NoHitsMaterialSurface) { Fixture f(125_MeV, makeMaterialSurface()); + // call.actor: pre propagation + f.state.stage = Acts::PropagatorStage::prePropagation; + f.actor.act(f.state, f.stepper, f.navigator, f.result, + Acts::getDummyLogger()); + // call.actor: no surface sel. -> no hit, material -> one secondary + f.state.stage = Acts::PropagatorStage::postStep; f.actor.act(f.state, f.stepper, f.navigator, f.result, Acts::getDummyLogger()); BOOST_CHECK(f.result.isAlive); @@ -483,6 +510,7 @@ BOOST_AUTO_TEST_CASE(NoHitsMaterialSurface) { tol); // call.actor again: still no hit, one more secondary + f.state.stage = Acts::PropagatorStage::postStep; f.actor.act(f.state, f.stepper, f.navigator, f.result, Acts::getDummyLogger()); BOOST_CHECK(f.result.isAlive); @@ -523,7 +551,13 @@ BOOST_AUTO_TEST_CASE(Decay) { // inverse Lorentz factor for proper time dilation: 1/gamma = m/E const auto gammaInv = f.m / f.e; + // call.actor: pre propagation + f.state.stage = Acts::PropagatorStage::prePropagation; + f.actor.act(f.state, f.stepper, f.navigator, f.result, + Acts::getDummyLogger()); + // first step w/ defaults leaves particle alive + f.state.stage = Acts::PropagatorStage::postStep; f.actor.act(f.state, f.stepper, f.navigator, f.result, Acts::getDummyLogger()); BOOST_CHECK(f.result.isAlive); @@ -536,6 +570,7 @@ BOOST_AUTO_TEST_CASE(Decay) { BOOST_CHECK_EQUAL(f.result.particle.properTime(), 0_ns); // second step w/ defaults increases proper time + f.state.stage = Acts::PropagatorStage::postStep; f.state.stepping.time += 1_ns; f.actor.act(f.state, f.stepper, f.navigator, f.result, Acts::getDummyLogger()); @@ -549,6 +584,7 @@ BOOST_AUTO_TEST_CASE(Decay) { CHECK_CLOSE_REL(f.result.particle.properTime(), gammaInv * 1_ns, tol); // third step w/ proper time limit decays the particle + f.state.stage = Acts::PropagatorStage::postStep; f.state.stepping.time += 1_ns; f.result.properTimeLimit = f.result.particle.properTime() + gammaInv * 0.5_ns; f.actor.act(f.state, f.stepper, f.navigator, f.result,