diff --git a/PWGHF/D2H/Tasks/taskCd.cxx b/PWGHF/D2H/Tasks/taskCd.cxx index 8282a7eecf5..9e46327cccc 100644 --- a/PWGHF/D2H/Tasks/taskCd.cxx +++ b/PWGHF/D2H/Tasks/taskCd.cxx @@ -257,6 +257,12 @@ struct HfTaskCd { using CollisionsWithEvSelFT0M = soa::Join; using CollisionsMcWithEvSelFT0M = soa::Join; + struct : ConfigurableGroup { + std::string prefix = "mixing"; + Configurable enabled{"enabled", false, "Analyze mixed-event input in a separate run"}; + Configurable type{"type", 0, "QA label matching builder mixing.type: 0 d, 1 K, 2 pion"}; + } mixing; + using CdCandidates = soa::Filtered>; using CdCandidatesMc = soa::Filtered>; using McParticles3ProngMatched = soa::Join; @@ -327,6 +333,30 @@ struct HfTaskCd { } const bool isData = doprocessDataStd || doprocessDataStdWithFT0C || doprocessDataStdWithFT0M; + if (mixing.enabled) { + constexpr int LastMixingType{2}; + if (!isData || mixing.type < 0 || mixing.type > LastMixingType) { + LOGP(fatal, "Mixed-event taskCd requires a data process and a valid mixing.type."); + } + registry.add("Mixed/hMassVsPtVsType", "Mixed candidates after selector;M(dK#pi);p_{T};mixing type", HistType::kTH3F, + {{400, 2.4, 4.4}, {binsPt, "p_{T}"}, {3, -0.5, 2.5}}); + } + if (isData && (fillCandLiteTree || fillCandFullTree)) { + registry.add("Data/hTreeCutFlow", "Cumulative write preselection per mass hypothesis;Stage;Hypothesis", HistType::kTH2D, {{4, -0.5, 3.5}, {2, -0.5, 1.5}}); + registry.add("Data/hDeuteronTofStatus", "TOF status before write preselection;Status;Hypothesis", HistType::kTH2D, {{5, -0.5, 4.5}, {2, -0.5, 1.5}}); + const std::array cutLabels{"Before tree cuts", "After TOF (or disabled)", "After minimum d DCA", "After DCA ordering (or disabled)"}; + const std::array tofLabels{"No TOF match", "Nonfinite nSigma", "Missing nSigma sentinel", "Finite nSigma within cut", "Finite nSigma outside cut"}; + for (size_t i = 0; i < cutLabels.size(); ++i) { + registry.get(HIST("Data/hTreeCutFlow"))->GetXaxis()->SetBinLabel(i + 1, cutLabels[i].c_str()); + } + for (size_t i = 0; i < tofLabels.size(); ++i) { + registry.get(HIST("Data/hDeuteronTofStatus"))->GetXaxis()->SetBinLabel(i + 1, tofLabels[i].c_str()); + } + for (const auto& histogram : {registry.get(HIST("Data/hTreeCutFlow")), registry.get(HIST("Data/hDeuteronTofStatus"))}) { + histogram->GetYaxis()->SetBinLabel(1, "DeKPi"); + histogram->GetYaxis()->SetBinLabel(2, "PiKDe"); + } + } auto addHistogramsRec = [&](const std::string& histoName, const std::string& xAxisTitle, const std::string& yAxisTitle, const HistogramConfigSpec& configSpec) { if (!isData) { registry.add(("MC/reconstructed/signal/" + histoName + "RecSig").c_str(), ("3-prong candidates (matched);" + xAxisTitle + ";" + yAxisTitle).c_str(), configSpec); @@ -768,6 +798,14 @@ struct HfTaskCd { const auto chi2PCA = candidate.chi2PCA(); const auto cpa = candidate.cpa(); const auto cpaXY = candidate.cpaXY(); + if (mixing.enabled) { + if (candidate.isSelCdToDeKPi() >= selectionFlagCd) { + registry.fill(HIST("Mixed/hMassVsPtVsType"), HfHelper::invMassCdToDeKPi(candidate), candidate.pt(), mixing.type.value); + } + if (candidate.isSelCdToPiKDe() >= selectionFlagCd) { + registry.fill(HIST("Mixed/hMassVsPtVsType"), HfHelper::invMassCdToPiKDe(candidate), candidate.pt(), mixing.type.value); + } + } if (candidate.isSelCdToDeKPi() >= selectionFlagCd) { registry.fill(HIST("Data/hMass"), HfHelper::invMassCdToDeKPi(candidate)); registry.fill(HIST("Data/hMassVsPtVsNPvContributors"), HfHelper::invMassCdToDeKPi(candidate), pt, numPvContributors); @@ -874,16 +912,34 @@ struct HfTaskCd { registry.fill(HIST("Data/hNsigmaTPCKaVsP"), prong1.tpcInnerParam() * prong1.sign(), nSigmaTpcKa); registry.fill(HIST("Data/hNsigmaTOFKaVsP"), prong1.tpcInnerParam() * prong1.sign(), nSigmaTofKa); + // Diagnostic categories only: preserve the existing selection predicates below. + enum TofStatus { NoMatch, + Nonfinite, + Missing, + WithinCut, + OutsideCut }; + constexpr float MissingTofNSigma{-999.f}; + const int hypothesis = isDeKPi ? 0 : 1; + const int tofStatus = !deuteronProng.hasTOF() ? NoMatch : !std::isfinite(nSigmaTofDe) ? Nonfinite + : nSigmaTofDe <= MissingTofNSigma ? Missing + : std::abs(nSigmaTofDe) <= cfgMaxDeuteronTofPidPreselection ? WithinCut + : OutsideCut; + registry.fill(HIST("Data/hDeuteronTofStatus"), tofStatus, hypothesis); + registry.fill(HIST("Data/hTreeCutFlow"), 0, hypothesis); if (cfgUseTofPidForDeuteron && std::abs(nSigmaTofDe) > cfgMaxDeuteronTofPidPreselection) { return; } + registry.fill(HIST("Data/hTreeCutFlow"), 1, hypothesis); if (std::abs(dcaDeuteron) < cfgMinDeuteronDcaPreselection) { return; } + registry.fill(HIST("Data/hTreeCutFlow"), 2, hypothesis); if (cfgCutOnDeuteronDcaOrdering && (std::abs(dcaDeuteron) > std::abs(dcaKaon) || std::abs(dcaDeuteron) > std::abs(dcaPion))) { return; } + registry.fill(HIST("Data/hTreeCutFlow"), 3, hypothesis); + if (fillCandLiteTree) { rowCandCdLite( invMassCd, invMassLc, pt, eta, phi, ptProng0, ptProng1, ptProng2, diff --git a/PWGHF/TableProducer/candidateCreator3Prong.cxx b/PWGHF/TableProducer/candidateCreator3Prong.cxx index 1091c00a0a2..e29f4f77f5b 100644 --- a/PWGHF/TableProducer/candidateCreator3Prong.cxx +++ b/PWGHF/TableProducer/candidateCreator3Prong.cxx @@ -69,14 +69,21 @@ #include +#include #include +#include #include #include #include +#include +#include #include #include +#include #include #include +#include +#include #include #include @@ -140,6 +147,16 @@ struct HfCandidateCreator3Prong { Configurable createLc{"createLc", false, "enable Lc+/- candidate creation"}; Configurable createXic{"createXic", false, "enable Xic+/- candidate creation"}; Configurable createCharmNuclei{"createCharmNuclei", false, "enable createCharmNuclei candidate creation"}; + struct : ConfigurableGroup { + std::string prefix = "mixing"; + Configurable enabled{"enabled", false, "Enable Cd event mixing in the existing no-PV-refit DCA data process"}; + Configurable type{"type", 0, "Replace 0: deuteron, 1: kaon, 2: pion"}; + Configurable depth{"depth", 5, "Previous distinct events per pool within a dataframe"}; + Configurable zBinWidth{"zBinWidth", 1.f, "PV z pool width (cm); no track translation"}; + Configurable pvMultBinWidth{"pvMultBinWidth", 20, "PV track multiplicity pool width; used only without centrality"}; + Configurable centralityBinWidth{"centralityBinWidth", 10.f, "Centrality pool width (%); replaces PV multiplicity binning in FT0C/FT0M modes"}; + } mixing; + // KF Configurable applyTopoConstraint{"applyTopoConstraint", false, "apply origin from PV hypothesis for created candidate, works only in KF mode"}; Configurable applyInvMassConstraint{"applyInvMassConstraint", false, "apply particle type hypothesis to recalculate created candidate's momentum, works only in KF mode"}; @@ -178,6 +195,31 @@ struct HfCandidateCreator3Prong { if ((std::accumulate(doprocessDF.begin(), doprocessDF.end(), 0) + std::accumulate(doprocessKF.begin(), doprocessKF.end(), 0)) != 1) { LOGP(fatal, "One and only one process function must be enabled at a time."); } + if (mixing.enabled) { + constexpr int LastMixingType{2}; + if (!(doprocessNoPvRefitWithDCAFitterN || doprocessNoPvRefitWithDCAFitterNCentFT0C || doprocessNoPvRefitWithDCAFitterNCentFT0M) || + !createCharmNuclei || createDplus || createDs || createLc || createXic) { + LOGP(fatal, "Cd mixing requires a no-PV-refit DCA data process and only createCharmNuclei enabled."); + } + if (mixing.type < 0 || mixing.type > LastMixingType || mixing.depth <= 0 || + !std::isfinite(mixing.zBinWidth.value) || mixing.zBinWidth <= 0.f || (doprocessNoPvRefitWithDCAFitterN && mixing.pvMultBinWidth <= 0) || + !std::isfinite(mixing.centralityBinWidth.value) || mixing.centralityBinWidth <= 0.f) { + LOGP(fatal, "Invalid mixing configuration."); + } + registry.add("Mixing/hCounter", "Cd mixing;stage;Entries", HistType::kTH1D, {{5, 0., 5.}}); + registry.get(HIST("Mixing/hCounter"))->GetXaxis()->SetBinLabel(1, "seed candidates"); + registry.get(HIST("Mixing/hCounter"))->GetXaxis()->SetBinLabel(2, "accepted events"); + registry.get(HIST("Mixing/hCounter"))->GetXaxis()->SetBinLabel(3, "candidate pairs"); + registry.get(HIST("Mixing/hCounter"))->GetXaxis()->SetBinLabel(4, "unique fit attempts"); + registry.get(HIST("Mixing/hCounter"))->GetXaxis()->SetBinLabel(5, "written candidates"); + const AxisSpec partnerAxis{mixing.depth.value + 1, -0.5, mixing.depth.value + 0.5}; + registry.add("Mixing/hPoolSize", "Buffered events before insertion;Pool size;Event-charge entries", HistType::kTH1D, {partnerAxis}); + registry.add("Mixing/hNPartners", "Partners with at least one usable candidate pair;Partners;Event-charge entries", HistType::kTH1D, {partnerAxis}); + registry.add("Mixing/hEventPairDeltaPVZ", "One entry per usable event pair and charge;z_{current}-z_{previous} (cm);Event pairs", HistType::kTH1D, {{200, -10., 10.}}); + registry.add("Mixing/hEventPairCentrality", "One entry per usable event pair and charge;Current centrality (%);Previous centrality (%)", HistType::kTH2D, {{102, -1.5, 100.5}, {102, -1.5, 100.5}}); + registry.add("Mixing/hPoolOccupancy", "Occupancy at end of dataframe;Buffered events;Pools", HistType::kTH1D, {partnerAxis}); + registry.add("Mixing/hDeltaPVZ", "Mixed event PV difference;z_{A}-z_{B} (cm);Candidate pairs", HistType::kTH1F, {{200, -10., 10.}}); + } std::array processesCollisions = {doprocessCollisions, doprocessCollisionsCentFT0C, doprocessCollisionsCentFT0M, doprocessCollisionsUpc}; const int nProcessesCollisions = std::accumulate(processesCollisions.begin(), processesCollisions.end(), 0); @@ -715,6 +757,180 @@ struct HfCandidateCreator3Prong { /// /// /////////////////////////////////// + // In-memory adapter only: no new AOD table and no changes to the ordinary fitter. + template + struct HfMixed3ProngSeed { + TSeed source; + TTracks const* tracks{nullptr}; + std::array ids{}; + uint8_t flag{0}; + template + auto collision_as() const + { + return source.template collision_as(); + } + template + auto prong0_as() const + { + return tracks->rawIteratorAt(ids[0]); + } + template + auto prong1_as() const + { + return tracks->rawIteratorAt(ids[1]); + } + template + auto prong2_as() const + { + return tracks->rawIteratorAt(ids[2]); + } + int64_t prong0Id() const { return ids[0]; } + int64_t prong1Id() const { return ids[1]; } + int64_t prong2Id() const { return ids[2]; } + uint8_t hfflag() const { return flag; } + }; + + template + void runCreator3ProngMixedWithDCAFitterN(TCollisions const& collisions, TCandidates const& candidates, TTracks const& tracks, TBCs const& bcs, uint8_t channelFlag) + { + constexpr int LastProng{2}; + using Seed = std::decay_t; + using Event = std::vector; + using PoolKey = std::tuple; // run, charge, PV z, PV multiplicity OR centrality (unused bin is zero) + std::map events; + std::map> pools; + std::set> seedKeys; + std::map, Seed> mixedCandidates; // reference collision and ordered prong track IDs + + // These are skim candidates, before final topology/PID/BDT selection. + // Enumerate both d hypotheses for each source below, including cross-prong exchanges. + for (const auto& seed : candidates) { + if (!(seed.hfflag() & channelFlag)) { + continue; + } + std::array key{seed.collisionId(), seed.prong0Id(), seed.prong1Id(), seed.prong2Id()}; + if (seedKeys.insert(key).second) { + events[seed.collisionId()].push_back(seed); + registry.fill(HIST("Mixing/hCounter"), 0.5); + } + } + + for (const auto& [collisionId, event] : events) { + auto collision = event.front().template collision_as(); + float centrality{-1.f}; + if (hfEvSel.getHfCollisionRejectionMask(collision, centrality, ccdb, registry) != 0) { + continue; + } + if (!std::isfinite(collision.posZ())) { + continue; + } + registry.fill(HIST("Mixing/hCounter"), 1.5); + const int eventRun = collision.template bc_as().runNumber(); + const int zBin = static_cast(std::floor(collision.posZ() / mixing.zBinWidth)); + int multBin = 0; + if constexpr (CentEstimator == CentralityEstimator::None) { + multBin = collision.multNTracksPV() / mixing.pvMultBinWidth; + } + int centralityBin = 0; + if constexpr (CentEstimator != CentralityEstimator::None) { + if (!std::isfinite(centrality) || centrality < 0.f) { + continue; + } + centralityBin = static_cast(std::floor(centrality / mixing.centralityBinWidth)); + } + + // Separate matter/antimatter within each accepted event before adding it to pools. + for (const int& charge : {-1, 1}) { + Event seeds; + for (const auto& seed : event) { + auto t0 = seed.template prong0_as(); + auto t1 = seed.template prong1_as(); + auto t2 = seed.template prong2_as(); + if (t0.sign() == charge && t1.sign() == -charge && t2.sign() == charge && + t0.globalIndex() != t2.globalIndex()) { + seeds.push_back(seed); + } + } + if (seeds.empty()) { + continue; + } + auto& pool = pools[PoolKey{eventRun, charge, zBin, multBin, centralityBin}]; + registry.fill(HIST("Mixing/hPoolSize"), pool.size()); + int nPartners = 0; + for (const auto& previousEvent : pool) { + bool hasUsablePair = false; + for (const auto& seedA : seeds) { + for (const auto& seedB : previousEvent) { + if (seedA.collisionId() == seedB.collisionId()) { + continue; + } + std::array idsA{seedA.prong0Id(), seedA.prong1Id(), seedA.prong2Id()}; + std::array idsB{seedB.prong0Id(), seedB.prong1Id(), seedB.prong2Id()}; + bool sharesTrack = false; + for (const auto& id : idsA) { + sharesTrack |= std::find(idsB.begin(), idsB.end(), id) != idsB.end(); + } + if (sharesTrack) { + continue; + } + hasUsablePair = true; + registry.fill(HIST("Mixing/hCounter"), 2.5); + auto collisionB = seedB.template collision_as(); + registry.fill(HIST("Mixing/hDeltaPVZ"), collision.posZ() - collisionB.posZ()); + + auto buildMixed = [&](auto const& referenceCollision, auto const& referenceSeed, + std::array ids, std::array const& donorIds, + int deuteronProng, int donorDeuteronProng) { + const int replacedProng = mixing.type == 0 ? deuteronProng : (mixing.type == 1 ? 1 : LastProng - deuteronProng); + const int donorProng = mixing.type == 0 ? donorDeuteronProng : (mixing.type == 1 ? 1 : LastProng - donorDeuteronProng); + ids[replacedProng] = donorIds[donorProng]; + // Preserve prong positions; the selector evaluates both mass hypotheses. + std::array key{referenceCollision.globalIndex(), ids[0], ids[1], ids[2]}; + mixedCandidates.try_emplace(key, referenceSeed); + }; + for (const int& deuteronProngA : {0, LastProng}) { + for (const int& deuteronProngB : {0, LastProng}) { + buildMixed(collision, seedA, idsA, idsB, deuteronProngA, deuteronProngB); + buildMixed(collisionB, seedB, idsB, idsA, deuteronProngB, deuteronProngA); + } + } + } + } + if (hasUsablePair) { + ++nPartners; + auto previousCollision = previousEvent.front().template collision_as(); + registry.fill(HIST("Mixing/hEventPairDeltaPVZ"), collision.posZ() - previousCollision.posZ()); + float previousCentrality = -1.f; + if constexpr (CentEstimator == CentralityEstimator::FT0C) { + previousCentrality = previousCollision.centFT0C(); + } else if constexpr (CentEstimator == CentralityEstimator::FT0M) { + previousCentrality = previousCollision.centFT0M(); + } + registry.fill(HIST("Mixing/hEventPairCentrality"), centrality, previousCentrality); + } + } + registry.fill(HIST("Mixing/hNPartners"), nPartners); + pool.push_back(std::move(seeds)); + if (pool.size() > static_cast(mixing.depth.value)) { + pool.pop_front(); + } + } + } + for (const auto& entry : pools) { + registry.fill(HIST("Mixing/hPoolOccupancy"), entry.second.size()); + } + // The ordinary reconstruction function consumes these adapters without any changes. + std::vector> mixedSeeds; + mixedSeeds.reserve(mixedCandidates.size()); + for (const auto& [key, source] : mixedCandidates) { + mixedSeeds.push_back({source, &tracks, {key[1], key[2], key[3]}, channelFlag}); + registry.fill(HIST("Mixing/hCounter"), 3.5); + } + const auto firstIndex = rowCandidateBase.lastIndex(); + runCreator3ProngWithDCAFitterN(collisions, mixedSeeds, tracks, bcs); + registry.fill(HIST("Mixing/hCounter"), 4.5, rowCandidateBase.lastIndex() - firstIndex); + } + /// @brief process function using DCA fitter w/ PV refit and w/o centrality selections void processPvRefitWithDCAFitterN(soa::Join const& collisions, FilteredPvRefitHf3Prongs const& rowsTrackIndexProng3, @@ -731,7 +947,11 @@ struct HfCandidateCreator3Prong { TracksWCovExtraPidPiKaPrLightNuclei const& tracks, aod::BCsWithTimestamps const& bcWithTimeStamps) { - runCreator3ProngWithDCAFitterN(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps); + if (mixing.enabled) { + runCreator3ProngMixedWithDCAFitterN(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps, static_cast(BIT(DecayType::CdToDeKPi))); + } else { + runCreator3ProngWithDCAFitterN(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps); + } } PROCESS_SWITCH(HfCandidateCreator3Prong, processNoPvRefitWithDCAFitterN, "Run candidate creator using DCA fitter without PV refit and w/o centrality selections", true); @@ -777,7 +997,11 @@ struct HfCandidateCreator3Prong { TracksWCovExtraPidPiKaPrLightNuclei const& tracks, aod::BCsWithTimestamps const& bcWithTimeStamps) { - runCreator3ProngWithDCAFitterN(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps); + if (mixing.enabled) { + runCreator3ProngMixedWithDCAFitterN(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps, static_cast(BIT(DecayType::CdToDeKPi))); + } else { + runCreator3ProngWithDCAFitterN(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps); + } } PROCESS_SWITCH(HfCandidateCreator3Prong, processNoPvRefitWithDCAFitterNCentFT0C, "Run candidate creator using DCA fitter without PV refit and w/ centrality selection on FT0C", false); @@ -823,7 +1047,11 @@ struct HfCandidateCreator3Prong { TracksWCovExtraPidPiKaPrLightNuclei const& tracks, aod::BCsWithTimestamps const& bcWithTimeStamps) { - runCreator3ProngWithDCAFitterN(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps); + if (mixing.enabled) { + runCreator3ProngMixedWithDCAFitterN(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps, static_cast(BIT(DecayType::CdToDeKPi))); + } else { + runCreator3ProngWithDCAFitterN(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps); + } } PROCESS_SWITCH(HfCandidateCreator3Prong, processNoPvRefitWithDCAFitterNCentFT0M, "Run candidate creator using DCA fitter without PV refit and w/ centrality selection on FT0M", false);