From 783c4bfc86caecf1ca1a1d189c72c612aa9124d2 Mon Sep 17 00:00:00 2001 From: Alberto Caliva Date: Wed, 3 Jun 2026 11:52:46 +0200 Subject: [PATCH 1/2] [PWGLF] Refactoring and bug fixes --- .../Nuspex/coalescenceTreeProducer.cxx | 716 +++++------------- 1 file changed, 208 insertions(+), 508 deletions(-) diff --git a/PWGLF/TableProducer/Nuspex/coalescenceTreeProducer.cxx b/PWGLF/TableProducer/Nuspex/coalescenceTreeProducer.cxx index 4b478ad4d89..1dbc95139c6 100644 --- a/PWGLF/TableProducer/Nuspex/coalescenceTreeProducer.cxx +++ b/PWGLF/TableProducer/Nuspex/coalescenceTreeProducer.cxx @@ -10,8 +10,7 @@ // or submit itself to any jurisdiction. // /// \file coalescenceTreeProducer.cxx -/// \brief Tree producer for deuteron, triton, helium-3, and hypertriton -/// coalescence studies. +/// \brief Tree producer for deuteron, triton, helium-3, and hypertriton coalescence studies. /// /// The task loops over generated MC particles and builds bound-state /// candidates according to the selected species. For each candidate, @@ -55,125 +54,75 @@ #include #include +using namespace std; using namespace o2; +using namespace o2::soa; +using namespace o2::aod; using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::constants::physics; +using namespace o2::constants::math; + +// Lightweight particle container +struct ReducedParticle { + int64_t idPart = 0; + int pdgPart = 0; + float xPart = 0.f; + float yPart = 0.f; + float zPart = 0.f; + float tPart = 0.f; + float pxPart = 0.f; + float pyPart = 0.f; + float pzPart = 0.f; +}; struct CoalescenceTreeProducer { - // Supported bound-state species - enum BoundStateSpecies { - kDeuteron = 0, - kTriton = 1, - kHelium3 = 2, - kHypertriton = 3 - }; - + // Histogram registry HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; - // Analysis parameters - Configurable zVtx{"zVtx", 10.f, "Maximum |z vertex| in cm"}; + // Configurable analysis parameters + Configurable zVtxMax{"zVtxMax", 10.f, "Maximum |z vertex| in cm"}; Configurable etaMax{"etaMax", 1.5f, "Maximum |eta| for generated particles"}; Configurable pRhoMax{"pRhoMax", 1.0f, "Maximum Jacobi p_rho in GeV/c"}; Configurable pLambdaMax{"pLambdaMax", 1.0f, "Maximum Jacobi p_lambda in GeV/c"}; - Configurable boundStateSpecies{"boundStateSpecies", kHypertriton, "Species selector: 0 = deuteron, 1 = triton, 2 = helium3, 3 = hypertriton"}; - - Preslice mcParticlesPerMcCollision = aod::mcparticle::mcCollisionId; + // Output tree storing bound-state candidates OutputObj treeBoundState{"treeBoundState", OutputObjHandlingPolicy::AnalysisObject}; - int64_t eventID = 0; // Event ID - int64_t idB1 = 0, idB2 = 0, idB3 = 0; // MC particle IDs + // Variables to be stored in the tree + Long64_t eventID = 0; // Event identifier + Long64_t nMatterCandidatesPerEvent = 0; // Number of matter candidates per event + Long64_t nAntimatterCandidatesPerEvent = 0; // Number of antimatter candidates per event + Long64_t idB1 = 0, idB2 = 0, idB3 = 0; // Constituent identifiers + int pdgB1 = 0, pdgB2 = 0, pdgB3 = 0; // Constituent PDG codes - int pdgB1 = 0, pdgB2 = 0, pdgB3 = 0; - int chargeB1 = 0, chargeB2 = 0, chargeB3 = 0; - - // Space-time coordinates and momentum components + // Constituent space-time coordinates and momentum components float xB1 = 0.f, yB1 = 0.f, zB1 = 0.f, tB1 = 0.f, pxB1 = 0.f, pyB1 = 0.f, pzB1 = 0.f; float xB2 = 0.f, yB2 = 0.f, zB2 = 0.f, tB2 = 0.f, pxB2 = 0.f, pyB2 = 0.f, pzB2 = 0.f; float xB3 = 0.f, yB3 = 0.f, zB3 = 0.f, tB3 = 0.f, pxB3 = 0.f, pyB3 = 0.f, pzB3 = 0.f; - struct Particle { - int64_t id = 0; - int pdg = 0; - int charge = 0; - - float x = 0.f; - float y = 0.f; - float z = 0.f; - float t = 0.f; - - float px = 0.f; - float py = 0.f; - float pz = 0.f; - - float mass = 0.f; - - Particle() = default; - - template - explicit Particle(T const& p) - : id(p.globalIndex()), - pdg(p.pdgCode()), - charge(chargeFromPdg(pdg)), - x(p.vx()), - y(p.vy()), - z(p.vz()), - t(p.vt()), - px(p.px()), - py(p.py()), - pz(p.pz()), - mass(static_cast(massFromPdg(pdg))) - { - } - - static int chargeFromPdg(int pdg) - { - if (pdg == PDG_t::kProton) { - return 1; - } - if (pdg == PDG_t::kProtonBar) { - return -1; - } - return 0; - } - - static double massFromPdg(int pdg) - { - switch (std::abs(pdg)) { - case PDG_t::kProton: - return o2::constants::physics::MassProton; - case PDG_t::kNeutron: - return o2::constants::physics::MassNeutron; - case PDG_t::kLambda0: - return o2::constants::physics::MassLambda0; - default: - return -1.; - } - } - }; - - void init(InitContext&) + // Initialize output objects and histograms + void init(InitContext const&) { + // Event selection counter registry.add("eventCounter", "event Counter", HistType::kTH1F, {{5, 0, 5, "counter"}}); registry.get(HIST("eventCounter"))->GetXaxis()->SetBinLabel(1, "Before z-vertex cut"); registry.get(HIST("eventCounter"))->GetXaxis()->SetBinLabel(2, "After z-vertex cut"); - registry.get(HIST("eventCounter"))->GetXaxis()->SetBinLabel(3, "After non-empty lists"); + registry.get(HIST("eventCounter"))->GetXaxis()->SetBinLabel(3, "After non-empty constituent lists"); + registry.get(HIST("eventCounter"))->GetXaxis()->SetBinLabel(4, "After non-empty candidate lists"); + + // Output tree for bound-state candidates + treeBoundState.setObject(new TTree("BoundStateTree", "Tree for coalescence")); - // Tree for pairs/triplets. - // For deuterons only the first two baryons are stored. - // For three-body states, the third-baryon branches are also created. - treeBoundState.setObject(new TTree("BoundStateTree", "Tree for coalescence studies")); + // Event information treeBoundState->Branch("eventID", &eventID, "eventID/L"); + treeBoundState->Branch("nMatterCandidatesPerEvent", &nMatterCandidatesPerEvent, "nMatterCandidatesPerEvent/L"); + treeBoundState->Branch("nAntimatterCandidatesPerEvent", &nAntimatterCandidatesPerEvent, "nAntimatterCandidatesPerEvent/L"); + // Constituent 1: ID, PDG code, space-time coordinates, and momentum components treeBoundState->Branch("idB1", &idB1, "idB1/L"); - treeBoundState->Branch("idB2", &idB2, "idB2/L"); - treeBoundState->Branch("pdgB1", &pdgB1, "pdgB1/I"); - treeBoundState->Branch("pdgB2", &pdgB2, "pdgB2/I"); - - treeBoundState->Branch("chargeB1", &chargeB1, "chargeB1/I"); - treeBoundState->Branch("chargeB2", &chargeB2, "chargeB2/I"); - treeBoundState->Branch("xB1", &xB1, "xB1/F"); treeBoundState->Branch("yB1", &yB1, "yB1/F"); treeBoundState->Branch("zB1", &zB1, "zB1/F"); @@ -182,6 +131,9 @@ struct CoalescenceTreeProducer { treeBoundState->Branch("pyB1", &pyB1, "pyB1/F"); treeBoundState->Branch("pzB1", &pzB1, "pzB1/F"); + // Constituent 2: ID, PDG code, space-time coordinates, and momentum components + treeBoundState->Branch("idB2", &idB2, "idB2/L"); + treeBoundState->Branch("pdgB2", &pdgB2, "pdgB2/I"); treeBoundState->Branch("xB2", &xB2, "xB2/F"); treeBoundState->Branch("yB2", &yB2, "yB2/F"); treeBoundState->Branch("zB2", &zB2, "zB2/F"); @@ -190,12 +142,10 @@ struct CoalescenceTreeProducer { treeBoundState->Branch("pyB2", &pyB2, "pyB2/F"); treeBoundState->Branch("pzB2", &pzB2, "pzB2/F"); - if (static_cast(boundStateSpecies) != kDeuteron) { + // Constituent 3: ID, PDG code, space-time coordinates, and momentum components (three-body bound states only) + if (doprocessHypertriton) { treeBoundState->Branch("idB3", &idB3, "idB3/L"); - treeBoundState->Branch("pdgB3", &pdgB3, "pdgB3/I"); - treeBoundState->Branch("chargeB3", &chargeB3, "chargeB3/I"); - treeBoundState->Branch("xB3", &xB3, "xB3/F"); treeBoundState->Branch("yB3", &yB3, "yB3/F"); treeBoundState->Branch("zB3", &zB3, "zB3/F"); @@ -206,477 +156,227 @@ struct CoalescenceTreeProducer { } } - ROOT::Math::PxPyPzMVector makeFourVector(Particle const& p) - { - return ROOT::Math::PxPyPzMVector{p.px, p.py, p.pz, p.mass}; - } - - void storePair(Particle const& b1, Particle const& b2) - { - idB1 = b1.id; - idB2 = b2.id; - - pdgB1 = b1.pdg; - pdgB2 = b2.pdg; - - chargeB1 = b1.charge; - chargeB2 = b2.charge; - - xB1 = b1.x; - yB1 = b1.y; - zB1 = b1.z; - tB1 = b1.t; - pxB1 = b1.px; - pyB1 = b1.py; - pzB1 = b1.pz; - - xB2 = b2.x; - yB2 = b2.y; - zB2 = b2.z; - tB2 = b2.t; - pxB2 = b2.px; - pyB2 = b2.py; - pzB2 = b2.pz; - - treeBoundState->Fill(); - } - - void storeTriplet(Particle const& b1, Particle const& b2, Particle const& b3) + // Assign mass based on PDG code + static double massFromPdg(int pdg) { - idB1 = b1.id; - idB2 = b2.id; - idB3 = b3.id; - - pdgB1 = b1.pdg; - pdgB2 = b2.pdg; - pdgB3 = b3.pdg; - - chargeB1 = b1.charge; - chargeB2 = b2.charge; - chargeB3 = b3.charge; - - xB1 = b1.x; - yB1 = b1.y; - zB1 = b1.z; - tB1 = b1.t; - pxB1 = b1.px; - pyB1 = b1.py; - pzB1 = b1.pz; - - xB2 = b2.x; - yB2 = b2.y; - zB2 = b2.z; - tB2 = b2.t; - pxB2 = b2.px; - pyB2 = b2.py; - pzB2 = b2.pz; - - xB3 = b3.x; - yB3 = b3.y; - zB3 = b3.z; - tB3 = b3.t; - pxB3 = b3.px; - pyB3 = b3.py; - pzB3 = b3.pz; - - treeBoundState->Fill(); + switch (std::abs(pdg)) { + case PDG_t::kProton: + return o2::constants::physics::MassProton; + case PDG_t::kNeutron: + return o2::constants::physics::MassNeutron; + case PDG_t::kLambda0: + return o2::constants::physics::MassLambda0; + default: + return 0.0; + } } - bool passTwoBodySkim(Particle const& b1, Particle const& b2) + // Apply a momentum-space skim using Jacobi momenta in the candidate rest frame + template + bool passThreeBodySkim(const ReducedPart& b1, const ReducedPart& b2, const ReducedPart& b3) { - auto p4B1 = makeFourVector(b1); - auto p4B2 = makeFourVector(b2); - - auto candidate = p4B1 + p4B2; - - // BoostToCM() returns the boost vector needed to go from the lab frame - // to the candidate center-of-mass frame. - auto betaCandidate = candidate.BoostToCM(); - ROOT::Math::Boost boostToRest(betaCandidate); - - auto p4B1Star = boostToRest(p4B1); - auto p4B2Star = boostToRest(p4B2); - - /* - // Space-time propagation block. - // Not needed for the present skim because the selection is based only on - // the relative momentum pRho in the candidate rest frame. - // - // This block can be re-enabled later if one wants to apply coordinate-space - // or full phase-space coalescence selections, e.g. cuts on relative distance - // after propagation to a common time. - - ROOT::Math::XYZTVector x4B1{b1.x, b1.y, b1.z, b1.t}; - ROOT::Math::XYZTVector x4B2{b2.x, b2.y, b2.z, b2.t}; - - auto x4B1Star = boostToRest(x4B1); - auto x4B2Star = boostToRest(x4B2); - - const double tMax = std::max(x4B1Star.T(), x4B2Star.T()); - - ROOT::Math::XYZVector rB1{x4B1Star.X(), x4B1Star.Y(), x4B1Star.Z()}; - ROOT::Math::XYZVector rB2{x4B2Star.X(), x4B2Star.Y(), x4B2Star.Z()}; - - ROOT::Math::XYZVector vB1 = p4B1Star.Vect() / p4B1Star.E(); - ROOT::Math::XYZVector vB2 = p4B2Star.Vect() / p4B2Star.E(); - - rB1 += vB1 * (tMax - x4B1Star.T()); - rB2 += vB2 * (tMax - x4B2Star.T()); - */ - - const ROOT::Math::XYZVector pB1Star = p4B1Star.Vect(); - const ROOT::Math::XYZVector pB2Star = p4B2Star.Vect(); - - const double m12 = b1.mass + b2.mass; - - ROOT::Math::XYZVector pRho = (b2.mass * pB1Star - b1.mass * pB2Star) / m12; - - if (pRho.R() > pRhoMax) { + // Constituent masses from PDG codes + double m1 = massFromPdg(b1.pdgPart); + double m2 = massFromPdg(b2.pdgPart); + double m3 = massFromPdg(b3.pdgPart); + if (m1 <= 0.0 || m2 <= 0.0 || m3 <= 0.0) { return false; } - return true; - } - bool passThreeBodySkim(Particle const& b1, Particle const& b2, Particle const& b3) - { - auto p4B1 = makeFourVector(b1); - auto p4B2 = makeFourVector(b2); - auto p4B3 = makeFourVector(b3); + // Constituent four-momenta in the laboratory frame + auto p4B1 = ROOT::Math::PxPyPzMVector{b1.pxPart, b1.pyPart, b1.pzPart, m1}; + auto p4B2 = ROOT::Math::PxPyPzMVector{b2.pxPart, b2.pyPart, b2.pzPart, m2}; + auto p4B3 = ROOT::Math::PxPyPzMVector{b3.pxPart, b3.pyPart, b3.pzPart, m3}; + // Candidate four-momentum auto candidate = p4B1 + p4B2 + p4B3; - // BoostToCM() returns the boost vector needed to go from the lab frame - // to the candidate center-of-mass frame. + // Boost to the candidate rest frame auto betaCandidate = candidate.BoostToCM(); ROOT::Math::Boost boostToRest(betaCandidate); + // Constituent momenta in the candidate rest frame auto p4B1Star = boostToRest(p4B1); auto p4B2Star = boostToRest(p4B2); auto p4B3Star = boostToRest(p4B3); - - /* - // Space-time propagation block. - // Not needed for the present skim because the selection is based only on - // the Jacobi relative momenta pRho and pLambda in the candidate rest frame. - // - // This block can be re-enabled later if one wants to apply coordinate-space - // or full phase-space coalescence selections, e.g. cuts on relative distances - // after propagation to a common time. - - ROOT::Math::XYZTVector x4B1{b1.x, b1.y, b1.z, b1.t}; - ROOT::Math::XYZTVector x4B2{b2.x, b2.y, b2.z, b2.t}; - ROOT::Math::XYZTVector x4B3{b3.x, b3.y, b3.z, b3.t}; - - auto x4B1Star = boostToRest(x4B1); - auto x4B2Star = boostToRest(x4B2); - auto x4B3Star = boostToRest(x4B3); - - const double tMax = std::max({x4B1Star.T(), x4B2Star.T(), x4B3Star.T()}); - - ROOT::Math::XYZVector rB1{x4B1Star.X(), x4B1Star.Y(), x4B1Star.Z()}; - ROOT::Math::XYZVector rB2{x4B2Star.X(), x4B2Star.Y(), x4B2Star.Z()}; - ROOT::Math::XYZVector rB3{x4B3Star.X(), x4B3Star.Y(), x4B3Star.Z()}; - - ROOT::Math::XYZVector vB1 = p4B1Star.Vect() / p4B1Star.E(); - ROOT::Math::XYZVector vB2 = p4B2Star.Vect() / p4B2Star.E(); - ROOT::Math::XYZVector vB3 = p4B3Star.Vect() / p4B3Star.E(); - - rB1 += vB1 * (tMax - x4B1Star.T()); - rB2 += vB2 * (tMax - x4B2Star.T()); - rB3 += vB3 * (tMax - x4B3Star.T()); - */ - const ROOT::Math::XYZVector pB1Star = p4B1Star.Vect(); const ROOT::Math::XYZVector pB2Star = p4B2Star.Vect(); const ROOT::Math::XYZVector pB3Star = p4B3Star.Vect(); - const double m12 = b1.mass + b2.mass; - const double mTot = b1.mass + b2.mass + b3.mass; + // Reduced masses entering the Jacobi coordinates + const double m12 = m1 + m2; + const double mTot = m1 + m2 + m3; - ROOT::Math::XYZVector pRho = (b2.mass * pB1Star - b1.mass * pB2Star) / m12; - ROOT::Math::XYZVector pLambda = (b3.mass * (pB1Star + pB2Star) - m12 * pB3Star) / mTot; + // Jacobi momenta in the candidate rest frame + ROOT::Math::XYZVector pRho = (m2 * pB1Star - m1 * pB2Star) / m12; + ROOT::Math::XYZVector pLambda = (m3 * (pB1Star + pB2Star) - m12 * pB3Star) / mTot; - if (pRho.R() > pRhoMax) { - return false; - } - - if (pLambda.R() > pLambdaMax) { + // Reject candidates outside the momentum-space acceptance + if (pRho.R() > pRhoMax || pLambda.R() > pLambdaMax) { return false; } return true; } - void buildPairs(std::vector const& b1Candidates, - std::vector const& b2Candidates) - { - if (b1Candidates.empty() || b2Candidates.empty()) { - return; - } - - for (const auto& b1 : b1Candidates) { - for (const auto& b2 : b2Candidates) { - if (!passTwoBodySkim(b1, b2)) { - continue; - } + // Group MC particles by their associated MC collision + Preslice mcParticlesPerMcCollision = o2::aod::mcparticle::mcCollisionId; - storePair(b1, b2); - } - } - } - - void buildTriplets(std::vector const& b1Candidates, - std::vector const& b2Candidates, - std::vector const& b3Candidates, - bool identicalB1B2 = false, - bool identicalB2B3 = false) + // Process Hypertriton + void processHypertriton(aod::McCollisions const& collisions, aod::McParticles const& mcParticles) { - if (b1Candidates.empty() || b2Candidates.empty() || b3Candidates.empty()) { - return; - } - - if (identicalB1B2) { - for (size_t iB1 = 0; iB1 < b1Candidates.size(); ++iB1) { - for (size_t iB2 = iB1 + 1; iB2 < b2Candidates.size(); ++iB2) { - for (const auto& b3 : b3Candidates) { - const auto& b1 = b1Candidates[iB1]; - const auto& b2 = b2Candidates[iB2]; - - if (!passThreeBodySkim(b1, b2, b3)) { - continue; - } - - storeTriplet(b1, b2, b3); - } - } - } - return; - } + // Loop over MC collisions + for (const auto& collision : collisions) { - if (identicalB2B3) { - for (const auto& b1 : b1Candidates) { - for (size_t iB2 = 0; iB2 < b2Candidates.size(); ++iB2) { - for (size_t iB3 = iB2 + 1; iB3 < b3Candidates.size(); ++iB3) { - const auto& b2 = b2Candidates[iB2]; - const auto& b3 = b3Candidates[iB3]; - - if (!passThreeBodySkim(b1, b2, b3)) { - continue; - } + // Event counter before selections + registry.fill(HIST("eventCounter"), 0.5); - storeTriplet(b1, b2, b3); - } - } - } - return; - } + // Apply z-vertex selection + if (std::fabs(collision.posZ()) > zVtxMax) + continue; - for (const auto& b1 : b1Candidates) { - for (const auto& b2 : b2Candidates) { - for (const auto& b3 : b3Candidates) { - if (!passThreeBodySkim(b1, b2, b3)) { - continue; - } + // Event counter after z-vertex cut + registry.fill(HIST("eventCounter"), 1.5); - storeTriplet(b1, b2, b3); - } - } - } - } + // To be implemented: maybe add INEL>0 selection - template - void selectGeneratedParticles(McParticles const& mcParticlesThisMcColl, - std::vector& protons, - std::vector& antiProtons, - std::vector& neutrons, - std::vector& antiNeutrons, - std::vector& lambdas, - std::vector& antiLambdas) - { - for (const auto& particle : mcParticlesThisMcColl) { - if (!particle.isPhysicalPrimary()) { - continue; - } + // Get particles in this MC collision + const auto mcParticlesThisMcColl = mcParticles.sliceBy(mcParticlesPerMcCollision, collision.globalIndex()); - if (std::fabs(particle.eta()) > etaMax) { - continue; - } + // Event ID + eventID = collision.globalIndex(); - const int pdg = particle.pdgCode(); - - if (pdg == PDG_t::kProton) { - protons.emplace_back(particle); - } else if (pdg == PDG_t::kProtonBar) { - antiProtons.emplace_back(particle); - } else if (pdg == PDG_t::kNeutron) { - neutrons.emplace_back(particle); - } else if (pdg == PDG_t::kNeutronBar) { - antiNeutrons.emplace_back(particle); - } else if (pdg == PDG_t::kLambda0) { - lambdas.emplace_back(particle); - } else if (pdg == PDG_t::kLambda0Bar) { - antiLambdas.emplace_back(particle); - } - } - } + // Containers for candidate constituents + std::vector protons; + std::vector neutrons; + std::vector lambdas; - void fillDeuteron(std::vector const& protons, - std::vector const& antiProtons, - std::vector const& neutrons, - std::vector const& antiNeutrons) - { - // Deuteron: p + n - buildPairs(protons, neutrons); + // Loop over MC particles + for (const auto& particle : mcParticlesThisMcColl) { - // anti-deuteron: anti-p + anti-n - buildPairs(antiProtons, antiNeutrons); - } + // Select physical primary particles only + if (!particle.isPhysicalPrimary()) + continue; - void fillHelium3(std::vector const& protons, - std::vector const& antiProtons, - std::vector const& neutrons, - std::vector const& antiNeutrons) - { - // 3He: p + p + n - buildTriplets(protons, protons, neutrons, true, false); + // Select particles with |eta|< eta_max + if (std::fabs(particle.eta()) > etaMax) + continue; - // anti-3He: anti-p + anti-p + anti-n - buildTriplets(antiProtons, antiProtons, antiNeutrons, true, false); - } + // Store protons + if (std::abs(particle.pdgCode()) == PDG_t::kProton) { + protons.push_back({particle.globalIndex(), particle.pdgCode(), particle.vx(), particle.vy(), particle.vz(), particle.vt(), particle.px(), particle.py(), particle.pz()}); + } - void fillTriton(std::vector const& protons, - std::vector const& antiProtons, - std::vector const& neutrons, - std::vector const& antiNeutrons) - { - // Triton: p + n + n - buildTriplets(protons, neutrons, neutrons, false, true); + // Store neutrons + if (std::abs(particle.pdgCode()) == PDG_t::kNeutron) { + neutrons.push_back({particle.globalIndex(), particle.pdgCode(), particle.vx(), particle.vy(), particle.vz(), particle.vt(), particle.px(), particle.py(), particle.pz()}); + } - // anti-triton: anti-p + anti-n + anti-n - buildTriplets(antiProtons, antiNeutrons, antiNeutrons, false, true); - } + // Store lambdas + if (std::abs(particle.pdgCode()) == PDG_t::kLambda0) { + lambdas.push_back({particle.globalIndex(), particle.pdgCode(), particle.vx(), particle.vy(), particle.vz(), particle.vt(), particle.px(), particle.py(), particle.pz()}); + } + }// end of loop over MC particles - void fillHypertriton(std::vector const& protons, - std::vector const& antiProtons, - std::vector const& neutrons, - std::vector const& antiNeutrons, - std::vector const& lambdas, - std::vector const& antiLambdas) - { - // Hypertriton: p + n + Lambda - buildTriplets(protons, neutrons, lambdas); + // Reject events that do not contain at least one proton, one neutron, and one lambda + if (protons.empty() || neutrons.empty() || lambdas.empty()) + continue; - // anti-hypertriton: anti-p + anti-n + anti-Lambda - buildTriplets(antiProtons, antiNeutrons, antiLambdas); - } + // Event counter: events containing all three constituent species + registry.fill(HIST("eventCounter"), 2.5); - bool hasDeuteronCandidates(std::vector const& protons, - std::vector const& antiProtons, - std::vector const& neutrons, - std::vector const& antiNeutrons) const - { - return (!protons.empty() && !neutrons.empty()) || - (!antiProtons.empty() && !antiNeutrons.empty()); - } + // Count matter and antimatter candidates in the event + Long64_t nCandidatesMatter(0); + Long64_t nCandidatesAntiMatter(0); + for (const auto& proton : protons) { + for (const auto& neutron : neutrons) { + for (const auto& lambda : lambdas) { - bool hasHelium3Candidates(std::vector const& protons, - std::vector const& antiProtons, - std::vector const& neutrons, - std::vector const& antiNeutrons) const - { - constexpr std::size_t minimumSizeContainer = 2; - return (protons.size() >= minimumSizeContainer && !neutrons.empty()) || - (antiProtons.size() >= minimumSizeContainer && !antiNeutrons.empty()); - } + const bool isMatter = proton.pdgPart > 0 && neutron.pdgPart > 0 && lambda.pdgPart > 0; + const bool isAntimatter = proton.pdgPart < 0 && neutron.pdgPart < 0 && lambda.pdgPart < 0; - bool hasTritonCandidates(std::vector const& protons, - std::vector const& antiProtons, - std::vector const& neutrons, - std::vector const& antiNeutrons) const - { - constexpr std::size_t minimumSizeContainer = 2; - return (!protons.empty() && neutrons.size() >= minimumSizeContainer) || - (!antiProtons.empty() && antiNeutrons.size() >= minimumSizeContainer); - } + // Skip mixed-sign combinations + if (!isMatter && !isAntimatter) + continue; + const bool passSkim = passThreeBodySkim(proton, neutron, lambda); - bool hasHypertritonCandidates(std::vector const& protons, - std::vector const& antiProtons, - std::vector const& neutrons, - std::vector const& antiNeutrons, - std::vector const& lambdas, - std::vector const& antiLambdas) const - { - return (!protons.empty() && !neutrons.empty() && !lambdas.empty()) || - (!antiProtons.empty() && !antiNeutrons.empty() && !antiLambdas.empty()); - } + if (isMatter && passSkim) + nCandidatesMatter++; + if (isAntimatter && passSkim) + nCandidatesAntiMatter++; + } + } + } - void processCoalescence(aod::McCollision const& collision, aod::McParticles const& mcParticles) - { - registry.fill(HIST("eventCounter"), 0.5); + // Reject events with no accepted candidates + if (nCandidatesMatter == 0 && nCandidatesAntiMatter == 0) + continue; - if (std::fabs(collision.posZ()) > zVtx) { - return; - } + // Event counter: number of events with at least one candidate + registry.fill(HIST("eventCounter"), 3.5); - registry.fill(HIST("eventCounter"), 1.5); + // Store number of candidates per event + nMatterCandidatesPerEvent = nCandidatesMatter; + nAntimatterCandidatesPerEvent = nCandidatesAntiMatter; - eventID = collision.globalIndex(); + // Loop over accepted triplets and fill the tree + for (const auto& proton : protons) { + for (const auto& neutron : neutrons) { + for (const auto& lambda : lambdas) { - std::vector protons; - std::vector antiProtons; - std::vector neutrons; - std::vector antiNeutrons; - std::vector lambdas; - std::vector antiLambdas; + const bool isMatter = proton.pdgPart > 0 && neutron.pdgPart > 0 && lambda.pdgPart > 0; + const bool isAntimatter = proton.pdgPart < 0 && neutron.pdgPart < 0 && lambda.pdgPart < 0; - const auto mcParticlesThisMcColl = mcParticles.sliceBy(mcParticlesPerMcCollision, collision.globalIndex()); + // Skip mixed-sign combinations + if (!isMatter && !isAntimatter) + continue; - selectGeneratedParticles(mcParticlesThisMcColl, - protons, - antiProtons, - neutrons, - antiNeutrons, - lambdas, - antiLambdas); + // Apply momentum-space selection + const bool passSkim = passThreeBodySkim(proton, neutron, lambda); + if (!passSkim) + continue; - switch (boundStateSpecies) { - case kDeuteron: - if (!hasDeuteronCandidates(protons, antiProtons, neutrons, antiNeutrons)) { - return; - } - registry.fill(HIST("eventCounter"), 2.5); - fillDeuteron(protons, antiProtons, neutrons, antiNeutrons); - break; - case kHelium3: - if (!hasHelium3Candidates(protons, antiProtons, neutrons, antiNeutrons)) { - return; - } - registry.fill(HIST("eventCounter"), 2.5); - fillHelium3(protons, antiProtons, neutrons, antiNeutrons); - break; - case kTriton: - if (!hasTritonCandidates(protons, antiProtons, neutrons, antiNeutrons)) { - return; - } - registry.fill(HIST("eventCounter"), 2.5); - fillTriton(protons, antiProtons, neutrons, antiNeutrons); - break; - case kHypertriton: - if (!hasHypertritonCandidates(protons, antiProtons, neutrons, antiNeutrons, lambdas, antiLambdas)) { - return; + // Fill tree + idB1 = proton.idPart; + pdgB1 = proton.pdgPart; + xB1 = proton.xPart; + yB1 = proton.yPart; + zB1 = proton.zPart; + tB1 = proton.tPart; + pxB1 = proton.pxPart; + pyB1 = proton.pyPart; + pzB1 = proton.pzPart; + + idB2 = neutron.idPart; + pdgB2 = neutron.pdgPart; + xB2 = neutron.xPart; + yB2 = neutron.yPart; + zB2 = neutron.zPart; + tB2 = neutron.tPart; + pxB2 = neutron.pxPart; + pyB2 = neutron.pyPart; + pzB2 = neutron.pzPart; + + idB3 = lambda.idPart; + pdgB3 = lambda.pdgPart; + xB3 = lambda.xPart; + yB3 = lambda.yPart; + zB3 = lambda.zPart; + tB3 = lambda.tPart; + pxB3 = lambda.pxPart; + pyB3 = lambda.pyPart; + pzB3 = lambda.pzPart; + + treeBoundState->Fill(); + } } - registry.fill(HIST("eventCounter"), 2.5); - fillHypertriton(protons, antiProtons, neutrons, antiNeutrons, lambdas, antiLambdas); - break; - default: - LOGF(fatal, "Unknown boundStateSpecies=%d. Use 0 = deuteron, 1 = triton, 2 = helium3, 3 = hypertriton", static_cast(boundStateSpecies)); - } + } + }// end of loop over mc collisions } - PROCESS_SWITCH(CoalescenceTreeProducer, processCoalescence, "process coalescence", true); + PROCESS_SWITCH(CoalescenceTreeProducer, processHypertriton, "process hypertriton", true); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { - return WorkflowSpec{ - adaptAnalysisTask(cfgc)}; + return WorkflowSpec{adaptAnalysisTask(cfgc)}; } From 913296af5b17647671da4cbbac172f59eaff9243 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Wed, 3 Jun 2026 09:54:50 +0000 Subject: [PATCH 2/2] Please consider the following formatting changes --- PWGLF/TableProducer/Nuspex/coalescenceTreeProducer.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/PWGLF/TableProducer/Nuspex/coalescenceTreeProducer.cxx b/PWGLF/TableProducer/Nuspex/coalescenceTreeProducer.cxx index 1dbc95139c6..339fd71a655 100644 --- a/PWGLF/TableProducer/Nuspex/coalescenceTreeProducer.cxx +++ b/PWGLF/TableProducer/Nuspex/coalescenceTreeProducer.cxx @@ -95,7 +95,7 @@ struct CoalescenceTreeProducer { Long64_t nMatterCandidatesPerEvent = 0; // Number of matter candidates per event Long64_t nAntimatterCandidatesPerEvent = 0; // Number of antimatter candidates per event Long64_t idB1 = 0, idB2 = 0, idB3 = 0; // Constituent identifiers - int pdgB1 = 0, pdgB2 = 0, pdgB3 = 0; // Constituent PDG codes + int pdgB1 = 0, pdgB2 = 0, pdgB3 = 0; // Constituent PDG codes // Constituent space-time coordinates and momentum components float xB1 = 0.f, yB1 = 0.f, zB1 = 0.f, tB1 = 0.f, pxB1 = 0.f, pyB1 = 0.f, pzB1 = 0.f; @@ -276,7 +276,7 @@ struct CoalescenceTreeProducer { if (std::abs(particle.pdgCode()) == PDG_t::kLambda0) { lambdas.push_back({particle.globalIndex(), particle.pdgCode(), particle.vx(), particle.vy(), particle.vz(), particle.vt(), particle.px(), particle.py(), particle.pz()}); } - }// end of loop over MC particles + } // end of loop over MC particles // Reject events that do not contain at least one proton, one neutron, and one lambda if (protons.empty() || neutrons.empty() || lambdas.empty()) @@ -371,7 +371,7 @@ struct CoalescenceTreeProducer { } } } - }// end of loop over mc collisions + } // end of loop over mc collisions } PROCESS_SWITCH(CoalescenceTreeProducer, processHypertriton, "process hypertriton", true); };