From 91e94fcdbe061f49231bf05cd599eb32338b50d7 Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Thu, 22 Oct 2026 13:09:04 +0200 Subject: [PATCH 1/5] Add a tool that writes the Geant4 particle table This commit provides a new tool o2-sim-g4-particle-table that writes the Geant4 particle table to a text file. - One line per particle with a PDG code: pdg, Geant4 name, ctau in mm, short-lived flag, decay-table flag. - The Geant4 setup is the same as in o2-sim-g4-determine-unknown-pdg-properties. - Particles that o2-sim adds at run time are not in the table. Usage: o2-sim-g4-particle-table --output g4_particle_table.dat https://its.cern.ch/jira/browse/O2-5903 Co-Authored-By: Claude Opus 5.5 --- run/CMakeLists.txt | 6 +++ run/o2sim_g4_particle_table.cxx | 77 +++++++++++++++++++++++++++++++++ 2 files changed, 83 insertions(+) create mode 100644 run/o2sim_g4_particle_table.cxx diff --git a/run/CMakeLists.txt b/run/CMakeLists.txt index 8886bfcbaf1cc..5fda996400ed2 100644 --- a/run/CMakeLists.txt +++ b/run/CMakeLists.txt @@ -113,6 +113,12 @@ o2_add_executable(g4-determine-unknown-pdg-properties SOURCES g4DetermineUnknownPdgProperties.cxx PUBLIC_LINK_LIBRARIES O2::O2TrivialMC MC::Geant4 MC::Geant4VMC O2::SimConfig) +# Writes the Geant4 particle table to a text file. +o2_add_executable(g4-particle-table + COMPONENT_NAME sim + SOURCES o2sim_g4_particle_table.cxx + PUBLIC_LINK_LIBRARIES O2::O2TrivialMC MC::Geant4 MC::Geant4VMC O2::SimConfig) + # Reads a geometry file and a field and reports where the two do not fit together. # It needs neither the simulation engines nor any detector code -- the medium flags # it reasons about are GSTMED parameters stored in the geometry file itself -- so it diff --git a/run/o2sim_g4_particle_table.cxx b/run/o2sim_g4_particle_table.cxx new file mode 100644 index 0000000000000..fc205d92cf8a8 --- /dev/null +++ b/run/o2sim_g4_particle_table.cxx @@ -0,0 +1,77 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +// Writes the Geant4 particle table (one line per particle with a PDG code) to a text file. + +#include +#include +#include +#include + +#include "G4ParticleTable.hh" +#include "G4ParticleDefinition.hh" +#include "G4Version.hh" +#include "CLHEP/Units/PhysicalConstants.h" + +#include "TGeant4.h" +#include "TG4RunConfiguration.h" + +#include "SimConfig/G4Params.h" +#include "O2TrivialMC/O2TrivialMCApplication.h" + +int main(int argc, char** argv) +{ + if (argc != 3 || std::string(argv[1]) != "--output") { + std::cerr << "Usage: " << argv[0] << " --output \n"; + return 1; + } + + std::ofstream os{argv[2]}; + if (!os.is_open()) { + std::cerr << "Cannot open " << argv[2] << "\n"; + return 1; + } + + // We use our O2 VMC setup to make sure we are in line with our physics definition. + new o2::mc::O2TrivialMCApplication(); + auto& physicsSetup = ::o2::conf::G4Params::Instance().getPhysicsConfigString(); + auto runConfiguration = new TG4RunConfiguration("geomRoot", physicsSetup); + auto vmc = new TGeant4("TGeant4", "The Geant4 Monte Carlo", runConfiguration); + // needs to be initialised to have the particle table / physics entirely ready + vmc->Init(); + + os << "# Geant4 particle table, Geant4 version " << G4VERSION_NUMBER << "\n"; + os << "# physics setup: " << physicsSetup << "\n"; + os << "# pdg g4name ctau_mm shortlived hasDecayTable\n"; + + auto table = G4ParticleTable::GetParticleTable(); + auto it = table->GetIterator(); + it->reset(); + int nParticles = 0; + while ((*it)()) { + auto particle = it->value(); + int pdg = particle->GetPDGEncoding(); + if (pdg == 0) { + continue; + } + double lifetime = particle->GetPDGLifeTime(); + // lifetime is in G4 internal units, times c_light gives mm + double ctau = lifetime > 0. ? lifetime * CLHEP::c_light : 0.; + char buf[64]; + std::snprintf(buf, sizeof(buf), "%.6g", ctau); + os << pdg << " " << particle->GetParticleName() << " " << buf << " " + << (particle->IsShortLived() ? 1 : 0) << " " + << (particle->GetDecayTable() ? 1 : 0) << "\n"; + ++nParticles; + } + std::cout << "Wrote " << nParticles << " particles to " << argv[2] << "\n"; + return 0; +} From 3e0d632d189f79c4b7012103e07fcb0417817699 Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Thu, 22 Oct 2026 13:09:04 +0200 Subject: [PATCH 2/5] Add an event decayer that decays generator particles before transport This commit provides a class EventDecayer that decays unstable generator particles with Pythia8 before transport and keeps the full decay chain in the MC truth. - It works in place on a std::vector, so that a generator can call it directly. - Only primaries flagged for transport are decayed; the products are appended with their displaced vertices. - Decayed particles keep their generator status and get HepMC status 2. - By default it reads the DecayerPythia8 configs and never decays species above tau0Max (10 mm/c) or listed in skipPdg. - Particles left for transport are checked against O2DatabasePDG and the Geant4 particle table. - The settings come from the new param EventDecayer. https://its.cern.ch/jira/browse/O2-5903 Co-Authored-By: Claude Opus 5.5 --- Generators/CMakeLists.txt | 12 + Generators/include/Generators/EventDecayer.h | 90 +++ .../include/Generators/EventDecayerParam.h | 39 + Generators/src/EventDecayer.cxx | 403 ++++++++++ Generators/src/EventDecayerParam.cxx | 13 + Generators/src/GeneratorsLinkDef.h | 2 + Generators/test/test_EventDecayer.cxx | 698 ++++++++++++++++++ 7 files changed, 1257 insertions(+) create mode 100644 Generators/include/Generators/EventDecayer.h create mode 100644 Generators/include/Generators/EventDecayerParam.h create mode 100644 Generators/src/EventDecayer.cxx create mode 100644 Generators/src/EventDecayerParam.cxx create mode 100644 Generators/test/test_EventDecayer.cxx diff --git a/Generators/CMakeLists.txt b/Generators/CMakeLists.txt index 094510c2c61a5..52f4f5968698a 100644 --- a/Generators/CMakeLists.txt +++ b/Generators/CMakeLists.txt @@ -47,6 +47,8 @@ o2_add_library(Generators $<$:src/DecayerPythia8.cxx> $<$:src/GeneratorPythia8Param.cxx> $<$:src/DecayerPythia8Param.cxx> + $<$:src/EventDecayer.cxx> + $<$:src/EventDecayerParam.cxx> $<$:src/GeneratorHepMC.cxx> $<$:src/GeneratorHepMCParam.cxx> $<$:src/AODToHepMC.cxx> @@ -103,6 +105,7 @@ if(pythia_FOUND) include/Generators/DecayerPythia8.h include/Generators/GeneratorPythia8Param.h include/Generators/DecayerPythia8Param.h + include/Generators/EventDecayerParam.h include/Generators/GeneratorFactory.h) endif() @@ -147,6 +150,15 @@ if(doBuildSimulation) PUBLIC_LINK_LIBRARIES O2::Generators) endif() + if(pythia_FOUND) + o2_add_test(EventDecayer NAME test_Generator_test_EventDecayer + SOURCES test/test_EventDecayer.cxx + COMPONENT_NAME Generator + ENVIRONMENT O2_ROOT=${CMAKE_BINARY_DIR}/stage + LABELS generator + PUBLIC_LINK_LIBRARIES O2::Generators) + endif() + # o2_add_test(GeneratorPythia8Param NAME test_Generator_test_GeneratorPythia8Param # SOURCES test/test_GeneratorPythia8Param.cxx # COMPONENT_NAME Generator diff --git a/Generators/include/Generators/EventDecayer.h b/Generators/include/Generators/EventDecayer.h new file mode 100644 index 0000000000000..625db7ce8f582 --- /dev/null +++ b/Generators/include/Generators/EventDecayer.h @@ -0,0 +1,90 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef ALICEO2_EVENTGEN_EVENTDECAYER_H_ +#define ALICEO2_EVENTGEN_EVENTDECAYER_H_ + +#include "Pythia8/Pythia.h" +#include +#include +#include +#include +#include +#include + +namespace o2 +{ +namespace eventgen +{ + +/// Counts of what one or more EventDecayer::decay calls did, per PDG code +struct EventDecayerReport { + std::map decayed; // particles decayed by Pythia8 + std::map failed; // particles Pythia8 could not decay + std::map unknownInO2; // particles to transport unknown to O2DatabasePDG + std::map unknownInG4; // particles to transport unknown to Geant4 + std::map> lifetimeMismatch; // pdg -> {pythia ctau mm, g4 ctau mm} + void merge(EventDecayerReport const& other); + void print() const; +}; + +/// One entry of the Geant4 particle table +struct G4ParticleInfo { + std::string name; + double ctauMM; + bool shortLived; + bool hasDecayTable; +}; + +/// Reads the Geant4 particle table written by o2-sim-g4-particle-table +std::unordered_map loadG4ParticleTable(std::string const& path); + +/// Decays unstable generator particles with Pythia8 before transport, keeping the full chain +class EventDecayer +{ + public: + bool init(); + /// Sets the caller's units in cm, s, GeV/c and GeV, as the unit setters of Generator (default cm, s, GeV/c, GeV) + void setUnits(double positionUnit, double timeUnit, double momentumUnit = 1., double energyUnit = 1.) + { + mPositionUnit = positionUnit; + mTimeUnit = timeUnit; + mMomentumUnit = momentumUnit; + mEnergyUnit = energyUnit; + } + EventDecayerReport decay(std::vector& particles); + std::string configSummary() const; + bool vertexLimitsActive() const; + int seed() const { return mSeed; } + + private: + bool decayable(int pdg) const; + bool decayOne(std::vector& particles, int index); + void checkParticles(std::vector const& particles, EventDecayerReport& report) const; + void checkLifetimes(EventDecayerReport& report) const; + + mutable ::Pythia8::Pythia mPythia; // mutable since Pythia8 settings lookups are not const + std::vector mConfigs; // expanded config paths + std::unordered_map mG4Particles; // Geant4 particle table + bool mHasG4Table = false; + bool mVerbose = false; + double mPositionUnit = 1.; // caller's position unit in cm + double mTimeUnit = 1.; // caller's time unit in s + double mMomentumUnit = 1.; // caller's momentum unit in GeV/c + double mEnergyUnit = 1.; // caller's energy unit in GeV + double mTau0Max = 0.; // tau0 limit in mm/c in use after init, 0 when off + int mSeed = 0; // Pythia8 seed in use after init +}; + +} // namespace eventgen +} // namespace o2 + +#endif // ALICEO2_EVENTGEN_EVENTDECAYER_H_ diff --git a/Generators/include/Generators/EventDecayerParam.h b/Generators/include/Generators/EventDecayerParam.h new file mode 100644 index 0000000000000..5962405b5cd91 --- /dev/null +++ b/Generators/include/Generators/EventDecayerParam.h @@ -0,0 +1,39 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef ALICEO2_EVENTGEN_EVENTDECAYERPARAM_H_ +#define ALICEO2_EVENTGEN_EVENTDECAYERPARAM_H_ + +#include "CommonUtils/ConfigurableParam.h" +#include "CommonUtils/ConfigurableParamHelper.h" +#include + +namespace o2 +{ +namespace eventgen +{ + +/// Settings of the EventDecayer, which decays generator particles with Pythia8 before transport +struct EventDecayerParam : public o2::conf::ConfigurableParamHelper { + std::string config[8] = {"", "", "", "", "", "", "", ""}; // Pythia8 config files, an empty slot takes DecayerPythia8.config + double tau0Max = 10.; // species with a larger tau0 in mm/c stay undecayed (<= 0 switches it off) + std::string skipPdg = "115 225 335 10221 10331 20223 20333 102132 102134 9010221"; // PDG codes never decayed, matched on |pdg| + std::string g4ParticleTable = ""; // Geant4 particle table from o2-sim-g4-particle-table + int seed = 0; // Pythia8 random seed (0 derives it from gRandom) + bool verbose = false; // list every Pythia8 decay + bool showChanged = false; // show changed particle data at init + O2ParamDef(EventDecayerParam, "EventDecayer"); +}; + +} // end namespace eventgen +} // end namespace o2 + +#endif // ALICEO2_EVENTGEN_EVENTDECAYERPARAM_H_ diff --git a/Generators/src/EventDecayer.cxx b/Generators/src/EventDecayer.cxx new file mode 100644 index 0000000000000..2af5717c14da9 --- /dev/null +++ b/Generators/src/EventDecayer.cxx @@ -0,0 +1,403 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include "Generators/EventDecayer.h" +#include "Generators/EventDecayerParam.h" +#include "Generators/DecayerPythia8Param.h" +#include "SimulationDataFormat/MCGenProperties.h" +#include "SimulationDataFormat/MCUtils.h" +#include "SimulationDataFormat/O2DatabasePDG.h" +#include "SimulationDataFormat/ParticleStatus.h" +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace o2 +{ +namespace eventgen +{ + +namespace +{ +// speed of light in mm/s, to convert Pythia times (mm/c) to seconds +constexpr double kSpeedOfLightMMPerS = 2.99792458e11; + +// largest seed Pythia8 accepts for Random:seed +constexpr long kMaxSeed = 900000000; + +// mixed into a seed taken from gRandom, so that it differs from the one of GeneratorPythia8 +constexpr unsigned int kGRandomSeedXor = 0x5eed0d3c; + +// Pythia settings that cut decays on the decay vertex, not supported here +const char* const kVertexLimits[] = {"ParticleDecays:limitTau", "ParticleDecays:limitRadius", "ParticleDecays:limitCylinder"}; + +std::string expand(std::string const& path) +{ + TString expanded = path; + gSystem->ExpandPathName(expanded); + return expanded.Data(); +} + +// sets the Pythia8 seed as GeneratorPythia8::seedGenerator, a seed of 0 derives it from gRandom +bool setSeed(::Pythia8::Pythia& pythia, long seed) +{ + if (seed < 0) { + LOG(error) << "Failed to init \'EventDecayer\': negative seed " << seed; + return false; + } + if (seed == 0) { + seed = gRandom->TRandom::GetSeed() ^ kGRandomSeedXor; // from the seed set by o2-sim --seed + LOG(info) << "EventDecayer: using random seed from gRandom ^ " << kGRandomSeedXor << " % " << kMaxSeed + 1; + } + seed = seed % (kMaxSeed + 1); + if (!pythia.readString("Random:setSeed on") || !pythia.readString("Random:seed " + std::to_string(seed))) { + LOG(error) << "Failed to init \'EventDecayer\': cannot set the Pythia8 seed " << seed; + return false; + } + return true; +} + +// parses one skipPdg entry into |pdg|, false if it is not an integer +bool parsePdg(std::string const& code, int& pdg) +{ + try { + size_t used = 0; + pdg = std::abs(std::stoi(code, &used)); + return used == code.size(); + } catch (std::exception const&) { + return false; + } +} + +// switches off the decays of the species listed in skipPdg +bool skipDecays(::Pythia8::ParticleData& particleData, std::string const& skipPdg) +{ + std::istringstream skip(skipPdg); + std::string code; + while (skip >> code) { + int pdg = 0; + if (!parsePdg(code, pdg)) { + LOG(error) << "Failed to init \'EventDecayer\': skipPdg entry " << code << " is not a PDG code"; + return false; + } + if (particleData.isParticle(pdg)) { + particleData.mayDecay(pdg, false); + } else { + LOG(info) << "EventDecayer: skipPdg entry " << pdg << " is unknown to Pythia8, so never decayed anyway"; + } + } + return true; +} +} // namespace + +/*****************************************************************/ + +void EventDecayerReport::merge(EventDecayerReport const& other) +{ + auto add = [](std::map& to, std::map const& from) { + for (auto const& [pdg, count] : from) { + to[pdg] += count; + } + }; + add(decayed, other.decayed); + add(failed, other.failed); + add(unknownInO2, other.unknownInO2); + add(unknownInG4, other.unknownInG4); + for (auto const& [pdg, lifetimes] : other.lifetimeMismatch) { + lifetimeMismatch[pdg] = lifetimes; + } +} + +void EventDecayerReport::print() const +{ + for (auto const& [pdg, count] : decayed) { + LOG(info) << "EventDecayer: decayed pdg " << pdg << ": " << count; + } + for (auto const& [pdg, count] : failed) { + LOG(warn) << "EventDecayer: Pythia8 failed to decay pdg " << pdg << ": " << count; + } + for (auto const& [pdg, count] : unknownInO2) { + LOG(warn) << "EventDecayer: pdg " << pdg << " left for transport is unknown to O2DatabasePDG: " << count; + } + for (auto const& [pdg, count] : unknownInG4) { + LOG(warn) << "EventDecayer: pdg " << pdg << " left for transport is not in the Geant4 table (o2-sim may define it at runtime): " << count; + } + for (auto const& [pdg, lifetimes] : lifetimeMismatch) { + LOG(warn) << "EventDecayer: lifetime of pdg " << pdg << " differs between Pythia8 (ctau " << lifetimes.first + << " mm) and Geant4 (ctau " << lifetimes.second << " mm)"; + } +} + +/*****************************************************************/ + +std::unordered_map loadG4ParticleTable(std::string const& path) +{ + std::unordered_map table; + std::ifstream in(path); + if (!in) { + LOG(error) << "EventDecayer: cannot read Geant4 particle table " << path; + return table; + } + std::string line; + while (std::getline(in, line)) { + if (line.empty() || line[0] == '#') { + continue; + } + std::istringstream fields(line); + int pdg = 0; + G4ParticleInfo info; + if (!(fields >> pdg >> info.name >> info.ctauMM >> info.shortLived >> info.hasDecayTable)) { + LOG(warn) << "EventDecayer: skipping malformed line in " << path << ": " << line; + continue; + } + table[pdg] = info; + } + return table; +} + +/*****************************************************************/ + +bool EventDecayer::init() +{ + mPythia.readString("ProcessLevel:all off"); + + auto& param = EventDecayerParam::Instance(); + LOG(info) << "Init \'EventDecayer\' with following parameters"; + LOG(info) << param; + + // before the config files, so that they can override it + if (!setSeed(mPythia, param.seed)) { + return false; + } + + mConfigs.clear(); + auto& fallback = DecayerPythia8Param::Instance(); + for (int i = 0; i < 8; ++i) { + auto const& entry = param.config[i].empty() ? fallback.config[i] : param.config[i]; + if (entry.empty()) { + continue; + } + auto config = expand(entry); + LOG(info) << "Reading configuration from file: " << config; + if (!mPythia.readFile(config, true)) { + LOG(error) << "Failed to init \'EventDecayer\': problems with configuration file " << config; + return false; + } + mConfigs.push_back(config); + } + + for (auto name : kVertexLimits) { + if (mPythia.settings.flag(name)) { + LOG(error) << "EventDecayer does not support " << name << ", switching it off"; + mPythia.settings.flag(name, false); + } + } + if (param.tau0Max > 0.) { + mPythia.settings.flag("ParticleDecays:limitTau0", true); + mPythia.settings.parm("ParticleDecays:tau0Max", param.tau0Max); + } + + if (!skipDecays(mPythia.particleData, param.skipPdg)) { + return false; + } + + mVerbose = param.verbose; + mPythia.readString(param.showChanged ? "Init:showChangedParticleData on" : "Init:showChangedParticleData off"); + + if (!mPythia.init()) { + LOG(error) << "Failed to init \'EventDecayer\': init returned with error"; + return false; + } + mSeed = mPythia.settings.mode("Random:seed"); + LOG(info) << "EventDecayer: Pythia8 seed " << mSeed; + mTau0Max = mPythia.settings.flag("ParticleDecays:limitTau0") ? mPythia.settings.parm("ParticleDecays:tau0Max") : 0.; + if (mTau0Max > 0.) { + LOG(info) << "EventDecayer: species with tau0 above " << mTau0Max << " mm/c stay undecayed"; + } + + mG4Particles.clear(); + mHasG4Table = !param.g4ParticleTable.empty(); + if (mHasG4Table) { + auto path = expand(param.g4ParticleTable); + mG4Particles = loadG4ParticleTable(path); + if (mG4Particles.empty()) { + LOG(error) << "Failed to init \'EventDecayer\': no particles in Geant4 particle table " << path; + return false; + } + LOG(info) << "EventDecayer: read " << mG4Particles.size() << " particles from Geant4 particle table " << path; + } else { + LOG(info) << "EventDecayer: no Geant4 particle table given, skipping the Geant4 PDG and lifetime check"; + } + return true; +} + +/*****************************************************************/ + +std::string EventDecayer::configSummary() const +{ + std::string summary; + for (auto const& config : mConfigs) { + summary += (summary.empty() ? "" : ";") + config; + } + return summary; +} + +bool EventDecayer::vertexLimitsActive() const +{ + for (auto name : kVertexLimits) { + if (mPythia.settings.flag(name)) { + return true; + } + } + return false; +} + +/*****************************************************************/ + +bool EventDecayer::decayable(int pdg) const +{ + auto& particleData = mPythia.particleData; + return particleData.mayDecay(pdg) && particleData.canDecay(pdg) && !(mTau0Max > 0. && particleData.tau0(pdg) > mTau0Max); +} + +EventDecayerReport EventDecayer::decay(std::vector& particles) +{ + EventDecayerReport report; + const int nOriginal = particles.size(); + for (int i = 0; i < nOriginal; ++i) { + auto const& particle = particles[i]; + auto pdg = particle.GetPdgCode(); + if (!particle.TestBit(ParticleStatus::kToBeDone) || particle.GetUniqueID() != kPPrimary || !decayable(pdg)) { + continue; + } + if (decayOne(particles, i)) { + report.decayed[pdg]++; + } else { + report.failed[pdg]++; + } + } + + checkParticles(particles, report); + checkLifetimes(report); + return report; +} + +/*****************************************************************/ + +bool EventDecayer::decayOne(std::vector& particles, int index) +{ + // a copy, since appending to particles may reallocate + const TParticle mother = particles[index]; + auto& event = mPythia.event; + event.clear(); + // the caller's units to GeV/c and GeV + const double px = mother.Px() * mMomentumUnit; + const double py = mother.Py() * mMomentumUnit; + const double pz = mother.Pz() * mMomentumUnit; + const double e = mother.Energy() * mEnergyUnit; + const double m2 = e * e - px * px - py * py - pz * pz; + event.append(mother.GetPdgCode(), 1, 0, 0, px, py, pz, e, m2 > 0. ? std::sqrt(m2) : -std::sqrt(-m2)); + event[0].vProd(0., 0., 0., 0.); + event[0].tau(event[0].tau0() * mPythia.rndm.exp()); + if (!mPythia.moreDecays() || event[0].isFinal()) { + return false; + } + if (mVerbose) { + event.list(); + } + // partons (meMode 91-100) or a product that Pythia failed to decay make the whole decay fail + for (int k = 1; k < event.size(); ++k) { + auto const& entry = event[k]; + if (entry.isFinal() && (entry.isParton() || decayable(entry.id()))) { + return false; + } + } + + // Pythia entry 0 is the original particle, entry k > 0 goes to offset + k - 1 + const int offset = particles.size(); + auto toIndex = [index, offset](int k) { return k == 0 ? index : offset + k - 1; }; + // Pythia mm and mm/c to the caller's units + const double toPosition = 0.1 / mPositionUnit; + const double toTime = 1. / (kSpeedOfLightMMPerS * mTimeUnit); + + for (int k = 1; k < event.size(); ++k) { + auto const& entry = event[k]; + int mother1 = toIndex(entry.mother1()); + int mother2 = (entry.mother2() > 0 && entry.mother2() != entry.mother1()) ? toIndex(entry.mother2()) : -1; + int daughter1 = entry.daughter1() > 0 ? toIndex(entry.daughter1()) : -1; + int daughter2 = entry.daughter2() > 0 ? toIndex(entry.daughter2()) : daughter1; + TParticle product(entry.id(), 0, mother1, mother2, daughter1, daughter2, + entry.px() / mMomentumUnit, entry.py() / mMomentumUnit, entry.pz() / mMomentumUnit, entry.e() / mEnergyUnit, + mother.Vx() + entry.xProd() * toPosition, mother.Vy() + entry.yProd() * toPosition, + mother.Vz() + entry.zProd() * toPosition, mother.T() + entry.tProd() * toTime); + bool isFinal = entry.isFinal(); + o2::mcutils::MCGenHelper::encodeParticleStatusAndTracking(product, isFinal ? 1 : 2, std::abs(entry.status()), isFinal); + product.SetWeight(mother.GetWeight()); + product.SetUniqueID(kPPrimary); + particles.push_back(product); + } + + auto& decayed = particles[index]; + // the mother keeps its generator status code, only its HepMC status becomes decayed + o2::mcutils::MCGenHelper::encodeParticleStatusAndTracking(decayed, 2, o2::mcgenstatus::getGenStatusCode(decayed.GetStatusCode()), false); + decayed.SetFirstDaughter(toIndex(event[0].daughter1())); + decayed.SetLastDaughter(toIndex(event[0].daughter2() > 0 ? event[0].daughter2() : event[0].daughter1())); + return true; +} + +/*****************************************************************/ + +void EventDecayer::checkParticles(std::vector const& particles, EventDecayerReport& report) const +{ + auto pdgDB = o2::O2DatabasePDG::Instance(); + for (auto const& particle : particles) { + if (!particle.TestBit(ParticleStatus::kToBeDone)) { + continue; + } + auto pdg = particle.GetPdgCode(); + if (!pdgDB->GetParticle(pdg)) { + report.unknownInO2[pdg]++; + } + // ions (codes above 10^9) are created by Geant4 on demand, so they are not in the table + if (mHasG4Table && std::abs(pdg) < 1000000000 && mG4Particles.find(pdg) == mG4Particles.end()) { + report.unknownInG4[pdg]++; + } + } +} + +void EventDecayer::checkLifetimes(EventDecayerReport& report) const +{ + if (!mHasG4Table) { + return; + } + for (auto const& [pdg, count] : report.decayed) { + auto g4 = mG4Particles.find(pdg); + if (g4 == mG4Particles.end() || g4->second.ctauMM == 0.) { + continue; + } + double pythiaCtau = mPythia.particleData.tau0(pdg); + double g4Ctau = g4->second.ctauMM; + if (std::abs(pythiaCtau - g4Ctau) / g4Ctau > 0.01) { + report.lifetimeMismatch[pdg] = {pythiaCtau, g4Ctau}; + } + } +} + +} // namespace eventgen +} // namespace o2 diff --git a/Generators/src/EventDecayerParam.cxx b/Generators/src/EventDecayerParam.cxx new file mode 100644 index 0000000000000..cad1fefcff23d --- /dev/null +++ b/Generators/src/EventDecayerParam.cxx @@ -0,0 +1,13 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include "Generators/EventDecayerParam.h" +O2ParamImpl(o2::eventgen::EventDecayerParam); diff --git a/Generators/src/GeneratorsLinkDef.h b/Generators/src/GeneratorsLinkDef.h index 6269a311ee2e0..e8050312f51fd 100644 --- a/Generators/src/GeneratorsLinkDef.h +++ b/Generators/src/GeneratorsLinkDef.h @@ -59,6 +59,8 @@ #pragma link C++ class o2::eventgen::DecayerPythia8Param + ; #pragma link C++ class o2::conf::ConfigurableParamPromoter < o2::eventgen::GeneratorPythia8Param, o2::eventgen::Pythia8GenConfig> + ; #pragma link C++ class o2::conf::ConfigurableParamHelper < o2::eventgen::DecayerPythia8Param> + ; +#pragma link C++ class o2::eventgen::EventDecayerParam + ; +#pragma link C++ class o2::conf::ConfigurableParamHelper < o2::eventgen::EventDecayerParam> + ; #pragma link C++ class o2::eventgen::GeneratorFactory + ; #endif diff --git a/Generators/test/test_EventDecayer.cxx b/Generators/test/test_EventDecayer.cxx new file mode 100644 index 0000000000000..82121c4f31e16 --- /dev/null +++ b/Generators/test/test_EventDecayer.cxx @@ -0,0 +1,698 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#define BOOST_TEST_MODULE Test EventDecayer class +#define BOOST_TEST_MAIN +#define BOOST_TEST_DYN_LINK +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace o2::eventgen; + +namespace +{ + +// forces Xic+ -> Xi(1530)0 pi+ with the Geant4 lifetime +const std::string FORCE = "4232:mayDecay = on\n4232:onMode = off\n4232:oneChannel = 1 1.0 0 3324 211\n4232:tau0 = 0.1358\n"; + +// the default of EventDecayer.skipPdg +const std::string DEFAULTSKIP = "115 225 335 10221 10331 20223 20333 102132 102134 9010221"; + +// writes text to a fresh temporary file and returns its path +std::string writeCfg(std::string const& text) +{ + char name[] = "/tmp/test_EventDecayer_XXXXXX"; + int fd = mkstemp(name); + close(fd); + std::ofstream out(name); + out << text; + return name; +} + +// a primary generator particle moving along z +TParticle makeParticle(int pdg, double p, bool toBeDone, double vx = 0, double vy = 0, double vz = 0) +{ + bool ok = false; + double m = o2::O2DatabasePDG::Mass(pdg, ok); + double e = std::sqrt(p * p + m * m); + TParticle part(pdg, 1, -1, -1, -1, -1, 0., 0., p, e, vx, vy, vz, 0.); + o2::mcutils::MCGenHelper::encodeParticleStatusAndTracking(part, 1, 1, toBeDone); + part.SetUniqueID(kPPrimary); + return part; +} + +int hepmc(TParticle const& p) +{ + return o2::mcgenstatus::getHepMCStatusCode(p.GetStatusCode()); +} + +bool toBeDone(TParticle const& p) +{ + return p.TestBit(ParticleStatus::kToBeDone); +} + +// resets the parameter to the forcing configuration with nothing else set +void setParams(std::string const& g4table = "", std::string const& extraCfg = "", int seed = 0) +{ + static const std::string forcePath = writeCfg(FORCE); + o2::conf::ConfigurableParam::setValue("EventDecayer.config[0]", ""); + o2::conf::ConfigurableParam::setValue("EventDecayer.config[1]", forcePath); + o2::conf::ConfigurableParam::setValue("EventDecayer.config[2]", extraCfg); + o2::conf::ConfigurableParam::setValue("EventDecayer.g4ParticleTable", g4table); + o2::conf::ConfigurableParam::setValue("EventDecayer.seed", std::to_string(seed)); + o2::conf::ConfigurableParam::setValue("EventDecayer.tau0Max", "10."); + o2::conf::ConfigurableParam::setValue("EventDecayer.skipPdg", DEFAULTSKIP); +} + +// the output of decaying one Xic+ at (0.1,0.2,0.3) cm with p = 5 GeV +std::vector decayOneXic(EventDecayerReport& report) +{ + setParams(); + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + auto xic = makeParticle(4232, 5., true, 0.1, 0.2, 0.3); + // a distinctive generator status and a non-zero start time + o2::mcutils::MCGenHelper::encodeParticleStatusAndTracking(xic, 1, 85, true); + xic.SetProductionVertex(0.1, 0.2, 0.3, 1e-9); + std::vector particles{xic}; + report = decayer.decay(particles); + return particles; +} + +bool sameParticle(TParticle const& a, TParticle const& b) +{ + return a.GetPdgCode() == b.GetPdgCode() && a.Px() == b.Px() && a.Py() == b.Py() && a.Pz() == b.Pz() && + a.Energy() == b.Energy() && a.Vx() == b.Vx() && a.Vy() == b.Vy() && a.Vz() == b.Vz() && a.T() == b.T() && + a.GetStatusCode() == b.GetStatusCode() && a.TestBits(0xffffffff) == b.TestBits(0xffffffff) && + a.GetFirstMother() == b.GetFirstMother() && a.GetFirstDaughter() == b.GetFirstDaughter() && + a.GetLastDaughter() == b.GetLastDaughter() && a.GetUniqueID() == b.GetUniqueID(); +} + +// checks that a decay call left the particles as they were +void checkUnchanged(std::vector const& particles, std::vector const& before) +{ + BOOST_REQUIRE_EQUAL(particles.size(), before.size()); + for (size_t i = 0; i < before.size(); ++i) { + BOOST_CHECK(sameParticle(particles[i], before[i])); + } +} + +// a Xi- or Xi0, or their antiparticles +bool isXi(TParticle const& p) +{ + return std::abs(p.GetPdgCode()) == 3312 || std::abs(p.GetPdgCode()) == 3322; +} + +// 20000 Xic+ at p = 3 GeV decayed by one decayer, shared by the statistics cases +std::vector const& manyXic() +{ + static std::vector particles; + if (particles.empty()) { + setParams(); + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + for (int i = 0; i < 20000; ++i) { + particles.push_back(makeParticle(4232, 3., true)); + } + auto report = decayer.decay(particles); + BOOST_REQUIRE(report.decayed[4232] == 20000); + } + return particles; +} + +} // namespace + +BOOST_AUTO_TEST_CASE(xic_to_xistar_chain) +{ + EventDecayerReport report; + auto particles = decayOneXic(report); + + auto const& xic = particles[0]; + BOOST_CHECK(!toBeDone(xic)); + BOOST_CHECK_EQUAL(hepmc(xic), 2); + BOOST_CHECK_EQUAL(o2::mcgenstatus::getGenStatusCode(xic.GetStatusCode()), 85); + + int nXiStar = 0; + int iXiStar = -1; + int iPi = -1; + for (int i = 1; i < (int)particles.size(); ++i) { + if (particles[i].GetPdgCode() == 3324) { + ++nXiStar; + iXiStar = i; + } + if (particles[i].GetPdgCode() == 211 && particles[i].GetFirstMother() == 0) { + iPi = i; + } + } + BOOST_REQUIRE_EQUAL(nXiStar, 1); + BOOST_REQUIRE(iPi > 0); + + // the daughter range of every decayed particle holds exactly its direct products + for (int m = 0; m < (int)particles.size(); ++m) { + auto const& mother = particles[m]; + int nChildren = 0; + for (auto const& p : particles) { + nChildren += p.GetFirstMother() == m; + } + if (hepmc(mother) != 2) { + BOOST_CHECK_EQUAL(nChildren, 0); + continue; + } + BOOST_CHECK_EQUAL(mother.GetLastDaughter() - mother.GetFirstDaughter() + 1, nChildren); + for (int d = mother.GetFirstDaughter(); d <= mother.GetLastDaughter(); ++d) { + BOOST_CHECK_EQUAL(particles[d].GetFirstMother(), m); + } + } + BOOST_CHECK_EQUAL(xic.GetLastDaughter() - xic.GetFirstDaughter() + 1, 2); + + // products carry the Pythia decay status 91 and are flagged as primaries + for (int i = 1; i < (int)particles.size(); ++i) { + BOOST_CHECK_EQUAL(o2::mcgenstatus::getGenStatusCode(particles[i].GetStatusCode()), 91); + BOOST_CHECK_EQUAL(particles[i].GetUniqueID(), kPPrimary); + } + auto const& xiStar = particles[iXiStar]; + BOOST_CHECK_EQUAL(xiStar.GetFirstMother(), 0); + BOOST_CHECK(!toBeDone(xiStar)); + BOOST_CHECK_EQUAL(hepmc(xiStar), 2); + + std::vector dPdg; + for (int d = xiStar.GetFirstDaughter(); d >= 0 && d <= xiStar.GetLastDaughter(); ++d) { + dPdg.push_back(particles[d].GetPdgCode()); + BOOST_CHECK(toBeDone(particles[d])); + BOOST_CHECK_EQUAL(hepmc(particles[d]), 1); + BOOST_CHECK_EQUAL(particles[d].GetFirstMother(), iXiStar); + } + BOOST_REQUIRE_EQUAL(dPdg.size(), 2); + bool charged = (dPdg[0] == 3312 && dPdg[1] == 211) || (dPdg[0] == 211 && dPdg[1] == 3312); + bool neutral = (dPdg[0] == 3322 && dPdg[1] == 111) || (dPdg[0] == 111 && dPdg[1] == 3322); + BOOST_CHECK(charged || neutral); + + int first = xic.GetFirstDaughter(); + int last = xic.GetLastDaughter(); + BOOST_CHECK(first <= iXiStar && iXiStar <= last); + BOOST_CHECK(first <= iPi && iPi <= last); + + // all products start at the Xic decay point + double dx = xiStar.Vx() - 0.1; + double dy = xiStar.Vy() - 0.2; + double dz = xiStar.Vz() - 0.3; + BOOST_CHECK(std::sqrt(dx * dx + dy * dy + dz * dz) > 0.); + for (int i = 1; i < (int)particles.size(); ++i) { + BOOST_CHECK_CLOSE(particles[i].Vx(), xiStar.Vx(), 1e-6); + BOOST_CHECK_CLOSE(particles[i].Vy(), xiStar.Vy(), 1e-6); + BOOST_CHECK_CLOSE(particles[i].Vz(), xiStar.Vz(), 1e-6); + } + + // products start at the Xic time plus its flight time, |d| / (beta c) with c in cm/ns and T in s + double flightTime = std::sqrt(dx * dx + dy * dy + dz * dz) / (xic.P() / xic.Energy() * 29.9792458) * 1e-9; + for (int i = 1; i < (int)particles.size(); ++i) { + BOOST_CHECK_CLOSE(particles[i].T(), xic.T() + flightTime, 1e-6); + } + BOOST_CHECK_EQUAL(report.decayed[4232], 1); +} + +BOOST_AUTO_TEST_CASE(antiparticle_chain) +{ + setParams(); + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + std::vector particles{makeParticle(-4232, 5., true)}; + auto report = decayer.decay(particles); + int nAntiXiStar = 0; + for (auto const& p : particles) { + nAntiXiStar += p.GetPdgCode() == -3324; + } + BOOST_CHECK_EQUAL(nAntiXiStar, 1); + BOOST_CHECK_EQUAL(report.decayed[-4232], 1); +} + +BOOST_AUTO_TEST_CASE(proper_decay_length) +{ + auto const& particles = manyXic(); + double sum = 0.; + int n = 0; + for (int i = 0; i < 20000; ++i) { + auto const& mother = particles[i]; + auto const& daughter = particles[mother.GetFirstDaughter()]; + double dx = daughter.Vx() - mother.Vx(); + double dy = daughter.Vy() - mother.Vy(); + double dz = daughter.Vz() - mother.Vz(); + sum += std::sqrt(dx * dx + dy * dy + dz * dz) * mother.GetCalcMass() / mother.P(); + ++n; + } + BOOST_CHECK_CLOSE(sum / n, 0.01358, 2.); +} + +BOOST_AUTO_TEST_CASE(xistar_charge_split) +{ + auto const& particles = manyXic(); + int nXiStar = 0; + int nCharged = 0; + for (auto const& p : particles) { + if (p.GetPdgCode() != 3324) { + continue; + } + ++nXiStar; + for (int d = p.GetFirstDaughter(); d >= 0 && d <= p.GetLastDaughter(); ++d) { + if (particles[d].GetPdgCode() == 3312) { + ++nCharged; + } + } + } + BOOST_REQUIRE_EQUAL(nXiStar, 20000); + BOOST_CHECK_SMALL(double(nCharged) / nXiStar - 0.667, 0.02); +} + +BOOST_AUTO_TEST_CASE(strange_baryons_stay) +{ + setParams(); + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + std::vector particles; + for (int pdg : {3312, 3122, 310, 3334}) { + particles.push_back(makeParticle(pdg, 2., true)); + } + auto before = particles; + auto report = decayer.decay(particles); + checkUnchanged(particles, before); + BOOST_CHECK(report.decayed.empty()); +} + +BOOST_AUTO_TEST_CASE(nothing_to_decay) +{ + setParams(); + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + std::vector particles; + for (int pdg : {211, -321, 2212, 22}) { + particles.push_back(makeParticle(pdg, 1., true, 0.01, 0.02, 0.03)); + } + auto before = particles; + decayer.decay(particles); + checkUnchanged(particles, before); +} + +BOOST_AUTO_TEST_CASE(not_tobedone_untouched) +{ + setParams(); + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + std::vector particles{makeParticle(4232, 5., false)}; + auto before = particles; + auto report = decayer.decay(particles); + checkUnchanged(particles, before); + BOOST_CHECK(report.decayed.empty()); +} + +BOOST_AUTO_TEST_CASE(second_pass_is_noop) +{ + EventDecayerReport report; + auto particles = decayOneXic(report); + auto before = particles; + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + auto second = decayer.decay(particles); + checkUnchanged(particles, before); + BOOST_CHECK(second.decayed.empty()); +} + +BOOST_AUTO_TEST_CASE(weight_propagates) +{ + setParams(); + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + std::vector particles{makeParticle(4232, 5., true)}; + particles[0].SetWeight(2.5); + decayer.decay(particles); + BOOST_REQUIRE(particles.size() > 1); + for (size_t i = 1; i < particles.size(); ++i) { + BOOST_CHECK_EQUAL(particles[i].GetWeight(), 2.5); + } +} + +BOOST_AUTO_TEST_CASE(ions_and_unknown_to_pythia) +{ + setParams(); + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + std::vector particles{makeParticle(1000020040, 5., true), makeParticle(1010010030, 5., true)}; + auto before = particles; + auto report = decayer.decay(particles); + checkUnchanged(particles, before); + BOOST_CHECK(report.failed.find(1000020040) == report.failed.end()); + BOOST_CHECK(report.failed.find(1010010030) == report.failed.end()); +} + +BOOST_AUTO_TEST_CASE(impossible_decay) +{ + setParams(); + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + auto xic = makeParticle(4232, 5., true); + // a mass of 1 GeV, below the 3324 + 211 threshold + xic.SetMomentum(0., 0., 5., std::sqrt(25. + 1.0)); + xic.SetCalcMass(1.0); + std::vector particles{xic}; + auto before = particles; + auto report = decayer.decay(particles); + checkUnchanged(particles, before); + BOOST_CHECK_EQUAL(report.failed[4232], 1); + BOOST_CHECK(report.decayed.find(4232) == report.decayed.end()); +} + +BOOST_AUTO_TEST_CASE(vertex_limits_refused) +{ + for (auto limit : {"limitTau", "limitRadius", "limitCylinder"}) { + setParams("", writeCfg(std::string("ParticleDecays:") + limit + " = on\n")); + EventDecayer decayer; + BOOST_CHECK(decayer.init()); + BOOST_CHECK(!decayer.vertexLimitsActive()); + } + // the species-level cut of a config is kept when the tau0Max parameter is off + setParams("", writeCfg("ParticleDecays:limitTau0 = on\nParticleDecays:tau0Max = 0.05\n")); + o2::conf::ConfigurableParam::setValue("EventDecayer.tau0Max", "0"); + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + BOOST_CHECK(!decayer.vertexLimitsActive()); + std::vector particles{makeParticle(4232, 5., true)}; + auto report = decayer.decay(particles); + BOOST_CHECK_EQUAL(particles.size(), 1); + BOOST_CHECK(report.decayed.empty()); + BOOST_CHECK(report.failed.empty()); + setParams(); +} + +BOOST_AUTO_TEST_CASE(pdg_check) +{ + auto table = writeCfg("# test table\n111 pi0 0 0 0\n211 pi+ 7804.5 0 1\n3312 xi- 49.1 0 1\n3322 xi0 87.1 0 1\n"); + setParams(table); + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + std::vector particles{makeParticle(4232, 5., true)}; + auto report = decayer.decay(particles); + BOOST_CHECK_EQUAL(report.decayed[4232], 1); + BOOST_CHECK(report.unknownInG4.empty()); + + auto unknown = makeParticle(211, 1., true); + unknown.SetPdgCode(9900001); + std::vector more{unknown}; + report = decayer.decay(more); + BOOST_CHECK_EQUAL(report.unknownInO2[9900001], 1); + BOOST_CHECK_EQUAL(report.unknownInG4[9900001], 1); +} + +BOOST_AUTO_TEST_CASE(lifetime_mismatch) +{ + auto table = writeCfg("4232 xi_c+ 0.1300 0 1\n"); + setParams(table); + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + std::vector particles{makeParticle(4232, 5., true)}; + auto report = decayer.decay(particles); + BOOST_REQUIRE(report.lifetimeMismatch.count(4232) == 1); + BOOST_CHECK_CLOSE(report.lifetimeMismatch[4232].first, 0.1358, 1e-6); + BOOST_CHECK_CLOSE(report.lifetimeMismatch[4232].second, 0.1300, 1e-6); +} + +BOOST_AUTO_TEST_CASE(units_mm_and_mmc) +{ + // the same decay in cm and s, and in mm and mm/c (the Generator default units) + const double mmc = 3.3356410e-12; + std::vector inCm; + std::vector inMm; + for (bool mm : {false, true}) { + setParams("", "", 4711); + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + auto xic = makeParticle(4232, 5., true); + if (mm) { + decayer.setUnits(0.1, mmc); + xic.SetProductionVertex(1., 2., 3., 1e-9 / mmc); + } else { + xic.SetProductionVertex(0.1, 0.2, 0.3, 1e-9); + } + auto& particles = mm ? inMm : inCm; + particles.push_back(xic); + decayer.decay(particles); + } + setParams(); + BOOST_REQUIRE(inCm.size() > 1); + BOOST_REQUIRE_EQUAL(inCm.size(), inMm.size()); + for (size_t i = 1; i < inCm.size(); ++i) { + BOOST_CHECK_CLOSE(inMm[i].Vz(), 10. * inCm[i].Vz(), 1e-6); + BOOST_CHECK_CLOSE(inMm[i].T(), inCm[i].T() / mmc, 1e-6); + } + // in mm and mm/c the flight time is the decay length over beta + auto const& xic = inMm[0]; + auto const& d = inMm[xic.GetFirstDaughter()]; + double length = std::sqrt(std::pow(d.Vx() - 1., 2) + std::pow(d.Vy() - 2., 2) + std::pow(d.Vz() - 3., 2)); + BOOST_CHECK(length > 0.); + BOOST_CHECK_CLOSE(d.T() - xic.T(), length / (xic.P() / xic.Energy()), 1e-4); +} + +BOOST_AUTO_TEST_CASE(seeding) +{ + // decays one Xic with the given seed and returns the seed in use and the decay position + auto run = [](int seed) { + setParams("", "", seed); + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + std::vector particles{makeParticle(4232, 5., true)}; + decayer.decay(particles); + BOOST_REQUIRE(particles.size() > 1); + return std::make_pair(decayer.seed(), particles[1].Vz()); + }; + // a seed above the Pythia8 maximum is taken modulo 900000001 + auto large = run(1000000123); + auto reduced = run(100000122); + BOOST_CHECK_EQUAL(large.first, 100000122); + BOOST_CHECK_EQUAL(large.second, reduced.second); + // seed 0 derives the seed from gRandom, differently from GeneratorPythia8 which takes it unchanged + gRandom->SetSeed(777); + auto fromGRandom = run(0); + BOOST_CHECK(fromGRandom.first != 777); + BOOST_CHECK(fromGRandom.first >= 0 && fromGRandom.first <= 900000000); + auto explicitSeed = run(fromGRandom.first); + BOOST_CHECK_EQUAL(fromGRandom.second, explicitSeed.second); + gRandom->SetSeed(777); + BOOST_CHECK_EQUAL(run(0).first, fromGRandom.first); + gRandom->SetSeed(778); + BOOST_CHECK(run(0).first != fromGRandom.first); + BOOST_CHECK(run(1).second != run(2).second); + setParams(); +} + +BOOST_AUTO_TEST_CASE(partonic_decay_fails) +{ + // Xic+ -> s uu_1 with meMode 91, which Pythia8 leaves as partons + setParams("", writeCfg("4232:oneChannel = 1 1.0 91 3 2203\n")); + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + std::vector particles{makeParticle(4232, 5., true)}; + auto before = particles; + auto report = decayer.decay(particles); + checkUnchanged(particles, before); + BOOST_CHECK_EQUAL(report.failed[4232], 1); + BOOST_CHECK(report.decayed.empty()); +} + +BOOST_AUTO_TEST_CASE(failed_intermediate_decay_fails) +{ + // Xi(1530)0 -> Xi- p is closed, so Pythia8 cannot decay the intermediate + setParams("", writeCfg("3324:oneChannel = 1 1.0 0 3312 2212\n")); + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + std::vector particles{makeParticle(4232, 5., true)}; + auto before = particles; + auto report = decayer.decay(particles); + checkUnchanged(particles, before); + BOOST_CHECK_EQUAL(report.failed[4232], 1); + BOOST_CHECK(report.decayed.empty()); +} + +BOOST_AUTO_TEST_CASE(config_fallback_to_decayerpythia8) +{ + static const std::string forcePath = writeCfg(FORCE); + auto clear = []() { + for (int i = 0; i < 8; ++i) { + o2::conf::ConfigurableParam::setValue("EventDecayer.config[" + std::to_string(i) + "]", ""); + } + }; + // with every EventDecayer.config empty the DecayerPythia8 configs are read + setParams(); + clear(); + o2::conf::ConfigurableParam::setValue("DecayerPythia8.config[1]", forcePath); + { + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + auto summary = decayer.configSummary(); + BOOST_CHECK(summary.find("pythia8/decays/base.cfg;" + forcePath) != std::string::npos); + std::vector particles{makeParticle(4232, 5., true)}; + BOOST_CHECK_EQUAL(decayer.decay(particles).decayed[4232], 1); + } + // an explicitly set EventDecayer.config wins over the DecayerPythia8 one of the same slot + auto other = writeCfg("4232:mayDecay = off\n"); + o2::conf::ConfigurableParam::setValue("EventDecayer.config[1]", other); + { + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + auto summary = decayer.configSummary(); + BOOST_CHECK(summary.find("pythia8/decays/base.cfg;" + other) != std::string::npos); + BOOST_CHECK(summary.find(forcePath) == std::string::npos); + std::vector particles{makeParticle(4232, 5., true)}; + BOOST_CHECK(decayer.decay(particles).decayed.empty()); + } + o2::conf::ConfigurableParam::setValue("DecayerPythia8.config[1]", ""); + setParams(); +} + +BOOST_AUTO_TEST_CASE(tau0max_leaves_long_lived) +{ + // Xic+ (tau0 0.1358 mm) above tau0Max stays, and is not counted as failed + setParams(); + o2::conf::ConfigurableParam::setValue("EventDecayer.tau0Max", "0.1"); + { + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + std::vector particles{makeParticle(4232, 5., true)}; + auto before = particles; + auto report = decayer.decay(particles); + checkUnchanged(particles, before); + BOOST_CHECK(report.decayed.empty()); + BOOST_CHECK(report.failed.empty()); + } + // a decayable product above tau0Max (Xi-, Xi0) is left for transport + setParams("", writeCfg("3312:mayDecay = on\n3322:mayDecay = on\n")); + { + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + std::vector particles{makeParticle(4232, 5., true)}; + auto report = decayer.decay(particles); + BOOST_CHECK_EQUAL(report.decayed[4232], 1); + int nXi = 0; + for (auto const& p : particles) { + if (isXi(p)) { + ++nXi; + BOOST_CHECK(toBeDone(p)); + BOOST_CHECK_EQUAL(hepmc(p), 1); + } + } + BOOST_CHECK_EQUAL(nXi, 1); + } + // tau0Max = 0 switches the cut off, so the Xi decays as well + o2::conf::ConfigurableParam::setValue("EventDecayer.tau0Max", "0"); + { + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + std::vector particles{makeParticle(4232, 5., true)}; + decayer.decay(particles); + for (auto const& p : particles) { + if (isXi(p)) { + BOOST_CHECK_EQUAL(hepmc(p), 2); + } + } + } + setParams(); +} + +BOOST_AUTO_TEST_CASE(skip_pdg) +{ + // a listed code is never decayed, for the particle and its antiparticle + setParams(); + o2::conf::ConfigurableParam::setValue("EventDecayer.skipPdg", "421 4232"); + { + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + std::vector particles{makeParticle(4232, 5., true), makeParticle(-4232, 5., true)}; + auto before = particles; + auto report = decayer.decay(particles); + checkUnchanged(particles, before); + BOOST_CHECK(report.decayed.empty()); + BOOST_CHECK(report.failed.empty()); + } + // a listed product stays undecayed and is left for transport + o2::conf::ConfigurableParam::setValue("EventDecayer.skipPdg", "3324"); + { + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + std::vector particles{makeParticle(4232, 5., true)}; + auto report = decayer.decay(particles); + BOOST_CHECK_EQUAL(report.decayed[4232], 1); + BOOST_REQUIRE_EQUAL(particles.size(), 3); + for (size_t i = 1; i < particles.size(); ++i) { + BOOST_CHECK(toBeDone(particles[i])); + BOOST_CHECK_EQUAL(hepmc(particles[i]), 1); + } + } + // the default list keeps the Lambda(1520) for its Geant4 decay modes in O2 + setParams(); + { + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + std::vector particles{makeParticle(102134, 2., true), makeParticle(-102134, 2., true)}; + auto report = decayer.decay(particles); + BOOST_CHECK_EQUAL(particles.size(), 2); + BOOST_CHECK(report.decayed.empty()); + } + // with an empty list Pythia8 decays it + o2::conf::ConfigurableParam::setValue("EventDecayer.skipPdg", ""); + { + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + std::vector particles{makeParticle(102134, 2., true)}; + auto report = decayer.decay(particles); + BOOST_CHECK_EQUAL(report.decayed[102134], 1); + } + setParams(); +} + +BOOST_AUTO_TEST_CASE(units_momentum_and_energy) +{ + // the same decay in GeV and in MeV + std::vector inGeV; + std::vector inMeV; + for (bool mev : {false, true}) { + setParams("", "", 4711); + EventDecayer decayer; + BOOST_REQUIRE(decayer.init()); + auto xic = makeParticle(4232, 5., true); + if (mev) { + decayer.setUnits(1., 1., 1e-3, 1e-3); + xic.SetMomentum(1e3 * xic.Px(), 1e3 * xic.Py(), 1e3 * xic.Pz(), 1e3 * xic.Energy()); + } + auto& particles = mev ? inMeV : inGeV; + particles.push_back(xic); + decayer.decay(particles); + } + setParams(); + BOOST_REQUIRE(inGeV.size() > 1); + BOOST_REQUIRE_EQUAL(inGeV.size(), inMeV.size()); + for (size_t i = 1; i < inGeV.size(); ++i) { + BOOST_CHECK_EQUAL(inMeV[i].GetPdgCode(), inGeV[i].GetPdgCode()); + BOOST_CHECK_CLOSE(inMeV[i].Pz(), 1e3 * inGeV[i].Pz(), 1e-6); + BOOST_CHECK_CLOSE(inMeV[i].Energy(), 1e3 * inGeV[i].Energy(), 1e-6); + BOOST_CHECK_CLOSE(inMeV[i].Vz(), inGeV[i].Vz(), 1e-6); + } +} From 71da12c7fb20128cfc3ba41fa06a04fa3d52bec3 Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Thu, 22 Oct 2026 13:09:04 +0200 Subject: [PATCH 3/5] Decay generator particles before transport in Generator::ReadEvent This commit adds an optional decay step to Generator::ReadEvent, run after the trigger and before addTracks, and adds a unit test. - The new parameter EventDecayer.enable (default off) switches the step on. - The generator owns a lazily initialised EventDecayer and prints its report when destroyed. - The step uses the generator's position, time, momentum and energy units. - The step adds the header keys eventDecayer and eventDecayerSeed. - The step compiles out without Pythia8. Usage: o2-sim -g evtpool --configKeyValues "GeneratorEventPool.eventPoolPath=pool.root;EventDecayer.enable=true;EventDecayer.config[1]=my.cfg" https://its.cern.ch/jira/browse/O2-5903 Co-Authored-By: Claude Opus 5.5 --- Generators/CMakeLists.txt | 7 ++ .../include/Generators/EventDecayerParam.h | 1 + Generators/include/Generators/Generator.h | 7 ++ Generators/src/Generator.cxx | 58 ++++++++-- .../test/test_GeneratorEventDecayer.cxx | 101 ++++++++++++++++++ 5 files changed, 168 insertions(+), 6 deletions(-) create mode 100644 Generators/test/test_GeneratorEventDecayer.cxx diff --git a/Generators/CMakeLists.txt b/Generators/CMakeLists.txt index 52f4f5968698a..cf8afe3235712 100644 --- a/Generators/CMakeLists.txt +++ b/Generators/CMakeLists.txt @@ -157,6 +157,13 @@ if(doBuildSimulation) ENVIRONMENT O2_ROOT=${CMAKE_BINARY_DIR}/stage LABELS generator PUBLIC_LINK_LIBRARIES O2::Generators) + + o2_add_test(GeneratorEventDecayer NAME test_Generator_test_GeneratorEventDecayer + SOURCES test/test_GeneratorEventDecayer.cxx + COMPONENT_NAME Generator + ENVIRONMENT O2_ROOT=${CMAKE_BINARY_DIR}/stage + LABELS generator + PUBLIC_LINK_LIBRARIES O2::Generators) endif() # o2_add_test(GeneratorPythia8Param NAME test_Generator_test_GeneratorPythia8Param diff --git a/Generators/include/Generators/EventDecayerParam.h b/Generators/include/Generators/EventDecayerParam.h index 5962405b5cd91..832c445f357bf 100644 --- a/Generators/include/Generators/EventDecayerParam.h +++ b/Generators/include/Generators/EventDecayerParam.h @@ -23,6 +23,7 @@ namespace eventgen /// Settings of the EventDecayer, which decays generator particles with Pythia8 before transport struct EventDecayerParam : public o2::conf::ConfigurableParamHelper { + bool enable = false; // decay step in Generator::ReadEvent std::string config[8] = {"", "", "", "", "", "", "", ""}; // Pythia8 config files, an empty slot takes DecayerPythia8.config double tau0Max = 10.; // species with a larger tau0 in mm/c stay undecayed (<= 0 switches it off) std::string skipPdg = "115 225 335 10221 10331 20223 20333 102132 102134 9010221"; // PDG codes never decayed, matched on |pdg| diff --git a/Generators/include/Generators/Generator.h b/Generators/include/Generators/Generator.h index ab69ac5232f18..500a40d3579f7 100644 --- a/Generators/include/Generators/Generator.h +++ b/Generators/include/Generators/Generator.h @@ -35,6 +35,8 @@ namespace eventgen { class GenTPCLoopers; // Forward declaration +class EventDecayer; +struct EventDecayerReport; /*****************************************************************/ /*****************************************************************/ @@ -181,6 +183,11 @@ class Generator : public FairGenerator // Loopers generator instance o2::eventgen::GenTPCLoopers* mTPCLoopersGen = nullptr; + + // decay step after the trigger, created on first use when EventDecayer.enable is set + o2::eventgen::EventDecayer* mEventDecayer = nullptr; //! + o2::eventgen::EventDecayerReport* mEventDecayerReport = nullptr; //! + Bool_t decayParticles(o2::dataformats::MCEventHeader* header); #ifdef GENERATORS_WITH_TPCLOOPERS bool initTPCLoopersGen(); #endif diff --git a/Generators/src/Generator.cxx b/Generators/src/Generator.cxx index 1e2a1b7e36171..63b0fa36d5016 100644 --- a/Generators/src/Generator.cxx +++ b/Generators/src/Generator.cxx @@ -33,6 +33,10 @@ #include "Generators/TPCLoopers.h" #include "Generators/TPCLoopersParam.h" #endif +#ifdef GENERATORS_WITH_PYTHIA8 +#include "Generators/EventDecayer.h" +#include "Generators/EventDecayerParam.h" +#endif namespace o2 { @@ -117,6 +121,14 @@ Generator::~Generator() mTPCLoopersGen = nullptr; } #endif +#ifdef GENERATORS_WITH_PYTHIA8 + if (mEventDecayerReport) { + LOG(info) << "EventDecayer summary of generator " << GetName(); + mEventDecayerReport->print(); + } + delete mEventDecayer; + delete mEventDecayerReport; +#endif } /*****************************************************************/ @@ -335,6 +347,19 @@ Bool_t } } + auto header = primGen->GetEvent(); + auto o2header = dynamic_cast(header); + if (!header) { + LOG(fatal) << "MC event header is not a 'o2::dataformats::MCEventHeader' object"; + return kFALSE; + } + + /** decay particles before transport **/ + if (!decayParticles(o2header)) { + LOG(error) << "ReadEvent failed in decayParticles"; + return kFALSE; + } + /** add tracks **/ if (!addTracks(primGen)) { LOG(error) << "ReadEvent failed in addTracks"; @@ -342,12 +367,6 @@ Bool_t } /** update header **/ - auto header = primGen->GetEvent(); - auto o2header = dynamic_cast(header); - if (!header) { - LOG(fatal) << "MC event header is not a 'o2::dataformats::MCEventHeader' object"; - return kFALSE; - } updateHeader(o2header); updateSubGeneratorInformation(o2header); @@ -357,6 +376,33 @@ Bool_t /*****************************************************************/ +Bool_t + Generator::decayParticles(o2::dataformats::MCEventHeader* header) +{ +#ifdef GENERATORS_WITH_PYTHIA8 + if (!EventDecayerParam::Instance().enable) { + return kTRUE; + } + if (!mEventDecayer) { + mEventDecayer = new EventDecayer(); + if (!mEventDecayer->init()) { + LOG(fatal) << "Could not initialise the EventDecayer of generator " << GetName(); + return kFALSE; + } + mEventDecayerReport = new EventDecayerReport(); + } + mEventDecayer->setUnits(mPositionUnit, mTimeUnit, mMomentumUnit, mEnergyUnit); + mEventDecayerReport->merge(mEventDecayer->decay(mParticles)); + if (header) { + header->putInfo("eventDecayer", mEventDecayer->configSummary()); + header->putInfo("eventDecayerSeed", mEventDecayer->seed()); + } +#endif + return kTRUE; +} + +/*****************************************************************/ + Bool_t Generator::addTracks(FairPrimaryGenerator* primGen) { diff --git a/Generators/test/test_GeneratorEventDecayer.cxx b/Generators/test/test_GeneratorEventDecayer.cxx new file mode 100644 index 0000000000000..a91d82f09c6db --- /dev/null +++ b/Generators/test/test_GeneratorEventDecayer.cxx @@ -0,0 +1,101 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#define BOOST_TEST_MODULE Test decay step of Generator +#define BOOST_TEST_MAIN +#define BOOST_TEST_DYN_LINK +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace o2::eventgen; + +namespace +{ +// forces Xic+ -> Xi(1530)0 pi+ with the Geant4 lifetime +const std::string FORCE = "4232:mayDecay = on\n4232:onMode = off\n4232:oneChannel = 1 1.0 0 3324 211\n4232:tau0 = 0.1358\n"; + +std::string writeCfg(std::string const& text) +{ + char name[] = "/tmp/test_GeneratorEventDecayer_XXXXXX"; + int fd = mkstemp(name); + close(fd); + std::ofstream out(name); + out << text; + return name; +} + +// generates one boxgen event of five Xic+ with the decay step on or off +std::pair, o2::dataformats::MCEventHeader> generate(bool enable) +{ + static const std::string forcePath = writeCfg(FORCE); + o2::conf::ConfigurableParam::updateFromString("BoxGun.pdg=4232;BoxGun.number=5;EventDecayer.config[1]=" + forcePath + + ";EventDecayer.enable=" + (enable ? "true" : "false")); + GeneratorService service; + service.initService("boxgen", "", o2::eventgen::NoVertexOption()); + return service.generateEvent(); +} +} // namespace + +BOOST_AUTO_TEST_CASE(decay_step_on) +{ + auto [tracks, header] = generate(true); + int nXic = 0; + int nXiStar = 0; + for (int i = 0; i < (int)tracks.size(); ++i) { + auto const& t = tracks[i]; + if (t.GetPdgCode() == 4232) { + ++nXic; + BOOST_CHECK(!t.getToBeDone()); + BOOST_CHECK_EQUAL(o2::mcgenstatus::getHepMCStatusCode(t.getStatusCode()), 2); + } + if (t.GetPdgCode() == 3324) { + ++nXiStar; + int m = t.getMotherTrackId(); + BOOST_REQUIRE(m >= 0 && m < (int)tracks.size()); + BOOST_CHECK_EQUAL(tracks[m].GetPdgCode(), 4232); + BOOST_CHECK(!t.getToBeDone()); + // the Xic+ decay point in cm, displaced from the origin by about its decay length + double r = std::sqrt(t.Vx() * t.Vx() + t.Vy() * t.Vy() + t.Vz() * t.Vz()); + BOOST_CHECK(r > 0. && r < 10.); + for (int d = t.getFirstDaughterTrackId(); d >= 0 && d <= t.getLastDaughterTrackId(); ++d) { + BOOST_CHECK_EQUAL(tracks[d].getMotherTrackId(), i); + BOOST_CHECK(tracks[d].getToBeDone()); + } + } + } + BOOST_CHECK_EQUAL(nXic, 5); + BOOST_CHECK_EQUAL(nXiStar, 5); + bool valid = false; + auto summary = header.getInfo("eventDecayer", valid); + BOOST_CHECK(valid); + BOOST_CHECK(summary.find("pythia8/decays/base.cfg") != std::string::npos); +} + +BOOST_AUTO_TEST_CASE(decay_step_off) +{ + auto [tracks, header] = generate(false); + BOOST_REQUIRE_EQUAL(tracks.size(), 5); + for (auto const& t : tracks) { + BOOST_CHECK_EQUAL(t.GetPdgCode(), 4232); + BOOST_CHECK(t.getToBeDone()); + BOOST_CHECK_EQUAL(t.getFirstDaughterTrackId(), -1); + } + BOOST_CHECK(!header.hasInfo("eventDecayer")); +} From f2a6da169e30d2ebf9b7761c103b7efe1e137b6b Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Thu, 22 Oct 2026 13:09:04 +0200 Subject: [PATCH 4/5] Update charm and beauty lifetimes in the Pythia8 decay base config This commit updates the tau0 of eight weakly decaying charm and beauty hadrons in base.cfg to the PDG 2026 values. - Xi_c+ (4232): 0.132 -> 0.1358 mm. - Xi_c0 (4132): 0.0336 -> 0.04497 mm. - Omega_c0 (4332): 0.021 -> 0.08214 mm. - B_s0 (531): 0.439 -> 0.4542 mm, from the mean life 1/Gamma_s. - B_c+ (541): 0.138 -> 0.1529 mm. - Xi_b- (5132): 0.364 -> 0.4707 mm. - Xi_b0 (5232): 0.364 -> 0.4413 mm. - Omega_b- (5332): 0.364 -> 0.4917 mm. - Xi_c0 and Omega_c0 now differ from Geant4 by 1.3% and 2.2%, so the event decayer lists them in its lifetime report. - DecayerPythia8 reads the same file, but Geant4 takes only the momenta of the external decayer's products and places them at its own decay point (TG4ExtDecayer, G4Decay), so decays during transport do not change. https://pdglive.lbl.gov/Particle.action?init=0&node=S045 https://pdglive.lbl.gov/Particle.action?init=0&node=S048 https://pdglive.lbl.gov/Particle.action?init=0&node=S047 https://pdglive.lbl.gov/Particle.action?init=0&node=S086 https://pdglive.lbl.gov/Particle.action?init=0&node=S091 https://pdglive.lbl.gov/Particle.action?init=0&node=S069 https://pdglive.lbl.gov/Particle.action?init=0&node=S070 https://pdglive.lbl.gov/Particle.action?init=0&node=S063 https://its.cern.ch/jira/browse/O2-5903 Co-Authored-By: Claude Opus 5.5 --- Generators/share/pythia8/decays/base.cfg | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/Generators/share/pythia8/decays/base.cfg b/Generators/share/pythia8/decays/base.cfg index 5e19446b67316..ac8c9b7572d0e 100644 --- a/Generators/share/pythia8/decays/base.cfg +++ b/Generators/share/pythia8/decays/base.cfg @@ -32,3 +32,14 @@ Init:showChangedParticleData off 4232:mayDecay off # Xi_c+ 4332:mayDecay off # Omega_c0 5122:mayDecay off # Lambda_b0 + +### lifetimes (PDG 2026), Pythia defaults are older measurements + +531:tau0 = 0.4542 # B_s0 +541:tau0 = 0.1529 # B_c+ +4132:tau0 = 0.04497 # Xi_c0 +4232:tau0 = 0.1358 # Xi_c+ +4332:tau0 = 0.08214 # Omega_c0 +5132:tau0 = 0.4707 # Xi_b- +5232:tau0 = 0.4413 # Xi_b0 +5332:tau0 = 0.4917 # Omega_b- From fbd23f91b423f53eedff1703674bf097eb97b0ff Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Thu, 22 Oct 2026 14:29:33 +0200 Subject: [PATCH 5/5] Make the partonic-decay test independent of the Pythia8 version This rewrites the partonic decay test so that it passes with Pythia8 8.315 and 8.318. - Pythia8 8.318 hadronises the partons of a meMode 91 decay itself, so the decay succeeds, whereas 8.315 leaves partons and the decay fails. - The test now checks that no partonic particle is left to be transported, and the matching outcome in either case. https://its.cern.ch/jira/browse/O2-5903 Co-Authored-By: Claude Opus 5.5 --- Generators/test/test_EventDecayer.cxx | 26 +++++++++++++++++++++----- 1 file changed, 21 insertions(+), 5 deletions(-) diff --git a/Generators/test/test_EventDecayer.cxx b/Generators/test/test_EventDecayer.cxx index 82121c4f31e16..58f31ce50b422 100644 --- a/Generators/test/test_EventDecayer.cxx +++ b/Generators/test/test_EventDecayer.cxx @@ -504,18 +504,34 @@ BOOST_AUTO_TEST_CASE(seeding) setParams(); } -BOOST_AUTO_TEST_CASE(partonic_decay_fails) +// quarks, gluon, diquarks and the Pythia8 pseudo-particles 81-99 +bool isPartonic(int pdg) { - // Xic+ -> s uu_1 with meMode 91, which Pythia8 leaves as partons + int a = std::abs(pdg); + bool diquark = a > 1100 && a < 5600 && a / 1000 < 6 && (a / 10) % 10 == 0 && (a / 100) % 10 <= a / 1000 && (a / 100) % 10 > 0 && (a % 10) % 2 == 1; + return (a >= 1 && a <= 8) || a == 21 || (a >= 81 && a <= 99) || diquark; +} + +BOOST_AUTO_TEST_CASE(partonic_decay_never_transports_partons) +{ + // Xic+ -> s uu_1 with meMode 91; whether Pythia8 hadronises the partons depends on its version (8.318 does, 8.315 does not) setParams("", writeCfg("4232:oneChannel = 1 1.0 91 3 2203\n")); EventDecayer decayer; BOOST_REQUIRE(decayer.init()); std::vector particles{makeParticle(4232, 5., true)}; auto before = particles; auto report = decayer.decay(particles); - checkUnchanged(particles, before); - BOOST_CHECK_EQUAL(report.failed[4232], 1); - BOOST_CHECK(report.decayed.empty()); + for (auto const& p : particles) { + BOOST_CHECK(!(toBeDone(p) && isPartonic(p.GetPdgCode()))); + } + if (report.failed[4232] == 1) { + checkUnchanged(particles, before); + BOOST_CHECK(report.decayed.empty()); + } else { + BOOST_CHECK_EQUAL(report.decayed[4232], 1); + BOOST_CHECK(!toBeDone(particles[0])); + BOOST_CHECK(particles.size() > 1); + } } BOOST_AUTO_TEST_CASE(failed_intermediate_decay_fails)