From 13336cd05f7fd7f69599ac3fa81dd5bf0655465c Mon Sep 17 00:00:00 2001 From: creetz16 Date: Thu, 10 Sep 2026 15:53:15 +0200 Subject: [PATCH] Add centrality and gen daughter momenta --- PWGLF/DataModel/Vtx3BodyTables.h | 38 ++++++-- .../Nuspex/decay3bodybuilder.cxx | 42 ++++----- .../Nuspex/trackedHypertritonRecoTask.cxx | 87 +++++++++++++------ 3 files changed, 110 insertions(+), 57 deletions(-) diff --git a/PWGLF/DataModel/Vtx3BodyTables.h b/PWGLF/DataModel/Vtx3BodyTables.h index bb0dcc89fe1..fbfb258a0f1 100644 --- a/PWGLF/DataModel/Vtx3BodyTables.h +++ b/PWGLF/DataModel/Vtx3BodyTables.h @@ -17,6 +17,8 @@ #ifndef PWGLF_DATAMODEL_VTX3BODYTABLES_H_ #define PWGLF_DATAMODEL_VTX3BODYTABLES_H_ +#include "PWGLF/DataModel/LFHypernucleiTables.h" + #include "Common/Core/RecoDecay.h" #include @@ -128,12 +130,15 @@ DECLARE_SOA_COLUMN(GenCt, genCt, float); // generated Ct of the hy DECLARE_SOA_COLUMN(GenPhi, genPhi, float); // generated Phi of the hypertriton DECLARE_SOA_COLUMN(GenEta, genEta, float); // Eta of the hypertriton DECLARE_SOA_COLUMN(GenRap, genRap, float); // generated rapidity of the hypertriton -DECLARE_SOA_COLUMN(GenPPr, genPPr, float); //! generated momentum proton daughter particle -DECLARE_SOA_COLUMN(GenPPi, genPPi, float); //! generated momentum pion daughter particle -DECLARE_SOA_COLUMN(GenPDe, genPDe, float); //! generated momentum deuteron daughter particle -DECLARE_SOA_COLUMN(GenPtPr, genPtPr, float); //! generated transverse momentum proton daughter particle -DECLARE_SOA_COLUMN(GenPtPi, genPtPi, float); //! generated transverse momentum pion daughter particle -DECLARE_SOA_COLUMN(GenPtDe, genPtDe, float); //! generated transverse momentum deuteron daughter particle +DECLARE_SOA_COLUMN(GenPxPr, genPxPr, float); //! generated Px of proton daughter particle +DECLARE_SOA_COLUMN(GenPyPr, genPyPr, float); //! generated Py of proton daughter particle +DECLARE_SOA_COLUMN(GenPzPr, genPzPr, float); //! generated Pz of proton daughter particle +DECLARE_SOA_COLUMN(GenPxPi, genPxPi, float); //! generated Px of pion daughter particle +DECLARE_SOA_COLUMN(GenPyPi, genPyPi, float); //! generated Py of pion daughter particle +DECLARE_SOA_COLUMN(GenPzPi, genPzPi, float); //! generated Pz of pion daughter particle +DECLARE_SOA_COLUMN(GenPxDe, genPxDe, float); //! generated Px of deuteron daughter particle +DECLARE_SOA_COLUMN(GenPyDe, genPyDe, float); //! generated Py of deuteron daughter particle +DECLARE_SOA_COLUMN(GenPzDe, genPzDe, float); //! generated Pz of deuteron daughter particle DECLARE_SOA_COLUMN(MotherPdgCode, motherPdgCode, int); //! PDG code of the mother particle DECLARE_SOA_COLUMN(PrPdgCode, prPdgCode, int); //! MC particle proton PDG code DECLARE_SOA_COLUMN(PiPdgCode, piPdgCode, int); //! MC particle pion PDG code @@ -271,6 +276,22 @@ DECLARE_SOA_TABLE(Vtx3BodyTrackedInfo, "AOD", "VTX3BODYTR", //! vtx3body::ITSTrackDCAXYToSV, vtx3body::ITSTrackDCAZToSV); +// collision info table +DECLARE_SOA_TABLE(Vtx3BodyCollision, "AOD", "VTX3BODYCOLL", //! + o2::soa::Index<>, + hyperrec::CentralityFT0A, hyperrec::CentralityFT0C, hyperrec::CentralityFT0M, + hyperrec::TrackOccupancyInTimeRange, hyperrec::Ft0cOccupancyInTimeRange, + hyperrec::XPrimVtx, hyperrec::YPrimVtx, hyperrec::ZPrimVtx, + hyperrec::RunNumber); + +DECLARE_SOA_TABLE(McVtx3BodyCollision, "AOD", "MCVTX3BODYCOLL", //! + o2::soa::Index<>, + hyperrec::IsRecoMCCollision, + hyperrec::CentralityFT0A, hyperrec::CentralityFT0C, hyperrec::CentralityFT0M, + hyperrec::TrackOccupancyInTimeRange, hyperrec::Ft0cOccupancyInTimeRange, + hyperrec::XPrimVtx, hyperrec::YPrimVtx, hyperrec::ZPrimVtx, + hyperrec::RunNumber); + // MC candidate table for analysis DECLARE_SOA_TABLE(McVtx3BodyDatas, "AOD", "MC3BODYDATA", //! o2::soa::Index<>, @@ -303,8 +324,9 @@ DECLARE_SOA_TABLE(McVtx3BodyDatas, "AOD", "MC3BODYDATA", //! vtx3body::GenX, vtx3body::GenY, vtx3body::GenZ, vtx3body::GenCt, vtx3body::GenPhi, vtx3body::GenEta, vtx3body::GenRap, - vtx3body::GenPPr, vtx3body::GenPPi, vtx3body::GenPDe, - vtx3body::GenPtPr, vtx3body::GenPtPi, vtx3body::GenPtDe, + vtx3body::GenPxPr, vtx3body::GenPyPr, vtx3body::GenPzPr, + vtx3body::GenPxPi, vtx3body::GenPyPi, vtx3body::GenPzPi, + vtx3body::GenPxDe, vtx3body::GenPyDe, vtx3body::GenPzDe, vtx3body::IsReco, vtx3body::MotherLabel, vtx3body::MotherPdgCode, diff --git a/PWGLF/TableProducer/Nuspex/decay3bodybuilder.cxx b/PWGLF/TableProducer/Nuspex/decay3bodybuilder.cxx index b71c38eb678..2f96571ae61 100644 --- a/PWGLF/TableProducer/Nuspex/decay3bodybuilder.cxx +++ b/PWGLF/TableProducer/Nuspex/decay3bodybuilder.cxx @@ -227,12 +227,9 @@ struct decay3bodyBuilder { float genPhi; float genEta; float genRapidity; - float genMomProton; - float genMomPion; - float genMomDeuteron; - float genPtProton; - float genPtPion; - float genPtDeuteron; + std::array genMomProton{0.f}; + std::array genMomPion{0.f}; + std::array genMomDeuteron{0.f}; bool isReco; int motherLabel; int motherPdgCode; @@ -822,12 +819,9 @@ struct decay3bodyBuilder { this3BodyMCInfo.daughterPiPdgCode = mcTrackPion.pdgCode(); this3BodyMCInfo.daughterDePdgCode = mcTrackDeuteron.pdgCode(); this3BodyMCInfo.isDeuteronPrimary = mcTrackDeuteron.isPhysicalPrimary(); - this3BodyMCInfo.genMomProton = mcTrackProton.p(); - this3BodyMCInfo.genMomPion = mcTrackPion.p(); - this3BodyMCInfo.genMomDeuteron = mcTrackDeuteron.p(); - this3BodyMCInfo.genPtProton = mcTrackProton.pt(); - this3BodyMCInfo.genPtPion = mcTrackPion.pt(); - this3BodyMCInfo.genPtDeuteron = mcTrackDeuteron.pt(); + this3BodyMCInfo.genMomProton = {mcTrackProton.px(), mcTrackProton.py(), mcTrackProton.pz()}; + this3BodyMCInfo.genMomPion = {mcTrackPion.px(), mcTrackPion.py(), mcTrackPion.pz()}; + this3BodyMCInfo.genMomDeuteron = {mcTrackDeuteron.px(), mcTrackDeuteron.py(), mcTrackDeuteron.pz()}; // daughters are matched to MC, now we check if reco mother is true H3L/Anti-H3l and decayed via three-body decay this3BodyMCInfo.motherLabel = checkH3LTruth(mcTrackProton, mcTrackPion, mcTrackDeuteron); // returns global index of mother if true H3L/Anti-H3L mother decaying via three-body decay, otherwise negative value for background @@ -905,17 +899,14 @@ struct decay3bodyBuilder { // get daughters for (const auto& mcparticleDaughter : mcparticle.template daughters_as()) { if (std::abs(mcparticleDaughter.pdgCode()) == PDG_t::kProton) { // proton - this3BodyMCInfo.genMomProton = mcparticleDaughter.p(); - this3BodyMCInfo.genPtProton = mcparticleDaughter.pt(); + this3BodyMCInfo.genMomProton = {mcparticleDaughter.px(), mcparticleDaughter.py(), mcparticleDaughter.pz()}; this3BodyMCInfo.daughterPrPdgCode = mcparticleDaughter.pdgCode(); this3BodyMCInfo.genDecVtx = {mcparticleDaughter.vx(), mcparticleDaughter.vy(), mcparticleDaughter.vz()}; } else if (std::abs(mcparticleDaughter.pdgCode()) == PDG_t::kPiPlus) { // pion - this3BodyMCInfo.genMomPion = mcparticleDaughter.p(); - this3BodyMCInfo.genPtPion = mcparticleDaughter.pt(); + this3BodyMCInfo.genMomPion = {mcparticleDaughter.px(), mcparticleDaughter.py(), mcparticleDaughter.pz()}; this3BodyMCInfo.daughterPiPdgCode = mcparticleDaughter.pdgCode(); } else if (std::abs(mcparticleDaughter.pdgCode()) == o2::constants::physics::Pdg::kDeuteron) { // deuteron - this3BodyMCInfo.genMomDeuteron = mcparticleDaughter.p(); - this3BodyMCInfo.genPtDeuteron = mcparticleDaughter.pt(); + this3BodyMCInfo.genMomDeuteron = {mcparticleDaughter.px(), mcparticleDaughter.py(), mcparticleDaughter.pz()}; this3BodyMCInfo.daughterDePdgCode = mcparticleDaughter.pdgCode(); this3BodyMCInfo.isDeuteronPrimary = mcparticleDaughter.isPhysicalPrimary(); } @@ -957,8 +948,9 @@ struct decay3bodyBuilder { this3BodyMCInfo.genDecVtx[0], this3BodyMCInfo.genDecVtx[1], this3BodyMCInfo.genDecVtx[2], this3BodyMCInfo.genCt, mcparticle.phi(), mcparticle.eta(), mcparticle.y(), - this3BodyMCInfo.genMomProton, this3BodyMCInfo.genMomPion, this3BodyMCInfo.genMomDeuteron, - this3BodyMCInfo.genPtProton, this3BodyMCInfo.genPtPion, this3BodyMCInfo.genPtDeuteron, + this3BodyMCInfo.genMomProton[0], this3BodyMCInfo.genMomProton[1], this3BodyMCInfo.genMomProton[2], + this3BodyMCInfo.genMomPion[0], this3BodyMCInfo.genMomPion[1], this3BodyMCInfo.genMomPion[2], + this3BodyMCInfo.genMomDeuteron[0], this3BodyMCInfo.genMomDeuteron[1], this3BodyMCInfo.genMomDeuteron[2], this3BodyMCInfo.isReco, mcparticle.globalIndex(), // motherLabel mcparticle.pdgCode(), // motherPdgCode @@ -1153,8 +1145,9 @@ struct decay3bodyBuilder { this3BodyMCInfo.genDecVtx[0], this3BodyMCInfo.genDecVtx[1], this3BodyMCInfo.genDecVtx[2], this3BodyMCInfo.genCt, this3BodyMCInfo.genPhi, this3BodyMCInfo.genEta, this3BodyMCInfo.genRapidity, - this3BodyMCInfo.genMomProton, this3BodyMCInfo.genMomPion, this3BodyMCInfo.genMomDeuteron, - this3BodyMCInfo.genPtProton, this3BodyMCInfo.genPtPion, this3BodyMCInfo.genPtDeuteron, + this3BodyMCInfo.genMomProton[0], this3BodyMCInfo.genMomProton[1], this3BodyMCInfo.genMomProton[2], + this3BodyMCInfo.genMomPion[0], this3BodyMCInfo.genMomPion[1], this3BodyMCInfo.genMomPion[2], + this3BodyMCInfo.genMomDeuteron[0], this3BodyMCInfo.genMomDeuteron[1], this3BodyMCInfo.genMomDeuteron[2], this3BodyMCInfo.isReco, this3BodyMCInfo.motherLabel, this3BodyMCInfo.motherPdgCode, @@ -1280,8 +1273,9 @@ struct decay3bodyBuilder { mcInfo.genDecVtx[0] = -1., mcInfo.genDecVtx[1] = -1., mcInfo.genDecVtx[2] = -1.; mcInfo.genCt = -1.; mcInfo.genPhi = -1., mcInfo.genEta = -1., mcInfo.genRapidity = -1.; - mcInfo.genMomProton = -1., mcInfo.genMomPion = -1., mcInfo.genMomDeuteron = -1.; - mcInfo.genPtProton = -1., mcInfo.genPtPion = -1., mcInfo.genPtDeuteron = -1.; + mcInfo.genMomProton[0] = -1., mcInfo.genMomProton[1] = -1., mcInfo.genMomProton[2] = -1.; + mcInfo.genMomPion[0] = -1., mcInfo.genMomPion[1] = -1., mcInfo.genMomPion[2] = -1.; + mcInfo.genMomDeuteron[0] = -1., mcInfo.genMomDeuteron[1] = -1., mcInfo.genMomDeuteron[2] = -1.; mcInfo.isReco = false; mcInfo.motherPdgCode = 0; mcInfo.daughterPrPdgCode = -1, mcInfo.daughterPiPdgCode = -1, mcInfo.daughterDePdgCode = -1; diff --git a/PWGLF/TableProducer/Nuspex/trackedHypertritonRecoTask.cxx b/PWGLF/TableProducer/Nuspex/trackedHypertritonRecoTask.cxx index da2ca3ee1ee..c8bc8cdb20d 100644 --- a/PWGLF/TableProducer/Nuspex/trackedHypertritonRecoTask.cxx +++ b/PWGLF/TableProducer/Nuspex/trackedHypertritonRecoTask.cxx @@ -101,7 +101,9 @@ struct TrackedHypertritonRecoTask { Produces vtx3BodyDatas; Produces vtx3BodyCovs; Produces vtx3BodyTrackedInfo; + Produces vtx3BodyCollision; Produces mcVtx3BodyDatas; + Produces mcVtx3BodyCollision; Service ccdb{}; HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject}; @@ -238,12 +240,9 @@ struct TrackedHypertritonRecoTask { float genPhi = -1.f; float genEta = -1.f; float genRapidity = -1.f; - float genMomentumProton = -1.f; - float genMomentumPion = -1.f; - float genMomentumDeuteron = -1.f; - float genPtProton = -1.f; - float genPtPion = -1.f; - float genPtDeuteron = -1.f; + std::array genMomProton{-1.f, -1.f, -1.f}; + std::array genMomPion{-1.f, -1.f, -1.f}; + std::array genMomDeuteron{-1.f, -1.f, -1.f}; bool isReco = true; int motherLabel = -1; int motherPdgCode = 0; @@ -251,6 +250,7 @@ struct TrackedHypertritonRecoTask { int pionPdgCode = -1; int deuteronPdgCode = -1; bool isDeuteronPrimary = false; + bool isRecoMCCollision = false; bool survivedEventSelection = false; }; @@ -521,12 +521,9 @@ struct TrackedHypertritonRecoTask { info.pionPdgCode = mcPion.pdgCode(); info.deuteronPdgCode = mcDeuteron.pdgCode(); info.isDeuteronPrimary = mcDeuteron.isPhysicalPrimary(); - info.genMomentumProton = mcProton.p(); - info.genMomentumPion = mcPion.p(); - info.genMomentumDeuteron = mcDeuteron.p(); - info.genPtProton = mcProton.pt(); - info.genPtPion = mcPion.pt(); - info.genPtDeuteron = mcDeuteron.pt(); + info.genMomProton = {mcProton.px(), mcProton.py(), mcProton.pz()}; + info.genMomPion = {mcPion.px(), mcPion.py(), mcPion.pz()}; + info.genMomDeuteron = {mcDeuteron.px(), mcDeuteron.py(), mcDeuteron.pz()}; const int motherLabel = findCommonMother(mcProton, mcPion, mcDeuteron); if (motherLabel < 0) { @@ -548,7 +545,8 @@ struct TrackedHypertritonRecoTask { info.genPhi = mother.phi(); info.genEta = mother.eta(); info.genRapidity = mother.y(); - if (mother.mcCollisionId() >= 0 && mother.mcCollisionId() < static_cast(survivedMCEventSelection.size())) { + if (mother.mcCollisionId() >= 0 && mother.mcCollisionId() < static_cast(recoCollisionForMC.size())) { + info.isRecoMCCollision = recoCollisionForMC[mother.mcCollisionId()] >= 0; info.survivedEventSelection = survivedMCEventSelection[mother.mcCollisionId()]; } return info; @@ -790,7 +788,8 @@ struct TrackedHypertritonRecoTask { return dcaInfo; } - void fillThreeBodyTables() + template + void fillThreeBodyTables(TCollision const& collision) { const auto& candidate = builder3Body.decay3body; vtx3BodyDatas(static_cast(candidate.sign), @@ -816,6 +815,10 @@ struct TrackedHypertritonRecoTask { static_cast(candidate.pidForTrackingDeuteron)); vtx3BodyCovs(candidate.covProton.data(), candidate.covPion.data(), candidate.covDeuteron.data(), candidate.covariance.data()); vtx3BodyTrackedInfo(candidate.itsTrackDCAToSV[0], candidate.itsTrackDCAToSV[1]); + vtx3BodyCollision(collision.centFT0A(), collision.centFT0C(), collision.centFT0M(), + collision.trackOccupancyInTimeRange(), collision.ft0cOccupancyInTimeRange(), + collision.posX(), collision.posY(), collision.posZ(), + runNumber); } void fillThreeBodyMCTable(ThreeBodyMCInfo const& info) @@ -845,8 +848,9 @@ struct TrackedHypertritonRecoTask { info.genMomentum[0], info.genMomentum[1], info.genMomentum[2], info.genDecayVertex[0], info.genDecayVertex[1], info.genDecayVertex[2], info.genCt, info.genPhi, info.genEta, info.genRapidity, - info.genMomentumProton, info.genMomentumPion, info.genMomentumDeuteron, - info.genPtProton, info.genPtPion, info.genPtDeuteron, + info.genMomProton[0], info.genMomProton[1], info.genMomProton[2], + info.genMomPion[0], info.genMomPion[1], info.genMomPion[2], + info.genMomDeuteron[0], info.genMomDeuteron[1], info.genMomDeuteron[2], static_cast(info.isReco), info.motherLabel, info.motherPdgCode, info.protonPdgCode, info.pionPdgCode, info.deuteronPdgCode, info.isDeuteronPrimary, static_cast(info.survivedEventSelection)); @@ -879,8 +883,9 @@ struct TrackedHypertritonRecoTask { info.genMomentum[0], info.genMomentum[1], info.genMomentum[2], info.genDecayVertex[0], info.genDecayVertex[1], info.genDecayVertex[2], info.genCt, info.genPhi, info.genEta, info.genRapidity, - info.genMomentumProton, info.genMomentumPion, info.genMomentumDeuteron, - info.genPtProton, info.genPtPion, info.genPtDeuteron, + info.genMomProton[0], info.genMomProton[1], info.genMomProton[2], + info.genMomPion[0], info.genMomPion[1], info.genMomPion[2], + info.genMomDeuteron[0], info.genMomDeuteron[1], info.genMomDeuteron[2], 0, info.motherLabel, info.motherPdgCode, info.protonPdgCode, info.pionPdgCode, info.deuteronPdgCode, info.isDeuteronPrimary, static_cast(info.survivedEventSelection)); @@ -972,7 +977,7 @@ struct TrackedHypertritonRecoTask { } fillQAHistograms(trackProton, trackPion, trackDeuteron, true); - fillThreeBodyTables(); + fillThreeBodyTables(collision); } } } @@ -1077,6 +1082,12 @@ struct TrackedHypertritonRecoTask { } fillThreeBodyMCTable(mcInfo); + mcVtx3BodyCollision(true, // IsRecoMCCollision + collision.centFT0A(), collision.centFT0C(), collision.centFT0M(), + collision.trackOccupancyInTimeRange(), collision.ft0cOccupancyInTimeRange(), + collision.posX(), collision.posY(), collision.posZ(), + runNumber); + if (mcInfo.motherLabel >= 0) { reconstructedThreeBody[mcInfo.motherLabel] = true; } @@ -1117,19 +1128,16 @@ struct TrackedHypertritonRecoTask { hasPionTwoBody = true; hasPionThreeBody = true; threeBodyInfo.pionPdgCode = daughter.pdgCode(); - threeBodyInfo.genMomentumPion = daughter.p(); - threeBodyInfo.genPtPion = daughter.pt(); + threeBodyInfo.genMomPion = {daughter.px(), daughter.py(), daughter.pz()}; } else if (daughter.pdgCode() == sign * PDG_t::kProton) { hasProton = true; threeBodyInfo.protonPdgCode = daughter.pdgCode(); - threeBodyInfo.genMomentumProton = daughter.p(); - threeBodyInfo.genPtProton = daughter.pt(); + threeBodyInfo.genMomProton = {daughter.px(), daughter.py(), daughter.pz()}; threeBodyInfo.genDecayVertex = {daughter.vx(), daughter.vy(), daughter.vz()}; } else if (daughter.pdgCode() == sign * constants::physics::Pdg::kDeuteron) { hasDeuteron = true; threeBodyInfo.deuteronPdgCode = daughter.pdgCode(); - threeBodyInfo.genMomentumDeuteron = daughter.p(); - threeBodyInfo.genPtDeuteron = daughter.pt(); + threeBodyInfo.genMomDeuteron = {daughter.px(), daughter.py(), daughter.pz()}; threeBodyInfo.isDeuteronPrimary = daughter.isPhysicalPrimary(); } } @@ -1178,11 +1186,40 @@ struct TrackedHypertritonRecoTask { } if (hasProton && hasPionThreeBody && hasDeuteron && !reconstructedThreeBody[mother.globalIndex()]) { + float centralityFT0A = -1.f; + float centralityFT0C = -1.f; + float centralityFT0M = -1.f; + int trackOccupancyInTimeRange = -1; + float ft0cOccupancyInTimeRange = -1.f; + float primaryVertexX = -1.f; + float primaryVertexY = -1.f; + float primaryVertexZ = -1.f; + bool isRecoMCCollision = false; + if (mother.mcCollisionId() >= 0 && mother.mcCollisionId() < static_cast(recoCollisionForMC.size())) { + const int recoCollisionId = recoCollisionForMC[mother.mcCollisionId()]; + isRecoMCCollision = recoCollisionId >= 0; + if (isRecoMCCollision) { + const auto collision = collisions.rawIteratorAt(recoCollisionId); + centralityFT0A = collision.centFT0A(); + centralityFT0C = collision.centFT0C(); + centralityFT0M = collision.centFT0M(); + trackOccupancyInTimeRange = collision.trackOccupancyInTimeRange(); + ft0cOccupancyInTimeRange = collision.ft0cOccupancyInTimeRange(); + primaryVertexX = collision.posX(); + primaryVertexY = collision.posY(); + primaryVertexZ = collision.posZ(); + } + } threeBodyInfo.genCt = RecoDecay::sqrtSumOfSquares(threeBodyInfo.genDecayVertex[0] - mother.vx(), threeBodyInfo.genDecayVertex[1] - mother.vy(), threeBodyInfo.genDecayVertex[2] - mother.vz()) * constants::physics::MassHyperTriton / (mother.p() + 1.e-10f); fillGeneratedThreeBodyMCTable(threeBodyInfo); + mcVtx3BodyCollision(isRecoMCCollision, + centralityFT0A, centralityFT0C, centralityFT0M, + trackOccupancyInTimeRange, ft0cOccupancyInTimeRange, + primaryVertexX, primaryVertexY, primaryVertexZ, + runNumber); } } }