From aeeec7b9b1c779cb4703cba032db9f7fa00a0622 Mon Sep 17 00:00:00 2001
From: "s.zharko@gsi.de" <s.zharko@gsi.de>
Date: Thu, 13 Feb 2025 03:39:52 +0100
Subject: [PATCH] online: V0Finder implementation
---
algo/ca/core/data/CaTrack.h | 8 +-
algo/ca/core/utils/CaMonitorData.h | 4 +
algo/data/sts/Hit.h | 15 +-
algo/detectors/tof/Hit.h | 15 +-
algo/evselector/RecoEventSelectorMonitor.h | 2 +
algo/global/Reco.cxx | 7 +-
algo/global/RecoResults.h | 8 +-
algo/kfp/KfpV0Finder.cxx | 398 +++++++++++++++++++--
algo/kfp/KfpV0Finder.h | 123 ++++++-
algo/kfp/KfpV0FinderConfig.cxx | 6 +
algo/kfp/KfpV0FinderMonitor.h | 24 +-
11 files changed, 540 insertions(+), 70 deletions(-)
diff --git a/algo/ca/core/data/CaTrack.h b/algo/ca/core/data/CaTrack.h
index 29b7f3859d..cf46636107 100644
--- a/algo/ca/core/data/CaTrack.h
+++ b/algo/ca/core/data/CaTrack.h
@@ -27,6 +27,8 @@ namespace cbm::algo::ca
///
class Track {
public:
+ using TrackParam_t = cbm::algo::kf::TrackParamS;
+
friend class boost::serialization::access;
Track() = default;
@@ -43,9 +45,9 @@ namespace cbm::algo::ca
public:
int fNofHits{kfdefs::Undef<int>}; ///< Number of hits in track
- cbm::algo::kf::TrackParamS fParFirst; ///< Track parameters on the first station
- cbm::algo::kf::TrackParamS fParLast; ///< Track parameters on the last station
- cbm::algo::kf::TrackParamS fParPV; ///< Track parameters in the primary vertex
+ TrackParam_t fParFirst; ///< Track parameters on the first station
+ TrackParam_t fParLast; ///< Track parameters on the last station
+ TrackParam_t fParPV; ///< Track parameters in the primary vertex
};
} // namespace cbm::algo::ca
diff --git a/algo/ca/core/utils/CaMonitorData.h b/algo/ca/core/utils/CaMonitorData.h
index 09b3d0e117..5d93816c1f 100644
--- a/algo/ca/core/utils/CaMonitorData.h
+++ b/algo/ca/core/utils/CaMonitorData.h
@@ -82,6 +82,10 @@ namespace cbm::algo::ca
/// \brief Resets all the counters and timers
void Reset();
+ /// \brief Resets a particular counter
+ /// \param key Counter key
+ void ResetCounter(ECounterKey key) { faCounters[key] = 0; }
+
/// \brief Starts timer
/// \param key Timer key
void StartTimer(ETimerKey key) { faTimers[key].Start(); }
diff --git a/algo/data/sts/Hit.h b/algo/data/sts/Hit.h
index ab1a733f2c..2df0d8370f 100644
--- a/algo/data/sts/Hit.h
+++ b/algo/data/sts/Hit.h
@@ -26,14 +26,13 @@ namespace cbm::algo::sts
u32 fBackClusterId; ///< Index of back-side cluster, used by tracking to reduce combinatorics
// Interface for tracking
- real X() const { return fX; }
- real Y() const { return fY; }
- real Z() const { return fZ; }
- u32 Time() const { return fTime; }
-
- real Dx() const { return fDx; }
- real Dy() const { return fDy; }
- real TimeError() const { return fTimeError; }
+ double X() const { return fX; }
+ double Y() const { return fY; }
+ double Z() const { return fZ; }
+ double Time() const { return fTime; }
+ double Dx() const { return fDx; }
+ double Dy() const { return fDy; }
+ double TimeError() const { return fTimeError; }
private: // serialization
friend class boost::serialization::access;
diff --git a/algo/detectors/tof/Hit.h b/algo/detectors/tof/Hit.h
index a7608fc72e..8e3a1a0685 100644
--- a/algo/detectors/tof/Hit.h
+++ b/algo/detectors/tof/Hit.h
@@ -53,14 +53,13 @@ namespace cbm::algo::tof
// Interface for tracker
- real X() const { return hitPos.X(); }
- real Y() const { return hitPos.Y(); }
- real Z() const { return hitPos.Z(); }
- u32 Time() const { return hitTime; }
-
- real Dx() const { return hitPosErr.X(); }
- real Dy() const { return hitPosErr.Y(); }
- real TimeError() const { return hitTimeErr; }
+ double X() const { return hitPos.X(); }
+ double Y() const { return hitPos.Y(); }
+ double Z() const { return hitPos.Z(); }
+ double Time() const { return hitTime; }
+ double Dx() const { return hitPosErr.X(); }
+ double Dy() const { return hitPosErr.Y(); }
+ double TimeError() const { return hitTimeErr; }
// Interface end
diff --git a/algo/evselector/RecoEventSelectorMonitor.h b/algo/evselector/RecoEventSelectorMonitor.h
index 82c68535c7..de2eb4a27e 100644
--- a/algo/evselector/RecoEventSelectorMonitor.h
+++ b/algo/evselector/RecoEventSelectorMonitor.h
@@ -38,6 +38,7 @@ namespace cbm::algo::evselect
TofHitFinder,
TrdHitFinder,
TrackFinder,
+ V0Finder,
END
};
/* clang-format on */
@@ -65,6 +66,7 @@ namespace cbm::algo::evselect
SetTimerName(ETimer::TofHitFinder, "hit finding in TOF");
SetTimerName(ETimer::TrdHitFinder, "hit finding in TRD");
SetTimerName(ETimer::TrackFinder, "track finding");
+ SetTimerName(ETimer::V0Finder, "V0 finding");
}
private:
diff --git a/algo/global/Reco.cxx b/algo/global/Reco.cxx
index 37c37b5088..b2a7355599 100644
--- a/algo/global/Reco.cxx
+++ b/algo/global/Reco.cxx
@@ -597,9 +597,10 @@ bool Reco::ReconstructEvent(const DigiEvent& digiEvent)
}
//* V0-selector
- {
- auto triggers = fV0Finder->ProcessEvent(recoEvent);
- }
+ fEvSelectingMonitor.StartTimer(evselect::ETimer::V0Finder);
+ auto triggers = fV0Finder->ProcessEvent(recoEvent);
+ fEvSelectingMonitor.StopTimer(evselect::ETimer::V0Finder);
+
fEvSelectingMonitor.IncrementCounter(evselect::ECounter::EventsSelected);
return true;
diff --git a/algo/global/RecoResults.h b/algo/global/RecoResults.h
index 7795dda396..3bcee8be6a 100644
--- a/algo/global/RecoResults.h
+++ b/algo/global/RecoResults.h
@@ -27,6 +27,8 @@ namespace cbm::algo
/// @name RecoResults
/// @brief Structure to keep reconstructed results: digi-events, hits and tracks
struct RecoResults {
+ using HitId_t = std::pair<uint32_t, uint32_t>; // Hit ID by track
+
PODVector<CbmBmonDigi> bmonDigis;
PODVector<CbmStsDigi> stsDigis;
PODVector<CbmMuchDigi> muchDigis;
@@ -45,8 +47,8 @@ namespace cbm::algo
PartitionedVector<bmon::Hit> bmonHits;
ca::Vector<ca::Track> tracks;
- ca::Vector<std::vector<std::pair<uint32_t, uint32_t>>> trackStsHitIndices; // [trk][hit][(iPart, iHit)]
- ca::Vector<std::vector<std::pair<uint32_t, uint32_t>>> trackTofHitIndices; // [trk][hit][(iPart, iHit)]
- ca::Vector<std::vector<std::pair<uint32_t, uint32_t>>> trackTrdHitIndices; // [trk][hit][(iPart, iHit)]
+ ca::Vector<std::vector<HitId_t>> trackStsHitIndices; // [trk][hit][(iPart, iHit)]
+ ca::Vector<std::vector<HitId_t>> trackTofHitIndices; // [trk][hit][(iPart, iHit)]
+ ca::Vector<std::vector<HitId_t>> trackTrdHitIndices; // [trk][hit][(iPart, iHit)]
};
} // namespace cbm::algo
diff --git a/algo/kfp/KfpV0Finder.cxx b/algo/kfp/KfpV0Finder.cxx
index 5a5ddbc6e8..2488ba89e2 100644
--- a/algo/kfp/KfpV0Finder.cxx
+++ b/algo/kfp/KfpV0Finder.cxx
@@ -18,32 +18,99 @@
using cbm::algo::RecoResults;
using cbm::algo::kfp::V0Finder;
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+bool V0Finder::AssignMomentum(const PartitionedVector<tof::Hit>& tofHits,
+ const std::vector<RecoResults::HitId_t>& tofHitIds, double t0, ParticleInfo& particleInfo)
+{
+ if (tofHitIds.empty()) {
+ fEventMonitor.IncrementCounter(ECounter::TracksWoTofHits);
+ return false;
+ }
+
+ double beta{0.};
+ if constexpr (kUseAverageSpeed) {
+ for (const auto& hitId : tofHitIds) {
+ beta += EstimateBeta(tofHits[hitId.first][hitId.second], t0);
+ }
+ beta /= tofHitIds.size();
+ }
+ else {
+ const auto& hitId = tofHitIds.back();
+ beta = EstimateBeta(tofHits[hitId.first][hitId.second], t0);
+ }
+ if (beta < 0.) {
+ fEventMonitor.IncrementCounter(ECounter::TracksWNegativeTofHitTime);
+ return false;
+ }
+ else if (beta > 1.) {
+ fEventMonitor.IncrementCounter(ECounter::TracksWUnphysicalBeta);
+ return false;
+ }
+ double gamma{1. / sqrt(1. - beta * beta)};
+ particleInfo.fBeta = beta;
+ particleInfo.fQp = particleInfo.fQp / (gamma * beta * particleInfo.fMass);
+ return true;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void V0Finder::AssignPid(ParticleInfo& particleInfo)
+{
+ if (std::isnan(particleInfo.fDca)) {
+ particleInfo.fPdg = kUndefPdg;
+ fEventMonitor.IncrementCounter(ECounter::TracksWoPid);
+ }
+ else if (particleInfo.fDca > fMinPionDca) {
+ // pi-
+ particleInfo.fPdg = -211;
+ particleInfo.fMass = kPionMass;
+ particleInfo.fCharge = -1;
+ fEventMonitor.IncrementCounter(ECounter::PionsDca);
+ }
+ else if (particleInfo.fDca > fMinProtonDca) {
+ // proton
+ particleInfo.fPdg = 2212;
+ particleInfo.fMass = kProtonMass;
+ particleInfo.fCharge = 1;
+ fEventMonitor.IncrementCounter(ECounter::ProtonsDca);
+ }
+ else {
+ // primary
+ //particleInfo.fPdg = fPrimaryAssignedPdg;
+ particleInfo.fPdg = kUndefPdg;
+ fEventMonitor.IncrementCounter(ECounter::TracksWoPid);
+ }
+}
+
// ---------------------------------------------------------------------------------------------------------------------
//
-std::vector<double> V0Finder::CollectDca(const RecoResults& recoEvent) const
+void V0Finder::CollectDca(const RecoResults& recoEvent)
{
- std::vector<double> dca(recoEvent.trackStsHitIndices.size(), std::numeric_limits<double>::signaling_NaN());
const auto& stsHitIndices = recoEvent.trackStsHitIndices;
for (size_t iTrk = 0; iTrk < stsHitIndices.size(); ++iTrk) {
const auto& stsHitIndicesInTrack = stsHitIndices[iTrk];
if (stsHitIndicesInTrack.size() < 2) { // less then two sts hits
+ fEventMonitor.IncrementCounter(ECounter::TracksWoStsHits);
continue;
}
- auto [iPtFst, iHitFst] = stsHitIndicesInTrack[0];
- auto [iPtSnd, iHitSnd] = stsHitIndicesInTrack[1];
- dca[iTrk] = EstimateDca(recoEvent.stsHits[iPtFst][iHitFst], recoEvent.stsHits[iPtSnd][iHitSnd]);
+ auto& particleInfo = fvParticleInfo[iTrk];
+ auto [iPtFst, iHitFst] = stsHitIndicesInTrack[0];
+ auto [iPtSnd, iHitSnd] = stsHitIndicesInTrack[1];
+ fvParticleInfo[iTrk].fDca = EstimateDca(recoEvent.stsHits[iPtFst][iHitFst], recoEvent.stsHits[iPtSnd][iHitSnd]);
+ AssignPid(particleInfo);
}
- return dca;
}
// ---------------------------------------------------------------------------------------------------------------------
//
-std::vector<double> V0Finder::CollectT0(gsl::span<const bmon::Hit> bmonHits) const
+void V0Finder::CollectT0(gsl::span<const bmon::Hit> bmonHits)
{
- std::vector<double> t0;
- t0.reserve(bmonHits.size());
- std::transform(bmonHits.begin(), bmonHits.end(), std::back_inserter(t0), [&](const auto& h) { return h.GetTime(); });
- return t0;
+ fvT0s.clear();
+ fvT0s.reserve(bmonHits.size());
+ std::transform(bmonHits.begin(), bmonHits.end(), std::back_inserter(fvT0s),
+ [&](const auto& h) { return h.GetTime(); });
}
// ---------------------------------------------------------------------------------------------------------------------
@@ -58,38 +125,317 @@ double V0Finder::EstimateDca(const sts::Hit& fst, const sts::Hit& snd) const
// ---------------------------------------------------------------------------------------------------------------------
//
-void V0Finder::Init() {}
+double V0Finder::EstimateBeta(const tof::Hit& hit, double t0) const
+{
+ double t = hit.Time() - t0 - fTzeroOffset;
+ double x{hit.X() - fOrigin[0]};
+ double y{hit.Y() - fOrigin[1]};
+ double z{hit.Y() - fOrigin[2]};
+ double x2{x * x};
+ double y2{y * y};
+ double z2{z * z};
+ double r2{x2 + y2 + z2};
+ return std::sqrt(r2) / (t * kSpeedOfLight);
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+bool V0Finder::FindV0Candidates(const RecoResults& recoEvent, double t0)
+{
+ const auto& tracks = recoEvent.tracks;
+
+ // Reset temporary data structures
+ InitTrackParamVectors(tracks);
+ fpTopoReconstructor->Clear();
+ fvSelectedTrackIds.clear();
+ fvSelectedTrackIds.reserve(tracks.size());
+ fEventMonitor.ResetCounter(ECounter::TracksWoMomentum);
+ fEventMonitor.ResetCounter(ECounter::TracksSelected);
+ fEventMonitor.ResetCounter(ECounter::Pions);
+ fEventMonitor.ResetCounter(ECounter::Protons);
+ fEventMonitor.ResetCounter(ECounter::EventsLambdaCand);
+
+ // Preselect tracks
+ uint32_t nProtonCandidates{0};
+ uint32_t nPionCandidates{0};
+ uint32_t nSelectedTracks{0};
+ fEventMonitor.StartTimer(ETimer::PreselectTracks);
+ for (size_t iTrk = 0; iTrk < tracks.size(); ++iTrk) { // Over all tracks
+ auto& particleInfo = fvParticleInfo[iTrk];
+
+ // Cut tracks with undefined dca (by PID == -2)
+ // NOTE: if fPdg == kUndefPdg, beta and QP were not estimated for this track on previous iterations as well
+ if (particleInfo.fPdg == kUndefPdg) {
+ continue;
+ }
+
+ // Reset fields of the ParticleInfo, which could be filled on the previous iteration
+ particleInfo.fBeta = std::numeric_limits<double>::quiet_NaN();
+ particleInfo.fQp = std::numeric_limits<double>::quiet_NaN();
+ particleInfo.fbSelected = false;
+
+ // Assign momentum to tracks
+ if (!AssignMomentum(recoEvent.tofHits, recoEvent.trackTofHitIndices[iTrk], t0, particleInfo)) {
+ fEventMonitor.IncrementCounter(ECounter::TracksWoMomentum);
+ continue; // No momentum was assigned
+ }
+
+ // Select tracks
+ if (!SelectTrack(particleInfo)) {
+ continue;
+ }
+ fvSelectedTrackIds.push_back(iTrk);
+ particleInfo.fbSelected = true;
+
+ // Update track parameters
+ double qpVar{particleInfo.fQp * fQpAssignedUncertainty};
+ qpVar = qpVar * qpVar;
+ auto& trkParam{fvTrackParam[iTrk]};
+ trkParam.first.SetQp(particleInfo.fQp);
+ trkParam.first.SetC44(qpVar);
+ trkParam.second.SetQp(particleInfo.fQp);
+ trkParam.second.SetC44(qpVar);
+
+ switch (particleInfo.fPdg) {
+ case -211: ++nPionCandidates; break;
+ case 2212: ++nProtonCandidates; break;
+ }
+ ++nSelectedTracks;
+ }
+ fEventMonitor.StopTimer(ETimer::PreselectTracks);
+ if (!nPionCandidates || !nProtonCandidates) {
+ return false; // no Lambda can be found
+ }
+ fEventMonitor.IncrementCounter(ECounter::EventsLambdaCand);
+ fEventMonitor.IncrementCounter(ECounter::TracksSelected, nSelectedTracks);
+ fEventMonitor.IncrementCounter(ECounter::Pions, nPionCandidates);
+ fEventMonitor.IncrementCounter(ECounter::Protons, nProtonCandidates);
+
+ // Initialize and run the KFParticleFinder
+ fEventMonitor.StartTimer(ETimer::InitKfp);
+ KFPTrackVector kfpTracksFst;
+ KFPTrackVector kfpTracksLst;
+ kfpTracksFst.Resize(nSelectedTracks);
+ kfpTracksLst.Resize(nSelectedTracks);
+ fpTopoReconstructor->Init(kfpTracksFst, kfpTracksLst);
+ for (uint32_t iKfpTrk = 0; iKfpTrk < fvSelectedTrackIds.size(); ++iKfpTrk) { // Over selected tracks
+ uint32_t iCaTrk{fvSelectedTrackIds[iKfpTrk]};
+ const auto& trkParam{fvTrackParam[iCaTrk]};
+ const auto& particleInfo{fvParticleInfo[iCaTrk]};
+ SetKfpTrackParameters(kfpTracksFst, iKfpTrk, iCaTrk, trkParam.first, particleInfo);
+ SetKfpTrackParameters(kfpTracksLst, iKfpTrk, iCaTrk, trkParam.second, particleInfo);
+ }
+ fpTopoReconstructor->AddPV(MakeKfpPrimaryVertex(fOrigin));
+ fpTopoReconstructor->SortTracks();
+ fEventMonitor.StopTimer(ETimer::InitKfp);
+ fEventMonitor.StartTimer(ETimer::ExecKfp);
+ fpTopoReconstructor->ReconstructParticles();
+ fEventMonitor.StopTimer(ETimer::ExecKfp);
+ const auto& particles{fpTopoReconstructor->GetParticles()};
+
+ // Scan for Lambda-candidates
+ uint32_t nLambdaCandidates =
+ std::count_if(particles.begin(), particles.end(), [&](const auto& p) { return p.GetPDG() == 3122; });
+ fEventMonitor.IncrementCounter(ECounter::KfpLambdaCandidates, nLambdaCandidates);
+
+ return static_cast<bool>(nLambdaCandidates);
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void V0Finder::Init() { fpTopoReconstructor->SetTarget({float(fOrigin[0]), float(fOrigin[1]), float(fOrigin[2])}); }
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void V0Finder::InitTrackParamVectors(const ca::Vector<ca::Track>& tracks)
+{
+ fvTrackParam.clear();
+ fvTrackParam.reserve(tracks.size());
+ std::transform(tracks.begin(), tracks.end(), std::back_inserter(fvTrackParam),
+ [&](const auto& t) { return std::make_pair(t.fParFirst, t.fParLast); });
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+KFVertex V0Finder::MakeKfpPrimaryVertex(const std::array<float, 3>& r)
+{
+ KFVertex kfVertex;
+ kfVertex.X() = r[0];
+ kfVertex.Y() = r[1];
+ kfVertex.Z() = r[2];
+ for (int iC = 0; iC < 6; ++iC) {
+ kfVertex.Covariance(iC) = 0.f;
+ }
+ kfVertex.Chi2() = -100.f;
+ return kfVertex;
+}
// ---------------------------------------------------------------------------------------------------------------------
//
CbmEventTriggers V0Finder::Process(const RecoResults& recoEvent)
{
+ //L_(info) << "----------------------------- EVENT --------------";
CbmEventTriggers res;
fEventMonitor.Reset();
- fEventMonitor.StartTimer(ETimer::Event);
+ fEventMonitor.StartTimer(ETimer::ProcessEvent);
fEventMonitor.IncrementCounter(ECounter::EventsTotal);
+ fEventMonitor.IncrementCounter(ECounter::TracksTotal, recoEvent.tracks.size());
+
+ // ----- Initialize data-structures
+ fvParticleInfo.clear();
+ fvParticleInfo.resize(recoEvent.tracks.size());
// ----- Define T0
// So far we cannot preselect a hit from multiple ones, we will be using all of them iteratively to find lambdas
- auto vT0 = CollectT0(recoEvent.bmonHits[fBmonPartitionIndex]);
- if (vT0.empty()) {
+ fEventMonitor.StartTimer(ETimer::CollectT0);
+ CollectT0(recoEvent.bmonHits[fBmonPartitionIndex]);
+ if (fvT0s.empty()) {
fEventMonitor.IncrementCounter(ECounter::EventsWoTzero);
return res;
}
+ fEventMonitor.StopTimer(ETimer::CollectT0);
- //L_(info) << "------- Event: ";
-
- // ----- Estimate DCA of tracks
+ // ----- Estimate DCA of tracks and assign PID
// If a track has less then two STS hits, and undefined DCA value is stored
- std::vector<double> vDca = CollectDca(recoEvent);
+ fEventMonitor.StartTimer(ETimer::CollectDca);
+ CollectDca(recoEvent);
+ fEventMonitor.StopTimer(ETimer::CollectDca);
- // ***** DEBUG: BEGIN
- //for (auto dca: vDca) {
- // L_(info) << "dca=" << dca;
- //}
- // ***** DEBUG: END
+ // ----- Try to find lambdas for different T0
+ fSelectedT0 = std::numeric_limits<double>::quiet_NaN();
+ fEventMonitor.StartTimer(ETimer::FindV0Candidates);
+ for (double t0 : fvT0s) {
+ if (FindV0Candidates(recoEvent, t0)) {
+ fSelectedT0 = t0;
+ res.Set(CbmEventTriggers::ETrigger::Lambda);
+ fEventMonitor.IncrementCounter(ECounter::EventsLambdaCand);
+ break; // Lambda-candidates were found, there is no sense to scan further
+ }
+ }
+ fEventMonitor.StopTimer(ETimer::FindV0Candidates);
+ fEventMonitor.StopTimer(ETimer::ProcessEvent);
+ return res;
+}
+// ---------------------------------------------------------------------------------------------------------------------
+//
+bool V0Finder::SelectTrack(const ParticleInfo& particleInfo) const
+{
+ // Speed cut
+ if (particleInfo.fPdg == -211) {
+ if (particleInfo.fBeta < fMinBetaPion || particleInfo.fBeta > fMaxBetaPion) {
+ return false;
+ }
+ }
+ else if (particleInfo.fPdg == 2212) {
+ if (particleInfo.fBeta < fMinBetaProton || particleInfo.fBeta > fMaxBetaProton) {
+ return false;
+ }
+ }
+ return true;
+}
- fEventMonitor.StopTimer(ETimer::Event);
- return res;
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void V0Finder::SetKfpTrackParameters(KFPTrackVector& trackVector, uint32_t iKfpTrk, uint32_t iCaTrk,
+ const ca::Track::TrackParam_t& trkParam, const ParticleInfo& particleInfo) const
+{
+ // ----- Parameter definition
+ double tx{trkParam.GetTx()};
+ double ty{trkParam.GetTy()};
+ double qp{trkParam.GetQp()};
+ double p{particleInfo.fCharge / qp};
+ double p2{p * p};
+ double t2inv{1. / (1. + tx * tx + ty * ty)};
+ double pz{std::sqrt(t2inv * p2)};
+ double px{tx * pz};
+ double py{ty * pz};
+
+ trackVector.SetParameter(trkParam.GetX(), 0, iKfpTrk);
+ trackVector.SetParameter(trkParam.GetY(), 1, iKfpTrk);
+ trackVector.SetParameter(trkParam.GetZ(), 2, iKfpTrk);
+ trackVector.SetParameter(px, 3, iKfpTrk);
+ trackVector.SetParameter(py, 4, iKfpTrk);
+ trackVector.SetParameter(pz, 5, iKfpTrk);
+
+ // Jacobian matrix for (tx, ty, qp) -> (px, py, pz)
+ std::array<std::array<double, 3>, 3> Jp;
+ Jp[2][0] = -t2inv * px; // d(pz)/d(tx)
+ Jp[2][1] = -t2inv * py; // d(pz)/d(ty)
+ Jp[2][2] = -pz / qp; // d(pz)/d(qp)
+ Jp[0][0] = tx * Jp[2][0] + pz; // d(px)/d(tx)
+ Jp[0][1] = tx * Jp[2][1]; // d(px)/d(ty)
+ Jp[0][2] = tx * Jp[2][2]; // d(px)/d(qp)
+ Jp[1][0] = ty * Jp[2][0]; // d(py)/d(tx)
+ Jp[1][1] = ty * Jp[2][1] + pz; // d(py)/d(ty)
+ Jp[1][2] = ty * Jp[2][2]; // d(py)/d(qp)
+
+
+ // ----- Covariance matrix definition
+ // Position covariance
+ trackVector.SetCovariance(trkParam.C00(), 0, iKfpTrk); // var(x)
+ trackVector.SetCovariance(trkParam.C01(), 1, iKfpTrk); // cov(x, y)
+ trackVector.SetCovariance(trkParam.C11(), 2, iKfpTrk); // var(y)
+
+ // Momentum-position covariances
+ auto MomPosCovariance = [&](const int k, const int l) constexpr->double
+ {
+ double val{0.};
+ const auto& JpA = Jp[k];
+ for (int i = 0; i < 3; ++i) {
+ val += JpA[i] * trkParam.C(i + 2, l);
+ }
+ return val;
+ };
+ trackVector.SetCovariance(MomPosCovariance(0, 0), 6, iKfpTrk); // cov(x, px)
+ trackVector.SetCovariance(MomPosCovariance(0, 1), 7, iKfpTrk); // cov(y, px)
+ trackVector.SetCovariance(MomPosCovariance(1, 0), 10, iKfpTrk); // cov(x, py)
+ trackVector.SetCovariance(MomPosCovariance(1, 1), 11, iKfpTrk); // cov(y, py)
+ trackVector.SetCovariance(MomPosCovariance(2, 0), 15, iKfpTrk); // cov(x, pz)
+ trackVector.SetCovariance(MomPosCovariance(2, 1), 16, iKfpTrk); // cov(y, pz)
+
+ // Momentum covariances
+ auto MomentumCovariance = [&](const int k, const int l) constexpr->double
+ {
+ double val{0.};
+ const auto& JpA = Jp[k];
+ const auto& JpB = Jp[l];
+ for (int i = 0; i < 3; ++i) {
+ double factor{0.};
+ for (int j = 0; j < 3; ++j) {
+ factor += JpB[j] * trkParam.C(i + 2, j + 2);
+ }
+ val += JpA[i] * factor;
+ }
+ return val;
+ };
+ trackVector.SetCovariance(MomentumCovariance(0, 0), 9, iKfpTrk); // var(px)
+ trackVector.SetCovariance(MomentumCovariance(1, 0), 13, iKfpTrk); // cov(px, py)
+ trackVector.SetCovariance(MomentumCovariance(1, 1), 14, iKfpTrk); // var(py)
+ trackVector.SetCovariance(MomentumCovariance(2, 0), 18, iKfpTrk); // cov(px, pz)
+ trackVector.SetCovariance(MomentumCovariance(2, 1), 19, iKfpTrk); // cov(py, pz)
+ trackVector.SetCovariance(MomentumCovariance(2, 2), 20, iKfpTrk); // var(pz)
+
+ // Zero covariances (with z-coordinate)
+ trackVector.SetCovariance(0.f, 3, iKfpTrk); // cov(x,z)
+ trackVector.SetCovariance(0.f, 4, iKfpTrk); // cov(y,z)
+ trackVector.SetCovariance(0.f, 5, iKfpTrk); // var(z)
+ trackVector.SetCovariance(0.f, 8, iKfpTrk); // cov(z,px)
+ trackVector.SetCovariance(0.f, 12, iKfpTrk); // cov(z,py)
+ trackVector.SetCovariance(0.f, 17, iKfpTrk); // var(z,pz)
+
+ // ----- Other quantities
+ // Magnetic field (NOTE: zero fom mCBM)
+ // FIXME: Provide a proper initialization for full CBM
+ for (int iF = 0; iF < 10; ++iF) {
+ trackVector.SetFieldCoefficient(0.f, iF, iKfpTrk);
+ }
+
+ trackVector.SetId(iCaTrk, iKfpTrk);
+ trackVector.SetPDG(particleInfo.fPdg, iKfpTrk);
+ trackVector.SetQ(particleInfo.fCharge, iKfpTrk);
+ trackVector.SetNPixelHits(0, iKfpTrk);
+
+ // NOTE: 0 - primary tracks, -1 - secondary tracks. Here for now we assign ALL tracks as secondary
+ trackVector.SetPVIndex(-1, iKfpTrk);
}
diff --git a/algo/kfp/KfpV0Finder.h b/algo/kfp/KfpV0Finder.h
index ef8b1e22b7..0272b0ab89 100644
--- a/algo/kfp/KfpV0Finder.h
+++ b/algo/kfp/KfpV0Finder.h
@@ -24,6 +24,18 @@ namespace cbm::algo::kfp
/// \brief A V0-finding algorithm
class V0Finder {
public:
+ /// \struct ParticleInfo
+ /// \brief A structure to keep temporary PID information for tracks
+ struct ParticleInfo {
+ double fMass{std::numeric_limits<double>::quiet_NaN()}; //< Estimated mass of the particle [GeV/c2]
+ double fDca{std::numeric_limits<double>::quiet_NaN()}; //< DCA to the origin [cm]
+ double fBeta{std::numeric_limits<double>::quiet_NaN()}; //< Speed of the particle [c]
+ double fQp{std::numeric_limits<double>::quiet_NaN()}; //< q/p of the particle [GeV^-1]
+ int32_t fPdg{V0Finder::kUndefPdg}; //< PDG code
+ int32_t fCharge{0}; //< Charge of the particle
+ bool fbSelected{false}; //< The track was selected
+ };
+
/// \brief Default constructor
V0Finder() = default;
@@ -46,14 +58,30 @@ namespace cbm::algo::kfp
/// \param pdg A PDG code of the particle to be reconstructed
void AddDecayToReconstructionList(int pdg) { GetKFParticleFinder()->AddDecayToReconstructionList(pdg); }
+ /// \brief Gets monitor data
+ const V0FinderMonitorData_t& GetEventMonitor() const { return fEventMonitor; }
+
+ /// \brief Gets origin
+ const std::array<float, 3>& GetOrigin() const { return fOrigin; }
+
/// \brief Mutable access to the KfParticleFinder of the run topology reconstructor
KFParticleFinder* GetKFParticleFinder() { return fpTopoReconstructor->GetKFParticleFinder(); }
/// \brief Constant access to the KfParticleFinder of the run topology reconstructor
const KFParticleFinder* GetKFParticleFinder() const { return fpTopoReconstructor->GetKFParticleFinder(); }
- /// \brief Gets monitor data
- const V0FinderMonitorData_t& GetEventMonitor() const { return fEventMonitor; }
+ /// \brief Gets a vector of particle info
+ const std::vector<ParticleInfo> GetParticleInfo() const { return fvParticleInfo; }
+
+ /// \brief Gets selected t0
+ /// \note NaN, if Lambda-candidate was not found
+ double GetSelectedT0() const { return fSelectedT0; }
+
+ /// \brief Gets indices of selected tracks
+ const std::vector<uint32_t>& GetSelectedTrackIds() const { return fvSelectedTrackIds; }
+
+ /// \brief Gets found t0s
+ const std::vector<double>& GetT0s() const { return fvT0s; }
/// \brief Initializes the instance (called in the beginning of the run)
void Init();
@@ -77,7 +105,7 @@ namespace cbm::algo::kfp
/// \param y Y-coordinate of the origin [cm]
/// \param z Z-coordinate of the origin [cm]
// FIXME: for now origin is defined as the target center, later it can be changed
- void SetOrigin(double x, double y, double z) { fOrigin = {x, y, z}; }
+ void SetOrigin(double x, double y, double z) { fOrigin = {float(x), float(y), float(z)}; }
/// \brief Sets minimal pion DCA to primary vertex
/// \brief Sets pion velocity range
@@ -110,7 +138,6 @@ namespace cbm::algo::kfp
/// \param offset An offset [ns]
void SetTzeroOffset(double offset) { fTzeroOffset = offset; }
- //* Specific parameters for the KFParticleFinder
/// \brief Sets cut on the distance to the primary vertex from the decay vertex
/// \param cut Cut value [cm]
void SetLCut(float cut) { GetKFParticleFinder()->SetLCut(cut); }
@@ -128,17 +155,40 @@ namespace cbm::algo::kfp
void SetLdLCut2D(float cut) { GetKFParticleFinder()->SetLdLCut2D(cut); }
private:
+ /// \brief Assigns momentum based on the TOF measurement
+ /// \param[in] tofHits Tof hits container
+ /// \param[in] tofHitIds Tof hit indices, used by the track
+ /// \param[in] t0 A t0 value
+ /// \param[inout] pidInfo PID information for the track
+ /// \return true A physically reasonable momentum assigned
+ /// \return false Momentum was not assigned, because it was nonphysical
+ /* clang-format off */
+ bool AssignMomentum(const PartitionedVector<tof::Hit>& tofHits,
+ const std::vector<RecoResults::HitId_t>& tofHitIds,
+ double t0,
+ ParticleInfo& pidInfo);
+ /* clang-format on */
+
+ /// \brief Assigns PID info based on the estimated DCA
+ /// \param dca DCA of track to origin
+ /// \return (mass, charge, pid)
+ void AssignPid(ParticleInfo& info);
+
+ /// \brief Collects a vector of DCA
+ /// \param recoEvent Instance of a reconstructed event
+ void CollectDca(const RecoResults& recoEvent);
+
/// \brief Collects T0 values among the BMON hits
/// \param bmonHits A span of BMON hits
- /// \return A vector of T0-s
///
/// If multiple T0-s are found, the routine will run multiple times, until V0-candidates are found
- std::vector<double> CollectT0(gsl::span<const bmon::Hit> bmonHits) const;
+ void CollectT0(gsl::span<const bmon::Hit> bmonHits);
- /// \brief Collects a vector of DCA
- /// \param recoEvent Instance of a reconstructed event
- /// \return A vector of DCAs to origin
- std::vector<double> CollectDca(const RecoResults& recoEvent) const;
+ /// \brief Estimates speed of particle, using TOF measurement
+ /// \param tofHit A TOF hit
+ /// \param t0 An t0 value
+ /// \return Speed of particle [c]
+ double EstimateBeta(const tof::Hit& tofHit, double t0) const;
/// \brief Estimate DCA of a track to origin
/// \param fst first STS hit
@@ -146,29 +196,68 @@ namespace cbm::algo::kfp
/// \return dca [cm]
double EstimateDca(const sts::Hit& fst, const sts::Hit& snd) const;
- //* Physical constants
+ /// \brief Tries to find V0-candidates for a given t0
+ /// \param recoEvent Instance of a reconstructed event
+ /// \param t0 An estimated t0 value
+ /// \return true V0-candidates found
+ /// \return false V0-candidates not found
+ bool FindV0Candidates(const RecoResults& recoEvent, double t0);
+
+ /// \brief Initializes copies of track parameter vectors
+ /// \param tracks A container of tracks
+ void InitTrackParamVectors(const ca::Vector<ca::Track>& tracks);
+
+ /// \brief Makes a KF vertex
+ /// \param r coordinates of PV [cm]
+ static KFVertex MakeKfpPrimaryVertex(const std::array<float, 3>& r);
+
+ /// \brief Applies selection cut on the track
+ /// \param particleInfo Particle collected information
+ /// \return true Track is selected
+ /// \return false Track is rejected
+ bool SelectTrack(const ParticleInfo& particleInfo) const;
+
+ /// \brief Sets KFP track parameters
+ /// \param[inout] kfpTrkVector Reference to the KFP track vector
+ /// \param[in] iKfpTrk Index of the KFP track
+ /// \param[in] iCaTrk Index of the CA track
+ /// \param[in] trkParam Track parameters
+ /// \param[in] particleInfo Particle information
+ void SetKfpTrackParameters(KFPTrackVector& kfpTrkVector, uint32_t iKfpTrk, uint32_t iCaTrk,
+ const ca::Track::TrackParam_t& trkParam, const ParticleInfo& particleInfo) const;
+
+ //* Framework and physical constants
static constexpr double kPionMass{0.13957039}; ///< Pion mass [GeV/c2]
static constexpr double kProtonMass{0.938272088}; ///< Proton mass [GeV/c2]
static constexpr double kSpeedOfLight{29.9792458}; ///< Speed of light [cm/ns]
+ static constexpr int32_t kUndefPdg{-2}; ///< PDG for tracks, which PID cannot be inferred
+ static constexpr bool kUseAverageSpeed{false}; ///< If an average speed of tof hits is used
-
- V0FinderMonitorData_t fEventMonitor; ///< Monitor data instance (per event)
+ V0FinderMonitorData_t fEventMonitor; ///< Main monitor data instance
//* Different run-time cuts and flags (TODO: define in a config)
double fTzeroOffset{0.}; ///< Offset for T0
double fMinPionDca{1.5}; ///< Minimum DCA to PV for pions
double fMinProtonDca{0.5}; ///< Minimum DCA to PV for protons
double fQpAssignedUncertainty{0.1}; ///< Assigned relative uncertainty for q/p estimation
- int fPrimaryAssignedPdg{321}; ///< Assigned PDG hypothesis for primary particles
double fMinBetaProton{0.}; ///< Minimal proton velocity (beta) [c]
double fMaxBetaProton{1.}; ///< Maximal proton velocity (beta) [c]
double fMinBetaPion{0.}; ///< Minimal proton velocity (beta) [c]
double fMaxBetaPion{1.}; ///< Maximal proton velocity (beta) [c]
-
+ int fPrimaryAssignedPdg{321}; ///< Assigned PDG hypothesis for primary particles
//* Run-time variables, provided by framework
- int fBmonPartitionIndex{-1}; ///< Index of selected partition in BMON hit vector
- std::array<double, 3> fOrigin{0., 0., 0.}; ///< Coordinates of origin [cm]
+ int fBmonPartitionIndex{-1}; ///< Index of selected partition in BMON hit vector
+
+ //* Temporary data arrays (NOTE: keep them here for QA)
+ std::array<float, 3> fOrigin{0.f, 0.f, 0.f}; ///< Coordinates of origin [cm]
+ std::vector<double> fvT0s; ///< Found t0s [ns] (in event)
+ std::vector<ParticleInfo> fvParticleInfo; ///< PID info of tracks (in event)
+ std::vector<uint32_t> fvSelectedTrackIds; ///< IDs of selected tracks (in event)
+ double fSelectedT0{std::numeric_limits<double>::quiet_NaN()}; ///< A t0 value selected by the lambda-finder
+
+ /// \brief A copy of track parameters (first, last)
+ std::vector<std::pair<ca::Track::TrackParam_t, ca::Track::TrackParam_t>> fvTrackParam;
//* Auxilary variables
diff --git a/algo/kfp/KfpV0FinderConfig.cxx b/algo/kfp/KfpV0FinderConfig.cxx
index 54e08ec323..d74e4d8b2c 100644
--- a/algo/kfp/KfpV0FinderConfig.cxx
+++ b/algo/kfp/KfpV0FinderConfig.cxx
@@ -9,6 +9,7 @@
#include "KfpV0FinderConfig.h"
+#include <iomanip>
#include <sstream>
using cbm::algo::kfp::Cuts;
@@ -58,12 +59,17 @@ std::string Cuts::ToString() const
//
std::string V0FinderConfig::ToString() const
{
+ using std::dec;
+ using std::hex;
+ using std::setfill;
+ using std::setw;
std::stringstream msg;
msg << "\n-------------- KFP V0-finder Configuration: ------------------------------------------------";
msg << "\nGOAL V0 (PDG): " << reconstructPdg;
msg << '\n' << cuts.ToString();
msg << "\nOTHER PARAMETERS:";
msg << "\n\tt0-offset: " << tZeroOffset << " [ns]";
+ msg << "\n\tReference BMON: 0x" << hex << setw(8) << setfill('0') << bmonAddress << dec;
msg << "\n\tq/p relative uncertainty: " << qpAssignedUncertainty;
msg << "\n\tassigned PDG of primary: " << primaryAssignedPdg;
msg << "\n--------------------------------------------------------------------------------------------";
diff --git a/algo/kfp/KfpV0FinderMonitor.h b/algo/kfp/KfpV0FinderMonitor.h
index 1f63bb246a..5bb4311b66 100644
--- a/algo/kfp/KfpV0FinderMonitor.h
+++ b/algo/kfp/KfpV0FinderMonitor.h
@@ -26,6 +26,8 @@ namespace cbm::algo::kfp
TracksWoPid, ///< Tracks, which has undefined PID
TracksWoMomentum, ///< Tracks, which has no momentum
TracksWUnphysicalBeta, ///< Tracks with beta > 1
+ PionsDca, ///< Number of raw pion-candidates
+ ProtonsDca, ///< Number of raw proton-candidates
Pions, ///< Number of pion-candidates
Protons, ///< Number of proton-candidates
EventsTotal, ///< Total number of events
@@ -38,11 +40,20 @@ namespace cbm::algo::kfp
/// \enum ETimer
/// \brief Timer keys for the V0FinderMonitor
+ /* clang-format off */
enum class ETimer
{
- Event, ///< Timer to process a single event
+ ProcessEvent, ///< Processing of a single event
+ CollectT0, ///< Collecting T0s
+ CollectDca, ///< Estimating DCAs
+ FindV0Candidates, ///< V0-finder procedure for a given t0
+ PrepareContainers, ///< Prepare data containers
+ PreselectTracks, ///< Track preselection
+ InitKfp, ///< Init KFParticleFinder inside the event
+ ExecKfp, ///< Run KFParticleFinder inside the event
END
};
+ /* clang-format on */
/// \brief Specification of ca::MonitorData for the V0Finder
using V0FinderMonitorData_t = ca::MonitorData<ECounter, ETimer>;
@@ -63,6 +74,8 @@ namespace cbm::algo::kfp
SetCounterName(ECounter::TracksWoPid, "tracks w/o PID");
SetCounterName(ECounter::TracksWoMomentum, "tracks w/o momentum");
SetCounterName(ECounter::TracksWUnphysicalBeta, "tracks w/ beta > 1");
+ SetCounterName(ECounter::PionsDca, "raw pion candidates");
+ SetCounterName(ECounter::ProtonsDca, "raw proton candidates");
SetCounterName(ECounter::Pions, "pion candidates");
SetCounterName(ECounter::Protons, "proton candidates");
SetCounterName(ECounter::EventsTotal, "all events");
@@ -71,7 +84,14 @@ namespace cbm::algo::kfp
SetCounterName(ECounter::KfpEventsLambdaCand, "events w/ lambda candidates");
SetCounterName(ECounter::KfpLambdaCandidates, "lambda candidates");
- SetTimerName(ETimer::Event, "event processing");
+ SetTimerName(ETimer::ProcessEvent, "event processing");
+ SetTimerName(ETimer::CollectT0, "t0 container preparation");
+ SetTimerName(ETimer::CollectDca, "DCA container preparation");
+ SetTimerName(ETimer::FindV0Candidates, "V0-candidates finding");
+ SetTimerName(ETimer::PrepareContainers, "Container initialization");
+ SetTimerName(ETimer::PreselectTracks, "Track preselection");
+ SetTimerName(ETimer::InitKfp, "KFParticleFinder initialization");
+ SetTimerName(ETimer::ExecKfp, "KFParticleFinder execution");
}
private:
--
GitLab