From d9d0affec7209b421b46dfaa57979865fa6fc233 Mon Sep 17 00:00:00 2001 From: Francesco Mazzaschi Date: Wed, 23 Apr 2025 10:20:10 +0200 Subject: [PATCH 1/6] Add header file for coalescence afterburner w/ Pythia --- Generators/CMakeLists.txt | 1 + .../Generators/NuclCoalescencePythia8.h | 151 ++++++++++++++++++ 2 files changed, 152 insertions(+) create mode 100644 Generators/include/Generators/NuclCoalescencePythia8.h diff --git a/Generators/CMakeLists.txt b/Generators/CMakeLists.txt index 02caa63df0d43..69f933aef0d3b 100644 --- a/Generators/CMakeLists.txt +++ b/Generators/CMakeLists.txt @@ -86,6 +86,7 @@ set(headers include/Generators/GeneratorService.h include/Generators/BoxGenerator.h include/Generators/FlowMapper.h + include/Generators/NuclCoalescencePythia8.h ) if(pythia_FOUND) diff --git a/Generators/include/Generators/NuclCoalescencePythia8.h b/Generators/include/Generators/NuclCoalescencePythia8.h new file mode 100644 index 0000000000000..c24eddd9c2769 --- /dev/null +++ b/Generators/include/Generators/NuclCoalescencePythia8.h @@ -0,0 +1,151 @@ +#ifndef O2_COALESCENCE_H +#define O2_COALESCENCE_H + +#include "Pythia8/Pythia.h" +#include "fairlogger/Logger.h" +#include "TParticlePDG.h" +#include "TDatabasePDG.h" +#include "TSystem.h" +#include "TMath.h" +#include +#include +#include +#include +using namespace Pythia8; + +namespace o2 +{ +/// Coalescence afterburner for Pythia8 +/// Utility to compute naive coalescence afterburner as done in PRL 126, 101101 (2021) + +enum NucleiBits { + kDeuteron = 0, + kTriton = 1, + kHe3 = 2, + kHyperTriton = 3, + kHe4 = 4, +}; + +std::vector pdgList = {10010010, 1000010030, 1000020030, 1010010030, 1000020040}; +std::vector massList = {1.875612, 2.80892113298, 2.808391, 2.991134, 3.727379}; + +bool coalPythia8(Pythia8::Event& event, int charge, int pdgCode, float mass, double coalescenceRadius, int iD1, int iD2, int iD3 = -1, int iD4 = -1) +{ + std::vector nucleonIDs = std::vector{iD1, iD2}; + // add A=3 and A=4 nuclei if enabled + if (iD3 > 0) { + nucleonIDs.push_back(iD3); + } + if (iD4 > 0) { + nucleonIDs.push_back(iD4); + } + // coalescence afterburner + Pythia8::Vec4 p; + for (auto nID : nucleonIDs) { + if (event[nID].status() < 0) { + return false; + } + p += event[nID].p(); + } + bool isCoalescence = true; + for (auto nID : nucleonIDs) { + auto pN = event[nID].p(); + pN.bstback(p); + if (pN.pAbs() > coalescenceRadius) { + isCoalescence = false; + break; + } + } + if (!isCoalescence) { + return false; + } + p.e(std::hypot(p.pAbs(), mass)); + /// In order to avoid the transport of the mother particles, but to still keep them in the stack, we set the status to negative and we mark the nucleus status as 94 (decay product) + event.append((charge * 2 - 1) * pdgCode, 94, 0, 0, 0, 0, 0, 0, p.px(), p.py(), p.pz(), p.e(), mass); + for (auto nID : nucleonIDs) { + event[nID].statusNeg(); + event[nID].daughter1(event.size() - 1); + } + fmt::printf(">> Adding a %i with p = %f, %f, %f, E = %f\n", (charge * 2 - 1) * pdgCode, p.px(), p.py(), p.pz(), p.e()); + return true; +} + +bool CoalAfterburner(Pythia8::Event& event, std::vector inputPdgList = {}, double coalMomentum = 0.4, int firstDauID = -1, int lastDauId = -1) +{ + const double coalescenceRadius{0.5 * 1.122462 * coalMomentum}; + + // if coalescence from a heavy hadron, loop only between firstDauID and lastDauID + int loopStart = firstDauID > 0 ? firstDauID : 0; + int loopEnd = lastDauId > 0 ? lastDauId : event.size(); + // fill the nuclear mask + uint8_t nuclearMask = 0; + for (auto nuclPdg : inputPdgList) { + if (nuclPdg == pdgList[NucleiBits::kDeuteron]) { + nuclearMask |= (1 << kDeuteron); + } else if (nuclPdg == pdgList[NucleiBits::kTriton]) { + nuclearMask |= (1 << kTriton); + } else if (nuclPdg == pdgList[NucleiBits::kHe3]) { + nuclearMask |= (1 << kHe3); + } else if (nuclPdg == pdgList[NucleiBits::kHyperTriton]) { + nuclearMask |= (1 << kHyperTriton); + } else if (nuclPdg == pdgList[NucleiBits::kHe4]) { + nuclearMask |= (1 << kHe4); + } else { + LOG(fatal) << "Unknown pdg code for coalescence generator: " << nuclPdg; + return false; + } + } + // fill nucleon pools + std::vector protons[2], neutrons[2], lambdas[2]; + for (auto iPart{loopStart}; iPart < loopEnd; ++iPart) { + if (std::abs(event[iPart].y()) > 1.) // skip particles with y > 1 + { + continue; + } + if (std::abs(event[iPart].id()) == 2212) { + protons[event[iPart].id() > 0].push_back(iPart); + } else if (std::abs(event[iPart].id()) == 2112) { + neutrons[event[iPart].id() > 0].push_back(iPart); + } else if (std::abs(event[iPart].id()) == 3122 && (nuclearMask & (1 << kHyperTriton))) { + lambdas[event[iPart].id() > 0].push_back(iPart); + } + } + // run coalescence + bool coalHappened = false; + for (int iC{0}; iC < 2; ++iC) { + for (int iP{0}; iP < protons[iC].size(); ++iP) { + for (int iN{0}; iN < neutrons[iC].size(); ++iN) { + if (nuclearMask & (1 << kDeuteron)) { + coalHappened |= coalPythia8(event, iC, pdgList[kDeuteron], massList[kDeuteron], coalescenceRadius, protons[iC][iP], neutrons[iC][iN]); + } + if (nuclearMask & (1 << kTriton)) { + for (int iN2{iN + 1}; iN2 < neutrons[iC].size(); ++iN2) { + coalHappened |= coalPythia8(event, iC, pdgList[kTriton], massList[kTriton], coalescenceRadius, protons[iC][iP], neutrons[iC][iN], neutrons[iC][iN2]); + } + } + if (nuclearMask & (1 << kHe3)) { + for (int iP2{iP + 1}; iP2 < protons[iC].size(); ++iP2) { + coalHappened |= coalPythia8(event, iC, pdgList[kHe3], massList[kHe3], coalescenceRadius, protons[iC][iP], protons[iC][iP2], neutrons[iC][iN]); + } + } + if (nuclearMask & (1 << kHyperTriton)) { + for (int iL{0}; iL < lambdas[iC].size(); ++iL) { + coalHappened |= coalPythia8(event, iC, pdgList[kHyperTriton], massList[kHyperTriton], coalescenceRadius, protons[iC][iP], neutrons[iC][iN], lambdas[iC][iL]); + } + } + if (nuclearMask & (1 << kHe4)) { + for (int iP2{iP + 1}; iP2 < protons[iC].size(); ++iP2) { + for (int iN2{iN + 1}; iN2 < neutrons[iC].size(); ++iN2) { + coalHappened |= coalPythia8(event, iC, pdgList[kHe4], massList[kHe4], coalescenceRadius, protons[iC][iP], protons[iC][iP2], neutrons[iC][iN], neutrons[iC][iN2]); + } + } + } + } + } + } + return coalHappened; +} + +} // namespace o2 + +#endif // O2_COALESCENCE_H \ No newline at end of file From 21ec6c2620ee91b213e03b9c15132ca524fa4ff7 Mon Sep 17 00:00:00 2001 From: Francesco Mazzaschi Date: Wed, 23 Apr 2025 10:50:52 +0200 Subject: [PATCH 2/6] Fix loop + copyright --- .../include/Generators/NuclCoalescencePythia8.h | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/Generators/include/Generators/NuclCoalescencePythia8.h b/Generators/include/Generators/NuclCoalescencePythia8.h index c24eddd9c2769..c53dd18e0b20d 100644 --- a/Generators/include/Generators/NuclCoalescencePythia8.h +++ b/Generators/include/Generators/NuclCoalescencePythia8.h @@ -1,3 +1,14 @@ +// 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 O2_COALESCENCE_H #define O2_COALESCENCE_H @@ -73,10 +84,9 @@ bool coalPythia8(Pythia8::Event& event, int charge, int pdgCode, float mass, dou bool CoalAfterburner(Pythia8::Event& event, std::vector inputPdgList = {}, double coalMomentum = 0.4, int firstDauID = -1, int lastDauId = -1) { const double coalescenceRadius{0.5 * 1.122462 * coalMomentum}; - // if coalescence from a heavy hadron, loop only between firstDauID and lastDauID int loopStart = firstDauID > 0 ? firstDauID : 0; - int loopEnd = lastDauId > 0 ? lastDauId : event.size(); + int loopEnd = lastDauId > 0 ? lastDauId : event.size() - 1; // fill the nuclear mask uint8_t nuclearMask = 0; for (auto nuclPdg : inputPdgList) { @@ -97,7 +107,7 @@ bool CoalAfterburner(Pythia8::Event& event, std::vector inputPdgLi } // fill nucleon pools std::vector protons[2], neutrons[2], lambdas[2]; - for (auto iPart{loopStart}; iPart < loopEnd; ++iPart) { + for (auto iPart{loopStart}; iPart <= loopEnd; ++iPart) { if (std::abs(event[iPart].y()) > 1.) // skip particles with y > 1 { continue; From cc9bc73d9f93d6d88f384d5bd9d5da26ec30e480 Mon Sep 17 00:00:00 2001 From: Francesco Mazzaschi Date: Thu, 24 Apr 2025 13:59:29 +0200 Subject: [PATCH 3/6] Add trivial coalescence + feed-down index reassignment --- .../Generators/NuclCoalescencePythia8.h | 64 ++++++++++++++----- 1 file changed, 47 insertions(+), 17 deletions(-) diff --git a/Generators/include/Generators/NuclCoalescencePythia8.h b/Generators/include/Generators/NuclCoalescencePythia8.h index c53dd18e0b20d..4bdc9c1976b9a 100644 --- a/Generators/include/Generators/NuclCoalescencePythia8.h +++ b/Generators/include/Generators/NuclCoalescencePythia8.h @@ -40,7 +40,7 @@ enum NucleiBits { std::vector pdgList = {10010010, 1000010030, 1000020030, 1010010030, 1000020040}; std::vector massList = {1.875612, 2.80892113298, 2.808391, 2.991134, 3.727379}; -bool coalPythia8(Pythia8::Event& event, int charge, int pdgCode, float mass, double coalescenceRadius, int iD1, int iD2, int iD3 = -1, int iD4 = -1) +bool coalPythia8(Pythia8::Event& event, int charge, int pdgCode, float mass, bool trivialCoal, double coalescenceRadius, bool nuclFromDecay, int iD1, int iD2, int iD3 = -1, int iD4 = -1) { std::vector nucleonIDs = std::vector{iD1, iD2}; // add A=3 and A=4 nuclei if enabled @@ -50,10 +50,10 @@ bool coalPythia8(Pythia8::Event& event, int charge, int pdgCode, float mass, dou if (iD4 > 0) { nucleonIDs.push_back(iD4); } - // coalescence afterburner Pythia8::Vec4 p; for (auto nID : nucleonIDs) { if (event[nID].status() < 0) { + // nucleon already used in coalescence return false; } p += event[nID].p(); @@ -62,7 +62,8 @@ bool coalPythia8(Pythia8::Event& event, int charge, int pdgCode, float mass, dou for (auto nID : nucleonIDs) { auto pN = event[nID].p(); pN.bstback(p); - if (pN.pAbs() > coalescenceRadius) { + // trivial coal does not check the distance of the nucleons + if (pN.pAbs() > coalescenceRadius && !trivialCoal) { isCoalescence = false; break; } @@ -71,22 +72,49 @@ bool coalPythia8(Pythia8::Event& event, int charge, int pdgCode, float mass, dou return false; } p.e(std::hypot(p.pAbs(), mass)); - /// In order to avoid the transport of the mother particles, but to still keep them in the stack, we set the status to negative and we mark the nucleus status as 94 (decay product) - event.append((charge * 2 - 1) * pdgCode, 94, 0, 0, 0, 0, 0, 0, p.px(), p.py(), p.pz(), p.e(), mass); - for (auto nID : nucleonIDs) { - event[nID].statusNeg(); - event[nID].daughter1(event.size() - 1); + + if (!nuclFromDecay) { + /// keep the original nucleons with negative status, store the mother with status 94 + event.append((charge * 2 - 1) * pdgCode, 94, 0, 0, 0, 0, 0, 0, p.px(), p.py(), p.pz(), p.e(), mass); + for (auto nID : nucleonIDs) { + event[nID].statusNeg(); + event[nID].daughter1(event.size() - 1); + } + } else { + // first nucleon will be replaced by the nucleus, the others will be removed + bool swap = true; + int nRemoved = 0; + for (auto iPart{0}; iPart < event.size(); ++iPart) { + for (auto nID : nucleonIDs) { + if (iPart == nID && swap) { + // replace the nucleon with the nucleus + LOG(debug) << "Replacing nucleon with index " << iPart << " and pdg code " << event[iPart].id() << " with nucleus with pdg code " << (charge * 2 - 1) * pdgCode; + event[iPart].id((charge * 2 - 1) * pdgCode); + event[iPart].status(94); + event[iPart].px(p.px()); + event[iPart].py(p.py()); + event[iPart].pz(p.pz()); + event[iPart].e(std::hypot(p.pAbs(), mass)); + event[iPart].m(mass); + swap = false; + } else if (iPart == nID - nRemoved && !swap) { + LOG(debug) << "Removing nucleon with index " << iPart << " and pdg code " << event[iPart].id(); + event.remove(iPart, iPart, true); + nRemoved++; + } + } + } } - fmt::printf(">> Adding a %i with p = %f, %f, %f, E = %f\n", (charge * 2 - 1) * pdgCode, p.px(), p.py(), p.pz(), p.e()); + LOG(debug) "Adding a " << (charge * 2 - 1) * pdgCode << " with p = " << p.px() << ", " << p.py() << ", " << p.pz() << ", E = " << p.e(); return true; } -bool CoalAfterburner(Pythia8::Event& event, std::vector inputPdgList = {}, double coalMomentum = 0.4, int firstDauID = -1, int lastDauId = -1) +bool CoalAfterburner(Pythia8::Event& event, std::vector inputPdgList = {}, bool trivialCoal = false, double coalMomentum = 0.4, int firstDauID = -1, int lastDauId = -1) { const double coalescenceRadius{0.5 * 1.122462 * coalMomentum}; // if coalescence from a heavy hadron, loop only between firstDauID and lastDauID - int loopStart = firstDauID > 0 ? firstDauID : 0; - int loopEnd = lastDauId > 0 ? lastDauId : event.size() - 1; + int loopStart = firstDauID > -1 ? firstDauID : 0; + int loopEnd = lastDauId > -1 ? lastDauId : event.size() - 1; // fill the nuclear mask uint8_t nuclearMask = 0; for (auto nuclPdg : inputPdgList) { @@ -121,32 +149,34 @@ bool CoalAfterburner(Pythia8::Event& event, std::vector inputPdgLi } } // run coalescence + bool nuclFromDecay = firstDauID > -1; bool coalHappened = false; + for (int iC{0}; iC < 2; ++iC) { for (int iP{0}; iP < protons[iC].size(); ++iP) { for (int iN{0}; iN < neutrons[iC].size(); ++iN) { if (nuclearMask & (1 << kDeuteron)) { - coalHappened |= coalPythia8(event, iC, pdgList[kDeuteron], massList[kDeuteron], coalescenceRadius, protons[iC][iP], neutrons[iC][iN]); + coalHappened |= coalPythia8(event, iC, pdgList[kDeuteron], massList[kDeuteron], trivialCoal, coalescenceRadius, nuclFromDecay, protons[iC][iP], neutrons[iC][iN]); } if (nuclearMask & (1 << kTriton)) { for (int iN2{iN + 1}; iN2 < neutrons[iC].size(); ++iN2) { - coalHappened |= coalPythia8(event, iC, pdgList[kTriton], massList[kTriton], coalescenceRadius, protons[iC][iP], neutrons[iC][iN], neutrons[iC][iN2]); + coalHappened |= coalPythia8(event, iC, pdgList[kTriton], massList[kTriton], trivialCoal, coalescenceRadius, nuclFromDecay, protons[iC][iP], neutrons[iC][iN], neutrons[iC][iN2]); } } if (nuclearMask & (1 << kHe3)) { for (int iP2{iP + 1}; iP2 < protons[iC].size(); ++iP2) { - coalHappened |= coalPythia8(event, iC, pdgList[kHe3], massList[kHe3], coalescenceRadius, protons[iC][iP], protons[iC][iP2], neutrons[iC][iN]); + coalHappened |= coalPythia8(event, iC, pdgList[kHe3], massList[kHe3], trivialCoal, coalescenceRadius, nuclFromDecay, protons[iC][iP], protons[iC][iP2], neutrons[iC][iN]); } } if (nuclearMask & (1 << kHyperTriton)) { for (int iL{0}; iL < lambdas[iC].size(); ++iL) { - coalHappened |= coalPythia8(event, iC, pdgList[kHyperTriton], massList[kHyperTriton], coalescenceRadius, protons[iC][iP], neutrons[iC][iN], lambdas[iC][iL]); + coalHappened |= coalPythia8(event, iC, pdgList[kHyperTriton], massList[kHyperTriton], trivialCoal, coalescenceRadius, nuclFromDecay, protons[iC][iP], neutrons[iC][iN], lambdas[iC][iL]); } } if (nuclearMask & (1 << kHe4)) { for (int iP2{iP + 1}; iP2 < protons[iC].size(); ++iP2) { for (int iN2{iN + 1}; iN2 < neutrons[iC].size(); ++iN2) { - coalHappened |= coalPythia8(event, iC, pdgList[kHe4], massList[kHe4], coalescenceRadius, protons[iC][iP], protons[iC][iP2], neutrons[iC][iN], neutrons[iC][iN2]); + coalHappened |= coalPythia8(event, iC, pdgList[kHe4], massList[kHe4], trivialCoal, coalescenceRadius, nuclFromDecay, protons[iC][iP], protons[iC][iP2], neutrons[iC][iN], neutrons[iC][iN2]); } } } From 9cf8f6aafc7c61468e74e1d86b383f1ea1fe4f5a Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Thu, 24 Apr 2025 12:00:38 +0000 Subject: [PATCH 4/6] Please consider the following formatting changes --- Generators/include/Generators/NuclCoalescencePythia8.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Generators/include/Generators/NuclCoalescencePythia8.h b/Generators/include/Generators/NuclCoalescencePythia8.h index 4bdc9c1976b9a..8734122a52e77 100644 --- a/Generators/include/Generators/NuclCoalescencePythia8.h +++ b/Generators/include/Generators/NuclCoalescencePythia8.h @@ -105,7 +105,8 @@ bool coalPythia8(Pythia8::Event& event, int charge, int pdgCode, float mass, boo } } } - LOG(debug) "Adding a " << (charge * 2 - 1) * pdgCode << " with p = " << p.px() << ", " << p.py() << ", " << p.pz() << ", E = " << p.e(); + LOG(debug) + "Adding a " << (charge * 2 - 1) * pdgCode << " with p = " << p.px() << ", " << p.py() << ", " << p.pz() << ", E = " << p.e(); return true; } From 7d73f8725b1849a5bf86149a99313f5e82f97f8a Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Thu, 24 Apr 2025 12:12:26 +0000 Subject: [PATCH 5/6] Please consider the following formatting changes --- Generators/include/Generators/NuclCoalescencePythia8.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Generators/include/Generators/NuclCoalescencePythia8.h b/Generators/include/Generators/NuclCoalescencePythia8.h index 8734122a52e77..d7498a52756e9 100644 --- a/Generators/include/Generators/NuclCoalescencePythia8.h +++ b/Generators/include/Generators/NuclCoalescencePythia8.h @@ -106,7 +106,7 @@ bool coalPythia8(Pythia8::Event& event, int charge, int pdgCode, float mass, boo } } LOG(debug) - "Adding a " << (charge * 2 - 1) * pdgCode << " with p = " << p.px() << ", " << p.py() << ", " << p.pz() << ", E = " << p.e(); + "Adding a " << (charge * 2 - 1) * pdgCode << " with p = " << p.px() << ", " << p.py() << ", " << p.pz() << ", E = " << p.e(); return true; } From b2e89b03c8dd624cc8340968393eabb7ef90663b Mon Sep 17 00:00:00 2001 From: Francesco Mazzaschi Date: Thu, 24 Apr 2025 16:37:13 +0200 Subject: [PATCH 6/6] Fix log printout --- Generators/include/Generators/NuclCoalescencePythia8.h | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Generators/include/Generators/NuclCoalescencePythia8.h b/Generators/include/Generators/NuclCoalescencePythia8.h index d7498a52756e9..dfd63624eb759 100644 --- a/Generators/include/Generators/NuclCoalescencePythia8.h +++ b/Generators/include/Generators/NuclCoalescencePythia8.h @@ -105,8 +105,7 @@ bool coalPythia8(Pythia8::Event& event, int charge, int pdgCode, float mass, boo } } } - LOG(debug) - "Adding a " << (charge * 2 - 1) * pdgCode << " with p = " << p.px() << ", " << p.py() << ", " << p.pz() << ", E = " << p.e(); + LOG(debug) << "Adding a " << (charge * 2 - 1) * pdgCode << " with p = " << p.px() << ", " << p.py() << ", " << p.pz() << ", E = " << p.e(); return true; }