Skip to content

Add a unit test for isFromRadDecay and fix its ancestry walk #15675

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
sawenzel wants to merge 1 commit into
base: dev
Choose a base branch
from
+68 −32
Open

Add a unit test for isFromRadDecay and fix its ancestry walk #15675

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
56 changes: 24 additions & 32 deletions Detectors/Base/include/DetectorsBase/Stack.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,8 +24,7 @@
#include "SimulationDataFormat/ParticleStatus.h"
#include "Rtypes.h"
#include "TParticle.h"
#include "TVirtualMC.h"
#include "TMCProcess.h"

#include <map>
#include <memory>
#include <stack>
Expand Down Expand Up @@ -211,7 +210,10 @@ class Stack : public FairGenericStack
/// query if a track is a direct **or** indirect daughter of a parentID
/// if trackid is same as parentid it returns true
bool isTrackDaughterOf(int /*trackid*/, int /*parentid*/) const;
bool isFromRadDecay(const int id);
/// query if a track originates, directly or indirectly, from a radioactive decay
/// only meaningful during transport, before selectTracks() remaps mother indices
bool isFromRadDecay(int trackid) const;

bool isCurrentTrackDaughterOf(int parentid) const;

// returns the index of the currently transported primary
Expand Down Expand Up @@ -349,39 +351,29 @@ inline int Stack::getMotherTrackId(int trackid) const
return mParticles[entryinParticles].getMotherTrackId();
}

inline bool Stack::isFromRadDecay(const int id)
inline bool Stack::isFromRadDecay(int trackid) const
{
// Check whether particle originates directly or indirectly from radioactive decay
// Check whether particle originates directly or indirectly from radioactive decay.
// Walks up the mother chain until a primary is reached. Only meaningful during
// transport, since selectTracks() later rewrites the mother indices in mParticles.
//
if (id < 0 || id >= static_cast<int>(mTrackIDtoParticlesEntry.size())) {
return false;
}
const auto entry = mTrackIDtoParticlesEntry[id];
if (entry < 0 || entry >= static_cast<int>(mParticles.size()))
return false;
auto part = (mParticles[entry]);

// primary particle ?
if (part.getProcess() == 0)
return false;
// particle directly from radioactive decay ?
if (part.getProcess() == kPRadDecay) {
return true;
}

// search in particle history
auto imo = mTrackIDtoParticlesEntry[part.getMotherTrackId()];
auto isRad = false;
while (imo > 0) {
auto mother = (mParticles[imo]);
if (mother.getProcess() == kPRadDecay) {
isRad = true;
break;
// Note that primaries are not kept in mParticles and that mTrackIDtoParticlesEntry
// is meaningless for them, so the chain has to be terminated on the trackID itself.
for (int id = trackid; id >= mNumberOfPrimaryParticles;) {
if (id >= static_cast<int>(mTrackIDtoParticlesEntry.size())) {
return false;
}
const auto entry = mTrackIDtoParticlesEntry[id];
if (entry < 0 || entry >= static_cast<int>(mParticles.size())) {
return false;
}
const auto& part = mParticles[entry];
if (part.getProcess() == kPRadDecay) {
return true;
}
part = mother;
imo = mTrackIDtoParticlesEntry[mother.getMotherTrackId()];
id = part.getMotherTrackId();
}
return isRad;
return false;
}

inline bool Stack::isCurrentTrackDaughterOf(int parentid) const
Expand Down
44 changes: 44 additions & 0 deletions Detectors/Base/test/testStack.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,3 +44,47 @@ BOOST_AUTO_TEST_CASE(Stack_test)
BOOST_CHECK(inst->getPrimaries().size() == 2);
}
}

// convenience wrapper to push a track and return the assigned trackID
static int pushTrack(o2::data::Stack& st, int parentId, TMCProcess proc)
{
int trackId;
st.PushTrack(1, parentId, 0, 0., 0., 0., 10., 5., 5., 5., 0.1, 0., 0., 0., proc, trackId, 1., 1);
return trackId;
}

// unit test for the radioactive-decay ancestry query
BOOST_AUTO_TEST_CASE(Stack_isFromRadDecay_test)
{
o2::data::Stack st;

// two primaries; note that primaries do not enter mParticles, only secondaries do
const auto prim0 = pushTrack(st, -1, kPPrimary);
const auto prim1 = pushTrack(st, -1, kPPrimary);

// a radioactive decay product of the second primary, and its descendants.
// this is deliberately the *first* secondary of the primary, so that it lands
// in the first entry of the particle buffer
const auto radDecay = pushTrack(st, prim1, kPRadDecay);
const auto radChild = pushTrack(st, radDecay, kPHadronic);
const auto radGrandChild = pushTrack(st, radChild, kPHadronic);

// a plain secondary of the second primary: no radioactive decay anywhere in its history
const auto ordinary = pushTrack(st, prim1, kPHadronic);

// primaries can never come from a radioactive decay
BOOST_CHECK(!st.isFromRadDecay(prim0));
BOOST_CHECK(!st.isFromRadDecay(prim1));

// a secondary whose ancestry ends in a primary must terminate the search with false
BOOST_CHECK(!st.isFromRadDecay(ordinary));

// directly and indirectly from a radioactive decay
BOOST_CHECK(st.isFromRadDecay(radDecay));
BOOST_CHECK(st.isFromRadDecay(radChild));
BOOST_CHECK(st.isFromRadDecay(radGrandChild));

// out-of-range track IDs are rejected rather than looked up
BOOST_CHECK(!st.isFromRadDecay(-1));
BOOST_CHECK(!st.isFromRadDecay(1000000000));
}