MCTrack: Generalize hit bit property

So far, every possible sensitive detector for ALICE Run3
had a dedicated bit assigned inside MCTracks. The number of such
bit was fixed to be 22 maximally.
In this way we could store information if a particular Monte Carlo track
caused energy deposits in detectors.

Since more and more detectors are being added to the code
(ALICE3, ITS3, R&D), we'll soon have more than 22 DetIDs in our code
which will conflict with the space available in MCTrack.

Since enlarging the MCTrack class in memory is to be avoided, this commit
proposes a solution based on an intermediate lookup table which dynamically
maps the detector ID to a particular bit. The lookup table is setup during
geometry construction and stored in MCEventHeaders for offline use.

This solution should scale much further than before. The limit is now 22 simultaneous
sensitive detectors in the geometry and not 22 total sensitive detectors in the code.

Author: sawenzel

DataFormats/simulation/include/SimulationDataFormat/MCEventHeader.h

@@ -86,6 +86,9 @@ class MCEventHeader : public FairMCEventHeader
 
   MCEventStats& getMCEventStats() { return mEventStats; }
 
+  void setDetId2HitBitLUT(std::vector<int> const& v) { mDetId2HitBitIndex = v; }
+  std::vector<int> const& getDetId2HitBitLUT() const { return mDetId2HitBitIndex; }
+
   /// create a standalone ROOT file/tree with only the MCHeader branch
   static void extractFileFromKinematics(std::string_view kinefilename, std::string_view targetfilename);
 
@@ -96,9 +99,10 @@ class MCEventHeader : public FairMCEventHeader
   // store a few global properties that this event
   // had in the current simulation (which can be used quick filtering/searching)
   MCEventStats mEventStats{};
+  std::vector<int> mDetId2HitBitIndex;
   o2::utils::RootSerializableKeyValueStore mEventInfo;
 
-  ClassDefOverride(MCEventHeader, 3);
+  ClassDefOverride(MCEventHeader, 4);
 
 }; /** class MCEventHeader **/
 

DataFormats/simulation/include/SimulationDataFormat/MCTrack.h

@@ -46,6 +46,8 @@ template <class _T>
 class MCTrackT
 {
  public:
+  static constexpr int NHITBITS = 22; // do not modify this
+
   ///  Default constructor
   MCTrackT();
 
@@ -159,24 +161,40 @@ class MCTrackT
   void SetSecondMotherTrackId(Int_t id) { mSecondMotherTrackId = id; }
   void SetFirstDaughterTrackId(Int_t id) { mFirstDaughterTrackId = id; }
   void SetLastDaughterTrackId(Int_t id) { mLastDaughterTrackId = id; }
+
   // set bit indicating that this track
-  // left a hit in detector with id iDet
-  void setHit(Int_t iDet)
+  // left a hit in detector that corresponds to bit iDetBit. This bit is found in
+  // a lookup table.
+  void setHit(Int_t iDetBit)
   {
-    assert(0 <= iDet && iDet < o2::detectors::DetID::nDetectors);
+    assert(0 <= iDetBit && iDetBit < o2::detectors::DetID::nDetectors);
     auto prop = ((PropEncoding)mProp);
-    prop.hitmask |= 1 << iDet;
+    prop.hitmask |= 1 << iDetBit;
     mProp = prop.i;
   }
 
-  // did detector iDet see this track?
-  bool leftTrace(Int_t iDet) const { return (((PropEncoding)mProp).hitmask & (1 << iDet)) > 0; }
+  bool leftTraceGivenBitField(int bit) const
+  {
+    return (((PropEncoding)mProp).hitmask & (1 << bit)) > 0;
+  }
+
+  /// Returns of detector with DetID iDet has seen this track.
+  /// Needs lookup table detIDToBit mapping detectorIDs to bitfields (which is persistified in MCEventHeaders).
+  bool leftTrace(Int_t iDet, std::vector<int> const& detIDtoBit) const
+  {
+    auto bit = detIDtoBit[iDet];
+    if (bit != -1) {
+      return leftTraceGivenBitField(bit);
+    }
+    return false;
+  }
+
   // determine how many detectors "saw" this track
   int getNumDet() const
   {
     int count = 0;
-    for (auto i = o2::detectors::DetID::First; i < o2::detectors::DetID::nDetectors; ++i) {
-      if (leftTrace(i)) {
+    for (auto i = 0; i < NHITBITS; ++i) {
+      if (leftTraceGivenBitField(i)) {
         count++;
       }
     }
@@ -261,7 +279,7 @@ class MCTrackT
     struct {
       int storage : 1;  // encoding whether to store this track to the output
       unsigned int process : 6; // encoding process that created this track (enough to store TMCProcess from ROOT)
-      int hitmask : 22;         // encoding hits per detector
+      int hitmask : NHITBITS;   // encoding hits per detector
       static_assert(o2::detectors::DetID::nDetectors <= 22); // ensure that all known detectors can be encoded here by a bit
       int reserved1 : 1;        // bit reserved for possible future purposes
       int inhibited : 1; // whether tracking of this was inhibited

DataFormats/simulation/test/MCTrack.cxx

@@ -24,10 +24,17 @@ using namespace o2;
 BOOST_AUTO_TEST_CASE(MCTrack_test)
 {
   MCTrack track;
+
+  // auxiliary lookup table needed to fetch and set hit properties
+  std::vector<int> hitLUT(o2::detectors::DetID::nDetectors, -1);
+  // in this test we have a single fictional detector 1, which we map to
+  // the first bit
+  hitLUT[1] = 0;
+
   // check properties on default constructed object
   BOOST_CHECK(track.getStore() == false);
   for (auto i = o2::detectors::DetID::First; i < o2::detectors::DetID::nDetectors; ++i) {
-    BOOST_CHECK(track.leftTrace(i) == false);
+    BOOST_CHECK(track.leftTrace(i, hitLUT) == false);
   }
   BOOST_CHECK(track.getNumDet() == 0);
   BOOST_CHECK(track.hasHits() == false);
@@ -41,10 +48,10 @@ BOOST_AUTO_TEST_CASE(MCTrack_test)
   BOOST_CHECK(track.getStore() == true);
 
   // set hit for first detector
-  BOOST_CHECK(track.leftTrace(1) == false);
-  track.setHit(1);
+  BOOST_CHECK(track.leftTrace(1, hitLUT) == false);
+  track.setHit(hitLUT[1]);
   BOOST_CHECK(track.hasHits() == true);
-  BOOST_CHECK(track.leftTrace(1) == true);
+  BOOST_CHECK(track.leftTrace(1, hitLUT) == true);
   BOOST_CHECK(track.getNumDet() == 1);
 
   // check process encoding

Detectors/Base/include/DetectorsBase/Detector.h

@@ -216,6 +216,10 @@ class Detector : public FairDetector
   // which is required for media creation
   static void initFieldTrackingParams(int& mode, float& maxfield);
 
+  /// set the DetID to HitBitIndex mapping. Succeeds if not already set.
+  static void setDetId2HitBitIndex(std::vector<int> const& v) { Detector::sDetId2HitBitIndex = v; }
+  static std::vector<int> const& getDetId2HitBitIndex() { return Detector::sDetId2HitBitIndex; }
+
  protected:
   Detector(const Detector& origin);
 
@@ -233,6 +237,7 @@ class Detector : public FairDetector
 
   static Float_t mDensityFactor; //! factor that is multiplied to all material densities (ONLY for
   // systematic studies)
+  static std::vector<int> sDetId2HitBitIndex; //! global lookup table keeping mapping of DetID to index in hit bit field (used in MCTrack)
 
   ClassDefOverride(Detector, 1); // Base class for ALICE Modules
 };
@@ -738,6 +743,7 @@ class DetImpl : public o2::base::Detector
   int mInitialized = false;
 
   char* mHitCollectorBufferPtr = nullptr; //! pointer to hit (collector) buffer location (strictly internal)
+
   ClassDefOverride(DetImpl, 0);
 };
 } // namespace base

Detectors/Base/src/Detector.cxx

@@ -30,6 +30,7 @@ using namespace o2::base;
 using namespace o2::detectors;
 
 Float_t Detector::mDensityFactor = 1.0;
+std::vector<int> o2::base::Detector::sDetId2HitBitIndex{}; // initialize empty vector
 
 Detector::Detector() : FairDetector(), mMapMaterial(), mMapMedium() {}
 Detector::Detector(const char* name, Bool_t Active)

Detectors/Base/src/Stack.cxx

@@ -638,14 +638,17 @@ void Stack::Print(Option_t* option) const
 
 void Stack::addHit(int iDet)
 {
+  // translate detector ID to bit id
+  const auto bitID = o2::base::Detector::getDetId2HitBitIndex()[iDet];
+
   if (mIndexOfCurrentTrack < mNumberOfPrimaryParticles) {
     auto& part = mTracks->at(mIndexOfCurrentTrack);
-    part.setHit(iDet);
+    part.setHit(bitID);
 
   } else {
     Int_t iTrack = mTrackIDtoParticlesEntry[mIndexOfCurrentTrack];
     auto& part = mParticles[iTrack];
-    part.setHit(iDet);
+    part.setHit(bitID);
   }
   mCurrentParticle.SetBit(ParticleStatus::kHasHits, 1);
   mHitCounter++;
@@ -654,7 +657,9 @@ void Stack::addHit(int iDet, Int_t iTrack)
 {
   mHitCounter++;
   auto& part = mParticles[iTrack];
-  part.setHit(iDet);
+  // fetch the bit encoding for hits
+  const auto bitID = o2::base::Detector::getDetId2HitBitIndex()[iDet];
+  part.setHit(bitID);
   mCurrentParticle.SetBit(ParticleStatus::kHasHits, 1);
 }
 

Steer/src/O2MCApplication.cxx

@@ -197,6 +197,7 @@ void O2MCApplicationBase::finishEventCommon()
 
   auto header = static_cast<o2::dataformats::MCEventHeader*>(fMCEventHeader);
   header->getMCEventStats().setNSteps(mStepCounter);
+  header->setDetId2HitBitLUT(o2::base::Detector::getDetId2HitBitIndex());
 
   static_cast<o2::data::Stack*>(GetStack())->updateEventStats();
 }

macro/build_geometry.C

@@ -179,103 +179,120 @@ void build_geometry(FairRunSim* run = nullptr)
     run->AddModule(new o2::passive::FrameStructure("FRAME", "Frame"));
   }
 
+  std::vector<int> detId2RunningId = std::vector<int>(o2::detectors::DetID::nDetectors, -1); // a mapping of detectorId to a dense runtime index
+  // used for instance to set bits in the hit structure of MCTracks; -1 means that there is no bit associated
+
+  auto addReadoutDetector = [&detId2RunningId, &run](o2::base::Detector* detector) {
+    static int runningid = 0; // this is static for constant lambda interfaces --> use fixed type and not auto in the lambda!
+    run->AddModule(detector);
+    if (detector->IsActive()) {
+      auto detID = detector->GetDetId();
+      detId2RunningId[detID] = runningid;
+      LOG(info) << " DETID " << detID << " vs " << detector->GetDetId() << " mapped to hit bit index " << runningid;
+      runningid++;
+    }
+  };
+
   if (isActivated("TOF")) {
     // TOF
-    run->AddModule(new o2::tof::Detector(isReadout("TOF")));
+    addReadoutDetector(new o2::tof::Detector(isReadout("TOF")));
   }
 
   if (isActivated("TRD")) {
     // TRD
-    run->AddModule(new o2::trd::Detector(isReadout("TRD")));
+    addReadoutDetector(new o2::trd::Detector(isReadout("TRD")));
   }
 
   if (isActivated("TPC")) {
     // tpc
-    run->AddModule(new o2::tpc::Detector(isReadout("TPC")));
+    addReadoutDetector(new o2::tpc::Detector(isReadout("TPC")));
   }
 #ifdef ENABLE_UPGRADES
   if (isActivated("IT3")) {
     // IT3
-    run->AddModule(new o2::its::Detector(isReadout("IT3"), "IT3"));
+    addReadoutDetector(new o2::its::Detector(isReadout("IT3"), "IT3"));
   }
 
   if (isActivated("TRK")) {
     // ALICE 3 TRK
-    run->AddModule(new o2::trk::Detector(isReadout("TRK")));
+    addReadoutDetector(new o2::trk::Detector(isReadout("TRK")));
   }
 
   if (isActivated("FT3")) {
     // ALICE 3 FT3
-    run->AddModule(new o2::ft3::Detector(isReadout("FT3")));
+    addReadoutDetector(new o2::ft3::Detector(isReadout("FT3")));
   }
 
   if (isActivated("FCT")) {
     // ALICE 3 FCT
-    run->AddModule(new o2::fct::Detector(isReadout("FCT")));
+    addReadoutDetector(new o2::fct::Detector(isReadout("FCT")));
   }
 #endif
 
   if (isActivated("ITS")) {
     // its
-    run->AddModule(new o2::its::Detector(isReadout("ITS")));
+    addReadoutDetector(new o2::its::Detector(isReadout("ITS")));
   }
 
   if (isActivated("MFT")) {
     // mft
-    run->AddModule(new o2::mft::Detector(isReadout("MFT")));
+    addReadoutDetector(new o2::mft::Detector(isReadout("MFT")));
   }
 
   if (isActivated("MCH")) {
     // mch
-    run->AddModule(new o2::mch::Detector(isReadout("MCH")));
+    addReadoutDetector(new o2::mch::Detector(isReadout("MCH")));
   }
 
   if (isActivated("MID")) {
     // mid
-    run->AddModule(new o2::mid::Detector(isReadout("MID")));
+    addReadoutDetector(new o2::mid::Detector(isReadout("MID")));
   }
 
   if (isActivated("EMC")) {
     // emcal
-    run->AddModule(new o2::emcal::Detector(isReadout("EMC")));
+    addReadoutDetector(new o2::emcal::Detector(isReadout("EMC")));
   }
 
   if (isActivated("PHS")) {
     // phos
-    run->AddModule(new o2::phos::Detector(isReadout("PHS")));
+    addReadoutDetector(new o2::phos::Detector(isReadout("PHS")));
   }
 
   if (isActivated("CPV")) {
     // cpv
-    run->AddModule(new o2::cpv::Detector(isReadout("CPV")));
+    addReadoutDetector(new o2::cpv::Detector(isReadout("CPV")));
   }
 
   if (isActivated("FT0")) {
     // FIT-T0
-    run->AddModule(new o2::ft0::Detector(isReadout("FT0")));
+    addReadoutDetector(new o2::ft0::Detector(isReadout("FT0")));
   }
 
   if (isActivated("FV0")) {
     // FIT-V0
-    run->AddModule(new o2::fv0::Detector(isReadout("FV0")));
+    addReadoutDetector(new o2::fv0::Detector(isReadout("FV0")));
   }
 
   if (isActivated("FDD")) {
     // FIT-FDD
-    run->AddModule(new o2::fdd::Detector(isReadout("FDD")));
+    addReadoutDetector(new o2::fdd::Detector(isReadout("FDD")));
   }
 
   if (isActivated("HMP")) {
     // HMP
-    run->AddModule(new o2::hmpid::Detector(isReadout("HMP")));
+    addReadoutDetector(new o2::hmpid::Detector(isReadout("HMP")));
   }
 
   if (isActivated("ZDC")) {
     // ZDC
-    run->AddModule(new o2::zdc::Detector(isReadout("ZDC")));
+    addReadoutDetector(new o2::zdc::Detector(isReadout("ZDC")));
   }
 
   if (geomonly) {
     run->Init();
   }
+
+  // register the DetId2HitIndex lookup with the detector class by copying the vector
+  o2::base::Detector::setDetId2HitBitIndex(detId2RunningId);
 }