From 4307e998c4b801e8849463dc389ea123d4d7cde8 Mon Sep 17 00:00:00 2001 From: Henrik Fribert Date: Mon, 14 Sep 2026 18:16:07 +0200 Subject: [PATCH] SigmaPlus improved selections and tables --- PWGLF/DataModel/LFKinkDecayTables.h | 73 ++- .../Strangeness/sigmaHadCorr.cxx | 40 +- .../Strangeness/sigmaplusbuilder.cxx | 589 ++++++++++++++---- 3 files changed, 547 insertions(+), 155 deletions(-) diff --git a/PWGLF/DataModel/LFKinkDecayTables.h b/PWGLF/DataModel/LFKinkDecayTables.h index 4443727a151..a8a5cffa8a3 100644 --- a/PWGLF/DataModel/LFKinkDecayTables.h +++ b/PWGLF/DataModel/LFKinkDecayTables.h @@ -163,13 +163,12 @@ DECLARE_SOA_TABLE(SlimKinkCandsMC, "AOD", "SLIMKINKCANDSMC", namespace sigmapluscand { -DECLARE_SOA_COLUMN(XDecVtx, xDecVtx, float); //! Decay vertex of the candidate (x direction) -DECLARE_SOA_COLUMN(YDecVtx, yDecVtx, float); //! Decay vertex of the candidate (y direction) -DECLARE_SOA_COLUMN(ZDecVtx, zDecVtx, float); //! Decay vertex of the candidate (z direction) -DECLARE_SOA_COLUMN(Radius, radius, float); //! Decay radius of the candidate (cm) -DECLARE_SOA_COLUMN(FlightDistance, flightDistance, float); //! Flight distance of the candidate (PV to decay vertex, cm) -DECLARE_SOA_COLUMN(DcaProtonGamma, dcaProtonGamma, float); //! DCA between proton and photon at the fitted vertex (cm) -DECLARE_SOA_COLUMN(Chi2, chi2, float); //! chi2 of the proton-photon vertex fit +DECLARE_SOA_COLUMN(XDecVtx, xDecVtx, float); //! Decay vertex of the candidate (x direction) +DECLARE_SOA_COLUMN(YDecVtx, yDecVtx, float); //! Decay vertex of the candidate (y direction) +DECLARE_SOA_COLUMN(ZDecVtx, zDecVtx, float); //! Decay vertex of the candidate (z direction) +DECLARE_SOA_COLUMN(FlightDistance, flightDistance, float); //! 3D distance from the PV to the decay vertex (cm) +DECLARE_SOA_COLUMN(TransDecayRadius, transDecayRadius, float); //! Transverse decay radius of the candidate (cm), stored directly in the slim tables (no XDecVtx/YDecVtx there) +DECLARE_SOA_COLUMN(DcaProtonGamma, dcaProtonGamma, float); //! DCA between proton and photon at the fitted vertex (cm) DECLARE_SOA_COLUMN(PxProton, pxProton, float); //! Px of the proton DECLARE_SOA_COLUMN(PyProton, pyProton, float); //! Py of the proton @@ -194,6 +193,11 @@ DECLARE_SOA_COLUMN(PhotonOpeningAngle, photonOpeningAngle, float); //! Opening DECLARE_SOA_COLUMN(PhotonPointingAngle, photonPointingAngle, float); //! Angle between the photon momentum and the line from its conversion point to the candidate decay vertex (rad) DECLARE_SOA_COLUMN(PhotonDcaToPV, photonDcaToPV, float); //! DCA of the photon's flight line to the primary vertex (cm) +DECLARE_SOA_COLUMN(RootCenter, rootCenter, float); //! -coefB/(2*coefA) of the missing-photon quadratic solve: negative flags an unphysical phase-space point +DECLARE_SOA_COLUMN(AntiSigmaPointingAngle, antiSigmaPointingAngle, float); //! Angle between the field-unbent proton momentum (from its original reference point) and the PV->decay-vertex direction +DECLARE_SOA_COLUMN(CandDcaToPV, candDcaToPV, float); //! DCA of the candidate's total reconstructed momentum line to the PV (cm) + +DECLARE_SOA_COLUMN(ProtonSign, protonSign, int); //! Charge sign of the proton track (= sign of the whole candidate, since the photon is neutral) DECLARE_SOA_COLUMN(ProtonItsNCls, protonItsNCls, uint8_t); //! Number of ITS clusters of the proton track DECLARE_SOA_COLUMN(ProtonTpcNCls, protonTpcNCls, int16_t); //! Number of found TPC clusters of the proton track DECLARE_SOA_COLUMN(ProtonDcaXY, protonDcaXY, float); //! DCA of the proton track to the primary vertex, xy (cm) @@ -205,18 +209,15 @@ DECLARE_SOA_COLUMN(PhotonNegItsNCls, photonNegItsNCls, uint8_t); //! Number of I DECLARE_SOA_COLUMN(PhotonNegTpcNCls, photonNegTpcNCls, int16_t); //! Number of found TPC clusters of the photon's negative daughter // MC columns -DECLARE_SOA_COLUMN(IsSignal, isSignal, bool); //! True if proton and photon are MC-truth matched to the same Sigma+ +DECLARE_SOA_COLUMN(CollisionIdCheck, collisionIdCheck, bool); //! True if the proton's collision ID matches the reconstructed collision ID -DECLARE_SOA_COLUMN(ProtonPdgCode, protonPdgCode, int); //! PDG code of the proton's MC particle -DECLARE_SOA_COLUMN(ProtonMotherPdgCode, protonMotherPdgCode, int); //! PDG code of the proton's MC mother -DECLARE_SOA_COLUMN(GammaPdgCode, gammaPdgCode, int); //! PDG code of the measured photon's MC particle -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(IsSignal, isSignal, bool); //! True if the proton and photon share the same true Sigma+ mother 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(DecayRadiusMC, decayRadiusMC, float); //! MC-truth Sigma+ decay radius +DECLARE_SOA_COLUMN(MassMC, massMC, float); //! MC-truth invariant mass of the Sigma+ mother 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 @@ -229,6 +230,12 @@ DECLARE_SOA_COLUMN(PzGammaMC, pzGammaMC, float); //! MC-truth momentum o // DYNAMIC COLUMNS +DECLARE_SOA_DYNAMIC_COLUMN(Radius, radius, //! Transverse decay radius of the candidate (cm), from the beam axis + [](float xDecVtx, float yDecVtx) -> float { return std::hypot(xDecVtx, yDecVtx); }); + +DECLARE_SOA_DYNAMIC_COLUMN(Chi2, chi2, //! chi2 of the proton-photon vertex fit; dcaProtonGamma = sqrt(chi2) by construction in the builder + [](float dcaProtonGamma) -> float { return dcaProtonGamma * dcaProtonGamma; }); + DECLARE_SOA_DYNAMIC_COLUMN(PxSigmaPlus, pxSigmaPlus, //! Px of the Sigma+ candidate [](float pxProton, float pxGamma1, float pxGamma2) -> float { return pxProton + pxGamma1 + pxGamma2; }); @@ -258,9 +265,8 @@ DECLARE_SOA_DYNAMIC_COLUMN(YSigmaPlusMC, ySigmaPlusMC, //! True rapidity of the } // namespace sigmapluscand DECLARE_SOA_TABLE(SigmaPlusCands, "AOD", "SIGMAPLUSCANDS", - sigmapluscand::XDecVtx, sigmapluscand::YDecVtx, sigmapluscand::ZDecVtx, - sigmapluscand::Radius, sigmapluscand::FlightDistance, - sigmapluscand::DcaProtonGamma, sigmapluscand::Chi2, + sigmapluscand::XDecVtx, sigmapluscand::YDecVtx, sigmapluscand::ZDecVtx, sigmapluscand::FlightDistance, + sigmapluscand::DcaProtonGamma, sigmapluscand::PxProton, sigmapluscand::PyProton, sigmapluscand::PzProton, sigmapluscand::PxGamma1, sigmapluscand::PyGamma1, sigmapluscand::PzGamma1, sigmapluscand::PxGamma2, sigmapluscand::PyGamma2, sigmapluscand::PzGamma2, @@ -268,10 +274,14 @@ DECLARE_SOA_TABLE(SigmaPlusCands, "AOD", "SIGMAPLUSCANDS", sigmapluscand::NSigmaTPCElPos, sigmapluscand::NSigmaTPCElNeg, sigmapluscand::PhotonMass, sigmapluscand::PhotonAlpha, sigmapluscand::PhotonQt, sigmapluscand::PhotonConvRadius, sigmapluscand::PhotonOpeningAngle, sigmapluscand::PhotonPointingAngle, sigmapluscand::PhotonDcaToPV, + sigmapluscand::RootCenter, sigmapluscand::AntiSigmaPointingAngle, sigmapluscand::CandDcaToPV, + sigmapluscand::ProtonSign, sigmapluscand::ProtonItsNCls, sigmapluscand::ProtonTpcNCls, sigmapluscand::ProtonDcaXY, sigmapluscand::ProtonDcaZ, sigmapluscand::PhotonPosItsNCls, sigmapluscand::PhotonPosTpcNCls, sigmapluscand::PhotonNegItsNCls, sigmapluscand::PhotonNegTpcNCls, // dynamic columns + sigmapluscand::Radius, + sigmapluscand::Chi2, sigmapluscand::PxSigmaPlus, sigmapluscand::PySigmaPlus, sigmapluscand::PzSigmaPlus, @@ -279,9 +289,8 @@ DECLARE_SOA_TABLE(SigmaPlusCands, "AOD", "SIGMAPLUSCANDS", sigmapluscand::MassSigmaPlus); DECLARE_SOA_TABLE(SigmaPlusCandsMC, "AOD", "SIGMAPLUSMC", - sigmapluscand::XDecVtx, sigmapluscand::YDecVtx, sigmapluscand::ZDecVtx, - sigmapluscand::Radius, sigmapluscand::FlightDistance, - sigmapluscand::DcaProtonGamma, sigmapluscand::Chi2, + sigmapluscand::XDecVtx, sigmapluscand::YDecVtx, sigmapluscand::ZDecVtx, sigmapluscand::FlightDistance, + sigmapluscand::DcaProtonGamma, sigmapluscand::PxProton, sigmapluscand::PyProton, sigmapluscand::PzProton, sigmapluscand::PxGamma1, sigmapluscand::PyGamma1, sigmapluscand::PzGamma1, sigmapluscand::PxGamma2, sigmapluscand::PyGamma2, sigmapluscand::PzGamma2, @@ -289,17 +298,21 @@ DECLARE_SOA_TABLE(SigmaPlusCandsMC, "AOD", "SIGMAPLUSMC", sigmapluscand::NSigmaTPCElPos, sigmapluscand::NSigmaTPCElNeg, sigmapluscand::PhotonMass, sigmapluscand::PhotonAlpha, sigmapluscand::PhotonQt, sigmapluscand::PhotonConvRadius, sigmapluscand::PhotonOpeningAngle, sigmapluscand::PhotonPointingAngle, sigmapluscand::PhotonDcaToPV, + sigmapluscand::RootCenter, sigmapluscand::AntiSigmaPointingAngle, sigmapluscand::CandDcaToPV, + sigmapluscand::ProtonSign, sigmapluscand::ProtonItsNCls, sigmapluscand::ProtonTpcNCls, sigmapluscand::ProtonDcaXY, sigmapluscand::ProtonDcaZ, sigmapluscand::PhotonPosItsNCls, sigmapluscand::PhotonPosTpcNCls, sigmapluscand::PhotonNegItsNCls, sigmapluscand::PhotonNegTpcNCls, + sigmapluscand::CollisionIdCheck, sigmapluscand::IsSignal, - sigmapluscand::ProtonPdgCode, sigmapluscand::ProtonMotherPdgCode, - sigmapluscand::GammaPdgCode, sigmapluscand::GammaMotherPdgCode, sigmapluscand::GammaGMotherPdgCode, sigmapluscand::XDecVtxMC, sigmapluscand::YDecVtxMC, sigmapluscand::ZDecVtxMC, sigmapluscand::PxProtonMC, sigmapluscand::PyProtonMC, sigmapluscand::PzProtonMC, sigmapluscand::PxGammaMC, sigmapluscand::PyGammaMC, sigmapluscand::PzGammaMC, sigmapluscand::PxSigmaPlusMC, sigmapluscand::PySigmaPlusMC, sigmapluscand::PzSigmaPlusMC, + sigmapluscand::DecayRadiusMC, sigmapluscand::MassMC, // dynamic columns + sigmapluscand::Radius, + sigmapluscand::Chi2, sigmapluscand::PxSigmaPlus, sigmapluscand::PySigmaPlus, sigmapluscand::PzSigmaPlus, @@ -309,7 +322,11 @@ DECLARE_SOA_TABLE(SigmaPlusCandsMC, "AOD", "SIGMAPLUSMC", sigmapluscand::YSigmaPlusMC); DECLARE_SOA_TABLE(SlimSigmaPlusCands, "AOD", "SLIMSIGMAPLUS", - sigmapluscand::Radius, sigmapluscand::DcaProtonGamma, + sigmapluscand::TransDecayRadius, + sigmapluscand::CandDcaToPV, + sigmapluscand::DcaProtonGamma, + sigmapluscand::ProtonSign, + sigmapluscand::ProtonDcaXY, sigmapluscand::ProtonDcaZ, sigmapluscand::PxProton, sigmapluscand::PyProton, sigmapluscand::PzProton, sigmapluscand::PxGamma1, sigmapluscand::PyGamma1, sigmapluscand::PzGamma1, sigmapluscand::PxGamma2, sigmapluscand::PyGamma2, sigmapluscand::PzGamma2, @@ -325,16 +342,20 @@ DECLARE_SOA_TABLE(SlimSigmaPlusCands, "AOD", "SLIMSIGMAPLUS", sigmapluscand::MassSigmaPlus); DECLARE_SOA_TABLE(SlimSigmaPlusCandsMC, "AOD", "SLIMSIGMAPLUSMC", - sigmapluscand::Radius, sigmapluscand::DcaProtonGamma, + sigmapluscand::TransDecayRadius, + sigmapluscand::CandDcaToPV, + sigmapluscand::DcaProtonGamma, + sigmapluscand::ProtonSign, + sigmapluscand::ProtonDcaXY, sigmapluscand::ProtonDcaZ, sigmapluscand::PxProton, sigmapluscand::PyProton, sigmapluscand::PzProton, sigmapluscand::PxGamma1, sigmapluscand::PyGamma1, sigmapluscand::PzGamma1, sigmapluscand::PxGamma2, sigmapluscand::PyGamma2, sigmapluscand::PzGamma2, sigmapluscand::NSigmaTPCProton, sigmapluscand::NSigmaTOFProton, sigmapluscand::NSigmaTPCElPos, sigmapluscand::NSigmaTPCElNeg, sigmapluscand::PhotonMass, - sigmapluscand::ProtonPdgCode, sigmapluscand::ProtonMotherPdgCode, - sigmapluscand::GammaPdgCode, sigmapluscand::GammaMotherPdgCode, sigmapluscand::GammaGMotherPdgCode, - sigmapluscand::DecayRadiusMC, + sigmapluscand::CollisionIdCheck, + sigmapluscand::IsSignal, + sigmapluscand::DecayRadiusMC, sigmapluscand::MassMC, sigmapluscand::PxSigmaPlusMC, sigmapluscand::PySigmaPlusMC, sigmapluscand::PzSigmaPlusMC, // dynamic columns diff --git a/PWGLF/TableProducer/Strangeness/sigmaHadCorr.cxx b/PWGLF/TableProducer/Strangeness/sigmaHadCorr.cxx index 58acb86cfb7..6f2a8a2336c 100644 --- a/PWGLF/TableProducer/Strangeness/sigmaHadCorr.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaHadCorr.cxx @@ -113,6 +113,8 @@ struct SigmaHadCorr { Configurable cutMaxKStar{"cutMaxKStar", 1.5, "Maximum k* for Sigma-hadron pairs (GeV/c)"}; Configurable minPtSigma{"minPtSigma", 1.f, "Minimum pT for Sigma candidates (GeV/c)"}; + Configurable cutRapMotherMC{"cutRapMotherMC", 1.0f, "Rapidity cut for generated mother Sigma in MC"}; + Configurable cutPtGenMC{"cutPtGenMC", 0.5f, "Minimum pT for generated Sigma particles in MC"}; Configurable useRecalculatedSigmaMomentum{"useRecalculatedSigmaMomentum", true, "If true, compute k* using Sigma momentum recalculated from daughter kinematics"}; Configurable cutDCAtoPVSigma{"cutDCAtoPVSigma", 0.1f, "Max DCA to primary vertex for Sigma candidates (cm)"}; Configurable cutSigmaRadius{"cutSigmaRadius", 20.f, "Minimum radius for Sigma candidates (cm)"}; @@ -237,29 +239,29 @@ struct SigmaHadCorr { float versorZ = sigmaPz / pMother; float eChDau = std::sqrt(massChargedDau * massChargedDau + sigmaDauPx * sigmaDauPx + sigmaDauPy * sigmaDauPy + sigmaDauPz * sigmaDauPz); float a = versorX * sigmaDauPx + versorY * sigmaDauPy + versorZ * sigmaDauPz; - float K = massSigma * massSigma + massChargedDau * massChargedDau - massNeutralDau * massNeutralDau; - 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) < epsilon) { + float coefK = massSigma * massSigma + massChargedDau * massChargedDau - massNeutralDau * massNeutralDau; + float coefA = 4.f * (eChDau * eChDau - a * a); + float coefB = -4.f * a * coefK; + float coefC = 4.f * eChDau * eChDau * massSigma * massSigma - coefK * coefK; + if (std::abs(coefA) < epsilon) { return -999.f; } - float D = B * B - 4.f * A * C; - if (D < 0.f) { + float discriminant = coefB * coefB - 4.f * coefA * coefC; + if (discriminant < 0.f) { return -999.f; } - float sqrtD = std::sqrt(D); - float P1 = (-B + sqrtD) / (2.f * A); - float P2 = (-B - sqrtD) / (2.f * A); - if (P2 < 0.f && P1 < 0.f) { + float sqrtDiscriminant = std::sqrt(discriminant); + float root1 = (-coefB + sqrtDiscriminant) / (2.f * coefA); + float root2 = (-coefB - sqrtDiscriminant) / (2.f * coefA); + if (root2 < 0.f && root1 < 0.f) { return -999.f; } - if (P2 < 0.f) { - return P1; + if (root2 < 0.f) { + return root1; } - float p1Diff = std::abs(P1 - pMother); - float p2Diff = std::abs(P2 - pMother); - return (p1Diff < p2Diff) ? P1 : P2; + float p1Diff = std::abs(root1 - pMother); + float p2Diff = std::abs(root2 - pMother); + return (p1Diff < p2Diff) ? root1 : root2; } std::array getSigmaMomentumForKstar(float sigmaPx, float sigmaPy, float sigmaPz, float sigmaDauPx, float sigmaDauPy, float sigmaDauPz) @@ -794,6 +796,12 @@ struct SigmaHadCorr { if (!isValidMother) { continue; } + if (std::abs(mcPart.y()) > cutRapMotherMC) { + continue; + } + if (mcPart.pt() < cutPtGenMC) { + continue; + } bool hasChargedDaughter = false; std::array genDecVtx{-999.f, -999.f, -999.f}; int daugPdgCode = 0; diff --git a/PWGLF/TableProducer/Strangeness/sigmaplusbuilder.cxx b/PWGLF/TableProducer/Strangeness/sigmaplusbuilder.cxx index 1b521798bf3..7969042b519 100644 --- a/PWGLF/TableProducer/Strangeness/sigmaplusbuilder.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaplusbuilder.cxx @@ -18,11 +18,13 @@ #include "Common/Core/RecoDecay.h" #include "Common/Core/trackUtilities.h" +#include "Common/DataModel/EventSelection.h" #include "Common/DataModel/PIDResponseTOF.h" #include "Common/DataModel/PIDResponseTPC.h" #include "Common/DataModel/TrackSelectionTables.h" #include +#include #include #include #include @@ -56,37 +58,61 @@ using namespace o2::framework; using TracksFull = soa::Join; using TracksFullMC = soa::Join; -using CollisionsFull = aod::Collisions; -using CollisionsFullMC = soa::Join; +using CollisionsFull = soa::Join; +using CollisionsFullMC = soa::Join; struct Sigmaplusbuilder { + // event selection + Configurable cutZVertex{"cutZVertex", 10.0f, "Accepted z-vertex range (cm)"}; + // photon (PCM) selection Configurable photonMaxMass{"photonMaxMass", 0.20, "Max photon mass (GeV/c^2)"}; Configurable photonMinRapidity{"photonMinRapidity", -0.8, "Min photon rapidity"}; Configurable photonMaxRapidity{"photonMaxRapidity", 0.8, "Max photon rapidity"}; + Configurable cutRapMotherMC{"cutRapMotherMC", 1.0f, "Rapidity cut for generated mother Sigma+ in MC"}; + Configurable cutPtGenMC{"cutPtGenMC", 0.5f, "Minimum pT for generated Sigma+ in MC"}; Configurable photonDauEtaMin{"photonDauEtaMin", -0.8, "Min eta of photon daughter tracks"}; Configurable photonDauEtaMax{"photonDauEtaMax", 0.8, "Max eta of photon daughter tracks"}; Configurable photonMinRadius{"photonMinRadius", 3.0, "Min photon conversion radius (cm)"}; Configurable photonMaxRadius{"photonMaxRadius", 115., "Max photon conversion radius (cm)"}; Configurable photonMinV0cospa{"photonMinV0cospa", 0.80, "Min V0 CosPA"}; Configurable photonMaxDCAV0Dau{"photonMaxDCAV0Dau", 3.5, "Max DCA between photon daughters (cm)"}; + Configurable photonMaxOpeningAngle{"photonMaxOpeningAngle", 0.4, "Max opening angle between the photon's e+/e- daughter momenta (rad)"}; + Configurable photonMaxDeltaTheta{"photonMaxDeltaTheta", 0.15, "Max |theta_pos - theta_neg| of the photon's daughter tracks (rad)"}; Configurable photonMaxQt{"photonMaxQt", 0.15, "Max Armenteros qT for photons (GeV/c)"}; Configurable photonMaxAlpha{"photonMaxAlpha", 1.0, "Max |Armenteros alpha| for photons"}; Configurable photonDauMinTPCNSigmaEl{"photonDauMinTPCNSigmaEl", -5., "Min TPC nSigma_el of the photon daughters"}; Configurable photonDauMaxTPCNSigmaEl{"photonDauMaxTPCNSigmaEl", 5., "Max TPC nSigma_el of the photon daughters"}; + Configurable photonDauMinTpcNCls{"photonDauMinTpcNCls", 30, "Min number of found TPC clusters for the photon (V0) daughter tracks"}; // proton selection Configurable protonMinPt{"protonMinPt", 0.3, "Minimum proton pT (GeV/c)"}; Configurable protonMaxEta{"protonMaxEta", 0.9, "Maximum |eta| for proton track"}; + Configurable protonMinTpcNCls{"protonMinTpcNCls", 80, "Min number of found TPC clusters for the proton track"}; Configurable protonMaxTPCNSigma{"protonMaxTPCNSigma", 4, "Max |TPC nSigma_pr| for proton"}; Configurable protonMaxTOFNSigma{"protonMaxTOFNSigma", 4, "Max |TOF nSigma_pr| for proton, if TOF available"}; Configurable protonPtMinRequireTOF{"protonPtMinRequireTOF", 0.75, "Above this pT, require TOF PID for proton"}; + Configurable protonRequireTofHit{"protonRequireTofHit", false, "Above protonPtMinRequireTOF, reject a proton with no TOF hit at all"}; + Configurable protonMinDcaToPV{"protonMinDcaToPV", 0.005, "Min DCAxy of the proton track to the PV (cm)"}; + Configurable protonMaxDcaToPV{"protonMaxDcaToPV", 5.0, "Max DCAxy of the proton track to the PV (cm)"}; // proton-photon candidate selection - Configurable candMaxDcaProtonGamma{"candMaxDcaProtonGamma", 5.0, "Max DCA between proton and photon at the fitted vertex (cm)"}; + Configurable candMaxDcaProtonGamma{"candMaxDcaProtonGamma", 0.5, "Max DCA between proton and photon at the fitted vertex (cm)"}; + Configurable candMaxDcaToPV{"candMaxDcaToPV", 0.1, "Max DCA of the candidate's total (reconstructed) momentum line to the PV (cm)"}; + Configurable candRejectNegRootCenter{"candRejectNegRootCenter", true, "Reject candidates with rootCenter<0"}; + Configurable candMaxRootCenter{"candMaxRootCenter", 15, "Max rootCenter=-coefB/(2*coefA) (GeV/c)"}; + Configurable candMinAntiSigmaPointingAngle{"candMinAntiSigmaPointingAngle", 0.0, "Min AntiSigmaPointingAngle (rad)"}; + Configurable candMaxAntiSigmaPointingAngle{"candMaxAntiSigmaPointingAngle", 0.3, "Max AntiSigmaPointingAngle (rad)"}; + Configurable candMaxSigmaMass{"candMaxSigmaMass", 1.35, "Max reconstructed Sigma+ candidate mass (GeV/c^2)"}; + Configurable candMaxRapidity{"candMaxRapidity", 0.9, "Max |rapidity| of the reconstructed Sigma+ candidate"}; Configurable candMinRadius{"candMinRadius", 1.0, "Min candidate decay radius (cm)"}; Configurable candMaxRadius{"candMaxRadius", 100., "Max candidate decay radius (cm)"}; + Configurable candMinFlightDistance{"candMinFlightDistance", 0.0, "Min 3D distance from PV to candidate decay vertex (cm)"}; + Configurable candMaxFlightDistance{"candMaxFlightDistance", 250.0, "Max 3D distance from PV to candidate decay vertex (cm)"}; + Configurable candMaxPhotonOpeningAngle{"candMaxPhotonOpeningAngle", 3.15, "Max photon opening angle (rad), recomputed from the daughters' raw track momenta"}; + Configurable candMaxPhotonPointingAngle{"candMaxPhotonPointingAngle", 0.5, "Max angle between the photon's fitted momentum and the decay-vertex-to-conversion-point line (rad)"}; + Configurable candMaxPhotonDcaToPV{"candMaxPhotonDcaToPV", 250.0, "Max DCA of the photon's flight line to the PV (cm)"}; // missing-photon discriminant retry // resolution-shaped function and tuned parameters from Run 2 used @@ -143,7 +169,7 @@ struct Sigmaplusbuilder { const AxisSpec axisVertexZ{100, -15., 15., "vrtx_{Z} (cm)"}; - const AxisSpec axisPhotonSel{11, -0.5, 10.5, "selection step"}; + const AxisSpec axisPhotonSel{14, -0.5, 13.5, "selection step"}; const AxisSpec axisPhotonMass{200, 0., 0.3, "m_{#gamma} (GeV/c^{2})"}; const AxisSpec axisPhotonPt{100, 0., 5., "#it{p}_{T,#gamma} (GeV/c)"}; const AxisSpec axisPhotonRadius{200, 0., 200., "R_{conv} (cm)"}; @@ -151,10 +177,14 @@ struct Sigmaplusbuilder { const AxisSpec axisQt{100, 0., 0.3, "q_{T,AP} (GeV/c)"}; const AxisSpec axisConvXY{200, -100., 100., "conv. point (cm)"}; const AxisSpec axisNSigmaEl{100, -10., 10., "n#sigma_{el}"}; + const AxisSpec axisTpcNCls{160, -0.5, 159.5, "TPC clusters"}; + const AxisSpec axisPhotonOpeningAngle{180, 0., 3.15, "opening angle (rad)"}; + const AxisSpec axisPhotonDeltaTheta{200, -1., 1., "#Delta#theta (rad)"}; - const AxisSpec axisProtonSel{5, -0.5, 4.5, "selection step"}; + const AxisSpec axisProtonSel{7, -0.5, 6.5, "selection step"}; const AxisSpec axisProtonPt{100, 0., 5., "#it{p}_{T,p} (GeV/c)"}; const AxisSpec axisNSigma{100, -5., 5., "n#sigma"}; + const AxisSpec axisProtonDcaToPV{1000, 0., 5., "DCA_{xy,p} to PV (cm)"}; histos.add("hVertexZ", "hVertexZ", kTH1F, {axisVertexZ}); @@ -163,14 +193,17 @@ struct Sigmaplusbuilder { hPhotonSel->GetXaxis()->SetBinLabel(1, "All"); hPhotonSel->GetXaxis()->SetBinLabel(2, "Neg eta"); hPhotonSel->GetXaxis()->SetBinLabel(3, "Pos eta"); - hPhotonSel->GetXaxis()->SetBinLabel(4, "TPC nSigma_{el}"); - 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, "Mass"); + hPhotonSel->GetXaxis()->SetBinLabel(4, "TPC clusters"); + hPhotonSel->GetXaxis()->SetBinLabel(5, "TPC nSigma_{el}"); + hPhotonSel->GetXaxis()->SetBinLabel(6, "DCA daughters"); + hPhotonSel->GetXaxis()->SetBinLabel(7, "Radius"); + hPhotonSel->GetXaxis()->SetBinLabel(8, "Opening angle"); + hPhotonSel->GetXaxis()->SetBinLabel(9, "Delta theta"); + hPhotonSel->GetXaxis()->SetBinLabel(10, "CosPA"); + hPhotonSel->GetXaxis()->SetBinLabel(11, "Rapidity"); + hPhotonSel->GetXaxis()->SetBinLabel(12, "Qt"); + hPhotonSel->GetXaxis()->SetBinLabel(13, "Alpha"); + hPhotonSel->GetXaxis()->SetBinLabel(14, "Mass"); histos.add("Photon/hMass", "Photon/hMass", kTH1F, {axisPhotonMass}); histos.add("Photon/hPt", "Photon/hPt", kTH1F, {axisPhotonPt}); @@ -179,42 +212,74 @@ struct Sigmaplusbuilder { histos.add("Photon/h2ConvPointXY", "Photon/h2ConvPointXY", kTH2F, {axisConvXY, axisConvXY}); histos.add("Photon/h2TPCNSigmaElPosVsPt", "Photon/h2TPCNSigmaElPosVsPt", kTH2F, {axisPhotonPt, axisNSigmaEl}); histos.add("Photon/h2TPCNSigmaElNegVsPt", "Photon/h2TPCNSigmaElNegVsPt", kTH2F, {axisPhotonPt, axisNSigmaEl}); + histos.add("Photon/h2TPCNClsPosVsPt", "Photon/h2TPCNClsPosVsPt", kTH2F, {axisPhotonPt, axisTpcNCls}); + histos.add("Photon/h2TPCNClsNegVsPt", "Photon/h2TPCNClsNegVsPt", kTH2F, {axisPhotonPt, axisTpcNCls}); + histos.add("Photon/hOpeningAngle", "Photon/hOpeningAngle", kTH1F, {axisPhotonOpeningAngle}); + histos.add("Photon/hDeltaTheta", "Photon/hDeltaTheta", kTH1F, {axisPhotonDeltaTheta}); histos.add("Proton/hSelectionCounter", "Proton/hSelectionCounter", kTH1F, {axisProtonSel}); auto hProtonSel = histos.get(HIST("Proton/hSelectionCounter")); hProtonSel->GetXaxis()->SetBinLabel(1, "All"); hProtonSel->GetXaxis()->SetBinLabel(2, "Pt"); hProtonSel->GetXaxis()->SetBinLabel(3, "Eta"); - hProtonSel->GetXaxis()->SetBinLabel(4, "TPC nSigma"); - hProtonSel->GetXaxis()->SetBinLabel(5, "TOF nSigma"); + hProtonSel->GetXaxis()->SetBinLabel(4, "TPC clusters"); + hProtonSel->GetXaxis()->SetBinLabel(5, "TPC nSigma"); + hProtonSel->GetXaxis()->SetBinLabel(6, "TOF nSigma"); + hProtonSel->GetXaxis()->SetBinLabel(7, "DCA to PV"); histos.add("Proton/hPt", "Proton/hPt", kTH1F, {axisProtonPt}); histos.add("Proton/h2TPCNSigmaVsPt", "Proton/h2TPCNSigmaVsPt", kTH2F, {axisProtonPt, axisNSigma}); histos.add("Proton/h2TOFNSigmaVsPt", "Proton/h2TOFNSigmaVsPt", kTH2F, {axisProtonPt, axisNSigma}); + histos.add("Proton/h2TPCNClsVsPt", "Proton/h2TPCNClsVsPt", kTH2F, {axisProtonPt, axisTpcNCls}); + histos.add("Proton/h2DcaToPVVsPt", "Proton/h2DcaToPVVsPt", kTH2F, {axisProtonPt, axisProtonDcaToPV}); - const AxisSpec axisCandSel{7, -0.5, 6.5, "selection step"}; + const AxisSpec axisCandSel{15, -0.5, 14.5, "selection step"}; + const AxisSpec axisCandPhotonDcaToPV{250, 0., 250., "DCA_{#gamma-line} to PV (cm)"}; const AxisSpec axisDca{100, 0., 10., "DCA(p,#gamma) (cm)"}; + const AxisSpec axisDcaToPV{500, 0., 0.5, "DCA_{cand} to PV (cm)"}; const AxisSpec axisCandRadius{200, 0., 200., "R_{dec} (cm)"}; + const AxisSpec axisFlightDistance{250, 0., 250., "|SV-PV| (cm)"}; const AxisSpec axisDiscriminant{200, -1., 1., "discriminant (GeV^{4}/#it{c}^{4})"}; + const AxisSpec axisRootCenter{400, -10., 10., "-b/(2a) (GeV/#it{c})"}; + const AxisSpec axisAntiSigmaPA{180, 0., 3.15, "AntiPA (rad)"}; const AxisSpec axisMassSigma{200, 1.0, 1.4, "m_{p#gamma#gamma} (GeV/#it{c}^{2})"}; + const AxisSpec axisRapidity{200, -2., 2., "y_{#Sigma^{+}}"}; const AxisSpec axisSigmaPt{100, 0., 6., "#it{p}_{T,#Sigma^{+}} (GeV/#it{c})"}; const AxisSpec axisMomentum{100, 0., 10., "#it{p} (GeV/#it{c})"}; histos.add("Candidate/hSelectionCounter", "Candidate/hSelectionCounter", kTH1F, {axisCandSel}); auto hCandSel = histos.get(HIST("Candidate/hSelectionCounter")); hCandSel->GetXaxis()->SetBinLabel(1, "All pairs"); - hCandSel->GetXaxis()->SetBinLabel(2, "Vertex fit"); - hCandSel->GetXaxis()->SetBinLabel(3, "DCA(p,#gamma)"); - hCandSel->GetXaxis()->SetBinLabel(4, "Radius"); - hCandSel->GetXaxis()->SetBinLabel(5, "Real root"); - hCandSel->GetXaxis()->SetBinLabel(6, "Valid root"); - hCandSel->GetXaxis()->SetBinLabel(7, "Filled"); + hCandSel->GetXaxis()->SetBinLabel(2, "Autocorrelation"); + hCandSel->GetXaxis()->SetBinLabel(3, "Vertex fit"); + hCandSel->GetXaxis()->SetBinLabel(4, "DCA(p,#gamma)"); + hCandSel->GetXaxis()->SetBinLabel(5, "Radius"); + hCandSel->GetXaxis()->SetBinLabel(6, "Flight distance"); + hCandSel->GetXaxis()->SetBinLabel(7, "Real root"); + hCandSel->GetXaxis()->SetBinLabel(8, "Valid root"); + hCandSel->GetXaxis()->SetBinLabel(9, "DCA to PV"); + hCandSel->GetXaxis()->SetBinLabel(10, "Mass"); + hCandSel->GetXaxis()->SetBinLabel(11, "Rapidity"); + hCandSel->GetXaxis()->SetBinLabel(12, "Photon opening angle"); + hCandSel->GetXaxis()->SetBinLabel(13, "Photon pointing angle"); + hCandSel->GetXaxis()->SetBinLabel(14, "Photon DCA to PV"); + hCandSel->GetXaxis()->SetBinLabel(15, "Filled"); histos.add("Candidate/hDcaProtonGamma", "Candidate/hDcaProtonGamma", kTH1F, {axisDca}); + histos.add("Candidate/hDcaToPV", "Candidate/hDcaToPV", kTH1F, {axisDcaToPV}); histos.add("Candidate/hRadius", "Candidate/hRadius", kTH1F, {axisCandRadius}); + histos.add("Candidate/hFlightDistance", "Candidate/hFlightDistance", kTH1F, {axisFlightDistance}); + histos.add("Candidate/hPhotonOpeningAngle", "Candidate/hPhotonOpeningAngle", kTH1F, {axisPhotonOpeningAngle}); + histos.add("Candidate/hPhotonPointingAngle", "Candidate/hPhotonPointingAngle", kTH1F, {axisPhotonOpeningAngle}); + histos.add("Candidate/hPhotonDcaToPV", "Candidate/hPhotonDcaToPV", kTH1F, {axisCandPhotonDcaToPV}); histos.add("Candidate/hDiscriminant", "Candidate/hDiscriminant", kTH1F, {axisDiscriminant}); + histos.add("Candidate/hRootCenter", "Candidate/hRootCenter", kTH1F, {axisRootCenter}); histos.add("Candidate/hMassSigmaPlus", "Candidate/hMassSigmaPlus", kTH1F, {axisMassSigma}); + histos.add("Candidate/hRapidity", "Candidate/hRapidity", kTH1F, {axisRapidity}); histos.add("Candidate/h2MassVsPt", "Candidate/h2MassVsPt", kTH2F, {axisSigmaPt, axisMassSigma}); + histos.add("Candidate/h2MassVsRootCenter", "Candidate/h2MassVsRootCenter", kTH2F, {axisRootCenter, axisMassSigma}); + histos.add("Candidate/hAntiSigmaPointingAngle", "Candidate/hAntiSigmaPointingAngle", kTH1F, {axisAntiSigmaPA}); + histos.add("Candidate/h2MassVsAntiSigmaPointingAngle", "Candidate/h2MassVsAntiSigmaPointingAngle", kTH2F, {axisAntiSigmaPA, axisMassSigma}); const AxisSpec axisDiscrIter{discrRetryMaxIter + 2, -0.5, discrRetryMaxIter + 1.5, "discriminant retry iteration"}; histos.add("Candidate/hDiscriminantRetryIter", "Candidate/hDiscriminantRetryIter", kTH1F, {axisDiscrIter}); @@ -225,14 +290,22 @@ struct Sigmaplusbuilder { hPhotonSelSignal->GetXaxis()->SetBinLabel(1, "All"); hPhotonSelSignal->GetXaxis()->SetBinLabel(2, "Neg eta"); hPhotonSelSignal->GetXaxis()->SetBinLabel(3, "Pos eta"); - hPhotonSelSignal->GetXaxis()->SetBinLabel(4, "TPC nSigma_{el}"); - 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, "Mass"); + hPhotonSelSignal->GetXaxis()->SetBinLabel(4, "TPC clusters"); + hPhotonSelSignal->GetXaxis()->SetBinLabel(5, "TPC nSigma_{el}"); + hPhotonSelSignal->GetXaxis()->SetBinLabel(6, "DCA daughters"); + hPhotonSelSignal->GetXaxis()->SetBinLabel(7, "Radius"); + hPhotonSelSignal->GetXaxis()->SetBinLabel(8, "Opening angle"); + hPhotonSelSignal->GetXaxis()->SetBinLabel(9, "Delta theta"); + hPhotonSelSignal->GetXaxis()->SetBinLabel(10, "CosPA"); + hPhotonSelSignal->GetXaxis()->SetBinLabel(11, "Rapidity"); + hPhotonSelSignal->GetXaxis()->SetBinLabel(12, "Qt"); + hPhotonSelSignal->GetXaxis()->SetBinLabel(13, "Alpha"); + hPhotonSelSignal->GetXaxis()->SetBinLabel(14, "Mass"); + + histos.add("Photon/True/h2TPCNSigmaElPosVsPt", "Photon/True/h2TPCNSigmaElPosVsPt", kTH2F, {axisPhotonPt, axisNSigmaEl}); + histos.add("Photon/True/h2TPCNSigmaElNegVsPt", "Photon/True/h2TPCNSigmaElNegVsPt", kTH2F, {axisPhotonPt, axisNSigmaEl}); + histos.add("Photon/True/h2TPCNClsPosVsPt", "Photon/True/h2TPCNClsPosVsPt", kTH2F, {axisPhotonPt, axisTpcNCls}); + histos.add("Photon/True/h2TPCNClsNegVsPt", "Photon/True/h2TPCNClsNegVsPt", kTH2F, {axisPhotonPt, axisTpcNCls}); } if (doprocessMc) { @@ -241,23 +314,50 @@ struct Sigmaplusbuilder { histos.add("Photon/True/hRadius", "Photon/True/hRadius", kTH1F, {axisPhotonRadius}); histos.add("Photon/True/h2ArmenterosPodolanski", "Photon/True/h2ArmenterosPodolanski", kTH2F, {axisAlpha, axisQt}); histos.add("Photon/True/h2ConvPointXY", "Photon/True/h2ConvPointXY", kTH2F, {axisConvXY, axisConvXY}); + histos.add("Photon/True/hOpeningAngle", "Photon/True/hOpeningAngle", kTH1F, {axisPhotonOpeningAngle}); + histos.add("Photon/True/hDeltaTheta", "Photon/True/hDeltaTheta", kTH1F, {axisPhotonDeltaTheta}); histos.add("Proton/True/hPt", "Proton/True/hPt", kTH1F, {axisProtonPt}); histos.add("Proton/True/h2TPCNSigmaVsPt", "Proton/True/h2TPCNSigmaVsPt", kTH2F, {axisProtonPt, axisNSigma}); histos.add("Proton/True/h2TOFNSigmaVsPt", "Proton/True/h2TOFNSigmaVsPt", kTH2F, {axisProtonPt, axisNSigma}); + histos.add("Proton/True/h2TPCNClsVsPt", "Proton/True/h2TPCNClsVsPt", kTH2F, {axisProtonPt, axisTpcNCls}); + histos.add("Proton/True/h2DcaToPVVsPt", "Proton/True/h2DcaToPVVsPt", kTH2F, {axisProtonPt, axisProtonDcaToPV}); + + histos.add("Proton/True/hSelectionCounter", "Proton/True/hSelectionCounter", kTH1F, {axisProtonSel}); + auto hProtonSelSignal = histos.get(HIST("Proton/True/hSelectionCounter")); + hProtonSelSignal->GetXaxis()->SetBinLabel(1, "All"); + hProtonSelSignal->GetXaxis()->SetBinLabel(2, "Pt"); + hProtonSelSignal->GetXaxis()->SetBinLabel(3, "Eta"); + hProtonSelSignal->GetXaxis()->SetBinLabel(4, "TPC clusters"); + hProtonSelSignal->GetXaxis()->SetBinLabel(5, "TPC nSigma"); + hProtonSelSignal->GetXaxis()->SetBinLabel(6, "TOF nSigma"); + hProtonSelSignal->GetXaxis()->SetBinLabel(7, "DCA to PV"); histos.add("Candidate/True/hSelectionCounter", "Candidate/True/hSelectionCounter", kTH1F, {axisCandSel}); auto hCandSelSignal = histos.get(HIST("Candidate/True/hSelectionCounter")); hCandSelSignal->GetXaxis()->SetBinLabel(1, "All pairs"); - hCandSelSignal->GetXaxis()->SetBinLabel(2, "Vertex fit"); - hCandSelSignal->GetXaxis()->SetBinLabel(3, "DCA(p,#gamma)"); - hCandSelSignal->GetXaxis()->SetBinLabel(4, "Radius"); - hCandSelSignal->GetXaxis()->SetBinLabel(5, "Real root"); - hCandSelSignal->GetXaxis()->SetBinLabel(6, "Valid root"); - hCandSelSignal->GetXaxis()->SetBinLabel(7, "Filled"); + hCandSelSignal->GetXaxis()->SetBinLabel(2, "Autocorrelation"); + hCandSelSignal->GetXaxis()->SetBinLabel(3, "Vertex fit"); + hCandSelSignal->GetXaxis()->SetBinLabel(4, "DCA(p,#gamma)"); + hCandSelSignal->GetXaxis()->SetBinLabel(5, "Radius"); + hCandSelSignal->GetXaxis()->SetBinLabel(6, "Flight distance"); + hCandSelSignal->GetXaxis()->SetBinLabel(7, "Real root"); + hCandSelSignal->GetXaxis()->SetBinLabel(8, "Valid root"); + hCandSelSignal->GetXaxis()->SetBinLabel(9, "DCA to PV"); + hCandSelSignal->GetXaxis()->SetBinLabel(10, "Mass"); + hCandSelSignal->GetXaxis()->SetBinLabel(11, "Rapidity"); + hCandSelSignal->GetXaxis()->SetBinLabel(12, "Photon opening angle"); + hCandSelSignal->GetXaxis()->SetBinLabel(13, "Photon pointing angle"); + hCandSelSignal->GetXaxis()->SetBinLabel(14, "Photon DCA to PV"); + hCandSelSignal->GetXaxis()->SetBinLabel(15, "Filled"); histos.add("Candidate/True/hDcaProtonGamma", "Candidate/True/hDcaProtonGamma", kTH1F, {axisDca}); + histos.add("Candidate/True/hDcaToPV", "Candidate/True/hDcaToPV", kTH1F, {axisDcaToPV}); histos.add("Candidate/True/hRadius", "Candidate/True/hRadius", kTH1F, {axisCandRadius}); + histos.add("Candidate/True/hFlightDistance", "Candidate/True/hFlightDistance", kTH1F, {axisFlightDistance}); + histos.add("Candidate/True/hPhotonOpeningAngle", "Candidate/True/hPhotonOpeningAngle", kTH1F, {axisPhotonOpeningAngle}); + histos.add("Candidate/True/hPhotonPointingAngle", "Candidate/True/hPhotonPointingAngle", kTH1F, {axisPhotonOpeningAngle}); + histos.add("Candidate/True/hPhotonDcaToPV", "Candidate/True/hPhotonDcaToPV", kTH1F, {axisCandPhotonDcaToPV}); const AxisSpec axisVtxRes{200, 0., 20., "|vtx_{fit} - vtx_{MC}| (cm)"}; histos.add("Candidate/True/hVertexResFromMcTruth", "Candidate/True/hVertexResFromMcTruth", kTH1F, {axisVtxRes}); const AxisSpec axisMomRes{200, -1., 1., "(p_{fit} - p_{MC}) / p_{MC}"}; @@ -267,8 +367,13 @@ struct Sigmaplusbuilder { histos.add("Candidate/True/hProtonFlightAngle", "Candidate/True/hProtonFlightAngle", kTH1F, {axisProtonFlightAngle}); histos.add("Candidate/True/hDiscriminant", "Candidate/True/hDiscriminant", kTH1F, {axisDiscriminant}); histos.add("Candidate/True/hDiscriminantRetryIter", "Candidate/True/hDiscriminantRetryIter", kTH1F, {axisDiscrIter}); + histos.add("Candidate/True/hRootCenter", "Candidate/True/hRootCenter", kTH1F, {axisRootCenter}); histos.add("Candidate/True/hMassSigmaPlus", "Candidate/True/hMassSigmaPlus", kTH1F, {axisMassSigma}); + histos.add("Candidate/True/hRapidity", "Candidate/True/hRapidity", kTH1F, {axisRapidity}); histos.add("Candidate/True/h2MassVsPt", "Candidate/True/h2MassVsPt", kTH2F, {axisSigmaPt, axisMassSigma}); + histos.add("Candidate/True/h2MassVsRootCenter", "Candidate/True/h2MassVsRootCenter", kTH2F, {axisRootCenter, axisMassSigma}); + histos.add("Candidate/True/hAntiSigmaPointingAngle", "Candidate/True/hAntiSigmaPointingAngle", kTH1F, {axisAntiSigmaPA}); + histos.add("Candidate/True/h2MassVsAntiSigmaPointingAngle", "Candidate/True/h2MassVsAntiSigmaPointingAngle", kTH2F, {axisAntiSigmaPA, axisMassSigma}); histos.add("MC/hGenSigmaPlusPt", "MC/hGenSigmaPlusPt", kTH1F, {axisSigmaPt}); @@ -385,28 +490,76 @@ struct Sigmaplusbuilder { } fillPhotonStep(2); + histos.fill(HIST("Photon/h2TPCNClsPosVsPt"), posTrack.pt(), posTrack.tpcNClsFound()); + histos.fill(HIST("Photon/h2TPCNClsNegVsPt"), negTrack.pt(), negTrack.tpcNClsFound()); + if constexpr (IsMC) { + if (isSignal) { + histos.fill(HIST("Photon/True/h2TPCNClsPosVsPt"), posTrack.pt(), posTrack.tpcNClsFound()); + histos.fill(HIST("Photon/True/h2TPCNClsNegVsPt"), negTrack.pt(), negTrack.tpcNClsFound()); + } + } + if (posTrack.tpcNClsFound() < photonDauMinTpcNCls || negTrack.tpcNClsFound() < photonDauMinTpcNCls) { + continue; + } + fillPhotonStep(3); + histos.fill(HIST("Photon/h2TPCNSigmaElPosVsPt"), posTrack.pt(), posTrack.tpcNSigmaEl()); histos.fill(HIST("Photon/h2TPCNSigmaElNegVsPt"), negTrack.pt(), negTrack.tpcNSigmaEl()); + if constexpr (IsMC) { + if (isSignal) { + histos.fill(HIST("Photon/True/h2TPCNSigmaElPosVsPt"), posTrack.pt(), posTrack.tpcNSigmaEl()); + histos.fill(HIST("Photon/True/h2TPCNSigmaElNegVsPt"), negTrack.pt(), negTrack.tpcNSigmaEl()); + } + } if (posTrack.tpcNSigmaEl() < photonDauMinTPCNSigmaEl || posTrack.tpcNSigmaEl() > photonDauMaxTPCNSigmaEl || negTrack.tpcNSigmaEl() < photonDauMinTPCNSigmaEl || negTrack.tpcNSigmaEl() > photonDauMaxTPCNSigmaEl) { continue; } - fillPhotonStep(3); + fillPhotonStep(4); if (v0.dcaV0daughters() > photonMaxDCAV0Dau) { continue; } - fillPhotonStep(4); + fillPhotonStep(5); std::array secVtx{v0.x(), v0.y(), v0.z()}; float radius = v0.v0radius(); if (radius < photonMinRadius || radius > photonMaxRadius) { continue; } - fillPhotonStep(5); + fillPhotonStep(6); std::array pNeg{v0.pxneg(), v0.pyneg(), v0.pzneg()}; std::array pPos{v0.pxpos(), v0.pypos(), v0.pzpos()}; + + // opening angle between the daughter momenta + float photonOpeningAngleV0 = std::acos(std::clamp(dot3(pPos, pNeg) / std::sqrt(dot3(pPos, pPos) * dot3(pNeg, pNeg)), -1.f, 1.f)); + histos.fill(HIST("Photon/hOpeningAngle"), photonOpeningAngleV0); + if constexpr (IsMC) { + if (isSignal) { + histos.fill(HIST("Photon/True/hOpeningAngle"), photonOpeningAngleV0); + } + } + if (photonOpeningAngleV0 > photonMaxOpeningAngle) { + continue; + } + fillPhotonStep(7); + + // delta theta between the daughter tracks' own polar angles + float posTheta = 2.f * std::atan(std::exp(-posTrack.eta())); + float negTheta = 2.f * std::atan(std::exp(-negTrack.eta())); + float photonDeltaTheta = posTheta - negTheta; + histos.fill(HIST("Photon/hDeltaTheta"), photonDeltaTheta); + if constexpr (IsMC) { + if (isSignal) { + histos.fill(HIST("Photon/True/hDeltaTheta"), photonDeltaTheta); + } + } + if (std::abs(photonDeltaTheta) > photonMaxDeltaTheta) { + continue; + } + fillPhotonStep(8); + std::array pGamma{pNeg[0] + pPos[0], pNeg[1] + pPos[1], pNeg[2] + pPos[2]}; float gammaP = std::sqrt(dot3(pGamma, pGamma)); @@ -416,31 +569,31 @@ struct Sigmaplusbuilder { if (cosPA < photonMinV0cospa) { continue; } - fillPhotonStep(6); + fillPhotonStep(9); float photonY = RecoDecay::y(pGamma, o2::constants::physics::MassGamma); if (photonY < photonMinRapidity || photonY > photonMaxRapidity) { continue; } - fillPhotonStep(7); + fillPhotonStep(10); float qtarm = v0.qtarm(); float alpha = v0.alpha(); if (qtarm > photonMaxQt) { continue; } - fillPhotonStep(8); + fillPhotonStep(11); if (std::abs(alpha) > photonMaxAlpha) { continue; } - fillPhotonStep(9); + fillPhotonStep(12); float mGamma = v0.mGamma(); if (mGamma > photonMaxMass) { continue; } - fillPhotonStep(10); + fillPhotonStep(13); histos.fill(HIST("Photon/hMass"), mGamma); histos.fill(HIST("Photon/hPt"), std::hypot(pGamma[0], pGamma[1])); @@ -455,38 +608,91 @@ struct Sigmaplusbuilder { } // proton candidate selection - template + template bool selectProton(const TTrack& track) { - histos.fill(HIST("Proton/hSelectionCounter"), 0); + bool isSignal = false; + if constexpr (IsMC) { + if (track.has_mcParticle()) { + auto mcProton = track.template mcParticle_as(); + isSignal = findSigmaPlusMotherOfProton(mcProton) >= 0; + } + } + auto fillProtonStep = [&](int step) { + histos.fill(HIST("Proton/hSelectionCounter"), step); + if constexpr (IsMC) { + if (isSignal) { + histos.fill(HIST("Proton/True/hSelectionCounter"), step); + } + } + }; + fillProtonStep(0); if (track.pt() < protonMinPt) { return false; } - histos.fill(HIST("Proton/hSelectionCounter"), 1); + fillProtonStep(1); if (std::abs(track.eta()) > protonMaxEta) { return false; } - histos.fill(HIST("Proton/hSelectionCounter"), 2); + fillProtonStep(2); + + histos.fill(HIST("Proton/h2TPCNClsVsPt"), track.pt(), track.tpcNClsFound()); + if constexpr (IsMC) { + if (isSignal) { + histos.fill(HIST("Proton/True/h2TPCNClsVsPt"), track.pt(), track.tpcNClsFound()); + } + } + if (track.tpcNClsFound() < protonMinTpcNCls) { + return false; + } + fillProtonStep(3); if (std::abs(track.tpcNSigmaPr()) > protonMaxTPCNSigma) { return false; } - histos.fill(HIST("Proton/hSelectionCounter"), 3); + fillProtonStep(4); + // A track with no TOF hit at all is tolerated even above protonPtMinRequireTOF - only an existing-but-failing + // TOF nSigma gets rejected. Set protonRequireTofHit=true to instead mandate a TOF hit above the threshold. if (track.pt() > protonPtMinRequireTOF) { - if (!track.hasTOF() || std::abs(track.tofNSigmaPr()) > protonMaxTOFNSigma) { + if (track.hasTOF()) { + if (std::abs(track.tofNSigmaPr()) > protonMaxTOFNSigma) { + return false; + } + } else if (protonRequireTofHit) { return false; } } - histos.fill(HIST("Proton/hSelectionCounter"), 4); + fillProtonStep(5); + + // proton DCA to PV + histos.fill(HIST("Proton/h2DcaToPVVsPt"), track.pt(), std::abs(track.dcaXY())); + if constexpr (IsMC) { + if (isSignal) { + histos.fill(HIST("Proton/True/h2DcaToPVVsPt"), track.pt(), std::abs(track.dcaXY())); + } + } + if (std::abs(track.dcaXY()) < protonMinDcaToPV || std::abs(track.dcaXY()) > protonMaxDcaToPV) { + return false; + } + fillProtonStep(6); histos.fill(HIST("Proton/hPt"), track.pt()); histos.fill(HIST("Proton/h2TPCNSigmaVsPt"), track.pt(), track.tpcNSigmaPr()); if (track.hasTOF()) { histos.fill(HIST("Proton/h2TOFNSigmaVsPt"), track.pt(), track.tofNSigmaPr()); } + if constexpr (IsMC) { + if (isSignal) { + histos.fill(HIST("Proton/True/hPt"), track.pt()); + histos.fill(HIST("Proton/True/h2TPCNSigmaVsPt"), track.pt(), track.tpcNSigmaPr()); + if (track.hasTOF()) { + histos.fill(HIST("Proton/True/h2TOFNSigmaVsPt"), track.pt(), track.tofNSigmaPr()); + } + } + } return true; } @@ -625,24 +831,23 @@ struct Sigmaplusbuilder { auto negTrack = photon.negTrack; bool isSignal = false; + bool collisionIdCheck = false; // MC only: true if the proton's true MC collision matches the reconstructed collision 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 + float decayRadiusMC = -999.f; // MC-truth decay radius + float massMC = -999.f; // MC-truth Sigma+ mass int matchedSigmaId = -1; - int protonPdgCode = 0; - int protonMotherPdgCode = 0; - int gammaPdgCode = 0; - int gammaMotherPdgCode = 0; - int gammaGMotherPdgCode = 0; if constexpr (IsMC) { if (protonTrack.has_mcParticle()) { auto mcProton = protonTrack.template mcParticle_as(); - protonPdgCode = mcProton.pdgCode(); - auto const& protonMothers = mcProton.template mothers_as(); - if (!protonMothers.empty()) { - protonMotherPdgCode = protonMothers.front().pdgCode(); + + auto protonCollision = protonTrack.template collision_as(); + if (protonCollision.has_mcCollision()) { + collisionIdCheck = protonCollision.mcCollision().globalIndex() == mcProton.mcCollisionId(); } + auto const& protonMothers = mcProton.template mothers_as(); int protonSigmaIdx = findSigmaPlusMotherOfProton(mcProton); @@ -653,19 +858,7 @@ struct Sigmaplusbuilder { auto const& posMothers = mcPos.template mothers_as(); if (!posMothers.empty()) { auto mcGamma = posMothers.front(); - gammaPdgCode = mcGamma.pdgCode(); mcTrueMomGamma = {mcGamma.px(), mcGamma.py(), mcGamma.pz()}; - - auto const& gammaMothers = mcGamma.template mothers_as(); - if (!gammaMothers.empty()) { - auto mcPi0 = gammaMothers.front(); - gammaMotherPdgCode = mcPi0.pdgCode(); - - auto const& pi0Mothers = mcPi0.template mothers_as(); - if (!pi0Mothers.empty()) { - gammaGMotherPdgCode = pi0Mothers.front().pdgCode(); - } - } } int photonSigmaIdx = findSigmaPlusMotherOfPhoton(mcPos, mcNeg); @@ -675,8 +868,11 @@ 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(); + auto mcSigmaPlusMother = protonMothers.front(); + mcTrueMomSigmaPlus = {mcSigmaPlusMother.px(), mcSigmaPlusMother.py(), mcSigmaPlusMother.pz()}; + matchedSigmaId = mcSigmaPlusMother.globalIndex(); + decayRadiusMC = std::hypot(mcTrueVtx[0] - mcSigmaPlusMother.vx(), mcTrueVtx[1] - mcSigmaPlusMother.vy()); + massMC = std::sqrt(mcSigmaPlusMother.e() * mcSigmaPlusMother.e() - mcSigmaPlusMother.p() * mcSigmaPlusMother.p()); } } } @@ -691,7 +887,15 @@ struct Sigmaplusbuilder { }; fillCandStep(0); // all pairs + // Reject the pair if the proton track is one of the photon's own e+/e- daughter tracks. + if (protonTrack.globalIndex() == posTrack.globalIndex() || protonTrack.globalIndex() == negTrack.globalIndex()) { + return -1; + } + fillCandStep(1); // autocorrelation + auto protonTrackParCov = getTrackParCov(protonTrack); + std::array protonOrigPos{}; + protonTrackParCov.getXYZGlo(protonOrigPos); std::array zeroCov{}; auto photonTrackParCov = o2::track::TrackParCov({photon.x, photon.y, photon.z}, {photon.px, photon.py, photon.pz}, zeroCov, 0, true); @@ -707,7 +911,7 @@ struct Sigmaplusbuilder { if (nCand == 0 || !fitter.propagateTracksToVertex()) { return -1; } - fillCandStep(1); // Vertex fit + fillCandStep(2); // Vertex fit float fitChi2 = fitter.getChi2AtPCACandidate(); float dcaProtonGamma = std::sqrt(fitChi2); @@ -720,7 +924,7 @@ struct Sigmaplusbuilder { if (dcaProtonGamma > candMaxDcaProtonGamma) { return -1; } - fillCandStep(2); // DCA(p,gamma) + fillCandStep(3); // DCA(p,gamma) std::array secVtx = fitter.getPCACandidatePos(); float radius = std::hypot(secVtx[0], secVtx[1]); @@ -735,13 +939,24 @@ struct Sigmaplusbuilder { if (radius < candMinRadius || radius > candMaxRadius) { return -1; } - fillCandStep(3); // radius + fillCandStep(4); // radius // flight direction n and the decay-plane basis n, eIn, eOut std::array flightVec{secVtx[0] - pv[0], secVtx[1] - pv[1], secVtx[2] - pv[2]}; std::array nHat = normalize3(flightVec); float flightDistance = std::sqrt(dot3(flightVec, flightVec)); + histos.fill(HIST("Candidate/hFlightDistance"), flightDistance); + if constexpr (IsMC) { + if (isSignal) { + histos.fill(HIST("Candidate/True/hFlightDistance"), flightDistance); + } + } + if (flightDistance < candMinFlightDistance || flightDistance > candMaxFlightDistance) { + return -1; + } + fillCandStep(5); // flight distance + std::array pProton{}; std::array pGamma1{}; fitter.getTrack(0).getPxPyPzGlo(pProton); @@ -832,12 +1047,27 @@ struct Sigmaplusbuilder { if (discriminant < 0.f) { return -1; } - fillCandStep(4); // real root + fillCandStep(6); // real root // two roots from the quadratic, among both we keep the mass closest to the nominal Sigma+ mass float sqrtDisc = std::sqrt(discriminant); std::array roots{(-coefB + sqrtDisc) / (2.f * coefA), (-coefB - sqrtDisc) / (2.f * coefA)}; + // perturbation-free center of the two roots + float rootCenter = -coefB / (2.f * coefA); + histos.fill(HIST("Candidate/hRootCenter"), rootCenter); + if constexpr (IsMC) { + if (isSignal) { + histos.fill(HIST("Candidate/True/hRootCenter"), rootCenter); + } + } + if (candRejectNegRootCenter && rootCenter < 0.f) { + return -1; + } + if (rootCenter > candMaxRootCenter) { + return -1; + } + bool haveCandidate = false; float bestMass = -999.f; std::array bestMomGamma2{}; @@ -864,54 +1094,164 @@ struct Sigmaplusbuilder { if (!haveCandidate) { return -1; } - fillCandStep(5); // valid root + fillCandStep(7); // valid root std::array pSigma{pProton[0] + pGamma1[0] + bestMomGamma2[0], pProton[1] + pGamma1[1] + bestMomGamma2[1], pProton[2] + pGamma1[2] + bestMomGamma2[2]}; float ptSigma = std::hypot(pSigma[0], pSigma[1]); + // candidate DCA to PV + o2::track::TrackPar sigmaTrackPar({secVtx[0], secVtx[1], secVtx[2]}, {pSigma[0], pSigma[1], pSigma[2]}, protonTrack.sign(), true); + std::array dcaSigmaToPv{}; + o2::base::Propagator::Instance()->propagateToDCA(o2::math_utils::Point3D{pv[0], pv[1], pv[2]}, sigmaTrackPar, mBz, 2.f, + o2::base::Propagator::MatCorrType::USEMatCorrNONE, &dcaSigmaToPv); + float candDcaToPV = std::abs(dcaSigmaToPv[0]); + histos.fill(HIST("Candidate/hDcaToPV"), candDcaToPV); + if constexpr (IsMC) { + if (isSignal) { + histos.fill(HIST("Candidate/True/hDcaToPV"), candDcaToPV); + } + } + if (candDcaToPV > candMaxDcaToPV) { + return -1; + } + fillCandStep(8); // DCA to PV + + // AntiSigmaPointingAngle - Fake pointing angle + std::array protonPath{secVtx[0] - protonOrigPos[0], secVtx[1] - protonOrigPos[1], secVtx[2] - protonOrigPos[2]}; + float protonPathLength = std::sqrt(dot3(protonPath, protonPath)); + float svRadiusFromPv = std::hypot(secVtx[0] - pv[0], secVtx[1] - pv[1]); + float protonOrigRadiusFromPv = std::hypot(protonOrigPos[0] - pv[0], protonOrigPos[1] - pv[1]); + float propDir = (svRadiusFromPv < protonOrigRadiusFromPv) ? -1.f : 1.f; + float qProton = (protonTrack.sign() < 0) ? -1.f : 1.f; + float protonP = std::sqrt(dot3(pProton, pProton)); + float rCurve = protonP * 1000.f / (0.2998f * std::abs(mBz)); + float bzSign = (mBz < 0) ? -1.f : 1.f; + float alphaRot = -propDir * qProton * bzSign * o2::constants::math::PI; + if (protonPathLength / (2.f * rCurve) < 1.f) { + alphaRot = -2.f * propDir * qProton * bzSign * std::asin(protonPathLength / (2.f * rCurve)); + } + std::array pProtonRot{ + pProton[0] * std::cos(alphaRot) - pProton[1] * std::sin(alphaRot), + pProton[0] * std::sin(alphaRot) + pProton[1] * std::cos(alphaRot), + pProton[2]}; + float alphaVtxRot = -propDir * qProton * bzSign * o2::constants::math::PI * 0.5f; + if (flightDistance / (2.f * rCurve) < 1.f) { + alphaVtxRot = -propDir * qProton * bzSign * std::asin(flightDistance / (2.f * rCurve)); + } + std::array sigmaVertexVec{secVtx[0] - pv[0], secVtx[1] - pv[1], secVtx[2] - pv[2]}; + std::array sigmaVertexVecRot{ + sigmaVertexVec[0] * std::cos(alphaVtxRot) - sigmaVertexVec[1] * std::sin(alphaVtxRot), + sigmaVertexVec[0] * std::sin(alphaVtxRot) + sigmaVertexVec[1] * std::cos(alphaVtxRot), + sigmaVertexVec[2]}; + float antiSigmaPointingAngle = std::acos(std::clamp(dot3(pProtonRot, sigmaVertexVecRot) / std::sqrt(dot3(pProtonRot, pProtonRot) * dot3(sigmaVertexVecRot, sigmaVertexVecRot)), -1.f, 1.f)); + histos.fill(HIST("Candidate/hAntiSigmaPointingAngle"), antiSigmaPointingAngle); + if constexpr (IsMC) { + if (isSignal) { + histos.fill(HIST("Candidate/True/hAntiSigmaPointingAngle"), antiSigmaPointingAngle); + } + } + // window cut + if (antiSigmaPointingAngle < candMinAntiSigmaPointingAngle || antiSigmaPointingAngle > candMaxAntiSigmaPointingAngle) { + return -1; + } + histos.fill(HIST("Candidate/hMassSigmaPlus"), bestMass); histos.fill(HIST("Candidate/h2MassVsPt"), ptSigma, bestMass); + histos.fill(HIST("Candidate/h2MassVsRootCenter"), rootCenter, bestMass); + histos.fill(HIST("Candidate/h2MassVsAntiSigmaPointingAngle"), antiSigmaPointingAngle, bestMass); if constexpr (IsMC) { if (isSignal) { histos.fill(HIST("Candidate/True/hMassSigmaPlus"), bestMass); histos.fill(HIST("Candidate/True/h2MassVsPt"), ptSigma, bestMass); + histos.fill(HIST("Candidate/True/h2MassVsAntiSigmaPointingAngle"), antiSigmaPointingAngle, bestMass); + histos.fill(HIST("Candidate/True/h2MassVsRootCenter"), rootCenter, bestMass); + } + } + + if (bestMass > candMaxSigmaMass) { + return -1; + } + fillCandStep(9); // mass + + float candRapidity = RecoDecay::y(pSigma, o2::constants::physics::MassSigmaPlus); + histos.fill(HIST("Candidate/hRapidity"), candRapidity); + if constexpr (IsMC) { + if (isSignal) { + histos.fill(HIST("Candidate/True/hRapidity"), candRapidity); } } - fillCandStep(6); // filled + if (std::abs(candRapidity) > candMaxRapidity) { + return -1; + } + fillCandStep(10); // rapidity - // photon (V0) opening angle: angle between the e+/e- daughter momenta at their own reference point + // photon (V0) opening angle, recomputed here from the daughters' own raw track momenta std::array pPosDau{posTrack.px(), posTrack.py(), posTrack.pz()}; std::array pNegDau{negTrack.px(), negTrack.py(), negTrack.pz()}; float photonOpeningAngle = std::acos(std::clamp(dot3(pPosDau, pNegDau) / std::sqrt(dot3(pPosDau, pPosDau) * dot3(pNegDau, pNegDau)), -1.f, 1.f)); + histos.fill(HIST("Candidate/hPhotonOpeningAngle"), photonOpeningAngle); + if constexpr (IsMC) { + if (isSignal) { + histos.fill(HIST("Candidate/True/hPhotonOpeningAngle"), photonOpeningAngle); + } + } + if (photonOpeningAngle > candMaxPhotonOpeningAngle) { + return -1; + } + fillCandStep(11); // photon opening angle - // 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}; - float photonPointingAngle = std::acos(std::clamp(dot3(pGamma1, convToDecVtx) / std::sqrt(dot3(pGamma1, pGamma1) * dot3(convToDecVtx, convToDecVtx)), -1.f, 1.f)); + // photon pointing angle + std::array decVtxToConv{photon.x - secVtx[0], photon.y - secVtx[1], photon.z - secVtx[2]}; + float photonPointingAngle = std::acos(std::clamp(dot3(pGamma1, decVtxToConv) / std::sqrt(dot3(pGamma1, pGamma1) * dot3(decVtxToConv, decVtxToConv)), -1.f, 1.f)); + histos.fill(HIST("Candidate/hPhotonPointingAngle"), photonPointingAngle); + if constexpr (IsMC) { + if (isSignal) { + histos.fill(HIST("Candidate/True/hPhotonPointingAngle"), photonPointingAngle); + } + } + if (photonPointingAngle > candMaxPhotonPointingAngle) { + return -1; + } + fillCandStep(12); // photon pointing angle - // photon DCA to PV: distance from the PV to the line through the conversion point along the photon momentum direction + // photon DCA to PV 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)); + histos.fill(HIST("Candidate/hPhotonDcaToPV"), photonDcaToPV); + if constexpr (IsMC) { + if (isSignal) { + histos.fill(HIST("Candidate/True/hPhotonDcaToPV"), photonDcaToPV); + } + } + if (photonDcaToPV > candMaxPhotonDcaToPV) { + return -1; + } + fillCandStep(13); // photon DCA to PV + fillCandStep(14); // filled if constexpr (IsMC) { if (fillSlimTables) { - slimSigmaPlusCandsMC(radius, dcaProtonGamma, + slimSigmaPlusCandsMC(radius, + candDcaToPV, dcaProtonGamma, + protonTrack.sign(), + protonTrack.dcaXY(), protonTrack.dcaZ(), pProton[0], pProton[1], pProton[2], pGamma1[0], pGamma1[1], pGamma1[2], bestMomGamma2[0], bestMomGamma2[1], bestMomGamma2[2], protonTrack.tpcNSigmaPr(), protonTrack.tofNSigmaPr(), posTrack.tpcNSigmaEl(), negTrack.tpcNSigmaEl(), photon.mGamma, - protonPdgCode, protonMotherPdgCode, - gammaPdgCode, gammaMotherPdgCode, gammaGMotherPdgCode, - std::hypot(mcTrueVtx[0], mcTrueVtx[1]), + collisionIdCheck, + isSignal, + decayRadiusMC, massMC, mcTrueMomSigmaPlus[0], mcTrueMomSigmaPlus[1], mcTrueMomSigmaPlus[2]); return matchedSigmaId; } sigmaPlusCandsMC(secVtx[0], secVtx[1], secVtx[2], - radius, flightDistance, dcaProtonGamma, fitChi2, + flightDistance, dcaProtonGamma, pProton[0], pProton[1], pProton[2], pGamma1[0], pGamma1[1], pGamma1[2], bestMomGamma2[0], bestMomGamma2[1], bestMomGamma2[2], @@ -919,18 +1259,23 @@ struct Sigmaplusbuilder { posTrack.tpcNSigmaEl(), negTrack.tpcNSigmaEl(), photon.mGamma, photon.alpha, photon.qtarm, photon.radius, photonOpeningAngle, photonPointingAngle, photonDcaToPV, + rootCenter, antiSigmaPointingAngle, candDcaToPV, + protonTrack.sign(), protonTrack.itsNCls(), protonTrack.tpcNClsFound(), protonTrack.dcaXY(), protonTrack.dcaZ(), posTrack.itsNCls(), posTrack.tpcNClsFound(), negTrack.itsNCls(), negTrack.tpcNClsFound(), + collisionIdCheck, isSignal, - protonPdgCode, protonMotherPdgCode, - gammaPdgCode, gammaMotherPdgCode, gammaGMotherPdgCode, mcTrueVtx[0], mcTrueVtx[1], mcTrueVtx[2], mcTrueMomProton[0], mcTrueMomProton[1], mcTrueMomProton[2], mcTrueMomGamma[0], mcTrueMomGamma[1], mcTrueMomGamma[2], - mcTrueMomSigmaPlus[0], mcTrueMomSigmaPlus[1], mcTrueMomSigmaPlus[2]); + mcTrueMomSigmaPlus[0], mcTrueMomSigmaPlus[1], mcTrueMomSigmaPlus[2], + decayRadiusMC, massMC); } else { if (fillSlimTables) { - slimSigmaPlusCands(radius, dcaProtonGamma, + slimSigmaPlusCands(radius, + candDcaToPV, dcaProtonGamma, + protonTrack.sign(), + protonTrack.dcaXY(), protonTrack.dcaZ(), pProton[0], pProton[1], pProton[2], pGamma1[0], pGamma1[1], pGamma1[2], bestMomGamma2[0], bestMomGamma2[1], bestMomGamma2[2], @@ -940,7 +1285,7 @@ struct Sigmaplusbuilder { return matchedSigmaId; } sigmaPlusCands(secVtx[0], secVtx[1], secVtx[2], - radius, flightDistance, dcaProtonGamma, fitChi2, + flightDistance, dcaProtonGamma, pProton[0], pProton[1], pProton[2], pGamma1[0], pGamma1[1], pGamma1[2], bestMomGamma2[0], bestMomGamma2[1], bestMomGamma2[2], @@ -948,6 +1293,8 @@ struct Sigmaplusbuilder { posTrack.tpcNSigmaEl(), negTrack.tpcNSigmaEl(), photon.mGamma, photon.alpha, photon.qtarm, photon.radius, photonOpeningAngle, photonPointingAngle, photonDcaToPV, + rootCenter, antiSigmaPointingAngle, candDcaToPV, + protonTrack.sign(), protonTrack.itsNCls(), protonTrack.tpcNClsFound(), protonTrack.dcaXY(), protonTrack.dcaZ(), posTrack.itsNCls(), posTrack.tpcNClsFound(), negTrack.itsNCls(), negTrack.tpcNClsFound()); } @@ -970,6 +1317,9 @@ struct Sigmaplusbuilder { void processData(CollisionsFull const& collisions, aod::V0Datas const& v0s, TracksFull const& tracks, aod::BCs const&) { for (const auto& collision : collisions) { + if (std::abs(collision.posZ()) > cutZVertex || !collision.sel8()) { + continue; + } initCCDB(collision.bc_as()); histos.fill(HIST("hVertexZ"), collision.posZ()); std::array pv{collision.posX(), collision.posY(), collision.posZ()}; @@ -981,7 +1331,7 @@ struct Sigmaplusbuilder { std::vector acceptedProtons; for (const auto& track : tracksThisCollision) { - if (selectProton(track)) { + if (selectProton(track)) { acceptedProtons.push_back(track); } } @@ -995,11 +1345,14 @@ struct Sigmaplusbuilder { } 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, aod::V0Datas const& v0s, TracksFullMC const& tracks, aod::BCs const&, aod::McParticles const& mcParticles, aod::McCollisions const&) { std::vector matchedSigmaPlusMcIds; // Sigma+ MC indices that got at least one signal candidate for (const auto& collision : collisions) { + if (std::abs(collision.posZ()) > cutZVertex || !collision.sel8()) { + continue; + } initCCDB(collision.bc_as()); histos.fill(HIST("hVertexZ"), collision.posZ()); std::array pv{collision.posX(), collision.posY(), collision.posZ()}; @@ -1044,7 +1397,7 @@ struct Sigmaplusbuilder { std::vector acceptedProtons; for (const auto& track : tracksThisCollision) { - if (selectProton(track)) { + if (selectProton(track)) { acceptedProtons.push_back(track); histos.fill(HIST("MC/hProtonTruthQA"), 0); @@ -1053,12 +1406,6 @@ struct Sigmaplusbuilder { auto mcProton = track.template mcParticle_as(); if (findSigmaPlusMotherOfProton(mcProton) >= 0) { histos.fill(HIST("MC/hProtonTruthQA"), 2); - - histos.fill(HIST("Proton/True/hPt"), track.pt()); - histos.fill(HIST("Proton/True/h2TPCNSigmaVsPt"), track.pt(), track.tpcNSigmaPr()); - if (track.hasTOF()) { - histos.fill(HIST("Proton/True/h2TOFNSigmaVsPt"), track.pt(), track.tofNSigmaPr()); - } } } } @@ -1079,6 +1426,12 @@ struct Sigmaplusbuilder { if (!isSigmaPlusToProtonPi0(mcPart)) { continue; } + if (std::abs(mcPart.y()) > cutRapMotherMC) { + continue; + } + if (mcPart.pt() < cutPtGenMC) { + continue; + } histos.fill(HIST("MC/hGenSigmaPlusPt"), mcPart.pt()); bool wasReconstructed = std::find(matchedSigmaPlusMcIds.begin(), matchedSigmaPlusMcIds.end(), mcPart.globalIndex()) != matchedSigmaPlusMcIds.end(); @@ -1098,24 +1451,29 @@ struct Sigmaplusbuilder { break; } } + float genDecayRadiusMC = std::hypot(genDecVtx[0] - mcPart.vx(), genDecVtx[1] - mcPart.vy()); + float genMassMC = std::sqrt(mcPart.e() * mcPart.e() - mcPart.p() * mcPart.p()); if (fillSlimTables) { - slimSigmaPlusCandsMC(-999.f, -999.f, + slimSigmaPlusCandsMC(-999.f, + -999.f, -999.f, + 0, + -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, - pdgProton, mcPart.pdgCode(), - 0, 0, 0, - std::hypot(genDecVtx[0], genDecVtx[1]), + false, + false, + genDecayRadiusMC, genMassMC, mcPart.px(), mcPart.py(), mcPart.pz()); continue; } 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, @@ -1123,15 +1481,17 @@ struct Sigmaplusbuilder { -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, 0, -999, -999.f, -999.f, 0, -999, 0, -999, - true, - pdgProton, mcPart.pdgCode(), - 0, 0, 0, + false, + false, genDecVtx[0], genDecVtx[1], genDecVtx[2], genMomProton[0], genMomProton[1], genMomProton[2], -999.f, -999.f, -999.f, - mcPart.px(), mcPart.py(), mcPart.pz()); + mcPart.px(), mcPart.py(), mcPart.pz(), + genDecayRadiusMC, genMassMC); } } PROCESS_SWITCH(Sigmaplusbuilder, processMc, "Process MC", false); @@ -1140,6 +1500,9 @@ struct Sigmaplusbuilder { { constexpr int MinDauTpcCls = 90; for (const auto& collision : collisions) { + if (std::abs(collision.posZ()) > cutZVertex || !collision.sel8()) { + continue; + } initCCDB(collision.bc_as()); auto tracksThisCollision = tracks.sliceBy(tracksPerCollisionMC, collision.globalIndex()); auto v0sThisCollision = v0s.sliceBy(v0PerCollision, collision.globalIndex());