diff --git a/PWGCF/Femto/FemtoNuclei/TableProducer/HadNucleiFemto.cxx b/PWGCF/Femto/FemtoNuclei/TableProducer/HadNucleiFemto.cxx index 53a26284222..063ba4dd841 100644 --- a/PWGCF/Femto/FemtoNuclei/TableProducer/HadNucleiFemto.cxx +++ b/PWGCF/Femto/FemtoNuclei/TableProducer/HadNucleiFemto.cxx @@ -79,7 +79,7 @@ using std::array; using CollBracket = o2::math_utils::Bracket; using HyperCandidates = aod::DataHypCandsWColl; -using HyperCandidatesMC = aod::MCHypCands; +using HyperCandidatesMC = aod::MCHypCandsWColl; using CollisionsFull = soa::Join; using CollisionsFullMC = soa::Join; using HadHyperCollisionsFull = soa::Join; @@ -448,14 +448,14 @@ struct HadNucleiFemto { float phiPi{0.f}; bool isMatter{false}; - float pt() const { return std::hypot(momentum[0], momentum[1]); } - float eta() const + [[nodiscard]] float pt() const { return std::hypot(momentum[0], momentum[1]); } + [[nodiscard]] float eta() const { const float transverseMomentum = pt(); return transverseMomentum > 0.f ? std::asinh(momentum[2] / transverseMomentum) : 999.f; } - float phi() const { return std::atan2(momentum[1], momentum[0]); } - int8_t sign() const { return isMatter ? 1 : -1; } + [[nodiscard]] float phi() const { return std::atan2(momentum[1], momentum[0]); } + [[nodiscard]] int8_t sign() const { return isMatter ? 1 : -1; } }; struct HadHyperHadron { @@ -470,10 +470,10 @@ struct HadNucleiFemto { float phiValue{0.f}; int8_t signValue{0}; - float pt() const { return std::hypot(momentum[0], momentum[1]); } - float eta() const { return etaValue; } - float phi() const { return phiValue; } - int8_t sign() const { return signValue; } + [[nodiscard]] float pt() const { return std::hypot(momentum[0], momentum[1]); } + [[nodiscard]] float eta() const { return etaValue; } + [[nodiscard]] float phi() const { return phiValue; } + [[nodiscard]] int8_t sign() const { return signValue; } }; struct HadHyperEvent { diff --git a/PWGLF/DataModel/LFHypernucleiTables.h b/PWGLF/DataModel/LFHypernucleiTables.h index a4162c55f75..c737e867f82 100644 --- a/PWGLF/DataModel/LFHypernucleiTables.h +++ b/PWGLF/DataModel/LFHypernucleiTables.h @@ -134,6 +134,36 @@ DECLARE_SOA_TABLE(DataHypCandsFlow, "AOD", "HYPCANDSFLOW", hyperrec::Flags, hyperrec::TrackedClSize); DECLARE_SOA_TABLE(MCHypCands, "AOD", "MCHYPCANDS", + o2::soa::Index<>, + hyperrec::CentralityFT0A, hyperrec::CentralityFT0C, hyperrec::CentralityFT0M, + hyperrec::TrackOccupancyInTimeRange, hyperrec::Ft0cOccupancyInTimeRange, + hyperrec::XPrimVtx, hyperrec::YPrimVtx, hyperrec::ZPrimVtx, + + hyperrec::RunNumber, hyperrec::IsMatter, + hyperrec::PtHe3, hyperrec::PhiHe3, hyperrec::EtaHe3, + hyperrec::PtPi, hyperrec::PhiPi, hyperrec::EtaPi, + hyperrec::XDecVtx, hyperrec::YDecVtx, hyperrec::ZDecVtx, + hyperrec::DcaV0Daug, hyperrec::DcaHe, hyperrec::DcaPi, + hyperrec::NSigmaHe, hyperrec::NTPCclusHe, hyperrec::NTPCclusPi, hyperrec::NTPCpidClusHe, hyperrec::NTPCpidClusPi, hyperrec::NTPCCrossedRowsHe, hyperrec::NTPCCrossedRowsPi, + hyperrec::TpcMomHe, hyperrec::TpcMomPi, hyperrec::TpcSignalHe, hyperrec::TpcSignalPi, hyperrec::TpcChi2He, hyperrec::ItsChi2He, hyperrec::ItsChi2Pi, + hyperrec::TofMass, + hyperrec::ItsClusterSizesHe, hyperrec::ItsClusterSizesPi, + hyperrec::Flags, hyperrec::TrackedClSize, + hyperrec::GenPt, + hyperrec::GenPhi, + hyperrec::GenEta, + hyperrec::GenPtHe3, + hyperrec::GenXDecVtx, + hyperrec::GenYDecVtx, + hyperrec::GenZDecVtx, + hyperrec::IsReco, + hyperrec::IsFakeHeOnITSLayer, + hyperrec::IsSignal, + hyperrec::IsRecoMCCollision, + hyperrec::IsSurvEvSel, + hyperrec::IsTwoBodyDecay, aod::mcparticle::StatusCode); + +DECLARE_SOA_TABLE(MCHypCandsWColl, "AOD", "MCHYPCANDSWCOLL", o2::soa::Index<>, hyperrec::CollisionId, hyperrec::HeTrackId, hyperrec::PiTrackId, hyperrec::CentralityFT0A, hyperrec::CentralityFT0C, hyperrec::CentralityFT0M, @@ -185,6 +215,7 @@ DECLARE_SOA_TABLE(DataHypCandsWColl, "AOD", "HYPCANDSWCOLL", using DataHypCand = DataHypCands::iterator; using DataHypCandFlow = DataHypCandsFlow::iterator; using MCHypCand = MCHypCands::iterator; +using MCHypCandWColl = MCHypCandsWColl::iterator; using DataHypCandWColl = DataHypCandsWColl::iterator; namespace hyperkink diff --git a/PWGLF/DataModel/LFPiHypertritonFemtoTables.h b/PWGLF/DataModel/LFPiHypertritonFemtoTables.h new file mode 100644 index 00000000000..ba6416885c0 --- /dev/null +++ b/PWGLF/DataModel/LFPiHypertritonFemtoTables.h @@ -0,0 +1,239 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// +/// \file LFPiHypertritonFemtoTables.h +/// \brief Slim tables for pi-hypertriton femto pairs +/// + +#ifndef PWGLF_DATAMODEL_LFPIHYPERTRITONFEMTOTABLES_H_ +#define PWGLF_DATAMODEL_LFPIHYPERTRITONFEMTOTABLES_H_ + +#include + +#include + +namespace o2::aod +{ +namespace pihypertritonfemto +{ +constexpr uint8_t ClosePairAngular = 1 << 0; +constexpr uint8_t ClosePairDistance = 1 << 1; +constexpr uint8_t ClosePairAngularUnavailable = 1 << 2; +constexpr uint8_t ClosePairDistanceUnavailable = 1 << 3; +DECLARE_SOA_COLUMN(IsMixed, isMixed, bool); +// Offline CPR bits, OR-ed across same-sign daughters: 1=angular rejection, 2=distance rejection, +// 4=angular unavailable, 8=distance unavailable. Zero means all applicable checks passed. +// Opposite-sign comparisons require no CPR. All flagged pairs remain in the table. +DECLARE_SOA_COLUMN(IsClosePairRejected, isClosePairRejected, uint8_t); +DECLARE_SOA_COLUMN(MixingDepth, mixingDepth, int); +DECLARE_SOA_COLUMN(PosZ, posZ, float); +DECLARE_SOA_COLUMN(CentFT0C, centFT0C, float); +DECLARE_SOA_COLUMN(TrackOccupancy, trackOccupancy, int); +DECLARE_SOA_COLUMN(Ft0cOccupancy, ft0cOccupancy, float); +DECLARE_SOA_COLUMN(MultFT0C, multFT0C, float); +DECLARE_SOA_COLUMN(XPrimVtx, xPrimVtx, float); +DECLARE_SOA_COLUMN(YPrimVtx, yPrimVtx, float); +DECLARE_SOA_COLUMN(ZPrimVtx, zPrimVtx, float); +DECLARE_SOA_COLUMN(HypIsMatter, hypIsMatter, bool); +DECLARE_SOA_COLUMN(HypPtHe3, hypPtHe3, float); +DECLARE_SOA_COLUMN(HypEtaHe3, hypEtaHe3, float); +DECLARE_SOA_COLUMN(HypPhiHe3, hypPhiHe3, float); +DECLARE_SOA_COLUMN(HypPtPi, hypPtPi, float); +DECLARE_SOA_COLUMN(HypEtaPi, hypEtaPi, float); +DECLARE_SOA_COLUMN(HypPhiPi, hypPhiPi, float); +DECLARE_SOA_COLUMN(HypDcaV0Daug, hypDcaV0Daug, float); +DECLARE_SOA_COLUMN(HypDcaHe, hypDcaHe, float); +DECLARE_SOA_COLUMN(HypDcaPi, hypDcaPi, float); +DECLARE_SOA_COLUMN(HypNSigmaHe, hypNSigmaHe, float); +DECLARE_SOA_COLUMN(HypNTPCCrossedRowsHe, hypNTPCCrossedRowsHe, uint8_t); +DECLARE_SOA_COLUMN(HypNTPCCrossedRowsPi, hypNTPCCrossedRowsPi, uint8_t); +DECLARE_SOA_COLUMN(HypTpcMomHe, hypTpcMomHe, float); +DECLARE_SOA_COLUMN(HypTpcMomPi, hypTpcMomPi, float); +DECLARE_SOA_COLUMN(HypTpcSignalHe, hypTpcSignalHe, uint16_t); +DECLARE_SOA_COLUMN(HypTpcSignalPi, hypTpcSignalPi, uint16_t); +DECLARE_SOA_COLUMN(HypItsClusterSizesHe, hypItsClusterSizesHe, uint32_t); +DECLARE_SOA_COLUMN(HypItsClusterSizesPi, hypItsClusterSizesPi, uint32_t); +DECLARE_SOA_COLUMN(HypXDecVtx, hypXDecVtx, float); +DECLARE_SOA_COLUMN(HypYDecVtx, hypYDecVtx, float); +DECLARE_SOA_COLUMN(HypZDecVtx, hypZDecVtx, float); +DECLARE_SOA_COLUMN(HadPt, hadPt, float); +DECLARE_SOA_COLUMN(HadEta, hadEta, float); +DECLARE_SOA_COLUMN(HadPhi, hadPhi, float); +DECLARE_SOA_COLUMN(HadSign, hadSign, int8_t); +DECLARE_SOA_COLUMN(HadDcaXY, hadDcaXY, float); +DECLARE_SOA_COLUMN(HadTpcNClsCrossedRows, hadTpcNClsCrossedRows, uint8_t); +DECLARE_SOA_COLUMN(HadTpcNClsPID, hadTpcNClsPID, uint8_t); +DECLARE_SOA_COLUMN(HadTpcChi2NCl, hadTpcChi2NCl, float); +DECLARE_SOA_COLUMN(HadItsClusterSizes, hadItsClusterSizes, uint32_t); +DECLARE_SOA_COLUMN(HadItsChi2NCl, hadItsChi2NCl, float); +DECLARE_SOA_COLUMN(HadHasTOF, hadHasTOF, bool); +DECLARE_SOA_COLUMN(HadTpcNSigmaPi, hadTpcNSigmaPi, float); +DECLARE_SOA_COLUMN(HadTofNSigmaPi, hadTofNSigmaPi, float); +DECLARE_SOA_COLUMN(HypGenPt, hypGenPt, float); +DECLARE_SOA_COLUMN(HypGenEta, hypGenEta, float); +DECLARE_SOA_COLUMN(HypGenPhi, hypGenPhi, float); +DECLARE_SOA_COLUMN(HypGenXDecVtx, hypGenXDecVtx, float); +DECLARE_SOA_COLUMN(HypGenYDecVtx, hypGenYDecVtx, float); +DECLARE_SOA_COLUMN(HypGenZDecVtx, hypGenZDecVtx, float); +DECLARE_SOA_COLUMN(HypIsReco, hypIsReco, bool); +DECLARE_SOA_COLUMN(HypIsSignal, hypIsSignal, bool); +DECLARE_SOA_COLUMN(HypIsRecoMCCollision, hypIsRecoMCCollision, bool); +DECLARE_SOA_COLUMN(HypIsSurvEvSel, hypIsSurvEvSel, bool); +DECLARE_SOA_COLUMN(HypIsTwoBodyDecay, hypIsTwoBodyDecay, bool); +DECLARE_SOA_COLUMN(HypStatusCode, hypStatusCode, int16_t); +DECLARE_SOA_COLUMN(HeGenPt, heGenPt, float); +DECLARE_SOA_COLUMN(HeIsPhysicalPrimary, heIsPhysicalPrimary, bool); +DECLARE_SOA_COLUMN(DecayPiIsPhysicalPrimary, decayPiIsPhysicalPrimary, bool); +DECLARE_SOA_COLUMN(HadRecoPt, hadRecoPt, float); +DECLARE_SOA_COLUMN(HadRecoEta, hadRecoEta, float); +DECLARE_SOA_COLUMN(HadRecoPhi, hadRecoPhi, float); +DECLARE_SOA_COLUMN(HadGenPt, hadGenPt, float); +DECLARE_SOA_COLUMN(HadGenEta, hadGenEta, float); +DECLARE_SOA_COLUMN(HadGenPhi, hadGenPhi, float); +DECLARE_SOA_COLUMN(HadIsPhysicalPrimary, hadIsPhysicalPrimary, bool); +DECLARE_SOA_COLUMN(HadProcess, hadProcess, int16_t); +DECLARE_SOA_COLUMN(SameMCCollision, sameMCCollision, bool); +DECLARE_SOA_COLUMN(MatchesHypRecoMCCollision, matchesHypRecoMCCollision, bool); +DECLARE_SOA_COLUMN(MatchesPairRecoMCCollision, matchesPairRecoMCCollision, bool); +DECLARE_SOA_COLUMN(IsTruthSelfCorrelation, isTruthSelfCorrelation, bool); +DECLARE_SOA_COLUMN(IsTruePrimaryHadHyperPair, isTruePrimaryHadHyperPair, bool); +DECLARE_SOA_COLUMN(HadIsTruePion, hadIsTruePion, bool); // Truth species, independent of SE/ME and primary status. +} // namespace pihypertritonfemto + +DECLARE_SOA_TABLE(PiHypertritonFemtoTable, "AOD", "PIHYPFEMTO", + pihypertritonfemto::IsMixed, + pihypertritonfemto::IsClosePairRejected, + pihypertritonfemto::MixingDepth, + pihypertritonfemto::PosZ, + pihypertritonfemto::CentFT0C, + pihypertritonfemto::TrackOccupancy, + pihypertritonfemto::Ft0cOccupancy, + pihypertritonfemto::MultFT0C, + pihypertritonfemto::XPrimVtx, + pihypertritonfemto::YPrimVtx, + pihypertritonfemto::ZPrimVtx, + pihypertritonfemto::HypIsMatter, + pihypertritonfemto::HypPtHe3, + pihypertritonfemto::HypEtaHe3, + pihypertritonfemto::HypPhiHe3, + pihypertritonfemto::HypPtPi, + pihypertritonfemto::HypEtaPi, + pihypertritonfemto::HypPhiPi, + pihypertritonfemto::HypDcaV0Daug, + pihypertritonfemto::HypDcaHe, + pihypertritonfemto::HypDcaPi, + pihypertritonfemto::HypNSigmaHe, + pihypertritonfemto::HypNTPCCrossedRowsHe, + pihypertritonfemto::HypNTPCCrossedRowsPi, + pihypertritonfemto::HypTpcMomHe, + pihypertritonfemto::HypTpcMomPi, + pihypertritonfemto::HypTpcSignalHe, + pihypertritonfemto::HypTpcSignalPi, + pihypertritonfemto::HypItsClusterSizesHe, + pihypertritonfemto::HypItsClusterSizesPi, + pihypertritonfemto::HypXDecVtx, + pihypertritonfemto::HypYDecVtx, + pihypertritonfemto::HypZDecVtx, + pihypertritonfemto::HadPt, + pihypertritonfemto::HadEta, + pihypertritonfemto::HadPhi, + pihypertritonfemto::HadSign, + pihypertritonfemto::HadDcaXY, + pihypertritonfemto::HadTpcNClsCrossedRows, + pihypertritonfemto::HadTpcNClsPID, + pihypertritonfemto::HadTpcChi2NCl, + pihypertritonfemto::HadItsClusterSizes, + pihypertritonfemto::HadItsChi2NCl, + pihypertritonfemto::HadHasTOF, + pihypertritonfemto::HadTpcNSigmaPi, + pihypertritonfemto::HadTofNSigmaPi); + +DECLARE_SOA_TABLE(PiHypertritonFemtoTableMC, "AOD", "PIHYPFEMTOMC", + pihypertritonfemto::IsMixed, + pihypertritonfemto::IsClosePairRejected, + pihypertritonfemto::MixingDepth, + pihypertritonfemto::PosZ, + pihypertritonfemto::CentFT0C, + pihypertritonfemto::TrackOccupancy, + pihypertritonfemto::Ft0cOccupancy, + pihypertritonfemto::MultFT0C, + pihypertritonfemto::XPrimVtx, + pihypertritonfemto::YPrimVtx, + pihypertritonfemto::ZPrimVtx, + pihypertritonfemto::HypIsMatter, + pihypertritonfemto::HypPtHe3, + pihypertritonfemto::HypEtaHe3, + pihypertritonfemto::HypPhiHe3, + pihypertritonfemto::HypPtPi, + pihypertritonfemto::HypEtaPi, + pihypertritonfemto::HypPhiPi, + pihypertritonfemto::HypDcaV0Daug, + pihypertritonfemto::HypDcaHe, + pihypertritonfemto::HypDcaPi, + pihypertritonfemto::HypNSigmaHe, + pihypertritonfemto::HypNTPCCrossedRowsHe, + pihypertritonfemto::HypNTPCCrossedRowsPi, + pihypertritonfemto::HypTpcMomHe, + pihypertritonfemto::HypTpcMomPi, + pihypertritonfemto::HypTpcSignalHe, + pihypertritonfemto::HypTpcSignalPi, + pihypertritonfemto::HypItsClusterSizesHe, + pihypertritonfemto::HypItsClusterSizesPi, + pihypertritonfemto::HypXDecVtx, + pihypertritonfemto::HypYDecVtx, + pihypertritonfemto::HypZDecVtx, + pihypertritonfemto::HadPt, + pihypertritonfemto::HadEta, + pihypertritonfemto::HadPhi, + pihypertritonfemto::HadSign, + pihypertritonfemto::HadDcaXY, + pihypertritonfemto::HadTpcNClsCrossedRows, + pihypertritonfemto::HadTpcNClsPID, + pihypertritonfemto::HadTpcChi2NCl, + pihypertritonfemto::HadItsClusterSizes, + pihypertritonfemto::HadItsChi2NCl, + pihypertritonfemto::HadHasTOF, + pihypertritonfemto::HadTpcNSigmaPi, + pihypertritonfemto::HadTofNSigmaPi, + pihypertritonfemto::HypGenPt, + pihypertritonfemto::HypGenEta, + pihypertritonfemto::HypGenPhi, + pihypertritonfemto::HypGenXDecVtx, + pihypertritonfemto::HypGenYDecVtx, + pihypertritonfemto::HypGenZDecVtx, + pihypertritonfemto::HypIsReco, + pihypertritonfemto::HypIsSignal, + pihypertritonfemto::HypIsRecoMCCollision, + pihypertritonfemto::HypIsSurvEvSel, + pihypertritonfemto::HypIsTwoBodyDecay, + pihypertritonfemto::HypStatusCode, + pihypertritonfemto::HeGenPt, + pihypertritonfemto::HeIsPhysicalPrimary, + pihypertritonfemto::DecayPiIsPhysicalPrimary, + pihypertritonfemto::HadRecoPt, + pihypertritonfemto::HadRecoEta, + pihypertritonfemto::HadRecoPhi, + pihypertritonfemto::HadGenPt, + pihypertritonfemto::HadGenEta, + pihypertritonfemto::HadGenPhi, + pihypertritonfemto::HadIsPhysicalPrimary, + pihypertritonfemto::HadProcess, + pihypertritonfemto::SameMCCollision, + pihypertritonfemto::MatchesHypRecoMCCollision, + pihypertritonfemto::MatchesPairRecoMCCollision, + pihypertritonfemto::IsTruthSelfCorrelation, + pihypertritonfemto::IsTruePrimaryHadHyperPair, + pihypertritonfemto::HadIsTruePion); + +} // namespace o2::aod + +#endif // PWGLF_DATAMODEL_LFPIHYPERTRITONFEMTOTABLES_H_ diff --git a/PWGLF/TableProducer/Nuspex/CMakeLists.txt b/PWGLF/TableProducer/Nuspex/CMakeLists.txt index eeadf523aac..ef89cf84878 100644 --- a/PWGLF/TableProducer/Nuspex/CMakeLists.txt +++ b/PWGLF/TableProducer/Nuspex/CMakeLists.txt @@ -49,6 +49,11 @@ o2physics_add_dpl_workflow(he3hadronfemto PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::EventFilteringUtils COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(pi-hypertriton-femto + SOURCES piHypertritonFemto.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter O2Physics::EventFilteringUtils + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(nuclei-spectra SOURCES nucleiSpectra.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DetectorsBase O2Physics::EventFilteringUtils diff --git a/PWGLF/TableProducer/Nuspex/hyperRecoTask.cxx b/PWGLF/TableProducer/Nuspex/hyperRecoTask.cxx index 8b9a782ea3c..1bf6fd11035 100644 --- a/PWGLF/TableProducer/Nuspex/hyperRecoTask.cxx +++ b/PWGLF/TableProducer/Nuspex/hyperRecoTask.cxx @@ -58,7 +58,6 @@ #include #include #include -#include #include #include @@ -89,16 +88,16 @@ static const std::vector particleName{"He3"}; } // namespace struct hyperCandidate { - float recoPtHe3() const { return std::hypot(momHe3[0], momHe3[1]); } - float recoPhiHe3() const { return std::atan2(momHe3[1], momHe3[0]); } - float recoEtaHe3() const { return std::asinh(momHe3[2] / recoPtHe3()); } - float recoPtPi() const { return std::hypot(momPi[0], momPi[1]); } - float recoPhiPi() const { return std::atan2(momPi[1], momPi[0]); } - float recoEtaPi() const { return std::asinh(momPi[2] / recoPtPi()); } - float genPt() const { return std::hypot(gMom[0], gMom[1]); } - float genPtHe3() const { return std::hypot(gMomHe3[0], gMomHe3[1]); } - float genPhi() const { return std::atan2(gMom[1], gMom[0]); } - float genEta() const { return std::asinh(gMom[2] / genPt()); } + [[nodiscard]] float recoPtHe3() const { return std::hypot(momHe3[0], momHe3[1]); } + [[nodiscard]] float recoPhiHe3() const { return std::atan2(momHe3[1], momHe3[0]); } + [[nodiscard]] float recoEtaHe3() const { return std::asinh(momHe3[2] / recoPtHe3()); } + [[nodiscard]] float recoPtPi() const { return std::hypot(momPi[0], momPi[1]); } + [[nodiscard]] float recoPhiPi() const { return std::atan2(momPi[1], momPi[0]); } + [[nodiscard]] float recoEtaPi() const { return std::asinh(momPi[2] / recoPtPi()); } + [[nodiscard]] float genPt() const { return std::hypot(gMom[0], gMom[1]); } + [[nodiscard]] float genPtHe3() const { return std::hypot(gMomHe3[0], gMomHe3[1]); } + [[nodiscard]] float genPhi() const { return std::atan2(gMom[1], gMom[0]); } + [[nodiscard]] float genEta() const { return std::asinh(gMom[2] / genPt()); } int v0ID = -1; int heTrackID = -1; @@ -150,6 +149,7 @@ struct HyperRecoTask { Produces outputDataTable; Produces outputDataTableWithFlow; Produces outputMCTable; + Produces outputMCTableWithCollID; Produces outputDataTableWithCollID; Service ccdb; Zorro zorro; @@ -297,7 +297,7 @@ struct HyperRecoTask { hEventsZorro->GetXaxis()->SetBinLabel(1, "Zorro before evsel"); hEventsZorro->GetXaxis()->SetBinLabel(2, "Zorro after evsel"); - if (doprocessMC || doprocessMCTracked) { + if (doprocessMC || doprocessMCWithCollID || doprocessMCTracked) { const auto hDecayChannel = qaRegistry.add("hDecayChannel", ";Decay channel; ", HistType::kTH1D, {{2, -0.5, 1.5}}); hDecayChannel->GetXaxis()->SetBinLabel(1, "2-body"); hDecayChannel->GetXaxis()->SetBinLabel(2, "3-body"); @@ -481,8 +481,9 @@ struct HyperRecoTask { if (collision.has_mcCollision()) { recoCollisionIds[collision.mcCollisionId()] = collision.globalIndex(); } - if (!collision.selection_bit(aod::evsel::kIsTriggerTVX) || !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) + if (!collision.selection_bit(aod::evsel::kIsTriggerTVX) || !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) { continue; + } qaRegistry.fill(HIST("hEvents"), 1.); @@ -535,8 +536,9 @@ struct HyperRecoTask { hypCand.clusterSizeITSPi = piTrack.itsClusterSizes(); bool heliumPID = heTrack.pidForTracking() == o2::track::PID::Helium3 || heTrack.pidForTracking() == o2::track::PID::Alpha; hypCand.momHe3TPC = (heliumPID && cfgCompensatePIDinTracking) ? heTrack.tpcInnerParam() / 2 : heTrack.tpcInnerParam(); - if (hypCand.momHe3TPC < tpcRigidityMinHe) + if (hypCand.momHe3TPC < tpcRigidityMinHe) { return; + } hypCand.momPiTPC = piTrack.tpcInnerParam(); qaRegistry.fill(HIST("hDeDxTot"), hypCand.momHe3TPC * heTrack.sign(), heTrack.tpcSignal()); qaRegistry.fill(HIST("hDeDxTot"), hypCand.momPiTPC * piTrack.sign(), piTrack.tpcSignal()); @@ -578,17 +580,21 @@ struct HyperRecoTask { hypMom[i] = hypCand.momHe3[i] + hypCand.momPi[i]; } float hypPt = std::hypot(hypMom[0], hypMom[1]); - if (hypPt < ptMin) + if (hypPt < ptMin) { return; + } float massH3L = std::sqrt(h3lE * h3lE - hypMom[0] * hypMom[0] - hypMom[1] * hypMom[1] - hypMom[2] * hypMom[2]); float massH4L = std::sqrt(h4lE * h4lE - hypMom[0] * hypMom[0] - hypMom[1] * hypMom[1] - hypMom[2] * hypMom[2]); bool isHypMass = false; - if (massH3L > o2::constants::physics::MassHyperTriton - massWidth && massH3L < o2::constants::physics::MassHyperTriton + massWidth) + if (massH3L > o2::constants::physics::MassHyperTriton - massWidth && massH3L < o2::constants::physics::MassHyperTriton + massWidth) { isHypMass = true; - if (massH4L > o2::constants::physics::MassHyperhydrog4 - massWidth && massH4L < o2::constants::physics::MassHyperhydrog4 + massWidth) + } + if (massH4L > o2::constants::physics::MassHyperhydrog4 - massWidth && massH4L < o2::constants::physics::MassHyperhydrog4 + massWidth) { isHypMass = true; - if (!isHypMass) + } + if (!isHypMass) { return; + } qaRegistry.fill(HIST("hH3LMassBefSel"), massH3L); qaRegistry.fill(HIST("hH4LMassBefSel"), massH4L); @@ -670,8 +676,9 @@ struct HyperRecoTask { auto posTrack = tracks.rawIteratorAt(v0.posTrackId()); auto negTrack = tracks.rawIteratorAt(v0.negTrackId()); - if (std::abs(posTrack.eta()) > etaMax || std::abs(negTrack.eta()) > etaMax) + if (std::abs(posTrack.eta()) > etaMax || std::abs(negTrack.eta()) > etaMax) { continue; + } // temporary fix: tpcInnhRecoCentralityColvsFT0CmultiplicityerParam() returns the momentum in all the software tags before: https://github.com/AliceO2Group/AliceO2/pull/12521 auto nSigmaTPCpos = computeNSigmaHe3(posTrack); @@ -710,24 +717,30 @@ struct HyperRecoTask { for (const auto& track : tracks) { - if (std::abs(track.eta()) > etaMax) + if (std::abs(track.eta()) > etaMax) { continue; + } - if (!track.hasITS()) + if (!track.hasITS()) { continue; + } auto nSigmaHe = computeNSigmaHe3(track); bool isHe = nSigmaHe > -1 * nSigmaMaxHe; int pdgHypo = isHe ? heDauPdg : PDG_t::kPiPlus; // LOG(info) << "ncls found: " << track.tpcNClsFound(); - if (isHe && track.tpcNClsFound() < nTPCClusMinHe) + if (isHe && track.tpcNClsFound() < nTPCClusMinHe) { continue; - if (!isHe && track.tpcNClsFound() < nTPCClusMinPi) + } + if (!isHe && track.tpcNClsFound() < nTPCClusMinPi) { continue; - if (isHe && track.tpcNClsCrossedRows() < nTPCCrossedRowsMinHe) + } + if (isHe && track.tpcNClsCrossedRows() < nTPCCrossedRowsMinHe) { continue; - if (!isHe && track.tpcNClsCrossedRows() < nTPCCrossedRowsMinPi) + } + if (!isHe && track.tpcNClsCrossedRows() < nTPCCrossedRowsMinPi) { continue; + } svCreator.appendTrackCand(track, collisions, pdgHypo, ambiguousTracks, bcs); } @@ -755,12 +768,15 @@ struct HyperRecoTask { if (mcTrackHe.has_mothers() && mcTrackPi.has_mothers()) { for (const auto& heMother : mcTrackHe.mothers_as()) { for (const auto& piMother : mcTrackPi.mothers_as()) { - if (heMother.globalIndex() != piMother.globalIndex()) + if (heMother.globalIndex() != piMother.globalIndex()) { continue; - if (std::abs(mcTrackHe.pdgCode()) != heDauPdg || std::abs(mcTrackPi.pdgCode()) != PDG_t::kPiPlus) + } + if (std::abs(mcTrackHe.pdgCode()) != heDauPdg || std::abs(mcTrackPi.pdgCode()) != PDG_t::kPiPlus) { continue; - if (std::abs(heMother.pdgCode()) != hyperPdg) + } + if (std::abs(heMother.pdgCode()) != hyperPdg) { continue; + } auto primVtx = std::array{heMother.vx(), heMother.vy(), heMother.vz()}; auto secVtx = std::array{mcTrackHe.vx(), mcTrackHe.vy(), mcTrackHe.vz()}; @@ -782,6 +798,143 @@ struct HyperRecoTask { } } + template + void fillMCCandidate(int64_t collisionId, int64_t heTrackId, int64_t piTrackId, Columns const&... columns) + { + if constexpr (withCollID) { + outputMCTableWithCollID(collisionId, heTrackId, piTrackId, columns...); + } else { + outputMCTable(columns...); + } + } + + template + void runMC(CollisionsFullMC const& collisions, aod::McCollisions const& mcCollisions, aod::V0s const& V0s, TracksFull const& tracks, aod::AmbiguousTracks const& ambiTracks, aod::BCsWithTimestamps const& bcs, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const& particlesMC) + { + filledMothers.clear(); + recoCollisionIds.clear(); + recoCollisionIds.resize(mcCollisions.size(), -1); + isSurvEvSelCollision.clear(); + isSurvEvSelCollision.resize(mcCollisions.size(), false); + goodCollision.clear(); + goodCollision.resize(collisions.size(), false); + hyperCandidates.clear(); + + selectGoodCollisionsMC(collisions); + useCustomVertexer ? fillCustomV0s(collisions, tracks, ambiTracks, bcs) : fillV0s(collisions, tracks, V0s); + fillMCinfo(trackLabelsMC, particlesMC); + for (const auto& hypCand : hyperCandidates) { + auto collision = collisions.rawIteratorAt(hypCand.collisionID); + if (!hypCand.isSignal && mcSignalOnly) { + continue; + } + int chargeFactor = -1 + 2 * (hypCand.pdgCode > 0); + float trackedHypClSize = !trackedClSize.empty() ? trackedClSize[hypCand.v0ID] : 0; + fillMCCandidate(hypCand.collisionID, hypCand.heTrackID, hypCand.piTrackID, + collision.centFT0A(), collision.centFT0C(), collision.centFT0M(), + collision.trackOccupancyInTimeRange(), collision.ft0cOccupancyInTimeRange(), + collision.posX(), collision.posY(), collision.posZ(), + mRunNumber, hypCand.isMatter, + hypCand.recoPtHe3(), hypCand.recoPhiHe3(), hypCand.recoEtaHe3(), + hypCand.recoPtPi(), hypCand.recoPhiPi(), hypCand.recoEtaPi(), + hypCand.decVtx[0], hypCand.decVtx[1], hypCand.decVtx[2], + hypCand.dcaV0dau, hypCand.he3DCAXY, hypCand.piDCAXY, + hypCand.nSigmaHe3, hypCand.nTPCClustersHe3, hypCand.nTPCClustersPi, hypCand.nTPCpidClusHe3, hypCand.nTPCpidClusPi, hypCand.nTPCCrossedRowsHe3, hypCand.nTPCCrossedRowsPi, + hypCand.momHe3TPC, hypCand.momPiTPC, hypCand.tpcSignalHe3, hypCand.tpcSignalPi, hypCand.tpcChi2He3, hypCand.itsChi2He3, hypCand.itsChi2Pi, + hypCand.massTOFHe3, + hypCand.clusterSizeITSHe3, hypCand.clusterSizeITSPi, hypCand.flags, trackedHypClSize, + chargeFactor * hypCand.genPt(), hypCand.genPhi(), hypCand.genEta(), hypCand.genPtHe3(), + hypCand.gDecVtx[0], hypCand.gDecVtx[1], hypCand.gDecVtx[2], + hypCand.isReco, hypCand.isFakeHeOnITSLayer, hypCand.isSignal, hypCand.isRecoMCCollision, hypCand.isSurvEvSelection, 1, 0); + } + + // now we fill only the signal candidates that were not reconstructed + for (const auto& mcPart : particlesMC) { + + if (std::abs(mcPart.pdgCode()) != hyperPdg) { + continue; + } + std::array secVtx{0.f, 0.f, 0.f}; + std::array lastDaugVtx{0.f, 0.f, 0.f}; + std::array primVtx = {mcPart.vx(), mcPart.vy(), mcPart.vz()}; + std::array momMother = {mcPart.px(), mcPart.py(), mcPart.pz()}; + std::array momHe3{0.f, 0.f, 0.f}; + bool isHeFound = false; + int mcProcess = {0}; + for (const auto& mcDaught : mcPart.daughters_as()) { + if (mcDaught.pdgCode() != PDG_t::kElectron) { // we do not care about delta electrons + lastDaugVtx = {mcDaught.vx(), mcDaught.vy(), mcDaught.vz()}; + mcProcess = mcDaught.getProcess(); + } + if (std::abs(mcDaught.pdgCode()) == heDauPdg) { + secVtx = lastDaugVtx; + momHe3 = {mcDaught.px(), mcDaught.py(), mcDaught.pz()}; + isHeFound = true; + } + } + if (mcPart.pdgCode() > 0) { + qaRegistry.fill(HIST("hIsMatterGen"), 0.); + } else { + qaRegistry.fill(HIST("hIsMatterGen"), 1.); + } + if (!isHeFound) { + qaRegistry.fill(HIST("hDecayChannel"), 1.); + } + qaRegistry.fill(HIST("hDecayChannel"), 0.); + if (mcPart.pdgCode() > 0) { + qaRegistry.fill(HIST("hIsMatterGenTwoBody"), 0.); + } else { + qaRegistry.fill(HIST("hIsMatterGenTwoBody"), 1.); + } + if (std::find(filledMothers.begin(), filledMothers.end(), mcPart.globalIndex()) != std::end(filledMothers)) { + continue; + } + hyperCandidate hypCand; + hypCand.pdgCode = mcPart.pdgCode(); + hypCand.isRecoMCCollision = recoCollisionIds[mcPart.mcCollisionId()] >= 0; + hypCand.isSurvEvSelection = isSurvEvSelCollision[mcPart.mcCollisionId()]; + int chargeFactor = -1 + 2 * (hypCand.pdgCode > 0); + for (std::size_t i = 0; i < hypCand.gDecVtx.size(); i++) { + hypCand.gDecVtx[i] = (isHeFound ? secVtx[i] : lastDaugVtx[i]) - primVtx[i]; + hypCand.gMom[i] = momMother[i]; + hypCand.gMomHe3[i] = momHe3[i]; + } + hypCand.heTrackID = -1; + hypCand.piTrackID = -1; + hypCand.isSignal = true; + + float centFT0A = -1, centFT0C = -1, centFT0M = -1; + int trackOccupancyInTimeRange = -1; + float ft0cOccupancyInTimeRange = -1.f; + int recoCollisionId = -1; + if (hypCand.isRecoMCCollision) { + recoCollisionId = recoCollisionIds[mcPart.mcCollisionId()]; + auto recoCollision = collisions.rawIteratorAt(recoCollisionId); + centFT0A = recoCollision.centFT0A(); + centFT0C = recoCollision.centFT0C(); + centFT0M = recoCollision.centFT0M(); + trackOccupancyInTimeRange = recoCollision.trackOccupancyInTimeRange(); + ft0cOccupancyInTimeRange = recoCollision.ft0cOccupancyInTimeRange(); + } + + fillMCCandidate(recoCollisionId, -1, -1, + centFT0A, centFT0C, centFT0M, + trackOccupancyInTimeRange, ft0cOccupancyInTimeRange, + primVtx[0], primVtx[1], primVtx[2], + mRunNumber, 0, + -1, -1, -1, + -1, -1, -1, + -1, -1, -1, + -1, -1, -1, + -1, -1, -1, -1, -1, -1, -1, + -1, -1, -1, -1, 0, 0, 0, 0, + -1, -1, -1, false, + chargeFactor * hypCand.genPt(), hypCand.genPhi(), hypCand.genEta(), hypCand.genPtHe3(), + hypCand.gDecVtx[0], hypCand.gDecVtx[1], hypCand.gDecVtx[2], + hypCand.isReco, -1, hypCand.isSignal, hypCand.isRecoMCCollision, hypCand.isSurvEvSelection, isHeFound, mcProcess); + } + } + void processDataTracked(CollisionsFull const& collisions, aod::V0s const& V0s, aod::TrackedV0s const& tV0s, TracksFull const& tracks, aod::AmbiguousTracks const& ambiTracks, aod::BCsWithTimestamps const& bcs) { trackedClSize.clear(); @@ -893,128 +1046,16 @@ struct HyperRecoTask { void processMC(CollisionsFullMC const& collisions, aod::McCollisions const& mcCollisions, aod::V0s const& V0s, TracksFull const& tracks, aod::AmbiguousTracks const& ambiTracks, aod::BCsWithTimestamps const& bcs, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const& particlesMC) { - filledMothers.clear(); - recoCollisionIds.clear(); - recoCollisionIds.resize(mcCollisions.size(), -1); - isSurvEvSelCollision.clear(); - isSurvEvSelCollision.resize(mcCollisions.size(), false); - goodCollision.clear(); - goodCollision.resize(collisions.size(), false); - hyperCandidates.clear(); - - selectGoodCollisionsMC(collisions); - useCustomVertexer ? fillCustomV0s(collisions, tracks, ambiTracks, bcs) : fillV0s(collisions, tracks, V0s); - fillMCinfo(trackLabelsMC, particlesMC); - for (const auto& hypCand : hyperCandidates) { - auto collision = collisions.rawIteratorAt(hypCand.collisionID); - if (!hypCand.isSignal && mcSignalOnly) - continue; - int chargeFactor = -1 + 2 * (hypCand.pdgCode > 0); - float trackedHypClSize = !trackedClSize.empty() ? trackedClSize[hypCand.v0ID] : 0; - outputMCTable(hypCand.collisionID, hypCand.heTrackID, hypCand.piTrackID, - collision.centFT0A(), collision.centFT0C(), collision.centFT0M(), - collision.trackOccupancyInTimeRange(), collision.ft0cOccupancyInTimeRange(), - collision.posX(), collision.posY(), collision.posZ(), - mRunNumber, hypCand.isMatter, - hypCand.recoPtHe3(), hypCand.recoPhiHe3(), hypCand.recoEtaHe3(), - hypCand.recoPtPi(), hypCand.recoPhiPi(), hypCand.recoEtaPi(), - hypCand.decVtx[0], hypCand.decVtx[1], hypCand.decVtx[2], - hypCand.dcaV0dau, hypCand.he3DCAXY, hypCand.piDCAXY, - hypCand.nSigmaHe3, hypCand.nTPCClustersHe3, hypCand.nTPCClustersPi, hypCand.nTPCpidClusHe3, hypCand.nTPCpidClusPi, hypCand.nTPCCrossedRowsHe3, hypCand.nTPCCrossedRowsPi, - hypCand.momHe3TPC, hypCand.momPiTPC, hypCand.tpcSignalHe3, hypCand.tpcSignalPi, hypCand.tpcChi2He3, hypCand.itsChi2He3, hypCand.itsChi2Pi, - hypCand.massTOFHe3, - hypCand.clusterSizeITSHe3, hypCand.clusterSizeITSPi, hypCand.flags, trackedHypClSize, - chargeFactor * hypCand.genPt(), hypCand.genPhi(), hypCand.genEta(), hypCand.genPtHe3(), - hypCand.gDecVtx[0], hypCand.gDecVtx[1], hypCand.gDecVtx[2], - hypCand.isReco, hypCand.isFakeHeOnITSLayer, hypCand.isSignal, hypCand.isRecoMCCollision, hypCand.isSurvEvSelection, 1, 0); - } - - // now we fill only the signal candidates that were not reconstructed - for (const auto& mcPart : particlesMC) { - - if (std::abs(mcPart.pdgCode()) != hyperPdg) { - continue; - } - std::array secVtx{0.f, 0.f, 0.f}; - std::array lastDaugVtx{0.f, 0.f, 0.f}; - std::array primVtx = {mcPart.vx(), mcPart.vy(), mcPart.vz()}; - std::array momMother = {mcPart.px(), mcPart.py(), mcPart.pz()}; - std::array momHe3{0.f, 0.f, 0.f}; - bool isHeFound = false; - int mcProcess = {0}; - for (const auto& mcDaught : mcPart.daughters_as()) { - if (mcDaught.pdgCode() != PDG_t::kElectron) { // we do not care about delta electrons - lastDaugVtx = {mcDaught.vx(), mcDaught.vy(), mcDaught.vz()}; - mcProcess = mcDaught.getProcess(); - } - if (std::abs(mcDaught.pdgCode()) == heDauPdg) { - secVtx = lastDaugVtx; - momHe3 = {mcDaught.px(), mcDaught.py(), mcDaught.pz()}; - isHeFound = true; - } - } - if (mcPart.pdgCode() > 0) { - qaRegistry.fill(HIST("hIsMatterGen"), 0.); - } else { - qaRegistry.fill(HIST("hIsMatterGen"), 1.); - } - if (!isHeFound) { - qaRegistry.fill(HIST("hDecayChannel"), 1.); - } - qaRegistry.fill(HIST("hDecayChannel"), 0.); - if (mcPart.pdgCode() > 0) { - qaRegistry.fill(HIST("hIsMatterGenTwoBody"), 0.); - } else { - qaRegistry.fill(HIST("hIsMatterGenTwoBody"), 1.); - } - if (std::find(filledMothers.begin(), filledMothers.end(), mcPart.globalIndex()) != std::end(filledMothers)) { - continue; - } - hyperCandidate hypCand; - hypCand.pdgCode = mcPart.pdgCode(); - hypCand.isRecoMCCollision = recoCollisionIds[mcPart.mcCollisionId()] >= 0; - hypCand.isSurvEvSelection = isSurvEvSelCollision[mcPart.mcCollisionId()]; - int chargeFactor = -1 + 2 * (hypCand.pdgCode > 0); - for (std::size_t i = 0; i < hypCand.gDecVtx.size(); i++) { - hypCand.gDecVtx[i] = (isHeFound ? secVtx[i] : lastDaugVtx[i]) - primVtx[i]; - hypCand.gMom[i] = momMother[i]; - hypCand.gMomHe3[i] = momHe3[i]; - } - hypCand.heTrackID = -1; - hypCand.piTrackID = -1; - hypCand.isSignal = true; - - float centFT0A = -1, centFT0C = -1, centFT0M = -1; - int trackOccupancyInTimeRange = -1; - float ft0cOccupancyInTimeRange = -1.f; - if (hypCand.isRecoMCCollision) { - auto recoCollision = collisions.rawIteratorAt(recoCollisionIds[mcPart.mcCollisionId()]); - centFT0A = recoCollision.centFT0A(); - centFT0C = recoCollision.centFT0C(); - centFT0M = recoCollision.centFT0M(); - trackOccupancyInTimeRange = recoCollision.trackOccupancyInTimeRange(); - ft0cOccupancyInTimeRange = recoCollision.ft0cOccupancyInTimeRange(); - } - - outputMCTable(-1, -1, -1, - centFT0A, centFT0C, centFT0M, - trackOccupancyInTimeRange, ft0cOccupancyInTimeRange, - primVtx[0], primVtx[1], primVtx[2], - mRunNumber, 0, - -1, -1, -1, - -1, -1, -1, - -1, -1, -1, - -1, -1, -1, - -1, -1, -1, -1, -1, -1, -1, - -1, -1, -1, -1, 0, 0, 0, 0, - -1, -1, -1, false, - chargeFactor * hypCand.genPt(), hypCand.genPhi(), hypCand.genEta(), hypCand.genPtHe3(), - hypCand.gDecVtx[0], hypCand.gDecVtx[1], hypCand.gDecVtx[2], - hypCand.isReco, -1, hypCand.isSignal, hypCand.isRecoMCCollision, hypCand.isSurvEvSelection, isHeFound, mcProcess); - } + runMC(collisions, mcCollisions, V0s, tracks, ambiTracks, bcs, trackLabelsMC, particlesMC); } PROCESS_SWITCH(HyperRecoTask, processMC, "MC analysis", false); + void processMCWithCollID(CollisionsFullMC const& collisions, aod::McCollisions const& mcCollisions, aod::V0s const& V0s, TracksFull const& tracks, aod::AmbiguousTracks const& ambiTracks, aod::BCsWithTimestamps const& bcs, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const& particlesMC) + { + runMC(collisions, mcCollisions, V0s, tracks, ambiTracks, bcs, trackLabelsMC, particlesMC); + } + PROCESS_SWITCH(HyperRecoTask, processMCWithCollID, "MC analysis with collision ID", false); + void processMCTracked(CollisionsFullMC const& collisions, aod::McCollisions const& mcCollisions, aod::V0s const& V0s, aod::TrackedV0s const& tV0s, TracksFull const& tracks, aod::AmbiguousTracks const& ambiTracks, aod::BCsWithTimestamps const& bcs, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const& particlesMC) { trackedClSize.clear(); @@ -1029,17 +1070,21 @@ struct HyperRecoTask { template bool passEvtSel(const CollType& collision) { - if (!collision.sel8()) + if (!collision.sel8()) { return false; + } - if ((std::abs(collision.posZ())) > MaxAbsVertexZ) + if ((std::abs(collision.posZ())) > MaxAbsVertexZ) { return false; + } - if (cfgEvSelkNoSameBunchPileup && !collision.selection_bit(aod::evsel::kNoSameBunchPileup)) + if (cfgEvSelkNoSameBunchPileup && !collision.selection_bit(aod::evsel::kNoSameBunchPileup)) { return false; + } - if (cfgEvSelkIsGoodZvtxFT0vsPV && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) + if (cfgEvSelkIsGoodZvtxFT0vsPV && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { return false; + } return true; } @@ -1101,10 +1146,12 @@ struct HyperRecoTask { ROOT::Math::PxPyPzMVector daugh1, daugh2, mother; for (const auto& genParticle : GenParticles) { - if (std::abs(genParticle.y()) > 1) + if (std::abs(genParticle.y()) > 1) { continue; - if (std::abs(genParticle.pdgCode()) != hyperPdg) + } + if (std::abs(genParticle.pdgCode()) != hyperPdg) { continue; + } auto daughters = genParticle.daughters_as(); @@ -1127,8 +1174,9 @@ struct HyperRecoTask { } } // Check pairs to avoid wrong charge associations - if (!((dauHe3 && dauPiMinus) || !(dauAntiHe3 && dauPiPos))) + if (!((dauHe3 && dauPiMinus) || !(dauAntiHe3 && dauPiPos))) { continue; + } mother = daugh1 + daugh2; diff --git a/PWGLF/TableProducer/Nuspex/piHypertritonFemto.cxx b/PWGLF/TableProducer/Nuspex/piHypertritonFemto.cxx new file mode 100644 index 00000000000..ce1433b873b --- /dev/null +++ b/PWGLF/TableProducer/Nuspex/piHypertritonFemto.cxx @@ -0,0 +1,1222 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +// + +/// \file piHypertritonFemto.cxx +/// \brief LF pi-hypertriton femto pair table producer Shifted from PWGCF +/// \author zhengqing and meiyi +/// \date 2026-09-17 + +#include "PWGLF/DataModel/LFHypernucleiTables.h" +#include "PWGLF/DataModel/LFPiHypertritonFemtoTables.h" + +#include "Common/CCDB/EventSelectionParams.h" +#include "Common/Core/RecoDecay.h" +#include "Common/Core/Zorro.h" +#include "Common/Core/ZorroSummary.h" +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/PIDResponseTOF.h" +#include "Common/DataModel/PIDResponseTPC.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +using HyperCandidates = aod::DataHypCandsWColl; +using HyperCandidatesMC = aod::MCHypCandsWColl; +using HadHyperCollisionsFull = soa::Join; +using HadHyperCollisionsFullMC = soa::Join; +using TrackCandidates = soa::Join; +using TrackCandidatesMC = soa::Join; + +namespace +{ +constexpr std::array InvalidMomentum{-1.f, -1.f, -1.f}; +constexpr int HyperTritonPDG = o2::constants::physics::Pdg::kHyperTriton; +// FemtoDream average-phi-star reference radii, in cm. +constexpr std::array CPRTPCRadii{85.f, 105.f, 125.f, 145.f, 165.f, 185.f, 205.f, 225.f, 245.f}; +using PairLorentzVector = ROOT::Math::LorentzVector>; +} // namespace + +struct PiHypertritonFemto { + // First aggregate constructor argument, supplied by defineDataProcessing. + std::shared_ptr> mEndOfStreamCallback; + Produces mOutputDataTable; + Produces mOutputMCTable; + + struct : o2::framework::ConfigurableGroup { + // cppcheck-suppress unusedStructMember + std::string prefix{"eventMixing"}; + Configurable settingCutVertex{"settingCutVertex", 10.0f, "Accepted z-vertex range"}; + Configurable settingNoMixedEvents{"settingNoMixedEvents", 5, "Number of mixed events per event"}; + } eventMixing; + + struct : o2::framework::ConfigurableGroup { + // cppcheck-suppress unusedStructMember + std::string prefix{"eventSelection"}; + Configurable disableITSROFCut{"disableITSROFCut", false, "Disable the ITS readout-frame border cut for data, as in hyperRecoTask"}; + Configurable cfgEvSelkNoSameBunchPileup{"cfgEvSelkNoSameBunchPileup", false, "Reject collisions sharing the same found-by-T0 bunch crossing"}; + Configurable cfgEvSelkIsGoodZvtxFT0vsPV{"cfgEvSelkIsGoodZvtxFT0vsPV", false, "Require consistent primary vertex z positions from tracking and FT0 timing"}; + } eventSelection; + + struct : o2::framework::ConfigurableGroup { + // cppcheck-suppress unusedStructMember + std::string prefix{"pionTrack"}; + Configurable settingCutEta{"settingCutEta", 0.8f, "Maximum pion track |eta|"}; + Configurable settingPtMin{"settingPtMin", 0.14f, "Minimum pion pT"}; + Configurable settingPtMax{"settingPtMax", 2.5f, "Maximum pion pT"}; + Configurable settingITSInnerBarrelMin{"settingITSInnerBarrelMin", 3, "Minimum pion ITS inner-barrel clusters"}; + Configurable settingITSNClsMin{"settingITSNClsMin", 7, "Minimum pion ITS clusters"}; + Configurable settingTPCNClsFoundMin{"settingTPCNClsFoundMin", 80, "Minimum pion found TPC clusters"}; + Configurable settingTPCCrossedRowsMin{"settingTPCCrossedRowsMin", 90, "Minimum pion crossed TPC rows"}; + Configurable settingDCAxyOffset{"settingDCAxyOffset", 0.004f, "Pion DCAxy offset"}; + Configurable settingDCAxyPtCoeff{"settingDCAxyPtCoeff", 0.013f, "Pion DCAxy 1/pT coefficient"}; + Configurable settingDCAzOffset{"settingDCAzOffset", 0.004f, "Pion DCAz offset"}; + Configurable settingDCAzPtCoeff{"settingDCAzPtCoeff", 0.013f, "Pion DCAz 1/pT coefficient"}; + } pionTrack; + + struct : o2::framework::ConfigurableGroup { + // cppcheck-suppress unusedStructMember + std::string prefix{"pionPid"}; + Configurable settingMomCombMin{"settingMomCombMin", 0.5f, "Minimum momentum to use combined TPC+TOF pion PID"}; + Configurable settingTPCNsigMax{"settingTPCNsigMax", 3.0f, "Maximum pion TPC n-sigma below the TOF threshold"}; + Configurable settingCombNsigMax{"settingCombNsigMax", 3.0f, "Maximum combined pion TPC+TOF n-sigma"}; + Configurable settingReqSingleNsig{"settingReqSingleNsig", false, "Also require individual TPC and TOF n-sigma cuts in combined PID"}; + } pionPid; + + struct : o2::framework::ConfigurableGroup { + // cppcheck-suppress unusedStructMember + std::string prefix{"CPR"}; + Configurable settingClosePairDeltaPhiMax{"settingClosePairDeltaPhiMax", 0.01f, "Average delta-phi-star ellipse semiaxis for the offline close-pair flag"}; + Configurable settingClosePairDeltaEtaMax{"settingClosePairDeltaEtaMax", 0.01f, "Delta-eta ellipse semiaxis for the offline close-pair flag"}; + Configurable settingClosePairDistanceMax{"settingClosePairDistanceMax", 8.f, "Average 3D TPC separation below which a same-sign daughter pair is tagged (cm)"}; + } CPR; + + struct : o2::framework::ConfigurableGroup { + // cppcheck-suppress unusedStructMember + std::string prefix{"mc"}; + Configurable settingRequireSel8{"settingRequireSel8", false, "Additionally require sel8 for reconstructed MC; TVX and timeframe border cuts always apply"}; + Configurable settingRequireRecoMCCollisionMatch{"settingRequireRecoMCCollisionMatch", true, "Require reconstructed collision MC labels"}; + } mc; + + struct : o2::framework::ConfigurableGroup { + // cppcheck-suppress unusedStructMember + std::string prefix{"hypertriton"}; + Configurable settingHypMassMin{"settingHypMassMin", 2.94f, "Minimum hypertriton invariant mass"}; + Configurable settingHypMassMax{"settingHypMassMax", 3.10f, "Maximum hypertriton invariant mass"}; + } hypertriton; + + struct : o2::framework::ConfigurableGroup { + // cppcheck-suppress unusedStructMember + std::string prefix{"output"}; + Configurable settingFillTable{"settingFillTable", false, "Enable output table filling"}; + } output; + + struct : o2::framework::ConfigurableGroup { + // cppcheck-suppress unusedStructMember + std::string prefix{"hadHyper"}; + Configurable enableMixing{"enableMixing", true, "Build mixed-event pion-hypertriton pairs"}; + Configurable maxOutputKstar{"maxOutputKstar", -1.f, "Maximum pair k* (GeV/c); negative saves all selected pairs"}; + } hadHyper; + + struct : o2::framework::ConfigurableGroup { + // cppcheck-suppress unusedStructMember + std::string prefix{"zorro"}; + Configurable settingSkimmedProcessing{"settingSkimmedProcessing", false, "Skimmed dataset processing"}; + Configurable settingTriggerMask{"settingTriggerMask", "fPiHypertritonFemto", "Zorro trigger mask"}; + } zorro; + + struct : o2::framework::ConfigurableGroup { + // cppcheck-suppress unusedStructMember + std::string prefix{"ccdb"}; + Configurable settingCcdburl{"settingCcdburl", "http://alice-ccdb.cern.ch", "URL of the CCDB repository used by Zorro and the magnetic field"}; + Configurable settingGrpmagPath{"settingGrpmagPath", "GLO/Config/GRPMagField", "CCDB path of the Run 3 magnetic field"}; + } ccdb; + + Preslice mPerCol = aod::track::collisionId; + Preslice mPerColMC = aod::track::collisionId; + PresliceUnsorted hypPerCol = o2::aod::hyperrec::collisionId; + PresliceUnsorted hypPerColMC = o2::aod::hyperrec::collisionId; + + ConfigurableAxis axisVertex{"axisVertex", {30, -10, 10}, "Binning for mixed-event vertex z"}; + ConfigurableAxis axisCentrality{"axisCentrality", {40, 0, 100}, "Binning for mixed-event centrality"}; + ConfigurableAxis axisPairMultiplicity{"axisPairMultiplicity", {1501, -0.5, 1500.5}, "Binning for accepted pairs per hypertriton candidate (SE and ME)"}; + using BinningType = ColumnBinningPolicy; + + using HadHyperTrackInfo = std::tuple; + using HadHyperEventInfo = std::tuple; + using HadHyperCandidateInfo = std::tuple; + using HadHyperDataInfo = decltype(std::tuple_cat(std::declval>(), + std::declval(), + std::declval(), + std::declval())); + using HadHyperMCInfo = std::tuple; + + struct HadHyperParticleTruth { + int64_t particleId{-1}; + int64_t collisionId{-1}; + int32_t pdgCode{0}; + float pt{-1.f}; + float eta{-999.f}; + float phi{-999.f}; + bool isPhysicalPrimary{false}; + int16_t statusCode{0}; + int16_t process{0}; + std::array momentum{InvalidMomentum}; + }; + + struct ClosePairTrack { + float eta{0.f}; + int8_t sign{0}; + std::array phiStar{}; + std::array valid{}; + std::array, CPRTPCRadii.size()> positions{}; // Relative to the source collision PV. + std::array positionValid{}; + }; + + struct ClosePairResult { + float deltaEta{0.f}; + float deltaPhiStar{0.f}; + int validRadii{0}; + float averageDistance{0.f}; + int validPositionRadii{0}; + bool angularRejected{false}; + bool distanceRejected{false}; + bool angularUnavailable{false}; + bool distanceUnavailable{false}; + }; + + struct ClosePairQA { + std::shared_ptr before; + std::shared_ptr after; + std::shared_ptr validRadii; + std::shared_ptr distanceBefore; + std::shared_ptr distanceAfter; + std::shared_ptr validPositionRadii; + std::shared_ptr decisions; + }; + + struct HadHyperCandidate { + HadHyperCandidateInfo info{}; + HadHyperParticleTruth hyperTruth{}; + HadHyperParticleTruth heTruth{}; + HadHyperParticleTruth decayPionTruth{}; + int64_t heTrackId{-1}; + int64_t pionTrackId{-1}; + int16_t statusCode{0}; + bool isReco{true}; + bool isSignal{false}; + bool isRecoMCCollision{false}; + bool isSurvEvSel{false}; + bool isTwoBodyDecay{false}; + float genPt{-1.f}; + float genEta{-1.f}; + float genPhi{-1.f}; + float genPtHe3{-1.f}; + std::array genDecVtx{InvalidMomentum}; + std::array momentum{}; + std::array genMomentum{InvalidMomentum}; + float mass{0.f}; + bool isMatter{false}; + std::array daughterClosePairTracks{}; // He3, decay pion. + uint64_t mixedPairCount{0}; // Accumulated over the complete residence in the mixing pool. + + [[nodiscard]] float pt() const { return std::hypot(momentum[0], momentum[1]); } + [[nodiscard]] float eta() const + { + const float transverseMomentum = pt(); + return transverseMomentum > 0.f ? std::asinh(momentum[2] / transverseMomentum) : 999.f; + } + [[nodiscard]] float phi() const { return std::atan2(momentum[1], momentum[0]); } + [[nodiscard]] int8_t sign() const { return isMatter ? 1 : -1; } + }; + + struct HadHyperHadron { + HadHyperTrackInfo info{}; + HadHyperParticleTruth truth{}; + std::vector motherIds; + int64_t sourceId{-1}; + std::array momentum{}; + int8_t signValue{0}; + ClosePairTrack closePairTrack{}; + + [[nodiscard]] float eta() const { return std::get<1>(info); } + [[nodiscard]] float phi() const { return std::get<2>(info); } + [[nodiscard]] int8_t sign() const { return signValue; } + }; + + struct HadHyperEvent { + HadHyperEventInfo info{}; + uint64_t sourceFrameId{0}; // Task-local scope for track and MC indices stored in the mixing pool. + int64_t mcCollisionId{-1}; + bool hasMCCollision{false}; + float centrality{0.f}; + std::vector candidates; + std::vector hadrons; + }; + + std::unordered_map> mMixingPools; + std::unordered_map> mMCMixingPools; + int mMixingRunNumber{-1}; + int mMCMixingRunNumber{-1}; + uint64_t mSourceFrameId{0}; + int mRunNumber{-1}; + float mMagneticFieldTesla{0.f}; + std::array, 2> mClosePairQA{}; // SE/ME, He3/decay pion. + Service mCcdb{}; + Zorro mZorro; + OutputObj mZorroSummary{"zorroSummary"}; + + HistogramRegistry mQaRegistry{ + "QA", + {{"hEvents", "Event selection cut-flow;Selection step;Events", {HistType::kTH1F, {{11, -0.5, 10.5}}}}, + {"hVtxZ", "Vertex distribution in Z;Z (cm);Entries", {HistType::kTH1F, {{400, -20.0, 20.0}}}}, + {"hNcontributor", "Number of primary vertex contributors;contributors;Entries", {HistType::kTH1F, {{2000, 0.0f, 2000.0f}}}}, + {"hPionTPCNSigmaPreselection", "Pion TPC n#sigma before PID;signed #it{p}_{TPC};n#sigma_{TPC}^{#pi}", {HistType::kTH2F, {{200, -5.0f, 5.0f}, {400, -10.0f, 10.0f}}}}, + {"hPionTOFNSigmaPreselection", "Pion TOF n#sigma before PID;signed #it{p}_{T};n#sigma_{TOF}^{#pi}", {HistType::kTH2F, {{280, -7.0f, 7.0f}, {400, -10.0f, 10.0f}}}}, + {"hPionCombNSigmaPreselection", "Pion combined n#sigma before PID;signed #it{p}_{T};n#sigma_{comb}^{#pi}", {HistType::kTH2F, {{280, -7.0f, 7.0f}, {100, 0.0f, 5.0f}}}}, + {"hPionTPCNSigma", "Selected pion TPC n#sigma;signed #it{p}_{T};n#sigma_{TPC}^{#pi}", {HistType::kTH2F, {{280, -7.0f, 7.0f}, {200, -5.0f, 5.0f}}}}, + {"hPionTOFNSigma", "Selected pion TOF n#sigma;signed #it{p}_{T};n#sigma_{TOF}^{#pi}", {HistType::kTH2F, {{280, -7.0f, 7.0f}, {200, -5.0f, 5.0f}}}}, + {"hPionCombNSigma", "Selected pion combined n#sigma;signed #it{p}_{T};n#sigma_{comb}^{#pi}", {HistType::kTH2F, {{280, -7.0f, 7.0f}, {100, 0.0f, 5.0f}}}}, + {"hPionTPC", "Selected pion TPC d#it{E}/d#it{x};signed #it{p}_{TPC};TPC signal", {HistType::kTH2F, {{200, -5.0f, 5.0f}, {100, 0.0f, 2000.0f}}}}, + {"hHypHe3TPCNSigma", "Selected hypertriton candidates;#it{p}_{T}^{^{3}He};n#sigma_{TPC}^{^{3}He}", {HistType::kTH2F, {{280, -7.0f, 7.0f}, {200, -5.0f, 5.0f}}}}, + {"hHypHe3TPCMomPreselection", "Hypertriton candidate He3 TPC momentum before mass cut;#it{p}_{TPC};Entries", {HistType::kTH1F, {{120, -3.0f, 3.0f}}}}, + {"hHypHe3TPCMom", "Hypertriton candidate He3 TPC momentum after mass cut;#it{p}_{TPC};Entries", {HistType::kTH1F, {{120, -3.0f, 3.0f}}}}}, + OutputObjHandlingPolicy::AnalysisObject, + false, + true}; + + HistogramRegistry hadHyperRegistry{ + "hadHyperRegistry", + {{"hSE", "Raw same-event pairs;k*;Entries", {HistType::kTH1D, {{300, 0., 3.}}}}, + {"hME", "Raw mixed-event pairs;k*;Entries", {HistType::kTH1D, {{300, 0., 3.}}}}, + {"hMass", "Unique selected candidates;mass He3-pi;Entries", {HistType::kTH1D, {{160, 2.94, 3.10}}}}, + {"hSelfPairs", "Rejected daughter reuse;reason (0=He3,1=pion,2=truth);k*", {HistType::kTH2F, {{3, -0.5, 2.5}, {300, 0., 3.}}}}, + {"hSameEventSelfPairs", "Same-event rejected daughter reuse;reason (0=He3,1=pion,2=truth);k*", {HistType::kTH2F, {{3, -0.5, 2.5}, {300, 0., 3.}}}}, + {"hMixEventDeltaPosZVsCent", "Mixed-event #Delta z_{vtx} vs hypertriton centrality;hypertriton CentFT0C;#Delta z_{vtx}", {HistType::kTH2F, {{100, 0., 100.}, {120, -30., 30.}}}}, + {"hMixEventDeltaCentFT0CVsCent", "Mixed-event #Delta CentFT0C vs hypertriton centrality;hypertriton CentFT0C;#Delta CentFT0C", {HistType::kTH2F, {{100, 0., 100.}, {200, -100., 100.}}}}, + {"hMixingDepth", "Available partner events;events;anchor events", {HistType::kTH1F, {{101, -0.5, 100.5}}}}, + {"hPoolFlow", "Mixing pool selection;0=selected,1=in range;events", {HistType::kTH1F, {{2, -0.5, 1.5}}}}, + {"hDaughterHeTPC", "He3 daughter QA per candidate;TPC rigidity;TPC signal", {HistType::kTH2F, {{100, 0., 10.}, {300, 0., 1500.}}}}, + {"hDaughterPiTPC", "Pion daughter QA per candidate;TPC rigidity;TPC signal", {HistType::kTH2F, {{100, 0., 10.}, {300, 0., 1500.}}}}, + {"hPionTPC", "Pair pion QA per selected track;TPC rigidity;TPC signal", {HistType::kTH2F, {{100, 0., 10.}, {300, 0., 1500.}}}}, + {"MC/SE/hKstarRecVsGenHyperReco", "Pion-hypertriton k* response for reconstructed true hypertritons;generated k* (GeV/#it{c});reconstructed k* (GeV/#it{c})", {HistType::kTH2F, {{300, 0., 3.}, {300, 0., 3.}}}}, + {"MC/SE/hKstarResolutionHyperReco", "Pion-hypertriton k* resolution for reconstructed true hypertritons;reconstructed k* (GeV/#it{c});k*_{reco}-k*_{gen} (GeV/#it{c})", {HistType::kTH2F, {{300, 0., 3.}, {200, -0.2, 0.2}}}}, + {"MC/SE/hPrimaryPionVsKstarDen", "Truth-matched selected pion pair-weighted denominator;k* (GeV/#it{c});Entries", {HistType::kTH1D, {{300, 0., 3.}}}}, + {"MC/SE/hPrimaryPionVsKstarNum", "Physical-primary selected pion pair-weighted numerator;k* (GeV/#it{c});Entries", {HistType::kTH1D, {{300, 0., 3.}}}}, + {"MC/SE/hPrimaryPionVsCentDen", "Truth-matched selected pion pair-weighted denominator;hypertriton CentFT0C;Entries", {HistType::kTH1D, {{100, 0., 100.}}}}, + {"MC/SE/hPrimaryPionVsCentNum", "Physical-primary selected pion pair-weighted numerator;hypertriton CentFT0C;Entries", {HistType::kTH1D, {{100, 0., 100.}}}}}}; + + void init(o2::framework::InitContext&); + template + void initCCDB(const aod::BCsWithTimestamps::iterator& bc); + template + bool selectCollision(const Tcollision& collision, const aod::BCsWithTimestamps& bcs); + template + bool selectPionTrack(const Ttrack& candidate) const; + template + bool selectionPIDPion(const Ttrack& candidate); + template + bool selectHyperCandidate(const Tcandidate& candidate); + template + void processPairs(const Tcollisions& collisions, const Ttracks& tracks, const Tcandidates& candidates, + const aod::BCsWithTimestamps& bcs, + std::unordered_map>& mixingPools, + int& mixingRunNumber); + + template + float computeHyperCandidateMass(const Tcandidate& candidate) const + { + const std::array heMomentum{candidate.ptHe3() * std::cos(candidate.phiHe3()), candidate.ptHe3() * std::sin(candidate.phiHe3()), candidate.ptHe3() * std::sinh(candidate.etaHe3())}; + const std::array pionMomentum{candidate.ptPi() * std::cos(candidate.phiPi()), candidate.ptPi() * std::sin(candidate.phiPi()), candidate.ptPi() * std::sinh(candidate.etaPi())}; + return RecoDecay::m(std::array{heMomentum, pionMomentum}, + std::array{static_cast(o2::constants::physics::MassHelium3), + static_cast(o2::constants::physics::MassPiPlus)}); + } + + float wrapDeltaPhi(float dphi) const { return std::atan2(std::sin(dphi), std::cos(dphi)); } + + template + ClosePairTrack makeClosePairTrack(const Ttrack& track, const std::array& primaryVertex) const + { + ClosePairTrack result; + result.eta = track.eta(); + result.sign = track.sign(); + if (!std::isfinite(result.eta) || !std::isfinite(track.phi()) || !std::isfinite(track.signed1Pt()) || track.signed1Pt() == 0.f) { + return result; + } + // FemtoDream helix approximation, using the original daughter track, not the mother. + // AO2D signed1Pt is signed inverse rigidity: He3's |Z|=2 is already included. + // Do not multiply it by two again. Cache each source event's field-dependent values. + constexpr float CurvatureFactor = 0.3f * 0.01f / 2.f; // B in T, radius in cm. + for (std::size_t i = 0; i < CPRTPCRadii.size(); ++i) { + const float argument = CurvatureFactor * mMagneticFieldTesla * CPRTPCRadii[i] * track.signed1Pt(); + if (std::isfinite(argument) && std::abs(argument) < 1.f) { + result.phiStar[i] = track.phi() - std::asin(argument); + result.valid[i] = true; + } + } + // Uniform-field extrapolation from the measured track state, retaining daughter displacement. + // This is not material-aware transport. getTrackPar keeps q/pT (including He3 rigidity). + const auto trackPar = getTrackPar(track); + constexpr float TeslaToKilogauss = 10.f; + constexpr float TPCHalfLength = 250.f; // cm + const float fieldKilogauss = TeslaToKilogauss * mMagneticFieldTesla; + for (std::size_t i = 0; i < CPRTPCRadii.size(); ++i) { + float localX = 0.f; + if (!trackPar.getXatLabR(CPRTPCRadii[i], localX, fieldKilogauss, o2::track::DirOutward)) { + continue; + } + bool propagated = false; + const auto point = trackPar.getXYZGloAt(localX, fieldKilogauss, propagated); + if (!propagated || !std::isfinite(point.X()) || !std::isfinite(point.Y()) || !std::isfinite(point.Z()) || std::abs(point.Z()) > TPCHalfLength) { + continue; + } + // Apply the same source-PV convention to SE and ME; in SE the common PV cancels. + result.positions[i] = {point.X() - primaryVertex[0], point.Y() - primaryVertex[1], point.Z() - primaryVertex[2]}; + result.positionValid[i] = true; + } + return result; + } + + ClosePairResult evaluateClosePair(const ClosePairTrack& pion, const ClosePairTrack& daughter) const + { + ClosePairResult result; + result.deltaEta = pion.eta - daughter.eta; + for (std::size_t i = 0; i < CPRTPCRadii.size(); ++i) { + if (pion.valid[i] && daughter.valid[i]) { + result.deltaPhiStar += wrapDeltaPhi(pion.phiStar[i] - daughter.phiStar[i]); + ++result.validRadii; + } + if (pion.positionValid[i] && daughter.positionValid[i]) { + const auto& pionPosition = pion.positions[i]; + const auto& daughterPosition = daughter.positions[i]; + result.averageDistance += std::hypot(pionPosition[0] - daughterPosition[0], + pionPosition[1] - daughterPosition[1], + pionPosition[2] - daughterPosition[2]); + ++result.validPositionRadii; + } + } + const bool sameSign = pion.sign * daughter.sign > 0; + result.angularUnavailable = sameSign && result.validRadii == 0; + result.distanceUnavailable = sameSign && result.validPositionRadii == 0; + // Undefined comparisons set separate status bits, never a rejection bit. + if (result.validRadii > 0) { + result.deltaPhiStar /= result.validRadii; + const float scaledEta = result.deltaEta / CPR.settingClosePairDeltaEtaMax.value; + const float scaledPhi = result.deltaPhiStar / CPR.settingClosePairDeltaPhiMax.value; + result.angularRejected = sameSign && scaledEta * scaledEta + scaledPhi * scaledPhi < 1.f; + } + if (result.validPositionRadii > 0) { + result.averageDistance /= result.validPositionRadii; + result.distanceRejected = sameSign && result.averageDistance < CPR.settingClosePairDistanceMax.value; + } + return result; + } + + uint8_t flagClosePair(const HadHyperCandidate& candidate, const HadHyperHadron& pion, bool mixed) + { + const std::array results{ + evaluateClosePair(pion.closePairTrack, candidate.daughterClosePairTracks[0]), + evaluateClosePair(pion.closePairTrack, candidate.daughterClosePairTracks[1])}; + uint8_t rejectionFlags = 0; + for (const auto& result : results) { + if (result.angularRejected) { + rejectionFlags |= aod::pihypertritonfemto::ClosePairAngular; + } + if (result.distanceRejected) { + rejectionFlags |= aod::pihypertritonfemto::ClosePairDistance; + } + if (result.angularUnavailable) { + rejectionFlags |= aod::pihypertritonfemto::ClosePairAngularUnavailable; + } + if (result.distanceUnavailable) { + rejectionFlags |= aod::pihypertritonfemto::ClosePairDistanceUnavailable; + } + } + for (std::size_t i = 0; i < results.size(); ++i) { + const auto& result = results[i]; + const auto& qa = mClosePairQA[mixed ? 1 : 0][i]; + qa.validRadii->Fill(result.validRadii); + qa.validPositionRadii->Fill(result.validPositionRadii); + if (result.validRadii > 0 && result.validPositionRadii > 0) { + qa.decisions->Fill(result.angularRejected, result.distanceRejected); + } + if (result.validPositionRadii > 0) { + qa.distanceBefore->Fill(result.averageDistance); + if (rejectionFlags == 0) { + qa.distanceAfter->Fill(result.averageDistance); + } + } + if (result.validRadii > 0) { + qa.before->Fill(result.deltaEta, result.deltaPhiStar); + if (rejectionFlags == 0) { + qa.after->Fill(result.deltaEta, result.deltaPhiStar); + } + } + } + return rejectionFlags; + } + + float computePairKstar(const std::array& momPion, const std::array& momHyper) const + { + const PairLorentzVector vecPion(momPion[0], momPion[1], momPion[2], o2::constants::physics::MassPiPlus); + const PairLorentzVector vecHyper(momHyper[0], momHyper[1], momHyper[2], o2::constants::physics::MassHyperTriton); + const PairLorentzVector trackSum = vecPion + vecHyper; + const float beta = trackSum.Beta(); + const float betax = beta * std::cos(trackSum.Phi()) * std::sin(trackSum.Theta()); + const float betay = beta * std::sin(trackSum.Phi()) * std::sin(trackSum.Theta()); + const float betaz = beta * std::cos(trackSum.Theta()); + PairLorentzVector pionCMS(vecPion); + PairLorentzVector hyperCMS(vecHyper); + const ROOT::Math::Boost boostPRF = ROOT::Math::Boost(-betax, -betay, -betaz); + pionCMS = boostPRF(pionCMS); + hyperCMS = boostPRF(hyperCMS); + return 0.5f * (pionCMS - hyperCMS).P(); + } + + template + HadHyperTrackInfo makePionInfo(const Ttrack& track) + { + constexpr float InvalidPID = -999.f; + return {track.pt(), track.eta(), track.phi(), static_cast(track.sign()), + track.dcaXY(), track.tpcNClsCrossedRows(), track.tpcNClsPID(), track.tpcChi2NCl(), + track.itsClusterSizes(), track.itsChi2NCl(), track.hasTOF(), + track.tpcNSigmaPi(), track.hasTOF() ? track.tofNSigmaPi() : InvalidPID}; + } + + template + HadHyperEventInfo makeEventInfo(const Tcollision& collision) + { + return {collision.posZ(), collision.centFT0C(), + collision.trackOccupancyInTimeRange(), collision.ft0cOccupancyInTimeRange(), + collision.multFT0C(), + collision.posX(), collision.posY(), collision.posZ()}; + } + + template + HadHyperParticleTruth makeParticleTruth(const Tparticle& particle) + { + return {particle.globalIndex(), particle.mcCollisionId(), particle.pdgCode(), particle.pt(), particle.eta(), particle.phi(), particle.isPhysicalPrimary(), static_cast(particle.statusCode()), static_cast(particle.getProcess()), {particle.px(), particle.py(), particle.pz()}}; + } + + template + HadHyperHadron makePion(const Ttrack& track, const std::array& primaryVertex) + { + HadHyperHadron pion{makePionInfo(track), {}, {}, track.globalIndex(), {track.px(), track.py(), track.pz()}, static_cast(track.sign())}; + pion.closePairTrack = makeClosePairTrack(track, primaryVertex); + if constexpr (isMC) { + if (track.has_mcParticle()) { + const auto particle = track.template mcParticle_as(); + pion.truth = makeParticleTruth(particle); + if (particle.has_mothers()) { + for (const auto& mother : particle.template mothers_as()) { + pion.motherIds.push_back(mother.globalIndex()); + } + } + } + } + return pion; + } + + template + HadHyperCandidate makeHyperCandidate(const Tcandidate& candidate, const Ttracks& tracks, const std::array& primaryVertex) + { + const auto he = tracks.rawIteratorAt(candidate.heTrackId() - tracks.offset()); + const auto decayPion = tracks.rawIteratorAt(candidate.piTrackId() - tracks.offset()); + const std::array heMomentum{candidate.ptHe3() * std::cos(candidate.phiHe3()), candidate.ptHe3() * std::sin(candidate.phiHe3()), candidate.ptHe3() * std::sinh(candidate.etaHe3())}; + const std::array decayPionMomentum{candidate.ptPi() * std::cos(candidate.phiPi()), candidate.ptPi() * std::sin(candidate.phiPi()), candidate.ptPi() * std::sinh(candidate.etaPi())}; + + HadHyperCandidate result; + for (std::size_t i = 0; i < result.momentum.size(); ++i) { + result.momentum[i] = heMomentum[i] + decayPionMomentum[i]; + } + result.mass = computeHyperCandidateMass(candidate); + result.heTrackId = candidate.heTrackId(); + result.pionTrackId = candidate.piTrackId(); + result.isMatter = candidate.isMatter(); + result.daughterClosePairTracks = {makeClosePairTrack(he, primaryVertex), makeClosePairTrack(decayPion, primaryVertex)}; + + if constexpr (isMC) { + if (he.has_mcParticle()) { + result.heTruth = makeParticleTruth(he.template mcParticle_as()); + } + if (decayPion.has_mcParticle()) { + result.decayPionTruth = makeParticleTruth(decayPion.template mcParticle_as()); + } + if (he.has_mcParticle() && decayPion.has_mcParticle()) { + const auto heParticle = he.template mcParticle_as(); + const auto decayPionParticle = decayPion.template mcParticle_as(); + if (heParticle.has_mothers() && decayPionParticle.has_mothers()) { + for (const auto& heMother : heParticle.template mothers_as()) { + for (const auto& pionMother : decayPionParticle.template mothers_as()) { + if (heMother.globalIndex() == pionMother.globalIndex() && std::abs(heMother.pdgCode()) == HyperTritonPDG) { + result.hyperTruth = makeParticleTruth(heMother); + } + } + } + } + } + } + + result.info = {candidate.isMatter(), candidate.ptHe3(), candidate.etaHe3(), candidate.phiHe3(), + candidate.ptPi(), candidate.etaPi(), candidate.phiPi(), + candidate.dcaV0Daug(), candidate.dcaHe(), candidate.dcaPi(), + candidate.nSigmaHe(), + candidate.nTPCCrossedRowsHe(), candidate.nTPCCrossedRowsPi(), + candidate.tpcMomHe(), candidate.tpcMomPi(), candidate.tpcSignalHe(), candidate.tpcSignalPi(), + candidate.itsClusterSizesHe(), candidate.itsClusterSizesPi(), + candidate.xDecVtx(), candidate.yDecVtx(), candidate.zDecVtx()}; + + if constexpr (isMC) { + result.statusCode = static_cast(candidate.statusCode()); + result.isReco = candidate.isReco(); + result.isSignal = candidate.isSignal(); + result.isRecoMCCollision = candidate.isRecoMCCollision(); + result.isSurvEvSel = candidate.isSurvEvSel(); + result.isTwoBodyDecay = candidate.isTwoBodyDecay(); + result.genPt = candidate.genPt(); + result.genEta = candidate.genEta(); + result.genPhi = candidate.genPhi(); + result.genPtHe3 = candidate.genPtHe3(); + result.genDecVtx = {candidate.genXDecVtx(), candidate.genYDecVtx(), candidate.genZDecVtx()}; + const float absGenPt = std::abs(candidate.genPt()); + result.genMomentum = {absGenPt * std::cos(candidate.genPhi()), absGenPt * std::sin(candidate.genPhi()), absGenPt * std::sinh(candidate.genEta())}; + } + + hadHyperRegistry.fill(HIST("hMass"), result.mass); + hadHyperRegistry.fill(HIST("hDaughterHeTPC"), he.tpcInnerParam(), he.tpcSignal()); + hadHyperRegistry.fill(HIST("hDaughterPiTPC"), decayPion.tpcInnerParam(), decayPion.tpcSignal()); + return result; + } + + HadHyperDataInfo makeDataInfo(const HadHyperCandidate& candidate, const HadHyperHadron& pion, const HadHyperEvent& hyperEvent, bool mixed, uint8_t closePairRejected, int mixingDepth) + { + return std::tuple_cat(std::make_tuple(mixed, closePairRejected, mixingDepth), hyperEvent.info, candidate.info, pion.info); + } + + bool hasTruthMother(const HadHyperHadron& pion, int64_t motherId) const + { + return motherId >= 0 && std::find(pion.motherIds.begin(), pion.motherIds.end(), motherId) != pion.motherIds.end(); + } + + bool isRecoSelfCorrelation(const HadHyperCandidate& candidate, const HadHyperHadron& pion, + const HadHyperEvent& hyperEvent, const HadHyperEvent& pionEvent) const + { + return hyperEvent.sourceFrameId == pionEvent.sourceFrameId && pion.sourceId >= 0 && + (pion.sourceId == candidate.heTrackId || pion.sourceId == candidate.pionTrackId); + } + + bool isTruthSelfCorrelation(const HadHyperCandidate& candidate, const HadHyperHadron& pion, + const HadHyperEvent& hyperEvent, const HadHyperEvent& pionEvent) const + { + return hyperEvent.sourceFrameId == pionEvent.sourceFrameId && pion.truth.particleId >= 0 && + (pion.truth.particleId == candidate.heTruth.particleId || + pion.truth.particleId == candidate.decayPionTruth.particleId || + hasTruthMother(pion, candidate.hyperTruth.particleId)); + } + + HadHyperMCInfo makeMCInfo(const HadHyperCandidate& candidate, const HadHyperHadron& pion, + const HadHyperEvent& hyperEvent, const HadHyperEvent& pionEvent) const + { + const bool sameMCCollision = hyperEvent.sourceFrameId == pionEvent.sourceFrameId && + candidate.hyperTruth.collisionId >= 0 && pion.truth.collisionId == candidate.hyperTruth.collisionId; + const bool matchesHypRecoMCCollision = hyperEvent.hasMCCollision && candidate.hyperTruth.collisionId == hyperEvent.mcCollisionId; + const bool matchesPairRecoMCCollision = pionEvent.hasMCCollision && pion.truth.collisionId == pionEvent.mcCollisionId; + const bool truthSelfCorrelation = isTruthSelfCorrelation(candidate, pion, hyperEvent, pionEvent); + const bool hadIsTruePion = pion.truth.particleId >= 0 && std::abs(pion.truth.pdgCode) == static_cast(PDG_t::kPiPlus); + const bool isTruePrimaryPiHyperPair = candidate.isSignal && + std::abs(candidate.hyperTruth.pdgCode) == HyperTritonPDG && + hadIsTruePion && + pion.truth.isPhysicalPrimary && sameMCCollision && !truthSelfCorrelation; + return {candidate.genPt, candidate.genEta, candidate.genPhi, + candidate.genDecVtx[0], candidate.genDecVtx[1], candidate.genDecVtx[2], + candidate.isReco, candidate.isSignal, candidate.isRecoMCCollision, candidate.isSurvEvSel, candidate.isTwoBodyDecay, candidate.statusCode, + candidate.heTruth.pt, candidate.heTruth.isPhysicalPrimary, candidate.decayPionTruth.isPhysicalPrimary, + std::get<0>(pion.info), std::get<1>(pion.info), std::get<2>(pion.info), + pion.truth.pt, pion.truth.eta, pion.truth.phi, pion.truth.isPhysicalPrimary, pion.truth.process, + sameMCCollision, matchesHypRecoMCCollision, matchesPairRecoMCCollision, truthSelfCorrelation, isTruePrimaryPiHyperPair, hadIsTruePion}; + } + + void fillMCQA(const HadHyperCandidate& candidate, const HadHyperHadron& pion, const HadHyperEvent& hyperEvent, float kstar, bool mixed) + { + if (!candidate.isReco || !candidate.isSignal || std::abs(candidate.hyperTruth.pdgCode) != HyperTritonPDG || std::abs(pion.truth.pdgCode) != static_cast(PDG_t::kPiPlus)) { + return; + } + + const float kstarMC = computePairKstar(pion.truth.momentum, candidate.genMomentum); + const auto fill = [&](const auto& prefix) { + if (std::isfinite(kstarMC)) { + hadHyperRegistry.fill(prefix + HIST("hKstarRecVsGenHyperReco"), kstarMC, kstar); + hadHyperRegistry.fill(prefix + HIST("hKstarResolutionHyperReco"), kstar, kstar - kstarMC); + } + + hadHyperRegistry.fill(prefix + HIST("hPrimaryPionVsKstarDen"), kstar); + hadHyperRegistry.fill(prefix + HIST("hPrimaryPionVsCentDen"), hyperEvent.centrality); + if (pion.truth.isPhysicalPrimary) { + hadHyperRegistry.fill(prefix + HIST("hPrimaryPionVsKstarNum"), kstar); + hadHyperRegistry.fill(prefix + HIST("hPrimaryPionVsCentNum"), hyperEvent.centrality); + } + }; + if (mixed) { + fill(HIST("MC/ME/")); + } else { + fill(HIST("MC/SE/")); + } + } + + template + bool fillPair(const HadHyperCandidate& candidate, const HadHyperHadron& pion, const HadHyperEvent& hyperEvent, const HadHyperEvent& pionEvent, bool mixed, int mixingDepth) + { + const float kstar = computePairKstar(pion.momentum, candidate.momentum); + if (isRecoSelfCorrelation(candidate, pion, hyperEvent, pionEvent)) { + const int reason = pion.sourceId == candidate.heTrackId ? 0 : 1; + hadHyperRegistry.fill(HIST("hSelfPairs"), reason, kstar); + if (!mixed) { + hadHyperRegistry.fill(HIST("hSameEventSelfPairs"), reason, kstar); + } + return false; + } + if constexpr (isMC) { + if (isTruthSelfCorrelation(candidate, pion, hyperEvent, pionEvent)) { + hadHyperRegistry.fill(HIST("hSelfPairs"), 2, kstar); + if (!mixed) { + hadHyperRegistry.fill(HIST("hSameEventSelfPairs"), 2, kstar); + } + return false; + } + } + + if (!std::isfinite(kstar)) { + return false; + } + // CPR is an offline flag only: do not remove rows or change raw-pair/MC QA counts. + const uint8_t closePairRejected = flagClosePair(candidate, pion, mixed); + + if (mixed) { + hadHyperRegistry.fill(HIST("hME"), kstar); + } else { + hadHyperRegistry.fill(HIST("hSE"), kstar); + } + if constexpr (isMC) { + fillMCQA(candidate, pion, hyperEvent, kstar, mixed); + } + + if (!output.settingFillTable || (hadHyper.maxOutputKstar.value > 0.f && kstar >= hadHyper.maxOutputKstar.value)) { + return true; + } + + const auto dataInfo = makeDataInfo(candidate, pion, hyperEvent, mixed, closePairRejected, mixingDepth); + if constexpr (isMC) { + std::apply([this](const auto&... columns) { mOutputMCTable(columns...); }, std::tuple_cat(dataInfo, makeMCInfo(candidate, pion, hyperEvent, pionEvent))); + } else { + std::apply([this](const auto&... columns) { mOutputDataTable(columns...); }, dataInfo); + } + return true; + } + + template + bool hasValidHyperDaughterIndices(const Ttracks& tracks, int64_t heTrackId, int64_t pionTrackId) const + { + const auto first = static_cast(tracks.offset()); + const auto last = first + static_cast(tracks.size()); + return heTrackId >= first && heTrackId < last && pionTrackId >= first && pionTrackId < last; + } + + template + void collectPions(HadHyperEvent& event, const Ttracks& eventTracks) + { + const std::array primaryVertex{std::get<5>(event.info), std::get<6>(event.info), std::get<7>(event.info)}; + for (const auto& track : eventTracks) { + if (!selectPionTrack(track) || !selectionPIDPion(track)) { + continue; + } + event.hadrons.push_back(makePion(track, primaryVertex)); + hadHyperRegistry.fill(HIST("hPionTPC"), track.tpcInnerParam(), track.tpcSignal()); + } + } + + template + void collectHyperCandidates(HadHyperEvent& event, const Tcandidates& eventCandidates, const Ttracks& tracks) + { + const std::array primaryVertex{std::get<5>(event.info), std::get<6>(event.info), std::get<7>(event.info)}; + for (const auto& candidate : eventCandidates) { + if constexpr (isMC) { + if (!candidate.isReco()) { + continue; + } + } + if (!hasValidHyperDaughterIndices(tracks, candidate.heTrackId(), candidate.piTrackId())) { + if constexpr (isMC) { + continue; + } else { + LOG(fatal) << "Hypertriton daughter indices must reference the input Tracks table"; + } + } + if (!selectHyperCandidate(candidate)) { + continue; + } + event.candidates.push_back(makeHyperCandidate(candidate, tracks, primaryVertex)); + } + } + + template + HadHyperEvent buildEvent(const Tcollision& collision, const Ttracks& tracks, const Tcandidates& candidates) + { + HadHyperEvent event; + event.info = makeEventInfo(collision); + event.centrality = collision.centFT0C(); + if constexpr (isMC) { + event.hasMCCollision = collision.has_mcCollision(); + event.mcCollisionId = collision.has_mcCollision() ? collision.mcCollisionId() : -1; + } + + if constexpr (isMC) { + const auto eventTracks = tracks.sliceBy(mPerColMC, collision.globalIndex()); + collectPions(event, eventTracks); + const auto eventCandidates = candidates.sliceBy(hypPerColMC, collision.globalIndex()); + collectHyperCandidates(event, eventCandidates, tracks); + } else { + const auto eventTracks = tracks.sliceBy(mPerCol, collision.globalIndex()); + collectPions(event, eventTracks); + const auto eventCandidates = candidates.sliceBy(hypPerCol, collision.globalIndex()); + collectHyperCandidates(event, eventCandidates, tracks); + } + return event; + } + + template + void fillSameEventPairs(const HadHyperEvent& event) + { + for (const auto& candidate : event.candidates) { + int acceptedPairs = 0; + for (const auto& pion : event.hadrons) { + if (fillPair(candidate, pion, event, event, false, 0)) { + ++acceptedPairs; + } + } + hadHyperRegistry.fill(HIST("hCandidatePairMultiplicitySE"), acceptedPairs, event.centrality); + } + } + + template + void fillMixedEventPairs(HadHyperEvent& currentEvent, std::deque& pool) + { + const int depth = static_cast(pool.size()); + hadHyperRegistry.fill(HIST("hMixingDepth"), depth); + for (std::size_t partnerIndex = 0; partnerIndex < pool.size(); ++partnerIndex) { + auto& partner = pool[partnerIndex]; + const float currentPosZ = std::get<0>(currentEvent.info); + const float partnerPosZ = std::get<0>(partner.info); + hadHyperRegistry.fill(HIST("hMixEventDeltaPosZVsCent"), currentEvent.centrality, currentPosZ - partnerPosZ); + hadHyperRegistry.fill(HIST("hMixEventDeltaCentFT0CVsCent"), currentEvent.centrality, currentEvent.centrality - partner.centrality); + hadHyperRegistry.fill(HIST("hMixEventDeltaPosZVsCent"), partner.centrality, partnerPosZ - currentPosZ); + hadHyperRegistry.fill(HIST("hMixEventDeltaCentFT0CVsCent"), partner.centrality, partner.centrality - currentEvent.centrality); + for (std::size_t candidateIndex = 0; candidateIndex < currentEvent.candidates.size(); ++candidateIndex) { + auto& candidate = currentEvent.candidates[candidateIndex]; + for (const auto& pion : partner.hadrons) { + if (fillPair(candidate, pion, currentEvent, partner, true, depth)) { + ++candidate.mixedPairCount; + } + } + } + for (std::size_t candidateIndex = 0; candidateIndex < partner.candidates.size(); ++candidateIndex) { + auto& candidate = partner.candidates[candidateIndex]; + for (const auto& pion : currentEvent.hadrons) { + if (fillPair(candidate, pion, partner, currentEvent, true, depth)) { + ++candidate.mixedPairCount; + } + } + } + } + } + + void flushMixedEventMultiplicity(const HadHyperEvent& event) + { + // One entry per candidate, including zero partners, after all its actual ME pairings. + for (const auto& candidate : event.candidates) { + hadHyperRegistry.fill(HIST("hCandidatePairMultiplicityME"), candidate.mixedPairCount, event.centrality); + } + } + + void flushMixingPools(std::unordered_map>& pools) + { + for (const auto& [bin, pool] : pools) { + for (const auto& event : pool) { + flushMixedEventMultiplicity(event); + } + } + pools.clear(); + } + + void storeEventInPool(std::deque& pool, HadHyperEvent&& event) + { + const int requestedMixingDepth = eventMixing.settingNoMixedEvents.value; + if (requestedMixingDepth <= 0) { + return; + } + const auto mixingDepth = static_cast(requestedMixingDepth); + if (pool.size() >= mixingDepth) { + flushMixedEventMultiplicity(pool.front()); + pool.pop_front(); + } + pool.push_back(std::move(event)); + } + + void endOfStream(o2::framework::EndOfStreamContext&) + { + flushMixingPools(mMixingPools); + flushMixingPools(mMCMixingPools); + } + + void processHyper(const HadHyperCollisionsFull& collisions, const TrackCandidates& tracks, + const HyperCandidates& candidates, const aod::BCsWithTimestamps& bcs) + { + processPairs(collisions, tracks, candidates, bcs, mMixingPools, mMixingRunNumber); + } + PROCESS_SWITCH(PiHypertritonFemto, processHyper, "Process pion-hypertriton same-event and mixed-event pairs", false); + + void processMCHyper(const HadHyperCollisionsFullMC& collisions, const TrackCandidatesMC& tracks, + const HyperCandidatesMC& candidates, const aod::McParticles&, const aod::BCsWithTimestamps& bcs) + { + processPairs(collisions, tracks, candidates, bcs, mMCMixingPools, mMCMixingRunNumber); + } + PROCESS_SWITCH(PiHypertritonFemto, processMCHyper, "Process MC pion-hypertriton same-event and mixed-event pairs", false); +}; + +void PiHypertritonFemto::init(o2::framework::InitContext&) +{ + *mEndOfStreamCallback = [this](EndOfStreamContext& context) { endOfStream(context); }; + hadHyperRegistry.addClone("MC/SE/", "MC/ME/"); + const AxisSpec pairMultiplicityAxis{axisPairMultiplicity, "pairs"}; + hadHyperRegistry.add("hCandidatePairMultiplicitySE", "Accepted same-event pion pairs per hypertriton candidate;pairs;centrality", + HistType::kTH2F, {pairMultiplicityAxis, {10, 0., 100.}}); + hadHyperRegistry.add("hCandidatePairMultiplicityME", "Total mixed-event pion pairs per hypertriton candidate over its pool lifetime;pairs;centrality", + HistType::kTH2F, {pairMultiplicityAxis, {10, 0., 100.}}); + + const std::array eventNames{"SE", "ME"}; + const std::array daughterNames{"He3", "DecayPi"}; + for (std::size_t eventIndex = 0; eventIndex < eventNames.size(); ++eventIndex) { + for (std::size_t daughterIndex = 0; daughterIndex < daughterNames.size(); ++daughterIndex) { + const std::string prefix = "CPR/" + eventNames[eventIndex] + "/" + daughterNames[daughterIndex]; + auto& qa = mClosePairQA[eventIndex][daughterIndex]; + qa.before = hadHyperRegistry.add((prefix + "/hBefore").c_str(), "All pairs;#Delta#eta;#LT#Delta#varphi*#GT", + HistType::kTH2F, {{200, -0.1, 0.1}, {200, -0.1, 0.1}}); + qa.after = hadHyperRegistry.add((prefix + "/hAfter").c_str(), "Pairs with offline CPR flags zero;#Delta#eta;#LT#Delta#varphi*#GT", + HistType::kTH2F, {{200, -0.1, 0.1}, {200, -0.1, 0.1}}); + qa.validRadii = hadHyperRegistry.add((prefix + "/hValidRadii").c_str(), "Common valid TPC radii (0=undefined CPR);radii;Comparisons", + HistType::kTH1F, {{10, -0.5, 9.5}}); + qa.distanceBefore = hadHyperRegistry.add((prefix + "/hDistanceBefore").c_str(), "All pairs;Average TPC separation (cm);Comparisons", + HistType::kTH1F, {{500, 0., 100.}}); + qa.distanceAfter = hadHyperRegistry.add((prefix + "/hDistanceAfter").c_str(), "Pairs with offline CPR flags zero;Average TPC separation (cm);Comparisons", + HistType::kTH1F, {{500, 0., 100.}}); + qa.validPositionRadii = hadHyperRegistry.add((prefix + "/hValidPositionRadii").c_str(), "Common extrapolated TPC radii (0=undefined distance);radii;Comparisons", + HistType::kTH1F, {{10, -0.5, 9.5}}); + qa.decisions = hadHyperRegistry.add((prefix + "/hDecisions").c_str(), "Valid angular and distance comparisons;Angular flag;Distance flag", + HistType::kTH2F, {{2, -0.5, 1.5}, {2, -0.5, 1.5}}); + } + } + + const bool processHyperPairs = doprocessHyper || doprocessMCHyper; + if (processHyperPairs && hadHyper.maxOutputKstar.value == 0.f) { + LOG(fatal) << "Pion-hypertriton mode requires a nonzero output k* range"; + } + if (processHyperPairs && hadHyper.enableMixing.value && eventMixing.settingNoMixedEvents.value <= 0) { + LOG(fatal) << "Pion-hypertriton mixed-event mode requires a positive mixing depth"; + } + if (processHyperPairs && hypertriton.settingHypMassMin.value >= hypertriton.settingHypMassMax.value) { + LOG(fatal) << "Pion-hypertriton mode requires settingHypMassMin < settingHypMassMax"; + } + if (!std::isfinite(CPR.settingClosePairDeltaEtaMax.value) || !std::isfinite(CPR.settingClosePairDeltaPhiMax.value) || + CPR.settingClosePairDeltaEtaMax.value <= 0.f || CPR.settingClosePairDeltaPhiMax.value <= 0.f) { + LOG(fatal) << "Close-pair rejection requires positive delta-eta and delta-phi limits"; + } + if (!std::isfinite(CPR.settingClosePairDistanceMax.value) || CPR.settingClosePairDistanceMax.value <= 0.f) { + LOG(fatal) << "Close-pair tagging requires a finite positive distance threshold"; + } + mZorroSummary.setObject(mZorro.getZorroSummary()); + mRunNumber = -1; + mCcdb->setURL(ccdb.settingCcdburl); + mCcdb->setCaching(true); + mCcdb->setLocalObjectValidityChecking(); + mCcdb->setFatalWhenNull(false); + + const std::array eventsLabels = { + "All", + "kNoITSROFrameBorder (data, optional)", + "kIsTriggerTVX", + "kNoTimeFrameBorder", + "z_{vtx}", + "kNoSameBunchPileup (optional)", + "kIsGoodZvtxFT0vsPV (optional)", + "sel8 (MC, optional)", + "Zorro (data, optional)", + "MC collision label (MC, optional)", + "Mixing bin (if enabled)"}; + for (std::size_t i = 0; i < eventsLabels.size(); i++) { + mQaRegistry.get(HIST("hEvents"))->GetXaxis()->SetBinLabel(i + 1, eventsLabels[i].c_str()); + } +} + +template +void PiHypertritonFemto::initCCDB(const aod::BCsWithTimestamps::iterator& bc) +{ + if (mRunNumber == bc.runNumber()) { + return; + } + auto* magneticField = mCcdb->getForTimeStamp(ccdb.settingGrpmagPath.value, bc.timestamp()); + if (!magneticField) { + LOG(fatal) << "Cannot load GRPMagField for daughter close-pair tagging at timestamp " << bc.timestamp(); + return; + } + constexpr float KilogaussToTesla = 0.1f; + mMagneticFieldTesla = KilogaussToTesla * magneticField->getNominalL3Field(); + if constexpr (!isMC) { + if (zorro.settingSkimmedProcessing.value) { + mZorro.initCCDB(mCcdb.service, bc.runNumber(), bc.timestamp(), zorro.settingTriggerMask.value.c_str()); + mZorro.populateHistRegistry(mQaRegistry, bc.runNumber()); + } + } + mRunNumber = bc.runNumber(); +} + +template +bool PiHypertritonFemto::selectCollision(const Tcollision& collision, const aod::BCsWithTimestamps&) +{ + auto bc = collision.template bc_as(); + initCCDB(bc); + mQaRegistry.fill(HIST("hEvents"), 0); + + // Match hyperRecoTask: the ITS readout-frame border cut applies only to data. + if constexpr (!isMC) { + if (!eventSelection.disableITSROFCut.value && !collision.selection_bit(aod::evsel::kNoITSROFrameBorder)) { + return false; + } + } + // Disabled or inapplicable cuts pass through, keeping the cut-flow cumulative. + mQaRegistry.fill(HIST("hEvents"), 1); + if (!collision.selection_bit(aod::evsel::kIsTriggerTVX)) { + return false; + } + mQaRegistry.fill(HIST("hEvents"), 2); + if (!collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) { + return false; + } + mQaRegistry.fill(HIST("hEvents"), 3); + if (std::abs(collision.posZ()) > eventMixing.settingCutVertex.value) { + return false; + } + mQaRegistry.fill(HIST("hEvents"), 4); + if (eventSelection.cfgEvSelkNoSameBunchPileup.value && !collision.selection_bit(aod::evsel::kNoSameBunchPileup)) { + return false; + } + mQaRegistry.fill(HIST("hEvents"), 5); + if (eventSelection.cfgEvSelkIsGoodZvtxFT0vsPV.value && !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) { + return false; + } + mQaRegistry.fill(HIST("hEvents"), 6); + if constexpr (isMC) { + if (mc.settingRequireSel8.value && !collision.sel8()) { + return false; + } + } + + mQaRegistry.fill(HIST("hEvents"), 7); + if constexpr (!isMC) { + if (zorro.settingSkimmedProcessing.value) { + if (!mZorro.isSelected(bc.globalBC())) { + return false; + } + } + } + mQaRegistry.fill(HIST("hEvents"), 8); + mQaRegistry.fill(HIST("hNcontributor"), collision.numContrib()); + mQaRegistry.fill(HIST("hVtxZ"), collision.posZ()); + return true; +} + +template +bool PiHypertritonFemto::selectPionTrack(const Ttrack& candidate) const +{ + const float absPt = std::abs(candidate.pt()); + if (std::abs(candidate.eta()) > pionTrack.settingCutEta.value || absPt < pionTrack.settingPtMin.value || absPt > pionTrack.settingPtMax.value || absPt <= 0.f) { + return false; + } + if (candidate.itsNClsInnerBarrel() < pionTrack.settingITSInnerBarrelMin.value || + candidate.itsNCls() < pionTrack.settingITSNClsMin.value || + candidate.tpcNClsFound() < pionTrack.settingTPCNClsFoundMin.value || + candidate.tpcNClsCrossedRows() < pionTrack.settingTPCCrossedRowsMin.value) { + return false; + } + const float pionDCAxyMax = pionTrack.settingDCAxyOffset.value + pionTrack.settingDCAxyPtCoeff.value / absPt; + const float pionDCAzMax = pionTrack.settingDCAzOffset.value + pionTrack.settingDCAzPtCoeff.value / absPt; + return !(std::abs(candidate.dcaXY()) > pionDCAxyMax || std::abs(candidate.dcaZ()) > pionDCAzMax); +} + +template +bool PiHypertritonFemto::selectionPIDPion(const Ttrack& candidate) +{ + const float tpcNSigmaPi = candidate.tpcNSigmaPi(); + const float absP = std::abs(candidate.p()); + mQaRegistry.fill(HIST("hPionTPCNSigmaPreselection"), candidate.sign() * candidate.tpcInnerParam(), tpcNSigmaPi); + if (absP <= pionPid.settingMomCombMin.value) { + if (std::abs(tpcNSigmaPi) > pionPid.settingTPCNsigMax.value) { + return false; + } + mQaRegistry.fill(HIST("hPionTPCNSigma"), candidate.sign() * candidate.pt(), tpcNSigmaPi); + mQaRegistry.fill(HIST("hPionTPC"), candidate.sign() * candidate.tpcInnerParam(), candidate.tpcSignal()); + return true; + } + if (!candidate.hasTOF()) { + return false; + } + const float tofNSigmaPi = candidate.tofNSigmaPi(); + const float combNsigma = std::sqrt(tofNSigmaPi * tofNSigmaPi + tpcNSigmaPi * tpcNSigmaPi); + mQaRegistry.fill(HIST("hPionTOFNSigmaPreselection"), candidate.sign() * candidate.pt(), tofNSigmaPi); + mQaRegistry.fill(HIST("hPionCombNSigmaPreselection"), candidate.sign() * candidate.pt(), combNsigma); + if (combNsigma > pionPid.settingCombNsigMax.value) { + return false; + } + if (pionPid.settingReqSingleNsig.value && (std::abs(tpcNSigmaPi) > pionPid.settingCombNsigMax.value || std::abs(tofNSigmaPi) > pionPid.settingCombNsigMax.value)) { + return false; + } + mQaRegistry.fill(HIST("hPionTPCNSigma"), candidate.sign() * candidate.pt(), tpcNSigmaPi); + mQaRegistry.fill(HIST("hPionTOFNSigma"), candidate.sign() * candidate.pt(), tofNSigmaPi); + mQaRegistry.fill(HIST("hPionCombNSigma"), candidate.sign() * candidate.pt(), combNsigma); + mQaRegistry.fill(HIST("hPionTPC"), candidate.sign() * candidate.tpcInnerParam(), candidate.tpcSignal()); + return true; +} + +template +bool PiHypertritonFemto::selectHyperCandidate(const Tcandidate& candidate) +{ + mQaRegistry.fill(HIST("hHypHe3TPCMomPreselection"), candidate.tpcMomHe()); + const float mass = computeHyperCandidateMass(candidate); + if (!std::isfinite(mass) || mass < hypertriton.settingHypMassMin.value || mass > hypertriton.settingHypMassMax.value) { + return false; + } + mQaRegistry.fill(HIST("hHypHe3TPCMom"), candidate.tpcMomHe()); + mQaRegistry.fill(HIST("hHypHe3TPCNSigma"), candidate.ptHe3(), candidate.nSigmaHe()); + return true; +} + +template +void PiHypertritonFemto::processPairs(const Tcollisions& collisions, const Ttracks& tracks, const Tcandidates& candidates, + const aod::BCsWithTimestamps& bcs, + std::unordered_map>& mixingPools, + int& mixingRunNumber) +{ + // Each invocation receives a new input-table scope; row indices can repeat in later frames. + const uint64_t sourceFrameId = ++mSourceFrameId; + const BinningType configuredBinningPolicy{{axisVertex, axisCentrality}, true}; + for (const auto& collision : collisions) { + if (!selectCollision(collision, bcs)) { + continue; + } + if constexpr (isMC) { + if (mc.settingRequireRecoMCCollisionMatch.value && !collision.has_mcCollision()) { + continue; + } + } + mQaRegistry.fill(HIST("hEvents"), 9); + hadHyperRegistry.fill(HIST("hPoolFlow"), 0); + int poolBin = -1; + if (hadHyper.enableMixing.value) { + poolBin = configuredBinningPolicy.getBin(std::make_tuple(collision.posZ(), collision.centFT0C())); + if (poolBin < 0) { + continue; + } + hadHyperRegistry.fill(HIST("hPoolFlow"), 1); + } + mQaRegistry.fill(HIST("hEvents"), 10); + auto event = buildEvent(collision, tracks, candidates); + event.sourceFrameId = sourceFrameId; + fillSameEventPairs(event); + if (!hadHyper.enableMixing.value) { + continue; + } + const auto bc = collision.template bc_as(); + if (mixingRunNumber != bc.runNumber()) { + flushMixingPools(mixingPools); + mixingRunNumber = bc.runNumber(); + } + auto& pool = mixingPools[poolBin]; + fillMixedEventPairs(event, pool); + storeEventInPool(pool, std::move(event)); + } +} + +WorkflowSpec defineDataProcessing(const ConfigContext& cfgc) +{ + auto finishMixing = std::make_shared>(); + auto task = adaptAnalysisTask(cfgc, finishMixing); + auto initializeTask = task.algorithm.onInit; + task.algorithm.onInit = [initializeTask, finishMixing](InitContext& context) { + // CallbackService preserves registration order. Register before adaptAnalysisTask's + // initializer so the pools are flushed BEFORE its histogram snapshot/cleanup callback. + context.services().get().set( + [finishMixing](EndOfStreamContext& eosContext) { + if (*finishMixing) { + (*finishMixing)(eosContext); + } + }); + return initializeTask(context); + }; + return WorkflowSpec{std::move(task)}; +} diff --git a/PWGLF/TableProducer/Nuspex/trackedHypertritonRecoTask.cxx b/PWGLF/TableProducer/Nuspex/trackedHypertritonRecoTask.cxx index d63be90ba45..4297ebc4eba 100644 --- a/PWGLF/TableProducer/Nuspex/trackedHypertritonRecoTask.cxx +++ b/PWGLF/TableProducer/Nuspex/trackedHypertritonRecoTask.cxx @@ -1033,8 +1033,7 @@ struct TrackedHypertritonRecoTask { continue; } buildTwoBody(heTrack, piTrack, collision, trackedV0.itsClsSize(), [&](auto... values) { - mcHypCands(collision.globalIndex(), heTrack.globalIndex(), piTrack.globalIndex(), - values..., + mcHypCands(values..., mcInfo.genPt, mcInfo.genPhi, mcInfo.genEta, mcInfo.genPtHe3, mcInfo.genDecayVertex[0], mcInfo.genDecayVertex[1], mcInfo.genDecayVertex[2], true, mcInfo.fakeHeITSLayerMap, mcInfo.isSignal, @@ -1170,8 +1169,7 @@ struct TrackedHypertritonRecoTask { primaryVertexZ = collision.posZ(); } } - mcHypCands(-1, -1, -1, - centralityFT0A, centralityFT0C, centralityFT0M, + mcHypCands(centralityFT0A, centralityFT0C, centralityFT0M, trackOccupancyInTimeRange, ft0cOccupancyInTimeRange, primaryVertexX, primaryVertexY, primaryVertexZ, runNumber, mother.pdgCode() > 0,