diff --git a/PWGCF/Femto/FemtoNuclei/DataModel/HadronNucleiTables.h b/PWGCF/Femto/FemtoNuclei/DataModel/HadronNucleiTables.h index c83fca57f9b..161be5544ab 100644 --- a/PWGCF/Femto/FemtoNuclei/DataModel/HadronNucleiTables.h +++ b/PWGCF/Femto/FemtoNuclei/DataModel/HadronNucleiTables.h @@ -71,8 +71,8 @@ DECLARE_SOA_COLUMN(MassTOFNu, massTOFNu, float); DECLARE_SOA_COLUMN(MassTOFHad, massTOFHad, float); DECLARE_SOA_COLUMN(PidTrkNu, pidTrkNu, uint32_t); DECLARE_SOA_COLUMN(PidTrkHad, pidTrkHad, uint32_t); -DECLARE_SOA_COLUMN(TrackIDHad, trackIDHad, int); -DECLARE_SOA_COLUMN(TrackIDNu, trackIDNu, int); +DECLARE_SOA_COLUMN(TrackIDHad, trackIDHad, int64_t); +DECLARE_SOA_COLUMN(TrackIDNu, trackIDNu, int64_t); DECLARE_SOA_COLUMN(ItsClusterSizeNu, itsClusterSizeNu, uint32_t); DECLARE_SOA_COLUMN(ItsClusterSizeHad, itsClusterSizeHad, uint32_t); @@ -82,6 +82,23 @@ DECLARE_SOA_COLUMN(SharedClustersHad, sharedClustersHad, uint8_t); DECLARE_SOA_COLUMN(DeltaEta, deltaEta, float); DECLARE_SOA_COLUMN(DeltaPhi, deltaPhi, float); +DECLARE_SOA_COLUMN(RunNumber, runNumber, int32_t); +DECLARE_SOA_COLUMN(MagneticField, magneticField, float); +DECLARE_SOA_COLUMN(CollisionIdNu, collisionIdNu, int64_t); +DECLARE_SOA_COLUMN(CollisionIdHad, collisionIdHad, int64_t); +DECLARE_SOA_COLUMN(NClsFindableTPCNu, nClsFindableTPCNu, uint8_t); +DECLARE_SOA_COLUMN(NClsFindableTPCHad, nClsFindableTPCHad, uint8_t); +DECLARE_SOA_COLUMN(FractionSharedTPCNu, fractionSharedTPCNu, float); +DECLARE_SOA_COLUMN(FractionSharedTPCHad, fractionSharedTPCHad, float); +DECLARE_SOA_COLUMN(NClsITSNu, nClsITSNu, uint8_t); +DECLARE_SOA_COLUMN(NClsITSHad, nClsITSHad, uint8_t); +DECLARE_SOA_COLUMN(NClsITSInnerBarrelNu, nClsITSInnerBarrelNu, uint8_t); +DECLARE_SOA_COLUMN(NClsITSInnerBarrelHad, nClsITSInnerBarrelHad, uint8_t); +DECLARE_SOA_COLUMN(Chi2ITSNu, chi2ITSNu, float); +DECLARE_SOA_COLUMN(Chi2ITSHad, chi2ITSHad, float); +DECLARE_SOA_COLUMN(NSigmaTPCNuDe, nSigmaTPCNuDe, float); +DECLARE_SOA_COLUMN(NSigmaTPCNuPr, nSigmaTPCNuPr, float); +DECLARE_SOA_COLUMN(NSigmaTPCNuPi, nSigmaTPCNuPi, float); // Reconstructed-MC pair information. The signed generated pT follows the // convention used by PtNu/PtHad: particles are positive and antiparticles @@ -154,7 +171,26 @@ DECLARE_SOA_TABLE(HadronNucleiTable, "AOD", "HADNUCLEITABLE", hadron_nuclei_tables::NSigmaTOFNu, hadron_nuclei_tables::NSigmaITSNu, hadron_nuclei_tables::NSigmaTOFHad, - hadron_nuclei_tables::NSigmaITSHad) + hadron_nuclei_tables::NSigmaITSHad, + hadron_nuclei_tables::RunNumber, + hadron_nuclei_tables::MagneticField, + hadron_nuclei_tables::TrackIDNu, + hadron_nuclei_tables::TrackIDHad, + hadron_nuclei_tables::CollisionIdNu, + hadron_nuclei_tables::CollisionIdHad, + hadron_nuclei_tables::NClsFindableTPCNu, + hadron_nuclei_tables::NClsFindableTPCHad, + hadron_nuclei_tables::FractionSharedTPCNu, + hadron_nuclei_tables::FractionSharedTPCHad, + hadron_nuclei_tables::NClsITSNu, + hadron_nuclei_tables::NClsITSHad, + hadron_nuclei_tables::NClsITSInnerBarrelNu, + hadron_nuclei_tables::NClsITSInnerBarrelHad, + hadron_nuclei_tables::Chi2ITSNu, + hadron_nuclei_tables::Chi2ITSHad, + hadron_nuclei_tables::NSigmaTPCNuDe, + hadron_nuclei_tables::NSigmaTPCNuPr, + hadron_nuclei_tables::NSigmaTPCNuPi) DECLARE_SOA_TABLE(HadronNucleiTableMC, "AOD", "HADNUCLEIMC", hadron_nuclei_tables::PtNuMC, hadron_nuclei_tables::EtaNuMC, diff --git a/PWGCF/Femto/FemtoNuclei/TableProducer/HadNucleiFemto.cxx b/PWGCF/Femto/FemtoNuclei/TableProducer/HadNucleiFemto.cxx index 53a26284222..fdc18ea18d2 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 { + Primary = 0, + WeakDecay, + Material, + NDcaOrigins +}; + +enum CollisionAssociation { + CorrectCollision = 0, + WrongCollision, + NCollisionAssociations +}; + +// The base DCA origin stays primary/weak/material. This independent parent +// axis permits optional feed-down splits without changing the base templates. +enum ParentCategory { + NoSpecialParent = 0, + LambdaParent, + SigmaParent, + K0ShortParent, + K0LongParent, + ChargedKaonParent, + HypertritonParent, + Lithium4Parent, + OtherDecayParent, + MissingDecayParent, + NParentCategories +}; + +enum PurityCategory { + AllSelected = 0, + CorrectSpecies, + MisidentifiedSpecies, + NoMCLabel, + CorrectCollisionSpecies, + WrongCollisionSpecies, + CorrectPrimary, + CorrectWeakDecay, + CorrectMaterial, + NPurityCategories +}; + } // namespace struct HadNucandidate { @@ -141,11 +196,20 @@ struct HadNucandidate { uint8_t nTPCClustersHad = 0u; uint8_t nTPCCrossedRowsNu = 0u; uint8_t nTPCCrossedRowsHad = 0u; + uint8_t nTPCFindableNu = 0u; + uint8_t nTPCFindableHad = 0u; uint8_t sharedClustersNu = 0u; uint8_t sharedClustersHad = 0u; + float fractionSharedTPCNu = 0.f; + float fractionSharedTPCHad = 0.f; float chi2TPCNu = -10.f; float chi2TPCHad = -10.f; + float chi2ITSNu = -10.f; + float chi2ITSHad = -10.f; float nSigmaNu = -10.f; + float nSigmaTPCNuDe = -10.f; + float nSigmaTPCNuPr = -10.f; + float nSigmaTPCNuPi = -10.f; float nSigmaHad = -10.f; float nSigmaTOFNu = -10.f; float nSigmaITSNu = -10.f; @@ -168,12 +232,16 @@ struct HadNucandidate { uint8_t nClsItsNu = 0u; uint8_t nClsItsHad = 0u; + uint8_t nClsItsInnerBarrelNu = 0u; + uint8_t nClsItsInnerBarrelHad = 0u; bool isBkgUS = false; // unlike sign bool isBkgEM = false; // event mixing - int trackIDNu = -1; - int trackIDHad = -1; + int64_t trackIDNu = -1; + int64_t trackIDHad = -1; + int64_t collisionIDNu = -1; + int64_t collisionIDHad = -1; float deltaEta = -99.f; float deltaPhi = -99.f; @@ -202,10 +270,17 @@ struct BufferedTrack { float tpcInnerParam{-99.f}; uint8_t tpcNClsFound{0u}; uint8_t tpcNClsCrossedRows{0u}; + uint8_t tpcNClsFindable{0u}; uint8_t tpcNClsShared{0u}; + float tpcFractionSharedCls{0.f}; uint8_t itsNCls{0u}; + uint8_t itsNClsInnerBarrel{0u}; float tpcChi2NCl{-10.f}; + float itsChi2NCl{-10.f}; float nSigmaTPC{-10.f}; + float nSigmaTPCDe{-10.f}; + float nSigmaTPCPr{-10.f}; + float nSigmaTPCPi{-10.f}; float nSigmaTOF{-10.f}; float nSigmaITS{-10.f}; float nSigmaTPCHadPi{-10.f}; @@ -218,6 +293,7 @@ struct BufferedTrack { uint32_t pidForTracking{0xFFFFFu}; uint32_t itsClusterSizes{0u}; int64_t trackId{-1}; + int64_t collisionId{-1}; }; struct BufferedCollision { @@ -252,6 +328,7 @@ struct HadNucleiFemto { // Event selection and mixing configuration Configurable settingCutVertex{"settingCutVertex", 10.0f, "Accepted z-vertex range"}; Configurable settingNoMixedEvents{"settingNoMixedEvents", 5, "Number of mixed events per event"}; + Configurable settingRequireBothSpeciesForMixing{"settingRequireBothSpeciesForMixing", false, "Use only events containing at least one selected nucleus and one selected hadron for mixed-event pairing"}; Configurable settingEnableBkgUS{"settingEnableBkgUS", false, "Enable US background"}; Configurable settingSaveUSandLS{"settingSaveUSandLS", true, "Save All Pairs"}; } eventMixing; @@ -338,6 +415,7 @@ struct HadNucleiFemto { std::string prefix{"CPR"}; // Close pair rejection controls Configurable settingEnableClosePairRejection{"settingEnableClosePairRejection", false, "Enable close pair rejection for nucleus-hadron track pairs"}; + Configurable settingApplyClosePairRejection{"settingApplyClosePairRejection", true, "Reject close pairs; disable to calculate and store CPR variables without rejecting pairs"}; Configurable settingClosePairDeltaPhiMax{"settingClosePairDeltaPhiMax", 0.01f, "Maximum delta phi star for close pair rejection"}; Configurable settingClosePairDeltaEtaMax{"settingClosePairDeltaEtaMax", 0.01f, "Maximum delta eta for close pair rejection"}; Configurable settingClosePairRadiusMode{"settingClosePairRadiusMode", 1, "Close pair rejection mode: 0 = PV, 1 = average phi star, 2 = specific TPC radius"}; @@ -354,6 +432,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 +686,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}, {NDcaOrigins, -0.5f, static_cast(NDcaOrigins) - 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}, {NDcaOrigins, -0.5f, static_cast(NDcaOrigins) - 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}, {NCollisionAssociations, -0.5f, static_cast(NCollisionAssociations) - 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}, {NCollisionAssociations, -0.5f, static_cast(NCollisionAssociations) - 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}, {NDcaOrigins, -0.5f, static_cast(NDcaOrigins) - 0.5f}, {NCollisionAssociations, -0.5f, static_cast(NCollisionAssociations) - 0.5f}, {NParentCategories, -0.5f, static_cast(NParentCategories) - 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}, {NDcaOrigins, -0.5f, static_cast(NDcaOrigins) - 0.5f}, {NCollisionAssociations, -0.5f, static_cast(NCollisionAssociations) - 0.5f}, {NParentCategories, -0.5f, static_cast(NParentCategories) - 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}, {NDcaOrigins, -0.5f, static_cast(NDcaOrigins) - 0.5f}, {NCollisionAssociations, -0.5f, static_cast(NCollisionAssociations) - 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}, {NDcaOrigins, -0.5f, static_cast(NDcaOrigins) - 0.5f}, {NCollisionAssociations, -0.5f, static_cast(NCollisionAssociations) - 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}, {NPurityCategories, -0.5f, static_cast(NPurityCategories) - 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}, {NPurityCategories, -0.5f, static_cast(NPurityCategories) - 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 +854,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 < NDcaOrigins; ++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 < NCollisionAssociations; ++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 < NDcaOrigins; ++i) { + histogram->GetAxis(4)->SetBinLabel(i + 1, originLabels[i]); + } + for (int i = 0; i < NCollisionAssociations; ++i) { + histogram->GetAxis(5)->SetBinLabel(i + 1, collisionLabels[i]); + } + for (int i = 0; i < NParentCategories; ++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 < NDcaOrigins; ++i) { + histogram->GetAxis(4)->SetBinLabel(i + 1, originLabels[i]); + } + for (int i = 0; i < NCollisionAssociations; ++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 < NPurityCategories; ++i) { + histogram->GetYaxis()->SetBinLabel(i + 1, purityLabels[i]); + } + } } template @@ -1054,6 +1232,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); @@ -1182,6 +1414,12 @@ struct HadNucleiFemto { if (fillQA) { mQaRegistry.fill(HIST("h2CPRBefore"), deltaEta, deltaPhi); } + if (!CPR.settingApplyClosePairRejection.value) { + if (fillQA) { + mQaRegistry.fill(HIST("h2CPRAfter"), deltaEta, deltaPhi); + } + return false; + } const bool isRejected = std::pow(deltaPhi, 2.f) / std::pow(CPR.settingClosePairDeltaPhiMax.value, 2.f) + std::pow(deltaEta, 2.f) / std::pow(CPR.settingClosePairDeltaEtaMax.value, 2.f) < 1.f; @@ -1618,6 +1856,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 Primary; + } + if (particle.getProcess() == TMCProcess::kPDecay) { + return WeakDecay; + } + // 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 Material; + } + + int purityOriginCategory(int origin) const + { + switch (origin) { + case Primary: + return CorrectPrimary; + case WeakDecay: + return CorrectWeakDecay; + case Material: + return CorrectMaterial; + default: + // classifyDcaOrigin deliberately has only the three Giorgio classes. + return CorrectMaterial; + } + } + + template + int classifyParentCategory(const Tparticle& particle, int origin) const + { + if (!particle.has_mothers()) { + return origin == WeakDecay ? MissingDecayParent : NoSpecialParent; + } + + const int fallbackCategory = origin == WeakDecay ? OtherDecayParent : NoSpecialParent; + 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 HypertritonParent; + } + if (motherPdg == Lithium4PDG) { + return Lithium4Parent; + } + + if (origin != WeakDecay) { + continue; + } + if (motherPdg == PDG_t::kLambda0) { + return LambdaParent; + } + if (motherPdg == PDG_t::kSigmaMinus || motherPdg == PDG_t::kSigma0 || motherPdg == PDG_t::kSigmaPlus) { + return SigmaParent; + } + if (motherPdg == PDG_t::kK0Short) { + return K0ShortParent; + } + if (motherPdg == PDG_t::kK0Long) { + return K0LongParent; + } + if (motherPdg == PDG_t::kKPlus) { + return ChargedKaonParent; + } + } + return fallbackCategory; + } + template float getNucleusTPCNSigma(const Ttrack& candidate) { @@ -1800,7 +2160,12 @@ struct HadNucleiFemto { hadNucand.nTPCClustersHad = trackHad.tpcNClsFound(); hadNucand.nTPCCrossedRowsNu = trackDe.tpcNClsCrossedRows(); hadNucand.nTPCCrossedRowsHad = trackHad.tpcNClsCrossedRows(); + hadNucand.nTPCFindableNu = trackDe.tpcNClsFindable(); + hadNucand.nTPCFindableHad = trackHad.tpcNClsFindable(); hadNucand.nSigmaNu = getNucleusTPCNSigma(trackDe); + hadNucand.nSigmaTPCNuDe = trackDe.tpcNSigmaDe(); + hadNucand.nSigmaTPCNuPr = trackDe.tpcNSigmaPr(); + hadNucand.nSigmaTPCNuPi = trackDe.tpcNSigmaPi(); hadNucand.nSigmaHad = getHadronTPCNSigma(trackHad); hadNucand.nSigmaTOFNu = getNucleusTOFNSigma(trackDe); hadNucand.nSigmaITSNu = getNucleusITSNSigma(trackDe); @@ -1815,6 +2180,8 @@ struct HadNucleiFemto { hadNucand.chi2TPCNu = trackDe.tpcChi2NCl(); hadNucand.chi2TPCHad = trackHad.tpcChi2NCl(); + hadNucand.chi2ITSNu = trackDe.itsChi2NCl(); + hadNucand.chi2ITSHad = trackHad.itsChi2NCl(); hadNucand.pidTrkNu = trackDe.pidForTracking(); hadNucand.pidTrkHad = trackHad.pidForTracking(); @@ -1824,9 +2191,13 @@ struct HadNucleiFemto { hadNucand.nClsItsNu = trackDe.itsNCls(); hadNucand.nClsItsHad = trackHad.itsNCls(); + hadNucand.nClsItsInnerBarrelNu = trackDe.itsNClsInnerBarrel(); + hadNucand.nClsItsInnerBarrelHad = trackHad.itsNClsInnerBarrel(); hadNucand.sharedClustersNu = trackDe.tpcNClsShared(); hadNucand.sharedClustersHad = trackHad.tpcNClsShared(); + hadNucand.fractionSharedTPCNu = trackDe.tpcFractionSharedCls(); + hadNucand.fractionSharedTPCHad = trackHad.tpcFractionSharedCls(); hadNucand.isBkgUS = trackDe.sign() * trackHad.sign() < 0; hadNucand.isBkgEM = isMixedEvent; @@ -1835,6 +2206,8 @@ struct HadNucleiFemto { hadNucand.trackIDNu = trackDe.globalIndex(); hadNucand.trackIDHad = trackHad.globalIndex(); + hadNucand.collisionIDNu = trackDe.collisionId(); + hadNucand.collisionIDHad = trackHad.collisionId(); if (trackDe.hasTOF()) { float beta = o2::pid::tof::Beta::GetBeta(trackDe); @@ -1926,7 +2299,7 @@ struct HadNucleiFemto { } template - BufferedTrack makeBufferedTrack(const Ttrack& track, bool isNucleus) + BufferedTrack makeBufferedTrack(const Ttrack& track, bool isNucleus, int64_t collisionId) { BufferedTrack bufferedTrack; bufferedTrack.momentum = {track.px(), track.py(), track.pz()}; @@ -1940,15 +2313,23 @@ struct HadNucleiFemto { bufferedTrack.tpcInnerParam = isNucleus && useHelium3Nucleus() ? correctedTPCInnerParamHe3(track) : track.tpcInnerParam(); bufferedTrack.tpcNClsFound = track.tpcNClsFound(); bufferedTrack.tpcNClsCrossedRows = track.tpcNClsCrossedRows(); + bufferedTrack.tpcNClsFindable = track.tpcNClsFindable(); bufferedTrack.tpcNClsShared = track.tpcNClsShared(); + bufferedTrack.tpcFractionSharedCls = track.tpcFractionSharedCls(); bufferedTrack.itsNCls = track.itsNCls(); + bufferedTrack.itsNClsInnerBarrel = track.itsNClsInnerBarrel(); bufferedTrack.tpcChi2NCl = track.tpcChi2NCl(); + bufferedTrack.itsChi2NCl = track.itsChi2NCl(); bufferedTrack.pidForTracking = track.pidForTracking(); bufferedTrack.itsClusterSizes = track.itsClusterSizes(); bufferedTrack.trackId = track.globalIndex(); + bufferedTrack.collisionId = collisionId; if (isNucleus) { bufferedTrack.nSigmaTPC = getNucleusTPCNSigma(track); + bufferedTrack.nSigmaTPCDe = track.tpcNSigmaDe(); + bufferedTrack.nSigmaTPCPr = track.tpcNSigmaPr(); + bufferedTrack.nSigmaTPCPi = track.tpcNSigmaPi(); bufferedTrack.nSigmaTOF = getNucleusTOFNSigma(track); bufferedTrack.nSigmaITS = getNucleusITSNSigma(track); } else { @@ -2012,13 +2393,24 @@ struct HadNucleiFemto { hadNucand.nTPCClustersHad = hadron.tpcNClsFound; hadNucand.nTPCCrossedRowsNu = nucleus.tpcNClsCrossedRows; hadNucand.nTPCCrossedRowsHad = hadron.tpcNClsCrossedRows; + hadNucand.nTPCFindableNu = nucleus.tpcNClsFindable; + hadNucand.nTPCFindableHad = hadron.tpcNClsFindable; hadNucand.sharedClustersNu = nucleus.tpcNClsShared; hadNucand.sharedClustersHad = hadron.tpcNClsShared; + hadNucand.fractionSharedTPCNu = nucleus.tpcFractionSharedCls; + hadNucand.fractionSharedTPCHad = hadron.tpcFractionSharedCls; hadNucand.nClsItsNu = nucleus.itsNCls; hadNucand.nClsItsHad = hadron.itsNCls; + hadNucand.nClsItsInnerBarrelNu = nucleus.itsNClsInnerBarrel; + hadNucand.nClsItsInnerBarrelHad = hadron.itsNClsInnerBarrel; hadNucand.chi2TPCNu = nucleus.tpcChi2NCl; hadNucand.chi2TPCHad = hadron.tpcChi2NCl; + hadNucand.chi2ITSNu = nucleus.itsChi2NCl; + hadNucand.chi2ITSHad = hadron.itsChi2NCl; hadNucand.nSigmaNu = nucleus.nSigmaTPC; + hadNucand.nSigmaTPCNuDe = nucleus.nSigmaTPCDe; + hadNucand.nSigmaTPCNuPr = nucleus.nSigmaTPCPr; + hadNucand.nSigmaTPCNuPi = nucleus.nSigmaTPCPi; hadNucand.nSigmaHad = hadron.nSigmaTPC; hadNucand.nSigmaTOFNu = nucleus.nSigmaTOF; hadNucand.nSigmaITSNu = nucleus.nSigmaITS; @@ -2036,8 +2428,10 @@ struct HadNucleiFemto { hadNucand.pidTrkHad = hadron.pidForTracking; hadNucand.itsClSizeNu = nucleus.itsClusterSizes; hadNucand.itsClSizeHad = hadron.itsClusterSizes; - hadNucand.trackIDNu = static_cast(nucleus.trackId); - hadNucand.trackIDHad = static_cast(hadron.trackId); + hadNucand.trackIDNu = nucleus.trackId; + hadNucand.trackIDHad = hadron.trackId; + hadNucand.collisionIDNu = nucleus.collisionId; + hadNucand.collisionIDHad = hadron.collisionId; hadNucand.isBkgUS = nucleus.sign() * hadron.sign() < 0; hadNucand.isBkgEM = true; @@ -2094,7 +2488,26 @@ struct HadNucleiFemto { hadNucand.nSigmaTOFNu, hadNucand.nSigmaITSNu, hadNucand.nSigmaTOFHad, - hadNucand.nSigmaITSHad); + hadNucand.nSigmaITSHad, + mRunNumber, + mDbz, + hadNucand.trackIDNu, + hadNucand.trackIDHad, + hadNucand.collisionIDNu, + hadNucand.collisionIDHad, + hadNucand.nTPCFindableNu, + hadNucand.nTPCFindableHad, + hadNucand.fractionSharedTPCNu, + hadNucand.fractionSharedTPCHad, + hadNucand.nClsItsNu, + hadNucand.nClsItsHad, + hadNucand.nClsItsInnerBarrelNu, + hadNucand.nClsItsInnerBarrelHad, + hadNucand.chi2ITSNu, + hadNucand.chi2ITSHad, + hadNucand.nSigmaTPCNuDe, + hadNucand.nSigmaTPCNuPr, + hadNucand.nSigmaTPCNuPi); } template @@ -2838,16 +3251,21 @@ struct HadNucleiFemto { tracksThisCollision.bindExternalIndices(&tracks); for (const auto& track : tracksThisCollision) { if (selectTrackNu(track) && selectionPIDNu(track)) { - currentCollision.nuclei.push_back(makeBufferedTrack(track, /*isNucleus*/ true)); + currentCollision.nuclei.push_back(makeBufferedTrack(track, /*isNucleus*/ true, currentCollision.eventId)); } if (selectTrackHadron(track) && selectionPIDHadron(track)) { - currentCollision.hadrons.push_back(makeBufferedTrack(track, /*isNucleus*/ false)); + currentCollision.hadrons.push_back(makeBufferedTrack(track, /*isNucleus*/ false, currentCollision.eventId)); } } mQaRegistry.fill(HIST("hMixedNucleiPerEvent"), currentCollision.nuclei.size()); mQaRegistry.fill(HIST("hMixedHadronsPerEvent"), currentCollision.hadrons.size()); + if (eventMixing.settingRequireBothSpeciesForMixing.value && + (currentCollision.nuclei.empty() || currentCollision.hadrons.empty())) { + continue; + } + const int poolBin = configuredBinningPolicy.getBin(std::make_tuple(collision.posZ(), collision.centFT0C())); auto& pool = mMixingPools[poolBin]; mQaRegistry.fill(HIST("hMixingPoolOccupancy"), pool.size()); @@ -2882,6 +3300,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(CorrectCollision) : static_cast(WrongCollision); + 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 ? CorrectCollision : WrongCollision; + 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(AllSelected); + if (!correctSpecies) { + fillCategory(MisidentifiedSpecies); + return; + } + fillCategory(CorrectSpecies); + if (!matchesRecoCollision) { + fillCategory(WrongCollisionSpecies); + return; + } + fillCategory(CorrectCollisionSpecies); + 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(AllSelected), centrality); + mQaRegistry.fill(HIST("purityMC/hNucleus"), signedPt, static_cast(NoMCLabel), centrality); + } else { + mQaRegistry.fill(HIST("purityMC/hHadron"), signedPt, static_cast(AllSelected), centrality); + mQaRegistry.fill(HIST("purityMC/hHadron"), signedPt, static_cast(NoMCLabel), 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) {