From 0e3f140d6e685276cdebd89418b86d9140fbd7bd Mon Sep 17 00:00:00 2001 From: blacw Date: Thu, 10 Sep 2026 15:32:50 +0800 Subject: [PATCH 1/3] add DCA fit process --- .../TableProducer/HadNucleiFemto.cxx | 558 +++++++++++++++++- 1 file changed, 556 insertions(+), 2 deletions(-) diff --git a/PWGCF/Femto/FemtoNuclei/TableProducer/HadNucleiFemto.cxx b/PWGCF/Femto/FemtoNuclei/TableProducer/HadNucleiFemto.cxx index 53a26284222..07ceba0aabc 100644 --- a/PWGCF/Femto/FemtoNuclei/TableProducer/HadNucleiFemto.cxx +++ b/PWGCF/Femto/FemtoNuclei/TableProducer/HadNucleiFemto.cxx @@ -10,8 +10,8 @@ // or submit itself to any jurisdiction. // -/// \file HadNucleiFemto.cxx -/// \brief Analysis task for Nuclei-Hadron femto analysis +/// \file HadNucleiFemtoDcaPurity.cxx +/// \brief Nuclei-hadron femtoscopy task with DCA-fraction and purity inputs /// \author CMY /// \date 2025-04-10 @@ -57,6 +57,9 @@ #include #include #include +#include +#include +#include #include #include @@ -98,6 +101,16 @@ constexpr int DeuteronPDG = o2::constants::physics::Pdg::kDeuteron; constexpr int TritonPDG = o2::constants::physics::Pdg::kTriton; constexpr int He3PDG = o2::constants::physics::Pdg::kHelium3; constexpr int HyperTritonPDG = o2::constants::physics::Pdg::kHyperTriton; +constexpr int Lithium4PDG = o2::constants::physics::Pdg::kLithium4; +// Store an exact integer PDG code in sparse histograms as +// sign * (high * kMotherPdgChunkBase + low). Splitting the code avoids the +// loss of integer precision that a single float coordinate has near 10^9. +constexpr int MotherPdgChunkBase = 10000; // o2-linter: disable=pdg/explicit-code (encoding base, not a PDG code) +constexpr int MotherPdgHighMax = 120000; // o2-linter: disable=pdg/explicit-code (encoding-axis limit, not a PDG code) +constexpr int HadronDcaFitBins = 2400; // 0.002 cm/bin in [-2.4, 2.4] cm +constexpr float HadronDcaFitAxisMax = 2.4f; +constexpr int NucleusDcaFitBins = 2000; // 0.001 cm/bin in [-1, 1] cm +constexpr float NucleusDcaFitAxisMax = 1.f; constexpr float He3TPCChi2NClMin = 0.5f; using PairLorentzVector = ROOT::Math::LorentzVector>; @@ -108,6 +121,48 @@ enum Selections { kAll }; +enum DcaOrigin { + kPrimary = 0, + kWeakDecay, + kMaterial, + kNDcaOrigins +}; + +enum CollisionAssociation { + kCorrectCollision = 0, + kWrongCollision, + kNCollisionAssociations +}; + +// The base DCA origin stays primary/weak/material. This independent parent +// axis permits optional feed-down splits without changing the base templates. +enum ParentCategory { + kNoSpecialParent = 0, + kLambdaParent, + kSigmaParent, + kK0ShortParent, + kK0LongParent, + kChargedKaonParent, + kHypertritonParent, + kLithium4Parent, + kOtherDecayParent, + kMissingDecayParent, + kNParentCategories +}; + +enum PurityCategory { + kAllSelected = 0, + kCorrectSpecies, + kMisidentifiedSpecies, + kNoMCLabel, + kCorrectCollisionSpecies, + kWrongCollisionSpecies, + kCorrectPrimary, + kCorrectWeakDecay, + kCorrectMaterial, + kNPurityCategories +}; + } // namespace struct HadNucandidate { @@ -354,6 +409,14 @@ struct HadNucleiFemto { Configurable settingRequirePhysicalPrimaries{"settingRequirePhysicalPrimaries", false, "Store only pairs in which both truth particles are physical primaries"}; } mc; + struct : o2::framework::ConfigurableGroup { + // cppcheck-suppress unusedStructMember + std::string prefix{"fractionPurity"}; + Configurable settingHadronDcaFitAbsMax{"settingHadronDcaFitAbsMax", 2.4f, "Maximum absolute pion DCAxy and DCAz stored for the fraction fit"}; + Configurable settingNucleusDcaFitAbsMax{"settingNucleusDcaFitAbsMax", 1.0f, "Maximum absolute nucleus DCAxy and DCAz stored for the fraction fit"}; + Configurable settingRequireRecoMCCollisionMatch{"settingRequireRecoMCCollisionMatch", true, "For the base DCA templates, require the truth particle to belong to the reconstructed collision MC label; detailed templates always store both association classes"}; + } fractionPurity; + struct : o2::framework::ConfigurableGroup { // cppcheck-suppress unusedStructMember std::string prefix{"hypertriton"}; @@ -600,6 +663,40 @@ struct HadNucleiFemto { {"MC/hPtNuRecVsGen", "Reconstructed versus generated signed nucleus pT;generated pT (GeV/c);reconstructed pT (GeV/c)", {HistType::kTH2F, {{280, -7.f, 7.f}, {280, -7.f, 7.f}}}}, {"MC/hPtHadRecVsGen", "Reconstructed versus generated signed pion pT;generated pT (GeV/c);reconstructed pT (GeV/c)", {HistType::kTH2F, {{280, -7.f, 7.f}, {280, -7.f, 7.f}}}}, + // Base inputs for Giorgio-style offline DCA fits. The origin axis is: + // 0 primary, 1 weak decay, 2 material (every non-primary, non-decay process). + {"fraction/hDcaDataHad", "Pion DCA-fit input;signed reconstructed p_{T} (GeV/c);DCA_{xy} (cm);DCA_{z} (cm);centrality", {HistType::kTHnSparseF, {{280, -7.f, 7.f}, {HadronDcaFitBins, -HadronDcaFitAxisMax, HadronDcaFitAxisMax}, {HadronDcaFitBins, -HadronDcaFitAxisMax, HadronDcaFitAxisMax}, {40, 0.f, 100.f}}}}, + {"fraction/hDcaDataNu", "Nucleus DCA-fit input;signed physical reconstructed p_{T} (GeV/c);DCA_{xy} (cm);DCA_{z} (cm);centrality", {HistType::kTHnSparseF, {{280, -7.f, 7.f}, {NucleusDcaFitBins, -NucleusDcaFitAxisMax, NucleusDcaFitAxisMax}, {NucleusDcaFitBins, -NucleusDcaFitAxisMax, NucleusDcaFitAxisMax}, {40, 0.f, 100.f}}}}, + {"fraction/hDcaTemplateHad", "Truth pion DCA templates;signed reconstructed p_{T} (GeV/c);DCA_{xy} (cm);DCA_{z} (cm);centrality;origin", {HistType::kTHnSparseF, {{280, -7.f, 7.f}, {HadronDcaFitBins, -HadronDcaFitAxisMax, HadronDcaFitAxisMax}, {HadronDcaFitBins, -HadronDcaFitAxisMax, HadronDcaFitAxisMax}, {40, 0.f, 100.f}, {kNDcaOrigins, -0.5f, static_cast(kNDcaOrigins) - 0.5f}}}}, + {"fraction/hDcaTemplateNu", "Truth nucleus DCA templates;signed physical reconstructed p_{T} (GeV/c);DCA_{xy} (cm);DCA_{z} (cm);centrality;origin", {HistType::kTHnSparseF, {{280, -7.f, 7.f}, {NucleusDcaFitBins, -NucleusDcaFitAxisMax, NucleusDcaFitAxisMax}, {NucleusDcaFitBins, -NucleusDcaFitAxisMax, NucleusDcaFitAxisMax}, {40, 0.f, 100.f}, {kNDcaOrigins, -0.5f, static_cast(kNDcaOrigins) - 0.5f}}}}, + + // DCA shapes of selected candidates that are not the requested truth + // species. They are kept separate from primary/weak/material so that + // their normalisation can be fixed or constrained by the PID purity. + {"fraction/hDcaMisidentifiedHad", "Misidentified pion candidates;signed reconstructed p_{T} (GeV/c);DCA_{xy} (cm);DCA_{z} (cm);centrality;collision association", {HistType::kTHnSparseF, {{280, -7.f, 7.f}, {HadronDcaFitBins, -HadronDcaFitAxisMax, HadronDcaFitAxisMax}, {HadronDcaFitBins, -HadronDcaFitAxisMax, HadronDcaFitAxisMax}, {40, 0.f, 100.f}, {kNCollisionAssociations, -0.5f, static_cast(kNCollisionAssociations) - 0.5f}}}}, + {"fraction/hDcaMisidentifiedNu", "Misidentified nucleus candidates;signed physical reconstructed p_{T} (GeV/c);DCA_{xy} (cm);DCA_{z} (cm);centrality;collision association", {HistType::kTHnSparseF, {{280, -7.f, 7.f}, {NucleusDcaFitBins, -NucleusDcaFitAxisMax, NucleusDcaFitAxisMax}, {NucleusDcaFitBins, -NucleusDcaFitAxisMax, NucleusDcaFitAxisMax}, {40, 0.f, 100.f}, {kNCollisionAssociations, -0.5f, static_cast(kNCollisionAssociations) - 0.5f}}}}, + {"fraction/hDcaNoMCLabelHad", "Selected pion candidates without an MC label;signed reconstructed p_{T} (GeV/c);DCA_{xy} (cm);DCA_{z} (cm);centrality", {HistType::kTHnSparseF, {{280, -7.f, 7.f}, {HadronDcaFitBins, -HadronDcaFitAxisMax, HadronDcaFitAxisMax}, {HadronDcaFitBins, -HadronDcaFitAxisMax, HadronDcaFitAxisMax}, {40, 0.f, 100.f}}}}, + {"fraction/hDcaNoMCLabelNu", "Selected nucleus candidates without an MC label;signed physical reconstructed p_{T} (GeV/c);DCA_{xy} (cm);DCA_{z} (cm);centrality", {HistType::kTHnSparseF, {{280, -7.f, 7.f}, {NucleusDcaFitBins, -NucleusDcaFitAxisMax, NucleusDcaFitAxisMax}, {NucleusDcaFitBins, -NucleusDcaFitAxisMax, NucleusDcaFitAxisMax}, {40, 0.f, 100.f}}}}, + + // Superset templates for optional offline refinements. They retain the + // same base origin plus collision association, direct-parent category, + // and the particle production radius. Projecting away the extra axes + // recovers the base shapes; selecting them enables species-specific fits. + {"fraction/hDcaTemplateDetailHad", "Detailed truth pion DCA templates;signed reconstructed p_{T} (GeV/c);DCA_{xy} (cm);DCA_{z} (cm);centrality;origin;collision association;parent category;production R_{xy} (cm)", {HistType::kTHnSparseF, {{280, -7.f, 7.f}, {HadronDcaFitBins, -HadronDcaFitAxisMax, HadronDcaFitAxisMax}, {HadronDcaFitBins, -HadronDcaFitAxisMax, HadronDcaFitAxisMax}, {40, 0.f, 100.f}, {kNDcaOrigins, -0.5f, static_cast(kNDcaOrigins) - 0.5f}, {kNCollisionAssociations, -0.5f, static_cast(kNCollisionAssociations) - 0.5f}, {kNParentCategories, -0.5f, static_cast(kNParentCategories) - 0.5f}, {200, 0.f, 100.f}}}}, + {"fraction/hDcaTemplateDetailNu", "Detailed truth nucleus DCA templates;signed physical reconstructed p_{T} (GeV/c);DCA_{xy} (cm);DCA_{z} (cm);centrality;origin;collision association;parent category;production R_{xy} (cm)", {HistType::kTHnSparseF, {{280, -7.f, 7.f}, {NucleusDcaFitBins, -NucleusDcaFitAxisMax, NucleusDcaFitAxisMax}, {NucleusDcaFitBins, -NucleusDcaFitAxisMax, NucleusDcaFitAxisMax}, {40, 0.f, 100.f}, {kNDcaOrigins, -0.5f, static_cast(kNDcaOrigins) - 0.5f}, {kNCollisionAssociations, -0.5f, static_cast(kNCollisionAssociations) - 0.5f}, {kNParentCategories, -0.5f, static_cast(kNParentCategories) - 0.5f}, {200, 0.f, 100.f}}}}, + + // Mother-resolved templates. There is one entry per direct mother (or a + // single zero-code entry when no mother is stored). Reconstruct the exact + // PDG code offline as sign * (PDG high * 10000 + PDG low). + {"fraction/hDcaMotherPdgHad", "Pion DCA by exact direct-mother PDG; signed reconstructed p_{T} (GeV/c);DCA_{xy} (cm);DCA_{z} (cm);centrality;origin;collision association;mother PDG sign;mother PDG high;mother PDG low", {HistType::kTHnSparseF, {{280, -7.f, 7.f}, {HadronDcaFitBins, -HadronDcaFitAxisMax, HadronDcaFitAxisMax}, {HadronDcaFitBins, -HadronDcaFitAxisMax, HadronDcaFitAxisMax}, {40, 0.f, 100.f}, {kNDcaOrigins, -0.5f, static_cast(kNDcaOrigins) - 0.5f}, {kNCollisionAssociations, -0.5f, static_cast(kNCollisionAssociations) - 0.5f}, {3, -1.5f, 1.5f}, {MotherPdgHighMax + 1, -0.5f, static_cast(MotherPdgHighMax) + 0.5f}, {MotherPdgChunkBase, -0.5f, static_cast(MotherPdgChunkBase) - 0.5f}}}}, + {"fraction/hDcaMotherPdgNu", "Nucleus DCA by exact direct-mother PDG; signed physical reconstructed p_{T} (GeV/c);DCA_{xy} (cm);DCA_{z} (cm);centrality;origin;collision association;mother PDG sign;mother PDG high;mother PDG low", {HistType::kTHnSparseF, {{280, -7.f, 7.f}, {NucleusDcaFitBins, -NucleusDcaFitAxisMax, NucleusDcaFitAxisMax}, {NucleusDcaFitBins, -NucleusDcaFitAxisMax, NucleusDcaFitAxisMax}, {40, 0.f, 100.f}, {kNDcaOrigins, -0.5f, static_cast(kNDcaOrigins) - 0.5f}, {kNCollisionAssociations, -0.5f, static_cast(kNCollisionAssociations) - 0.5f}, {3, -1.5f, 1.5f}, {MotherPdgHighMax + 1, -0.5f, static_cast(MotherPdgHighMax) + 0.5f}, {MotherPdgChunkBase, -0.5f, static_cast(MotherPdgChunkBase) - 0.5f}}}}, + + // Hierarchical MC truth purity counters: all = correct species + mis-ID + // + no label; correct species = correct collision + wrong collision; + // correct collision = primary + weak + material. + {"purityMC/hHadron", "Selected pion truth composition;signed reconstructed p_{T} (GeV/c);category;centrality", {HistType::kTH3F, {{280, -7.f, 7.f}, {kNPurityCategories, -0.5f, static_cast(kNPurityCategories) - 0.5f}, {40, 0.f, 100.f}}}}, + {"purityMC/hNucleus", "Selected nucleus truth composition;signed physical reconstructed p_{T} (GeV/c);category;centrality", {HistType::kTH3F, {{280, -7.f, 7.f}, {kNPurityCategories, -0.5f, static_cast(kNPurityCategories) - 0.5f}, {40, 0.f, 100.f}}}}, + // dE/dx {"h2dEdxNucandidates", "dEdx distribution; #it{p} (GeV/#it{c}); dE/dx (a.u.)", {HistType::kTH2F, {{200, -5.0f, 5.0f}, {100, 0.0f, 2000.0f}}}}, {"h2dEdxHadcandidates", "dEdx distribution; #it{p} (GeV/#it{c}); dE/dx (a.u.)", {HistType::kTH2F, {{200, -5.0f, 5.0f}, {100, 0.0f, 2000.0f}}}}, @@ -734,6 +831,64 @@ struct HadNucleiFemto { for (size_t i = 0; i < mixedEventLabels.size(); i++) { mQaRegistry.get(HIST("hMixedEventSelections"))->GetXaxis()->SetBinLabel(i + 1, mixedEventLabels[i].c_str()); } + + const std::array originLabels = {"Primary", "WeakDecay", "Material"}; + const std::array collisionLabels = {"CorrectCollision", "WrongCollision"}; + const std::array parentLabels = {"NoSpecialParent", "Lambda", "Sigma", "K0Short", "K0Long", "ChargedKaon", "Hypertriton", "Lithium4", "OtherDecay", "MissingDecayMother"}; + const std::array purityLabels = {"AllSelected", "CorrectSpecies", "MisidentifiedSpecies", "NoMCLabel", "CorrectCollisionSpecies", "WrongCollisionSpecies", "CorrectPrimary", "CorrectWeakDecay", "CorrectMaterial"}; + + const auto baseDcaHad = mQaRegistry.get(HIST("fraction/hDcaTemplateHad")); + const auto baseDcaNu = mQaRegistry.get(HIST("fraction/hDcaTemplateNu")); + const std::array baseDcaHistograms = {baseDcaHad.get(), baseDcaNu.get()}; + for (const auto& histogram : baseDcaHistograms) { + for (int i = 0; i < kNDcaOrigins; ++i) { + histogram->GetAxis(4)->SetBinLabel(i + 1, originLabels[i]); + } + } + const auto misidentifiedDcaHad = mQaRegistry.get(HIST("fraction/hDcaMisidentifiedHad")); + const auto misidentifiedDcaNu = mQaRegistry.get(HIST("fraction/hDcaMisidentifiedNu")); + const std::array misidentifiedDcaHistograms = {misidentifiedDcaHad.get(), misidentifiedDcaNu.get()}; + for (const auto& histogram : misidentifiedDcaHistograms) { + for (int i = 0; i < kNCollisionAssociations; ++i) { + histogram->GetAxis(4)->SetBinLabel(i + 1, collisionLabels[i]); + } + } + const auto detailedDcaHad = mQaRegistry.get(HIST("fraction/hDcaTemplateDetailHad")); + const auto detailedDcaNu = mQaRegistry.get(HIST("fraction/hDcaTemplateDetailNu")); + const std::array detailedDcaHistograms = {detailedDcaHad.get(), detailedDcaNu.get()}; + for (const auto& histogram : detailedDcaHistograms) { + for (int i = 0; i < kNDcaOrigins; ++i) { + histogram->GetAxis(4)->SetBinLabel(i + 1, originLabels[i]); + } + for (int i = 0; i < kNCollisionAssociations; ++i) { + histogram->GetAxis(5)->SetBinLabel(i + 1, collisionLabels[i]); + } + for (int i = 0; i < kNParentCategories; ++i) { + histogram->GetAxis(6)->SetBinLabel(i + 1, parentLabels[i]); + } + } + const auto motherPdgHad = mQaRegistry.get(HIST("fraction/hDcaMotherPdgHad")); + const auto motherPdgNu = mQaRegistry.get(HIST("fraction/hDcaMotherPdgNu")); + const std::array motherPdgHistograms = {motherPdgHad.get(), motherPdgNu.get()}; + for (const auto& histogram : motherPdgHistograms) { + for (int i = 0; i < kNDcaOrigins; ++i) { + histogram->GetAxis(4)->SetBinLabel(i + 1, originLabels[i]); + } + for (int i = 0; i < kNCollisionAssociations; ++i) { + histogram->GetAxis(5)->SetBinLabel(i + 1, collisionLabels[i]); + } + histogram->GetAxis(6)->SetBinLabel(1, "Negative"); + histogram->GetAxis(6)->SetBinLabel(2, "NoMother"); + histogram->GetAxis(6)->SetBinLabel(3, "Positive"); + } + const auto purityHad = mQaRegistry.get(HIST("purityMC/hHadron")); + const auto purityNu = mQaRegistry.get(HIST("purityMC/hNucleus")); + const std::array purityHistograms = {purityHad.get(), purityNu.get()}; + for (const auto& histogram : purityHistograms) { + for (int i = 0; i < kNPurityCategories; ++i) { + histogram->GetYaxis()->SetBinLabel(i + 1, purityLabels[i]); + } + } } template @@ -1054,6 +1209,60 @@ struct HadNucleiFemto { return false; } + // DCA-template inputs must retain the tails. These predicates reproduce the + // nominal track-quality selections while deliberately omitting only DCA. + template + bool selectPionTrackForDcaFit(const Ttrack& candidate) const + { + const float absPt = std::abs(candidate.pt()); + return std::abs(candidate.eta()) <= trackCut.settingCutEta.value && + absPt >= hadronPid.settingHadptMin.value && + absPt <= hadronPid.settingHadptMax.value && + candidate.itsNClsInnerBarrel() >= hadronPid.settingPionITSInnerBarrelMin.value && + candidate.itsNCls() >= hadronPid.settingPionITSNClsMin.value && + candidate.tpcNClsFound() >= hadronPid.settingPionTPCNClsFoundMin.value && + candidate.tpcNClsCrossedRows() >= hadronPid.settingPionTPCCrossedRowsMin.value; + } + + template + bool selectTritonTrackForDcaFit(const Ttrack& candidate) const + { + constexpr float maxAbsEta = 0.8f; + constexpr int minTPCCrossedRows = 70; + constexpr float maxTPCChi2NCl = 5.f; + constexpr float maxTPCFractionSharedCls = 0.3f; + constexpr int minITSNCls = 5; + constexpr float maxITSChi2NCl = 10.f; + return std::abs(candidate.eta()) < maxAbsEta && + candidate.tpcNClsCrossedRows() >= minTPCCrossedRows && + candidate.tpcChi2NCl() < maxTPCChi2NCl && + candidate.tpcFractionSharedCls() < maxTPCFractionSharedCls && + candidate.itsNCls() >= minITSNCls && + candidate.itsChi2NCl() < maxITSChi2NCl; + } + + template + bool selectHadronTrackForDcaFit(const Ttrack& candidate) + { + if (species.settingHadPDGCode.value == static_cast(PDG_t::kPiPlus)) { + return selectPionTrackForDcaFit(candidate); + } + // The fraction feature is intended for pion-nucleus configurations. Keep + // the nominal selection for other supported hadrons instead of silently + // changing their established cuts. + return selectTrackHadron(candidate); + } + + template + bool selectNucleusTrackForDcaFit(const Ttrack& candidate) + { + if (useTritonNucleus()) { + return selectTritonTrackForDcaFit(candidate); + } + // The deuteron and helium-3 track-quality predicates contain no DCA cut. + return selectTrackNu(candidate); + } + void fillNucleusTrackSelection(const Selections selection) { mQaRegistry.fill(HIST("hTrackSelNu"), selection); @@ -1618,6 +1827,128 @@ struct HadNucleiFemto { return false; } + template + bool selectionPIDDeForDcaFit(const Ttrack& candidate) + { + const float absPt = std::abs(candidate.pt()); + const float absTPCInnerParam = std::abs(candidate.tpcInnerParam()); + if (absTPCInnerParam < deuteronPid.settingCutPinMinDe.value || + absPt < deuteronPid.settingCutDeptMin.value || + absPt > deuteronPid.settingCutDeptMax.value) { + return false; + } + + const float tpcNSigmaDe = output.settingUseBBcomputeDeNsigma.value ? computeNSigmaDe(candidate) : candidate.tpcNSigmaDe(); + if (absTPCInnerParam > deuteronPid.settingCutPinMinTOFITSDe.value) { + if (!candidate.hasTOF()) { + return false; + } + const float tofNSigmaDe = candidate.tofNSigmaDe(); + const float combinedNSigma = std::hypot(tpcNSigmaDe, tofNSigmaDe); + if (combinedNSigma > deuteronPid.settingCutNsigmaTOFTPCDe.value) { + return false; + } + return !deuteronPid.settingReqSingleNsig.value || + (std::abs(tpcNSigmaDe) <= deuteronPid.settingCutNsigmaTOFTPCDe.value && + std::abs(tofNSigmaDe) <= deuteronPid.settingCutNsigmaTOFTPCDe.value); + } + + if (std::abs(tpcNSigmaDe) > deuteronPid.settingCutNsigmaTPCDe.value) { + return false; + } + o2::aod::ITSResponse itsResponse; + const float itsNSigmaDe = itsResponse.nSigmaITS(candidate.itsClusterSizes(), candidate.p(), candidate.eta()); + return std::abs(itsNSigmaDe) <= deuteronPid.settingCutNsigmaITSDe.value; + } + + template + bool selectionPIDNuForDcaFit(const Ttrack& candidate) + { + if (useDeuteronNucleus()) { + return selectionPIDDeForDcaFit(candidate); + } + // The helium-3 and triton PID predicates do not contain DCA selections. + return selectionPIDNu(candidate); + } + + template + float signedPhysicalPt(const Ttrack& track, bool isNucleus) const + { + const float chargeFactor = isNucleus ? nucleusChargeFactor() : 1.f; + return track.sign() * chargeFactor * std::abs(track.pt()); + } + + template + int classifyDcaOrigin(const Tparticle& particle) const + { + if (particle.isPhysicalPrimary()) { + return kPrimary; + } + if (particle.getProcess() == TMCProcess::kPDecay) { + return kWeakDecay; + } + // Following the DCA-template convention used by Giorgio: after excluding + // physical primaries and decay daughters, every remaining transport + // process is treated as a secondary produced in detector material. + return kMaterial; + } + + int purityOriginCategory(int origin) const + { + switch (origin) { + case kPrimary: + return kCorrectPrimary; + case kWeakDecay: + return kCorrectWeakDecay; + case kMaterial: + return kCorrectMaterial; + default: + // classifyDcaOrigin deliberately has only the three Giorgio classes. + return kCorrectMaterial; + } + } + + template + int classifyParentCategory(const Tparticle& particle, int origin) const + { + if (!particle.has_mothers()) { + return origin == kWeakDecay ? kMissingDecayParent : kNoSpecialParent; + } + + const int fallbackCategory = origin == kWeakDecay ? kOtherDecayParent : kNoSpecialParent; + for (const auto& mother : particle.template mothers_as()) { + const int motherPdg = std::abs(mother.pdgCode()); + + // Keep the nuclear feed-down parents identifiable for all origin labels. + if (motherPdg == HyperTritonPDG) { + return kHypertritonParent; + } + if (motherPdg == Lithium4PDG) { + return kLithium4Parent; + } + + if (origin != kWeakDecay) { + continue; + } + if (motherPdg == PDG_t::kLambda0) { + return kLambdaParent; + } + if (motherPdg == PDG_t::kSigmaMinus || motherPdg == PDG_t::kSigma0 || motherPdg == PDG_t::kSigmaPlus) { + return kSigmaParent; + } + if (motherPdg == PDG_t::kK0Short) { + return kK0ShortParent; + } + if (motherPdg == PDG_t::kK0Long) { + return kK0LongParent; + } + if (motherPdg == PDG_t::kKPlus) { + return kChargedKaonParent; + } + } + return fallbackCategory; + } + template float getNucleusTPCNSigma(const Ttrack& candidate) { @@ -2882,6 +3213,229 @@ struct HadNucleiFemto { } PROCESS_SWITCH(HadNucleiFemto, processMixedEvent, "Process Mixed event", false); + // Produce the data distributions fitted offline with the MC templates. All + // nominal quality and PID selections are applied; only DCA is relaxed. + void processDcaFractionData(const CollisionsFull& collisions, const TrackCandidates& tracks, const aod::BCsWithTimestamps& bcs) + { + for (const auto& collision : collisions) { + if (!selectCollision(collision, bcs)) { + continue; + } + + const uint64_t collIdx = collision.globalIndex(); + auto tracksThisCollision = tracks.sliceBy(mPerCol, collIdx); + tracksThisCollision.bindExternalIndices(&tracks); + for (const auto& track : tracksThisCollision) { + if (selectHadronTrackForDcaFit(track) && selectionPIDHadron(track) && + std::abs(track.dcaXY()) <= fractionPurity.settingHadronDcaFitAbsMax.value && + std::abs(track.dcaZ()) <= fractionPurity.settingHadronDcaFitAbsMax.value) { + mQaRegistry.fill(HIST("fraction/hDcaDataHad"), signedPhysicalPt(track, false), track.dcaXY(), track.dcaZ(), collision.centFT0C()); + } + + if (selectNucleusTrackForDcaFit(track) && selectionPIDNuForDcaFit(track) && + std::abs(track.dcaXY()) <= fractionPurity.settingNucleusDcaFitAbsMax.value && + std::abs(track.dcaZ()) <= fractionPurity.settingNucleusDcaFitAbsMax.value) { + mQaRegistry.fill(HIST("fraction/hDcaDataNu"), signedPhysicalPt(track, true), track.dcaXY(), track.dcaZ(), collision.centFT0C()); + } + } + } + } + PROCESS_SWITCH(HadNucleiFemto, processDcaFractionData, "Produce data DCA-fraction fit inputs", false); + + template + bool truthBelongsToRecoCollision(const Ttrack&, const Tparticle& particle, const Tcollision& collision) const + { + return collision.has_mcCollision() && particle.mcCollisionId() == collision.mcCollisionId(); + } + + template + void fillMCDcaMisidentified(const Ttrack& track, float centrality, bool isNucleus, bool matchesRecoCollision) + { + const float signedPt = signedPhysicalPt(track, isNucleus); + const float collisionAssociation = matchesRecoCollision ? static_cast(kCorrectCollision) : static_cast(kWrongCollision); + if (isNucleus) { + mQaRegistry.fill(HIST("fraction/hDcaMisidentifiedNu"), signedPt, track.dcaXY(), track.dcaZ(), centrality, collisionAssociation); + } else { + mQaRegistry.fill(HIST("fraction/hDcaMisidentifiedHad"), signedPt, track.dcaXY(), track.dcaZ(), centrality, collisionAssociation); + } + } + + template + void fillMCDcaNoLabel(const Ttrack& track, float centrality, bool isNucleus) + { + const float signedPt = signedPhysicalPt(track, isNucleus); + if (isNucleus) { + mQaRegistry.fill(HIST("fraction/hDcaNoMCLabelNu"), signedPt, track.dcaXY(), track.dcaZ(), centrality); + } else { + mQaRegistry.fill(HIST("fraction/hDcaNoMCLabelHad"), signedPt, track.dcaXY(), track.dcaZ(), centrality); + } + } + + template + void fillMCDcaTemplate(const Ttrack& track, const Tparticle& particle, float centrality, bool isNucleus, bool matchesRecoCollision) + { + const int origin = classifyDcaOrigin(particle); + const int collisionAssociation = matchesRecoCollision ? kCorrectCollision : kWrongCollision; + const int parentCategory = classifyParentCategory(particle, origin); + const float productionRadius = std::hypot(particle.vx(), particle.vy()); + const float signedPt = signedPhysicalPt(track, isNucleus); + + // The detailed histogram is always filled for a truth-PDG-matched track, + // including wrong-collision associations, so tighter choices can be made + // offline without rerunning the table producer. + if (isNucleus) { + mQaRegistry.fill(HIST("fraction/hDcaTemplateDetailNu"), signedPt, track.dcaXY(), track.dcaZ(), centrality, static_cast(origin), static_cast(collisionAssociation), static_cast(parentCategory), productionRadius); + } else { + mQaRegistry.fill(HIST("fraction/hDcaTemplateDetailHad"), signedPt, track.dcaXY(), track.dcaZ(), centrality, static_cast(origin), static_cast(collisionAssociation), static_cast(parentCategory), productionRadius); + } + + const auto fillMotherPdg = [&](int motherPdg) { + const int motherSign = (motherPdg > 0) - (motherPdg < 0); + const int absoluteMotherPdg = std::abs(motherPdg); + const int motherPdgHigh = absoluteMotherPdg / MotherPdgChunkBase; + const int motherPdgLow = absoluteMotherPdg % MotherPdgChunkBase; + if (motherPdgHigh > MotherPdgHighMax) { + LOG(warning) << "Direct-mother PDG " << motherPdg << " exceeds the configured exact-PDG histogram range"; + return; + } + if (isNucleus) { + mQaRegistry.fill(HIST("fraction/hDcaMotherPdgNu"), signedPt, track.dcaXY(), track.dcaZ(), centrality, static_cast(origin), static_cast(collisionAssociation), static_cast(motherSign), static_cast(motherPdgHigh), static_cast(motherPdgLow)); + } else { + mQaRegistry.fill(HIST("fraction/hDcaMotherPdgHad"), signedPt, track.dcaXY(), track.dcaZ(), centrality, static_cast(origin), static_cast(collisionAssociation), static_cast(motherSign), static_cast(motherPdgHigh), static_cast(motherPdgLow)); + } + }; + + if (particle.has_mothers()) { + for (const auto& mother : particle.template mothers_as()) { + fillMotherPdg(mother.pdgCode()); + } + } else { + fillMotherPdg(0); + } + + // This is the stable three-component Giorgio template. The collision-match + // configurable preserves the previous strict/relaxed behavior. + if (fractionPurity.settingRequireRecoMCCollisionMatch.value && !matchesRecoCollision) { + return; + } + if (isNucleus) { + mQaRegistry.fill(HIST("fraction/hDcaTemplateNu"), signedPt, track.dcaXY(), track.dcaZ(), centrality, static_cast(origin)); + } else { + mQaRegistry.fill(HIST("fraction/hDcaTemplateHad"), signedPt, track.dcaXY(), track.dcaZ(), centrality, static_cast(origin)); + } + } + + template + void fillMCPurityComposition(const Ttrack& track, const Tparticle& particle, float centrality, bool isNucleus, bool correctSpecies, bool matchesRecoCollision) + { + const float signedPt = signedPhysicalPt(track, isNucleus); + const auto fillCategory = [&](int category) { + if (isNucleus) { + mQaRegistry.fill(HIST("purityMC/hNucleus"), signedPt, static_cast(category), centrality); + } else { + mQaRegistry.fill(HIST("purityMC/hHadron"), signedPt, static_cast(category), centrality); + } + }; + fillCategory(kAllSelected); + if (!correctSpecies) { + fillCategory(kMisidentifiedSpecies); + return; + } + fillCategory(kCorrectSpecies); + if (!matchesRecoCollision) { + fillCategory(kWrongCollisionSpecies); + return; + } + fillCategory(kCorrectCollisionSpecies); + const int origin = classifyDcaOrigin(particle); + fillCategory(purityOriginCategory(origin)); + } + + template + void fillMCNoLabelPurity(const Ttrack& track, float centrality, bool isNucleus) + { + const float signedPt = signedPhysicalPt(track, isNucleus); + if (isNucleus) { + mQaRegistry.fill(HIST("purityMC/hNucleus"), signedPt, static_cast(kAllSelected), centrality); + mQaRegistry.fill(HIST("purityMC/hNucleus"), signedPt, static_cast(kNoMCLabel), centrality); + } else { + mQaRegistry.fill(HIST("purityMC/hHadron"), signedPt, static_cast(kAllSelected), centrality); + mQaRegistry.fill(HIST("purityMC/hHadron"), signedPt, static_cast(kNoMCLabel), centrality); + } + } + + // MC DCA templates use the relaxed-DCA selection. MC purity uses the full + // nominal candidate selection, including its DCA requirement. + void processDcaFractionPurityMC(const CollisionsFullMC& collisions, const TrackCandidatesMC& tracks, const aod::McParticles&, const aod::BCsWithTimestamps& bcs) + { + for (const auto& collision : collisions) { + if (!selectCollision(collision, bcs)) { + continue; + } + + const uint64_t collIdx = collision.globalIndex(); + auto tracksThisCollision = tracks.sliceBy(mPerColMC, collIdx); + tracksThisCollision.bindExternalIndices(&tracks); + for (const auto& track : tracksThisCollision) { + const bool selectedHadronForFraction = selectHadronTrackForDcaFit(track) && selectionPIDHadron(track) && + std::abs(track.dcaXY()) <= fractionPurity.settingHadronDcaFitAbsMax.value && + std::abs(track.dcaZ()) <= fractionPurity.settingHadronDcaFitAbsMax.value; + const bool selectedNucleusForFraction = selectNucleusTrackForDcaFit(track) && selectionPIDNuForDcaFit(track) && + std::abs(track.dcaXY()) <= fractionPurity.settingNucleusDcaFitAbsMax.value && + std::abs(track.dcaZ()) <= fractionPurity.settingNucleusDcaFitAbsMax.value; + + const bool selectedHadronForPurity = selectTrackHadron(track) && selectionPIDHadron(track); + const bool selectedNucleusForPurity = selectTrackNu(track) && selectionPIDNu(track); + + if (!track.has_mcParticle()) { + if (selectedHadronForFraction) { + fillMCDcaNoLabel(track, collision.centFT0C(), false); + } + if (selectedNucleusForFraction) { + fillMCDcaNoLabel(track, collision.centFT0C(), true); + } + if (selectedHadronForPurity) { + fillMCNoLabelPurity(track, collision.centFT0C(), false); + } + if (selectedNucleusForPurity) { + fillMCNoLabelPurity(track, collision.centFT0C(), true); + } + continue; + } + + const auto particle = track.template mcParticle_as(); + const bool matchesRecoCollision = truthBelongsToRecoCollision(track, particle, collision); + const int expectedHadronPdg = track.sign() >= 0 ? std::abs(species.settingHadPDGCode.value) : -std::abs(species.settingHadPDGCode.value); + const int expectedNucleusPdg = track.sign() >= 0 ? std::abs(species.settingNuPDGCode.value) : -std::abs(species.settingNuPDGCode.value); + const bool correctHadronSpecies = particle.pdgCode() == expectedHadronPdg; + const bool correctNucleusSpecies = particle.pdgCode() == expectedNucleusPdg; + + if (selectedHadronForPurity) { + fillMCPurityComposition(track, particle, collision.centFT0C(), false, correctHadronSpecies, matchesRecoCollision); + } + if (selectedNucleusForPurity) { + fillMCPurityComposition(track, particle, collision.centFT0C(), true, correctNucleusSpecies, matchesRecoCollision); + } + + if (selectedHadronForFraction) { + if (correctHadronSpecies) { + fillMCDcaTemplate(track, particle, collision.centFT0C(), false, matchesRecoCollision); + } else { + fillMCDcaMisidentified(track, collision.centFT0C(), false, matchesRecoCollision); + } + } + if (selectedNucleusForFraction) { + if (correctNucleusSpecies) { + fillMCDcaTemplate(track, particle, collision.centFT0C(), true, matchesRecoCollision); + } else { + fillMCDcaMisidentified(track, collision.centFT0C(), true, matchesRecoCollision); + } + } + } + } + } + PROCESS_SWITCH(HadNucleiFemto, processDcaFractionPurityMC, "Produce MC DCA templates and truth-purity counters", false); + void processPurity(const CollisionsFull& collisions, const TrackCandidates& tracks, const aod::BCsWithTimestamps& bcs) { for (const auto& collision : collisions) { From ac9d30b6a92c5b61e9536a1b3c7d0ee549ffbf03 Mon Sep 17 00:00:00 2001 From: Meiyi Chen Date: Thu, 10 Sep 2026 15:59:27 +0800 Subject: [PATCH 2/3] Update HadNucleiFemto.cxx --- PWGCF/Femto/FemtoNuclei/TableProducer/HadNucleiFemto.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGCF/Femto/FemtoNuclei/TableProducer/HadNucleiFemto.cxx b/PWGCF/Femto/FemtoNuclei/TableProducer/HadNucleiFemto.cxx index 07ceba0aabc..71063e1fbd6 100644 --- a/PWGCF/Femto/FemtoNuclei/TableProducer/HadNucleiFemto.cxx +++ b/PWGCF/Femto/FemtoNuclei/TableProducer/HadNucleiFemto.cxx @@ -10,8 +10,8 @@ // or submit itself to any jurisdiction. // -/// \file HadNucleiFemtoDcaPurity.cxx -/// \brief Nuclei-hadron femtoscopy task with DCA-fraction and purity inputs +/// \file HadNucleiFemto.cxx +/// \brief Analysis task for Nuclei-Hadron femto analysis /// \author CMY /// \date 2025-04-10 From db417f552d82f32c23054d4930fb3ca49e4de4c3 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Thu, 10 Sep 2026 08:01:08 +0000 Subject: [PATCH 3/3] Please consider the following formatting changes --- PWGCF/Femto/FemtoNuclei/TableProducer/HadNucleiFemto.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGCF/Femto/FemtoNuclei/TableProducer/HadNucleiFemto.cxx b/PWGCF/Femto/FemtoNuclei/TableProducer/HadNucleiFemto.cxx index 71063e1fbd6..a523b068168 100644 --- a/PWGCF/Femto/FemtoNuclei/TableProducer/HadNucleiFemto.cxx +++ b/PWGCF/Femto/FemtoNuclei/TableProducer/HadNucleiFemto.cxx @@ -107,7 +107,7 @@ constexpr int Lithium4PDG = o2::constants::physics::Pdg::kLithium4; // loss of integer precision that a single float coordinate has near 10^9. constexpr int MotherPdgChunkBase = 10000; // o2-linter: disable=pdg/explicit-code (encoding base, not a PDG code) constexpr int MotherPdgHighMax = 120000; // o2-linter: disable=pdg/explicit-code (encoding-axis limit, not a PDG code) -constexpr int HadronDcaFitBins = 2400; // 0.002 cm/bin in [-2.4, 2.4] cm +constexpr int HadronDcaFitBins = 2400; // 0.002 cm/bin in [-2.4, 2.4] cm constexpr float HadronDcaFitAxisMax = 2.4f; constexpr int NucleusDcaFitBins = 2000; // 0.001 cm/bin in [-1, 1] cm constexpr float NucleusDcaFitAxisMax = 1.f;