From 4255ce2ac6366bdf6c4161196f0a75eb11057703 Mon Sep 17 00:00:00 2001 From: marcellocosti Date: Mon, 1 Jun 2026 13:13:01 +0200 Subject: [PATCH 1/5] MC matching for kinkBuilder + SP process function for L1405 --- PWGLF/DataModel/LFLambda1405Table.h | 16 + PWGLF/TableProducer/Common/kinkBuilder.cxx | 676 ++++++++++---- PWGLF/Tasks/Resonances/lambda1405analysis.cxx | 838 ++++++++++++++---- PWGLF/Utils/svPoolCreator.h | 35 +- 4 files changed, 1217 insertions(+), 348 deletions(-) diff --git a/PWGLF/DataModel/LFLambda1405Table.h b/PWGLF/DataModel/LFLambda1405Table.h index be61dc083bc..cb472f5fc58 100644 --- a/PWGLF/DataModel/LFLambda1405Table.h +++ b/PWGLF/DataModel/LFLambda1405Table.h @@ -29,6 +29,7 @@ namespace lambda1405 DECLARE_SOA_COLUMN(Px, px, float); //! Px of the candidate DECLARE_SOA_COLUMN(Py, py, float); //! Py of the candidate DECLARE_SOA_COLUMN(Pz, pz, float); //! Pz of the candidate +DECLARE_SOA_COLUMN(Pt, pt, float); //! Pt of the candidate DECLARE_SOA_COLUMN(Mass, mass, float); //! Invariant mass of the candidate DECLARE_SOA_COLUMN(MassXi1530, massXi1530, float); //! Invariant mass of the Xi(1530) candidate DECLARE_SOA_COLUMN(SigmaMinusMass, sigmaMinusMass, float); //! Invariant mass of the Sigma- candidate @@ -47,6 +48,10 @@ DECLARE_SOA_COLUMN(DCAKinkDauToPV, dcaKinkDauToPV, float); //! DCA of the kink DECLARE_SOA_COLUMN(NSigmaTPCPiDau, nSigmaTPCPiDau, float); //! Number of sigmas for the lambda1405 pion daughter in TPC DECLARE_SOA_COLUMN(NSigmaTOFPiDau, nSigmaTOFPiDau, float); //! Number of sigmas for the lambda1405 pion daughter in TOF +// Flow columns +DECLARE_SOA_COLUMN(ScalarProd, scalarProd, float); //! Scalar product of the candidate +DECLARE_SOA_COLUMN(Centrality, centrality, float); //! Centrality of the candidate + // MC Columns 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 @@ -67,6 +72,17 @@ DECLARE_SOA_TABLE(Lambda1405Cands, "AOD", "LAMBDA1405", lambda1405::DCAKinkDauToPV, lambda1405::NSigmaTPCPiDau, lambda1405::NSigmaTOFPiDau); +DECLARE_SOA_TABLE(Lambda1405Flow, "AOD", "LAMBDA1405FLOW", + o2::soa::Index<>, + lambda1405::Pt, + lambda1405::Mass, lambda1405::SigmaMinusMass, lambda1405::SigmaPlusMass, + lambda1405::AlphaAPSigma, lambda1405::QtAPSigma, + lambda1405::NSigmaTPCPiKink, lambda1405::NSigmaTOFPiKink, + lambda1405::NSigmaTPCPrKink, lambda1405::NSigmaTOFPrKink, + lambda1405::DCAKinkDauToPV, + lambda1405::NSigmaTPCPiDau, lambda1405::NSigmaTOFPiDau, + lambda1405::ScalarProd, lambda1405::Centrality); + DECLARE_SOA_TABLE(Lambda1405CandsMC, "AOD", "MCLAMBDA1405", o2::soa::Index<>, lambda1405::Px, lambda1405::Py, lambda1405::Pz, diff --git a/PWGLF/TableProducer/Common/kinkBuilder.cxx b/PWGLF/TableProducer/Common/kinkBuilder.cxx index cd2f4c89cf1..330e0f6315b 100644 --- a/PWGLF/TableProducer/Common/kinkBuilder.cxx +++ b/PWGLF/TableProducer/Common/kinkBuilder.cxx @@ -16,6 +16,7 @@ #include "PWGLF/DataModel/LFKinkDecayTables.h" #include "PWGLF/Utils/svPoolCreator.h" +#include "Common/Core/RecoDecay.h" #include "Common/Core/trackUtilities.h" #include @@ -40,6 +41,7 @@ #include #include +#include #include #include @@ -50,11 +52,20 @@ #include using namespace o2; +using namespace o2::aod; using namespace o2::framework; using namespace o2::framework::expressions; +using namespace o2::constants::physics; using std::array; using VBracket = o2::math_utils::Bracket; using TracksFull = soa::Join; +using TracksFullMc = soa::Join; +using McRecoCollisionsCent = soa::Join; + +enum KinkDecayType { kSigmaMinusToPiMinusNeutron = 0, + kSigmaPlusToPiPlusNeutron, + kSigmaPlusToProtonPi0, + kNMatchedDecays }; namespace { @@ -62,12 +73,29 @@ constexpr std::array LayerRadii{2.33959f, 3.14076f, 3.91924f, 19.6213f constexpr double betheBlochDefault[1][6]{{-1.e32, -1.e32, -1.e32, -1.e32, -1.e32, -1.e32}}; static const std::vector betheBlochParNames{"p0", "p1", "p2", "p3", "p4", "resolution"}; static const std::vector particleNames{"Daughter"}; +constexpr int DepthMcMatchMax = 1; // Max depth for MC matching std::shared_ptr h2ClsMapPtMoth; std::shared_ptr h2ClsMapPtDaug; std::shared_ptr h2DeDxDaugSel; std::shared_ptr h2KinkAnglePt; std::shared_ptr h2MothMassPt; +std::shared_ptr hSelMotherQA; +std::shared_ptr hSelDaugQA; +std::shared_ptr hSelKinkedTrackQA; +std::shared_ptr hMothDaughSignsInit; +std::shared_ptr hMothDaughSignsFinal; +std::shared_ptr h2ItsClsMothBeforeSel; +std::shared_ptr h2ItsClsDaugBeforeSel; +std::shared_ptr hZDiff; +std::shared_ptr hPhiDiff; +std::shared_ptr hDCAMothToPV; +std::shared_ptr hDCADaugToPV; +std::shared_ptr hMothDecRad2; +std::shared_ptr hGenCandidates; +std::shared_ptr hRecCandidates; +std::array, kNMatchedDecays> hGenPtKinkAngle; +std::array, kNMatchedDecays> hRecPtKinkAngle; } // namespace struct kinkCandidate { @@ -137,6 +165,7 @@ struct kinkBuilder { Configurable customVertexerTimeMargin{"customVertexerTimeMargin", 800, "Time margin for custom vertexer (ns)"}; Configurable skipAmbiTracks{"skipAmbiTracks", false, "Skip ambiguous tracks"}; Configurable unlikeSignBkg{"unlikeSignBkg", false, "Use unlike sign background"}; + Configurable skipBkgCands{"skipBkgCands", true, "Skip bkg candidates in MC process"}; // CCDB options Configurable ccdbPath{"ccdbPath", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; @@ -149,6 +178,11 @@ struct kinkBuilder { ConfigurableAxis rigidityBins{"rigidityBins", {200, -10.f, 10.f}, "Binning for rigidity #it{p}^{TPC}/#it{z}"}; ConfigurableAxis dedxBins{"dedxBins", {1000, 0.f, 1000.f}, "Binning for dE/dx"}; ConfigurableAxis nSigmaBins{"nSigmaBins", {200, -5.f, 5.f}, "Binning for n sigma"}; + ConfigurableAxis zDiffBins{"zDiffBins", {300, 0.f, 30.f}, "Binning for z difference"}; + ConfigurableAxis phiDiffBins{"phiDiffBins", {360, 0.f, 360.f}, "Binning for phi difference"}; + ConfigurableAxis dcaMothToPVBins{"dcaMothToPVBins", {120, 0.f, 60.f}, "Binning for DCA of mother to PV"}; + ConfigurableAxis dcaDaugToPVBins{"dcaDaugToPVBins", {200, 0.f, 20.f}, "Binning for DCA of daughter to PV"}; + ConfigurableAxis mothDecRad2Bins{"mothDecRad2Bins", {100, 0, 100}, "Binning for mother decay radius squared"}; // std vector of candidates std::vector kinkCandidates; @@ -161,6 +195,7 @@ struct kinkBuilder { // mother and daughter tracks' properties (absolute charge and mass) int charge = 1; + std::vector mothMatchPdgCodes{}; float mothMass = o2::constants::physics::MassSigmaMinus; float chargedDauMass = o2::constants::physics::MassPiMinus; float neutDauMass = o2::constants::physics::MassNeutron; @@ -169,6 +204,7 @@ struct kinkBuilder { { if (hypoMoth == kSigmaMinus) { charge = 1; + mothMatchPdgCodes = {PDG_t::kSigmaMinus, PDG_t::kSigmaPlus}; mothMass = o2::constants::physics::MassSigmaMinus; chargedDauMass = o2::constants::physics::MassPiMinus; neutDauMass = o2::constants::physics::MassNeutron; @@ -209,6 +245,7 @@ struct kinkBuilder { const AxisSpec itsClusterMapAxis(128, 0, 127, "ITS cluster map"); const AxisSpec rigidityAxis{rigidityBins, "#it{p}^{TPC}/#it{z}"}; const AxisSpec ptAxis{rigidityBins, "#it{p}_{T} (GeV/#it{c})"}; + const AxisSpec absPtAxis{1000, 0, 10, "#it{p}_{T} (GeV/#it{c})"}; const AxisSpec kinkAngleAxis{100, 0, 180, "#theta_{kink} (deg)"}; const AxisSpec dedxAxis{dedxBins, "d#it{E}/d#it{x}"}; @@ -226,208 +263,350 @@ struct kinkBuilder { h2MothMassPt = qaRegistry.add("h2MothMassPt", "; p_{T} (GeV/#it{c}); m (GeV/#it{c}^{2})", HistType::kTH2F, {ptAxis, massAxis}); h2ClsMapPtMoth = qaRegistry.add("h2ClsMapPtMoth", "; p_{T} (GeV/#it{c}); ITS cluster map", HistType::kTH2F, {ptAxis, itsClusterMapAxis}); h2ClsMapPtDaug = qaRegistry.add("h2ClsMapPtDaug", "; p_{T} (GeV/#it{c}); ITS cluster map", HistType::kTH2F, {ptAxis, itsClusterMapAxis}); + h2ItsClsMothBeforeSel = qaRegistry.add("h2ItsClsMothBeforeSel", "; Tot Cls; IB clusters", HistType::kTH2F, {{8, -0.5, 7.5}, {4, -0.5, 3.5}}); + h2ItsClsDaugBeforeSel = qaRegistry.add("h2ItsClsDaugBeforeSel", "; Tot Cls; IB clusters", HistType::kTH2F, {{8, -0.5, 7.5}, {4, -0.5, 3.5}}); + + // QA + hSelMotherQA = qaRegistry.add("hSelMotherQA", ";;Q_{Mother}", {HistType::kTH2F, {{9, -0.5, 8.5}, {2, -0.5, 1.5}}}); + hSelMotherQA->GetYaxis()->SetBinLabel(1, "-"); + hSelMotherQA->GetYaxis()->SetBinLabel(2, "+"); + hSelMotherQA->GetXaxis()->SetBinLabel(1, "All"); + hSelMotherQA->GetXaxis()->SetBinLabel(2, "Has Collision"); + hSelMotherQA->GetXaxis()->SetBinLabel(3, "Has ITS"); + hSelMotherQA->GetXaxis()->SetBinLabel(4, "Has NO TPC"); + hSelMotherQA->GetXaxis()->SetBinLabel(5, "Has NO TOF"); + hSelMotherQA->GetXaxis()->SetBinLabel(6, "ITS cls sel"); + hSelMotherQA->GetXaxis()->SetBinLabel(7, "ITS IB cls sel"); + hSelMotherQA->GetXaxis()->SetBinLabel(8, "ITS chi2 N cls sel"); + hSelMotherQA->GetXaxis()->SetBinLabel(9, "Pt sel"); + + // QA + hSelDaugQA = qaRegistry.add("hSelDaugQA", ";;Q_{Daug}", {HistType::kTH2F, {{8, -0.5, 7.5}, {2, -0.5, 1.5}}}); + hSelDaugQA->GetYaxis()->SetBinLabel(1, "-"); + hSelDaugQA->GetYaxis()->SetBinLabel(2, "+"); + hSelDaugQA->GetXaxis()->SetBinLabel(1, "All"); + hSelDaugQA->GetXaxis()->SetBinLabel(2, "Has ITS"); + hSelDaugQA->GetXaxis()->SetBinLabel(3, "Has TPC"); + hSelDaugQA->GetXaxis()->SetBinLabel(4, "Has TOF"); + hSelDaugQA->GetXaxis()->SetBinLabel(5, "ITS tot cls sel"); + hSelDaugQA->GetXaxis()->SetBinLabel(6, "ITS IB cls sel"); + hSelDaugQA->GetXaxis()->SetBinLabel(7, "TPC Crossed rows vs findable sel"); + hSelDaugQA->GetXaxis()->SetBinLabel(8, "TPC cls sel"); + + hSelKinkedTrackQA = qaRegistry.add("hSelKinkedTrackQA", ";;(Q_{Mother}, Q_{Daughter})", {HistType::kTH2F, {{11, -0.5, 10.5}, {4, -0.5, 3.5}}}); + hSelKinkedTrackQA->GetYaxis()->SetBinLabel(1, "(+,+)"); + hSelKinkedTrackQA->GetYaxis()->SetBinLabel(2, "(-,-)"); + hSelKinkedTrackQA->GetYaxis()->SetBinLabel(3, "(+,-)"); + hSelKinkedTrackQA->GetYaxis()->SetBinLabel(4, "(-,+)"); + hSelKinkedTrackQA->GetXaxis()->SetBinLabel(1, "All"); + hSelKinkedTrackQA->GetXaxis()->SetBinLabel(2, "MaxDCAMothToPV"); + hSelKinkedTrackQA->GetXaxis()->SetBinLabel(3, "MaxZDiff"); + hSelKinkedTrackQA->GetXaxis()->SetBinLabel(4, "MaxPhiDiff"); + hSelKinkedTrackQA->GetXaxis()->SetBinLabel(5, "MinDCADaugToPV"); + hSelKinkedTrackQA->GetXaxis()->SetBinLabel(6, "UpdateMothTrackUsePV"); + hSelKinkedTrackQA->GetXaxis()->SetBinLabel(7, "DCAFitter"); + hSelKinkedTrackQA->GetXaxis()->SetBinLabel(8, "PropagateTracksToVertex"); + hSelKinkedTrackQA->GetXaxis()->SetBinLabel(9, "MothDecRad2InsideL4"); + hSelKinkedTrackQA->GetXaxis()->SetBinLabel(10, "MothDaugLastLayers"); + hSelKinkedTrackQA->GetXaxis()->SetBinLabel(11, "MothLastLayerRadiusCheck"); + hMothDaughSignsInit = qaRegistry.add("hMothDaughSignsInit", "; Sign Mother; Sign Daughter", {HistType::kTH2F, {{3, -1.5, 1.5}, {3, -1.5, 1.5}}}); + hMothDaughSignsFinal = qaRegistry.add("hMothDaughSignsFinal", "; Sign Mother; Sign Daughter", {HistType::kTH2F, {{3, -1.5, 1.5}, {3, -1.5, 1.5}}}); + hZDiff = qaRegistry.add("hZDiff", "; #Delta z (#mu m);(Q_{Mother}, Q_{Daughter})", HistType::kTH2F, {zDiffBins, {4, -0.5, 3.5}}); + hZDiff->GetYaxis()->SetBinLabel(1, "(+,+)"); + hZDiff->GetYaxis()->SetBinLabel(2, "(-,-)"); + hZDiff->GetYaxis()->SetBinLabel(3, "(+,-)"); + hZDiff->GetYaxis()->SetBinLabel(4, "(-,+)"); + hPhiDiff = qaRegistry.add("hPhiDiff", "; #Delta #phi (rad);(Q_{Mother}, Q_{Daughter})", HistType::kTH2F, {phiDiffBins, {4, -0.5, 3.5}}); + hPhiDiff->GetYaxis()->SetBinLabel(1, "(+,+)"); + hPhiDiff->GetYaxis()->SetBinLabel(2, "(-,-)"); + hPhiDiff->GetYaxis()->SetBinLabel(3, "(+,-)"); + hPhiDiff->GetYaxis()->SetBinLabel(4, "(-,+)"); + hDCAMothToPV = qaRegistry.add("hDCAMothToPV", "; DCA moth to PV;(Q_{Mother}, Q_{Daughter})", HistType::kTH2F, {dcaMothToPVBins, {4, -0.5, 3.5}}); + hDCAMothToPV->GetYaxis()->SetBinLabel(1, "(+,+)"); + hDCAMothToPV->GetYaxis()->SetBinLabel(2, "(-,-)"); + hDCAMothToPV->GetYaxis()->SetBinLabel(3, "(+,-)"); + hDCAMothToPV->GetYaxis()->SetBinLabel(4, "(-,+)"); + hDCADaugToPV = qaRegistry.add("hDCADaugToPV", "; DCA daug to PV;(Q_{Mother}, Q_{Daughter})", HistType::kTH2F, {dcaDaugToPVBins, {4, -0.5, 3.5}}); + hDCADaugToPV->GetYaxis()->SetBinLabel(1, "(+,+)"); + hDCADaugToPV->GetYaxis()->SetBinLabel(2, "(-,-)"); + hDCADaugToPV->GetYaxis()->SetBinLabel(3, "(+,-)"); + hDCADaugToPV->GetYaxis()->SetBinLabel(4, "(-,+)"); + hMothDecRad2 = qaRegistry.add("hMothDecRad2", "; Mother Dec. Radius Squared (cm^{2});(Q_{Mother}, Q_{Daughter})", HistType::kTH2F, {mothDecRad2Bins, {4, -0.5, 3.5}}); + hMothDecRad2->GetYaxis()->SetBinLabel(1, "(+,+)"); + hMothDecRad2->GetYaxis()->SetBinLabel(2, "(-,-)"); + hMothDecRad2->GetYaxis()->SetBinLabel(3, "(+,-)"); + hMothDecRad2->GetYaxis()->SetBinLabel(4, "(-,+)"); for (int i = 0; i < 5; i++) { mBBparamsDaug[i] = cfgBetheBlochParams->get("Daughter", Form("p%i", i)); } mBBparamsDaug[5] = cfgBetheBlochParams->get("Daughter", "resolution"); - } - template - bool selectMothTrack(const T& candidate) - { - if (candidate.has_collision() && candidate.hasITS() && !candidate.hasTPC() && !candidate.hasTOF() && candidate.itsNCls() < 6 && - candidate.itsNClsInnerBarrel() == 3 && candidate.itsChi2NCl() < 36 && candidate.pt() > minPtMoth) { - return true; + if (doprocessMc) { + if (skipBkgCands) { + hRecCandidates = qaRegistry.add("hRecCandidates", ";Counts;", {HistType::kTH2F, {{kNMatchedDecays, -0.5, kNMatchedDecays-0.5}, absPtAxis}}); + hRecCandidates->GetXaxis()->SetBinLabel(1, "#Sigma^{-} #rightarrow n#pi^{-}"); + hRecCandidates->GetXaxis()->SetBinLabel(2, "#Sigma^{+} #rightarrow n#pi^{+}"); + hRecCandidates->GetXaxis()->SetBinLabel(3, "#Sigma^{+} #rightarrow p#pi^{0}"); + } + hGenCandidates = qaRegistry.add("hGenCandidates", ";Counts;", {HistType::kTH2F, {{kNMatchedDecays, -0.5, kNMatchedDecays-0.5}, absPtAxis}}); + hGenCandidates->GetXaxis()->SetBinLabel(1, "#Sigma^{-} #rightarrow n#pi^{-}"); + hGenCandidates->GetXaxis()->SetBinLabel(2, "#Sigma^{+} #rightarrow n#pi^{+}"); + hGenCandidates->GetXaxis()->SetBinLabel(3, "#Sigma^{+} #rightarrow p#pi^{0}"); + hRecPtKinkAngle[kSigmaMinusToPiMinusNeutron] = qaRegistry.add("hRecPtKinkAngleSigmaMinusToPiMinusNeutron", ";Counts;", {HistType::kTH2F, {absPtAxis, kinkAngleAxis}}); + hRecPtKinkAngle[kSigmaPlusToPiPlusNeutron] = qaRegistry.add("hRecPtKinkAngleSigmaPlusToPiPlusNeutron", ";Counts;", {HistType::kTH2F, {absPtAxis, kinkAngleAxis}}); + hRecPtKinkAngle[kSigmaPlusToProtonPi0] = qaRegistry.add("hRecPtKinkAngleSigmaPlusToProtonPi0", ";Counts;", {HistType::kTH2F, {absPtAxis, kinkAngleAxis}}); + hGenPtKinkAngle[kSigmaMinusToPiMinusNeutron] = qaRegistry.add("hGenPtKinkAngleSigmaMinusToPiMinusNeutron", ";Counts;", {HistType::kTH2F, {absPtAxis, kinkAngleAxis}}); + hGenPtKinkAngle[kSigmaPlusToPiPlusNeutron] = qaRegistry.add("hGenPtKinkAngleSigmaPlusToPiPlusNeutron", ";Counts;", {HistType::kTH2F, {absPtAxis, kinkAngleAxis}}); + hGenPtKinkAngle[kSigmaPlusToProtonPi0] = qaRegistry.add("hGenPtKinkAngleSigmaPlusToProtonPi0", ";Counts;", {HistType::kTH2F, {absPtAxis, kinkAngleAxis}}); } - return false; } template - bool selectDaugTrack(const T& candidate) + bool selectMothTrack(const T& candidate) { - if (!candidate.hasTPC() || !candidate.hasITS()) { + bool isPositive = candidate.sign() > 0; + hSelMotherQA->Fill(0.f, isPositive); + + if (!candidate.has_collision()) return false; - } + hSelMotherQA->Fill(1.f, isPositive); - if (askTOFforDaug && !candidate.hasTOF()) { + if (!candidate.hasITS()) return false; - } + hSelMotherQA->Fill(2.f, isPositive); - bool isGoodTPCCand = false; - if (candidate.itsNClsInnerBarrel() == 0 && candidate.itsNCls() < 4 && - candidate.tpcNClsCrossedRows() > 0.8 * candidate.tpcNClsFindable() && candidate.tpcNClsFound() > nTPCClusMinDaug) { - isGoodTPCCand = true; - } + if (candidate.hasTPC()) + return false; + hSelMotherQA->Fill(3.f, isPositive); - if (!isGoodTPCCand) { + if (candidate.hasTOF()) return false; - } + hSelMotherQA->Fill(4.f, isPositive); + + h2ItsClsMothBeforeSel->Fill(candidate.itsNCls(), candidate.itsNClsInnerBarrel()); + if (candidate.itsNCls() >= 6) + return false; + hSelMotherQA->Fill(5.f, isPositive); + + if (candidate.itsNClsInnerBarrel() != 3) + return false; + hSelMotherQA->Fill(6.f, isPositive); + + if (candidate.itsChi2NCl() >= 36) + return false; + hSelMotherQA->Fill(7.f, isPositive); + + if (candidate.pt() <= minPtMoth) + return false; + hSelMotherQA->Fill(8.f, isPositive); return true; } - template - void fillCandidateData(const Tcolls& collisions, const Ttracks& tracks, aod::AmbiguousTracks const& ambiguousTracks, aod::BCs const& bcs) + template + bool selectDaugTrack(const T& candidate) { - svCreator.clearPools(); - svCreator.fillBC2Coll(collisions, bcs); + bool isPositive = candidate.sign() > 0; + hSelDaugQA->Fill(0.f, isPositive); - for (const auto& track : tracks) { - if (std::abs(track.eta()) > etaMax) - continue; + if (!candidate.hasITS()) + return false; + hSelDaugQA->Fill(1.f, isPositive); - bool isDaug = selectDaugTrack(track); - bool isMoth = selectMothTrack(track); + if (!candidate.hasTPC()) + return false; + hSelDaugQA->Fill(2.f, isPositive); - if (!isDaug && !isMoth) - continue; + if (askTOFforDaug && !candidate.hasTOF()) + return false; + hSelDaugQA->Fill(3.f, isPositive); - int pdgHypo = isMoth ? 1 : 0; - svCreator.appendTrackCand(track, collisions, pdgHypo, ambiguousTracks, bcs); - } - auto& kinkPool = svCreator.getSVCandPool(collisions, !unlikeSignBkg); + h2ItsClsDaugBeforeSel->Fill(candidate.itsNCls(), candidate.itsNClsInnerBarrel()); + if (candidate.itsNClsInnerBarrel() != 0) + return false; + hSelDaugQA->Fill(4.f, isPositive); + + if (candidate.itsNCls() >= 4) + return false; + hSelDaugQA->Fill(5.f, isPositive); + + if (candidate.tpcNClsCrossedRows() <= 0.8 * candidate.tpcNClsFindable()) + return false; + hSelDaugQA->Fill(6.f, isPositive); + + if (candidate.tpcNClsFound() <= nTPCClusMinDaug) + return false; + hSelDaugQA->Fill(7.f, isPositive); - for (const auto& svCand : kinkPool) { - kinkCandidate kinkCand; + return true; + } - auto trackMoth = tracks.rawIteratorAt(svCand.tr0Idx); - auto trackDaug = tracks.rawIteratorAt(svCand.tr1Idx); + template + void buildKinkCand(const TSVCand& svCand, const TTracks& tracks, bool& isAccepted, const TColls& collisions) + { + isAccepted = false; + kinkCandidate kinkCand; - auto const& collision = trackMoth.template collision_as(); - auto const& bc = collision.template bc_as(); - initCCDB(bc); + auto trackMoth = tracks.rawIteratorAt(svCand.tr0Idx); + auto trackDaug = tracks.rawIteratorAt(svCand.tr1Idx); - o2::dataformats::VertexBase primaryVertex; - primaryVertex.setPos({collision.posX(), collision.posY(), collision.posZ()}); - primaryVertex.setCov(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ()); - kinkCand.primVtx = {primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}; + // Fill Selections QA histo + int chargeCombSvCand = 2*unlikeSignBkg + (trackMoth.sign() == -1 ? 1 : 0); + hSelKinkedTrackQA->Fill(0.f, chargeCombSvCand); // all candidates bin + hMothDaughSignsInit->Fill(trackMoth.sign(), trackDaug.sign()); - o2::track::TrackParCov trackParCovMoth = getTrackParCov(trackMoth); - o2::track::TrackParCov trackParCovMothPV{trackParCovMoth}; - o2::base::Propagator::Instance()->PropagateToXBxByBz(trackParCovMoth, LayerRadii[trackMoth.itsNCls() - 1]); + auto const& collision = trackMoth.template collision_as(); + auto const& bc = collision.template bc_as(); + initCCDB(bc); - std::array dcaInfoMoth; - o2::base::Propagator::Instance()->propagateToDCABxByBz({primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, trackParCovMothPV, 2.f, static_cast(cfgMaterialCorrection.value), &dcaInfoMoth); + o2::dataformats::VertexBase primaryVertex; + primaryVertex.setPos({collision.posX(), collision.posY(), collision.posZ()}); + primaryVertex.setCov(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ()); + kinkCand.primVtx = {primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}; - if (std::abs(dcaInfoMoth[0]) > maxDCAMothToPV) { - continue; - } + o2::track::TrackParCov trackParCovMoth = getTrackParCov(trackMoth); + o2::track::TrackParCov trackParCovMothPV{trackParCovMoth}; + o2::base::Propagator::Instance()->PropagateToXBxByBz(trackParCovMoth, LayerRadii[trackMoth.itsNCls() - 1]); - o2::track::TrackParCov trackParCovDaug = getTrackParCov(trackDaug); + std::array dcaInfoMoth; + o2::base::Propagator::Instance()->propagateToDCABxByBz({primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, trackParCovMothPV, 2.f, static_cast(cfgMaterialCorrection.value), &dcaInfoMoth); - // check if the kink daughter is close to the mother - if (std::abs(trackParCovMoth.getZ() - trackParCovDaug.getZ()) > maxZDiff) { - continue; - } + hDCAMothToPV->Fill(std::abs(dcaInfoMoth[0]), chargeCombSvCand); + if (std::abs(dcaInfoMoth[0]) > maxDCAMothToPV) { + return; + } + hSelKinkedTrackQA->Fill(1.f, chargeCombSvCand); // MaxDCAMothToPV cut bin - if ((std::abs(trackParCovMoth.getPhi() - trackParCovDaug.getPhi()) * radToDeg) > maxPhiDiff) { - continue; - } + o2::track::TrackParCov trackParCovDaug = getTrackParCov(trackDaug); - // propagate to PV - std::array dcaInfoDaug; - o2::base::Propagator::Instance()->propagateToDCABxByBz({primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, trackParCovDaug, 2.f, static_cast(cfgMaterialCorrection.value), &dcaInfoDaug); - if (std::abs(dcaInfoDaug[0]) < minDCADaugToPV) { - continue; - } + // check if the kink daughter is close to the mother + hZDiff->Fill(std::abs(trackParCovMoth.getZ() - trackParCovDaug.getZ()), chargeCombSvCand); + if (std::abs(trackParCovMoth.getZ() - trackParCovDaug.getZ()) > maxZDiff) { + return; + } + hSelKinkedTrackQA->Fill(2.f, chargeCombSvCand); // MaxZDiff cut bin - if (updateMothTrackUsePV) { - // update the mother track parameters using the primary vertex - trackParCovMoth = trackParCovMothPV; - if (!trackParCovMoth.update(primaryVertex)) { - continue; - } - } + hPhiDiff->Fill(std::abs(trackParCovMoth.getPhi() - trackParCovDaug.getPhi()) * radToDeg, chargeCombSvCand); + if ((std::abs(trackParCovMoth.getPhi() - trackParCovDaug.getPhi()) * radToDeg) > maxPhiDiff) { + return; + } + hSelKinkedTrackQA->Fill(3.f, chargeCombSvCand); // MaxPhiDiff cut bin - int nCand = 0; - try { - nCand = fitter.process(trackParCovMoth, trackParCovDaug); - } catch (...) { - LOG(error) << "Exception caught in DCA fitter process call!"; - continue; - } - if (nCand == 0) { - continue; - } + // propagate to PV + std::array dcaInfoDaug; + o2::base::Propagator::Instance()->propagateToDCABxByBz({primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, trackParCovDaug, 2.f, static_cast(cfgMaterialCorrection.value), &dcaInfoDaug); + hDCADaugToPV->Fill(std::abs(dcaInfoDaug[0]), chargeCombSvCand); + if (std::abs(dcaInfoDaug[0]) < minDCADaugToPV) { + return; + } + hSelKinkedTrackQA->Fill(4.f, chargeCombSvCand); // MinDCADaugToPV cut bin - if (!fitter.propagateTracksToVertex()) { - continue; + if (updateMothTrackUsePV) { + // update the mother track parameters using the primary vertex + trackParCovMoth = trackParCovMothPV; + if (!trackParCovMoth.update(primaryVertex)) { + return; } + } + hSelKinkedTrackQA->Fill(5.f, chargeCombSvCand); // UpdateMothTrackUsePV cut bin - auto propMothTrack = fitter.getTrack(0); - auto propDaugTrack = fitter.getTrack(1); - kinkCand.decVtx = fitter.getPCACandidatePos(); + int nCand = 0; + try { + nCand = fitter.process(trackParCovMoth, trackParCovDaug); + } catch (...) { + LOG(error) << "Exception caught in DCA fitter process call!"; + return; + } + if (nCand == 0) { + return; + } + hSelKinkedTrackQA->Fill(6.f, chargeCombSvCand); // DCAFitter cut bin - // cut on decay radius to 17 cm - float decRad2 = kinkCand.decVtx[0] * kinkCand.decVtx[0] + kinkCand.decVtx[1] * kinkCand.decVtx[1]; - if (doSVRadiusCut && decRad2 < LayerRadii[3] * LayerRadii[3]) { - continue; - } + if (!fitter.propagateTracksToVertex()) { + return; + } + hSelKinkedTrackQA->Fill(7.f, chargeCombSvCand); // PropagateTracksToVertex cut bin - // get last layer hitted by the mother and the first layer hitted by the daughter - int lastLayerMoth = 0, firstLayerDaug = 0; - for (int i = 0; i < 7; i++) { - if (trackMoth.itsClusterMap() & (1 << i)) { - lastLayerMoth = i; - } - } + auto propMothTrack = fitter.getTrack(0); + auto propDaugTrack = fitter.getTrack(1); + kinkCand.decVtx = fitter.getPCACandidatePos(); - for (int i = 0; i < 7; i++) { - if (trackDaug.itsClusterMap() & (1 << i)) { - firstLayerDaug = i; - break; - } - } + // cut on decay radius to 17 cm + float decRad2 = kinkCand.decVtx[0] * kinkCand.decVtx[0] + kinkCand.decVtx[1] * kinkCand.decVtx[1]; + hMothDecRad2->Fill(decRad2, chargeCombSvCand); + if (doSVRadiusCut && decRad2 < LayerRadii[3] * LayerRadii[3]) { + return; + } + hSelKinkedTrackQA->Fill(8.f, chargeCombSvCand); // MothDecRad2InsideL4 cut bin - if (doSVRadiusCut && lastLayerMoth >= firstLayerDaug) { - continue; + // get last layer hitted by the mother and the first layer hitted by the daughter + int lastLayerMoth = 0, firstLayerDaug = 0; + for (int i = 0; i < 7; i++) { + if (trackMoth.itsClusterMap() & (1 << i)) { + lastLayerMoth = i; } + } - if (doSVRadiusCut && decRad2 < LayerRadii[lastLayerMoth] * LayerRadii[lastLayerMoth]) { - continue; + for (int i = 0; i < 7; i++) { + if (trackDaug.itsClusterMap() & (1 << i)) { + firstLayerDaug = i; + break; } + } - for (int i = 0; i < 3; i++) { - kinkCand.decVtx[i] -= kinkCand.primVtx[i]; - } + if (doSVRadiusCut && lastLayerMoth >= firstLayerDaug) { + return; + } + hSelKinkedTrackQA->Fill(9.f, chargeCombSvCand); // MothDaugLastLayers cut bin - propMothTrack.getPxPyPzGlo(kinkCand.momMoth); - propDaugTrack.getPxPyPzGlo(kinkCand.momDaug); - for (int i = 0; i < 3; i++) { - kinkCand.momMoth[i] *= charge; - kinkCand.momDaug[i] *= charge; - } - float pMoth = propMothTrack.getP() * charge; - float pDaug = propDaugTrack.getP() * charge; - float spKink = kinkCand.momMoth[0] * kinkCand.momDaug[0] + kinkCand.momMoth[1] * kinkCand.momDaug[1] + kinkCand.momMoth[2] * kinkCand.momDaug[2]; - kinkCand.kinkAngle = std::acos(spKink / (pMoth * pDaug)); - - std::array neutDauMom{0.f, 0.f, 0.f}; - for (int i = 0; i < 3; i++) { - neutDauMom[i] = kinkCand.momMoth[i] - kinkCand.momDaug[i]; - } + if (doSVRadiusCut && decRad2 < LayerRadii[lastLayerMoth] * LayerRadii[lastLayerMoth]) { + return; + } + hSelKinkedTrackQA->Fill(10.f, chargeCombSvCand); // MothLastLayerRadiusCheck cut bin - float chargedDauE = std::sqrt(pDaug * pDaug + chargedDauMass * chargedDauMass); - float neutE = std::sqrt(neutDauMom[0] * neutDauMom[0] + neutDauMom[1] * neutDauMom[1] + neutDauMom[2] * neutDauMom[2] + neutDauMass * neutDauMass); - float invMass = std::sqrt((chargedDauE + neutE) * (chargedDauE + neutE) - (pMoth * pMoth)); - - h2DeDxDaugSel->Fill(trackDaug.tpcInnerParam() * trackDaug.sign(), trackDaug.tpcSignal()); - h2KinkAnglePt->Fill(trackMoth.pt() * charge * trackMoth.sign(), kinkCand.kinkAngle * radToDeg); - h2MothMassPt->Fill(trackMoth.pt() * charge * trackMoth.sign(), invMass); - h2ClsMapPtMoth->Fill(trackMoth.pt() * charge * trackMoth.sign(), trackMoth.itsClusterMap()); - h2ClsMapPtDaug->Fill(trackDaug.pt() * charge * trackDaug.sign(), trackDaug.itsClusterMap()); - - kinkCand.collisionID = collision.globalIndex(); - kinkCand.mothTrackID = trackMoth.globalIndex(); - kinkCand.daugTrackID = trackDaug.globalIndex(); - - kinkCand.dcaXYmoth = dcaInfoMoth[0]; - kinkCand.mothSign = trackMoth.sign(); - kinkCand.dcaXYdaug = dcaInfoDaug[0]; - kinkCand.dcaKinkTopo = std::sqrt(fitter.getChi2AtPCACandidate()); - kinkCandidates.push_back(kinkCand); + for (int i = 0; i < 3; i++) { + kinkCand.decVtx[i] -= kinkCand.primVtx[i]; + } + + propMothTrack.getPxPyPzGlo(kinkCand.momMoth); + propDaugTrack.getPxPyPzGlo(kinkCand.momDaug); + for (int i = 0; i < 3; i++) { + kinkCand.momMoth[i] *= charge; + kinkCand.momDaug[i] *= charge; + } + float pMoth = propMothTrack.getP() * charge; + float pDaug = propDaugTrack.getP() * charge; + float spKink = kinkCand.momMoth[0] * kinkCand.momDaug[0] + kinkCand.momMoth[1] * kinkCand.momDaug[1] + kinkCand.momMoth[2] * kinkCand.momDaug[2]; + kinkCand.kinkAngle = std::acos(spKink / (pMoth * pDaug)); + + std::array neutDauMom{0.f, 0.f, 0.f}; + for (int i = 0; i < 3; i++) { + neutDauMom[i] = kinkCand.momMoth[i] - kinkCand.momDaug[i]; } + + float chargedDauE = std::sqrt(pDaug * pDaug + chargedDauMass * chargedDauMass); + float neutE = std::sqrt(neutDauMom[0] * neutDauMom[0] + neutDauMom[1] * neutDauMom[1] + neutDauMom[2] * neutDauMom[2] + neutDauMass * neutDauMass); + float invMass = std::sqrt((chargedDauE + neutE) * (chargedDauE + neutE) - (pMoth * pMoth)); + + h2DeDxDaugSel->Fill(trackDaug.tpcInnerParam() * trackDaug.sign(), trackDaug.tpcSignal()); + h2KinkAnglePt->Fill(trackMoth.pt() * charge * trackMoth.sign(), kinkCand.kinkAngle * radToDeg); + h2MothMassPt->Fill(trackMoth.pt() * charge * trackMoth.sign(), invMass); + h2ClsMapPtMoth->Fill(trackMoth.pt() * charge * trackMoth.sign(), trackMoth.itsClusterMap()); + h2ClsMapPtDaug->Fill(trackDaug.pt() * charge * trackDaug.sign(), trackDaug.itsClusterMap()); + hMothDaughSignsFinal->Fill(trackMoth.sign(), trackDaug.sign()); + + kinkCand.collisionID = collision.globalIndex(); + kinkCand.mothTrackID = trackMoth.globalIndex(); + kinkCand.daugTrackID = trackDaug.globalIndex(); + + kinkCand.dcaXYmoth = dcaInfoMoth[0]; + kinkCand.mothSign = trackMoth.sign(); + kinkCand.dcaXYdaug = dcaInfoDaug[0]; + kinkCand.dcaKinkTopo = std::sqrt(fitter.getChi2AtPCACandidate()); + kinkCandidates.push_back(kinkCand); + isAccepted = true; + return; } void initCCDB(aod::BCs::iterator const& bc) @@ -451,11 +630,185 @@ struct kinkBuilder { LOG(info) << "Task initialized for run " << mRunNumber << " with magnetic field " << mBz << " kZG"; } - void process(aod::Collisions const& collisions, TracksFull const& tracks, aod::AmbiguousTracks const& ambiTracks, aod::BCs const& bcs) + template + void buildSvPool(const TColls& collisions, const TTracks& tracks, const TAmbiTracks& ambiguousTracks, const aod::BCs& bcs) { + svCreator.clearPools(); + svCreator.fillBC2Coll(collisions, bcs); + + for (const auto& track : tracks) { + if (std::abs(track.eta()) > etaMax) { + continue; + } + + bool isDaug = selectDaugTrack(track); + bool isMoth = selectMothTrack(track); + if (!isDaug && !isMoth) { + continue; + } + int pdgHypo = isMoth ? 1 : 0; + svCreator.appendTrackCand(track, collisions, pdgHypo, ambiguousTracks, bcs); + } + } + + void processData(aod::Collisions const& collisions, TracksFull const& tracks, aod::AmbiguousTracks const& ambiTracks, aod::BCs const& bcs) { + kinkCandidates.clear(); + // LOG(info) << "[kink] "; + // LOG(info) << "[kink] *****************************************"; + // LOG(info) << "[kink] Processing " << tracks.size() << " tracks and " << collisions.size() << " collisions"; + + buildSvPool(collisions, tracks, ambiTracks, bcs); + auto& kinkPool = svCreator.getSVCandPool(collisions, !unlikeSignBkg); + // LOG(info) << "[kink] Number of kink candidates in the pool: " << kinkPool.size(); + bool isAccepted = false; + for (const auto& svCand : kinkPool) { + buildKinkCand(svCand, tracks, isAccepted, collisions); + } + + // sort kinkCandidates by collisionID to allow joining with collision table + std::sort(kinkCandidates.begin(), kinkCandidates.end(), [](const kinkCandidate& a, const kinkCandidate& b) { return a.collisionID < b.collisionID; }); + for (const auto& kinkCand : kinkCandidates) { + if (fillDebugTable) { + outputDataTableUB(kinkCand.decVtx[0], kinkCand.decVtx[1], kinkCand.decVtx[2], + kinkCand.mothSign, kinkCand.momMoth[0], kinkCand.momMoth[1], kinkCand.momMoth[2], + kinkCand.momDaug[0], kinkCand.momDaug[1], kinkCand.momDaug[2], + kinkCand.dcaXYmoth, kinkCand.dcaXYdaug, kinkCand.dcaKinkTopo); + } else { + outputDataTable(kinkCand.collisionID, kinkCand.mothTrackID, kinkCand.daugTrackID, + kinkCand.decVtx[0], kinkCand.decVtx[1], kinkCand.decVtx[2], + kinkCand.mothSign, kinkCand.momMoth[0], kinkCand.momMoth[1], kinkCand.momMoth[2], + kinkCand.momDaug[0], kinkCand.momDaug[1], kinkCand.momDaug[2], + kinkCand.dcaXYmoth, kinkCand.dcaXYdaug, kinkCand.dcaKinkTopo); + } + } + } + PROCESS_SWITCH(kinkBuilder, processData, "Data processing", false); + + template + int matchKinkDecay(const TMother& motherPart, const aod::McParticles& mcParticles) { + int pdgMother = motherPart.pdgCode(); + int8_t sign = 0; + int pdgCodeNeutralDaug{-1}, pdgCodeChargedDaug{-1}; + std::array finState = {-1, -1}; + switch (std::abs(pdgMother)) { + case PDG_t::kSigmaMinus: { + // Swap the sign of the neutral decay products in case of anti-particles + pdgCodeNeutralDaug = (pdgMother > 0) ? +PDG_t::kNeutron : -PDG_t::kNeutron; + pdgCodeChargedDaug = (pdgMother > 0) ? +PDG_t::kPiMinus : +PDG_t::kPiPlus; + finState = {pdgCodeChargedDaug, pdgCodeNeutralDaug}; // Both decay channels have the same neutral daughter + // if constexpr (checkKinkDaugPdg) { + // LOG(info) << "[kink] Matching Sigma- decay with daughter PDG " << pdgDaug << " and neutral daughter PDG " << pdgCodeNeutralDaug; + // } + if (RecoDecay::isMatchedMCGen(mcParticles, motherPart, pdgMother, finState, true, &sign, DepthMcMatchMax)) { + // if constexpr (checkKinkDaugPdg) { + // LOG(info) << "[kink] Matched Sigma- to pi- neutron decay"; + // } + return kSigmaMinusToPiMinusNeutron; + } + break; + } + case PDG_t::kSigmaPlus: { + // Swap the sign of the neutral decay products in case of anti-particles + pdgCodeNeutralDaug = (pdgMother > 0) ? +PDG_t::kNeutron : -PDG_t::kNeutron; + pdgCodeChargedDaug = (pdgMother > 0) ? +PDG_t::kPiPlus : +PDG_t::kPiMinus; + finState = {pdgCodeChargedDaug, pdgCodeNeutralDaug}; + // if constexpr (checkKinkDaugPdg) { + // LOG(info) << "[kink] Matching Sigma+ decay with daughter PDG " << pdgDaug << " and neutral daughter PDG " << pdgCodeNeutralDaug; + // } + if (RecoDecay::isMatchedMCGen(mcParticles, motherPart, pdgMother, finState, true, &sign, DepthMcMatchMax)) { + // if constexpr (checkKinkDaugPdg) { + // LOG(info) << "[kink] Matched Sigma+ to pi+ neutron decay"; + // } + return kSigmaPlusToPiPlusNeutron; + } + + // Swap the sign of the neutral decay products in case of anti-particles + pdgCodeNeutralDaug = (pdgMother > 0) ? +PDG_t::kPi0 : -PDG_t::kPi0; + pdgCodeChargedDaug = (pdgMother > 0) ? +PDG_t::kProton : -PDG_t::kProton; + finState = {pdgCodeChargedDaug, pdgCodeNeutralDaug}; + // if constexpr (checkKinkDaugPdg) { + // LOG(info) << "[kink] Matching Sigma+ decay with daughter PDG " << pdgDaug << " and neutral daughter PDG " << pdgCodeNeutralDaug; + // } + if (RecoDecay::isMatchedMCGen(mcParticles, motherPart, pdgMother, finState, true, &sign, DepthMcMatchMax)) { + // if constexpr (checkKinkDaugPdg) { + // LOG(info) << "[kink] Matched Sigma+ to proton pi0 decay"; + // } + return kSigmaPlusToProtonPi0; + } + break; + } + default: + LOG(warning) << "No matching function implemented for the selected mother particle hypothesis. Returning -1."; + return -1; + } + + return -1; + } + + void processMc(aod::McCollisions const& mcGenCollisions, + McRecoCollisionsCent const& mcRecoCollisions, + aod::McParticles const& mcParticles, + TracksFullMc const& tracksMc, + aod::AmbiguousTracks const& ambiTracksMc, + aod::BCs const& bcs) + { kinkCandidates.clear(); - fillCandidateData(collisions, tracks, ambiTracks, bcs); + // LOG(info) << "[kink] "; + // LOG(info) << "[kink] *****************************************"; + // LOG(info) << "[kink] Processing " << tracksMc.size() << " tracks and " << mcRecoCollisions.size() << " collisions"; + + buildSvPool(mcRecoCollisions, tracksMc, ambiTracksMc, bcs); + auto& kinkPool = svCreator.getSVCandPool(mcRecoCollisions, !unlikeSignBkg); + // LOG(info) << "[kink] Number of kink candidates in the pool: " << kinkPool.size(); + int countSameMother = 0; + for (const auto& svCand : kinkPool) { + // Perform matching of the kink candidate + auto trackMoth = tracksMc.rawIteratorAt(svCand.tr0Idx); + auto genMothPart = trackMoth.mcParticle_as(); + auto trackDaug = tracksMc.rawIteratorAt(svCand.tr1Idx); + auto genDaugPart = trackDaug.mcParticle_as(); + int decayChannel{-1}; + if (skipBkgCands) { + + // Check mother PDG first + if (std::find(mothMatchPdgCodes.begin(), mothMatchPdgCodes.end(), std::abs(genMothPart.pdgCode())) == mothMatchPdgCodes.end()) { + continue; + } + + // Check that the daughter was generated from the selected mother + int genDaugMothIdx = RecoDecay::getMother(mcParticles, genDaugPart, genMothPart.pdgCode(), true); + if (genDaugMothIdx != trackMoth.mcParticleId()) { + continue; // Skip candidates where the daughter is not coming from the selected mother + } + // LOG(info) << "[kink] "; + // LOG(info) << "[kink] -----------------"; + // LOG(info) << "[kink] MC particle offset: " << mcParticles.offset(); + // LOG(info) << "[kink] [Rec] Matching kink candidate with mother PDG " << genMothPart.pdgCode() << " and daughter PDG " << genDaugPart.pdgCode(); + // LOG(info) << "[kink] Mother global index: " << trackMoth.globalIndex() << ", trackMoth.mcParticleId(): " << trackMoth.mcParticleId() << ", PDG code: " << genMothPart.pdgCode(); + // for (const auto& mcMother : genDaugPart.mothers_as()) { + // LOG(info) << "[kink] Daughter's mother PDG: " << mcMother.pdgCode() << ", global index: " << mcMother.globalIndex(); + // } + // std::vector arrDaughIdxs = {}; + // RecoDecay::getDaughters(genMothPart, &arrDaughIdxs, std::array{0}, 1); + // for (auto iProng = 0u; iProng < arrDaughIdxs.size(); ++iProng) { + // auto daughI = mcParticles.rawIteratorAt(arrDaughIdxs[iProng]); + // LOG(info) << "[kink] Daughter PDG: " << daughI.pdgCode() << ", global index: " << daughI.globalIndex(); + // } + decayChannel = matchKinkDecay(genMothPart, mcParticles); + // LOG(info) << "[kink] Decay channel matched: " << decayChannel; + if (decayChannel < 0) { + continue; // Skip candidates that do not match the decay channels of interest + } + hRecCandidates->Fill(decayChannel, trackMoth.pt()); // Decay channel match bin + } + bool isAccepted = false; + buildKinkCand(svCand, tracksMc, isAccepted, mcRecoCollisions); + // If candidate passed all selections and was built, fill histogram + if (isAccepted && skipBkgCands) { + hRecPtKinkAngle[decayChannel]->Fill(trackMoth.pt(), std::abs(trackMoth.phi() - trackDaug.phi()) * radToDeg); + } + } // sort kinkCandidates by collisionID to allow joining with collision table std::sort(kinkCandidates.begin(), kinkCandidates.end(), [](const kinkCandidate& a, const kinkCandidate& b) { return a.collisionID < b.collisionID; }); @@ -474,7 +827,28 @@ struct kinkBuilder { kinkCand.dcaXYmoth, kinkCand.dcaXYdaug, kinkCand.dcaKinkTopo); } } + + // Generated kinks + for (const auto& mcPart : mcParticles) { + if (std::find(mothMatchPdgCodes.begin(), mothMatchPdgCodes.end(), std::abs(mcPart.pdgCode())) == mothMatchPdgCodes.end()) { + continue; // Skip if mother particle does not match the PDG codes of interest + } + int decayChannel = matchKinkDecay(mcPart, mcParticles); // Dummy daughter PDG code since we are not checking it in this case + if (decayChannel < 0) { + continue; // Skip if no matching decay channel is found + } + hGenCandidates->Fill(decayChannel, mcPart.pt()); + std::vector arrDaughIdxs = {}; + RecoDecay::getDaughters(mcPart, &arrDaughIdxs, std::array{0}, DepthMcMatchMax); + for (auto iProng = 0u; iProng < arrDaughIdxs.size(); ++iProng) { + auto daughI = mcParticles.rawIteratorAt(arrDaughIdxs[iProng]); + if (std::abs(daughI.pdgCode()) == PDG_t::kPiPlus || std::abs(daughI.pdgCode()) == PDG_t::kProton) { + hGenPtKinkAngle[decayChannel]->Fill(mcPart.pt(), std::abs(mcPart.phi() - daughI.phi()) * radToDeg); + } + } + } } + PROCESS_SWITCH(kinkBuilder, processMc, "MC processing", false); }; WorkflowSpec diff --git a/PWGLF/Tasks/Resonances/lambda1405analysis.cxx b/PWGLF/Tasks/Resonances/lambda1405analysis.cxx index 7c12c2b9042..0511d02bcc9 100644 --- a/PWGLF/Tasks/Resonances/lambda1405analysis.cxx +++ b/PWGLF/Tasks/Resonances/lambda1405analysis.cxx @@ -17,9 +17,11 @@ #include "PWGLF/DataModel/LFLambda1405Table.h" #include "Common/Core/RecoDecay.h" +#include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/PIDResponseTOF.h" #include "Common/DataModel/PIDResponseTPC.h" +#include "Common/DataModel/Qvectors.h" #include #include @@ -43,19 +45,28 @@ using namespace o2::framework::expressions; using TracksFull = soa::Join; using CollisionsFull = soa::Join; +using CollisionsFullWCentQVecs = soa::Join; using CollisionsFullMC = soa::Join; +enum L1405DecayType { kSigmaMinusPiToPiPiNeutron = 0, + kSigmaPlusPiToPiPiNeutron, + kSigmaPlusPiToPiPiProton, + kNL1405Decays }; + struct lambda1405candidate { // Columns for Lambda(1405) candidate - float mass = -1; // Invariant mass of the Lambda(1405) candidate + float massL1405 = -1; // Invariant mass of the Lambda(1405) candidate float massXi1530 = -1; // Invariant mass of the Xi(1530) candidate float px = -1; // Px of the Lambda(1405) candidate float py = -1; // Py of the Lambda(1405) candidate float pz = -1; // Pz of the Lambda(1405) candidate - float pt() const { return std::sqrt(px * px + py * py); } // pT of the Lambda(1405 candidate + float pt() const { return std::sqrt(px * px + py * py); } // pT of the Lambda(1405) candidate + float phi = -1; // Phi of the Lambda(1405) candidate bool isSigmaPlus = false; // True if compatible with Sigma+ bool isSigmaMinus = false; // True if compatible with Sigma- + bool hasPiKink = false; // True if the Sigma candidate has a kink topology compatible with a pion + bool hasPrKink = false; // True if the Sigma candidate has a kink topology compatible with a proton float sigmaMinusMass = -1; // Invariant mass of the Sigma- candidate float sigmaPlusMass = -1; // Invariant mass of the Sigma+ candidate float xiMinusMass = -1; // Invariant mass of the Xi- candidate @@ -64,46 +75,85 @@ struct lambda1405candidate { float sigmaAlphaAP = -1; // Alpha of the Sigma float sigmaQtAP = -1; // qT of the Sigma float kinkPt = -1; // pT of the kink daughter - float kinkTPCNSigmaPi = -1; // Number of sigmas for the pion candidate from Sigma kink in TPC - float kinkTOFNSigmaPi = -1; // Number of sigmas for the pion candidate from Sigma kink in TOF - float kinkTPCNSigmaPr = -1; // Number of sigmas for the proton candidate from Sigma kink in TPC - float kinkTOFNSigmaPr = -1; // Number of sigmas for the proton candidate from Sigma kink in TOF - float dcaKinkDauToPV = -1; // DCA of the kink daughter to the primary vertex + float kinkPiNSigTpc = -1; // Number of sigmas for the pion candidate from Sigma kink in Tpc + float kinkPiNSigTof = -1; // Number of sigmas for the pion candidate from Sigma kink in Tof + float kinkPrNSigTpc = -1; // Number of sigmas for the proton candidate from Sigma kink in Tpc + float kinkPrNSigTof = -1; // Number of sigmas for the proton candidate from Sigma kink in Tof + float kinkDcaDauToPv = -1; // DCA of the kink daughter to the primary vertex float sigmaRadius = -1; // Radius of the Sigma decay vertex float piPt = -1; // pT of the pion daughter - float nSigmaTPCPi = -1; // Number of sigmas for the pion candidate - float nSigmaTOFPi = -1; // Number of sigmas for the pion candidate using TOF - int kinkDauID = 0; // ID of the pion from Sigma decay in MC - int sigmaID = 0; // ID of the Sigma candidate in MC - int piID = 0; // ID of the pion candidate in MC + float bachPiNSigTpc = -1; // Number of sigmas for the pion candidate + float bachPiNSigTof = -1; // Number of sigmas for the pion candidate using Tof + int kinkDauId = 0; // Id of the pion from Sigma decay in MC + int sigmaId = 0; // Id of the Sigma candidate in MC + int piId = 0; // Id of the pion candidate in MC + + float scalarProd = -1; // Scalar product for flow analysis }; struct lambda1405analysis { int lambda1405PdgCode = 102132; // PDG code for Lambda(1405) lambda1405candidate lambda1405Cand; // Lambda(1405) candidate structure Produces outputDataTable; // Output table for Lambda(1405) candidates + Produces outputDataFlowTable; // Output table for Lambda(1405) flow analysis Produces outputDataTableMC; // Output table for Lambda(1405) candidates in MC // Histograms are defined with HistogramRegistry HistogramRegistry rEventSelection{"eventSelection", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; HistogramRegistry rLambda1405{"lambda1405", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + HistogramRegistry rSigmaMinus{"sigmaMinus", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + HistogramRegistry rSigmaPlus{"sigmaPlus", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + HistogramRegistry rSelections{"selections", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; // Configurable for event selection Configurable cutzvertex{"cutZVertex", 10.0f, "Accepted z-vertex range (cm)"}; Configurable cutEtaDaught{"cutEtaDaughter", 0.8f, "Eta cut for daughter tracks"}; - Configurable cutDCAtoPVSigma{"cutDCAtoPVSigma", 0.1f, "Max DCA to primary vertex for Sigma candidates (cm)"}; - Configurable cutDCAtoPVPiFromSigma{"cutDCAtoPVPiFromSigma", 2., "Min DCA to primary vertex for pion from Sigma candidates (cm)"}; + Configurable cutDCAtoPvSigma{"cutDCAtoPvSigma", 0.1f, "Max DCA to primary vertex for Sigma candidates (cm)"}; + Configurable cutDCAtoPvPiFromSigma{"cutDCAtoPvPiFromSigma", 2., "Min DCA to primary vertex for pion from Sigma candidates (cm)"}; Configurable cutUpperMass{"cutUpperMass", 1.6f, "Upper mass cut for Lambda(1405) candidates (GeV/c^2)"}; Configurable cutSigmaRadius{"cutSigmaRadius", 20.f, "Minimum radius for Sigma candidates (cm)"}; Configurable cutSigmaMass{"cutSigmaMass", 0.1, "Sigma mass window (MeV/c^2)"}; - Configurable cutNITSClusPi{"cutNITSClusPi", 5, "Minimum number of ITS clusters for pion candidate"}; - Configurable cutNTPCClusPi{"cutNTPCClusPi", 90, "Minimum number of TPC clusters for pion candidate"}; - Configurable cutNSigmaTPC{"cutNSigmaTPC", 3, "NSigmaTPCPion"}; - Configurable cutNSigmaTOF{"cutNSigmaTOF", 3, "NSigmaTOFPion"}; + Configurable cutNITSClusKink{"cutNITSClusKink", 3, "Minimum number of ITS clusters for pion candidate"}; + Configurable cutNTpcClusPi{"cutNTpcClusPi", 90, "Minimum number of Tpc clusters for pion candidate"}; + Configurable cutNSigTpc{"cutNSigTpc", 3, "NSigTpcPion"}; + Configurable cutNSigTof{"cutNSigTof", 3, "NSigTofPion"}; + Configurable cutSigmaQtAPMin{"cutSigmaQtAPMin", 0.17, "Lower limit for Sigma qT"}; + Configurable cutSigmaQtAPMax{"cutSigmaQtAPMax", 0.2, "Upper limit for Sigma qT"}; + Configurable cutSigmaAlphaAPMin{"cutSigmaAlphaAPMin", -0.9, "Lower limit for Sigma qT"}; + Configurable cutSigmaAlphaAPMax{"cutSigmaAlphaAPMax", -0.4, "Upper limit for Sigma qT"}; + + // Configurables for flow analysis + Configurable centralityMin{"centralityMin", 0., "Minimum centrality accepted in SP/EP computation (not applied in resolution process)"}; + Configurable centralityMax{"centralityMax", 100., "Maximum centrality accepted in SP/EP computation (not applied in resolution process)"}; + Configurable fillFlowTree{"fillFlowTree", false, "If true, fill the output tree with Lambda(1405) candidates"}; + + // Downsample for table filling + Configurable downSampleFactor{"downSampleFactor", 1., "Fraction of candidates to keep in TTree"}; + Configurable ptDownSampleMax{"ptDownSampleMax", 10., "Maximum pt for the application of the downsampling factor"}; Configurable fillOutputTree{"fillOutputTree", true, "If true, fill the output tree with Lambda(1405) candidates"}; Configurable doLSBkg{"doLikeSignBkg", false, "Use like-sign background"}; - Configurable useTOF{"useTOF", false, "Use TOF for PID for pion candidates"}; + Configurable useTof{"useTof", false, "Use Tof for PId for pion candidates"}; + + // Configurable axes + ConfigurableAxis axisCent{"axisCent", {10000, 0., 100.}, ""}; + ConfigurableAxis axisPtL1405{"axisPtL1405", {100, 0., 10.}, ""}; + ConfigurableAxis axisPCompsL1405{"axisPCompsL1405", {100, -10., 10.}, ""}; + ConfigurableAxis axisPtKinkDaug{"axisPtKinkDaug", {100, -5., 5.}, ""}; + ConfigurableAxis axisPtResolution{"axisPtResolution", {100, -0.5, 0.5}, ""}; + ConfigurableAxis axisMassL1405{"axisMassL1405", {100, 1.3, 1.8}, ""}; + ConfigurableAxis axisXi1530Mass{"axisXi1530Mass", {100, 1.4, 1.6}, ""}; + ConfigurableAxis axisMassResolution{"axisMassResolution", {100, -0.1, 0.1}, ""}; + ConfigurableAxis axisNSig{"axisNSig", {100, -5., 5.}, ""}; + ConfigurableAxis axisSigmaMass{"axisSigmaMass", {100, 1.1, 1.4}, ""}; + ConfigurableAxis axisXiMass{"axisXiMass", {100, 1.2, 1.5}, ""}; + ConfigurableAxis axisVertexZ{"axisVertexZ", {100, -15., 15.}, ""}; + ConfigurableAxis axisQtAP{"axisQtAP", {100, 0., 0.3}, ""}; + ConfigurableAxis axisAlphaAP{"axisAlphaAP", {200, -1., 1.}, ""}; + ConfigurableAxis axisSigmaRadius{"axisSigmaRadius", {100, 0., 100.}, ""}; + ConfigurableAxis axisDcaSigmaToPv{"axisDcaSigmaToPv", {120, 0., 0.05}, ""}; + ConfigurableAxis axisDcaKinkToPv{"axisDcaKinkToPv", {200, 0., 20.}, ""}; + ConfigurableAxis axisScalarProd{"axisScalarProd", {200, -4., 4.}, ""}; Preslice mKinkPerCol = aod::track::collisionId; Preslice mPerColTracks = aod::track::collisionId; @@ -111,40 +161,158 @@ struct lambda1405analysis { void init(InitContext const&) { // Axes - const AxisSpec ptAxis{100, -10, 10, "#it{p}_{T} (GeV/#it{c})"}; - const AxisSpec ptPiAxis{100, -5, 5, "#it{p}_{T}^{#pi} (GeV/#it{c})"}; - const AxisSpec ptResolutionAxis{100, -0.5, 0.5, "#it{p}_{T}^{rec} - #it{p}_{T}^{gen} (GeV/#it{c})"}; - const AxisSpec massAxis{100, 1.3, 1.6, "m (GeV/#it{c}^{2})"}; - const AxisSpec massResolutionAxis{100, -0.1, 0.1, "m_{rec} - m_{gen} (GeV/#it{c}^{2})"}; - const AxisSpec nSigmaPiAxis{100, -5, 5, "n#sigma_{#pi}"}; - const AxisSpec sigmaMassAxis{100, 1.1, 1.4, "m (GeV/#it{c}^{2})"}; - const AxisSpec vertexZAxis{100, -15., 15., "vrtx_{Z} [cm]"}; + const AxisSpec ptAxis{axisPtL1405, "#it{p}_{T} (GeV/#it{c})"}; + const AxisSpec pCompAxisL1405{axisPCompsL1405, "#it{p}_{comp} (GeV/#it{c})"}; + const AxisSpec ptKinkDaugAxis{axisPtKinkDaug, "#it{p}_{T}^{#pi} (GeV/#it{c})"}; + const AxisSpec ptResolutionAxis{axisPtResolution, "#it{p}_{T}^{rec} - #it{p}_{T}^{gen} (GeV/#it{c})"}; + const AxisSpec lambda1405MassAxis{axisMassL1405, "m (GeV/#it{c}^{2})"}; + const AxisSpec xi1530MassAxis{axisXi1530Mass, "m (GeV/#it{c}^{2})"}; + const AxisSpec massResolutionAxis{axisMassResolution, "m_{rec} - m_{gen} (GeV/#it{c}^{2})"}; + const AxisSpec nSigmaAxis{axisNSig, "n#sigma_{#pi}"}; + const AxisSpec sigmaMassAxis{axisSigmaMass, "m (GeV/#it{c}^{2})"}; + const AxisSpec xiMassAxis{axisXiMass, "m (GeV/#it{c}^{2})"}; + const AxisSpec vertexZAxis{axisVertexZ, "vrtx_{Z} [cm]"}; + const AxisSpec qtAxis{axisQtAP, "q_{T, AP}"}; + const AxisSpec alphaAxis{axisAlphaAP, "#alpha_{AP}"}; + const AxisSpec sigmaRadiusAxis{axisSigmaRadius, "#Sigma radius (cm)"}; + const AxisSpec centAxis{axisCent, "Centrality"}; + const AxisSpec dcaSigmaToPvBinsAxis{axisDcaSigmaToPv, "DCA of mother to Pv"}; + const AxisSpec dcaKinkToPvBinsAxis{axisDcaKinkToPv, "DCA of kink to Pv"}; + const AxisSpec scalarProdAxis{axisScalarProd, "SP"}; + rEventSelection.add("hVertexZRec", "hVertexZRec", {HistType::kTH1F, {vertexZAxis}}); - // lambda1405 to sigmaminus - rLambda1405.add("h2PtMass_0", "h2PtMass_0", {HistType::kTH2F, {ptAxis, massAxis}}); - rLambda1405.add("h2PtMassSigmaBeforeCuts_0", "h2PtMassSigmaBeforeCuts_0", {HistType::kTH2F, {ptAxis, sigmaMassAxis}}); - rLambda1405.add("h2PtMassSigma_0", "h2PtMassSigma_0", {HistType::kTH2F, {ptAxis, sigmaMassAxis}}); - rLambda1405.add("h2SigmaMassVsMass_0", "h2SigmaMassVsMass_0", {HistType::kTH2F, {massAxis, sigmaMassAxis}}); - rLambda1405.add("h2PtPiNSigma_0", "h2PtPiNSigma_0", {HistType::kTH2F, {ptPiAxis, nSigmaPiAxis}}); - rLambda1405.add("h2PtPiNSigmaTOF_0", "h2PtPiNSigmaTOF_0", {HistType::kTH2F, {ptPiAxis, nSigmaPiAxis}}); - // lambda1405 to sigmaplus - rLambda1405.add("h2PtMass_1", "h2PtMass_1", {HistType::kTH2F, {ptAxis, massAxis}}); - rLambda1405.add("h2PtMassSigmaBeforeCuts_1", "h2PtMassSigmaBeforeCuts_1", {HistType::kTH2F, {ptAxis, sigmaMassAxis}}); - rLambda1405.add("h2PtMassSigma_1", "h2PtMassSigma_1", {HistType::kTH2F, {ptAxis, sigmaMassAxis}}); - rLambda1405.add("h2SigmaMassVsMass_1", "h2SigmaMassVsMass_1", {HistType::kTH2F, {massAxis, sigmaMassAxis}}); - rLambda1405.add("h2PtPiNSigma_1", "h2PtPiNSigma_1", {HistType::kTH2F, {ptPiAxis, nSigmaPiAxis}}); - rLambda1405.add("h2PtPiNSigmaTOF_1", "h2PtPiNSigmaTOF_1", {HistType::kTH2F, {ptPiAxis, nSigmaPiAxis}}); + // Sigma- candidate properties + rSigmaMinus.add("hSigmaMinusMass", "hSigmaMinusMass", {HistType::kTH1D, {sigmaMassAxis}}); + rSigmaMinus.add("hSigmaPlusMass", "hSigmaPlusMass", {HistType::kTH1D, {sigmaMassAxis}}); + rSigmaMinus.add("hSigmaMinusPt", "hSigmaMinusPt", {HistType::kTH1D, {ptAxis}}); + rSigmaMinus.add("hMassXiMinusSigmaMinus", "hMassXiMinusSigmaMinus", {HistType::kTH1D, {xiMassAxis}}); + rSigmaMinus.add("h2PtMassSigmaMinusBeforeCuts", "h2PtMassSigmaMinusBeforeCuts", {HistType::kTH2F, {ptAxis, sigmaMassAxis}}); + rSigmaMinus.add("h2PtPiKinkNSigBeforeCutsSigmaMinus", "h2PtPiKinkNSigBeforeCutsSigmaMinus", {HistType::kTH2F, {ptAxis, nSigmaAxis}}); + rSigmaMinus.add("h2dPtMassSigmaMinus", "h2dPtMassSigmaMinus", {HistType::kTH2F, {ptAxis, sigmaMassAxis}}); + rSigmaMinus.add("h2SigmaMinusMassVsLambdaMass", "h2SigmaMinusMassVsLambdaMass", {HistType::kTH2F, {lambda1405MassAxis, sigmaMassAxis}}); + rSigmaMinus.add("h2KinkPiPtNSigTofSigmaMinus", "h2KinkPiPtNSigTofSigmaMinus", {HistType::kTH2F, {ptKinkDaugAxis, nSigmaAxis}}); + rSigmaMinus.add("h2KinkPrPtNSigTofSigmaMinus", "h2KinkPrPtNSigTofSigmaMinus", {HistType::kTH2F, {ptKinkDaugAxis, nSigmaAxis}}); + rSigmaMinus.add("hSigmaMinusArmPod", "hSigmaMinusArmPod", {HistType::kTH2D, {alphaAxis, qtAxis}}); + rSigmaMinus.add("hSigmaMinusRadius", "hSigmaMinusRadius", {HistType::kTH1D, {sigmaRadiusAxis}}); + rSigmaMinus.add("hSigmaMinusDcaToPv", "hSigmaMinusDcaToPv", {HistType::kTH1D, {dcaSigmaToPvBinsAxis}}); + rSigmaMinus.add("hSigmaMinusKinkPt", "hSigmaMinusKinkPt", {HistType::kTH1D, {ptKinkDaugAxis}}); + rSigmaMinus.add("hSigmaMinusKinkTpcNSigPi", "hSigmaMinusKinkTpcNSigPi", {HistType::kTH1D, {nSigmaAxis}}); + rSigmaMinus.add("hSigmaMinusKinkTofNSigPi", "hSigmaMinusKinkTofNSigPi", {HistType::kTH1D, {nSigmaAxis}}); + rSigmaMinus.add("hSigmaMinusDcaKinkDauToPv", "hSigmaMinusDcaKinkDauToPv", {HistType::kTH1D, {dcaKinkToPvBinsAxis}}); + + // Sigma+ candidate properties + rSigmaPlus.add("hSigmaMinusMass", "hSigmaMinusMass", {HistType::kTH1D, {sigmaMassAxis}}); + rSigmaPlus.add("hSigmaPlusMass", "hSigmaPlusMass", {HistType::kTH1D, {sigmaMassAxis}}); + rSigmaPlus.add("hSigmaPlusPt", "hSigmaPlusPt", {HistType::kTH1D, {ptAxis}}); + rSigmaPlus.add("h2PtMassSigmaPlusBeforeCuts", "h2PtMassSigmaPlusBeforeCuts", {HistType::kTH2F, {ptAxis, sigmaMassAxis}}); + rSigmaPlus.add("h2PtPiKinkNSigBeforeCutsSigmaPlus", "h2PtPiKinkNSigBeforeCutsSigmaPlus", {HistType::kTH2F, {ptAxis, nSigmaAxis}}); + rSigmaPlus.add("h2PtPrKinkNSigBeforeCutsSigmaPlus", "h2PtPrKinkNSigBeforeCutsSigmaPlus", {HistType::kTH2F, {ptAxis, nSigmaAxis}}); + rSigmaPlus.add("h2dPtMassSigmaPlus", "h2dPtMassSigmaPlus", {HistType::kTH2F, {ptAxis, sigmaMassAxis}}); + rSigmaPlus.add("h2SigmaPlusMassVsLambdaMass", "h2SigmaPlusMassVsLambdaMass", {HistType::kTH2F, {lambda1405MassAxis, sigmaMassAxis}}); + rSigmaPlus.add("h2KinkPrPtNSigTofSigmaPlus", "h2KinkPrPtNSigTofSigmaPlus", {HistType::kTH2F, {ptKinkDaugAxis, nSigmaAxis}}); + rSigmaPlus.add("hSigmaPlusArmPod", "hSigmaPlusArmPod", {HistType::kTH2D, {alphaAxis, qtAxis}}); + rSigmaPlus.add("hSigmaPlusRadius", "hSigmaPlusRadius", {HistType::kTH1D, {sigmaRadiusAxis}}); + rSigmaPlus.add("hSigmaPlusDcaToPv", "hSigmaPlusDcaToPv", {HistType::kTH1D, {dcaSigmaToPvBinsAxis}}); + rSigmaPlus.add("hSigmaPlusKinkPt", "hSigmaPlusKinkPt", {HistType::kTH1D, {ptKinkDaugAxis}}); + rSigmaPlus.add("hSigmaPlusKinkTpcNSigPi", "hSigmaPlusKinkTpcNSigPi", {HistType::kTH1D, {nSigmaAxis}}); + rSigmaPlus.add("hSigmaPlusKinkTofNSigPi", "hSigmaPlusKinkTofNSigPi", {HistType::kTH1D, {nSigmaAxis}}); + rSigmaPlus.add("hSigmaPlusKinkTpcNSigPr", "hSigmaPlusKinkTpcNSigPr", {HistType::kTH1D, {nSigmaAxis}}); + rSigmaPlus.add("hSigmaPlusKinkTofNSigPr", "hSigmaPlusKinkTofNSigPr", {HistType::kTH1D, {nSigmaAxis}}); + rSigmaPlus.add("hSigmaPlusDcaKinkDauToPv", "hSigmaPlusDcaKinkDauToPv", {HistType::kTH1D, {dcaKinkToPvBinsAxis}}); + rSigmaPlus.add("hMassXiMinusSigmaPlus", "hMassXiMinusSigmaPlus", {HistType::kTH1D, {xiMassAxis}}); + + // Mass QA + rLambda1405.add("hMassL1405", "hMassL1405", {HistType::kTH1D, {lambda1405MassAxis}}); + rLambda1405.add("hMassXi1530", "hMassXi1530", {HistType::kTH1D, {xi1530MassAxis}}); + // Kinematic distributions + rLambda1405.add("hPx", "hPx;#it{p}_x;Counts", {HistType::kTH1D, {pCompAxisL1405}}); + rLambda1405.add("hPy", "hPy;#it{p}_y;Counts", {HistType::kTH1D, {pCompAxisL1405}}); + rLambda1405.add("hPz", "hPz;#it{p}_z;Counts", {HistType::kTH1D, {pCompAxisL1405}}); + rLambda1405.add("hPt", "hPt", {HistType::kTH1D, {ptAxis}}); + rLambda1405.add("hPhi", "hPhi", {HistType::kTH1D, {{128, -o2::constants::math::PI, o2::constants::math::PI}}}); + // Pion daughter properties + rLambda1405.add("hBachPiPt", "hBachPiPt", {HistType::kTH1D, {ptKinkDaugAxis}}); + rLambda1405.add("hBachPiNSigTpc", "hBachPiNSigTpc", {HistType::kTH1D, {nSigmaAxis}}); + rLambda1405.add("hBachPiNSigTof", "hBachPtNSigTof", {HistType::kTH1D, {nSigmaAxis}}); + rLambda1405.add("h2BachPiPtNSigTof", "h2BachPiPtNSigTof", {HistType::kTH2F, {ptKinkDaugAxis, nSigmaAxis}}); + rLambda1405.add("h2BachPiPtNSigTpc", "h2BachPiPtNSigTpc", {HistType::kTH2F, {ptKinkDaugAxis, nSigmaAxis}}); + + // Selection QA + rSelections.add("hSelectionsL1405", "hSelectionsL1405", {HistType::kTH1D, {{5, -0.f, 4.5f}}}); + rSelections.get(HIST("hSelectionsL1405"))->GetXaxis()->SetBinLabel(1, "All"); + rSelections.get(HIST("hSelectionsL1405"))->GetXaxis()->SetBinLabel(2, "Passed Sigma sel"); + rSelections.get(HIST("hSelectionsL1405"))->GetXaxis()->SetBinLabel(3, "Passed Bach PID sel"); + rSelections.get(HIST("hSelectionsL1405"))->GetXaxis()->SetBinLabel(4, "Upper mass sel"); + rSelections.get(HIST("hSelectionsL1405"))->GetXaxis()->SetBinLabel(5, "Accepted"); + rSelections.add("hSelectionsSigmaPlus", "hSelectionsSigmaPlus", {HistType::kTH1D, {{8, -0.f, 7.5f}}}); + rSelections.get(HIST("hSelectionsSigmaPlus"))->GetXaxis()->SetBinLabel(1, "All"); + rSelections.get(HIST("hSelectionsSigmaPlus"))->GetXaxis()->SetBinLabel(2, "Passed kink sel"); + rSelections.get(HIST("hSelectionsSigmaPlus"))->GetXaxis()->SetBinLabel(3, "Passed mass sel"); + rSelections.get(HIST("hSelectionsSigmaPlus"))->GetXaxis()->SetBinLabel(4, "Passed cutDCAtoPvSigma"); + rSelections.get(HIST("hSelectionsSigmaPlus"))->GetXaxis()->SetBinLabel(5, "Passed cutDCAtoPvPiFromSigma"); + rSelections.get(HIST("hSelectionsSigmaPlus"))->GetXaxis()->SetBinLabel(6, "Passed cutSigmaRadius"); + rSelections.get(HIST("hSelectionsSigmaPlus"))->GetXaxis()->SetBinLabel(7, "Passed cutSigmaQtAP"); + rSelections.get(HIST("hSelectionsSigmaPlus"))->GetXaxis()->SetBinLabel(8, "Passed cutSigmaAlphaAP"); + rSelections.add("hSelectionsSigmaMinus", "hSelectionsSigmaMinus", {HistType::kTH1D, {{8, -0.f, 7.5f}}}); + rSelections.get(HIST("hSelectionsSigmaMinus"))->GetXaxis()->SetBinLabel(1, "All"); + rSelections.get(HIST("hSelectionsSigmaMinus"))->GetXaxis()->SetBinLabel(2, "Passed kink sel"); + rSelections.get(HIST("hSelectionsSigmaMinus"))->GetXaxis()->SetBinLabel(3, "Passed mass sel"); + rSelections.get(HIST("hSelectionsSigmaMinus"))->GetXaxis()->SetBinLabel(4, "Passed cutDCAtoPvSigma"); + rSelections.get(HIST("hSelectionsSigmaMinus"))->GetXaxis()->SetBinLabel(5, "Passed cutDCAtoPvPiFromSigma"); + rSelections.get(HIST("hSelectionsSigmaMinus"))->GetXaxis()->SetBinLabel(6, "Passed cutSigmaRadius"); + rSelections.get(HIST("hSelectionsSigmaMinus"))->GetXaxis()->SetBinLabel(7, "Passed cutSigmaQtAP"); + rSelections.get(HIST("hSelectionsSigmaMinus"))->GetXaxis()->SetBinLabel(8, "Passed cutSigmaAlphaAP"); + rSelections.add("hSelectionsBachPi", "hSelectionsBachPi", {HistType::kTH1D, {{6, -0.f, 5.5f}}}); + rSelections.get(HIST("hSelectionsBachPi"))->GetXaxis()->SetBinLabel(1, "All"); + rSelections.get(HIST("hSelectionsBachPi"))->GetXaxis()->SetBinLabel(2, "Sign sel"); + rSelections.get(HIST("hSelectionsBachPi"))->GetXaxis()->SetBinLabel(3, "Nsigma Tpc sel"); + rSelections.get(HIST("hSelectionsBachPi"))->GetXaxis()->SetBinLabel(4, "Tpc clusters sel"); + rSelections.get(HIST("hSelectionsBachPi"))->GetXaxis()->SetBinLabel(5, "#eta sel"); + rSelections.get(HIST("hSelectionsBachPi"))->GetXaxis()->SetBinLabel(6, "PID sel"); + rSelections.add("hSelectionsKinkPi", "hSelectionsKinkPi", {HistType::kTH1D, {{7, -0.f, 6.5f}}}); + rSelections.get(HIST("hSelectionsKinkPi"))->GetXaxis()->SetBinLabel(1, "All"); + rSelections.get(HIST("hSelectionsKinkPi"))->GetXaxis()->SetBinLabel(2, "Tpc N#sigma PID sel"); + rSelections.get(HIST("hSelectionsKinkPi"))->GetXaxis()->SetBinLabel(3, "Tpc N clusters sel"); + rSelections.get(HIST("hSelectionsKinkPi"))->GetXaxis()->SetBinLabel(4, "#eta sel"); + rSelections.get(HIST("hSelectionsKinkPi"))->GetXaxis()->SetBinLabel(5, "ITS clusters sel"); + rSelections.get(HIST("hSelectionsKinkPi"))->GetXaxis()->SetBinLabel(6, "has Tof sel"); + rSelections.get(HIST("hSelectionsKinkPi"))->GetXaxis()->SetBinLabel(7, "N#sigma Tof sel"); + rSelections.add("hSelectionsKinkPr", "hSelectionsKinkPr", {HistType::kTH1D, {{7, -0.f, 6.5f}}}); + rSelections.get(HIST("hSelectionsKinkPr"))->GetXaxis()->SetBinLabel(1, "All"); + rSelections.get(HIST("hSelectionsKinkPr"))->GetXaxis()->SetBinLabel(2, "Tpc N#sigma PID sel"); + rSelections.get(HIST("hSelectionsKinkPr"))->GetXaxis()->SetBinLabel(3, "Tpc N clusters sel"); + rSelections.get(HIST("hSelectionsKinkPr"))->GetXaxis()->SetBinLabel(4, "#eta sel"); + rSelections.get(HIST("hSelectionsKinkPr"))->GetXaxis()->SetBinLabel(5, "ITS clusters sel"); + rSelections.get(HIST("hSelectionsKinkPr"))->GetXaxis()->SetBinLabel(6, "has Tof sel"); + rSelections.get(HIST("hSelectionsKinkPr"))->GetXaxis()->SetBinLabel(7, "N#sigma Tof sel"); + + if (doprocessDataWCentQVecs) { + rLambda1405.add("hScalarProd", "hScalarProd", {HistType::kTH1D, {scalarProdAxis}}); + std::vector axesFlow = {lambda1405MassAxis, ptAxis, centAxis, scalarProdAxis}; + rLambda1405.add("hSparseFlowL1405", "THn for SP", HistType::kTHnSparseF, axesFlow); + } if (doprocessMC) { // Add MC histograms if needed, to sigmaminus - rLambda1405.add("h2MassResolution_0", "h2MassResolution_0", {HistType::kTH2F, {massAxis, massResolutionAxis}}); - rLambda1405.add("h2PtResolution_0", "h2PtResolution_0", {HistType::kTH2F, {ptAxis, ptResolutionAxis}}); - rLambda1405.add("h2PtMassMC_0", "h2PtMassMC_0", {HistType::kTH2F, {ptAxis, massAxis}}); + rLambda1405.add("hRecoL1405", "hRecoL1405;;Counts", {HistType::kTH2F, {{6, -0.5, 5.5}, ptAxis}}); + rLambda1405.get(HIST("hRecoL1405"))->GetXaxis()->SetBinLabel(1, "Reconstructed #Lambda(1405)"); + rLambda1405.get(HIST("hRecoL1405"))->GetXaxis()->SetBinLabel(2, "Has MC particle"); + rLambda1405.get(HIST("hRecoL1405"))->GetXaxis()->SetBinLabel(3, "Has #Sigma daug"); + rLambda1405.get(HIST("hRecoL1405"))->GetXaxis()->SetBinLabel(4, "Has bach #pi"); + rLambda1405.get(HIST("hRecoL1405"))->GetXaxis()->SetBinLabel(5, "Has mothers"); + rLambda1405.get(HIST("hRecoL1405"))->GetXaxis()->SetBinLabel(6, "Has same mother"); + rLambda1405.add("hGenL1405", "hGenL1405;;Counts", {HistType::kTH2F, {{3, -0.5, 2.5}, ptAxis}}); + rLambda1405.get(HIST("hGenL1405"))->GetXaxis()->SetBinLabel(1, "#Lambda(1405) #rightarrow #Sigma^{-} #pi^{+} #rightarrow n #pi^{-} #pi^{+}"); + rLambda1405.get(HIST("hGenL1405"))->GetXaxis()->SetBinLabel(2, "#Lambda(1405) #rightarrow #Sigma^{+} #pi^{-} #rightarrow n #pi^{+} #pi^{-}"); + rLambda1405.get(HIST("hGenL1405"))->GetXaxis()->SetBinLabel(3, "#Lambda(1405) #rightarrow #Sigma^{+} #pi^{-} #rightarrow p #pi^{0} #pi^{-}"); + rLambda1405.add("h2MassResolutionFromSigmaMinus", "h2MassResolutionFromSigmaMinus", {HistType::kTH2F, {lambda1405MassAxis, massResolutionAxis}}); + rLambda1405.add("h2PtResolutionFromSigmaMinus", "h2PtResolutionFromSigmaMinus", {HistType::kTH2F, {ptAxis, ptResolutionAxis}}); // Add MC histograms if needed, to sigmaplus - rLambda1405.add("h2MassResolution_1", "h2MassResolution_1", {HistType::kTH2F, {massAxis, massResolutionAxis}}); - rLambda1405.add("h2PtResolution_1", "h2PtResolution_1", {HistType::kTH2F, {ptAxis, ptResolutionAxis}}); - rLambda1405.add("h2PtMassMC_1", "h2PtMassMC_1", {HistType::kTH2F, {ptAxis, massAxis}}); + rLambda1405.add("h2MassResolutionFromSigmaPlus", "h2MassResolutionFromSigmaPlus", {HistType::kTH2F, {lambda1405MassAxis, massResolutionAxis}}); + rLambda1405.add("h2PtResolutionFromSigmaPlus", "h2PtResolutionFromSigmaPlus", {HistType::kTH2F, {ptAxis, ptResolutionAxis}}); + // rLambda1405.add("h2PtMassMCWithSigmaDaug", "h2PtMassMCWithSigmaDaug", {HistType::kTH2F, {ptAxis, lambda1405MassAxis}}); + // rLambda1405.add("h2PtMassMCNoSigmaDaug", "h2PtMassMCNoSigmaDaug", {HistType::kTH2F, {ptAxis, lambda1405MassAxis}}); } } @@ -164,139 +332,330 @@ struct lambda1405analysis { return std::sqrt(p2A - dp * dp / p2V0); } - template - bool selectPiTrack(const Ttrack& candidate, bool piFromSigma) + template + bool selectPiBach(const TTrack& candidate) { - if (std::abs(candidate.tpcNSigmaPi()) > cutNSigmaTPC || candidate.tpcNClsFound() < cutNTPCClusPi || std::abs(candidate.eta()) > cutEtaDaught) { + if (std::abs(candidate.tpcNSigmaPi()) > cutNSigTpc) { return false; } - if (piFromSigma) { - return true; + rSelections.fill(HIST("hSelectionsBachPi"), 2); // Nsigma Tpc + + if (candidate.tpcNClsFound() < cutNTpcClusPi) { + return false; } + rSelections.fill(HIST("hSelectionsBachPi"), 3); // Tpc clusters - if (candidate.itsNCls() < cutNITSClusPi) { + if (std::abs(std::abs(candidate.eta()) > cutEtaDaught)) { return false; } + rSelections.fill(HIST("hSelectionsBachPi"), 4); // Eta selection + + return true; + } - if (useTOF && !candidate.hasTOF()) { + template + bool selectPiKink(const TTrack& candidate) + { + rSelections.fill(HIST("hSelectionsKinkPi"), 0); // All pion kink candidates + + if (std::abs(candidate.tpcNSigmaPi()) > cutNSigTpc) { + return false; + } + rSelections.fill(HIST("hSelectionsKinkPi"), 1); // Nsigma Tpc + + if (candidate.tpcNClsFound() < cutNTpcClusPi) { return false; } + rSelections.fill(HIST("hSelectionsKinkPi"), 2); // Tpc clusters - if (useTOF && std::abs(candidate.tofNSigmaPi()) > cutNSigmaTOF) { + if (std::abs(std::abs(candidate.eta()) > cutEtaDaught)) { return false; } + rSelections.fill(HIST("hSelectionsKinkPi"), 3); // Eta selection + + if (candidate.itsNCls() < cutNITSClusKink) { + return false; + } + rSelections.fill(HIST("hSelectionsKinkPi"), 4); // ITS clusters + + if (useTof && !candidate.hasTOF()) { + return false; + } + rSelections.fill(HIST("hSelectionsKinkPi"), 5); // has Tof + + if (useTof && std::abs(candidate.tofNSigmaPi()) > cutNSigTof) { + return false; + } + rSelections.fill(HIST("hSelectionsKinkPi"), 6); // Nsigma Tof return true; // Track is selected } - template - bool selectProTrack(const Ttrack& candidate, bool prFromSigma) + template + bool selectPrKink(const TTrack& candidate) { - if (std::abs(candidate.tpcNSigmaPr()) > cutNSigmaTPC || candidate.tpcNClsFound() < cutNTPCClusPi || std::abs(candidate.eta()) > cutEtaDaught) { + rSelections.fill(HIST("hSelectionsKinkPr"), 0); // All proton kink candidates + + if (std::abs(candidate.tpcNSigmaPr()) > cutNSigTpc) { + return false; + } + rSelections.fill(HIST("hSelectionsKinkPr"), 1); // Nsigma Tpc + + if (candidate.tpcNClsFound() < cutNTpcClusPi) { return false; } - if (prFromSigma) { - return true; + rSelections.fill(HIST("hSelectionsKinkPr"), 2); // Tpc clusters + + if (std::abs(candidate.eta()) > cutEtaDaught) { + return false; } - if (candidate.itsNCls() < cutNITSClusPi) { + rSelections.fill(HIST("hSelectionsKinkPr"), 3); // eta selection + + if (candidate.itsNCls() < cutNITSClusKink) { return false; } - if (useTOF && !candidate.hasTOF()) { + rSelections.fill(HIST("hSelectionsKinkPr"), 4); // ITS clusters + + if (useTof && !candidate.hasTOF()) { return false; } - if (useTOF && std::abs(candidate.tofNSigmaPr()) > cutNSigmaTOF) { + rSelections.fill(HIST("hSelectionsKinkPr"), 5); // has Tof + + if (useTof && std::abs(candidate.tofNSigmaPr()) > cutNSigTof) { return false; } + rSelections.fill(HIST("hSelectionsKinkPr"), 6); // Nsigma Tof + return true; // Track is selected } - bool selectCandidate(aod::KinkCands::iterator const& sigmaCand, TracksFull const& tracks) + template + void fillHistosSigma(const lambda1405candidate& lambda1405Cand, const TCand& sigmaCand, const TTrack& kinkDauTrack) { + + if (sigmaCand.mothSign() > 0) { + rSigmaPlus.fill(HIST("hSigmaPlusMass"), sigmaCand.mSigmaPlus()); + rSigmaPlus.fill(HIST("h2dPtMassSigmaPlus"), sigmaCand.ptMoth(), sigmaCand.mSigmaPlus()); + rSigmaPlus.fill(HIST("hMassXiMinusSigmaPlus"), sigmaCand.mXiMinus()); + rSigmaPlus.fill(HIST("hSigmaPlusArmPod"), lambda1405Cand.sigmaAlphaAP, lambda1405Cand.sigmaQtAP); + rSigmaPlus.fill(HIST("hSigmaPlusPt"), sigmaCand.ptMoth()); + rSigmaPlus.fill(HIST("hSigmaPlusRadius"), lambda1405Cand.sigmaRadius); + rSigmaPlus.fill(HIST("hSigmaPlusDcaToPv"), sigmaCand.dcaMothPv()); + rSigmaPlus.fill(HIST("hSigmaPlusDcaKinkDauToPv"), sigmaCand.dcaDaugPv()); + // Fill QA histos for kink daughter + rSigmaPlus.fill(HIST("hSigmaPlusKinkPt"), kinkDauTrack.pt()); + rSigmaPlus.fill(HIST("hSigmaPlusKinkTpcNSigPi"), kinkDauTrack.tpcNSigmaPi()); + rSigmaPlus.fill(HIST("hSigmaPlusKinkTofNSigPi"), kinkDauTrack.tofNSigmaPi()); + rSigmaPlus.fill(HIST("hSigmaPlusKinkTpcNSigPr"), kinkDauTrack.tpcNSigmaPr()); + rSigmaPlus.fill(HIST("hSigmaPlusKinkTofNSigPr"), kinkDauTrack.tofNSigmaPr()); + } else { + rSigmaMinus.fill(HIST("hSigmaMinusMass"), sigmaCand.mSigmaMinus()); + rSigmaMinus.fill(HIST("h2dPtMassSigmaMinus"), sigmaCand.ptMoth(), sigmaCand.mSigmaMinus()); + rSigmaMinus.fill(HIST("hMassXiMinusSigmaMinus"), sigmaCand.mXiMinus()); + rSigmaMinus.fill(HIST("hSigmaMinusArmPod"), lambda1405Cand.sigmaAlphaAP, lambda1405Cand.sigmaQtAP); + rSigmaMinus.fill(HIST("hSigmaMinusPt"), sigmaCand.ptMoth()); + rSigmaMinus.fill(HIST("hSigmaMinusRadius"), lambda1405Cand.sigmaRadius); + rSigmaMinus.fill(HIST("hSigmaMinusDcaToPv"), sigmaCand.dcaMothPv()); + rSigmaMinus.fill(HIST("hSigmaMinusDcaKinkDauToPv"), sigmaCand.dcaDaugPv()); + // Fill QA histos for kink daughter + rSigmaMinus.fill(HIST("hSigmaMinusKinkPt"), kinkDauTrack.pt()); + rSigmaMinus.fill(HIST("hSigmaMinusKinkTpcNSigPi"), kinkDauTrack.tpcNSigmaPi()); + rSigmaMinus.fill(HIST("hSigmaMinusKinkTofNSigPi"), kinkDauTrack.tofNSigmaPi()); + } + } + + template + void fillHistosLambda1405(const lambda1405candidate& cand, const TTrack& piTrack) { + + // Fill QA histos for Lambda(1405) candidate + rLambda1405.fill(HIST("hMassL1405"), cand.massL1405); + rLambda1405.fill(HIST("hMassXi1530"), cand.massXi1530); + rLambda1405.fill(HIST("hPx"), cand.px); + rLambda1405.fill(HIST("hPy"), cand.py); + rLambda1405.fill(HIST("hPz"), cand.pz); + rLambda1405.fill(HIST("hPt"), cand.pt()); + rLambda1405.fill(HIST("hPhi"), cand.phi); + + // Bachelor Pi + rLambda1405.fill(HIST("hBachPiPt"), piTrack.pt() * (cand.isSigmaPlus ? -1 : 1)); // Invert pt for Sigma+ to have the correct charge correlation + rLambda1405.fill(HIST("hBachPiNSigTpc"), piTrack.tpcNSigmaPi()); + rLambda1405.fill(HIST("h2BachPiPtNSigTpc"), piTrack.pt(), piTrack.tpcNSigmaPi()); + rLambda1405.fill(HIST("hBachPiNSigTof"), piTrack.tofNSigmaPi()); + rLambda1405.fill(HIST("h2BachPiPtNSigTof"), piTrack.pt(), piTrack.tofNSigmaPi()); + } + + void constructCollCandidates(aod::KinkCands::iterator const& sigmaCand, TracksFull const& tracks, std::vector& selectedCandidates) + { + rSelections.fill(HIST("hSelectionsL1405"), 0); // All candidates + + if (sigmaCand.mothSign() < 0) { + rSelections.fill(HIST("hSelectionsSigmaMinus"), 0); // All Sigma- candidates + } else { + rSelections.fill(HIST("hSelectionsSigmaPlus"), 0); // All Sigma+ candidates + } + auto kinkDauTrack = sigmaCand.trackDaug_as(); - bool isPiKink = selectPiTrack(kinkDauTrack, true); - bool isProKink = selectProTrack(kinkDauTrack, true); - if (!isPiKink && !isProKink) { - return false; + bool isPiKink = selectPiKink(kinkDauTrack); + bool isPrKink = selectPrKink(kinkDauTrack); + if (!isPiKink && !isPrKink) { + return; + } + lambda1405Cand.hasPiKink = isPiKink; + lambda1405Cand.hasPrKink = isPrKink; + + if (sigmaCand.mothSign() < 0) { + rSelections.fill(HIST("hSelectionsSigmaMinus"), 1); // Passed kink sel + } else { + rSelections.fill(HIST("hSelectionsSigmaPlus"), 1); // Passed kink sel } - if (isPiKink) { - rLambda1405.fill(HIST("h2PtMassSigmaBeforeCuts_0"), sigmaCand.mothSign() * sigmaCand.ptMoth(), sigmaCand.mSigmaMinus()); - rLambda1405.fill(HIST("h2PtPiNSigma_0"), sigmaCand.mothSign() * kinkDauTrack.pt(), kinkDauTrack.tpcNSigmaPi()); + // Sigma- or AntiSigma+ candidates + if (isPiKink && sigmaCand.mothSign() < 0) { + rSigmaMinus.fill(HIST("h2PtMassSigmaMinusBeforeCuts"), sigmaCand.mothSign() * sigmaCand.ptMoth(), sigmaCand.mSigmaMinus()); + rSigmaMinus.fill(HIST("h2PtPiKinkNSigBeforeCutsSigmaMinus"), sigmaCand.mothSign() * kinkDauTrack.pt(), kinkDauTrack.tpcNSigmaPi()); } - if (isProKink) { - rLambda1405.fill(HIST("h2PtMassSigmaBeforeCuts_1"), sigmaCand.mothSign() * sigmaCand.ptMoth(), sigmaCand.mSigmaPlus()); - rLambda1405.fill(HIST("h2PtPiNSigma_1"), sigmaCand.mothSign() * kinkDauTrack.pt(), kinkDauTrack.tpcNSigmaPr()); + if (isPiKink && sigmaCand.mothSign() > 0) { + rSigmaPlus.fill(HIST("h2PtMassSigmaPlusBeforeCuts"), sigmaCand.mothSign() * sigmaCand.ptMoth(), sigmaCand.mSigmaPlus()); + rSigmaPlus.fill(HIST("h2PtPiKinkNSigBeforeCutsSigmaPlus"), sigmaCand.mothSign() * kinkDauTrack.pt(), kinkDauTrack.tpcNSigmaPi()); + } + // Only Sigma+ can have a proton as kink daughter + if (isPrKink) { + rSigmaPlus.fill(HIST("h2PtMassSigmaPlusBeforeCuts"), sigmaCand.mothSign() * sigmaCand.ptMoth(), sigmaCand.mSigmaPlus()); + rSigmaPlus.fill(HIST("h2PtPrKinkNSigBeforeCutsSigmaPlus"), sigmaCand.mothSign() * kinkDauTrack.pt(), kinkDauTrack.tpcNSigmaPr()); } - lambda1405Cand.isSigmaPlus = isProKink && (sigmaCand.mSigmaPlus() > o2::constants::physics::MassSigmaPlus - cutSigmaMass && sigmaCand.mSigmaPlus() < o2::constants::physics::MassSigmaPlus + cutSigmaMass); + lambda1405Cand.isSigmaPlus = isPrKink && (sigmaCand.mSigmaPlus() > o2::constants::physics::MassSigmaPlus - cutSigmaMass && sigmaCand.mSigmaPlus() < o2::constants::physics::MassSigmaPlus + cutSigmaMass); lambda1405Cand.isSigmaMinus = isPiKink && (sigmaCand.mSigmaMinus() > o2::constants::physics::MassSigmaMinus - cutSigmaMass && sigmaCand.mSigmaMinus() < o2::constants::physics::MassSigmaMinus + cutSigmaMass); if (!lambda1405Cand.isSigmaPlus && !lambda1405Cand.isSigmaMinus) { - return false; + return; + } + if (sigmaCand.mothSign() < 0) { + rSelections.fill(HIST("hSelectionsSigmaMinus"), 2); // Passed mass sel + } else { + rSelections.fill(HIST("hSelectionsSigmaPlus"), 2); // Passed mass sel } + float sigmaRad = std::hypot(sigmaCand.xDecVtx(), sigmaCand.yDecVtx()); - if (std::abs(sigmaCand.dcaMothPv()) > cutDCAtoPVSigma || std::abs(sigmaCand.dcaDaugPv()) < cutDCAtoPVPiFromSigma || sigmaRad < cutSigmaRadius) { - return false; + if (std::abs(sigmaCand.dcaMothPv()) > cutDCAtoPvSigma) { + return; + } + if (sigmaCand.mothSign() < 0) { + rSelections.fill(HIST("hSelectionsSigmaMinus"), 3); // Passed cutDCAtoPvSigma + } else { + rSelections.fill(HIST("hSelectionsSigmaPlus"), 3); // Passed cutDCAtoPvSigma + } + + if (std::abs(sigmaCand.dcaDaugPv()) < cutDCAtoPvPiFromSigma) { + return; + } + if (sigmaCand.mothSign() < 0) { + rSelections.fill(HIST("hSelectionsSigmaMinus"), 4); // cutDCAtoPvPiFromSigma + } else { + rSelections.fill(HIST("hSelectionsSigmaPlus"), 4); // cutDCAtoPvPiFromSigma + } + + if (sigmaRad < cutSigmaRadius) { + return; + } + if (sigmaCand.mothSign() < 0) { + rSelections.fill(HIST("hSelectionsSigmaMinus"), 5); // Passed mass sel + } else { + rSelections.fill(HIST("hSelectionsSigmaPlus"), 5); // Passed mass sel + } + + auto kinkDauMom = std::array{sigmaCand.pxDaug(), sigmaCand.pyDaug(), sigmaCand.pzDaug()}; + auto sigmaMom = std::array{sigmaCand.pxMoth(), sigmaCand.pyMoth(), sigmaCand.pzMoth()}; + // Sigma properties + lambda1405Cand.sigmaId = sigmaCand.globalIndex(); + lambda1405Cand.sigmaMinusMass = sigmaCand.mSigmaMinus(); + lambda1405Cand.sigmaPlusMass = sigmaCand.mSigmaPlus(); + lambda1405Cand.xiMinusMass = sigmaCand.mXiMinus(); + lambda1405Cand.sigmaSign = sigmaCand.mothSign(); + lambda1405Cand.sigmaAlphaAP = alphaAP(sigmaMom, kinkDauMom); + lambda1405Cand.sigmaQtAP = qtAP(sigmaMom, kinkDauMom); + lambda1405Cand.sigmaPt = sigmaCand.ptMoth(); + lambda1405Cand.sigmaRadius = sigmaRad; + lambda1405Cand.kinkDcaDauToPv = sigmaCand.dcaDaugPv(); + + if (lambda1405Cand.sigmaQtAP < cutSigmaQtAPMin || lambda1405Cand.sigmaQtAP > cutSigmaQtAPMax) { + return; + } + if (sigmaCand.mothSign() < 0) { + rSelections.fill(HIST("hSelectionsSigmaMinus"), 6); // Passed mass sel + } else { + rSelections.fill(HIST("hSelectionsSigmaPlus"), 6); // Passed mass sel } + if (lambda1405Cand.sigmaAlphaAP < cutSigmaAlphaAPMin || lambda1405Cand.sigmaAlphaAP > cutSigmaAlphaAPMax) { + return; + } + if (sigmaCand.mothSign() < 0) { + rSelections.fill(HIST("hSelectionsSigmaMinus"), 7); // Passed mass sel + } else { + rSelections.fill(HIST("hSelectionsSigmaPlus"), 7); // Passed mass sel + } + + // Kink daughter properties + lambda1405Cand.kinkDauId = kinkDauTrack.globalIndex(); + lambda1405Cand.kinkPt = kinkDauTrack.pt(); + lambda1405Cand.kinkPiNSigTpc = kinkDauTrack.tpcNSigmaPi(); + lambda1405Cand.kinkPiNSigTof = kinkDauTrack.tofNSigmaPi(); + lambda1405Cand.kinkPrNSigTpc = kinkDauTrack.tpcNSigmaPr(); + lambda1405Cand.kinkPrNSigTof = kinkDauTrack.tofNSigmaPr(); + + fillHistosSigma(lambda1405Cand, sigmaCand, kinkDauTrack); + rSelections.fill(HIST("hSelectionsL1405"), 1); // Passed Sigma sel + for (const auto& piTrack : tracks) { - if (!doLSBkg) { - if (piTrack.sign() == sigmaCand.mothSign()) { - continue; - } - } else { - if (piTrack.sign() != sigmaCand.mothSign()) { - continue; - } - } + rSelections.fill(HIST("hSelectionsBachPi"), 0); // All bachelors - if (!selectPiTrack(piTrack, false)) { + bool isUnlikeSign = (piTrack.sign() != sigmaCand.mothSign()); + bool acceptPair = doLSBkg ? !isUnlikeSign : isUnlikeSign; + if (!acceptPair) { continue; } + rSelections.fill(HIST("hSelectionsBachPi"), 1); - auto kinkDauMom = std::array{sigmaCand.pxDaug(), sigmaCand.pyDaug(), sigmaCand.pzDaug()}; - auto sigmaMom = std::array{sigmaCand.pxMoth(), sigmaCand.pyMoth(), sigmaCand.pzMoth()}; + if (!selectPiBach(piTrack)) { + continue; + } + rSelections.fill(HIST("hSelectionsBachPi"), 5); // PID sel + rSelections.fill(HIST("hSelectionsL1405"), 2); // Bach Pi selection + auto piMom = std::array{piTrack.px(), piTrack.py(), piTrack.pz()}; - float invMass = RecoDecay::m(std::array{sigmaMom, piMom}, std::array{o2::constants::physics::MassSigmaMinus, o2::constants::physics::MassPiPlus}); - float invMassXiPi = RecoDecay::m(std::array{sigmaMom, kinkDauMom}, std::array{o2::constants::physics::MassXiMinus, o2::constants::physics::MassPiPlus}); + float invMass{-1.f}; + if (lambda1405Cand.isSigmaMinus) { + invMass = RecoDecay::m(std::array{sigmaMom, piMom}, std::array{o2::constants::physics::MassSigmaMinus, o2::constants::physics::MassPiPlus}); + } else if (lambda1405Cand.isSigmaPlus) { + invMass = RecoDecay::m(std::array{sigmaMom, piMom}, std::array{o2::constants::physics::MassSigmaPlus, o2::constants::physics::MassPiMinus}); + } if (invMass > cutUpperMass) { continue; } + rSelections.fill(HIST("hSelectionsL1405"), 3); // Upper mass selection - lambda1405Cand.kinkDauID = kinkDauTrack.globalIndex(); - lambda1405Cand.sigmaID = sigmaCand.globalIndex(); - lambda1405Cand.piID = piTrack.globalIndex(); + // Daughter Pi properties + lambda1405Cand.piId = piTrack.globalIndex(); + lambda1405Cand.piPt = piTrack.pt(); + lambda1405Cand.bachPiNSigTpc = piTrack.tpcNSigmaPi(); + if (useTof) { + lambda1405Cand.bachPiNSigTof = piTrack.tofNSigmaPi(); + } else { + lambda1405Cand.bachPiNSigTof = -999; // Not used if Tof is not enabled + } + // Lambda(1405) candidate properties + lambda1405Cand.massL1405 = invMass; + lambda1405Cand.massXi1530 = RecoDecay::m(std::array{sigmaMom, piMom}, std::array{o2::constants::physics::MassXiMinus, o2::constants::physics::MassPiPlus}); lambda1405Cand.px = sigmaMom[0] + piMom[0]; lambda1405Cand.py = sigmaMom[1] + piMom[1]; lambda1405Cand.pz = sigmaMom[2] + piMom[2]; - lambda1405Cand.mass = invMass; - lambda1405Cand.massXi1530 = invMassXiPi; - - lambda1405Cand.sigmaMinusMass = sigmaCand.mSigmaMinus(); - lambda1405Cand.sigmaPlusMass = sigmaCand.mSigmaPlus(); - lambda1405Cand.xiMinusMass = sigmaCand.mXiMinus(); - lambda1405Cand.sigmaSign = sigmaCand.mothSign(); - lambda1405Cand.sigmaAlphaAP = alphaAP(sigmaMom, kinkDauMom); - lambda1405Cand.sigmaQtAP = qtAP(sigmaMom, kinkDauMom); - lambda1405Cand.sigmaPt = sigmaCand.ptMoth(); - lambda1405Cand.sigmaRadius = sigmaRad; - lambda1405Cand.kinkPt = kinkDauTrack.pt(); - lambda1405Cand.kinkTPCNSigmaPi = kinkDauTrack.tpcNSigmaPi(); - lambda1405Cand.kinkTOFNSigmaPi = kinkDauTrack.tofNSigmaPi(); - lambda1405Cand.kinkTPCNSigmaPr = kinkDauTrack.tpcNSigmaPr(); - lambda1405Cand.kinkTOFNSigmaPr = kinkDauTrack.tofNSigmaPr(); - lambda1405Cand.dcaKinkDauToPV = sigmaCand.dcaDaugPv(); - - lambda1405Cand.piPt = piTrack.pt(); - lambda1405Cand.nSigmaTPCPi = piTrack.tpcNSigmaPi(); - if (useTOF) { - lambda1405Cand.nSigmaTOFPi = piTrack.tofNSigmaPi(); - } else { - lambda1405Cand.nSigmaTOFPi = -999; // Not used if TOF is not enabled - } - return true; // Candidate is selected + lambda1405Cand.phi = std::atan2(lambda1405Cand.py, lambda1405Cand.px); + lambda1405Cand.scalarProd = -1; + fillHistosLambda1405(lambda1405Cand, piTrack); + rSelections.fill(HIST("hSelectionsL1405"), 4); // Accepted + selectedCandidates.push_back(lambda1405Cand); } - return false; // No valid pion track found } template @@ -326,53 +685,182 @@ struct lambda1405analysis { } rEventSelection.fill(HIST("hVertexZRec"), collision.posZ()); for (const auto& sigmaCand : kinkCands) { - if (selectCandidate(sigmaCand, tracks)) { + std::vector selectedCandidates; + constructCollCandidates(sigmaCand, tracks, selectedCandidates); + for (const auto& lambda1405Cand : selectedCandidates) { if (lambda1405Cand.isSigmaMinus) { - rLambda1405.fill(HIST("h2PtMass_0"), lambda1405Cand.sigmaSign * lambda1405Cand.pt(), lambda1405Cand.mass); - rLambda1405.fill(HIST("h2PtMassSigma_0"), lambda1405Cand.sigmaSign * lambda1405Cand.sigmaPt, lambda1405Cand.sigmaMinusMass); - rLambda1405.fill(HIST("h2SigmaMassVsMass_0"), lambda1405Cand.mass, lambda1405Cand.sigmaMinusMass); - rLambda1405.fill(HIST("h2PtPiNSigmaTOF_0"), lambda1405Cand.sigmaSign * lambda1405Cand.piPt, lambda1405Cand.nSigmaTOFPi); + rSigmaMinus.fill(HIST("h2SigmaMinusMassVsLambdaMass"), lambda1405Cand.sigmaSign * lambda1405Cand.sigmaPt, lambda1405Cand.sigmaMinusMass); + if (lambda1405Cand.hasPiKink) { + rSigmaMinus.fill(HIST("h2KinkPiPtNSigTofSigmaMinus"), lambda1405Cand.sigmaSign * lambda1405Cand.piPt, lambda1405Cand.bachPiNSigTof); + } + if (lambda1405Cand.hasPrKink) { + rSigmaMinus.fill(HIST("h2KinkPrPtNSigTofSigmaMinus"), lambda1405Cand.sigmaSign * lambda1405Cand.piPt, lambda1405Cand.bachPiNSigTof); + } } if (lambda1405Cand.isSigmaPlus) { - rLambda1405.fill(HIST("h2PtMass_1"), lambda1405Cand.sigmaSign * lambda1405Cand.pt(), lambda1405Cand.mass); - rLambda1405.fill(HIST("h2PtMassSigma_1"), lambda1405Cand.sigmaSign * lambda1405Cand.sigmaPt, lambda1405Cand.sigmaPlusMass); - rLambda1405.fill(HIST("h2SigmaMassVsMass_1"), lambda1405Cand.mass, lambda1405Cand.sigmaPlusMass); - rLambda1405.fill(HIST("h2PtPiNSigmaTOF_1"), lambda1405Cand.sigmaSign * lambda1405Cand.piPt, lambda1405Cand.nSigmaTOFPi); + rSigmaPlus.fill(HIST("h2SigmaPlusMassVsLambdaMass"), lambda1405Cand.sigmaSign * lambda1405Cand.sigmaPt, lambda1405Cand.sigmaPlusMass); + rSigmaPlus.fill(HIST("h2KinkPrPtNSigTofSigmaPlus"), lambda1405Cand.sigmaSign * lambda1405Cand.piPt, lambda1405Cand.bachPiNSigTof); } if (fillOutputTree) { + float const ptCand = lambda1405Cand.pt(); + if (downSampleFactor < 1.) { + float const pseudoRndm = ptCand * 1000. - static_cast(ptCand * 1000); + if (ptCand < ptDownSampleMax && pseudoRndm >= downSampleFactor) { + continue; + } + } outputDataTable(lambda1405Cand.px, lambda1405Cand.py, lambda1405Cand.pz, - lambda1405Cand.mass, lambda1405Cand.massXi1530, + lambda1405Cand.massL1405, lambda1405Cand.massXi1530, lambda1405Cand.sigmaMinusMass, lambda1405Cand.sigmaPlusMass, lambda1405Cand.xiMinusMass, lambda1405Cand.sigmaPt, lambda1405Cand.sigmaAlphaAP, lambda1405Cand.sigmaQtAP, lambda1405Cand.sigmaRadius, lambda1405Cand.kinkPt, - lambda1405Cand.kinkTPCNSigmaPi, lambda1405Cand.kinkTOFNSigmaPi, - lambda1405Cand.kinkTPCNSigmaPr, lambda1405Cand.kinkTOFNSigmaPr, - lambda1405Cand.dcaKinkDauToPV, - lambda1405Cand.nSigmaTPCPi, lambda1405Cand.nSigmaTOFPi); + lambda1405Cand.kinkPiNSigTpc, lambda1405Cand.kinkPiNSigTof, + lambda1405Cand.kinkPrNSigTpc, lambda1405Cand.kinkPrNSigTof, + lambda1405Cand.kinkDcaDauToPv, + lambda1405Cand.bachPiNSigTpc, lambda1405Cand.bachPiNSigTof); } } } } PROCESS_SWITCH(lambda1405analysis, processData, "Data processing", true); + void processDataWCentQVecs(CollisionsFullWCentQVecs::iterator const& collision, aod::KinkCands const& kinkCands, TracksFull const& tracks) + { + LOG(info) << "Processing collision with centrality " << collision.centFT0C() << " and Q vector (" << collision.qvecFT0CRe() << ", " << collision.qvecFT0CIm() << ")"; + if (collision.centFT0C() < centralityMin || collision.centFT0C() > centralityMax) { + LOG(info) << "Skipping collision due to centrality cut"; + return; + } + if (std::abs(collision.posZ()) > cutzvertex || !collision.sel8()) { + LOG(info) << "Skipping collision due to vertex cut"; + return; + } + rEventSelection.fill(HIST("hVertexZRec"), collision.posZ()); + for (const auto& sigmaCand : kinkCands) { + LOG(info) << "Processing Sigma candidate with mother sign " << sigmaCand.mothSign() << " and pt " << sigmaCand.ptMoth(); + std::vector selectedCandidates; + constructCollCandidates(sigmaCand, tracks, selectedCandidates); + for (auto& lambda1405Cand : selectedCandidates) { + if (lambda1405Cand.isSigmaMinus) { + rSigmaMinus.fill(HIST("h2SigmaMinusMassVsLambdaMass"), lambda1405Cand.sigmaSign * lambda1405Cand.sigmaPt, lambda1405Cand.sigmaMinusMass); + if (lambda1405Cand.hasPiKink) { + rSigmaMinus.fill(HIST("h2KinkPiPtNSigTofSigmaMinus"), lambda1405Cand.sigmaSign * lambda1405Cand.piPt, lambda1405Cand.bachPiNSigTof); + } + if (lambda1405Cand.hasPrKink) { + rSigmaMinus.fill(HIST("h2KinkPrPtNSigTofSigmaMinus"), lambda1405Cand.sigmaSign * lambda1405Cand.piPt, lambda1405Cand.bachPiNSigTof); + } + } + if (lambda1405Cand.isSigmaPlus) { + rSigmaPlus.fill(HIST("h2SigmaPlusMassVsLambdaMass"), lambda1405Cand.sigmaSign * lambda1405Cand.sigmaPt, lambda1405Cand.sigmaPlusMass); + rSigmaPlus.fill(HIST("h2KinkPrPtNSigTofSigmaPlus"), lambda1405Cand.sigmaSign * lambda1405Cand.piPt, lambda1405Cand.bachPiNSigTof); + } + LOG(info) << "Filled Sigma histograms for candidate with Sigma pt " << lambda1405Cand.sigmaPt << " and mass " << (lambda1405Cand.isSigmaMinus ? lambda1405Cand.sigmaMinusMass : lambda1405Cand.sigmaPlusMass); + float const xQVec = collision.qvecFT0CRe(); + float const yQVec = collision.qvecFT0CIm(); + float const cos2Phi = std::cos(2 * lambda1405Cand.phi); + float const sin2Phi = std::sin(2 * lambda1405Cand.phi); + lambda1405Cand.scalarProd = cos2Phi * xQVec + sin2Phi * yQVec; + float const ptCand = lambda1405Cand.pt(); + rLambda1405.fill(HIST("hSparseFlowL1405"), ptCand, lambda1405Cand.massL1405, collision.centFT0C(), lambda1405Cand.scalarProd); + LOG(info) << "Filled flow sparse for candidate with pt " << ptCand << ", mass " << lambda1405Cand.massL1405 << ", centrality " << collision.centFT0C() << " and scalar product " << lambda1405Cand.scalarProd; + if (fillFlowTree) { + if (downSampleFactor < 1.) { + float const pseudoRndm = ptCand * 1000. - static_cast(ptCand * 1000); + if (ptCand < ptDownSampleMax && pseudoRndm >= downSampleFactor) { + continue; + } + } + outputDataFlowTable(ptCand, lambda1405Cand.massL1405, + lambda1405Cand.sigmaMinusMass, lambda1405Cand.sigmaPlusMass, + lambda1405Cand.sigmaAlphaAP, lambda1405Cand.sigmaQtAP, + lambda1405Cand.kinkPiNSigTpc, lambda1405Cand.kinkPiNSigTof, + lambda1405Cand.kinkPrNSigTpc, lambda1405Cand.kinkPrNSigTof, + lambda1405Cand.kinkDcaDauToPv, + lambda1405Cand.bachPiNSigTpc, lambda1405Cand.bachPiNSigTof, + lambda1405Cand.scalarProd, collision.centFT0C()); + } + LOG(info) << "Filled flow tree for candidate with pt " << ptCand; + } + } + } + PROCESS_SWITCH(lambda1405analysis, processDataWCentQVecs, "Data processing with centrality and Q vectors info", false); + + template + int matchGenDecay(const TMother& motherPart, const aod::McParticles& mcParticles) { + LOG(info) << "Matching MC decay for particle with PDG code " << motherPart.pdgCode() << " and index " << motherPart.globalIndex(); + int pdgMother = motherPart.pdgCode(); + int8_t sign = 0; + + int MaxDepthFinState = 2; // Maximum depth to look for the decay chain + + // Match L(1405) --> n pi- pi+ final state + std::array finalState; + if (pdgMother > 0) { // Change sign of neutral decay products + finalState = {PDG_t::kNeutron, PDG_t::kPiMinus, PDG_t::kPiPlus}; + } else { + finalState = {-PDG_t::kNeutron, PDG_t::kPiMinus, PDG_t::kPiPlus}; + } + + if (RecoDecay::isMatchedMCGen(mcParticles, motherPart, lambda1405PdgCode, finalState, true, &sign, MaxDepthFinState)) { + // Match the intermediate Sigma + std::vector arrDaughIdxs = {}; + RecoDecay::getDaughters(motherPart, &arrDaughIdxs, std::array{0}, 1); + for (auto iProng = 0u; iProng < arrDaughIdxs.size(); ++iProng) { + auto daughI = mcParticles.rawIteratorAt(arrDaughIdxs[iProng]); + if (std::abs(daughI.pdgCode()) == PDG_t::kSigmaPlus) { + return kSigmaPlusPiToPiPiNeutron; + } + } + return kSigmaMinusPiToPiPiNeutron; + } + + // Match L(1405) --> p pi0 pi+ final state, only possible for Sigma+ + if (pdgMother > 0) { // Change sign of neutral decay products + finalState = {PDG_t::kProton, PDG_t::kPi0, PDG_t::kPiMinus}; + } else { + finalState = {PDG_t::kProton, -PDG_t::kPi0, PDG_t::kPiMinus}; + } + + if (RecoDecay::isMatchedMCGen(mcParticles, motherPart, lambda1405PdgCode, finalState, true, &sign, MaxDepthFinState)) { + // Match the intermediate Sigma + return kSigmaPlusPiToPiPiProton; + } + + return -1; + } + void processMC(CollisionsFullMC const& collisions, aod::KinkCands const& kinkCands, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const& particlesMC, TracksFull const& tracks) { for (const auto& collision : collisions) { - if (std::abs(collision.posZ()) > cutzvertex || !collision.sel8()) { + if (std::abs(collision.posZ()) > cutzvertex) { // || !collision.sel8()) { continue; } rEventSelection.fill(HIST("hVertexZRec"), collision.posZ()); auto sigmaCandsPerCol = kinkCands.sliceBy(mKinkPerCol, collision.globalIndex()); auto tracksPerCol = tracks.sliceBy(mPerColTracks, collision.globalIndex()); + if (sigmaCandsPerCol.size() != 0) { + LOG(info) << "Processing collision with " << sigmaCandsPerCol.size() << " Sigma kink candidates and " << tracksPerCol.size() << " tracks"; + } for (const auto& sigmaCand : sigmaCandsPerCol) { - if (selectCandidate(sigmaCand, tracksPerCol)) { + std::vector selectedCandidates; + constructCollCandidates(sigmaCand, tracksPerCol, selectedCandidates); + LOG(info) << "Selected " << selectedCandidates.size() << " Lambda(1405) candidates in this collision"; + for (const auto& lambda1405Cand : selectedCandidates) { + rLambda1405.fill(HIST("hRecoL1405"), 0., lambda1405Cand.pt()); // All reconstructed + // Do MC association - auto mcLabPiKink = trackLabelsMC.rawIteratorAt(lambda1405Cand.kinkDauID); - auto mcLabSigma = trackLabelsMC.rawIteratorAt(lambda1405Cand.sigmaID); - auto mcLabPi = trackLabelsMC.rawIteratorAt(lambda1405Cand.piID); + auto mcLabPiKink = trackLabelsMC.rawIteratorAt(lambda1405Cand.kinkDauId); + auto mcLabSigma = trackLabelsMC.rawIteratorAt(lambda1405Cand.sigmaId); + auto mcLabPi = trackLabelsMC.rawIteratorAt(lambda1405Cand.piId); if (!mcLabSigma.has_mcParticle() || mcLabPiKink.has_mcParticle() || mcLabPi.has_mcParticle()) { + LOG(info) << "Skipping candidate due to missing MC association: " + << "Sigma MC assoc: " << mcLabSigma.has_mcParticle() << ", " + << "Kink daughter MC assoc: " << mcLabPiKink.has_mcParticle() << ", " + << "Bachelor Pi MC assoc: " << mcLabPi.has_mcParticle(); continue; // Skip if no valid MC association } + rLambda1405.fill(HIST("hRecoL1405"), 1., lambda1405Cand.pt()); // All with associated MC particle + auto mcTrackKink = mcLabPiKink.mcParticle_as(); auto mcTrackSigma = mcLabSigma.mcParticle_as(); auto mcTrackPi = mcLabPi.mcParticle_as(); @@ -382,16 +870,28 @@ struct lambda1405analysis { bool isSigmaPlusToPrKink = checkSigmaKinkMC(mcTrackSigma, mcTrackKink, 3222, 2212, particlesMC); if (!isSigmaMinusKink && !isSigmaPlusToPiKink && !isSigmaPlusToPrKink) { + LOG(info) << "Skipping candidate due to failed MC kink decay check: " + << "isSigmaMinusKink: " << isSigmaMinusKink << ", " + << "isSigmaPlusToPiKink: " << isSigmaPlusToPiKink << ", " + << "isSigmaPlusToPrKink: " << isSigmaPlusToPrKink; continue; // Skip if not a valid Sigma kink decay } + rLambda1405.fill(HIST("hRecoL1405"), 2., lambda1405Cand.pt()); // Has kink decay in MC if (std::abs(mcTrackPi.pdgCode()) != 211) { + LOG(info) << "Skipping candidate due to bachelor Pi not being a pion in MC: " + << "Bachelor Pi PDG code: " << mcTrackPi.pdgCode(); continue; // Skip if not a valid pion candidate } + rLambda1405.fill(HIST("hRecoL1405"), 3., lambda1405Cand.pt()); // Has bach pi if (!mcTrackSigma.has_mothers() || !mcTrackPi.has_mothers()) { + LOG(info) << "Skipping candidate due to missing mothers in MC: " + << "Sigma has mothers: " << mcTrackSigma.has_mothers() << ", " + << "Pi has mothers: " << mcTrackPi.has_mothers(); continue; // Skip if no mothers found } + rLambda1405.fill(HIST("hRecoL1405"), 4., lambda1405Cand.pt()); // Has mothers for Sigma and Pi // check that labpi and labsigma have the same mother (a lambda1405 candidate) int lambda1405Id = -1; @@ -399,42 +899,40 @@ struct lambda1405analysis { for (const auto& sigmaMother : mcTrackSigma.mothers_as()) { if (piMother.globalIndex() == sigmaMother.globalIndex() && std::abs(piMother.pdgCode()) == lambda1405PdgCode) { lambda1405Id = piMother.globalIndex(); + LOG(info) << "Found common mother for Sigma and Pi with global index: " << lambda1405Id; break; // Found the mother, exit loop } } } if (lambda1405Id == -1) { + LOG(info) << "Skipping candidate due to Sigma and pion not sharing the same lambda1405 candidate"; continue; // Skip if the Sigma and pion do not share the same lambda1405 candidate } + rLambda1405.fill(HIST("hRecoL1405"), 4., lambda1405Cand.pt()); // Has same mother + auto lambda1405Mother = particlesMC.rawIteratorAt(lambda1405Id); float lambda1405Mass = std::sqrt(lambda1405Mother.e() * lambda1405Mother.e() - lambda1405Mother.p() * lambda1405Mother.p()); if (lambda1405Cand.isSigmaMinus) { - rLambda1405.fill(HIST("h2PtMass_0"), lambda1405Cand.sigmaSign * lambda1405Cand.pt(), lambda1405Cand.mass); - rLambda1405.fill(HIST("h2PtMassSigma_0"), lambda1405Cand.sigmaSign * lambda1405Cand.sigmaPt, lambda1405Cand.sigmaMinusMass); - rLambda1405.fill(HIST("h2SigmaMassVsMass_0"), lambda1405Cand.mass, lambda1405Cand.sigmaMinusMass); - rLambda1405.fill(HIST("h2PtPiNSigma_0"), lambda1405Cand.piPt, lambda1405Cand.nSigmaTPCPi); - rLambda1405.fill(HIST("h2MassResolution_0"), lambda1405Mass, lambda1405Mass - lambda1405Cand.mass); - rLambda1405.fill(HIST("h2PtResolution_0"), lambda1405Cand.pt(), lambda1405Cand.pt() - lambda1405Mother.pt()); + rSigmaMinus.fill(HIST("h2SigmaMinusMassVsLambdaMass"), lambda1405Cand.massL1405, lambda1405Cand.sigmaMinusMass); + rLambda1405.fill(HIST("h2MassResolutionFromSigmaMinus"), lambda1405Mass, lambda1405Mass - lambda1405Cand.massL1405); + rLambda1405.fill(HIST("h2PtResolutionFromSigmaMinus"), lambda1405Cand.pt(), lambda1405Cand.pt() - lambda1405Mother.pt()); } if (lambda1405Cand.isSigmaPlus) { - rLambda1405.fill(HIST("h2PtMass_1"), lambda1405Cand.sigmaSign * lambda1405Cand.pt(), lambda1405Cand.mass); - rLambda1405.fill(HIST("h2PtMassSigma_1"), lambda1405Cand.sigmaSign * lambda1405Cand.sigmaPt, lambda1405Cand.sigmaPlusMass); - rLambda1405.fill(HIST("h2SigmaMassVsMass_1"), lambda1405Cand.mass, lambda1405Cand.sigmaPlusMass); - rLambda1405.fill(HIST("h2PtPiNSigma_1"), lambda1405Cand.piPt, lambda1405Cand.nSigmaTPCPi); - rLambda1405.fill(HIST("h2MassResolution_1"), lambda1405Mass, lambda1405Mass - lambda1405Cand.mass); - rLambda1405.fill(HIST("h2PtResolution_1"), lambda1405Cand.pt(), lambda1405Cand.pt() - lambda1405Mother.pt()); + rSigmaPlus.fill(HIST("h2SigmaPlusMassVsLambdaMass"), lambda1405Cand.massL1405, lambda1405Cand.sigmaPlusMass); + rLambda1405.fill(HIST("h2MassResolutionFromSigmaPlus"), lambda1405Mass, lambda1405Mass - lambda1405Cand.massL1405); + rLambda1405.fill(HIST("h2PtResolutionFromSigmaPlus"), lambda1405Cand.pt(), lambda1405Cand.pt() - lambda1405Mother.pt()); } if (fillOutputTree) { outputDataTableMC(lambda1405Cand.px, lambda1405Cand.py, lambda1405Cand.pz, - lambda1405Cand.mass, lambda1405Cand.massXi1530, + lambda1405Cand.massL1405, lambda1405Cand.massXi1530, lambda1405Cand.sigmaMinusMass, lambda1405Cand.sigmaPlusMass, lambda1405Cand.xiMinusMass, lambda1405Cand.sigmaPt, lambda1405Cand.sigmaAlphaAP, lambda1405Cand.sigmaQtAP, lambda1405Cand.sigmaRadius, lambda1405Cand.kinkPt, - lambda1405Cand.kinkTPCNSigmaPi, lambda1405Cand.kinkTOFNSigmaPi, - lambda1405Cand.kinkTPCNSigmaPr, lambda1405Cand.kinkTOFNSigmaPr, - lambda1405Cand.dcaKinkDauToPV, - lambda1405Cand.nSigmaTPCPi, lambda1405Cand.nSigmaTOFPi, + lambda1405Cand.kinkPiNSigTpc, lambda1405Cand.kinkPiNSigTof, + lambda1405Cand.kinkPrNSigTpc, lambda1405Cand.kinkPrNSigTof, + lambda1405Cand.kinkDcaDauToPv, + lambda1405Cand.bachPiNSigTpc, lambda1405Cand.bachPiNSigTof, lambda1405Mother.pt(), lambda1405Mass, mcTrackSigma.pdgCode(), mcTrackKink.pdgCode()); } } @@ -446,27 +944,11 @@ struct lambda1405analysis { if (std::abs(mcPart.pdgCode()) != lambda1405PdgCode) { continue; // Only consider Lambda(1405) candidates } - - if (!mcPart.has_daughters()) { - continue; // Skip if no daughters - } - - // Check if the Lambda(1405) has a Sigma daughter - bool hasSigmaDaughter = false; - int dauPdgCode = 0; - for (const auto& daughter : mcPart.daughters_as()) { - if (std::abs(daughter.pdgCode()) == 3122 || std::abs(daughter.pdgCode()) == 3222) { // Sigma PDG code - hasSigmaDaughter = true; - dauPdgCode = daughter.pdgCode(); - break; // Found a Sigma daughter, exit loop - } + int decayChannel = matchGenDecay(mcPart, particlesMC); + if (decayChannel == -1) { + continue; // Skip if it doesn't match the decay channels of interest } - if (!hasSigmaDaughter) { - continue; // Skip if no Sigma daughter found - } - - float mcMass = std::sqrt(mcPart.e() * mcPart.e() - mcPart.p() * mcPart.p()); - dauPdgCode ? rLambda1405.fill(HIST("h2PtMassMC_0"), mcPart.pt(), mcMass) : rLambda1405.fill(HIST("h2PtMassMC_1"), mcPart.pt(), mcMass); + rLambda1405.fill(HIST("hGenL1405"), decayChannel, mcPart.pt()); } } PROCESS_SWITCH(lambda1405analysis, processMC, "MC processing", false); diff --git a/PWGLF/Utils/svPoolCreator.h b/PWGLF/Utils/svPoolCreator.h index 68a621fd977..bcb62b95f61 100644 --- a/PWGLF/Utils/svPoolCreator.h +++ b/PWGLF/Utils/svPoolCreator.h @@ -87,8 +87,8 @@ class svPoolCreator } } - template - void appendTrackCand(const T& trackCand, const C& collisions, int pdgHypo, o2::aod::AmbiguousTracks const& ambiTracks, BC const&) + template + void appendTrackCand(const T& trackCand, const C& collisions, int pdgHypo, const AT& ambiTracks, BC const&) { if (pdgHypo != track0Pdg && pdgHypo != track1Pdg) { LOGP(debug, "Wrong pdg hypothesis"); @@ -106,11 +106,11 @@ class svPoolCreator if (ambTrack.trackId() != trackCand.globalIndex()) { continue; } - if (!ambTrack.has_bc() || ambTrack.bc_as().size() == 0) { + if (!ambTrack.has_bc() || ambTrack.template bc_as().size() == 0) { globalBC = BcInvalid; break; } - globalBC = ambTrack.bc_as().begin().globalBC(); + globalBC = ambTrack.template bc_as().begin().globalBC(); break; } } else { @@ -170,7 +170,6 @@ class svPoolCreator } else if (TESTBIT(trackCand.flags(), o2::aod::track::TrackTimeResIsRange)) { thresholdTime = std::sqrt(sigmaTimeRes2); thresholdTime += timeMarginNS; - } else { thresholdTime = 4. * std::sqrt(sigmaTimeRes2); thresholdTime += timeMarginNS; @@ -197,6 +196,7 @@ class svPoolCreator // is Sorting Needed ? TBD } + template std::vector& getSVCandPool(const C& collisions, bool combineLikeSign = false) { @@ -208,37 +208,34 @@ class svPoolCreator mVtxTrack0[i].clear(); mVtxTrack0[i].resize(collisions.size(), -1); } - - for (int pn = 0; pn < 2; pn++) { - auto& vtxFirstT = mVtxTrack0[pn]; - const auto& signTrack0Pool = track0Pool[pn]; + for (int iCharge = 0; iCharge < 2; iCharge++) { + auto& vtxFirstT = mVtxTrack0[iCharge]; + const auto& signTrack0Pool = track0Pool[iCharge]; for (unsigned i = 0; i < signTrack0Pool.size(); i++) { const auto& track0Seed = signTrack0Pool[i]; - LOG(debug) << "Processsing track0 with index: " << track0Seed.Idxtr << " min bracket: " << track0Seed.collBracket.getMin() << " max bracket: " << track0Seed.collBracket.getMax(); for (int j{track0Seed.collBracket.getMin()}; j <= track0Seed.collBracket.getMax(); ++j) { if (vtxFirstT[j] == -1) { vtxFirstT[j] = i; } } } - int track1sign = combineLikeSign ? pn : 1 - pn; + int track1sign = combineLikeSign ? iCharge : 1 - iCharge; auto& signTrack1 = track1Pool[track1sign]; - for (unsigned itp = 0; itp < signTrack1.size(); itp++) { - auto& track1Seed = signTrack1[itp]; - LOG(debug) << "Processing track1 with index: " << track1Seed.Idxtr << " min bracket: " << track1Seed.collBracket.getMin() << " max bracket: " << track1Seed.collBracket.getMax(); - int firsOverlapIdx = -1; + for (unsigned iTrack1 = 0; iTrack1 < signTrack1.size(); iTrack1++) { + auto& track1Seed = signTrack1[iTrack1]; + int firstOverlapIdx = -1; for (int j{track1Seed.collBracket.getMin()}; j <= track1Seed.collBracket.getMax(); ++j) { LOG(debug) << "Checking vtxFirstT at position " << j << " with value " << vtxFirstT[j]; if (vtxFirstT[j] != -1) { - firsOverlapIdx = vtxFirstT[j]; + firstOverlapIdx = vtxFirstT[j]; break; } } - if (firsOverlapIdx < 0) { + if (firstOverlapIdx < 0) { continue; } - for (unsigned itn = firsOverlapIdx; itn < signTrack0Pool.size(); itn++) { - auto& track0Seed = signTrack0Pool[itn]; + for (unsigned iTrack0 = firstOverlapIdx; iTrack0 < signTrack0Pool.size(); iTrack0++) { + auto& track0Seed = signTrack0Pool[iTrack0]; if (track0Seed.collBracket.getMin() > track1Seed.collBracket.getMax()) { break; From ffc68474b3fbbf7a75cec3bba3754771d9691d7c Mon Sep 17 00:00:00 2001 From: marcellocosti Date: Tue, 2 Jun 2026 05:46:49 +0200 Subject: [PATCH 2/5] Process functions with centrality selection --- PWGLF/TableProducer/Common/kinkBuilder.cxx | 158 ++++----- PWGLF/Tasks/Resonances/lambda1405analysis.cxx | 327 +++++++++++------- 2 files changed, 282 insertions(+), 203 deletions(-) diff --git a/PWGLF/TableProducer/Common/kinkBuilder.cxx b/PWGLF/TableProducer/Common/kinkBuilder.cxx index 330e0f6315b..42afc9f37bc 100644 --- a/PWGLF/TableProducer/Common/kinkBuilder.cxx +++ b/PWGLF/TableProducer/Common/kinkBuilder.cxx @@ -16,6 +16,7 @@ #include "PWGLF/DataModel/LFKinkDecayTables.h" #include "PWGLF/Utils/svPoolCreator.h" +#include "Common/DataModel/Centrality.h" #include "Common/Core/RecoDecay.h" #include "Common/Core/trackUtilities.h" @@ -60,7 +61,9 @@ using std::array; using VBracket = o2::math_utils::Bracket; using TracksFull = soa::Join; using TracksFullMc = soa::Join; -using McRecoCollisionsCent = soa::Join; +using CollisionsCent = soa::Join; +using McRecoCollisions = soa::Join; +using McRecoCollisionsCent = soa::Join; enum KinkDecayType { kSigmaMinusToPiMinusNeutron = 0, kSigmaPlusToPiPlusNeutron, @@ -167,13 +170,15 @@ struct kinkBuilder { Configurable unlikeSignBkg{"unlikeSignBkg", false, "Use unlike sign background"}; Configurable skipBkgCands{"skipBkgCands", true, "Skip bkg candidates in MC process"}; + // Centrality selection + Configurable centralityMin{"centralityMin", 0., "Minimum centrality accepted in SP/EP computation (not applied in resolution process)"}; + Configurable centralityMax{"centralityMax", 100., "Maximum centrality accepted in SP/EP computation (not applied in resolution process)"}; + // CCDB options 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"}; - // PDG codes - // histogram axes ConfigurableAxis rigidityBins{"rigidityBins", {200, -10.f, 10.f}, "Binning for rigidity #it{p}^{TPC}/#it{z}"}; ConfigurableAxis dedxBins{"dedxBins", {1000, 0.f, 1000.f}, "Binning for dE/dx"}; @@ -184,6 +189,12 @@ struct kinkBuilder { ConfigurableAxis dcaDaugToPVBins{"dcaDaugToPVBins", {200, 0.f, 20.f}, "Binning for DCA of daughter to PV"}; ConfigurableAxis mothDecRad2Bins{"mothDecRad2Bins", {100, 0, 100}, "Binning for mother decay radius squared"}; + // Filter collisions with centrality information (for the process with centrality selection) + using CollisionsCentSel = soa::Filtered; + using McRecoCollisionsCentSel = soa::Filtered; + + Filter filterCentrality = aod::cent::centFT0C >= centralityMin && aod::cent::centFT0C <= centralityMax; + // std vector of candidates std::vector kinkCandidates; @@ -258,16 +269,16 @@ struct kinkBuilder { massAxis = AxisSpec{100, 3.85, 4.25, "m (GeV/#it{c}^{2})"}; } - h2DeDxDaugSel = qaRegistry.add("h2DeDxDaugSel", ";p_{TPC}/z (GeV/#it{c}); dE/dx", HistType::kTH2F, {rigidityAxis, dedxAxis}); - h2KinkAnglePt = qaRegistry.add("h2KinkAnglePt", "; p_{T} (GeV/#it{c}); #theta_{kink} (deg)", HistType::kTH2F, {ptAxis, kinkAngleAxis}); - h2MothMassPt = qaRegistry.add("h2MothMassPt", "; p_{T} (GeV/#it{c}); m (GeV/#it{c}^{2})", HistType::kTH2F, {ptAxis, massAxis}); - h2ClsMapPtMoth = qaRegistry.add("h2ClsMapPtMoth", "; p_{T} (GeV/#it{c}); ITS cluster map", HistType::kTH2F, {ptAxis, itsClusterMapAxis}); - h2ClsMapPtDaug = qaRegistry.add("h2ClsMapPtDaug", "; p_{T} (GeV/#it{c}); ITS cluster map", HistType::kTH2F, {ptAxis, itsClusterMapAxis}); - h2ItsClsMothBeforeSel = qaRegistry.add("h2ItsClsMothBeforeSel", "; Tot Cls; IB clusters", HistType::kTH2F, {{8, -0.5, 7.5}, {4, -0.5, 3.5}}); - h2ItsClsDaugBeforeSel = qaRegistry.add("h2ItsClsDaugBeforeSel", "; Tot Cls; IB clusters", HistType::kTH2F, {{8, -0.5, 7.5}, {4, -0.5, 3.5}}); + h2DeDxDaugSel = qaRegistry.add("h2DeDxDaugSel", "h2DeDxDaugSel; p_{TPC}/z (GeV/#it{c}); dE/dx", HistType::kTH2F, {rigidityAxis, dedxAxis}); + h2KinkAnglePt = qaRegistry.add("h2KinkAnglePt", "h2KinkAnglePt; p_{T} (GeV/#it{c}); #theta_{kink} (deg)", HistType::kTH2F, {ptAxis, kinkAngleAxis}); + h2MothMassPt = qaRegistry.add("h2MothMassPt", "h2MothMassPt; p_{T} (GeV/#it{c}); m (GeV/#it{c}^{2})", HistType::kTH2F, {ptAxis, massAxis}); + h2ClsMapPtMoth = qaRegistry.add("h2ClsMapPtMoth", "h2ClsMapPtMoth; p_{T} (GeV/#it{c}); ITS cluster map", HistType::kTH2F, {ptAxis, itsClusterMapAxis}); + h2ClsMapPtDaug = qaRegistry.add("h2ClsMapPtDaug", "h2ClsMapPtDaug; p_{T} (GeV/#it{c}); ITS cluster map", HistType::kTH2F, {ptAxis, itsClusterMapAxis}); + h2ItsClsMothBeforeSel = qaRegistry.add("h2ItsClsMothBeforeSel", "h2ItsClsMothBeforeSel; Tot Cls; IB clusters", HistType::kTH2F, {{8, -0.5, 7.5}, {4, -0.5, 3.5}}); + h2ItsClsDaugBeforeSel = qaRegistry.add("h2ItsClsDaugBeforeSel", "h2ItsClsDaugBeforeSel; Tot Cls; IB clusters", HistType::kTH2F, {{8, -0.5, 7.5}, {4, -0.5, 3.5}}); // QA - hSelMotherQA = qaRegistry.add("hSelMotherQA", ";;Q_{Mother}", {HistType::kTH2F, {{9, -0.5, 8.5}, {2, -0.5, 1.5}}}); + hSelMotherQA = qaRegistry.add("hSelMotherQA", "hSelMotherQA;;Q_{Mother}", {HistType::kTH2F, {{9, -0.5, 8.5}, {2, -0.5, 1.5}}}); hSelMotherQA->GetYaxis()->SetBinLabel(1, "-"); hSelMotherQA->GetYaxis()->SetBinLabel(2, "+"); hSelMotherQA->GetXaxis()->SetBinLabel(1, "All"); @@ -281,7 +292,7 @@ struct kinkBuilder { hSelMotherQA->GetXaxis()->SetBinLabel(9, "Pt sel"); // QA - hSelDaugQA = qaRegistry.add("hSelDaugQA", ";;Q_{Daug}", {HistType::kTH2F, {{8, -0.5, 7.5}, {2, -0.5, 1.5}}}); + hSelDaugQA = qaRegistry.add("hSelDaugQA", "hSelDaugQA;;Q_{Daug}", {HistType::kTH2F, {{8, -0.5, 7.5}, {2, -0.5, 1.5}}}); hSelDaugQA->GetYaxis()->SetBinLabel(1, "-"); hSelDaugQA->GetYaxis()->SetBinLabel(2, "+"); hSelDaugQA->GetXaxis()->SetBinLabel(1, "All"); @@ -293,7 +304,7 @@ struct kinkBuilder { hSelDaugQA->GetXaxis()->SetBinLabel(7, "TPC Crossed rows vs findable sel"); hSelDaugQA->GetXaxis()->SetBinLabel(8, "TPC cls sel"); - hSelKinkedTrackQA = qaRegistry.add("hSelKinkedTrackQA", ";;(Q_{Mother}, Q_{Daughter})", {HistType::kTH2F, {{11, -0.5, 10.5}, {4, -0.5, 3.5}}}); + hSelKinkedTrackQA = qaRegistry.add("hSelKinkedTrackQA", "hSelKinkedTrackQA;;(Q_{Mother}, Q_{Daughter})", {HistType::kTH2F, {{11, -0.5, 10.5}, {4, -0.5, 3.5}}}); hSelKinkedTrackQA->GetYaxis()->SetBinLabel(1, "(+,+)"); hSelKinkedTrackQA->GetYaxis()->SetBinLabel(2, "(-,-)"); hSelKinkedTrackQA->GetYaxis()->SetBinLabel(3, "(+,-)"); @@ -309,29 +320,29 @@ struct kinkBuilder { hSelKinkedTrackQA->GetXaxis()->SetBinLabel(9, "MothDecRad2InsideL4"); hSelKinkedTrackQA->GetXaxis()->SetBinLabel(10, "MothDaugLastLayers"); hSelKinkedTrackQA->GetXaxis()->SetBinLabel(11, "MothLastLayerRadiusCheck"); - hMothDaughSignsInit = qaRegistry.add("hMothDaughSignsInit", "; Sign Mother; Sign Daughter", {HistType::kTH2F, {{3, -1.5, 1.5}, {3, -1.5, 1.5}}}); - hMothDaughSignsFinal = qaRegistry.add("hMothDaughSignsFinal", "; Sign Mother; Sign Daughter", {HistType::kTH2F, {{3, -1.5, 1.5}, {3, -1.5, 1.5}}}); - hZDiff = qaRegistry.add("hZDiff", "; #Delta z (#mu m);(Q_{Mother}, Q_{Daughter})", HistType::kTH2F, {zDiffBins, {4, -0.5, 3.5}}); + hMothDaughSignsInit = qaRegistry.add("hMothDaughSignsInit", "hMothDaughSignsInit; Sign Mother; Sign Daughter", {HistType::kTH2F, {{3, -1.5, 1.5}, {3, -1.5, 1.5}}}); + hMothDaughSignsFinal = qaRegistry.add("hMothDaughSignsFinal", "hMothDaughSignsFinal; Sign Mother; Sign Daughter", {HistType::kTH2F, {{3, -1.5, 1.5}, {3, -1.5, 1.5}}}); + hZDiff = qaRegistry.add("hZDiff", "hZDiff; #Delta z (#mu m);(Q_{Mother}, Q_{Daughter})", HistType::kTH2F, {zDiffBins, {4, -0.5, 3.5}}); hZDiff->GetYaxis()->SetBinLabel(1, "(+,+)"); hZDiff->GetYaxis()->SetBinLabel(2, "(-,-)"); hZDiff->GetYaxis()->SetBinLabel(3, "(+,-)"); hZDiff->GetYaxis()->SetBinLabel(4, "(-,+)"); - hPhiDiff = qaRegistry.add("hPhiDiff", "; #Delta #phi (rad);(Q_{Mother}, Q_{Daughter})", HistType::kTH2F, {phiDiffBins, {4, -0.5, 3.5}}); + hPhiDiff = qaRegistry.add("hPhiDiff", "hPhiDiff; #Delta #phi (rad);(Q_{Mother}, Q_{Daughter})", HistType::kTH2F, {phiDiffBins, {4, -0.5, 3.5}}); hPhiDiff->GetYaxis()->SetBinLabel(1, "(+,+)"); hPhiDiff->GetYaxis()->SetBinLabel(2, "(-,-)"); hPhiDiff->GetYaxis()->SetBinLabel(3, "(+,-)"); hPhiDiff->GetYaxis()->SetBinLabel(4, "(-,+)"); - hDCAMothToPV = qaRegistry.add("hDCAMothToPV", "; DCA moth to PV;(Q_{Mother}, Q_{Daughter})", HistType::kTH2F, {dcaMothToPVBins, {4, -0.5, 3.5}}); + hDCAMothToPV = qaRegistry.add("hDCAMothToPV", "hDCAMothToPV; DCA moth to PV;(Q_{Mother}, Q_{Daughter})", HistType::kTH2F, {dcaMothToPVBins, {4, -0.5, 3.5}}); hDCAMothToPV->GetYaxis()->SetBinLabel(1, "(+,+)"); hDCAMothToPV->GetYaxis()->SetBinLabel(2, "(-,-)"); hDCAMothToPV->GetYaxis()->SetBinLabel(3, "(+,-)"); hDCAMothToPV->GetYaxis()->SetBinLabel(4, "(-,+)"); - hDCADaugToPV = qaRegistry.add("hDCADaugToPV", "; DCA daug to PV;(Q_{Mother}, Q_{Daughter})", HistType::kTH2F, {dcaDaugToPVBins, {4, -0.5, 3.5}}); + hDCADaugToPV = qaRegistry.add("hDCADaugToPV", "hDCADaugToPV; DCA daug to PV;(Q_{Mother}, Q_{Daughter})", HistType::kTH2F, {dcaDaugToPVBins, {4, -0.5, 3.5}}); hDCADaugToPV->GetYaxis()->SetBinLabel(1, "(+,+)"); hDCADaugToPV->GetYaxis()->SetBinLabel(2, "(-,-)"); hDCADaugToPV->GetYaxis()->SetBinLabel(3, "(+,-)"); hDCADaugToPV->GetYaxis()->SetBinLabel(4, "(-,+)"); - hMothDecRad2 = qaRegistry.add("hMothDecRad2", "; Mother Dec. Radius Squared (cm^{2});(Q_{Mother}, Q_{Daughter})", HistType::kTH2F, {mothDecRad2Bins, {4, -0.5, 3.5}}); + hMothDecRad2 = qaRegistry.add("hMothDecRad2", "hMothDecRad2; Mother Dec. Radius Squared (cm^{2});(Q_{Mother}, Q_{Daughter})", HistType::kTH2F, {mothDecRad2Bins, {4, -0.5, 3.5}}); hMothDecRad2->GetYaxis()->SetBinLabel(1, "(+,+)"); hMothDecRad2->GetYaxis()->SetBinLabel(2, "(-,-)"); hMothDecRad2->GetYaxis()->SetBinLabel(3, "(+,-)"); @@ -342,23 +353,23 @@ struct kinkBuilder { } mBBparamsDaug[5] = cfgBetheBlochParams->get("Daughter", "resolution"); - if (doprocessMc) { + if (doprocessMc || doprocessMcWCent) { if (skipBkgCands) { - hRecCandidates = qaRegistry.add("hRecCandidates", ";Counts;", {HistType::kTH2F, {{kNMatchedDecays, -0.5, kNMatchedDecays-0.5}, absPtAxis}}); + hRecCandidates = qaRegistry.add("hRecCandidates", "hRecCandidates;Counts;", {HistType::kTH2F, {{kNMatchedDecays, -0.5, static_cast(kNMatchedDecays)-0.5}, absPtAxis}}); hRecCandidates->GetXaxis()->SetBinLabel(1, "#Sigma^{-} #rightarrow n#pi^{-}"); hRecCandidates->GetXaxis()->SetBinLabel(2, "#Sigma^{+} #rightarrow n#pi^{+}"); hRecCandidates->GetXaxis()->SetBinLabel(3, "#Sigma^{+} #rightarrow p#pi^{0}"); } - hGenCandidates = qaRegistry.add("hGenCandidates", ";Counts;", {HistType::kTH2F, {{kNMatchedDecays, -0.5, kNMatchedDecays-0.5}, absPtAxis}}); + hGenCandidates = qaRegistry.add("hGenCandidates", "hGenCandidates;Counts;", {HistType::kTH2F, {{kNMatchedDecays, -0.5, static_cast(kNMatchedDecays)-0.5}, absPtAxis}}); hGenCandidates->GetXaxis()->SetBinLabel(1, "#Sigma^{-} #rightarrow n#pi^{-}"); hGenCandidates->GetXaxis()->SetBinLabel(2, "#Sigma^{+} #rightarrow n#pi^{+}"); hGenCandidates->GetXaxis()->SetBinLabel(3, "#Sigma^{+} #rightarrow p#pi^{0}"); - hRecPtKinkAngle[kSigmaMinusToPiMinusNeutron] = qaRegistry.add("hRecPtKinkAngleSigmaMinusToPiMinusNeutron", ";Counts;", {HistType::kTH2F, {absPtAxis, kinkAngleAxis}}); - hRecPtKinkAngle[kSigmaPlusToPiPlusNeutron] = qaRegistry.add("hRecPtKinkAngleSigmaPlusToPiPlusNeutron", ";Counts;", {HistType::kTH2F, {absPtAxis, kinkAngleAxis}}); - hRecPtKinkAngle[kSigmaPlusToProtonPi0] = qaRegistry.add("hRecPtKinkAngleSigmaPlusToProtonPi0", ";Counts;", {HistType::kTH2F, {absPtAxis, kinkAngleAxis}}); - hGenPtKinkAngle[kSigmaMinusToPiMinusNeutron] = qaRegistry.add("hGenPtKinkAngleSigmaMinusToPiMinusNeutron", ";Counts;", {HistType::kTH2F, {absPtAxis, kinkAngleAxis}}); - hGenPtKinkAngle[kSigmaPlusToPiPlusNeutron] = qaRegistry.add("hGenPtKinkAngleSigmaPlusToPiPlusNeutron", ";Counts;", {HistType::kTH2F, {absPtAxis, kinkAngleAxis}}); - hGenPtKinkAngle[kSigmaPlusToProtonPi0] = qaRegistry.add("hGenPtKinkAngleSigmaPlusToProtonPi0", ";Counts;", {HistType::kTH2F, {absPtAxis, kinkAngleAxis}}); + hRecPtKinkAngle[kSigmaMinusToPiMinusNeutron] = qaRegistry.add("hRecPtKinkAngleSigmaMinusToPiMinusNeutron", "Rec Pt vs KinkAngle #Sigma^{-} #rightarrow #pi^{-}n;Counts;", {HistType::kTH2F, {absPtAxis, kinkAngleAxis}}); + hRecPtKinkAngle[kSigmaPlusToPiPlusNeutron] = qaRegistry.add("hRecPtKinkAngleSigmaPlusToPiPlusNeutron", "Rec Pt vs KinkAngle #Sigma^{+} #rightarrow #pi^{+}n;Counts;", {HistType::kTH2F, {absPtAxis, kinkAngleAxis}}); + hRecPtKinkAngle[kSigmaPlusToProtonPi0] = qaRegistry.add("hRecPtKinkAngleSigmaPlusToProtonPi0", "Rec Pt vs KinkAngle #Sigma^{+} #rightarrow p#pi^{0};Counts;", {HistType::kTH2F, {absPtAxis, kinkAngleAxis}}); + hGenPtKinkAngle[kSigmaMinusToPiMinusNeutron] = qaRegistry.add("hGenPtKinkAngleSigmaMinusToPiMinusNeutron", "Gen Pt vs KinkAngle #Sigma^{-} #rightarrow #pi^{-}n;Counts;", {HistType::kTH2F, {absPtAxis, kinkAngleAxis}}); + hGenPtKinkAngle[kSigmaPlusToPiPlusNeutron] = qaRegistry.add("hGenPtKinkAngleSigmaPlusToPiPlusNeutron", "Gen Pt vs KinkAngle #Sigma^{+} #rightarrow #pi^{+}n;Counts;", {HistType::kTH2F, {absPtAxis, kinkAngleAxis}}); + hGenPtKinkAngle[kSigmaPlusToProtonPi0] = qaRegistry.add("hGenPtKinkAngleSigmaPlusToProtonPi0", "Gen Pt vs KinkAngle #Sigma^{+} #rightarrow p#pi^{0};Counts;", {HistType::kTH2F, {absPtAxis, kinkAngleAxis}}); } } @@ -650,8 +661,8 @@ struct kinkBuilder { } } - void processData(aod::Collisions const& collisions, TracksFull const& tracks, aod::AmbiguousTracks const& ambiTracks, aod::BCs const& bcs) - { + template + void fillOutputsData(const TColls& collisions, const TTracks& tracks, const aod::AmbiguousTracks& ambiTracks, const aod::BCs& bcs) { kinkCandidates.clear(); // LOG(info) << "[kink] "; // LOG(info) << "[kink] *****************************************"; @@ -683,7 +694,18 @@ struct kinkBuilder { } } } - PROCESS_SWITCH(kinkBuilder, processData, "Data processing", false); + + void processData(aod::Collisions const& collisions, TracksFull const& tracks, aod::AmbiguousTracks const& ambiTracks, aod::BCs const& bcs) + { + fillOutputsData(collisions, tracks, ambiTracks, bcs); + } + PROCESS_SWITCH(kinkBuilder, processData, "Data processing", true); + + void processDataWCentSel(CollisionsCentSel const& collisions, TracksFull const& tracks, aod::AmbiguousTracks const& ambiTracks, aod::BCs const& bcs) + { + fillOutputsData(collisions, tracks, ambiTracks, bcs); + } + PROCESS_SWITCH(kinkBuilder, processDataWCentSel, "Data processing with centrality selection", false); template int matchKinkDecay(const TMother& motherPart, const aod::McParticles& mcParticles) { @@ -697,13 +719,7 @@ struct kinkBuilder { pdgCodeNeutralDaug = (pdgMother > 0) ? +PDG_t::kNeutron : -PDG_t::kNeutron; pdgCodeChargedDaug = (pdgMother > 0) ? +PDG_t::kPiMinus : +PDG_t::kPiPlus; finState = {pdgCodeChargedDaug, pdgCodeNeutralDaug}; // Both decay channels have the same neutral daughter - // if constexpr (checkKinkDaugPdg) { - // LOG(info) << "[kink] Matching Sigma- decay with daughter PDG " << pdgDaug << " and neutral daughter PDG " << pdgCodeNeutralDaug; - // } if (RecoDecay::isMatchedMCGen(mcParticles, motherPart, pdgMother, finState, true, &sign, DepthMcMatchMax)) { - // if constexpr (checkKinkDaugPdg) { - // LOG(info) << "[kink] Matched Sigma- to pi- neutron decay"; - // } return kSigmaMinusToPiMinusNeutron; } break; @@ -713,13 +729,7 @@ struct kinkBuilder { pdgCodeNeutralDaug = (pdgMother > 0) ? +PDG_t::kNeutron : -PDG_t::kNeutron; pdgCodeChargedDaug = (pdgMother > 0) ? +PDG_t::kPiPlus : +PDG_t::kPiMinus; finState = {pdgCodeChargedDaug, pdgCodeNeutralDaug}; - // if constexpr (checkKinkDaugPdg) { - // LOG(info) << "[kink] Matching Sigma+ decay with daughter PDG " << pdgDaug << " and neutral daughter PDG " << pdgCodeNeutralDaug; - // } if (RecoDecay::isMatchedMCGen(mcParticles, motherPart, pdgMother, finState, true, &sign, DepthMcMatchMax)) { - // if constexpr (checkKinkDaugPdg) { - // LOG(info) << "[kink] Matched Sigma+ to pi+ neutron decay"; - // } return kSigmaPlusToPiPlusNeutron; } @@ -727,13 +737,7 @@ struct kinkBuilder { pdgCodeNeutralDaug = (pdgMother > 0) ? +PDG_t::kPi0 : -PDG_t::kPi0; pdgCodeChargedDaug = (pdgMother > 0) ? +PDG_t::kProton : -PDG_t::kProton; finState = {pdgCodeChargedDaug, pdgCodeNeutralDaug}; - // if constexpr (checkKinkDaugPdg) { - // LOG(info) << "[kink] Matching Sigma+ decay with daughter PDG " << pdgDaug << " and neutral daughter PDG " << pdgCodeNeutralDaug; - // } if (RecoDecay::isMatchedMCGen(mcParticles, motherPart, pdgMother, finState, true, &sign, DepthMcMatchMax)) { - // if constexpr (checkKinkDaugPdg) { - // LOG(info) << "[kink] Matched Sigma+ to proton pi0 decay"; - // } return kSigmaPlusToProtonPi0; } break; @@ -746,28 +750,18 @@ struct kinkBuilder { return -1; } - void processMc(aod::McCollisions const& mcGenCollisions, - McRecoCollisionsCent const& mcRecoCollisions, - aod::McParticles const& mcParticles, - TracksFullMc const& tracksMc, - aod::AmbiguousTracks const& ambiTracksMc, - aod::BCs const& bcs) - { + template + void fillOutputsMc(const TColls& mcRecoCollisions, const TTracks& tracksMc, const aod::AmbiguousTracks& ambiTracksMc, const aod::BCs& bcs, const aod::McParticles& mcParticles) { kinkCandidates.clear(); - // LOG(info) << "[kink] "; - // LOG(info) << "[kink] *****************************************"; - // LOG(info) << "[kink] Processing " << tracksMc.size() << " tracks and " << mcRecoCollisions.size() << " collisions"; buildSvPool(mcRecoCollisions, tracksMc, ambiTracksMc, bcs); auto& kinkPool = svCreator.getSVCandPool(mcRecoCollisions, !unlikeSignBkg); - // LOG(info) << "[kink] Number of kink candidates in the pool: " << kinkPool.size(); - int countSameMother = 0; for (const auto& svCand : kinkPool) { // Perform matching of the kink candidate auto trackMoth = tracksMc.rawIteratorAt(svCand.tr0Idx); - auto genMothPart = trackMoth.mcParticle_as(); + auto genMothPart = trackMoth.template mcParticle_as(); auto trackDaug = tracksMc.rawIteratorAt(svCand.tr1Idx); - auto genDaugPart = trackDaug.mcParticle_as(); + auto genDaugPart = trackDaug.template mcParticle_as(); int decayChannel{-1}; if (skipBkgCands) { @@ -781,22 +775,7 @@ struct kinkBuilder { if (genDaugMothIdx != trackMoth.mcParticleId()) { continue; // Skip candidates where the daughter is not coming from the selected mother } - // LOG(info) << "[kink] "; - // LOG(info) << "[kink] -----------------"; - // LOG(info) << "[kink] MC particle offset: " << mcParticles.offset(); - // LOG(info) << "[kink] [Rec] Matching kink candidate with mother PDG " << genMothPart.pdgCode() << " and daughter PDG " << genDaugPart.pdgCode(); - // LOG(info) << "[kink] Mother global index: " << trackMoth.globalIndex() << ", trackMoth.mcParticleId(): " << trackMoth.mcParticleId() << ", PDG code: " << genMothPart.pdgCode(); - // for (const auto& mcMother : genDaugPart.mothers_as()) { - // LOG(info) << "[kink] Daughter's mother PDG: " << mcMother.pdgCode() << ", global index: " << mcMother.globalIndex(); - // } - // std::vector arrDaughIdxs = {}; - // RecoDecay::getDaughters(genMothPart, &arrDaughIdxs, std::array{0}, 1); - // for (auto iProng = 0u; iProng < arrDaughIdxs.size(); ++iProng) { - // auto daughI = mcParticles.rawIteratorAt(arrDaughIdxs[iProng]); - // LOG(info) << "[kink] Daughter PDG: " << daughI.pdgCode() << ", global index: " << daughI.globalIndex(); - // } decayChannel = matchKinkDecay(genMothPart, mcParticles); - // LOG(info) << "[kink] Decay channel matched: " << decayChannel; if (decayChannel < 0) { continue; // Skip candidates that do not match the decay channels of interest } @@ -848,7 +827,28 @@ struct kinkBuilder { } } } + + void processMc(aod::McCollisions const& mcGenCollisions, + McRecoCollisions const& mcRecoCollisions, + aod::McParticles const& mcParticles, + TracksFullMc const& tracksMc, + aod::AmbiguousTracks const& ambiTracksMc, + aod::BCs const& bcs) + { + fillOutputsMc(mcRecoCollisions, tracksMc, ambiTracksMc, bcs, mcParticles); + } PROCESS_SWITCH(kinkBuilder, processMc, "MC processing", false); + + void processMcWCent(aod::McCollisions const& mcGenCollisions, + McRecoCollisionsCentSel const& mcRecoCollisions, + aod::McParticles const& mcParticles, + TracksFullMc const& tracksMc, + aod::AmbiguousTracks const& ambiTracksMc, + aod::BCs const& bcs) + { + fillOutputsMc(mcRecoCollisions, tracksMc, ambiTracksMc, bcs, mcParticles); + } + PROCESS_SWITCH(kinkBuilder, processMcWCent, "MC processing with centrality selection", false); }; WorkflowSpec diff --git a/PWGLF/Tasks/Resonances/lambda1405analysis.cxx b/PWGLF/Tasks/Resonances/lambda1405analysis.cxx index 0511d02bcc9..17fd745b794 100644 --- a/PWGLF/Tasks/Resonances/lambda1405analysis.cxx +++ b/PWGLF/Tasks/Resonances/lambda1405analysis.cxx @@ -46,7 +46,8 @@ using namespace o2::framework::expressions; using TracksFull = soa::Join; using CollisionsFull = soa::Join; using CollisionsFullWCentQVecs = soa::Join; -using CollisionsFullMC = soa::Join; +using CollisionsFullMc = soa::Join; +using CollisionsFullMcWCent = soa::Join; enum L1405DecayType { kSigmaMinusPiToPiPiNeutron = 0, kSigmaPlusPiToPiPiNeutron, @@ -117,10 +118,12 @@ struct lambda1405analysis { Configurable cutNTpcClusPi{"cutNTpcClusPi", 90, "Minimum number of Tpc clusters for pion candidate"}; Configurable cutNSigTpc{"cutNSigTpc", 3, "NSigTpcPion"}; Configurable cutNSigTof{"cutNSigTof", 3, "NSigTofPion"}; - Configurable cutSigmaQtAPMin{"cutSigmaQtAPMin", 0.17, "Lower limit for Sigma qT"}; - Configurable cutSigmaQtAPMax{"cutSigmaQtAPMax", 0.2, "Upper limit for Sigma qT"}; - Configurable cutSigmaAlphaAPMin{"cutSigmaAlphaAPMin", -0.9, "Lower limit for Sigma qT"}; - Configurable cutSigmaAlphaAPMax{"cutSigmaAlphaAPMax", -0.4, "Upper limit for Sigma qT"}; + Configurable cutSigmaQtAPMin{"cutSigmaQtAPMin", "0.17", "Functional form of minimum qT(alpha) selection"}; + Configurable cutSigmaQtAPMax{"cutSigmaQtAPMax", "0.2", "Functional form of minimum qT(alpha) selection"}; + Configurable cutMinBachPiPtVsL1405Pt{"cutMinBachPiPtVsL1405Pt", "0", "Functional form of minimum qT(alpha) selection"}; + Configurable cutMaxBachPiPtVsL1405Pt{"cutMaxBachPiPtVsL1405Pt", "10", "Functional form of minimum qT(alpha) selection"}; + Configurable cutMinSigmaPtVsL1405Pt{"cutMinSigmaPtVsL1405Pt", "0", "Functional form of minimum qT(alpha) selection"}; + Configurable cutMaxSigmaPtVsL1405Pt{"cutMaxSigmaPtVsL1405Pt", "10", "Functional form of minimum qT(alpha) selection"}; // Configurables for flow analysis Configurable centralityMin{"centralityMin", 0., "Minimum centrality accepted in SP/EP computation (not applied in resolution process)"}; @@ -135,6 +138,18 @@ struct lambda1405analysis { Configurable doLSBkg{"doLikeSignBkg", false, "Use like-sign background"}; Configurable useTof{"useTof", false, "Use Tof for PId for pion candidates"}; + TF1* funcMinQtAlphaAP = nullptr; + TF1* funcMaxQtAlphaAP = nullptr; + TF1* funcMinBachPiPtVsL1405Pt = nullptr; + TF1* funcMaxBachPiPtVsL1405Pt = nullptr; + TF1* funcMinSigmaPtVsL1405Pt = nullptr; + TF1* funcMaxSigmaPtVsL1405Pt = nullptr; + + using CollisionsCentSel = soa::Filtered; + using McRecoCollisionsCentSel = soa::Filtered; + + Filter filterCentrality = aod::cent::centFT0C >= centralityMin && aod::cent::centFT0C <= centralityMax; + // Configurable axes ConfigurableAxis axisCent{"axisCent", {10000, 0., 100.}, ""}; ConfigurableAxis axisPtL1405{"axisPtL1405", {100, 0., 10.}, ""}; @@ -239,13 +254,15 @@ struct lambda1405analysis { rLambda1405.add("h2BachPiPtNSigTpc", "h2BachPiPtNSigTpc", {HistType::kTH2F, {ptKinkDaugAxis, nSigmaAxis}}); // Selection QA - rSelections.add("hSelectionsL1405", "hSelectionsL1405", {HistType::kTH1D, {{5, -0.f, 4.5f}}}); + rSelections.add("hSelectionsL1405", "hSelectionsL1405", {HistType::kTH1D, {{7, -0.f, 6.5f}}}); rSelections.get(HIST("hSelectionsL1405"))->GetXaxis()->SetBinLabel(1, "All"); rSelections.get(HIST("hSelectionsL1405"))->GetXaxis()->SetBinLabel(2, "Passed Sigma sel"); rSelections.get(HIST("hSelectionsL1405"))->GetXaxis()->SetBinLabel(3, "Passed Bach PID sel"); rSelections.get(HIST("hSelectionsL1405"))->GetXaxis()->SetBinLabel(4, "Upper mass sel"); - rSelections.get(HIST("hSelectionsL1405"))->GetXaxis()->SetBinLabel(5, "Accepted"); - rSelections.add("hSelectionsSigmaPlus", "hSelectionsSigmaPlus", {HistType::kTH1D, {{8, -0.f, 7.5f}}}); + rSelections.get(HIST("hSelectionsL1405"))->GetXaxis()->SetBinLabel(5, "Correl. #Sigma #it{p}_{T} sel"); + rSelections.get(HIST("hSelectionsL1405"))->GetXaxis()->SetBinLabel(6, "Correl. bach. #pi #it{p}_{T} sel"); + rSelections.get(HIST("hSelectionsL1405"))->GetXaxis()->SetBinLabel(7, "Accepted"); + rSelections.add("hSelectionsSigmaPlus", "hSelectionsSigmaPlus", {HistType::kTH1D, {{7, -0.f, 6.5f}}}); rSelections.get(HIST("hSelectionsSigmaPlus"))->GetXaxis()->SetBinLabel(1, "All"); rSelections.get(HIST("hSelectionsSigmaPlus"))->GetXaxis()->SetBinLabel(2, "Passed kink sel"); rSelections.get(HIST("hSelectionsSigmaPlus"))->GetXaxis()->SetBinLabel(3, "Passed mass sel"); @@ -253,8 +270,7 @@ struct lambda1405analysis { rSelections.get(HIST("hSelectionsSigmaPlus"))->GetXaxis()->SetBinLabel(5, "Passed cutDCAtoPvPiFromSigma"); rSelections.get(HIST("hSelectionsSigmaPlus"))->GetXaxis()->SetBinLabel(6, "Passed cutSigmaRadius"); rSelections.get(HIST("hSelectionsSigmaPlus"))->GetXaxis()->SetBinLabel(7, "Passed cutSigmaQtAP"); - rSelections.get(HIST("hSelectionsSigmaPlus"))->GetXaxis()->SetBinLabel(8, "Passed cutSigmaAlphaAP"); - rSelections.add("hSelectionsSigmaMinus", "hSelectionsSigmaMinus", {HistType::kTH1D, {{8, -0.f, 7.5f}}}); + rSelections.add("hSelectionsSigmaMinus", "hSelectionsSigmaMinus", {HistType::kTH1D, {{7, -0.f, 6.5f}}}); rSelections.get(HIST("hSelectionsSigmaMinus"))->GetXaxis()->SetBinLabel(1, "All"); rSelections.get(HIST("hSelectionsSigmaMinus"))->GetXaxis()->SetBinLabel(2, "Passed kink sel"); rSelections.get(HIST("hSelectionsSigmaMinus"))->GetXaxis()->SetBinLabel(3, "Passed mass sel"); @@ -262,7 +278,6 @@ struct lambda1405analysis { rSelections.get(HIST("hSelectionsSigmaMinus"))->GetXaxis()->SetBinLabel(5, "Passed cutDCAtoPvPiFromSigma"); rSelections.get(HIST("hSelectionsSigmaMinus"))->GetXaxis()->SetBinLabel(6, "Passed cutSigmaRadius"); rSelections.get(HIST("hSelectionsSigmaMinus"))->GetXaxis()->SetBinLabel(7, "Passed cutSigmaQtAP"); - rSelections.get(HIST("hSelectionsSigmaMinus"))->GetXaxis()->SetBinLabel(8, "Passed cutSigmaAlphaAP"); rSelections.add("hSelectionsBachPi", "hSelectionsBachPi", {HistType::kTH1D, {{6, -0.f, 5.5f}}}); rSelections.get(HIST("hSelectionsBachPi"))->GetXaxis()->SetBinLabel(1, "All"); rSelections.get(HIST("hSelectionsBachPi"))->GetXaxis()->SetBinLabel(2, "Sign sel"); @@ -293,10 +308,10 @@ struct lambda1405analysis { rLambda1405.add("hSparseFlowL1405", "THn for SP", HistType::kTHnSparseF, axesFlow); } - if (doprocessMC) { + if (doprocessMc || doprocessMcWCentSel) { // Add MC histograms if needed, to sigmaminus rLambda1405.add("hRecoL1405", "hRecoL1405;;Counts", {HistType::kTH2F, {{6, -0.5, 5.5}, ptAxis}}); - rLambda1405.get(HIST("hRecoL1405"))->GetXaxis()->SetBinLabel(1, "Reconstructed #Lambda(1405)"); + rLambda1405.get(HIST("hRecoL1405"))->GetXaxis()->SetBinLabel(1, "Reconstructed #Sigma (1405)"); rLambda1405.get(HIST("hRecoL1405"))->GetXaxis()->SetBinLabel(2, "Has MC particle"); rLambda1405.get(HIST("hRecoL1405"))->GetXaxis()->SetBinLabel(3, "Has #Sigma daug"); rLambda1405.get(HIST("hRecoL1405"))->GetXaxis()->SetBinLabel(4, "Has bach #pi"); @@ -306,6 +321,17 @@ struct lambda1405analysis { rLambda1405.get(HIST("hGenL1405"))->GetXaxis()->SetBinLabel(1, "#Lambda(1405) #rightarrow #Sigma^{-} #pi^{+} #rightarrow n #pi^{-} #pi^{+}"); rLambda1405.get(HIST("hGenL1405"))->GetXaxis()->SetBinLabel(2, "#Lambda(1405) #rightarrow #Sigma^{+} #pi^{-} #rightarrow n #pi^{+} #pi^{-}"); rLambda1405.get(HIST("hGenL1405"))->GetXaxis()->SetBinLabel(3, "#Lambda(1405) #rightarrow #Sigma^{+} #pi^{-} #rightarrow p #pi^{0} #pi^{-}"); + rLambda1405.add("h2GenSigmaMinusArmPod", "h2GenSigmaMinusArmPod", {HistType::kTH2F, {alphaAxis, qtAxis}}); + rLambda1405.add("h2GenSigmaPlusArmPod", "h2GenSigmaPlusArmPod", {HistType::kTH2F, {alphaAxis, qtAxis}}); + rLambda1405.add("h2GenPtVsBachPtSigmaMinusPiToPiPiNeutron", "h2GenPtVsBachPtSigmaMinusPiToPiPiNeutron;#Lambda(1405) #it{p}_{T} (GeV/c); Bach #pi #it{p}_{T} (GeV/c)", {HistType::kTH2F, {ptAxis, ptAxis}}); + rLambda1405.add("h2GenPtVsBachPtSigmaPlusPiToPiPiN", "h2GenPtVsBachPtSigmaPlusPiToPiPiN;#Lambda(1405) #it{p}_{T} (GeV/c); Bach #pi #it{p}_{T} (GeV/c)", {HistType::kTH2F, {ptAxis, ptAxis}}); + rLambda1405.add("h2GenPtVsBachPtSigmaPlusPiToPiPiP", "h2GenPtVsBachPtSigmaPlusPiToPiPiP;#Lambda(1405) #it{p}_{T} (GeV/c); Bach #pi #it{p}_{T} (GeV/c)", {HistType::kTH2F, {ptAxis, ptAxis}}); + rLambda1405.add("h2GenPtVsSigmaPtSigmaMinusPiToPiPiNeutron", "h2GenPtVsBachPtSigmaMinusPiToPiPiNeutron;#Lambda(1405) #it{p}_{T} (GeV/c); #Sigma #it{p}_{T} (GeV/c)", {HistType::kTH2F, {ptAxis, ptAxis}}); + rLambda1405.add("h2GenPtVsSigmaPtSigmaPlusPiToPiPiN", "h2GenPtVsBachPtSigmaPlusPiToPiPiN;#Lambda(1405) #it{p}_{T} (GeV/c); #Sigma #it{p}_{T} (GeV/c)", {HistType::kTH2F, {ptAxis, ptAxis}}); + rLambda1405.add("h2GenPtVsSigmaPtSigmaPlusPiToPiPiP", "h2GenPtVsBachPtSigmaPlusPiToPiPiP;#Lambda(1405) #it{p}_{T} (GeV/c); #Sigma #it{p}_{T} (GeV/c)", {HistType::kTH2F, {ptAxis, ptAxis}}); + rLambda1405.add("h2GenSigmaPtVsKinkPtSigmaMinusPiToPiPiNeutron", "h2GenSigmaPtVsKinkPtSigmaMinusPiToPiPiNeutron;#Sigma #it{p}_{T} (GeV/c); Kink #it{p}_{T} (GeV/c)", {HistType::kTH2F, {ptAxis, ptAxis}}); + rLambda1405.add("h2GenSigmaPtVsKinkPtSigmaPlusPiToPiPiN", "h2GenSigmaPtVsKinkPtSigmaPlusPiTo Pi PiNeutron;#Sigma #it{p}_{T} (GeV/c); Kink #it{p}_{T} (GeV/c)", {HistType::kTH2F, {ptAxis, ptAxis}}); + rLambda1405.add("h2GenSigmaPtVsKinkPtSigmaPlusPiToPiPiP", "h2GenSigmaPtVsKinkPtSigmaPlusPiToPiPiP;#Sigma #it{p}_{T} (GeV/c); Kink #it{p}_{T} (GeV/c)", {HistType::kTH2F, {ptAxis, ptAxis}}); rLambda1405.add("h2MassResolutionFromSigmaMinus", "h2MassResolutionFromSigmaMinus", {HistType::kTH2F, {lambda1405MassAxis, massResolutionAxis}}); rLambda1405.add("h2PtResolutionFromSigmaMinus", "h2PtResolutionFromSigmaMinus", {HistType::kTH2F, {ptAxis, ptResolutionAxis}}); // Add MC histograms if needed, to sigmaplus @@ -314,6 +340,20 @@ struct lambda1405analysis { // rLambda1405.add("h2PtMassMCWithSigmaDaug", "h2PtMassMCWithSigmaDaug", {HistType::kTH2F, {ptAxis, lambda1405MassAxis}}); // rLambda1405.add("h2PtMassMCNoSigmaDaug", "h2PtMassMCNoSigmaDaug", {HistType::kTH2F, {ptAxis, lambda1405MassAxis}}); } + + // Functional selections + funcMinQtAlphaAP = new TF1("funcMinQtAlphaAP", Form("%s", cutSigmaQtAPMin.value.data()), -1, 1); + LOGF(info, "funcMinQtAlphaAP: %s", Form("%s", cutSigmaQtAPMin.value.data())); + funcMaxQtAlphaAP = new TF1("funcMaxQtAlphaAP", Form("%s", cutSigmaQtAPMax.value.data()), -1, 1); + LOGF(info, "funcMaxQtAlphaAP: %s", Form("%s", cutSigmaQtAPMax.value.data())); + funcMinBachPiPtVsL1405Pt = new TF1("funcMinBachPiPtVsL1405Pt", Form("%s", cutMinBachPiPtVsL1405Pt.value.data()), 0., 100); + LOGF(info, "funcMinBachPiPtVsL1405Pt: %s", Form("%s", cutMinBachPiPtVsL1405Pt.value.data())); + funcMaxBachPiPtVsL1405Pt = new TF1("funcMaxBachPiPtVsL1405Pt", Form("%s", cutMaxBachPiPtVsL1405Pt.value.data()), 0., 100); + LOGF(info, "funcMaxBachPiPtVsL1405Pt: %s", Form("%s", cutMaxBachPiPtVsL1405Pt.value.data())); + funcMinSigmaPtVsL1405Pt = new TF1("funcMinSigmaPtVsL1405Pt", Form("%s", cutMinSigmaPtVsL1405Pt.value.data()), 0., 100); + LOGF(info, "funcMinSigmaPtVsL1405Pt: %s", Form("%s", cutMinSigmaPtVsL1405Pt.value.data())); + funcMaxSigmaPtVsL1405Pt = new TF1("funcMaxSigmaPtVsL1405Pt", Form("%s", cutMaxSigmaPtVsL1405Pt.value.data()), 0., 100); + LOGF(info, "funcMaxSigmaPtVsL1405Pt: %s", Form("%s", cutMaxSigmaPtVsL1405Pt.value.data())); } float alphaAP(const std::array& momMother, const std::array& momKink) @@ -577,7 +617,8 @@ struct lambda1405analysis { lambda1405Cand.sigmaRadius = sigmaRad; lambda1405Cand.kinkDcaDauToPv = sigmaCand.dcaDaugPv(); - if (lambda1405Cand.sigmaQtAP < cutSigmaQtAPMin || lambda1405Cand.sigmaQtAP > cutSigmaQtAPMax) { + if (lambda1405Cand.sigmaQtAP < funcMinQtAlphaAP->Eval(lambda1405Cand.sigmaAlphaAP) || + lambda1405Cand.sigmaQtAP > funcMaxQtAlphaAP->Eval(lambda1405Cand.sigmaAlphaAP)) { return; } if (sigmaCand.mothSign() < 0) { @@ -586,15 +627,6 @@ struct lambda1405analysis { rSelections.fill(HIST("hSelectionsSigmaPlus"), 6); // Passed mass sel } - if (lambda1405Cand.sigmaAlphaAP < cutSigmaAlphaAPMin || lambda1405Cand.sigmaAlphaAP > cutSigmaAlphaAPMax) { - return; - } - if (sigmaCand.mothSign() < 0) { - rSelections.fill(HIST("hSelectionsSigmaMinus"), 7); // Passed mass sel - } else { - rSelections.fill(HIST("hSelectionsSigmaPlus"), 7); // Passed mass sel - } - // Kink daughter properties lambda1405Cand.kinkDauId = kinkDauTrack.globalIndex(); lambda1405Cand.kinkPt = kinkDauTrack.pt(); @@ -652,8 +684,22 @@ struct lambda1405analysis { lambda1405Cand.pz = sigmaMom[2] + piMom[2]; lambda1405Cand.phi = std::atan2(lambda1405Cand.py, lambda1405Cand.px); lambda1405Cand.scalarProd = -1; - fillHistosLambda1405(lambda1405Cand, piTrack); + + // Check correlations between transverse momenta of L1405, sigma and bachelor pi + if (std::hypot(sigmaCand.pxMoth(), sigmaCand.pyMoth()) < funcMinSigmaPtVsL1405Pt->Eval(lambda1405Cand.pt()) || + std::hypot(sigmaCand.pxMoth(), sigmaCand.pyMoth()) > funcMaxSigmaPtVsL1405Pt->Eval(lambda1405Cand.pt())) { + continue; + } rSelections.fill(HIST("hSelectionsL1405"), 4); // Accepted + + if (piTrack.pt() < funcMinBachPiPtVsL1405Pt->Eval(lambda1405Cand.pt()) || + piTrack.pt() > funcMaxBachPiPtVsL1405Pt->Eval(lambda1405Cand.pt())) { + continue; + } + rSelections.fill(HIST("hSelectionsL1405"), 5); // Accepted + + fillHistosLambda1405(lambda1405Cand, piTrack); + rSelections.fill(HIST("hSelectionsL1405"), 6); // Accepted selectedCandidates.push_back(lambda1405Cand); } } @@ -678,16 +724,22 @@ struct lambda1405analysis { return isKinkFromSigma; // Return true if the kink comes from the Sigma } - void processData(CollisionsFull::iterator const& collision, aod::KinkCands const& kinkCands, TracksFull const& tracks) - { + template + void fillOutputData(const TCollision& collision, const TCand& sigmaCands, const TTrack& kinkDauTrack) { + if constexpr (requires{collision.centFT0C();}) { + if (collision.centFT0C() < centralityMin || collision.centFT0C() > centralityMax) { + LOG(info) << "Skipping collision due to centrality cut"; + return; + } + } if (std::abs(collision.posZ()) > cutzvertex || !collision.sel8()) { return; } rEventSelection.fill(HIST("hVertexZRec"), collision.posZ()); - for (const auto& sigmaCand : kinkCands) { + for (const auto& sigmaCand : sigmaCands) { std::vector selectedCandidates; - constructCollCandidates(sigmaCand, tracks, selectedCandidates); - for (const auto& lambda1405Cand : selectedCandidates) { + constructCollCandidates(sigmaCand, kinkDauTrack, selectedCandidates); + for (auto& lambda1405Cand : selectedCandidates) { if (lambda1405Cand.isSigmaMinus) { rSigmaMinus.fill(HIST("h2SigmaMinusMassVsLambdaMass"), lambda1405Cand.sigmaSign * lambda1405Cand.sigmaPt, lambda1405Cand.sigmaMinusMass); if (lambda1405Cand.hasPiKink) { @@ -719,99 +771,82 @@ struct lambda1405analysis { lambda1405Cand.kinkDcaDauToPv, lambda1405Cand.bachPiNSigTpc, lambda1405Cand.bachPiNSigTof); } + if constexpr (requires{collision.qvecFT0CRe();}) { + float const xQVec = collision.qvecFT0CRe(); + float const yQVec = collision.qvecFT0CIm(); + float const cos2Phi = std::cos(2 * lambda1405Cand.phi); + float const sin2Phi = std::sin(2 * lambda1405Cand.phi); + lambda1405Cand.scalarProd = cos2Phi * xQVec + sin2Phi * yQVec; + float const ptCand = lambda1405Cand.pt(); + rLambda1405.fill(HIST("hSparseFlowL1405"), ptCand, lambda1405Cand.massL1405, collision.centFT0C(), lambda1405Cand.scalarProd); + LOG(info) << "Filled flow sparse for candidate with pt " << ptCand << ", mass " << lambda1405Cand.massL1405 << ", centrality " << collision.centFT0C() << " and scalar product " << lambda1405Cand.scalarProd; + if (fillFlowTree) { + if (downSampleFactor < 1.) { + float const pseudoRndm = ptCand * 1000. - static_cast(ptCand * 1000); + if (ptCand < ptDownSampleMax && pseudoRndm >= downSampleFactor) { + continue; + } + } + outputDataFlowTable(ptCand, lambda1405Cand.massL1405, + lambda1405Cand.sigmaMinusMass, lambda1405Cand.sigmaPlusMass, + lambda1405Cand.sigmaAlphaAP, lambda1405Cand.sigmaQtAP, + lambda1405Cand.kinkPiNSigTpc, lambda1405Cand.kinkPiNSigTof, + lambda1405Cand.kinkPrNSigTpc, lambda1405Cand.kinkPrNSigTof, + lambda1405Cand.kinkDcaDauToPv, + lambda1405Cand.bachPiNSigTpc, lambda1405Cand.bachPiNSigTof, + lambda1405Cand.scalarProd, collision.centFT0C()); + } + } } } } + + void processData(CollisionsFull::iterator const& collision, aod::KinkCands const& kinkCands, TracksFull const& tracks) + { + fillOutputData(collision, kinkCands, tracks); + } PROCESS_SWITCH(lambda1405analysis, processData, "Data processing", true); - void processDataWCentQVecs(CollisionsFullWCentQVecs::iterator const& collision, aod::KinkCands const& kinkCands, TracksFull const& tracks) + void processDataWCentQVecs(CollisionsCentSel::iterator const& collision, aod::KinkCands const& kinkCands, TracksFull const& tracks) { - LOG(info) << "Processing collision with centrality " << collision.centFT0C() << " and Q vector (" << collision.qvecFT0CRe() << ", " << collision.qvecFT0CIm() << ")"; - if (collision.centFT0C() < centralityMin || collision.centFT0C() > centralityMax) { - LOG(info) << "Skipping collision due to centrality cut"; - return; - } - if (std::abs(collision.posZ()) > cutzvertex || !collision.sel8()) { - LOG(info) << "Skipping collision due to vertex cut"; - return; - } - rEventSelection.fill(HIST("hVertexZRec"), collision.posZ()); - for (const auto& sigmaCand : kinkCands) { - LOG(info) << "Processing Sigma candidate with mother sign " << sigmaCand.mothSign() << " and pt " << sigmaCand.ptMoth(); - std::vector selectedCandidates; - constructCollCandidates(sigmaCand, tracks, selectedCandidates); - for (auto& lambda1405Cand : selectedCandidates) { - if (lambda1405Cand.isSigmaMinus) { - rSigmaMinus.fill(HIST("h2SigmaMinusMassVsLambdaMass"), lambda1405Cand.sigmaSign * lambda1405Cand.sigmaPt, lambda1405Cand.sigmaMinusMass); - if (lambda1405Cand.hasPiKink) { - rSigmaMinus.fill(HIST("h2KinkPiPtNSigTofSigmaMinus"), lambda1405Cand.sigmaSign * lambda1405Cand.piPt, lambda1405Cand.bachPiNSigTof); - } - if (lambda1405Cand.hasPrKink) { - rSigmaMinus.fill(HIST("h2KinkPrPtNSigTofSigmaMinus"), lambda1405Cand.sigmaSign * lambda1405Cand.piPt, lambda1405Cand.bachPiNSigTof); - } - } - if (lambda1405Cand.isSigmaPlus) { - rSigmaPlus.fill(HIST("h2SigmaPlusMassVsLambdaMass"), lambda1405Cand.sigmaSign * lambda1405Cand.sigmaPt, lambda1405Cand.sigmaPlusMass); - rSigmaPlus.fill(HIST("h2KinkPrPtNSigTofSigmaPlus"), lambda1405Cand.sigmaSign * lambda1405Cand.piPt, lambda1405Cand.bachPiNSigTof); - } - LOG(info) << "Filled Sigma histograms for candidate with Sigma pt " << lambda1405Cand.sigmaPt << " and mass " << (lambda1405Cand.isSigmaMinus ? lambda1405Cand.sigmaMinusMass : lambda1405Cand.sigmaPlusMass); - float const xQVec = collision.qvecFT0CRe(); - float const yQVec = collision.qvecFT0CIm(); - float const cos2Phi = std::cos(2 * lambda1405Cand.phi); - float const sin2Phi = std::sin(2 * lambda1405Cand.phi); - lambda1405Cand.scalarProd = cos2Phi * xQVec + sin2Phi * yQVec; - float const ptCand = lambda1405Cand.pt(); - rLambda1405.fill(HIST("hSparseFlowL1405"), ptCand, lambda1405Cand.massL1405, collision.centFT0C(), lambda1405Cand.scalarProd); - LOG(info) << "Filled flow sparse for candidate with pt " << ptCand << ", mass " << lambda1405Cand.massL1405 << ", centrality " << collision.centFT0C() << " and scalar product " << lambda1405Cand.scalarProd; - if (fillFlowTree) { - if (downSampleFactor < 1.) { - float const pseudoRndm = ptCand * 1000. - static_cast(ptCand * 1000); - if (ptCand < ptDownSampleMax && pseudoRndm >= downSampleFactor) { - continue; - } - } - outputDataFlowTable(ptCand, lambda1405Cand.massL1405, - lambda1405Cand.sigmaMinusMass, lambda1405Cand.sigmaPlusMass, - lambda1405Cand.sigmaAlphaAP, lambda1405Cand.sigmaQtAP, - lambda1405Cand.kinkPiNSigTpc, lambda1405Cand.kinkPiNSigTof, - lambda1405Cand.kinkPrNSigTpc, lambda1405Cand.kinkPrNSigTof, - lambda1405Cand.kinkDcaDauToPv, - lambda1405Cand.bachPiNSigTpc, lambda1405Cand.bachPiNSigTof, - lambda1405Cand.scalarProd, collision.centFT0C()); - } - LOG(info) << "Filled flow tree for candidate with pt " << ptCand; - } - } + fillOutputData(collision, kinkCands, tracks); } PROCESS_SWITCH(lambda1405analysis, processDataWCentQVecs, "Data processing with centrality and Q vectors info", false); template - int matchGenDecay(const TMother& motherPart, const aod::McParticles& mcParticles) { - LOG(info) << "Matching MC decay for particle with PDG code " << motherPart.pdgCode() << " and index " << motherPart.globalIndex(); + int matchGenDecay(const TMother& motherPart, const aod::McParticles& mcParticles, std::array& daugsIdxs) { int pdgMother = motherPart.pdgCode(); int8_t sign = 0; int MaxDepthFinState = 2; // Maximum depth to look for the decay chain - // Match L(1405) --> n pi- pi+ final state + std::vector arrL1405Daugs = {}; std::array finalState; + int decayChannel = -1; + + // Match L(1405) --> n pi- pi+ final state if (pdgMother > 0) { // Change sign of neutral decay products finalState = {PDG_t::kNeutron, PDG_t::kPiMinus, PDG_t::kPiPlus}; } else { finalState = {-PDG_t::kNeutron, PDG_t::kPiMinus, PDG_t::kPiPlus}; } - - if (RecoDecay::isMatchedMCGen(mcParticles, motherPart, lambda1405PdgCode, finalState, true, &sign, MaxDepthFinState)) { + if (RecoDecay::isMatchedMCGen(mcParticles, motherPart, lambda1405PdgCode, finalState, true, &sign, MaxDepthFinState)) { // Match the intermediate Sigma - std::vector arrDaughIdxs = {}; - RecoDecay::getDaughters(motherPart, &arrDaughIdxs, std::array{0}, 1); - for (auto iProng = 0u; iProng < arrDaughIdxs.size(); ++iProng) { - auto daughI = mcParticles.rawIteratorAt(arrDaughIdxs[iProng]); + RecoDecay::getDaughters(motherPart, &arrL1405Daugs, std::array{0}, 1); + for (auto iProng = 0u; iProng < arrL1405Daugs.size(); ++iProng) { + auto daughI = mcParticles.rawIteratorAt(arrL1405Daugs[iProng]); if (std::abs(daughI.pdgCode()) == PDG_t::kSigmaPlus) { - return kSigmaPlusPiToPiPiNeutron; + decayChannel = kSigmaPlusPiToPiPiNeutron; + daugsIdxs[1] = arrL1405Daugs[iProng]; + } + if (std::abs(daughI.pdgCode()) == PDG_t::kSigmaMinus) { + decayChannel = kSigmaMinusPiToPiPiNeutron; + daugsIdxs[1] = arrL1405Daugs[iProng]; + } + if (std::abs(daughI.pdgCode()) == PDG_t::kPiPlus) { + daugsIdxs[0] = arrL1405Daugs[iProng]; } } - return kSigmaMinusPiToPiPiNeutron; } // Match L(1405) --> p pi0 pi+ final state, only possible for Sigma+ @@ -821,26 +856,48 @@ struct lambda1405analysis { finalState = {PDG_t::kProton, -PDG_t::kPi0, PDG_t::kPiMinus}; } - if (RecoDecay::isMatchedMCGen(mcParticles, motherPart, lambda1405PdgCode, finalState, true, &sign, MaxDepthFinState)) { + if (RecoDecay::isMatchedMCGen(mcParticles, motherPart, lambda1405PdgCode, finalState, true, &sign, MaxDepthFinState)) { // Match the intermediate Sigma - return kSigmaPlusPiToPiPiProton; + RecoDecay::getDaughters(motherPart, &arrL1405Daugs, std::array{0}, 1); + for (auto iProng = 0u; iProng < arrL1405Daugs.size(); ++iProng) { + auto daughI = mcParticles.rawIteratorAt(arrL1405Daugs[iProng]); + if (std::abs(daughI.pdgCode()) == PDG_t::kSigmaPlus) { + decayChannel = kSigmaPlusPiToPiPiProton; + daugsIdxs[1] = arrL1405Daugs[iProng]; + } + if (std::abs(daughI.pdgCode()) == PDG_t::kPiPlus) { + daugsIdxs[0] = arrL1405Daugs[iProng]; + } + } } - return -1; + std::vector arrSigmaDaugs = {}; + auto sigmaDaug = mcParticles.rawIteratorAt(daugsIdxs[1]); + RecoDecay::getDaughters(sigmaDaug, &arrSigmaDaugs, std::array{0}, 1); + for (auto iProng = 0u; iProng < arrSigmaDaugs.size(); ++iProng) { + auto daughSigma = mcParticles.rawIteratorAt(arrSigmaDaugs[iProng]); + if (std::abs(daughSigma.pdgCode()) == PDG_t::kPiPlus || std::abs(daughSigma.pdgCode()) == PDG_t::kProton) { + daugsIdxs[2] = arrSigmaDaugs[iProng]; + } + } + + return decayChannel; } - void processMC(CollisionsFullMC const& collisions, aod::KinkCands const& kinkCands, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const& particlesMC, TracksFull const& tracks) - { + template + void fillOutputMc(const TCollision& collisions, const aod::KinkCands& sigmaCands, const aod::McTrackLabels& trackLabelsMC, const TTrack& tracks, const aod::McParticles& particlesMC) { for (const auto& collision : collisions) { + if constexpr (requires{collision.centFT0C();}) { + if (collision.centFT0C() < centralityMin || collision.centFT0C() > centralityMax) { + return; + } + } if (std::abs(collision.posZ()) > cutzvertex) { // || !collision.sel8()) { continue; } rEventSelection.fill(HIST("hVertexZRec"), collision.posZ()); - auto sigmaCandsPerCol = kinkCands.sliceBy(mKinkPerCol, collision.globalIndex()); + auto sigmaCandsPerCol = sigmaCands.sliceBy(mKinkPerCol, collision.globalIndex()); auto tracksPerCol = tracks.sliceBy(mPerColTracks, collision.globalIndex()); - if (sigmaCandsPerCol.size() != 0) { - LOG(info) << "Processing collision with " << sigmaCandsPerCol.size() << " Sigma kink candidates and " << tracksPerCol.size() << " tracks"; - } for (const auto& sigmaCand : sigmaCandsPerCol) { std::vector selectedCandidates; constructCollCandidates(sigmaCand, tracksPerCol, selectedCandidates); @@ -852,11 +909,7 @@ struct lambda1405analysis { auto mcLabPiKink = trackLabelsMC.rawIteratorAt(lambda1405Cand.kinkDauId); auto mcLabSigma = trackLabelsMC.rawIteratorAt(lambda1405Cand.sigmaId); auto mcLabPi = trackLabelsMC.rawIteratorAt(lambda1405Cand.piId); - if (!mcLabSigma.has_mcParticle() || mcLabPiKink.has_mcParticle() || mcLabPi.has_mcParticle()) { - LOG(info) << "Skipping candidate due to missing MC association: " - << "Sigma MC assoc: " << mcLabSigma.has_mcParticle() << ", " - << "Kink daughter MC assoc: " << mcLabPiKink.has_mcParticle() << ", " - << "Bachelor Pi MC assoc: " << mcLabPi.has_mcParticle(); + if (!mcLabSigma.has_mcParticle() || !mcLabPiKink.has_mcParticle() || !mcLabPi.has_mcParticle()) { continue; // Skip if no valid MC association } rLambda1405.fill(HIST("hRecoL1405"), 1., lambda1405Cand.pt()); // All with associated MC particle @@ -870,25 +923,16 @@ struct lambda1405analysis { bool isSigmaPlusToPrKink = checkSigmaKinkMC(mcTrackSigma, mcTrackKink, 3222, 2212, particlesMC); if (!isSigmaMinusKink && !isSigmaPlusToPiKink && !isSigmaPlusToPrKink) { - LOG(info) << "Skipping candidate due to failed MC kink decay check: " - << "isSigmaMinusKink: " << isSigmaMinusKink << ", " - << "isSigmaPlusToPiKink: " << isSigmaPlusToPiKink << ", " - << "isSigmaPlusToPrKink: " << isSigmaPlusToPrKink; continue; // Skip if not a valid Sigma kink decay } rLambda1405.fill(HIST("hRecoL1405"), 2., lambda1405Cand.pt()); // Has kink decay in MC if (std::abs(mcTrackPi.pdgCode()) != 211) { - LOG(info) << "Skipping candidate due to bachelor Pi not being a pion in MC: " - << "Bachelor Pi PDG code: " << mcTrackPi.pdgCode(); continue; // Skip if not a valid pion candidate } rLambda1405.fill(HIST("hRecoL1405"), 3., lambda1405Cand.pt()); // Has bach pi if (!mcTrackSigma.has_mothers() || !mcTrackPi.has_mothers()) { - LOG(info) << "Skipping candidate due to missing mothers in MC: " - << "Sigma has mothers: " << mcTrackSigma.has_mothers() << ", " - << "Pi has mothers: " << mcTrackPi.has_mothers(); continue; // Skip if no mothers found } rLambda1405.fill(HIST("hRecoL1405"), 4., lambda1405Cand.pt()); // Has mothers for Sigma and Pi @@ -899,13 +943,11 @@ struct lambda1405analysis { for (const auto& sigmaMother : mcTrackSigma.mothers_as()) { if (piMother.globalIndex() == sigmaMother.globalIndex() && std::abs(piMother.pdgCode()) == lambda1405PdgCode) { lambda1405Id = piMother.globalIndex(); - LOG(info) << "Found common mother for Sigma and Pi with global index: " << lambda1405Id; break; // Found the mother, exit loop } } } if (lambda1405Id == -1) { - LOG(info) << "Skipping candidate due to Sigma and pion not sharing the same lambda1405 candidate"; continue; // Skip if the Sigma and pion do not share the same lambda1405 candidate } rLambda1405.fill(HIST("hRecoL1405"), 4., lambda1405Cand.pt()); // Has same mother @@ -944,14 +986,51 @@ struct lambda1405analysis { if (std::abs(mcPart.pdgCode()) != lambda1405PdgCode) { continue; // Only consider Lambda(1405) candidates } - int decayChannel = matchGenDecay(mcPart, particlesMC); + std::array idxsSigmaKinkBachPi{-1, -1, -1}; + int decayChannel = matchGenDecay(mcPart, particlesMC, idxsSigmaKinkBachPi); if (decayChannel == -1) { continue; // Skip if it doesn't match the decay channels of interest } rLambda1405.fill(HIST("hGenL1405"), decayChannel, mcPart.pt()); + + auto bachPi = particlesMC.rawIteratorAt(idxsSigmaKinkBachPi[0]); + auto sigmaDaug = particlesMC.rawIteratorAt(idxsSigmaKinkBachPi[1]); + auto sigmaKinkDaug = particlesMC.rawIteratorAt(idxsSigmaKinkBachPi[2]); + // Generated Armenteros-Podolanski variables + float genSigmaAlphaAP = alphaAP({sigmaDaug.px(), sigmaDaug.py(), sigmaDaug.pz()}, {sigmaKinkDaug.px(), sigmaKinkDaug.py(), sigmaKinkDaug.pz()}); + float genSigmaQtAP = qtAP({sigmaDaug.px(), sigmaDaug.py(), sigmaDaug.pz()}, {sigmaKinkDaug.px(), sigmaKinkDaug.py(), sigmaKinkDaug.pz()}); + if (decayChannel == kSigmaMinusPiToPiPiNeutron) { + rLambda1405.fill(HIST("h2GenSigmaMinusArmPod"), genSigmaAlphaAP, genSigmaQtAP); + rLambda1405.fill(HIST("h2GenPtVsBachPtSigmaMinusPiToPiPiNeutron"), mcPart.pt(), bachPi.pt()); + rLambda1405.fill(HIST("h2GenPtVsSigmaPtSigmaMinusPiToPiPiNeutron"), mcPart.pt(), sigmaDaug.pt()); + rLambda1405.fill(HIST("h2GenSigmaPtVsKinkPtSigmaMinusPiToPiPiNeutron"), sigmaDaug.pt(), sigmaKinkDaug.pt()); + } + if (decayChannel == kSigmaPlusPiToPiPiNeutron) { + rLambda1405.fill(HIST("h2GenSigmaPlusArmPod"), genSigmaAlphaAP, genSigmaQtAP); + rLambda1405.fill(HIST("h2GenPtVsBachPtSigmaPlusPiToPiPiN"), mcPart.pt(), bachPi.pt()); + rLambda1405.fill(HIST("h2GenPtVsSigmaPtSigmaPlusPiToPiPiN"), mcPart.pt(), sigmaDaug.pt()); + rLambda1405.fill(HIST("h2GenSigmaPtVsKinkPtSigmaPlusPiToPiPiN"), sigmaDaug.pt(), sigmaKinkDaug.pt()); + } + if (decayChannel == kSigmaPlusPiToPiPiProton) { + rLambda1405.fill(HIST("h2GenSigmaPlusArmPod"), genSigmaAlphaAP, genSigmaQtAP); + rLambda1405.fill(HIST("h2GenPtVsBachPtSigmaPlusPiToPiPiP"), mcPart.pt(), bachPi.pt()); + rLambda1405.fill(HIST("h2GenPtVsSigmaPtSigmaPlusPiToPiPiP"), mcPart.pt(), sigmaDaug.pt()); + rLambda1405.fill(HIST("h2GenSigmaPtVsKinkPtSigmaPlusPiToPiPiP"), sigmaDaug.pt(), sigmaKinkDaug.pt()); + } } } - PROCESS_SWITCH(lambda1405analysis, processMC, "MC processing", false); + + void processMc(CollisionsFullMc const& collisions, aod::KinkCands const& kinkCands, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const& particlesMC, TracksFull const& tracks) + { + fillOutputMc(collisions, kinkCands, trackLabelsMC, tracks, particlesMC); + } + PROCESS_SWITCH(lambda1405analysis, processMc, "MC processing", false); + + void processMcWCentSel(McRecoCollisionsCentSel const& collisions, aod::KinkCands const& kinkCands, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const& particlesMC, TracksFull const& tracks) + { + fillOutputMc(collisions, kinkCands, trackLabelsMC, tracks, particlesMC); + } + PROCESS_SWITCH(lambda1405analysis, processMcWCentSel, "MC processing with centrality selection", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { From 71ae571855060a038515c70fb83ae2170705b24e Mon Sep 17 00:00:00 2001 From: marcellocosti Date: Tue, 2 Jun 2026 05:52:12 +0200 Subject: [PATCH 3/5] Delete comment --- PWGLF/TableProducer/Common/kinkBuilder.cxx | 4 ---- 1 file changed, 4 deletions(-) diff --git a/PWGLF/TableProducer/Common/kinkBuilder.cxx b/PWGLF/TableProducer/Common/kinkBuilder.cxx index 42afc9f37bc..f091dfe4f49 100644 --- a/PWGLF/TableProducer/Common/kinkBuilder.cxx +++ b/PWGLF/TableProducer/Common/kinkBuilder.cxx @@ -664,13 +664,9 @@ struct kinkBuilder { template void fillOutputsData(const TColls& collisions, const TTracks& tracks, const aod::AmbiguousTracks& ambiTracks, const aod::BCs& bcs) { kinkCandidates.clear(); - // LOG(info) << "[kink] "; - // LOG(info) << "[kink] *****************************************"; - // LOG(info) << "[kink] Processing " << tracks.size() << " tracks and " << collisions.size() << " collisions"; buildSvPool(collisions, tracks, ambiTracks, bcs); auto& kinkPool = svCreator.getSVCandPool(collisions, !unlikeSignBkg); - // LOG(info) << "[kink] Number of kink candidates in the pool: " << kinkPool.size(); bool isAccepted = false; for (const auto& svCand : kinkPool) { buildKinkCand(svCand, tracks, isAccepted, collisions); From 0de9f2a666fb1d67e65ce915b8b5ef8cdf09c9bc Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Tue, 2 Jun 2026 03:52:44 +0000 Subject: [PATCH 4/5] Please consider the following formatting changes --- PWGLF/DataModel/LFLambda1405Table.h | 8 +- PWGLF/TableProducer/Common/kinkBuilder.cxx | 44 ++++---- PWGLF/Tasks/Resonances/lambda1405analysis.cxx | 102 +++++++++--------- 3 files changed, 81 insertions(+), 73 deletions(-) diff --git a/PWGLF/DataModel/LFLambda1405Table.h b/PWGLF/DataModel/LFLambda1405Table.h index cb472f5fc58..31721f634e3 100644 --- a/PWGLF/DataModel/LFLambda1405Table.h +++ b/PWGLF/DataModel/LFLambda1405Table.h @@ -49,8 +49,8 @@ DECLARE_SOA_COLUMN(NSigmaTPCPiDau, nSigmaTPCPiDau, float); //! Number of sigma DECLARE_SOA_COLUMN(NSigmaTOFPiDau, nSigmaTOFPiDau, float); //! Number of sigmas for the lambda1405 pion daughter in TOF // Flow columns -DECLARE_SOA_COLUMN(ScalarProd, scalarProd, float); //! Scalar product of the candidate -DECLARE_SOA_COLUMN(Centrality, centrality, float); //! Centrality of the candidate +DECLARE_SOA_COLUMN(ScalarProd, scalarProd, float); //! Scalar product of the candidate +DECLARE_SOA_COLUMN(Centrality, centrality, float); //! Centrality of the candidate // MC Columns DECLARE_SOA_COLUMN(PtMC, ptMC, float); //! pT of the candidate in MC @@ -73,8 +73,8 @@ DECLARE_SOA_TABLE(Lambda1405Cands, "AOD", "LAMBDA1405", lambda1405::NSigmaTPCPiDau, lambda1405::NSigmaTOFPiDau); DECLARE_SOA_TABLE(Lambda1405Flow, "AOD", "LAMBDA1405FLOW", - o2::soa::Index<>, - lambda1405::Pt, + o2::soa::Index<>, + lambda1405::Pt, lambda1405::Mass, lambda1405::SigmaMinusMass, lambda1405::SigmaPlusMass, lambda1405::AlphaAPSigma, lambda1405::QtAPSigma, lambda1405::NSigmaTPCPiKink, lambda1405::NSigmaTOFPiKink, diff --git a/PWGLF/TableProducer/Common/kinkBuilder.cxx b/PWGLF/TableProducer/Common/kinkBuilder.cxx index f091dfe4f49..496f986d158 100644 --- a/PWGLF/TableProducer/Common/kinkBuilder.cxx +++ b/PWGLF/TableProducer/Common/kinkBuilder.cxx @@ -16,9 +16,9 @@ #include "PWGLF/DataModel/LFKinkDecayTables.h" #include "PWGLF/Utils/svPoolCreator.h" -#include "Common/DataModel/Centrality.h" #include "Common/Core/RecoDecay.h" #include "Common/Core/trackUtilities.h" +#include "Common/DataModel/Centrality.h" #include #include @@ -322,7 +322,7 @@ struct kinkBuilder { hSelKinkedTrackQA->GetXaxis()->SetBinLabel(11, "MothLastLayerRadiusCheck"); hMothDaughSignsInit = qaRegistry.add("hMothDaughSignsInit", "hMothDaughSignsInit; Sign Mother; Sign Daughter", {HistType::kTH2F, {{3, -1.5, 1.5}, {3, -1.5, 1.5}}}); hMothDaughSignsFinal = qaRegistry.add("hMothDaughSignsFinal", "hMothDaughSignsFinal; Sign Mother; Sign Daughter", {HistType::kTH2F, {{3, -1.5, 1.5}, {3, -1.5, 1.5}}}); - hZDiff = qaRegistry.add("hZDiff", "hZDiff; #Delta z (#mu m);(Q_{Mother}, Q_{Daughter})", HistType::kTH2F, {zDiffBins, {4, -0.5, 3.5}}); + hZDiff = qaRegistry.add("hZDiff", "hZDiff; #Delta z (#mu m);(Q_{Mother}, Q_{Daughter})", HistType::kTH2F, {zDiffBins, {4, -0.5, 3.5}}); hZDiff->GetYaxis()->SetBinLabel(1, "(+,+)"); hZDiff->GetYaxis()->SetBinLabel(2, "(-,-)"); hZDiff->GetYaxis()->SetBinLabel(3, "(+,-)"); @@ -355,12 +355,12 @@ struct kinkBuilder { if (doprocessMc || doprocessMcWCent) { if (skipBkgCands) { - hRecCandidates = qaRegistry.add("hRecCandidates", "hRecCandidates;Counts;", {HistType::kTH2F, {{kNMatchedDecays, -0.5, static_cast(kNMatchedDecays)-0.5}, absPtAxis}}); + hRecCandidates = qaRegistry.add("hRecCandidates", "hRecCandidates;Counts;", {HistType::kTH2F, {{kNMatchedDecays, -0.5, static_cast(kNMatchedDecays) - 0.5}, absPtAxis}}); hRecCandidates->GetXaxis()->SetBinLabel(1, "#Sigma^{-} #rightarrow n#pi^{-}"); hRecCandidates->GetXaxis()->SetBinLabel(2, "#Sigma^{+} #rightarrow n#pi^{+}"); hRecCandidates->GetXaxis()->SetBinLabel(3, "#Sigma^{+} #rightarrow p#pi^{0}"); } - hGenCandidates = qaRegistry.add("hGenCandidates", "hGenCandidates;Counts;", {HistType::kTH2F, {{kNMatchedDecays, -0.5, static_cast(kNMatchedDecays)-0.5}, absPtAxis}}); + hGenCandidates = qaRegistry.add("hGenCandidates", "hGenCandidates;Counts;", {HistType::kTH2F, {{kNMatchedDecays, -0.5, static_cast(kNMatchedDecays) - 0.5}, absPtAxis}}); hGenCandidates->GetXaxis()->SetBinLabel(1, "#Sigma^{-} #rightarrow n#pi^{-}"); hGenCandidates->GetXaxis()->SetBinLabel(2, "#Sigma^{+} #rightarrow n#pi^{+}"); hGenCandidates->GetXaxis()->SetBinLabel(3, "#Sigma^{+} #rightarrow p#pi^{0}"); @@ -437,15 +437,15 @@ struct kinkBuilder { if (candidate.itsNClsInnerBarrel() != 0) return false; hSelDaugQA->Fill(4.f, isPositive); - + if (candidate.itsNCls() >= 4) return false; hSelDaugQA->Fill(5.f, isPositive); - + if (candidate.tpcNClsCrossedRows() <= 0.8 * candidate.tpcNClsFindable()) return false; hSelDaugQA->Fill(6.f, isPositive); - + if (candidate.tpcNClsFound() <= nTPCClusMinDaug) return false; hSelDaugQA->Fill(7.f, isPositive); @@ -463,7 +463,7 @@ struct kinkBuilder { auto trackDaug = tracks.rawIteratorAt(svCand.tr1Idx); // Fill Selections QA histo - int chargeCombSvCand = 2*unlikeSignBkg + (trackMoth.sign() == -1 ? 1 : 0); + int chargeCombSvCand = 2 * unlikeSignBkg + (trackMoth.sign() == -1 ? 1 : 0); hSelKinkedTrackQA->Fill(0.f, chargeCombSvCand); // all candidates bin hMothDaughSignsInit->Fill(trackMoth.sign(), trackDaug.sign()); @@ -641,8 +641,9 @@ struct kinkBuilder { LOG(info) << "Task initialized for run " << mRunNumber << " with magnetic field " << mBz << " kZG"; } - template - void buildSvPool(const TColls& collisions, const TTracks& tracks, const TAmbiTracks& ambiguousTracks, const aod::BCs& bcs) { + template + void buildSvPool(const TColls& collisions, const TTracks& tracks, const TAmbiTracks& ambiguousTracks, const aod::BCs& bcs) + { svCreator.clearPools(); svCreator.fillBC2Coll(collisions, bcs); @@ -661,8 +662,9 @@ struct kinkBuilder { } } - template - void fillOutputsData(const TColls& collisions, const TTracks& tracks, const aod::AmbiguousTracks& ambiTracks, const aod::BCs& bcs) { + template + void fillOutputsData(const TColls& collisions, const TTracks& tracks, const aod::AmbiguousTracks& ambiTracks, const aod::BCs& bcs) + { kinkCandidates.clear(); buildSvPool(collisions, tracks, ambiTracks, bcs); @@ -703,15 +705,16 @@ struct kinkBuilder { } PROCESS_SWITCH(kinkBuilder, processDataWCentSel, "Data processing with centrality selection", false); - template - int matchKinkDecay(const TMother& motherPart, const aod::McParticles& mcParticles) { + template + int matchKinkDecay(const TMother& motherPart, const aod::McParticles& mcParticles) + { int pdgMother = motherPart.pdgCode(); int8_t sign = 0; int pdgCodeNeutralDaug{-1}, pdgCodeChargedDaug{-1}; std::array finState = {-1, -1}; switch (std::abs(pdgMother)) { case PDG_t::kSigmaMinus: { - // Swap the sign of the neutral decay products in case of anti-particles + // Swap the sign of the neutral decay products in case of anti-particles pdgCodeNeutralDaug = (pdgMother > 0) ? +PDG_t::kNeutron : -PDG_t::kNeutron; pdgCodeChargedDaug = (pdgMother > 0) ? +PDG_t::kPiMinus : +PDG_t::kPiPlus; finState = {pdgCodeChargedDaug, pdgCodeNeutralDaug}; // Both decay channels have the same neutral daughter @@ -723,7 +726,7 @@ struct kinkBuilder { case PDG_t::kSigmaPlus: { // Swap the sign of the neutral decay products in case of anti-particles pdgCodeNeutralDaug = (pdgMother > 0) ? +PDG_t::kNeutron : -PDG_t::kNeutron; - pdgCodeChargedDaug = (pdgMother > 0) ? +PDG_t::kPiPlus : +PDG_t::kPiMinus; + pdgCodeChargedDaug = (pdgMother > 0) ? +PDG_t::kPiPlus : +PDG_t::kPiMinus; finState = {pdgCodeChargedDaug, pdgCodeNeutralDaug}; if (RecoDecay::isMatchedMCGen(mcParticles, motherPart, pdgMother, finState, true, &sign, DepthMcMatchMax)) { return kSigmaPlusToPiPlusNeutron; @@ -746,8 +749,9 @@ struct kinkBuilder { return -1; } - template - void fillOutputsMc(const TColls& mcRecoCollisions, const TTracks& tracksMc, const aod::AmbiguousTracks& ambiTracksMc, const aod::BCs& bcs, const aod::McParticles& mcParticles) { + template + void fillOutputsMc(const TColls& mcRecoCollisions, const TTracks& tracksMc, const aod::AmbiguousTracks& ambiTracksMc, const aod::BCs& bcs, const aod::McParticles& mcParticles) + { kinkCandidates.clear(); buildSvPool(mcRecoCollisions, tracksMc, ambiTracksMc, bcs); @@ -828,7 +832,7 @@ struct kinkBuilder { McRecoCollisions const& mcRecoCollisions, aod::McParticles const& mcParticles, TracksFullMc const& tracksMc, - aod::AmbiguousTracks const& ambiTracksMc, + aod::AmbiguousTracks const& ambiTracksMc, aod::BCs const& bcs) { fillOutputsMc(mcRecoCollisions, tracksMc, ambiTracksMc, bcs, mcParticles); @@ -839,7 +843,7 @@ struct kinkBuilder { McRecoCollisionsCentSel const& mcRecoCollisions, aod::McParticles const& mcParticles, TracksFullMc const& tracksMc, - aod::AmbiguousTracks const& ambiTracksMc, + aod::AmbiguousTracks const& ambiTracksMc, aod::BCs const& bcs) { fillOutputsMc(mcRecoCollisions, tracksMc, ambiTracksMc, bcs, mcParticles); diff --git a/PWGLF/Tasks/Resonances/lambda1405analysis.cxx b/PWGLF/Tasks/Resonances/lambda1405analysis.cxx index 17fd745b794..dd4b51076b7 100644 --- a/PWGLF/Tasks/Resonances/lambda1405analysis.cxx +++ b/PWGLF/Tasks/Resonances/lambda1405analysis.cxx @@ -76,21 +76,21 @@ struct lambda1405candidate { float sigmaAlphaAP = -1; // Alpha of the Sigma float sigmaQtAP = -1; // qT of the Sigma float kinkPt = -1; // pT of the kink daughter - float kinkPiNSigTpc = -1; // Number of sigmas for the pion candidate from Sigma kink in Tpc - float kinkPiNSigTof = -1; // Number of sigmas for the pion candidate from Sigma kink in Tof - float kinkPrNSigTpc = -1; // Number of sigmas for the proton candidate from Sigma kink in Tpc - float kinkPrNSigTof = -1; // Number of sigmas for the proton candidate from Sigma kink in Tof + float kinkPiNSigTpc = -1; // Number of sigmas for the pion candidate from Sigma kink in Tpc + float kinkPiNSigTof = -1; // Number of sigmas for the pion candidate from Sigma kink in Tof + float kinkPrNSigTpc = -1; // Number of sigmas for the proton candidate from Sigma kink in Tpc + float kinkPrNSigTof = -1; // Number of sigmas for the proton candidate from Sigma kink in Tof float kinkDcaDauToPv = -1; // DCA of the kink daughter to the primary vertex float sigmaRadius = -1; // Radius of the Sigma decay vertex float piPt = -1; // pT of the pion daughter float bachPiNSigTpc = -1; // Number of sigmas for the pion candidate float bachPiNSigTof = -1; // Number of sigmas for the pion candidate using Tof - int kinkDauId = 0; // Id of the pion from Sigma decay in MC - int sigmaId = 0; // Id of the Sigma candidate in MC - int piId = 0; // Id of the pion candidate in MC + int kinkDauId = 0; // Id of the pion from Sigma decay in MC + int sigmaId = 0; // Id of the Sigma candidate in MC + int piId = 0; // Id of the pion candidate in MC - float scalarProd = -1; // Scalar product for flow analysis + float scalarProd = -1; // Scalar product for flow analysis }; struct lambda1405analysis { @@ -470,7 +470,8 @@ struct lambda1405analysis { } template - void fillHistosSigma(const lambda1405candidate& lambda1405Cand, const TCand& sigmaCand, const TTrack& kinkDauTrack) { + void fillHistosSigma(const lambda1405candidate& lambda1405Cand, const TCand& sigmaCand, const TTrack& kinkDauTrack) + { if (sigmaCand.mothSign() > 0) { rSigmaPlus.fill(HIST("hSigmaPlusMass"), sigmaCand.mSigmaPlus()); @@ -606,18 +607,18 @@ struct lambda1405analysis { auto kinkDauMom = std::array{sigmaCand.pxDaug(), sigmaCand.pyDaug(), sigmaCand.pzDaug()}; auto sigmaMom = std::array{sigmaCand.pxMoth(), sigmaCand.pyMoth(), sigmaCand.pzMoth()}; // Sigma properties - lambda1405Cand.sigmaId = sigmaCand.globalIndex(); + lambda1405Cand.sigmaId = sigmaCand.globalIndex(); lambda1405Cand.sigmaMinusMass = sigmaCand.mSigmaMinus(); - lambda1405Cand.sigmaPlusMass = sigmaCand.mSigmaPlus(); - lambda1405Cand.xiMinusMass = sigmaCand.mXiMinus(); - lambda1405Cand.sigmaSign = sigmaCand.mothSign(); - lambda1405Cand.sigmaAlphaAP = alphaAP(sigmaMom, kinkDauMom); - lambda1405Cand.sigmaQtAP = qtAP(sigmaMom, kinkDauMom); - lambda1405Cand.sigmaPt = sigmaCand.ptMoth(); - lambda1405Cand.sigmaRadius = sigmaRad; + lambda1405Cand.sigmaPlusMass = sigmaCand.mSigmaPlus(); + lambda1405Cand.xiMinusMass = sigmaCand.mXiMinus(); + lambda1405Cand.sigmaSign = sigmaCand.mothSign(); + lambda1405Cand.sigmaAlphaAP = alphaAP(sigmaMom, kinkDauMom); + lambda1405Cand.sigmaQtAP = qtAP(sigmaMom, kinkDauMom); + lambda1405Cand.sigmaPt = sigmaCand.ptMoth(); + lambda1405Cand.sigmaRadius = sigmaRad; lambda1405Cand.kinkDcaDauToPv = sigmaCand.dcaDaugPv(); - if (lambda1405Cand.sigmaQtAP < funcMinQtAlphaAP->Eval(lambda1405Cand.sigmaAlphaAP) || + if (lambda1405Cand.sigmaQtAP < funcMinQtAlphaAP->Eval(lambda1405Cand.sigmaAlphaAP) || lambda1405Cand.sigmaQtAP > funcMaxQtAlphaAP->Eval(lambda1405Cand.sigmaAlphaAP)) { return; } @@ -628,8 +629,8 @@ struct lambda1405analysis { } // Kink daughter properties - lambda1405Cand.kinkDauId = kinkDauTrack.globalIndex(); - lambda1405Cand.kinkPt = kinkDauTrack.pt(); + lambda1405Cand.kinkDauId = kinkDauTrack.globalIndex(); + lambda1405Cand.kinkPt = kinkDauTrack.pt(); lambda1405Cand.kinkPiNSigTpc = kinkDauTrack.tpcNSigmaPi(); lambda1405Cand.kinkPiNSigTof = kinkDauTrack.tofNSigmaPi(); lambda1405Cand.kinkPrNSigTpc = kinkDauTrack.tpcNSigmaPr(); @@ -652,8 +653,8 @@ struct lambda1405analysis { continue; } rSelections.fill(HIST("hSelectionsBachPi"), 5); // PID sel - rSelections.fill(HIST("hSelectionsL1405"), 2); // Bach Pi selection - + rSelections.fill(HIST("hSelectionsL1405"), 2); // Bach Pi selection + auto piMom = std::array{piTrack.px(), piTrack.py(), piTrack.pz()}; float invMass{-1.f}; if (lambda1405Cand.isSigmaMinus) { @@ -664,7 +665,7 @@ struct lambda1405analysis { if (invMass > cutUpperMass) { continue; } - rSelections.fill(HIST("hSelectionsL1405"), 3); // Upper mass selection + rSelections.fill(HIST("hSelectionsL1405"), 3); // Upper mass selection // Daughter Pi properties lambda1405Cand.piId = piTrack.globalIndex(); @@ -686,20 +687,20 @@ struct lambda1405analysis { lambda1405Cand.scalarProd = -1; // Check correlations between transverse momenta of L1405, sigma and bachelor pi - if (std::hypot(sigmaCand.pxMoth(), sigmaCand.pyMoth()) < funcMinSigmaPtVsL1405Pt->Eval(lambda1405Cand.pt()) || + if (std::hypot(sigmaCand.pxMoth(), sigmaCand.pyMoth()) < funcMinSigmaPtVsL1405Pt->Eval(lambda1405Cand.pt()) || std::hypot(sigmaCand.pxMoth(), sigmaCand.pyMoth()) > funcMaxSigmaPtVsL1405Pt->Eval(lambda1405Cand.pt())) { continue; } - rSelections.fill(HIST("hSelectionsL1405"), 4); // Accepted - - if (piTrack.pt() < funcMinBachPiPtVsL1405Pt->Eval(lambda1405Cand.pt()) || - piTrack.pt() > funcMaxBachPiPtVsL1405Pt->Eval(lambda1405Cand.pt())) { + rSelections.fill(HIST("hSelectionsL1405"), 4); // Accepted + + if (piTrack.pt() < funcMinBachPiPtVsL1405Pt->Eval(lambda1405Cand.pt()) || + piTrack.pt() > funcMaxBachPiPtVsL1405Pt->Eval(lambda1405Cand.pt())) { continue; } - rSelections.fill(HIST("hSelectionsL1405"), 5); // Accepted + rSelections.fill(HIST("hSelectionsL1405"), 5); // Accepted fillHistosLambda1405(lambda1405Cand, piTrack); - rSelections.fill(HIST("hSelectionsL1405"), 6); // Accepted + rSelections.fill(HIST("hSelectionsL1405"), 6); // Accepted selectedCandidates.push_back(lambda1405Cand); } } @@ -724,9 +725,10 @@ struct lambda1405analysis { return isKinkFromSigma; // Return true if the kink comes from the Sigma } - template - void fillOutputData(const TCollision& collision, const TCand& sigmaCands, const TTrack& kinkDauTrack) { - if constexpr (requires{collision.centFT0C();}) { + template + void fillOutputData(const TCollision& collision, const TCand& sigmaCands, const TTrack& kinkDauTrack) + { + if constexpr (requires { collision.centFT0C(); }) { if (collision.centFT0C() < centralityMin || collision.centFT0C() > centralityMax) { LOG(info) << "Skipping collision due to centrality cut"; return; @@ -771,7 +773,7 @@ struct lambda1405analysis { lambda1405Cand.kinkDcaDauToPv, lambda1405Cand.bachPiNSigTpc, lambda1405Cand.bachPiNSigTof); } - if constexpr (requires{collision.qvecFT0CRe();}) { + if constexpr (requires { collision.qvecFT0CRe(); }) { float const xQVec = collision.qvecFT0CRe(); float const yQVec = collision.qvecFT0CIm(); float const cos2Phi = std::cos(2 * lambda1405Cand.phi); @@ -813,8 +815,9 @@ struct lambda1405analysis { } PROCESS_SWITCH(lambda1405analysis, processDataWCentQVecs, "Data processing with centrality and Q vectors info", false); - template - int matchGenDecay(const TMother& motherPart, const aod::McParticles& mcParticles, std::array& daugsIdxs) { + template + int matchGenDecay(const TMother& motherPart, const aod::McParticles& mcParticles, std::array& daugsIdxs) + { int pdgMother = motherPart.pdgCode(); int8_t sign = 0; @@ -825,7 +828,7 @@ struct lambda1405analysis { int decayChannel = -1; // Match L(1405) --> n pi- pi+ final state - if (pdgMother > 0) { // Change sign of neutral decay products + if (pdgMother > 0) { // Change sign of neutral decay products finalState = {PDG_t::kNeutron, PDG_t::kPiMinus, PDG_t::kPiPlus}; } else { finalState = {-PDG_t::kNeutron, PDG_t::kPiMinus, PDG_t::kPiPlus}; @@ -843,14 +846,14 @@ struct lambda1405analysis { decayChannel = kSigmaMinusPiToPiPiNeutron; daugsIdxs[1] = arrL1405Daugs[iProng]; } - if (std::abs(daughI.pdgCode()) == PDG_t::kPiPlus) { + if (std::abs(daughI.pdgCode()) == PDG_t::kPiPlus) { daugsIdxs[0] = arrL1405Daugs[iProng]; } } } // Match L(1405) --> p pi0 pi+ final state, only possible for Sigma+ - if (pdgMother > 0) { // Change sign of neutral decay products + if (pdgMother > 0) { // Change sign of neutral decay products finalState = {PDG_t::kProton, PDG_t::kPi0, PDG_t::kPiMinus}; } else { finalState = {PDG_t::kProton, -PDG_t::kPi0, PDG_t::kPiMinus}; @@ -884,10 +887,11 @@ struct lambda1405analysis { return decayChannel; } - template - void fillOutputMc(const TCollision& collisions, const aod::KinkCands& sigmaCands, const aod::McTrackLabels& trackLabelsMC, const TTrack& tracks, const aod::McParticles& particlesMC) { + template + void fillOutputMc(const TCollision& collisions, const aod::KinkCands& sigmaCands, const aod::McTrackLabels& trackLabelsMC, const TTrack& tracks, const aod::McParticles& particlesMC) + { for (const auto& collision : collisions) { - if constexpr (requires{collision.centFT0C();}) { + if constexpr (requires { collision.centFT0C(); }) { if (collision.centFT0C() < centralityMin || collision.centFT0C() > centralityMax) { return; } @@ -903,7 +907,7 @@ struct lambda1405analysis { constructCollCandidates(sigmaCand, tracksPerCol, selectedCandidates); LOG(info) << "Selected " << selectedCandidates.size() << " Lambda(1405) candidates in this collision"; for (const auto& lambda1405Cand : selectedCandidates) { - rLambda1405.fill(HIST("hRecoL1405"), 0., lambda1405Cand.pt()); // All reconstructed + rLambda1405.fill(HIST("hRecoL1405"), 0., lambda1405Cand.pt()); // All reconstructed // Do MC association auto mcLabPiKink = trackLabelsMC.rawIteratorAt(lambda1405Cand.kinkDauId); @@ -912,7 +916,7 @@ struct lambda1405analysis { if (!mcLabSigma.has_mcParticle() || !mcLabPiKink.has_mcParticle() || !mcLabPi.has_mcParticle()) { continue; // Skip if no valid MC association } - rLambda1405.fill(HIST("hRecoL1405"), 1., lambda1405Cand.pt()); // All with associated MC particle + rLambda1405.fill(HIST("hRecoL1405"), 1., lambda1405Cand.pt()); // All with associated MC particle auto mcTrackKink = mcLabPiKink.mcParticle_as(); auto mcTrackSigma = mcLabSigma.mcParticle_as(); @@ -925,17 +929,17 @@ struct lambda1405analysis { if (!isSigmaMinusKink && !isSigmaPlusToPiKink && !isSigmaPlusToPrKink) { continue; // Skip if not a valid Sigma kink decay } - rLambda1405.fill(HIST("hRecoL1405"), 2., lambda1405Cand.pt()); // Has kink decay in MC + rLambda1405.fill(HIST("hRecoL1405"), 2., lambda1405Cand.pt()); // Has kink decay in MC if (std::abs(mcTrackPi.pdgCode()) != 211) { continue; // Skip if not a valid pion candidate } - rLambda1405.fill(HIST("hRecoL1405"), 3., lambda1405Cand.pt()); // Has bach pi + rLambda1405.fill(HIST("hRecoL1405"), 3., lambda1405Cand.pt()); // Has bach pi if (!mcTrackSigma.has_mothers() || !mcTrackPi.has_mothers()) { continue; // Skip if no mothers found } - rLambda1405.fill(HIST("hRecoL1405"), 4., lambda1405Cand.pt()); // Has mothers for Sigma and Pi + rLambda1405.fill(HIST("hRecoL1405"), 4., lambda1405Cand.pt()); // Has mothers for Sigma and Pi // check that labpi and labsigma have the same mother (a lambda1405 candidate) int lambda1405Id = -1; @@ -950,7 +954,7 @@ struct lambda1405analysis { if (lambda1405Id == -1) { continue; // Skip if the Sigma and pion do not share the same lambda1405 candidate } - rLambda1405.fill(HIST("hRecoL1405"), 4., lambda1405Cand.pt()); // Has same mother + rLambda1405.fill(HIST("hRecoL1405"), 4., lambda1405Cand.pt()); // Has same mother auto lambda1405Mother = particlesMC.rawIteratorAt(lambda1405Id); float lambda1405Mass = std::sqrt(lambda1405Mother.e() * lambda1405Mother.e() - lambda1405Mother.p() * lambda1405Mother.p()); @@ -998,7 +1002,7 @@ struct lambda1405analysis { auto sigmaKinkDaug = particlesMC.rawIteratorAt(idxsSigmaKinkBachPi[2]); // Generated Armenteros-Podolanski variables float genSigmaAlphaAP = alphaAP({sigmaDaug.px(), sigmaDaug.py(), sigmaDaug.pz()}, {sigmaKinkDaug.px(), sigmaKinkDaug.py(), sigmaKinkDaug.pz()}); - float genSigmaQtAP = qtAP({sigmaDaug.px(), sigmaDaug.py(), sigmaDaug.pz()}, {sigmaKinkDaug.px(), sigmaKinkDaug.py(), sigmaKinkDaug.pz()}); + float genSigmaQtAP = qtAP({sigmaDaug.px(), sigmaDaug.py(), sigmaDaug.pz()}, {sigmaKinkDaug.px(), sigmaKinkDaug.py(), sigmaKinkDaug.pz()}); if (decayChannel == kSigmaMinusPiToPiPiNeutron) { rLambda1405.fill(HIST("h2GenSigmaMinusArmPod"), genSigmaAlphaAP, genSigmaQtAP); rLambda1405.fill(HIST("h2GenPtVsBachPtSigmaMinusPiToPiPiNeutron"), mcPart.pt(), bachPi.pt()); From c8b6d61a8359695b83e08b34205d95aa513080f4 Mon Sep 17 00:00:00 2001 From: marcellocosti Date: Tue, 2 Jun 2026 13:48:53 +0200 Subject: [PATCH 5/5] Remove unused pars --- PWGLF/TableProducer/Common/kinkBuilder.cxx | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/PWGLF/TableProducer/Common/kinkBuilder.cxx b/PWGLF/TableProducer/Common/kinkBuilder.cxx index 496f986d158..48cdcdf1f7f 100644 --- a/PWGLF/TableProducer/Common/kinkBuilder.cxx +++ b/PWGLF/TableProducer/Common/kinkBuilder.cxx @@ -454,7 +454,7 @@ struct kinkBuilder { } template - void buildKinkCand(const TSVCand& svCand, const TTracks& tracks, bool& isAccepted, const TColls& collisions) + void buildKinkCand(const TSVCand& svCand, const TTracks& tracks, bool& isAccepted, const TColls&) { isAccepted = false; kinkCandidate kinkCand; @@ -828,8 +828,7 @@ struct kinkBuilder { } } - void processMc(aod::McCollisions const& mcGenCollisions, - McRecoCollisions const& mcRecoCollisions, + void processMc(McRecoCollisions const& mcRecoCollisions, aod::McParticles const& mcParticles, TracksFullMc const& tracksMc, aod::AmbiguousTracks const& ambiTracksMc, @@ -839,8 +838,7 @@ struct kinkBuilder { } PROCESS_SWITCH(kinkBuilder, processMc, "MC processing", false); - void processMcWCent(aod::McCollisions const& mcGenCollisions, - McRecoCollisionsCentSel const& mcRecoCollisions, + void processMcWCent(McRecoCollisionsCentSel const& mcRecoCollisions, aod::McParticles const& mcParticles, TracksFullMc const& tracksMc, aod::AmbiguousTracks const& ambiTracksMc,