From 3e5efc85f0cc297ec0334e055a276fa31b476414 Mon Sep 17 00:00:00 2001 From: Alberto Caliva Date: Sat, 31 Oct 2026 16:31:59 +0100 Subject: [PATCH] [PWGLF] Add hybrid PYTHIA8 generator for antideuteron coalescence in jets --- .../PWGLF/hybrid/hybrid_coalescence.json | 24 ++++++ .../GeneratorAntideuteronCoalescenceJet.ini | 6 ++ .../GeneratorAntideuteronCoalescenceJet.C | 85 +++++++++++++++++++ .../pythia8/generator_pythia8_coalescence.C | 70 +++++++++++---- 4 files changed, 169 insertions(+), 16 deletions(-) create mode 100644 MC/config/PWGLF/hybrid/hybrid_coalescence.json create mode 100644 MC/config/PWGLF/ini/GeneratorAntideuteronCoalescenceJet.ini create mode 100644 MC/config/PWGLF/ini/tests/GeneratorAntideuteronCoalescenceJet.C diff --git a/MC/config/PWGLF/hybrid/hybrid_coalescence.json b/MC/config/PWGLF/hybrid/hybrid_coalescence.json new file mode 100644 index 000000000..32c3f54ae --- /dev/null +++ b/MC/config/PWGLF/hybrid/hybrid_coalescence.json @@ -0,0 +1,24 @@ + +{ + "mode": "sequential", + "generators": [ + { + "name": "pythia8", + "config": { + "config": "${O2_ROOT}/share/Generators/egconfig/pythia8_inel.cfg", + "hooksFileName": "", + "hooksFuncName": "", + "includePartonEvent": false, + "particleFilter": "", + "verbose": 0 + } + }, + { + "name": "external", + "config": { + "iniFile": "${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/ini/GeneratorAntideuteronCoalescenceJet.ini" + } + } + ], + "fractions": [4, 1] +} diff --git a/MC/config/PWGLF/ini/GeneratorAntideuteronCoalescenceJet.ini b/MC/config/PWGLF/ini/GeneratorAntideuteronCoalescenceJet.ini new file mode 100644 index 000000000..a33512a27 --- /dev/null +++ b/MC/config/PWGLF/ini/GeneratorAntideuteronCoalescenceJet.ini @@ -0,0 +1,6 @@ +[GeneratorExternal] +fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/generator_pythia8_coalescence.C +funcName = generateCoalescence({-1000010020}, 1, 0.4) + +[GeneratorPythia8] +config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/pythia8/generator/pythia8_jet.cfg \ No newline at end of file diff --git a/MC/config/PWGLF/ini/tests/GeneratorAntideuteronCoalescenceJet.C b/MC/config/PWGLF/ini/tests/GeneratorAntideuteronCoalescenceJet.C new file mode 100644 index 000000000..c9f8171c0 --- /dev/null +++ b/MC/config/PWGLF/ini/tests/GeneratorAntideuteronCoalescenceJet.C @@ -0,0 +1,85 @@ + +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include "SimulationDataFormat/MCTrack.h" +#include "TFile.h" +#include "TTree.h" + +#include +#include +#include +#include +#endif + +int External() +{ + std::string path{"o2sim_Kine.root"}; + + std::vector injectedPDGs = { + -1000010020 // anti-deuteron + }; + + TFile file(path.c_str(), "READ"); + if (file.IsZombie()) { + std::cerr << "Cannot open ROOT file " << path << "\n"; + return 1; + } + + auto* tree = (TTree*)file.Get("o2sim"); + if (!tree) { + std::cerr << "Cannot find tree o2sim in file " << path << "\n"; + return 1; + } + + std::vector* tracks = nullptr; + tree->SetBranchAddress("MCTrack", &tracks); + + std::vector nSignal(injectedPDGs.size(), 0); + std::vector nEventsWithSignal(injectedPDGs.size(), 0); + + const auto nEvents = tree->GetEntries(); + + for (Long64_t i = 0; i < nEvents; ++i) { + tree->GetEntry(i); + + std::vector found(injectedPDGs.size(), false); + + for (const auto& track : *tracks) { + const int pdg = track.GetPdgCode(); + + auto it = std::find( + injectedPDGs.begin(), + injectedPDGs.end(), + pdg); + + if (it != injectedPDGs.end()) { + const auto index = std::distance(injectedPDGs.begin(), it); + ++nSignal[index]; + found[index] = true; + } + } + + for (size_t j = 0; j < injectedPDGs.size(); ++j) { + if (found[j]) { + ++nEventsWithSignal[j]; + } + } + } + + std::cout << "--------------------------------\n"; + std::cout << "# Events: " << nEvents << "\n"; + + for (size_t i = 0; i < injectedPDGs.size(); ++i) { + std::cout << "PDG: " << injectedPDGs[i] << "\n"; + std::cout << "# Generated particles: " << nSignal[i] << "\n"; + std::cout << "# Events containing particle: " + << nEventsWithSignal[i] << "\n"; + + if (nSignal[i] == 0) { + std::cerr << "No generated particles with PDG " + << injectedPDGs[i] << "\n"; + return 1; + } + } + + return 0; +} diff --git a/MC/config/PWGLF/pythia8/generator_pythia8_coalescence.C b/MC/config/PWGLF/pythia8/generator_pythia8_coalescence.C index 790710476..966376b53 100644 --- a/MC/config/PWGLF/pythia8/generator_pythia8_coalescence.C +++ b/MC/config/PWGLF/pythia8/generator_pythia8_coalescence.C @@ -9,6 +9,7 @@ #include "TDatabasePDG.h" #include "TSystem.h" #include "TMath.h" +#include #include #include #include @@ -18,20 +19,26 @@ using namespace Pythia8; R__ADD_INCLUDE_PATH($O2DPG_MC_CONFIG_ROOT) #include "MC/config/common/external/generator/CoalescencePythia8.h" -/// First version of the simple coalescence generator based PYTHIA8 +/// First version of the simple coalescence generator based on PYTHIA8 class GeneratorPythia8Coalescence : public o2::eventgen::GeneratorPythia8 { public: /// Constructor - GeneratorPythia8Coalescence(std::vector pdgList, int input_trigger_ratio = 1, double coal_momentum = 0.4) + GeneratorPythia8Coalescence(std::vector pdgList, int input_trigger_ratio = 1, double coal_momentum = 0.4) : o2::eventgen::GeneratorPythia8() { fmt::printf(">> Coalescence generator %d\n", input_trigger_ratio); mInverseTriggerRatio = input_trigger_ratio; mCoalMomentum = coal_momentum; mPdgList = pdgList; + + // CoalescencePythia8 accepts unsigned PDG codes + for (auto pdg : mPdgList) { + mAbsPdgList.push_back(static_cast(std::abs(pdg))); + } } + /// Destructor ~GeneratorPythia8Coalescence() = default; @@ -45,55 +52,86 @@ public: protected: bool generateEvent() override { - fmt::printf(">> Generating event %d\n", mGeneratedEvents); + fmt::printf(">> Generating event %llu\n", + static_cast(mGeneratedEvents)); + // Simple straightforward check to alternate generators if (mGeneratedEvents % mInverseTriggerRatio == 0) { - fmt::printf(">> Generating coalescence event %d\n", mGeneratedEvents); + fmt::printf(">> Generating coalescence event %llu\n", + static_cast(mGeneratedEvents)); + bool genOk = false; int localCounter{0}; + while (!genOk) { if (GeneratorPythia8::generateEvent()) { - genOk = CoalescencePythia8(mPythia.event, mPdgList, mCoalMomentum); + // Store the event size before coalescence + const int originalSize = mPythia.event.size(); + + // Run coalescence with absolute PDG codes + bool coalOk = CoalescencePythia8( + mPythia.event, mAbsPdgList, false, mCoalMomentum); + + if (coalOk && mPythia.event.size() > originalSize) { + for (int i = originalSize; i < mPythia.event.size(); ++i) { + const int producedPdg = mPythia.event[i].id(); + if (std::find(mPdgList.begin(), mPdgList.end(), producedPdg) != mPdgList.end()) { + genOk = true; + break; + } + } + } } + localCounter++; } - fmt::printf(">> Coalescence successful after %i generations\n", localCounter); + + fmt::printf(">> Coalescence successful after %i generations\n", + localCounter); std::cout << std::endl << std::endl; notifySubGenerator(1); } else { - fmt::printf(">> Generating minimum-bias event %d\n", mGeneratedEvents); + fmt::printf(">> Generating minimum-bias event %llu\n", + static_cast(mGeneratedEvents)); + // Generate minimum-bias event bool genOk = false; while (!genOk) { genOk = GeneratorPythia8::generateEvent(); } + notifySubGenerator(0); } + mGeneratedEvents++; return true; } - private: - std::vector mPdgList; /// list of pdg codes to be generated - float mCoalMomentum = 0.4; /// coalescence momentum - uint64_t mGeneratedEvents = 0; /// number of events generated so far - int mInverseTriggerRatio = 1; /// injection gap + std::vector mPdgList; /// Signed PDG codes + std::vector mAbsPdgList; /// Absolute PDG codes for coalescence + float mCoalMomentum = 0.4; /// Coalescence momentum + uint64_t mGeneratedEvents = 0; /// Number of generated events + int mInverseTriggerRatio = 1; /// Injection gap }; -///___________________________________________________________ -FairGenerator *generateCoalescence(std::vector pdgList, int input_trigger_ratio, double coal_momentum = 0.4) +/// Factory function +FairGenerator* generateCoalescence(std::vector pdgList, + int input_trigger_ratio, + double coal_momentum = 0.4) { - auto myGen = new GeneratorPythia8Coalescence(pdgList, input_trigger_ratio, coal_momentum); + auto myGen = new GeneratorPythia8Coalescence( + pdgList, input_trigger_ratio, coal_momentum); + auto seed = (gRandom->TRandom::GetSeed() % 900000000); myGen->readString("Random:setSeed on"); myGen->readString("Random:seed " + std::to_string(seed)); + return myGen; } -///___________________________________________________________ \ No newline at end of file