diff --git a/PWGLF/DataModel/LFDoubleOmegaTables.h b/PWGLF/DataModel/LFDoubleOmegaTables.h index 16c03e444f7..83cee30b84a 100644 --- a/PWGLF/DataModel/LFDoubleOmegaTables.h +++ b/PWGLF/DataModel/LFDoubleOmegaTables.h @@ -62,6 +62,8 @@ DECLARE_SOA_COLUMN(GenDecayLength, genDecayLength, float); DECLARE_SOA_COLUMN(PdgDoubleOmega, pdgDoubleOmega, int); DECLARE_SOA_COLUMN(GenMotherId, genMotherId, int64_t); // MC particle index within the input dataframe. DECLARE_SOA_COLUMN(IsReco, isReco, bool); + +DECLARE_SOA_COLUMN(MixingType, mixingType, uint8_t); // 1: Omega, 2: Lambda, 3: kaon from another collision; 4: Omega, Lambda and kaon from three different collisions. } // namespace DoubleOmegaTables #define DOUBLE_OMEGA_RECO_COLUMNS \ @@ -107,8 +109,116 @@ DECLARE_SOA_TABLE(DoubleOmegaTableMC, "AOD", "DBLOMEGAMCTABLE", DoubleOmegaTables::GenMotherId, DoubleOmegaTables::IsReco); +DECLARE_SOA_TABLE(DoubleOmegaMixedTable, "AOD", "DBLOMEGAMIXTABLE", + DOUBLE_OMEGA_RECO_COLUMNS, + DoubleOmegaTables::MixingType); + #undef DOUBLE_OMEGA_RECO_COLUMNS +// Derived data: Omega, Lambda and kaon candidates of the collisions with at least one Omega, one +// Lambda and one kaon candidate not sharing tracks with it. The tables only index each other, so +// they are self-contained and dataframes can be merged. +namespace domega +{ +DECLARE_SOA_COLUMN(MultFT0M, multFT0M, float); //! FT0M multiplicity of the collision +DECLARE_SOA_COLUMN(Bz, bz, float); //! Nominal magnetic field (kG) +} // namespace domega + +DECLARE_SOA_TABLE(DOmegaColls, "AOD", "DOMEGACOLL", //! Collisions with double-Omega building blocks + o2::soa::Index<>, + collision::PosX, collision::PosY, collision::PosZ, + domega::MultFT0M, + bc::RunNumber, + domega::Bz); + +namespace domega +{ +DECLARE_SOA_INDEX_COLUMN_FULL(DOmegaColl, dOmegaColl, int, DOmegaColls, ""); //! Collision index + +// Track parametrisation covariance, in the order expected by getTrackParCov. +DECLARE_SOA_COLUMN(CYY, cYY, float); +DECLARE_SOA_COLUMN(CZY, cZY, float); +DECLARE_SOA_COLUMN(CZZ, cZZ, float); +DECLARE_SOA_COLUMN(CSnpY, cSnpY, float); +DECLARE_SOA_COLUMN(CSnpZ, cSnpZ, float); +DECLARE_SOA_COLUMN(CSnpSnp, cSnpSnp, float); +DECLARE_SOA_COLUMN(CTglY, cTglY, float); +DECLARE_SOA_COLUMN(CTglZ, cTglZ, float); +DECLARE_SOA_COLUMN(CTglSnp, cTglSnp, float); +DECLARE_SOA_COLUMN(CTglTgl, cTglTgl, float); +DECLARE_SOA_COLUMN(C1PtY, c1PtY, float); +DECLARE_SOA_COLUMN(C1PtZ, c1PtZ, float); +DECLARE_SOA_COLUMN(C1PtSnp, c1PtSnp, float); +DECLARE_SOA_COLUMN(C1PtTgl, c1PtTgl, float); +DECLARE_SOA_COLUMN(C1Pt21Pt2, c1Pt21Pt2, float); + +DECLARE_SOA_COLUMN(TPCNSigmaKa, tpcNSigmaKa, float); //! -999 when TPC is unavailable +DECLARE_SOA_COLUMN(TOFNSigmaKa, tofNSigmaKa, float); //! -999 when TOF is unavailable +DECLARE_SOA_COLUMN(DetectorMap, detectorMap, uint8_t); //! ITS, TPC, TOF bits +DECLARE_SOA_COLUMN(DcaXYToPV, dcaXYToPV, float); //! Transverse DCA to the primary vertex of its collision +DECLARE_SOA_COLUMN(IsKaon, isKaon, bool); //! Selected as a direct-kaon candidate +} // namespace domega + +#define DOMEGA_TRACK_PAR_COV_COLUMNS \ + track::X, track::Alpha, track::Y, track::Z, track::Snp, track::Tgl, track::Signed1Pt, \ + domega::CYY, domega::CZY, domega::CZZ, domega::CSnpY, domega::CSnpZ, domega::CSnpSnp, \ + domega::CTglY, domega::CTglZ, domega::CTglSnp, domega::CTglTgl, \ + domega::C1PtY, domega::C1PtZ, domega::C1PtSnp, domega::C1PtTgl, domega::C1Pt21Pt2 + +DECLARE_SOA_TABLE(DOmegaTracks, "AOD", "DOMEGATRACK", //! Direct-kaon candidates and Omega/Lambda daughters, at the innermost update + o2::soa::Index<>, + domega::DOmegaCollId, + DOMEGA_TRACK_PAR_COV_COLUMNS, + domega::TPCNSigmaKa, + domega::TOFNSigmaKa, + domega::DetectorMap, + domega::DcaXYToPV, + domega::IsKaon, + track::Sign); + +namespace domega +{ +DECLARE_SOA_INDEX_COLUMN_FULL(PosTrack, posTrack, int, DOmegaTracks, "_Pos"); //! Positive daughter +DECLARE_SOA_INDEX_COLUMN_FULL(NegTrack, negTrack, int, DOmegaTracks, "_Neg"); //! Negative daughter +DECLARE_SOA_INDEX_COLUMN_FULL(Bachelor, bachelor, int, DOmegaTracks, "_Bach"); //! Bachelor kaon of the Omega + +DECLARE_SOA_COLUMN(DecayX, decayX, float); +DECLARE_SOA_COLUMN(DecayY, decayY, float); +DECLARE_SOA_COLUMN(DecayZ, decayZ, float); +DECLARE_SOA_COLUMN(Px, px, float); +DECLARE_SOA_COLUMN(Py, py, float); +DECLARE_SOA_COLUMN(Pz, pz, float); +DECLARE_SOA_COLUMN(MassOmega, massOmega, float); +DECLARE_SOA_COLUMN(MassXi, massXi, float); +DECLARE_SOA_COLUMN(MassLambda, massLambda, float); +DECLARE_SOA_COLUMN(DcaDaughters, dcaDaughters, float); //! DCA between the daughters +DECLARE_SOA_COLUMN(Sign, sign, int8_t); //! Charge of the Omega (Omega table) or of the Omega the Lambda can be paired with (Lambda table) +} // namespace domega + +DECLARE_SOA_TABLE(DOmegaOmegas, "AOD", "DOMEGAOMEGA", //! Omega candidates; track parameters at the Omega decay vertex + o2::soa::Index<>, + domega::DOmegaCollId, + domega::PosTrackId, domega::NegTrackId, domega::BachelorId, + DOMEGA_TRACK_PAR_COV_COLUMNS, + domega::DecayX, domega::DecayY, domega::DecayZ, + domega::Px, domega::Py, domega::Pz, + domega::MassOmega, domega::MassXi, + domega::DcaDaughters, + domega::Sign); + +DECLARE_SOA_TABLE(DOmegaLambdas, "AOD", "DOMEGALAMBDA", //! Lambda candidates; neutral track parameters at the Lambda decay vertex + o2::soa::Index<>, + domega::DOmegaCollId, + domega::PosTrackId, domega::NegTrackId, + DOMEGA_TRACK_PAR_COV_COLUMNS, + domega::DecayX, domega::DecayY, domega::DecayZ, + domega::Px, domega::Py, domega::Pz, + domega::MassLambda, + domega::DcaDaughters, + domega::Sign); + +#undef DOMEGA_TRACK_PAR_COV_COLUMNS + } // namespace o2::aod #endif // PWGLF_DATAMODEL_LFDOUBLEOMEGATABLES_H_ diff --git a/PWGLF/TableProducer/Strangeness/CMakeLists.txt b/PWGLF/TableProducer/Strangeness/CMakeLists.txt index a61a01e6a56..40ff199253f 100644 --- a/PWGLF/TableProducer/Strangeness/CMakeLists.txt +++ b/PWGLF/TableProducer/Strangeness/CMakeLists.txt @@ -61,6 +61,11 @@ o2physics_add_dpl_workflow(double-omega-tree-creator PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::EventFilteringUtils COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(double-omega-builder + SOURCES doubleOmegaBuilder.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::EventFilteringUtils + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(lambdakzerofinder SOURCES lambdakzerofinder.cxx PUBLIC_LINK_LIBRARIES O2::DCAFitter O2Physics::AnalysisCore diff --git a/PWGLF/TableProducer/Strangeness/doubleOmegaBuilder.cxx b/PWGLF/TableProducer/Strangeness/doubleOmegaBuilder.cxx new file mode 100644 index 00000000000..459997de976 --- /dev/null +++ b/PWGLF/TableProducer/Strangeness/doubleOmegaBuilder.cxx @@ -0,0 +1,321 @@ +// 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 doubleOmegaBuilder.cxx +/// \brief Builds the Omega, Lambda and direct-kaon candidates of the double-Omega analysis and stores +/// them in derived tables, consumed on the fly (or from saved derived data) by doubleOmegaTreeCreator +/// for the same-event candidates and the event mixing. + +#include "PWGLF/DataModel/LFDoubleOmegaTables.h" +#include "PWGLF/Utils/doubleOmegaBuilderHelper.h" + +#include "Common/CCDB/EventSelectionParams.h" +#include "Common/Core/Zorro.h" +#include "Common/Core/ZorroSummary.h" +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/PIDResponseTOF.h" +#include "Common/DataModel/PIDResponseTPC.h" + +#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::pwglf::doubleomega; + +using CollisionsTable = soa::Join; +using TracksFull = soa::Join; + +struct doubleOmegaBuilder { + Produces dOmegaColls; + Produces dOmegaTracks; + Produces dOmegaOmegas; + Produces dOmegaLambdas; + Service ccdb{}; + + enum CollisionStep : uint8_t { + kAllCollisions, + kSelectedCollisions, + kWithOmega, + kStoredCollisions, + kNCollisionSteps + }; + + Preslice tracksPerCollision = aod::track::collisionId; + Preslice cascadesPerCollision = aod::cascade::collisionId; + Preslice v0sPerCollision = aod::v0::collisionId; + + DoubleOmegaBuilderHelper helper; + int mRunNumber = 0; + float mBz = 0.f; + + Zorro zorro; + OutputObj zorroSummary{"zorroSummary"}; + + Configurable cfgSkimmedProcessing{"cfgSkimmedProcessing", false, "Skimmed dataset processing"}; + Configurable cfgApplyEventSelection{"cfgApplyEventSelection", true, "Apply the standard collision event selection"}; + Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable cfgMaterialCorrection{"cfgMaterialCorrection", static_cast(o2::base::Propagator::MatCorrType::USEMatCorrLUT), "Material correction for the V0/cascade fits"}; + + ConfigurableAxis zVtxAxis{"zVtxBins", {100, -20.f, 20.f}, "Binning for the vertex z in cm"}; + ConfigurableAxis massOmegaAxis{"massOmegaAxis", {400, o2::constants::physics::MassOmegaMinus - 0.05f, o2::constants::physics::MassOmegaMinus + 0.05f}, "binning for the Omega invariant-mass"}; + ConfigurableAxis massXiAxis{"massXiAxis", {400, o2::constants::physics::MassXiMinus - 0.05f, o2::constants::physics::MassXiMinus + 0.05f}, "binning for the Xi invariant-mass"}; + ConfigurableAxis massLambdaAxis{"massLambdaAxis", {400, o2::constants::physics::MassLambda0 - 0.05f, o2::constants::physics::MassLambda0 + 0.05f}, "binning for the Lambda invariant-mass"}; + ConfigurableAxis momAxis{"momAxisFine", {5.e2, 0.f, 5.f}, "momentum axis binning"}; + + // Production selections: loose, they can be tightened when the derived tables are read. + Configurable zVtxMax{"zVtxMax", 10.0f, "maximum z position of the primary vertex"}; + Configurable etaMax{"etaMax", 0.9f, "maximum eta"}; + Configurable minNCrossedRows{"minNCrossedRows", 100, "Minimum number of crossed TPC rows"}; + Configurable minNITSClus{"minNITSClus", 0., "Minimum number of ITS clusters"}; + Configurable minNTPCClus{"minNTPCClus", 80, "Minimum number of TPC clusters"}; + Configurable maxNSharedTPCClus{"maxNSharedTPCClus", 5, "Maximum number of shared TPC clusters"}; + Configurable nSigmaTPCCut{"nSigmaTPCCut", 3.f, "Number of sigmas for the TPC PID"}; + Configurable dcaBachToPV{"dcaBachToPV", 0.05f, "Minimum DCA of a cascade bachelor to the primary vertex"}; + Configurable dcaV0DauToPV{"dcaV0DauToPV", 0.05f, "Minimum DCA of Lambda daughters to the primary vertex"}; + Configurable dcaV0Bach{"dcaV0Bach", 1.f, "Maximum DCA between the V0 and cascade bachelor"}; + Configurable dcaLambdaDaughters{"dcaLambdaDaughters", 1.f, "Maximum DCA between Lambda daughters"}; + Configurable mOmegaWindow{"mOmegaWindow", 0.01f, "Omega mass window"}; + Configurable mLambdaWindow{"mLambdaWindow", 0.01f, "Lambda mass window"}; + + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + + void init(InitContext const&) + { + ccdb->setURL(ccdburl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setFatalWhenNull(false); + + auto& selections = helper.selections; + selections.etaMax = etaMax; + selections.minNITSClus = minNITSClus; + selections.minNTPCClus = minNTPCClus; + selections.minNCrossedRows = minNCrossedRows; + selections.maxNSharedTPCClus = maxNSharedTPCClus; + selections.nSigmaTPCCut = nSigmaTPCCut; + selections.dcaBachToPV = dcaBachToPV; + selections.dcaV0DauToPV = dcaV0DauToPV; + selections.dcaV0Bach = dcaV0Bach; + selections.dcaLambdaDaughters = dcaLambdaDaughters; + selections.mOmegaWindow = mOmegaWindow; + selections.mLambdaWindow = mLambdaWindow; + helper.initFitters(); + helper.setMatCorrType(static_cast(cfgMaterialCorrection.value)); + + zorroSummary.setObject(zorro.getZorroSummary()); + + histos.add("QA/zVtx", ";#it{z}_{vtx} (cm);Entries", HistType::kTH1F, {zVtxAxis}); + histos.add("QA/massXi", ";#it{p}_{T} (GeV/#it{c});#it{M}(#Lambda + #pi) (GeV/#it{c}^{2});Entries", HistType::kTH2F, {momAxis, massXiAxis}); + histos.add("QA/massOmega", ";#it{p}_{T} (GeV/#it{c});#it{M}(#Lambda + K) (GeV/#it{c}^{2});Entries", HistType::kTH2F, {momAxis, massOmegaAxis}); + histos.add("QA/massLambda", ";#it{p}_{T} (GeV/#it{c});#it{M}(p + #pi) (GeV/#it{c}^{2});Entries", HistType::kTH2F, {momAxis, massLambdaAxis}); + histos.add("QA/collisions", ";;Collisions", HistType::kTH1F, {{kNCollisionSteps, -0.5f, static_cast(kNCollisionSteps) - 0.5f}}); + auto collisions = histos.get(HIST("QA/collisions")); + collisions->GetXaxis()->SetBinLabel(kAllCollisions + 1, "All"); + collisions->GetXaxis()->SetBinLabel(kSelectedCollisions + 1, "Event selection"); + collisions->GetXaxis()->SetBinLabel(kWithOmega + 1, "With #Omega"); + collisions->GetXaxis()->SetBinLabel(kStoredCollisions + 1, "With #Omega, #Lambda and K (stored)"); + } + + template + void initCCDB(Bc const& bc) + { + if (mRunNumber == bc.runNumber()) { + return; + } + auto* grpmag = ccdb->getForTimeStamp("GLO/Config/GRPMagField", bc.timestamp()); + if (!grpmag) { + LOG(fatal) << "Could not retrieve GRPMagField for timestamp " << bc.timestamp(); + } + o2::base::Propagator::initFieldFromGRP(grpmag); + mBz = o2::base::Propagator::Instance()->getNominalBz(); + helper.setBz(mBz); + + if (static_cast(cfgMaterialCorrection.value) == o2::base::Propagator::MatCorrType::USEMatCorrLUT) { + auto* lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->getForTimeStamp("GLO/Param/MatLUT", bc.timestamp())); + if (!lut) { + LOG(fatal) << "Could not retrieve material LUT for timestamp " << bc.timestamp(); + } + o2::base::Propagator::Instance()->setMatLUT(lut); + } + + mRunNumber = bc.runNumber(); + if (cfgSkimmedProcessing) { + zorro.initCCDB(ccdb.service, bc.runNumber(), bc.timestamp(), "fDoubleOmega,fOmegaXi"); + zorro.populateHistRegistry(histos, bc.runNumber()); + } + } + + template + bool acceptCollision(C const& collision) + { + auto bc = collision.template bc_as(); + initCCDB(bc); + + if (cfgApplyEventSelection && + (!collision.sel8() || + std::abs(collision.posZ()) > zVtxMax || + !collision.selection_bit(aod::evsel::kNoITSROFrameBorder) || + !collision.selection_bit(aod::evsel::kNoTimeFrameBorder))) { + return false; + } + if (cfgSkimmedProcessing) { + zorro.isSelected(bc.globalBC()); + } + histos.fill(HIST("QA/zVtx"), collision.posZ()); + return true; + } + + void fillCandidateQA(CollisionCandidates const& candidates) + { + for (const auto& omega : candidates.omegas) { + const float omegaPt = std::hypot(omega.px, omega.py); + histos.fill(HIST("QA/massXi"), omegaPt, omega.massXi); + histos.fill(HIST("QA/massOmega"), omegaPt, omega.massOmega); + } + for (const auto& lambdas : candidates.lambdas) { + for (const auto& lambda : lambdas) { + histos.fill(HIST("QA/massLambda"), std::hypot(lambda.candidate.px, lambda.candidate.py), lambda.candidate.mass); + } + } + } + + // Writes the candidates of a collision. Tracks are written once per collision: the direct-kaon + // candidates first, then the Omega and Lambda daughters not written yet. + template + void writeCollision(C const& collision, CollisionCandidates const& candidates, T const& tracks) + { + dOmegaColls(collision.posX(), collision.posY(), collision.posZ(), collision.multFT0M(), mRunNumber, mBz); + const auto collisionIndex = static_cast(dOmegaColls.lastIndex()); + + std::unordered_map trackRows; // Track index -> DOmegaTracks row. + const auto writeTrack = [&](int64_t trackId, bool isKaon) { + const auto [row, inserted] = trackRows.try_emplace(trackId, -1); + if (!inserted) { + return row->second; + } + const auto track = tracks.rawIteratorAt(trackId); + float dcaXY = -999.f; + helper.dcaXYToPV(collision, track, dcaXY); + const auto trackParCov = getTrackParCov(track); + const auto& cov = trackParCov.getCov(); + dOmegaTracks(collisionIndex, + trackParCov.getX(), trackParCov.getAlpha(), trackParCov.getY(), trackParCov.getZ(), trackParCov.getSnp(), trackParCov.getTgl(), trackParCov.getQ2Pt(), + cov[0], cov[1], cov[2], cov[3], cov[4], cov[5], cov[6], cov[7], cov[8], cov[9], cov[10], cov[11], cov[12], cov[13], cov[14], + track.hasTPC() ? track.tpcNSigmaKa() : -999.f, + track.hasTOF() ? track.tofNSigmaKa() : -999.f, + static_cast(DoubleOmegaBuilderHelper::detectorMap(track)), + dcaXY, + isKaon); + row->second = static_cast(dOmegaTracks.lastIndex()); + return row->second; + }; + + for (const auto& kaons : candidates.kaons) { + for (const auto& kaon : kaons) { + writeTrack(kaon.trackId, true); + } + } + for (const auto& omega : candidates.omegas) { + const auto posTrackRow = writeTrack(omega.posTrackId, false); + const auto negTrackRow = writeTrack(omega.negTrackId, false); + const auto bachelorRow = writeTrack(omega.bachelorId, false); + const auto& parent = omega.parentTrack; + const auto& cov = parent.getCov(); + dOmegaOmegas(collisionIndex, posTrackRow, negTrackRow, bachelorRow, + parent.getX(), parent.getAlpha(), parent.getY(), parent.getZ(), parent.getSnp(), parent.getTgl(), parent.getQ2Pt(), + cov[0], cov[1], cov[2], cov[3], cov[4], cov[5], cov[6], cov[7], cov[8], cov[9], cov[10], cov[11], cov[12], cov[13], cov[14], + omega.x, omega.y, omega.z, + omega.px, omega.py, omega.pz, + omega.massOmega, omega.massXi, + omega.dcaDaughters, + omega.sign); + } + for (const int8_t charge : {int8_t{-1}, int8_t{1}}) { + for (const auto& lambda : candidates.lambdas[DoubleOmegaBuilderHelper::chargeIndex(charge)]) { + const auto posTrackRow = writeTrack(lambda.posTrackId, false); + const auto negTrackRow = writeTrack(lambda.negTrackId, false); + const auto& parent = lambda.parentTrack; + const auto& cov = parent.getCov(); + dOmegaLambdas(collisionIndex, posTrackRow, negTrackRow, + parent.getX(), parent.getAlpha(), parent.getY(), parent.getZ(), parent.getSnp(), parent.getTgl(), parent.getQ2Pt(), + cov[0], cov[1], cov[2], cov[3], cov[4], cov[5], cov[6], cov[7], cov[8], cov[9], cov[10], cov[11], cov[12], cov[13], cov[14], + lambda.decayVertex[0], lambda.decayVertex[1], lambda.decayVertex[2], + lambda.candidate.px, lambda.candidate.py, lambda.candidate.pz, + lambda.candidate.mass, + lambda.dcaDaughters, + charge); + } + } + } + + void processData(CollisionsTable const& collisions, + TracksFull const& tracks, + aod::V0s const& v0s, + aod::Cascades const& cascades, + aod::BCsWithTimestamps const&) + { + for (const auto& collision : collisions) { + histos.fill(HIST("QA/collisions"), kAllCollisions); + if (!acceptCollision(collision)) { + continue; + } + histos.fill(HIST("QA/collisions"), kSelectedCollisions); + auto tracksThisCollision = tracks.sliceBy(tracksPerCollision, collision.globalIndex()); + auto cascadesThisCollision = cascades.sliceBy(cascadesPerCollision, collision.globalIndex()); + auto v0sThisCollision = v0s.sliceBy(v0sPerCollision, collision.globalIndex()); + v0sThisCollision.bindExternalIndices(&tracks); + cascadesThisCollision.bindExternalIndices(&tracks); + cascadesThisCollision.bindExternalIndices(&v0s); + const auto candidates = helper.buildCollisionCandidates(collision, tracksThisCollision, cascadesThisCollision, v0sThisCollision, true); + fillCandidateQA(candidates); + if (candidates.omegas.empty()) { + continue; + } + histos.fill(HIST("QA/collisions"), kWithOmega); + if (!DoubleOmegaBuilderHelper::hasOmegaLambdaKaon(candidates)) { + continue; + } + histos.fill(HIST("QA/collisions"), kStoredCollisions); + writeCollision(collision, candidates, tracks); + } + } + PROCESS_SWITCH(doubleOmegaBuilder, processData, "Build the Omega, Lambda and kaon candidates and store them in derived tables", true); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} diff --git a/PWGLF/TableProducer/Strangeness/doubleOmegaTreeCreator.cxx b/PWGLF/TableProducer/Strangeness/doubleOmegaTreeCreator.cxx index 6ce665dee2d..8622f2e626a 100644 --- a/PWGLF/TableProducer/Strangeness/doubleOmegaTreeCreator.cxx +++ b/PWGLF/TableProducer/Strangeness/doubleOmegaTreeCreator.cxx @@ -11,6 +11,7 @@ #include "PWGLF/DataModel/LFDoubleOmegaTables.h" #include "PWGLF/DataModel/LFStrangenessTables.h" +#include "PWGLF/Utils/doubleOmegaBuilderHelper.h" #include "Common/CCDB/EventSelectionParams.h" #include "Common/Core/RecoDecay.h" @@ -23,14 +24,16 @@ #include "Common/DataModel/PIDResponseTPC.h" #include +#include #include -#include #include #include #include +#include #include #include #include +#include #include #include #include @@ -43,20 +46,29 @@ #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 namespace o2::pwglf::doubleomega; using CollisionsTable = soa::Join; using Collisions = CollisionsTable::iterator; @@ -67,38 +79,6 @@ using FullCascades = aod::CascDataExt; using TracksFull = soa::Join; using TracksFullMC = soa::Join; -struct DoubleOmegaCandidate { - float pt = -999.f; - float eta = -999.f; - float phi = -999.f; - float x = -999.f; - float y = -999.f; - float z = -999.f; - float pvX = -999.f; - float pvY = -999.f; - float pvZ = -999.f; - float cosPAOmega = -999.f; - float cosPADirectLambda = -999.f; - float cosPADoubleOmega = -999.f; - float cosPAOmegaToPV = -999.f; - float cosPADirectLambdaToPV = -999.f; - float dcaXYOmegaToPV = -999.f; - float dcaXYDirectLambdaToPV = -999.f; - float dcaXYDirectKaonToPV = -999.f; - float dcaXYOmegaKaonToPV = -999.f; - float dcaXYOmegaProtonToPV = -999.f; - float dcaXYOmegaPionToPV = -999.f; - float dcaXYDirectProtonToPV = -999.f; - float dcaXYDirectPionToPV = -999.f; - float tpcNSigmaDirectKaon = -999.f; - float tofNSigmaDirectKaon = -999.f; - uint32_t daughterDetectorMap = 0; - float mass = -999.f; - float massOmega = -999.f; - float massXi = -999.f; - int8_t sign = 0; -}; - struct DoubleOmegaMCInfo { int64_t motherId = -1; float pt = -999.f; @@ -114,47 +94,22 @@ struct DoubleOmegaMCInfo { int8_t sign = 0; }; -struct LambdaCandidate { - float px = 0.f; - float py = 0.f; - float pz = 0.f; - float mass = -999.f; -}; - -struct BuiltLambdaCandidate { - LambdaCandidate candidate; - o2::track::TrackParCov parentTrack; - std::array decayVertex{}; - int64_t posTrackId = -1; - int64_t negTrackId = -1; - uint32_t daughterDetectorMap = 0; // Proton, pion; three detector bits each. - float dcaXYProtonToPV = -999.f; - float dcaXYPionToPV = -999.f; -}; - -struct BuiltOmegaCandidate { - float px = 0.f; - float py = 0.f; - float pz = 0.f; - float x = 0.f; - float y = 0.f; - float z = 0.f; - float massOmega = -999.f; - float massXi = -999.f; - int8_t sign = 0; - int64_t posTrackId = -1; - int64_t negTrackId = -1; - int64_t bachelorId = -1; - o2::track::TrackParCov parentTrack; - uint32_t daughterDetectorMap = 0; // Kaon, proton, pion; three detector bits each. - float dcaXYKaonToPV = -999.f; - float dcaXYProtonToPV = -999.f; - float dcaXYPionToPV = -999.f; +// Same-event double-Omega candidate used as a seed of the candidate-level event mixing. Seeds are +// kept in a pool across dataframes, so they store copies of their daughters, with the collision and +// track indices made unique across dataframes (see uniqueIndex). +struct MixingSeed { + int64_t collisionId = -1; + std::array primaryVertex{}; + float bz = 0.f; + BuiltOmegaCandidate omega; + BuiltLambdaCandidate lambda; + BuiltKaonCandidate kaon; }; struct doubleOmegaTreeCreator { Produces doubleOmegaTable; Produces doubleOmegaTableMC; + Produces doubleOmegaMixedTable; Service ccdb{}; static constexpr int kDoubleOmegaPdg = 1060020020; @@ -169,16 +124,26 @@ struct doubleOmegaTreeCreator { kFindableCascadeAndV0SelectedTracks, kNFindabilitySteps }; + enum MixingType : uint8_t { + kSameEvent = 0, + kMixedOmega, // Lambda and kaon from the same collision, Omega from another one + kMixedLambda, // Omega and kaon from the same collision, Lambda from another one + kMixedKaon, // Omega and Lambda from the same collision, kaon from another one + kMixedAll, // Omega, Lambda and kaon from three different collisions + kNMixingTypes + }; std::vector reconstructedDoubleOmegaIds; Preslice tracksPerCollision = aod::track::collisionId; Preslice rawCascadesPerCollision = aod::cascade::collisionId; Preslice rawV0sPerCollision = aod::v0::collisionId; + Preslice derivedTracksPerCollision = aod::domega::dOmegaCollId; + Preslice derivedOmegasPerCollision = aod::domega::dOmegaCollId; + Preslice derivedLambdasPerCollision = aod::domega::dOmegaCollId; int mRunNumber = 0; float mBz = 0.f; - o2::vertexing::DCAFitterN<2> fitter; - o2::vertexing::DCAFitterN<3> fitter3Body; + DoubleOmegaBuilderHelper helper; Zorro zorro; OutputObj zorroSummary{"zorroSummary"}; @@ -221,32 +186,25 @@ struct doubleOmegaTreeCreator { Configurable minCosPADoubleOmega{"minCosPADoubleOmega", -1.f, "Minimum double-Omega cosPA relative to the primary vertex"}; Configurable minDoubleOmegaDecayRadius{"minDoubleOmegaDecayRadius", 1.f, "Minimum double-Omega transverse decay radius in cm"}; - HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; - - template - static uint32_t detectorMap(Track const& track) - { - return static_cast(track.hasITS()) | - (static_cast(track.hasTPC()) << 1) | - (static_cast(track.hasTOF()) << 2); - } + struct : ConfigurableGroup { + std::string prefix = "mixingOpts"; + Configurable nMixingCollisions{"nMixingCollisions", 5, "Number of neighbouring collisions in the same bin each collision is mixed with"}; + Configurable mixingTypes{"mixingTypes", 0b1111, "Bitmask of the mixed backgrounds: 1 Omega, 2 Lambda, 4 kaon from another collision, 8 all from different collisions"}; + Configurable maxDeltaPVz{"maxDeltaPVz", 1.f, "Maximum difference of the primary-vertex z among the mixed collisions (cm)"}; + ConfigurableAxis binsVtxZ{"binsVtxZ", {VARIABLE_WIDTH, -10.f, -8.f, -6.f, -4.f, -2.f, 0.f, 2.f, 4.f, 6.f, 8.f, 10.f}, "Mixing bins - primary-vertex z (cm)"}; + ConfigurableAxis binsMultiplicity{"binsMultiplicity", {VARIABLE_WIDTH, 0.f, 20.f, 50.f, 100.f, 200.f, 500.f, 1.e5f}, "Mixing bins - FT0M multiplicity"}; + // Candidate-level mixing (processDerivedMixingCandidates) + Configurable nMixingCandidates{"nMixingCandidates", 10, "Number of same-event candidates from other collisions each candidate is mixed with, per (R, phi, charge) bin"}; + Configurable poolDepth{"poolDepth", 50, "Number of same-event candidates kept in the mixing pool per (R, phi, charge) bin, across dataframes"}; + Configurable monitoringInterval{"monitoringInterval", 50, "Log the mixing pool size, memory and timing every this many dataframes (0: never)"}; + ConfigurableAxis binsRadius{"binsRadius", {VARIABLE_WIDTH, 0.f, 0.25f, 0.5f, 1.f, 2.f, 4.f, 10.f, 40.f}, "Mixing bins - double-Omega decay radius (cm)"}; + ConfigurableAxis binsPhi{"binsPhi", {VARIABLE_WIDTH, -o2::constants::math::PI, -2.f * o2::constants::math::PI / 3.f, -o2::constants::math::PI / 3.f, 0.f, o2::constants::math::PI / 3.f, 2.f * o2::constants::math::PI / 3.f, o2::constants::math::PI}, "Mixing bins - double-Omega momentum azimuth (rad)"}; + } mixingOpts; + + using DerivedMixingBinning = ColumnBinningPolicy; + static constexpr float kBzTolerance = 1.e-3f; // kG - template - bool selectTrack(T const& track) - { - if (std::abs(track.eta()) > etaMax) { - return false; - } - if (track.itsNCls() < minNITSClus || - track.tpcNClsFound() < minNTPCClus || - track.tpcNClsCrossedRows() < minNCrossedRows || - track.tpcNClsCrossedRows() < 0.8f * track.tpcNClsFindable() || - track.tpcNClsShared() > maxNSharedTPCClus) { - LOG(debug) << "Track failed selection: itsNCls=" << static_cast(track.itsNCls()) << ", tpcNClsFound=" << track.tpcNClsFound() << ", tpcNClsCrossedRows=" << static_cast(track.tpcNClsCrossedRows()) << ", tpcNClsFindable=" << static_cast(track.tpcNClsFindable()) << ", tpcNClsShared=" << static_cast(track.tpcNClsShared()); - return false; - } - return true; - } + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; template void initCCDB(Bc const& bc) @@ -260,8 +218,7 @@ struct doubleOmegaTreeCreator { } o2::base::Propagator::initFieldFromGRP(grpmag); mBz = o2::base::Propagator::Instance()->getNominalBz(); - fitter.setBz(mBz); - fitter3Body.setBz(mBz); + helper.setBz(mBz); if (static_cast(cfgMaterialCorrection.value) == o2::base::Propagator::MatCorrType::USEMatCorrLUT) { auto* lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->getForTimeStamp("GLO/Param/MatLUT", bc.timestamp())); @@ -270,8 +227,7 @@ struct doubleOmegaTreeCreator { } o2::base::Propagator::Instance()->setMatLUT(lut); } - fitter.setMatCorrType(static_cast(cfgMaterialCorrection.value)); - fitter3Body.setMatCorrType(static_cast(cfgMaterialCorrection.value)); + helper.setMatCorrType(static_cast(cfgMaterialCorrection.value)); mRunNumber = bc.runNumber(); if (cfgSkimmedProcessing) { @@ -332,10 +288,23 @@ struct doubleOmegaTreeCreator { printValue(minCosPADirectLambda); printValue(minCosPADoubleOmega); printValue(minDoubleOmegaDecayRadius); + printValue(mixingOpts.nMixingCollisions); + printValue(mixingOpts.mixingTypes); + printValue(mixingOpts.maxDeltaPVz); + printAxis(mixingOpts.binsVtxZ); + printAxis(mixingOpts.binsMultiplicity); + printValue(mixingOpts.nMixingCandidates); + printValue(mixingOpts.poolDepth); + printValue(mixingOpts.monitoringInterval); + printAxis(mixingOpts.binsRadius); + printAxis(mixingOpts.binsPhi); printValue(doprocessData); printValue(doprocessDataFromCascades); printValue(doprocessMC); printValue(doprocessFindableTracks); + printValue(doprocessDerivedData); + printValue(doprocessDerivedMixing); + printValue(doprocessDerivedMixingCandidates); } void init(InitContext const&) @@ -347,19 +316,29 @@ struct doubleOmegaTreeCreator { ccdb->setLocalObjectValidityChecking(); ccdb->setFatalWhenNull(false); - fitter.setPropagateToPCA(true); - fitter.setMaxR(200.); - fitter.setMaxDZIni(4.); - fitter.setMinParamChange(1.e-3); - fitter.setMinRelChi2Change(0.9); - fitter.setUseAbsDCA(true); - - fitter3Body.setPropagateToPCA(true); - fitter3Body.setMaxR(200.); - fitter3Body.setMaxDZIni(4.); - fitter3Body.setMinParamChange(1.e-3); - fitter3Body.setMinRelChi2Change(0.9); - fitter3Body.setUseAbsDCA(true); + auto& selections = helper.selections; + selections.etaMax = etaMax; + selections.minNITSClus = minNITSClus; + selections.minNTPCClus = minNTPCClus; + selections.minNCrossedRows = minNCrossedRows; + selections.maxNSharedTPCClus = maxNSharedTPCClus; + selections.nSigmaTPCCut = nSigmaTPCCut; + selections.dcaBachToPV = dcaBachToPV; + selections.dcaKaonToPV = dcaKaonToPV; + selections.dcaOmegaToPV = dcaOmegaToPV; + selections.dcaDirectLambdaToPV = dcaDirectLambdaToPV; + selections.dcaV0DauToPV = dcaV0DauToPV; + selections.dcaV0Bach = dcaV0Bach; + selections.dcaLambdaDaughters = dcaLambdaDaughters; + selections.mOmegaWindow = mOmegaWindow; + selections.mLambdaWindow = mLambdaWindow; + selections.maxDoubleOmegaMass = maxDoubleOmegaMass; + selections.minCosPAOmega = minCosPAOmega; + selections.minCosPADirectLambda = minCosPADirectLambda; + selections.minCosPADoubleOmega = minCosPADoubleOmega; + selections.minDoubleOmegaDecayRadius = minDoubleOmegaDecayRadius; + helper.initFitters(); + helper.setMatCorrType(static_cast(cfgMaterialCorrection.value)); zorroSummary.setObject(zorro.getZorroSummary()); @@ -377,302 +356,29 @@ struct doubleOmegaTreeCreator { generatedAndFindable->GetXaxis()->SetBinLabel(kFindableSelectedTracks + 1, "Findable, selected tracks"); generatedAndFindable->GetXaxis()->SetBinLabel(kFindableCascadeAndV0 + 1, "Findable cascade + V0"); generatedAndFindable->GetXaxis()->SetBinLabel(kFindableCascadeAndV0SelectedTracks + 1, "Findable cascade + V0, selected tracks"); - } - - template - bool isSelectedKaon(Track const& track, int8_t charge) - { - return track.sign() == charge && - selectTrack(track) && - std::abs(track.tpcNSigmaKa()) <= nSigmaTPCCut; - } - - static float invariantMass2Body(std::array const& momentum1, float mass1, - std::array const& momentum2, float mass2) - { - const float momentum1Squared = momentum1[0] * momentum1[0] + momentum1[1] * momentum1[1] + momentum1[2] * momentum1[2]; - const float momentum2Squared = momentum2[0] * momentum2[0] + momentum2[1] * momentum2[1] + momentum2[2] * momentum2[2]; - const float energy = std::sqrt(momentum1Squared + mass1 * mass1) + std::sqrt(momentum2Squared + mass2 * mass2); - const std::array totalMomentum{ - momentum1[0] + momentum2[0], momentum1[1] + momentum2[1], momentum1[2] + momentum2[2]}; - const float massSquared = energy * energy - totalMomentum[0] * totalMomentum[0] - totalMomentum[1] * totalMomentum[1] - totalMomentum[2] * totalMomentum[2]; - return std::sqrt(std::max(0.f, massSquared)); - } - - // Transverse DCA of a track to the collision primary vertex, computed by propagating - // a copy of the track parameters. Returns false if the propagation fails. - template - bool dcaXYToPV(C const& collision, Track const& track, float& dcaXY) const - { - auto trackParCov = getTrackParCov(track); - const o2::math_utils::Point3D primaryVertex{collision.posX(), collision.posY(), collision.posZ()}; - const auto matCorr = static_cast(cfgMaterialCorrection.value); - std::array dca{}; - if (!o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertex, trackParCov, 2.f, matCorr, &dca)) { - return false; - } - dcaXY = dca[0]; - return true; - } - - template - bool buildLambda(C const& collision, Track const& posTrack, Track const& negTrack, int8_t charge, BuiltLambdaCandidate& builtLambda) - { - if (posTrack.sign() <= 0 || negTrack.sign() >= 0 || - !selectTrack(posTrack) || !selectTrack(negTrack)) { - return false; - } - - const bool isMatter = charge < 0; - const float protonNSigma = isMatter ? posTrack.tpcNSigmaPr() : negTrack.tpcNSigmaPr(); - const float pionNSigma = isMatter ? negTrack.tpcNSigmaPi() : posTrack.tpcNSigmaPi(); - if (std::abs(protonNSigma) > nSigmaTPCCut || std::abs(pionNSigma) > nSigmaTPCCut) { - return false; - } - float dcaXYPos = -999.f; - float dcaXYNeg = -999.f; - if (!dcaXYToPV(collision, posTrack, dcaXYPos) || !dcaXYToPV(collision, negTrack, dcaXYNeg)) { - LOG(debug) << "Propagation to the PV failed for Lambda daughters " << posTrack.globalIndex() << ", " << negTrack.globalIndex(); - return false; - } - if (std::abs(dcaXYPos) < dcaV0DauToPV || std::abs(dcaXYNeg) < dcaV0DauToPV) { - LOG(debug) << "Lambda dau don't pass selections: " << "nsigmaPR: " << protonNSigma << ", pionNSigma: " << pionNSigma << ", dcaXYPos: " << dcaXYPos << ", dcaXYNeg: " << dcaXYNeg; - return false; - } - - auto posTrackParCov = getTrackParCov(posTrack); - auto negTrackParCov = getTrackParCov(negTrack); - int nCandidates = 0; - try { - nCandidates = fitter.process(posTrackParCov, negTrackParCov); - } catch (...) { - LOG(error) << "Exception while fitting Lambda daughter tracks " << posTrack.globalIndex() << ", " << negTrack.globalIndex(); - return false; - } - if (nCandidates == 0) { - return false; - } - - std::array posMomentum{}; - std::array negMomentum{}; - fitter.getTrack(0).getPxPyPzGlo(posMomentum); - fitter.getTrack(1).getPxPyPzGlo(negMomentum); - const std::array lambdaMomentum{ - posMomentum[0] + negMomentum[0], - posMomentum[1] + negMomentum[1], - posMomentum[2] + negMomentum[2]}; - const auto& fittedVertex = fitter.getPCACandidate(); - const std::array decayVertex{ - static_cast(fittedVertex[0]), static_cast(fittedVertex[1]), static_cast(fittedVertex[2])}; - const float mass = isMatter ? invariantMass2Body(posMomentum, o2::constants::physics::MassProton, - negMomentum, o2::constants::physics::MassPionCharged) - : invariantMass2Body(posMomentum, o2::constants::physics::MassPionCharged, - negMomentum, o2::constants::physics::MassProton); - const float dcaDaughters = std::sqrt(std::abs(fitter.getChi2AtPCACandidate())); - const float eta = etaFromMomentum(lambdaMomentum[0], lambdaMomentum[1], lambdaMomentum[2]); - if (std::abs(eta) > etaMax || - dcaDaughters > dcaLambdaDaughters || - std::abs(mass - o2::constants::physics::MassLambda0) > mLambdaWindow) { - return false; - } - - builtLambda.candidate = { - .px = lambdaMomentum[0], - .py = lambdaMomentum[1], - .pz = lambdaMomentum[2], - .mass = mass}; - builtLambda.parentTrack = fitter.createParentTrackParCov(0); - builtLambda.decayVertex = decayVertex; - builtLambda.posTrackId = posTrack.globalIndex(); - builtLambda.negTrackId = negTrack.globalIndex(); - // Keep particle-role ordering for matter and antimatter alike. - builtLambda.daughterDetectorMap = isMatter ? detectorMap(posTrack) | (detectorMap(negTrack) << 3) - : detectorMap(negTrack) | (detectorMap(posTrack) << 3); - builtLambda.dcaXYProtonToPV = isMatter ? dcaXYPos : dcaXYNeg; - builtLambda.dcaXYPionToPV = isMatter ? dcaXYNeg : dcaXYPos; - return true; - } - - template - bool buildOmega(C const& collision, Track const& posTrack, Track const& negTrack, Track const& bachelor, BuiltOmegaCandidate& builtOmega) - { - if (bachelor.sign() == 0 || - !selectTrack(bachelor) || - std::abs(bachelor.tpcNSigmaKa()) > nSigmaTPCCut) { - return false; - } - float dcaXYBachelor = -999.f; - if (!dcaXYToPV(collision, bachelor, dcaXYBachelor) || std::abs(dcaXYBachelor) < dcaBachToPV) { - return false; - } - - BuiltLambdaCandidate lambda; - if (!buildLambda(collision, posTrack, negTrack, bachelor.sign(), lambda)) { - return false; - } - auto v0TrackParCov = lambda.parentTrack; - auto bachelorTrackParCov = getTrackParCov(bachelor); - int nCandidates = 0; - try { - nCandidates = fitter.process(v0TrackParCov, bachelorTrackParCov); - } catch (...) { - LOG(error) << "Exception while fitting Omega daughter tracks " << posTrack.globalIndex() << ", " << negTrack.globalIndex() << ", " << bachelor.globalIndex(); - return false; - } - if (nCandidates == 0) { - return false; + if (doprocessDerivedMixingCandidates) { + histos.add("Mixing/hSeedRadiusPhi", ";Double-#Omega decay radius (cm);Double-#Omega #varphi (rad);Same-event candidates", HistType::kTH2F, {mixingOpts.binsRadius, mixingOpts.binsPhi}); + // Resource monitoring of the candidate-level mixing, one entry per dataframe. + histos.add("Mixing/hProcessingTime", ";Processing time per dataframe (ms);Dataframes", HistType::kTH1F, {{500, 0.f, 5000.f}}); + histos.add("Mixing/hPoolSeeds", ";Seeds in the mixing pool;Dataframes", HistType::kTH1F, {{500, 0.f, 50000.f}}); + histos.add("Mixing/hPoolMemory", ";Estimated mixing-pool memory (MB);Dataframes", HistType::kTH1F, {{500, 0.f, 50.f}}); + histos.add("Mixing/hProcessingTimeVsPoolSeeds", ";Seeds in the mixing pool;Processing time per dataframe (ms)", HistType::kTH2F, {{100, 0.f, 50000.f}, {100, 0.f, 5000.f}}); } - - std::array lambdaMomentum{}; - std::array bachelorMomentum{}; - fitter.getTrack(0).getPxPyPzGlo(lambdaMomentum); - fitter.getTrack(1).getPxPyPzGlo(bachelorMomentum); - const std::array omegaMomentum{ - lambdaMomentum[0] + bachelorMomentum[0], - lambdaMomentum[1] + bachelorMomentum[1], - lambdaMomentum[2] + bachelorMomentum[2]}; - const auto& fittedVertex = fitter.getPCACandidate(); - const std::array decayVertex{ - static_cast(fittedVertex[0]), static_cast(fittedVertex[1]), static_cast(fittedVertex[2])}; - const float massOmega = invariantMass2Body(lambdaMomentum, o2::constants::physics::MassLambda0, - bachelorMomentum, o2::constants::physics::MassKaonCharged); - const float massXi = invariantMass2Body(lambdaMomentum, o2::constants::physics::MassLambda0, - bachelorMomentum, o2::constants::physics::MassPionCharged); - const float dcaDaughters = std::sqrt(std::abs(fitter.getChi2AtPCACandidate())); - const float eta = etaFromMomentum(omegaMomentum[0], omegaMomentum[1], omegaMomentum[2]); - if (dcaDaughters > dcaV0Bach || - std::abs(eta) > etaMax || - std::abs(massOmega - o2::constants::physics::MassOmegaMinus) > mOmegaWindow) { - return false; + if (doprocessDerivedMixing || doprocessDerivedMixingCandidates) { + histos.add("Mixing/hMixedCollisions", ";;Entries", HistType::kTH1F, {{4, -0.5f, 3.5f}}); + auto mixedCollisions = histos.get(HIST("Mixing/hMixedCollisions")); + mixedCollisions->GetXaxis()->SetBinLabel(1, "Pairs"); + mixedCollisions->GetXaxis()->SetBinLabel(2, "Pairs, #Delta#it{z}_{PV}, same #it{B}"); + mixedCollisions->GetXaxis()->SetBinLabel(3, "Triples"); + mixedCollisions->GetXaxis()->SetBinLabel(4, "Triples, #Delta#it{z}_{PV}, same #it{B}"); + histos.add("Mixing/hMixedCandidates", ";;Entries", HistType::kTH1F, {{kNMixingTypes - 1, 0.5f, static_cast(kNMixingTypes) - 0.5f}}); + auto mixedCandidates = histos.get(HIST("Mixing/hMixedCandidates")); + mixedCandidates->GetXaxis()->SetBinLabel(kMixedOmega, "#Omega mixed"); + mixedCandidates->GetXaxis()->SetBinLabel(kMixedLambda, "#Lambda mixed"); + mixedCandidates->GetXaxis()->SetBinLabel(kMixedKaon, "K mixed"); + mixedCandidates->GetXaxis()->SetBinLabel(kMixedAll, "All mixed"); } - - builtOmega = { - omegaMomentum[0], - omegaMomentum[1], - omegaMomentum[2], - decayVertex[0], - decayVertex[1], - decayVertex[2], - massOmega, - massXi, - static_cast(bachelor.sign()), - posTrack.globalIndex(), - negTrack.globalIndex(), - bachelor.globalIndex(), - fitter.createParentTrackParCov(), - detectorMap(bachelor) | (lambda.daughterDetectorMap << 3), - dcaXYBachelor, - lambda.dcaXYProtonToPV, - lambda.dcaXYPionToPV}; - return true; - } - - static float etaFromMomentum(float px, float py, float pz) - { - const float pt = std::hypot(px, py); - return pt > 0.f ? std::asinh(pz / pt) : 0.f; - } - - template - bool buildDoubleOmega(C const& collision, BuiltOmegaCandidate const& omega, BuiltLambdaCandidate const& directLambda, Kaon const& kaon, DoubleOmegaCandidate& cand) - { - auto omegaTrackParCov = omega.parentTrack; - auto lambdaTrackParCov = directLambda.parentTrack; - auto kaonTrackParCov = getTrackParCov(kaon); - int nCandidates = 0; - try { - nCandidates = fitter3Body.process(omegaTrackParCov, lambdaTrackParCov, kaonTrackParCov); - } catch (...) { - LOG(error) << "Exception while fitting the double-Omega candidate"; - return false; - } - if (nCandidates == 0) { - return false; - } - - auto omegaTrackAtVertex = fitter3Body.getTrack(0); - auto lambdaTrackAtVertex = fitter3Body.getTrack(1); - auto kaonTrackAtVertex = fitter3Body.getTrack(2); - std::array omegaMomentum{}; - std::array lambdaMomentum{}; - std::array kaonMomentum{}; - omegaTrackAtVertex.getPxPyPzGlo(omegaMomentum); - lambdaTrackAtVertex.getPxPyPzGlo(lambdaMomentum); - kaonTrackAtVertex.getPxPyPzGlo(kaonMomentum); - - const std::array totalMomentum{ - omegaMomentum[0] + lambdaMomentum[0] + kaonMomentum[0], - omegaMomentum[1] + lambdaMomentum[1] + kaonMomentum[1], - omegaMomentum[2] + lambdaMomentum[2] + kaonMomentum[2]}; - const auto& fittedVertex = fitter3Body.getPCACandidate(); - const std::array decayVertex{ - static_cast(fittedVertex[0]), static_cast(fittedVertex[1]), static_cast(fittedVertex[2])}; - const std::array primaryVertex{collision.posX(), collision.posY(), collision.posZ()}; - const o2::math_utils::Point3D primaryVertexPoint{collision.posX(), collision.posY(), collision.posZ()}; - const std::array omegaDecayVertex{omega.x, omega.y, omega.z}; - - std::array dcaOmega{}; - std::array dcaDirectLambda{}; - std::array dcaDirectKaon{}; - const auto matCorr = static_cast(cfgMaterialCorrection.value); - if (!o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertexPoint, omegaTrackAtVertex, 2.f, matCorr, &dcaOmega) || - !o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertexPoint, lambdaTrackAtVertex, 2.f, matCorr, &dcaDirectLambda) || - !o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertexPoint, kaonTrackAtVertex, 2.f, matCorr, &dcaDirectKaon)) { - return false; - } - - const auto momentumSquared = [](std::array const& momentum) { - return momentum[0] * momentum[0] + momentum[1] * momentum[1] + momentum[2] * momentum[2]; - }; - const float energy = std::sqrt(momentumSquared(omegaMomentum) + o2::constants::physics::MassOmegaMinus * o2::constants::physics::MassOmegaMinus) + - std::sqrt(momentumSquared(lambdaMomentum) + o2::constants::physics::MassLambda0 * o2::constants::physics::MassLambda0) + - std::sqrt(momentumSquared(kaonMomentum) + o2::constants::physics::MassKaonCharged * o2::constants::physics::MassKaonCharged); - const float massSquared = energy * energy - totalMomentum[0] * totalMomentum[0] - totalMomentum[1] * totalMomentum[1] - totalMomentum[2] * totalMomentum[2]; - - cand.pt = std::hypot(totalMomentum[0], totalMomentum[1]); - cand.eta = etaFromMomentum(totalMomentum[0], totalMomentum[1], totalMomentum[2]); - cand.phi = std::atan2(totalMomentum[1], totalMomentum[0]); - cand.x = decayVertex[0]; - cand.y = decayVertex[1]; - cand.z = decayVertex[2]; - cand.pvX = primaryVertex[0]; - cand.pvY = primaryVertex[1]; - cand.pvZ = primaryVertex[2]; - cand.cosPAOmega = RecoDecay::cpa(decayVertex, omegaDecayVertex, omegaMomentum); - cand.cosPADirectLambda = RecoDecay::cpa(decayVertex, directLambda.decayVertex, lambdaMomentum); - cand.cosPADoubleOmega = RecoDecay::cpa(primaryVertex, decayVertex, totalMomentum); - // PV pointing angles use the daughter momenta at their own decay vertices. - cand.cosPAOmegaToPV = RecoDecay::cpa(primaryVertex, omegaDecayVertex, std::array{omega.px, omega.py, omega.pz}); - cand.cosPADirectLambdaToPV = RecoDecay::cpa(primaryVertex, directLambda.decayVertex, std::array{directLambda.candidate.px, directLambda.candidate.py, directLambda.candidate.pz}); - cand.dcaXYOmegaToPV = dcaOmega[0]; - cand.dcaXYDirectLambdaToPV = dcaDirectLambda[0]; - cand.dcaXYDirectKaonToPV = dcaDirectKaon[0]; - cand.dcaXYOmegaKaonToPV = omega.dcaXYKaonToPV; - cand.dcaXYOmegaProtonToPV = omega.dcaXYProtonToPV; - cand.dcaXYOmegaPionToPV = omega.dcaXYPionToPV; - cand.dcaXYDirectProtonToPV = directLambda.dcaXYProtonToPV; - cand.dcaXYDirectPionToPV = directLambda.dcaXYPionToPV; - cand.tpcNSigmaDirectKaon = kaon.hasTPC() ? kaon.tpcNSigmaKa() : -999.f; - cand.tofNSigmaDirectKaon = kaon.hasTOF() ? kaon.tofNSigmaKa() : -999.f; - cand.daughterDetectorMap = omega.daughterDetectorMap | - (directLambda.daughterDetectorMap << 9) | - (detectorMap(kaon) << 15); - cand.mass = std::sqrt(std::max(0.f, massSquared)); - if (cand.cosPAOmega < minCosPAOmega || - cand.cosPADirectLambda < minCosPADirectLambda || - cand.cosPADoubleOmega < minCosPADoubleOmega || - std::hypot(decayVertex[0], decayVertex[1]) < minDoubleOmegaDecayRadius || - std::abs(cand.dcaXYOmegaToPV) < dcaOmegaToPV || - std::abs(cand.dcaXYDirectLambdaToPV) < dcaDirectLambdaToPV || - std::abs(cand.dcaXYDirectKaonToPV) < dcaKaonToPV || - cand.mass > maxDoubleOmegaMass) { - return false; - } - cand.massOmega = omega.massOmega; - cand.massXi = omega.massXi; - cand.sign = omega.sign; - return true; } template @@ -882,6 +588,39 @@ struct doubleOmegaTreeCreator { cand.massXi); } + void writeMixedCandidate(DoubleOmegaCandidate const& cand, MixingType mixingType) + { + doubleOmegaMixedTable(cand.pt * cand.sign, + cand.eta, + cand.phi, + cand.x, + cand.y, + cand.z, + cand.pvX, + cand.pvY, + cand.pvZ, + cand.cosPAOmega, + cand.cosPADirectLambda, + cand.cosPADoubleOmega, + cand.cosPAOmegaToPV, + cand.cosPADirectLambdaToPV, + cand.dcaXYOmegaToPV, + cand.dcaXYDirectLambdaToPV, + cand.dcaXYDirectKaonToPV, + cand.dcaXYOmegaKaonToPV, + cand.dcaXYOmegaProtonToPV, + cand.dcaXYOmegaPionToPV, + cand.dcaXYDirectProtonToPV, + cand.dcaXYDirectPionToPV, + cand.tpcNSigmaDirectKaon, + cand.tofNSigmaDirectKaon, + cand.daughterDetectorMap, + cand.mass, + cand.massOmega, + cand.massXi, + mixingType); + } + void writeMCCandidate(DoubleOmegaCandidate const& cand, DoubleOmegaMCInfo const& mcInfo, bool isReco) { doubleOmegaTableMC(cand.pt * cand.sign, @@ -942,74 +681,34 @@ struct doubleOmegaTreeCreator { LOG(debug) << "collision " << collision.globalIndex() << " (direct kaon), numContrib = " << collision.numContrib(); } - template - void fillFromRawTables(C const& collision, T const& tracks, Cascades const& cascades, V0s const& v0s, - aod::McParticles const* mcParticles = nullptr) + void fillCandidateQA(CollisionCandidates const& candidates) { - for (const auto& rawCascade : cascades) { - const auto omegaV0 = rawCascade.template v0_as(); - const auto omegaPosTrack = omegaV0.template posTrack_as(); - const auto omegaNegTrack = omegaV0.template negTrack_as(); - const auto omegaBachelorTrack = rawCascade.template bachelor_as(); - BuiltOmegaCandidate omega; - if (!buildOmega(collision, omegaPosTrack, omegaNegTrack, omegaBachelorTrack, omega)) { - continue; - } + for (const auto& omega : candidates.omegas) { const float omegaPt = std::hypot(omega.px, omega.py); histos.fill(HIST("QA/massXi"), omegaPt, omega.massXi); histos.fill(HIST("QA/massOmega"), omegaPt, omega.massOmega); - const std::array omegaTrackIds{omega.posTrackId, omega.negTrackId, omega.bachelorId}; - - for (const auto& rawV0 : v0s) { - const auto lambdaPosTrack = rawV0.template posTrack_as(); - const auto lambdaNegTrack = rawV0.template negTrack_as(); - BuiltLambdaCandidate lambda; - if (!buildLambda(collision, lambdaPosTrack, lambdaNegTrack, omega.sign, lambda)) { - continue; - } - if (std::find(omegaTrackIds.begin(), omegaTrackIds.end(), lambda.posTrackId) != omegaTrackIds.end() || - std::find(omegaTrackIds.begin(), omegaTrackIds.end(), lambda.negTrackId) != omegaTrackIds.end()) { - continue; - } + } + for (const auto& lambdas : candidates.lambdas) { + for (const auto& lambda : lambdas) { histos.fill(HIST("QA/massLambda"), std::hypot(lambda.candidate.px, lambda.candidate.py), lambda.candidate.mass); - - for (const auto& kaon : tracks) { - if (!isSelectedKaon(kaon, omega.sign) || - kaon.globalIndex() == lambda.posTrackId || - kaon.globalIndex() == lambda.negTrackId || - std::find(omegaTrackIds.begin(), omegaTrackIds.end(), kaon.globalIndex()) != omegaTrackIds.end()) { - continue; - } - DoubleOmegaCandidate cand; - if (!buildDoubleOmega(collision, omega, lambda, kaon, cand)) { - continue; - } - if constexpr (isMC) { - if (mcParticles == nullptr) { - continue; - } - DoubleOmegaMCInfo mcInfo; - if (getMCInfo(mcInfo, - omegaPosTrack, omegaNegTrack, omegaBachelorTrack, - lambdaPosTrack, lambdaNegTrack, kaon, - omega.sign, *mcParticles)) { - reconstructedDoubleOmegaIds.push_back(mcInfo.motherId); - writeMCCandidate(cand, mcInfo, true); - } - } else { - writeDataCandidate(cand); - } - } } } } + template + CollisionCandidates buildCollisionCandidates(C const& collision, T const& tracks, Cascades const& cascades, V0s const& v0s, bool requireOmega) + { + auto candidates = helper.buildCollisionCandidates(collision, tracks, cascades, v0s, requireOmega); + fillCandidateQA(candidates); + return candidates; + } + template void fillFromCascades(C const& collision, T const&, Cascades const& cascades) { for (const auto& omegaRow : cascades) { BuiltOmegaCandidate omega; - if (!buildOmega(collision, omegaRow.template posTrack_as(), omegaRow.template negTrack_as(), omegaRow.template bachelor_as(), omega)) { + if (!helper.buildOmega(collision, omegaRow.template posTrack_as(), omegaRow.template negTrack_as(), omegaRow.template bachelor_as(), omega)) { continue; } const float omegaPt = std::hypot(omega.px, omega.py); @@ -1037,16 +736,16 @@ struct doubleOmegaTreeCreator { } auto kaon = lambdaKaonSource.template bachelor_as(); - if (!isSelectedKaon(kaon, omega.sign)) { + if (!helper.isSelectedKaon(kaon, omega.sign)) { continue; } BuiltLambdaCandidate lambda; - if (!buildLambda(collision, lambdaKaonSource.template posTrack_as(), lambdaKaonSource.template negTrack_as(), omega.sign, lambda)) { + if (!helper.buildLambda(collision, lambdaKaonSource.template posTrack_as(), lambdaKaonSource.template negTrack_as(), omega.sign, lambda)) { continue; } histos.fill(HIST("QA/massLambda"), std::hypot(lambda.candidate.px, lambda.candidate.py), lambda.candidate.mass); DoubleOmegaCandidate cand; - if (buildDoubleOmega(collision, omega, lambda, kaon, cand)) { + if (helper.buildDoubleOmega({collision.posX(), collision.posY(), collision.posZ()}, omega, lambda, DoubleOmegaBuilderHelper::buildKaon(kaon), cand)) { writeDataCandidate(cand); } } @@ -1078,7 +777,7 @@ struct doubleOmegaTreeCreator { const auto lambdaPosTrack = tracks.rawIteratorAt(daughterTrackIds[isMatter ? 3 : 4]); const auto lambdaNegTrack = tracks.rawIteratorAt(daughterTrackIds[isMatter ? 4 : 3]); const auto kaon = tracks.rawIteratorAt(daughterTrackIds[5]); - if (!isSelectedKaon(kaon, charge)) { + if (!helper.isSelectedKaon(kaon, charge)) { return false; } const int64_t collisionId = selectFindableCollision(kaon); @@ -1090,19 +789,19 @@ struct doubleOmegaTreeCreator { initCCDB(bc); BuiltOmegaCandidate omega; - if (!buildOmega(collision, omegaPosTrack, omegaNegTrack, bachelor, omega)) { + if (!helper.buildOmega(collision, omegaPosTrack, omegaNegTrack, bachelor, omega)) { return false; } const float omegaPt = std::hypot(omega.px, omega.py); histos.fill(HIST("QA/massXi"), omegaPt, omega.massXi); histos.fill(HIST("QA/massOmega"), omegaPt, omega.massOmega); BuiltLambdaCandidate lambda; - if (!buildLambda(collision, lambdaPosTrack, lambdaNegTrack, charge, lambda)) { + if (!helper.buildLambda(collision, lambdaPosTrack, lambdaNegTrack, charge, lambda)) { return false; } histos.fill(HIST("QA/massLambda"), std::hypot(lambda.candidate.px, lambda.candidate.py), lambda.candidate.mass); DoubleOmegaCandidate cand; - if (!buildDoubleOmega(collision, omega, lambda, kaon, cand)) { + if (!helper.buildDoubleOmega({collision.posX(), collision.posY(), collision.posZ()}, omega, lambda, DoubleOmegaBuilderHelper::buildKaon(kaon), cand)) { return false; } logVertexComparison(collision, cand, omega, mcInfo); @@ -1130,15 +829,54 @@ struct doubleOmegaTreeCreator { return true; } + void mix(CollisionCandidates const& reference, CollisionCandidates const& omegaSource, + CollisionCandidates const& lambdaSource, CollisionCandidates const& kaonSource, MixingType mixingType) + { + if (!(mixingOpts.mixingTypes & (1 << (mixingType - 1)))) { + return; + } + helper.combine(reference, omegaSource, lambdaSource, kaonSource, [&](DoubleOmegaCandidate const& cand, auto const&...) { + histos.fill(HIST("Mixing/hMixedCandidates"), mixingType); + writeMixedCandidate(cand, mixingType); + }); + } + + // Two candidates from one collision and one from the other, in both directions. The collision + // providing two candidates is the reference one. + void mixPair(CollisionCandidates const& candidates0, CollisionCandidates const& candidates1) + { + mix(candidates0, candidates1, candidates0, candidates0, kMixedOmega); + mix(candidates1, candidates0, candidates1, candidates1, kMixedOmega); + mix(candidates0, candidates0, candidates1, candidates0, kMixedLambda); + mix(candidates1, candidates1, candidates0, candidates1, kMixedLambda); + mix(candidates0, candidates0, candidates0, candidates1, kMixedKaon); + mix(candidates1, candidates1, candidates1, candidates0, kMixedKaon); + } + + // All the assignments of the Omega, Lambda and kaon to three different collisions. The collision + // providing the Omega is the reference one. + void mixTriple(CollisionCandidates const& candidates0, CollisionCandidates const& candidates1, CollisionCandidates const& candidates2) + { + const std::array sources{&candidates0, &candidates1, &candidates2}; + std::array roles{0, 1, 2}; // Sources of the Omega, Lambda and kaon. + do { + const auto& omegaSource = *sources[roles[0]]; + mix(omegaSource, omegaSource, *sources[roles[1]], *sources[roles[2]], kMixedAll); + } while (std::next_permutation(roles.begin(), roles.end())); + } + void processData(Collisions const& collision, TracksFull const& tracks, aod::V0s const& v0s, aod::Cascades const& cascades, aod::BCsWithTimestamps const& bcs) { - if (acceptCollision(collision, bcs)) { - fillFromRawTables(collision, tracks, cascades, v0s); + if (!acceptCollision(collision, bcs)) { + return; } + const auto candidates = buildCollisionCandidates(collision, tracks, cascades, v0s, true); + helper.combine(candidates, candidates, candidates, candidates, + [&](DoubleOmegaCandidate const& cand, auto const&...) { writeDataCandidate(cand); }); } PROCESS_SWITCH(doubleOmegaTreeCreator, processData, "Reconstruct Omega + Lambda + kaon from raw V0/cascade indices", true); @@ -1174,7 +912,18 @@ struct doubleOmegaTreeCreator { v0sThisCollision.bindExternalIndices(&tracks); cascadesThisCollision.bindExternalIndices(&tracks); cascadesThisCollision.bindExternalIndices(&v0s); - fillFromRawTables(collision, tracksThisCollision, cascadesThisCollision, v0sThisCollision, &mcParticles); + const auto candidates = buildCollisionCandidates(collision, tracksThisCollision, cascadesThisCollision, v0sThisCollision, true); + helper.combine(candidates, candidates, candidates, candidates, + [&](DoubleOmegaCandidate const& cand, BuiltOmegaCandidate const& omega, BuiltLambdaCandidate const& lambda, BuiltKaonCandidate const& kaon) { + DoubleOmegaMCInfo mcInfo; + if (getMCInfo(mcInfo, + tracks.rawIteratorAt(omega.posTrackId), tracks.rawIteratorAt(omega.negTrackId), tracks.rawIteratorAt(omega.bachelorId), + tracks.rawIteratorAt(lambda.posTrackId), tracks.rawIteratorAt(lambda.negTrackId), tracks.rawIteratorAt(kaon.trackId), + omega.sign, mcParticles)) { + reconstructedDoubleOmegaIds.push_back(mcInfo.motherId); + writeMCCandidate(cand, mcInfo, true); + } + }); } for (const auto& mcParticle : mcParticles) { DoubleOmegaMCInfo mcInfo; @@ -1278,7 +1027,7 @@ struct doubleOmegaTreeCreator { if (daughterTrackIds[daughterIndex] < 0) { daughterTrackIds[daughterIndex] = track.globalIndex(); } - if (selectTrack(track)) { + if (helper.selectTrack(track)) { selectedReconstructedDaughters[daughterIndex] = true; if (selectedTrackIds[daughterIndex] < 0) { selectedTrackIds[daughterIndex] = track.globalIndex(); @@ -1384,6 +1133,384 @@ struct doubleOmegaTreeCreator { } } PROCESS_SWITCH(doubleOmegaTreeCreator, processFindableTracks, "Count findable double-Omega decays and reconstruct selected daughter tracks", false); + + // Bin index of value on a ConfigurableAxis (variable or fixed width), -1 outside the axis. + static int findBin(std::vector const& axis, double value) + { + if (axis.size() < 3) { + return -1; + } + if (axis[0] == VARIABLE_WIDTH) { + const auto edge = std::upper_bound(axis.begin() + 1, axis.end(), value); + if (edge == axis.begin() + 1 || edge == axis.end()) { + return -1; + } + return static_cast(std::distance(axis.begin() + 2, edge)); + } + const int nBins = static_cast(axis[0]); + if (value < axis[1] || value >= axis[2]) { + return -1; + } + return static_cast((value - axis[1]) / (axis[2] - axis[1]) * nBins); + } + + static int nBins(std::vector const& axis) + { + return axis[0] == VARIABLE_WIDTH ? static_cast(axis.size()) - 2 : static_cast(axis[0]); + } + + // Daughter track indices identifying a combination. Seeds sharing a daughter would produce the same + // combination more than once: each combination is written only once. + using MixedCombination = std::array; + + // Mixing pool of same-event candidates per (radius, phi, charge) bin, kept across dataframes. + std::unordered_map> mSeedPool; + std::size_t mPoolSeeds = 0; + int64_t mDataframeCounter = 0; + double mProcessingTimeSum = 0.; + + // Indices are only unique within a dataframe: the dataframe counter goes in the upper bits. + int64_t uniqueIndex(int64_t index) const + { + return index < 0 ? index : (mDataframeCounter << 32) | index; + } + + // The track identifiers read from the derived tables (see trackFingerprint) already identify a track + // in any dataframe; the collision index is made unique across dataframes. + MixingSeed makeSeed(int64_t collisionId, std::array const& primaryVertex, BuiltOmegaCandidate const& omega, + BuiltLambdaCandidate const& lambda, BuiltKaonCandidate const& kaon) const + { + return {uniqueIndex(collisionId), primaryVertex, mBz, omega, lambda, kaon}; + } + + // Builds and writes one mixed candidate, with the primary vertex of the reference seed. + void mixSeedDaughters(MixingSeed const& reference, BuiltOmegaCandidate const& omega, BuiltLambdaCandidate const& lambda, + BuiltKaonCandidate const& kaon, MixingType mixingType, std::set& written) + { + if (!(mixingOpts.mixingTypes & (1 << (mixingType - 1)))) { + return; + } + const MixedCombination combination{omega.posTrackId, omega.negTrackId, omega.bachelorId, lambda.posTrackId, lambda.negTrackId, kaon.trackId}; + if (!written.insert(combination).second || DoubleOmegaBuilderHelper::sharesTracks(omega, lambda, kaon)) { + return; + } + DoubleOmegaCandidate cand; + if (!helper.buildDoubleOmega(reference.primaryVertex, omega, lambda, kaon, cand)) { + return; + } + histos.fill(HIST("Mixing/hMixedCandidates"), mixingType); + writeMixedCandidate(cand, mixingType); + } + + // One daughter of each seed replaced by the corresponding daughter of the other seed. The seed + // keeping two of its daughters provides the primary vertex. + void mixSeedPair(MixingSeed const& seed0, MixingSeed const& seed1, std::set& written) + { + mixSeedDaughters(seed0, seed1.omega, seed0.lambda, seed0.kaon, kMixedOmega, written); + mixSeedDaughters(seed1, seed0.omega, seed1.lambda, seed1.kaon, kMixedOmega, written); + mixSeedDaughters(seed0, seed0.omega, seed1.lambda, seed0.kaon, kMixedLambda, written); + mixSeedDaughters(seed1, seed1.omega, seed0.lambda, seed1.kaon, kMixedLambda, written); + mixSeedDaughters(seed0, seed0.omega, seed0.lambda, seed1.kaon, kMixedKaon, written); + mixSeedDaughters(seed1, seed1.omega, seed1.lambda, seed0.kaon, kMixedKaon, written); + } + + // Two daughters of a seed replaced by those of two other seeds: Omega, Lambda and kaon from three + // different collisions, in all the assignments. The seed providing the Omega is the reference. + void mixSeedTriple(MixingSeed const& seed0, MixingSeed const& seed1, MixingSeed const& seed2, std::set& written) + { + const std::array seeds{&seed0, &seed1, &seed2}; + std::array roles{0, 1, 2}; // Seeds providing the Omega, Lambda and kaon. + do { + const auto& omegaSeed = *seeds[roles[0]]; + mixSeedDaughters(omegaSeed, omegaSeed.omega, seeds[roles[1]]->lambda, seeds[roles[2]]->kaon, kMixedAll, written); + } while (std::next_permutation(roles.begin(), roles.end())); + } + + // Seeds can be mixed if they come from different collisions with close primary vertices, recorded + // with the same magnetic field. + bool canMix(MixingSeed const& seed0, MixingSeed const& seed1) const + { + return seed0.collisionId != seed1.collisionId && + std::abs(seed0.bz - seed1.bz) < 1.e-3f && + std::abs(seed0.primaryVertex[2] - seed1.primaryVertex[2]) <= mixingOpts.maxDeltaPVz.value; + } + + // Mixes a new seed with the most recent compatible seeds of its pool bin: pairs with each of them, + // triples with each pair of them coming from different collisions. + void mixWithPool(MixingSeed const& seed, std::deque const& pool, std::set& written) + { + const auto nPartners = static_cast(mixingOpts.nMixingCandidates.value); + std::vector partners; + for (auto partner = pool.rbegin(); partner != pool.rend() && partners.size() < nPartners; ++partner) { + histos.fill(HIST("Mixing/hMixedCollisions"), 0); + if (!canMix(seed, *partner)) { + continue; + } + histos.fill(HIST("Mixing/hMixedCollisions"), 1); + partners.push_back(&*partner); + mixSeedPair(seed, *partner, written); + } + if (!(mixingOpts.mixingTypes & (1 << (kMixedAll - 1)))) { + return; + } + for (std::size_t j = 0; j < partners.size(); ++j) { + for (std::size_t k = j + 1; k < partners.size(); ++k) { + histos.fill(HIST("Mixing/hMixedCollisions"), 2); + if (!canMix(*partners[j], *partners[k])) { + continue; + } + histos.fill(HIST("Mixing/hMixedCollisions"), 3); + mixSeedTriple(seed, *partners[j], *partners[k], written); + } + } + } + + // Same-event candidates of a collision become seeds, binned in decay radius, momentum azimuth and + // charge (mixingOpts.binsRadius, mixingOpts.binsPhi). The seeds are mixed with the pool first and + // added to it afterwards, so seeds of the same collision are never mixed with each other. + void mixCollisionSeeds(int64_t collisionId, CollisionCandidates const& candidates, std::set& written) + { + const int nPhiBins = nBins(mixingOpts.binsPhi.value); + const auto poolDepth = static_cast(mixingOpts.poolDepth.value); + std::vector> newSeeds; + helper.combine(candidates, candidates, candidates, candidates, + [&](DoubleOmegaCandidate const& cand, BuiltOmegaCandidate const& omega, BuiltLambdaCandidate const& lambda, BuiltKaonCandidate const& kaon) { + const float radius = std::hypot(cand.x, cand.y); + const int radiusBin = findBin(mixingOpts.binsRadius.value, radius); + const int phiBin = findBin(mixingOpts.binsPhi.value, cand.phi); + if (radiusBin < 0 || phiBin < 0) { + return; + } + histos.fill(HIST("Mixing/hSeedRadiusPhi"), radius, cand.phi); + const int key = (radiusBin * nPhiBins + phiBin) * 2 + DoubleOmegaBuilderHelper::chargeIndex(cand.sign); + newSeeds.emplace_back(key, makeSeed(collisionId, candidates.primaryVertex, omega, lambda, kaon)); + }); + for (const auto& [key, seed] : newSeeds) { + mixWithPool(seed, mSeedPool[key], written); + } + for (auto& [key, seed] : newSeeds) { + auto& pool = mSeedPool[key]; + pool.push_back(std::move(seed)); + ++mPoolSeeds; + if (pool.size() > poolDepth) { + pool.pop_front(); + --mPoolSeeds; + } + } + } + + void monitorMixing(std::chrono::steady_clock::time_point const& start) + { + const double processingTime = std::chrono::duration(std::chrono::steady_clock::now() - start).count(); + mProcessingTimeSum += processingTime; + // Lower bound: deque and hash-map bookkeeping are not included. + const double poolMemory = static_cast(mPoolSeeds * sizeof(MixingSeed)) / (1024. * 1024.); + histos.fill(HIST("Mixing/hProcessingTime"), processingTime); + histos.fill(HIST("Mixing/hPoolSeeds"), mPoolSeeds); + histos.fill(HIST("Mixing/hPoolMemory"), poolMemory); + histos.fill(HIST("Mixing/hProcessingTimeVsPoolSeeds"), mPoolSeeds, processingTime); + if (mixingOpts.monitoringInterval > 0 && mDataframeCounter % mixingOpts.monitoringInterval == 0) { + LOGP(info, "Candidate mixing, dataframe {}: {} seeds in {} pool bins ({:.2f} MB), {:.1f} ms in this dataframe, {:.1f} ms on average", + mDataframeCounter, mPoolSeeds, mSeedPool.size(), poolMemory, processingTime, mProcessingTimeSum / mDataframeCounter); + } + } + + // --------------------------------------------------------------------------------------------- + // Derived data produced by doubleOmegaBuilder, read on the fly or from saved (merged) derived data. + + // The derived tables store a track once per collision, so the same track can appear in the rows of + // different collisions (V0s and cascades can use tracks associated to another collision). Tracks are + // identified across collisions by their parameters, which are bitwise identical for the same track. + template + static int64_t trackFingerprint(Row const& row) + { + std::size_t hash = 0; + for (const float value : {row.x(), row.alpha(), row.y(), row.z(), row.snp(), row.tgl(), row.signed1Pt()}) { + hash ^= std::hash{}(value) + 0x9e3779b97f4a7c15ULL + (hash << 6) + (hash >> 2); + } + return static_cast(hash); + } + + // Derived collisions store the run number and the magnetic field: the propagator, the material LUT + // and the fitters are configured again only when the field changes. + void initCCDBForDerived(int runNumber, float bz) + { + if (mRunNumber != 0 && std::abs(bz - mBz) < kBzTolerance) { + return; + } + auto* grpmag = ccdb->getForRun("GLO/Config/GRPMagField", runNumber); + if (!grpmag) { + LOG(fatal) << "Could not retrieve GRPMagField for run " << runNumber; + } + o2::base::Propagator::initFieldFromGRP(grpmag); + mBz = o2::base::Propagator::Instance()->getNominalBz(); + helper.setBz(mBz); + if (static_cast(cfgMaterialCorrection.value) == o2::base::Propagator::MatCorrType::USEMatCorrLUT) { + auto* lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->getForRun("GLO/Param/MatLUT", runNumber)); + if (!lut) { + LOG(fatal) << "Could not retrieve material LUT for run " << runNumber; + } + o2::base::Propagator::Instance()->setMatLUT(lut); + } + mRunNumber = runNumber; + } + + // Reads the candidates of a derived collision back, applying again the selections that can be + // tightened with respect to the production ones. + template + CollisionCandidates readCollisionCandidates(C const& collision, aod::DOmegaTracks const& tracks, + aod::DOmegaOmegas const& omegas, aod::DOmegaLambdas const& lambdas) + { + CollisionCandidates candidates; + candidates.primaryVertex = {collision.posX(), collision.posY(), collision.posZ()}; + + for (const auto& row : omegas.sliceBy(derivedOmegasPerCollision, collision.globalIndex())) { + const auto posTrack = tracks.rawIteratorAt(row.posTrackId()); + const auto negTrack = tracks.rawIteratorAt(row.negTrackId()); + const auto bachelor = tracks.rawIteratorAt(row.bachelorId()); + const bool isMatter = row.sign() < 0; + const auto& proton = isMatter ? posTrack : negTrack; + const auto& pion = isMatter ? negTrack : posTrack; + const uint32_t lambdaDetectorMap = static_cast(proton.detectorMap()) | (static_cast(pion.detectorMap()) << 3); + const BuiltOmegaCandidate omega{row.px(), row.py(), row.pz(), + row.decayX(), row.decayY(), row.decayZ(), + row.massOmega(), row.massXi(), + row.sign(), + trackFingerprint(posTrack), trackFingerprint(negTrack), trackFingerprint(bachelor), + getTrackParCov(row), + static_cast(bachelor.detectorMap()) | (lambdaDetectorMap << 3), + bachelor.dcaXYToPV(), proton.dcaXYToPV(), pion.dcaXYToPV(), + row.dcaDaughters()}; + if (helper.isSelectedOmega(omega)) { + candidates.omegas.push_back(omega); + } + } + + for (const auto& row : lambdas.sliceBy(derivedLambdasPerCollision, collision.globalIndex())) { + const auto posTrack = tracks.rawIteratorAt(row.posTrackId()); + const auto negTrack = tracks.rawIteratorAt(row.negTrackId()); + const bool isMatter = row.sign() < 0; + const auto& proton = isMatter ? posTrack : negTrack; + const auto& pion = isMatter ? negTrack : posTrack; + BuiltLambdaCandidate lambda; + lambda.candidate = {.px = row.px(), .py = row.py(), .pz = row.pz(), .mass = row.massLambda()}; + lambda.parentTrack = getTrackParCov(row); + lambda.parentTrack.setAbsCharge(0); // Neutral, as created by the V0 fit. + lambda.decayVertex = {row.decayX(), row.decayY(), row.decayZ()}; + lambda.posTrackId = trackFingerprint(posTrack); + lambda.negTrackId = trackFingerprint(negTrack); + lambda.daughterDetectorMap = static_cast(proton.detectorMap()) | (static_cast(pion.detectorMap()) << 3); + lambda.dcaXYProtonToPV = proton.dcaXYToPV(); + lambda.dcaXYPionToPV = pion.dcaXYToPV(); + lambda.dcaDaughters = row.dcaDaughters(); + if (helper.isSelectedLambda(lambda)) { + candidates.lambdas[DoubleOmegaBuilderHelper::chargeIndex(row.sign())].push_back(lambda); + } + } + + for (const auto& row : tracks.sliceBy(derivedTracksPerCollision, collision.globalIndex())) { + if (!row.isKaon()) { + continue; + } + const BuiltKaonCandidate kaon{getTrackParCov(row), + trackFingerprint(row), + static_cast(row.sign()), + row.tpcNSigmaKa(), + row.tofNSigmaKa(), + row.detectorMap()}; + if (helper.isSelectedKaon(kaon)) { + candidates.kaons[DoubleOmegaBuilderHelper::chargeIndex(kaon.sign)].push_back(kaon); + } + } + fillCandidateQA(candidates); + return candidates; + } + + void processDerivedData(aod::DOmegaColls const& collisions, + aod::DOmegaTracks const& tracks, + aod::DOmegaOmegas const& omegas, + aod::DOmegaLambdas const& lambdas) + { + for (const auto& collision : collisions) { + initCCDBForDerived(collision.runNumber(), collision.bz()); + histos.fill(HIST("QA/zVtx"), collision.posZ()); + const auto candidates = readCollisionCandidates(collision, tracks, omegas, lambdas); + helper.combine(candidates, candidates, candidates, candidates, + [&](DoubleOmegaCandidate const& cand, auto const&...) { writeDataCandidate(cand); }); + } + } + PROCESS_SWITCH(doubleOmegaTreeCreator, processDerivedData, "Same-event Omega + Lambda + kaon from the derived tables of doubleOmegaBuilder", false); + + // Mixed-event background from the derived tables. Saved derived data can be merged into large + // dataframes, which provides many more mixing partners than the input AO2D dataframes. Only + // collisions recorded with the same magnetic field are mixed. + void processDerivedMixing(aod::DOmegaColls const& collisions, + aod::DOmegaTracks const& tracks, + aod::DOmegaOmegas const& omegas, + aod::DOmegaLambdas const& lambdas) + { + std::unordered_map candidatesPerCollision; + for (const auto& collision : collisions) { + candidatesPerCollision.emplace(collision.globalIndex(), readCollisionCandidates(collision, tracks, omegas, lambdas)); + } + + DerivedMixingBinning binning{{mixingOpts.binsVtxZ, mixingOpts.binsMultiplicity}, true}; + for (const auto& [collision0, collision1] : selfPairCombinations(binning, mixingOpts.nMixingCollisions, -1, collisions)) { + histos.fill(HIST("Mixing/hMixedCollisions"), 0); + if (std::abs(collision0.posZ() - collision1.posZ()) > mixingOpts.maxDeltaPVz || + std::abs(collision0.bz() - collision1.bz()) > kBzTolerance) { + continue; + } + histos.fill(HIST("Mixing/hMixedCollisions"), 1); + initCCDBForDerived(collision0.runNumber(), collision0.bz()); + mixPair(candidatesPerCollision.at(collision0.globalIndex()), candidatesPerCollision.at(collision1.globalIndex())); + } + + if (!(mixingOpts.mixingTypes & (1 << (kMixedAll - 1)))) { + return; + } + for (const auto& [collision0, collision1, collision2] : selfTripleCombinations(binning, mixingOpts.nMixingCollisions, -1, collisions)) { + histos.fill(HIST("Mixing/hMixedCollisions"), 2); + const auto [minPVz, maxPVz] = std::minmax({collision0.posZ(), collision1.posZ(), collision2.posZ()}); + if (maxPVz - minPVz > mixingOpts.maxDeltaPVz || + std::abs(collision0.bz() - collision1.bz()) > kBzTolerance || + std::abs(collision0.bz() - collision2.bz()) > kBzTolerance) { + continue; + } + histos.fill(HIST("Mixing/hMixedCollisions"), 3); + initCCDBForDerived(collision0.runNumber(), collision0.bz()); + mixTriple(candidatesPerCollision.at(collision0.globalIndex()), candidatesPerCollision.at(collision1.globalIndex()), + candidatesPerCollision.at(collision2.globalIndex())); + } + } + PROCESS_SWITCH(doubleOmegaTreeCreator, processDerivedMixing, "Mixed-event Omega + Lambda + kaon background from the derived tables of doubleOmegaBuilder", false); + + // Candidate-level mixed-event background from the derived tables of doubleOmegaBuilder, as in + // decay3bodybuilder. The same-event candidates (seeds) are binned in decay radius, momentum azimuth + // and charge, so seeds in the same bin have their decay vertex in the same region of the detector and + // point in the same direction. Mixed candidates are built by swapping one (or two) daughters between + // seeds of different collisions, which keeps the vertex topology without moving any track. Seeds are + // pooled across dataframes. + void processDerivedMixingCandidates(aod::DOmegaColls const& collisions, + aod::DOmegaTracks const& tracks, + aod::DOmegaOmegas const& omegas, + aod::DOmegaLambdas const& lambdas) + { + const auto start = std::chrono::steady_clock::now(); + ++mDataframeCounter; + std::set written; // Every combination has a daughter from this dataframe. + + for (const auto& collision : collisions) { + // Seeds are mixed only with seeds recorded with the same field (canMix), so the field of the + // current collision is the right one for all its mixed candidates. + initCCDBForDerived(collision.runNumber(), collision.bz()); + histos.fill(HIST("QA/zVtx"), collision.posZ()); + const auto candidates = readCollisionCandidates(collision, tracks, omegas, lambdas); + mixCollisionSeeds(collision.globalIndex(), candidates, written); + } + monitorMixing(start); + } + PROCESS_SWITCH(doubleOmegaTreeCreator, processDerivedMixingCandidates, "Candidate-level mixed-event background from the derived tables of doubleOmegaBuilder, seeds binned in decay radius and azimuth", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGLF/Utils/doubleOmegaBuilderHelper.h b/PWGLF/Utils/doubleOmegaBuilderHelper.h new file mode 100644 index 00000000000..ac071658516 --- /dev/null +++ b/PWGLF/Utils/doubleOmegaBuilderHelper.h @@ -0,0 +1,628 @@ +// 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 doubleOmegaBuilderHelper.h +/// \brief Building blocks of the double-Omega (Omega + Lambda + kaon) reconstruction, shared by the +/// candidate builder (derived data producer) and the tree creator. + +#ifndef PWGLF_UTILS_DOUBLEOMEGABUILDERHELPER_H_ +#define PWGLF_UTILS_DOUBLEOMEGABUILDERHELPER_H_ + +#include "Common/Core/RecoDecay.h" +#include "Common/Core/trackUtilities.h" + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace o2::pwglf::doubleomega +{ + +struct DoubleOmegaCandidate { + float pt = -999.f; + float eta = -999.f; + float phi = -999.f; + float x = -999.f; + float y = -999.f; + float z = -999.f; + float pvX = -999.f; + float pvY = -999.f; + float pvZ = -999.f; + float cosPAOmega = -999.f; + float cosPADirectLambda = -999.f; + float cosPADoubleOmega = -999.f; + float cosPAOmegaToPV = -999.f; + float cosPADirectLambdaToPV = -999.f; + float dcaXYOmegaToPV = -999.f; + float dcaXYDirectLambdaToPV = -999.f; + float dcaXYDirectKaonToPV = -999.f; + float dcaXYOmegaKaonToPV = -999.f; + float dcaXYOmegaProtonToPV = -999.f; + float dcaXYOmegaPionToPV = -999.f; + float dcaXYDirectProtonToPV = -999.f; + float dcaXYDirectPionToPV = -999.f; + float tpcNSigmaDirectKaon = -999.f; + float tofNSigmaDirectKaon = -999.f; + uint32_t daughterDetectorMap = 0; + float mass = -999.f; + float massOmega = -999.f; + float massXi = -999.f; + int8_t sign = 0; +}; + +struct LambdaCandidate { + float px = 0.f; + float py = 0.f; + float pz = 0.f; + float mass = -999.f; +}; + +struct BuiltLambdaCandidate { + LambdaCandidate candidate; + o2::track::TrackParCov parentTrack; + std::array decayVertex{}; + int64_t posTrackId = -1; + int64_t negTrackId = -1; + uint32_t daughterDetectorMap = 0; // Proton, pion; three detector bits each. + float dcaXYProtonToPV = -999.f; + float dcaXYPionToPV = -999.f; + float dcaDaughters = -999.f; +}; + +struct BuiltOmegaCandidate { + float px = 0.f; + float py = 0.f; + float pz = 0.f; + float x = 0.f; + float y = 0.f; + float z = 0.f; + float massOmega = -999.f; + float massXi = -999.f; + int8_t sign = 0; + int64_t posTrackId = -1; + int64_t negTrackId = -1; + int64_t bachelorId = -1; + o2::track::TrackParCov parentTrack; + uint32_t daughterDetectorMap = 0; // Kaon, proton, pion; three detector bits each. + float dcaXYKaonToPV = -999.f; + float dcaXYProtonToPV = -999.f; + float dcaXYPionToPV = -999.f; + float dcaDaughters = -999.f; +}; + +struct BuiltKaonCandidate { + o2::track::TrackParCov track; + int64_t trackId = -1; + int8_t sign = 0; + float tpcNSigma = -999.f; // -999 when TPC is unavailable. + float tofNSigma = -999.f; // -999 when TOF is unavailable. + uint32_t detectorMap = 0; +}; + +// Omega, Lambda and kaon candidates of one collision, built once and paired afterwards. +// Lambdas and kaons are stored per charge of the Omega they can be paired with (see chargeIndex). +struct CollisionCandidates { + std::array primaryVertex{}; + std::vector omegas; + std::array, 2> lambdas; + std::array, 2> kaons; +}; + +struct DoubleOmegaSelections { + float etaMax = 0.9f; + float minNITSClus = 0.f; + float minNTPCClus = 80.f; + float minNCrossedRows = 100.f; + float maxNSharedTPCClus = 5.f; + float nSigmaTPCCut = 3.f; + float dcaBachToPV = 0.05f; + float dcaKaonToPV = 0.05f; + float dcaOmegaToPV = 0.f; + float dcaDirectLambdaToPV = 0.f; + float dcaV0DauToPV = 0.05f; + float dcaV0Bach = 1.f; + float dcaLambdaDaughters = 1.f; + float mOmegaWindow = 0.01f; + float mLambdaWindow = 0.01f; + float maxDoubleOmegaMass = 3.6f; + float minCosPAOmega = -1.f; + float minCosPADirectLambda = -1.f; + float minCosPADoubleOmega = -1.f; + float minDoubleOmegaDecayRadius = 1.f; +}; + +class DoubleOmegaBuilderHelper +{ + public: + DoubleOmegaSelections selections; + o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrLUT; + o2::vertexing::DCAFitterN<2> fitter; + o2::vertexing::DCAFitterN<3> fitter3Body; + + void initFitters() + { + configureFitter(fitter); + configureFitter(fitter3Body); + } + + void setBz(float bz) + { + fitter.setBz(bz); + fitter3Body.setBz(bz); + } + + void setMatCorrType(o2::base::Propagator::MatCorrType matCorrType) + { + matCorr = matCorrType; + fitter.setMatCorrType(matCorrType); + fitter3Body.setMatCorrType(matCorrType); + } + + static int chargeIndex(int8_t charge) + { + return charge < 0 ? 0 : 1; + } + + template + static uint32_t detectorMap(Track const& track) + { + return static_cast(track.hasITS()) | + (static_cast(track.hasTPC()) << 1) | + (static_cast(track.hasTOF()) << 2); + } + + template + static BuiltKaonCandidate buildKaon(Track const& track) + { + return {getTrackParCov(track), + track.globalIndex(), + static_cast(track.sign()), + track.hasTPC() ? track.tpcNSigmaKa() : -999.f, + track.hasTOF() ? track.tofNSigmaKa() : -999.f, + detectorMap(track)}; + } + + static bool sharesTracks(BuiltOmegaCandidate const& omega, BuiltLambdaCandidate const& lambda) + { + const std::array omegaTrackIds{omega.posTrackId, omega.negTrackId, omega.bachelorId}; + return std::find(omegaTrackIds.begin(), omegaTrackIds.end(), lambda.posTrackId) != omegaTrackIds.end() || + std::find(omegaTrackIds.begin(), omegaTrackIds.end(), lambda.negTrackId) != omegaTrackIds.end(); + } + + static bool sharesTracks(BuiltOmegaCandidate const& omega, BuiltKaonCandidate const& kaon) + { + return kaon.trackId == omega.posTrackId || kaon.trackId == omega.negTrackId || kaon.trackId == omega.bachelorId; + } + + static bool sharesTracks(BuiltOmegaCandidate const& omega, BuiltLambdaCandidate const& lambda, BuiltKaonCandidate const& kaon) + { + return sharesTracks(omega, lambda) || sharesTracks(omega, kaon) || + kaon.trackId == lambda.posTrackId || kaon.trackId == lambda.negTrackId; + } + + template + bool selectTrack(T const& track) const + { + if (std::abs(track.eta()) > selections.etaMax) { + return false; + } + if (track.itsNCls() < selections.minNITSClus || + track.tpcNClsFound() < selections.minNTPCClus || + track.tpcNClsCrossedRows() < selections.minNCrossedRows || + track.tpcNClsCrossedRows() < 0.8f * track.tpcNClsFindable() || + track.tpcNClsShared() > selections.maxNSharedTPCClus) { + LOG(debug) << "Track failed selection: itsNCls=" << static_cast(track.itsNCls()) << ", tpcNClsFound=" << track.tpcNClsFound() << ", tpcNClsCrossedRows=" << static_cast(track.tpcNClsCrossedRows()) << ", tpcNClsFindable=" << static_cast(track.tpcNClsFindable()) << ", tpcNClsShared=" << static_cast(track.tpcNClsShared()); + return false; + } + return true; + } + + template + bool isSelectedKaon(Track const& track, int8_t charge) const + { + return track.sign() == charge && isSelectedKaon(track); + } + + template + bool isSelectedKaon(Track const& track) const + { + return selectTrack(track) && + std::abs(track.tpcNSigmaKa()) <= selections.nSigmaTPCCut; + } + + // Transverse DCA of a track to the collision primary vertex, computed by propagating + // a copy of the track parameters. Returns false if the propagation fails. + template + bool dcaXYToPV(C const& collision, Track const& track, float& dcaXY) const + { + auto trackParCov = getTrackParCov(track); + const o2::math_utils::Point3D primaryVertex{collision.posX(), collision.posY(), collision.posZ()}; + std::array dca{}; + if (!o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertex, trackParCov, 2.f, matCorr, &dca)) { + return false; + } + dcaXY = dca[0]; + return true; + } + + template + bool buildLambda(C const& collision, Track const& posTrack, Track const& negTrack, int8_t charge, BuiltLambdaCandidate& builtLambda) + { + if (posTrack.sign() <= 0 || negTrack.sign() >= 0 || + !selectTrack(posTrack) || !selectTrack(negTrack)) { + return false; + } + + const bool isMatter = charge < 0; + const float protonNSigma = isMatter ? posTrack.tpcNSigmaPr() : negTrack.tpcNSigmaPr(); + const float pionNSigma = isMatter ? negTrack.tpcNSigmaPi() : posTrack.tpcNSigmaPi(); + if (std::abs(protonNSigma) > selections.nSigmaTPCCut || std::abs(pionNSigma) > selections.nSigmaTPCCut) { + return false; + } + float dcaXYPos = -999.f; + float dcaXYNeg = -999.f; + if (!dcaXYToPV(collision, posTrack, dcaXYPos) || !dcaXYToPV(collision, negTrack, dcaXYNeg)) { + LOG(debug) << "Propagation to the PV failed for Lambda daughters " << posTrack.globalIndex() << ", " << negTrack.globalIndex(); + return false; + } + if (std::abs(dcaXYPos) < selections.dcaV0DauToPV || std::abs(dcaXYNeg) < selections.dcaV0DauToPV) { + LOG(debug) << "Lambda dau don't pass selections: " << "nsigmaPR: " << protonNSigma << ", pionNSigma: " << pionNSigma << ", dcaXYPos: " << dcaXYPos << ", dcaXYNeg: " << dcaXYNeg; + return false; + } + + auto posTrackParCov = getTrackParCov(posTrack); + auto negTrackParCov = getTrackParCov(negTrack); + int nCandidates = 0; + try { + nCandidates = fitter.process(posTrackParCov, negTrackParCov); + } catch (...) { + LOG(error) << "Exception while fitting Lambda daughter tracks " << posTrack.globalIndex() << ", " << negTrack.globalIndex(); + return false; + } + if (nCandidates == 0) { + return false; + } + + std::array posMomentum{}; + std::array negMomentum{}; + fitter.getTrack(0).getPxPyPzGlo(posMomentum); + fitter.getTrack(1).getPxPyPzGlo(negMomentum); + const std::array lambdaMomentum{ + posMomentum[0] + negMomentum[0], + posMomentum[1] + negMomentum[1], + posMomentum[2] + negMomentum[2]}; + const auto& fittedVertex = fitter.getPCACandidate(); + const std::array decayVertex{ + static_cast(fittedVertex[0]), static_cast(fittedVertex[1]), static_cast(fittedVertex[2])}; + const float mass = isMatter ? RecoDecay::m(std::array{posMomentum, negMomentum}, + std::array{o2::constants::physics::MassProton, o2::constants::physics::MassPionCharged}) + : RecoDecay::m(std::array{posMomentum, negMomentum}, + std::array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassProton}); + const float dcaDaughters = std::sqrt(std::abs(fitter.getChi2AtPCACandidate())); + const float eta = RecoDecay::eta(lambdaMomentum); + if (std::abs(eta) > selections.etaMax || + dcaDaughters > selections.dcaLambdaDaughters || + std::abs(mass - o2::constants::physics::MassLambda0) > selections.mLambdaWindow) { + return false; + } + + builtLambda.candidate = { + .px = lambdaMomentum[0], + .py = lambdaMomentum[1], + .pz = lambdaMomentum[2], + .mass = mass}; + builtLambda.parentTrack = fitter.createParentTrackParCov(0); + builtLambda.decayVertex = decayVertex; + builtLambda.posTrackId = posTrack.globalIndex(); + builtLambda.negTrackId = negTrack.globalIndex(); + // Keep particle-role ordering for matter and antimatter alike. + builtLambda.daughterDetectorMap = isMatter ? detectorMap(posTrack) | (detectorMap(negTrack) << 3) + : detectorMap(negTrack) | (detectorMap(posTrack) << 3); + builtLambda.dcaXYProtonToPV = isMatter ? dcaXYPos : dcaXYNeg; + builtLambda.dcaXYPionToPV = isMatter ? dcaXYNeg : dcaXYPos; + builtLambda.dcaDaughters = dcaDaughters; + return true; + } + + template + bool buildOmega(C const& collision, Track const& posTrack, Track const& negTrack, Track const& bachelor, BuiltOmegaCandidate& builtOmega) + { + if (bachelor.sign() == 0 || + !selectTrack(bachelor) || + std::abs(bachelor.tpcNSigmaKa()) > selections.nSigmaTPCCut) { + return false; + } + float dcaXYBachelor = -999.f; + if (!dcaXYToPV(collision, bachelor, dcaXYBachelor) || std::abs(dcaXYBachelor) < selections.dcaBachToPV) { + return false; + } + + BuiltLambdaCandidate lambda; + if (!buildLambda(collision, posTrack, negTrack, bachelor.sign(), lambda)) { + return false; + } + + auto v0TrackParCov = lambda.parentTrack; + auto bachelorTrackParCov = getTrackParCov(bachelor); + int nCandidates = 0; + try { + nCandidates = fitter.process(v0TrackParCov, bachelorTrackParCov); + } catch (...) { + LOG(error) << "Exception while fitting Omega daughter tracks " << posTrack.globalIndex() << ", " << negTrack.globalIndex() << ", " << bachelor.globalIndex(); + return false; + } + if (nCandidates == 0) { + return false; + } + + std::array lambdaMomentum{}; + std::array bachelorMomentum{}; + fitter.getTrack(0).getPxPyPzGlo(lambdaMomentum); + fitter.getTrack(1).getPxPyPzGlo(bachelorMomentum); + const std::array omegaMomentum{ + lambdaMomentum[0] + bachelorMomentum[0], + lambdaMomentum[1] + bachelorMomentum[1], + lambdaMomentum[2] + bachelorMomentum[2]}; + const auto& fittedVertex = fitter.getPCACandidate(); + const std::array decayVertex{ + static_cast(fittedVertex[0]), static_cast(fittedVertex[1]), static_cast(fittedVertex[2])}; + const float massOmega = RecoDecay::m(std::array{lambdaMomentum, bachelorMomentum}, + std::array{o2::constants::physics::MassLambda0, o2::constants::physics::MassKaonCharged}); + const float massXi = RecoDecay::m(std::array{lambdaMomentum, bachelorMomentum}, + std::array{o2::constants::physics::MassLambda0, o2::constants::physics::MassPionCharged}); + const float dcaDaughters = std::sqrt(std::abs(fitter.getChi2AtPCACandidate())); + const float eta = RecoDecay::eta(omegaMomentum); + if (dcaDaughters > selections.dcaV0Bach || + std::abs(eta) > selections.etaMax || + std::abs(massOmega - o2::constants::physics::MassOmegaMinus) > selections.mOmegaWindow) { + return false; + } + + builtOmega = { + omegaMomentum[0], + omegaMomentum[1], + omegaMomentum[2], + decayVertex[0], + decayVertex[1], + decayVertex[2], + massOmega, + massXi, + static_cast(bachelor.sign()), + posTrack.globalIndex(), + negTrack.globalIndex(), + bachelor.globalIndex(), + fitter.createParentTrackParCov(), + detectorMap(bachelor) | (lambda.daughterDetectorMap << 3), + dcaXYBachelor, + lambda.dcaXYProtonToPV, + lambda.dcaXYPionToPV, + dcaDaughters}; + return true; + } + + // Selections of the Omega, Lambda and kaon candidates that can be applied again on their stored + // properties (e.g. when reading them back from derived data), to tighten the production selections. + bool isSelectedOmega(BuiltOmegaCandidate const& omega) const + { + return std::abs(RecoDecay::eta(std::array{omega.px, omega.py, omega.pz})) <= selections.etaMax && + std::abs(omega.massOmega - o2::constants::physics::MassOmegaMinus) <= selections.mOmegaWindow && + omega.dcaDaughters <= selections.dcaV0Bach && + std::abs(omega.dcaXYKaonToPV) >= selections.dcaBachToPV && + std::abs(omega.dcaXYProtonToPV) >= selections.dcaV0DauToPV && + std::abs(omega.dcaXYPionToPV) >= selections.dcaV0DauToPV; + } + + bool isSelectedLambda(BuiltLambdaCandidate const& lambda) const + { + return std::abs(RecoDecay::eta(std::array{lambda.candidate.px, lambda.candidate.py, lambda.candidate.pz})) <= selections.etaMax && + std::abs(lambda.candidate.mass - o2::constants::physics::MassLambda0) <= selections.mLambdaWindow && + lambda.dcaDaughters <= selections.dcaLambdaDaughters && + std::abs(lambda.dcaXYProtonToPV) >= selections.dcaV0DauToPV && + std::abs(lambda.dcaXYPionToPV) >= selections.dcaV0DauToPV; + } + + bool isSelectedKaon(BuiltKaonCandidate const& kaon) const + { + return kaon.tpcNSigma > -999.f && std::abs(kaon.tpcNSigma) <= selections.nSigmaTPCCut; + } + + bool buildDoubleOmega(std::array const& primaryVertex, BuiltOmegaCandidate const& omega, BuiltLambdaCandidate const& directLambda, BuiltKaonCandidate const& kaon, DoubleOmegaCandidate& cand) + { + auto omegaTrackParCov = omega.parentTrack; + auto lambdaTrackParCov = directLambda.parentTrack; + auto kaonTrackParCov = kaon.track; + int nCandidates = 0; + try { + nCandidates = fitter3Body.process(omegaTrackParCov, lambdaTrackParCov, kaonTrackParCov); + } catch (...) { + LOG(error) << "Exception while fitting the double-Omega candidate"; + return false; + } + if (nCandidates == 0) { + return false; + } + + auto omegaTrackAtVertex = fitter3Body.getTrack(0); + auto lambdaTrackAtVertex = fitter3Body.getTrack(1); + auto kaonTrackAtVertex = fitter3Body.getTrack(2); + std::array omegaMomentum{}; + std::array lambdaMomentum{}; + std::array kaonMomentum{}; + omegaTrackAtVertex.getPxPyPzGlo(omegaMomentum); + lambdaTrackAtVertex.getPxPyPzGlo(lambdaMomentum); + kaonTrackAtVertex.getPxPyPzGlo(kaonMomentum); + + const std::array totalMomentum{ + omegaMomentum[0] + lambdaMomentum[0] + kaonMomentum[0], + omegaMomentum[1] + lambdaMomentum[1] + kaonMomentum[1], + omegaMomentum[2] + lambdaMomentum[2] + kaonMomentum[2]}; + const auto& fittedVertex = fitter3Body.getPCACandidate(); + const std::array decayVertex{ + static_cast(fittedVertex[0]), static_cast(fittedVertex[1]), static_cast(fittedVertex[2])}; + const o2::math_utils::Point3D primaryVertexPoint{primaryVertex[0], primaryVertex[1], primaryVertex[2]}; + const std::array omegaDecayVertex{omega.x, omega.y, omega.z}; + + std::array dcaOmega{}; + std::array dcaDirectLambda{}; + std::array dcaDirectKaon{}; + if (!o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertexPoint, omegaTrackAtVertex, 2.f, matCorr, &dcaOmega) || + !o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertexPoint, lambdaTrackAtVertex, 2.f, matCorr, &dcaDirectLambda) || + !o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertexPoint, kaonTrackAtVertex, 2.f, matCorr, &dcaDirectKaon)) { + return false; + } + + const float massSquared = RecoDecay::m2(std::array{omegaMomentum, lambdaMomentum, kaonMomentum}, + std::array{o2::constants::physics::MassOmegaMinus, o2::constants::physics::MassLambda0, o2::constants::physics::MassKaonCharged}); + + cand.pt = std::hypot(totalMomentum[0], totalMomentum[1]); + cand.eta = RecoDecay::eta(totalMomentum); + cand.phi = std::atan2(totalMomentum[1], totalMomentum[0]); + cand.x = decayVertex[0]; + cand.y = decayVertex[1]; + cand.z = decayVertex[2]; + cand.pvX = primaryVertex[0]; + cand.pvY = primaryVertex[1]; + cand.pvZ = primaryVertex[2]; + cand.cosPAOmega = RecoDecay::cpa(decayVertex, omegaDecayVertex, omegaMomentum); + cand.cosPADirectLambda = RecoDecay::cpa(decayVertex, directLambda.decayVertex, lambdaMomentum); + cand.cosPADoubleOmega = RecoDecay::cpa(primaryVertex, decayVertex, totalMomentum); + // PV pointing angles use the daughter momenta at their own decay vertices. + cand.cosPAOmegaToPV = RecoDecay::cpa(primaryVertex, omegaDecayVertex, std::array{omega.px, omega.py, omega.pz}); + cand.cosPADirectLambdaToPV = RecoDecay::cpa(primaryVertex, directLambda.decayVertex, std::array{directLambda.candidate.px, directLambda.candidate.py, directLambda.candidate.pz}); + cand.dcaXYOmegaToPV = dcaOmega[0]; + cand.dcaXYDirectLambdaToPV = dcaDirectLambda[0]; + cand.dcaXYDirectKaonToPV = dcaDirectKaon[0]; + cand.dcaXYOmegaKaonToPV = omega.dcaXYKaonToPV; + cand.dcaXYOmegaProtonToPV = omega.dcaXYProtonToPV; + cand.dcaXYOmegaPionToPV = omega.dcaXYPionToPV; + cand.dcaXYDirectProtonToPV = directLambda.dcaXYProtonToPV; + cand.dcaXYDirectPionToPV = directLambda.dcaXYPionToPV; + cand.tpcNSigmaDirectKaon = kaon.tpcNSigma; + cand.tofNSigmaDirectKaon = kaon.tofNSigma; + cand.daughterDetectorMap = omega.daughterDetectorMap | + (directLambda.daughterDetectorMap << 9) | + (kaon.detectorMap << 15); + cand.mass = std::sqrt(std::max(0.f, massSquared)); + if (cand.cosPAOmega < selections.minCosPAOmega || + cand.cosPADirectLambda < selections.minCosPADirectLambda || + cand.cosPADoubleOmega < selections.minCosPADoubleOmega || + std::hypot(decayVertex[0], decayVertex[1]) < selections.minDoubleOmegaDecayRadius || + std::abs(cand.dcaXYOmegaToPV) < selections.dcaOmegaToPV || + std::abs(cand.dcaXYDirectLambdaToPV) < selections.dcaDirectLambdaToPV || + std::abs(cand.dcaXYDirectKaonToPV) < selections.dcaKaonToPV || + cand.mass > selections.maxDoubleOmegaMass) { + return false; + } + cand.massOmega = omega.massOmega; + cand.massXi = omega.massXi; + cand.sign = omega.sign; + return true; + } + + // Builds the Omega, Lambda and kaon candidates of a collision. If requireOmega is set, Lambdas and + // kaons are built only if the collision has at least one Omega. + template + CollisionCandidates buildCollisionCandidates(C const& collision, T const& tracks, Cascades const& cascades, V0s const& v0s, bool requireOmega) + { + CollisionCandidates candidates; + candidates.primaryVertex = {collision.posX(), collision.posY(), collision.posZ()}; + for (const auto& rawCascade : cascades) { + const auto omegaV0 = rawCascade.template v0_as(); + BuiltOmegaCandidate omega; + if (buildOmega(collision, omegaV0.template posTrack_as(), omegaV0.template negTrack_as(), rawCascade.template bachelor_as(), omega)) { + candidates.omegas.push_back(omega); + } + } + if (requireOmega && candidates.omegas.empty()) { + return candidates; + } + + for (const int8_t charge : {int8_t{-1}, int8_t{1}}) { + for (const auto& rawV0 : v0s) { + BuiltLambdaCandidate lambda; + if (buildLambda(collision, rawV0.template posTrack_as(), rawV0.template negTrack_as(), charge, lambda)) { + candidates.lambdas[chargeIndex(charge)].push_back(lambda); + } + } + for (const auto& track : tracks) { + if (isSelectedKaon(track, charge)) { + candidates.kaons[chargeIndex(charge)].push_back(buildKaon(track)); + } + } + } + return candidates; + } + + // Pairs the Omegas, Lambdas and kaons of (possibly different) collisions. The double-Omega candidate + // is built with respect to the primary vertex of the reference collision. + template + void combine(CollisionCandidates const& reference, CollisionCandidates const& omegaSource, + CollisionCandidates const& lambdaSource, CollisionCandidates const& kaonSource, F&& onCandidate) + { + for (const auto& omega : omegaSource.omegas) { + for (const auto& lambda : lambdaSource.lambdas[chargeIndex(omega.sign)]) { + for (const auto& kaon : kaonSource.kaons[chargeIndex(omega.sign)]) { + // Tracks can also be shared by V0s and cascades associated to different collisions. + if (sharesTracks(omega, lambda, kaon)) { + continue; + } + DoubleOmegaCandidate cand; + if (!buildDoubleOmega(reference.primaryVertex, omega, lambda, kaon, cand)) { + continue; + } + onCandidate(cand, omega, lambda, kaon); + } + } + } + } + + // True if at least one Omega has a Lambda and a kaon of compatible charge that do not share its tracks. + static bool hasOmegaLambdaKaon(CollisionCandidates const& candidates) + { + for (const auto& omega : candidates.omegas) { + const auto& lambdas = candidates.lambdas[chargeIndex(omega.sign)]; + const auto& kaons = candidates.kaons[chargeIndex(omega.sign)]; + const bool hasLambda = std::any_of(lambdas.begin(), lambdas.end(), [&](BuiltLambdaCandidate const& lambda) { return !sharesTracks(omega, lambda); }); + const bool hasKaon = std::any_of(kaons.begin(), kaons.end(), [&](BuiltKaonCandidate const& kaon) { return !sharesTracks(omega, kaon); }); + if (hasLambda && hasKaon) { + return true; + } + } + return false; + } + + private: + template + static void configureFitter(Fitter& dcaFitter) + { + dcaFitter.setPropagateToPCA(true); + dcaFitter.setMaxR(200.); + dcaFitter.setMaxDZIni(4.); + dcaFitter.setMinParamChange(1.e-3); + dcaFitter.setMinRelChi2Change(0.9); + dcaFitter.setUseAbsDCA(true); + } +}; + +} // namespace o2::pwglf::doubleomega + +#endif // PWGLF_UTILS_DOUBLEOMEGABUILDERHELPER_H_