Give every timeframe its own slot in the collision context

This fixes a problem in the timeframe index structure of the collision
context and adds a unit test.

- getTimeFrameBoundaries closed only one timeframe per collision, so a
  timeframe without collisions was left out of the index structure
  entirely and the collisions after it were assigned to the wrong
  timeframe.
- The number of extracted per-timeframe contexts was therefore the number
  of non-empty timeframes, not the number of timeframes asked for, and
  the last tf<N>/collisioncontext.root could be missing.
- The scan now closes every timeframe a collision skips over and pads the
  result to the number of timeframes the caller asks for, so entry i always
  describes orbits [start + i*orbitsPerTF, start + (i+1)*orbitsPerTF).
- applyMaxCollisionFilter keeps an empty timeframe empty when it re-indexes,
  and extractSingleTimeframe returns a valid empty context for it.
- o2-steer-colcontexttool passes the number of timeframes it asked for,
  reports timeframes that came out empty together with the mean number of
  collisions per timeframe implied by the interaction rate, and refuses to
  continue when --noEmptyTF was requested.

https://its.cern.ch/jira/browse/O2-7132

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

Author: sawenzel

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

Steer/src/CollisionContextTool.cxx

@@ -20,6 +20,7 @@
 #include "SimulationDataFormat/DigitizationContext.h"
 #include "SimConfig/InteractionDiamondParam.h"
 #include "DataFormatsFT0/EventsPerBc.h"
+#include "CommonConstants/LHCConstants.h"
 #include <cmath>
 #include <TRandom.h>
 #include <numeric>
@@ -56,7 +57,8 @@ struct Options {
   uint32_t firstBC = 0;    // first bunch crossing (relative to firstOrbit) of the first interaction;
   int orbitsPerTF = 256; // number of orbits per timeframe --> used to calculate start orbit for collisions
   bool useexistingkinematics = false;
-  bool noEmptyTF = false; // prevent empty timeframes; the first interaction will be shifted backwards to fall within the range given by Options.orbits
+  bool noEmptyTF = false;     // prevent empty timeframes; the first interaction will be shifted backwards to fall within the range given by Options.orbits
+  bool failOnEmptyTF = false; // stop rather than continue when a timeframe holds no collision
   int maxCollsPerTF = -1; // the maximal number of hadronic collisions per TF (can be used to constrain number of collisions per timeframe to some maximal value)
   std::string configKeyValues = ""; // string to init config key values
   long timestamp = -1;              // timestamp for CCDB queries
@@ -238,7 +240,8 @@ bool parseOptions(int argc, char* argv[], Options& optvalues)
     "timeframeID", bpo::value<int>(&optvalues.tfid)->default_value(0), "Timeframe id of the first timeframe int this context. Allows to generate contexts for different start orbits")(
     "first-orbit", bpo::value<double>(&optvalues.firstFractionalOrbit)->default_value(0), "First (fractional) orbit in the run (HBFUtils.firstOrbit + BC from decimal)")(
     "maxCollsPerTF", bpo::value<int>(&optvalues.maxCollsPerTF)->default_value(-1), "Maximal number of MC collisions to put into one timeframe. By default no constraint.")(
-    "noEmptyTF", bpo::bool_switch(&optvalues.noEmptyTF), "Enforce to have at least one collision")(
+    "noEmptyTF", bpo::bool_switch(&optvalues.noEmptyTF), "Shift the first collision backwards so that it falls within the sampled orbit range")(
+    "failOnEmptyTF", bpo::bool_switch(&optvalues.failOnEmptyTF), "Stop instead of continuing when one of the timeframes asked for ends up without a collision")(
     "configKeyValues", bpo::value<std::string>(&optvalues.configKeyValues)->default_value(""), "Semicolon separated key=value strings (e.g.: 'TPC.gasDensity=1;...')")(
     "with-vertices", bpo::value<std::string>(&optvalues.vertexModeString)->default_value("kNoVertex"), "Assign vertices to collisions. Argument is the vertex mode. Defaults to no vertexing applied")(
     "timestamp", bpo::value<long>(&optvalues.timestamp)->default_value(-1L), "Timestamp for CCDB queries / anchoring")(
@@ -664,7 +667,10 @@ int main(int argc, char* argv[])
   }
   LOG(info) << "-------- DENSE CONTEXT ------->>";
 
-  auto timeframeindices = digicontext.calcTimeframeIndices(orbitstart, options.orbitsPerTF, options.orbitsEarly);
+  // the number of timeframes we were asked for; passing it makes sure that a timeframe without
+  // collisions keeps its own slot instead of shifting every later timeframe down by one
+  long const num_timeframes_asked = usetimeframelength ? (orbits_total / options.orbitsPerTF) : -1;
+  auto timeframeindices = digicontext.calcTimeframeIndices(orbitstart, options.orbitsPerTF, options.orbitsEarly, num_timeframes_asked);
   LOG(info) << "Fixed " << timeframeindices.size() << " timeframes ";
   for (auto p : timeframeindices) {
     LOG(info) << std::get<0>(p) << " " << std::get<1>(p) << " " << std::get<2>(p);
@@ -688,6 +694,45 @@ int main(int argc, char* argv[])
 
   auto numTimeFrames = timeframeindices.size(); // digicontext.finalizeTimeframeStructure(orbitstart, options.orbitsPerTF, options.orbitsEarly);
 
+  // report - and, if asked, refuse - timeframes without a single collision. A timeframe with no
+  // collision cannot be simulated, and the rest of the MC workflow expects one collision context
+  // file per timeframe, so this has to be visible here and not five hours later in the simulation.
+  {
+    std::vector<int> empty_timeframes;
+    auto const first_real_tf = options.orbitsEarly > 0. ? 1 : 0;
+    for (int tf_id = first_real_tf; tf_id < (int)numTimeFrames; ++tf_id) {
+      if (std::get<0>(timeframeindices[tf_id]) > std::get<1>(timeframeindices[tf_id])) {
+        empty_timeframes.push_back(tf_id - first_real_tf + 1);
+      }
+    }
+    if (!empty_timeframes.empty()) {
+      std::stringstream tflist;
+      for (auto tf : empty_timeframes) {
+        tflist << " tf" << tf;
+      }
+      // the mean number of collisions in one timeframe, from the rate we were given
+      auto const tf_length_s = options.orbitsPerTF * o2::constants::lhc::LHCOrbitMUS * 1e-6;
+      double rate = 0.;
+      for (auto& p : ispecs) {
+        rate = std::max(rate, (double)p.interactionRate);
+      }
+      auto const mu_per_tf = rate * tf_length_s;
+      LOG(warn) << empty_timeframes.size() << " of " << (numTimeFrames - first_real_tf)
+                << " timeframes contain no collision:" << tflist.str();
+      LOG(warn) << "with interaction rate " << rate << " Hz and " << options.orbitsPerTF
+                << " orbits per timeframe there are only " << mu_per_tf
+                << " collisions per timeframe on average, so a fraction " << std::exp(-mu_per_tf)
+                << " of the timeframes comes out empty";
+      if (mu_per_tf > 0.) {
+        LOG(warn) << "use at least " << (int)std::ceil(8. / (rate * o2::constants::lhc::LHCOrbitMUS * 1e-6))
+                  << " orbits per timeframe to keep that fraction below 1 per mille";
+      }
+      if (options.failOnEmptyTF) {
+        LOG(fatal) << "--failOnEmptyTF was requested and timeframes without collisions were produced; refusing to continue";
+      }
+    }
+  }
+
   if (options.vertexMode != o2::conf::VertexMode::kNoVertex) {
     switch (options.vertexMode) {
       case o2::conf::VertexMode::kCCDB: {