From d4d86bca9be7f1f233918073ba83831acca10dbc Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Mon, 10 Aug 2026 13:59:59 +0200 Subject: [PATCH] Make isFromRadDecay walk the mother chain on trackIDs The new query returned the wrong answer for two common cases, both of which come down to how the stack stores its particles during transport. Primaries never enter mParticles - only secondaries do - and mTrackIDtoParticlesEntry is written for every pushed track, so for a primary it points at whichever secondary happened to be next in the buffer. Asking about a primary therefore inspected an unrelated secondary, and returned true whenever that secondary came from a radioactive decay. Since the buffer is emptied after every primary, its first entry is the first secondary of the current primary, so the `imo > 0` loop guard also skipped exactly that particle and lost any radioactive decay recorded there, together with all of its descendants. Walking the chain on trackIDs instead of on buffer entries removes both problems: primaries are the first mNumberOfPrimaryParticles trackIDs, so that single comparison ends the search without a lookup, and every remaining step is a genuine secondary. The method becomes const and takes the parameter by value like its neighbours, the two includes are dropped again because MCTrack.h already provides TMCProcess, and the declaration gains the note that the answer is only meaningful during transport, since selectTracks() rewrites the mother indices afterwards. The accompanying test builds a small stack with a radioactive decay as the first secondary of its primary and checks the direct, indirect and negative cases; it fails on all four of the affected checks before this change. --- Detectors/Base/include/DetectorsBase/Stack.h | 56 +++++++++----------- Detectors/Base/test/testStack.cxx | 44 +++++++++++++++ 2 files changed, 68 insertions(+), 32 deletions(-) diff --git a/Detectors/Base/include/DetectorsBase/Stack.h b/Detectors/Base/include/DetectorsBase/Stack.h index 3845f10c4f5b4..479981a65477a 100644 --- a/Detectors/Base/include/DetectorsBase/Stack.h +++ b/Detectors/Base/include/DetectorsBase/Stack.h @@ -24,8 +24,7 @@ #include "SimulationDataFormat/ParticleStatus.h" #include "Rtypes.h" #include "TParticle.h" -#include "TVirtualMC.h" -#include "TMCProcess.h" + #include #include #include @@ -211,7 +210,10 @@ class Stack : public FairGenericStack /// query if a track is a direct **or** indirect daughter of a parentID /// if trackid is same as parentid it returns true bool isTrackDaughterOf(int /*trackid*/, int /*parentid*/) const; - bool isFromRadDecay(const int id); + /// query if a track originates, directly or indirectly, from a radioactive decay + /// only meaningful during transport, before selectTracks() remaps mother indices + bool isFromRadDecay(int trackid) const; + bool isCurrentTrackDaughterOf(int parentid) const; // returns the index of the currently transported primary @@ -349,39 +351,29 @@ inline int Stack::getMotherTrackId(int trackid) const return mParticles[entryinParticles].getMotherTrackId(); } -inline bool Stack::isFromRadDecay(const int id) +inline bool Stack::isFromRadDecay(int trackid) const { - // Check whether particle originates directly or indirectly from radioactive decay + // Check whether particle originates directly or indirectly from radioactive decay. + // Walks up the mother chain until a primary is reached. Only meaningful during + // transport, since selectTracks() later rewrites the mother indices in mParticles. // - if (id < 0 || id >= static_cast(mTrackIDtoParticlesEntry.size())) { - return false; - } - const auto entry = mTrackIDtoParticlesEntry[id]; - if (entry < 0 || entry >= static_cast(mParticles.size())) - return false; - auto part = (mParticles[entry]); - - // primary particle ? - if (part.getProcess() == 0) - return false; - // particle directly from radioactive decay ? - if (part.getProcess() == kPRadDecay) { - return true; - } - - // search in particle history - auto imo = mTrackIDtoParticlesEntry[part.getMotherTrackId()]; - auto isRad = false; - while (imo > 0) { - auto mother = (mParticles[imo]); - if (mother.getProcess() == kPRadDecay) { - isRad = true; - break; + // Note that primaries are not kept in mParticles and that mTrackIDtoParticlesEntry + // is meaningless for them, so the chain has to be terminated on the trackID itself. + for (int id = trackid; id >= mNumberOfPrimaryParticles;) { + if (id >= static_cast(mTrackIDtoParticlesEntry.size())) { + return false; + } + const auto entry = mTrackIDtoParticlesEntry[id]; + if (entry < 0 || entry >= static_cast(mParticles.size())) { + return false; + } + const auto& part = mParticles[entry]; + if (part.getProcess() == kPRadDecay) { + return true; } - part = mother; - imo = mTrackIDtoParticlesEntry[mother.getMotherTrackId()]; + id = part.getMotherTrackId(); } - return isRad; + return false; } inline bool Stack::isCurrentTrackDaughterOf(int parentid) const diff --git a/Detectors/Base/test/testStack.cxx b/Detectors/Base/test/testStack.cxx index 150fb9515c5f1..6d5349711e437 100644 --- a/Detectors/Base/test/testStack.cxx +++ b/Detectors/Base/test/testStack.cxx @@ -44,3 +44,47 @@ BOOST_AUTO_TEST_CASE(Stack_test) BOOST_CHECK(inst->getPrimaries().size() == 2); } } + +// convenience wrapper to push a track and return the assigned trackID +static int pushTrack(o2::data::Stack& st, int parentId, TMCProcess proc) +{ + int trackId; + st.PushTrack(1, parentId, 0, 0., 0., 0., 10., 5., 5., 5., 0.1, 0., 0., 0., proc, trackId, 1., 1); + return trackId; +} + +// unit test for the radioactive-decay ancestry query +BOOST_AUTO_TEST_CASE(Stack_isFromRadDecay_test) +{ + o2::data::Stack st; + + // two primaries; note that primaries do not enter mParticles, only secondaries do + const auto prim0 = pushTrack(st, -1, kPPrimary); + const auto prim1 = pushTrack(st, -1, kPPrimary); + + // a radioactive decay product of the second primary, and its descendants. + // this is deliberately the *first* secondary of the primary, so that it lands + // in the first entry of the particle buffer + const auto radDecay = pushTrack(st, prim1, kPRadDecay); + const auto radChild = pushTrack(st, radDecay, kPHadronic); + const auto radGrandChild = pushTrack(st, radChild, kPHadronic); + + // a plain secondary of the second primary: no radioactive decay anywhere in its history + const auto ordinary = pushTrack(st, prim1, kPHadronic); + + // primaries can never come from a radioactive decay + BOOST_CHECK(!st.isFromRadDecay(prim0)); + BOOST_CHECK(!st.isFromRadDecay(prim1)); + + // a secondary whose ancestry ends in a primary must terminate the search with false + BOOST_CHECK(!st.isFromRadDecay(ordinary)); + + // directly and indirectly from a radioactive decay + BOOST_CHECK(st.isFromRadDecay(radDecay)); + BOOST_CHECK(st.isFromRadDecay(radChild)); + BOOST_CHECK(st.isFromRadDecay(radGrandChild)); + + // out-of-range track IDs are rejected rather than looked up + BOOST_CHECK(!st.isFromRadDecay(-1)); + BOOST_CHECK(!st.isFromRadDecay(1000000000)); +}