diff --git a/DataFormats/simulation/include/SimulationDataFormat/MCTrack.h b/DataFormats/simulation/include/SimulationDataFormat/MCTrack.h index b4cf26b11f82e..bb505b3b2e493 100644 --- a/DataFormats/simulation/include/SimulationDataFormat/MCTrack.h +++ b/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 }; diff --git a/Detectors/Base/include/DetectorsBase/Stack.h b/Detectors/Base/include/DetectorsBase/Stack.h index 479981a65477a..ed470337de4f5 100644 --- a/Detectors/Base/include/DetectorsBase/Stack.h +++ b/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; diff --git a/Detectors/Base/src/Stack.cxx b/Detectors/Base/src/Stack.cxx index a00c21c0589b9..a215ae27facc7 100644 --- a/Detectors/Base/src/Stack.cxx +++ b/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++; diff --git a/Detectors/Base/test/testStack.cxx b/Detectors/Base/test/testStack.cxx index f6d32d3bf7157..50caf7f244876 100644 --- a/Detectors/Base/test/testStack.cxx +++ b/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()); + } + } +} diff --git a/Detectors/EMCAL/simulation/src/Detector.cxx b/Detectors/EMCAL/simulation/src/Detector.cxx index e5163dbda8288..6dc4aded8a68d 100644 --- a/Detectors/EMCAL/simulation/src/Detector.cxx +++ b/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()); } } diff --git a/Detectors/FIT/FT0/simulation/src/Detector.cxx b/Detectors/FIT/FT0/simulation/src/Detector.cxx index 3b7bff8d51633..86b538581c2f8 100644 --- a/Detectors/FIT/FT0/simulation/src/Detector.cxx +++ b/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()); } diff --git a/Detectors/FOCAL/simulation/src/Detector.cxx b/Detectors/FOCAL/simulation/src/Detector.cxx index 464e57c07e676..9c89baa8f271b 100644 --- a/Detectors/FOCAL/simulation/src/Detector.cxx +++ b/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; } diff --git a/Detectors/PHOS/simulation/src/Detector.cxx b/Detectors/PHOS/simulation/src/Detector.cxx index a8ff99b3805f1..5e28a27230ef0 100644 --- a/Detectors/PHOS/simulation/src/Detector.cxx +++ b/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( diff --git a/Detectors/Upgrades/ALICE3/ECal/simulation/src/Detector.cxx b/Detectors/Upgrades/ALICE3/ECal/simulation/src/Detector.cxx index f0de8aa4022a6..1e049eaa741e1 100644 --- a/Detectors/Upgrades/ALICE3/ECal/simulation/src/Detector.cxx +++ b/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; } diff --git a/Detectors/ZDC/simulation/src/Detector.cxx b/Detectors/ZDC/simulation/src/Detector.cxx index c4466092c17e0..7685890b61426 100644 --- a/Detectors/ZDC/simulation/src/Detector.cxx +++ b/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()); }