diff --git a/PWGLF/DataModel/LFKinkDecayTables.h b/PWGLF/DataModel/LFKinkDecayTables.h index 0b9137d64b1..a8a5cffa8a3 100644 --- a/PWGLF/DataModel/LFKinkDecayTables.h +++ b/PWGLF/DataModel/LFKinkDecayTables.h @@ -25,6 +25,7 @@ #include #include +#include namespace o2::aod { @@ -162,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 @@ -193,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) @@ -204,26 +209,33 @@ 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(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(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(CollisionIdCheck, collisionIdCheck, bool); //! True if the proton's collision ID matches the reconstructed collision ID + +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 +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 +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; }); @@ -244,12 +256,17 @@ 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", - 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, @@ -257,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, @@ -268,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, @@ -278,14 +298,41 @@ 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, + sigmapluscand::PtSigmaPlus, + sigmapluscand::MassSigmaPlus, + sigmapluscand::PtSigmaPlusMC, + sigmapluscand::YSigmaPlusMC); + +DECLARE_SOA_TABLE(SlimSigmaPlusCands, "AOD", "SLIMSIGMAPLUS", + 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, // dynamic columns sigmapluscand::PxSigmaPlus, @@ -294,6 +341,32 @@ DECLARE_SOA_TABLE(SigmaPlusCandsMC, "AOD", "SIGMAPLUSMC", sigmapluscand::PtSigmaPlus, sigmapluscand::MassSigmaPlus); +DECLARE_SOA_TABLE(SlimSigmaPlusCandsMC, "AOD", "SLIMSIGMAPLUSMC", + 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::CollisionIdCheck, + sigmapluscand::IsSignal, + sigmapluscand::DecayRadiusMC, sigmapluscand::MassMC, + sigmapluscand::PxSigmaPlusMC, sigmapluscand::PySigmaPlusMC, sigmapluscand::PzSigmaPlusMC, + + // dynamic columns + sigmapluscand::PxSigmaPlus, + sigmapluscand::PySigmaPlus, + sigmapluscand::PzSigmaPlus, + sigmapluscand::PtSigmaPlus, + sigmapluscand::MassSigmaPlus, + sigmapluscand::PtSigmaPlusMC, + sigmapluscand::YSigmaPlusMC); + } // namespace o2::aod #endif // PWGLF_DATAMODEL_LFKINKDECAYTABLES_H_ diff --git a/PWGLF/TableProducer/Strangeness/sigmaHadCorr.cxx b/PWGLF/TableProducer/Strangeness/sigmaHadCorr.cxx index dc7df2f567a..6f2a8a2336c 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,22 +97,24 @@ 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+"}; 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)"}; @@ -119,7 +122,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 +139,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 +225,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; @@ -235,29 +239,29 @@ struct sigmaHadCorrTask { 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) < 1e-6f) { + 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) @@ -309,8 +313,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 +336,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 +408,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 +430,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 +500,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 +529,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 +562,7 @@ struct sigmaHadCorrTask { continue; } - sigmaHadCand candidate; + SigmaHadCand candidate; candidate.sigmaCharge = sigmaCand.mothSign(); candidate.sigmaPx = sigmaCand.pxMoth(); candidate.sigmaPy = sigmaCand.pyMoth(); @@ -575,33 +600,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 +652,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 +663,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 +685,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 +707,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 +787,64 @@ 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; + } + 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; + 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 +883,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 +931,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..7969042b519 100644 --- a/PWGLF/TableProducer/Strangeness/sigmaplusbuilder.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaplusbuilder.cxx @@ -18,19 +18,25 @@ #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 #include #include +#include #include #include #include +#include +#include +#include #include #include #include @@ -52,36 +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 photonMaxTPCNSigmaEl{"photonMaxTPCNSigmaEl", 15, "Max |TPC nSigma_el| for photon daughters"}; + 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 @@ -93,13 +124,17 @@ struct Sigmaplusbuilder { Configurable discrRetryPhiPar0{"discrRetryPhiPar0", 0.000129845, "par0 of the delta-phi resolution function"}; Configurable discrRetryPhiPar1{"discrRetryPhiPar1", 0.00199688, "par1 of the delta-phi resolution function"}; + Configurable fillSlimTables{"fillSlimTables", false, "write the slim candidate tables instead of the full ones"}; + Configurable ccdbPath{"ccdbPath", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; Configurable grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; Produces sigmaPlusCands; Produces sigmaPlusCandsMC; + Produces slimSigmaPlusCands; + Produces slimSigmaPlusCandsMC; - Service ccdb; + Service ccdb{}; o2::vertexing::DCAFitterN<2> fitter; int mRunNumber = 0; float mBz = 0; @@ -123,6 +158,7 @@ struct Sigmaplusbuilder { fitter.setMinParamChange(1e-3); fitter.setMinRelChi2Change(0.9); fitter.setMaxDZIni(1e9); + fitter.setMaxDXYIni(1e9); fitter.setMaxChi2(1e9); fitter.setUseAbsDCA(true); @@ -133,116 +169,195 @@ 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)"}; const AxisSpec axisAlpha{100, -1., 1., "#alpha_{AP}"}; 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}); 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(2, "Neg eta"); + hPhotonSel->GetXaxis()->SetBinLabel(3, "Pos eta"); + 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, "CosPA"); - hPhotonSel->GetXaxis()->SetBinLabel(9, "Qt"); - hPhotonSel->GetXaxis()->SetBinLabel(10, "Alpha"); - hPhotonSel->GetXaxis()->SetBinLabel(11, "TPC nSigma el"); + 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}); histos.add("Photon/hRadius", "Photon/hRadius", kTH1F, {axisPhotonRadius}); histos.add("Photon/h2ArmenterosPodolanski", "Photon/h2ArmenterosPodolanski", kTH2F, {axisAlpha, axisQt}); 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}); - if (doprocessMc) { + if (doprocessMc || doprocessFindable) { 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(2, "Neg eta"); + hPhotonSelSignal->GetXaxis()->SetBinLabel(3, "Pos eta"); + 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, "CosPA"); - hPhotonSelSignal->GetXaxis()->SetBinLabel(9, "Qt"); - hPhotonSelSignal->GetXaxis()->SetBinLabel(10, "Alpha"); - hPhotonSelSignal->GetXaxis()->SetBinLabel(11, "TPC nSigma el"); + 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) { histos.add("Photon/True/hMass", "Photon/True/hMass", kTH1F, {axisPhotonMass}); histos.add("Photon/True/hPt", "Photon/True/hPt", kTH1F, {axisPhotonPt}); 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}"}; @@ -252,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}); @@ -277,11 +397,11 @@ struct Sigmaplusbuilder { if (doprocessFindable) { const AxisSpec axisDetectorPresence{3, -0.5, 2.5, "detector"}; const AxisSpec axisDuplicateTrack{2, -0.5, 1.5, "track"}; - const AxisSpec axisV0Presence{2, -0.5, 1.5, "V0 match"}; const AxisSpec axisTPCClusters{160, -0.5, 159.5, "TPC clusters"}; const AxisSpec axisPhotonMomResolution{200, -1., 1., "(#it{p}_{reco,#gamma} - #it{p}_{MC,#gamma})/#it{p}_{MC,#gamma}"}; const AxisSpec axisPhotonCosPA{200, -1., 1., "cosPA_{#gamma}"}; const AxisSpec axisPairMass{200, 0., 0.1, "m_{e^{+}e^{-}} (GeV/#it{c}^{2})"}; + const AxisSpec axisV0Presence{2, -0.5, 1.5, "V0 match"}; histos.add("Findable/hSigmaPlusPt", "Findable/hSigmaPlusPt", kTH1F, {axisSigmaPt}); histos.add("Findable/h2ProtonPtVsSigmaPlusPt", "Findable/h2ProtonPtVsSigmaPlusPt", kTH2F, {axisSigmaPt, axisMomentum}); histos.add("Findable/hElectronPt", "Findable/hElectronPt", kTH1F, {axisMomentum}); @@ -291,14 +411,19 @@ 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}); histos.add("Findable/hElectronDetectorPresence", "Findable/hElectronDetectorPresence", kTH1F, {axisDetectorPresence}); histos.add("Findable/hPositronDetectorPresence", "Findable/hPositronDetectorPresence", kTH1F, {axisDetectorPresence}); + 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, "passed photon selection"); auto hDuplicateConversionTrackCounter = histos.get(HIST("Findable/hDuplicateConversionTrackCounter")); hDuplicateConversionTrackCounter->GetXaxis()->SetBinLabel(1, "e^{-}"); hDuplicateConversionTrackCounter->GetXaxis()->SetBinLabel(2, "e^{+}"); @@ -313,125 +438,261 @@ 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 candidates + template + std::vector> findPhotonsFromV0s(const TV0s& v0s, const TTracks&, 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& v0 : v0s) { + auto posTrack = v0.template posTrack_as(); + auto negTrack = v0.template negTrack_as(); + + 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; + } } - } - auto fillPhotonStep = [&](int step) { - histos.fill(HIST("Photon/hSelectionCounter"), step); + auto fillPhotonStep = [&](int step) { + histos.fill(HIST("Photon/hSelectionCounter"), step); + if constexpr (IsMC) { + if (isSignal) { + histos.fill(HIST("Photon/True/hSelectionCounter"), step); + } + } + }; + fillPhotonStep(0); + + if (negTrack.eta() < photonDauEtaMin || negTrack.eta() > photonDauEtaMax) { + continue; + } + fillPhotonStep(1); + + if (posTrack.eta() < photonDauEtaMin || posTrack.eta() > photonDauEtaMax) { + continue; + } + 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/hSelectionCounter"), step); + histos.fill(HIST("Photon/True/h2TPCNClsPosVsPt"), posTrack.pt(), posTrack.tpcNClsFound()); + histos.fill(HIST("Photon/True/h2TPCNClsNegVsPt"), negTrack.pt(), negTrack.tpcNClsFound()); } } - }; - fillPhotonStep(0); + if (posTrack.tpcNClsFound() < photonDauMinTpcNCls || negTrack.tpcNClsFound() < photonDauMinTpcNCls) { + continue; + } + fillPhotonStep(3); - if (v0.mGamma() < 0 || v0.mGamma() > photonMaxMass) { - return false; - } - fillPhotonStep(1); + 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(4); - 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 (v0.dcaV0daughters() > photonMaxDCAV0Dau) { + continue; + } + fillPhotonStep(5); - if (v0.negativeeta() < photonDauEtaMin || v0.negativeeta() > photonDauEtaMax) { - return false; - } - fillPhotonStep(3); + std::array secVtx{v0.x(), v0.y(), v0.z()}; + float radius = v0.v0radius(); + if (radius < photonMinRadius || radius > photonMaxRadius) { + continue; + } + fillPhotonStep(6); - if (v0.positiveeta() < photonDauEtaMin || v0.positiveeta() > photonDauEtaMax) { - return false; - } - fillPhotonStep(4); + std::array pNeg{v0.pxneg(), v0.pyneg(), v0.pzneg()}; + std::array pPos{v0.pxpos(), v0.pypos(), v0.pzpos()}; - if (std::abs(v0.dcaV0daughters()) > photonMaxDCAV0Dau) { - return false; - } - fillPhotonStep(5); + // 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); - if (v0.v0radius() < photonMinRadius || v0.v0radius() > photonMaxRadius) { - return false; - } - fillPhotonStep(6); + // 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); - if (v0.v0cosPA() < photonMinV0cospa) { - return false; - } - fillPhotonStep(7); + std::array pGamma{pNeg[0] + pPos[0], pNeg[1] + pPos[1], pNeg[2] + pPos[2]}; + float gammaP = std::sqrt(dot3(pGamma, pGamma)); - if (v0.qtarm() > photonMaxQt) { - return false; - } - fillPhotonStep(8); + 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(9); - if (std::abs(v0.alpha()) > photonMaxAlpha) { - return false; - } - fillPhotonStep(9); + float photonY = RecoDecay::y(pGamma, o2::constants::physics::MassGamma); + if (photonY < photonMinRapidity || photonY > photonMaxRapidity) { + continue; + } + fillPhotonStep(10); - if (std::abs(posTrack.tpcNSigmaEl()) > photonMaxTPCNSigmaEl || std::abs(negTrack.tpcNSigmaEl()) > photonMaxTPCNSigmaEl) { - return false; - } - fillPhotonStep(10); + float qtarm = v0.qtarm(); + float alpha = v0.alpha(); + if (qtarm > photonMaxQt) { + continue; + } + fillPhotonStep(11); - 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()); + if (std::abs(alpha) > photonMaxAlpha) { + continue; + } + fillPhotonStep(12); - return true; + float mGamma = v0.mGamma(); + if (mGamma > photonMaxMass) { + continue; + } + fillPhotonStep(13); + + 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]); + + photons.push_back({secVtx[0], secVtx[1], secVtx[2], pGamma[0], pGamma[1], pGamma[2], mGamma, alpha, qtarm, radius, negTrack, posTrack}); + } + + return photons; } // 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; } @@ -562,29 +823,31 @@ 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 - int protonPdgCode = 0; - int protonMotherPdgCode = 0; - int gammaPdgCode = 0; - int gammaMotherPdgCode = 0; - int gammaGMotherPdgCode = 0; + 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; 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); @@ -595,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); @@ -617,6 +868,11 @@ struct Sigmaplusbuilder { if (isSignal) { mcTrueVtx = {mcProton.vx(), mcProton.vy(), mcProton.vz()}; mcTrueMomProton = {mcProton.px(), mcProton.py(), mcProton.pz()}; + 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()); } } } @@ -631,10 +887,18 @@ 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); + 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,12 +906,12 @@ struct Sigmaplusbuilder { try { nCand = fitter.process(protonTrackParCov, photonTrackParCov); } catch (...) { - return; + return -1; } if (nCand == 0 || !fitter.propagateTracksToVertex()) { - return; + return -1; } - fillCandStep(1); // Vertex fit + fillCandStep(2); // Vertex fit float fitChi2 = fitter.getChi2AtPCACandidate(); float dcaProtonGamma = std::sqrt(fitChi2); @@ -658,9 +922,9 @@ struct Sigmaplusbuilder { } } if (dcaProtonGamma > candMaxDcaProtonGamma) { - return; + 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]); @@ -673,17 +937,28 @@ struct Sigmaplusbuilder { } } if (radius < candMinRadius || radius > candMaxRadius) { - return; + 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)); - std::array pProton; - std::array pGamma1; + 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); fitter.getTrack(1).getPxPyPzGlo(pGamma1); @@ -770,14 +1045,29 @@ struct Sigmaplusbuilder { } } if (discriminant < 0.f) { - return; + 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{}; @@ -802,70 +1092,213 @@ struct Sigmaplusbuilder { } } if (!haveCandidate) { - return; + 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); } } - fillCandStep(6); // filled - // photon (V0) opening angle: angle between the e+/e- daughter momenta at their own reference point + 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); + } + } + if (std::abs(candRapidity) > candMaxRapidity) { + return -1; + } + fillCandStep(10); // rapidity + + // 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 - std::array convPoint{photon.x(), photon.y(), photon.z()}; - std::array photonDir = normalize3({photon.px(), photon.py(), photon.pz()}); + // 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, + 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, + 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], 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, + 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]); + mcTrueMomGamma[0], mcTrueMomGamma[1], mcTrueMomGamma[2], + mcTrueMomSigmaPlus[0], mcTrueMomSigmaPlus[1], mcTrueMomSigmaPlus[2], + decayRadiusMC, massMC); } else { + if (fillSlimTables) { + 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], + protonTrack.tpcNSigmaPr(), protonTrack.tofNSigmaPr(), + posTrack.tpcNSigmaEl(), negTrack.tpcNSigmaEl(), + photon.mGamma); + 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], 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, + rootCenter, antiSigmaPointingAngle, candDcaToPV, + protonTrack.sign(), 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) @@ -884,75 +1317,77 @@ 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()}; + auto tracksThisCollision = tracks.sliceBy(tracksPerCollision, collision.globalIndex()); 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 = findPhotonsFromV0s(v0sThisCollision, tracksThisCollision, pv); + std::vector acceptedProtons; for (const auto& track : tracksThisCollision) { - if (selectProton(track)) { + if (selectProton(track)) { acceptedProtons.push_back(track); } } 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, 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()}; + auto tracksThisCollision = tracks.sliceBy(tracksPerCollisionMC, collision.globalIndex()); 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 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 acceptedPhotons = findPhotonsFromV0s(v0sThisCollision, 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"), 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,10 +1395,9 @@ struct Sigmaplusbuilder { } } - auto tracksThisCollision = tracks.sliceBy(tracksPerCollisionMC, collision.globalIndex()); std::vector acceptedProtons; for (const auto& track : tracksThisCollision) { - if (selectProton(track)) { + if (selectProton(track)) { acceptedProtons.push_back(track); histos.fill(HIST("MC/hProtonTruthQA"), 0); @@ -972,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()); - } } } } @@ -985,25 +1413,101 @@ 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; + } + 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(); + 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; + } + } + 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, -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, + 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, + -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, + 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(), + genDecayRadiusMC, genMassMC); } } PROCESS_SWITCH(Sigmaplusbuilder, processMc, "Process MC", false); - void processFindable(CollisionsFullMC const& collisions, aod::V0Datas const& v0s, TracksFullMC const& tracks, aod::McParticles const&) + void processFindable(CollisionsFullMC const& collisions, aod::V0Datas const& v0s, TracksFullMC const& tracks, aod::McParticles const&, aod::BCs const&) { + 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()); + std::array pv{collision.posX(), collision.posY(), collision.posZ()}; + auto acceptedPhotons = findPhotonsFromV0s(v0sThisCollision, tracksThisCollision, pv); for (const auto& protonTrack : tracksThisCollision) { if (!protonTrack.has_mcParticle()) { @@ -1035,7 +1539,7 @@ struct Sigmaplusbuilder { continue; } - if (electronTrack.tpcNClsFound() < 90 || electronTrack.sign() > 0) { + if (electronTrack.tpcNClsFound() < MinDauTpcCls || electronTrack.sign() > 0) { continue; } @@ -1055,7 +1559,7 @@ struct Sigmaplusbuilder { continue; } - if (positronTrack.tpcNClsFound() < 90 || positronTrack.sign() < 0) { + if (positronTrack.tpcNClsFound() < MinDauTpcCls || positronTrack.sign() < 0) { continue; } @@ -1121,10 +1625,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 +1659,15 @@ 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); + 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]));