diff --git a/PWGLF/DataModel/LFLambda1405Table.h b/PWGLF/DataModel/LFLambda1405Table.h index be61dc083bc..31721f634e3 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..48cdcdf1f7f 100644 --- a/PWGLF/TableProducer/Common/kinkBuilder.cxx +++ b/PWGLF/TableProducer/Common/kinkBuilder.cxx @@ -16,7 +16,9 @@ #include "PWGLF/DataModel/LFKinkDecayTables.h" #include "PWGLF/Utils/svPoolCreator.h" +#include "Common/Core/RecoDecay.h" #include "Common/Core/trackUtilities.h" +#include "Common/DataModel/Centrality.h" #include #include @@ -40,6 +42,7 @@ #include #include +#include #include #include @@ -50,11 +53,22 @@ #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 CollisionsCent = soa::Join; +using McRecoCollisions = soa::Join; +using McRecoCollisionsCent = soa::Join; + +enum KinkDecayType { kSigmaMinusToPiMinusNeutron = 0, + kSigmaPlusToPiPlusNeutron, + kSigmaPlusToProtonPi0, + kNMatchedDecays }; namespace { @@ -62,12 +76,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,18 +168,32 @@ 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"}; + + // 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"}; 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"}; + + // 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; @@ -161,6 +206,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 +215,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 +256,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}"}; @@ -221,213 +269,355 @@ 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}); + 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", "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", "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", "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", "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", "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", "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", "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", "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 || doprocessMcWCent) { + if (skipBkgCands) { + 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->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", "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}}); } - 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); - for (const auto& svCand : kinkPool) { - kinkCandidate kinkCand; + if (candidate.itsNCls() >= 4) + return false; + hSelDaugQA->Fill(5.f, isPositive); - auto trackMoth = tracks.rawIteratorAt(svCand.tr0Idx); - auto trackDaug = tracks.rawIteratorAt(svCand.tr1Idx); + if (candidate.tpcNClsCrossedRows() <= 0.8 * candidate.tpcNClsFindable()) + return false; + hSelDaugQA->Fill(6.f, isPositive); - auto const& collision = trackMoth.template collision_as(); - auto const& bc = collision.template bc_as(); - initCCDB(bc); + if (candidate.tpcNClsFound() <= nTPCClusMinDaug) + return false; + hSelDaugQA->Fill(7.f, isPositive); - 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()}; + return true; + } - o2::track::TrackParCov trackParCovMoth = getTrackParCov(trackMoth); - o2::track::TrackParCov trackParCovMothPV{trackParCovMoth}; - o2::base::Propagator::Instance()->PropagateToXBxByBz(trackParCovMoth, LayerRadii[trackMoth.itsNCls() - 1]); + template + void buildKinkCand(const TSVCand& svCand, const TTracks& tracks, bool& isAccepted, const TColls&) + { + isAccepted = false; + kinkCandidate kinkCand; - std::array dcaInfoMoth; - o2::base::Propagator::Instance()->propagateToDCABxByBz({primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, trackParCovMothPV, 2.f, static_cast(cfgMaterialCorrection.value), &dcaInfoMoth); + auto trackMoth = tracks.rawIteratorAt(svCand.tr0Idx); + auto trackDaug = tracks.rawIteratorAt(svCand.tr1Idx); - if (std::abs(dcaInfoMoth[0]) > maxDCAMothToPV) { - continue; - } + // 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 trackParCovDaug = getTrackParCov(trackDaug); + auto const& collision = trackMoth.template collision_as(); + auto const& bc = collision.template bc_as(); + initCCDB(bc); - // check if the kink daughter is close to the mother - if (std::abs(trackParCovMoth.getZ() - trackParCovDaug.getZ()) > maxZDiff) { - continue; - } + 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(trackParCovMoth.getPhi() - trackParCovDaug.getPhi()) * radToDeg) > maxPhiDiff) { - continue; - } + o2::track::TrackParCov trackParCovMoth = getTrackParCov(trackMoth); + o2::track::TrackParCov trackParCovMothPV{trackParCovMoth}; + o2::base::Propagator::Instance()->PropagateToXBxByBz(trackParCovMoth, LayerRadii[trackMoth.itsNCls() - 1]); - // 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; - } + std::array dcaInfoMoth; + o2::base::Propagator::Instance()->propagateToDCABxByBz({primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, trackParCovMothPV, 2.f, static_cast(cfgMaterialCorrection.value), &dcaInfoMoth); - if (updateMothTrackUsePV) { - // update the mother track parameters using the primary vertex - trackParCovMoth = trackParCovMothPV; - if (!trackParCovMoth.update(primaryVertex)) { - continue; - } - } + hDCAMothToPV->Fill(std::abs(dcaInfoMoth[0]), chargeCombSvCand); + if (std::abs(dcaInfoMoth[0]) > maxDCAMothToPV) { + return; + } + hSelKinkedTrackQA->Fill(1.f, chargeCombSvCand); // MaxDCAMothToPV 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; - } + o2::track::TrackParCov trackParCovDaug = getTrackParCov(trackDaug); - if (!fitter.propagateTracksToVertex()) { - 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 - auto propMothTrack = fitter.getTrack(0); - auto propDaugTrack = fitter.getTrack(1); - kinkCand.decVtx = fitter.getPCACandidatePos(); + 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 - // 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; - } + // 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 - // 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 (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 - for (int i = 0; i < 7; i++) { - if (trackDaug.itsClusterMap() & (1 << i)) { - firstLayerDaug = i; - break; - } - } + 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 - if (doSVRadiusCut && lastLayerMoth >= firstLayerDaug) { - continue; - } + if (!fitter.propagateTracksToVertex()) { + return; + } + hSelKinkedTrackQA->Fill(7.f, chargeCombSvCand); // PropagateTracksToVertex cut bin - if (doSVRadiusCut && decRad2 < LayerRadii[lastLayerMoth] * LayerRadii[lastLayerMoth]) { - continue; - } + auto propMothTrack = fitter.getTrack(0); + auto propDaugTrack = fitter.getTrack(1); + kinkCand.decVtx = fitter.getPCACandidatePos(); - for (int i = 0; i < 3; i++) { - kinkCand.decVtx[i] -= kinkCand.primVtx[i]; - } + // 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 - propMothTrack.getPxPyPzGlo(kinkCand.momMoth); - propDaugTrack.getPxPyPzGlo(kinkCand.momDaug); - for (int i = 0; i < 3; i++) { - kinkCand.momMoth[i] *= charge; - kinkCand.momDaug[i] *= charge; + // 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; } - 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]; + } + + for (int i = 0; i < 7; i++) { + if (trackDaug.itsClusterMap() & (1 << i)) { + firstLayerDaug = i; + break; } + } + + if (doSVRadiusCut && lastLayerMoth >= firstLayerDaug) { + return; + } + hSelKinkedTrackQA->Fill(9.f, chargeCombSvCand); // MothDaugLastLayers 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); + if (doSVRadiusCut && decRad2 < LayerRadii[lastLayerMoth] * LayerRadii[lastLayerMoth]) { + return; + } + hSelKinkedTrackQA->Fill(10.f, chargeCombSvCand); // MothLastLayerRadiusCheck cut bin + + 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 +641,153 @@ 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); + } + } + + template + void fillOutputsData(const TColls& collisions, const TTracks& tracks, const aod::AmbiguousTracks& ambiTracks, const aod::BCs& bcs) + { + kinkCandidates.clear(); + + buildSvPool(collisions, tracks, ambiTracks, bcs); + auto& kinkPool = svCreator.getSVCandPool(collisions, !unlikeSignBkg); + 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); + } + } + } + + 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) + { + 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 (RecoDecay::isMatchedMCGen(mcParticles, motherPart, pdgMother, finState, true, &sign, DepthMcMatchMax)) { + 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 (RecoDecay::isMatchedMCGen(mcParticles, motherPart, pdgMother, finState, true, &sign, DepthMcMatchMax)) { + 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 (RecoDecay::isMatchedMCGen(mcParticles, motherPart, pdgMother, finState, true, &sign, DepthMcMatchMax)) { + return kSigmaPlusToProtonPi0; + } + break; + } + default: + LOG(warning) << "No matching function implemented for the selected mother particle hypothesis. Returning -1."; + return -1; + } + + return -1; + } + template + void fillOutputsMc(const TColls& mcRecoCollisions, const TTracks& tracksMc, const aod::AmbiguousTracks& ambiTracksMc, const aod::BCs& bcs, const aod::McParticles& mcParticles) + { kinkCandidates.clear(); - fillCandidateData(collisions, tracks, ambiTracks, bcs); + + buildSvPool(mcRecoCollisions, tracksMc, ambiTracksMc, bcs); + auto& kinkPool = svCreator.getSVCandPool(mcRecoCollisions, !unlikeSignBkg); + for (const auto& svCand : kinkPool) { + // Perform matching of the kink candidate + auto trackMoth = tracksMc.rawIteratorAt(svCand.tr0Idx); + auto genMothPart = trackMoth.template mcParticle_as(); + auto trackDaug = tracksMc.rawIteratorAt(svCand.tr1Idx); + auto genDaugPart = trackDaug.template 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 + } + decayChannel = matchKinkDecay(genMothPart, mcParticles); + 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 +806,47 @@ 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); + } + } + } + } + + void processMc(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(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 7c12c2b9042..dd4b51076b7 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,29 @@ using namespace o2::framework::expressions; using TracksFull = soa::Join; using CollisionsFull = soa::Join; -using CollisionsFullMC = soa::Join; +using CollisionsFullWCentQVecs = soa::Join; +using CollisionsFullMc = soa::Join; +using CollisionsFullMcWCent = 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 +76,99 @@ 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", "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)"}; + 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"}; + + 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.}, ""}; + 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,41 +176,184 @@ 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}}); - - if (doprocessMC) { + // 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, {{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, "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"); + 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.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"); + 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.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 || doprocessMcWCentSel) { // 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 #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"); + 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("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 - 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}}); } + + // 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) @@ -164,139 +372,337 @@ 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; } - if (prFromSigma) { - return true; + rSelections.fill(HIST("hSelectionsKinkPr"), 1); // Nsigma Tpc + + if (candidate.tpcNClsFound() < cutNTpcClusPi) { + return false; + } + 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 (isPiKink) { - rLambda1405.fill(HIST("h2PtMassSigmaBeforeCuts_0"), sigmaCand.mothSign() * sigmaCand.ptMoth(), sigmaCand.mSigmaMinus()); - rLambda1405.fill(HIST("h2PtPiNSigma_0"), sigmaCand.mothSign() * kinkDauTrack.pt(), kinkDauTrack.tpcNSigmaPi()); + if (sigmaCand.mothSign() < 0) { + rSelections.fill(HIST("hSelectionsSigmaMinus"), 1); // Passed kink sel + } else { + rSelections.fill(HIST("hSelectionsSigmaPlus"), 1); // Passed kink sel } - if (isProKink) { - rLambda1405.fill(HIST("h2PtMassSigmaBeforeCuts_1"), sigmaCand.mothSign() * sigmaCand.ptMoth(), sigmaCand.mSigmaPlus()); - rLambda1405.fill(HIST("h2PtPiNSigma_1"), sigmaCand.mothSign() * kinkDauTrack.pt(), kinkDauTrack.tpcNSigmaPr()); + + // 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 (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 < funcMinQtAlphaAP->Eval(lambda1405Cand.sigmaAlphaAP) || + lambda1405Cand.sigmaQtAP > funcMaxQtAlphaAP->Eval(lambda1405Cand.sigmaAlphaAP)) { + return; } + if (sigmaCand.mothSign() < 0) { + rSelections.fill(HIST("hSelectionsSigmaMinus"), 6); // Passed mass sel + } else { + rSelections.fill(HIST("hSelectionsSigmaPlus"), 6); // 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 + + bool isUnlikeSign = (piTrack.sign() != sigmaCand.mothSign()); + bool acceptPair = doLSBkg ? !isUnlikeSign : isUnlikeSign; + if (!acceptPair) { + continue; } + rSelections.fill(HIST("hSelectionsBachPi"), 1); - if (!selectPiTrack(piTrack, false)) { + if (!selectPiBach(piTrack)) { continue; } + rSelections.fill(HIST("hSelectionsBachPi"), 5); // PID sel + rSelections.fill(HIST("hSelectionsL1405"), 2); // Bach Pi selection - auto kinkDauMom = std::array{sigmaCand.pxDaug(), sigmaCand.pyDaug(), sigmaCand.pzDaug()}; - auto sigmaMom = std::array{sigmaCand.pxMoth(), sigmaCand.pyMoth(), sigmaCand.pzMoth()}; 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.phi = std::atan2(lambda1405Cand.py, lambda1405Cand.px); + lambda1405Cand.scalarProd = -1; - 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 + // 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; } - return true; // Candidate is selected + 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); } - return false; // No valid pion track found } template @@ -319,60 +725,199 @@ 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) { - if (selectCandidate(sigmaCand, tracks)) { + for (const auto& sigmaCand : sigmaCands) { + std::vector selectedCandidates; + constructCollCandidates(sigmaCand, kinkDauTrack, selectedCandidates); + for (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); + } + 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 processMC(CollisionsFullMC const& collisions, aod::KinkCands const& kinkCands, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const& particlesMC, TracksFull const& tracks) + void processDataWCentQVecs(CollisionsCentSel::iterator const& collision, aod::KinkCands const& kinkCands, TracksFull const& tracks) + { + 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, std::array& daugsIdxs) + { + int pdgMother = motherPart.pdgCode(); + int8_t sign = 0; + + int MaxDepthFinState = 2; // Maximum depth to look for the decay chain + + 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)) { + // Match the intermediate Sigma + 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 = 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]; + } + } + } + + // 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 + 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]; + } + } + } + + 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; + } + + 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 (std::abs(collision.posZ()) > cutzvertex || !collision.sel8()) { + 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()); 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); - if (!mcLabSigma.has_mcParticle() || mcLabPiKink.has_mcParticle() || mcLabPi.has_mcParticle()) { + 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()) { 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(); @@ -384,14 +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 if (std::abs(mcTrackPi.pdgCode()) != 211) { 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()) { 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; @@ -406,35 +954,31 @@ 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 + 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,30 +990,51 @@ struct lambda1405analysis { if (std::abs(mcPart.pdgCode()) != lambda1405PdgCode) { continue; // Only consider Lambda(1405) candidates } - - if (!mcPart.has_daughters()) { - continue; // Skip if no daughters + 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()); - // 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 - } + 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 (!hasSigmaDaughter) { - continue; // Skip if no Sigma daughter found + 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()); } - - 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); } } - 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) { 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;