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)); +}