Fix mcstack per event indexmap

DataFormats/simulation/include/SimulationDataFormat/MCCompLabel.h

@@ -52,7 +52,14 @@ class MCCompLabel
   // mask for all used fields
   static constexpr uint64_t maskFull = (ul0x1 << (nbitsTrackID + nbitsEvID + nbitsSrcID)) - 1;
 
-  MCCompLabel(int trackID, int evID, int srcID, bool fake = false) { set(trackID, evID, srcID, fake); }
+  MCCompLabel(int trackID, int evID, int srcID, bool fake = false)
+  {
+    // a negative trackID means no MC particle is attached to this signal;
+    // the label stays unset rather than encoding a track that does not exist
+    if (trackID >= 0) {
+      set(trackID, evID, srcID, fake);
+    }
+  }
   MCCompLabel(bool noise = false)
   {
     if (noise) {

DataFormats/simulation/test/testMCCompLabel.cxx

@@ -53,3 +53,10 @@ BOOST_AUTO_TEST_CASE(MCCompLabel_test)
   MCCompLabel dummy;
   BOOST_CHECK(dummy.isEmpty() && !dummy.isNoise() && dummy.isFake() && !dummy.isValid());
 }
+
+// A hit whose track was pruned carries trackID -1 and has no MC particle
+BOOST_AUTO_TEST_CASE(MCCompLabel_no_particle_test)
+{
+  MCCompLabel unmapped(-1, 200, 10);
+  BOOST_CHECK(!unmapped.isValid());
+}

Detectors/Base/include/DetectorsBase/Detector.h

@@ -163,6 +163,25 @@ class Detector : public FairDetector
   // FIXME: make private friend of stack?
   virtual void updateHitTrackIndices(std::map<int, int> const&) = 0;
 
+  // the index a hit's track ended up at after the stack filtered the event;
+  // a pruned track has no entry in the mapping and gets the invalid index -1
+  static int updatedTrackIndex(std::map<int, int> const& indexmapping, int trackID)
+  {
+    const auto iter = indexmapping.find(trackID);
+    return iter != indexmapping.end() ? iter->second : -1;
+  }
+
+  // the index a hit or a track reference gets when the sub-events of one event
+  // are concatenated; an index already flagged invalid carries no track to
+  // shift and stays invalid
+  static int offsetTrackIndex(int trackID, int nprimaries, int primaryOffset, int secondaryOffset)
+  {
+    if (trackID < 0) {
+      return trackID;
+    }
+    return trackID + (trackID < nprimaries ? primaryOffset : secondaryOffset);
+  }
+
   // interfaces to attach properly encoded hit information to a FairMQ message
   // and to decode it
   virtual void attachHits(fair::mq::Channel&, fair::mq::Parts&) = 0;
@@ -323,8 +342,7 @@ class DetImpl : public o2::base::Detector
                    // them via a probe integer until we get a nullptr
     while (auto hits = static_cast<Det*>(this)->Det::getHits(probe++)) {
       for (auto& hit : *hits) {
-        auto iter = indexmapping.find(hit.GetTrackID());
-        hit.SetTrackID(iter->second);
+        hit.SetTrackID(updatedTrackIndex(indexmapping, hit.GetTrackID()));
       }
     }
   }
@@ -391,10 +409,7 @@ class DetImpl : public o2::base::Detector
           if (incomingdata) {
             // fix the trackIDs for this data
             for (auto& hit : *incomingdata) {
-              const auto oldID = hit.GetTrackID();
-              // offset depends on whether the trackis a primary or secondary
-              Int_t offset = (oldID < nprim) ? idelta0 : idelta1;
-              hit.SetTrackID(oldID + offset);
+              hit.SetTrackID(offsetTrackIndex(hit.GetTrackID(), nprim, idelta0, idelta1));
             }
             // this could be further generalized by using a policy for T
             std::copy(incomingdata->begin(), incomingdata->end(), std::back_inserter(*targetdata));
@@ -458,10 +473,7 @@ class DetImpl : public o2::base::Detector
         if (incomingdata) {
           // fix the trackIDs for this data
           for (auto& hit : *incomingdata) {
-            const auto oldID = hit.GetTrackID();
-            // offset depends on whether the trackis a primary or secondary
-            int offset = (oldID < nprim) ? idelta0 : idelta1;
-            hit.SetTrackID(oldID + offset);
+            hit.SetTrackID(offsetTrackIndex(hit.GetTrackID(), nprim, idelta0, idelta1));
           }
           // this could be further generalized by using a policy for T
           std::copy(incomingdata->begin(), incomingdata->end(), std::back_inserter(*targetdata));

Detectors/Base/src/Stack.cxx

@@ -193,8 +193,8 @@
                       TMCProcess proc2)
 {
   //  printf("Pushing %s toBeDone %5d parentId %5d pdgCode %5d is %5d entries %5d \n",
   //	 proc == kPPrimary ? "Primary:   " : "Secondary: ",
   //	 toBeDone, parentId, pdgCode, is, mNumberOfEntriesInParticles);
 
   //
   // This method is called
@@ -605,6 +605,7 @@
   mPrimaryParticles.clear();
   mTrackRefs->clear();
   mTrackIDtoParticlesEntry.clear();
+  mIndexMap.clear();
   mHitCounter = 0;
 }
 

Detectors/Base/test/testStack.cxx

@@ -13,9 +13,15 @@
 #define BOOST_TEST_MAIN
 #define BOOST_TEST_DYN_LINK
 #include <boost/test/unit_test.hpp>
+#include "DetectorsBase/Detector.h"
 #include "DetectorsBase/Stack.h"
+#include "SimulationDataFormat/BaseHits.h"
 #include "TFile.h"
 #include "TMCProcess.h"
+#include "TRefArray.h"
+#include <map>
+#include <string>
+#include <vector>
 
 using namespace o2;
 
@@ -88,3 +94,89 @@ BOOST_AUTO_TEST_CASE(Stack_isFromRadDecay_test)
   BOOST_CHECK(!st.isFromRadDecay(-1));
   BOOST_CHECK(!st.isFromRadDecay(1000000000));
 }
+
+namespace
+{
+// A test detector to exercise hit creation and its interaction with the MCStack
+class TestDetector : public o2::base::Detector
+{
+ public:
+  // the name is turned into a DetID, so it has to be one of the real detectors
+  TestDetector() : o2::base::Detector("ITS", true) {}
+
+  void updateHitTrackIndices(std::map<int, int> const& indexmapping) override
+  {
+    for (auto& hit : mHits) {
+      hit.SetTrackID(updatedTrackIndex(indexmapping, hit.GetTrackID()));
+    }
+  }
+
+  std::vector<o2::BaseHit> mHits;
+
+  // rest of the interface, unused here
+  std::string getHitBranchNames(int) const override { return {}; }
+  void attachHits(fair::mq::Channel&, fair::mq::Parts&) override {}
+  void fillHitBranch(TTree&, fair::mq::Parts&, int&) override {}
+  void collectHits(int, fair::mq::Parts&, int&) override {}
+  void mergeHitEntriesAndFlush(int, TTree&, std::vector<int> const&, std::vector<int> const&,
+                               std::vector<int> const&) override {}
+  void mergeHitEntries(TTree&, TTree&, std::vector<int> const&, std::vector<int> const&,
+                       std::vector<int> const&) override {}
+  void InitializeO2Detector() override {}
+  void initializeLate() override {}
+  Bool_t ProcessHits(FairVolume* = nullptr) override { return kFALSE; }
+  void Register() override {}
+  void Reset() override {}
+  void ConstructGeometry() override {}
+};
+
+// Transport one primary with n secondaries, so that the stack builds its mapping
+void transportOnePrimary(o2::data::Stack& st, int nsecondaries)
+{
+  int ntr = 0;
+  st.PushTrack(1, -1, 11, 0., 0., 1., 1., 0., 0., 0., 0., 0., 0., 0., kPPrimary, ntr, 1., 1);
+  st.SetCurrentTrack(0);
+  for (int i = 0; i < nsecondaries; ++i) {
+    st.PushTrack(1, 0, 11, 0., 0., 0.1, 0.1, 0., 0., 0., 0., 0., 0., 0., kPHadronic, ntr, 1., 1);
+  }
+  st.FinishPrimary();
+}
+} // namespace
+
+// A pruned track has no entry in the mapping
+BOOST_AUTO_TEST_CASE(Unmapped_trackID_yields_invalid_index)
+{
+  const std::map<int, int> indexmapping{{0, 0}, {1, 1}};
+
+  BOOST_CHECK_EQUAL(o2::base::Detector::updatedTrackIndex(indexmapping, 1), 1);
+  BOOST_CHECK_EQUAL(o2::base::Detector::updatedTrackIndex(indexmapping, 99), -1);
+}
+
+// The mapping is per event and must not survive Reset()
+BOOST_AUTO_TEST_CASE(Stack_does_not_reuse_index_map_of_previous_event)
+{
+  TestDetector det;
+  TRefArray detlist;
+  detlist.Add(&det);
+
+  o2::data::Stack st;
+  transportOnePrimary(st, 20); // event 1: trackIDs 0 to 20
+  st.UpdateTrackIndex(&detlist);
+  st.Reset();
+
+  transportOnePrimary(st, 1); // event 2: trackIDs 0 and 1 only
+  det.mHits.emplace_back(15); // only valid in event 1
+  st.UpdateTrackIndex(&detlist);
+
+  BOOST_CHECK_EQUAL(det.mHits[0].GetTrackID(), -1);
+}
+
+// An invalid index must not be offset when sub-events are merged
+BOOST_AUTO_TEST_CASE(Offsetting_keeps_an_invalid_index_invalid)
+{
+  const int nprimaries = 5, primaryOffset = 10, secondaryOffset = 100;
+
+  BOOST_CHECK_EQUAL(o2::base::Detector::offsetTrackIndex(3, nprimaries, primaryOffset, secondaryOffset), 13);
+  BOOST_CHECK_EQUAL(o2::base::Detector::offsetTrackIndex(7, nprimaries, primaryOffset, secondaryOffset), 107);
+  BOOST_CHECK_EQUAL(o2::base::Detector::offsetTrackIndex(-1, nprimaries, primaryOffset, secondaryOffset), -1);
+}

run/O2HitMerger.h

@@ -627,9 +627,7 @@ class O2HitMerger : public fair::mq::Device
 
   void updateTrackIdWithOffset(TrackReference& ref, Int_t nprim, Int_t idelta0, Int_t idelta1)
   {
-    Int_t cId = ref.getTrackID();
-    Int_t ioffset = (cId < nprim) ? idelta0 : idelta1;
-    ref.setTrackID(cId + ioffset);
+    ref.setTrackID(o2::base::Detector::offsetTrackIndex(ref.getTrackID(), nprim, idelta0, idelta1));
   }
 
   void initHitTreeAndOutFile(std::string prefix, int detID)