Fix isFromRadDecay ancestry walk and add a unit test

The query added in #15470 gave two wrong answers. Primaries returned true,
because mTrackIDtoParticlesEntry is meaningless for them - primaries never
enter mParticles - so the lookup landed on an unrelated secondary. Descendants
of a radioactive decay returned false, because the `imo > 0` guard skipped
buffer entry 0, which after FinishPrimary() is the first secondary of the
current primary.

Walking trackIDs and stopping at mNumberOfPrimaryParticles removes both, and
the method becomes const, binds a reference instead of copying the MCTrack, and
drops the two includes MCTrack.h already provides. The new test case fails on
all four affected checks without this change.

Author: sawenzel

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 <map>
 #include <memory>
 #include <stack>
@@ -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<int>(mTrackIDtoParticlesEntry.size())) {
-    return false;
-  }
-  const auto entry = mTrackIDtoParticlesEntry[id];
-  if (entry < 0 || entry >= static_cast<int>(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<int>(mTrackIDtoParticlesEntry.size())) {
+      return false;
+    }
+    const auto entry = mTrackIDtoParticlesEntry[id];
+    if (entry < 0 || entry >= static_cast<int>(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

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