Fix instabilities when we have empty timeframes due low interaction rate

DataFormats/simulation/CMakeLists.txt

@@ -55,6 +55,11 @@ o2_target_root_dictionary(
 # * src/SimulationDataLinkDef.h
 # * and not src/SimulationDataFormatLinkDef.h
 
+o2_add_test(DigitizationContext
+            SOURCES test/testDigitizationContext.cxx
+            COMPONENT_NAME SimulationDataFormat
+            PUBLIC_LINK_LIBRARIES O2::SimulationDataFormat)
+
 o2_add_test(InteractionSampler
             SOURCES test/testInteractionSampler.cxx
             COMPONENT_NAME SimulationDataFormat

DataFormats/simulation/include/SimulationDataFormat/DigitizationContext.h

@@ -135,7 +135,11 @@ class DigitizationContext
   void applyMaxCollisionFilter(std::vector<std::tuple<int, int, int>>& timeframeindices, long startOrbit, long orbitsPerTF, int maxColl, double orbitsEarly = 0.);
 
   /// get timeframe structure --> index markers where timeframe starts/ends/is_influenced_by
-  std::vector<std::tuple<int, int, int>> calcTimeframeIndices(long startOrbit, long orbitsPerTF, double orbitsEarly = 0.) const;
+  /// One entry is produced per timeframe, including timeframes which contain no collision at all.
+  /// nTimeframes is the number of timeframes the caller asked for; when given, the result has exactly
+  /// that many entries, so that a timeframe without collisions keeps its own slot instead of shifting
+  /// all later timeframes down by one.
+  std::vector<std::tuple<int, int, int>> calcTimeframeIndices(long startOrbit, long orbitsPerTF, double orbitsEarly = 0., long nTimeframes = -1) const;
 
   // Sample and fix interaction vertices (according to some distribution). Makes sure that same event ids
   // have to have same vertex, as well as event ids associated to same collision.

DataFormats/simulation/src/DigitizationContext.cxx

@@ -389,20 +389,33 @@ void DigitizationContext::fillQED(std::string_view QEDprefix, std::vector<o2::In
 namespace
 {
 // a common helper for timeframe structure
-std::vector<std::pair<int, int>> getTimeFrameBoundaries(std::vector<o2::InteractionTimeRecord> const& irecords, long startOrbit, long orbitsPerTF)
+// One entry is produced per timeframe. A timeframe without collisions gets an empty range
+// (first > second) rather than being left out, so that entry i always describes the timeframe
+// covering orbits [startOrbit + i * orbitsPerTF, startOrbit + (i+1) * orbitsPerTF).
+// nTimeframes, when positive, is the number of timeframes the caller asked for; the result is
+// padded with empty timeframes (or truncated) to exactly that length.
+std::vector<std::pair<int, int>> getTimeFrameBoundaries(std::vector<o2::InteractionTimeRecord> const& irecords, long startOrbit, long orbitsPerTF, long nTimeframes = -1)
 {
   std::vector<std::pair<int, int>> result;
 
+  auto pad_and_return = [&result, nTimeframes](int index) {
+    if (nTimeframes > 0) {
+      while ((long)result.size() < nTimeframes) {
+        result.emplace_back(std::pair<int, int>(index, index - 1)); // an empty timeframe
+      }
+      result.resize(nTimeframes);
+    }
+    return result;
+  };
+
   // the goal is to determine timeframe boundaries inside the interaction record vectors
-  // determine if we can do anything
   if (irecords.size() == 0) {
-    // nothing to do
-    return result;
+    return pad_and_return(0);
   }
 
   if (irecords.back().orbit < startOrbit) {
     LOG(error) << "start orbit larger than last collision entry";
-    return result;
+    return pad_and_return((int)irecords.size());
   }
 
   // skip to the first index falling within our constrained
@@ -413,10 +426,13 @@ std::vector<std::pair<int, int>> getTimeFrameBoundaries(std::vector<o2::Interact
 
   // now we can start (2 pointer approach)
   auto right = left;
-  int timeframe_count = 1;
+  long timeframe_count = 1;
   while (right < irecords.size()) {
-    if (irecords[right].orbit >= startOrbit + timeframe_count * orbitsPerTF) {
-      // we finished one timeframe
+    // a collision may lie several timeframes ahead of the previous one; close every timeframe it
+    // skips over, as an empty one, so that the collision ends up in the timeframe it belongs to.
+    // (A plain "if" here closed only one timeframe per collision, which both dropped the empty
+    //  timeframes and mis-assigned the collisions after them.)
+    while (irecords[right].orbit >= startOrbit + timeframe_count * orbitsPerTF) {
       result.emplace_back(std::pair<int, int>(left, right - 1));
       timeframe_count++;
       left = right;
@@ -425,17 +441,18 @@ std::vector<std::pair<int, int>> getTimeFrameBoundaries(std::vector<o2::Interact
   }
   // finished last timeframe
   result.emplace_back(std::pair<int, int>(left, right - 1));
-  return result;
+  return pad_and_return((int)irecords.size());
 }
 
 // a common helper for timeframe structure - includes indices for orbits-early (orbits from last timeframe still affecting current one)
 std::vector<std::tuple<int, int, int>> getTimeFrameBoundaries(std::vector<o2::InteractionTimeRecord> const& irecords,
                                                               long startOrbit,
                                                               long orbitsPerTF,
-                                                              float orbitsEarly)
+                                                              float orbitsEarly,
+                                                              long nTimeframes = -1)
 {
   // we could actually use the other method first ... then do another pass to fix the early-index ... or impact index
-  auto true_indices = getTimeFrameBoundaries(irecords, startOrbit, orbitsPerTF);
+  auto true_indices = getTimeFrameBoundaries(irecords, startOrbit, orbitsPerTF, nTimeframes);
 
   std::vector<std::tuple<int, int, int>> indices_with_early{};
   for (int ti = 0; ti < true_indices.size(); ++ti) {
@@ -447,7 +464,7 @@ std::vector<std::tuple<int, int, int>> getTimeFrameBoundaries(std::vector<o2::In
 
     // from the second timeframe on we can determine the index in the previous timeframe
     // which matches our criterion
-    if (orbitsEarly > 0. && ti > 0) {
+    if (orbitsEarly > 0. && ti > 0 && tf_range.first <= tf_range.second) {
       auto& prev_tf_range = true_indices[ti - 1];
       // in this range search the smallest index which precedes
       // timeframe ti by not more than "orbitsEarly" orbits
@@ -518,7 +535,8 @@ void DigitizationContext::applyMaxCollisionFilter(std::vector<std::tuple<int, in
 
     LOG(info) << "timeframe indices " << previndex << " : " << firstindex << " : " << lastindex;
 
-    int collCount = 0; // counting collisions within timeframe
+    int collCount = 0;                            // counting collisions within timeframe
+    const size_t nrecords_before = newrecords.size(); // to detect a timeframe that stays empty
     // copy to new structure
     for (int index = previndex >= 0 ? previndex : firstindex; index <= lastindex; ++index) {
       if (collCount >= maxColl) {
@@ -571,6 +589,14 @@ void DigitizationContext::applyMaxCollisionFilter(std::vector<std::tuple<int, in
     } // ends one timeframe
 
     // correct the timeframe indices
+    if (newrecords.size() == nrecords_before) {
+      // this timeframe received no collision at all; give it an empty range at the current
+      // position so that it keeps its slot and the timeframes after it are not shifted
+      std::get<0>(tf_indices) = (int)newrecords.size();
+      std::get<1>(tf_indices) = (int)newrecords.size() - 1;
+      std::get<2>(tf_indices) = -1;
+      continue;
+    }
     if (indices_old_to_new.find(firstindex) != indices_old_to_new.end()) {
       std::get<0>(tf_indices) = indices_old_to_new[firstindex]; // start
     }
@@ -588,9 +614,9 @@ void DigitizationContext::applyMaxCollisionFilter(std::vector<std::tuple<int, in
   mEventParts = newparts;
 }
 
-std::vector<std::tuple<int, int, int>> DigitizationContext::calcTimeframeIndices(long startOrbit, long orbitsPerTF, double orbitsEarly) const
+std::vector<std::tuple<int, int, int>> DigitizationContext::calcTimeframeIndices(long startOrbit, long orbitsPerTF, double orbitsEarly, long nTimeframes) const
 {
-  auto timeframeindices = getTimeFrameBoundaries(mEventRecords, startOrbit, orbitsPerTF, orbitsEarly);
+  auto timeframeindices = getTimeFrameBoundaries(mEventRecords, startOrbit, orbitsPerTF, orbitsEarly, nTimeframes);
   return timeframeindices;
 }
 
@@ -710,6 +736,11 @@ DigitizationContext DigitizationContext::extractSingleTimeframe(int timeframeid,
     if (earlyindex >= 0) {
       startindex = earlyindex;
     }
+    if (endindex < startindex) {
+      // a timeframe without any collision: return a valid but empty context rather than
+      // copying a negative range
+      endindex = startindex;
+    }
     std::copy(mEventRecords.begin() + startindex, mEventRecords.begin() + endindex, std::back_inserter(r.mEventRecords));
     std::copy(mEventParts.begin() + startindex, mEventParts.begin() + endindex, std::back_inserter(r.mEventParts));
     if (mInteractionVertices.size() >= endindex) {

DataFormats/simulation/test/testDigitizationContext.cxx

@@ -0,0 +1,127 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+#define BOOST_TEST_MODULE Test DigitizationContext class
+#define BOOST_TEST_MAIN
+#define BOOST_TEST_DYN_LINK
+
+#include <boost/test/unit_test.hpp>
+#include "SimulationDataFormat/DigitizationContext.h"
+#include <vector>
+
+namespace o2
+{
+
+// build a context whose collisions sit at the given orbits (one collision each, source 0)
+steer::DigitizationContext makeContext(std::vector<long> const& orbits)
+{
+  steer::DigitizationContext ctx;
+  auto& records = ctx.getEventRecords();
+  auto& parts = ctx.getEventParts();
+  int entry = 0;
+  for (auto o : orbits) {
+    records.emplace_back(o2::InteractionTimeRecord(o2::InteractionRecord(0, o), 0.));
+    parts.push_back({steer::EventPart(0, entry++)});
+  }
+  ctx.setNCollisions(records.size());
+  ctx.setMaxNumberParts(1);
+  return ctx;
+}
+
+// The timeframe index structure must have one entry per timeframe asked for, and entry i must
+// describe exactly the collisions falling into orbits [start + i*orbitsPerTF, start + (i+1)*orbitsPerTF).
+BOOST_AUTO_TEST_CASE(TimeframeIndicesAreSlotAligned)
+{
+  long const orbitsPerTF = 6;
+  long const start = 0;
+  long const nTF = 5; // orbits 0..29
+
+  // timeframe 1 (orbits 6..11) and timeframe 4 (orbits 24..29) hold no collision
+  std::vector<long> orbits{0, 3, 5, 12, 14, 17, 18, 21};
+  auto ctx = makeContext(orbits);
+
+  auto indices = ctx.calcTimeframeIndices(start, orbitsPerTF, 0., nTF);
+  BOOST_CHECK_EQUAL(indices.size(), (size_t)nTF);
+
+  for (int tf = 0; tf < nTF; ++tf) {
+    auto first = std::get<0>(indices[tf]);
+    auto last = std::get<1>(indices[tf]);
+    long const lo = start + tf * orbitsPerTF;
+    long const hi = lo + orbitsPerTF;
+    // count what should be in this timeframe
+    int expected = 0;
+    for (auto o : orbits) {
+      if (o >= lo && o < hi) {
+        expected++;
+      }
+    }
+    BOOST_CHECK_EQUAL(last - first + 1, expected);
+    for (int i = first; i <= last; ++i) {
+      BOOST_CHECK(orbits[i] >= lo);
+      BOOST_CHECK(orbits[i] < hi);
+    }
+  }
+}
+
+// A timeframe without collisions must survive extraction as a valid, empty context
+BOOST_AUTO_TEST_CASE(EmptyTimeframeExtracts)
+{
+  long const orbitsPerTF = 6;
+  long const nTF = 3;
+  auto ctx = makeContext({0, 2, 13});         // timeframe 1 (orbits 6..11) is empty
+  auto indices = ctx.calcTimeframeIndices(0, orbitsPerTF, 0., nTF);
+  BOOST_CHECK_EQUAL(indices.size(), (size_t)nTF);
+
+  auto tf0 = ctx.extractSingleTimeframe(0, indices, {});
+  auto tf1 = ctx.extractSingleTimeframe(1, indices, {});
+  auto tf2 = ctx.extractSingleTimeframe(2, indices, {});
+  BOOST_CHECK_EQUAL(tf0.getEventRecords().size(), (size_t)2);
+  BOOST_CHECK_EQUAL(tf1.getEventRecords().size(), (size_t)0);
+  BOOST_CHECK_EQUAL(tf2.getEventRecords().size(), (size_t)1);
+  BOOST_CHECK_EQUAL(tf2.getEventRecords()[0].orbit, 13);
+}
+
+// The trailing timeframes of the requested range must be present even when the last collision
+// falls well before the end of the range
+BOOST_AUTO_TEST_CASE(TrailingTimeframesArePresent)
+{
+  long const orbitsPerTF = 6;
+  long const nTF = 9; // this is what an 8-timeframe anchored MC job with orbitsEarly asks for
+  auto ctx = makeContext({1, 2, 7});
+  auto indices = ctx.calcTimeframeIndices(0, orbitsPerTF, 0., nTF);
+  BOOST_CHECK_EQUAL(indices.size(), (size_t)nTF);
+  for (int tf = 2; tf < nTF; ++tf) {
+    BOOST_CHECK(std::get<0>(indices[tf]) > std::get<1>(indices[tf])); // empty, but present
+  }
+}
+
+// applyMaxCollisionFilter must not shift timeframes when one of them is empty
+BOOST_AUTO_TEST_CASE(MaxCollisionFilterKeepsSlots)
+{
+  long const orbitsPerTF = 6;
+  long const nTF = 4;
+  //  tf0: orbits 0,1,2   tf1: empty   tf2: orbits 12,13   tf3: orbit 19
+  auto ctx = makeContext({0, 1, 2, 12, 13, 19});
+  auto indices = ctx.calcTimeframeIndices(0, orbitsPerTF, 0., nTF);
+  ctx.applyMaxCollisionFilter(indices, 0, orbitsPerTF, 2, 0.); // keep at most 2 per timeframe
+
+  BOOST_CHECK_EQUAL(indices.size(), (size_t)nTF);
+  BOOST_CHECK_EQUAL(std::get<1>(indices[0]) - std::get<0>(indices[0]) + 1, 2); // capped
+  BOOST_CHECK(std::get<0>(indices[1]) > std::get<1>(indices[1]));              // still empty
+  BOOST_CHECK_EQUAL(std::get<1>(indices[2]) - std::get<0>(indices[2]) + 1, 2);
+  BOOST_CHECK_EQUAL(std::get<1>(indices[3]) - std::get<0>(indices[3]) + 1, 1);
+
+  auto tf2 = ctx.extractSingleTimeframe(2, indices, {});
+  BOOST_CHECK_EQUAL(tf2.getEventRecords().size(), (size_t)2);
+  BOOST_CHECK_EQUAL(tf2.getEventRecords()[0].orbit, 12);
+}
+
+} // namespace o2

Detectors/CPV/workflow/src/ClusterReaderSpec.cxx

@@ -41,7 +41,16 @@ void ClusterReader::init(InitContext& ic)
 void ClusterReader::run(ProcessingContext& pc)
 {
   auto ent = mTree->GetReadEntry() + 1;
-  assert(ent < mTree->GetEntries()); // this should not happen
+  if (ent >= mTree->GetEntries()) {
+    // A timeframe holds no collision at all whenever the interaction rate is low enough, and
+    // the tree then has no entry to read. End the stream instead of reading past the end and
+    // publishing branch addresses that GetEntry has not filled. This was an assert, which is
+    // compiled out of every production build since ENABLE_CASSERT defaults to OFF.
+    LOG(info) << "no entry to read, ending the stream";
+    pc.services().get<ControlService>().endOfStream();
+    pc.services().get<ControlService>().readyToQuit(QuitRequest::Me);
+    return;
+  }
   mTree->GetEntry(ent);
   LOG(info) << "Pushing " << mClusters.size() << " Clusters in " << mTRs.size() << " TriggerRecords at entry " << ent;
   pc.outputs().snapshot(Output{mOrigin, "CLUSTERS", 0}, mClusters);

Detectors/CPV/workflow/src/DigitReaderSpec.cxx

@@ -41,7 +41,16 @@ void DigitReader::init(InitContext& ic)
 void DigitReader::run(ProcessingContext& pc)
 {
   auto ent = mTree->GetReadEntry() + 1;
-  assert(ent < mTree->GetEntries()); // this should not happen
+  if (ent >= mTree->GetEntries()) {
+    // A timeframe holds no collision at all whenever the interaction rate is low enough, and
+    // the tree then has no entry to read. End the stream instead of reading past the end and
+    // publishing branch addresses that GetEntry has not filled. This was an assert, which is
+    // compiled out of every production build since ENABLE_CASSERT defaults to OFF.
+    LOG(info) << "no entry to read, ending the stream";
+    pc.services().get<ControlService>().endOfStream();
+    pc.services().get<ControlService>().readyToQuit(QuitRequest::Me);
+    return;
+  }
   mTree->GetEntry(ent);
   LOG(info) << "Pushing " << mDigits.size() << " Digits in " << mTRs.size() << " TriggerRecords at entry " << ent;
   pc.outputs().snapshot(Output{mOrigin, "DIGITS", 0}, mDigits);

Detectors/CTP/workflowIO/src/DigitReaderSpec.cxx

@@ -86,7 +86,16 @@ void DigitReader::run(ProcessingContext& pc)
   auto ent = mTree->GetReadEntry();
   if (!mUseIRFrames) {
     ent++;
-    assert(ent < mTree->GetEntries()); // this should not happen
+    if (ent >= mTree->GetEntries()) {
+      // A timeframe holds no collision at all whenever the interaction rate is low enough, and
+      // the tree then has no entry to read. End the stream instead of reading past the end and
+      // publishing branch addresses that GetEntry has not filled. This was an assert, which is
+      // compiled out of every production build since ENABLE_CASSERT defaults to OFF.
+      LOG(info) << "no entry to read, ending the stream";
+      pc.services().get<ControlService>().endOfStream();
+      pc.services().get<ControlService>().readyToQuit(QuitRequest::Me);
+      return;
+    }
     mTree->GetEntry(ent);
     LOG(info) << "DigitReader pushes " << mDigits.size() << " digits at entry " << ent;
     pc.outputs().snapshot(Output{"CTP", "DIGITS", 0}, mDigits);

Detectors/FIT/FDD/workflow/src/DigitReaderSpec.cxx

@@ -77,7 +77,16 @@ void DigitReader::run(ProcessingContext& pc)
     }
   }
   auto ent = mTree->GetReadEntry() + 1;
-  assert(ent < mTree->GetEntries()); // this should not happen
+  if (ent >= mTree->GetEntries()) {
+    // A timeframe holds no collision at all whenever the interaction rate is low enough, and
+    // the tree then has no entry to read. End the stream instead of reading past the end and
+    // publishing branch addresses that GetEntry has not filled. This was an assert, which is
+    // compiled out of every production build since ENABLE_CASSERT defaults to OFF.
+    LOG(info) << "no entry to read, ending the stream";
+    pc.services().get<ControlService>().endOfStream();
+    pc.services().get<ControlService>().readyToQuit(QuitRequest::Me);
+    return;
+  }
   mTree->GetEntry(ent);
 
   LOG(info) << "FDD DigitReader pushes " << digitsBC->size() << " digits";

Detectors/FIT/FDD/workflow/src/RecPointReaderSpec.cxx

@@ -45,7 +45,16 @@ void RecPointReader::init(InitContext& ic)
 void RecPointReader::run(ProcessingContext& pc)
 {
   auto ent = mTree->GetReadEntry() + 1;
-  assert(ent < mTree->GetEntries()); // this should not happen
+  if (ent >= mTree->GetEntries()) {
+    // A timeframe holds no collision at all whenever the interaction rate is low enough, and
+    // the tree then has no entry to read. End the stream instead of reading past the end and
+    // publishing branch addresses that GetEntry has not filled. This was an assert, which is
+    // compiled out of every production build since ENABLE_CASSERT defaults to OFF.
+    LOG(info) << "no entry to read, ending the stream";
+    pc.services().get<ControlService>().endOfStream();
+    pc.services().get<ControlService>().readyToQuit(QuitRequest::Me);
+    return;
+  }
   mTree->GetEntry(ent);
 
   LOG(info) << "FDD RecPointReader pushes " << mRecPoints->size() << " recpoints with " << mChannelData->size() << " channels at entry " << ent;

Detectors/FIT/FT0/workflow/src/DigitReaderSpec.cxx

@@ -61,7 +61,16 @@ void DigitReader::run(ProcessingContext& pc)
     mTree->SetBranchAddress("FT0DIGITSMCTR", &plabels);
   }
   auto ent = mTree->GetReadEntry() + 1;
-  assert(ent < mTree->GetEntries()); // this should not happen
+  if (ent >= mTree->GetEntries()) {
+    // A timeframe holds no collision at all whenever the interaction rate is low enough, and
+    // the tree then has no entry to read. End the stream instead of reading past the end and
+    // publishing branch addresses that GetEntry has not filled. This was an assert, which is
+    // compiled out of every production build since ENABLE_CASSERT defaults to OFF.
+    LOG(info) << "no entry to read, ending the stream";
+    pc.services().get<ControlService>().endOfStream();
+    pc.services().get<ControlService>().readyToQuit(QuitRequest::Me);
+    return;
+  }
   mTree->GetEntry(ent);
   LOG(debug) << "FT0DigitReader pushed " << channels.size() << " channels in " << digits.size() << " digits";
   pc.outputs().snapshot(Output{"FT0", "DIGITSBC", 0}, digits);