diff --git a/PWGLF/DataModel/LFKinkDecayTables.h b/PWGLF/DataModel/LFKinkDecayTables.h index 0b9137d64b1..797e6c9a546 100644 --- a/PWGLF/DataModel/LFKinkDecayTables.h +++ b/PWGLF/DataModel/LFKinkDecayTables.h @@ -25,6 +25,7 @@ #include #include +#include namespace o2::aod { @@ -212,15 +213,18 @@ DECLARE_SOA_COLUMN(GammaPdgCode, gammaPdgCode, int); //! PDG code DECLARE_SOA_COLUMN(GammaMotherPdgCode, gammaMotherPdgCode, int); //! PDG code of the photon's MC mother (expected: pi0) DECLARE_SOA_COLUMN(GammaGMotherPdgCode, gammaGMotherPdgCode, int); //! PDG code of the photon's MC grandmother (expected: Sigma+) -DECLARE_SOA_COLUMN(XDecVtxMC, xDecVtxMC, float); //! MC-truth Sigma+ decay vertex (x direction) -DECLARE_SOA_COLUMN(YDecVtxMC, yDecVtxMC, float); //! MC-truth Sigma+ decay vertex (y direction) -DECLARE_SOA_COLUMN(ZDecVtxMC, zDecVtxMC, float); //! MC-truth Sigma+ decay vertex (z direction) -DECLARE_SOA_COLUMN(PxProtonMC, pxProtonMC, float); //! MC-truth proton Px -DECLARE_SOA_COLUMN(PyProtonMC, pyProtonMC, float); //! MC-truth proton Py -DECLARE_SOA_COLUMN(PzProtonMC, pzProtonMC, float); //! MC-truth proton Pz -DECLARE_SOA_COLUMN(PxGammaMC, pxGammaMC, float); //! MC-truth momentum of the measured photon (Px) -DECLARE_SOA_COLUMN(PyGammaMC, pyGammaMC, float); //! MC-truth momentum of the measured photon (Py) -DECLARE_SOA_COLUMN(PzGammaMC, pzGammaMC, float); //! MC-truth momentum of the measured photon (Pz) +DECLARE_SOA_COLUMN(XDecVtxMC, xDecVtxMC, float); //! MC-truth Sigma+ decay vertex (x direction) +DECLARE_SOA_COLUMN(YDecVtxMC, yDecVtxMC, float); //! MC-truth Sigma+ decay vertex (y direction) +DECLARE_SOA_COLUMN(ZDecVtxMC, zDecVtxMC, float); //! MC-truth Sigma+ decay vertex (z direction) +DECLARE_SOA_COLUMN(PxSigmaPlusMC, pxSigmaPlusMC, float); //! MC-truth Sigma+ mother Px +DECLARE_SOA_COLUMN(PySigmaPlusMC, pySigmaPlusMC, float); //! MC-truth Sigma+ mother Py +DECLARE_SOA_COLUMN(PzSigmaPlusMC, pzSigmaPlusMC, float); //! MC-truth Sigma+ mother Pz +DECLARE_SOA_COLUMN(PxProtonMC, pxProtonMC, float); //! MC-truth proton Px +DECLARE_SOA_COLUMN(PyProtonMC, pyProtonMC, float); //! MC-truth proton Py +DECLARE_SOA_COLUMN(PzProtonMC, pzProtonMC, float); //! MC-truth proton Pz +DECLARE_SOA_COLUMN(PxGammaMC, pxGammaMC, float); //! MC-truth momentum of the measured photon (Px) +DECLARE_SOA_COLUMN(PyGammaMC, pyGammaMC, float); //! MC-truth momentum of the measured photon (Py) +DECLARE_SOA_COLUMN(PzGammaMC, pzGammaMC, float); //! MC-truth momentum of the measured photon (Pz) // DYNAMIC COLUMNS @@ -244,6 +248,12 @@ DECLARE_SOA_DYNAMIC_COLUMN(MassSigmaPlus, massSigmaPlus, //! Invariant mass of t std::array{pxGamma2, pyGamma2, pzGamma2}}, std::array{o2::constants::physics::MassProton, o2::constants::physics::MassGamma, o2::constants::physics::MassGamma}); }); +DECLARE_SOA_DYNAMIC_COLUMN(PtSigmaPlusMC, ptSigmaPlusMC, //! True pT of the Sigma+ mother + [](float pxSigmaPlusMC, float pySigmaPlusMC) -> float { return std::hypot(pxSigmaPlusMC, pySigmaPlusMC); }); + +DECLARE_SOA_DYNAMIC_COLUMN(YSigmaPlusMC, ySigmaPlusMC, //! True rapidity of the Sigma+ mother + [](float pxSigmaPlusMC, float pySigmaPlusMC, float pzSigmaPlusMC) -> float { return RecoDecay::y(std::array{pxSigmaPlusMC, pySigmaPlusMC, pzSigmaPlusMC}, o2::constants::physics::MassSigmaPlus); }); + } // namespace sigmapluscand DECLARE_SOA_TABLE(SigmaPlusCands, "AOD", "SIGMAPLUSCANDS", @@ -286,13 +296,16 @@ DECLARE_SOA_TABLE(SigmaPlusCandsMC, "AOD", "SIGMAPLUSMC", sigmapluscand::XDecVtxMC, sigmapluscand::YDecVtxMC, sigmapluscand::ZDecVtxMC, sigmapluscand::PxProtonMC, sigmapluscand::PyProtonMC, sigmapluscand::PzProtonMC, sigmapluscand::PxGammaMC, sigmapluscand::PyGammaMC, sigmapluscand::PzGammaMC, + sigmapluscand::PxSigmaPlusMC, sigmapluscand::PySigmaPlusMC, sigmapluscand::PzSigmaPlusMC, // dynamic columns sigmapluscand::PxSigmaPlus, sigmapluscand::PySigmaPlus, sigmapluscand::PzSigmaPlus, sigmapluscand::PtSigmaPlus, - sigmapluscand::MassSigmaPlus); + sigmapluscand::MassSigmaPlus, + sigmapluscand::PtSigmaPlusMC, + sigmapluscand::YSigmaPlusMC); } // namespace o2::aod diff --git a/PWGLF/TableProducer/Strangeness/sigmaHadCorr.cxx b/PWGLF/TableProducer/Strangeness/sigmaHadCorr.cxx index dc7df2f567a..58acb86cfb7 100644 --- a/PWGLF/TableProducer/Strangeness/sigmaHadCorr.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaHadCorr.cxx @@ -37,7 +37,8 @@ #include #include -#include +#include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) +#include #include #include @@ -56,7 +57,7 @@ using TracksFullMC = soa::Join; using CollisionsFullMC = soa::Join; -struct sigmaHadCand { +struct SigmaHadCand { float ptHad() const { @@ -96,16 +97,16 @@ struct sigmaHadCand { float multiplicity = -1.f; }; -struct sigmaHadCorrTask { +struct SigmaHadCorr { - std::vector sigmaHadCandidates; // Vector to store Sigma-hadron candidates Produces outputDataTable; // Output table for Sigma-hadron candidates Produces outputDataTableMC; // Output table for Sigma-hadron candidates in MC + Produces outputKinkCandsMC; // Single-Sigma-level MC truth record, filled before hadron pairing // Histograms are defined with HistogramRegistry HistogramRegistry rEventSelection{"eventSelection", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; HistogramRegistry rSigmaHad{"sigmaHad", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; // Configurable for event selection - Configurable cutzvertex{"cutZVertex", 10.0f, "Accepted z-vertex range (cm)"}; + Configurable cutZVertex{"cutZVertex", 10.0f, "Accepted z-vertex range (cm)"}; Configurable doSigmaPion{"doSigmaPion", false, "If true, pair Sigma with pions instead of protons"}; Configurable doSigmaMinus{"doSigmaMinus", true, "If true, pair Sigma- candidates, else Sigma+"}; @@ -119,7 +120,7 @@ struct sigmaHadCorrTask { Configurable alphaAPCut{"alphaAPCut", 0., "Alpha AP cut for Sigma candidates"}; Configurable qtAPCutLow{"qtAPCutLow", 0.15, "Lower qT AP cut for Sigma candidates (GeV/c)"}; Configurable qtAPCutHigh{"qtAPCutHigh", 0.2, "Upper qT AP cut for Sigma candidates (GeV/c)"}; - Configurable cutEtaDaught{"cutEtaDaughter", 0.8f, "Eta cut for daughter tracks"}; + Configurable cutEtaDaughter{"cutEtaDaughter", 0.8f, "Eta cut for daughter tracks"}; Configurable ptMinTOFKinkDau{"ptMinTOFKinkDau", 0.75f, "Minimum pT to require TOF for kink daughter PID (GeV/c)"}; Configurable applyTOFPIDKinkDaughter{"applyTOFPIDKinkDaughter", false, "If true, apply TOF PID cut to the kink daughter track"}; @@ -136,8 +137,8 @@ struct sigmaHadCorrTask { 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"}; + 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"}; Configurable nEvtMixingBkg{"nEvtMixingBkg", 5, "Number of events to mix for background reconstruction"}; Preslice kinkCandsPerCollisionPreslice = aod::kinkcand::collisionId; @@ -222,12 +223,13 @@ struct sigmaHadCorrTask { { // Sigma- -> n + pi- (charged daughter = pion, neutral daughter = neutron) // Sigma+ -> p + pi0 (charged daughter = proton, neutral daughter = pi0) + const float epsilon = 1e-6f; float massChargedDau = doSigmaMinus ? o2::constants::physics::MassPionCharged : o2::constants::physics::MassProton; float massNeutralDau = doSigmaMinus ? o2::constants::physics::MassNeutron : o2::constants::physics::MassPionNeutral; float massSigma = doSigmaMinus ? o2::constants::physics::MassSigmaMinus : o2::constants::physics::MassSigmaPlus; float pMother = std::sqrt(sigmaPx * sigmaPx + sigmaPy * sigmaPy + sigmaPz * sigmaPz); - if (pMother < 1e-6f) { + if (pMother < epsilon) { return -999.f; } float versorX = sigmaPx / pMother; @@ -239,7 +241,7 @@ struct sigmaHadCorrTask { float A = 4.f * (eChDau * eChDau - a * a); float B = -4.f * a * K; float C = 4.f * eChDau * eChDau * massSigma * massSigma - K * K; - if (std::abs(A) < 1e-6f) { + if (std::abs(A) < epsilon) { return -999.f; } float D = B * B - 4.f * A * C; @@ -309,8 +311,9 @@ struct sigmaHadCorrTask { template int findFirstPartonicMotherPDG(const TMcParticle& mcParticle, const TMcParticles& mcParticles) { + const int notFoundPdg = -999; if (!mcParticle.has_mothers()) { - return -999; + return notFoundPdg; } auto motherIds = mcParticle.mothersIds(); const int defaultMotherSize = 2; @@ -331,10 +334,10 @@ struct sigmaHadCorrTask { return mother.pdgCode(); } int found = findFirstPartonicMotherPDG(mother, mcParticles); - if (found != -999) + if (found != notFoundPdg) return found; } - return -999; + return notFoundPdg; } // Walk up the decay chain iteratively and return the PDG of the last quark or gluon before beam remnants @@ -403,24 +406,19 @@ struct sigmaHadCorrTask { return doSigmaPion ? track.tofNSigmaPi() : track.tofNSigmaPr(); } - TLorentzVector trackSum, PartOneCMS, PartTwoCMS, trackRelK; float getKStar(float sigmaPx, float sigmaPy, float sigmaPz, float pxHad, float pyHad, float pzHad) { - TLorentzVector part1; // Sigma - TLorentzVector part2; // Hadron track (proton/pion) - part1.SetXYZM(sigmaPx, sigmaPy, sigmaPz, getSigmaMassForKstar()); - part2.SetXYZM(pxHad, pyHad, pzHad, getHadTrackMass()); - trackSum = part1 + part2; + ROOT::Math::PxPyPzMVector part1(sigmaPx, sigmaPy, sigmaPz, getSigmaMassForKstar()); // Sigma + ROOT::Math::PxPyPzMVector part2(pxHad, pyHad, pzHad, getHadTrackMass()); // Hadron track (proton/pion) + ROOT::Math::PxPyPzMVector trackSum = part1 + part2; const float beta = trackSum.Beta(); const float betax = beta * std::cos(trackSum.Phi()) * std::sin(trackSum.Theta()); const float betay = beta * std::sin(trackSum.Phi()) * std::sin(trackSum.Theta()); const float betaz = beta * std::cos(trackSum.Theta()); - PartOneCMS.SetXYZM(part1.Px(), part1.Py(), part1.Pz(), part1.M()); - PartTwoCMS.SetXYZM(part2.Px(), part2.Py(), part2.Pz(), part2.M()); const ROOT::Math::Boost boostPRF = ROOT::Math::Boost(-betax, -betay, -betaz); - PartOneCMS = boostPRF(PartOneCMS); - PartTwoCMS = boostPRF(PartTwoCMS); - trackRelK = PartOneCMS - PartTwoCMS; + ROOT::Math::PxPyPzMVector partOneCMS = boostPRF(part1); + ROOT::Math::PxPyPzMVector partTwoCMS = boostPRF(part2); + ROOT::Math::PxPyPzMVector trackRelK = partOneCMS - partTwoCMS; return 0.5 * trackRelK.P(); } @@ -430,7 +428,7 @@ struct sigmaHadCorrTask { if (candidate.pt() < ptMinHad) { return false; } - if (std::abs(getTPCNSigmaHad(candidate)) > cutNSigmaTPC || candidate.tpcNClsFound() < cutNTPCClusHad || std::abs(candidate.eta()) > cutEtaDaught) { + if (std::abs(getTPCNSigmaHad(candidate)) > cutNSigmaTPC || candidate.tpcNClsFound() < cutNTPCClusHad || std::abs(candidate.eta()) > cutEtaDaughter) { return false; } @@ -500,9 +498,10 @@ struct sigmaHadCorrTask { return true; } - template - void fillTreeAndHistograms(aod::KinkCands const& kinkCands, Ttrack const& tracksDauSigma, Ttrack const& tracks, Tcollision const& collision, bool isMC) + template + std::vector fillTreeAndHistograms(aod::KinkCands const& kinkCands, Ttrack const& tracksDauSigma, Ttrack const& tracks, Tcollision const& collision) { + std::vector sigmaHadCandidates; for (const auto& sigmaCand : kinkCands) { auto kinkDauTrack = tracksDauSigma.rawIteratorAt(sigmaCand.trackDaugId()); if (!selectSigma(sigmaCand, kinkDauTrack)) { @@ -528,6 +527,30 @@ struct sigmaHadCorrTask { rSigmaHad.fill(HIST("QA/hSigmaPtRecal"), sigmaPtRecal); rSigmaHad.fill(HIST("QA/h2InvMassVsPtSigma"), sigmaPtRecal, sigmaMassForQa); + // single-Sigma-level MC truth record, filled once per accepted candidate, independent of hadron pairing + if constexpr (IsMC) { + auto mothTrack = tracksDauSigma.rawIteratorAt(sigmaCand.trackMothId()); + if (mothTrack.has_mcParticle() && kinkDauTrack.has_mcParticle()) { + auto mcMoth = mothTrack.template mcParticle_as(); + auto mcDaug = kinkDauTrack.template mcParticle_as(); + float massMC = std::sqrt(mcMoth.e() * mcMoth.e() - mcMoth.p() * mcMoth.p()); + float decayRadiusMC = std::hypot(mcDaug.vx() - mcMoth.vx(), mcDaug.vy() - mcMoth.vy()); + bool collisionIdCheck = false; + if (collision.has_mcCollision()) { + collisionIdCheck = collision.mcCollision().globalIndex() == mcDaug.mcCollisionId(); + } + outputKinkCandsMC(sigmaCand.xDecVtx(), sigmaCand.yDecVtx(), sigmaCand.zDecVtx(), + sigmaCand.pxMoth(), sigmaCand.pyMoth(), sigmaCand.pzMoth(), + sigmaCand.pxDaug(), sigmaCand.pyDaug(), sigmaCand.pzDaug(), + sigmaCand.dcaMothPv(), sigmaCand.dcaDaugPv(), sigmaCand.dcaKinkTopo(), + sigmaCand.mothSign(), + kinkDauTrack.tpcNSigmaPi(), kinkDauTrack.tpcNSigmaPr(), -999.f, + kinkDauTrack.tofNSigmaPi(), kinkDauTrack.tofNSigmaPr(), -999.f, + mcMoth.pdgCode(), mcDaug.pdgCode(), + mcMoth.pt(), mcMoth.pz(), massMC, decayRadiusMC, collisionIdCheck); + } + } + for (const auto& hadTrack : tracks) { if (hadTrack.globalIndex() == sigmaCand.trackDaugId()) { continue; @@ -537,7 +560,7 @@ struct sigmaHadCorrTask { continue; } - sigmaHadCand candidate; + SigmaHadCand candidate; candidate.sigmaCharge = sigmaCand.mothSign(); candidate.sigmaPx = sigmaCand.pxMoth(); candidate.sigmaPy = sigmaCand.pyMoth(); @@ -575,33 +598,35 @@ struct sigmaHadCorrTask { if (hadTrack.hasTOF()) { rSigmaHad.fill(HIST("QA/h2TOFNSigmaHadVsPtHad"), candidate.ptHad(), candidate.nSigmaTOFHad); } - if (fillSparseInvMassKstar && !isMC) { - rSigmaHad.fill(HIST("hSparseSigmaHad"), - candidate.sigmaMass, - kStar, - candidate.sigmaCharge, - candidate.chargeHad, - candidate.sigmaDecRadius, - candidate.sigmaCosPA, - sigmaPtRecal); + if constexpr (!IsMC) { + if (fillSparseInvMassKstar) { + rSigmaHad.fill(HIST("hSparseSigmaHad"), + candidate.sigmaMass, + kStar, + candidate.sigmaCharge, + candidate.chargeHad, + candidate.sigmaDecRadius, + candidate.sigmaCosPA, + sigmaPtRecal); + } } sigmaHadCandidates.push_back(candidate); } } + return sigmaHadCandidates; } void processSameEvent(CollisionsFull const& collisions, aod::KinkCands const& kinkCands, TracksFull const& tracks) { for (auto const& collision : collisions) { - sigmaHadCandidates.clear(); - auto kinkCands_c = kinkCands.sliceBy(kinkCandsPerCollisionPreslice, collision.globalIndex()); - auto tracks_c = tracks.sliceBy(tracksPerCollisionPreslice, collision.globalIndex()); - if (std::abs(collision.posZ()) > cutzvertex || !collision.sel8()) { + auto kinkCandsC = kinkCands.sliceBy(kinkCandsPerCollisionPreslice, collision.globalIndex()); + auto tracksC = tracks.sliceBy(tracksPerCollisionPreslice, collision.globalIndex()); + if (std::abs(collision.posZ()) > cutZVertex || !collision.sel8()) { continue; } rEventSelection.fill(HIST("hVertexZRec"), collision.posZ()); - fillTreeAndHistograms(kinkCands_c, tracks, tracks_c, collision, false); + auto sigmaHadCandidates = fillTreeAndHistograms(kinkCandsC, tracks, tracksC, collision); if (fillOutputTree) { // Fill output table for (const auto& candidate : sigmaHadCandidates) { @@ -625,7 +650,7 @@ struct sigmaHadCorrTask { } } } - PROCESS_SWITCH(sigmaHadCorrTask, processSameEvent, "Process Same event", true); + PROCESS_SWITCH(SigmaHadCorr, processSameEvent, "Process Same event", true); // Processing Event Mixing SliceCache cache; @@ -636,17 +661,16 @@ struct sigmaHadCorrTask { { if (useMultNTracksPV.value) { for (auto const& [collision1, collision2] : - selfCombinations(BinningTypeMultNTracksPV{{CfgVtxBins, CfgMultBins}, true}, nEvtMixingBkg, -1, collisions, collisions)) { + selfCombinations(BinningTypeMultNTracksPV{{cfgVtxBins, cfgMultBins}, true}, nEvtMixingBkg, -1, collisions, collisions)) { if (collision1.index() == collision2.index()) continue; - sigmaHadCandidates.clear(); - if (std::abs(collision1.posZ()) > cutzvertex || !collision1.sel8()) + if (std::abs(collision1.posZ()) > cutZVertex || !collision1.sel8()) continue; - if (std::abs(collision2.posZ()) > cutzvertex || !collision2.sel8()) + 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); + auto kinkCandsC1 = kinkCands.sliceBy(kinkCandsPerCollisionPreslice, collision1.globalIndex()); + auto tracksC2 = tracks.sliceBy(tracksPerCollisionPreslice, collision2.globalIndex()); + auto sigmaHadCandidates = fillTreeAndHistograms(kinkCandsC1, tracks, tracksC2, collision1); if (fillOutputTree) { for (const auto& candidate : sigmaHadCandidates) { outputDataTable(candidate.sigmaCharge, candidate.sigmaPx, candidate.sigmaPy, candidate.sigmaPz, @@ -659,17 +683,16 @@ struct sigmaHadCorrTask { } } else { for (auto const& [collision1, collision2] : - selfCombinations(BinningTypeNumContrib{{CfgVtxBins, CfgMultBins}, true}, nEvtMixingBkg, -1, collisions, collisions)) { + selfCombinations(BinningTypeNumContrib{{cfgVtxBins, cfgMultBins}, true}, nEvtMixingBkg, -1, collisions, collisions)) { if (collision1.index() == collision2.index()) continue; - sigmaHadCandidates.clear(); - if (std::abs(collision1.posZ()) > cutzvertex || !collision1.sel8()) + if (std::abs(collision1.posZ()) > cutZVertex || !collision1.sel8()) continue; - if (std::abs(collision2.posZ()) > cutzvertex || !collision2.sel8()) + 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); + auto kinkCandsC1 = kinkCands.sliceBy(kinkCandsPerCollisionPreslice, collision1.globalIndex()); + auto tracksC2 = tracks.sliceBy(tracksPerCollisionPreslice, collision2.globalIndex()); + auto sigmaHadCandidates = fillTreeAndHistograms(kinkCandsC1, tracks, tracksC2, collision1); if (fillOutputTree) { for (const auto& candidate : sigmaHadCandidates) { outputDataTable(candidate.sigmaCharge, candidate.sigmaPx, candidate.sigmaPy, candidate.sigmaPz, @@ -682,21 +705,20 @@ struct sigmaHadCorrTask { } } } - PROCESS_SWITCH(sigmaHadCorrTask, processMixedEvent, "Process Mixed event", false); + PROCESS_SWITCH(SigmaHadCorr, processMixedEvent, "Process Mixed event", false); - void processSameEventMC(CollisionsFullMC const& collisions, aod::KinkCands const& kinkCands, TracksFullMC const& tracks, aod::McParticles const& mcParticles) + void processSameEventMC(CollisionsFullMC const& collisions, aod::KinkCands const& kinkCands, TracksFullMC const& tracks, aod::McParticles const& mcParticles, aod::McCollisions const&) { for (auto const& collision : collisions) { - sigmaHadCandidates.clear(); - auto kinkCands_c = kinkCands.sliceBy(kinkCandsPerCollisionPreslice, collision.globalIndex()); - auto tracks_c = tracks.sliceBy(tracksMCPerCollisionPreslice, collision.globalIndex()); + auto kinkCandsC = kinkCands.sliceBy(kinkCandsPerCollisionPreslice, collision.globalIndex()); + auto tracksC = tracks.sliceBy(tracksMCPerCollisionPreslice, collision.globalIndex()); - if (std::abs(collision.posZ()) > cutzvertex || !collision.sel8()) { + if (std::abs(collision.posZ()) > cutZVertex || !collision.sel8()) { continue; } rEventSelection.fill(HIST("hVertexZRec"), collision.posZ()); - fillTreeAndHistograms(kinkCands_c, tracks, tracks_c, collision, true); + auto sigmaHadCandidates = fillTreeAndHistograms(kinkCandsC, tracks, tracksC, collision); for (const auto& candidate : sigmaHadCandidates) { auto mcLabelSigma = tracks.rawIteratorAt(candidate.sigmaID); auto mcLabelSigmaDau = tracks.rawIteratorAt(candidate.kinkDauID); @@ -763,24 +785,58 @@ struct sigmaHadCorrTask { } } } + + // all generated Sigma -> chargedDau + neutralDau decays + int pdgChargedDauAbs = doSigmaMinus ? PDG_t::kPiPlus : PDG_t::kProton; + for (const auto& mcPart : mcParticles) { + int pdgMothAbs = std::abs(mcPart.pdgCode()); + bool isValidMother = doSigmaMinus ? (pdgMothAbs == PDG_t::kSigmaMinus || pdgMothAbs == PDG_t::kSigmaPlus) : (pdgMothAbs == PDG_t::kSigmaPlus); + if (!isValidMother) { + continue; + } + bool hasChargedDaughter = false; + std::array genDecVtx{-999.f, -999.f, -999.f}; + int daugPdgCode = 0; + for (const auto& daughter : mcPart.daughters_as()) { + if (std::abs(daughter.pdgCode()) == pdgChargedDauAbs) { + hasChargedDaughter = true; + genDecVtx = {daughter.vx(), daughter.vy(), daughter.vz()}; + daugPdgCode = daughter.pdgCode(); + break; + } + } + if (!hasChargedDaughter) { + continue; + } + float massMC = std::sqrt(mcPart.e() * mcPart.e() - mcPart.p() * mcPart.p()); + float decayRadiusMC = std::hypot(genDecVtx[0] - mcPart.vx(), genDecVtx[1] - mcPart.vy()); + outputKinkCandsMC(-999.f, -999.f, -999.f, + -999.f, -999.f, -999.f, + -999.f, -999.f, -999.f, + -999.f, -999.f, -999.f, + mcPart.pdgCode() > 0 ? 1 : -1, + -999.f, -999.f, -999.f, + -999.f, -999.f, -999.f, + mcPart.pdgCode(), daugPdgCode, + mcPart.pt(), mcPart.pz(), massMC, decayRadiusMC, false); + } } - PROCESS_SWITCH(sigmaHadCorrTask, processSameEventMC, "Process Same event MC", false); + PROCESS_SWITCH(SigmaHadCorr, processSameEventMC, "Process Same event MC", false); - void processMixedEventMC(const CollisionsFullMC& collisions, const aod::KinkCands& kinkCands, const TracksFullMC& tracks, const aod::McParticles& mcParticles) + void processMixedEventMC(const CollisionsFullMC& collisions, const aod::KinkCands& kinkCands, const TracksFullMC& tracks, const aod::McParticles& mcParticles, aod::McCollisions const&) { if (useMultNTracksPV.value) { for (auto const& [collision1, collision2] : - selfCombinations(BinningTypeMultNTracksPV{{CfgVtxBins, CfgMultBins}, true}, nEvtMixingBkg, -1, collisions, collisions)) { + selfCombinations(BinningTypeMultNTracksPV{{cfgVtxBins, cfgMultBins}, true}, nEvtMixingBkg, -1, collisions, collisions)) { if (collision1.index() == collision2.index()) continue; - sigmaHadCandidates.clear(); - if (std::abs(collision1.posZ()) > cutzvertex || !collision1.sel8()) + if (std::abs(collision1.posZ()) > cutZVertex || !collision1.sel8()) continue; - if (std::abs(collision2.posZ()) > cutzvertex || !collision2.sel8()) + 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); + auto kinkCandsC1 = kinkCands.sliceBy(kinkCandsPerCollisionPreslice, collision1.globalIndex()); + auto tracksC2 = tracks.sliceBy(tracksPerCollisionPreslice, collision2.globalIndex()); + auto sigmaHadCandidates = fillTreeAndHistograms(kinkCandsC1, tracks, tracksC2, collision1); for (const auto& candidate : sigmaHadCandidates) { auto mcLabelSigma = tracks.rawIteratorAt(candidate.sigmaID); auto mcLabelSigmaDau = tracks.rawIteratorAt(candidate.kinkDauID); @@ -819,17 +875,16 @@ struct sigmaHadCorrTask { } } else { for (auto const& [collision1, collision2] : - selfCombinations(BinningTypeNumContrib{{CfgVtxBins, CfgMultBins}, true}, nEvtMixingBkg, -1, collisions, collisions)) { + selfCombinations(BinningTypeNumContrib{{cfgVtxBins, cfgMultBins}, true}, nEvtMixingBkg, -1, collisions, collisions)) { if (collision1.index() == collision2.index()) continue; - sigmaHadCandidates.clear(); - if (std::abs(collision1.posZ()) > cutzvertex || !collision1.sel8()) + if (std::abs(collision1.posZ()) > cutZVertex || !collision1.sel8()) continue; - if (std::abs(collision2.posZ()) > cutzvertex || !collision2.sel8()) + 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); + auto kinkCandsC1 = kinkCands.sliceBy(kinkCandsPerCollisionPreslice, collision1.globalIndex()); + auto tracksC2 = tracks.sliceBy(tracksPerCollisionPreslice, collision2.globalIndex()); + auto sigmaHadCandidates = fillTreeAndHistograms(kinkCandsC1, tracks, tracksC2, collision1); for (const auto& candidate : sigmaHadCandidates) { auto mcLabelSigma = tracks.rawIteratorAt(candidate.sigmaID); auto mcLabelSigmaDau = tracks.rawIteratorAt(candidate.kinkDauID); @@ -868,10 +923,10 @@ struct sigmaHadCorrTask { } } } - PROCESS_SWITCH(sigmaHadCorrTask, processMixedEventMC, "Process Mixed event MC", false); + PROCESS_SWITCH(SigmaHadCorr, processMixedEventMC, "Process Mixed event MC", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ - adaptAnalysisTask(cfgc)}; + adaptAnalysisTask(cfgc)}; } diff --git a/PWGLF/TableProducer/Strangeness/sigmaplusbuilder.cxx b/PWGLF/TableProducer/Strangeness/sigmaplusbuilder.cxx index 739e7859c0f..96f70efd7e2 100644 --- a/PWGLF/TableProducer/Strangeness/sigmaplusbuilder.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaplusbuilder.cxx @@ -28,9 +28,13 @@ #include #include #include +#include #include #include #include +#include +#include +#include #include #include #include @@ -69,7 +73,6 @@ struct Sigmaplusbuilder { Configurable photonMaxDCAV0Dau{"photonMaxDCAV0Dau", 3.5, "Max DCA between photon daughters (cm)"}; Configurable photonMaxQt{"photonMaxQt", 0.15, "Max Armenteros qT for photons (GeV/c)"}; Configurable photonMaxAlpha{"photonMaxAlpha", 1.0, "Max |Armenteros alpha| for photons"}; - Configurable photonMaxTPCNSigmaEl{"photonMaxTPCNSigmaEl", 15, "Max |TPC nSigma_el| for photon daughters"}; // proton selection Configurable protonMinPt{"protonMinPt", 0.3, "Minimum proton pT (GeV/c)"}; @@ -99,7 +102,7 @@ struct Sigmaplusbuilder { Produces sigmaPlusCands; Produces sigmaPlusCandsMC; - Service ccdb; + Service ccdb{}; o2::vertexing::DCAFitterN<2> fitter; int mRunNumber = 0; float mBz = 0; @@ -123,6 +126,7 @@ struct Sigmaplusbuilder { fitter.setMinParamChange(1e-3); fitter.setMinRelChi2Change(0.9); fitter.setMaxDZIni(1e9); + fitter.setMaxDXYIni(1e9); fitter.setMaxChi2(1e9); fitter.setUseAbsDCA(true); @@ -150,16 +154,16 @@ struct Sigmaplusbuilder { histos.add("Photon/hSelectionCounter", "Photon/hSelectionCounter", kTH1F, {axisPhotonSel}); auto hPhotonSel = histos.get(HIST("Photon/hSelectionCounter")); hPhotonSel->GetXaxis()->SetBinLabel(1, "All"); - hPhotonSel->GetXaxis()->SetBinLabel(2, "Mass"); - hPhotonSel->GetXaxis()->SetBinLabel(3, "Rapidity"); - hPhotonSel->GetXaxis()->SetBinLabel(4, "Neg eta"); - hPhotonSel->GetXaxis()->SetBinLabel(5, "Pos eta"); - hPhotonSel->GetXaxis()->SetBinLabel(6, "DCA daughters"); - hPhotonSel->GetXaxis()->SetBinLabel(7, "Radius"); - hPhotonSel->GetXaxis()->SetBinLabel(8, "CosPA"); + hPhotonSel->GetXaxis()->SetBinLabel(2, "Neg eta"); + hPhotonSel->GetXaxis()->SetBinLabel(3, "Pos eta"); + hPhotonSel->GetXaxis()->SetBinLabel(4, "Vertex fit"); + hPhotonSel->GetXaxis()->SetBinLabel(5, "DCA daughters"); + hPhotonSel->GetXaxis()->SetBinLabel(6, "Radius"); + hPhotonSel->GetXaxis()->SetBinLabel(7, "CosPA"); + hPhotonSel->GetXaxis()->SetBinLabel(8, "Rapidity"); hPhotonSel->GetXaxis()->SetBinLabel(9, "Qt"); hPhotonSel->GetXaxis()->SetBinLabel(10, "Alpha"); - hPhotonSel->GetXaxis()->SetBinLabel(11, "TPC nSigma el"); + hPhotonSel->GetXaxis()->SetBinLabel(11, "Mass"); histos.add("Photon/hMass", "Photon/hMass", kTH1F, {axisPhotonMass}); histos.add("Photon/hPt", "Photon/hPt", kTH1F, {axisPhotonPt}); @@ -210,16 +214,16 @@ struct Sigmaplusbuilder { histos.add("Photon/True/hSelectionCounter", "Photon/True/hSelectionCounter", kTH1F, {axisPhotonSel}); auto hPhotonSelSignal = histos.get(HIST("Photon/True/hSelectionCounter")); hPhotonSelSignal->GetXaxis()->SetBinLabel(1, "All"); - hPhotonSelSignal->GetXaxis()->SetBinLabel(2, "Mass"); - hPhotonSelSignal->GetXaxis()->SetBinLabel(3, "Rapidity"); - hPhotonSelSignal->GetXaxis()->SetBinLabel(4, "Neg eta"); - hPhotonSelSignal->GetXaxis()->SetBinLabel(5, "Pos eta"); - hPhotonSelSignal->GetXaxis()->SetBinLabel(6, "DCA daughters"); - hPhotonSelSignal->GetXaxis()->SetBinLabel(7, "Radius"); - hPhotonSelSignal->GetXaxis()->SetBinLabel(8, "CosPA"); + hPhotonSelSignal->GetXaxis()->SetBinLabel(2, "Neg eta"); + hPhotonSelSignal->GetXaxis()->SetBinLabel(3, "Pos eta"); + hPhotonSelSignal->GetXaxis()->SetBinLabel(4, "Vertex fit"); + hPhotonSelSignal->GetXaxis()->SetBinLabel(5, "DCA daughters"); + hPhotonSelSignal->GetXaxis()->SetBinLabel(6, "Radius"); + hPhotonSelSignal->GetXaxis()->SetBinLabel(7, "CosPA"); + hPhotonSelSignal->GetXaxis()->SetBinLabel(8, "Rapidity"); hPhotonSelSignal->GetXaxis()->SetBinLabel(9, "Qt"); hPhotonSelSignal->GetXaxis()->SetBinLabel(10, "Alpha"); - hPhotonSelSignal->GetXaxis()->SetBinLabel(11, "TPC nSigma el"); + hPhotonSelSignal->GetXaxis()->SetBinLabel(11, "Mass"); histos.add("Photon/True/hMass", "Photon/True/hMass", kTH1F, {axisPhotonMass}); histos.add("Photon/True/hPt", "Photon/True/hPt", kTH1F, {axisPhotonPt}); @@ -291,6 +295,7 @@ struct Sigmaplusbuilder { histos.add("Findable/hPhotonMomentumResolution", "Findable/hPhotonMomentumResolution", kTH1F, {axisPhotonMomResolution}); histos.add("Findable/hPhotonCosPA", "Findable/hPhotonCosPA", kTH1F, {axisPhotonCosPA}); histos.add("Findable/hConversionPairV0Presence", "Findable/hConversionPairV0Presence", kTH1F, {axisV0Presence}); + histos.add("Findable/hPhotonSearchPresence", "Findable/hPhotonSearchPresence", kTH1F, {axisV0Presence}); histos.add("Findable/hDuplicateConversionTrackCounter", "Findable/hDuplicateConversionTrackCounter", kTH1F, {axisDuplicateTrack}); histos.add("Findable/hDuplicateElectronTPCNClsFound", "Findable/hDuplicateElectronTPCNClsFound", kTH1F, {axisTPCClusters}); histos.add("Findable/hDuplicatePositronTPCNClsFound", "Findable/hDuplicatePositronTPCNClsFound", kTH1F, {axisTPCClusters}); @@ -299,6 +304,9 @@ struct Sigmaplusbuilder { auto hConversionPairV0Presence = histos.get(HIST("Findable/hConversionPairV0Presence")); hConversionPairV0Presence->GetXaxis()->SetBinLabel(1, "valid pair"); hConversionPairV0Presence->GetXaxis()->SetBinLabel(2, "in V0"); + auto hPhotonSearchPresence = histos.get(HIST("Findable/hPhotonSearchPresence")); + hPhotonSearchPresence->GetXaxis()->SetBinLabel(1, "valid pair"); + hPhotonSearchPresence->GetXaxis()->SetBinLabel(2, "found by findPhotons"); auto hDuplicateConversionTrackCounter = histos.get(HIST("Findable/hDuplicateConversionTrackCounter")); hDuplicateConversionTrackCounter->GetXaxis()->SetBinLabel(1, "e^{-}"); hDuplicateConversionTrackCounter->GetXaxis()->SetBinLabel(2, "e^{+}"); @@ -313,90 +321,145 @@ struct Sigmaplusbuilder { } } - // photon (PCM) candidate selection - template - bool selectPhoton(const TV0& v0) + // photon (electron/positron conversion pair) candidate + template + struct PhotonCand { + float x = 0.f, y = 0.f, z = 0.f; + float px = 0.f, py = 0.f, pz = 0.f; + float mGamma = 0.f; + float alpha = 0.f; + float qtarm = 0.f; + float radius = 0.f; + TTrack negTrack; + TTrack posTrack; + }; + + // photon (PCM) candidate search: pair up all neg/pos tracks in the collision and fit a common vertex ourselves + template + std::vector> findPhotons(const TTracks& tracks, const std::array& pv) { - auto posTrack = v0.template posTrack_as(); - auto negTrack = v0.template negTrack_as(); + std::vector> photons; - bool isSignal = false; - if constexpr (IsMC) { - if (posTrack.has_mcParticle() && negTrack.has_mcParticle()) { - auto mcPos = posTrack.template mcParticle_as(); - auto mcNeg = negTrack.template mcParticle_as(); - isSignal = findSigmaPlusMotherOfPhoton(mcPos, mcNeg) >= 0; + for (const auto& negTrack : tracks) { + if (negTrack.sign() > 0) { + continue; } - } + for (const auto& posTrack : tracks) { + if (posTrack.sign() < 0) { + continue; + } - auto fillPhotonStep = [&](int step) { - histos.fill(HIST("Photon/hSelectionCounter"), step); - if constexpr (IsMC) { - if (isSignal) { - histos.fill(HIST("Photon/True/hSelectionCounter"), step); + bool isSignal = false; + if constexpr (IsMC) { + if (posTrack.has_mcParticle() && negTrack.has_mcParticle()) { + auto mcPos = posTrack.template mcParticle_as(); + auto mcNeg = negTrack.template mcParticle_as(); + isSignal = findSigmaPlusMotherOfPhoton(mcPos, mcNeg) >= 0; + } } - } - }; - fillPhotonStep(0); - if (v0.mGamma() < 0 || v0.mGamma() > photonMaxMass) { - return false; - } - fillPhotonStep(1); + auto fillPhotonStep = [&](int step) { + histos.fill(HIST("Photon/hSelectionCounter"), step); + if constexpr (IsMC) { + if (isSignal) { + histos.fill(HIST("Photon/True/hSelectionCounter"), step); + } + } + }; + fillPhotonStep(0); - float photonY = RecoDecay::y(std::array{v0.px(), v0.py(), v0.pz()}, o2::constants::physics::MassGamma); - if (photonY < photonMinRapidity || photonY > photonMaxRapidity) { - return false; - } - fillPhotonStep(2); + if (negTrack.eta() < photonDauEtaMin || negTrack.eta() > photonDauEtaMax) { + continue; + } + fillPhotonStep(1); - if (v0.negativeeta() < photonDauEtaMin || v0.negativeeta() > photonDauEtaMax) { - return false; - } - fillPhotonStep(3); + if (posTrack.eta() < photonDauEtaMin || posTrack.eta() > photonDauEtaMax) { + continue; + } + fillPhotonStep(2); + + auto negTrackParCov = getTrackParCov(negTrack); + auto posTrackParCov = getTrackParCov(posTrack); + int nCand = 0; + try { + nCand = fitter.process(negTrackParCov, posTrackParCov); + } catch (...) { + continue; + } + if (nCand == 0 || !fitter.propagateTracksToVertex()) { + continue; + } + fillPhotonStep(3); - if (v0.positiveeta() < photonDauEtaMin || v0.positiveeta() > photonDauEtaMax) { - return false; - } - fillPhotonStep(4); + float dcaDaughters = std::sqrt(fitter.getChi2AtPCACandidate()); + if (dcaDaughters > photonMaxDCAV0Dau) { + continue; + } + fillPhotonStep(4); - if (std::abs(v0.dcaV0daughters()) > photonMaxDCAV0Dau) { - return false; - } - fillPhotonStep(5); + std::array secVtx = fitter.getPCACandidatePos(); + float radius = std::hypot(secVtx[0], secVtx[1]); + if (radius < photonMinRadius || radius > photonMaxRadius) { + continue; + } + fillPhotonStep(5); + + std::array pNeg{}; + std::array pPos{}; + fitter.getTrack(0).getPxPyPzGlo(pNeg); + fitter.getTrack(1).getPxPyPzGlo(pPos); + std::array pGamma{pNeg[0] + pPos[0], pNeg[1] + pPos[1], pNeg[2] + pPos[2]}; + float gammaP = std::sqrt(dot3(pGamma, pGamma)); + + std::array flightVec{secVtx[0] - pv[0], secVtx[1] - pv[1], secVtx[2] - pv[2]}; + float flightNorm = std::sqrt(dot3(flightVec, flightVec)); + float cosPA = dot3(flightVec, pGamma) / (flightNorm * gammaP); + if (cosPA < photonMinV0cospa) { + continue; + } + fillPhotonStep(6); - if (v0.v0radius() < photonMinRadius || v0.v0radius() > photonMaxRadius) { - return false; - } - fillPhotonStep(6); + float photonY = RecoDecay::y(pGamma, o2::constants::physics::MassGamma); + if (photonY < photonMinRapidity || photonY > photonMaxRapidity) { + continue; + } + fillPhotonStep(7); - if (v0.v0cosPA() < photonMinV0cospa) { - return false; - } - fillPhotonStep(7); + float qLNeg = dot3(pNeg, pGamma) / gammaP; + float qLPos = dot3(pPos, pGamma) / gammaP; + float alpha = (qLPos - qLNeg) / (qLPos + qLNeg); + float qtarm = std::sqrt(std::max(dot3(pNeg, pNeg) - qLNeg * qLNeg, 0.f)); + if (qtarm > photonMaxQt) { + continue; + } + fillPhotonStep(8); - if (v0.qtarm() > photonMaxQt) { - return false; - } - fillPhotonStep(8); + if (std::abs(alpha) > photonMaxAlpha) { + continue; + } + fillPhotonStep(9); - if (std::abs(v0.alpha()) > photonMaxAlpha) { - return false; - } - fillPhotonStep(9); + constexpr float MassEl = o2::constants::physics::MassElectron; + float eNeg = std::sqrt(dot3(pNeg, pNeg) + MassEl * MassEl); + float ePos = std::sqrt(dot3(pPos, pPos) + MassEl * MassEl); + float mGamma2 = (eNeg + ePos) * (eNeg + ePos) - gammaP * gammaP; + if (mGamma2 < 0.f || std::sqrt(mGamma2) > photonMaxMass) { + continue; + } + float mGamma = std::sqrt(mGamma2); + fillPhotonStep(10); - if (std::abs(posTrack.tpcNSigmaEl()) > photonMaxTPCNSigmaEl || std::abs(negTrack.tpcNSigmaEl()) > photonMaxTPCNSigmaEl) { - return false; - } - fillPhotonStep(10); + histos.fill(HIST("Photon/hMass"), mGamma); + histos.fill(HIST("Photon/hPt"), std::hypot(pGamma[0], pGamma[1])); + histos.fill(HIST("Photon/hRadius"), radius); + histos.fill(HIST("Photon/h2ArmenterosPodolanski"), alpha, qtarm); + histos.fill(HIST("Photon/h2ConvPointXY"), secVtx[0], secVtx[1]); - histos.fill(HIST("Photon/hMass"), v0.mGamma()); - histos.fill(HIST("Photon/hPt"), v0.pt()); - histos.fill(HIST("Photon/hRadius"), v0.v0radius()); - histos.fill(HIST("Photon/h2ArmenterosPodolanski"), v0.alpha(), v0.qtarm()); - histos.fill(HIST("Photon/h2ConvPointXY"), v0.x(), v0.y()); + photons.push_back({secVtx[0], secVtx[1], secVtx[2], pGamma[0], pGamma[1], pGamma[2], mGamma, alpha, qtarm, radius, negTrack, posTrack}); + } + } - return true; + return photons; } // proton candidate selection @@ -562,16 +625,19 @@ struct Sigmaplusbuilder { } // Build a Sigma+ -> p pi0 candidate from a proton track and a PCM photon - template - void buildSigmaPlusCandidate(const TTrack& protonTrack, const TV0& photon, const std::array& pv) + // Returns the MC index of the matched true Sigma+ mother if a signal candidate was built, -1 otherwise + template + int buildSigmaPlusCandidate(const TTrack& protonTrack, const PhotonCand& photon, const std::array& pv) { - auto posTrack = photon.template posTrack_as(); - auto negTrack = photon.template negTrack_as(); + auto posTrack = photon.posTrack; + auto negTrack = photon.negTrack; bool isSignal = false; - std::array mcTrueVtx{}; // Sigma+ decay vertex - std::array mcTrueMomProton{}; // true MC proton momentum - std::array mcTrueMomGamma{}; // true MC momentum of the measured photon + std::array mcTrueVtx{}; // Sigma+ decay vertex + std::array mcTrueMomProton{}; // true MC proton momentum + std::array mcTrueMomGamma{}; // true MC momentum of the measured photon + std::array mcTrueMomSigmaPlus{}; // true MC momentum of the Sigma+ mother + int matchedSigmaId = -1; int protonPdgCode = 0; int protonMotherPdgCode = 0; int gammaPdgCode = 0; @@ -617,6 +683,8 @@ struct Sigmaplusbuilder { if (isSignal) { mcTrueVtx = {mcProton.vx(), mcProton.vy(), mcProton.vz()}; mcTrueMomProton = {mcProton.px(), mcProton.py(), mcProton.pz()}; + mcTrueMomSigmaPlus = {protonMothers.front().px(), protonMothers.front().py(), protonMothers.front().pz()}; + matchedSigmaId = protonMothers.front().globalIndex(); } } } @@ -634,7 +702,7 @@ struct Sigmaplusbuilder { auto protonTrackParCov = getTrackParCov(protonTrack); std::array zeroCov{}; - auto photonTrackParCov = o2::track::TrackParCov({photon.x(), photon.y(), photon.z()}, {photon.px(), photon.py(), photon.pz()}, zeroCov, 0, true); + auto photonTrackParCov = o2::track::TrackParCov({photon.x, photon.y, photon.z}, {photon.px, photon.py, photon.pz}, zeroCov, 0, true); photonTrackParCov.setAbsCharge(0); photonTrackParCov.setPID(o2::track::PID::Photon); @@ -642,10 +710,10 @@ struct Sigmaplusbuilder { try { nCand = fitter.process(protonTrackParCov, photonTrackParCov); } catch (...) { - return; + return -1; } if (nCand == 0 || !fitter.propagateTracksToVertex()) { - return; + return -1; } fillCandStep(1); // Vertex fit @@ -658,7 +726,7 @@ struct Sigmaplusbuilder { } } if (dcaProtonGamma > candMaxDcaProtonGamma) { - return; + return -1; } fillCandStep(2); // DCA(p,gamma) @@ -673,7 +741,7 @@ struct Sigmaplusbuilder { } } if (radius < candMinRadius || radius > candMaxRadius) { - return; + return -1; } fillCandStep(3); // radius @@ -682,8 +750,8 @@ struct Sigmaplusbuilder { std::array nHat = normalize3(flightVec); float flightDistance = std::sqrt(dot3(flightVec, flightVec)); - std::array pProton; - std::array pGamma1; + std::array pProton{}; + std::array pGamma1{}; fitter.getTrack(0).getPxPyPzGlo(pProton); fitter.getTrack(1).getPxPyPzGlo(pGamma1); @@ -770,7 +838,7 @@ struct Sigmaplusbuilder { } } if (discriminant < 0.f) { - return; + return -1; } fillCandStep(4); // real root @@ -802,7 +870,7 @@ struct Sigmaplusbuilder { } } if (!haveCandidate) { - return; + return -1; } fillCandStep(5); // valid root @@ -825,12 +893,12 @@ struct Sigmaplusbuilder { float photonOpeningAngle = std::acos(std::clamp(dot3(pPosDau, pNegDau) / std::sqrt(dot3(pPosDau, pPosDau) * dot3(pNegDau, pNegDau)), -1.f, 1.f)); // photon pointing angle: angle between the fitted photon momentum and the line from its conversion point to the p-gamma decay vertex - std::array convToDecVtx{secVtx[0] - photon.x(), secVtx[1] - photon.y(), secVtx[2] - photon.z()}; + std::array convToDecVtx{secVtx[0] - photon.x, secVtx[1] - photon.y, secVtx[2] - photon.z}; float photonPointingAngle = std::acos(std::clamp(dot3(pGamma1, convToDecVtx) / std::sqrt(dot3(pGamma1, pGamma1) * dot3(convToDecVtx, convToDecVtx)), -1.f, 1.f)); // photon DCA to PV: distance from the PV to the line through the conversion point along the photon momentum direction - std::array convPoint{photon.x(), photon.y(), photon.z()}; - std::array photonDir = normalize3({photon.px(), photon.py(), photon.pz()}); + std::array convPoint{photon.x, photon.y, photon.z}; + std::array photonDir = normalize3({photon.px, photon.py, photon.pz}); std::array pvToConv{pv[0] - convPoint[0], pv[1] - convPoint[1], pv[2] - convPoint[2]}; std::array pvToConvCrossDir = cross3(pvToConv, photonDir); float photonDcaToPV = std::sqrt(dot3(pvToConvCrossDir, pvToConvCrossDir)); @@ -843,7 +911,7 @@ struct Sigmaplusbuilder { bestMomGamma2[0], bestMomGamma2[1], bestMomGamma2[2], protonTrack.tpcNSigmaPr(), protonTrack.tofNSigmaPr(), posTrack.tpcNSigmaEl(), negTrack.tpcNSigmaEl(), - photon.mGamma(), photon.alpha(), photon.qtarm(), photon.v0radius(), + photon.mGamma, photon.alpha, photon.qtarm, photon.radius, photonOpeningAngle, photonPointingAngle, photonDcaToPV, protonTrack.itsNCls(), protonTrack.tpcNClsFound(), protonTrack.dcaXY(), protonTrack.dcaZ(), posTrack.itsNCls(), posTrack.tpcNClsFound(), negTrack.itsNCls(), negTrack.tpcNClsFound(), @@ -852,7 +920,8 @@ struct Sigmaplusbuilder { gammaPdgCode, gammaMotherPdgCode, gammaGMotherPdgCode, mcTrueVtx[0], mcTrueVtx[1], mcTrueVtx[2], mcTrueMomProton[0], mcTrueMomProton[1], mcTrueMomProton[2], - mcTrueMomGamma[0], mcTrueMomGamma[1], mcTrueMomGamma[2]); + mcTrueMomGamma[0], mcTrueMomGamma[1], mcTrueMomGamma[2], + mcTrueMomSigmaPlus[0], mcTrueMomSigmaPlus[1], mcTrueMomSigmaPlus[2]); } else { sigmaPlusCands(secVtx[0], secVtx[1], secVtx[2], radius, flightDistance, dcaProtonGamma, fitChi2, @@ -861,11 +930,12 @@ struct Sigmaplusbuilder { bestMomGamma2[0], bestMomGamma2[1], bestMomGamma2[2], protonTrack.tpcNSigmaPr(), protonTrack.tofNSigmaPr(), posTrack.tpcNSigmaEl(), negTrack.tpcNSigmaEl(), - photon.mGamma(), photon.alpha(), photon.qtarm(), photon.v0radius(), + photon.mGamma, photon.alpha, photon.qtarm, photon.radius, photonOpeningAngle, photonPointingAngle, photonDcaToPV, protonTrack.itsNCls(), protonTrack.tpcNClsFound(), protonTrack.dcaXY(), protonTrack.dcaZ(), posTrack.itsNCls(), posTrack.tpcNClsFound(), negTrack.itsNCls(), negTrack.tpcNClsFound()); } + return matchedSigmaId; } void initCCDB(aod::BCs::iterator const& bc) @@ -881,22 +951,17 @@ struct Sigmaplusbuilder { LOG(info) << "Task initialized for run " << mRunNumber << " with magnetic field " << mBz << " kZG"; } - void processData(CollisionsFull const& collisions, aod::V0Datas const& v0s, TracksFull const& tracks, aod::BCs const&) + void processData(CollisionsFull const& collisions, TracksFull const& tracks, aod::BCs const&) { for (const auto& collision : collisions) { initCCDB(collision.bc_as()); histos.fill(HIST("hVertexZ"), collision.posZ()); std::array pv{collision.posX(), collision.posY(), collision.posZ()}; - auto v0sThisCollision = v0s.sliceBy(v0PerCollision, collision.globalIndex()); - std::vector acceptedPhotons; - for (const auto& v0 : v0sThisCollision) { - if (selectPhoton(v0)) { - acceptedPhotons.push_back(v0); - } - } - auto tracksThisCollision = tracks.sliceBy(tracksPerCollision, collision.globalIndex()); + + auto acceptedPhotons = findPhotons(tracksThisCollision, pv); + std::vector acceptedProtons; for (const auto& track : tracksThisCollision) { if (selectProton(track)) { @@ -906,53 +971,52 @@ struct Sigmaplusbuilder { for (const auto& photon : acceptedPhotons) { for (const auto& proton : acceptedProtons) { - buildSigmaPlusCandidate(proton, photon, pv); + buildSigmaPlusCandidate(proton, photon, pv); } } } } PROCESS_SWITCH(Sigmaplusbuilder, processData, "Process data", true); - void processMc(CollisionsFullMC const& collisions, aod::V0Datas const& v0s, TracksFullMC const& tracks, aod::BCs const&, aod::McParticles const& mcParticles) + void processMc(CollisionsFullMC const& collisions, TracksFullMC const& tracks, aod::BCs const&, aod::McParticles const& mcParticles) { + std::vector matchedSigmaPlusMcIds; // Sigma+ MC indices that got at least one signal candidate + for (const auto& collision : collisions) { initCCDB(collision.bc_as()); histos.fill(HIST("hVertexZ"), collision.posZ()); std::array pv{collision.posX(), collision.posY(), collision.posZ()}; - auto v0sThisCollision = v0s.sliceBy(v0PerCollision, collision.globalIndex()); - std::vector acceptedPhotons; - for (const auto& v0 : v0sThisCollision) { - if (selectPhoton(v0)) { - acceptedPhotons.push_back(v0); - - histos.fill(HIST("MC/hPhotonTruthQA"), 0); - auto posTrack = v0.template posTrack_as(); - auto negTrack = v0.template negTrack_as(); - if (posTrack.has_mcParticle() && negTrack.has_mcParticle()) { - histos.fill(HIST("MC/hPhotonTruthQA"), 1); - auto mcPos = posTrack.template mcParticle_as(); - auto mcNeg = negTrack.template mcParticle_as(); + auto tracksThisCollision = tracks.sliceBy(tracksPerCollisionMC, collision.globalIndex()); + + auto acceptedPhotons = findPhotons(tracksThisCollision, pv); + for (const auto& photon : acceptedPhotons) { + histos.fill(HIST("MC/hPhotonTruthQA"), 0); + auto posTrack = photon.posTrack; + auto negTrack = photon.negTrack; + if (posTrack.has_mcParticle() && negTrack.has_mcParticle()) { + histos.fill(HIST("MC/hPhotonTruthQA"), 1); + auto mcPos = posTrack.template mcParticle_as(); + auto mcNeg = negTrack.template mcParticle_as(); - auto const& posMothers = mcPos.template mothers_as(); - auto const& negMothers = mcNeg.template mothers_as(); - if (!posMothers.empty() && !negMothers.empty()) { - auto mcGamma = posMothers.front(); - if (mcGamma.globalIndex() == negMothers.front().globalIndex() && mcGamma.pdgCode() == PDG_t::kGamma) { - histos.fill(HIST("MC/hPhotonTruthQA"), 2); - auto const& pi0Mothers = mcGamma.template mothers_as(); - if (!pi0Mothers.empty() && std::abs(pi0Mothers.front().pdgCode()) == PDG_t::kPi0) { - histos.fill(HIST("MC/hPhotonTruthQA"), 3); - auto const& sigmaMothers = pi0Mothers.front().template mothers_as(); - if (!sigmaMothers.empty() && std::abs(sigmaMothers.front().pdgCode()) == PDG_t::kSigmaPlus) { - histos.fill(HIST("MC/hPhotonTruthQA"), 4); - - histos.fill(HIST("Photon/True/hMass"), v0.mGamma()); - histos.fill(HIST("Photon/True/hPt"), v0.pt()); - histos.fill(HIST("Photon/True/hRadius"), v0.v0radius()); - histos.fill(HIST("Photon/True/h2ArmenterosPodolanski"), v0.alpha(), v0.qtarm()); - histos.fill(HIST("Photon/True/h2ConvPointXY"), v0.x(), v0.y()); - } + auto const& posMothers = mcPos.template mothers_as(); + auto const& negMothers = mcNeg.template mothers_as(); + if (!posMothers.empty() && !negMothers.empty()) { + auto mcGamma = posMothers.front(); + if (mcGamma.globalIndex() == negMothers.front().globalIndex() && mcGamma.pdgCode() == PDG_t::kGamma) { + histos.fill(HIST("MC/hPhotonTruthQA"), 2); + auto const& pi0Mothers = mcGamma.template mothers_as(); + if (!pi0Mothers.empty() && std::abs(pi0Mothers.front().pdgCode()) == PDG_t::kPi0) { + histos.fill(HIST("MC/hPhotonTruthQA"), 3); + auto const& sigmaMothers = pi0Mothers.front().template mothers_as(); + if (!sigmaMothers.empty() && std::abs(sigmaMothers.front().pdgCode()) == PDG_t::kSigmaPlus) { + histos.fill(HIST("MC/hPhotonTruthQA"), 4); + + histos.fill(HIST("Photon/True/hMass"), photon.mGamma); + histos.fill(HIST("Photon/True/hPt"), std::hypot(photon.px, photon.py)); + histos.fill(HIST("Photon/True/hRadius"), photon.radius); + histos.fill(HIST("Photon/True/h2ArmenterosPodolanski"), photon.alpha, photon.qtarm); + histos.fill(HIST("Photon/True/h2ConvPointXY"), photon.x, photon.y); } } } @@ -960,7 +1024,6 @@ struct Sigmaplusbuilder { } } - auto tracksThisCollision = tracks.sliceBy(tracksPerCollisionMC, collision.globalIndex()); std::vector acceptedProtons; for (const auto& track : tracksThisCollision) { if (selectProton(track)) { @@ -985,25 +1048,69 @@ struct Sigmaplusbuilder { for (const auto& photon : acceptedPhotons) { for (const auto& proton : acceptedProtons) { - buildSigmaPlusCandidate(proton, photon, pv); + int matchedId = buildSigmaPlusCandidate(proton, photon, pv); + if (matchedId >= 0) { + matchedSigmaPlusMcIds.push_back(matchedId); + } } } } // all generated Sigma+ -> p pi0 decays, regardless of reconstruction for (const auto& mcPart : mcParticles) { - if (isSigmaPlusToProtonPi0(mcPart)) { - histos.fill(HIST("MC/hGenSigmaPlusPt"), mcPart.pt()); + if (!isSigmaPlusToProtonPi0(mcPart)) { + continue; } + histos.fill(HIST("MC/hGenSigmaPlusPt"), mcPart.pt()); + + bool wasReconstructed = std::find(matchedSigmaPlusMcIds.begin(), matchedSigmaPlusMcIds.end(), mcPart.globalIndex()) != matchedSigmaPlusMcIds.end(); + if (wasReconstructed) { + continue; + } + + // this true Sigma+ never made it into any signal candidate: still record its truth info, + // with the reconstructed-side columns set to -999 + int pdgProton = mcPart.pdgCode() > 0 ? PDG_t::kProton : PDG_t::kProtonBar; + std::array genDecVtx{-999.f, -999.f, -999.f}; + std::array genMomProton{-999.f, -999.f, -999.f}; + for (const auto& daughter : mcPart.template daughters_as()) { + if (daughter.pdgCode() == pdgProton) { + genDecVtx = {daughter.vx(), daughter.vy(), daughter.vz()}; + genMomProton = {daughter.px(), daughter.py(), daughter.pz()}; + break; + } + } + + sigmaPlusCandsMC(-999.f, -999.f, -999.f, + -999.f, -999.f, -999.f, -999.f, + -999.f, -999.f, -999.f, + -999.f, -999.f, -999.f, + -999.f, -999.f, -999.f, + -999.f, -999.f, + -999.f, -999.f, + -999.f, -999.f, -999.f, -999.f, + -999.f, -999.f, -999.f, + 0, -999, -999.f, -999.f, + 0, -999, 0, -999, + true, + pdgProton, mcPart.pdgCode(), + 0, 0, 0, + genDecVtx[0], genDecVtx[1], genDecVtx[2], + genMomProton[0], genMomProton[1], genMomProton[2], + -999.f, -999.f, -999.f, + mcPart.px(), mcPart.py(), mcPart.pz()); } } PROCESS_SWITCH(Sigmaplusbuilder, processMc, "Process MC", false); void processFindable(CollisionsFullMC const& collisions, aod::V0Datas const& v0s, TracksFullMC const& tracks, aod::McParticles const&) { + constexpr int MinDauTpcCls = 90; for (const auto& collision : collisions) { auto tracksThisCollision = tracks.sliceBy(tracksPerCollisionMC, collision.globalIndex()); auto v0sThisCollision = v0s.sliceBy(v0PerCollision, collision.globalIndex()); + std::array pv{collision.posX(), collision.posY(), collision.posZ()}; + auto acceptedPhotons = findPhotons(tracksThisCollision, pv); for (const auto& protonTrack : tracksThisCollision) { if (!protonTrack.has_mcParticle()) { @@ -1035,7 +1142,7 @@ struct Sigmaplusbuilder { continue; } - if (electronTrack.tpcNClsFound() < 90 || electronTrack.sign() > 0) { + if (electronTrack.tpcNClsFound() < MinDauTpcCls || electronTrack.sign() > 0) { continue; } @@ -1055,7 +1162,7 @@ struct Sigmaplusbuilder { continue; } - if (positronTrack.tpcNClsFound() < 90 || positronTrack.sign() < 0) { + if (positronTrack.tpcNClsFound() < MinDauTpcCls || positronTrack.sign() < 0) { continue; } @@ -1121,10 +1228,10 @@ struct Sigmaplusbuilder { std::array recoGammaMom{electronTrack.px() + positronTrack.px(), electronTrack.py() + positronTrack.py(), electronTrack.pz() + positronTrack.pz()}; float recoGammaP = std::sqrt(dot3(recoGammaMom, recoGammaMom)); - constexpr float electronMass = o2::constants::physics::MassElectron; + constexpr float ElectronMass = o2::constants::physics::MassElectron; float electronP2 = electronTrack.px() * electronTrack.px() + electronTrack.py() * electronTrack.py() + electronTrack.pz() * electronTrack.pz(); float positronP2 = positronTrack.px() * positronTrack.px() + positronTrack.py() * positronTrack.py() + positronTrack.pz() * positronTrack.pz(); - float pairEnergy = std::sqrt(electronP2 + electronMass * electronMass) + std::sqrt(positronP2 + electronMass * electronMass); + float pairEnergy = std::sqrt(electronP2 + ElectronMass * ElectronMass) + std::sqrt(positronP2 + ElectronMass * ElectronMass); float pairMass2 = pairEnergy * pairEnergy - recoGammaP * recoGammaP; histos.fill(HIST("Findable/hElectronPositronMass"), std::sqrt(std::max(pairMass2, 0.f))); @@ -1155,6 +1262,16 @@ struct Sigmaplusbuilder { break; } } + + histos.fill(HIST("Findable/hPhotonSearchPresence"), 0); + for (const auto& photon : acceptedPhotons) { + bool sameChargeMatched = photon.posTrack.globalIndex() == positronTrack.globalIndex() && photon.negTrack.globalIndex() == electronTrack.globalIndex(); + if (sameChargeMatched) { + histos.fill(HIST("Findable/hPhotonSearchPresence"), 1); + histos.fill(HIST("Findable/hSigmaPlusPtFoundByPhotonSearch"), mcSigmaPlus.pt()); + break; + } + } histos.fill(HIST("Findable/hElectronPt"), electronTrack.pt()); histos.fill(HIST("Findable/hPositronPt"), positronTrack.pt()); histos.fill(HIST("Findable/hPhotonConversionRadius"), std::hypot(electronVertex[0], electronVertex[1]));