diff --git a/Generators/CMakeLists.txt b/Generators/CMakeLists.txt index 094510c2c61a5..cf8afe3235712 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,22 @@ 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) + + 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 # 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..832c445f357bf --- /dev/null +++ b/Generators/include/Generators/EventDecayerParam.h @@ -0,0 +1,40 @@ +// 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 { + 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| + 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/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/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- 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/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/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..58f31ce50b422 --- /dev/null +++ b/Generators/test/test_EventDecayer.cxx @@ -0,0 +1,714 @@ +// 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(); +} + +// quarks, gluon, diquarks and the Pythia8 pseudo-particles 81-99 +bool isPartonic(int pdg) +{ + 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); + 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) +{ + // 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); + } +} 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")); +} 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; +}