Skip to content

Commit

Permalink
fix: Handle Fatras immediate abort edge case (#3744)
Browse files Browse the repository at this point in the history 
Currently Fatras might abort before reaching the pointer where we initialize the particle. This leads to desynchronization between initial and final particles which will then lead to warning/error logs downstream

```
ParticleSele   WARNING   No final particle found for 4|11|169|0|0
```

The initial call of the actor is now caught be checking the propagation stage which then initializes the result object gracefully.

Discovered in #3742
  • Loading branch information
@andiwand
andiwand authored Oct 17, 2024
1 parent 057d33d commit c862d2d
Show file tree
Hide file tree
Showing 3 changed files with 57 additions and 20 deletions.
10 changes: 9 additions & 1 deletion Fatras/include/ActsFatras/Kernel/Simulation.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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.
Expand Down Expand Up @@ -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.
Expand Down
31 changes: 12 additions & 19 deletions Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand All @@ -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) &&
Expand Down
36 changes: 36 additions & 0 deletions Tests/UnitTests/Fatras/Kernel/SimulationActorTests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -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);
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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);
Expand All @@ -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);
Expand Down Expand Up @@ -455,7 +476,13 @@ BOOST_AUTO_TEST_CASE(NoHitsEmptySurface) {
BOOST_AUTO_TEST_CASE(NoHitsMaterialSurface) {
Fixture<NoSurface> 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);
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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);
Expand All @@ -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());
Expand All @@ -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,
Expand Down

0 comments on commit c862d2d

Please sign in to comment.