diff --git a/PWGLF/DataModel/LFKinkDecayTables.h b/PWGLF/DataModel/LFKinkDecayTables.h index 9533bd66116..d8df2eb865f 100644 --- a/PWGLF/DataModel/LFKinkDecayTables.h +++ b/PWGLF/DataModel/LFKinkDecayTables.h @@ -47,6 +47,21 @@ DECLARE_SOA_COLUMN(DcaMothPv, dcaMothPv, float); //! DCA of the mother to th DECLARE_SOA_COLUMN(DcaDaugPv, dcaDaugPv, float); //! DCA of the daughter kink to the primary vertex DECLARE_SOA_COLUMN(DcaKinkTopo, dcaKinkTopo, float); //! DCA of the kink topology +DECLARE_SOA_COLUMN(NSigmaTPCPi, nSigmaTPCPi, float); //! Number of sigmas for the pion candidate from Sigma kink in TPC +DECLARE_SOA_COLUMN(NSigmaTPCPr, nSigmaTPCPr, float); //! Number of sigmas for the proton candidate from Sigma kink in TPC +DECLARE_SOA_COLUMN(NSigmaTPCKa, nSigmaTPCKa, float); //! Number of sigmas for the kaon candidate from Sigma kink in TPC +DECLARE_SOA_COLUMN(NSigmaTOFPi, nSigmaTOFPi, float); //! Number of sigmas for the pion candidate from Sigma kink in TOF +DECLARE_SOA_COLUMN(NSigmaTOFPr, nSigmaTOFPr, float); //! Number of sigmas for the proton candidate from Sigma kink in TOF +DECLARE_SOA_COLUMN(NSigmaTOFKa, nSigmaTOFKa, float); //! Number of sigmas for the kaon candidate from Sigma kink in TOF + +// MC Columns +DECLARE_SOA_COLUMN(MothPdgCode, mothPdgCode, int); //! PDG code of the Sigma daughter +DECLARE_SOA_COLUMN(DaugPdgCode, daugPdgCode, int); //! PDG code of the kink daughter +DECLARE_SOA_COLUMN(PtMC, ptMC, float); //! pT of the candidate in MC +DECLARE_SOA_COLUMN(MassMC, massMC, float); //! Invariant mass of the candidate in MC +DECLARE_SOA_COLUMN(DecayRadiusMC, decayRadiusMC, float); //! Decay radius of the candidate in MC +DECLARE_SOA_COLUMN(CollisionIdCheck, collisionIdCheck, bool); //! Check if mcDaughter collision ID matches the reconstructed collision ID + // DYNAMIC COLUMNS DECLARE_SOA_DYNAMIC_COLUMN(PxDaugNeut, pxDaugNeut, //! Px of the daughter neutral particle @@ -120,6 +135,26 @@ DECLARE_SOA_TABLE(KinkCandsUnbound, "AOD", "UBKINKCANDS", kinkcand::MSigmaPlus, kinkcand::MXiMinus); +DECLARE_SOA_TABLE(SlimKinkCands, "AOD", "SLIMKINKCANDS", + kinkcand::XDecVtx, kinkcand::YDecVtx, kinkcand::ZDecVtx, + kinkcand::PxMoth, kinkcand::PyMoth, kinkcand::PzMoth, + kinkcand::PxDaug, kinkcand::PyDaug, kinkcand::PzDaug, + kinkcand::DcaMothPv, kinkcand::DcaDaugPv, kinkcand::DcaKinkTopo, + kinkcand::MothSign, + kinkcand::NSigmaTPCPi, kinkcand::NSigmaTPCPr, kinkcand::NSigmaTPCKa, + kinkcand::NSigmaTOFPi, kinkcand::NSigmaTOFPr, kinkcand::NSigmaTOFKa); + +DECLARE_SOA_TABLE(SlimKinkCandsMC, "AOD", "SLIMKINKCANDSMC", + kinkcand::XDecVtx, kinkcand::YDecVtx, kinkcand::ZDecVtx, + kinkcand::PxMoth, kinkcand::PyMoth, kinkcand::PzMoth, + kinkcand::PxDaug, kinkcand::PyDaug, kinkcand::PzDaug, + kinkcand::DcaMothPv, kinkcand::DcaDaugPv, kinkcand::DcaKinkTopo, + kinkcand::MothSign, + kinkcand::NSigmaTPCPi, kinkcand::NSigmaTPCPr, kinkcand::NSigmaTPCKa, + kinkcand::NSigmaTOFPi, kinkcand::NSigmaTOFPr, kinkcand::NSigmaTOFKa, + kinkcand::MothPdgCode, kinkcand::DaugPdgCode, + kinkcand::PtMC, kinkcand::MassMC, kinkcand::DecayRadiusMC, kinkcand::CollisionIdCheck); + } // namespace o2::aod #endif // PWGLF_DATAMODEL_LFKINKDECAYTABLES_H_ diff --git a/PWGLF/TableProducer/Strangeness/CMakeLists.txt b/PWGLF/TableProducer/Strangeness/CMakeLists.txt index 55bf4550601..5b1080cab74 100644 --- a/PWGLF/TableProducer/Strangeness/CMakeLists.txt +++ b/PWGLF/TableProducer/Strangeness/CMakeLists.txt @@ -91,6 +91,11 @@ o2physics_add_dpl_workflow(lambdakzerospawner PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(sigmaminus-task + SOURCES sigmaminustask.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(strange-tree-creator SOURCES strangeTreeCreator.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore diff --git a/PWGLF/Tasks/Strangeness/sigmaminustask.cxx b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx similarity index 54% rename from PWGLF/Tasks/Strangeness/sigmaminustask.cxx rename to PWGLF/TableProducer/Strangeness/sigmaminustask.cxx index ac588fb8597..d60109009a5 100644 --- a/PWGLF/Tasks/Strangeness/sigmaminustask.cxx +++ b/PWGLF/TableProducer/Strangeness/sigmaminustask.cxx @@ -26,11 +26,18 @@ using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; -using TracksFull = soa::Join; +using TracksFull = soa::Join; using CollisionsFull = soa::Join; using CollisionsFullMC = soa::Join; struct sigmaminustask { + + // Output Tables + Produces outputDataTable; + Produces outputDataTableMC; + // Histograms are defined with HistogramRegistry HistogramRegistry rEventSelection{"eventSelection", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; HistogramRegistry rSigmaMinus{"sigmaminus", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; @@ -38,19 +45,25 @@ 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 fillOutputTree{"fillOutputTree", true, "If true, fill the output tree with Kink candidates"}; Preslice mPerCol = aod::track::collisionId; void init(InitContext const&) { // Axes - const AxisSpec ptAxis{50, -10, 10, "#it{p}_{T} (GeV/#it{c})"}; + const AxisSpec ptAxis{100, -10, 10, "#it{p}_{T} (GeV/#it{c})"}; const AxisSpec nSigmaPiAxis{100, -5, 5, "n#sigma_{#pi}"}; 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 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})"}; + // Event selection rEventSelection.add("hVertexZRec", "hVertexZRec", {HistType::kTH1F, {vertexZAxis}}); // Sigma-minus reconstruction @@ -62,6 +75,11 @@ struct sigmaminustask { // 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}}); } } @@ -71,14 +89,27 @@ struct sigmaminustask { return; } rEventSelection.fill(HIST("hVertexZRec"), collision.posZ()); + for (const auto& kinkCand : KinkCands) { auto dauTrack = kinkCand.trackDaug_as(); - if (abs(dauTrack.tpcNSigmaPi()) > cutNSigmaPi) { + + if (std::abs(dauTrack.tpcNSigmaPi()) > cutNSigmaPi) { 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()); + + if (fillOutputTree) { + outputDataTable(kinkCand.xDecVtx(), kinkCand.yDecVtx(), kinkCand.zDecVtx(), + kinkCand.pxMoth(), kinkCand.pyMoth(), kinkCand.pzMoth(), + kinkCand.pxDaug(), kinkCand.pyDaug(), kinkCand.pzDaug(), + kinkCand.dcaMothPv(), kinkCand.dcaDaugPv(), kinkCand.dcaKinkTopo(), + kinkCand.mothSign(), + dauTrack.tpcNSigmaPi(), dauTrack.tpcNSigmaPr(), dauTrack.tpcNSigmaKa(), + dauTrack.tofNSigmaPi(), dauTrack.tofNSigmaPr(), dauTrack.tofNSigmaKa()); + } } } PROCESS_SWITCH(sigmaminustask, processData, "Data processing", true); @@ -92,20 +123,23 @@ struct sigmaminustask { rEventSelection.fill(HIST("hVertexZRec"), collision.posZ()); auto kinkCandPerColl = KinkCands.sliceBy(mPerCol, collision.globalIndex()); + for (const auto& kinkCand : kinkCandPerColl) { + auto dauTrack = kinkCand.trackDaug_as(); auto mothTrack = kinkCand.trackMoth_as(); if (dauTrack.sign() != mothTrack.sign()) { LOG(info) << "Skipping kink candidate with opposite sign daughter and mother: " << kinkCand.globalIndex(); continue; // Skip if the daughter has the opposite sign as the mother } - if (abs(dauTrack.tpcNSigmaPi()) > cutNSigmaPi) { + if (std::abs(dauTrack.tpcNSigmaPi()) > cutNSigmaPi) { 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()); + // do MC association auto mcLabSigma = trackLabelsMC.rawIteratorAt(mothTrack.globalIndex()); auto mcLabPiDau = trackLabelsMC.rawIteratorAt(dauTrack.globalIndex()); @@ -119,36 +153,93 @@ struct sigmaminustask { if (piMother.globalIndex() != mcTrackSigma.globalIndex()) { continue; } - if (std::abs(mcTrackSigma.pdgCode()) != 3112 || std::abs(mcTrackPiDau.pdgCode()) != 211) { + if (std::abs(mcTrackSigma.pdgCode()) != 3112) { + continue; + } + if (std::abs(mcTrackPiDau.pdgCode()) != 211 && std::abs(mcTrackPiDau.pdgCode()) != 2212) { continue; } + + float MotherMassMC = std::sqrt(piMother.e() * piMother.e() - piMother.p() * piMother.p()); + float MotherpTMC = piMother.pt(); + float deltaXMother = mcTrackPiDau.vx() - piMother.vx(); + float deltaYMother = mcTrackPiDau.vy() - piMother.vy(); + float decayRadiusMC = std::sqrt(deltaXMother * deltaXMother + deltaYMother * deltaYMother); + + // Check coherence of MCcollision Id for daughter MCparticle and reconstructed collision + auto mcCollision = mcTrackPiDau.template mcCollision_as(); + bool mcCollisionIdCheck = collision.mcCollisionId() == mcCollision.globalIndex(); + 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); + } + + // fill the output table with Mc information + if (fillOutputTree) { + outputDataTableMC(kinkCand.xDecVtx(), kinkCand.yDecVtx(), kinkCand.zDecVtx(), + kinkCand.pxMoth(), kinkCand.pyMoth(), kinkCand.pzMoth(), + kinkCand.pxDaug(), kinkCand.pyDaug(), kinkCand.pzDaug(), + kinkCand.dcaMothPv(), kinkCand.dcaDaugPv(), kinkCand.dcaKinkTopo(), + kinkCand.mothSign(), + dauTrack.tpcNSigmaPi(), dauTrack.tpcNSigmaPr(), dauTrack.tpcNSigmaKa(), + dauTrack.tofNSigmaPi(), dauTrack.tofNSigmaPr(), dauTrack.tofNSigmaKa(), + mcTrackSigma.pdgCode(), mcTrackPiDau.pdgCode(), + MotherpTMC, MotherMassMC, decayRadiusMC, mcCollisionIdCheck); + } } - } - } - } + } // MC association and selection + } // kink cand loop + } // collision loop + + // Loop over all generated particles to fill MC histograms for (const auto& mcPart : particlesMC) { - if (std::abs(mcPart.pdgCode()) != 3112 || std::abs(mcPart.y()) > 0.5) { + if (std::abs(mcPart.pdgCode()) != 3112 || std::abs(mcPart.y()) > cutEtaMotherMC) { // only sigma mothers and rapidity cut continue; } if (!mcPart.has_daughters()) { continue; // Skip if no daughters } bool hasSigmaDaughter = false; + int daug_pdg = 0; + std::array secVtx; + std::array momDaug; for (const auto& daughter : mcPart.daughters_as()) { - if (std::abs(daughter.pdgCode()) == 211) { // Sigma PDG code + if (std::abs(daughter.pdgCode()) == 211 || std::abs(daughter.pdgCode()) == 2212) { // Pi or proton daughter hasSigmaDaughter = true; - break; // Found a pi daughter, exit loop + secVtx = {daughter.vx(), daughter.vy(), daughter.vz()}; + momDaug = {daughter.px(), daughter.py(), daughter.pz()}; + daug_pdg = daughter.pdgCode(); + break; // Found a daughter, exit loop } } if (!hasSigmaDaughter) { - continue; // Skip if no pi daughter found + continue; // Skip if no pi/proton daughter found } float mcMass = std::sqrt(mcPart.e() * mcPart.e() - mcPart.p() * mcPart.p()); + float mcDecayRadius = std::sqrt((secVtx[0] - mcPart.vx()) * (secVtx[0] - mcPart.vx()) + (secVtx[1] - mcPart.vy()) * (secVtx[1] - mcPart.vy())); int sigmaSign = mcPart.pdgCode() > 0 ? 1 : -1; // Determine the sign of the Sigma rSigmaMinus.fill(HIST("h2MassPtMCGen"), sigmaSign * mcPart.pt(), mcMass); + + // Fill output table with non reconstructed MC candidates + if (fillOutputTree) { + outputDataTableMC(-999, -999, -999, + -999, -999, -999, + -999, -999, -999, + -999, -999, -999, + sigmaSign, + -999, -999, -999, + -999, -999, -999, + mcPart.pdgCode(), daug_pdg, + mcPart.pt(), mcMass, mcDecayRadius, false); + } } } + PROCESS_SWITCH(sigmaminustask, processMC, "MC processing", false); }; diff --git a/PWGLF/Tasks/Strangeness/CMakeLists.txt b/PWGLF/Tasks/Strangeness/CMakeLists.txt index a6bc1f60334..58e53693fb8 100644 --- a/PWGLF/Tasks/Strangeness/CMakeLists.txt +++ b/PWGLF/Tasks/Strangeness/CMakeLists.txt @@ -29,11 +29,6 @@ o2physics_add_dpl_workflow(cascadeanalysis PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) -o2physics_add_dpl_workflow(sigmaminus-task - SOURCES sigmaminustask.cxx - PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore - COMPONENT_NAME Analysis) - o2physics_add_dpl_workflow(cascadeanalysismc SOURCES cascadeanalysisMC.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore