diff --git a/PWGLF/DataModel/LFSigmaHadTables.h b/PWGLF/DataModel/LFSigmaHadTables.h index 77c91c30349..984aa9005a5 100644 --- a/PWGLF/DataModel/LFSigmaHadTables.h +++ b/PWGLF/DataModel/LFSigmaHadTables.h @@ -36,14 +36,17 @@ DECLARE_SOA_COLUMN(PyHad, pyHad, float); //! Py of the hadron cand DECLARE_SOA_COLUMN(PzHad, pzHad, float); //! Pz of the hadron candidate DECLARE_SOA_COLUMN(NSigmaTPCHad, nSigmaTPCHad, float); //! Number of sigmas for the hadron candidate from Sigma kink in TPC DECLARE_SOA_COLUMN(NSigmaTOFHad, nSigmaTOFHad, float); //! Number of sigmas for the hadron candidate from Sigma kink in TOF +DECLARE_SOA_COLUMN(Multiplicity, multiplicity, float); //! Collision multiplicity // MC Columns -DECLARE_SOA_COLUMN(SigmaPDG, sigmaPDG, int); //! PDG code of the Sigma daughter -DECLARE_SOA_COLUMN(DaughterPDG, daughterPDG, int); //! PDG code of the kink daughter -DECLARE_SOA_COLUMN(HadPDG, hadPDG, int); //! PDG code of the hadron candidate -DECLARE_SOA_COLUMN(SigmaGenPt, sigmaGenPt, float); //! Generated pT of the Sigma candidate -DECLARE_SOA_COLUMN(HadGenPt, hadGenPt, float); //! Generated pT of the hadron candidate -DECLARE_SOA_COLUMN(GenKStar, genKStar, float); //! Generated k* of the Sigma-hadron pair +DECLARE_SOA_COLUMN(SigmaPDG, sigmaPDG, int); //! PDG code of the Sigma candidate +DECLARE_SOA_COLUMN(DaughterPDG, daughterPDG, int); //! PDG code of the kink daughter +DECLARE_SOA_COLUMN(HadPDG, hadPDG, int); //! PDG code of the hadron candidate +DECLARE_SOA_COLUMN(SigmaMotherPDG, sigmaMotherPDG, int); //! PDG code of the direct mother of the Sigma +DECLARE_SOA_COLUMN(SigmaPartonicMotherPDG, sigmaPartonicMotherPDG, int); //! PDG code of the partonic (quark/gluon) ancestor of the Sigma +DECLARE_SOA_COLUMN(SigmaGenPt, sigmaGenPt, float); //! Generated pT of the Sigma candidate +DECLARE_SOA_COLUMN(HadGenPt, hadGenPt, float); //! Generated pT of the hadron candidate +DECLARE_SOA_COLUMN(GenKStar, genKStar, float); //! Generated k* of the Sigma-hadron pair } // namespace sigmaproton @@ -52,15 +55,16 @@ DECLARE_SOA_TABLE(SigmaProtonCands, "AOD", "SIGMAPROTONCANDS", sigmaproton::ChargeSigma, kinkcand::PxMoth, kinkcand::PyMoth, kinkcand::PzMoth, kinkcand::PxDaug, kinkcand::PyDaug, kinkcand::PzDaug, sigmaproton::SigmaDecRad, sigmaproton::SigmaCosPA, sigmaproton::ChargeHad, sigmaproton::PxHad, sigmaproton::PyHad, sigmaproton::PzHad, - sigmaproton::NSigmaTPCHad, sigmaproton::NSigmaTOFHad); + sigmaproton::NSigmaTPCHad, sigmaproton::NSigmaTOFHad, sigmaproton::Multiplicity); DECLARE_SOA_TABLE(SigmaProtonMCCands, "AOD", "SIGMAPROTONMCCANDS", o2::soa::Index<>, sigmaproton::ChargeSigma, kinkcand::PxMoth, kinkcand::PyMoth, kinkcand::PzMoth, kinkcand::PxDaug, kinkcand::PyDaug, kinkcand::PzDaug, sigmaproton::SigmaDecRad, sigmaproton::SigmaCosPA, sigmaproton::ChargeHad, sigmaproton::PxHad, sigmaproton::PyHad, sigmaproton::PzHad, - sigmaproton::NSigmaTPCHad, sigmaproton::NSigmaTOFHad, + sigmaproton::NSigmaTPCHad, sigmaproton::NSigmaTOFHad, sigmaproton::Multiplicity, sigmaproton::SigmaPDG, sigmaproton::DaughterPDG, sigmaproton::HadPDG, + sigmaproton::SigmaMotherPDG, sigmaproton::SigmaPartonicMotherPDG, sigmaproton::SigmaGenPt, sigmaproton::HadGenPt, sigmaproton::GenKStar); } // namespace o2::aod diff --git a/PWGLF/TableProducer/Strangeness/sigmaHadCorr.cxx b/PWGLF/TableProducer/Strangeness/sigmaHadCorr.cxx index 3eaee6351b3..1884caad6a7 100644 --- a/PWGLF/TableProducer/Strangeness/sigmaHadCorr.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaHadCorr.cxx @@ -17,6 +17,7 @@ #include "PWGLF/DataModel/LFSigmaHadTables.h" #include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" #include "Common/DataModel/PIDResponseTOF.h" #include "Common/DataModel/PIDResponseTPC.h" @@ -37,6 +38,7 @@ #include #include +#include #include #include @@ -50,8 +52,8 @@ using namespace o2::framework::expressions; using TracksFull = soa::Join; using TracksFullMC = soa::Join; -using CollisionsFull = soa::Join; -using CollisionsFullMC = soa::Join; +using CollisionsFull = soa::Join; +using CollisionsFullMC = soa::Join; struct sigmaHadCand { @@ -86,6 +88,11 @@ struct sigmaHadCand { int kinkDauID = -1; // ID of the pion from Sigma decay in MC int sigmaID = -1; // ID of the Sigma candidate in MC int hadID = -1; // ID of the hadron candidate in MC + + int sigmaMotherPDG = -999; // PDG of the direct mother of the Sigma in MC + int sigmaPartonicMotherPDG = -999; // PDG of the first or last partonic ancestor of the Sigma in MC + + float multiplicity = -1.f; }; struct sigmaHadCorrTask { @@ -124,6 +131,9 @@ struct sigmaHadCorrTask { Configurable fillOutputTree{"fillOutputTree", true, "If true, fill the output tree with Sigma-hadron candidates"}; Configurable fillSparseInvMassKstar{"fillSparseInvMassKstar", false, "If true, fill THn with invmass, k*, sigma charge, proton charge, sigma decay radius, cosPA, sigma pt"}; + Configurable saveMultiplicity{"saveMultiplicity", false, "If true, save collision multiplicity in the output table"}; + Configurable useMultNTracksPV{"useMultNTracksPV", false, "If true, use multNTracksPV for multiplicity and mixing bins; if false, use numContrib"}; + Configurable findLastPartonicMother{"findLastPartonicMother", true, "If true, store the initial hard-scattering parton (last partonic mother). If false, store the last parton before hadronization (first partonic mother)"}; ConfigurableAxis CfgVtxBins{"CfgVtxBins", {10, -10, 10}, "Mixing bins - z-vertex"}; ConfigurableAxis CfgMultBins{"CfgMultBins", {VARIABLE_WIDTH, 0.0, 40.0, 80.0, 500.0}, "Mixing bins - number of contributor"}; @@ -278,6 +288,103 @@ struct sigmaHadCorrTask { return doSigmaPion ? o2::constants::physics::MassPionCharged : o2::constants::physics::MassProton; } + // Section of functions to retrieve MC truth information for the Sigma candidate (adapted from PWGCF/Femto) + // Get the PDG code of the direct mother of the Sigma candidate + template + int getSigmaMotherPDG(const TMcParticle& mcParticle, const TMcParticles& mcParticles) + { + if (!mcParticle.has_mothers()) { + return -999; + } + auto motherIds = mcParticle.mothersIds(); + int mcMotherIndex = motherIds.front(); + if (mcMotherIndex < 0 || mcMotherIndex >= static_cast(mcParticles.size())) { + return -999; + } + return mcParticles.iteratorAt(mcMotherIndex).pdgCode(); + } + + // Walk up the decay chain and return the PDG of the first quark or gluon ancestor found + template + int findFirstPartonicMotherPDG(const TMcParticle& mcParticle, const TMcParticles& mcParticles) + { + if (!mcParticle.has_mothers()) { + return -999; + } + auto motherIds = mcParticle.mothersIds(); + const int defaultMotherSize = 2; + std::vector allMotherIds; + if (motherIds.size() == defaultMotherSize && motherIds[1] > motherIds[0]) { + for (int i = motherIds[0]; i <= motherIds[1]; i++) + allMotherIds.push_back(i); + } else { + for (const int& id : motherIds) + allMotherIds.push_back(id); + } + for (const int& i : allMotherIds) { + if (i < 0 || i >= static_cast(mcParticles.size())) + continue; + const auto& mother = mcParticles.iteratorAt(i); + int pdgAbs = std::abs(mother.pdgCode()); + if (pdgAbs <= PDG_t::kTop || pdgAbs == PDG_t::kGluon) { + return mother.pdgCode(); + } + int found = findFirstPartonicMotherPDG(mother, mcParticles); + if (found != -999) + return found; + } + return -999; + } + + // Walk up the decay chain iteratively and return the PDG of the last quark or gluon before beam remnants + template + int findLastPartonicMotherPDG(const TMcParticle& mcParticle, const TMcParticles& mcParticles) + { + int lastPartonPDG = -999; + int64_t currentIndex = mcParticle.globalIndex(); + while (currentIndex >= 0 && currentIndex < static_cast(mcParticles.size())) { + const auto& current = mcParticles.iteratorAt(currentIndex); + if (!current.has_mothers()) + break; + auto motherIds = current.mothersIds(); + int nextIndex = -1; + const int defaultMotherSize = 2; + if (motherIds.size() == defaultMotherSize && motherIds[1] > motherIds[0]) { + nextIndex = motherIds[0]; + } else { + for (const int& id : motherIds) { + if (id >= 0 && id < static_cast(mcParticles.size())) { + nextIndex = id; + break; + } + } + } + if (nextIndex < 0 || nextIndex >= static_cast(mcParticles.size())) + break; + const auto& mother = mcParticles.iteratorAt(nextIndex); + int pdgAbs = std::abs(mother.pdgCode()); + int status = std::abs(o2::mcgenstatus::getGenStatusCode(mother.statusCode())); + bool isParton = (pdgAbs <= PDG_t::kTop || pdgAbs == PDG_t::kGluon); + const int isBeamParticleLowerLimit = 11; + const int isBeamParticleUpperLimit = 19; + bool isBeam = (status >= isBeamParticleLowerLimit && status <= isBeamParticleUpperLimit); + if (isBeam) + return lastPartonPDG; + if (isParton) + lastPartonPDG = mother.pdgCode(); + currentIndex = nextIndex; + } + return lastPartonPDG; + } + + int getPartonicMotherPDG(const auto& mcParticle, const auto& mcParticles) + { + if (findLastPartonicMother.value) { + return findLastPartonicMotherPDG(mcParticle, mcParticles); + } + return findFirstPartonicMotherPDG(mcParticle, mcParticles); + } + float getSigmaMassForKstar() { return doSigmaMinus ? o2::constants::physics::MassSigmaMinus : o2::constants::physics::MassSigmaPlus; @@ -455,6 +562,7 @@ struct sigmaHadCorrTask { candidate.sigmaID = sigmaCand.trackMothId(); candidate.kinkDauID = sigmaCand.trackDaugId(); candidate.hadID = hadTrack.globalIndex(); + candidate.multiplicity = saveMultiplicity.value ? (useMultNTracksPV.value ? collision.multNTracksPV() : collision.numContrib()) : -1.f; float kStar = getKStar(sigmaPRecal[0], sigmaPRecal[1], sigmaPRecal[2], candidate.pxHad, candidate.pyHad, candidate.pzHad); if (kStar > cutMaxKStar) { @@ -510,7 +618,8 @@ struct sigmaHadCorrTask { candidate.pyHad, candidate.pzHad, candidate.nSigmaTPCHad, - candidate.nSigmaTOFHad); + candidate.nSigmaTOFHad, + candidate.multiplicity); } } } @@ -519,54 +628,64 @@ struct sigmaHadCorrTask { // Processing Event Mixing SliceCache cache; - using BinningType = ColumnBinningPolicy; - BinningType colBinning{{CfgVtxBins, CfgMultBins}, true}; + using BinningTypeNumContrib = ColumnBinningPolicy; + using BinningTypeMultNTracksPV = ColumnBinningPolicy; + BinningTypeNumContrib colBinningNumContrib{{CfgVtxBins, CfgMultBins}, true}; + BinningTypeMultNTracksPV colBinningPVMult{{CfgVtxBins, CfgMultBins}, true}; void processMixedEvent(const CollisionsFull& collisions, const aod::KinkCands& kinkCands, const TracksFull& tracks) { - for (auto const& [collision1, collision2] : - selfCombinations(colBinning, nEvtMixingBkg, -1, collisions, collisions)) { - if (collision1.index() == collision2.index()) - continue; - - sigmaHadCandidates.clear(); - if (std::abs(collision1.posZ()) > cutzvertex || !collision1.sel8()) { - continue; - } - if (std::abs(collision2.posZ()) > cutzvertex || !collision2.sel8()) { - continue; + if (useMultNTracksPV.value) { + for (auto const& [collision1, collision2] : + selfCombinations(colBinningPVMult, nEvtMixingBkg, -1, collisions, collisions)) { + if (collision1.index() == collision2.index()) + continue; + sigmaHadCandidates.clear(); + if (std::abs(collision1.posZ()) > cutzvertex || !collision1.sel8()) + continue; + if (std::abs(collision2.posZ()) > cutzvertex || !collision2.sel8()) + continue; + auto kinkCands_c1 = kinkCands.sliceBy(kinkCandsPerCollisionPreslice, collision1.globalIndex()); + auto tracks_c2 = tracks.sliceBy(tracksPerCollisionPreslice, collision2.globalIndex()); + fillTreeAndHistograms(kinkCands_c1, tracks, tracks_c2, collision1, false); + if (fillOutputTree) { + for (const auto& candidate : sigmaHadCandidates) { + outputDataTable(candidate.sigmaCharge, candidate.sigmaPx, candidate.sigmaPy, candidate.sigmaPz, + candidate.sigmaDauPx, candidate.sigmaDauPy, candidate.sigmaDauPz, + candidate.sigmaDecRadius, candidate.sigmaCosPA, candidate.chargeHad, + candidate.pxHad, candidate.pyHad, candidate.pzHad, + candidate.nSigmaTPCHad, candidate.nSigmaTOFHad, candidate.multiplicity); + } + } } - auto kinkCands_c1 = kinkCands.sliceBy(kinkCandsPerCollisionPreslice, collision1.globalIndex()); - auto tracks_c1 = tracks.sliceBy(tracksPerCollisionPreslice, collision1.globalIndex()); - auto tracks_c2 = tracks.sliceBy(tracksPerCollisionPreslice, collision2.globalIndex()); - fillTreeAndHistograms(kinkCands_c1, tracks, tracks_c2, collision1, false); - - if (fillOutputTree) { - // Fill output table - for (const auto& candidate : sigmaHadCandidates) { - outputDataTable(candidate.sigmaCharge, - candidate.sigmaPx, - candidate.sigmaPy, - candidate.sigmaPz, - candidate.sigmaDauPx, - candidate.sigmaDauPy, - candidate.sigmaDauPz, - candidate.sigmaDecRadius, - candidate.sigmaCosPA, - candidate.chargeHad, - candidate.pxHad, - candidate.pyHad, - candidate.pzHad, - candidate.nSigmaTPCHad, - candidate.nSigmaTOFHad); + } else { + for (auto const& [collision1, collision2] : + selfCombinations(colBinningNumContrib, nEvtMixingBkg, -1, collisions, collisions)) { + if (collision1.index() == collision2.index()) + continue; + sigmaHadCandidates.clear(); + if (std::abs(collision1.posZ()) > cutzvertex || !collision1.sel8()) + continue; + if (std::abs(collision2.posZ()) > cutzvertex || !collision2.sel8()) + continue; + auto kinkCands_c1 = kinkCands.sliceBy(kinkCandsPerCollisionPreslice, collision1.globalIndex()); + auto tracks_c2 = tracks.sliceBy(tracksPerCollisionPreslice, collision2.globalIndex()); + fillTreeAndHistograms(kinkCands_c1, tracks, tracks_c2, collision1, false); + if (fillOutputTree) { + for (const auto& candidate : sigmaHadCandidates) { + outputDataTable(candidate.sigmaCharge, candidate.sigmaPx, candidate.sigmaPy, candidate.sigmaPz, + candidate.sigmaDauPx, candidate.sigmaDauPy, candidate.sigmaDauPz, + candidate.sigmaDecRadius, candidate.sigmaCosPA, candidate.chargeHad, + candidate.pxHad, candidate.pyHad, candidate.pzHad, + candidate.nSigmaTPCHad, candidate.nSigmaTOFHad, candidate.multiplicity); + } } } } - LOG(debug) << "Processing mixed event"; } PROCESS_SWITCH(sigmaHadCorrTask, processMixedEvent, "Process Mixed event", false); - void processSameEventMC(CollisionsFullMC const& collisions, aod::KinkCands const& kinkCands, TracksFullMC const& tracks, aod::McParticles const&) + void processSameEventMC(CollisionsFullMC const& collisions, aod::KinkCands const& kinkCands, TracksFullMC const& tracks, aod::McParticles const& mcParticles) { for (auto const& collision : collisions) { @@ -594,6 +713,8 @@ struct sigmaHadCorrTask { auto pdgSigma = mcPartSigma.pdgCode(); auto pdgSigmaDau = mcLabelSigmaDau.has_mcParticle() ? mcPartSigmaDau.pdgCode() : -999; auto pdgHad = mcLabelHad.has_mcParticle() ? mcPartHad.pdgCode() : -999; + auto sigmaMotherPDG = getSigmaMotherPDG(mcPartSigma, mcParticles); + auto sigmaPartonicMotherPDG = getPartonicMotherPDG(mcPartSigma, mcParticles); float sigmaPtGen = std::hypot(mcPartSigma.px(), mcPartSigma.py()); float hadPtGen = std::hypot(mcPartHad.px(), mcPartHad.py()); @@ -631,9 +752,12 @@ struct sigmaHadCorrTask { candidate.pzHad, candidate.nSigmaTPCHad, candidate.nSigmaTOFHad, + candidate.multiplicity, pdgSigma, pdgSigmaDau, pdgHad, + sigmaMotherPDG, + sigmaPartonicMotherPDG, sigmaPtGen, hadPtGen, kStarGen); @@ -643,86 +767,107 @@ struct sigmaHadCorrTask { } PROCESS_SWITCH(sigmaHadCorrTask, processSameEventMC, "Process Same event MC", false); - void processMixedEventMC(const CollisionsFullMC& collisions, const aod::KinkCands& kinkCands, const TracksFullMC& tracks, const aod::McParticles&) + void processMixedEventMC(const CollisionsFullMC& collisions, const aod::KinkCands& kinkCands, const TracksFullMC& tracks, const aod::McParticles& mcParticles) { - for (auto const& [collision1, collision2] : - selfCombinations(colBinning, nEvtMixingBkg, -1, collisions, collisions)) { - if (collision1.index() == collision2.index()) - continue; - - sigmaHadCandidates.clear(); - if (std::abs(collision1.posZ()) > cutzvertex || !collision1.sel8()) { - continue; - } - if (std::abs(collision2.posZ()) > cutzvertex || !collision2.sel8()) { - continue; - } - auto kinkCands_c1 = kinkCands.sliceBy(kinkCandsPerCollisionPreslice, collision1.globalIndex()); - auto tracks_c1 = tracks.sliceBy(tracksPerCollisionPreslice, collision1.globalIndex()); - auto tracks_c2 = tracks.sliceBy(tracksPerCollisionPreslice, collision2.globalIndex()); - fillTreeAndHistograms(kinkCands_c1, tracks, tracks_c2, collision1, true); - - for (const auto& candidate : sigmaHadCandidates) { - auto mcLabelSigma = tracks.rawIteratorAt(candidate.sigmaID); - auto mcLabelSigmaDau = tracks.rawIteratorAt(candidate.kinkDauID); - auto mcLabelHad = tracks.rawIteratorAt(candidate.hadID); - - if (!mcLabelSigma.has_mcParticle() || !mcLabelSigmaDau.has_mcParticle() || !mcLabelHad.has_mcParticle()) { - continue; // Skip candidates where MC truth is not available - } - - auto mcPartSigma = mcLabelSigma.mcParticle_as(); - auto mcPartSigmaDau = mcLabelSigmaDau.mcParticle_as(); - auto mcPartHad = mcLabelHad.mcParticle_as(); - auto pdgSigma = mcPartSigma.pdgCode(); - auto pdgSigmaDau = mcLabelSigmaDau.has_mcParticle() ? mcPartSigmaDau.pdgCode() : -999; - auto pdgHad = mcLabelHad.has_mcParticle() ? mcPartHad.pdgCode() : -999; - float sigmaPtGen = std::hypot(mcPartSigma.px(), mcPartSigma.py()); - float hadPtGen = std::hypot(mcPartHad.px(), mcPartHad.py()); - float kStarGen = getKStar(mcPartSigma.px(), mcPartSigma.py(), mcPartSigma.pz(), mcPartHad.px(), mcPartHad.py(), mcPartHad.pz()); - - if (fillSparseInvMassKstar) { - auto sigmaMomForKstar = getSigmaMomentumForKstar(candidate.sigmaPx, candidate.sigmaPy, candidate.sigmaPz, - candidate.sigmaDauPx, candidate.sigmaDauPy, candidate.sigmaDauPz); - float kStarRec = getKStar(sigmaMomForKstar[0], sigmaMomForKstar[1], sigmaMomForKstar[2], candidate.pxHad, candidate.pyHad, candidate.pzHad); - float sigmaPtUsed = std::hypot(sigmaMomForKstar[0], sigmaMomForKstar[1]); - rSigmaHad.fill(HIST("hSparseSigmaHadMC"), - candidate.sigmaMass, - kStarRec, - candidate.sigmaCharge, - candidate.chargeHad, - candidate.sigmaDecRadius, - candidate.sigmaCosPA, - sigmaPtUsed, - kStarGen); + if (useMultNTracksPV.value) { + for (auto const& [collision1, collision2] : + selfCombinations(colBinningPVMult, nEvtMixingBkg, -1, collisions, collisions)) { + if (collision1.index() == collision2.index()) + continue; + sigmaHadCandidates.clear(); + if (std::abs(collision1.posZ()) > cutzvertex || !collision1.sel8()) + continue; + if (std::abs(collision2.posZ()) > cutzvertex || !collision2.sel8()) + continue; + auto kinkCands_c1 = kinkCands.sliceBy(kinkCandsPerCollisionPreslice, collision1.globalIndex()); + auto tracks_c2 = tracks.sliceBy(tracksPerCollisionPreslice, collision2.globalIndex()); + fillTreeAndHistograms(kinkCands_c1, tracks, tracks_c2, collision1, true); + for (const auto& candidate : sigmaHadCandidates) { + auto mcLabelSigma = tracks.rawIteratorAt(candidate.sigmaID); + auto mcLabelSigmaDau = tracks.rawIteratorAt(candidate.kinkDauID); + auto mcLabelHad = tracks.rawIteratorAt(candidate.hadID); + if (!mcLabelSigma.has_mcParticle() || !mcLabelSigmaDau.has_mcParticle() || !mcLabelHad.has_mcParticle()) + continue; + auto mcPartSigma = mcLabelSigma.mcParticle_as(); + auto mcPartSigmaDau = mcLabelSigmaDau.mcParticle_as(); + auto mcPartHad = mcLabelHad.mcParticle_as(); + auto pdgSigma = mcPartSigma.pdgCode(); + auto pdgSigmaDau = mcLabelSigmaDau.has_mcParticle() ? mcPartSigmaDau.pdgCode() : -999; + auto pdgHad = mcLabelHad.has_mcParticle() ? mcPartHad.pdgCode() : -999; + auto sigmaMotherPDG = getSigmaMotherPDG(mcPartSigma, mcParticles); + auto sigmaPartonicMotherPDG = getPartonicMotherPDG(mcPartSigma, mcParticles); + float sigmaPtGen = std::hypot(mcPartSigma.px(), mcPartSigma.py()); + float hadPtGen = std::hypot(mcPartHad.px(), mcPartHad.py()); + float kStarGen = getKStar(mcPartSigma.px(), mcPartSigma.py(), mcPartSigma.pz(), mcPartHad.px(), mcPartHad.py(), mcPartHad.pz()); + if (fillSparseInvMassKstar) { + auto sigmaMomForKstar = getSigmaMomentumForKstar(candidate.sigmaPx, candidate.sigmaPy, candidate.sigmaPz, + candidate.sigmaDauPx, candidate.sigmaDauPy, candidate.sigmaDauPz); + float kStarRec = getKStar(sigmaMomForKstar[0], sigmaMomForKstar[1], sigmaMomForKstar[2], candidate.pxHad, candidate.pyHad, candidate.pzHad); + float sigmaPtUsed = std::hypot(sigmaMomForKstar[0], sigmaMomForKstar[1]); + rSigmaHad.fill(HIST("hSparseSigmaHadMC"), candidate.sigmaMass, kStarRec, candidate.sigmaCharge, + candidate.chargeHad, candidate.sigmaDecRadius, candidate.sigmaCosPA, sigmaPtUsed, kStarGen); + } + if (fillOutputTree) { + outputDataTableMC(candidate.sigmaCharge, candidate.sigmaPx, candidate.sigmaPy, candidate.sigmaPz, + candidate.sigmaDauPx, candidate.sigmaDauPy, candidate.sigmaDauPz, + candidate.sigmaDecRadius, candidate.sigmaCosPA, candidate.chargeHad, + candidate.pxHad, candidate.pyHad, candidate.pzHad, + candidate.nSigmaTPCHad, candidate.nSigmaTOFHad, candidate.multiplicity, + pdgSigma, pdgSigmaDau, pdgHad, sigmaMotherPDG, sigmaPartonicMotherPDG, + sigmaPtGen, hadPtGen, kStarGen); + } } - - if (fillOutputTree) { - outputDataTableMC(candidate.sigmaCharge, - candidate.sigmaPx, - candidate.sigmaPy, - candidate.sigmaPz, - candidate.sigmaDauPx, - candidate.sigmaDauPy, - candidate.sigmaDauPz, - candidate.sigmaDecRadius, - candidate.sigmaCosPA, - candidate.chargeHad, - candidate.pxHad, - candidate.pyHad, - candidate.pzHad, - candidate.nSigmaTPCHad, - candidate.nSigmaTOFHad, - pdgSigma, - pdgSigmaDau, - pdgHad, - sigmaPtGen, - hadPtGen, - kStarGen); + } + } else { + for (auto const& [collision1, collision2] : + selfCombinations(colBinningNumContrib, nEvtMixingBkg, -1, collisions, collisions)) { + if (collision1.index() == collision2.index()) + continue; + sigmaHadCandidates.clear(); + if (std::abs(collision1.posZ()) > cutzvertex || !collision1.sel8()) + continue; + if (std::abs(collision2.posZ()) > cutzvertex || !collision2.sel8()) + continue; + auto kinkCands_c1 = kinkCands.sliceBy(kinkCandsPerCollisionPreslice, collision1.globalIndex()); + auto tracks_c2 = tracks.sliceBy(tracksPerCollisionPreslice, collision2.globalIndex()); + fillTreeAndHistograms(kinkCands_c1, tracks, tracks_c2, collision1, true); + for (const auto& candidate : sigmaHadCandidates) { + auto mcLabelSigma = tracks.rawIteratorAt(candidate.sigmaID); + auto mcLabelSigmaDau = tracks.rawIteratorAt(candidate.kinkDauID); + auto mcLabelHad = tracks.rawIteratorAt(candidate.hadID); + if (!mcLabelSigma.has_mcParticle() || !mcLabelSigmaDau.has_mcParticle() || !mcLabelHad.has_mcParticle()) + continue; + auto mcPartSigma = mcLabelSigma.mcParticle_as(); + auto mcPartSigmaDau = mcLabelSigmaDau.mcParticle_as(); + auto mcPartHad = mcLabelHad.mcParticle_as(); + auto pdgSigma = mcPartSigma.pdgCode(); + auto pdgSigmaDau = mcLabelSigmaDau.has_mcParticle() ? mcPartSigmaDau.pdgCode() : -999; + auto pdgHad = mcLabelHad.has_mcParticle() ? mcPartHad.pdgCode() : -999; + auto sigmaMotherPDG = getSigmaMotherPDG(mcPartSigma, mcParticles); + auto sigmaPartonicMotherPDG = getPartonicMotherPDG(mcPartSigma, mcParticles); + float sigmaPtGen = std::hypot(mcPartSigma.px(), mcPartSigma.py()); + float hadPtGen = std::hypot(mcPartHad.px(), mcPartHad.py()); + float kStarGen = getKStar(mcPartSigma.px(), mcPartSigma.py(), mcPartSigma.pz(), mcPartHad.px(), mcPartHad.py(), mcPartHad.pz()); + if (fillSparseInvMassKstar) { + auto sigmaMomForKstar = getSigmaMomentumForKstar(candidate.sigmaPx, candidate.sigmaPy, candidate.sigmaPz, + candidate.sigmaDauPx, candidate.sigmaDauPy, candidate.sigmaDauPz); + float kStarRec = getKStar(sigmaMomForKstar[0], sigmaMomForKstar[1], sigmaMomForKstar[2], candidate.pxHad, candidate.pyHad, candidate.pzHad); + float sigmaPtUsed = std::hypot(sigmaMomForKstar[0], sigmaMomForKstar[1]); + rSigmaHad.fill(HIST("hSparseSigmaHadMC"), candidate.sigmaMass, kStarRec, candidate.sigmaCharge, + candidate.chargeHad, candidate.sigmaDecRadius, candidate.sigmaCosPA, sigmaPtUsed, kStarGen); + } + if (fillOutputTree) { + outputDataTableMC(candidate.sigmaCharge, candidate.sigmaPx, candidate.sigmaPy, candidate.sigmaPz, + candidate.sigmaDauPx, candidate.sigmaDauPy, candidate.sigmaDauPz, + candidate.sigmaDecRadius, candidate.sigmaCosPA, candidate.chargeHad, + candidate.pxHad, candidate.pyHad, candidate.pzHad, + candidate.nSigmaTPCHad, candidate.nSigmaTOFHad, candidate.multiplicity, + pdgSigma, pdgSigmaDau, pdgHad, sigmaMotherPDG, sigmaPartonicMotherPDG, + sigmaPtGen, hadPtGen, kStarGen); + } } } } - LOG(debug) << "Processing mixed event MC"; } PROCESS_SWITCH(sigmaHadCorrTask, processMixedEventMC, "Process Mixed event MC", false); };