Skip to content

Commit 9d2fc17

Browse files
skundu692skundu692
authored
[PWGLF] Improved DCA fitter for secondary reconstruction and add new PID selection (#9583)
1 parent 815f6df commit 9d2fc17

File tree

1 file changed

+109
-47
lines changed

1 file changed

+109
-47
lines changed

‎PWGLF/Tasks/Resonances/highmasslambda.cxx‎

Lines changed: 109 additions & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -170,7 +170,7 @@ struct highmasslambda {
170170
// std::vector<double> ptLambdaBinning = {2.0, 3.0, 4.0, 5.0, 6.0};
171171

172172
std::vector<double> occupancyBinning = {-0.5, 500.0, 1000.0, 1500.0, 2000.0, 3000.0, 4000.0, 5000.0, 50000.0};
173-
AxisSpec resAxis = {1600, -16, 16, "Res"};
173+
AxisSpec resAxis = {1600, -30, 30, "Res"};
174174
AxisSpec phiAxis = {500, -6.28, 6.28, "phi"};
175175
AxisSpec centAxis = {8, 0, 80, "V0M (%)"};
176176
const AxisSpec thnAxisInvMass{configThnAxisInvMass, "#it{M} (GeV/#it{c}^{2})"};
@@ -212,9 +212,9 @@ struct highmasslambda {
212212
histos.add("hImpactPar1", "hImpactPar1", kTH1F, {{500, 0.0f, 0.1f}});
213213
histos.add("hCPA", "hCPA", kTH1F, {{220, -1.1f, 1.1f}});
214214

215-
histos.add("hSparseV2SASameEvent_V2", "hSparseV2SASameEvent_V2", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisV2, thnAxisDCA});
216-
histos.add("hSparseV2SASameEventRotational_V2", "hSparseV2SASameEventRotational", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisV2, thnAxisDCA});
217-
histos.add("hSparseV2SAMixedEvent_V2", "hSparseV2SAMixedEvent_V2", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisV2, thnAxisDCA});
215+
histos.add("hSparseV2SASameEvent_V2", "hSparseV2SASameEvent_V2", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisV2, thnAxisDCA, thnAxisPtProton});
216+
histos.add("hSparseV2SASameEventRotational_V2", "hSparseV2SASameEventRotational", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisV2, thnAxisDCA, thnAxisPtProton});
217+
histos.add("hSparseV2SAMixedEvent_V2", "hSparseV2SAMixedEvent_V2", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisV2, thnAxisDCA, thnAxisPtProton});
218218

219219
histos.add("hSparseV2SASameEvent_V2_SVX", "hSparseV2SASameEvent_V2_SVX", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisV2, thnAxisDecayLength, thnAxisCPA});
220220
histos.add("hSparseV2SASameEventRotational_V2_SVX", "hSparseV2SASameEventRotational_SVX", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisV2, thnAxisDecayLength, thnAxisCPA});
@@ -307,6 +307,27 @@ struct highmasslambda {
307307
return false;
308308
}
309309

310+
// TPC TOF
311+
template <typename T>
312+
bool selectionPID3(const T& candidate)
313+
{
314+
if (candidate.hasTOF() && candidate.tpcInnerParam() > 0.6 && TMath::Abs(candidate.tpcNSigmaPr()) < nsigmaCutTPC && TMath::Abs(candidate.tofNSigmaPr()) < nsigmaCutTOF) {
315+
return true;
316+
}
317+
if (!candidate.hasTOF()) {
318+
if (candidate.tpcInnerParam() < 1.0 && TMath::Abs(candidate.tpcNSigmaPr()) < nsigmaCutTPC) {
319+
return true;
320+
}
321+
if (candidate.tpcInnerParam() >= 1.0 && candidate.tpcInnerParam() < 1.6 && candidate.tpcNSigmaPr() < nsigmaCutTPC && candidate.tpcNSigmaPr() > -1.5) {
322+
return true;
323+
}
324+
if (candidate.tpcInnerParam() >= 1.6 && TMath::Abs(candidate.tpcNSigmaPr()) < nsigmaCutTPC) {
325+
return true;
326+
}
327+
}
328+
return false;
329+
}
330+
310331
template <typename Collision, typename V0>
311332
bool SelectionV0(Collision const& collision, V0 const& candidate)
312333
{
@@ -448,6 +469,9 @@ struct highmasslambda {
448469
if (PIDstrategy == 1 && !selectionPID2(track1)) {
449470
continue;
450471
}
472+
if (PIDstrategy == 2 && !selectionPID3(track1)) {
473+
continue;
474+
}
451475
if (track1.p() < 1.0 && !(itsResponse.nSigmaITS<o2::track::PID::Proton>(track1) > -2.5 && itsResponse.nSigmaITS<o2::track::PID::Proton>(track1) < 2.5)) {
452476
continue;
453477
}
@@ -490,7 +514,7 @@ struct highmasslambda {
490514
auto phiminuspsi = GetPhiInRange(Lambdac.Phi() - psiFT0C);
491515
v2 = TMath::Cos(2.0 * phiminuspsi) * QFT0C;
492516
if (Lambdac.M() > cMinLambdaMass && Lambdac.M() <= cMaxLambdaMass && TMath::Abs(Lambdac.Rapidity()) < confRapidity && Lambdac.Pt() > 2.0 && Lambdac.Pt() <= 6.0) {
493-
histos.fill(HIST("hSparseV2SASameEvent_V2"), Lambdac.M(), Lambdac.Pt(), v2, TMath::Abs(track1.dcaXY()));
517+
histos.fill(HIST("hSparseV2SASameEvent_V2"), Lambdac.M(), Lambdac.Pt(), v2, TMath::Abs(track1.dcaXY()), Proton.Pt());
494518
}
495519
if (fillRotation) {
496520
for (int nrotbkg = 0; nrotbkg < nBkgRotations; nrotbkg++) {
@@ -506,7 +530,7 @@ struct highmasslambda {
506530
auto phiminuspsiRot = GetPhiInRange(LambdacRot.Phi() - psiFT0C);
507531
v2Rot = TMath::Cos(2.0 * phiminuspsiRot) * QFT0C;
508532
if (LambdacRot.M() > cMinLambdaMass && LambdacRot.M() <= cMaxLambdaMass && TMath::Abs(LambdacRot.Rapidity()) < confRapidity && LambdacRot.Pt() > 2.0 && LambdacRot.Pt() <= 6.0) {
509-
histos.fill(HIST("hSparseV2SASameEventRotational_V2"), LambdacRot.M(), LambdacRot.Pt(), v2Rot, TMath::Abs(track1.dcaXY()));
533+
histos.fill(HIST("hSparseV2SASameEventRotational_V2"), LambdacRot.M(), LambdacRot.Pt(), v2Rot, TMath::Abs(track1.dcaXY()), Proton.Pt());
510534
}
511535
}
512536
}
@@ -559,6 +583,9 @@ struct highmasslambda {
559583
if (PIDstrategy == 1 && !selectionPID2(track1)) {
560584
continue;
561585
}
586+
if (PIDstrategy == 2 && !selectionPID3(track1)) {
587+
continue;
588+
}
562589
if (track1.p() < 1.0 && !(itsResponse.nSigmaITS<o2::track::PID::Proton>(track1) > -2.5 && itsResponse.nSigmaITS<o2::track::PID::Proton>(track1) < 2.5)) {
563590
continue;
564591
}
@@ -586,7 +613,7 @@ struct highmasslambda {
586613
auto phiminuspsi = GetPhiInRange(Lambdac.Phi() - psiFT0C);
587614
v2 = TMath::Cos(2.0 * phiminuspsi) * QFT0C;
588615
if (occupancy1 < cfgOccupancyCut && occupancy2 < cfgOccupancyCut && Lambdac.M() > cMinLambdaMass && Lambdac.M() <= cMaxLambdaMass && TMath::Abs(Lambdac.Rapidity()) < confRapidity && Lambdac.Pt() > 2.0 && Lambdac.Pt() <= 6.0) {
589-
histos.fill(HIST("hSparseV2SAMixedEvent_V2"), Lambdac.M(), Lambdac.Pt(), v2, TMath::Abs(track1.dcaXY()));
616+
histos.fill(HIST("hSparseV2SAMixedEvent_V2"), Lambdac.M(), Lambdac.Pt(), v2, TMath::Abs(track1.dcaXY()), Proton.Pt());
590617
}
591618
}
592619
}
@@ -655,6 +682,9 @@ struct highmasslambda {
655682
if (PIDstrategy == 1 && !selectionPID2(track1)) {
656683
continue;
657684
}
685+
if (PIDstrategy == 2 && !selectionPID3(track1)) {
686+
continue;
687+
}
658688
if (track1.p() < 1.0 && !(itsResponse.nSigmaITS<o2::track::PID::Proton>(track1) > -2.5 && itsResponse.nSigmaITS<o2::track::PID::Proton>(track1) < 2.5)) {
659689
continue;
660690
}
@@ -667,6 +697,7 @@ struct highmasslambda {
667697
}
668698
auto track1ID = track1.globalIndex();
669699
auto trackParCovBach = getTrackParCov(track1);
700+
// auto trackParCovBach = getTrackParCov(bach);
670701

671702
for (auto v0 : V0s) {
672703
if (!SelectionV0(collision, v0)) {
@@ -693,53 +724,84 @@ struct highmasslambda {
693724
histos.fill(HIST("hInvMassKs0"), v0.mK0Short());
694725
}
695726
firstprimarytrack = firstprimarytrack + 1;
696-
float v0x, v0y, v0z, v0px, v0py, v0pz;
697-
float posTrackX, negTrackX;
698-
o2::track::TrackParCov trackParCovV0DaughPos;
699-
o2::track::TrackParCov trackParCovV0DaughNeg;
700-
trackParCovV0DaughPos = getTrackParCov(postrack); // check that aod::TracksWCov does not need TracksDCA!
701-
trackParCovV0DaughNeg = getTrackParCov(negtrack); // check that aod::TracksWCov does not need TracksDCA!
702-
posTrackX = v0.posX();
703-
negTrackX = v0.negX();
704-
v0x = v0.x();
705-
v0y = v0.y();
706-
v0z = v0.z();
707-
const std::array<float, 3> vertexV0 = {v0x, v0y, v0z};
708-
709-
v0px = v0.px();
710-
v0py = v0.py();
711-
v0pz = v0.pz();
712-
const std::array<float, 3> momentumV0 = {v0px, v0py, v0pz};
713-
714-
std::array<float, 6> covV0Pos = {0.};
727+
// LOGF(info, "Before dca fitter");
728+
std::array<float, 3> pVecV0 = {0., 0., 0.};
729+
std::array<float, 3> pVecBach = {0., 0., 0.};
730+
// std::array<float, 3> pVecCand = {0., 0., 0.};
731+
const std::array<float, 3> vertexV0 = {v0.x(), v0.y(), v0.z()};
732+
const std::array<float, 3> momentumV0 = {v0.px(), v0.py(), v0.pz()};
733+
// we build the neutral track to then build the cascade
734+
std::array<float, 21> covV = {0.};
735+
constexpr int MomInd[6] = {9, 13, 14, 18, 19, 20}; // cov matrix elements for momentum component
715736
for (int i = 0; i < 6; i++) {
716-
covV0Pos[i] = v0.positionCovMat()[i];
737+
covV[MomInd[i]] = v0.momentumCovMat()[i];
738+
covV[i] = v0.positionCovMat()[i];
717739
}
718-
trackParCovV0DaughPos.propagateTo(posTrackX, bz); // propagate the track to the X closest to the V0 vertex
719-
trackParCovV0DaughNeg.propagateTo(negTrackX, bz); // propagate the track to the X closest to the V0 vertex
740+
auto trackV0 = o2::track::TrackParCov(vertexV0, momentumV0, covV, 0, true);
741+
trackV0.setAbsCharge(0);
742+
trackV0.setPID(o2::track::PID::K0);
720743

721-
// we build the neutral track to then build the cascade
722-
// auto trackV0 = o2::dataformats::V0(vertexV0, momentumV0, {0, 0, 0, 0, 0, 0}, trackParCovV0DaughPos, trackParCovV0DaughNeg); // build the V0 track (indices for v0 daughters set to 0 for now)
723-
auto trackV0 = o2::dataformats::V0(vertexV0, momentumV0, covV0Pos, trackParCovV0DaughPos, trackParCovV0DaughNeg); // build the V0 track (indices for v0 daughters set to 0 for now)
724-
std::array<float, 3> pVecV0 = {0., 0., 0.};
725-
std::array<float, 3> pVecBach = {0., 0., 0.};
726-
std::array<float, 3> pVecCand = {0., 0., 0.};
744+
int nCand2 = 0;
727745
try {
728-
if (df.process(trackV0, trackParCovBach) == 0) {
729-
continue;
730-
} else {
731-
}
732-
} catch (const std::runtime_error& error) {
746+
nCand2 = df.process(trackV0, trackParCovBach);
747+
} catch (...) {
733748
continue;
734749
}
735-
df.propagateTracksToVertex(); // propagate the bachelor and V0 to the Lambdac vertex
736-
trackV0.getPxPyPzGlo(pVecV0); // momentum of D0 at the Lambdac vertex
737-
trackParCovBach.getPxPyPzGlo(pVecBach); // momentum of proton at the Lambdac vertex
738-
pVecCand = RecoDecay::pVec(pVecV0, pVecBach);
739-
const auto& secondaryVertex = df.getPCACandidate();
740750

741-
// get track impact parameters
742-
// This modifies track momenta!
751+
if (nCand2 == 0) {
752+
continue;
753+
}
754+
df.propagateTracksToVertex(); // propagate the bach and V0 to the Lc vertex
755+
df.getTrack(0).getPxPyPzGlo(pVecV0); // take the momentum at the Lc vertex
756+
df.getTrack(1).getPxPyPzGlo(pVecBach);
757+
// LOGF(info, "after dca fitter");
758+
759+
/*
760+
float v0x, v0y, v0z, v0px, v0py, v0pz;
761+
float posTrackX, negTrackX;
762+
o2::track::TrackParCov trackParCovV0DaughPos;
763+
o2::track::TrackParCov trackParCovV0DaughNeg;
764+
trackParCovV0DaughPos = getTrackParCov(postrack); // check that aod::TracksWCov does not need TracksDCA!
765+
trackParCovV0DaughNeg = getTrackParCov(negtrack); // check that aod::TracksWCov does not need TracksDCA!
766+
posTrackX = v0.posX();
767+
negTrackX = v0.negX();
768+
v0x = v0.x();
769+
v0y = v0.y();
770+
v0z = v0.z();
771+
const std::array<float, 3> vertexV0 = {v0x, v0y, v0z};
772+
773+
v0px = v0.px();
774+
v0py = v0.py();
775+
v0pz = v0.pz();
776+
const std::array<float, 3> momentumV0 = {v0px, v0py, v0pz};
777+
778+
std::array<float, 6> covV0Pos = {0.};
779+
for (int i = 0; i < 6; i++) {
780+
covV0Pos[i] = v0.positionCovMat()[i];
781+
}
782+
trackParCovV0DaughPos.propagateTo(posTrackX, bz); // propagate the track to the X closest to the V0 vertex
783+
trackParCovV0DaughNeg.propagateTo(negTrackX, bz); // propagate the track to the X closest to the V0 vertex
784+
785+
// we build the neutral track to then build the cascade
786+
// auto trackV0 = o2::dataformats::V0(vertexV0, momentumV0, {0, 0, 0, 0, 0, 0}, trackParCovV0DaughPos, trackParCovV0DaughNeg); // build the V0 track (indices for v0 daughters set to 0 for now)
787+
auto trackV0 = o2::dataformats::V0(vertexV0, momentumV0, covV0Pos, trackParCovV0DaughPos, trackParCovV0DaughNeg); // build the V0 track (indices for v0 daughters set to 0 for now)
788+
std::array<float, 3> pVecV0 = {0., 0., 0.};
789+
std::array<float, 3> pVecBach = {0., 0., 0.};
790+
std::array<float, 3> pVecCand = {0., 0., 0.};
791+
try {
792+
if (df.process(trackV0, trackParCovBach) == 0) {
793+
continue;
794+
} else {
795+
}
796+
} catch (const std::runtime_error& error) {
797+
continue;
798+
}
799+
df.propagateTracksToVertex(); // propagate the bachelor and V0 to the Lambdac vertex
800+
trackV0.getPxPyPzGlo(pVecV0); // momentum of D0 at the Lambdac vertex
801+
trackParCovBach.getPxPyPzGlo(pVecBach); // momentum of proton at the Lambdac vertex
802+
*/
803+
804+
const auto& secondaryVertex = df.getPCACandidate();
743805
auto primaryVertex = getPrimaryVertex(collision);
744806
o2::dataformats::DCA impactParameter0;
745807
o2::dataformats::DCA impactParameter1;

0 commit comments

Comments
 (0)