Add contributesToAHit flag to MCTrack (#15996)

* Add contributesToAHit flag to MCTrack

This adds a per-track flag for tracks that deposit into a hit created by another track, and a unit test.

- MCTrack uses the previously reserved bit in mProp for the flag, so the storage layout is unchanged.
- Stack::modifiesHit(iDet) sets the flag on the current track without touching the hit mask, kHasHits or the hit counter.
- Kinematics pruning and output size are therefore unchanged.
- EMCAL, PHOS, FOCAL, ZDC, FT0 and ALICE 3 ECal call modifiesHit where a track adds to an existing hit.

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

Author: sawenzel

DataFormats/simulation/include/SimulationDataFormat/MCTrack.h

@@ -250,6 +250,15 @@ class MCTrackT
 
   bool isTransported() const { return getToBeDone() && !getInhibited(); };
 
+  /// mark that this track deposited into a hit created by another track
+  void setContributesToAHit(bool f = true)
+  {
+    auto prop = ((PropEncoding)mProp);
+    prop.contributesToAHit = f;
+    mProp = prop.i;
+  }
+  bool getContributesToAHit() const { return ((PropEncoding)mProp).contributesToAHit; }
+
   /// get the string representation of the production process
   const char* getProdProcessAsString() const;
 
@@ -286,7 +295,7 @@ class MCTrackT
       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 : NHITBITS;   // encoding hits per detector
-      int reserved1 : 1;        // bit reserved for possible future purposes
+      int contributesToAHit : 1; // whether this track deposited into a hit created by another track
       int inhibited : 1; // whether tracking of this was inhibited
       int toBeDone : 1; // whether this (still) needs tracking --> we might more complete information to cover full ParticleStatus space
     };

Detectors/Base/include/DetectorsBase/Stack.h

@@ -172,6 +172,9 @@ class Stack : public FairGenericStack
   /// Increment number of hits for the current track in a given detector
   /// \param iDet  Detector unique identifier
   void addHit(int iDet);
+  /// Mark the current track as contributing to a hit created by another track
+  /// \param iDet  Detector unique identifier
+  void modifiesHit(int iDet);
 
   TClonesArray* GetListOfParticles() override;
 

Detectors/Base/src/Stack.cxx

@@ -193,8 +193,8 @@ void Stack::PushTrack(Int_t toBeDone, Int_t parentId, Int_t pdgCode, Double_t px
                       TMCProcess proc2)
 {
   //  printf("Pushing %s toBeDone %5d parentId %5d pdgCode %5d is %5d entries %5d \n",
-  //	 proc == kPPrimary ? "Primary:   " : "Secondary: ",
-  //	 toBeDone, parentId, pdgCode, is, mNumberOfEntriesInParticles);
+  //   proc == kPPrimary ? "Primary:   " : "Secondary: ",
+  //   toBeDone, parentId, pdgCode, is, mNumberOfEntriesInParticles);
 
   //
   // This method is called
@@ -653,6 +653,16 @@ void Stack::addHit(int iDet)
   mCurrentParticle.SetBit(ParticleStatus::kHasHits, 1);
   mHitCounter++;
 }
+void Stack::modifiesHit(int iDet)
+{
+  // flag only; hit mask, kHasHits and hit counter stay untouched
+  if (mIndexOfCurrentTrack < mNumberOfPrimaryParticles) {
+    mTracks->at(mIndexOfCurrentTrack).setContributesToAHit();
+  } else {
+    Int_t iTrack = mTrackIDtoParticlesEntry[mIndexOfCurrentTrack];
+    mParticles[iTrack].setContributesToAHit();
+  }
+}
 void Stack::addHit(int iDet, Int_t iTrack)
 {
   mHitCounter++;

Detectors/Base/test/testStack.cxx

@@ -180,3 +180,42 @@ BOOST_AUTO_TEST_CASE(Offsetting_keeps_an_invalid_index_invalid)
   BOOST_CHECK_EQUAL(o2::base::Detector::offsetTrackIndex(7, nprimaries, primaryOffset, secondaryOffset), 107);
   BOOST_CHECK_EQUAL(o2::base::Detector::offsetTrackIndex(-1, nprimaries, primaryOffset, secondaryOffset), -1);
 }
+
+// modifiesHit flags the current track without touching hit mask or pruning
+BOOST_AUTO_TEST_CASE(Stack_modifiesHit_test)
+{
+  const int iDet = o2::detectors::DetID::EMC;
+  for (bool prune : {false, true}) {
+    o2::data::Stack st;
+    st.pruneKinematics(prune);
+    const auto prim = pushTrack(st, -1, kPPrimary);
+    st.SetCurrentTrack(prim);
+    // use grandchildren, since direct daughters of primaries are always kept
+    const auto sec = pushTrack(st, prim, kPHadronic);
+    const auto flagged = pushTrack(st, sec, kPHadronic);
+    const auto other = pushTrack(st, sec, kPHadronic);
+
+    st.modifiesHit(iDet);
+    st.SetCurrentTrack(flagged);
+    st.modifiesHit(iDet);
+    BOOST_CHECK(!st.currentTrackLeftHit());
+    st.SetCurrentTrack(other);
+    st.FinishPrimary();
+
+    auto const& tracks = *st.getMCTracks();
+    if (prune) {
+      // a flagged track without hits is still pruned
+      BOOST_CHECK_EQUAL(tracks.size(), 2);
+    } else {
+      BOOST_REQUIRE_EQUAL(tracks.size(), 4);
+      BOOST_CHECK(!tracks[sec].getContributesToAHit());
+      BOOST_CHECK(tracks[flagged].getContributesToAHit());
+      BOOST_CHECK(!tracks[other].getContributesToAHit());
+    }
+    BOOST_CHECK(tracks[prim].getContributesToAHit());
+    for (auto const& t : tracks) {
+      BOOST_CHECK_EQUAL(t.getHitMask(), 0);
+      BOOST_CHECK(!t.hasHits());
+    }
+  }
+}

Detectors/EMCAL/simulation/src/Detector.cxx

@@ -313,6 +313,7 @@ Bool_t Detector::ProcessHits(FairVolume* v)
     } else {
       LOG(debug3) << "Adding energy to the current hit" << std::endl;
       currenthit->SetEnergyLoss(currenthit->GetEnergyLoss() + lightyield);
+      o2stack->modifiesHit(GetDetId());
     }
   }
 

Detectors/FIT/FT0/simulation/src/Detector.cxx

@@ -1180,9 +1180,12 @@ Bool_t Detector::ProcessHits(FairVolume* v)
 o2::ft0::HitType* Detector::AddHit(float x, float y, float z, float time, float energy, Int_t trackId, Int_t detId)
 {
   mHits->emplace_back(x, y, z, time, energy, trackId, detId);
+  auto stack = (o2::data::Stack*)fMC->GetStack();
   if (energy == 10) {
-    auto stack = (o2::data::Stack*)fMC->GetStack();
     stack->addHit(GetDetId());
+  } else {
+    // optical photon contributing to a hit recorded for its parent
+    stack->modifiesHit(GetDetId());
   }
   return &(mHits->back());
 }

Detectors/FOCAL/simulation/src/Detector.cxx

@@ -1335,6 +1335,7 @@ bool Detector::ProcessHitsEPad(FairVolume* v)
   } else {
     LOG(debug3) << "Adding energy to the current hit";
     currenthit->SetEnergyLoss(currenthit->GetEnergyLoss() + eloss);
+    o2stack->modifiesHit(GetDetId());
   }
   return true;
 }
@@ -1384,6 +1385,7 @@ bool Detector::ProcessHitsHCAL(FairVolume* v)
   } else {
     LOG(debug3) << "Adding energy to the current hit";
     currenthit->SetEnergyLoss(currenthit->GetEnergyLoss() + eloss);
+    o2stack->modifiesHit(GetDetId());
   }
   return true;
 }
@@ -1433,6 +1435,7 @@ bool Detector::ProcessHitsEPix(FairVolume* v)
   } else {
     LOG(debug3) << "Adding energy to the current hit";
     currenthit->SetEnergyLoss(currenthit->GetEnergyLoss() + eloss);
+    o2stack->modifiesHit(GetDetId());
   }
   return true;
 }

Detectors/PHOS/simulation/src/Detector.cxx

@@ -169,6 +169,9 @@ Bool_t Detector::ProcessHits(FairVolume* v)
   if (lostenergy < DBL_EPSILON && !isNewPartile) {
     return false; // do not create hits with zero energy deposition
   }
+  if (!isNewPartile) { // track deposits into a hit of its super-parent
+    stack->modifiesHit(GetDetId());
+  }
 
   Int_t moduleNumber;
   fMC->CurrentVolOffID(

Detectors/Upgrades/ALICE3/ECal/simulation/src/Detector.cxx

@@ -383,6 +383,7 @@ bool Detector::ProcessHits(FairVolume* vol)
     stack->addHit(GetDetId());
   } else {
     hit->SetEnergyLoss(hit->GetEnergyLoss() + eloss);
+    stack->modifiesHit(GetDetId());
   }
   return true;
 }

Detectors/ZDC/simulation/src/Detector.cxx

@@ -527,6 +527,7 @@ Bool_t Detector::ProcessHits(FairVolume* v)
       curHit.SetEnergyLoss(incenloss);
       curHit.setPMCLightYield(curHit.getPMCLightYield() + nPMC);
       curHit.setPMQLightYield(curHit.getPMQLightYield() + nPMQ);
+      stack->modifiesHit(GetDetId());
       //printf("  >>> Hit updated in vol %d %d  for track %d (mother: %d) \t light %1.0f %1.0f \n",
       //detector, sector, trackn, stack->GetCurrentTrack()->GetFirstMother(), curHit.getPMCLightYield(), curHit.getPMQLightYield());
     }