Skip to content

[ALICE3] Fix indices bug in otf decayer #16587

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

Merged
jesgum merged 9 commits into from
Jun 12, 2026
+58 −39
Merged

[ALICE3] Fix indices bug in otf decayer #16587

Show file tree
Hide file tree
Changes from 7 commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
572bd82
Fix indices bug in otf decayer
jesgum Jun 9, 2026
fac0b57
update
jesgum Jun 9, 2026
da02d0c
remove unused
jesgum Jun 9, 2026
bc39fbc
fix offset
jesgum Jun 9, 2026
2d572a0
fix const& in loops
jesgum Jun 11, 2026
c490e90
anton's fixes
jesgum Jun 12, 2026
41354cb
fix for particles that already decayed in generator
jesgum Jun 12, 2026
9b46d38
remove unused
jesgum Jun 12, 2026
21c7915
update
jesgum Jun 12, 2026
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
17 changes: 14 additions & 3 deletions ALICE3/Core/OTFParticle.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,9 @@ class OTFParticle
if (particle.has_mothers()) {
mIndicesMother = {particle.mothersIds().front(), particle.mothersIds().back()};
}
if (particle.has_daughters()) {
mIndicesDaughter = {particle.daughtersIds().front(), particle.daughtersIds().back()};
}
if constexpr (requires { particle.decayerBits(); }) {
mBits = particle.decayerBits();
} else {
Expand Down Expand Up @@ -88,6 +91,14 @@ class OTFParticle
mPz = pz;
mE = e;
}
void setIndexOffset(const std::size_t offset)
{
static constexpr int NotFound = -1;
mIndicesMother[0] = (mIndicesMother[0] >= 0) ? mIndicesMother[0] + static_cast<int>(offset) : NotFound;
mIndicesMother[1] = (mIndicesMother[1] >= 0) ? mIndicesMother[1] + static_cast<int>(offset) : NotFound;
mIndicesDaughter[0] = (mIndicesDaughter[0] >= 0) ? mIndicesDaughter[0] + static_cast<int>(offset) : NotFound;
mIndicesDaughter[1] = (mIndicesDaughter[1] >= 0) ? mIndicesDaughter[1] + static_cast<int>(offset) : NotFound;
}

// Getters
int pdgCode() const { return mPdgCode; }
Expand Down Expand Up @@ -147,8 +158,8 @@ class OTFParticle
std::span<const int> getMotherSpan() const { return hasMothers() ? std::span<const int>(mIndicesMother.data(), 2) : std::span<const int>(); }

// Checks
bool hasDaughters() const { return (mIndicesDaughter[0] > 0); }
bool hasMothers() const { return (mIndicesMother[0] > 0); }
bool hasDaughters() const { return (mIndicesDaughter[0] >= 0); }
bool hasMothers() const { return (mIndicesMother[0] >= 0); }
bool hasNaN() const
{
return std::isnan(mPx) || std::isnan(mPy) || std::isnan(mPz) || std::isnan(mE) ||
Expand All @@ -171,7 +182,7 @@ class OTFParticle

private:
int mPdgCode{}, mGlobalIndex{-1};
int mCollisionId{};
int mCollisionId{-1};
float mVx{}, mVy{}, mVz{}, mVt{};
float mPx{}, mPy{}, mPz{}, mE{};
bool mIsAlive{}, mIsFromMcParticles{false};
Expand Down
63 changes: 36 additions & 27 deletions ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,8 +49,8 @@
static constexpr int NumDecays = 7;
static constexpr int NumParameters = 1;
static constexpr int DefaultParameters[NumDecays][NumParameters]{{1}, {1}, {1}, {1}, {1}, {1}, {1}};
static const std::vector<std::string> ParameterNames{"enable"};
static const std::vector<std::string> ParticleNames{"K0s",
static const std::vector<std::string> parameterNames{"enable"};
static const std::vector<std::string> particleNames{"K0s",
"Lambda",
"Anti-Lambda",
"Xi",
Expand Down Expand Up @@ -83,12 +83,12 @@
Configurable<int> seed{"seed", 0, "Set seed for particle decayer"};
Configurable<float> magneticField{"magneticField", 20., "Magnetic field (kG)"};
Configurable<LabeledArray<int>> enabledDecays{"enabledDecays",
{DefaultParameters[0], NumDecays, NumParameters, ParticleNames, ParameterNames},
{DefaultParameters[0], NumDecays, NumParameters, particleNames, parameterNames},
"Enable option for particle to be decayed: 0 - no, 1 - yes"};

std::size_t indexOffset = 0;
std::size_t particlesInDataframe = 0;
static constexpr float PicoToNano = 1.e-3f;
int mCollisionId{-1};

std::vector<int> mEnabledDecays;
void init(o2::framework::InitContext&)
{
Expand All @@ -98,7 +98,7 @@
decayer.setSeed(seed);
decayer.setBField(magneticField);
for (int i = 0; i < NumDecays; ++i) {
if (enabledDecays->get(ParticleNames[i].c_str(), "enable")) {
if (enabledDecays->get(particleNames[i].c_str(), "enable")) {
LOG(info) << " --- Decay enabled: " << pdgCodes[i];
mEnabledDecays.push_back(pdgCodes[i]);
}
Expand All @@ -121,6 +121,10 @@
std::vector<o2::upgrade::OTFParticle> allParticles;
void decayParticles(const int start, const int stop)
{
if (start >= stop) {
return;
}

int ndau = 0;
for (int i = start; i < stop; ++i) {
o2::upgrade::OTFParticle& particle = allParticles[i];
Expand All @@ -135,6 +139,10 @@

particle.setBitOff(o2::upgrade::DecayerBits::IsAlive);
std::vector<o2::upgrade::OTFParticle> decayStack = decayer.decayParticle(pdgDB, particle);
if (decayStack.empty()) {
continue;
}

const float decayRadius = decayer.getDecayRadius();
const float trackVelocity = o2::upgrade::computeParticleVelocity(particle.p(), pdgDB->GetParticle(particle.pdgCode())->Mass());
const int charge = pdgDB->GetParticle(particle.pdgCode())->Charge() / 3;
Expand All @@ -151,20 +159,16 @@
}

const float trackTimeNS = trackLength / trackVelocity * PicoToNano;
particle.setIndicesDaughter(allParticles.size(), allParticles.size() + (decayStack.size() - 1));
for (o2::upgrade::OTFParticle daughter : decayStack) {
daughter.setIndicesMother(i, i);
daughter.setCollisionId(mCollisionId);
particle.setIndicesDaughter(particlesInDataframe - indexOffset + allParticles.size(), particlesInDataframe - indexOffset + allParticles.size() + (decayStack.size() - 1));
for (auto& daughter : decayStack) {

Check failure on line 163 in ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.
daughter.setIndicesMother(particlesInDataframe - indexOffset + i, particlesInDataframe - indexOffset + i);
daughter.setCollisionId(particle.collisionId());
daughter.setBitOn(o2::upgrade::DecayerBits::IsAlive);
daughter.setBitOff(o2::upgrade::DecayerBits::IsPrimary);
daughter.setProductionTime(trackTimeNS);
allParticles.push_back(daughter);
ndau++;
}
}

if (start >= stop) {
return;
ndau += decayStack.size();
}

decayParticles(stop, stop + ndau);
Expand All @@ -173,18 +177,16 @@
void process(aod::McCollisions_001From<aod::Hash<"TMP"_h>>::iterator const& collision, aod::McParticles_001From<aod::Hash<"TMP"_h>> const& mcParticles)
{
allParticles.clear();
allParticles.reserve(mcParticles.size() * 2);
if (collision.globalIndex() == 0) {
particlesInDataframe = 0;
indexOffset = 0;
}

// Reproduce collision table to have AOD origin
mCollisionId = collision.globalIndex();
tableMcCollisions(collision.bcId(),
collision.generatorsID(),
collision.posX(),
collision.posY(),
collision.posZ(),
collision.t(),
collision.weight(),
collision.impactParameter(),
collision.eventPlaneAngle());
tableMcCollisions(collision.bcId(), collision.generatorsID(),
collision.posX(), collision.posY(), collision.posZ(), collision.t(),
collision.weight(), collision.impactParameter(), collision.eventPlaneAngle());

// First we copy the particles from the table into a vector that is extendable
for (const auto& particle : mcParticles) {
Expand All @@ -195,19 +197,26 @@
decayParticles(0, allParticles.size());

// Fill output table
for (const auto& otfParticle : allParticles) {
int id = particlesInDataframe;
for (auto& otfParticle : allParticles) {

Check failure on line 201 in ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.
otfParticle.setIndexOffset(indexOffset);
if (otfParticle.hasNaN()) {
histos.fill(HIST("hNaNBookkeeping"), 1);
} else {
histos.fill(HIST("hNaNBookkeeping"), 0);
}

tableOTFDecayerBits(otfParticle.getBitsValue());
tableMcParticles(otfParticle.collisionId(), otfParticle.pdgCode(), otfParticle.statusCode(), otfParticle.flags(),
tableMcParticles(tableMcCollisions.lastIndex(), otfParticle.pdgCode(), otfParticle.statusCode(), otfParticle.flags(),
otfParticle.getMotherSpan(), otfParticle.getDaughters().data(), otfParticle.weight(),
otfParticle.px(), otfParticle.py(), otfParticle.pz(), otfParticle.e(),
otfParticle.vx(), otfParticle.vy(), otfParticle.vz(), otfParticle.vt());
}

// Particles for later collisions in df's needs to have thier mother
// and daughter indices adjusted since their global index will be
// shifted due to the appending of decay products
indexOffset += (allParticles.size() - mcParticles.size());
}
};

Expand Down
17 changes: 8 additions & 9 deletions ALICE3/Tasks/alice3DecayerQa.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -132,7 +132,6 @@ struct Alice3DecayerQa {

void process(const aod::McCollision& collision, const aod::McParticles& particles)
{
LOG(info) << particles.size();
// Group with collision
auto trueElectronsGrouped = trueElectrons->sliceByCached(aod::mcparticle::mcCollisionId, collision.globalIndex(), cache);
auto trueMuonsGrouped = trueMuons->sliceByCached(aod::mcparticle::mcCollisionId, collision.globalIndex(), cache);
Expand Down Expand Up @@ -164,8 +163,8 @@ struct Alice3DecayerQa {
histos.fill(HIST("K0S/hHasDecayed"), 1);
auto daughters = particle.daughtersIds();
if (daughters.size() == NV0Daughters) {
auto dau0 = particles.rawIteratorAt(daughters.front());
auto dau1 = particles.rawIteratorAt(daughters.back());
auto dau0 = particles.rawIteratorAt(daughters.front() - particles.offset());
auto dau1 = particles.rawIteratorAt(daughters.back() - particles.offset());

// K0S -> pi+ pi-
const bool k0sDecay = (dau0.pdgCode() == PDG_t::kPiPlus && dau1.pdgCode() == PDG_t::kPiMinus) ||
Expand All @@ -188,8 +187,8 @@ struct Alice3DecayerQa {
histos.fill(HIST("Lambda/hHasDecayed"), 1);
auto daughters = particle.daughtersIds();
if (daughters.size() == NV0Daughters) {
auto dau0 = particles.rawIteratorAt(daughters[0]);
auto dau1 = particles.rawIteratorAt(daughters[1]);
auto dau0 = particles.rawIteratorAt(daughters[0] - particles.offset());
auto dau1 = particles.rawIteratorAt(daughters[1] - particles.offset());

// Lambda -> p pi-
const bool lambdaDecay = (dau0.pdgCode() == PDG_t::kProton && dau1.pdgCode() == PDG_t::kPiMinus) ||
Expand All @@ -212,8 +211,8 @@ struct Alice3DecayerQa {
histos.fill(HIST("XiMinus/hHasDecayed"), 1);
auto daughters = particle.daughtersIds();
if (daughters.size() == NCascadeDaughters) {
auto dau0 = particles.rawIteratorAt(daughters.front());
auto dau1 = particles.rawIteratorAt(daughters.back());
auto dau0 = particles.rawIteratorAt(daughters.front() - particles.offset());
auto dau1 = particles.rawIteratorAt(daughters.back() - particles.offset());

// Xi- -> Lambda pi-
const bool xiDecay = (dau0.pdgCode() == PDG_t::kLambda0 && dau1.pdgCode() == PDG_t::kPiMinus) ||
Expand All @@ -229,8 +228,8 @@ struct Alice3DecayerQa {
if (v0.has_daughters()) {
auto v0daughters = v0.daughtersIds();
if (v0daughters.size() == NV0Daughters) {
auto v0dau0 = particles.rawIteratorAt(v0daughters.front());
auto v0dau1 = particles.rawIteratorAt(v0daughters.back());
auto v0dau0 = particles.rawIteratorAt(v0daughters.front() - particles.offset());
auto v0dau1 = particles.rawIteratorAt(v0daughters.back() - particles.offset());
const bool lambdaDecay = (v0dau0.pdgCode() == PDG_t::kProton && v0dau1.pdgCode() == PDG_t::kPiMinus) ||
(v0dau0.pdgCode() == PDG_t::kPiMinus && v0dau1.pdgCode() == PDG_t::kProton);
if (lambdaDecay) {
Expand Down
Loading