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.

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