Fix in indexing for Vertex-Track matching

DataFormats/Reconstruction/include/ReconstructionDataFormats/VtxTrackRef.h

@@ -27,17 +27,26 @@ namespace o2
 namespace dataformats
 {
 
+/* Class to refer in start and number of contributors in the container of with consecutively filled conributors.
+     The contributors are suppossed to be sorted according to their sources.
+     Note: the only way to fill the references is to fill them all in increasing order and set the end!
+     VtxTrackIndex ref;
+     for (int i=0;i<VtxTrackIndex::Source::NSources;i++) {
+       ref.setFirstEntryOfSource(i, idxI); // idxI must be >= idxI-1 (if it is =, then source i has not entries
+     }
+     ref.setEnd(idxLast + 1); // i.e. idxLast+1 = idx0 + TotalNumberOfEntries
+  */
+
 class VtxTrackRef : public RangeReference<int, int>
 {
  public:
-  VtxTrackRef(int ent, int n) : RangeReference(ent, n)
+  VtxTrackRef() : RangeReference(-1, 0)
   {
-    auto end = ent + n;
     for (int i = VtxTrackIndex::Source::NSources - 1; i--;) {
-      mFirstEntrySource[i] = end; // only 1st source (base reference) is filled at constructor level
+      mFirstEntrySource[i] = -1; // only 1st source (base reference) is filled at constructor level
     }
   }
-  using RangeReference<int, int>::RangeReference;
+
   void print() const;
   std::string asString() const;
 
@@ -64,7 +73,12 @@ class VtxTrackRef : public RangeReference<int, int>
     }
   }
 
+  // set the last +1 element index and finalize all references
+  void setEnd(int end);
+
  private:
+  using RangeReference<int, int>::RangeReference;
+
   std::array<int, VtxTrackIndex::Source::NSources - 1> mFirstEntrySource{0};
 
   ClassDefNV(VtxTrackRef, 1);

DataFormats/Reconstruction/src/VtxTrackRef.cxx

@@ -17,6 +17,7 @@
 #include <fmt/printf.h>
 #include <iostream>
 #include <bitset>
+#include <climits>
 
 using namespace o2::dataformats;
 
@@ -29,14 +30,32 @@ std::string VtxTrackRef::asString() const
   return str;
 }
 
+// set the last +1 element index and finalize all references
+void VtxTrackRef::print() const
+{
+  LOG(INFO) << asString();
+}
+
+// set the last +1 element index and check consistency
+void VtxTrackRef::setEnd(int end)
+{
+  if (end <= 0) {
+    return; // empty
+  }
+  setEntries(end - getFirstEntry());
+  for (int i = VtxTrackIndex::NSources - 1; i--;) {
+    if (getFirstEntryOfSource(i) < 0) {
+      throw std::runtime_error(fmt::format("1st entry for source {:d} was not set", i));
+    }
+    if (getEntriesOfSource(i) < 0) {
+      throw std::runtime_error(fmt::format("Source {:d} has negative number of entrie", getEntriesOfSource(i)));
+    }
+  }
+}
+
 std::ostream& o2::dataformats::operator<<(std::ostream& os, const o2::dataformats::VtxTrackRef& v)
 {
   // stream itself
   os << v.asString();
   return os;
 }
-
-void VtxTrackRef::print() const
-{
-  LOG(INFO) << asString();
-}

Detectors/Vertexing/include/DetectorsVertexing/VertexTrackMatcher.h

@@ -42,7 +42,7 @@ class VertexTrackMatcher
   using TrackITS = o2::its::TrackITS;
   using ITSROFR = o2::itsmft::ROFRecord;
   using TrackTPC = o2::tpc::TrackTPC;
-  using TmpMap = std::unordered_map<int, std::vector<GIndex>>;
+  using TmpMap = std::vector<std::vector<GIndex>>;
   using TimeEst = o2::dataformats::TimeStampWithError<float, float>;
   using TBracket = o2::math_utils::Bracketf_t;
 

Detectors/Vertexing/src/VertexTrackMatcher.cxx

@@ -49,21 +49,30 @@ void VertexTrackMatcher::process(const gsl::span<const PVertex>& vertices,
                                  std::vector<VRef>& vtxRefs)
 {
 
-  TmpMap tmpMap;
-  tmpMap.reserve(vertices.size());
-
+  int nv = vertices.size();
+  TmpMap tmpMap(nv);
+  // 1st register vertex contributors
   // TPC/ITS and TPC tracks are not sorted in time, do this and exclude indiced of tracks used in the vertex fit
   std::vector<int> idTPCITS(tpcitsTracks.size()); // indices of TPCITS tracks sorted in time
+  std::iota(idTPCITS.begin(), idTPCITS.end(), 0);
+  for (int iv = 0; iv < nv; iv++) {
+    int idMin = v2tfitRefs[iv].getFirstEntry(), idMax = idMin + v2tfitRefs[iv].getEntries();
+    auto& vtxIds = tmpMap[iv]; // global IDs of contibuting tracks
+    vtxIds.reserve(v2tfitRefs[iv].getEntries());
+    for (int id = idMin; id < idMax; id++) {
+      auto gid = v2tfitIDs[id];
+      vtxIds.push_back(gid);
+      // flag already accounted tracks
+      idTPCITS[gid.getIndex()] = -1; // RS Attention: this will not work once not only ITSTPC contributes to vertex, FIXME!!!
+    }
+  }
+
   std::vector<int> idVtxIRMin(vertices.size());   // indices of vertices sorted in IRmin
   std::vector<int> flgITS(itsTracks.size(), 0);
   std::vector<int> idTPC(tpcTracks.size()); // indices of TPC tracks sorted in min time
   std::iota(idTPC.begin(), idTPC.end(), 0);
-  std::iota(idTPCITS.begin(), idTPCITS.end(), 0);
   std::iota(idVtxIRMin.begin(), idVtxIRMin.end(), 0);
-  for (auto id : v2tfitIDs) { // flag matched ITS-TPC tracks used in the vertex fit, we exclude them from this association
-    idTPCITS[id.getIndex()] = -1;
-  }
-  for (const auto& tpcits : tpcitsTracks) { // flag standalone ITS and TPC tracks used in global matched, we exclude them from association to vertex
+  for (const auto& tpcits : tpcitsTracks) { // flag standalone ITS and TPC tracks used in global matches, we exclude them from association to vertex
     flgITS[tpcits.getRefITS()] = -1;
     idTPC[tpcits.getRefTPC()] = -1;
   }
@@ -87,22 +96,21 @@ void VertexTrackMatcher::process(const gsl::span<const PVertex>& vertices,
     return tI < tJ;
   });
 
-  // Important: do this in the same order in which VtxTrackIndex::Source enums are defined!
   attachTPCITS(tmpMap, tpcitsTracks, idTPCITS, vertices);
   attachITS(tmpMap, itsTracks, itsROFR, flgITS, vertices, idVtxIRMin);
   attachTPC(tmpMap, tpcTimes, idTPC, vertices, idVtxIRMin);
 
   // build vector of global indices
   trackIndex.clear();
   vtxRefs.clear();
+  // just to reuse these vectors for ambiguous attachment counting
   memset(idTPCITS.data(), 0, sizeof(int) * idTPCITS.size());
   memset(idTPC.data(), 0, sizeof(int) * idTPC.size());
   memset(flgITS.data(), 0, sizeof(int) * flgITS.size());
   std::array<std::vector<int>*, GIndex::NSources> vptr;
   vptr[GIndex::ITSTPC] = &idTPCITS;
   vptr[GIndex::ITS] = &flgITS;
   vptr[GIndex::TPC] = &idTPC;
-  int nv = vertices.size();
   // flag tracks attached to >1 vertex
   for (int iv = 0; iv < nv; iv++) {
     const auto& trvec = tmpMap[iv];
@@ -112,40 +120,29 @@ void VertexTrackMatcher::process(const gsl::span<const PVertex>& vertices,
   }
 
   for (int iv = 0; iv < nv; iv++) {
-    int srcStart[GIndex::NSources + 1];
-    int entry = trackIndex.size();
-    srcStart[GIndex::ITSTPC] = entry;
-    for (int is = 1; is < GIndex::NSources; is++) {
-      srcStart[is] = -1;
-    }
-
-    // 1st: attach indices of global tracks used in vertex fit
-    int idMin = v2tfitRefs[iv].getFirstEntry(), idMax = idMin + v2tfitRefs[iv].getEntries();
-    for (int id = idMin; id < idMax; id++) {
-      trackIndex.push_back(v2tfitIDs[id]);
-    }
+    auto& trvec = tmpMap[iv];
+    // sort entries in each vertex track indices list according to the source
+    std::sort(trvec.begin(), trvec.end(), [](GIndex a, GIndex b) { return a.getSource() < b.getSource(); });
 
-    // 2nd: attach non-contributing tracks
-    const auto& trvec = tmpMap[iv];
+    auto entry0 = trackIndex.size();   // start of entries for this vertex
+    auto& vr = vtxRefs.emplace_back(); //entry0, 0);
+    int oldSrc = -1;
     for (const auto gid0 : trvec) {
       int src = gid0.getSource();
-      if (srcStart[src] == -1) {
-        srcStart[src] = trackIndex.size();
+      while (oldSrc < src) {
+        oldSrc++;
+        vr.setFirstEntryOfSource(oldSrc, trackIndex.size()); // register start of new source
       }
       auto& gid = trackIndex.emplace_back(gid0);
       if ((*vptr[src])[gid.getIndex()] > 1) {
         gid.setAmbiguous();
       }
     }
-
-    auto& vr = vtxRefs.emplace_back(entry, trackIndex.size() - entry);
-    srcStart[GIndex::NSources] = trackIndex.size();
-    for (int is = GIndex::NSources; is > 0; is--) {
-      if (srcStart[is] == -1) { // in case the source did not contribute
-        srcStart[is] = srcStart[is + 1];
-      }
-      vr.setFirstEntryOfSource(is, srcStart[is]);
+    while (++oldSrc < GIndex::NSources) {
+      vr.setFirstEntryOfSource(oldSrc, trackIndex.size());
     }
+    vr.setEnd(trackIndex.size());
+    LOG(INFO) << "Vertxex " << iv << " Tracks " << vr;
   }
 }