TPC quasi-triggered mode with --TPCtriggered option

For every processed sector the the TPCDigitRootWriter will recieve a vector of
RangeReference objects, telling how to group the digits and labels in the output
tree entries.
In default continuous RO mode this vector contains just 1 RangeReference, referring to
the all entries in the digits/labels container, therefore it will be dumped as a single entry.
In triggered mode, this vector will contain as many elements as MC interactions processed
each RangeReference will refer to the digits/labels of particular event, which will be stored
in separate entries.

Each RangeReference contains the index of of the 1st digit entry and N digits for the group it
refers to.

Author: shahor02

DataFormats/common/CMakeLists.txt

@@ -10,6 +10,7 @@ set(SRCS
 Set(HEADERS
   include/${MODULE_NAME}/TimeStamp.h
   include/${MODULE_NAME}/EvIndex.h
+  include/${MODULE_NAME}/RangeReference.h
   include/${MODULE_NAME}/InteractionRecord.h
   include/${MODULE_NAME}/BunchFilling.h  
 )

DataFormats/common/include/CommonDataFormat/RangeReference.h

@@ -0,0 +1,54 @@
+// Copyright CERN and copyright holders of ALICE O2. This software is
+// distributed under the terms of the GNU General Public License v3 (GPL
+// Version 3), copied verbatim in the file "COPYING".
+//
+// See http://alice-o2.web.cern.ch/license for full licensing information.
+//
+// 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.
+
+/// \file RangeReference.h
+/// \brief Class to refered to the 1st entry and N elements of some group in the continuous container
+/// \author ruben.shahoyan@cern.ch
+
+#ifndef ALICEO2_RANGEREFERENCE_H
+#define ALICEO2_RANGEREFERENCE_H
+
+#include <Rtypes.h>
+
+namespace o2
+{
+namespace dataformats
+{
+// Composed range reference
+
+template <typename FirstEntry = int, typename NElem = int>
+class RangeReference
+{
+ public:
+  RangeReference(FirstEntry ent, NElem n) { set(ent, n); }
+  RangeReference(const RangeReference<FirstEntry, NElem>& src) = default;
+  RangeReference() = default;
+  ~RangeReference() = default;
+  void set(FirstEntry ent, NElem n)
+  {
+    mFirstEntry = ent;
+    mEntries = n;
+  }
+  FirstEntry getFirstEntry() const { return mFirstEntry; }
+  NElem getEntries() const { return mEntries; }
+  void setFirstEntry(FirstEntry ent) { mFirstEntry = ent; }
+  void setEntries(NElem n) { mEntries = n; }
+  void changeEntriesBy(NElem inc) { mEntries += inc; }
+
+ private:
+  FirstEntry mFirstEntry; ///< 1st entry of the group
+  NElem mEntries = 0;     ///< number of entries
+
+  ClassDefNV(RangeReference, 1);
+};
+} // namespace dataformats
+} // namespace o2
+
+#endif

DataFormats/common/src/CommonDataFormatLinkDef.h

@@ -30,6 +30,9 @@
 #pragma link C++ class o2::dataformats::TimeStampWithError < int, int > +;
 
 #pragma link C++ class o2::dataformats::EvIndex < int, int > +;
+#pragma link C++ class o2::dataformats::RangeReference < int, int > +;
+#pragma link C++ class o2::dataformats::RangeReference < o2::dataformats::EvIndex < int, int >, int > +;
+
 #pragma link C++ class o2::InteractionRecord + ;
 #pragma link C++ class o2::BunchFilling + ;
 

Detectors/TPC/simulation/include/TPCSimulation/DigitizerTask.h

@@ -63,6 +63,9 @@ class DigitizerTask : public FairTask
   /// \param isContinuous - false for triggered readout, true for continuous readout
   void setContinuousReadout(bool isContinuous);
 
+  /// query if the r/o mode is continuous
+  bool isContinuousReadout() const { return mIsContinuousReadout; }
+
   /// Enable the use of space-charge distortions
   /// \param distortionType select the type of space-charge distortions (constant or realistic)
   /// \param hisInitialSCDensity optional space-charge density histogram to use at the beginning of the simulation

Detectors/TPC/simulation/include/TPCSimulation/HitDriftFilter.h

@@ -80,6 +80,69 @@ void getHits(TChain& chain, const Collection& eventrecords, std::vector<std::vec
   }
 }
 
+inline void getHits(std::vector<TChain*> const& chains, const o2::steer::RunContext& runcontext,
+                    std::vector<std::vector<o2::TPC::HitGroup>*>& hitvectors,
+                    std::vector<o2::TPC::TPCHitGroupID>& hitids, const std::string_view branchname,
+                    int collision)
+{
+  // Collect hits of single collision in the run context
+  hitids.clear();
+
+  // helper to fetch right branch
+  auto fetchBranch = [chains, branchname](int source) -> TBranch* {
+    if (!chains[source]) {
+      return nullptr;
+    }
+    auto& chain = *chains[source];
+    auto br = chain.GetBranch(branchname.data());
+    if (!br) {
+      return nullptr;
+    }
+    return br;
+  };
+
+  // number of collision given by runcontext
+  const auto& eventrecords = runcontext.getEventRecords();
+  auto ncollisions = eventrecords.size();
+  if (collision >= ncollisions) {
+    LOG(ERROR) << "runContext contains" << ncollisions << " collisions, " << collision << " is requested";
+    return;
+  }
+  const auto& parts = runcontext.getEventParts();
+  const auto maxpartspercollision = runcontext.getMaxNumberParts();
+
+  const auto& theseparts = parts[collision];
+  for (int p = 0; p < theseparts.size(); ++p) {
+    // compute where to store the hits
+    const auto eventID = theseparts[p].entryID;
+    const auto source = theseparts[p].sourceID;
+    const auto storeentry = eventID * maxpartspercollision + source;
+
+    // retrieve correct branch and fill the vector
+    // This needs to be done only once for any entry per chain
+    // retrieve source
+    if (hitvectors[storeentry] == nullptr) {
+      // TODO: instead of trying to fetch all the time .. do this outside and cache
+      auto br = fetchBranch(source);
+      br->SetAddress(&hitvectors[storeentry]);
+      br->GetEntry(eventID);
+    }
+
+    int groupid = -1;
+    auto groups = hitvectors[storeentry];
+    for (auto& singlegroup : *groups) {
+      if (singlegroup.getSize() == 0) {
+        // there are not hits in this group .. so continue
+        // TODO: figure out why such a group would exist??
+        continue;
+      }
+      groupid++;
+      // need to record index of the group
+      hitids.emplace_back(storeentry, collision, eventID, groupid, source);
+    } // end loop over entries
+  }
+}
+
 inline void getHits(std::vector<TChain*> const& chains, const o2::steer::RunContext& runcontext,
                     std::vector<std::vector<o2::TPC::HitGroup>*>& hitvectors,
                     std::vector<o2::TPC::TPCHitGroupID>& hitids, const std::string_view branchname, float tmin /*NS*/,

Detectors/TPC/simulation/src/Digitizer.cxx

@@ -59,7 +59,6 @@ DigitContainer* Digitizer::Process(const Sector& sector, const std::vector<o2::T
   if (!mIsContinuous) {
     eventTime = 0.f;
   }
-
   /// TODO: if eventtime-lastUpdate>=one space-charge time slice
   ///  1) Propagate current space-charge density
   ///  2) recalculate distortion lookup tables with updated space-charge density
@@ -79,7 +78,8 @@ DigitContainer* Digitizer::Process2(const Sector& sector, const std::vector<std:
     const auto hitvector = hits[id.storeindex];
     auto& group = (*hitvector)[id.groupID];
     auto& MCrecord = interactRecords[id.entry];
-    ProcessHitGroup(group, sector, MCrecord.timeNS * 0.001f, id.entry, id.sourceID);
+    float evTime = mIsContinuous ? MCrecord.timeNS * 0.001f : 0.f;
+    ProcessHitGroup(group, sector, evTime, id.entry, id.sourceID);
   }
 
   return mDigitContainer;

Steer/DigitizerWorkflow/src/ITSMFTDigitizerSpec.cxx

@@ -339,8 +339,7 @@ DataProcessorSpec getITSDigitizerSpec(int channel)
                               { "simFile", VariantType::String, "o2sim.root", { "Sim (background) input filename" } },
                               { "simFileS", VariantType::String, "", { "Sim (signal) input filename" } },
                               { "simFileQED", VariantType::String, "", { "Sim (QED) input filename" } },
-			      { (detStr + "triggered").c_str(), VariantType::Bool, false,
-			      { "Impose triggered RO mode (default: continuous)" } } } };
+                              { (detStr + "triggered").c_str(), VariantType::Bool, false, { "Impose triggered RO mode (default: continuous)" } } } };
 }
 
 DataProcessorSpec getMFTDigitizerSpec(int channel)
@@ -361,8 +360,7 @@ DataProcessorSpec getMFTDigitizerSpec(int channel)
                               { "simFile", VariantType::String, "o2sim.root", { "Sim (background) input filename" } },
                               { "simFileS", VariantType::String, "", { "Sim (signal) input filename" } },
                               { "simFileQED", VariantType::String, "", { "Sim (QED) input filename" } },
-                              { (detStr + "triggered").c_str(), VariantType::Bool, false,
-			      { "Impose triggered RO mode (default: continuous)" } } } };
+                              { (detStr + "triggered").c_str(), VariantType::Bool, false, { "Impose triggered RO mode (default: continuous)" } } } };
 }
 
 } // end namespace ITSMFT

Steer/DigitizerWorkflow/src/TPCDigitRootWriterSpec.cxx

@@ -15,6 +15,7 @@
 
 #include "TPCDigitRootWriterSpec.h"
 #include "DataFormatsTPC/TPCSectorHeader.h"
+#include "CommonDataFormat/RangeReference.h"
 #include "Framework/CallbackService.h"
 #include "Framework/ControlService.h"
 #include "TPCBase/Sector.h"
@@ -29,9 +30,11 @@
 #include <sstream>
 #include <string>
 #include <vector>
+#include <utility>
 
 using namespace o2::framework;
 using SubSpecificationType = o2::framework::DataAllocator::SubSpecificationType;
+using DigiGroupRef = o2::dataformats::RangeReference<int, int>;
 
 namespace o2
 {
@@ -59,31 +62,36 @@ DataProcessorSpec getTPCDigitRootWriterSpec(int numberofsourcedevices)
 {
   // assign input names to each channel
   auto digitchannelname = std::make_shared<std::vector<std::string>>();
+  auto triggerchannelname = std::make_shared<std::vector<std::string>>();
   auto labelchannelname = std::make_shared<std::vector<std::string>>();
   for (int i = 0; i < numberofsourcedevices; ++i) {
     {
       std::stringstream ss;
       ss << "digitinput" << i;
       digitchannelname->push_back(ss.str());
     }
+    {
+      std::stringstream ss;
+      ss << "triggerinput" << i;
+      triggerchannelname->push_back(ss.str());
+    }
     {
       std::stringstream ss;
       ss << "labelinput" << i;
       labelchannelname->push_back(ss.str());
     }
   }
 
-  auto initFunction = [numberofsourcedevices, digitchannelname, labelchannelname](InitContext& ic) {
+  auto initFunction = [numberofsourcedevices, digitchannelname, labelchannelname, triggerchannelname](InitContext& ic) {
     // get the option from the init context
     auto filename = ic.options().get<std::string>("tpc-digit-outfile");
     auto treename = ic.options().get<std::string>("treename");
 
     auto outputfile = std::make_shared<TFile>(filename.c_str(), "RECREATE");
     auto outputtree = std::make_shared<TTree>(treename.c_str(), treename.c_str());
 
-    // container for incoming digits + label
-    auto digits = std::make_shared<std::vector<o2::TPC::Digit>>();
-    auto labels = std::make_shared<o2::dataformats::MCTruthContainer<o2::MCCompLabel>>();
+    // container for cashed grouping of digits
+    auto trigP2Sect = std::make_shared<std::array<std::vector<DigiGroupRef>, 36>>();
 
     // the callback to be set as hook at stop of processing for the framework
     auto finishWriting = [outputfile, outputtree]() {
@@ -116,8 +124,8 @@ DataProcessorSpec getTPCDigitRootWriterSpec(int numberofsourcedevices)
     // using by-copy capture of the worker instance shared pointer
     // the shared pointer makes sure to clean up the instance when the processing
     // function gets out of scope
-    auto processingFct = [outputfile, outputtree, digits, digitchannelname, labelchannelname,
-                          numberofsourcedevices](ProcessingContext& pc) {
+    auto processingFct = [outputfile, outputtree, trigP2Sect, digitchannelname, labelchannelname,
+                          triggerchannelname, numberofsourcedevices](ProcessingContext& pc) {
       static bool finished = false;
       if (finished) {
         // avoid being executed again when marked as finished;
@@ -130,9 +138,11 @@ DataProcessorSpec getTPCDigitRootWriterSpec(int numberofsourcedevices)
       // need to record which channel has completed in order to decide when we can shutdown
       static std::vector<bool> digitsdone;
       static std::vector<bool> labelsdone;
+      static std::vector<bool> triggersdone;
       if (invocation == 1) {
         digitsdone.resize(numberofsourcedevices, false);
         labelsdone.resize(numberofsourcedevices, false);
+        triggersdone.resize(numberofsourcedevices, false);
       }
 
       // find out if all source devices (channels) are done
@@ -157,6 +167,19 @@ DataProcessorSpec getTPCDigitRootWriterSpec(int numberofsourcedevices)
       for (int d = 0; d < numberofsourcedevices; ++d) {
         const auto dname = digitchannelname->operator[](d);
         const auto lname = labelchannelname->operator[](d);
+        const auto tname = triggerchannelname->operator[](d);
+        if (pc.inputs().isValid(tname.c_str())) {
+          sector = extractSector(tname.c_str());
+          LOG(INFO) << "HAVE TRIGGER DATA FOR SECTOR " << sector << " ON CHANNEL " << d;
+          if (sector == -1) {
+            triggersdone[d] = true;
+          } else {
+            auto triggers = pc.inputs().get<std::vector<DigiGroupRef>>(tname.c_str());
+            (*trigP2Sect.get())[sector] = std::move(triggers);
+            const auto& trigS = (*trigP2Sect.get())[sector];
+            LOG(INFO) << "GOT Triggers of sector " << sector << " | SIZE " << trigS.size();
+          }
+        }
 
         // probe which channel has data and of what kind
         // a) probe for digits:
@@ -168,14 +191,39 @@ DataProcessorSpec getTPCDigitRootWriterSpec(int numberofsourcedevices)
           } else {
             // have to do work ...
             // the digits
-            auto indata = pc.inputs().get<std::vector<o2::TPC::Digit>>(dname.c_str());
-            LOG(INFO) << "DIGIT SIZE " << indata.size();
-            *digits.get() = std::move(indata);
+            auto digiData = pc.inputs().get<std::vector<o2::TPC::Digit>>(dname.c_str());
+            LOG(INFO) << "DIGIT SIZE " << digiData.size();
+            const auto& trigS = (*trigP2Sect.get())[sector];
+            if (!trigS.size()) {
+              LOG(FATAL) << "Digits for sector " << sector << " are received w/o info on grouping in triggers";
+            } else { // check consistency of Ndigits with that of expected from the trigger
+              int nExp = trigS.back().getFirstEntry() + trigS.back().getEntries() - trigS.front().getFirstEntry();
+              if (nExp != digiData.size()) {
+                LOG(ERROR) << "Number of digits " << digiData.size() << " is inconsistent with expectation " << nExp
+                           << " from digits grouping for sector " << sector;
+              }
+            }
+
             {
-              // connect this to a particular branch
-              auto br = getOrMakeBranch(*outputtree.get(), "TPCDigit", sector, digits.get());
-              br->Fill();
-              br->ResetAddress();
+              if (trigS.size() == 1) { // just 1 entry (continous mode?), use digits directly
+                // connect this to a particular branch
+                auto digP = &digiData;
+                auto br = getOrMakeBranch(*outputtree.get(), "TPCDigit", sector, digP);
+                br->Fill();
+                br->ResetAddress();
+              } else {                                // triggered mode (>1 entrie will be written)
+                std::vector<o2::TPC::Digit> digGroup; // group of digits related to single trigger
+                auto digGroupPtr = &digGroup;
+                auto br = getOrMakeBranch(*outputtree.get(), "TPCDigit", sector, digGroupPtr);
+                for (auto grp : trigS) {
+                  digGroup.clear();
+                  for (int i = 0; i < grp.getEntries(); i++) {
+                    digGroup.emplace_back(digiData[grp.getFirstEntry() + i]); // fetch digits of given trigger
+                  }
+                  br->Fill();
+                }
+                br->ResetAddress();
+              }
             }
           }
         } // end digit case
@@ -191,15 +239,36 @@ DataProcessorSpec getTPCDigitRootWriterSpec(int numberofsourcedevices)
             auto labeldata = pc.inputs().get<o2::dataformats::MCTruthContainer<o2::MCCompLabel>*>(lname.c_str());
             auto labeldataRaw = labeldata.get();
             LOG(INFO) << "MCTRUTH ELEMENTS " << labeldataRaw->getIndexedSize()
-                      << " WITH " << labeldataRaw->getNElements() << " LABELS ";
-            if (labeldataRaw->getIndexedSize() != digits->size()) {
-              LOG(WARNING) << "Inconsistent number of indexed (label) slots "
-                           << labeldataRaw->getIndexedSize() << " versus digits " << digits->size();
+                      << " WITH " << labeldataRaw->getNElements() << " LABELS";
+            const auto& trigS = (*trigP2Sect.get())[sector];
+            if (!trigS.size()) {
+              LOG(FATAL) << "MCTruth for sector " << sector << " are received w/o info on grouping in triggers";
+            } else {
+              int nExp = trigS.back().getFirstEntry() + trigS.back().getEntries() - trigS.front().getFirstEntry();
+              if (nExp != labeldataRaw->getIndexedSize()) {
+                LOG(ERROR) << "Number of indexed (label) slots " << labeldataRaw->getIndexedSize()
+                           << " is inconsistent with expectation " << nExp
+                           << " from digits grouping for sector " << sector;
+              }
             }
             {
-              auto br = getOrMakeBranch(*outputtree.get(), "TPCDigitMCTruth", sector, &labeldataRaw);
-              br->Fill();
-              br->ResetAddress();
+              if (trigS.size() == 1) { // just 1 entry (continous mode?), use labels directly
+                auto br = getOrMakeBranch(*outputtree.get(), "TPCDigitMCTruth", sector, &labeldataRaw);
+                br->Fill();
+                br->ResetAddress();
+              } else {
+                o2::dataformats::MCTruthContainer<o2::MCCompLabel> lblGroup; // labels for group of digits related to single trigger
+                auto lblGroupPtr = &lblGroup;
+                auto br = getOrMakeBranch(*outputtree.get(), "TPCDigitMCTruth", sector, &lblGroupPtr);
+                for (auto grp : trigS) {
+                  lblGroup.clear();
+                  for (int i = 0; i < grp.getEntries(); i++) {
+                    auto lbls = labeldataRaw->getLabels(grp.getFirstEntry() + i);
+                    lblGroup.addElements(i, lbls);
+                  }
+                  br->Fill();
+                }
+              }
             }
           }
         } // end label case
@@ -231,6 +300,8 @@ DataProcessorSpec getTPCDigitRootWriterSpec(int numberofsourcedevices)
   for (int d = 0; d < numberofsourcedevices; ++d) {
     inputs.emplace_back(InputSpec{ (*digitchannelname.get())[d].c_str(), "TPC", "DIGITS",
                                    static_cast<SubSpecificationType>(d), Lifetime::Timeframe }); // digit input
+    inputs.emplace_back(InputSpec{ (*triggerchannelname.get())[d].c_str(), "TPC", "DIGTRIGGERS",
+                                   static_cast<SubSpecificationType>(d), Lifetime::Timeframe }); // groupping in triggers
     inputs.emplace_back(InputSpec{ (*labelchannelname.get())[d].c_str(), "TPC", "DIGITSMCTR",
                                    static_cast<SubSpecificationType>(d), Lifetime::Timeframe });
   }