From 17a568724fb11c1168b94bb5f7e57c6b0865686d Mon Sep 17 00:00:00 2001 From: sandeep dudi Date: Fri, 21 Nov 2025 11:05:03 +0100 Subject: [PATCH] Improvement in Kink vertex finder --- PWGLF/Tasks/Nuspex/spectraKinkPiKa.cxx | 143 ++++++++++++++++--------- 1 file changed, 94 insertions(+), 49 deletions(-) diff --git a/PWGLF/Tasks/Nuspex/spectraKinkPiKa.cxx b/PWGLF/Tasks/Nuspex/spectraKinkPiKa.cxx index a6ed41677e7..aa96ef670d4 100644 --- a/PWGLF/Tasks/Nuspex/spectraKinkPiKa.cxx +++ b/PWGLF/Tasks/Nuspex/spectraKinkPiKa.cxx @@ -67,7 +67,7 @@ using CollisionsFull = soa::Join; namespace { -constexpr std::array LayerRadii{2.33959f, 3.14076f, 3.91924f, 19.6213f, 24.5597f, 34.388f, 39.3329f}; +// constexpr std::array LayerRadii{2.33959f, 3.14076f, 3.91924f, 19.6213f, 24.5597f, 34.388f, 39.3329f}; 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"}; @@ -96,7 +96,7 @@ struct KinkBuilder { Service ccdb; // Selection criteria Configurable maxDCAMothToPV{"maxDCAMothToPV", 0.2, "Max DCA of the mother to the PV"}; - Configurable minDCADaugToPV{"minDCADaugToPV", 0., "Min DCA of the daughter to the PV"}; + Configurable minDCADaugToPV{"minDCADaugToPV", 0.1, "Min DCA of the daughter to the PV"}; Configurable minPtMoth{"minPtMoth", 0.15, "Minimum pT of the hypercandidate"}; Configurable maxZDiff{"maxZDiff", 20., "Max z difference between the kink daughter and the mother"}; Configurable maxPhiDiff{"maxPhiDiff", 100, "Max phi difference between the kink daughter and the mother"}; @@ -107,6 +107,10 @@ struct KinkBuilder { Configurable itsChi2cut{"itsChi2cut", 36, "mother itsChi2 cut"}; Configurable askTOFforDaug{"askTOFforDaug", false, "If true, ask for TOF signal"}; Configurable kaontopologhy{"kaontopologhy", true, "If true, selected mother have both ITS+TPC "}; + Configurable vertexfinding{"vertexfinding", false, "If true, find the vextex in TPC and applied cut of z and phi"}; + Configurable rMin{"rMin", 120., "min radius for kink vertex"}; + Configurable rMax{"rMax", 200., "max radius for kink vertex"}; + Configurable rStep{"rStep", 2., "step size for scan radius in tpc"}; o2::vertexing::DCAFitterN<2> fitter; o2::base::MatLayerCylSet* lut = nullptr; @@ -214,10 +218,15 @@ struct KinkBuilder { svCreator.fillBC2Coll(collisions, bcs); for (const auto& track : tracks) { bool isDaug = selectDaugTrack(track); + bool isMoth = selectMothTrack(track); if (!isDaug && !isMoth) continue; + if (isDaug && track.isPVContributor()) + continue; + if (isMoth && !track.isPVContributor()) + continue; if (isDaug && std::abs(track.eta()) > etaMaxDaug) continue; if (isMoth && std::abs(track.eta()) > etaMaxMoth) @@ -245,25 +254,65 @@ struct KinkBuilder { o2::track::TrackParCov trackParCovMoth = getTrackParCov(trackMoth); o2::track::TrackParCov trackParCovMothPV{trackParCovMoth}; - o2::base::Propagator::Instance()->PropagateToXBxByBz(trackParCovMoth, LayerRadii[trackMoth.itsNCls() - 1]); std::array dcaInfoMoth; o2::base::Propagator::Instance()->propagateToDCABxByBz({primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, trackParCovMothPV, 2.f, static_cast(cfgMaterialCorrection.value), &dcaInfoMoth); - o2::track::TrackParCov trackParCovDaug = getTrackParCov(trackDaug); + // propagate to PV + std::array dcaInfoDaug; + o2::base::Propagator::Instance()->propagateToDCABxByBz({primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, trackParCovDaug, 2.f, static_cast(cfgMaterialCorrection.value), &dcaInfoDaug); - // check if the kink daughter is close to the mother - if (std::abs(trackParCovMoth.getZ() - trackParCovDaug.getZ()) > maxZDiff) { + if (std::abs(dcaInfoMoth[1]) > maxDCAMothToPV) { continue; } - if ((std::abs(trackParCovMoth.getPhi() - trackParCovDaug.getPhi()) * radToDeg) > maxPhiDiff) { + if (std::abs(dcaInfoDaug[1]) < minDCADaugToPV) { 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); - if (std::abs(dcaInfoDaug[0]) < minDCADaugToPV) { + + if (vertexfinding) { + float bestR = -1; + float bestDeltaPhi = 999; + float bestDeltaZ = 999; + float bestCost = 1e12; + // make local copies (don’t modify originals) + auto mothTmp0 = trackParCovMoth; + auto daugTmp0 = trackParCovDaug; + const float minr = rMin; + const float maxr = rMax; + const float rs = rStep; + for (float R = minr; R <= maxr; R += rs) { + auto mothTmp = mothTmp0; + auto daugTmp = daugTmp0; + if (!o2::base::Propagator::Instance()->PropagateToXBxByBz(mothTmp, R)) + continue; + if (!o2::base::Propagator::Instance()->PropagateToXBxByBz(daugTmp, R)) + continue; + float dphi = std::abs(mothTmp.getPhi() - daugTmp.getPhi()); + if (dphi > M_PI) + dphi = 2 * M_PI - dphi; // wrap φ + float dz = std::abs(mothTmp.getZ() - daugTmp.getZ()); + float cost = dphi * dphi + dz * dz; // <-- correct metric + if (cost < bestCost) { + bestCost = cost; + bestDeltaPhi = dphi; + bestDeltaZ = dz; + bestR = R; + } + } + if (bestR < 0) + continue; + if (bestDeltaZ > maxZDiff) + continue; + if (bestDeltaPhi * radToDeg > maxPhiDiff) + continue; + } + /* // check if the kink daughter is close to the mother + if (std::abs(trackParCovMoth.getZ() - trackParCovDaug.getZ()) > maxZDiff) { + continue; + } + if ((std::abs(trackParCovMoth.getPhi() - trackParCovDaug.getPhi()) * radToDeg) > maxPhiDiff) { continue; } + */ int nCand = 0; try { nCand = fitter.process(trackParCovMoth, trackParCovDaug); @@ -367,7 +416,7 @@ struct SpectraKinkPiKa { Configurable tpcChi2Cut{"tpcChi2Cut", 4.0, "tpcChi2Cut"}; Configurable minqt{"minqt", 0.12, "min qt for kaon"}; Configurable maxqt{"maxqt", 0.3, "max qt for kaon"}; - + Configurable minPtMothmc{"minPtMothmc", 0.15, "Minimum pT of the mother"}; Configurable centestimator{"centestimator", 0, "Select multiplicity estimator: 0 - FT0C, 1 - FT0A, 2 - FT0M, 3 - FV0A, 4 - PVTracks"}; Configurable pid{"pid", 321, ""}; Configurable dpid{"dpid", 13, ""}; @@ -400,14 +449,13 @@ struct SpectraKinkPiKa { // Event selection rEventSelection.add("hVertexZRec", "hVertexZRec", {HistType::kTH1F, {{200, -20.0, 20.0}}}); rEventSelection.add("hMultiplicity", "hMultiplicity", {HistType::kTH1F, {multAxis}}); - - rpiKkink.add("h2_dau_pt_vs_eta_rec", "pt_vs_eta_dau", {HistType::kTH3F, {ptAxis, etaAxis, multAxis}}); - rpiKkink.add("h2_moth_pt_vs_eta_rec", "pt_vs_eta_moth", {HistType::kTH3F, {ptAxis, etaAxis, multAxis}}); - rpiKkink.add("h2_pt_moth_vs_dau_rec", "pt_moth_vs_dau", {HistType::kTH2F, {ptAxis, ptAxis}}); - - rpiKkink.add("h2_qt_vs_pt", "qt_pt", {HistType::kTH2F, {qtAxis, ptAxis}}); + if (additionalhist) { + rpiKkink.add("h2_dau_pt_vs_eta_rec", "pt_vs_eta_dau", {HistType::kTH3F, {ptAxis, etaAxis, multAxis}}); + rpiKkink.add("h2_moth_pt_vs_eta_rec", "pt_vs_eta_moth", {HistType::kTH3F, {ptAxis, etaAxis, multAxis}}); + rpiKkink.add("h2_pt_moth_vs_dau_rec", "pt_moth_vs_dau", {HistType::kTH2F, {ptAxis, ptAxis}}); + rpiKkink.add("h2_qt_vs_pt", "qt_pt", {HistType::kTH2F, {qtAxis, ptAxis}}); + } rpiKkink.add("h2_kink_angle", "kink angle", {HistType::kTH2F, {ptAxis, kinkAxis}}); - // inv mass rpiKkink.add("h2_kaon_data", "h2_kaon_data", HistType::kTHnSparseF, {massAxis, ptAxis, qtAxis, multAxis}, true); rpiKkink.add("h1_tracks_data", "track_cut_data", {HistType::kTH1F, {{17, 0.5, 17.5}}}); @@ -454,19 +502,19 @@ struct SpectraKinkPiKa { } if (doprocessMC) { - rpiKkink.add("h2_qt", "qt", {HistType::kTH1F, {qtAxis}}); rpiKkink.add("h2_kaon_pt_vs_rap_rec_full", "pt_vs_rap_kaon", {HistType::kTH2F, {ptAxis, etaAxis}}); - rpiKkink.add("h2_kaon_pt_vs_rap_rec_full1", "pt_vs_rap_kaon1", {HistType::kTH2F, {ptAxis, etaAxis}}); + // rpiKkink.add("h2_kaon_pt_vs_qt_rec_full1", "pt_vs_qt_kaon1", {HistType::kTH2F, {ptAxis, qtAxis}}); + rpiKkink.add("h2_kaon_pt_vs_qt_rec_full1", "pt_vs_qt_kaon1", {HistType::kTH1F, {qtAxis}}); - rpiKkink.add("h2_moth_ptrapqt_egen", "pt_vs_rap_qt_egen", {HistType::kTH3F, {ptAxis, etaAxis, qtAxis}}); - rpiKkink.add("h2_moth_ptrapqt_mugen", "pt_vs_rap_qt_mugen", {HistType::kTH3F, {ptAxis, etaAxis, qtAxis}}); - rpiKkink.add("h2_moth_ptrapqt_pigen", "pt_vs_rap_qt_pigen", {HistType::kTH3F, {ptAxis, etaAxis, qtAxis}}); + rpiKkink.add("h2_moth_ptrapqt_egen", "pt_vs_rap_qt_egen", {HistType::kTH2F, {ptAxis, qtAxis}}); + rpiKkink.add("h2_moth_ptrapqt_mugen", "pt_vs_rap_qt_mugen", {HistType::kTH2F, {ptAxis, qtAxis}}); + rpiKkink.add("h2_moth_ptrapqt_pigen", "pt_vs_rap_qt_pigen", {HistType::kTH2F, {ptAxis, qtAxis}}); rpiKkink.add("h2_dau_pt_vs_eta_gen", "pt_vs_eta_dau", {HistType::kTH2F, {ptAxis, etaAxis}}); rpiKkink.add("h2_moth_pt_vs_eta_gen", "pt_vs_eta_moth", {HistType::kTH2F, {ptAxis, etaAxis}}); rpiKkink.add("h2_moth_pt_vs_rap_genall", "pt_vs_rap_moth", {HistType::kTH2F, {ptAxis, etaAxis}}); rpiKkink.add("h2_pt_moth_vs_dau_gen", "pt_moth_vs_dau", {HistType::kTH2F, {ptAxis, ptAxis}}); - rpiKkink.add("h1_tracks", "track_cut", {HistType::kTH1F, {{17, 0.5, 17.5}}}); + rpiKkink.add("h1_tracks", "track_cut", {HistType::kTH1F, {{18, 0.5, 18.5}}}); rpiKkink.add("h1_tracks_gen", "track_cut_gen", {HistType::kTH1F, {{15, 0.5, 15.5}}}); rpiKkink.add("h2_qt_gen", "qt", {HistType::kTH1F, {qtAxis}}); rpiKkink.add("h2_qt_rec", "qt", {HistType::kTH1F, {qtAxis}}); @@ -576,10 +624,8 @@ struct SpectraKinkPiKa { continue; rpiKkink.fill(HIST("h1_tracks_data"), 6.0); if (qa) { - rpiKkink.fill(HIST("tpc_dedx"), v0.P(), mothTrack.tpcSignal()); rpiKkink.fill(HIST("tpc_nsigma_kaon"), v0.Pt(), mothTrack.tpcNSigmaKa()); - rpiKkink.fill(HIST("tr_chi2nclM"), mothTrack.tpcChi2NCl()); rpiKkink.fill(HIST("tr_chi2nclD"), dauTrack.tpcChi2NCl()); rpiKkink.fill(HIST("tr_tpcnclfindM"), mothTrack.tpcNClsFound()); @@ -685,7 +731,6 @@ struct SpectraKinkPiKa { if (!collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder) || !collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) { continue; } - float multiplicity{-1}; const int kCentFT0C = 0; const int kCentFT0A = 1; @@ -743,7 +788,6 @@ struct SpectraKinkPiKa { if (qa) { rpiKkink.fill(HIST("tpc_dedx"), v0.P(), mothTrack.tpcSignal()); rpiKkink.fill(HIST("tpc_nsigma_kaon"), v0.Pt(), mothTrack.tpcNSigmaKa()); - rpiKkink.fill(HIST("tr_chi2nclM"), mothTrack.tpcChi2NCl()); rpiKkink.fill(HIST("tr_chi2nclD"), dauTrack.tpcChi2NCl()); rpiKkink.fill(HIST("tr_tpcnclfindM"), mothTrack.tpcNClsFound()); @@ -803,9 +847,6 @@ struct SpectraKinkPiKa { // Compute transverse component TVector3 motherDir(v0.Px(), v0.Py(), v0.Pz()); double ptd = pdlab.Perp(motherDir); // or p_d_lab.Mag() * sin(theta) - - rpiKkink.fill(HIST("h2_qt"), ptd); - double mass = computeMotherMass(v0, v1); rpiKkink.fill(HIST("h2_kaon_mc_rec"), mass, v0.Pt(), ptd, multiplicity); @@ -818,31 +859,30 @@ struct SpectraKinkPiKa { if (!mcTrackMoth.isPhysicalPrimary()) continue; rpiKkink.fill(HIST("h1_tracks"), 13.0); - if (std::abs(mcTrackMoth.pdgCode()) != pid) continue; - rpiKkink.fill(HIST("h2_kaon_pt_vs_rap_rec_full1"), v0.Pt(), v0.Rapidity()); + rpiKkink.fill(HIST("h1_tracks"), 14.0); + // rpiKkink.fill(HIST("h2_kaon_pt_vs_qt_rec_full1"), v0.Pt(), ptd); + rpiKkink.fill(HIST("h2_kaon_pt_vs_qt_rec_full1"), ptd); + if (mcLabDau.has_mcParticle()) { auto mcTrackDau = mcLabDau.mcParticle_as(); if (!mcTrackDau.has_mothers()) continue; - rpiKkink.fill(HIST("h1_tracks"), 14.0); + rpiKkink.fill(HIST("h1_tracks"), 15.0); const int process = 4; if (mcTrackDau.getProcess() != process) continue; - - rpiKkink.fill(HIST("h1_tracks"), 15.0); - + rpiKkink.fill(HIST("h1_tracks"), 16.0); for (const auto& piMother : mcTrackDau.mothers_as()) { if (piMother.globalIndex() != mcTrackMoth.globalIndex()) { continue; } - // std::cout< rapCut) { continue; } + if (std::abs(v0.Pt()) < minPtMothmc) { + continue; + } rpiKkink.fill(HIST("h2_moth_pt_vs_rap_genall"), v0.Pt(), v0.Rapidity()); rpiKkink.fill(HIST("h1_tracks_gen"), 3.0); if (!mcPart.has_daughters()) { @@ -904,24 +947,26 @@ struct SpectraKinkPiKa { double ptd = 0; const int process = 4; for (const auto& daughter : mcPart.daughters_as()) { + v1.SetCoordinates(daughter.px(), daughter.py(), daughter.pz(), o2::constants::physics::MassMuon); if (daughter.getProcess() != process) continue; - v1.SetCoordinates(daughter.px(), daughter.py(), daughter.pz(), o2::constants::physics::MassMuon); ptd = computeQt(v0, v1); if (std::abs(daughter.pdgCode()) == kElectron) { eld++; - } else if (std::abs(daughter.pdgCode()) == kMuonPlus) { + } else if (std::abs(daughter.pdgCode()) == dpid) { muond++; } else if (std::abs(daughter.pdgCode()) == kPiPlus) { piond++; } } - if (eld >= 1) - rpiKkink.fill(HIST("h2_moth_ptrapqt_egen"), v0.Pt(), v0.Rapidity(), ptd); - if (muond >= 1) - rpiKkink.fill(HIST("h2_moth_ptrapqt_mugen"), v0.Pt(), v0.Rapidity(), ptd); - if (piond >= 1) - rpiKkink.fill(HIST("h2_moth_ptrapqt_pigen"), v0.Pt(), v0.Rapidity(), ptd); + if (additionalhist) { + if (eld >= 1) + rpiKkink.fill(HIST("h2_moth_ptrapqt_egen"), v0.Pt(), ptd); + if (muond >= 1) + rpiKkink.fill(HIST("h2_moth_ptrapqt_mugen"), v0.Pt(), ptd); + if (piond >= 1) + rpiKkink.fill(HIST("h2_moth_ptrapqt_pigen"), v0.Pt(), ptd); + } rpiKkink.fill(HIST("h1_tracks_gen"), 5.0); float pMoth = v0.P(); float pDaug = v1.P();