diff --git a/PWGLF/DataModel/LFKinkDecayTables.h b/PWGLF/DataModel/LFKinkDecayTables.h index a8a5cffa8a3..db5e3825d43 100644 --- a/PWGLF/DataModel/LFKinkDecayTables.h +++ b/PWGLF/DataModel/LFKinkDecayTables.h @@ -191,11 +191,11 @@ DECLARE_SOA_COLUMN(PhotonQt, photonQt, float); //! Armente DECLARE_SOA_COLUMN(PhotonConvRadius, photonConvRadius, float); //! Conversion radius of the measured photon (cm) DECLARE_SOA_COLUMN(PhotonOpeningAngle, photonOpeningAngle, float); //! Opening angle between the photon's e+e- daughters (rad) 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(FlightDirTilt, flightDirTilt, float); //! Tilt of the PV->decay-vertex direction needed for the missing photon to have a real solution (rad), 0 if none needed 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 @@ -208,10 +208,26 @@ DECLARE_SOA_COLUMN(PhotonPosTpcNCls, photonPosTpcNCls, int16_t); //! Number of f DECLARE_SOA_COLUMN(PhotonNegItsNCls, photonNegItsNCls, uint8_t); //! Number of ITS clusters of the photon's negative daughter DECLARE_SOA_COLUMN(PhotonNegTpcNCls, photonNegTpcNCls, int16_t); //! Number of found TPC clusters of the photon's negative daughter +DECLARE_SOA_COLUMN(PhotonDcaDau, photonDcaDau, float); //! DCA between the photon's e+e- daughters at the conversion point (cm) +DECLARE_SOA_COLUMN(PhotonCosPAToPV, photonCosPAToPV, float); //! Cosine of the angle between the photon momentum and the PV->conversion-point direction +DECLARE_SOA_COLUMN(PhotonDcaXYToPV, photonDcaXYToPV, float); //! Transverse DCA of the photon's flight line to the PV (cm) +DECLARE_SOA_COLUMN(PhotonDcaZToPV, photonDcaZToPV, float); //! Longitudinal distance of the photon's flight line to the PV at its transverse DCA (cm) +DECLARE_SOA_COLUMN(PhotonPsiPair, photonPsiPair, float); //! psi_pair of the e+e- legs, evaluated after propagating them outward, 999 if the propagation failed +DECLARE_SOA_COLUMN(PhotonPosDcaXY, photonPosDcaXY, float); //! DCAxy of the photon's positive daughter to the PV of its own collision (cm) +DECLARE_SOA_COLUMN(PhotonPosDcaZ, photonPosDcaZ, float); //! DCAz of the photon's positive daughter to the PV of its own collision (cm) +DECLARE_SOA_COLUMN(PhotonNegDcaXY, photonNegDcaXY, float); //! DCAxy of the photon's negative daughter to the PV of its own collision (cm) +DECLARE_SOA_COLUMN(PhotonNegDcaZ, photonNegDcaZ, float); //! DCAz of the photon's negative daughter to the PV of its own collision (cm) +DECLARE_SOA_COLUMN(PhotonPosTpcNClsFindable, photonPosTpcNClsFindable, int16_t); //! Number of findable TPC clusters of the photon's positive daughter +DECLARE_SOA_COLUMN(PhotonNegTpcNClsFindable, photonNegTpcNClsFindable, int16_t); //! Number of findable TPC clusters of the photon's negative daughter +DECLARE_SOA_COLUMN(PhotonPosTpcChi2NCl, photonPosTpcChi2NCl, float); //! TPC chi2 per cluster of the photon's positive daughter +DECLARE_SOA_COLUMN(PhotonNegTpcChi2NCl, photonNegTpcChi2NCl, float); //! TPC chi2 per cluster of the photon's negative daughter + // MC columns 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(IsSignal, isSignal, bool); //! True if the proton and photon share the same true Sigma+ mother +DECLARE_SOA_COLUMN(IsProtonFromSigma, isProtonFromSigma, bool); //! True if the proton is the daughter of a true Sigma+ -> p pi0 +DECLARE_SOA_COLUMN(IsPhotonFromSigma, isPhotonFromSigma, bool); //! True if the photon is from the pi0 of a true Sigma+ -> p pi0 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) @@ -273,11 +289,14 @@ DECLARE_SOA_TABLE(SigmaPlusCands, "AOD", "SIGMAPLUSCANDS", sigmapluscand::NSigmaTPCProton, sigmapluscand::NSigmaTOFProton, 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::PhotonOpeningAngle, sigmapluscand::PhotonPointingAngle, + sigmapluscand::RootCenter, sigmapluscand::AntiSigmaPointingAngle, sigmapluscand::CandDcaToPV, sigmapluscand::FlightDirTilt, sigmapluscand::ProtonSign, sigmapluscand::ProtonItsNCls, sigmapluscand::ProtonTpcNCls, sigmapluscand::ProtonDcaXY, sigmapluscand::ProtonDcaZ, sigmapluscand::PhotonPosItsNCls, sigmapluscand::PhotonPosTpcNCls, sigmapluscand::PhotonNegItsNCls, sigmapluscand::PhotonNegTpcNCls, + sigmapluscand::PhotonDcaDau, sigmapluscand::PhotonCosPAToPV, sigmapluscand::PhotonDcaXYToPV, sigmapluscand::PhotonDcaZToPV, sigmapluscand::PhotonPsiPair, + sigmapluscand::PhotonPosDcaXY, sigmapluscand::PhotonPosDcaZ, sigmapluscand::PhotonNegDcaXY, sigmapluscand::PhotonNegDcaZ, + sigmapluscand::PhotonPosTpcNClsFindable, sigmapluscand::PhotonNegTpcNClsFindable, sigmapluscand::PhotonPosTpcChi2NCl, sigmapluscand::PhotonNegTpcChi2NCl, // dynamic columns sigmapluscand::Radius, @@ -297,13 +316,16 @@ DECLARE_SOA_TABLE(SigmaPlusCandsMC, "AOD", "SIGMAPLUSMC", sigmapluscand::NSigmaTPCProton, sigmapluscand::NSigmaTOFProton, 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::PhotonOpeningAngle, sigmapluscand::PhotonPointingAngle, + sigmapluscand::RootCenter, sigmapluscand::AntiSigmaPointingAngle, sigmapluscand::CandDcaToPV, sigmapluscand::FlightDirTilt, sigmapluscand::ProtonSign, sigmapluscand::ProtonItsNCls, sigmapluscand::ProtonTpcNCls, sigmapluscand::ProtonDcaXY, sigmapluscand::ProtonDcaZ, sigmapluscand::PhotonPosItsNCls, sigmapluscand::PhotonPosTpcNCls, sigmapluscand::PhotonNegItsNCls, sigmapluscand::PhotonNegTpcNCls, + sigmapluscand::PhotonDcaDau, sigmapluscand::PhotonCosPAToPV, sigmapluscand::PhotonDcaXYToPV, sigmapluscand::PhotonDcaZToPV, sigmapluscand::PhotonPsiPair, + sigmapluscand::PhotonPosDcaXY, sigmapluscand::PhotonPosDcaZ, sigmapluscand::PhotonNegDcaXY, sigmapluscand::PhotonNegDcaZ, + sigmapluscand::PhotonPosTpcNClsFindable, sigmapluscand::PhotonNegTpcNClsFindable, sigmapluscand::PhotonPosTpcChi2NCl, sigmapluscand::PhotonNegTpcChi2NCl, sigmapluscand::CollisionIdCheck, - sigmapluscand::IsSignal, + sigmapluscand::IsSignal, sigmapluscand::IsProtonFromSigma, sigmapluscand::IsPhotonFromSigma, sigmapluscand::XDecVtxMC, sigmapluscand::YDecVtxMC, sigmapluscand::ZDecVtxMC, sigmapluscand::PxProtonMC, sigmapluscand::PyProtonMC, sigmapluscand::PzProtonMC, sigmapluscand::PxGammaMC, sigmapluscand::PyGammaMC, sigmapluscand::PzGammaMC, @@ -323,7 +345,7 @@ DECLARE_SOA_TABLE(SigmaPlusCandsMC, "AOD", "SIGMAPLUSMC", DECLARE_SOA_TABLE(SlimSigmaPlusCands, "AOD", "SLIMSIGMAPLUS", sigmapluscand::TransDecayRadius, - sigmapluscand::CandDcaToPV, + sigmapluscand::CandDcaToPV, sigmapluscand::FlightDirTilt, sigmapluscand::DcaProtonGamma, sigmapluscand::ProtonSign, sigmapluscand::ProtonDcaXY, sigmapluscand::ProtonDcaZ, @@ -333,6 +355,9 @@ DECLARE_SOA_TABLE(SlimSigmaPlusCands, "AOD", "SLIMSIGMAPLUS", sigmapluscand::NSigmaTPCProton, sigmapluscand::NSigmaTOFProton, sigmapluscand::NSigmaTPCElPos, sigmapluscand::NSigmaTPCElNeg, sigmapluscand::PhotonMass, + sigmapluscand::PhotonDcaDau, sigmapluscand::PhotonCosPAToPV, sigmapluscand::PhotonDcaXYToPV, sigmapluscand::PhotonDcaZToPV, sigmapluscand::PhotonPsiPair, + sigmapluscand::PhotonPosDcaXY, sigmapluscand::PhotonPosDcaZ, sigmapluscand::PhotonNegDcaXY, sigmapluscand::PhotonNegDcaZ, + sigmapluscand::PhotonPosTpcNClsFindable, sigmapluscand::PhotonNegTpcNClsFindable, sigmapluscand::PhotonPosTpcChi2NCl, sigmapluscand::PhotonNegTpcChi2NCl, // dynamic columns sigmapluscand::PxSigmaPlus, @@ -343,7 +368,7 @@ DECLARE_SOA_TABLE(SlimSigmaPlusCands, "AOD", "SLIMSIGMAPLUS", DECLARE_SOA_TABLE(SlimSigmaPlusCandsMC, "AOD", "SLIMSIGMAPLUSMC", sigmapluscand::TransDecayRadius, - sigmapluscand::CandDcaToPV, + sigmapluscand::CandDcaToPV, sigmapluscand::FlightDirTilt, sigmapluscand::DcaProtonGamma, sigmapluscand::ProtonSign, sigmapluscand::ProtonDcaXY, sigmapluscand::ProtonDcaZ, @@ -353,8 +378,11 @@ DECLARE_SOA_TABLE(SlimSigmaPlusCandsMC, "AOD", "SLIMSIGMAPLUSMC", sigmapluscand::NSigmaTPCProton, sigmapluscand::NSigmaTOFProton, sigmapluscand::NSigmaTPCElPos, sigmapluscand::NSigmaTPCElNeg, sigmapluscand::PhotonMass, + sigmapluscand::PhotonDcaDau, sigmapluscand::PhotonCosPAToPV, sigmapluscand::PhotonDcaXYToPV, sigmapluscand::PhotonDcaZToPV, sigmapluscand::PhotonPsiPair, + sigmapluscand::PhotonPosDcaXY, sigmapluscand::PhotonPosDcaZ, sigmapluscand::PhotonNegDcaXY, sigmapluscand::PhotonNegDcaZ, + sigmapluscand::PhotonPosTpcNClsFindable, sigmapluscand::PhotonNegTpcNClsFindable, sigmapluscand::PhotonPosTpcChi2NCl, sigmapluscand::PhotonNegTpcChi2NCl, sigmapluscand::CollisionIdCheck, - sigmapluscand::IsSignal, + sigmapluscand::IsSignal, sigmapluscand::IsProtonFromSigma, sigmapluscand::IsPhotonFromSigma, sigmapluscand::DecayRadiusMC, sigmapluscand::MassMC, sigmapluscand::PxSigmaPlusMC, sigmapluscand::PySigmaPlusMC, sigmapluscand::PzSigmaPlusMC, diff --git a/PWGLF/TableProducer/Strangeness/sigmaplusbuilder.cxx b/PWGLF/TableProducer/Strangeness/sigmaplusbuilder.cxx index 7969042b519..e9df7b7e324 100644 --- a/PWGLF/TableProducer/Strangeness/sigmaplusbuilder.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaplusbuilder.cxx @@ -13,10 +13,14 @@ /// \brief Task for Sigma+ -> p + pi0 reconstruction, pi0 reconstructed via one converted photon (PCM) /// \author Henrik Fribert (TUM) +#include "PWGEM/Dilepton/Utils/PairUtilities.h" +#include "PWGEM/PhotonMeson/Utils/PCMUtilities.h" #include "PWGLF/DataModel/LFKinkDecayTables.h" #include "PWGLF/DataModel/LFStrangenessTables.h" +#include "PWGLF/Utils/svPoolCreator.h" #include "Common/Core/RecoDecay.h" +#include "Common/Core/TPCVDriftManager.h" #include "Common/Core/trackUtilities.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/PIDResponseTOF.h" @@ -24,32 +28,42 @@ #include "Common/DataModel/TrackSelectionTables.h" #include +#include #include #include #include #include +#include #include #include #include #include #include +#include #include #include #include #include #include +#include #include #include -#include +#include +#include +#include #include #include -#include #include #include #include +#include +#include +#include #include +#include +#include #include using namespace o2; @@ -67,27 +81,33 @@ struct Sigmaplusbuilder { 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 photonMaxMass{"photonMaxMass", 0.05, "Max photon mass (GeV/c^2)"}; + Configurable photonMinRapidity{"photonMinRapidity", -1.0, "Min photon rapidity"}; + Configurable photonMaxRapidity{"photonMaxRapidity", 1.0, "Max photon rapidity"}; + Configurable photonDauEtaMin{"photonDauEtaMin", -1.0, "Min eta of photon daughter tracks"}; + Configurable photonDauEtaMax{"photonDauEtaMax", 1.0, "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 photonMinV0cospa{"photonMinV0cospa", 0.9, "Min photon cosine of pointing angle to the PV"}; + Configurable photonMaxDCAV0Dau{"photonMaxDCAV0Dau", 1.5, "Max DCA between photon daughters (cm)"}; Configurable photonMaxOpeningAngle{"photonMaxOpeningAngle", 0.4, "Max opening angle between the photon's e+/e- daughter momenta (rad)"}; Configurable photonMaxDeltaTheta{"photonMaxDeltaTheta", 0.15, "Max |theta_pos - theta_neg| of the photon's daughter tracks (rad)"}; - Configurable photonMaxQt{"photonMaxQt", 0.15, "Max Armenteros qT for photons (GeV/c)"}; - Configurable photonMaxAlpha{"photonMaxAlpha", 1.0, "Max |Armenteros alpha| for photons"}; - Configurable photonDauMinTPCNSigmaEl{"photonDauMinTPCNSigmaEl", -5., "Min TPC nSigma_el of the photon daughters"}; - Configurable photonDauMaxTPCNSigmaEl{"photonDauMaxTPCNSigmaEl", 5., "Max TPC nSigma_el of the photon daughters"}; - Configurable photonDauMinTpcNCls{"photonDauMinTpcNCls", 30, "Min number of found TPC clusters for the photon (V0) daughter tracks"}; + Configurable photonMaxQt{"photonMaxQt", 0.05, "Max Armenteros qT for photons (GeV/c)"}; + Configurable photonMaxAlpha{"photonMaxAlpha", 0.95, "Max |Armenteros alpha| for photons"}; + Configurable photonDauMinTPCNSigmaEl{"photonDauMinTPCNSigmaEl", -3., "Min TPC nSigma_el of the photon daughters"}; + Configurable photonDauMaxTPCNSigmaEl{"photonDauMaxTPCNSigmaEl", 3., "Max TPC nSigma_el of the photon daughters"}; + Configurable photonDauMaxPt{"photonDauMaxPt", 0.5, "Max pT of the photon daughters (GeV/c)"}; + Configurable photonDauMinTpcNCls{"photonDauMinTpcNCls", 30, "Min number of found TPC clusters for the photon daughter tracks"}; + + // photon daughters paired by this task instead of taken from V0Datas (V0Datas assumes photons from primary vertex) + Configurable useCustomVertexer{"useCustomVertexer", false, "Pair the photon's e+/e- daughter tracks in this task instead of reading V0Datas"}; + Configurable photonSkipAmbiTracks{"photonSkipAmbiTracks", false, "Skip ambiguous tracks when pairing photon daughters"}; + Configurable photonPoolTimeMarginNS{"photonPoolTimeMarginNS", 800., "Time margin (ns) added to a daughter track's time range when matching it to collisions"}; + Configurable photonMaxDXYIni{"photonMaxDXYIni", 4., "Max xy distance (cm) between the two daughter tracks at the start of the photon vertex fit"}; + Configurable photonMaxCircleTouchDist{"photonMaxCircleTouchDist", 4., "Max |circle centre distance - (R1 + R2)| (cm) of the two daughter tracks, applied before the fit"}; // proton selection - Configurable protonMinPt{"protonMinPt", 0.3, "Minimum proton pT (GeV/c)"}; + Configurable protonMinPt{"protonMinPt", 0.4, "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"}; @@ -98,36 +118,39 @@ struct Sigmaplusbuilder { Configurable protonMaxDcaToPV{"protonMaxDcaToPV", 5.0, "Max DCAxy of the proton track to the PV (cm)"}; // proton-photon candidate selection - 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 candVertexProtonWeight{"candVertexProtonWeight", 1.0, "Sigma+ vertex between the proton's (1) and the photon's (0) point of closest approach"}; + Configurable candMaxDcaProtonGamma{"candMaxDcaProtonGamma", 1.5, "Max DCA between proton and photon at the fitted vertex (cm)"}; + Configurable candMaxDcaToPV{"candMaxDcaToPV", 20., "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 candMaxAntiSigmaPointingAngle{"candMaxAntiSigmaPointingAngle", 10., "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 candMaxPhotonOpeningAngle{"candMaxPhotonOpeningAngle", 3.15, "Max photon opening angle (rad), recomputed from the daughters' track momenta"}; + Configurable candMaxPhotonPointingAngle{"candMaxPhotonPointingAngle", 0.5, "Max angle between the photon 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)"}; + Configurable candDeduplicatePhotons{"candDeduplicatePhotons", true, "Per timeframe, write only the candidate with the smallest DCA to PV among candidates sharing the same photon"}; - // missing-photon discriminant retry - // resolution-shaped function and tuned parameters from Run 2 used - Configurable discrRetryMaxIter{"discrRetryMaxIter", 10, "Max attempts to recover a negative discriminant by perturbing the flight direction"}; - Configurable discrRetryThetaRange{"discrRetryThetaRange", 0.02, "Max |delta-theta| perturbation of the flight direction per retry (rad)"}; - Configurable discrRetryThetaPar0{"discrRetryThetaPar0", 0.000133299, "par0 of the delta-theta resolution function"}; - Configurable discrRetryThetaPar1{"discrRetryThetaPar1", 0.0016761, "par1 of the delta-theta resolution function"}; - Configurable discrRetryPhiRange{"discrRetryPhiRange", 0.02, "Max |delta-phi| perturbation of the flight direction per retry (rad)"}; - Configurable discrRetryPhiPar0{"discrRetryPhiPar0", 0.000129845, "par0 of the delta-phi resolution function"}; - Configurable discrRetryPhiPar1{"discrRetryPhiPar1", 0.00199688, "par1 of the delta-phi resolution function"}; + // tilt search: if the measured flight direction gives no real solution for the missing photon, + // tilted flight directions are tried starting with the smallest tilt + Configurable candMaxTilt{"candMaxTilt", 0.3, "Max tilt of the flight direction searched for a real root (rad)"}; + Configurable candTiltStep{"candTiltStep", 0.005, "Tilt step between the rings of tried flight directions (rad)"}; + Configurable candTiltNAzimuth{"candTiltNAzimuth", 36, "Flight directions tried per ring"}; + + // MC + 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 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"}; + Configurable lutPath{"lutPath", "GLO/Param/MatLUT", "CCDB path of the material budget LUT"}; Produces sigmaPlusCands; Produces sigmaPlusCandsMC; @@ -138,8 +161,79 @@ struct Sigmaplusbuilder { o2::vertexing::DCAFitterN<2> fitter; int mRunNumber = 0; float mBz = 0; - TF1 mThetaResoFunc; - TF1 mPhiResoFunc; + o2::base::MatLayerCylSet* mLut = nullptr; + o2::aod::common::TPCVDriftManager mVDriftMgr; + + // daughter pairing (electrons and positrons with the range of collisions each is compatible with) + svPoolCreator svPhotonPoolCreator{PDG_t::kElectron, PDG_t::kPositron}; + std::vector mElectronPool; + std::vector mPositronPool; + std::vector mGoodCollision; // collisions passing the event selection + + // values of a candidate, kept until the end of the timeframe and then written to the tables + struct SigmaPlusCandidate { + uint64_t photonId = 0; + bool isSignal = false; + int photonLegsWithoutIts = 0; + int matchedSigmaId = -1; + std::array decVtx{}; + float radius = 0.f; + float flightDistance = 0.f; + float dcaProtonGamma = 0.f; + float dcaToPV = 0.f; + float tiltAngle = 0.f; + float rootCenter = 0.f; + float antiSigmaPointingAngle = 0.f; + std::array pProton{}; + std::array pGamma1{}; + std::array pGamma2{}; + float protonTpcNSigma = 0.f; + float protonTofNSigma = 0.f; + int protonSign = 0; + uint8_t protonItsNCls = 0; + int16_t protonTpcNCls = 0; + float protonDcaXY = 0.f; + float protonDcaZ = 0.f; + float photonMass = 0.f; + float photonAlpha = 0.f; + float photonQt = 0.f; + float photonRadius = 0.f; + float photonOpeningAngle = 0.f; + float photonPointingAngle = 0.f; + float photonDcaDau = 0.f; + float photonCosPAToPV = 0.f; + float photonDcaXYToPV = 0.f; + float photonDcaZToPV = 0.f; + float photonPsiPair = 0.f; + float posTpcNSigmaEl = 0.f; + float negTpcNSigmaEl = 0.f; + uint8_t posItsNCls = 0; + uint8_t negItsNCls = 0; + int16_t posTpcNCls = 0; + int16_t negTpcNCls = 0; + int16_t posTpcNClsFindable = 0; + int16_t negTpcNClsFindable = 0; + float posTpcChi2NCl = 0.f; + float negTpcChi2NCl = 0.f; + float posDcaXY = 0.f; + float posDcaZ = 0.f; + float negDcaXY = 0.f; + float negDcaZ = 0.f; + // MC + bool collisionIdCheck = false; + bool protonIsSignal = false; + bool photonIsSignal = false; + std::array decVtxMC{}; + std::array pProtonMC{}; + std::array pGammaMC{}; + std::array pSigmaPlusMC{}; + float decayRadiusMC = -999.f; + float massMC = -999.f; + }; + std::vector mCandidatesOfTimeframe; + + // MC: generated Sigma+ -> p pi0 in acceptance of current timeframe and whether a true candidate of it was written + std::unordered_map mGenSigmaWritten; HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; @@ -162,10 +256,11 @@ struct Sigmaplusbuilder { fitter.setMaxChi2(1e9); fitter.setUseAbsDCA(true); - mThetaResoFunc = TF1("thetaResoFunc", "[0]/(abs(x)+[1])", -discrRetryThetaRange, discrRetryThetaRange); - mThetaResoFunc.SetParameters(discrRetryThetaPar0.value, discrRetryThetaPar1.value); - mPhiResoFunc = TF1("phiResoFunc", "[0]/(abs(x)+[1])", -discrRetryPhiRange, discrRetryPhiRange); - mPhiResoFunc.SetParameters(discrRetryPhiPar0.value, discrRetryPhiPar1.value); + mVDriftMgr.init(&ccdb->instance()); + svPhotonPoolCreator.setTimeMargin(photonPoolTimeMarginNS); + if (photonSkipAmbiTracks) { + svPhotonPoolCreator.setSkipAmbiTracks(); + } const AxisSpec axisVertexZ{100, -15., 15., "vrtx_{Z} (cm)"}; @@ -180,221 +275,102 @@ struct Sigmaplusbuilder { 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 axisPhotonDcaV0Dau{200, 0., 10., "DCA between photon daughters (cm)"}; + const AxisSpec axisCollisionsPerPhoton{100, 0.5, 100.5, "collisions a photon is kept in"}; + const AxisSpec axisTpcTimeRangeNColl{200, 0.5, 200.5, "collisions in the TPC time range of a TPC-only daughter"}; 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, "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, "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 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{15, -0.5, 14.5, "selection step"}; + const AxisSpec axisPhotonLegType{3, -0.5, 2.5, "photon daughters without ITS"}; 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 axisTilt{200, 0., 0.5, "flight direction tilt for a real root (rad)"}; 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, "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}); + const AxisSpec axisCandidatesPerPhoton{50, 0.5, 50.5, "candidates sharing one photon"}; + + const AxisSpec axisSigmaSel{6, -0.5, 5.5, "step"}; + + const std::vector photonSteps{"All", "Neg eta", "Pos eta", "TPC clusters", "TPC nSigma_{el}, p_{T}", "DCA daughters", "Radius", + "Opening angle", "Delta theta", "CosPA", "Rapidity", "Qt", "Alpha", "Mass"}; + const std::vector protonSteps{"All", "Pt", "Eta", "TPC clusters", "TPC nSigma", "TOF nSigma", "DCA to PV"}; + const std::vector candSteps{"All pairs", "Autocorrelation", "Vertex fit", "DCA(p,#gamma)", "Radius", "Flight distance", "Real root", + "Valid root", "DCA to PV", "Mass", "Rapidity", "Photon opening angle", "Photon pointing angle", "Photon DCA to PV", "Written"}; + const std::vector sigmaSteps{"Generated", "Proton reconstructed", "Proton selected", "Photon daughters reconstructed", "Photon selected", "Candidate written"}; + auto setBinLabels = [](TAxis* axis, const std::vector& labels) { + for (size_t i = 0; i < labels.size(); ++i) { + axis->SetBinLabel(i + 1, labels[i].c_str()); + } + }; + histos.add("hVertexZ", "hVertexZ", kTH1F, {axisVertexZ}); + + setBinLabels(histos.add("Photon/Inclusive/hSelectionCounter", "hSelectionCounter", kTH1D, {axisPhotonSel})->GetXaxis(), photonSteps); + histos.add("Photon/Inclusive/h2TPCNClsPosVsPt", "h2TPCNClsPosVsPt", kTH2F, {axisPhotonPt, axisTpcNCls}); + histos.add("Photon/Inclusive/h2TPCNClsNegVsPt", "h2TPCNClsNegVsPt", kTH2F, {axisPhotonPt, axisTpcNCls}); + histos.add("Photon/Inclusive/h2TPCNSigmaElPosVsPt", "h2TPCNSigmaElPosVsPt", kTH2F, {axisPhotonPt, axisNSigmaEl}); + histos.add("Photon/Inclusive/h2TPCNSigmaElNegVsPt", "h2TPCNSigmaElNegVsPt", kTH2F, {axisPhotonPt, axisNSigmaEl}); + histos.add("Photon/Inclusive/hDcaV0Daughters", "hDcaV0Daughters", kTH1F, {axisPhotonDcaV0Dau}); + histos.add("Photon/Inclusive/hOpeningAngle", "hOpeningAngle", kTH1F, {axisPhotonOpeningAngle}); + histos.add("Photon/Inclusive/hDeltaTheta", "hDeltaTheta", kTH1F, {axisPhotonDeltaTheta}); + histos.add("Photon/Inclusive/hMass", "hMass", kTH1F, {axisPhotonMass}); + histos.add("Photon/Inclusive/hPt", "hPt", kTH1F, {axisPhotonPt}); + histos.add("Photon/Inclusive/hRadius", "hRadius", kTH1F, {axisPhotonRadius}); + histos.add("Photon/Inclusive/h2ArmenterosPodolanski", "h2ArmenterosPodolanski", kTH2F, {axisAlpha, axisQt}); + histos.add("Photon/Inclusive/h2ConvPointXY", "h2ConvPointXY", kTH2F, {axisConvXY, axisConvXY}); + histos.add("Photon/Inclusive/hCollisionsPerPhoton", "hCollisionsPerPhoton", kTH1F, {axisCollisionsPerPhoton}); + + setBinLabels(histos.add("Proton/Inclusive/hSelectionCounter", "hSelectionCounter", kTH1D, {axisProtonSel})->GetXaxis(), protonSteps); + histos.add("Proton/Inclusive/hPt", "hPt", kTH1F, {axisProtonPt}); + histos.add("Proton/Inclusive/h2TPCNSigmaVsPt", "h2TPCNSigmaVsPt", kTH2F, {axisProtonPt, axisNSigma}); + histos.add("Proton/Inclusive/h2TOFNSigmaVsPt", "h2TOFNSigmaVsPt", kTH2F, {axisProtonPt, axisNSigma}); + histos.add("Proton/Inclusive/h2TPCNClsVsPt", "h2TPCNClsVsPt", kTH2F, {axisProtonPt, axisTpcNCls}); + histos.add("Proton/Inclusive/h2DcaToPVVsPt", "h2DcaToPVVsPt", kTH2F, {axisProtonPt, axisProtonDcaToPV}); + + setBinLabels(histos.add("Candidate/Inclusive/hSelectionCounter", "hSelectionCounter", kTH1D, {axisCandSel})->GetXaxis(), candSteps); + setBinLabels(histos.add("Candidate/Inclusive/h2SelectionCounterVsLegType", "h2SelectionCounterVsLegType", kTH2D, {axisCandSel, axisPhotonLegType})->GetXaxis(), candSteps); + histos.add("Candidate/Inclusive/hDcaProtonGamma", "hDcaProtonGamma", kTH1F, {axisDca}); + histos.add("Candidate/Inclusive/hRadius", "hRadius", kTH1F, {axisCandRadius}); + histos.add("Candidate/Inclusive/hFlightDistance", "hFlightDistance", kTH1F, {axisFlightDistance}); + histos.add("Candidate/Inclusive/hTiltAngle", "hTiltAngle", kTH1F, {axisTilt}); + histos.add("Candidate/Inclusive/hRootCenter", "hRootCenter", kTH1F, {axisRootCenter}); + histos.add("Candidate/Inclusive/hDcaToPV", "hDcaToPV", kTH1F, {axisDcaToPV}); + histos.add("Candidate/Inclusive/hAntiSigmaPointingAngle", "hAntiSigmaPointingAngle", kTH1F, {axisAntiSigmaPA}); + histos.add("Candidate/Inclusive/hMassSigmaPlus", "hMassSigmaPlus", kTH1F, {axisMassSigma}); + histos.add("Candidate/Inclusive/h2MassVsPt", "h2MassVsPt", kTH2F, {axisSigmaPt, axisMassSigma}); + histos.add("Candidate/Inclusive/h2MassVsTilt", "h2MassVsTilt", kTH2F, {axisTilt, axisMassSigma}); + histos.add("Candidate/Inclusive/hRapidity", "hRapidity", kTH1F, {axisRapidity}); + histos.add("Candidate/Inclusive/hPhotonOpeningAngle", "hPhotonOpeningAngle", kTH1F, {axisPhotonOpeningAngle}); + histos.add("Candidate/Inclusive/hPhotonPointingAngle", "hPhotonPointingAngle", kTH1F, {axisPhotonOpeningAngle}); + histos.add("Candidate/Inclusive/hPhotonDcaToPV", "hPhotonDcaToPV", kTH1F, {axisCandPhotonDcaToPV}); 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, "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, "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}); + histos.addClone("Photon/Inclusive/", "Photon/True/"); } - 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, "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}"}; - histos.add("Candidate/True/hProtonMomResFromMcTruth", "Candidate/True/hProtonMomResFromMcTruth", kTH1F, {axisMomRes}); - histos.add("Candidate/True/hPhotonMomResFromMcTruth", "Candidate/True/hPhotonMomResFromMcTruth", kTH1F, {axisMomRes}); - const AxisSpec axisProtonFlightAngle{180, 0., 3.15, "angle(p_{proton}, n) (rad)"}; - 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}); + histos.addClone("Proton/Inclusive/", "Proton/True/"); + histos.addClone("Candidate/Inclusive/", "Candidate/True/"); + } - const AxisSpec axisProtonTruthQA{3, -0.5, 2.5, "step"}; - histos.add("MC/hProtonTruthQA", "MC/hProtonTruthQA", kTH1F, {axisProtonTruthQA}); - auto hProtonTruthQA = histos.get(HIST("MC/hProtonTruthQA")); - hProtonTruthQA->GetXaxis()->SetBinLabel(1, "Accepted"); - hProtonTruthQA->GetXaxis()->SetBinLabel(2, "Has mcParticle"); - hProtonTruthQA->GetXaxis()->SetBinLabel(3, "Sigma+ mother"); + histos.add("Photon/Inclusive/hTpcTimeRangeNColl", "hTpcTimeRangeNColl", kTH1F, {axisTpcTimeRangeNColl}); + histos.add("Candidate/Inclusive/hCandidatesPerPhoton", "hCandidatesPerPhoton", kTH1F, {axisCandidatesPerPhoton}); - const AxisSpec axisPhotonTruthQA{5, -0.5, 4.5, "step"}; - histos.add("MC/hPhotonTruthQA", "MC/hPhotonTruthQA", kTH1F, {axisPhotonTruthQA}); - auto hPhotonTruthQA = histos.get(HIST("MC/hPhotonTruthQA")); - hPhotonTruthQA->GetXaxis()->SetBinLabel(1, "Accepted"); - hPhotonTruthQA->GetXaxis()->SetBinLabel(2, "Both legs have mcParticle"); - hPhotonTruthQA->GetXaxis()->SetBinLabel(3, "Shared real gamma mother"); - hPhotonTruthQA->GetXaxis()->SetBinLabel(4, "Gamma's mother is pi0"); - hPhotonTruthQA->GetXaxis()->SetBinLabel(5, "Pi0's mother is Sigma+"); + if (doprocessMc) { + histos.add("MC/hGenSigmaPlusPt", "MC/hGenSigmaPlusPt", kTH1F, {axisSigmaPt}); + setBinLabels(histos.add("MC/hSigmaPlusCounter", "MC/hSigmaPlusCounter", kTH1D, {axisSigmaSel})->GetXaxis(), sigmaSteps); } if (doprocessFindable) { + const AxisSpec axisMomentum{100, 0., 10., "#it{p} (GeV/#it{c})"}; const AxisSpec axisDetectorPresence{3, -0.5, 2.5, "detector"}; const AxisSpec axisDuplicateTrack{2, -0.5, 1.5, "track"}; const AxisSpec axisTPCClusters{160, -0.5, 159.5, "TPC clusters"}; @@ -418,23 +394,11 @@ struct Sigmaplusbuilder { 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^{+}"); - auto hElectronDetectorPresence = histos.get(HIST("Findable/hElectronDetectorPresence")); - auto hPositronDetectorPresence = histos.get(HIST("Findable/hPositronDetectorPresence")); - hElectronDetectorPresence->GetXaxis()->SetBinLabel(1, "ITS"); - hElectronDetectorPresence->GetXaxis()->SetBinLabel(2, "TPC"); - hElectronDetectorPresence->GetXaxis()->SetBinLabel(3, "TOF"); - hPositronDetectorPresence->GetXaxis()->SetBinLabel(1, "ITS"); - hPositronDetectorPresence->GetXaxis()->SetBinLabel(2, "TPC"); - hPositronDetectorPresence->GetXaxis()->SetBinLabel(3, "TOF"); + setBinLabels(histos.get(HIST("Findable/hConversionPairV0Presence"))->GetXaxis(), {"valid pair", "in V0"}); + setBinLabels(histos.get(HIST("Findable/hPhotonSearchPresence"))->GetXaxis(), {"valid pair", "passed photon selection"}); + setBinLabels(histos.get(HIST("Findable/hDuplicateConversionTrackCounter"))->GetXaxis(), {"e^{-}", "e^{+}"}); + setBinLabels(histos.get(HIST("Findable/hElectronDetectorPresence"))->GetXaxis(), {"ITS", "TPC", "TOF"}); + setBinLabels(histos.get(HIST("Findable/hPositronDetectorPresence"))->GetXaxis(), {"ITS", "TPC", "TOF"}); } } @@ -447,11 +411,48 @@ struct Sigmaplusbuilder { float alpha = 0.f; float qtarm = 0.f; float radius = 0.f; + float dcaDau = 0.f; + float cosPAToPV = 0.f; TTrack negTrack; TTrack posTrack; }; - // photon candidates + // MC counter of generated Sigma+ -> p pi0 in acceptance + void fillSigmaCounter(int sigmaId, int step) + { + if (mGenSigmaWritten.contains(sigmaId)) { + histos.fill(HIST("MC/hSigmaPlusCounter"), step); + } + } + + template + void fillPhotonHist(const TName& name, bool isSignal, Ts... values) + { + histos.fill(HIST("Photon/Inclusive/") + name, values...); + if (isSignal) { + histos.fill(HIST("Photon/True/") + name, values...); + } + } + + template + void fillProtonHist(const TName& name, bool isSignal, Ts... values) + { + histos.fill(HIST("Proton/Inclusive/") + name, values...); + if (isSignal) { + histos.fill(HIST("Proton/True/") + name, values...); + } + } + + template + void fillCandHist(const TName& name, bool isSignal, Ts... values) + { + histos.fill(HIST("Candidate/Inclusive/") + name, values...); + if (isSignal) { + histos.fill(HIST("Candidate/True/") + name, values...); + } + } + + // photon candidates from V0Datas template std::vector> findPhotonsFromV0s(const TV0s& v0s, const TTracks&, const std::array& pv) { @@ -461,22 +462,16 @@ struct Sigmaplusbuilder { auto posTrack = v0.template posTrack_as(); auto negTrack = v0.template negTrack_as(); - bool isSignal = false; + int sigmaId = -1; 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; + sigmaId = findSigmaPlusMotherOfPhoton(posTrack.template mcParticle_as(), negTrack.template mcParticle_as()); } } + bool isSignal = sigmaId >= 0; auto fillPhotonStep = [&](int step) { - histos.fill(HIST("Photon/hSelectionCounter"), step); - if constexpr (IsMC) { - if (isSignal) { - histos.fill(HIST("Photon/True/hSelectionCounter"), step); - } - } + fillPhotonHist(HIST("hSelectionCounter"), isSignal, step); }; fillPhotonStep(0); @@ -490,33 +485,25 @@ struct Sigmaplusbuilder { } fillPhotonStep(2); - histos.fill(HIST("Photon/h2TPCNClsPosVsPt"), posTrack.pt(), posTrack.tpcNClsFound()); - histos.fill(HIST("Photon/h2TPCNClsNegVsPt"), negTrack.pt(), negTrack.tpcNClsFound()); - if constexpr (IsMC) { - if (isSignal) { - histos.fill(HIST("Photon/True/h2TPCNClsPosVsPt"), posTrack.pt(), posTrack.tpcNClsFound()); - histos.fill(HIST("Photon/True/h2TPCNClsNegVsPt"), negTrack.pt(), negTrack.tpcNClsFound()); - } - } + fillPhotonHist(HIST("h2TPCNClsPosVsPt"), isSignal, posTrack.pt(), posTrack.tpcNClsFound()); + fillPhotonHist(HIST("h2TPCNClsNegVsPt"), isSignal, negTrack.pt(), negTrack.tpcNClsFound()); if (posTrack.tpcNClsFound() < photonDauMinTpcNCls || negTrack.tpcNClsFound() < photonDauMinTpcNCls) { continue; } fillPhotonStep(3); - histos.fill(HIST("Photon/h2TPCNSigmaElPosVsPt"), posTrack.pt(), posTrack.tpcNSigmaEl()); - histos.fill(HIST("Photon/h2TPCNSigmaElNegVsPt"), negTrack.pt(), negTrack.tpcNSigmaEl()); - if constexpr (IsMC) { - if (isSignal) { - histos.fill(HIST("Photon/True/h2TPCNSigmaElPosVsPt"), posTrack.pt(), posTrack.tpcNSigmaEl()); - histos.fill(HIST("Photon/True/h2TPCNSigmaElNegVsPt"), negTrack.pt(), negTrack.tpcNSigmaEl()); - } - } + fillPhotonHist(HIST("h2TPCNSigmaElPosVsPt"), isSignal, posTrack.pt(), posTrack.tpcNSigmaEl()); + fillPhotonHist(HIST("h2TPCNSigmaElNegVsPt"), isSignal, negTrack.pt(), negTrack.tpcNSigmaEl()); if (posTrack.tpcNSigmaEl() < photonDauMinTPCNSigmaEl || posTrack.tpcNSigmaEl() > photonDauMaxTPCNSigmaEl || negTrack.tpcNSigmaEl() < photonDauMinTPCNSigmaEl || negTrack.tpcNSigmaEl() > photonDauMaxTPCNSigmaEl) { continue; } + if (posTrack.pt() > photonDauMaxPt || negTrack.pt() > photonDauMaxPt) { + continue; + } fillPhotonStep(4); + fillPhotonHist(HIST("hDcaV0Daughters"), isSignal, v0.dcaV0daughters()); if (v0.dcaV0daughters() > photonMaxDCAV0Dau) { continue; } @@ -532,40 +519,23 @@ struct Sigmaplusbuilder { std::array pNeg{v0.pxneg(), v0.pyneg(), v0.pzneg()}; std::array pPos{v0.pxpos(), v0.pypos(), v0.pzpos()}; - // opening angle between the daughter momenta - float photonOpeningAngleV0 = std::acos(std::clamp(dot3(pPos, pNeg) / std::sqrt(dot3(pPos, pPos) * dot3(pNeg, pNeg)), -1.f, 1.f)); - histos.fill(HIST("Photon/hOpeningAngle"), photonOpeningAngleV0); - if constexpr (IsMC) { - if (isSignal) { - histos.fill(HIST("Photon/True/hOpeningAngle"), photonOpeningAngleV0); - } - } - if (photonOpeningAngleV0 > photonMaxOpeningAngle) { + float photonOpeningAngle = o2::aod::pwgem::dilepton::utils::pairutil::getOpeningAngle(pPos[0], pPos[1], pPos[2], pNeg[0], pNeg[1], pNeg[2]); + fillPhotonHist(HIST("hOpeningAngle"), isSignal, photonOpeningAngle); + if (photonOpeningAngle > photonMaxOpeningAngle) { continue; } fillPhotonStep(7); - // delta theta between the daughter tracks' own polar angles - float posTheta = 2.f * std::atan(std::exp(-posTrack.eta())); - float negTheta = 2.f * std::atan(std::exp(-negTrack.eta())); - float photonDeltaTheta = posTheta - negTheta; - histos.fill(HIST("Photon/hDeltaTheta"), photonDeltaTheta); - if constexpr (IsMC) { - if (isSignal) { - histos.fill(HIST("Photon/True/hDeltaTheta"), photonDeltaTheta); - } - } + // difference of the daughter tracks' polar angles + float photonDeltaTheta = 2.f * std::atan(std::exp(-posTrack.eta())) - 2.f * std::atan(std::exp(-negTrack.eta())); + fillPhotonHist(HIST("hDeltaTheta"), isSignal, photonDeltaTheta); if (std::abs(photonDeltaTheta) > photonMaxDeltaTheta) { continue; } fillPhotonStep(8); std::array pGamma{pNeg[0] + pPos[0], pNeg[1] + pPos[1], pNeg[2] + pPos[2]}; - float gammaP = std::sqrt(dot3(pGamma, pGamma)); - - 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); + float cosPA = RecoDecay::cpa(pv, secVtx, pGamma); if (cosPA < photonMinV0cospa) { continue; } @@ -594,19 +564,428 @@ struct Sigmaplusbuilder { continue; } fillPhotonStep(13); + fillSigmaCounter(sigmaId, 4); - 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]); + fillPhotonHist(HIST("hMass"), isSignal, mGamma); + fillPhotonHist(HIST("hPt"), isSignal, std::hypot(pGamma[0], pGamma[1])); + fillPhotonHist(HIST("hRadius"), isSignal, radius); + fillPhotonHist(HIST("h2ArmenterosPodolanski"), isSignal, alpha, qtarm); + fillPhotonHist(HIST("h2ConvPointXY"), isSignal, 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}); + photons.push_back({secVtx[0], secVtx[1], secVtx[2], pGamma[0], pGamma[1], pGamma[2], mGamma, alpha, qtarm, radius, v0.dcaV0daughters(), cosPA, negTrack, posTrack}); } return photons; } + static uint64_t photonId(int64_t posTrackId, int64_t negTrackId) + { + return (static_cast(posTrackId) << 32) | static_cast(negTrackId); + } + + // psi_pair of the photon daughters (similar to PsiPair in PWGEM) + template + float photonPsiPair(const TTrack& posTrack, const TTrack& negTrack, float convRadius) + { + for (const float& offsetR : {60.f, 30.f, 10.f}) { + auto posTrackPar = getTrackParCov(posTrack); + auto negTrackPar = getTrackParCov(negTrack); + posTrackPar.setPID(o2::track::PID::Electron); + negTrackPar.setPID(o2::track::PID::Electron); + if (!o2::base::Propagator::Instance()->propagateToR(posTrackPar, convRadius + offsetR) || !o2::base::Propagator::Instance()->propagateToR(negTrackPar, convRadius + offsetR)) { + continue; + } + std::array pPos{}; + std::array pNeg{}; + posTrackPar.getPxPyPzGlo(pPos); + negTrackPar.getPxPyPzGlo(pNeg); + return o2::aod::pwgem::dilepton::utils::pairutil::getPsiPair(pPos[0], pPos[1], pPos[2], pNeg[0], pNeg[1], pNeg[2]); + } + return 999.f; + } + + template + void markGoodCollisions(const TCollisions& collisions) + { + mGoodCollision.assign(collisions.size(), false); + for (const auto& collision : collisions) { + if (std::abs(collision.posZ()) > cutZVertex || !collision.sel8()) { + continue; + } + mGoodCollision[collision.globalIndex()] = true; + } + } + + // photon daughter selection of the self-pairing (TPC histograms filled before cuts) + template + bool selectPhotonDaughter(const TTrack& track, bool isSignalLeg) + { + if (track.eta() < photonDauEtaMin || track.eta() > photonDauEtaMax) { + return false; + } + if (track.sign() > 0) { + fillPhotonHist(HIST("h2TPCNClsPosVsPt"), isSignalLeg, track.pt(), track.tpcNClsFound()); + } else { + fillPhotonHist(HIST("h2TPCNClsNegVsPt"), isSignalLeg, track.pt(), track.tpcNClsFound()); + } + if (track.tpcNClsFound() < photonDauMinTpcNCls) { + return false; + } + if (track.sign() > 0) { + fillPhotonHist(HIST("h2TPCNSigmaElPosVsPt"), isSignalLeg, track.pt(), track.tpcNSigmaEl()); + } else { + fillPhotonHist(HIST("h2TPCNSigmaElNegVsPt"), isSignalLeg, track.pt(), track.tpcNSigmaEl()); + } + if (track.tpcNSigmaEl() < photonDauMinTPCNSigmaEl || track.tpcNSigmaEl() > photonDauMaxTPCNSigmaEl) { + return false; + } + return track.pt() <= photonDauMaxPt; + } + + // collisions of a track with only TPC time (all inside time range [t - backward, t + forward] plus a margin) + template + TrackCand tpcOnlyTrackCand(const TTrack& track, const std::vector& collisionBcNS, const std::vector& collisionTimeNS) + { + o2::aod::track::extensions::TPCTimeErrEncoding timeEncoding{}; + timeEncoding.encoding.timeErr = track.trackTimeRes(); + double trackTimeNS = collisionBcNS[track.collisionId()] + track.trackTime(); // trackTime() is relative to the BC of its collision + double timeMin = trackTimeNS - timeEncoding.getDeltaTBwd() - photonPoolTimeMarginNS; + double timeMax = trackTimeNS + timeEncoding.getDeltaTFwd() + photonPoolTimeMarginNS; + int firstCollIdx = track.collisionId(); + int lastCollIdx = firstCollIdx; + for (int collIdx = 0; collIdx < static_cast(collisionTimeNS.size()); ++collIdx) { + if (collisionTimeNS[collIdx] >= timeMin && collisionTimeNS[collIdx] <= timeMax) { + firstCollIdx = std::min(firstCollIdx, collIdx); + lastCollIdx = std::max(lastCollIdx, collIdx); + } + } + histos.fill(HIST("Photon/Inclusive/hTpcTimeRangeNColl"), lastCollIdx - firstCollIdx + 1); + return TrackCand{.Idxtr = static_cast(track.globalIndex()), .collBracket = {firstCollIdx, lastCollIdx}}; + } + + // electron-positron pairs with |theta+ - theta-| <= photonMaxDeltaTheta that share at least one collision (theta should be similar) + std::vector findPairsByPolarAngle(const std::vector& thetaOfTrack, int nCollisions) + { + // electrons of each collision sorted by polar angle, so a positron only looks at those inside its theta window + std::vector>> electronsByCollision(nCollisions); + for (int iElectron = 0; iElectron < static_cast(mElectronPool.size()); ++iElectron) { + const auto& electron = mElectronPool[iElectron]; + for (int collIdx = electron.collBracket.getMin(); collIdx <= electron.collBracket.getMax(); ++collIdx) { + electronsByCollision[collIdx].push_back({thetaOfTrack[electron.Idxtr], iElectron}); + } + } + for (size_t collIdx = 0; collIdx < electronsByCollision.size(); ++collIdx) { + std::sort(electronsByCollision[collIdx].begin(), electronsByCollision[collIdx].end()); + } + + std::vector pairs; + std::vector lastPositronPaired(mElectronPool.size(), -1); // an electron in several collisions is paired once per positron + for (int iPositron = 0; iPositron < static_cast(mPositronPool.size()); ++iPositron) { + const auto& positron = mPositronPool[iPositron]; + float thetaMin = thetaOfTrack[positron.Idxtr] - photonMaxDeltaTheta; + float thetaMax = thetaOfTrack[positron.Idxtr] + photonMaxDeltaTheta; + for (int collIdx = positron.collBracket.getMin(); collIdx <= positron.collBracket.getMax(); ++collIdx) { + const auto& electronsOfCollision = electronsByCollision[collIdx]; + auto electronIt = std::lower_bound(electronsOfCollision.begin(), electronsOfCollision.end(), std::make_pair(thetaMin, -1)); + for (; electronIt != electronsOfCollision.end() && electronIt->first <= thetaMax; ++electronIt) { + int iElectron = electronIt->second; + if (lastPositronPaired[iElectron] == iPositron) { + continue; + } + lastPositronPaired[iElectron] = iPositron; + const auto& electron = mElectronPool[iElectron]; + pairs.push_back(SVCand{.tr0Idx = electron.Idxtr, .tr1Idx = positron.Idxtr, .collBracket = electron.collBracket.getOverlap(positron.collBracket)}); + } + } + } + return pairs; + } + + // cuts before the fit (similar polar angles, and the circles of the two tracks in xy touching at the conversion point) + template + bool passesPreFitCuts(const TTrack& posTrack, const TTrack& negTrack, float deltaTheta, bool isSignal) + { + fillPhotonHist(HIST("hDeltaTheta"), isSignal, deltaTheta); + if (std::abs(deltaTheta) > photonMaxDeltaTheta) { + return false; + } + float sna = 0.f; + float csa = 0.f; + o2::math_utils::CircleXYf_t posCircle; + o2::math_utils::CircleXYf_t negCircle; + getTrackParCov(posTrack).getCircleParams(mBz, posCircle, sna, csa); + getTrackParCov(negTrack).getCircleParams(mBz, negCircle, sna, csa); + float circleTouchDist = std::abs(std::hypot(posCircle.xC - negCircle.xC, posCircle.yC - negCircle.yC) - (posCircle.rC + negCircle.rC)); + return circleTouchDist <= photonMaxCircleTouchDist; + } + + static constexpr int LastDaughterCutStep = 4; + + // photon vertex fitted in one collision + struct PhotonFit { + float dca = 0.f; + std::array secVtx{}; + std::array pPos{}; + std::array pNeg{}; + float posTrackZ = 0.f; // z of the positive daughter at the fit, after moving it to the collision's time + }; + + // moves a track with only a TPC time to the time of the collision + template + bool moveToCollision(const TCollision& collision, const TTrack& track, o2::track::TrackParCov& trackParCov) + { + if (!(track.flags() & o2::aod::track::TrackTimeAsym)) { + return true; + } + return mVDriftMgr.moveTPCTrack(collision, track, trackParCov); + } + + // vertex fit of the photon daughters, in the collinear mode if a daughter has no ITS + bool fitPhoton(o2::track::TrackParCov posTrackParCov, o2::track::TrackParCov negTrackParCov, bool collinear, PhotonFit& fit) + { + posTrackParCov.setPID(o2::track::PID::Electron); + negTrackParCov.setPID(o2::track::PID::Electron); + + // photon fit settings, reset afterwards since the fitter is shared with the Sigma+ vertex fit + fitter.setMatCorrType(o2::base::Propagator::MatCorrType::USEMatCorrLUT); + fitter.setMaxDXYIni(photonMaxDXYIni); + fitter.setCollinear(collinear); + int nCand = 0; + try { + nCand = fitter.process(posTrackParCov, negTrackParCov); + } catch (...) { + nCand = 0; + } + fitter.setMatCorrType(o2::base::Propagator::MatCorrType::USEMatCorrNONE); + fitter.setMaxDXYIni(1e9); + fitter.setCollinear(false); + if (nCand == 0 || !fitter.propagateTracksToVertex()) { + return false; + } + + fit.dca = std::sqrt(fitter.getChi2AtPCACandidate()); + fit.secVtx = fitter.getPCACandidatePos(); + fitter.getTrack(0).getPxPyPzGlo(fit.pPos); + fitter.getTrack(1).getPxPyPzGlo(fit.pNeg); + fit.posTrackZ = posTrackParCov.getZ(); + return true; + } + + // a photon is kept in every compatible collision where it passes the cuts, and the proton of that collision then decides + template + std::vector>> findPhotonsSelfPaired(const TCollisions& collisions, const TTracks& tracks, aod::AmbiguousTracks const& ambiguousTracks, aod::BCsWithTimestamps const& bcs) + { + std::vector>> photonsByCollision(collisions.size()); + + svPhotonPoolCreator.clearPools(); + svPhotonPoolCreator.fillBC2Coll(collisions, bcs); + mElectronPool.clear(); + mPositronPool.clear(); + + // BC start and time of every collision, in ns relative to the BC of the first collision + std::vector collisionBcNS(collisions.size()); + std::vector collisionTimeNS(collisions.size()); + uint64_t firstGlobalBC = collisions.begin().template bc_as().globalBC(); + for (const auto& collision : collisions) { + int64_t bcDiff = static_cast(collision.template bc_as().globalBC()) - static_cast(firstGlobalBC); + collisionBcNS[collision.globalIndex()] = bcDiff * o2::constants::lhc::LHCBunchSpacingNS; + collisionTimeNS[collision.globalIndex()] = collisionBcNS[collision.globalIndex()] + collision.collisionTime(); + } + + // select the daughter tracks and find the collisions each is compatible with + std::vector thetaOfTrack(tracks.size(), 0.f); + std::vector mcPhotonOfTrack(tracks.size(), -1); // MC: Sigma+ photon the track comes from + std::vector mcSigmaOfTrack(tracks.size(), -1); // MC: Sigma+ the track comes from + for (const auto& track : tracks) { + int mcPhoton = -1; + int mcSigma = -1; + if constexpr (IsMC) { + if (track.has_mcParticle()) { + mcSigma = findSigmaPlusAncestorOfPhotonDaughter(track.template mcParticle_as(), mcPhoton); + } + } + if (!selectPhotonDaughter(track, mcSigma >= 0)) { + continue; + } + thetaOfTrack[track.globalIndex()] = 2.f * std::atan(std::exp(-track.eta())); + mcPhotonOfTrack[track.globalIndex()] = mcPhoton; + mcSigmaOfTrack[track.globalIndex()] = mcSigma; + + if (track.flags() & o2::aod::track::TrackTimeAsym) { + // track with only a TPC time + if (!track.has_collision()) { + continue; + } + if (track.sign() < 0) { + mElectronPool.push_back(tpcOnlyTrackCand(track, collisionBcNS, collisionTimeNS)); + } else { + mPositronPool.push_back(tpcOnlyTrackCand(track, collisionBcNS, collisionTimeNS)); + } + } else { + // track with an ITS time (matched to collisions by svPoolCreator) + svPhotonPoolCreator.appendTrackCand(track, collisions, track.sign() < 0 ? PDG_t::kElectron : PDG_t::kPositron, ambiguousTracks, bcs); + } + } + + // svPoolCreator pools are ordered dau0 pos, dau0 neg, dau1 pos, dau1 neg: electrons are in 1, positrons are in 2 + auto timeMatchedPools = svPhotonPoolCreator.getTrackCandPool(); + mElectronPool.insert(mElectronPool.end(), timeMatchedPools[1].begin(), timeMatchedPools[1].end()); + mPositronPool.insert(mPositronPool.end(), timeMatchedPools[2].begin(), timeMatchedPools[2].end()); + + auto pairs = findPairsByPolarAngle(thetaOfTrack, collisions.size()); + + // fit and select each pair in every good collision it is compatible with + uint64_t nTruePairs = 0; + for (const auto& pair : pairs) { + int electronIdx = pair.tr0Idx; + int positronIdx = pair.tr1Idx; + int sigmaId = -1; + if constexpr (IsMC) { + if (mcPhotonOfTrack[electronIdx] >= 0 && mcPhotonOfTrack[electronIdx] == mcPhotonOfTrack[positronIdx]) { + sigmaId = mcSigmaOfTrack[electronIdx]; + nTruePairs++; + } + } + bool isSignal = sigmaId >= 0; + + auto negTrack = tracks.rawIteratorAt(electronIdx); + auto posTrack = tracks.rawIteratorAt(positronIdx); + if (!passesPreFitCuts(posTrack, negTrack, thetaOfTrack[positronIdx] - thetaOfTrack[electronIdx], isSignal)) { + continue; + } + + // both daughters with only a TPC time on the same TPC side: moving to another collision shifts both by the same z, + // so the pair is fitted once and the vertex is only shifted in z for the other collisions + bool posTpcTime = posTrack.flags() & o2::aod::track::TrackTimeAsym; + bool negTpcTime = negTrack.flags() & o2::aod::track::TrackTimeAsym; + bool sameShiftForAllCollisions = posTpcTime && negTpcTime && ((posTrack.tgl() > 0.f) == (negTrack.tgl() > 0.f)); + bool collinearFit = !posTrack.hasITS() || !negTrack.hasITS(); + std::optional referenceFit; + + int maxStepReached = LastDaughterCutStep; // furthest selection step reached in any collision (for selection counter histograms) + auto passStep = [&](int step) { maxStepReached = std::max(maxStepReached, step); }; + bool anyFitOk = false; + float smallestDca = 1e9f; + bool openingAngleFilled = false; + int nCopies = 0; + + for (int collIdx = pair.collBracket.getMin(); collIdx <= pair.collBracket.getMax(); ++collIdx) { + if (!mGoodCollision[collIdx]) { + continue; + } + auto collision = collisions.rawIteratorAt(collIdx); + + PhotonFit fit; + if (sameShiftForAllCollisions && referenceFit) { + auto posTrackParCov = getTrackParCov(posTrack); + if (!moveToCollision(collision, posTrack, posTrackParCov)) { + continue; + } + fit = *referenceFit; + fit.secVtx[2] += posTrackParCov.getZ() - referenceFit->posTrackZ; + } else { + auto posTrackParCov = getTrackParCov(posTrack); + auto negTrackParCov = getTrackParCov(negTrack); + if (!moveToCollision(collision, posTrack, posTrackParCov) || !moveToCollision(collision, negTrack, negTrackParCov)) { + continue; + } + bool fitOk = fitPhoton(posTrackParCov, negTrackParCov, collinearFit, fit); + if (fitOk) { + anyFitOk = true; + smallestDca = std::min(smallestDca, fit.dca); + } + if (!fitOk || fit.dca > photonMaxDCAV0Dau) { + if (sameShiftForAllCollisions) { + break; // same result in every collision + } + continue; + } + if (sameShiftForAllCollisions) { + referenceFit = fit; + } + } + passStep(5); + + float radius = std::hypot(fit.secVtx[0], fit.secVtx[1]); + if (radius < photonMinRadius || radius > photonMaxRadius) { + continue; + } + passStep(6); + + float photonOpeningAngle = o2::aod::pwgem::dilepton::utils::pairutil::getOpeningAngle(fit.pPos[0], fit.pPos[1], fit.pPos[2], fit.pNeg[0], fit.pNeg[1], fit.pNeg[2]); + if (!openingAngleFilled) { // once per pair + openingAngleFilled = true; + fillPhotonHist(HIST("hOpeningAngle"), isSignal, photonOpeningAngle); + } + if (photonOpeningAngle > photonMaxOpeningAngle) { + continue; + } + passStep(8); // opening angle, and delta theta that is cut before the fit + + std::array pGamma{fit.pNeg[0] + fit.pPos[0], fit.pNeg[1] + fit.pPos[1], fit.pNeg[2] + fit.pPos[2]}; + float cosPA = RecoDecay::cpa(std::array{collision.posX(), collision.posY(), collision.posZ()}, fit.secVtx, pGamma); + if (cosPA < photonMinV0cospa) { + continue; + } + passStep(9); + + float photonY = RecoDecay::y(pGamma, o2::constants::physics::MassGamma); + if (photonY < photonMinRapidity || photonY > photonMaxRapidity) { + continue; + } + passStep(10); + + float qtarm = v0Qt(fit.pPos[0], fit.pPos[1], fit.pPos[2], fit.pNeg[0], fit.pNeg[1], fit.pNeg[2]); + float alpha = v0Alpha(fit.pPos[0], fit.pPos[1], fit.pPos[2], fit.pNeg[0], fit.pNeg[1], fit.pNeg[2]); + if (qtarm > photonMaxQt) { + continue; + } + passStep(11); + + if (std::abs(alpha) > photonMaxAlpha) { + continue; + } + passStep(12); + + float mGamma = RecoDecay::m(std::array{fit.pPos, fit.pNeg}, std::array{o2::constants::physics::MassElectron, o2::constants::physics::MassElectron}); + if (mGamma > photonMaxMass) { + continue; + } + passStep(13); + + fillPhotonHist(HIST("hMass"), isSignal, mGamma); + fillPhotonHist(HIST("hPt"), isSignal, std::hypot(pGamma[0], pGamma[1])); + fillPhotonHist(HIST("hRadius"), isSignal, radius); + fillPhotonHist(HIST("h2ArmenterosPodolanski"), isSignal, alpha, qtarm); + fillPhotonHist(HIST("h2ConvPointXY"), isSignal, fit.secVtx[0], fit.secVtx[1]); + + nCopies++; + photonsByCollision[collIdx].push_back({fit.secVtx[0], fit.secVtx[1], fit.secVtx[2], pGamma[0], pGamma[1], pGamma[2], mGamma, alpha, qtarm, radius, fit.dca, cosPA, negTrack, posTrack}); + } + + if (anyFitOk) { + fillPhotonHist(HIST("hDcaV0Daughters"), isSignal, smallestDca); + } + if (nCopies > 0) { + fillPhotonHist(HIST("hCollisionsPerPhoton"), isSignal, nCopies); + fillSigmaCounter(sigmaId, 4); + } + for (int step = 5; step <= maxStepReached; ++step) { + fillPhotonHist(HIST("hSelectionCounter"), isSignal, step); + } + } + + // steps 0-4 are the daughter-track cuts applied before the pairing, so they count all pairs + for (int step = 0; step <= LastDaughterCutStep; ++step) { + histos.fill(HIST("Photon/Inclusive/hSelectionCounter"), step, static_cast(pairs.size())); + if constexpr (IsMC) { + histos.fill(HIST("Photon/True/hSelectionCounter"), step, static_cast(nTruePairs)); + } + } + + return photonsByCollision; + } + // proton candidate selection template bool selectProton(const TTrack& track) @@ -614,17 +993,11 @@ struct Sigmaplusbuilder { bool isSignal = false; if constexpr (IsMC) { if (track.has_mcParticle()) { - auto mcProton = track.template mcParticle_as(); - isSignal = findSigmaPlusMotherOfProton(mcProton) >= 0; + isSignal = findSigmaPlusMotherOfProton(track.template mcParticle_as()) >= 0; } } auto fillProtonStep = [&](int step) { - histos.fill(HIST("Proton/hSelectionCounter"), step); - if constexpr (IsMC) { - if (isSignal) { - histos.fill(HIST("Proton/True/hSelectionCounter"), step); - } - } + fillProtonHist(HIST("hSelectionCounter"), isSignal, step); }; fillProtonStep(0); @@ -638,12 +1011,7 @@ struct Sigmaplusbuilder { } 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()); - } - } + fillProtonHist(HIST("h2TPCNClsVsPt"), isSignal, track.pt(), track.tpcNClsFound()); if (track.tpcNClsFound() < protonMinTpcNCls) { return false; } @@ -654,8 +1022,6 @@ struct Sigmaplusbuilder { } 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()) { if (std::abs(track.tofNSigmaPr()) > protonMaxTOFNSigma) { @@ -667,37 +1033,21 @@ struct Sigmaplusbuilder { } 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())); - } - } + fillProtonHist(HIST("h2DcaToPVVsPt"), isSignal, 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()); + fillProtonHist(HIST("hPt"), isSignal, track.pt()); + fillProtonHist(HIST("h2TPCNSigmaVsPt"), isSignal, 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()); - } - } + fillProtonHist(HIST("h2TOFNSigmaVsPt"), isSignal, track.pt(), track.tofNSigmaPr()); } return true; } - // small vector helpers used in the Sigma+ reconstruction static float dot3(const std::array& a, const std::array& b) { return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; @@ -726,31 +1076,38 @@ struct Sigmaplusbuilder { return mothers.front().globalIndex(); } + // MC: for an e+ or e- from Sigma+ -> pi0 -> photon -> e+e-, returns the Sigma+ index and sets photonIndex (-1 otherwise) template - int findSigmaPlusMotherOfPhoton(const TMcPart& mcPos, const TMcPart& mcNeg) + int findSigmaPlusAncestorOfPhotonDaughter(const TMcPart& mcDaughter, int& photonIndex) { - auto const& posMothers = mcPos.template mothers_as(); - auto const& negMothers = mcNeg.template mothers_as(); - if (posMothers.empty() || negMothers.empty()) { - return -1; - } - auto mcGamma = posMothers.front(); - if (mcGamma.globalIndex() != negMothers.front().globalIndex() || mcGamma.pdgCode() != PDG_t::kGamma) { + auto const& mothers = mcDaughter.template mothers_as(); + if (mothers.empty() || mothers.front().pdgCode() != PDG_t::kGamma) { return -1; } - - auto const& pi0Mothers = mcGamma.template mothers_as(); + auto const& pi0Mothers = mothers.front().template mothers_as(); if (pi0Mothers.empty() || std::abs(pi0Mothers.front().pdgCode()) != PDG_t::kPi0) { return -1; } - auto const& sigmaMothers = pi0Mothers.front().template mothers_as(); if (sigmaMothers.empty() || std::abs(sigmaMothers.front().pdgCode()) != PDG_t::kSigmaPlus) { return -1; } + photonIndex = mothers.front().globalIndex(); return sigmaMothers.front().globalIndex(); } + template + int findSigmaPlusMotherOfPhoton(const TMcPart& mcPos, const TMcPart& mcNeg) + { + int posPhoton = -1; + int negPhoton = -1; + int sigmaId = findSigmaPlusAncestorOfPhotonDaughter(mcPos, posPhoton); + if (sigmaId < 0 || findSigmaPlusAncestorOfPhotonDaughter(mcNeg, negPhoton) < 0 || posPhoton != negPhoton) { + return -1; + } + return sigmaId; + } + template bool isSigmaPlusToProtonPi0(const TMcPart& mcPart) { @@ -823,21 +1180,23 @@ struct Sigmaplusbuilder { } // Build a Sigma+ -> p pi0 candidate from a proton track and a PCM photon - // 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) + void buildSigmaPlusCandidate(const TTrack& protonTrack, const PhotonCand& photon, const std::array& pv) { auto posTrack = photon.posTrack; auto negTrack = photon.negTrack; + int photonLegsWithoutIts = (posTrack.hasITS() ? 0 : 1) + (negTrack.hasITS() ? 0 : 1); bool isSignal = false; - bool collisionIdCheck = false; // MC only: true if the proton's true MC collision matches the reconstructed collision - std::array mcTrueVtx{}; // Sigma+ decay vertex - std::array mcTrueMomProton{}; // true MC proton momentum - std::array mcTrueMomGamma{}; // true MC momentum of the measured photon - std::array mcTrueMomSigmaPlus{}; // true MC momentum of the Sigma+ mother - float decayRadiusMC = -999.f; // MC-truth decay radius - float massMC = -999.f; // MC-truth Sigma+ mass + bool protonIsSignal = false; + bool photonIsSignal = false; + bool collisionIdCheck = false; // MC only: true if the proton's MC collision matches the reconstructed collision + std::array mcTrueVtx{}; + std::array mcTrueMomProton{}; + std::array mcTrueMomGamma{}; + std::array mcTrueMomSigmaPlus{}; + float decayRadiusMC = -999.f; + float massMC = -999.f; int matchedSigmaId = -1; if constexpr (IsMC) { if (protonTrack.has_mcParticle()) { @@ -847,9 +1206,9 @@ struct Sigmaplusbuilder { if (protonCollision.has_mcCollision()) { collisionIdCheck = protonCollision.mcCollision().globalIndex() == mcProton.mcCollisionId(); } - auto const& protonMothers = mcProton.template mothers_as(); int protonSigmaIdx = findSigmaPlusMotherOfProton(mcProton); + protonIsSignal = protonSigmaIdx >= 0; if (posTrack.has_mcParticle() && negTrack.has_mcParticle()) { auto mcPos = posTrack.template mcParticle_as(); @@ -862,13 +1221,14 @@ struct Sigmaplusbuilder { } int photonSigmaIdx = findSigmaPlusMotherOfPhoton(mcPos, mcNeg); - isSignal = (protonSigmaIdx >= 0 && photonSigmaIdx >= 0 && photonSigmaIdx == protonSigmaIdx); + photonIsSignal = photonSigmaIdx >= 0; + isSignal = protonIsSignal && photonIsSignal && photonSigmaIdx == protonSigmaIdx; } if (isSignal) { + auto mcSigmaPlusMother = mcProton.template mothers_as().front(); 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()); @@ -878,18 +1238,14 @@ struct Sigmaplusbuilder { } auto fillCandStep = [&](int step) { - histos.fill(HIST("Candidate/hSelectionCounter"), step); - if constexpr (IsMC) { - if (isSignal) { - histos.fill(HIST("Candidate/True/hSelectionCounter"), step); - } - } + fillCandHist(HIST("hSelectionCounter"), isSignal, step); + fillCandHist(HIST("h2SelectionCounterVsLegType"), isSignal, step, photonLegsWithoutIts); }; fillCandStep(0); // all pairs - // Reject the pair if the proton track is one of the photon's own e+/e- daughter tracks. + // reject the pair if the proton track is one of the photon's own daughter tracks if (protonTrack.globalIndex() == posTrack.globalIndex() || protonTrack.globalIndex() == negTrack.globalIndex()) { - return -1; + return; } fillCandStep(1); // autocorrelation @@ -906,54 +1262,43 @@ struct Sigmaplusbuilder { try { nCand = fitter.process(protonTrackParCov, photonTrackParCov); } catch (...) { - return -1; + return; } if (nCand == 0 || !fitter.propagateTracksToVertex()) { - return -1; + return; } - fillCandStep(2); // Vertex fit + fillCandStep(2); // vertex fit + std::array protonAtPca{}; + std::array photonAtPca{}; + fitter.getTrack(0).getXYZGlo(protonAtPca); + fitter.getTrack(1).getXYZGlo(photonAtPca); - float fitChi2 = fitter.getChi2AtPCACandidate(); - float dcaProtonGamma = std::sqrt(fitChi2); - histos.fill(HIST("Candidate/hDcaProtonGamma"), dcaProtonGamma); - if constexpr (IsMC) { - if (isSignal) { - histos.fill(HIST("Candidate/True/hDcaProtonGamma"), dcaProtonGamma); - } - } + float dcaProtonGamma = std::sqrt(fitter.getChi2AtPCACandidate()); + fillCandHist(HIST("hDcaProtonGamma"), isSignal, dcaProtonGamma); if (dcaProtonGamma > candMaxDcaProtonGamma) { - return -1; + return; } fillCandStep(3); // DCA(p,gamma) - std::array secVtx = fitter.getPCACandidatePos(); - float radius = std::hypot(secVtx[0], secVtx[1]); - histos.fill(HIST("Candidate/hRadius"), radius); - if constexpr (IsMC) { - if (isSignal) { - histos.fill(HIST("Candidate/True/hRadius"), radius); - float vtxResFromMc = std::hypot(secVtx[0] - mcTrueVtx[0], secVtx[1] - mcTrueVtx[1], secVtx[2] - mcTrueVtx[2]); - histos.fill(HIST("Candidate/True/hVertexResFromMcTruth"), vtxResFromMc); - } + // decay vertex between the proton's and the photon's point of closest approach (relative weight given by candVertexProtonWeight) + std::array secVtx{}; + for (size_t i = 0; i < secVtx.size(); ++i) { + secVtx[i] = candVertexProtonWeight * protonAtPca[i] + (1.f - candVertexProtonWeight) * photonAtPca[i]; } + float radius = std::hypot(secVtx[0], secVtx[1]); + fillCandHist(HIST("hRadius"), isSignal, radius); if (radius < candMinRadius || radius > candMaxRadius) { - return -1; + return; } 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); + std::array nHatOrig = normalize3(flightVec); float flightDistance = std::sqrt(dot3(flightVec, flightVec)); - histos.fill(HIST("Candidate/hFlightDistance"), flightDistance); - if constexpr (IsMC) { - if (isSignal) { - histos.fill(HIST("Candidate/True/hFlightDistance"), flightDistance); - } - } + fillCandHist(HIST("hFlightDistance"), isSignal, flightDistance); if (flightDistance < candMinFlightDistance || flightDistance > candMaxFlightDistance) { - return -1; + return; } fillCandStep(5); // flight distance @@ -962,110 +1307,93 @@ struct Sigmaplusbuilder { fitter.getTrack(0).getPxPyPzGlo(pProton); fitter.getTrack(1).getPxPyPzGlo(pGamma1); - if constexpr (IsMC) { - if (isSignal) { - float trueProtonP = std::sqrt(dot3(mcTrueMomProton, mcTrueMomProton)); - float fitProtonP = std::sqrt(dot3(pProton, pProton)); - histos.fill(HIST("Candidate/True/hProtonMomResFromMcTruth"), (fitProtonP - trueProtonP) / trueProtonP); - - float trueGammaP = std::sqrt(dot3(mcTrueMomGamma, mcTrueMomGamma)); - float fitGammaP = std::sqrt(dot3(pGamma1, pGamma1)); - histos.fill(HIST("Candidate/True/hPhotonMomResFromMcTruth"), (fitGammaP - trueGammaP) / trueGammaP); - } - } - - float e1 = std::sqrt(dot3(pGamma1, pGamma1)); // |p_gamma1| + float e1 = std::sqrt(dot3(pGamma1, pGamma1)); float massPi0 = o2::constants::physics::MassPionNeutral; + constexpr float CoefADegenerateThreshold = 1e-6f; + std::array nHat{}; std::array eOut{}; std::array eIn{}; - float coefA = 0.f, coefB = 0.f, coefK = 0.f, coefL = 0.f; + float coefA = 0.f, coefB = 0.f; float pGamma2In = 0.f, pGamma2Out = 0.f, tPerp2 = 0.f; - float discriminant = -1.f; - std::array nHatOrig = nHat; // each retry perturbs around this - constexpr float CoefADegenerateThreshold = 1e-6f; // guards against dividing by a near-zero quadratic coefficient - - if constexpr (IsMC) { - if (isSignal) { - float cosProtonFlight = dot3(pProton, nHatOrig) / std::sqrt(dot3(pProton, pProton)); - float protonFlightAngle = std::acos(std::clamp(cosProtonFlight, -1.f, 1.f)); - histos.fill(HIST("Candidate/True/hProtonFlightAngle"), protonFlightAngle); - } - } - - int discrIter = 0; - for (; discrIter <= discrRetryMaxIter; ++discrIter) { - if (discrIter > 0) { - TVector3 nHatVec(nHatOrig[0], nHatOrig[1], nHatOrig[2]); - float dTheta = mThetaResoFunc.GetRandom(); - float dPhi = mPhiResoFunc.GetRandom(); - nHatVec.SetMagThetaPhi(1., nHatVec.Theta() + dTheta, nHatVec.Phi() + dPhi); - nHat = {static_cast(nHatVec.X()), static_cast(nHatVec.Y()), static_cast(nHatVec.Z())}; - } + auto solveForFlightDir = [&](const std::array& nHatUse) { + nHat = nHatUse; eOut = normalize3(cross3(pProton, nHat)); // normal to the decay plane - eIn = cross3(eOut, nHat); // in-plane, transverse to n + eIn = cross3(eOut, nHat); // in the decay plane, transverse to n - float pProtonIn = dot3(pProton, eIn); // proton has no out-of-plane component, by construction + float pProtonIn = dot3(pProton, eIn); float pGamma1Long = dot3(pGamma1, nHat); float pGamma1In = dot3(pGamma1, eIn); float pGamma1Out = dot3(pGamma1, eOut); - // in-plane transverse momentum cancels between p, gamma1, gamma2 pGamma2In = -(pProtonIn + pGamma1In); - // out-of-plane momentum cancels between gamma1 and gamma2 alone pGamma2Out = -pGamma1Out; tPerp2 = pGamma2In * pGamma2In + pGamma2Out * pGamma2Out; - // m_pi0^2 = 2*(|p_gamma1||p_gamma2| - p_gamma1.p_gamma2), solved for - // the remaining unknown (p_gamma2 along n) -> quadratic coefA*x^2+coefB*x+coefC=0 - float dot1 = pGamma1In * pGamma2In + pGamma1Out * pGamma2Out; - coefK = massPi0 * massPi0 + 2.f * dot1; - coefL = 2.f * pGamma1Long; - + // m_pi0^2 = 2*(|p_gamma1||p_gamma2| - p_gamma1.p_gamma2) as coefA*x^2 + coefB*x + coefC = 0 + float coefK = massPi0 * massPi0 + 2.f * (pGamma1In * pGamma2In + pGamma1Out * pGamma2Out); + float coefL = 2.f * pGamma1Long; coefA = 4.f * e1 * e1 - coefL * coefL; if (std::abs(coefA) < CoefADegenerateThreshold) { - continue; + return -1.f; } coefB = -2.f * coefK * coefL; float coefC = 4.f * e1 * e1 * tPerp2 - coefK * coefK; + return coefB * coefB - 4.f * coefA * coefC; + }; - discriminant = coefB * coefB - 4.f * coefA * coefC; - if (discriminant >= 0.f) { - break; + // if the measured flight direction gives no real solution, flight directions on rings of increasing + // tilt around it are tried (the first ring with a real solution is taken, and on it the direction with the largest discriminant) + float tiltAngle = 0.f; + float discriminant = solveForFlightDir(nHatOrig); + if (discriminant < 0.f) { + constexpr float MaxAbsNzForZHelperAxis = 0.9f; + std::array helperAxis = std::abs(nHatOrig[2]) < MaxAbsNzForZHelperAxis ? std::array{0.f, 0.f, 1.f} : std::array{1.f, 0.f, 0.f}; + std::array tiltAxisU = normalize3(cross3(nHatOrig, helperAxis)); + std::array tiltAxisV = cross3(nHatOrig, tiltAxisU); + int nRings = static_cast(std::lround(candMaxTilt / candTiltStep)); + for (int ring = 1; ring <= nRings && discriminant < 0.f; ++ring) { + float tilt = ring * candTiltStep; + float bestRingDisc = -1.f; + std::array bestRingDir{}; + for (int iAzimuth = 0; iAzimuth < candTiltNAzimuth; ++iAzimuth) { + float azimuth = o2::constants::math::TwoPI * iAzimuth / candTiltNAzimuth; + std::array tiltDir{}; + for (size_t i = 0; i < tiltDir.size(); ++i) { + tiltDir[i] = std::cos(tilt) * nHatOrig[i] + std::sin(tilt) * (std::cos(azimuth) * tiltAxisU[i] + std::sin(azimuth) * tiltAxisV[i]); + } + float disc = solveForFlightDir(tiltDir); + if (disc >= 0.f && disc > bestRingDisc) { + bestRingDisc = disc; + bestRingDir = tiltDir; + } + } + if (bestRingDisc >= 0.f) { + discriminant = solveForFlightDir(bestRingDir); + tiltAngle = tilt; + } } } - histos.fill(HIST("Candidate/hDiscriminantRetryIter"), discrIter); - histos.fill(HIST("Candidate/hDiscriminant"), discriminant); - if constexpr (IsMC) { - if (isSignal) { - histos.fill(HIST("Candidate/True/hDiscriminantRetryIter"), discrIter); - histos.fill(HIST("Candidate/True/hDiscriminant"), discriminant); - } - } + float tiltForHist = discriminant >= 0.f ? tiltAngle : 0.4999f; // no real solution + fillCandHist(HIST("hTiltAngle"), isSignal, tiltForHist); if (discriminant < 0.f) { - return -1; + return; } fillCandStep(6); // real root - // two roots from the quadratic, among both we keep the mass closest to the nominal Sigma+ mass + // of the two roots, keep the one giving 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); - } - } + fillCandHist(HIST("hRootCenter"), isSignal, rootCenter); if (candRejectNegRootCenter && rootCenter < 0.f) { - return -1; + return; } if (rootCenter > candMaxRootCenter) { - return -1; + return; } bool haveCandidate = false; @@ -1092,31 +1420,26 @@ struct Sigmaplusbuilder { } } if (!haveCandidate) { - return -1; + return; } 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); - } - } + fillCandHist(HIST("hDcaToPV"), isSignal, candDcaToPV); if (candDcaToPV > candMaxDcaToPV) { - return -1; + return; } fillCandStep(8); // DCA to PV - // AntiSigmaPointingAngle - Fake pointing angle + // angle between the proton momentum and the flight direction, both rotated back to undo the bending + // of the proton between its reference point and the decay vertex 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]); @@ -1138,196 +1461,295 @@ struct Sigmaplusbuilder { 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]}; + flightVec[0] * std::cos(alphaVtxRot) - flightVec[1] * std::sin(alphaVtxRot), + flightVec[0] * std::sin(alphaVtxRot) + flightVec[1] * std::cos(alphaVtxRot), + flightVec[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 + fillCandHist(HIST("hAntiSigmaPointingAngle"), isSignal, antiSigmaPointingAngle); 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); - } + return; } + fillCandHist(HIST("hMassSigmaPlus"), isSignal, bestMass); + fillCandHist(HIST("h2MassVsPt"), isSignal, ptSigma, bestMass); + fillCandHist(HIST("h2MassVsTilt"), isSignal, tiltAngle, bestMass); if (bestMass > candMaxSigmaMass) { - return -1; + return; } 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); - } - } + fillCandHist(HIST("hRapidity"), isSignal, candRapidity); if (std::abs(candRapidity) > candMaxRapidity) { - return -1; + return; } 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); - } - } + // photon opening angle from the daughters' track momenta + float photonOpeningAngle = o2::aod::pwgem::dilepton::utils::pairutil::getOpeningAngle(posTrack.px(), posTrack.py(), posTrack.pz(), negTrack.px(), negTrack.py(), negTrack.pz()); + fillCandHist(HIST("hPhotonOpeningAngle"), isSignal, photonOpeningAngle); if (photonOpeningAngle > candMaxPhotonOpeningAngle) { - return -1; + return; } fillCandStep(11); // photon opening angle - // photon pointing angle + // angle between the photon momentum and the line from the decay vertex to the conversion point 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); - } - } + fillCandHist(HIST("hPhotonPointingAngle"), isSignal, photonPointingAngle); if (photonPointingAngle > candMaxPhotonPointingAngle) { - return -1; + return; } fillCandStep(12); // photon pointing angle - // 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 pvToConv{pv[0] - photon.x, pv[1] - photon.y, pv[2] - photon.z}; 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); - } - } + fillCandHist(HIST("hPhotonDcaToPV"), isSignal, photonDcaToPV); if (photonDcaToPV > candMaxPhotonDcaToPV) { - return -1; + return; } fillCandStep(13); // photon DCA to PV - fillCandStep(14); // filled + float photonPt2 = photon.px * photon.px + photon.py * photon.py; + float tAtDcaXY = photonPt2 > 0.f ? ((pv[0] - photon.x) * photon.px + (pv[1] - photon.y) * photon.py) / photonPt2 : 0.f; + float photonDcaXYToPV = std::hypot(photon.x + tAtDcaXY * photon.px - pv[0], photon.y + tAtDcaXY * photon.py - pv[1]); + float photonDcaZToPV = photon.z + tAtDcaXY * photon.pz - pv[2]; + float psiPair = photonPsiPair(posTrack, negTrack, photon.radius); + + SigmaPlusCandidate cand; + cand.photonId = photonId(posTrack.globalIndex(), negTrack.globalIndex()); + cand.isSignal = isSignal; + cand.photonLegsWithoutIts = photonLegsWithoutIts; + cand.matchedSigmaId = matchedSigmaId; + cand.decVtx = secVtx; + cand.radius = radius; + cand.flightDistance = flightDistance; + cand.dcaProtonGamma = dcaProtonGamma; + cand.dcaToPV = candDcaToPV; + cand.tiltAngle = tiltAngle; + cand.rootCenter = rootCenter; + cand.antiSigmaPointingAngle = antiSigmaPointingAngle; + cand.pProton = pProton; + cand.pGamma1 = pGamma1; + cand.pGamma2 = bestMomGamma2; + cand.protonTpcNSigma = protonTrack.tpcNSigmaPr(); + cand.protonTofNSigma = protonTrack.tofNSigmaPr(); + cand.protonSign = protonTrack.sign(); + cand.protonItsNCls = protonTrack.itsNCls(); + cand.protonTpcNCls = protonTrack.tpcNClsFound(); + cand.protonDcaXY = protonTrack.dcaXY(); + cand.protonDcaZ = protonTrack.dcaZ(); + cand.photonMass = photon.mGamma; + cand.photonAlpha = photon.alpha; + cand.photonQt = photon.qtarm; + cand.photonRadius = photon.radius; + cand.photonOpeningAngle = photonOpeningAngle; + cand.photonPointingAngle = photonPointingAngle; + cand.photonDcaDau = photon.dcaDau; + cand.photonCosPAToPV = photon.cosPAToPV; + cand.photonDcaXYToPV = photonDcaXYToPV; + cand.photonDcaZToPV = photonDcaZToPV; + cand.photonPsiPair = psiPair; + cand.posTpcNSigmaEl = posTrack.tpcNSigmaEl(); + cand.negTpcNSigmaEl = negTrack.tpcNSigmaEl(); + cand.posItsNCls = posTrack.itsNCls(); + cand.negItsNCls = negTrack.itsNCls(); + cand.posTpcNCls = posTrack.tpcNClsFound(); + cand.negTpcNCls = negTrack.tpcNClsFound(); + cand.posTpcNClsFindable = posTrack.tpcNClsFindable(); + cand.negTpcNClsFindable = negTrack.tpcNClsFindable(); + cand.posTpcChi2NCl = posTrack.tpcChi2NCl(); + cand.negTpcChi2NCl = negTrack.tpcChi2NCl(); + cand.posDcaXY = posTrack.dcaXY(); + cand.posDcaZ = posTrack.dcaZ(); + cand.negDcaXY = negTrack.dcaXY(); + cand.negDcaZ = negTrack.dcaZ(); + cand.collisionIdCheck = collisionIdCheck; + cand.protonIsSignal = protonIsSignal; + cand.photonIsSignal = photonIsSignal; + cand.decVtxMC = mcTrueVtx; + cand.pProtonMC = mcTrueMomProton; + cand.pGammaMC = mcTrueMomGamma; + cand.pSigmaPlusMC = mcTrueMomSigmaPlus; + cand.decayRadiusMC = decayRadiusMC; + cand.massMC = massMC; + mCandidatesOfTimeframe.push_back(cand); + } + + template + void fillCandidateTables(const SigmaPlusCandidate& cand) + { 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; + slimSigmaPlusCandsMC(cand.radius, + cand.dcaToPV, cand.tiltAngle, cand.dcaProtonGamma, + cand.protonSign, + cand.protonDcaXY, cand.protonDcaZ, + cand.pProton[0], cand.pProton[1], cand.pProton[2], + cand.pGamma1[0], cand.pGamma1[1], cand.pGamma1[2], + cand.pGamma2[0], cand.pGamma2[1], cand.pGamma2[2], + cand.protonTpcNSigma, cand.protonTofNSigma, + cand.posTpcNSigmaEl, cand.negTpcNSigmaEl, + cand.photonMass, + cand.photonDcaDau, cand.photonCosPAToPV, cand.photonDcaXYToPV, cand.photonDcaZToPV, cand.photonPsiPair, + cand.posDcaXY, cand.posDcaZ, cand.negDcaXY, cand.negDcaZ, + cand.posTpcNClsFindable, cand.negTpcNClsFindable, cand.posTpcChi2NCl, cand.negTpcChi2NCl, + cand.collisionIdCheck, + cand.isSignal, cand.protonIsSignal, cand.photonIsSignal, + cand.decayRadiusMC, cand.massMC, + cand.pSigmaPlusMC[0], cand.pSigmaPlusMC[1], cand.pSigmaPlusMC[2]); + return; } - sigmaPlusCandsMC(secVtx[0], secVtx[1], secVtx[2], - 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.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, - mcTrueVtx[0], mcTrueVtx[1], mcTrueVtx[2], - mcTrueMomProton[0], mcTrueMomProton[1], mcTrueMomProton[2], - mcTrueMomGamma[0], mcTrueMomGamma[1], mcTrueMomGamma[2], - mcTrueMomSigmaPlus[0], mcTrueMomSigmaPlus[1], mcTrueMomSigmaPlus[2], - decayRadiusMC, massMC); + sigmaPlusCandsMC(cand.decVtx[0], cand.decVtx[1], cand.decVtx[2], + cand.flightDistance, cand.dcaProtonGamma, + cand.pProton[0], cand.pProton[1], cand.pProton[2], + cand.pGamma1[0], cand.pGamma1[1], cand.pGamma1[2], + cand.pGamma2[0], cand.pGamma2[1], cand.pGamma2[2], + cand.protonTpcNSigma, cand.protonTofNSigma, + cand.posTpcNSigmaEl, cand.negTpcNSigmaEl, + cand.photonMass, cand.photonAlpha, cand.photonQt, cand.photonRadius, + cand.photonOpeningAngle, cand.photonPointingAngle, + cand.rootCenter, cand.antiSigmaPointingAngle, cand.dcaToPV, cand.tiltAngle, + cand.protonSign, + cand.protonItsNCls, cand.protonTpcNCls, cand.protonDcaXY, cand.protonDcaZ, + cand.posItsNCls, cand.posTpcNCls, cand.negItsNCls, cand.negTpcNCls, + cand.photonDcaDau, cand.photonCosPAToPV, cand.photonDcaXYToPV, cand.photonDcaZToPV, cand.photonPsiPair, + cand.posDcaXY, cand.posDcaZ, cand.negDcaXY, cand.negDcaZ, + cand.posTpcNClsFindable, cand.negTpcNClsFindable, cand.posTpcChi2NCl, cand.negTpcChi2NCl, + cand.collisionIdCheck, + cand.isSignal, cand.protonIsSignal, cand.photonIsSignal, + cand.decVtxMC[0], cand.decVtxMC[1], cand.decVtxMC[2], + cand.pProtonMC[0], cand.pProtonMC[1], cand.pProtonMC[2], + cand.pGammaMC[0], cand.pGammaMC[1], cand.pGammaMC[2], + cand.pSigmaPlusMC[0], cand.pSigmaPlusMC[1], cand.pSigmaPlusMC[2], + cand.decayRadiusMC, cand.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; + slimSigmaPlusCands(cand.radius, + cand.dcaToPV, cand.tiltAngle, cand.dcaProtonGamma, + cand.protonSign, + cand.protonDcaXY, cand.protonDcaZ, + cand.pProton[0], cand.pProton[1], cand.pProton[2], + cand.pGamma1[0], cand.pGamma1[1], cand.pGamma1[2], + cand.pGamma2[0], cand.pGamma2[1], cand.pGamma2[2], + cand.protonTpcNSigma, cand.protonTofNSigma, + cand.posTpcNSigmaEl, cand.negTpcNSigmaEl, + cand.photonMass, + cand.photonDcaDau, cand.photonCosPAToPV, cand.photonDcaXYToPV, cand.photonDcaZToPV, cand.photonPsiPair, + cand.posDcaXY, cand.posDcaZ, cand.negDcaXY, cand.negDcaZ, + cand.posTpcNClsFindable, cand.negTpcNClsFindable, cand.posTpcChi2NCl, cand.negTpcChi2NCl); + return; + } + sigmaPlusCands(cand.decVtx[0], cand.decVtx[1], cand.decVtx[2], + cand.flightDistance, cand.dcaProtonGamma, + cand.pProton[0], cand.pProton[1], cand.pProton[2], + cand.pGamma1[0], cand.pGamma1[1], cand.pGamma1[2], + cand.pGamma2[0], cand.pGamma2[1], cand.pGamma2[2], + cand.protonTpcNSigma, cand.protonTofNSigma, + cand.posTpcNSigmaEl, cand.negTpcNSigmaEl, + cand.photonMass, cand.photonAlpha, cand.photonQt, cand.photonRadius, + cand.photonOpeningAngle, cand.photonPointingAngle, + cand.rootCenter, cand.antiSigmaPointingAngle, cand.dcaToPV, cand.tiltAngle, + cand.protonSign, + cand.protonItsNCls, cand.protonTpcNCls, cand.protonDcaXY, cand.protonDcaZ, + cand.posItsNCls, cand.posTpcNCls, cand.negItsNCls, cand.negTpcNCls, + cand.photonDcaDau, cand.photonCosPAToPV, cand.photonDcaXYToPV, cand.photonDcaZToPV, cand.photonPsiPair, + cand.posDcaXY, cand.posDcaZ, cand.negDcaXY, cand.negDcaZ, + cand.posTpcNClsFindable, cand.negTpcNClsFindable, cand.posTpcChi2NCl, cand.negTpcChi2NCl); + } + } + + // The same photon can be in candidates of several collisions (with candDeduplicatePhotons only the one with the smallest DCA to PV is written) + template + void deduplicateAndWriteCandidates() + { + std::unordered_map> candsByPhoton; + for (int iCand = 0; iCand < static_cast(mCandidatesOfTimeframe.size()); ++iCand) { + candsByPhoton[mCandidatesOfTimeframe[iCand].photonId].push_back(iCand); + } + std::vector keep(mCandidatesOfTimeframe.size(), !candDeduplicatePhotons); + for (const auto& [photon, cands] : candsByPhoton) { + histos.fill(HIST("Candidate/Inclusive/hCandidatesPerPhoton"), cands.size()); + if (candDeduplicatePhotons) { + int best = cands[0]; + for (const int& iCand : cands) { + if (mCandidatesOfTimeframe[iCand].dcaToPV < mCandidatesOfTimeframe[best].dcaToPV) { + best = iCand; + } + } + keep[best] = true; } - sigmaPlusCands(secVtx[0], secVtx[1], secVtx[2], - 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.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; + } + + for (int iCand = 0; iCand < static_cast(mCandidatesOfTimeframe.size()); ++iCand) { + if (!keep[iCand]) { + continue; + } + const auto& cand = mCandidatesOfTimeframe[iCand]; + fillCandHist(HIST("hSelectionCounter"), cand.isSignal, 14); + fillCandHist(HIST("h2SelectionCounterVsLegType"), cand.isSignal, 14, cand.photonLegsWithoutIts); + if constexpr (IsMC) { + if (cand.isSignal) { + fillSigmaCounter(cand.matchedSigmaId, 5); + auto genIt = mGenSigmaWritten.find(cand.matchedSigmaId); + if (genIt != mGenSigmaWritten.end()) { + genIt->second = true; + } + } + } + fillCandidateTables(cand); + } + mCandidatesOfTimeframe.clear(); } - void initCCDB(aod::BCs::iterator const& bc) + void initCCDB(aod::BCsWithTimestamps::iterator const& bc) { if (mRunNumber == bc.runNumber()) { return; } mRunNumber = bc.runNumber(); - o2::parameters::GRPMagField* grpmag = ccdb->getForRun(grpmagPath, mRunNumber); + auto* grpmag = ccdb->getForRun(grpmagPath, mRunNumber); o2::base::Propagator::initFieldFromGRP(grpmag); mBz = grpmag->getNominalL3Field(); fitter.setBz(mBz); + if (!mLut) { + mLut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->get(lutPath)); + } + o2::base::Propagator::Instance()->setMatLUT(mLut); LOG(info) << "Task initialized for run " << mRunNumber << " with magnetic field " << mBz << " kZG"; } - void processData(CollisionsFull const& collisions, aod::V0Datas const& v0s, TracksFull const& tracks, aod::BCs const&) + void processData(CollisionsFull const& collisions, aod::V0Datas const& v0s, TracksFull const& tracks, aod::AmbiguousTracks const& ambiguousTracks, aod::BCsWithTimestamps const& bcs) { + markGoodCollisions(collisions); + std::vector>> photonsByCollision; + if (useCustomVertexer && collisions.size() > 0) { + auto firstBC = collisions.begin().bc_as(); + initCCDB(firstBC); + mVDriftMgr.update(firstBC.timestamp()); + photonsByCollision = findPhotonsSelfPaired(collisions, tracks, ambiguousTracks, bcs); + } + for (const auto& collision : collisions) { if (std::abs(collision.posZ()) > cutZVertex || !collision.sel8()) { continue; } - initCCDB(collision.bc_as()); + 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()); - auto acceptedPhotons = findPhotonsFromV0s(v0sThisCollision, tracksThisCollision, pv); + std::vector> acceptedPhotons; + if (useCustomVertexer) { + acceptedPhotons = photonsByCollision[collision.globalIndex()]; + } else { + auto v0sThisCollision = v0s.sliceBy(v0PerCollision, collision.globalIndex()); + acceptedPhotons = findPhotonsFromV0s(v0sThisCollision, tracksThisCollision, pv); + } std::vector acceptedProtons; for (const auto& track : tracksThisCollision) { @@ -1342,105 +1764,102 @@ struct Sigmaplusbuilder { } } } + deduplicateAndWriteCandidates(); } 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, aod::McCollisions const&) + void processMc(CollisionsFullMC const& collisions, aod::V0Datas const& v0s, TracksFullMC const& tracks, aod::AmbiguousTracks const& ambiguousTracks, aod::BCsWithTimestamps const& bcs, aod::McParticles const& mcParticles, aod::McCollisions const&) { - std::vector matchedSigmaPlusMcIds; // Sigma+ MC indices that got at least one signal candidate + // generated Sigma+ -> p pi0 in acceptance + mGenSigmaWritten.clear(); + for (const auto& mcPart : mcParticles) { + if (!isSigmaPlusToProtonPi0(mcPart) || std::abs(mcPart.y()) > cutRapMotherMC || mcPart.pt() < cutPtGenMC) { + continue; + } + mGenSigmaWritten[mcPart.globalIndex()] = false; + histos.fill(HIST("MC/hGenSigmaPlusPt"), mcPart.pt()); + histos.fill(HIST("MC/hSigmaPlusCounter"), 0); + } + + constexpr int ElectronLeg = 1; + constexpr int PositronLeg = 2; + std::unordered_map> legsOfPhoton; + for (const auto& track : tracks) { + if (!track.has_mcParticle()) { + continue; + } + auto mcParticle = track.template mcParticle_as(); + if (std::abs(mcParticle.pdgCode()) == PDG_t::kProton) { + fillSigmaCounter(findSigmaPlusMotherOfProton(mcParticle), 1); + } else if (std::abs(mcParticle.pdgCode()) == PDG_t::kElectron) { + int photonIndex = -1; + std::array conversionVertex{}; + int sigmaId = findSigmaPlusMotherOfConversionElectron(mcParticle, photonIndex, conversionVertex); + if (sigmaId >= 0) { + auto& legs = legsOfPhoton[photonIndex]; + legs.first = sigmaId; + legs.second |= (mcParticle.pdgCode() == PDG_t::kElectron ? ElectronLeg : PositronLeg); + } + } + } + for (const auto& [photonIndex, legs] : legsOfPhoton) { + if (legs.second == (ElectronLeg | PositronLeg)) { + fillSigmaCounter(legs.first, 3); + } + } + + markGoodCollisions(collisions); + std::vector>> photonsByCollision; + if (useCustomVertexer && collisions.size() > 0) { + auto firstBC = collisions.begin().bc_as(); + initCCDB(firstBC); + mVDriftMgr.update(firstBC.timestamp()); + photonsByCollision = findPhotonsSelfPaired(collisions, tracks, ambiguousTracks, bcs); + } for (const auto& collision : collisions) { if (std::abs(collision.posZ()) > cutZVertex || !collision.sel8()) { continue; } - initCCDB(collision.bc_as()); + 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()); - 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); - } - } - } - } - } + std::vector> acceptedPhotons; + if (useCustomVertexer) { + acceptedPhotons = photonsByCollision[collision.globalIndex()]; + } else { + auto v0sThisCollision = v0s.sliceBy(v0PerCollision, collision.globalIndex()); + acceptedPhotons = findPhotonsFromV0s(v0sThisCollision, tracksThisCollision, pv); } std::vector acceptedProtons; for (const auto& track : tracksThisCollision) { if (selectProton(track)) { acceptedProtons.push_back(track); - - histos.fill(HIST("MC/hProtonTruthQA"), 0); if (track.has_mcParticle()) { - histos.fill(HIST("MC/hProtonTruthQA"), 1); - auto mcProton = track.template mcParticle_as(); - if (findSigmaPlusMotherOfProton(mcProton) >= 0) { - histos.fill(HIST("MC/hProtonTruthQA"), 2); - } + fillSigmaCounter(findSigmaPlusMotherOfProton(track.template mcParticle_as()), 2); } } } for (const auto& photon : acceptedPhotons) { for (const auto& proton : acceptedProtons) { - int matchedId = buildSigmaPlusCandidate(proton, photon, pv); - if (matchedId >= 0) { - matchedSigmaPlusMcIds.push_back(matchedId); - } + buildSigmaPlusCandidate(proton, photon, pv); } } } + deduplicateAndWriteCandidates(); - // all generated Sigma+ -> p pi0 decays, regardless of reconstruction + // a table row with the true values for each generated Sigma+ in acceptance without a written candidate for (const auto& mcPart : mcParticles) { - if (!isSigmaPlusToProtonPi0(mcPart)) { - continue; - } - if (std::abs(mcPart.y()) > cutRapMotherMC) { + auto genIt = mGenSigmaWritten.find(mcPart.globalIndex()); + if (genIt == mGenSigmaWritten.end() || genIt->second) { 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}; @@ -1456,7 +1875,7 @@ struct Sigmaplusbuilder { if (fillSlimTables) { slimSigmaPlusCandsMC(-999.f, - -999.f, -999.f, + -999.f, -999.f, -999.f, 0, -999.f, -999.f, -999.f, -999.f, -999.f, @@ -1465,8 +1884,11 @@ struct Sigmaplusbuilder { -999.f, -999.f, -999.f, -999.f, -999.f, + -999.f, -999.f, -999.f, -999.f, -999.f, + -999.f, -999.f, -999.f, -999.f, + 0, 0, -999.f, -999.f, false, - false, + false, false, false, genDecayRadiusMC, genMassMC, mcPart.px(), mcPart.py(), mcPart.pz()); continue; @@ -1480,13 +1902,16 @@ struct Sigmaplusbuilder { -999.f, -999.f, -999.f, -999.f, -999.f, -999.f, -999.f, -999.f, - -999.f, -999.f, -999.f, - -999.f, -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, + -999.f, -999.f, -999.f, -999.f, -999.f, + -999.f, -999.f, -999.f, -999.f, + 0, 0, -999.f, -999.f, false, - false, + false, false, false, genDecVtx[0], genDecVtx[1], genDecVtx[2], genMomProton[0], genMomProton[1], genMomProton[2], -999.f, -999.f, -999.f, @@ -1496,14 +1921,15 @@ struct Sigmaplusbuilder { } PROCESS_SWITCH(Sigmaplusbuilder, processMc, "Process MC", false); - void processFindable(CollisionsFullMC const& collisions, aod::V0Datas const& v0s, TracksFullMC const& tracks, aod::McParticles const&, aod::BCs const&) + void processFindable(CollisionsFullMC const& collisions, aod::V0Datas const& v0s, TracksFullMC const& tracks, aod::McParticles const&, aod::BCsWithTimestamps const&) { + mGenSigmaWritten.clear(); constexpr int MinDauTpcCls = 90; for (const auto& collision : collisions) { if (std::abs(collision.posZ()) > cutZVertex || !collision.sel8()) { continue; } - initCCDB(collision.bc_as()); + 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()};