Skip to content

Commit 8c0e826

Browse files
jesgumjesgum
authored
[ALICE3] Fix indices bug in otf decayer (#16587)
1 parent 21a8426 commit 8c0e826

3 files changed

Lines changed: 58 additions & 38 deletions

File tree

‎ALICE3/Core/OTFParticle.h‎

Lines changed: 14 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -58,6 +58,9 @@ class OTFParticle
5858
if (particle.has_mothers()) {
5959
mIndicesMother = {particle.mothersIds().front(), particle.mothersIds().back()};
6060
}
61+
if (particle.has_daughters()) {
62+
mIndicesDaughter = {particle.daughtersIds().front(), particle.daughtersIds().back()};
63+
}
6164
if constexpr (requires { particle.decayerBits(); }) {
6265
mBits = particle.decayerBits();
6366
} else {
@@ -88,6 +91,14 @@ class OTFParticle
8891
mPz = pz;
8992
mE = e;
9093
}
94+
void setIndexOffset(const std::size_t offset)
95+
{
96+
static constexpr int NotFound = -1;
97+
mIndicesMother[0] = (mIndicesMother[0] >= 0) ? mIndicesMother[0] + static_cast<int>(offset) : NotFound;
98+
mIndicesMother[1] = (mIndicesMother[1] >= 0) ? mIndicesMother[1] + static_cast<int>(offset) : NotFound;
99+
mIndicesDaughter[0] = (mIndicesDaughter[0] >= 0) ? mIndicesDaughter[0] + static_cast<int>(offset) : NotFound;
100+
mIndicesDaughter[1] = (mIndicesDaughter[1] >= 0) ? mIndicesDaughter[1] + static_cast<int>(offset) : NotFound;
101+
}
91102

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

149160
// Checks
150-
bool hasDaughters() const { return (mIndicesDaughter[0] > 0); }
151-
bool hasMothers() const { return (mIndicesMother[0] > 0); }
161+
bool hasDaughters() const { return (mIndicesDaughter[0] >= 0); }
162+
bool hasMothers() const { return (mIndicesMother[0] >= 0); }
152163
bool hasNaN() const
153164
{
154165
return std::isnan(mPx) || std::isnan(mPy) || std::isnan(mPz) || std::isnan(mE) ||
@@ -171,7 +182,7 @@ class OTFParticle
171182

172183
private:
173184
int mPdgCode{}, mGlobalIndex{-1};
174-
int mCollisionId{};
185+
int mCollisionId{-1};
175186
float mVx{}, mVy{}, mVz{}, mVt{};
176187
float mPx{}, mPy{}, mPz{}, mE{};
177188
bool mIsAlive{}, mIsFromMcParticles{false};

‎ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx‎

Lines changed: 36 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -49,8 +49,8 @@ using namespace o2::framework;
4949
static constexpr int NumDecays = 7;
5050
static constexpr int NumParameters = 1;
5151
static constexpr int DefaultParameters[NumDecays][NumParameters]{{1}, {1}, {1}, {1}, {1}, {1}, {1}};
52-
static const std::vector<std::string> ParameterNames{"enable"};
53-
static const std::vector<std::string> ParticleNames{"K0s",
52+
static const std::vector<std::string> parameterNames{"enable"};
53+
static const std::vector<std::string> particleNames{"K0s",
5454
"Lambda",
5555
"Anti-Lambda",
5656
"Xi",
@@ -83,12 +83,12 @@ struct OnTheFlyDecayer {
8383
Configurable<int> seed{"seed", 0, "Set seed for particle decayer"};
8484
Configurable<float> magneticField{"magneticField", 20., "Magnetic field (kG)"};
8585
Configurable<LabeledArray<int>> enabledDecays{"enabledDecays",
86-
{DefaultParameters[0], NumDecays, NumParameters, ParticleNames, ParameterNames},
86+
{DefaultParameters[0], NumDecays, NumParameters, particleNames, parameterNames},
8787
"Enable option for particle to be decayed: 0 - no, 1 - yes"};
8888

89+
std::size_t indexOffset = 0;
90+
std::size_t particlesInDataframe = 0;
8991
static constexpr float PicoToNano = 1.e-3f;
90-
int mCollisionId{-1};
91-
9292
std::vector<int> mEnabledDecays;
9393
void init(o2::framework::InitContext&)
9494
{
@@ -98,7 +98,7 @@ struct OnTheFlyDecayer {
9898
decayer.setSeed(seed);
9999
decayer.setBField(magneticField);
100100
for (int i = 0; i < NumDecays; ++i) {
101-
if (enabledDecays->get(ParticleNames[i].c_str(), "enable")) {
101+
if (enabledDecays->get(particleNames[i].c_str(), "enable")) {
102102
LOG(info) << " --- Decay enabled: " << pdgCodes[i];
103103
mEnabledDecays.push_back(pdgCodes[i]);
104104
}
@@ -121,6 +121,10 @@ struct OnTheFlyDecayer {
121121
std::vector<o2::upgrade::OTFParticle> allParticles;
122122
void decayParticles(const int start, const int stop)
123123
{
124+
if (start >= stop) {
125+
return;
126+
}
127+
124128
int ndau = 0;
125129
for (int i = start; i < stop; ++i) {
126130
o2::upgrade::OTFParticle& particle = allParticles[i];
@@ -135,6 +139,10 @@ struct OnTheFlyDecayer {
135139

136140
particle.setBitOff(o2::upgrade::DecayerBits::IsAlive);
137141
std::vector<o2::upgrade::OTFParticle> decayStack = decayer.decayParticle(pdgDB, particle);
142+
if (decayStack.empty()) {
143+
continue;
144+
}
145+
138146
const float decayRadius = decayer.getDecayRadius();
139147
const float trackVelocity = o2::upgrade::computeParticleVelocity(particle.p(), pdgDB->GetParticle(particle.pdgCode())->Mass());
140148
const int charge = pdgDB->GetParticle(particle.pdgCode())->Charge() / 3;
@@ -151,20 +159,16 @@ struct OnTheFlyDecayer {
151159
}
152160

153161
const float trackTimeNS = trackLength / trackVelocity * PicoToNano;
154-
particle.setIndicesDaughter(allParticles.size(), allParticles.size() + (decayStack.size() - 1));
155-
for (o2::upgrade::OTFParticle daughter : decayStack) {
156-
daughter.setIndicesMother(i, i);
157-
daughter.setCollisionId(mCollisionId);
162+
particle.setIndicesDaughter(particlesInDataframe - indexOffset + allParticles.size(), particlesInDataframe - indexOffset + allParticles.size() + (decayStack.size() - 1));
163+
for (auto& daughter : decayStack) {
164+
daughter.setIndicesMother(particlesInDataframe - indexOffset + i, particlesInDataframe - indexOffset + i);
165+
daughter.setCollisionId(particle.collisionId());
158166
daughter.setBitOn(o2::upgrade::DecayerBits::IsAlive);
159167
daughter.setBitOff(o2::upgrade::DecayerBits::IsPrimary);
160168
daughter.setProductionTime(trackTimeNS);
161169
allParticles.push_back(daughter);
162-
ndau++;
163170
}
164-
}
165-
166-
if (start >= stop) {
167-
return;
171+
ndau += decayStack.size();
168172
}
169173

170174
decayParticles(stop, stop + ndau);
@@ -173,18 +177,16 @@ struct OnTheFlyDecayer {
173177
void process(aod::McCollisions_001From<aod::Hash<"TMP"_h>>::iterator const& collision, aod::McParticles_001From<aod::Hash<"TMP"_h>> const& mcParticles)
174178
{
175179
allParticles.clear();
180+
allParticles.reserve(mcParticles.size() * 2);
181+
if (collision.globalIndex() == 0) {
182+
particlesInDataframe = 0;
183+
indexOffset = 0;
184+
}
176185

177186
// Reproduce collision table to have AOD origin
178-
mCollisionId = collision.globalIndex();
179-
tableMcCollisions(collision.bcId(),
180-
collision.generatorsID(),
181-
collision.posX(),
182-
collision.posY(),
183-
collision.posZ(),
184-
collision.t(),
185-
collision.weight(),
186-
collision.impactParameter(),
187-
collision.eventPlaneAngle());
187+
tableMcCollisions(collision.bcId(), collision.generatorsID(),
188+
collision.posX(), collision.posY(), collision.posZ(), collision.t(),
189+
collision.weight(), collision.impactParameter(), collision.eventPlaneAngle());
188190

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

197199
// Fill output table
198-
for (const auto& otfParticle : allParticles) {
200+
for (auto& otfParticle : allParticles) {
201+
otfParticle.setIndexOffset(indexOffset);
199202
if (otfParticle.hasNaN()) {
200203
histos.fill(HIST("hNaNBookkeeping"), 1);
201204
} else {
202205
histos.fill(HIST("hNaNBookkeeping"), 0);
203206
}
204207

205208
tableOTFDecayerBits(otfParticle.getBitsValue());
206-
tableMcParticles(otfParticle.collisionId(), otfParticle.pdgCode(), otfParticle.statusCode(), otfParticle.flags(),
209+
tableMcParticles(tableMcCollisions.lastIndex(), otfParticle.pdgCode(), otfParticle.statusCode(), otfParticle.flags(),
207210
otfParticle.getMotherSpan(), otfParticle.getDaughters().data(), otfParticle.weight(),
208211
otfParticle.px(), otfParticle.py(), otfParticle.pz(), otfParticle.e(),
209212
otfParticle.vx(), otfParticle.vy(), otfParticle.vz(), otfParticle.vt());
210213
}
214+
215+
// Particles for later collisions in df's needs to have thier mother
216+
// and daughter indices adjusted since their global index will be
217+
// shifted due to the appending of decay products
218+
indexOffset += (allParticles.size() - mcParticles.size());
219+
particlesInDataframe += allParticles.size();
211220
}
212221
};
213222

‎ALICE3/Tasks/alice3DecayerQa.cxx‎

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -163,8 +163,8 @@ struct Alice3DecayerQa {
163163
histos.fill(HIST("K0S/hHasDecayed"), 1);
164164
auto daughters = particle.daughtersIds();
165165
if (daughters.size() == NV0Daughters) {
166-
auto dau0 = particles.rawIteratorAt(daughters.front());
167-
auto dau1 = particles.rawIteratorAt(daughters.back());
166+
auto dau0 = particles.rawIteratorAt(daughters.front() - particles.offset());
167+
auto dau1 = particles.rawIteratorAt(daughters.back() - particles.offset());
168168

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

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

217217
// Xi- -> Lambda pi-
218218
const bool xiDecay = (dau0.pdgCode() == PDG_t::kLambda0 && dau1.pdgCode() == PDG_t::kPiMinus) ||
@@ -228,8 +228,8 @@ struct Alice3DecayerQa {
228228
if (v0.has_daughters()) {
229229
auto v0daughters = v0.daughtersIds();
230230
if (v0daughters.size() == NV0Daughters) {
231-
auto v0dau0 = particles.rawIteratorAt(v0daughters.front());
232-
auto v0dau1 = particles.rawIteratorAt(v0daughters.back());
231+
auto v0dau0 = particles.rawIteratorAt(v0daughters.front() - particles.offset());
232+
auto v0dau1 = particles.rawIteratorAt(v0daughters.back() - particles.offset());
233233
const bool lambdaDecay = (v0dau0.pdgCode() == PDG_t::kProton && v0dau1.pdgCode() == PDG_t::kPiMinus) ||
234234
(v0dau0.pdgCode() == PDG_t::kPiMinus && v0dau1.pdgCode() == PDG_t::kProton);
235235
if (lambdaDecay) {

0 commit comments

Comments
 (0)