From 79ef45abf79a83483492f5c76718f5311862075b Mon Sep 17 00:00:00 2001 From: jesgum Date: Fri, 2 Oct 2026 14:25:01 +0200 Subject: [PATCH 1/2] Add option to modify the decay table in the otf decayer --- ALICE3/Core/Decayer.h | 90 ++++++++++++++++++-- ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx | 7 ++ 2 files changed, 90 insertions(+), 7 deletions(-) diff --git a/ALICE3/Core/Decayer.h b/ALICE3/Core/Decayer.h index 8f88e4e311b..3a2059fde20 100644 --- a/ALICE3/Core/Decayer.h +++ b/ALICE3/Core/Decayer.h @@ -19,10 +19,13 @@ #ifndef ALICE3_CORE_DECAYER_H_ #define ALICE3_CORE_DECAYER_H_ +#include "ALICE3/Core/ConfigurationParser.h" #include "ALICE3/Core/OTFParticle.h" #include "ALICE3/Core/TrackUtilities.h" +#include #include +#include #include #include @@ -34,6 +37,10 @@ #include #include #include +#include +#include +#include +#include #include namespace o2::upgrade @@ -45,6 +52,53 @@ class Decayer // Default constructor Decayer() = default; + static constexpr double BranchingRatioTolerance = 0.01; + + struct DecayChannel { + double branchingRatio{}; + std::vector daughters; + }; + + template + void loadDecayTable(const std::string& path, const TDatabase& pdgDB, o2::ccdb::BasicCCDBManager* ccdb = nullptr) + { + std::string fileName = o2::fastsim::ConfigurationParser::accessFile(path, "./.ALICE3/Configuration/", ccdb); + std::vector mothers; + const std::map> config = o2::fastsim::ConfigurationParser::parseTEnvConfiguration(fileName, mothers); + + mDecayTable.clear(); + for (const auto& [motherName, channels] : config) { + const int motherPdg = std::stoi(motherName); + if (!pdgDB->GetParticle(motherPdg)) { + LOG(fatal) << "Decay table: unknown mother PDG code " << motherPdg; + } + for (const auto& [channelName, value] : channels) { + DecayChannel channel; + std::istringstream stream(value); + stream >> channel.branchingRatio; + int daughterPdg{}; + while (stream >> daughterPdg) { + if (!pdgDB->GetParticle(daughterPdg)) { + LOG(fatal) << "Decay table: unknown daughter PDG code " << daughterPdg << " in " << motherName << "." << channelName; + } + channel.daughters.push_back(daughterPdg); + } + if (!stream.eof() || channel.branchingRatio <= 0. || channel.daughters.size() < 2) { + LOG(fatal) << "Decay table: invalid channel " << motherName << "." << channelName << ": \"" << value << "\""; + } + mDecayTable[motherPdg].push_back(channel); + } + double brTotal = 0.; + for (const auto& channel : mDecayTable[motherPdg]) { + brTotal += channel.branchingRatio; + } + if (std::abs(brTotal - 1.) > BranchingRatioTolerance) { + LOG(fatal) << "Decay table: branching ratios of PDG " << motherPdg << " sum to " << brTotal << ", expected 1 within " << BranchingRatioTolerance; + } + LOG(info) << "Decay table: overriding " << mDecayTable[motherPdg].size() << " decay channel(s) of PDG " << motherPdg; + } + } + template std::vector decayParticle(const OTFParticle& particle, const TDatabase& pdgDB) { @@ -69,21 +123,22 @@ class Decayer py = particle.py() * std::cos(mTheta) + particle.px() * std::sin(mTheta); } - double brTotal = 0.; e = std::sqrt(mass * mass + px * px + py * py + particle.pz() * particle.pz()); - for (int ch = 0; ch < particleInfo->NDecayChannels(); ++ch) { - brTotal += particleInfo->DecayChannel(ch)->BranchingRatio(); + const std::vector channels = getDecayChannels(particle.pdgCode(), particleInfo); + + double brTotal = 0.; + for (const auto& channel : channels) { + brTotal += channel.branchingRatio; } double brSum = 0.; std::vector dauMasses; std::vector pdgCodesDaughters; const double randomChannel = mRand3.Uniform(0., brTotal); - for (int ch = 0; ch < particleInfo->NDecayChannels(); ++ch) { - brSum += particleInfo->DecayChannel(ch)->BranchingRatio(); + for (const auto& channel : channels) { + brSum += channel.branchingRatio; if (randomChannel < brSum) { - for (int dau = 0; dau < particleInfo->DecayChannel(ch)->NDaughters(); ++dau) { - const int pdgDau = particleInfo->DecayChannel(ch)->DaughterPdgCode(dau); + for (const int pdgDau : channel.daughters) { pdgCodesDaughters.push_back(pdgDau); const auto& dauInfo = pdgDB->GetParticle(pdgDau); dauMasses.push_back(dauInfo->Mass()); @@ -166,6 +221,27 @@ class Decayer [[nodiscard]] float getDecayRadius() const { return static_cast(std::hypot(mVx, mVy)); } private: + /// Decay channels from the loaded decay table if the mother is listed there, otherwise from TDatabasePDG + template + std::vector getDecayChannels(const int pdgCode, const TParticleInfo& particleInfo) const + { + if (const auto it = mDecayTable.find(pdgCode); it != mDecayTable.end()) { + return it->second; + } + + std::vector channels; + for (int ch = 0; ch < particleInfo->NDecayChannels(); ++ch) { + DecayChannel channel; + channel.branchingRatio = particleInfo->DecayChannel(ch)->BranchingRatio(); + for (int dau = 0; dau < particleInfo->DecayChannel(ch)->NDaughters(); ++dau) { + channel.daughters.push_back(particleInfo->DecayChannel(ch)->DaughterPdgCode(dau)); + } + channels.push_back(channel); + } + return channels; + } + + std::unordered_map> mDecayTable; double mBz{20.}; // kG double mVx{-1.}, mVy{-1.}, mVz{-1.}; double mTheta{}; diff --git a/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx b/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx index 2c8c4dea0c9..c6f587fdd3c 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx @@ -20,6 +20,7 @@ #include "ALICE3/Core/TrackUtilities.h" #include "ALICE3/DataModel/tracksAlice3.h" +#include #include #include #include @@ -92,10 +93,12 @@ struct OnTheFlyDecayer { o2::upgrade::Decayer decayer; Service pdgDB{}; + Service ccdb{}; HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; Configurable seed{"seed", 0, "Set seed for particle decayer"}; Configurable magneticField{"magneticField", 20., "Magnetic field (kG)"}; + Configurable decayTable{"decayTable", "", "Decay table overriding the TDatabasePDG decay channels of the listed particles (local path or ccdb:), empty to disable"}; Configurable> enabledDecays{"enabledDecays", {DefaultParameters[0].data(), NumDecays, NumParameters, particleNames, parameterNames}, "Enable option for particle to be decayed: 0 - no, 1 - yes"}; @@ -111,6 +114,10 @@ struct OnTheFlyDecayer { LOG(info) << " --- Using magnetic field: " << magneticField; decayer.setSeed(seed); decayer.setBField(magneticField); + if (!decayTable.value.empty()) { + LOG(info) << " --- Using decay table: " << decayTable.value; + decayer.loadDecayTable(decayTable.value, pdgDB, ccdb.operator->()); + } for (int i = 0; i < NumDecays; ++i) { if (enabledDecays->get(particleNames[i].c_str(), "enable") != 0) { LOG(info) << " --- Decay enabled: " << pdgCodes[i]; From c76e6db8109e09435108a89e1d6863483a41ed9f Mon Sep 17 00:00:00 2001 From: jesgum Date: Fri, 2 Oct 2026 14:40:24 +0200 Subject: [PATCH 2/2] u --- ALICE3/Core/Decayer.h | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/ALICE3/Core/Decayer.h b/ALICE3/Core/Decayer.h index 3a2059fde20..fb39a372cb2 100644 --- a/ALICE3/Core/Decayer.h +++ b/ALICE3/Core/Decayer.h @@ -157,16 +157,16 @@ class Decayer decay.Generate(); std::vector decayProducts; + decayProducts.reserve(dauMasses.size()); for (size_t i = 0; i < dauMasses.size(); ++i) { - o2::upgrade::OTFParticle particle; - TLorentzVector dau = *decay.GetDecay(i); - particle.setPDG(pdgCodesDaughters[i]); - particle.setVxVyVz(mVx, mVy, mVz); - particle.setPxPyPzE(dau.Px(), dau.Py(), dau.Pz(), dau.E()); - particle.setBitOn(o2::upgrade::DecayerBits::ProducedByDecayer); - decayProducts.push_back(particle); + o2::upgrade::OTFParticle daughter; + const TLorentzVector& dau = *decay.GetDecay(i); + daughter.setPDG(pdgCodesDaughters[i]); + daughter.setVxVyVz(mVx, mVy, mVz); + daughter.setPxPyPzE(dau.Px(), dau.Py(), dau.Pz(), dau.E()); + daughter.setBitOn(o2::upgrade::DecayerBits::ProducedByDecayer); + decayProducts.push_back(daughter); } - return decayProducts; }