Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

fix: Handle Fatras immediate abort edge case #3744

Merged
merged 3 commits into from
Oct 17, 2024
Merged

fix: Handle Fatras immediate abort edge case #3744

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
1ad4fb0
fix fatras
andiwand Oct 16, 2024
2c3980f
fix unit test
andiwand Oct 16, 2024
3e436af
Merge branch 'main' into fix-fatras-actor-edge-case
kodiakhq[bot] Oct 17, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 9 additions & 1 deletion Fatras/include/ActsFatras/Kernel/Simulation.hpp
View file
Open in desktop
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));
andiwand marked this conversation as resolved.
Show resolved Hide resolved

// 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
Open in desktop
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
Open in desktop
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
Loading