From cc29826cf262af7f712d3219f6ecfcbb34fe5d94 Mon Sep 17 00:00:00 2001 From: morgmatt Date: Thu, 31 Jul 2025 10:14:59 +0200 Subject: [PATCH 01/18] Added sigma plus to MC selection criteria, corrected name for rapidity cut --- PWGLF/TableProducer/Strangeness/sigmaminustask.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx index fbf5929effc..a6eb227ee02 100644 --- a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx @@ -45,7 +45,7 @@ struct sigmaminustask { // Configurable for event selection Configurable cutzvertex{"cutzvertex", 10.0f, "Accepted z-vertex range (cm)"}; Configurable cutNSigmaPi{"cutNSigmaPi", 4, "NSigmaTPCPion"}; - Configurable cutEtaMotherMC{"cutEtaMotherMC", 1.0f, "Eta cut for mother Sigma in MC"}; + Configurable cutRapMotherMC{"cutRapMotherMC", 1.0f, "Rapidity cut for mother Sigma in MC"}; Configurable fillOutputTree{"fillOutputTree", true, "If true, fill the output tree with Kink candidates"}; @@ -153,7 +153,7 @@ struct sigmaminustask { if (piMother.globalIndex() != mcTrackSigma.globalIndex()) { continue; } - if (std::abs(mcTrackSigma.pdgCode()) != 3112) { + if (std::abs(mcTrackSigma.pdgCode()) != 3112 && std::abs(mcTrackSigma.pdgCode()) != 3222) { continue; } if (std::abs(mcTrackPiDau.pdgCode()) != 211 && std::abs(mcTrackPiDau.pdgCode()) != 2212) { @@ -200,7 +200,7 @@ struct sigmaminustask { // Loop over all generated particles to fill MC histograms for (const auto& mcPart : particlesMC) { - if (std::abs(mcPart.pdgCode()) != 3112 || std::abs(mcPart.y()) > cutEtaMotherMC) { // only sigma mothers and rapidity cut + if ((std::abs(mcPart.pdgCode()) != 3112 && std::abs(mcPart.pdgCode()) != 3222) || std::abs(mcPart.y()) > cutRapMotherMC) { // only sigma mothers and rapidity cut continue; } if (!mcPart.has_daughters()) { From 88109a0beda636097391a306d433adba42acbe52 Mon Sep 17 00:00:00 2001 From: morgmatt Date: Fri, 1 Aug 2025 16:58:33 +0200 Subject: [PATCH 02/18] Added control histograms to the task: nSigmaTOF for reconstructed particles, BC Id checks, mcCollision Id checks --- .../Strangeness/sigmaminustask.cxx | 51 +++++++++++++------ 1 file changed, 36 insertions(+), 15 deletions(-) diff --git a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx index a6eb227ee02..90a811da461 100644 --- a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx @@ -55,7 +55,8 @@ struct sigmaminustask { { // Axes const AxisSpec ptAxis{100, -10, 10, "#it{p}_{T} (GeV/#it{c})"}; - const AxisSpec nSigmaPiAxis{100, -5, 5, "n#sigma_{#pi}"}; + const AxisSpec nSigmaPiAxis{1000, -1000, 1000, "n#sigma_{#pi}"}; + const AxisSpec nSigmaPrAxis{1000, -1000, 1000, "n#sigma_{p}"}; const AxisSpec sigmaMassAxis{100, 1.1, 1.4, "m (GeV/#it{c}^{2})"}; const AxisSpec xiMassAxis{100, 1.2, 1.6, "m_{#Xi} (GeV/#it{c}^{2})"}; const AxisSpec pdgAxis{10001, -5000, 5000, "PDG code"}; @@ -63,23 +64,30 @@ struct sigmaminustask { const AxisSpec ptResolutionAxis{100, -0.5, 0.5, "#it{p}_{T}^{rec} - #it{p}_{T}^{gen} (GeV/#it{c})"}; const AxisSpec massResolutionAxis{100, -0.1, 0.1, "m_{rec} - m_{gen} (GeV/#it{c}^{2})"}; - + const AxisSpec boolAxis{2, -0.5, 1.5, "Boolean value (0=false, 1=true)"}; // Event selection rEventSelection.add("hVertexZRec", "hVertexZRec", {HistType::kTH1F, {vertexZAxis}}); // Sigma-minus reconstruction rSigmaMinus.add("h2MassSigmaMinusPt", "h2MassSigmaMinusPt", {HistType::kTH2F, {ptAxis, sigmaMassAxis}}); rSigmaMinus.add("h2SigmaMassVsXiMass", "h2SigmaMassVsXiMass", {HistType::kTH2F, {xiMassAxis, sigmaMassAxis}}); - rSigmaMinus.add("h2NSigmaPiPt", "h2NSigmaPiPt", {HistType::kTH2F, {ptAxis, nSigmaPiAxis}}); + rSigmaMinus.add("h2NSigmaTPCPiPt", "h2NSigmaTPCPiPt", {HistType::kTH2F, {ptAxis, nSigmaPiAxis}}); if (doprocessMC) { // Add MC histograms if needed rSigmaMinus.add("h2MassPtMCRec", "h2MassPtMCRec", {HistType::kTH2F, {ptAxis, sigmaMassAxis}}); rSigmaMinus.add("h2MassPtMCGen", "h2MassPtMCGen", {HistType::kTH2F, {ptAxis, sigmaMassAxis}}); - rSigmaMinus.add("h2MassResolution_minus", "h2MassResolution_minus", {HistType::kTH2F, {sigmaMassAxis, massResolutionAxis}}); - rSigmaMinus.add("h2PtResolution_minus", "h2PtResolution_minus", {HistType::kTH2F, {ptAxis, ptResolutionAxis}}); - rSigmaMinus.add("h2MassResolution_plus", "h2MassResolution_plus", {HistType::kTH2F, {sigmaMassAxis, massResolutionAxis}}); - rSigmaMinus.add("h2PtResolution_plus", "h2PtResolution_plus", {HistType::kTH2F, {ptAxis, ptResolutionAxis}}); + rSigmaMinus.add("h2MassResolution", "h2MassResolution", {HistType::kTH2F, {ptAxis, massResolutionAxis}}); + rSigmaMinus.add("h2PtResolution", "h2PtResolution", {HistType::kTH2F, {ptAxis, ptResolutionAxis}}); + + rSigmaMinus.add("h2NSigmaTOFPiPt", "h2NSigmaTOFPiPt", {HistType::kTH2F, {ptAxis, nSigmaPiAxis}}); + rSigmaMinus.add("h2NSigmaTOFPrPt", "h2NSigmaTOFPrPt", {HistType::kTH2F, {ptAxis, nSigmaPrAxis}}); + + // BC ID comparison histograms + rSigmaMinus.add("hMcCollIdCoherence", "McCollId (coll == daug)", {HistType::kTH1F, {boolAxis}}); + rSigmaMinus.add("h2CollId_BCId", "(McCollId coherence) vs (EvSelBC == McBC)", {HistType::kTH2F, {boolAxis, boolAxis}}); + rSigmaMinus.add("h2BCId_comp1", "(BC == McBC) vs (BC == EvSelBC)", {HistType::kTH2F, {boolAxis, boolAxis}}); + rSigmaMinus.add("h2BCId_comp2", "(McBC == EvSelBC) vs (BC == EvSelBC)", {HistType::kTH2F, {boolAxis, boolAxis}}); } } @@ -114,7 +122,7 @@ struct sigmaminustask { } PROCESS_SWITCH(sigmaminustask, processData, "Data processing", true); - void processMC(CollisionsFullMC const& collisions, aod::KinkCands const& KinkCands, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const& particlesMC, aod::McCollisions const&, TracksFull const&) + void processMC(CollisionsFullMC const& collisions, aod::KinkCands const& KinkCands, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const& particlesMC, TracksFull const&, aod::McCollisions const&) { for (const auto& collision : collisions) { if (std::abs(collision.posZ()) > cutzvertex || !collision.sel8()) { @@ -136,9 +144,10 @@ struct sigmaminustask { continue; } + // histograms filled with all kink candidates rSigmaMinus.fill(HIST("h2MassSigmaMinusPt"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.mSigmaMinus()); rSigmaMinus.fill(HIST("h2SigmaMassVsXiMass"), kinkCand.mXiMinus(), kinkCand.mSigmaMinus()); - rSigmaMinus.fill(HIST("h2NSigmaPiPt"), kinkCand.mothSign() * kinkCand.ptMoth(), dauTrack.tpcNSigmaPi()); + rSigmaMinus.fill(HIST("h2NSigmaTPCPiPt"), kinkCand.mothSign() * kinkCand.ptMoth(), dauTrack.tpcNSigmaPi()); // do MC association auto mcLabSigma = trackLabelsMC.rawIteratorAt(mothTrack.globalIndex()); @@ -171,14 +180,25 @@ struct sigmaminustask { if (collision.has_mcCollision()) { mcCollisionIdCheck = collision.mcCollision().globalIndex() == mcTrackPiDau.mcCollisionId(); } + // Check bunch crossing ID coherence + auto mcCollision = mcTrackPiDau.template mcCollision_as(); + bool BCId_vs_MCBCId = collision.bcId() == mcCollision.bcId(); + bool BCId_vs_EvSel = collision.bcId() == collision.foundBCId(); + bool EvSel_vs_MCBCId = collision.foundBCId() == mcCollision.bcId(); + + rSigmaMinus.fill(HIST("hMcCollIdCoherence"), static_cast(mcCollisionIdCheck)); + rSigmaMinus.fill(HIST("h2CollId_BCId"), static_cast(mcCollisionIdCheck), static_cast(EvSel_vs_MCBCId)); + rSigmaMinus.fill(HIST("h2BCId_comp1"), static_cast(BCId_vs_MCBCId), static_cast(BCId_vs_EvSel)); + rSigmaMinus.fill(HIST("h2BCId_comp2"), static_cast(EvSel_vs_MCBCId), static_cast(BCId_vs_EvSel)); rSigmaMinus.fill(HIST("h2MassPtMCRec"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.mSigmaMinus()); - if (mcTrackSigma.pdgCode() > 0) { - rSigmaMinus.fill(HIST("h2MassResolution_plus"), kinkCand.mSigmaMinus(), kinkCand.mSigmaMinus() - MotherMassMC); - rSigmaMinus.fill(HIST("h2PtResolution_plus"), kinkCand.ptMoth(), kinkCand.ptMoth() - MotherpTMC); - } else { - rSigmaMinus.fill(HIST("h2MassResolution_minus"), kinkCand.mSigmaMinus(), kinkCand.mSigmaMinus() - MotherMassMC); - rSigmaMinus.fill(HIST("h2PtResolution_minus"), kinkCand.ptMoth(), kinkCand.ptMoth() - MotherpTMC); + rSigmaMinus.fill(HIST("h2MassResolution"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.mSigmaMinus() - MotherMassMC); + rSigmaMinus.fill(HIST("h2PtResolution"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.ptMoth() - MotherpTMC); + + if (std::abs(mcTrackPiDau.pdgCode()) == 211) { + rSigmaMinus.fill(HIST("h2NSigmaTOFPiPt"), kinkCand.mothSign() * kinkCand.ptMoth(), dauTrack.tofNSigmaPi()); + } else if (std::abs(mcTrackPiDau.pdgCode()) == 2212) { + rSigmaMinus.fill(HIST("h2NSigmaTOFPrPt"), kinkCand.mothSign() * kinkCand.ptMoth(), dauTrack.tofNSigmaPr()); } // fill the output table with Mc information @@ -250,3 +270,4 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) return WorkflowSpec{ adaptAnalysisTask(cfgc)}; } + From 5cfa04643891f61e2871993bedc6f5cff37746ce Mon Sep 17 00:00:00 2001 From: morgmatt Date: Fri, 1 Aug 2025 17:00:15 +0200 Subject: [PATCH 03/18] Added processFindable function to determine generated particles that are findable by the kinkbuilder --- .../Strangeness/sigmaminustask.cxx | 64 +++++++++++++++++++ 1 file changed, 64 insertions(+) diff --git a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx index 90a811da461..246517f8456 100644 --- a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx @@ -89,6 +89,11 @@ struct sigmaminustask { rSigmaMinus.add("h2BCId_comp1", "(BC == McBC) vs (BC == EvSelBC)", {HistType::kTH2F, {boolAxis, boolAxis}}); rSigmaMinus.add("h2BCId_comp2", "(McBC == EvSelBC) vs (BC == EvSelBC)", {HistType::kTH2F, {boolAxis, boolAxis}}); } + + if (doprocessFindable) { + // Add findable Sigma histograms + rSigmaMinus.add("h2MassPtFindable", "h2MassPtFindable", {HistType::kTH2F, {ptAxis, sigmaMassAxis}}); + } } void processData(CollisionsFull::iterator const& collision, aod::KinkCands const& KinkCands, TracksFull const&) @@ -263,6 +268,65 @@ struct sigmaminustask { } PROCESS_SWITCH(sigmaminustask, processMC, "MC processing", false); + + void processFindable(TracksFull const& tracks, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const&) + { + for (const auto& motherTrack : tracks) { + + // Filter: ITS but no TPC (mother candidates) + if (!motherTrack.hasITS() || motherTrack.hasTPC()) { + continue; + } + + // Check if mother is Sigma in MC + auto mcLabelMother = trackLabelsMC.rawIteratorAt(motherTrack.globalIndex()); + if (!mcLabelMother.has_mcParticle()) { + continue; + } + auto mcMother = mcLabelMother.mcParticle_as(); + if (std::abs(mcMother.pdgCode()) != 3112 && std::abs(mcMother.pdgCode()) != 3222) { + continue; + } + + for (const auto& daughterTrack : tracks) { + // Filter: ITS and TPC (daughter candidates) + if (!daughterTrack.hasITS() || !daughterTrack.hasTPC()) { + continue; + } + + auto mcLabelDaughter = trackLabelsMC.rawIteratorAt(daughterTrack.globalIndex()); + if (!mcLabelDaughter.has_mcParticle()) { + continue; + } + + // Check if daughter is pi/proton + auto mcDaughter = mcLabelDaughter.mcParticle_as(); + if (std::abs(mcDaughter.pdgCode()) != 211 && std::abs(mcDaughter.pdgCode()) != 2212) { + continue; + } + + // Verify the MC mother-daughter relationship + bool isValidPair = false; + if (mcDaughter.has_mothers()) { + for (const auto& mother : mcDaughter.mothers_as()) { + if (mother.globalIndex() == mcMother.globalIndex()) { + isValidPair = true; + break; + } + } + } + + if (isValidPair) { + // All requirements satisfied: fill histogram + float mcMass = std::sqrt(mcMother.e() * mcMother.e() - mcMother.p() * mcMother.p()); + int sigmaSign = mcMother.pdgCode() > 0 ? 1 : -1; + rSigmaMinus.fill(HIST("h2MassPtFindable"), sigmaSign * mcMother.pt(), mcMass); + } + } + } + } + + PROCESS_SWITCH(sigmaminustask, processFindable, "Findable Sigma processing", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) From 0e48b78fc3836c7cc70f33ff6415e05331ce1c25 Mon Sep 17 00:00:00 2001 From: morgmatt Date: Mon, 4 Aug 2025 11:54:55 +0200 Subject: [PATCH 04/18] Added DCA histograms, started implementation of the efficiency analysis for the kinkbuilder algorithm in processFindable function --- .../Strangeness/sigmaminustask.cxx | 88 ++++++++++++++----- 1 file changed, 66 insertions(+), 22 deletions(-) diff --git a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx index 246517f8456..013f72b77b7 100644 --- a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx @@ -55,22 +55,29 @@ struct sigmaminustask { { // Axes const AxisSpec ptAxis{100, -10, 10, "#it{p}_{T} (GeV/#it{c})"}; - const AxisSpec nSigmaPiAxis{1000, -1000, 1000, "n#sigma_{#pi}"}; - const AxisSpec nSigmaPrAxis{1000, -1000, 1000, "n#sigma_{p}"}; + const AxisSpec nSigmaPiAxis{60, -30, 30, "n#sigma_{#pi}"}; + const AxisSpec nSigmaPrAxis{60, -30, 30, "n#sigma_{p}"}; const AxisSpec sigmaMassAxis{100, 1.1, 1.4, "m (GeV/#it{c}^{2})"}; const AxisSpec xiMassAxis{100, 1.2, 1.6, "m_{#Xi} (GeV/#it{c}^{2})"}; const AxisSpec pdgAxis{10001, -5000, 5000, "PDG code"}; const AxisSpec vertexZAxis{100, -15., 15., "vrtx_{Z} [cm]"}; + const AxisSpec dcaMothAxis{100, 0, 1, "DCA [cm]"}; + const AxisSpec dcaDaugAxis{200, 0, 20, "DCA [cm]"}; const AxisSpec ptResolutionAxis{100, -0.5, 0.5, "#it{p}_{T}^{rec} - #it{p}_{T}^{gen} (GeV/#it{c})"}; const AxisSpec massResolutionAxis{100, -0.1, 0.1, "m_{rec} - m_{gen} (GeV/#it{c}^{2})"}; + const AxisSpec boolAxis{2, -0.5, 1.5, "Boolean value (0=false, 1=true)"}; + const AxisSpec filtersAxis{10, -0.5, 9.5, "Filter index"}; + // Event selection rEventSelection.add("hVertexZRec", "hVertexZRec", {HistType::kTH1F, {vertexZAxis}}); // Sigma-minus reconstruction rSigmaMinus.add("h2MassSigmaMinusPt", "h2MassSigmaMinusPt", {HistType::kTH2F, {ptAxis, sigmaMassAxis}}); rSigmaMinus.add("h2SigmaMassVsXiMass", "h2SigmaMassVsXiMass", {HistType::kTH2F, {xiMassAxis, sigmaMassAxis}}); rSigmaMinus.add("h2NSigmaTPCPiPt", "h2NSigmaTPCPiPt", {HistType::kTH2F, {ptAxis, nSigmaPiAxis}}); + rSigmaMinus.add("h2DCAMothPt", "h2DCAMothPt", {HistType::kTH2F, {ptAxis, dcaMothAxis}}); + rSigmaMinus.add("h2DCADaugPt", "h2DCADaugPt", {HistType::kTH2F, {ptAxis, dcaDaugAxis}}); if (doprocessMC) { // Add MC histograms if needed @@ -93,6 +100,7 @@ struct sigmaminustask { if (doprocessFindable) { // Add findable Sigma histograms rSigmaMinus.add("h2MassPtFindable", "h2MassPtFindable", {HistType::kTH2F, {ptAxis, sigmaMassAxis}}); + rSigmaMinus.add("hFilterIndex", "hFilterIndex", {HistType::kTH1F, {filtersAxis}}); } } @@ -113,6 +121,8 @@ struct sigmaminustask { rSigmaMinus.fill(HIST("h2MassSigmaMinusPt"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.mSigmaMinus()); rSigmaMinus.fill(HIST("h2SigmaMassVsXiMass"), kinkCand.mXiMinus(), kinkCand.mSigmaMinus()); rSigmaMinus.fill(HIST("h2NSigmaPiPt"), kinkCand.mothSign() * kinkCand.ptMoth(), dauTrack.tpcNSigmaPi()); + rSigmaMinus.fill(HIST("h2DCAMothPt"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.dcaMothPv()); + rSigmaMinus.fill(HIST("h2DCADaugPt"), kinkCand.mothSign() * kinkCand.ptDaug(), kinkCand.dcaDaugPv()); if (fillOutputTree) { outputDataTable(kinkCand.xDecVtx(), kinkCand.yDecVtx(), kinkCand.zDecVtx(), @@ -199,7 +209,9 @@ struct sigmaminustask { rSigmaMinus.fill(HIST("h2MassPtMCRec"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.mSigmaMinus()); rSigmaMinus.fill(HIST("h2MassResolution"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.mSigmaMinus() - MotherMassMC); rSigmaMinus.fill(HIST("h2PtResolution"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.ptMoth() - MotherpTMC); - + rSigmaMinus.fill(HIST("h2DCAMothPt"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.dcaMothPv()); + rSigmaMinus.fill(HIST("h2DCADaugPt"), kinkCand.mothSign() * kinkCand.ptDaug(), kinkCand.dcaDaugPv()); + if (std::abs(mcTrackPiDau.pdgCode()) == 211) { rSigmaMinus.fill(HIST("h2NSigmaTOFPiPt"), kinkCand.mothSign() * kinkCand.ptMoth(), dauTrack.tofNSigmaPi()); } else if (std::abs(mcTrackPiDau.pdgCode()) == 2212) { @@ -272,12 +284,6 @@ struct sigmaminustask { void processFindable(TracksFull const& tracks, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const&) { for (const auto& motherTrack : tracks) { - - // Filter: ITS but no TPC (mother candidates) - if (!motherTrack.hasITS() || motherTrack.hasTPC()) { - continue; - } - // Check if mother is Sigma in MC auto mcLabelMother = trackLabelsMC.rawIteratorAt(motherTrack.globalIndex()); if (!mcLabelMother.has_mcParticle()) { @@ -289,17 +295,11 @@ struct sigmaminustask { } for (const auto& daughterTrack : tracks) { - // Filter: ITS and TPC (daughter candidates) - if (!daughterTrack.hasITS() || !daughterTrack.hasTPC()) { - continue; - } - + // Check if daughter is pi/proton auto mcLabelDaughter = trackLabelsMC.rawIteratorAt(daughterTrack.globalIndex()); if (!mcLabelDaughter.has_mcParticle()) { continue; } - - // Check if daughter is pi/proton auto mcDaughter = mcLabelDaughter.mcParticle_as(); if (std::abs(mcDaughter.pdgCode()) != 211 && std::abs(mcDaughter.pdgCode()) != 2212) { continue; @@ -315,13 +315,58 @@ struct sigmaminustask { } } } + if (!isValidPair) { + continue; + } + + float mcMass = std::sqrt(mcMother.e() * mcMother.e() - mcMother.p() * mcMother.p()); + int sigmaSign = mcMother.pdgCode() > 0 ? 1 : -1; + rSigmaMinus.fill(HIST("h2MassPtFindable"), sigmaSign * mcMother.pt(), mcMass); + + // Define filter index and progressively apply kinkbuilder cuts to track pairs + int filterIndex = 0; + rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex); + + // 1 - tracks with right ITS, TPC, TOF signals + if (motherTrack.has_collision() && motherTrack.hasITS() && !motherTrack.hasTPC() && !motherTrack.hasTOF() && + daughterTrack.hasITS() && daughterTrack.hasTPC()) { + filterIndex += 1; + rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex); + } else { + continue; + } - if (isValidPair) { - // All requirements satisfied: fill histogram - float mcMass = std::sqrt(mcMother.e() * mcMother.e() - mcMother.p() * mcMother.p()); - int sigmaSign = mcMother.pdgCode() > 0 ? 1 : -1; - rSigmaMinus.fill(HIST("h2MassPtFindable"), sigmaSign * mcMother.pt(), mcMass); + // 2, 3 - mother track ITS properties + if (motherTrack.itsNCls() < 6 && + motherTrack.itsNClsInnerBarrel() == 3 && motherTrack.itsChi2NCl() < 36) { + filterIndex += 1; + rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex); + } + if (filterIndex > 1 && motherTrack.pt() > 0.5) { + filterIndex += 1; + rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex); + } else { + continue; // Skip if mother track does not pass ITS cuts } + + // 4 - daughter track ITS+TPC properties + if (daughterTrack.itsNClsInnerBarrel() == 0 && daughterTrack.itsNCls() < 4 && + daughterTrack.tpcNClsCrossedRows() > 0.8 * daughterTrack.tpcNClsFindable() && daughterTrack.tpcNClsFound() > 80) { + filterIndex += 1; + rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex); + } else { + continue; // Skip if daughter track does not pass ITS+TPC cuts + } + + // 5 - geometric cuts + if (std::abs(motherTrack.eta()) < 1.0 && std::abs(daughterTrack.eta()) < 1.0 && + std::abs(motherTrack.phi() - daughterTrack.phi()) < 100.0) { + filterIndex += 1; + rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex); + } else { + continue; // Skip if geometric cuts are not satisfied + } + } } } @@ -334,4 +379,3 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) return WorkflowSpec{ adaptAnalysisTask(cfgc)}; } - From dbd5cd5f7425fe7430db9ecf48c0e647e2dbe329 Mon Sep 17 00:00:00 2001 From: morgmatt Date: Mon, 4 Aug 2025 15:17:06 +0200 Subject: [PATCH 05/18] Added geometric cuts to processFindable analysis --- .../Strangeness/sigmaminustask.cxx | 31 ++++++++++++++----- 1 file changed, 23 insertions(+), 8 deletions(-) diff --git a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx index 013f72b77b7..6e176f8a588 100644 --- a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx @@ -51,6 +51,8 @@ struct sigmaminustask { Preslice mPerCol = aod::track::collisionId; + float radToDeg = o2::constants::math::Rad2Deg; + void init(InitContext const&) { // Axes @@ -281,7 +283,7 @@ struct sigmaminustask { PROCESS_SWITCH(sigmaminustask, processMC, "MC processing", false); - void processFindable(TracksFull const& tracks, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const&) + void processFindable(TracksFull const& tracks, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const&, CollisionsFullMC const&) { for (const auto& motherTrack : tracks) { // Check if mother is Sigma in MC @@ -321,7 +323,7 @@ struct sigmaminustask { float mcMass = std::sqrt(mcMother.e() * mcMother.e() - mcMother.p() * mcMother.p()); int sigmaSign = mcMother.pdgCode() > 0 ? 1 : -1; - rSigmaMinus.fill(HIST("h2MassPtFindable"), sigmaSign * mcMother.pt(), mcMass); + rSigmaMinus.fill(HIST("h2MassPtFindable"), sigmaSign * motherTrack.pt(), mcMass); // Define filter index and progressively apply kinkbuilder cuts to track pairs int filterIndex = 0; @@ -346,7 +348,7 @@ struct sigmaminustask { filterIndex += 1; rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex); } else { - continue; // Skip if mother track does not pass ITS cuts + continue; } // 4 - daughter track ITS+TPC properties @@ -355,16 +357,29 @@ struct sigmaminustask { filterIndex += 1; rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex); } else { - continue; // Skip if daughter track does not pass ITS+TPC cuts + continue; } - // 5 - geometric cuts - if (std::abs(motherTrack.eta()) < 1.0 && std::abs(daughterTrack.eta()) < 1.0 && - std::abs(motherTrack.phi() - daughterTrack.phi()) < 100.0) { + // 5 - geometric cuts: eta + if (std::abs(motherTrack.eta()) < 1.0 && std::abs(daughterTrack.eta()) < 1.0) { filterIndex += 1; rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex); } else { - continue; // Skip if geometric cuts are not satisfied + continue; + } + // 6 - geometric cuts: phi + if (std::abs(motherTrack.phi() - daughterTrack.phi()) * radToDeg < 100.0) { + filterIndex += 1; + rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex); + } else { + continue; + } + + // 7 - collision selection + auto collision = motherTrack.template collision_as(); + if (!(std::abs(collision.posZ()) > cutzvertex || !collision.sel8())) { + filterIndex += 1; + rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex); } } From d8f1e5781e3487c574f856acbb5b2cdb618d5795 Mon Sep 17 00:00:00 2001 From: morgmatt Date: Mon, 4 Aug 2025 18:48:12 +0200 Subject: [PATCH 06/18] Added histograms for decay radii, created new registry for findable histograms --- .../Strangeness/sigmaminustask.cxx | 51 ++++++++++++++----- 1 file changed, 39 insertions(+), 12 deletions(-) diff --git a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx index 6e176f8a588..660b2ffc29f 100644 --- a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx @@ -18,9 +18,11 @@ #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/PIDResponse.h" +#include "Common/Core/trackUtilities.h" #include "Framework/AnalysisTask.h" #include "Framework/runDataProcessing.h" #include "ReconstructionDataFormats/PID.h" +#include "ReconstructionDataFormats/Track.h" using namespace o2; using namespace o2::framework; @@ -41,6 +43,7 @@ struct sigmaminustask { // Histograms are defined with HistogramRegistry HistogramRegistry rEventSelection{"eventSelection", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; HistogramRegistry rSigmaMinus{"sigmaminus", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + HistogramRegistry rFindable{"findable", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; // Configurable for event selection Configurable cutzvertex{"cutzvertex", 10.0f, "Accepted z-vertex range (cm)"}; @@ -65,9 +68,11 @@ struct sigmaminustask { const AxisSpec vertexZAxis{100, -15., 15., "vrtx_{Z} [cm]"}; const AxisSpec dcaMothAxis{100, 0, 1, "DCA [cm]"}; const AxisSpec dcaDaugAxis{200, 0, 20, "DCA [cm]"}; + const AxisSpec radiusAxis{100, -1, 40, "Decay radius [cm]"}; - const AxisSpec ptResolutionAxis{100, -0.5, 0.5, "#it{p}_{T}^{rec} - #it{p}_{T}^{gen} (GeV/#it{c})"}; - const AxisSpec massResolutionAxis{100, -0.1, 0.1, "m_{rec} - m_{gen} (GeV/#it{c}^{2})"}; + const AxisSpec ptResolutionAxis{100, -2, 2, "(#it{p}_{T}^{rec} - #it{p}_{T}^{gen}) / #it{p}_{T}^{gen}"}; + const AxisSpec massResolutionAxis{100, -2, 2, "(m_{rec} - m_{gen}) / m_{gen}"}; + const AxisSpec radiusResolutionAxis{100, -2, 2, "(r_{rec} - r_{gen}) / r_{gen}"}; const AxisSpec boolAxis{2, -0.5, 1.5, "Boolean value (0=false, 1=true)"}; const AxisSpec filtersAxis{10, -0.5, 9.5, "Filter index"}; @@ -101,8 +106,11 @@ struct sigmaminustask { if (doprocessFindable) { // Add findable Sigma histograms - rSigmaMinus.add("h2MassPtFindable", "h2MassPtFindable", {HistType::kTH2F, {ptAxis, sigmaMassAxis}}); - rSigmaMinus.add("hFilterIndex", "hFilterIndex", {HistType::kTH1F, {filtersAxis}}); + rFindable.add("h2MassPtFindableAll", "h2MassPtFindableAll", {HistType::kTH2F, {ptAxis, sigmaMassAxis}}); + rFindable.add("hFilterIndex", "hFilterIndex", {HistType::kTH1F, {filtersAxis}}); + rFindable.add("h2MCRadiusFilterIndex", "h2RadiusFilterIndex", {HistType::kTH2F, {filtersAxis, radiusAxis}}); + rFindable.add("h2RecRadiusFilterIndex", "h2RecRadiusFilterIndex", {HistType::kTH2F, {filtersAxis, radiusAxis}}); + } } @@ -191,6 +199,7 @@ struct sigmaminustask { float deltaXMother = mcTrackPiDau.vx() - piMother.vx(); float deltaYMother = mcTrackPiDau.vy() - piMother.vy(); float decayRadiusMC = std::sqrt(deltaXMother * deltaXMother + deltaYMother * deltaYMother); + float decayRadiusRec = std::sqrt(kinkCand.xDecVtx() * kinkCand.xDecVtx() + kinkCand.yDecVtx() * kinkCand.yDecVtx()); // Check coherence of MCcollision Id for daughter MCparticle and reconstructed collision bool mcCollisionIdCheck = false; @@ -209,10 +218,11 @@ struct sigmaminustask { rSigmaMinus.fill(HIST("h2BCId_comp2"), static_cast(EvSel_vs_MCBCId), static_cast(BCId_vs_EvSel)); rSigmaMinus.fill(HIST("h2MassPtMCRec"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.mSigmaMinus()); - rSigmaMinus.fill(HIST("h2MassResolution"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.mSigmaMinus() - MotherMassMC); - rSigmaMinus.fill(HIST("h2PtResolution"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.ptMoth() - MotherpTMC); + rSigmaMinus.fill(HIST("h2MassResolution"), kinkCand.mothSign() * kinkCand.ptMoth(), (kinkCand.mSigmaMinus() - MotherMassMC) / MotherMassMC); + rSigmaMinus.fill(HIST("h2PtResolution"), kinkCand.mothSign() * kinkCand.ptMoth(), (kinkCand.ptMoth() - MotherpTMC) / MotherpTMC); + rSigmaMinus.fill(HIST("h2RadiusResolution"), kinkCand.mothSign() * kinkCand.ptMoth(), (decayRadiusRec - decayRadiusMC) / decayRadiusMC); rSigmaMinus.fill(HIST("h2DCAMothPt"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.dcaMothPv()); - rSigmaMinus.fill(HIST("h2DCADaugPt"), kinkCand.mothSign() * kinkCand.ptDaug(), kinkCand.dcaDaugPv()); + rSigmaMinus.fill(HIST("h2DCADaugPt"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.dcaDaugPv()); if (std::abs(mcTrackPiDau.pdgCode()) == 211) { rSigmaMinus.fill(HIST("h2NSigmaTOFPiPt"), kinkCand.mothSign() * kinkCand.ptMoth(), dauTrack.tofNSigmaPi()); @@ -283,7 +293,7 @@ struct sigmaminustask { PROCESS_SWITCH(sigmaminustask, processMC, "MC processing", false); - void processFindable(TracksFull const& tracks, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const&, CollisionsFullMC const&) + void processFindable(TracksFull const& tracks, aod::McTrackLabels const& trackLabelsMC, aod::KinkCands const& kinkCands, aod::McParticles const&, CollisionsFullMC const&) { for (const auto& motherTrack : tracks) { // Check if mother is Sigma in MC @@ -367,15 +377,23 @@ struct sigmaminustask { } else { continue; } - // 6 - geometric cuts: phi - if (std::abs(motherTrack.phi() - daughterTrack.phi()) * radToDeg < 100.0) { + // 6 - geometric cuts: phi difference + if (std::abs(motherTrack.phi() - daughterTrack.phi()) * radToDeg < 50.0) { filterIndex += 1; rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex); } else { continue; } - - // 7 - collision selection + // 7 - geometric cuts: z difference + o2::track::TrackParCov trackParCovMoth = getTrackParCov(motherTrack); + o2::track::TrackParCov trackParCovDaug = getTrackParCov(daughterTrack); + if (std::abs(trackParCovMoth.getZ() - trackParCovDaug.getZ()) < 20.0) { + filterIndex += 1; + rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex); + } else { + continue; + } + // 8 - collision selection auto collision = motherTrack.template collision_as(); if (!(std::abs(collision.posZ()) > cutzvertex || !collision.sel8())) { filterIndex += 1; @@ -394,3 +412,12 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) return WorkflowSpec{ adaptAnalysisTask(cfgc)}; } + +// Next steps: +// 0. Resolution histograms should have relative values, not absolute OK +// 1. New h2 with genRadius (recRadius) vs FilterIndex +// 2. Get recRadius through a map on kinkCands, put a negative value if the candidate is not reconstructed +// 2.1 Consider adding step in filters with the cuts on radius +// 2.2 Add h2 of radius resolution vs pt +// 3. Rewrite the findable method using maps to avoid the nested loop +// 4. For generated h2, avoid mass axis, use a bool axis to easily distinguish sigma minus and sigma plus \ No newline at end of file From 76134d5974bdeec6e7c1b65f81a15378453fdbca Mon Sep 17 00:00:00 2001 From: morgmatt Date: Tue, 5 Aug 2025 11:34:36 +0200 Subject: [PATCH 07/18] Removed nested loop from processFindable method, substituted with index map structure. Added decayradius detailed analysis of the filters --- .../Strangeness/sigmaminustask.cxx | 296 +++++++++++------- 1 file changed, 184 insertions(+), 112 deletions(-) diff --git a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx index 660b2ffc29f..ac27e665620 100644 --- a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx @@ -18,11 +18,9 @@ #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/PIDResponse.h" -#include "Common/Core/trackUtilities.h" #include "Framework/AnalysisTask.h" #include "Framework/runDataProcessing.h" #include "ReconstructionDataFormats/PID.h" -#include "ReconstructionDataFormats/Track.h" using namespace o2; using namespace o2::framework; @@ -66,13 +64,13 @@ struct sigmaminustask { const AxisSpec xiMassAxis{100, 1.2, 1.6, "m_{#Xi} (GeV/#it{c}^{2})"}; const AxisSpec pdgAxis{10001, -5000, 5000, "PDG code"}; const AxisSpec vertexZAxis{100, -15., 15., "vrtx_{Z} [cm]"}; - const AxisSpec dcaMothAxis{100, 0, 1, "DCA [cm]"}; + const AxisSpec dcaMothAxis{100, 0, 0.2, "DCA [cm]"}; const AxisSpec dcaDaugAxis{200, 0, 20, "DCA [cm]"}; const AxisSpec radiusAxis{100, -1, 40, "Decay radius [cm]"}; - const AxisSpec ptResolutionAxis{100, -2, 2, "(#it{p}_{T}^{rec} - #it{p}_{T}^{gen}) / #it{p}_{T}^{gen}"}; - const AxisSpec massResolutionAxis{100, -2, 2, "(m_{rec} - m_{gen}) / m_{gen}"}; - const AxisSpec radiusResolutionAxis{100, -2, 2, "(r_{rec} - r_{gen}) / r_{gen}"}; + const AxisSpec ptResolutionAxis{100, -1.0, 1.0, "(#it{p}_{T}^{rec} - #it{p}_{T}^{gen}) / #it{p}_{T}^{gen}"}; + const AxisSpec massResolutionAxis{100, -0.5, 0.5, "(m_{rec} - m_{gen}) / m_{gen}"}; + const AxisSpec radiusResolutionAxis{100, -0.5, 0.5, "(r_{rec} - r_{gen}) / r_{gen}"}; const AxisSpec boolAxis{2, -0.5, 1.5, "Boolean value (0=false, 1=true)"}; const AxisSpec filtersAxis{10, -0.5, 9.5, "Filter index"}; @@ -93,7 +91,8 @@ struct sigmaminustask { rSigmaMinus.add("h2MassResolution", "h2MassResolution", {HistType::kTH2F, {ptAxis, massResolutionAxis}}); rSigmaMinus.add("h2PtResolution", "h2PtResolution", {HistType::kTH2F, {ptAxis, ptResolutionAxis}}); - + rSigmaMinus.add("h2RadiusResolution", "h2RadiusResolution", {HistType::kTH2F, {ptAxis, radiusResolutionAxis}}); + rSigmaMinus.add("h2NSigmaTOFPiPt", "h2NSigmaTOFPiPt", {HistType::kTH2F, {ptAxis, nSigmaPiAxis}}); rSigmaMinus.add("h2NSigmaTOFPrPt", "h2NSigmaTOFPrPt", {HistType::kTH2F, {ptAxis, nSigmaPrAxis}}); @@ -108,9 +107,8 @@ struct sigmaminustask { // Add findable Sigma histograms rFindable.add("h2MassPtFindableAll", "h2MassPtFindableAll", {HistType::kTH2F, {ptAxis, sigmaMassAxis}}); rFindable.add("hFilterIndex", "hFilterIndex", {HistType::kTH1F, {filtersAxis}}); - rFindable.add("h2MCRadiusFilterIndex", "h2RadiusFilterIndex", {HistType::kTH2F, {filtersAxis, radiusAxis}}); + rFindable.add("h2MCRadiusFilterIndex", "h2MCRadiusFilterIndex", {HistType::kTH2F, {filtersAxis, radiusAxis}}); rFindable.add("h2RecRadiusFilterIndex", "h2RecRadiusFilterIndex", {HistType::kTH2F, {filtersAxis, radiusAxis}}); - } } @@ -293,114 +291,188 @@ struct sigmaminustask { PROCESS_SWITCH(sigmaminustask, processMC, "MC processing", false); - void processFindable(TracksFull const& tracks, aod::McTrackLabels const& trackLabelsMC, aod::KinkCands const& kinkCands, aod::McParticles const&, CollisionsFullMC const&) - { - for (const auto& motherTrack : tracks) { - // Check if mother is Sigma in MC - auto mcLabelMother = trackLabelsMC.rawIteratorAt(motherTrack.globalIndex()); - if (!mcLabelMother.has_mcParticle()) { + void processFindable(aod::KinkCands const& kinkCands, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const& particlesMC, + TracksFull const& tracks, CollisionsFullMC const& mcCollisions){ + // A - generated findable track pairs map: mcMother.globalIndex() -> (motherTrack.globalIndex(), daughterTrack.globalIndex()) + std::unordered_map> allCandsIndices; + + for (const auto& track : tracks) { + auto mcLabel = trackLabelsMC.rawIteratorAt(track.globalIndex()); + if (!mcLabel.has_mcParticle()) { continue; } - auto mcMother = mcLabelMother.mcParticle_as(); - if (std::abs(mcMother.pdgCode()) != 3112 && std::abs(mcMother.pdgCode()) != 3222) { + auto mcParticle = mcLabel.mcParticle_as(); + + if (mcParticle.has_daughters() && (std::abs(mcParticle.pdgCode()) == 3112 || std::abs(mcParticle.pdgCode()) == 3222)) { + allCandsIndices[mcParticle.globalIndex()] = {track.globalIndex(), -1}; + } + } + + for (const auto& track : tracks) { + auto mcLabel = trackLabelsMC.rawIteratorAt(track.globalIndex()); + if (!mcLabel.has_mcParticle()) { continue; } - - for (const auto& daughterTrack : tracks) { - // Check if daughter is pi/proton - auto mcLabelDaughter = trackLabelsMC.rawIteratorAt(daughterTrack.globalIndex()); - if (!mcLabelDaughter.has_mcParticle()) { - continue; - } - auto mcDaughter = mcLabelDaughter.mcParticle_as(); - if (std::abs(mcDaughter.pdgCode()) != 211 && std::abs(mcDaughter.pdgCode()) != 2212) { - continue; - } - - // Verify the MC mother-daughter relationship - bool isValidPair = false; - if (mcDaughter.has_mothers()) { - for (const auto& mother : mcDaughter.mothers_as()) { - if (mother.globalIndex() == mcMother.globalIndex()) { - isValidPair = true; - break; - } + auto mcParticle = mcLabel.mcParticle_as(); + + if (mcParticle.has_mothers() && (std::abs(mcParticle.pdgCode()) == 211 || std::abs(mcParticle.pdgCode()) == 2212)) { + for (const auto& mother : mcParticle.mothers_as()) { + auto it = allCandsIndices.find(mother.globalIndex()); + if (it != allCandsIndices.end()) { + it->second.second = track.globalIndex(); + break; } } - if (!isValidPair) { - continue; - } + } + } - float mcMass = std::sqrt(mcMother.e() * mcMother.e() - mcMother.p() * mcMother.p()); - int sigmaSign = mcMother.pdgCode() > 0 ? 1 : -1; - rSigmaMinus.fill(HIST("h2MassPtFindable"), sigmaSign * motherTrack.pt(), mcMass); - - // Define filter index and progressively apply kinkbuilder cuts to track pairs - int filterIndex = 0; - rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex); - - // 1 - tracks with right ITS, TPC, TOF signals - if (motherTrack.has_collision() && motherTrack.hasITS() && !motherTrack.hasTPC() && !motherTrack.hasTOF() && - daughterTrack.hasITS() && daughterTrack.hasTPC()) { - filterIndex += 1; - rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex); - } else { - continue; - } + // B - reconstructed kinkcands map: mcMother.globalIndex() -> kinkCand.globalIndex() + std::unordered_map findableToKinkCand; + for (const auto& kinkCand : kinkCands) { + auto motherTrack = kinkCand.trackMoth_as(); + auto daughterTrack = kinkCand.trackDaug_as(); + + auto mcLabMoth = trackLabelsMC.rawIteratorAt(motherTrack.globalIndex()); + auto mcLabDaug = trackLabelsMC.rawIteratorAt(daughterTrack.globalIndex()); + if (!mcLabMoth.has_mcParticle() || !mcLabDaug.has_mcParticle()) { + continue; + } + auto mcMother = mcLabMoth.mcParticle_as(); + auto mcDaughter = mcLabDaug.mcParticle_as(); + + if (std::abs(mcMother.pdgCode()) != 3112 && std::abs(mcMother.pdgCode()) != 3222) { + continue; + } + if (std::abs(mcDaughter.pdgCode()) != 211 && std::abs(mcDaughter.pdgCode()) != 2212) { + continue; + } + + auto findableIt = allCandsIndices.find(mcMother.globalIndex()); + if (findableIt != allCandsIndices.end() && + findableIt->second.first == motherTrack.globalIndex() && + findableIt->second.second == daughterTrack.globalIndex()) { + + findableToKinkCand[mcMother.globalIndex()] = kinkCand.globalIndex(); + } + } + + + // C - loop on valid pairs for findable analysis + for (const auto& [mcMotherIndex, trackIndices] : allCandsIndices) { + if (trackIndices.second == -1 || trackIndices.first == -1) { + continue; + } + // Retrieve mother and daughter tracks and mcParticles + auto motherTrack = tracks.rawIteratorAt(trackIndices.first); + auto daughterTrack = tracks.rawIteratorAt(trackIndices.second); + auto mcLabMoth = trackLabelsMC.rawIteratorAt(motherTrack.globalIndex()); + auto mcLabDaug = trackLabelsMC.rawIteratorAt(daughterTrack.globalIndex()); + if (!mcLabMoth.has_mcParticle() || !mcLabDaug.has_mcParticle()) { + continue; + } + auto mcMother = mcLabMoth.mcParticle_as(); + auto mcDaughter = mcLabDaug.mcParticle_as(); + + // Compute mass and radii + float mcMass = std::sqrt(mcMother.e() * mcMother.e() - mcMother.p() * mcMother.p()); + int sigmaSign = mcMother.pdgCode() > 0 ? 1 : -1; + float mcRadius = std::sqrt((mcMother.vx() - mcDaughter.vx()) * (mcMother.vx() - mcDaughter.vx()) + + (mcMother.vy() - mcDaughter.vy()) * (mcMother.vy() - mcDaughter.vy())); + float recRadius = -1.0; + if (findableToKinkCand.find(mcMother.globalIndex()) != findableToKinkCand.end()) { + auto kinkCand = kinkCands.rawIteratorAt(findableToKinkCand[mcMother.globalIndex()]); + recRadius = std::sqrt(kinkCand.xDecVtx() * kinkCand.xDecVtx() + kinkCand.yDecVtx() * kinkCand.yDecVtx()); + } + + rFindable.fill(HIST("h2MassPtFindableAll"), sigmaSign * mcMother.pt(), mcMass); + + // Define filter index and progressively apply kinkbuilder cuts to track pairs + int filterIndex = 0; + rFindable.fill(HIST("hFilterIndex"), filterIndex); + rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); + rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius); + + // 1 - tracks with right ITS, TPC, TOF signals + if (motherTrack.has_collision() && motherTrack.hasITS() && !motherTrack.hasTPC() && !motherTrack.hasTOF() && + daughterTrack.hasITS() && daughterTrack.hasTPC()) { + filterIndex += 1; + rFindable.fill(HIST("hFilterIndex"), filterIndex); + rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); + rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius); + } else { + continue; + } - // 2, 3 - mother track ITS properties - if (motherTrack.itsNCls() < 6 && + // 2, 3 - mother track ITS properties + if (motherTrack.itsNCls() < 6 && motherTrack.itsNClsInnerBarrel() == 3 && motherTrack.itsChi2NCl() < 36) { - filterIndex += 1; - rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex); - } - if (filterIndex > 1 && motherTrack.pt() > 0.5) { - filterIndex += 1; - rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex); - } else { - continue; - } + filterIndex += 1; + rFindable.fill(HIST("hFilterIndex"), filterIndex); + rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); + rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius); + } else { + continue; + } - // 4 - daughter track ITS+TPC properties - if (daughterTrack.itsNClsInnerBarrel() == 0 && daughterTrack.itsNCls() < 4 && - daughterTrack.tpcNClsCrossedRows() > 0.8 * daughterTrack.tpcNClsFindable() && daughterTrack.tpcNClsFound() > 80) { - filterIndex += 1; - rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex); - } else { - continue; - } - - // 5 - geometric cuts: eta - if (std::abs(motherTrack.eta()) < 1.0 && std::abs(daughterTrack.eta()) < 1.0) { - filterIndex += 1; - rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex); - } else { - continue; - } - // 6 - geometric cuts: phi difference - if (std::abs(motherTrack.phi() - daughterTrack.phi()) * radToDeg < 50.0) { - filterIndex += 1; - rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex); - } else { - continue; - } - // 7 - geometric cuts: z difference - o2::track::TrackParCov trackParCovMoth = getTrackParCov(motherTrack); - o2::track::TrackParCov trackParCovDaug = getTrackParCov(daughterTrack); - if (std::abs(trackParCovMoth.getZ() - trackParCovDaug.getZ()) < 20.0) { - filterIndex += 1; - rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex); - } else { - continue; - } - // 8 - collision selection - auto collision = motherTrack.template collision_as(); - if (!(std::abs(collision.posZ()) > cutzvertex || !collision.sel8())) { - filterIndex += 1; - rSigmaMinus.fill(HIST("hFilterIndex"), filterIndex); - } - + if (motherTrack.pt() > 0.5) { + filterIndex += 1; + rFindable.fill(HIST("hFilterIndex"), filterIndex); + rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); + rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius); + } else { + continue; + } + + // 4 - daughter track ITS+TPC properties + if (daughterTrack.itsNClsInnerBarrel() == 0 && daughterTrack.itsNCls() < 4 && + daughterTrack.tpcNClsCrossedRows() > 0.8 * daughterTrack.tpcNClsFindable() && daughterTrack.tpcNClsFound() > 80) { + filterIndex += 1; + rFindable.fill(HIST("hFilterIndex"), filterIndex); + rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); + rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius); + } else { + continue; + } + + // 5 - geometric cuts: eta + if (std::abs(motherTrack.eta()) < 1.0 && std::abs(daughterTrack.eta()) < 1.0) { + filterIndex += 1; + rFindable.fill(HIST("hFilterIndex"), filterIndex); + rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); + rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius); + } else { + continue; } + + // 6 - geometric cuts: phi difference + if (std::abs(motherTrack.phi() - daughterTrack.phi()) * radToDeg < 50.0) { + filterIndex += 1; + rFindable.fill(HIST("hFilterIndex"), filterIndex); + rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); + rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius); + } else { + continue; + } + + // 7 - radius cut + if (recRadius > 19.6213f) { + filterIndex += 1; + rFindable.fill(HIST("hFilterIndex"), filterIndex); + rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); + rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius); + } else { + continue; + } + + // 8 - collision selection + auto collision = motherTrack.template collision_as(); + if (!(std::abs(collision.posZ()) > cutzvertex || !collision.sel8())) { + filterIndex += 1; + rFindable.fill(HIST("hFilterIndex"), filterIndex); + rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); + rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius); + } + } } @@ -414,10 +486,10 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) } // Next steps: -// 0. Resolution histograms should have relative values, not absolute OK -// 1. New h2 with genRadius (recRadius) vs FilterIndex -// 2. Get recRadius through a map on kinkCands, put a negative value if the candidate is not reconstructed -// 2.1 Consider adding step in filters with the cuts on radius -// 2.2 Add h2 of radius resolution vs pt -// 3. Rewrite the findable method using maps to avoid the nested loop +// 0. Resolution histograms should have relative values, not absolute: OK +// 1. New h2 with genRadius (recRadius) vs FilterIndex: OK +// 2. Get recRadius through a map on kinkCands, put a negative value if the candidate is not reconstructed: OK +// 2.1 Consider adding step in filters with the cuts on radius: OK +// 2.2 Add h2 of radius resolution vs pt: OK +// 3. Rewrite the findable method using maps to avoid the nested loop: OK // 4. For generated h2, avoid mass axis, use a bool axis to easily distinguish sigma minus and sigma plus \ No newline at end of file From 7b7d06e4e61d0449eebd0777ec6ddfacf6752fc0 Mon Sep 17 00:00:00 2001 From: morgmatt Date: Tue, 5 Aug 2025 14:34:28 +0200 Subject: [PATCH 08/18] Added configurables for parameters that control the kinkBuilder cuts --- .../TableProducer/Strangeness/sigmaminustask.cxx | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx index ac27e665620..abbdc34d971 100644 --- a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx @@ -49,6 +49,12 @@ struct sigmaminustask { Configurable cutRapMotherMC{"cutRapMotherMC", 1.0f, "Rapidity cut for mother Sigma in MC"}; Configurable fillOutputTree{"fillOutputTree", true, "If true, fill the output tree with Kink candidates"}; + // Configurables for findable tracks (see kinkBuilder.cxx) + Configurable KBminPtMoth{"KBminPtMoth", 0.5, "Minimum pT of the mother"}; + Configurable KBmaxPhiDiff{"KBmaxPhiDiff", 100, "Max phi difference between the kink daughter and the mother"}; + Configurable KBetaMax{"KBetaMax", 1., "Max eta for both mother and daughter"}; + Configurable KBnTPCClusMinDaug{"KBnTPCClusMinDaug", 80, "Min daug NTPC clusters"}; + Configurable KBradiusCut{"KBradiusCut", 19.6213f, "Min reconstructed decay radius of the mother"}; Preslice mPerCol = aod::track::collisionId; @@ -414,7 +420,7 @@ struct sigmaminustask { continue; } - if (motherTrack.pt() > 0.5) { + if (motherTrack.pt() > KBminPtMoth) { filterIndex += 1; rFindable.fill(HIST("hFilterIndex"), filterIndex); rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); @@ -425,7 +431,7 @@ struct sigmaminustask { // 4 - daughter track ITS+TPC properties if (daughterTrack.itsNClsInnerBarrel() == 0 && daughterTrack.itsNCls() < 4 && - daughterTrack.tpcNClsCrossedRows() > 0.8 * daughterTrack.tpcNClsFindable() && daughterTrack.tpcNClsFound() > 80) { + daughterTrack.tpcNClsCrossedRows() > 0.8 * daughterTrack.tpcNClsFindable() && daughterTrack.tpcNClsFound() > KBnTPCClusMinDaug) { filterIndex += 1; rFindable.fill(HIST("hFilterIndex"), filterIndex); rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); @@ -435,7 +441,7 @@ struct sigmaminustask { } // 5 - geometric cuts: eta - if (std::abs(motherTrack.eta()) < 1.0 && std::abs(daughterTrack.eta()) < 1.0) { + if (std::abs(motherTrack.eta()) < KBetaMax && std::abs(daughterTrack.eta()) < KBetaMax) { filterIndex += 1; rFindable.fill(HIST("hFilterIndex"), filterIndex); rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); @@ -445,7 +451,7 @@ struct sigmaminustask { } // 6 - geometric cuts: phi difference - if (std::abs(motherTrack.phi() - daughterTrack.phi()) * radToDeg < 50.0) { + if (std::abs(motherTrack.phi() - daughterTrack.phi()) * radToDeg < KBmaxPhiDiff) { filterIndex += 1; rFindable.fill(HIST("hFilterIndex"), filterIndex); rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); @@ -455,7 +461,7 @@ struct sigmaminustask { } // 7 - radius cut - if (recRadius > 19.6213f) { + if (recRadius > KBradiusCut) { filterIndex += 1; rFindable.fill(HIST("hFilterIndex"), filterIndex); rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); From 33fbac89a46f3de9c03392b3389735da4a5c1031 Mon Sep 17 00:00:00 2001 From: morgmatt Date: Tue, 5 Aug 2025 15:21:38 +0200 Subject: [PATCH 09/18] Added helper function for filling findable histograms --- .../Strangeness/sigmaminustask.cxx | 71 ++++++++++--------- 1 file changed, 38 insertions(+), 33 deletions(-) diff --git a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx index abbdc34d971..5f5a438bb5d 100644 --- a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx @@ -115,6 +115,8 @@ struct sigmaminustask { rFindable.add("hFilterIndex", "hFilterIndex", {HistType::kTH1F, {filtersAxis}}); rFindable.add("h2MCRadiusFilterIndex", "h2MCRadiusFilterIndex", {HistType::kTH2F, {filtersAxis, radiusAxis}}); rFindable.add("h2RecRadiusFilterIndex", "h2RecRadiusFilterIndex", {HistType::kTH2F, {filtersAxis, radiusAxis}}); + rFindable.add("h2MCRadiusFilter_protonkink", "h2MCRadiusFilter_protonkink", {HistType::kTH2F, {filtersAxis, radiusAxis}}); + rFindable.add("h2MCRadiusFilter_pikink", "h2MCRadiusFilter_pikink", {HistType::kTH2F, {filtersAxis, radiusAxis}}); } } @@ -297,6 +299,18 @@ struct sigmaminustask { PROCESS_SWITCH(sigmaminustask, processMC, "MC processing", false); + void fillFindableHistograms(int filterIndex, float mcRadius, float recRadius, bool isPiDaughter) { + rFindable.fill(HIST("hFilterIndex"), filterIndex); + rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); + rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius); + + if (isPiDaughter) { + rFindable.fill(HIST("h2MCRadiusFilter_pikink"), filterIndex, mcRadius); + } else { + rFindable.fill(HIST("h2MCRadiusFilter_protonkink"), filterIndex, mcRadius); + } + } + void processFindable(aod::KinkCands const& kinkCands, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const& particlesMC, TracksFull const& tracks, CollisionsFullMC const& mcCollisions){ // A - generated findable track pairs map: mcMother.globalIndex() -> (motherTrack.globalIndex(), daughterTrack.globalIndex()) @@ -382,6 +396,7 @@ struct sigmaminustask { // Compute mass and radii float mcMass = std::sqrt(mcMother.e() * mcMother.e() - mcMother.p() * mcMother.p()); int sigmaSign = mcMother.pdgCode() > 0 ? 1 : -1; + bool isPiDaughter = std::abs(mcDaughter.pdgCode()) == 211; float mcRadius = std::sqrt((mcMother.vx() - mcDaughter.vx()) * (mcMother.vx() - mcDaughter.vx()) + (mcMother.vy() - mcDaughter.vy()) * (mcMother.vy() - mcDaughter.vy())); float recRadius = -1.0; @@ -394,17 +409,14 @@ struct sigmaminustask { // Define filter index and progressively apply kinkbuilder cuts to track pairs int filterIndex = 0; - rFindable.fill(HIST("hFilterIndex"), filterIndex); - rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); - rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius); + fillFindableHistograms(filterIndex, mcRadius, recRadius, isPiDaughter); // 1 - tracks with right ITS, TPC, TOF signals if (motherTrack.has_collision() && motherTrack.hasITS() && !motherTrack.hasTPC() && !motherTrack.hasTOF() && daughterTrack.hasITS() && daughterTrack.hasTPC()) { filterIndex += 1; - rFindable.fill(HIST("hFilterIndex"), filterIndex); - rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); - rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius); + fillFindableHistograms(filterIndex, mcRadius, recRadius, isPiDaughter); + } else { continue; } @@ -413,18 +425,16 @@ struct sigmaminustask { if (motherTrack.itsNCls() < 6 && motherTrack.itsNClsInnerBarrel() == 3 && motherTrack.itsChi2NCl() < 36) { filterIndex += 1; - rFindable.fill(HIST("hFilterIndex"), filterIndex); - rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); - rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius); + fillFindableHistograms(filterIndex, mcRadius, recRadius, isPiDaughter); + } else { continue; } if (motherTrack.pt() > KBminPtMoth) { filterIndex += 1; - rFindable.fill(HIST("hFilterIndex"), filterIndex); - rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); - rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius); + fillFindableHistograms(filterIndex, mcRadius, recRadius, isPiDaughter); + } else { continue; } @@ -433,9 +443,8 @@ struct sigmaminustask { if (daughterTrack.itsNClsInnerBarrel() == 0 && daughterTrack.itsNCls() < 4 && daughterTrack.tpcNClsCrossedRows() > 0.8 * daughterTrack.tpcNClsFindable() && daughterTrack.tpcNClsFound() > KBnTPCClusMinDaug) { filterIndex += 1; - rFindable.fill(HIST("hFilterIndex"), filterIndex); - rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); - rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius); + fillFindableHistograms(filterIndex, mcRadius, recRadius, isPiDaughter); + } else { continue; } @@ -443,9 +452,8 @@ struct sigmaminustask { // 5 - geometric cuts: eta if (std::abs(motherTrack.eta()) < KBetaMax && std::abs(daughterTrack.eta()) < KBetaMax) { filterIndex += 1; - rFindable.fill(HIST("hFilterIndex"), filterIndex); - rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); - rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius); + fillFindableHistograms(filterIndex, mcRadius, recRadius, isPiDaughter); + } else { continue; } @@ -453,9 +461,8 @@ struct sigmaminustask { // 6 - geometric cuts: phi difference if (std::abs(motherTrack.phi() - daughterTrack.phi()) * radToDeg < KBmaxPhiDiff) { filterIndex += 1; - rFindable.fill(HIST("hFilterIndex"), filterIndex); - rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); - rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius); + fillFindableHistograms(filterIndex, mcRadius, recRadius, isPiDaughter); + } else { continue; } @@ -463,9 +470,8 @@ struct sigmaminustask { // 7 - radius cut if (recRadius > KBradiusCut) { filterIndex += 1; - rFindable.fill(HIST("hFilterIndex"), filterIndex); - rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); - rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius); + fillFindableHistograms(filterIndex, mcRadius, recRadius, isPiDaughter); + } else { continue; } @@ -474,11 +480,16 @@ struct sigmaminustask { auto collision = motherTrack.template collision_as(); if (!(std::abs(collision.posZ()) > cutzvertex || !collision.sel8())) { filterIndex += 1; - rFindable.fill(HIST("hFilterIndex"), filterIndex); - rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); - rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius); + fillFindableHistograms(filterIndex, mcRadius, recRadius, isPiDaughter); + } + // 9 - TOF daughter presence + if (daughterTrack.hasTOF()) { + filterIndex += 1; + fillFindableHistograms(filterIndex, mcRadius, recRadius, isPiDaughter); + + } } } @@ -492,10 +503,4 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) } // Next steps: -// 0. Resolution histograms should have relative values, not absolute: OK -// 1. New h2 with genRadius (recRadius) vs FilterIndex: OK -// 2. Get recRadius through a map on kinkCands, put a negative value if the candidate is not reconstructed: OK -// 2.1 Consider adding step in filters with the cuts on radius: OK -// 2.2 Add h2 of radius resolution vs pt: OK -// 3. Rewrite the findable method using maps to avoid the nested loop: OK // 4. For generated h2, avoid mass axis, use a bool axis to easily distinguish sigma minus and sigma plus \ No newline at end of file From 312d4f3d0bd1c843b7dc3fa6d673ac93bf474647 Mon Sep 17 00:00:00 2001 From: morgmatt Date: Wed, 6 Aug 2025 12:28:10 +0200 Subject: [PATCH 10/18] Added findable histograms with daughter pt analysis --- .../Strangeness/sigmaminustask.cxx | 83 ++++++++++++------- 1 file changed, 54 insertions(+), 29 deletions(-) diff --git a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx index 5f5a438bb5d..2e8cb3e6895 100644 --- a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx @@ -64,6 +64,7 @@ struct sigmaminustask { { // Axes const AxisSpec ptAxis{100, -10, 10, "#it{p}_{T} (GeV/#it{c})"}; + const AxisSpec ptDaughterAxis{200, 0, 10, "Daughter #it{p}_{T} (GeV/#it{c})"}; const AxisSpec nSigmaPiAxis{60, -30, 30, "n#sigma_{#pi}"}; const AxisSpec nSigmaPrAxis{60, -30, 30, "n#sigma_{p}"}; const AxisSpec sigmaMassAxis{100, 1.1, 1.4, "m (GeV/#it{c}^{2})"}; @@ -78,7 +79,7 @@ struct sigmaminustask { const AxisSpec massResolutionAxis{100, -0.5, 0.5, "(m_{rec} - m_{gen}) / m_{gen}"}; const AxisSpec radiusResolutionAxis{100, -0.5, 0.5, "(r_{rec} - r_{gen}) / r_{gen}"}; - const AxisSpec boolAxis{2, -0.5, 1.5, "Boolean value (0=false, 1=true)"}; + const AxisSpec boolAxis{2, -0.5, 1.5, "Boolean value"}; const AxisSpec filtersAxis{10, -0.5, 9.5, "Filter index"}; // Event selection @@ -111,12 +112,30 @@ struct sigmaminustask { if (doprocessFindable) { // Add findable Sigma histograms - rFindable.add("h2MassPtFindableAll", "h2MassPtFindableAll", {HistType::kTH2F, {ptAxis, sigmaMassAxis}}); rFindable.add("hFilterIndex", "hFilterIndex", {HistType::kTH1F, {filtersAxis}}); + + auto hFilterIndex = rFindable.get(HIST("hFilterIndex")); + hFilterIndex->GetXaxis()->SetBinLabel(1, "Initial"); // 0: all generated pairs with tracks + hFilterIndex->GetXaxis()->SetBinLabel(2, "Det. signals"); // 1: correct ITS, TPC signals present + hFilterIndex->GetXaxis()->SetBinLabel(3, "Moth ITS"); // 2: cuts on ITS mother signal + hFilterIndex->GetXaxis()->SetBinLabel(4, "Moth p_{T}"); // 3: mother pT > KBminPtMoth + hFilterIndex->GetXaxis()->SetBinLabel(5, "Daug TPC"); // 4: cuts on TPC daughter signal + hFilterIndex->GetXaxis()->SetBinLabel(6, "#eta cuts"); // 5: eta < KBetaMax + hFilterIndex->GetXaxis()->SetBinLabel(7, "#Delta#phi"); // 6: delta phi < KBmaxPhiDiff + hFilterIndex->GetXaxis()->SetBinLabel(8, "Radius"); // 7: decay radius > KBradiusCut + hFilterIndex->GetXaxis()->SetBinLabel(9, "Collision"); // 8: collision sel8 cut + hFilterIndex->GetXaxis()->SetBinLabel(10, "TOF daug"); // 9: daughter has TOF + rFindable.add("h2MCRadiusFilterIndex", "h2MCRadiusFilterIndex", {HistType::kTH2F, {filtersAxis, radiusAxis}}); rFindable.add("h2RecRadiusFilterIndex", "h2RecRadiusFilterIndex", {HistType::kTH2F, {filtersAxis, radiusAxis}}); - rFindable.add("h2MCRadiusFilter_protonkink", "h2MCRadiusFilter_protonkink", {HistType::kTH2F, {filtersAxis, radiusAxis}}); - rFindable.add("h2MCRadiusFilter_pikink", "h2MCRadiusFilter_pikink", {HistType::kTH2F, {filtersAxis, radiusAxis}}); + + rFindable.add("h2MCRadiusFilter_plus_protonkink", "h2MCRadiusFilter_plus_protonkink", {HistType::kTH2F, {filtersAxis, radiusAxis}}); + rFindable.add("h2MCRadiusFilter_plus_pikink", "h2MCRadiusFilter_plus_pikink", {HistType::kTH2F, {filtersAxis, radiusAxis}}); + rFindable.add("h2MCRadiusFilter_minus_pikink", "h2MCRadiusFilter_minus_pikink", {HistType::kTH2F, {filtersAxis, radiusAxis}}); + + rFindable.add("h2PtFilter_plus_protonkink", "h2PtFilter_plus_protonkink", {HistType::kTH2F, {filtersAxis, ptDaughterAxis}}); + rFindable.add("h2PtFilter_plus_pikink", "h2PtFilter_plus_pikink", {HistType::kTH2F, {filtersAxis, ptDaughterAxis}}); + rFindable.add("h2PtFilter_minus_pikink", "h2PtFilter_minus_pikink", {HistType::kTH2F, {filtersAxis, ptDaughterAxis}}); } } @@ -299,20 +318,29 @@ struct sigmaminustask { PROCESS_SWITCH(sigmaminustask, processMC, "MC processing", false); - void fillFindableHistograms(int filterIndex, float mcRadius, float recRadius, bool isPiDaughter) { + void fillFindableHistograms(int filterIndex, float mcRadius, float recRadius, float ptDaughter, bool isSigmaMinus, bool isPiDaughter) { rFindable.fill(HIST("hFilterIndex"), filterIndex); rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius); if (isPiDaughter) { - rFindable.fill(HIST("h2MCRadiusFilter_pikink"), filterIndex, mcRadius); + if (isSigmaMinus) { + rFindable.fill(HIST("h2MCRadiusFilter_minus_pikink"), filterIndex, mcRadius); + rFindable.fill(HIST("h2PtFilter_minus_pikink"), filterIndex, ptDaughter); + } else { + rFindable.fill(HIST("h2MCRadiusFilter_plus_pikink"), filterIndex, mcRadius); + rFindable.fill(HIST("h2PtFilter_plus_pikink"), filterIndex, ptDaughter); + } } else { - rFindable.fill(HIST("h2MCRadiusFilter_protonkink"), filterIndex, mcRadius); + if (!isSigmaMinus) { + rFindable.fill(HIST("h2MCRadiusFilter_plus_protonkink"), filterIndex, mcRadius); + rFindable.fill(HIST("h2PtFilter_plus_protonkink"), filterIndex, ptDaughter); + } } } - void processFindable(aod::KinkCands const& kinkCands, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const& particlesMC, - TracksFull const& tracks, CollisionsFullMC const& mcCollisions){ + void processFindable(aod::KinkCands const& kinkCands, aod::McTrackLabels const& trackLabelsMC, + TracksFull const& tracks, aod::McParticles const&, CollisionsFullMC const&){ // A - generated findable track pairs map: mcMother.globalIndex() -> (motherTrack.globalIndex(), daughterTrack.globalIndex()) std::unordered_map> allCandsIndices; @@ -394,28 +422,28 @@ struct sigmaminustask { auto mcDaughter = mcLabDaug.mcParticle_as(); // Compute mass and radii - float mcMass = std::sqrt(mcMother.e() * mcMother.e() - mcMother.p() * mcMother.p()); - int sigmaSign = mcMother.pdgCode() > 0 ? 1 : -1; - bool isPiDaughter = std::abs(mcDaughter.pdgCode()) == 211; - float mcRadius = std::sqrt((mcMother.vx() - mcDaughter.vx()) * (mcMother.vx() - mcDaughter.vx()) + - (mcMother.vy() - mcDaughter.vy()) * (mcMother.vy() - mcDaughter.vy())); + bool isSigmaMinus = (std::abs(mcMother.pdgCode()) == 3112); + bool isPiDaughter = (std::abs(mcDaughter.pdgCode()) == 211); + + //float mcMass = std::sqrt(mcMother.e() * mcMother.e() - mcMother.p() * mcMother.p()); + //int sigmaSign = mcMother.pdgCode() > 0 ? 1 : -1; + float mcRadius = std::sqrt((mcMother.vx() - mcDaughter.vx()) * (mcMother.vx() - mcDaughter.vx()) + (mcMother.vy() - mcDaughter.vy()) * (mcMother.vy() - mcDaughter.vy())); float recRadius = -1.0; if (findableToKinkCand.find(mcMother.globalIndex()) != findableToKinkCand.end()) { auto kinkCand = kinkCands.rawIteratorAt(findableToKinkCand[mcMother.globalIndex()]); recRadius = std::sqrt(kinkCand.xDecVtx() * kinkCand.xDecVtx() + kinkCand.yDecVtx() * kinkCand.yDecVtx()); } - - rFindable.fill(HIST("h2MassPtFindableAll"), sigmaSign * mcMother.pt(), mcMass); + float recPtDaughter = daughterTrack.pt(); // Define filter index and progressively apply kinkbuilder cuts to track pairs int filterIndex = 0; - fillFindableHistograms(filterIndex, mcRadius, recRadius, isPiDaughter); + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtDaughter, isSigmaMinus, isPiDaughter); // 1 - tracks with right ITS, TPC, TOF signals if (motherTrack.has_collision() && motherTrack.hasITS() && !motherTrack.hasTPC() && !motherTrack.hasTOF() && daughterTrack.hasITS() && daughterTrack.hasTPC()) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, isPiDaughter); + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtDaughter, isSigmaMinus, isPiDaughter); } else { continue; @@ -425,7 +453,7 @@ struct sigmaminustask { if (motherTrack.itsNCls() < 6 && motherTrack.itsNClsInnerBarrel() == 3 && motherTrack.itsChi2NCl() < 36) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, isPiDaughter); + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtDaughter, isSigmaMinus, isPiDaughter); } else { continue; @@ -433,7 +461,7 @@ struct sigmaminustask { if (motherTrack.pt() > KBminPtMoth) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, isPiDaughter); + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtDaughter, isSigmaMinus, isPiDaughter); } else { continue; @@ -443,7 +471,7 @@ struct sigmaminustask { if (daughterTrack.itsNClsInnerBarrel() == 0 && daughterTrack.itsNCls() < 4 && daughterTrack.tpcNClsCrossedRows() > 0.8 * daughterTrack.tpcNClsFindable() && daughterTrack.tpcNClsFound() > KBnTPCClusMinDaug) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, isPiDaughter); + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtDaughter, isSigmaMinus, isPiDaughter); } else { continue; @@ -452,7 +480,7 @@ struct sigmaminustask { // 5 - geometric cuts: eta if (std::abs(motherTrack.eta()) < KBetaMax && std::abs(daughterTrack.eta()) < KBetaMax) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, isPiDaughter); + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtDaughter, isSigmaMinus, isPiDaughter); } else { continue; @@ -461,7 +489,7 @@ struct sigmaminustask { // 6 - geometric cuts: phi difference if (std::abs(motherTrack.phi() - daughterTrack.phi()) * radToDeg < KBmaxPhiDiff) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, isPiDaughter); + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtDaughter, isSigmaMinus, isPiDaughter); } else { continue; @@ -470,7 +498,7 @@ struct sigmaminustask { // 7 - radius cut if (recRadius > KBradiusCut) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, isPiDaughter); + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtDaughter, isSigmaMinus, isPiDaughter); } else { continue; @@ -480,14 +508,14 @@ struct sigmaminustask { auto collision = motherTrack.template collision_as(); if (!(std::abs(collision.posZ()) > cutzvertex || !collision.sel8())) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, isPiDaughter); + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtDaughter, isSigmaMinus, isPiDaughter); } // 9 - TOF daughter presence if (daughterTrack.hasTOF()) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, isPiDaughter); + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtDaughter, isSigmaMinus, isPiDaughter); } } @@ -501,6 +529,3 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) return WorkflowSpec{ adaptAnalysisTask(cfgc)}; } - -// Next steps: -// 4. For generated h2, avoid mass axis, use a bool axis to easily distinguish sigma minus and sigma plus \ No newline at end of file From 9d233edca1855d3d3eeca5e6fec7167e6fa440c9 Mon Sep 17 00:00:00 2001 From: morgmatt Date: Wed, 6 Aug 2025 15:15:51 +0200 Subject: [PATCH 11/18] Added labels to findable histogram filter index axis --- .../Strangeness/sigmaminustask.cxx | 74 +++++++++---------- 1 file changed, 33 insertions(+), 41 deletions(-) diff --git a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx index 2e8cb3e6895..5ff82b381ea 100644 --- a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx @@ -64,7 +64,7 @@ struct sigmaminustask { { // Axes const AxisSpec ptAxis{100, -10, 10, "#it{p}_{T} (GeV/#it{c})"}; - const AxisSpec ptDaughterAxis{200, 0, 10, "Daughter #it{p}_{T} (GeV/#it{c})"}; + const AxisSpec ptUnsignedAxis{100, 0, 10, "#it{p}_{T} (GeV/#it{c})"}; const AxisSpec nSigmaPiAxis{60, -30, 30, "n#sigma_{#pi}"}; const AxisSpec nSigmaPrAxis{60, -30, 30, "n#sigma_{p}"}; const AxisSpec sigmaMassAxis{100, 1.1, 1.4, "m (GeV/#it{c}^{2})"}; @@ -111,31 +111,29 @@ struct sigmaminustask { } if (doprocessFindable) { + std::vector filterLabels = {"Initial", "ITS/TPC present", "Moth ITS", "Moth p_{T}", "Daug ITS+TPC", "#eta", "#Delta#phi", "Radius", "Collision", "Daug TOF"}; + // Add findable Sigma histograms rFindable.add("hFilterIndex", "hFilterIndex", {HistType::kTH1F, {filtersAxis}}); - - auto hFilterIndex = rFindable.get(HIST("hFilterIndex")); - hFilterIndex->GetXaxis()->SetBinLabel(1, "Initial"); // 0: all generated pairs with tracks - hFilterIndex->GetXaxis()->SetBinLabel(2, "Det. signals"); // 1: correct ITS, TPC signals present - hFilterIndex->GetXaxis()->SetBinLabel(3, "Moth ITS"); // 2: cuts on ITS mother signal - hFilterIndex->GetXaxis()->SetBinLabel(4, "Moth p_{T}"); // 3: mother pT > KBminPtMoth - hFilterIndex->GetXaxis()->SetBinLabel(5, "Daug TPC"); // 4: cuts on TPC daughter signal - hFilterIndex->GetXaxis()->SetBinLabel(6, "#eta cuts"); // 5: eta < KBetaMax - hFilterIndex->GetXaxis()->SetBinLabel(7, "#Delta#phi"); // 6: delta phi < KBmaxPhiDiff - hFilterIndex->GetXaxis()->SetBinLabel(8, "Radius"); // 7: decay radius > KBradiusCut - hFilterIndex->GetXaxis()->SetBinLabel(9, "Collision"); // 8: collision sel8 cut - hFilterIndex->GetXaxis()->SetBinLabel(10, "TOF daug"); // 9: daughter has TOF - rFindable.add("h2MCRadiusFilterIndex", "h2MCRadiusFilterIndex", {HistType::kTH2F, {filtersAxis, radiusAxis}}); rFindable.add("h2RecRadiusFilterIndex", "h2RecRadiusFilterIndex", {HistType::kTH2F, {filtersAxis, radiusAxis}}); + auto hFilterIndex = rFindable.get(HIST("hFilterIndex")); + auto h2MCRadiusFilterIndex = rFindable.get(HIST("h2MCRadiusFilterIndex")); + auto h2RecRadiusFilterIndex = rFindable.get(HIST("h2RecRadiusFilterIndex")); + for (size_t i = 0; i < filterLabels.size(); ++i) { + hFilterIndex->GetXaxis()->SetBinLabel(i + 1, filterLabels[i].c_str()); + h2MCRadiusFilterIndex->GetXaxis()->SetBinLabel(i + 1, filterLabels[i].c_str()); + h2RecRadiusFilterIndex->GetXaxis()->SetBinLabel(i + 1, filterLabels[i].c_str()); + } + rFindable.add("h2MCRadiusFilter_plus_protonkink", "h2MCRadiusFilter_plus_protonkink", {HistType::kTH2F, {filtersAxis, radiusAxis}}); rFindable.add("h2MCRadiusFilter_plus_pikink", "h2MCRadiusFilter_plus_pikink", {HistType::kTH2F, {filtersAxis, radiusAxis}}); rFindable.add("h2MCRadiusFilter_minus_pikink", "h2MCRadiusFilter_minus_pikink", {HistType::kTH2F, {filtersAxis, radiusAxis}}); - rFindable.add("h2PtFilter_plus_protonkink", "h2PtFilter_plus_protonkink", {HistType::kTH2F, {filtersAxis, ptDaughterAxis}}); - rFindable.add("h2PtFilter_plus_pikink", "h2PtFilter_plus_pikink", {HistType::kTH2F, {filtersAxis, ptDaughterAxis}}); - rFindable.add("h2PtFilter_minus_pikink", "h2PtFilter_minus_pikink", {HistType::kTH2F, {filtersAxis, ptDaughterAxis}}); + rFindable.add("h2PtFilter_plus_protonkink", "h2PtFilter_plus_protonkink", {HistType::kTH2F, {filtersAxis, ptUnsignedAxis}}); + rFindable.add("h2PtFilter_plus_pikink", "h2PtFilter_plus_pikink", {HistType::kTH2F, {filtersAxis, ptUnsignedAxis}}); + rFindable.add("h2PtFilter_minus_pikink", "h2PtFilter_minus_pikink", {HistType::kTH2F, {filtersAxis, ptUnsignedAxis}}); } } @@ -318,7 +316,7 @@ struct sigmaminustask { PROCESS_SWITCH(sigmaminustask, processMC, "MC processing", false); - void fillFindableHistograms(int filterIndex, float mcRadius, float recRadius, float ptDaughter, bool isSigmaMinus, bool isPiDaughter) { + void fillFindableHistograms(int filterIndex, float mcRadius, float recRadius, float pt, bool isSigmaMinus, bool isPiDaughter) { rFindable.fill(HIST("hFilterIndex"), filterIndex); rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius); @@ -326,15 +324,15 @@ struct sigmaminustask { if (isPiDaughter) { if (isSigmaMinus) { rFindable.fill(HIST("h2MCRadiusFilter_minus_pikink"), filterIndex, mcRadius); - rFindable.fill(HIST("h2PtFilter_minus_pikink"), filterIndex, ptDaughter); + rFindable.fill(HIST("h2PtFilter_minus_pikink"), filterIndex, pt); } else { rFindable.fill(HIST("h2MCRadiusFilter_plus_pikink"), filterIndex, mcRadius); - rFindable.fill(HIST("h2PtFilter_plus_pikink"), filterIndex, ptDaughter); - } + rFindable.fill(HIST("h2PtFilter_plus_pikink"), filterIndex, pt); + } } else { if (!isSigmaMinus) { rFindable.fill(HIST("h2MCRadiusFilter_plus_protonkink"), filterIndex, mcRadius); - rFindable.fill(HIST("h2PtFilter_plus_protonkink"), filterIndex, ptDaughter); + rFindable.fill(HIST("h2PtFilter_plus_protonkink"), filterIndex, pt); } } } @@ -434,17 +432,17 @@ struct sigmaminustask { recRadius = std::sqrt(kinkCand.xDecVtx() * kinkCand.xDecVtx() + kinkCand.yDecVtx() * kinkCand.yDecVtx()); } float recPtDaughter = daughterTrack.pt(); + float recPtMother = motherTrack.pt(); // Define filter index and progressively apply kinkbuilder cuts to track pairs int filterIndex = 0; - fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtDaughter, isSigmaMinus, isPiDaughter); + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, isSigmaMinus, isPiDaughter); // 1 - tracks with right ITS, TPC, TOF signals if (motherTrack.has_collision() && motherTrack.hasITS() && !motherTrack.hasTPC() && !motherTrack.hasTOF() && daughterTrack.hasITS() && daughterTrack.hasTPC()) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtDaughter, isSigmaMinus, isPiDaughter); - + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, isSigmaMinus, isPiDaughter); } else { continue; } @@ -453,16 +451,14 @@ struct sigmaminustask { if (motherTrack.itsNCls() < 6 && motherTrack.itsNClsInnerBarrel() == 3 && motherTrack.itsChi2NCl() < 36) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtDaughter, isSigmaMinus, isPiDaughter); - + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, isSigmaMinus, isPiDaughter); } else { continue; } if (motherTrack.pt() > KBminPtMoth) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtDaughter, isSigmaMinus, isPiDaughter); - + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, isSigmaMinus, isPiDaughter); } else { continue; } @@ -471,8 +467,7 @@ struct sigmaminustask { if (daughterTrack.itsNClsInnerBarrel() == 0 && daughterTrack.itsNCls() < 4 && daughterTrack.tpcNClsCrossedRows() > 0.8 * daughterTrack.tpcNClsFindable() && daughterTrack.tpcNClsFound() > KBnTPCClusMinDaug) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtDaughter, isSigmaMinus, isPiDaughter); - + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, isSigmaMinus, isPiDaughter); } else { continue; } @@ -480,8 +475,7 @@ struct sigmaminustask { // 5 - geometric cuts: eta if (std::abs(motherTrack.eta()) < KBetaMax && std::abs(daughterTrack.eta()) < KBetaMax) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtDaughter, isSigmaMinus, isPiDaughter); - + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, isSigmaMinus, isPiDaughter); } else { continue; } @@ -489,8 +483,7 @@ struct sigmaminustask { // 6 - geometric cuts: phi difference if (std::abs(motherTrack.phi() - daughterTrack.phi()) * radToDeg < KBmaxPhiDiff) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtDaughter, isSigmaMinus, isPiDaughter); - + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, isSigmaMinus, isPiDaughter); } else { continue; } @@ -498,8 +491,7 @@ struct sigmaminustask { // 7 - radius cut if (recRadius > KBradiusCut) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtDaughter, isSigmaMinus, isPiDaughter); - + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, isSigmaMinus, isPiDaughter); } else { continue; } @@ -508,15 +500,15 @@ struct sigmaminustask { auto collision = motherTrack.template collision_as(); if (!(std::abs(collision.posZ()) > cutzvertex || !collision.sel8())) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtDaughter, isSigmaMinus, isPiDaughter); - + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, isSigmaMinus, isPiDaughter); + } else { + continue; } // 9 - TOF daughter presence if (daughterTrack.hasTOF()) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtDaughter, isSigmaMinus, isPiDaughter); - + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, isSigmaMinus, isPiDaughter); } } } From 539046d0c4bbe4af725bce88d279d44767b9d25f Mon Sep 17 00:00:00 2001 From: morgmatt Date: Thu, 7 Aug 2025 09:59:01 +0200 Subject: [PATCH 12/18] Added qt cuts in dataProcess to reduce size of slimkinkcands table for slim derived data usage. Added fake cluster analysis for findable tracks --- .../Strangeness/sigmaminustask.cxx | 88 +++++++++++++++---- 1 file changed, 73 insertions(+), 15 deletions(-) diff --git a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx index 5ff82b381ea..b504e91b22a 100644 --- a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx @@ -47,14 +47,19 @@ struct sigmaminustask { Configurable cutzvertex{"cutzvertex", 10.0f, "Accepted z-vertex range (cm)"}; Configurable cutNSigmaPi{"cutNSigmaPi", 4, "NSigmaTPCPion"}; Configurable cutRapMotherMC{"cutRapMotherMC", 1.0f, "Rapidity cut for mother Sigma in MC"}; + Configurable cutMinQtAP{"cutMinQtAP", 0.15f, "Minimum Qt for Armenteros-Podolanski cut"}; + Configurable cutMaxQtAP{"cutMaxQtAP", 0.20f, "Maximum Qt for Armenteros-Podolanski cut"}; Configurable fillOutputTree{"fillOutputTree", true, "If true, fill the output tree with Kink candidates"}; + // Configurables for findable tracks (see kinkBuilder.cxx) Configurable KBminPtMoth{"KBminPtMoth", 0.5, "Minimum pT of the mother"}; Configurable KBmaxPhiDiff{"KBmaxPhiDiff", 100, "Max phi difference between the kink daughter and the mother"}; Configurable KBetaMax{"KBetaMax", 1., "Max eta for both mother and daughter"}; Configurable KBnTPCClusMinDaug{"KBnTPCClusMinDaug", 80, "Min daug NTPC clusters"}; Configurable KBradiusCut{"KBradiusCut", 19.6213f, "Min reconstructed decay radius of the mother"}; + Configurable KBdcaMothPv{"KBdcaMothPv", 0.1f, "Max DCA of the mother to PV"}; + Configurable KBdcaDaugPv{"KBdcaDaugPv", 0.1f, "Max DCA of the daughter to PV"}; Preslice mPerCol = aod::track::collisionId; @@ -81,6 +86,7 @@ struct sigmaminustask { const AxisSpec boolAxis{2, -0.5, 1.5, "Boolean value"}; const AxisSpec filtersAxis{10, -0.5, 9.5, "Filter index"}; + const AxisSpec fakeITSAxis{9, -2.5, 6.5, "Fake ITS cluster layer"}; // Event selection rEventSelection.add("hVertexZRec", "hVertexZRec", {HistType::kTH1F, {vertexZAxis}}); @@ -114,6 +120,7 @@ struct sigmaminustask { std::vector filterLabels = {"Initial", "ITS/TPC present", "Moth ITS", "Moth p_{T}", "Daug ITS+TPC", "#eta", "#Delta#phi", "Radius", "Collision", "Daug TOF"}; // Add findable Sigma histograms + rFindable.add("hfakeITSfindable", "hfakeITSfindable", {HistType::kTH1F, {fakeITSAxis}}); rFindable.add("hFilterIndex", "hFilterIndex", {HistType::kTH1F, {filtersAxis}}); rFindable.add("h2MCRadiusFilterIndex", "h2MCRadiusFilterIndex", {HistType::kTH2F, {filtersAxis, radiusAxis}}); rFindable.add("h2RecRadiusFilterIndex", "h2RecRadiusFilterIndex", {HistType::kTH2F, {filtersAxis, radiusAxis}}); @@ -134,9 +141,30 @@ struct sigmaminustask { rFindable.add("h2PtFilter_plus_protonkink", "h2PtFilter_plus_protonkink", {HistType::kTH2F, {filtersAxis, ptUnsignedAxis}}); rFindable.add("h2PtFilter_plus_pikink", "h2PtFilter_plus_pikink", {HistType::kTH2F, {filtersAxis, ptUnsignedAxis}}); rFindable.add("h2PtFilter_minus_pikink", "h2PtFilter_minus_pikink", {HistType::kTH2F, {filtersAxis, ptUnsignedAxis}}); + + rFindable.add("h2PtDaugFilter_plus_protonkink", "h2PtDaugFilter_plus_protonkink", {HistType::kTH2F, {filtersAxis, ptUnsignedAxis}}); + rFindable.add("h2PtDaugFilter_plus_pikink", "h2PtDaugFilter_plus_pikink", {HistType::kTH2F, {filtersAxis, ptUnsignedAxis}}); + rFindable.add("h2PtDaugFilter_minus_pikink", "h2PtDaugFilter_minus_pikink", {HistType::kTH2F, {filtersAxis, ptUnsignedAxis}}); + } } + float alphaAP(const std::array& momMother, const std::array& momKink) + { + std::array momMissing = {momMother[0] - momKink[0], momMother[1] - momKink[1], momMother[2] - momKink[2]}; + float lQlP = std::inner_product(momMother.begin(), momMother.end(), momKink.begin(), 0.f); + float lQlN = std::inner_product(momMother.begin(), momMother.end(), momMissing.begin(), 0.f); + return (lQlP - lQlN) / (lQlP + lQlN); + } + + float qtAP(const std::array& momMother, const std::array& momKink) + { + float dp = std::inner_product(momMother.begin(), momMother.end(), momKink.begin(), 0.f); + float p2V0 = std::inner_product(momMother.begin(), momMother.end(), momMother.begin(), 0.f); + float p2A = std::inner_product(momKink.begin(), momKink.end(), momKink.begin(), 0.f); + return std::sqrt(p2A - dp * dp / p2V0); + } + void processData(CollisionsFull::iterator const& collision, aod::KinkCands const& KinkCands, TracksFull const&) { if (std::abs(collision.posZ()) > cutzvertex || !collision.sel8()) { @@ -151,6 +179,11 @@ struct sigmaminustask { continue; } + float qt = qtAP(std::array{kinkCand.pxMoth(), kinkCand.pyMoth(), kinkCand.pzMoth()}, std::array{kinkCand.pxDaug(), kinkCand.pyDaug(), kinkCand.pzDaug()}); + if (qt < cutMinQtAP || qt > cutMaxQtAP) { + continue; + } + rSigmaMinus.fill(HIST("h2MassSigmaMinusPt"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.mSigmaMinus()); rSigmaMinus.fill(HIST("h2SigmaMassVsXiMass"), kinkCand.mXiMinus(), kinkCand.mSigmaMinus()); rSigmaMinus.fill(HIST("h2NSigmaPiPt"), kinkCand.mothSign() * kinkCand.ptMoth(), dauTrack.tpcNSigmaPi()); @@ -316,7 +349,8 @@ struct sigmaminustask { PROCESS_SWITCH(sigmaminustask, processMC, "MC processing", false); - void fillFindableHistograms(int filterIndex, float mcRadius, float recRadius, float pt, bool isSigmaMinus, bool isPiDaughter) { + void fillFindableHistograms(int filterIndex, float mcRadius, float recRadius, float ptMoth, float ptDaug, bool isSigmaMinus, bool isPiDaughter) + { rFindable.fill(HIST("hFilterIndex"), filterIndex); rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius); @@ -324,21 +358,25 @@ struct sigmaminustask { if (isPiDaughter) { if (isSigmaMinus) { rFindable.fill(HIST("h2MCRadiusFilter_minus_pikink"), filterIndex, mcRadius); - rFindable.fill(HIST("h2PtFilter_minus_pikink"), filterIndex, pt); + rFindable.fill(HIST("h2PtFilter_minus_pikink"), filterIndex, ptMoth); + rFindable.fill(HIST("h2PtDaugFilter_minus_pikink"), filterIndex, ptDaug); } else { rFindable.fill(HIST("h2MCRadiusFilter_plus_pikink"), filterIndex, mcRadius); - rFindable.fill(HIST("h2PtFilter_plus_pikink"), filterIndex, pt); + rFindable.fill(HIST("h2PtFilter_plus_pikink"), filterIndex, ptMoth); + rFindable.fill(HIST("h2PtDaugFilter_plus_pikink"), filterIndex, ptDaug); } } else { if (!isSigmaMinus) { rFindable.fill(HIST("h2MCRadiusFilter_plus_protonkink"), filterIndex, mcRadius); - rFindable.fill(HIST("h2PtFilter_plus_protonkink"), filterIndex, pt); + rFindable.fill(HIST("h2PtFilter_plus_protonkink"), filterIndex, ptMoth); + rFindable.fill(HIST("h2PtDaugFilter_plus_protonkink"), filterIndex, ptDaug); } } } void processFindable(aod::KinkCands const& kinkCands, aod::McTrackLabels const& trackLabelsMC, - TracksFull const& tracks, aod::McParticles const&, CollisionsFullMC const&){ + TracksFull const& tracks, aod::McParticles const&, CollisionsFullMC const&) + { // A - generated findable track pairs map: mcMother.globalIndex() -> (motherTrack.globalIndex(), daughterTrack.globalIndex()) std::unordered_map> allCandsIndices; @@ -434,15 +472,28 @@ struct sigmaminustask { float recPtDaughter = daughterTrack.pt(); float recPtMother = motherTrack.pt(); + // Check for detector mismatches in ITS mother tracks + auto mask_value = mcLabMoth.mcMask(); + int mismatchITS_index = -1; + + for (int i = 0; i < 7; ++i) { // ITS has layers 0-6, bit ON means mismatch + if ((mask_value & (1 << i)) != 0) { + mismatchITS_index = i; + break; + } + } + // Define filter index and progressively apply kinkbuilder cuts to track pairs int filterIndex = 0; - fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, isSigmaMinus, isPiDaughter); + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); // 1 - tracks with right ITS, TPC, TOF signals if (motherTrack.has_collision() && motherTrack.hasITS() && !motherTrack.hasTPC() && !motherTrack.hasTOF() && daughterTrack.hasITS() && daughterTrack.hasTPC()) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, isSigmaMinus, isPiDaughter); + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); + rFindable.fill(HIST("hfakeITSfindable"), mismatchITS_index); + } else { continue; } @@ -451,14 +502,14 @@ struct sigmaminustask { if (motherTrack.itsNCls() < 6 && motherTrack.itsNClsInnerBarrel() == 3 && motherTrack.itsChi2NCl() < 36) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, isSigmaMinus, isPiDaughter); + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); } else { continue; } if (motherTrack.pt() > KBminPtMoth) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, isSigmaMinus, isPiDaughter); + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); } else { continue; } @@ -467,7 +518,7 @@ struct sigmaminustask { if (daughterTrack.itsNClsInnerBarrel() == 0 && daughterTrack.itsNCls() < 4 && daughterTrack.tpcNClsCrossedRows() > 0.8 * daughterTrack.tpcNClsFindable() && daughterTrack.tpcNClsFound() > KBnTPCClusMinDaug) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, isSigmaMinus, isPiDaughter); + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); } else { continue; } @@ -475,7 +526,7 @@ struct sigmaminustask { // 5 - geometric cuts: eta if (std::abs(motherTrack.eta()) < KBetaMax && std::abs(daughterTrack.eta()) < KBetaMax) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, isSigmaMinus, isPiDaughter); + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); } else { continue; } @@ -483,7 +534,7 @@ struct sigmaminustask { // 6 - geometric cuts: phi difference if (std::abs(motherTrack.phi() - daughterTrack.phi()) * radToDeg < KBmaxPhiDiff) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, isSigmaMinus, isPiDaughter); + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); } else { continue; } @@ -491,7 +542,7 @@ struct sigmaminustask { // 7 - radius cut if (recRadius > KBradiusCut) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, isSigmaMinus, isPiDaughter); + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); } else { continue; } @@ -500,7 +551,7 @@ struct sigmaminustask { auto collision = motherTrack.template collision_as(); if (!(std::abs(collision.posZ()) > cutzvertex || !collision.sel8())) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, isSigmaMinus, isPiDaughter); + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); } else { continue; } @@ -508,7 +559,7 @@ struct sigmaminustask { // 9 - TOF daughter presence if (daughterTrack.hasTOF()) { filterIndex += 1; - fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, isSigmaMinus, isPiDaughter); + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); } } } @@ -521,3 +572,10 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) return WorkflowSpec{ adaptAnalysisTask(cfgc)}; } + +// TO DO: +// 2) Histogram of fake its tracks of mothers among findable not found candidates +// 2.1) -1 index for not fake, 0,1,2,3 etc. for fake ITS tracks, with index as layer of fake cluster +// 3) Add DCA mother and daughter to PV using track propagation from hyperhelium or hypertriton tasks +// 3.1) Do it for each findable, add cuts on DCA to the filterIndex histos to understand what's the most effective cut +// 4) Open PR \ No newline at end of file From b7b652a88455551dc8320b7189bfad82fc227fd5 Mon Sep 17 00:00:00 2001 From: morgmatt Date: Thu, 7 Aug 2025 12:32:29 +0200 Subject: [PATCH 13/18] Added dca computation and filters analysis to process findable. Added necessary track propagation dependencies --- .../Strangeness/sigmaminustask.cxx | 179 +++++++++++++----- 1 file changed, 128 insertions(+), 51 deletions(-) diff --git a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx index b504e91b22a..209c0fd6f65 100644 --- a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx @@ -15,12 +15,19 @@ #include "PWGLF/DataModel/LFKinkDecayTables.h" +#include "Common/Core/trackUtilities.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/PIDResponse.h" +#include "CCDB/BasicCCDBManager.h" +#include "DataFormatsParameters/GRPMagField.h" +#include "DataFormatsParameters/GRPObject.h" +#include "DetectorsBase/GeometryManager.h" +#include "DetectorsBase/Propagator.h" #include "Framework/AnalysisTask.h" #include "Framework/runDataProcessing.h" #include "ReconstructionDataFormats/PID.h" +#include "ReconstructionDataFormats/Track.h" using namespace o2; using namespace o2::framework; @@ -52,21 +59,43 @@ struct sigmaminustask { Configurable fillOutputTree{"fillOutputTree", true, "If true, fill the output tree with Kink candidates"}; - // Configurables for findable tracks (see kinkBuilder.cxx) - Configurable KBminPtMoth{"KBminPtMoth", 0.5, "Minimum pT of the mother"}; - Configurable KBmaxPhiDiff{"KBmaxPhiDiff", 100, "Max phi difference between the kink daughter and the mother"}; - Configurable KBetaMax{"KBetaMax", 1., "Max eta for both mother and daughter"}; - Configurable KBnTPCClusMinDaug{"KBnTPCClusMinDaug", 80, "Min daug NTPC clusters"}; - Configurable KBradiusCut{"KBradiusCut", 19.6213f, "Min reconstructed decay radius of the mother"}; - Configurable KBdcaMothPv{"KBdcaMothPv", 0.1f, "Max DCA of the mother to PV"}; - Configurable KBdcaDaugPv{"KBdcaDaugPv", 0.1f, "Max DCA of the daughter to PV"}; + // Configurables for findable tracks (kinkBuilder.cxx efficiency) + Configurable minPtMothKB{"minPtMothKB", 0.5f, "Minimum pT of the mother"}; + Configurable maxPhiDiffKB{"maxPhiDiffKB", 100.0f, "Max phi difference between the kink daughter and the mother"}; + Configurable maxZDiffKB{"maxZDiffKB", 20.0f, "Max z difference between the kink daughter and the mother"}; + Configurable etaMaxKB{"etaMaxKB", 1.0f, "Max eta for both mother and daughter"}; + Configurable nTPCClusMinDaugKB{"nTPCClusMinDaugKB", 80, "Min daug NTPC clusters"}; + Configurable radiusCutKB{"radiusCutKB", 19.6213f, "Min reconstructed decay radius of the mother"}; + Configurable maxDcaMothPvKB{"maxDcaMothPvKB", 0.1f, "Max DCA of the mother to PV"}; + Configurable minDcaDaugPvKB{"minDcaDaugPvKB", 0.1f, "Min DCA of the daughter to PV"}; Preslice mPerCol = aod::track::collisionId; + // Constants float radToDeg = o2::constants::math::Rad2Deg; + // Services CCDB + Service ccdb; + Configurable ccdbPath{"ccdbPath", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; + Configurable lutPath{"lutPath", "GLO/Param/MatLUT", "Path of the Lut parametrization"}; + Configurable cfgMaterialCorrection{"cfgMaterialCorrection", static_cast(o2::base::Propagator::MatCorrType::USEMatCorrLUT), "Type of material correction"}; + + // Runtime variables + int mRunNumber = 0; + float mBz = 0.0f; + o2::base::MatLayerCylSet* matLUT = nullptr; + void init(InitContext const&) { + // Initialize CCDB + ccdb->setURL(ccdbPath); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setFatalWhenNull(true); + + matLUT = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->get("GLO/Param/MatLUT")); + // Axes const AxisSpec ptAxis{100, -10, 10, "#it{p}_{T} (GeV/#it{c})"}; const AxisSpec ptUnsignedAxis{100, 0, 10, "#it{p}_{T} (GeV/#it{c})"}; @@ -76,7 +105,7 @@ struct sigmaminustask { const AxisSpec xiMassAxis{100, 1.2, 1.6, "m_{#Xi} (GeV/#it{c}^{2})"}; const AxisSpec pdgAxis{10001, -5000, 5000, "PDG code"}; const AxisSpec vertexZAxis{100, -15., 15., "vrtx_{Z} [cm]"}; - const AxisSpec dcaMothAxis{100, 0, 0.2, "DCA [cm]"}; + const AxisSpec dcaMothAxis{100, 0, 0.03, "DCA [cm]"}; const AxisSpec dcaDaugAxis{200, 0, 20, "DCA [cm]"}; const AxisSpec radiusAxis{100, -1, 40, "Decay radius [cm]"}; @@ -85,8 +114,8 @@ struct sigmaminustask { const AxisSpec radiusResolutionAxis{100, -0.5, 0.5, "(r_{rec} - r_{gen}) / r_{gen}"}; const AxisSpec boolAxis{2, -0.5, 1.5, "Boolean value"}; - const AxisSpec filtersAxis{10, -0.5, 9.5, "Filter index"}; - const AxisSpec fakeITSAxis{9, -2.5, 6.5, "Fake ITS cluster layer"}; + const AxisSpec filtersAxis{12, -0.5, 11.5, "Filter index"}; + const AxisSpec fakeITSAxis{8, -1.5, 6.5, "Fake ITS cluster layer"}; // Event selection rEventSelection.add("hVertexZRec", "hVertexZRec", {HistType::kTH1F, {vertexZAxis}}); @@ -117,8 +146,8 @@ struct sigmaminustask { } if (doprocessFindable) { - std::vector filterLabels = {"Initial", "ITS/TPC present", "Moth ITS", "Moth p_{T}", "Daug ITS+TPC", "#eta", "#Delta#phi", "Radius", "Collision", "Daug TOF"}; - + std::vector filterLabels = {"Initial", "ITS/TPC present", "ITS/TPC quality", "Moth p_{T}", "min #eta", "max #Delta#phi", "max #Delta Z", "max DCAmoth", "min DCAdaug", "min Radius", "sel8 coll", "Daug TOF"}; + // Add findable Sigma histograms rFindable.add("hfakeITSfindable", "hfakeITSfindable", {HistType::kTH1F, {fakeITSAxis}}); rFindable.add("hFilterIndex", "hFilterIndex", {HistType::kTH1F, {filtersAxis}}); @@ -146,7 +175,28 @@ struct sigmaminustask { rFindable.add("h2PtDaugFilter_plus_pikink", "h2PtDaugFilter_plus_pikink", {HistType::kTH2F, {filtersAxis, ptUnsignedAxis}}); rFindable.add("h2PtDaugFilter_minus_pikink", "h2PtDaugFilter_minus_pikink", {HistType::kTH2F, {filtersAxis, ptUnsignedAxis}}); + rFindable.add("h2DCAMothPt_protonkink", "h2DCAMothPt_protonkink", {HistType::kTH2F, {ptAxis, dcaMothAxis}}); + rFindable.add("h2DCADaugPt_protonkink", "h2DCADaugPt_protonkink", {HistType::kTH2F, {ptAxis, dcaDaugAxis}}); + rFindable.add("h2DCAMothPt_pikink", "h2DCAMothPt_pikink", {HistType::kTH2F, {ptAxis, dcaMothAxis}}); + rFindable.add("h2DCADaugPt_pikink", "h2DCADaugPt_pikink", {HistType::kTH2F, {ptAxis, dcaDaugAxis}}); + } + } + + void initCCDB(aod::BCs::iterator const& bc) + { + if (mRunNumber == bc.runNumber()) { + return; + } + mRunNumber = bc.runNumber(); + LOG(info) << "Initializing CCDB for run " << mRunNumber; + o2::parameters::GRPMagField* grpmag = ccdb->getForRun(grpmagPath, mRunNumber); + o2::base::Propagator::initFieldFromGRP(grpmag); + mBz = grpmag->getNominalL3Field(); + if (!matLUT) { + matLUT = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->get(lutPath)); } + o2::base::Propagator::Instance()->setMatLUT(matLUT); + LOG(info) << "Task initialized for run " << mRunNumber << " with magnetic field " << mBz << " kZG"; } float alphaAP(const std::array& momMother, const std::array& momKink) @@ -346,7 +396,6 @@ struct sigmaminustask { } } } - PROCESS_SWITCH(sigmaminustask, processMC, "MC processing", false); void fillFindableHistograms(int filterIndex, float mcRadius, float recRadius, float ptMoth, float ptDaug, bool isSigmaMinus, bool isPiDaughter) @@ -375,7 +424,7 @@ struct sigmaminustask { } void processFindable(aod::KinkCands const& kinkCands, aod::McTrackLabels const& trackLabelsMC, - TracksFull const& tracks, aod::McParticles const&, CollisionsFullMC const&) + TracksFull const& tracks, aod::McParticles const&, CollisionsFullMC const&, aod::BCs const&) { // A - generated findable track pairs map: mcMother.globalIndex() -> (motherTrack.globalIndex(), daughterTrack.globalIndex()) std::unordered_map> allCandsIndices; @@ -446,32 +495,28 @@ struct sigmaminustask { if (trackIndices.second == -1 || trackIndices.first == -1) { continue; } + // Retrieve mother and daughter tracks and mcParticles auto motherTrack = tracks.rawIteratorAt(trackIndices.first); auto daughterTrack = tracks.rawIteratorAt(trackIndices.second); auto mcLabMoth = trackLabelsMC.rawIteratorAt(motherTrack.globalIndex()); auto mcLabDaug = trackLabelsMC.rawIteratorAt(daughterTrack.globalIndex()); - if (!mcLabMoth.has_mcParticle() || !mcLabDaug.has_mcParticle()) { - continue; - } auto mcMother = mcLabMoth.mcParticle_as(); auto mcDaughter = mcLabDaug.mcParticle_as(); - // Compute mass and radii + // Compute useful quantities for histograms bool isSigmaMinus = (std::abs(mcMother.pdgCode()) == 3112); bool isPiDaughter = (std::abs(mcDaughter.pdgCode()) == 211); - - //float mcMass = std::sqrt(mcMother.e() * mcMother.e() - mcMother.p() * mcMother.p()); - //int sigmaSign = mcMother.pdgCode() > 0 ? 1 : -1; + int sigmaSign = mcMother.pdgCode() > 0 ? 1 : -1; + float recPtDaughter = daughterTrack.pt(); + float recPtMother = motherTrack.pt(); float mcRadius = std::sqrt((mcMother.vx() - mcDaughter.vx()) * (mcMother.vx() - mcDaughter.vx()) + (mcMother.vy() - mcDaughter.vy()) * (mcMother.vy() - mcDaughter.vy())); float recRadius = -1.0; if (findableToKinkCand.find(mcMother.globalIndex()) != findableToKinkCand.end()) { auto kinkCand = kinkCands.rawIteratorAt(findableToKinkCand[mcMother.globalIndex()]); recRadius = std::sqrt(kinkCand.xDecVtx() * kinkCand.xDecVtx() + kinkCand.yDecVtx() * kinkCand.yDecVtx()); } - float recPtDaughter = daughterTrack.pt(); - float recPtMother = motherTrack.pt(); - + // Check for detector mismatches in ITS mother tracks auto mask_value = mcLabMoth.mcMask(); int mismatchITS_index = -1; @@ -487,68 +532,107 @@ struct sigmaminustask { int filterIndex = 0; fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); - // 1 - tracks with right ITS, TPC, TOF signals + // 1 - tracks with right ITS, TPC, TOF signals if (motherTrack.has_collision() && motherTrack.hasITS() && !motherTrack.hasTPC() && !motherTrack.hasTOF() && daughterTrack.hasITS() && daughterTrack.hasTPC()) { filterIndex += 1; fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); rFindable.fill(HIST("hfakeITSfindable"), mismatchITS_index); - } else { continue; } - // 2, 3 - mother track ITS properties - if (motherTrack.itsNCls() < 6 && - motherTrack.itsNClsInnerBarrel() == 3 && motherTrack.itsChi2NCl() < 36) { + // 2 - moth+daug track quality cuts + bool motherGoodITS = motherTrack.hasITS() && motherTrack.itsNCls() < 6 && motherTrack.itsNClsInnerBarrel() == 3 && motherTrack.itsChi2NCl() < 36; + bool daughterGoodITSTPC = daughterTrack.hasITS() && daughterTrack.hasTPC() && daughterTrack.itsNClsInnerBarrel() == 0 && + daughterTrack.itsNCls() < 4 && daughterTrack.tpcNClsCrossedRows() > 0.8 * daughterTrack.tpcNClsFindable() && daughterTrack.tpcNClsFound() > nTPCClusMinDaugKB; + if (motherGoodITS && daughterGoodITSTPC) { filterIndex += 1; fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); } else { continue; } - if (motherTrack.pt() > KBminPtMoth) { + // 3 - mother track min pT + if (motherTrack.pt() > minPtMothKB) { filterIndex += 1; fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); } else { continue; } - // 4 - daughter track ITS+TPC properties - if (daughterTrack.itsNClsInnerBarrel() == 0 && daughterTrack.itsNCls() < 4 && - daughterTrack.tpcNClsCrossedRows() > 0.8 * daughterTrack.tpcNClsFindable() && daughterTrack.tpcNClsFound() > KBnTPCClusMinDaug) { + // 4 - geometric cuts: eta + if (std::abs(motherTrack.eta()) < etaMaxKB && std::abs(daughterTrack.eta()) < etaMaxKB) { filterIndex += 1; fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); } else { continue; } - // 5 - geometric cuts: eta - if (std::abs(motherTrack.eta()) < KBetaMax && std::abs(daughterTrack.eta()) < KBetaMax) { + // 5 - geometric cuts: phi difference + if (std::abs(motherTrack.phi() - daughterTrack.phi()) * radToDeg < maxPhiDiffKB) { filterIndex += 1; fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); } else { continue; } - // 6 - geometric cuts: phi difference - if (std::abs(motherTrack.phi() - daughterTrack.phi()) * radToDeg < KBmaxPhiDiff) { + // DCA calculation: initialization CCDB + auto collision = motherTrack.template collision_as(); + auto bc = collision.template bc_as(); + initCCDB(bc); + const o2::math_utils::Point3D collVtx{collision.posX(), collision.posY(), collision.posZ()}; + o2::track::TrackParCov trackParCovMoth = getTrackParCov(motherTrack); + o2::track::TrackParCov trackParCovDaug = getTrackParCov(daughterTrack); + + // get DCA to PV for mother and daughter tracks + std::array dcaInfoMoth; + o2::base::Propagator::Instance()->propagateToDCA(collVtx, trackParCovMoth, mBz, 2.f, static_cast(cfgMaterialCorrection.value), &dcaInfoMoth); + std::array dcaInfoDaug; + o2::base::Propagator::Instance()->propagateToDCA(collVtx, trackParCovDaug, mBz, 2.f, static_cast(cfgMaterialCorrection.value), &dcaInfoDaug); + float dcaXYMother = std::abs(dcaInfoMoth[0]); + float dcaXYDaughter = std::abs(dcaInfoDaug[0]); + if (isPiDaughter) { + rFindable.fill(HIST("h2DCAMothPt_pikink"), sigmaSign * recPtMother, dcaXYMother); + rFindable.fill(HIST("h2DCADaugPt_pikink"), sigmaSign * recPtDaughter, dcaXYDaughter); + } else { + rFindable.fill(HIST("h2DCAMothPt_protonkink"), sigmaSign * recPtMother, dcaXYMother); + rFindable.fill(HIST("h2DCADaugPt_protonkink"), sigmaSign * recPtDaughter, dcaXYDaughter); + } + + // 6 - max Z difference + if (std::abs(trackParCovMoth.getZ() - trackParCovDaug.getZ()) < maxZDiffKB) { + filterIndex += 1; + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); + } else { + continue; + } + + // 7 - DCA mother + if (dcaXYMother < maxDcaMothPvKB) { filterIndex += 1; fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); } else { continue; } - - // 7 - radius cut - if (recRadius > KBradiusCut) { + + // 8 - DCA daughter + if (dcaXYDaughter > minDcaDaugPvKB) { filterIndex += 1; fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); } else { continue; } - // 8 - collision selection - auto collision = motherTrack.template collision_as(); + // 9 - radius cut + if (recRadius > radiusCutKB) { + filterIndex += 1; + fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); + } else { + continue; + } + + // 10 - collision selection if (!(std::abs(collision.posZ()) > cutzvertex || !collision.sel8())) { filterIndex += 1; fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); @@ -556,7 +640,7 @@ struct sigmaminustask { continue; } - // 9 - TOF daughter presence + // 11 - TOF daughter presence if (daughterTrack.hasTOF()) { filterIndex += 1; fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); @@ -572,10 +656,3 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) return WorkflowSpec{ adaptAnalysisTask(cfgc)}; } - -// TO DO: -// 2) Histogram of fake its tracks of mothers among findable not found candidates -// 2.1) -1 index for not fake, 0,1,2,3 etc. for fake ITS tracks, with index as layer of fake cluster -// 3) Add DCA mother and daughter to PV using track propagation from hyperhelium or hypertriton tasks -// 3.1) Do it for each findable, add cuts on DCA to the filterIndex histos to understand what's the most effective cut -// 4) Open PR \ No newline at end of file From c71bd5cf081ebdcef0e8fc3cb9f235d7373a142f Mon Sep 17 00:00:00 2001 From: morgmatt Date: Thu, 7 Aug 2025 14:42:37 +0200 Subject: [PATCH 14/18] Removed one histogram for BC id comparison. Corrected label for filter indices axis --- PWGLF/TableProducer/Strangeness/sigmaminustask.cxx | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx index 209c0fd6f65..869dc2967b1 100644 --- a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx @@ -141,12 +141,11 @@ struct sigmaminustask { // BC ID comparison histograms rSigmaMinus.add("hMcCollIdCoherence", "McCollId (coll == daug)", {HistType::kTH1F, {boolAxis}}); rSigmaMinus.add("h2CollId_BCId", "(McCollId coherence) vs (EvSelBC == McBC)", {HistType::kTH2F, {boolAxis, boolAxis}}); - rSigmaMinus.add("h2BCId_comp1", "(BC == McBC) vs (BC == EvSelBC)", {HistType::kTH2F, {boolAxis, boolAxis}}); - rSigmaMinus.add("h2BCId_comp2", "(McBC == EvSelBC) vs (BC == EvSelBC)", {HistType::kTH2F, {boolAxis, boolAxis}}); + rSigmaMinus.add("h2BCId_comp", "(McBC == EvSelBC) vs (BC == EvSelBC)", {HistType::kTH2F, {boolAxis, boolAxis}}); } if (doprocessFindable) { - std::vector filterLabels = {"Initial", "ITS/TPC present", "ITS/TPC quality", "Moth p_{T}", "min #eta", "max #Delta#phi", "max #Delta Z", "max DCAmoth", "min DCAdaug", "min Radius", "sel8 coll", "Daug TOF"}; + std::vector filterLabels = {"Initial", "ITS/TPC present", "ITS/TPC quality", "Moth p_{T}", "max #eta", "max #Delta#phi", "max #Delta Z", "max DCAmoth", "min DCAdaug", "min Radius", "sel8 coll", "Daug TOF"}; // Add findable Sigma histograms rFindable.add("hfakeITSfindable", "hfakeITSfindable", {HistType::kTH1F, {fakeITSAxis}}); @@ -314,14 +313,13 @@ struct sigmaminustask { } // Check bunch crossing ID coherence auto mcCollision = mcTrackPiDau.template mcCollision_as(); - bool BCId_vs_MCBCId = collision.bcId() == mcCollision.bcId(); + //bool BCId_vs_MCBCId = collision.bcId() == mcCollision.bcId(); bool BCId_vs_EvSel = collision.bcId() == collision.foundBCId(); bool EvSel_vs_MCBCId = collision.foundBCId() == mcCollision.bcId(); rSigmaMinus.fill(HIST("hMcCollIdCoherence"), static_cast(mcCollisionIdCheck)); rSigmaMinus.fill(HIST("h2CollId_BCId"), static_cast(mcCollisionIdCheck), static_cast(EvSel_vs_MCBCId)); - rSigmaMinus.fill(HIST("h2BCId_comp1"), static_cast(BCId_vs_MCBCId), static_cast(BCId_vs_EvSel)); - rSigmaMinus.fill(HIST("h2BCId_comp2"), static_cast(EvSel_vs_MCBCId), static_cast(BCId_vs_EvSel)); + rSigmaMinus.fill(HIST("h2BCId_comp"), static_cast(EvSel_vs_MCBCId), static_cast(BCId_vs_EvSel)); rSigmaMinus.fill(HIST("h2MassPtMCRec"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.mSigmaMinus()); rSigmaMinus.fill(HIST("h2MassResolution"), kinkCand.mothSign() * kinkCand.ptMoth(), (kinkCand.mSigmaMinus() - MotherMassMC) / MotherMassMC); From 45fd85a1b71ae884dcbc8fc2cfdd37ed4119de31 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Thu, 7 Aug 2025 14:04:27 +0000 Subject: [PATCH 15/18] Please consider the following formatting changes --- .../Strangeness/sigmaminustask.cxx | 79 +++++++++---------- 1 file changed, 39 insertions(+), 40 deletions(-) diff --git a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx index 869dc2967b1..9927fce36cb 100644 --- a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx @@ -83,7 +83,7 @@ struct sigmaminustask { // Runtime variables int mRunNumber = 0; - float mBz = 0.0f; + float mBz = 0.0f; o2::base::MatLayerCylSet* matLUT = nullptr; void init(InitContext const&) @@ -157,15 +157,15 @@ struct sigmaminustask { auto h2MCRadiusFilterIndex = rFindable.get(HIST("h2MCRadiusFilterIndex")); auto h2RecRadiusFilterIndex = rFindable.get(HIST("h2RecRadiusFilterIndex")); for (size_t i = 0; i < filterLabels.size(); ++i) { - hFilterIndex->GetXaxis()->SetBinLabel(i + 1, filterLabels[i].c_str()); - h2MCRadiusFilterIndex->GetXaxis()->SetBinLabel(i + 1, filterLabels[i].c_str()); - h2RecRadiusFilterIndex->GetXaxis()->SetBinLabel(i + 1, filterLabels[i].c_str()); + hFilterIndex->GetXaxis()->SetBinLabel(i + 1, filterLabels[i].c_str()); + h2MCRadiusFilterIndex->GetXaxis()->SetBinLabel(i + 1, filterLabels[i].c_str()); + h2RecRadiusFilterIndex->GetXaxis()->SetBinLabel(i + 1, filterLabels[i].c_str()); } rFindable.add("h2MCRadiusFilter_plus_protonkink", "h2MCRadiusFilter_plus_protonkink", {HistType::kTH2F, {filtersAxis, radiusAxis}}); rFindable.add("h2MCRadiusFilter_plus_pikink", "h2MCRadiusFilter_plus_pikink", {HistType::kTH2F, {filtersAxis, radiusAxis}}); rFindable.add("h2MCRadiusFilter_minus_pikink", "h2MCRadiusFilter_minus_pikink", {HistType::kTH2F, {filtersAxis, radiusAxis}}); - + rFindable.add("h2PtFilter_plus_protonkink", "h2PtFilter_plus_protonkink", {HistType::kTH2F, {filtersAxis, ptUnsignedAxis}}); rFindable.add("h2PtFilter_plus_pikink", "h2PtFilter_plus_pikink", {HistType::kTH2F, {filtersAxis, ptUnsignedAxis}}); rFindable.add("h2PtFilter_minus_pikink", "h2PtFilter_minus_pikink", {HistType::kTH2F, {filtersAxis, ptUnsignedAxis}}); @@ -252,7 +252,7 @@ struct sigmaminustask { } PROCESS_SWITCH(sigmaminustask, processData, "Data processing", true); - void processMC(CollisionsFullMC const& collisions, aod::KinkCands const& KinkCands, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const& particlesMC, TracksFull const&, aod::McCollisions const&) + void processMC(CollisionsFullMC const& collisions, aod::KinkCands const& KinkCands, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const& particlesMC, TracksFull const&, aod::McCollisions const&) { for (const auto& collision : collisions) { if (std::abs(collision.posZ()) > cutzvertex || !collision.sel8()) { @@ -313,7 +313,7 @@ struct sigmaminustask { } // Check bunch crossing ID coherence auto mcCollision = mcTrackPiDau.template mcCollision_as(); - //bool BCId_vs_MCBCId = collision.bcId() == mcCollision.bcId(); + // bool BCId_vs_MCBCId = collision.bcId() == mcCollision.bcId(); bool BCId_vs_EvSel = collision.bcId() == collision.foundBCId(); bool EvSel_vs_MCBCId = collision.foundBCId() == mcCollision.bcId(); @@ -396,12 +396,12 @@ struct sigmaminustask { } PROCESS_SWITCH(sigmaminustask, processMC, "MC processing", false); - void fillFindableHistograms(int filterIndex, float mcRadius, float recRadius, float ptMoth, float ptDaug, bool isSigmaMinus, bool isPiDaughter) + void fillFindableHistograms(int filterIndex, float mcRadius, float recRadius, float ptMoth, float ptDaug, bool isSigmaMinus, bool isPiDaughter) { rFindable.fill(HIST("hFilterIndex"), filterIndex); rFindable.fill(HIST("h2MCRadiusFilterIndex"), filterIndex, mcRadius); rFindable.fill(HIST("h2RecRadiusFilterIndex"), filterIndex, recRadius); - + if (isPiDaughter) { if (isSigmaMinus) { rFindable.fill(HIST("h2MCRadiusFilter_minus_pikink"), filterIndex, mcRadius); @@ -414,15 +414,15 @@ struct sigmaminustask { } } else { if (!isSigmaMinus) { - rFindable.fill(HIST("h2MCRadiusFilter_plus_protonkink"), filterIndex, mcRadius); - rFindable.fill(HIST("h2PtFilter_plus_protonkink"), filterIndex, ptMoth); - rFindable.fill(HIST("h2PtDaugFilter_plus_protonkink"), filterIndex, ptDaug); + rFindable.fill(HIST("h2MCRadiusFilter_plus_protonkink"), filterIndex, mcRadius); + rFindable.fill(HIST("h2PtFilter_plus_protonkink"), filterIndex, ptMoth); + rFindable.fill(HIST("h2PtDaugFilter_plus_protonkink"), filterIndex, ptDaug); } } } - void processFindable(aod::KinkCands const& kinkCands, aod::McTrackLabels const& trackLabelsMC, - TracksFull const& tracks, aod::McParticles const&, CollisionsFullMC const&, aod::BCs const&) + void processFindable(aod::KinkCands const& kinkCands, aod::McTrackLabels const& trackLabelsMC, + TracksFull const& tracks, aod::McParticles const&, CollisionsFullMC const&, aod::BCs const&) { // A - generated findable track pairs map: mcMother.globalIndex() -> (motherTrack.globalIndex(), daughterTrack.globalIndex()) std::unordered_map> allCandsIndices; @@ -436,7 +436,7 @@ struct sigmaminustask { if (mcParticle.has_daughters() && (std::abs(mcParticle.pdgCode()) == 3112 || std::abs(mcParticle.pdgCode()) == 3222)) { allCandsIndices[mcParticle.globalIndex()] = {track.globalIndex(), -1}; - } + } } for (const auto& track : tracks) { @@ -472,26 +472,25 @@ struct sigmaminustask { auto mcDaughter = mcLabDaug.mcParticle_as(); if (std::abs(mcMother.pdgCode()) != 3112 && std::abs(mcMother.pdgCode()) != 3222) { - continue; + continue; } if (std::abs(mcDaughter.pdgCode()) != 211 && std::abs(mcDaughter.pdgCode()) != 2212) { - continue; + continue; } auto findableIt = allCandsIndices.find(mcMother.globalIndex()); - if (findableIt != allCandsIndices.end() && - findableIt->second.first == motherTrack.globalIndex() && - findableIt->second.second == daughterTrack.globalIndex()) { + if (findableIt != allCandsIndices.end() && + findableIt->second.first == motherTrack.globalIndex() && + findableIt->second.second == daughterTrack.globalIndex()) { findableToKinkCand[mcMother.globalIndex()] = kinkCand.globalIndex(); } } - // C - loop on valid pairs for findable analysis for (const auto& [mcMotherIndex, trackIndices] : allCandsIndices) { if (trackIndices.second == -1 || trackIndices.first == -1) { - continue; + continue; } // Retrieve mother and daughter tracks and mcParticles @@ -514,14 +513,14 @@ struct sigmaminustask { auto kinkCand = kinkCands.rawIteratorAt(findableToKinkCand[mcMother.globalIndex()]); recRadius = std::sqrt(kinkCand.xDecVtx() * kinkCand.xDecVtx() + kinkCand.yDecVtx() * kinkCand.yDecVtx()); } - + // Check for detector mismatches in ITS mother tracks auto mask_value = mcLabMoth.mcMask(); int mismatchITS_index = -1; - - for (int i = 0; i < 7; ++i) { // ITS has layers 0-6, bit ON means mismatch + + for (int i = 0; i < 7; ++i) { // ITS has layers 0-6, bit ON means mismatch if ((mask_value & (1 << i)) != 0) { - mismatchITS_index = i; + mismatchITS_index = i; break; } } @@ -530,25 +529,25 @@ struct sigmaminustask { int filterIndex = 0; fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); - // 1 - tracks with right ITS, TPC, TOF signals + // 1 - tracks with right ITS, TPC, TOF signals if (motherTrack.has_collision() && motherTrack.hasITS() && !motherTrack.hasTPC() && !motherTrack.hasTOF() && daughterTrack.hasITS() && daughterTrack.hasTPC()) { filterIndex += 1; fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); rFindable.fill(HIST("hfakeITSfindable"), mismatchITS_index); } else { - continue; + continue; } - - // 2 - moth+daug track quality cuts + + // 2 - moth+daug track quality cuts bool motherGoodITS = motherTrack.hasITS() && motherTrack.itsNCls() < 6 && motherTrack.itsNClsInnerBarrel() == 3 && motherTrack.itsChi2NCl() < 36; bool daughterGoodITSTPC = daughterTrack.hasITS() && daughterTrack.hasTPC() && daughterTrack.itsNClsInnerBarrel() == 0 && - daughterTrack.itsNCls() < 4 && daughterTrack.tpcNClsCrossedRows() > 0.8 * daughterTrack.tpcNClsFindable() && daughterTrack.tpcNClsFound() > nTPCClusMinDaugKB; + daughterTrack.itsNCls() < 4 && daughterTrack.tpcNClsCrossedRows() > 0.8 * daughterTrack.tpcNClsFindable() && daughterTrack.tpcNClsFound() > nTPCClusMinDaugKB; if (motherGoodITS && daughterGoodITSTPC) { filterIndex += 1; fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); } else { - continue; + continue; } // 3 - mother track min pT @@ -556,7 +555,7 @@ struct sigmaminustask { filterIndex += 1; fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); } else { - continue; + continue; } // 4 - geometric cuts: eta @@ -564,7 +563,7 @@ struct sigmaminustask { filterIndex += 1; fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); } else { - continue; + continue; } // 5 - geometric cuts: phi difference @@ -572,9 +571,9 @@ struct sigmaminustask { filterIndex += 1; fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); } else { - continue; + continue; } - + // DCA calculation: initialization CCDB auto collision = motherTrack.template collision_as(); auto bc = collision.template bc_as(); @@ -611,7 +610,7 @@ struct sigmaminustask { filterIndex += 1; fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); } else { - continue; + continue; } // 8 - DCA daughter @@ -619,7 +618,7 @@ struct sigmaminustask { filterIndex += 1; fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); } else { - continue; + continue; } // 9 - radius cut @@ -627,7 +626,7 @@ struct sigmaminustask { filterIndex += 1; fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); } else { - continue; + continue; } // 10 - collision selection @@ -635,7 +634,7 @@ struct sigmaminustask { filterIndex += 1; fillFindableHistograms(filterIndex, mcRadius, recRadius, recPtMother, recPtDaughter, isSigmaMinus, isPiDaughter); } else { - continue; + continue; } // 11 - TOF daughter presence From 32cf5219fdb1f71b8e24faf194487c836357a0c7 Mon Sep 17 00:00:00 2001 From: morgmatt Date: Thu, 7 Aug 2025 16:39:56 +0200 Subject: [PATCH 16/18] Corrected error in histogram name --- PWGLF/TableProducer/Strangeness/sigmaminustask.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx index 869dc2967b1..198599c3a82 100644 --- a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx @@ -235,7 +235,7 @@ struct sigmaminustask { rSigmaMinus.fill(HIST("h2MassSigmaMinusPt"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.mSigmaMinus()); rSigmaMinus.fill(HIST("h2SigmaMassVsXiMass"), kinkCand.mXiMinus(), kinkCand.mSigmaMinus()); - rSigmaMinus.fill(HIST("h2NSigmaPiPt"), kinkCand.mothSign() * kinkCand.ptMoth(), dauTrack.tpcNSigmaPi()); + rSigmaMinus.fill(HIST("h2NSigmaTPCPiPt"), kinkCand.mothSign() * kinkCand.ptMoth(), dauTrack.tpcNSigmaPi()); rSigmaMinus.fill(HIST("h2DCAMothPt"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.dcaMothPv()); rSigmaMinus.fill(HIST("h2DCADaugPt"), kinkCand.mothSign() * kinkCand.ptDaug(), kinkCand.dcaDaugPv()); From 5200c021a70f9e778056bcdda39868654e165e31 Mon Sep 17 00:00:00 2001 From: morgmatt Date: Fri, 8 Aug 2025 12:04:17 +0200 Subject: [PATCH 17/18] Added Armenteros Podolanski histogram and pointing angle cosine histogram --- .../Strangeness/sigmaminustask.cxx | 31 +++++++++++++++++-- 1 file changed, 29 insertions(+), 2 deletions(-) diff --git a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx index dddfd65ea52..e16666a18d8 100644 --- a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx @@ -108,6 +108,9 @@ struct sigmaminustask { const AxisSpec dcaMothAxis{100, 0, 0.03, "DCA [cm]"}; const AxisSpec dcaDaugAxis{200, 0, 20, "DCA [cm]"}; const AxisSpec radiusAxis{100, -1, 40, "Decay radius [cm]"}; + const AxisSpec alphaAPAxis{200, -1.0, 1.0, "#alpha_{AP}"}; + const AxisSpec qtAPAxis{200, 0.0, 0.5, "q_{T,AP}"}; + const AxisSpec cosPointingAngleAxis{100, -1.0, 1.0, "Cos#theta_{PA}"}; const AxisSpec ptResolutionAxis{100, -1.0, 1.0, "(#it{p}_{T}^{rec} - #it{p}_{T}^{gen}) / #it{p}_{T}^{gen}"}; const AxisSpec massResolutionAxis{100, -0.5, 0.5, "(m_{rec} - m_{gen}) / m_{gen}"}; @@ -125,6 +128,8 @@ struct sigmaminustask { rSigmaMinus.add("h2NSigmaTPCPiPt", "h2NSigmaTPCPiPt", {HistType::kTH2F, {ptAxis, nSigmaPiAxis}}); rSigmaMinus.add("h2DCAMothPt", "h2DCAMothPt", {HistType::kTH2F, {ptAxis, dcaMothAxis}}); rSigmaMinus.add("h2DCADaugPt", "h2DCADaugPt", {HistType::kTH2F, {ptAxis, dcaDaugAxis}}); + rSigmaMinus.add("h2ArmenterosPreCuts", "h2ArmenterosPreCuts", {HistType::kTH2F, {alphaAPAxis, qtAPAxis}}); + rSigmaMinus.add("h2CosPointingAnglePt", "h2CosPointingAnglePt", {HistType::kTH2F, {ptAxis, cosPointingAngleAxis}}); if (doprocessMC) { // Add MC histograms if needed @@ -214,6 +219,14 @@ struct sigmaminustask { return std::sqrt(p2A - dp * dp / p2V0); } + float cosPAngle(const std::array& momMother, const std::array& posMother, const std::array& posKink) + { + std::array vMother = {posKink[0] - posMother[0], posKink[1] - posMother[1], posKink[2] - posMother[2]}; + float pMother = std::sqrt(std::inner_product(momMother.begin(), momMother.end(), momMother.begin(), 0.f)); + float vMotherNorm = std::sqrt(std::inner_product(vMother.begin(), vMother.end(), vMother.begin(), 0.f)); + return (std::inner_product(momMother.begin(), momMother.end(), vMother.begin(), 0.f)) / (pMother * vMotherNorm); + } + void processData(CollisionsFull::iterator const& collision, aod::KinkCands const& KinkCands, TracksFull const&) { if (std::abs(collision.posZ()) > cutzvertex || !collision.sel8()) { @@ -228,10 +241,17 @@ struct sigmaminustask { continue; } - float qt = qtAP(std::array{kinkCand.pxMoth(), kinkCand.pyMoth(), kinkCand.pzMoth()}, std::array{kinkCand.pxDaug(), kinkCand.pyDaug(), kinkCand.pzDaug()}); - if (qt < cutMinQtAP || qt > cutMaxQtAP) { + float alphaAPValue = alphaAP(std::array{kinkCand.pxMoth(), kinkCand.pyMoth(), kinkCand.pzMoth()}, std::array{kinkCand.pxDaug(), kinkCand.pyDaug(), kinkCand.pzDaug()}); + float qtValue = qtAP(std::array{kinkCand.pxMoth(), kinkCand.pyMoth(), kinkCand.pzMoth()}, std::array{kinkCand.pxDaug(), kinkCand.pyDaug(), kinkCand.pzDaug()}); + rSigmaMinus.fill(HIST("h2ArmenterosPreCuts"), alphaAPValue, qtValue); + + if (qtValue < cutMinQtAP || qtValue > cutMaxQtAP) { continue; } + float cosPointingAngleRec = cosPAngle(std::array{kinkCand.pxMoth(), kinkCand.pyMoth(), kinkCand.pzMoth()}, + std::array{0.0f, 0.0f, 0.0f}, + std::array{kinkCand.xDecVtx(), kinkCand.yDecVtx(), kinkCand.zDecVtx()}); + rSigmaMinus.fill(HIST("h2CosPointingAnglePt"), kinkCand.mothSign() * kinkCand.ptMoth(), cosPointingAngleRec); rSigmaMinus.fill(HIST("h2MassSigmaMinusPt"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.mSigmaMinus()); rSigmaMinus.fill(HIST("h2SigmaMassVsXiMass"), kinkCand.mXiMinus(), kinkCand.mSigmaMinus()); @@ -275,6 +295,9 @@ struct sigmaminustask { } // histograms filled with all kink candidates + float alphaAPValue = alphaAP(std::array{kinkCand.pxMoth(), kinkCand.pyMoth(), kinkCand.pzMoth()}, std::array{kinkCand.pxDaug(), kinkCand.pyDaug(), kinkCand.pzDaug()}); + float qtValue = qtAP(std::array{kinkCand.pxMoth(), kinkCand.pyMoth(), kinkCand.pzMoth()}, std::array{kinkCand.pxDaug(), kinkCand.pyDaug(), kinkCand.pzDaug()}); + rSigmaMinus.fill(HIST("h2ArmenterosPreCuts"), alphaAPValue, qtValue); rSigmaMinus.fill(HIST("h2MassSigmaMinusPt"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.mSigmaMinus()); rSigmaMinus.fill(HIST("h2SigmaMassVsXiMass"), kinkCand.mXiMinus(), kinkCand.mSigmaMinus()); rSigmaMinus.fill(HIST("h2NSigmaTPCPiPt"), kinkCand.mothSign() * kinkCand.ptMoth(), dauTrack.tpcNSigmaPi()); @@ -305,6 +328,9 @@ struct sigmaminustask { float deltaYMother = mcTrackPiDau.vy() - piMother.vy(); float decayRadiusMC = std::sqrt(deltaXMother * deltaXMother + deltaYMother * deltaYMother); float decayRadiusRec = std::sqrt(kinkCand.xDecVtx() * kinkCand.xDecVtx() + kinkCand.yDecVtx() * kinkCand.yDecVtx()); + float cosPointingAngleRec = cosPAngle(std::array{kinkCand.pxMoth(), kinkCand.pyMoth(), kinkCand.pzMoth()}, + std::array{0.0f, 0.0f, 0.0f}, + std::array{kinkCand.xDecVtx(), kinkCand.yDecVtx(), kinkCand.zDecVtx()}); // Check coherence of MCcollision Id for daughter MCparticle and reconstructed collision bool mcCollisionIdCheck = false; @@ -327,6 +353,7 @@ struct sigmaminustask { rSigmaMinus.fill(HIST("h2RadiusResolution"), kinkCand.mothSign() * kinkCand.ptMoth(), (decayRadiusRec - decayRadiusMC) / decayRadiusMC); rSigmaMinus.fill(HIST("h2DCAMothPt"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.dcaMothPv()); rSigmaMinus.fill(HIST("h2DCADaugPt"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.dcaDaugPv()); + rSigmaMinus.fill(HIST("h2CosPointingAnglePt"), kinkCand.mothSign() * kinkCand.ptMoth(), cosPointingAngleRec); if (std::abs(mcTrackPiDau.pdgCode()) == 211) { rSigmaMinus.fill(HIST("h2NSigmaTOFPiPt"), kinkCand.mothSign() * kinkCand.ptMoth(), dauTrack.tofNSigmaPi()); From e49c53dd68dbed15aea7b503669288a7c990d691 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Fri, 8 Aug 2025 10:05:04 +0000 Subject: [PATCH 18/18] Please consider the following formatting changes --- PWGLF/TableProducer/Strangeness/sigmaminustask.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx index e16666a18d8..c9b2d58f08c 100644 --- a/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx @@ -249,8 +249,8 @@ struct sigmaminustask { continue; } float cosPointingAngleRec = cosPAngle(std::array{kinkCand.pxMoth(), kinkCand.pyMoth(), kinkCand.pzMoth()}, - std::array{0.0f, 0.0f, 0.0f}, - std::array{kinkCand.xDecVtx(), kinkCand.yDecVtx(), kinkCand.zDecVtx()}); + std::array{0.0f, 0.0f, 0.0f}, + std::array{kinkCand.xDecVtx(), kinkCand.yDecVtx(), kinkCand.zDecVtx()}); rSigmaMinus.fill(HIST("h2CosPointingAnglePt"), kinkCand.mothSign() * kinkCand.ptMoth(), cosPointingAngleRec); rSigmaMinus.fill(HIST("h2MassSigmaMinusPt"), kinkCand.mothSign() * kinkCand.ptMoth(), kinkCand.mSigmaMinus());