Merge 42a6068 into sapling-pr-archive-ktf

Author: ktf

CCDB/include/CCDB/CcdbObjectInfo.h

@@ -93,6 +93,16 @@ class CcdbObjectInfo
   [[nodiscard]] long getEndValidityTimestamp() const { return mEnd; }
   void setEndValidityTimestamp(long end) { mEnd = end; }
 
+  bool operator<(const CcdbObjectInfo& other) const
+  {
+    return mStart < other.mStart;
+  }
+
+  bool operator>(const CcdbObjectInfo& other) const
+  {
+    return mStart > other.mStart;
+  }
+
  private:
   std::string mObjType{};                 // object type (e.g. class)
   std::string mFileName{};                // file name in the CCDB
@@ -107,4 +117,17 @@ class CcdbObjectInfo
 
 } // namespace o2::ccdb
 
+namespace std
+{
+// defining std::hash for InteractionRecord to be used with std containers
+template <>
+struct hash<o2::ccdb::CcdbObjectInfo> {
+ public:
+  size_t operator()(const o2::ccdb::CcdbObjectInfo& info) const
+  {
+    return info.getStartValidityTimestamp();
+  }
+};
+} // namespace std
+
 #endif // O2_CCDB_CCDBOBJECTINFO_H_

DataFormats/simulation/include/SimulationDataFormat/DigitizationContext.h

@@ -115,8 +115,8 @@ class DigitizationContext
 
   /// retrieves collision context for a single timeframe-id (which may be needed by simulation)
   /// (Only copies collision context without QED information. This can be added to the result with the fillQED method
-  ///  in a second step. As a pre-condition, one should have called finalizeTimeframeStructure)
-  DigitizationContext extractSingleTimeframe(int timeframeid, std::vector<int> const& sources_to_offset);
+  ///  in a second step. Takes as input a timeframe indices collection)
+  DigitizationContext extractSingleTimeframe(int timeframeid, std::vector<std::tuple<int, int, int>> const& timeframeindices, std::vector<int> const& sources_to_offset);
 
   /// function reading the hits from a chain (previously initialized with initSimChains
   /// The hits pointer will be initialized (what to we do about ownership??)
@@ -130,12 +130,12 @@ class DigitizationContext
   /// returns the GRP object associated to this context
   o2::parameters::GRPObject const& getGRP() const;
 
-  // apply collision number cuts and potential relabeling of eventID
-  void applyMaxCollisionFilter(long startOrbit, long orbitsPerTF, int maxColl);
+  // apply collision number cuts and potential relabeling of eventID, (keeps collisions which fall into the orbitsEarly range for the next timeframe)
+  // needs a timeframe index structure (determined by calcTimeframeIndices), which is adjusted during the process to reflect the filtering
+  void applyMaxCollisionFilter(std::vector<std::tuple<int, int, int>>& timeframeindices, long startOrbit, long orbitsPerTF, int maxColl, double orbitsEarly = 0.);
 
-  /// finalize timeframe structure (fixes the indices in mTimeFrameStartIndex)
-  // returns the number of timeframes
-  int finalizeTimeframeStructure(long startOrbit, long orbitsPerTF);
+  /// 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;
 
   // 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.
@@ -176,17 +176,13 @@ class DigitizationContext
   // for each collision we record the constituents (which shall not exceed mMaxPartNumber)
   std::vector<std::vector<o2::steer::EventPart>> mEventParts;
 
-  // for each collision we may record/fix the interaction vertex (to be used in event generation)
+  // for each collisionstd::vector<std::tuple<int,int,int>> &timeframeindice we may record/fix the interaction vertex (to be used in event generation)
   std::vector<math_utils::Point3D<float>> mInteractionVertices;
 
   // the collision records **with** QED interleaved;
   std::vector<o2::InteractionTimeRecord> mEventRecordsWithQED;
   std::vector<std::vector<o2::steer::EventPart>> mEventPartsWithQED;
 
-  // timeframe structure
-  std::vector<std::pair<int, int>> mTimeFrameStartIndex;    // for each timeframe, the pair of start-index and end-index into mEventParts, mEventRecords
-  std::vector<std::pair<int, int>> mTimeFrameStartIndexQED; // for each timeframe, the pair of start-index and end-index into mEventParts, mEventRecords (QED version)
-
   o2::BunchFilling mBCFilling; // pattern of active BCs
 
   std::vector<std::string> mSimPrefixes;             // identifiers to the hit sim products; the key corresponds to the source ID of event record

DataFormats/simulation/src/DigitizationContext.cxx

@@ -49,6 +49,27 @@ void DigitizationContext::printCollisionSummary(bool withQED, int truncateOutput
     }
   } else {
     std::cout << "Number of Collisions " << mEventRecords.size() << "\n";
+    if (mEventPartsWithQED.size() > 0) {
+      auto num_qed_events = mEventPartsWithQED.size() - mEventRecords.size();
+      if (num_qed_events > 0) {
+        std::cout << "Number of QED events (but not shown) " << num_qed_events << "\n";
+        // find first and last QED collision so that we can give the range in orbits where these
+        // things are included
+        auto firstQEDcoll_iter = std::find_if(mEventPartsWithQED.begin(), mEventPartsWithQED.end(),
+                                              [](const std::vector<EventPart>& vec) {
+                                                return std::find_if(vec.begin(), vec.end(), [](EventPart const& p) { return p.sourceID == 99; }) != vec.end();
+                                              });
+
+        auto lastColl_iter = std::find_if(mEventPartsWithQED.rbegin(), mEventPartsWithQED.rend(),
+                                          [](const std::vector<EventPart>& vec) {
+                                            return std::find_if(vec.begin(), vec.end(), [](EventPart const& p) { return p.sourceID == 99; }) != vec.end();
+                                          });
+
+        auto firstindex = std::distance(mEventPartsWithQED.begin(), firstQEDcoll_iter);
+        auto lastindex = std::distance(mEventPartsWithQED.begin(), lastColl_iter.base()) - 1;
+        std::cout << "QED from: " << mEventRecordsWithQED[firstindex] << " ---> " << mEventRecordsWithQED[lastindex] << "\n";
+      }
+    }
     for (int i = 0; i < mEventRecords.size(); ++i) {
       if (truncateOutputTo >= 0 && i > truncateOutputTo) {
         std::cout << "--- Output truncated to " << truncateOutputTo << " ---\n";
@@ -380,32 +401,126 @@ std::vector<std::pair<int, int>> getTimeFrameBoundaries(std::vector<o2::Interact
   result.emplace_back(std::pair<int, int>(left, right - 1));
   return result;
 }
+
+// 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)
+{
+  // 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);
+
+  std::vector<std::tuple<int, int, int>> indices_with_early{};
+  for (int ti = 0; ti < true_indices.size(); ++ti) {
+    // for each timeframe we copy the true indices
+    auto& tf_range = true_indices[ti];
+
+    // init new index without fixing the early index yet
+    indices_with_early.push_back(std::make_tuple(tf_range.first, tf_range.second, -1));
+
+    // from the second timeframe on we can determine the index in the previous timeframe
+    // which matches our criterion
+    if (orbitsEarly > 0. && ti > 0) {
+      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
+      // (could probably use binary search, in case optimization becomes necessary)
+      int earlyOrbitIndex = prev_tf_range.second;
+
+      // this is the orbit of the ti-th timeframe start
+      auto orbit_timeframe_start = startOrbit + ti * orbitsPerTF;
+
+      auto orbit_timeframe_early_fractional = orbit_timeframe_start - orbitsEarly;
+      auto orbit_timeframe_early_integral = (uint32_t)(orbit_timeframe_early_fractional);
+
+      auto bc_early = (uint32_t)((orbit_timeframe_early_fractional - orbit_timeframe_early_integral) * o2::constants::lhc::LHCMaxBunches);
+
+      // this is the interaction record of the ti-th timeframe start
+      o2::InteractionRecord timeframe_start_record(0, orbit_timeframe_early_integral);
+      // this is the interaction record in some previous timeframe after which interactions could still
+      // influence the ti-th timeframe according to orbitsEarly
+      o2::InteractionRecord timeframe_early_record(bc_early, orbit_timeframe_early_integral);
+
+      auto differenceInBCNS_max = timeframe_start_record.differenceInBCNS(timeframe_early_record);
+
+      for (int j = prev_tf_range.second; j >= prev_tf_range.first; --j) {
+        // determine difference in timing in NS; compare that with the limit given by orbitsEarly
+        auto timediff_NS = timeframe_start_record.differenceInBCNS(irecords[j]);
+        if (timediff_NS < differenceInBCNS_max) {
+          earlyOrbitIndex = j;
+        } else {
+          break;
+        }
+      }
+      std::get<2>(indices_with_early.back()) = earlyOrbitIndex;
+    }
+  }
+  return indices_with_early;
+}
+
 } // namespace
 
-void DigitizationContext::applyMaxCollisionFilter(long startOrbit, long orbitsPerTF, int maxColl)
+void DigitizationContext::applyMaxCollisionFilter(std::vector<std::tuple<int, int, int>>& timeframeindices, long startOrbit, long orbitsPerTF, int maxColl, double orbitsEarly)
 {
   // the idea is to go through each timeframe and throw away collisions beyond a certain count
   // then the indices should be condensed
 
   std::vector<std::vector<o2::steer::EventPart>> newparts;
   std::vector<o2::InteractionTimeRecord> newrecords;
 
-  // get a timeframe boundary indexing
-  auto timeframeindices = getTimeFrameBoundaries(mEventRecords, startOrbit, orbitsPerTF);
-
   std::unordered_map<int, int> currMaxId;                           // the max id encountered for a source
   std::unordered_map<int, std::unordered_map<int, int>> reIndexMap; // for each source, a map of old to new index for the event parts
 
   if (maxColl == -1) {
     maxColl = mEventRecords.size();
   }
 
+  // the actual first actual timeframe
+  int first_timeframe = orbitsEarly > 0. ? 1 : 0;
+
+  // mapping of old to new indices
+  std::unordered_map<size_t, size_t> indices_old_to_new;
+
   // now we can go through the structure timeframe by timeframe
-  for (auto timeframe : timeframeindices) {
-    auto firstindex = timeframe.first;
-    auto lastindex = timeframe.second;
+  for (int tf_id = first_timeframe; tf_id < timeframeindices.size(); ++tf_id) {
+    auto& tf_indices = timeframeindices[tf_id];
+
+    auto firstindex = std::get<0>(tf_indices); // .first;
+    auto lastindex = std::get<1>(tf_indices);  // .second;
+    auto previndex = std::get<2>(tf_indices);
+
+    LOG(info) << "timeframe indices " << previndex << " : " << firstindex << " : " << lastindex;
+
+    int collCount = 0; // counting collisions within timeframe
     // copy to new structure
-    for (int index = firstindex; index <= std::min(lastindex, firstindex + maxColl - 1); ++index) {
+    for (int index = previndex >= 0 ? previndex : firstindex; index <= lastindex; ++index) {
+      if (collCount >= maxColl) {
+        continue;
+      }
+
+      // look if this index was already done?
+      // avoid duplicate entries in transformed records
+      if (indices_old_to_new.find(index) != indices_old_to_new.end()) {
+        continue;
+      }
+
+      // we put these events under a certain condition
+      bool keep = index < firstindex || collCount < maxColl;
+
+      if (!keep) {
+        continue;
+      }
+
+      if (index >= firstindex) {
+        collCount++;
+      }
+
+      // we must also make sure that we don't duplicate the records
+      // moreover some records are merely put as precoll of tf2 ---> so they shouldn't be part of tf1 in the final
+      // extraction, ouch !
+      // maybe we should combine the filter and individual tf extraction in one step !!
+      indices_old_to_new[index] = newrecords.size();
       newrecords.push_back(mEventRecords[index]);
       newparts.push_back(mEventParts[index]);
 
@@ -427,22 +542,35 @@ void DigitizationContext::applyMaxCollisionFilter(long startOrbit, long orbitsPe
           currMaxId[source] += 1;
         }
       }
+    } // ends one timeframe
+
+    // correct the timeframe indices
+    if (indices_old_to_new.find(firstindex) != indices_old_to_new.end()) {
+      std::get<0>(tf_indices) = indices_old_to_new[firstindex]; // start
+    }
+    if (indices_old_to_new.find(lastindex) != indices_old_to_new.end()) {
+      std::get<1>(tf_indices) = indices_old_to_new[lastindex]; // end;
+    } else {
+      std::get<1>(tf_indices) = newrecords.size(); // end;
+    }
+    if (indices_old_to_new.find(previndex) != indices_old_to_new.end()) {
+      std::get<2>(tf_indices) = indices_old_to_new[previndex]; // previous or "early" index
     }
   }
   // reassignment
   mEventRecords = newrecords;
   mEventParts = newparts;
 }
 
-int DigitizationContext::finalizeTimeframeStructure(long startOrbit, long orbitsPerTF)
+std::vector<std::tuple<int, int, int>> DigitizationContext::calcTimeframeIndices(long startOrbit, long orbitsPerTF, double orbitsEarly) const
 {
-  mTimeFrameStartIndex = getTimeFrameBoundaries(mEventRecords, startOrbit, orbitsPerTF);
-  LOG(info) << "Fixed " << mTimeFrameStartIndex.size() << " timeframes ";
-  for (auto p : mTimeFrameStartIndex) {
-    LOG(info) << p.first << " " << p.second;
+  auto timeframeindices = getTimeFrameBoundaries(mEventRecords, startOrbit, orbitsPerTF, orbitsEarly);
+  LOG(info) << "Fixed " << timeframeindices.size() << " timeframes ";
+  for (auto p : timeframeindices) {
+    LOG(info) << std::get<0>(p) << " " << std::get<1>(p) << " " << std::get<2>(p);
   }
 
-  return mTimeFrameStartIndex.size();
+  return timeframeindices;
 }
 
 std::unordered_map<int, int> DigitizationContext::getCollisionIndicesForSource(int source) const
@@ -529,21 +657,25 @@ void DigitizationContext::sampleInteractionVertices(o2::dataformats::MeanVertexO
   }
 }
 
-DigitizationContext DigitizationContext::extractSingleTimeframe(int timeframeid, std::vector<int> const& sources_to_offset)
+DigitizationContext DigitizationContext::extractSingleTimeframe(int timeframeid, std::vector<std::tuple<int, int, int>> const& timeframeindices, std::vector<int> const& sources_to_offset)
 {
   DigitizationContext r; // make a return object
-  if (mTimeFrameStartIndex.size() == 0) {
-    LOG(error) << "No timeframe structure determined; Returning empty object. Please call ::finalizeTimeframeStructure before calling this function";
+  if (timeframeindices.size() == 0) {
+    LOG(error) << "Timeframe index structure empty; Returning empty object.";
     return r;
   }
   r.mSimPrefixes = mSimPrefixes;
   r.mMuBC = mMuBC;
   try {
-    auto startend = mTimeFrameStartIndex.at(timeframeid);
+    auto tf_ranges = timeframeindices.at(timeframeid);
 
-    auto startindex = startend.first;
-    auto endindex = startend.second;
+    auto startindex = std::get<0>(tf_ranges);
+    auto endindex = std::get<1>(tf_ranges);
+    auto earlyindex = std::get<2>(tf_ranges);
 
+    if (earlyindex >= 0) {
+      startindex = earlyindex;
+    }
     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) {

Detectors/CTF/README.md

@@ -100,6 +100,11 @@ comma-separated list of detectors to skip
 ```
 max CTFs to process (<= 0 : infinite)
 
+```
+--max-tf-per-file arg (=-1)
+```
+max TFs to process from every CTF file (<= 0 : infinite)
+
 ```
 --loop arg (=0)
 ```

Detectors/CTF/workflow/include/CTFWorkflow/CTFReaderSpec.h

@@ -40,6 +40,7 @@ struct CTFReaderInp {
   int64_t delay_us = 0;
   int maxLoops = 0;
   int maxTFs = -1;
+  int maxTFsPerFile = -1;
   unsigned int subspec = 0;
   unsigned int decSSpecEMC = 0;
   int tfRateLimit = -999;

Detectors/CTF/workflow/src/CTFReaderSpec.cxx

@@ -148,6 +148,10 @@ void CTFReaderSpec::init(InitContext& ic)
   mUseLocalTFCounter = ic.options().get<bool>("local-tf-counter");
   mImposeRunStartMS = ic.options().get<int64_t>("impose-run-start-timstamp");
   mInput.checkTFLimitBeforeReading = ic.options().get<bool>("limit-tf-before-reading");
+  mInput.maxTFs = ic.options().get<int>("max-tf");
+  mInput.maxTFs = mInput.maxTFs > 0 ? mInput.maxTFs : 0x7fffffff;
+  mInput.maxTFsPerFile = ic.options().get<int>("max-tf-per-file");
+  mInput.maxTFsPerFile = mInput.maxTFsPerFile > 0 ? mInput.maxTFsPerFile : 0x7fffffff;
   mRunning = true;
   mFileFetcher = std::make_unique<o2::utils::FileFetcher>(mInput.inpdata, mInput.tffileRegex, mInput.remoteRegex, mInput.copyCmd);
   mFileFetcher->setMaxFilesInQueue(mInput.maxFileCache);
@@ -360,7 +364,7 @@ bool CTFReaderSpec::processTF(ProcessingContext& pc)
 void CTFReaderSpec::checkTreeEntries()
 {
   // check if the tree has entries left, if needed, close current tree/file
-  if (++mCurrTreeEntry >= mCTFTree->GetEntries()) { // this file is done, check if there are other files
+  if (++mCurrTreeEntry >= mCTFTree->GetEntries() || (mInput.maxTFsPerFile > 0 && mCurrTreeEntry >= mInput.maxTFsPerFile)) { // this file is done, check if there are other files
     mCTFTree.reset();
     mCTFFile->Close();
     mCTFFile.reset();
@@ -474,6 +478,9 @@ DataProcessorSpec getCTFReaderSpec(const CTFReaderInp& inp)
   options.emplace_back(ConfigParamSpec{"local-tf-counter", VariantType::Bool, false, {"reassign header.tfCounter from local TF counter"}});
   options.emplace_back(ConfigParamSpec{"fetch-failure-threshold", VariantType::Float, 0.f, {"Fail if too many failures( >0: fraction, <0: abs number, 0: no threshold)"}});
   options.emplace_back(ConfigParamSpec{"limit-tf-before-reading", VariantType::Bool, false, {"Check TF limiting before reading new TF, otherwhise before injecting it"}});
+  options.emplace_back(ConfigParamSpec{"max-tf", VariantType::Int, -1, {"max CTFs to process (<= 0 : infinite)"}});
+  options.emplace_back(ConfigParamSpec{"max-tf-per-file", VariantType::Int, -1, {"max TFs to process per ctf file (<= 0 : infinite)"}});
+
   if (!inp.metricChannel.empty()) {
     options.emplace_back(ConfigParamSpec{"channel-config", VariantType::String, inp.metricChannel, {"Out-of-band channel config for TF throttling"}});
   }

Detectors/CTF/workflow/src/ctf-reader-workflow.cxx

@@ -54,7 +54,6 @@ void customize(std::vector<o2::framework::ConfigParamSpec>& workflowOptions)
   options.push_back(ConfigParamSpec{"ctf-input", VariantType::String, "none", {"comma-separated list CTF input files"}});
   options.push_back(ConfigParamSpec{"onlyDet", VariantType::String, std::string{DetID::ALL}, {"comma-separated list of detectors to accept. Overrides skipDet"}});
   options.push_back(ConfigParamSpec{"skipDet", VariantType::String, std::string{DetID::NONE}, {"comma-separate list of detectors to skip"}});
-  options.push_back(ConfigParamSpec{"max-tf", VariantType::Int, -1, {"max CTFs to process (<= 0 : infinite)"}});
   options.push_back(ConfigParamSpec{"loop", VariantType::Int, 0, {"loop N times (infinite for N<0)"}});
   options.push_back(ConfigParamSpec{"delay", VariantType::Float, 0.f, {"delay in seconds between consecutive TFs sending"}});
   options.push_back(ConfigParamSpec{"copy-cmd", VariantType::String, "alien_cp ?src file://?dst", {"copy command for remote files or no-copy to avoid copying"}}); // Use "XrdSecPROTOCOL=sss,unix xrdcp -N root://eosaliceo2.cern.ch/?src ?dst" for direct EOS access
@@ -116,8 +115,6 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext)
   if (ctfInput.delay_us < 0) {
     ctfInput.delay_us = 0;
   }
-  int n = configcontext.options().get<int>("max-tf");
-  ctfInput.maxTFs = n > 0 ? n : 0x7fffffff;
 
   ctfInput.maxFileCache = std::max(1, configcontext.options().get<int>("max-cached-files"));