Skip to content

Commit 49b549d

Browse files
creetz16creetz16
committed
Add option to use KFParticle for two-body vertex fit
Add DCA of ITS matched track to SV to output table
1 parent 28f0ea4 commit 49b549d

3 files changed

Lines changed: 185 additions & 30 deletions

File tree

‎PWGLF/DataModel/Vtx3BodyTables.h‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -140,6 +140,10 @@ DECLARE_SOA_COLUMN(IsSurvEvSel, isSurvEvSel, int); //! flag if reco collisio
140140
DECLARE_SOA_COLUMN(IsReco, isreco, int); //! flag if candidate was reconstructed
141141
DECLARE_SOA_COLUMN(MotherLabel, motherLabel, int); //! label of the mother particle (signal: MC index, bkg: negative value)
142142

143+
// strangeness tracking information
144+
DECLARE_SOA_COLUMN(ITSTrackDCAXYToSVXY, itsTrackDcaXYToSv, float); //! DCAxy of ITS track to SV
145+
DECLARE_SOA_COLUMN(ITSTrackDCAZToSVXY, itsTrackDcaZToSv, float); //! DCAz of ITS track to SV
146+
143147
// Derived expressions
144148
// Momenta
145149
DECLARE_SOA_DYNAMIC_COLUMN(Pt, pt, //! 3 body pT in GeV/c
@@ -258,6 +262,11 @@ DECLARE_SOA_TABLE(Vtx3BodyCovs, "AOD", "VTX3BODYCOV", //!
258262
vtx3body::CovProton, vtx3body::CovPion, vtx3body::CovDeuteron,
259263
vtx3body::VtxCovMat);
260264

265+
// strangeness tracked extra table
266+
DECLARE_SOA_TABLE(Vtx3BodyTrackedInfo, "AOD", "VTX3BODYTR", //!
267+
vtx3body::ITSTrackDCAXYToSVXY,
268+
vtx3body::ITSTrackDCAZToSVXY);
269+
261270
// MC candidate table for analysis
262271
DECLARE_SOA_TABLE(McVtx3BodyDatas, "AOD", "MC3BODYDATA", //!
263272
o2::soa::Index<>,

‎PWGLF/TableProducer/Nuspex/trackedHypertritonRecoTask.cxx‎

Lines changed: 173 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -100,6 +100,7 @@ struct TrackedHypertritonRecoTask {
100100
Produces<aod::MCHypCands> mcHypCands;
101101
Produces<aod::Vtx3BodyDatas> vtx3BodyDatas;
102102
Produces<aod::Vtx3BodyCovs> vtx3BodyCovs;
103+
Produces<aod::Vtx3BodyTrackedInfo> vtx3BodyTrackedInfo;
103104
Produces<aod::McVtx3BodyDatas> mcVtx3BodyDatas;
104105

105106
Service<o2::ccdb::BasicCCDBManager> ccdb{};
@@ -124,7 +125,9 @@ struct TrackedHypertritonRecoTask {
124125

125126
struct : ConfigurableGroup {
126127
std::string prefix = "twoBody";
128+
Configurable<bool> useKFParticle{"useKFParticle", false, "Use KFParticle to build the two-body candidate"};
127129
Configurable<bool> useSelections{"useSelections", false, "Apply the two-body candidate selections"};
130+
Configurable<bool> kfSetTopologicalConstraint{"kfSetTopologicalConstraint", false, "Set topological vertex constraint in case of KFParticle reconstruction"};
128131
Configurable<float> maxEtaDaughters{"maxEtaDaughters", 1.e10f, "Maximum absolute daughter eta"};
129132
Configurable<float> minTPCNClsHe{"minTPCNClsHe", -1.f, "Minimum He TPC clusters"};
130133
Configurable<float> minTPCNClsPi{"minTPCNClsPi", -1.f, "Minimum pion TPC clusters"};
@@ -133,7 +136,7 @@ struct TrackedHypertritonRecoTask {
133136
Configurable<float> minTPCInnerParamHe{"minTPCInnerParamHe", -1.f, "Minimum He TPC rigidity"};
134137
Configurable<float> minPt{"minPt", -1.f, "Minimum candidate transverse momentum"};
135138
Configurable<float> massWindow{"massWindow", 1.e10f, "Half-width of the hypertriton mass window"};
136-
Configurable<float> maxDcaDaughters{"maxDcaDaughters", 1.e10f, "Maximum DCA between daughters"};
139+
Configurable<float> maxChi2{"maxChi2", 1.e10f, "KFParticle: Maximum SV chi2, DCA fitter: Maximum DCA between daughters"};
137140
Configurable<float> minCosPA{"minCosPA", -2.f, "Minimum cosine of the pointing angle"};
138141
Configurable<float> minDcaHeToPV{"minDcaHeToPV", -1.f, "Minimum absolute He DCA to PV"};
139142
Configurable<float> minDcaPiToPV{"minDcaPiToPV", -1.f, "Minimum absolute pion DCA to PV"};
@@ -179,6 +182,8 @@ struct TrackedHypertritonRecoTask {
179182
Configurable<float> maxCtau{"maxCtau", 100.f, "Maximum candidate c tau"};
180183
Configurable<float> minCosPA{"minCosPA", 0.9f, "Minimum candidate cosine of the pointing angle"};
181184
Configurable<float> maxChi2{"maxChi2", 100.f, "Maximum candidate chi2"};
185+
Configurable<float> maxITSDCAxytrackToSV{"maxITSDCAxytrackToSV", 10.0, "Maximum distance of ITS matched track to SV in xy"};
186+
Configurable<float> maxITSDCAztrackToSV{"maxITSDCAztrackToSV", 10.0, "Maximum distance of ITS matched track to SV in z"};
182187
} threeBody;
183188

184189
o2::vertexing::DCAFitterN<2> fitter2Body;
@@ -196,6 +201,22 @@ struct TrackedHypertritonRecoTask {
196201
std::vector<int> recoCollisionForMC;
197202
std::vector<bool> survivedMCEventSelection;
198203

204+
struct v0Candidate {
205+
// daughter properties
206+
std::array<float, 3> momHelium{};
207+
std::array<float, 3> momPion{};
208+
std::array<float, 3> posHelium{};
209+
std::array<float, 3> posPion{};
210+
// vertex properties
211+
float mass;
212+
float chi2;
213+
float cosPA;
214+
std::array<float, 3> decayVertex{};
215+
std::array<float, 3> momentum{};
216+
};
217+
218+
v0Candidate v0;
219+
199220
struct TwoBodyMCInfo {
200221
float genPt = -1.f;
201222
float genPhi = -1.f;
@@ -536,14 +557,101 @@ struct TrackedHypertritonRecoTask {
536557
return (track.tpcSignal() - expected) / (expected * bbParamsHe[5]);
537558
}
538559

539-
template <typename TTrack, typename TCollision, typename TFillCandidate>
540-
void buildTwoBody(TTrack const& heTrack, TTrack const& piTrack, TCollision const& collision, float trackedClSize, TFillCandidate const& fillCandidate)
560+
template <typename TCollision, typename TTrack, typename TTrackParCov>
561+
void fit2BodyWithKF(TCollision const& collision,
562+
TTrack const& trackHelium,
563+
TTrack const& trackPion,
564+
TTrackParCov const& trackHeliumCov,
565+
TTrackParCov const& trackPionCov)
566+
{
567+
// initialise KF primary vertex
568+
KFVertex kfpVertex = createKFPVertexFromCollision(collision);
569+
KFParticle kfpv(kfpVertex);
570+
571+
// create KFParticle objects
572+
KFParticle kfpHelium, kfpPion;
573+
// helium
574+
std::array<float, 3> xyz, pxpypz;
575+
float xyzpxpypz[6];
576+
trackHeliumCov.getPxPyPzGlo(pxpypz);
577+
trackHeliumCov.getXYZGlo(xyz);
578+
for (int i = 0; i < 3; ++i) {
579+
xyzpxpypz[i] = xyz[i];
580+
xyzpxpypz[i + 3] = pxpypz[i] * 2;
581+
}
582+
std::array<float, 21> cv{};
583+
trackHeliumCov.getCovXYZPxPyPzGlo(cv);
584+
KFParticle kfHelium;
585+
kfHelium.Create(xyzpxpypz, cv.data(), trackHelium.sign() * 2, constants::physics::MassHelium3);
586+
// pion
587+
kfpPion = createKFParticleFromTrackParCov(trackPionCov, trackPion.sign(), constants::physics::MassPionCharged);
588+
589+
// construct V0 vertex
590+
KFParticle KFV0;
591+
int nDaughtersV0 = 2;
592+
const KFParticle* DaughtersV0[2] = {&kfpHelium, &kfpPion};
593+
KFV0.SetConstructMethod(2);
594+
try {
595+
KFV0.Construct(DaughtersV0, nDaughtersV0);
596+
} catch (std::runtime_error& e) {
597+
LOG(debug) << "Failed to create V0 vertex." << e.what();
598+
return;
599+
}
600+
601+
// topological constraint
602+
if (twoBody.kfSetTopologicalConstraint) {
603+
KFV0.SetProductionVertex(kfpv);
604+
KFV0.TransportToDecayVertex();
605+
}
606+
607+
// get vertex position and momentum
608+
v0.decayVertex[0] = KFV0.GetX();
609+
v0.decayVertex[1] = KFV0.GetY();
610+
v0.decayVertex[2] = KFV0.GetZ();
611+
v0.momentum[0] = KFV0.GetPx();
612+
v0.momentum[1] = KFV0.GetPy();
613+
v0.momentum[2] = KFV0.GetPz();
614+
615+
// transport all daughter tracks to hypertriton vertex
616+
// float position[3];
617+
// for (int i; i < 3; i++) {
618+
// position[i] = v0.decayVertex[i];
619+
// }
620+
kfpHelium.TransportToPoint(v0.decayVertex.data());
621+
kfpPion.TransportToPoint(v0.decayVertex.data());
622+
623+
// daughter positions
624+
v0.posHelium[0] = kfpHelium.GetX();
625+
v0.posHelium[1] = kfpHelium.GetY();
626+
v0.posHelium[2] = kfpHelium.GetZ();
627+
v0.posPion[0] = kfpPion.GetX();
628+
v0.posPion[1] = kfpPion.GetY();
629+
v0.posPion[2] = kfpPion.GetZ();
630+
631+
// daughter momenta
632+
v0.momHelium[0] = kfpHelium.GetPx();
633+
v0.momHelium[1] = kfpHelium.GetPy();
634+
v0.momHelium[2] = kfpHelium.GetPz();
635+
v0.momPion[0] = kfpPion.GetPx();
636+
v0.momPion[1] = kfpPion.GetPy();
637+
v0.momPion[2] = kfpPion.GetPz();
638+
639+
// candidate mass
640+
float mass, massErr;
641+
KFV0.GetMass(mass, massErr);
642+
v0.mass = mass;
643+
644+
// vertex chi2
645+
v0.chi2 = KFV0.GetChi2() / KFV0.GetNDF();
646+
}
647+
648+
template <typename TTrackParCov>
649+
void fit2bodyWithDCAFitter(TTrackParCov const& trackHeliumCov,
650+
TTrackParCov const& trackPionCov)
541651
{
542-
auto heTrackCov = getTrackParCov(heTrack);
543-
auto piTrackCov = getTrackParCov(piTrack);
544652
int nCandidates = 0;
545653
try {
546-
nCandidates = fitter2Body.process(heTrackCov, piTrackCov);
654+
nCandidates = fitter2Body.process(trackHeliumCov, trackPionCov);
547655
} catch (...) {
548656
LOG(error) << "Exception while fitting a tracked two-body candidate";
549657
return;
@@ -552,36 +660,57 @@ struct TrackedHypertritonRecoTask {
552660
return;
553661
}
554662

555-
std::array<float, 3> heMomentum{};
556-
std::array<float, 3> piMomentum{};
557-
fitter2Body.getTrack(0).getPxPyPzGlo(heMomentum);
558-
fitter2Body.getTrack(1).getPxPyPzGlo(piMomentum);
559-
for (std::size_t i = 0; i < heMomentum.size(); ++i) {
560-
heMomentum[i] *= 2.f;
561-
}
562-
std::array<float, 3> momentum{heMomentum[0] + piMomentum[0], heMomentum[1] + piMomentum[1], heMomentum[2] + piMomentum[2]};
563-
if (twoBody.useSelections && std::hypot(momentum[0], momentum[1]) < twoBody.minPt) {
564-
return;
663+
// get daughter momenta
664+
fitter2Body.getTrack(0).getPxPyPzGlo(v0.momHelium);
665+
fitter2Body.getTrack(1).getPxPyPzGlo(v0.momPion);
666+
667+
for (std::size_t i = 0; i < v0.momHelium.size(); ++i) {
668+
v0.momHelium[i] *= 2.f;
565669
}
566670

567-
const float heMomentum2 = RecoDecay::sumOfSquares(heMomentum[0], heMomentum[1], heMomentum[2]);
568-
const float piMomentum2 = RecoDecay::sumOfSquares(piMomentum[0], piMomentum[1], piMomentum[2]);
569-
const float candidateMomentum2 = RecoDecay::sumOfSquares(momentum[0], momentum[1], momentum[2]);
671+
// compute candidate momentum
672+
v0.momentum = {v0.momHelium[0] + v0.momPion[0], v0.momHelium[1] + v0.momPion[1], v0.momHelium[2] + v0.momPion[2]};
673+
674+
// compute candidate mass
675+
const float heMomentum2 = RecoDecay::sumOfSquares(v0.momHelium[0], v0.momHelium[1], v0.momHelium[2]);
676+
const float piMomentum2 = RecoDecay::sumOfSquares(v0.momPion[0], v0.momPion[1], v0.momPion[2]);
677+
const float candidateMomentum2 = RecoDecay::sumOfSquares(v0.momentum[0], v0.momentum[1], v0.momentum[2]);
570678
const float heEnergy = std::sqrt(heMomentum2 + constants::physics::MassHelium3 * constants::physics::MassHelium3);
571679
const float piEnergy = std::sqrt(piMomentum2 + constants::physics::MassPionCharged * constants::physics::MassPionCharged);
572-
const float mass = std::sqrt((heEnergy + piEnergy) * (heEnergy + piEnergy) - candidateMomentum2);
573-
if (twoBody.useSelections && std::abs(mass - constants::physics::MassHyperTriton) > twoBody.massWindow) {
574-
return;
575-
}
680+
v0.mass = std::sqrt((heEnergy + piEnergy) * (heEnergy + piEnergy) - candidateMomentum2);
576681

682+
// get SV position
577683
const auto& secondaryVertex = fitter2Body.getPCACandidate();
684+
for (int i = 0; i < 3; i++) {
685+
v0.decayVertex[i] = secondaryVertex[i];
686+
}
687+
v0.chi2 = std::sqrt(fitter2Body.getChi2AtPCACandidate());
688+
}
689+
690+
template <typename TTrack, typename TCollision, typename TFillCandidate>
691+
void buildTwoBody(TTrack const& heTrack, TTrack const& piTrack, TCollision const& collision, float trackedClSize, TFillCandidate const& fillCandidate)
692+
{
578693
const std::array<float, 3> primaryVertex{collision.posX(), collision.posY(), collision.posZ()};
579-
const std::array<float, 3> decayVertex{static_cast<float>(secondaryVertex[0]), static_cast<float>(secondaryVertex[1]), static_cast<float>(secondaryVertex[2])};
580-
const float dcaDaughters = std::sqrt(fitter2Body.getChi2AtPCACandidate());
581-
if (twoBody.useSelections && (dcaDaughters > twoBody.maxDcaDaughters || RecoDecay::cpa(primaryVertex, decayVertex, momentum) < twoBody.minCosPA)) {
694+
695+
auto heTrackCov = getTrackParCov(heTrack);
696+
auto piTrackCov = getTrackParCov(piTrack);
697+
698+
if (twoBody.useKFParticle) {
699+
fit2BodyWithKF(collision, heTrack, piTrack, heTrackCov, piTrackCov);
700+
} else {
701+
fit2bodyWithDCAFitter(heTrackCov, piTrackCov);
702+
}
703+
704+
v0.cosPA = RecoDecay::cpa(primaryVertex, v0.decayVertex, v0.momentum);
705+
if (twoBody.useSelections &&
706+
(std::hypot(v0.momentum[0], v0.momentum[1]) < twoBody.minPt ||
707+
std::abs(v0.mass - constants::physics::MassHyperTriton) > twoBody.massWindow ||
708+
v0.chi2 > twoBody.maxChi2 ||
709+
v0.cosPA < twoBody.minCosPA)) {
582710
return;
583711
}
584712

713+
// Do propagation with Propagator including material interactions in all cases (KFParticle propagation would not include material)
585714
std::array<float, 2> dcaInfo{};
586715
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, heTrackCov, 2.f, fitter2Body.getMatCorrType(), &dcaInfo);
587716
const float dcaHe = dcaInfo[0];
@@ -605,10 +734,10 @@ struct TrackedHypertritonRecoTask {
605734
fillCandidate(collision.centFT0A(), collision.centFT0C(), collision.centFT0M(),
606735
collision.posX(), collision.posY(), collision.posZ(),
607736
runNumber, heTrack.sign() > 0,
608-
std::hypot(heMomentum[0], heMomentum[1]), std::atan2(heMomentum[1], heMomentum[0]), RecoDecay::eta(heMomentum),
609-
std::hypot(piMomentum[0], piMomentum[1]), std::atan2(piMomentum[1], piMomentum[0]), RecoDecay::eta(piMomentum),
610-
secondaryVertex[0], secondaryVertex[1], secondaryVertex[2],
611-
dcaDaughters, dcaHe, dcaPi,
737+
std::hypot(v0.momHelium[0], v0.momHelium[1]), std::atan2(v0.momHelium[1], v0.momHelium[0]), RecoDecay::eta(v0.momHelium),
738+
std::hypot(v0.momPion[0], v0.momPion[1]), std::atan2(v0.momPion[1], v0.momPion[0]), RecoDecay::eta(v0.momPion),
739+
v0.decayVertex[0], v0.decayVertex[1], v0.decayVertex[2],
740+
v0.chi2, dcaHe, dcaPi,
612741
nSigmaHe3(heTrack), heTrack.tpcNClsFound(), piTrack.tpcNClsFound(),
613742
static_cast<int16_t>(heTrack.tpcNClsFindable()) - heTrack.tpcNClsFindableMinusPID(),
614743
static_cast<int16_t>(piTrack.tpcNClsFindable()) - piTrack.tpcNClsFindableMinusPID(),
@@ -662,6 +791,7 @@ struct TrackedHypertritonRecoTask {
662791
static_cast<int>(candidate.tpcNCl[0]), static_cast<int>(candidate.tpcNCl[1]), static_cast<int>(candidate.tpcNCl[2]),
663792
static_cast<uint32_t>(candidate.pidForTrackingDeuteron));
664793
vtx3BodyCovs(candidate.covProton, candidate.covPion, candidate.covDeuteron, candidate.covariance);
794+
vtx3BodyTrackedInfo(candidate.itsTrackDCAToSV[0], candidate.itsTrackDCAToSV[1]);
665795
}
666796

667797
void fillThreeBodyMCTable(ThreeBodyMCInfo const& info)
@@ -776,12 +906,25 @@ struct TrackedHypertritonRecoTask {
776906
if (decay3Body.collisionId() < 0 || !goodCollision[decay3Body.collisionId()] || (skimmedProcessing && !zorroDecision[decay3Body.collisionId()][kTracked3Body])) {
777907
continue;
778908
}
909+
779910
const auto collision = decay3Body.collision_as<Collisions>();
780911
const auto trackPositive = decay3Body.track0_as<Tracks>();
781912
const auto trackNegative = decay3Body.track1_as<Tracks>();
782913
const auto trackDeuteron = decay3Body.track2_as<Tracks>();
783914
const auto trackProton = trackDeuteron.sign() > 0 ? trackPositive : trackNegative;
784915
const auto trackPion = trackDeuteron.sign() > 0 ? trackNegative : trackPositive;
916+
917+
// get DCA of ITS track to SV
918+
const auto itsTrack = tracked3Body.itsTrack_as<Tracks>();
919+
auto itsTrackParCov = getTrackParCov(itsTrack);
920+
std::array<float, 2> dcaInfoItsTrack{};
921+
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, itsTrackParCov, 2.f, fitter2Body.getMatCorrType(), &dcaInfoItsTrack);
922+
builder3Body.decay3body.itsTrackDCAToSV[0] = dcaInfoItsTrack[0];
923+
builder3Body.decay3body.itsTrackDCAToSV[1] = dcaInfoItsTrack[1];
924+
if (threeBody.useSelections && (builder3Body.decay3body.itsTrackDCAToSV[0] > threeBody.maxITSDCAxytrackToSV || builder3Body.decay3body.itsTrackDCAToSV[1] > threeBody.maxITSDCAztrackToSV)) {
925+
continue;
926+
}
927+
785928
if (builder3Body.buildDecay3BodyCandidate(collision, trackProton, trackPion, trackDeuteron,
786929
decay3Body.globalIndex(), deuteronTOFNSigma(collision, trackDeuteron), tracked3Body.itsClsSize(),
787930
threeBody.useKFParticle, threeBody.setTopologicalConstraint,

‎PWGLF/Utils/decay3bodyBuilderHelper.h‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -91,6 +91,9 @@ struct decay3bodyCandidate {
9191
float covPion[21] = {0.0f};
9292
float covDeuteron[21] = {0.0f};
9393
float covariance[21] = {0.0f};
94+
95+
// strangeness tracking
96+
std::array<float, 2> itsTrackDCAToSV;
9497
};
9598

9699
//_______________________________________________________________________

0 commit comments

Comments
 (0)