Skip to content

Commit 24d7b68

Browse files
sdudi123sandeep dudi
sdudi123
and
sandeep dudi
authored
[PWGLF] Improvement in Kink vertex finder (#13921)
Co-authored-by: sandeep dudi <sandeep.dudi@cern.ch>
1 parent 1a9840d commit 24d7b68

File tree

1 file changed

+94
-49
lines changed

1 file changed

+94
-49
lines changed

‎PWGLF/Tasks/Nuspex/spectraKinkPiKa.cxx‎

Lines changed: 94 additions & 49 deletions
Original file line numberDiff line numberDiff line change
@@ -67,7 +67,7 @@ using CollisionsFull = soa::Join<aod::Collisions, aod::EvSel, aod::FT0Mults, aod
6767
using CollisionsFullMC = soa::Join<aod::Collisions, aod::McCollisionLabels, aod::EvSel, aod::CentFT0Ms, aod::CentFT0As, aod::CentFT0Cs, aod::CentFV0As>;
6868
namespace
6969
{
70-
constexpr std::array<float, 7> LayerRadii{2.33959f, 3.14076f, 3.91924f, 19.6213f, 24.5597f, 34.388f, 39.3329f};
70+
// constexpr std::array<float, 7> LayerRadii{2.33959f, 3.14076f, 3.91924f, 19.6213f, 24.5597f, 34.388f, 39.3329f};
7171
constexpr double BetheBlochDefault[1][6]{{-1.e32, -1.e32, -1.e32, -1.e32, -1.e32, -1.e32}};
7272
static const std::vector<std::string> betheBlochParNames{"p0", "p1", "p2", "p3", "p4", "resolution"};
7373
static const std::vector<std::string> particleNames{"Daughter"};
@@ -96,7 +96,7 @@ struct KinkBuilder {
9696
Service<o2::ccdb::BasicCCDBManager> ccdb;
9797
// Selection criteria
9898
Configurable<float> maxDCAMothToPV{"maxDCAMothToPV", 0.2, "Max DCA of the mother to the PV"};
99-
Configurable<float> minDCADaugToPV{"minDCADaugToPV", 0., "Min DCA of the daughter to the PV"};
99+
Configurable<float> minDCADaugToPV{"minDCADaugToPV", 0.1, "Min DCA of the daughter to the PV"};
100100
Configurable<float> minPtMoth{"minPtMoth", 0.15, "Minimum pT of the hypercandidate"};
101101
Configurable<float> maxZDiff{"maxZDiff", 20., "Max z difference between the kink daughter and the mother"};
102102
Configurable<float> maxPhiDiff{"maxPhiDiff", 100, "Max phi difference between the kink daughter and the mother"};
@@ -107,6 +107,10 @@ struct KinkBuilder {
107107
Configurable<float> itsChi2cut{"itsChi2cut", 36, "mother itsChi2 cut"};
108108
Configurable<bool> askTOFforDaug{"askTOFforDaug", false, "If true, ask for TOF signal"};
109109
Configurable<bool> kaontopologhy{"kaontopologhy", true, "If true, selected mother have both ITS+TPC "};
110+
Configurable<bool> vertexfinding{"vertexfinding", false, "If true, find the vextex in TPC and applied cut of z and phi"};
111+
Configurable<float> rMin{"rMin", 120., "min radius for kink vertex"};
112+
Configurable<float> rMax{"rMax", 200., "max radius for kink vertex"};
113+
Configurable<float> rStep{"rStep", 2., "step size for scan radius in tpc"};
110114

111115
o2::vertexing::DCAFitterN<2> fitter;
112116
o2::base::MatLayerCylSet* lut = nullptr;
@@ -214,10 +218,15 @@ struct KinkBuilder {
214218
svCreator.fillBC2Coll(collisions, bcs);
215219
for (const auto& track : tracks) {
216220
bool isDaug = selectDaugTrack(track);
221+
217222
bool isMoth = selectMothTrack(track);
218223

219224
if (!isDaug && !isMoth)
220225
continue;
226+
if (isDaug && track.isPVContributor())
227+
continue;
228+
if (isMoth && !track.isPVContributor())
229+
continue;
221230
if (isDaug && std::abs(track.eta()) > etaMaxDaug)
222231
continue;
223232
if (isMoth && std::abs(track.eta()) > etaMaxMoth)
@@ -245,25 +254,65 @@ struct KinkBuilder {
245254

246255
o2::track::TrackParCov trackParCovMoth = getTrackParCov(trackMoth);
247256
o2::track::TrackParCov trackParCovMothPV{trackParCovMoth};
248-
o2::base::Propagator::Instance()->PropagateToXBxByBz(trackParCovMoth, LayerRadii[trackMoth.itsNCls() - 1]);
249257
std::array<float, 2> dcaInfoMoth;
250258
o2::base::Propagator::Instance()->propagateToDCABxByBz({primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, trackParCovMothPV, 2.f, static_cast<o2::base::Propagator::MatCorrType>(cfgMaterialCorrection.value), &dcaInfoMoth);
251-
252259
o2::track::TrackParCov trackParCovDaug = getTrackParCov(trackDaug);
260+
// propagate to PV
261+
std::array<float, 2> dcaInfoDaug;
262+
o2::base::Propagator::Instance()->propagateToDCABxByBz({primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, trackParCovDaug, 2.f, static_cast<o2::base::Propagator::MatCorrType>(cfgMaterialCorrection.value), &dcaInfoDaug);
253263

254-
// check if the kink daughter is close to the mother
255-
if (std::abs(trackParCovMoth.getZ() - trackParCovDaug.getZ()) > maxZDiff) {
264+
if (std::abs(dcaInfoMoth[1]) > maxDCAMothToPV) {
256265
continue;
257266
}
258-
if ((std::abs(trackParCovMoth.getPhi() - trackParCovDaug.getPhi()) * radToDeg) > maxPhiDiff) {
267+
if (std::abs(dcaInfoDaug[1]) < minDCADaugToPV) {
259268
continue;
260269
}
261-
// propagate to PV
262-
std::array<float, 2> dcaInfoDaug;
263-
o2::base::Propagator::Instance()->propagateToDCABxByBz({primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, trackParCovDaug, 2.f, static_cast<o2::base::Propagator::MatCorrType>(cfgMaterialCorrection.value), &dcaInfoDaug);
264-
if (std::abs(dcaInfoDaug[0]) < minDCADaugToPV) {
270+
271+
if (vertexfinding) {
272+
float bestR = -1;
273+
float bestDeltaPhi = 999;
274+
float bestDeltaZ = 999;
275+
float bestCost = 1e12;
276+
// make local copies (don’t modify originals)
277+
auto mothTmp0 = trackParCovMoth;
278+
auto daugTmp0 = trackParCovDaug;
279+
const float minr = rMin;
280+
const float maxr = rMax;
281+
const float rs = rStep;
282+
for (float R = minr; R <= maxr; R += rs) {
283+
auto mothTmp = mothTmp0;
284+
auto daugTmp = daugTmp0;
285+
if (!o2::base::Propagator::Instance()->PropagateToXBxByBz(mothTmp, R))
286+
continue;
287+
if (!o2::base::Propagator::Instance()->PropagateToXBxByBz(daugTmp, R))
288+
continue;
289+
float dphi = std::abs(mothTmp.getPhi() - daugTmp.getPhi());
290+
if (dphi > M_PI)
291+
dphi = 2 * M_PI - dphi; // wrap φ
292+
float dz = std::abs(mothTmp.getZ() - daugTmp.getZ());
293+
float cost = dphi * dphi + dz * dz; // <-- correct metric
294+
if (cost < bestCost) {
295+
bestCost = cost;
296+
bestDeltaPhi = dphi;
297+
bestDeltaZ = dz;
298+
bestR = R;
299+
}
300+
}
301+
if (bestR < 0)
302+
continue;
303+
if (bestDeltaZ > maxZDiff)
304+
continue;
305+
if (bestDeltaPhi * radToDeg > maxPhiDiff)
306+
continue;
307+
}
308+
/* // check if the kink daughter is close to the mother
309+
if (std::abs(trackParCovMoth.getZ() - trackParCovDaug.getZ()) > maxZDiff) {
310+
continue;
311+
}
312+
if ((std::abs(trackParCovMoth.getPhi() - trackParCovDaug.getPhi()) * radToDeg) > maxPhiDiff) {
265313
continue;
266314
}
315+
*/
267316
int nCand = 0;
268317
try {
269318
nCand = fitter.process(trackParCovMoth, trackParCovDaug);
@@ -367,7 +416,7 @@ struct SpectraKinkPiKa {
367416
Configurable<float> tpcChi2Cut{"tpcChi2Cut", 4.0, "tpcChi2Cut"};
368417
Configurable<float> minqt{"minqt", 0.12, "min qt for kaon"};
369418
Configurable<float> maxqt{"maxqt", 0.3, "max qt for kaon"};
370-
419+
Configurable<float> minPtMothmc{"minPtMothmc", 0.15, "Minimum pT of the mother"};
371420
Configurable<int> centestimator{"centestimator", 0, "Select multiplicity estimator: 0 - FT0C, 1 - FT0A, 2 - FT0M, 3 - FV0A, 4 - PVTracks"};
372421
Configurable<int> pid{"pid", 321, ""};
373422
Configurable<int> dpid{"dpid", 13, ""};
@@ -400,14 +449,13 @@ struct SpectraKinkPiKa {
400449
// Event selection
401450
rEventSelection.add("hVertexZRec", "hVertexZRec", {HistType::kTH1F, {{200, -20.0, 20.0}}});
402451
rEventSelection.add("hMultiplicity", "hMultiplicity", {HistType::kTH1F, {multAxis}});
403-
404-
rpiKkink.add("h2_dau_pt_vs_eta_rec", "pt_vs_eta_dau", {HistType::kTH3F, {ptAxis, etaAxis, multAxis}});
405-
rpiKkink.add("h2_moth_pt_vs_eta_rec", "pt_vs_eta_moth", {HistType::kTH3F, {ptAxis, etaAxis, multAxis}});
406-
rpiKkink.add("h2_pt_moth_vs_dau_rec", "pt_moth_vs_dau", {HistType::kTH2F, {ptAxis, ptAxis}});
407-
408-
rpiKkink.add("h2_qt_vs_pt", "qt_pt", {HistType::kTH2F, {qtAxis, ptAxis}});
452+
if (additionalhist) {
453+
rpiKkink.add("h2_dau_pt_vs_eta_rec", "pt_vs_eta_dau", {HistType::kTH3F, {ptAxis, etaAxis, multAxis}});
454+
rpiKkink.add("h2_moth_pt_vs_eta_rec", "pt_vs_eta_moth", {HistType::kTH3F, {ptAxis, etaAxis, multAxis}});
455+
rpiKkink.add("h2_pt_moth_vs_dau_rec", "pt_moth_vs_dau", {HistType::kTH2F, {ptAxis, ptAxis}});
456+
rpiKkink.add("h2_qt_vs_pt", "qt_pt", {HistType::kTH2F, {qtAxis, ptAxis}});
457+
}
409458
rpiKkink.add("h2_kink_angle", "kink angle", {HistType::kTH2F, {ptAxis, kinkAxis}});
410-
411459
// inv mass
412460
rpiKkink.add("h2_kaon_data", "h2_kaon_data", HistType::kTHnSparseF, {massAxis, ptAxis, qtAxis, multAxis}, true);
413461
rpiKkink.add("h1_tracks_data", "track_cut_data", {HistType::kTH1F, {{17, 0.5, 17.5}}});
@@ -454,19 +502,19 @@ struct SpectraKinkPiKa {
454502
}
455503
if (doprocessMC) {
456504

457-
rpiKkink.add("h2_qt", "qt", {HistType::kTH1F, {qtAxis}});
458505
rpiKkink.add("h2_kaon_pt_vs_rap_rec_full", "pt_vs_rap_kaon", {HistType::kTH2F, {ptAxis, etaAxis}});
459-
rpiKkink.add("h2_kaon_pt_vs_rap_rec_full1", "pt_vs_rap_kaon1", {HistType::kTH2F, {ptAxis, etaAxis}});
506+
// rpiKkink.add("h2_kaon_pt_vs_qt_rec_full1", "pt_vs_qt_kaon1", {HistType::kTH2F, {ptAxis, qtAxis}});
507+
rpiKkink.add("h2_kaon_pt_vs_qt_rec_full1", "pt_vs_qt_kaon1", {HistType::kTH1F, {qtAxis}});
460508

461-
rpiKkink.add("h2_moth_ptrapqt_egen", "pt_vs_rap_qt_egen", {HistType::kTH3F, {ptAxis, etaAxis, qtAxis}});
462-
rpiKkink.add("h2_moth_ptrapqt_mugen", "pt_vs_rap_qt_mugen", {HistType::kTH3F, {ptAxis, etaAxis, qtAxis}});
463-
rpiKkink.add("h2_moth_ptrapqt_pigen", "pt_vs_rap_qt_pigen", {HistType::kTH3F, {ptAxis, etaAxis, qtAxis}});
509+
rpiKkink.add("h2_moth_ptrapqt_egen", "pt_vs_rap_qt_egen", {HistType::kTH2F, {ptAxis, qtAxis}});
510+
rpiKkink.add("h2_moth_ptrapqt_mugen", "pt_vs_rap_qt_mugen", {HistType::kTH2F, {ptAxis, qtAxis}});
511+
rpiKkink.add("h2_moth_ptrapqt_pigen", "pt_vs_rap_qt_pigen", {HistType::kTH2F, {ptAxis, qtAxis}});
464512

465513
rpiKkink.add("h2_dau_pt_vs_eta_gen", "pt_vs_eta_dau", {HistType::kTH2F, {ptAxis, etaAxis}});
466514
rpiKkink.add("h2_moth_pt_vs_eta_gen", "pt_vs_eta_moth", {HistType::kTH2F, {ptAxis, etaAxis}});
467515
rpiKkink.add("h2_moth_pt_vs_rap_genall", "pt_vs_rap_moth", {HistType::kTH2F, {ptAxis, etaAxis}});
468516
rpiKkink.add("h2_pt_moth_vs_dau_gen", "pt_moth_vs_dau", {HistType::kTH2F, {ptAxis, ptAxis}});
469-
rpiKkink.add("h1_tracks", "track_cut", {HistType::kTH1F, {{17, 0.5, 17.5}}});
517+
rpiKkink.add("h1_tracks", "track_cut", {HistType::kTH1F, {{18, 0.5, 18.5}}});
470518
rpiKkink.add("h1_tracks_gen", "track_cut_gen", {HistType::kTH1F, {{15, 0.5, 15.5}}});
471519
rpiKkink.add("h2_qt_gen", "qt", {HistType::kTH1F, {qtAxis}});
472520
rpiKkink.add("h2_qt_rec", "qt", {HistType::kTH1F, {qtAxis}});
@@ -576,10 +624,8 @@ struct SpectraKinkPiKa {
576624
continue;
577625
rpiKkink.fill(HIST("h1_tracks_data"), 6.0);
578626
if (qa) {
579-
580627
rpiKkink.fill(HIST("tpc_dedx"), v0.P(), mothTrack.tpcSignal());
581628
rpiKkink.fill(HIST("tpc_nsigma_kaon"), v0.Pt(), mothTrack.tpcNSigmaKa());
582-
583629
rpiKkink.fill(HIST("tr_chi2nclM"), mothTrack.tpcChi2NCl());
584630
rpiKkink.fill(HIST("tr_chi2nclD"), dauTrack.tpcChi2NCl());
585631
rpiKkink.fill(HIST("tr_tpcnclfindM"), mothTrack.tpcNClsFound());
@@ -685,7 +731,6 @@ struct SpectraKinkPiKa {
685731
if (!collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder) || !collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) {
686732
continue;
687733
}
688-
689734
float multiplicity{-1};
690735
const int kCentFT0C = 0;
691736
const int kCentFT0A = 1;
@@ -743,7 +788,6 @@ struct SpectraKinkPiKa {
743788
if (qa) {
744789
rpiKkink.fill(HIST("tpc_dedx"), v0.P(), mothTrack.tpcSignal());
745790
rpiKkink.fill(HIST("tpc_nsigma_kaon"), v0.Pt(), mothTrack.tpcNSigmaKa());
746-
747791
rpiKkink.fill(HIST("tr_chi2nclM"), mothTrack.tpcChi2NCl());
748792
rpiKkink.fill(HIST("tr_chi2nclD"), dauTrack.tpcChi2NCl());
749793
rpiKkink.fill(HIST("tr_tpcnclfindM"), mothTrack.tpcNClsFound());
@@ -803,9 +847,6 @@ struct SpectraKinkPiKa {
803847
// Compute transverse component
804848
TVector3 motherDir(v0.Px(), v0.Py(), v0.Pz());
805849
double ptd = pdlab.Perp(motherDir); // or p_d_lab.Mag() * sin(theta)
806-
807-
rpiKkink.fill(HIST("h2_qt"), ptd);
808-
809850
double mass = computeMotherMass(v0, v1);
810851

811852
rpiKkink.fill(HIST("h2_kaon_mc_rec"), mass, v0.Pt(), ptd, multiplicity);
@@ -818,31 +859,30 @@ struct SpectraKinkPiKa {
818859
if (!mcTrackMoth.isPhysicalPrimary())
819860
continue;
820861
rpiKkink.fill(HIST("h1_tracks"), 13.0);
821-
822862
if (std::abs(mcTrackMoth.pdgCode()) != pid)
823863
continue;
824-
rpiKkink.fill(HIST("h2_kaon_pt_vs_rap_rec_full1"), v0.Pt(), v0.Rapidity());
864+
rpiKkink.fill(HIST("h1_tracks"), 14.0);
865+
// rpiKkink.fill(HIST("h2_kaon_pt_vs_qt_rec_full1"), v0.Pt(), ptd);
866+
rpiKkink.fill(HIST("h2_kaon_pt_vs_qt_rec_full1"), ptd);
867+
825868
if (mcLabDau.has_mcParticle()) {
826869
auto mcTrackDau = mcLabDau.mcParticle_as<aod::McParticles>();
827870
if (!mcTrackDau.has_mothers())
828871
continue;
829-
rpiKkink.fill(HIST("h1_tracks"), 14.0);
872+
rpiKkink.fill(HIST("h1_tracks"), 15.0);
830873
const int process = 4;
831874
if (mcTrackDau.getProcess() != process)
832875
continue;
833-
834-
rpiKkink.fill(HIST("h1_tracks"), 15.0);
835-
876+
rpiKkink.fill(HIST("h1_tracks"), 16.0);
836877
for (const auto& piMother : mcTrackDau.mothers_as<aod::McParticles>()) {
837878
if (piMother.globalIndex() != mcTrackMoth.globalIndex()) {
838879
continue;
839880
}
840-
// std::cout<<piMother.pdgCode()<<" check--- "<<mcTrackMoth.pdgCode()<<std::endl;
841-
rpiKkink.fill(HIST("h1_tracks"), 16.0);
881+
rpiKkink.fill(HIST("h1_tracks"), 17.0);
842882
if (dpidCut && std::abs(mcTrackDau.pdgCode()) != dpid) {
843883
continue;
844884
}
845-
rpiKkink.fill(HIST("h1_tracks"), 17.0);
885+
rpiKkink.fill(HIST("h1_tracks"), 18.0);
846886
if (additionalhist) {
847887
rpiKkink.fill(HIST("h2_moth_pt_vs_eta_rec_m"), v0.Pt(), v0.Eta(), multiplicity);
848888
rpiKkink.fill(HIST("h2_dau_pt_vs_eta_rec_m"), v1.Pt(), v1.Eta(), multiplicity);
@@ -892,6 +932,9 @@ struct SpectraKinkPiKa {
892932
if (std::abs(v0.Rapidity()) > rapCut) {
893933
continue;
894934
}
935+
if (std::abs(v0.Pt()) < minPtMothmc) {
936+
continue;
937+
}
895938
rpiKkink.fill(HIST("h2_moth_pt_vs_rap_genall"), v0.Pt(), v0.Rapidity());
896939
rpiKkink.fill(HIST("h1_tracks_gen"), 3.0);
897940
if (!mcPart.has_daughters()) {
@@ -904,24 +947,26 @@ struct SpectraKinkPiKa {
904947
double ptd = 0;
905948
const int process = 4;
906949
for (const auto& daughter : mcPart.daughters_as<aod::McParticles>()) {
950+
v1.SetCoordinates(daughter.px(), daughter.py(), daughter.pz(), o2::constants::physics::MassMuon);
907951
if (daughter.getProcess() != process)
908952
continue;
909-
v1.SetCoordinates(daughter.px(), daughter.py(), daughter.pz(), o2::constants::physics::MassMuon);
910953
ptd = computeQt(v0, v1);
911954
if (std::abs(daughter.pdgCode()) == kElectron) {
912955
eld++;
913-
} else if (std::abs(daughter.pdgCode()) == kMuonPlus) {
956+
} else if (std::abs(daughter.pdgCode()) == dpid) {
914957
muond++;
915958
} else if (std::abs(daughter.pdgCode()) == kPiPlus) {
916959
piond++;
917960
}
918961
}
919-
if (eld >= 1)
920-
rpiKkink.fill(HIST("h2_moth_ptrapqt_egen"), v0.Pt(), v0.Rapidity(), ptd);
921-
if (muond >= 1)
922-
rpiKkink.fill(HIST("h2_moth_ptrapqt_mugen"), v0.Pt(), v0.Rapidity(), ptd);
923-
if (piond >= 1)
924-
rpiKkink.fill(HIST("h2_moth_ptrapqt_pigen"), v0.Pt(), v0.Rapidity(), ptd);
962+
if (additionalhist) {
963+
if (eld >= 1)
964+
rpiKkink.fill(HIST("h2_moth_ptrapqt_egen"), v0.Pt(), ptd);
965+
if (muond >= 1)
966+
rpiKkink.fill(HIST("h2_moth_ptrapqt_mugen"), v0.Pt(), ptd);
967+
if (piond >= 1)
968+
rpiKkink.fill(HIST("h2_moth_ptrapqt_pigen"), v0.Pt(), ptd);
969+
}
925970
rpiKkink.fill(HIST("h1_tracks_gen"), 5.0);
926971
float pMoth = v0.P();
927972
float pDaug = v1.P();

0 commit comments

Comments
 (0)