diff --git a/algo/CMakeLists.txt b/algo/CMakeLists.txt
index a2a052fa4026f62246d4be605c6a50c2217bf65a..a70b0872c98516d229715c67c7f86b14d8102a83 100644
--- a/algo/CMakeLists.txt
+++ b/algo/CMakeLists.txt
@@ -52,10 +52,12 @@ if (CBM_ONLINE_STANDALONE)
add_subdirectory(../external external)
endif()
+add_subdirectory(kfp)
add_subdirectory(log)
add_subdirectory(data)
add_subdirectory(kf)
add_subdirectory(ca)
+#add_subdirectory(kfp) # For KFParticleOnline
# exclude unittests from being build inside the container
if (NOT CBM_ONLINE_STANDALONE)
@@ -160,7 +162,10 @@ set(SRCS
ca/TrackingSetup.cxx
ca/TrackingChain.cxx
ca/qa/CaQa.cxx
+ kfp/KfpV0Finder.cxx
+ kfp/KfpV0FinderChain.cxx
kfp/KfpV0FinderConfig.cxx
+ kfp/KfpV0FinderQa.cxx
)
set(BUILD_INFO_CXX ${CMAKE_CURRENT_BINARY_DIR}/base/BuildInfo.cxx)
@@ -223,6 +228,7 @@ target_link_libraries(Algo
external::fles_monitoring
cppzmq
poolstl
+ PRIVATE CbmKFParticleOnlineInterface
)
target_compile_definitions(Algo PUBLIC NO_ROOT)
xpu_attach(Algo ${DEVICE_SRCS})
@@ -296,6 +302,7 @@ if (NOT CBM_ONLINE_STANDALONE)
external::fles_monitoring
cppzmq
poolstl
+ PRIVATE CbmKFParticleOnlineInterface
)
xpu_attach(AlgoOffline ${DEVICE_SRCS})
diff --git a/algo/base/Definitions.h b/algo/base/Definitions.h
index 6c6d77ba8715786fada47c3483207254bb269b49..b1ff4c9ccc8e804a33a67fe0ab40d81d74dd616e 100644
--- a/algo/base/Definitions.h
+++ b/algo/base/Definitions.h
@@ -78,6 +78,7 @@ namespace cbm::algo
RecoTof,
RecoFsd,
Tracking,
+ V0Finder,
};
} // namespace cbm::algo
@@ -148,7 +149,8 @@ CBM_ENUM_DICT(cbm::algo::QaStep,
{"RecoMuch", cbm::algo::QaStep::RecoMuch},
{"RecoTof", cbm::algo::QaStep::RecoTof},
{"RecoFsd", cbm::algo::QaStep::RecoFsd},
- {"Tracking", cbm::algo::QaStep::Tracking}
+ {"Tracking", cbm::algo::QaStep::Tracking},
+ {"V0Finder", cbm::algo::QaStep::V0Finder}
);
#endif
diff --git a/algo/base/DigiData.cxx b/algo/base/DigiData.cxx
index 179c8ca31508e322b1bd1222b70bc417a9aafd7b..79fd5a4836022d2caae819e0715f3fcf0a1ecfe7 100644
--- a/algo/base/DigiData.cxx
+++ b/algo/base/DigiData.cxx
@@ -119,14 +119,16 @@ DigiEvent::DigiEvent(const CbmDigiEvent& storable)
: DigiData(storable.fData)
, fNumber(storable.fNumber)
, fTime(storable.fTime)
+ , fSelectionTriggers(storable.fSelectionTriggers)
{
}
CbmDigiEvent DigiEvent::ToStorable() const
{
return CbmDigiEvent{
- .fData = DigiData::ToStorable(),
- .fNumber = fNumber,
- .fTime = fTime,
+ .fData = DigiData::ToStorable(),
+ .fNumber = fNumber,
+ .fTime = fTime,
+ .fSelectionTriggers = fSelectionTriggers,
};
}
diff --git a/algo/base/DigiData.h b/algo/base/DigiData.h
index c19b685d996ea4fadc68b7c50e170f356e3c3a44..09b82d933478c57b65087b93e1d0a799120d287b 100644
--- a/algo/base/DigiData.h
+++ b/algo/base/DigiData.h
@@ -7,6 +7,7 @@
#include "CbmBmonDigi.h"
#include "CbmDigiData.h"
#include "CbmDigiEvent.h"
+#include "CbmEventTriggers.h"
#include "CbmFsdDigi.h"
#include "CbmMuchDigi.h"
#include "CbmPsdDigi.h"
@@ -80,6 +81,7 @@ namespace cbm::algo
// FIXME: Event number not set yet!
uint64_t fNumber = -1; ///< Event identifier
double fTime = 0; ///< Event trigger time [ns]
+ CbmEventTriggers fSelectionTriggers;
static std::vector<DigiEvent> FromCbmDigiEvents(const std::vector<CbmDigiEvent>& events);
static std::vector<CbmDigiEvent> ToCbmDigiEvents(const std::vector<DigiEvent>& events);
diff --git a/algo/base/Options.cxx b/algo/base/Options.cxx
index ea9ecace39cfe29f8fff9cec2161871e9f08c569..abce4f27d969904068407a2022446d9cc36946c4 100644
--- a/algo/base/Options.cxx
+++ b/algo/base/Options.cxx
@@ -93,6 +93,7 @@ Options::Options(int argc, char** argv)
("compress-archive", po::bool_switch(&fCompressArchive)->default_value(false), "Enable compression for output archives")
("steps", po::value(&fRecoSteps)->multitoken()->default_value({Step::Unpack, Step::DigiTrigger, Step::LocalReco, Step::Tracking})->value_name("<steps>"),
"space separated list of reconstruction steps (unpack, digitrigger, localreco, ...)")
+ ("event-reco", po::bool_switch(&fReconstructDigiEvents)->default_value(false), "runs digi event reconstruction (local reco, tracking, trigger)")
("systems,s", po::value(&fDetectors)->multitoken()->default_value({Subsystem::STS, Subsystem::TOF, Subsystem::BMON, Subsystem::MUCH, Subsystem::RICH, Subsystem::TRD, Subsystem::TRD2D})->value_name("<detectors>"),
"space separated list of detectors to process (sts, mvd, ...)")
("child-id,c", po::value(&fChildId)->default_value("00")->value_name("<id>"), "online process id on node")
diff --git a/algo/base/Options.h b/algo/base/Options.h
index 5a9d232fcac8b006b54cfc0ab3f4cdbe5b62a535..f1b508a0699f5839a0b1f5b0e764f0d6a29f8e34 100644
--- a/algo/base/Options.h
+++ b/algo/base/Options.h
@@ -68,6 +68,8 @@ namespace cbm::algo
bool Has(QaStep qastep) const;
+ bool ReconstructDigiEvents() const { return fReconstructDigiEvents; }
+
private: // members
std::string fParamsDir; // TODO: can we make this a std::path?
std::string fInputLocator;
@@ -97,6 +99,7 @@ namespace cbm::algo
uint64_t fRunId = 2391;
uint64_t fRunStartTime = 0;
bool fCollectAuxData = false;
+ bool fReconstructDigiEvents = false;
};
} // namespace cbm::algo
diff --git a/algo/ca/TrackingChain.cxx b/algo/ca/TrackingChain.cxx
index 1d481261a1af5143deb8128a969175122fcf94eb..374c23bb1e053d132dda601a4049746b910c90e7 100644
--- a/algo/ca/TrackingChain.cxx
+++ b/algo/ca/TrackingChain.cxx
@@ -40,8 +40,10 @@ using cbm::algo::ca::constants::clrs::GNb; // grin bald text
// ---------------------------------------------------------------------------------------------------------------------
//
-TrackingChain::TrackingChain(const std::unique_ptr<cbm::algo::qa::Manager>& qaManager, std::string_view name)
+TrackingChain::TrackingChain(ECbmRecoMode recoMode, const std::unique_ptr<cbm::algo::qa::Manager>& qaManager,
+ std::string_view name)
: fQa(Qa(qaManager, name))
+ , fRecoMode(recoMode)
{
}
@@ -116,8 +118,13 @@ void TrackingChain::Init()
// ------ Initialize CA framework
fCaMonitor.Reset();
- fCaFramework.SetNofThreads(Opts().NumOMPThreads() == std::nullopt ? openmp::GetMaxThreads()
- : *(Opts().NumOMPThreads()));
+ if (ECbmRecoMode::Timeslice == fRecoMode) {
+ fCaFramework.SetNofThreads(Opts().NumOMPThreads() == std::nullopt ? openmp::GetMaxThreads()
+ : *(Opts().NumOMPThreads()));
+ }
+ else {
+ fCaFramework.SetNofThreads(1);
+ }
fCaFramework.ReceiveParameters(std::move(parameters));
fCaFramework.Init(ca::TrackingMode::kMcbm);
@@ -126,6 +133,7 @@ void TrackingChain::Init()
fQa.RegisterParameters(&fCaFramework.GetParameters());
fQa.Init();
}
+ L_(info) << "TRACKING QA: " << fQa.IsActive();
}
// ---------------------------------------------------------------------------------------------------------------------
@@ -174,18 +182,26 @@ void TrackingChain::PrepareInput(Input_t recoResults)
faHitExternalIndices.clear();
faHitExternalIndices.reserve(nHitsTot);
if (fbDetUsed[EDetectorID::kSts]) {
+ fCaMonitorData.StartTimer(ca::ETimer::PrepareStsHits);
ReadHits<EDetectorID::kSts>(recoResults.stsHits);
+ fCaMonitorData.StopTimer(ca::ETimer::PrepareStsHits);
}
if (fbDetUsed[EDetectorID::kTrd]) {
+ fCaMonitorData.StartTimer(ca::ETimer::PrepareTrdHits);
ReadHits<EDetectorID::kTrd>(recoResults.trdHits);
+ fCaMonitorData.StopTimer(ca::ETimer::PrepareTrdHits);
}
if (fbDetUsed[EDetectorID::kTof]) {
+ fCaMonitorData.StartTimer(ca::ETimer::PrepareTofHits);
ReadHits<EDetectorID::kTof>(recoResults.tofHits);
+ fCaMonitorData.StopTimer(ca::ETimer::PrepareTofHits);
}
faHitExternalIndices.shrink_to_fit();
fCaDataManager.SetNhitKeys(fNofHitKeys);
L_(debug) << "Tracking chain: " << fCaDataManager.GetNofHits() << " hits will be passed to the ca::Framework";
+ fCaMonitorData.StartTimer(ca::ETimer::InputDataTransmission);
fCaFramework.ReceiveInputData(fCaDataManager.TakeInputData());
+ fCaMonitorData.StopTimer(ca::ETimer::InputDataTransmission);
}
// ---------------------------------------------------------------------------------------------------------------------
@@ -210,6 +226,7 @@ TrackingChain::Output_t TrackingChain::PrepareOutput()
int iHitInternal = fCaFramework.GetInputData().GetHit(fCaFramework.fRecoHits[trackFirstHit + iHit]).Id();
const auto [detID, iPartition, iPartHit] = faHitExternalIndices[iHitInternal];
switch (detID) {
+ // FIXME: store a global hit index instead of (partition, hit)
case ca::EDetectorID::kSts: output.stsHitIndices[iTrk].push_back(std::make_pair(iPartition, iPartHit)); break;
case ca::EDetectorID::kTof: output.tofHitIndices[iTrk].push_back(std::make_pair(iPartition, iPartHit)); break;
case ca::EDetectorID::kTrd: output.trdHitIndices[iTrk].push_back(std::make_pair(iPartition, iPartHit)); break;
@@ -223,12 +240,13 @@ TrackingChain::Output_t TrackingChain::PrepareOutput()
trackFirstHit += nHits;
}
- L_(info) << "TrackingChain: Timeslice contains " << fCaMonitorData.GetCounterValue(ca::ECounter::RecoTrack)
- << " tracks, with " << fCaMonitorData.GetCounterValue(ca::ECounter::RecoStsHit) << " sts hits, "
- << fCaMonitorData.GetCounterValue(ca::ECounter::RecoTofHit) << " tof hits, "
- << fCaMonitorData.GetCounterValue(ca::ECounter::RecoTrdHit) << " trd hits"
- << "; the FindTracks routine ran " << fCaMonitorData.GetTimer(ca::ETimer::FindTracks).GetTotal() << " s";
-
+ if (ECbmRecoMode::Timeslice == fRecoMode) {
+ L_(info) << "TrackingChain: Timeslice contains " << fCaMonitorData.GetCounterValue(ca::ECounter::RecoTrack)
+ << " tracks, with " << fCaMonitorData.GetCounterValue(ca::ECounter::RecoStsHit) << " sts hits, "
+ << fCaMonitorData.GetCounterValue(ca::ECounter::RecoTofHit) << " tof hits, "
+ << fCaMonitorData.GetCounterValue(ca::ECounter::RecoTrdHit) << " trd hits"
+ << "; the FindTracks routine ran " << fCaMonitorData.GetTimer(ca::ETimer::FindTracks).GetTotal() << " s";
+ }
// QA
if (fQa.IsActive()) {
@@ -327,6 +345,7 @@ void TrackingChain::ReadHits(PartitionedSpan<const ca::HitTypes_t::at<DetID>> hi
//int iHitExt = iOffset + iPartHit;
// ---- Fill ca::Hit
+ fCaMonitorData.StartTimer(ca::ETimer::CaHitCreation);
ca::Hit caHit;
if constexpr (IsSts) {
caHit.SetFrontKey(firstHitKey + hit.fFrontClusterId);
@@ -337,6 +356,15 @@ void TrackingChain::ReadHits(PartitionedSpan<const ca::HitTypes_t::at<DetID>> hi
caHit.SetFrontKey(firstHitKey + iPartHit);
caHit.SetBackKey(caHit.FrontKey());
}
+
+ if constexpr (IsTof) {
+ // Cut the BMon hits (if any)
+ if (hit.Z() < 0.1) {
+ // FIXME: Provide BMon addresses explicitly in all the parameter files
+ continue;
+ }
+ }
+
caHit.SetX(hit.X());
caHit.SetY(hit.Y());
caHit.SetZ(hit.Z());
@@ -366,11 +394,16 @@ void TrackingChain::ReadHits(PartitionedSpan<const ca::HitTypes_t::at<DetID>> hi
caHit.SetStation(iStActive);
caHit.SetId(fCaDataManager.GetNofHits());
if (caHit.Check()) {
- if ((caHit.T() < lastTime - 1000.) && (dataStream < 100000)) {
- dataStream++;
+ if (ECbmRecoMode::Timeslice == fRecoMode) {
+ if ((caHit.T() < lastTime - 1000.) && (dataStream < 100000)) {
+ dataStream++;
+ }
+ lastTime = caHit.T();
+ fCaDataManager.PushBackHit(caHit, dataStreamDet | dataStream);
+ }
+ else {
+ fCaDataManager.PushBackHit(caHit); // A single data stream in event-by-event mode
}
- lastTime = caHit.T();
- fCaDataManager.PushBackHit(caHit, dataStreamDet | dataStream);
faHitExternalIndices.push_back(std::make_tuple(DetID, iPartition, iPartHit));
if (fNofHitKeys <= caHit.FrontKey()) {
fNofHitKeys = caHit.FrontKey() + 1;
@@ -396,6 +429,7 @@ void TrackingChain::ReadHits(PartitionedSpan<const ca::HitTypes_t::at<DetID>> hi
fCaMonitorData.IncrementCounter(ca::ECounter::UndefinedTofHit);
}
}
+ fCaMonitorData.StopTimer(ca::ETimer::CaHitCreation);
//prevTime = caHit.T();
// ---- Update number of hit keys
} // iPartHit
diff --git a/algo/ca/TrackingChain.h b/algo/ca/TrackingChain.h
index 3f949d3605e87e1f58b5404886a63ffd44c636d1..39ba48e26fca9cea0cfa295408eb14345a856918 100644
--- a/algo/ca/TrackingChain.h
+++ b/algo/ca/TrackingChain.h
@@ -41,9 +41,11 @@ namespace cbm::algo
class TrackingChain : public SubChain {
public:
/// \brief Constructor from parameters
+ /// \param recoMode Reconstruction mode
/// \param pManager a QA-manager
/// \param name A name of the task (histograms directory)
- TrackingChain(const std::unique_ptr<cbm::algo::qa::Manager>& qaManager = nullptr, std::string_view name = "");
+ TrackingChain(ECbmRecoMode recoMode, const std::unique_ptr<cbm::algo::qa::Manager>& qaManager = nullptr,
+ std::string_view name = "");
/// \struct Input_t
/// \brief Input to the TrackingChain
@@ -92,6 +94,7 @@ namespace cbm::algo
/// \brief Provides action in the end of the run
void Finalize();
+
private:
// *********************
// ** Utility functions
@@ -132,6 +135,8 @@ namespace cbm::algo
/// \note Indexing: [globalHitID], value: (DetID, partitionID, hitID)
ca::Vector<std::tuple<ca::EDetectorID, uint32_t, uint32_t>> faHitExternalIndices{"faHitExternalIndices"};
+ ECbmRecoMode fRecoMode{ECbmRecoMode::Undefined}; ///< Reconstruction mode
+
static constexpr bool kDEBUG = false; ///< Debug mode
};
diff --git a/algo/ca/core/data/CaDataManager.cxx b/algo/ca/core/data/CaDataManager.cxx
index 7d237151d6b690ccdaa9cdc4439aba2ee0fb913b..4e74a2524c9aff2641deb285cce8262be38d310f 100644
--- a/algo/ca/core/data/CaDataManager.cxx
+++ b/algo/ca/core/data/CaDataManager.cxx
@@ -17,23 +17,6 @@
using cbm::algo::ca::DataManager;
using cbm::algo::ca::InputData;
-// ---------------------------------------------------------------------------------------------------------------------
-//
-bool DataManager::SendInputData(InputData& destination)
-{
- // Set boundary hit indexes
- InitData();
-
- // Check data before input
- if (CheckInputData<constants::control::InputDataQaLevel>()) {
- destination = std::move(fInputData);
- assert(this->GetNofHits() == 0);
- return true;
- }
- LOG(error) << "L1: Attempt to set up inconsistent input data";
- return false;
-}
-
// ---------------------------------------------------------------------------------------------------------------------
//
InputData&& DataManager::TakeInputData()
@@ -93,16 +76,21 @@ void DataManager::InitData()
{
// set the end pointers to data streams
// TODO: SZh 14.08.2023: Move the max allowed number of streams to the constants.h
- if (fInputData.fvStreamStartIndices.size() > 3000) {
- LOG(warning) << "ca::DataManager: unexpected order of input data: too many data streams!!! ";
- fInputData.fvStreamStartIndices.shrink(3000);
- }
+
int nStreams = fInputData.fvStreamStartIndices.size();
- fInputData.fvStreamStopIndices.reset(nStreams);
- for (int i = 0; i < nStreams - 1; i++) {
- fInputData.fvStreamStopIndices[i] = fInputData.fvStreamStartIndices[i + 1];
+ if (!nStreams) { // No data streams provided
+ fInputData.fvStreamStartIndices.push_back(0);
+ fInputData.fvStreamStopIndices.push_back(fInputData.fvHits.size());
}
- if (nStreams > 0) {
+ else {
+ if (nStreams > 3000) {
+ LOG(warning) << "ca::DataManager: unexpected order of input data: too many data streams!!! ";
+ fInputData.fvStreamStartIndices.shrink(3000);
+ }
+ fInputData.fvStreamStopIndices.reset(nStreams);
+ for (int i = 0; i < nStreams - 1; i++) {
+ fInputData.fvStreamStopIndices[i] = fInputData.fvStreamStartIndices[i + 1];
+ }
fInputData.fvStreamStopIndices[nStreams - 1] = fInputData.fvHits.size();
}
}
@@ -112,10 +100,12 @@ void DataManager::InitData()
void DataManager::WriteInputData(const std::string& fileName) const
{
// Check current data object for consistency
- if (!CheckInputData<constants::control::InputDataQaLevel>()) {
- LOG(error) << "ca::DataManager: input data writer: attempt to write invalid input data object to file \""
- << fileName << "\"";
- return;
+ if constexpr (constants::control::InputDataQaLevel > 0) {
+ if (!CheckInputData<constants::control::InputDataQaLevel>()) {
+ LOG(error) << "ca::DataManager: input data writer: attempt to write invalid input data object to file \""
+ << fileName << "\"";
+ return;
+ }
}
// Open output binary file
diff --git a/algo/ca/core/data/CaDataManager.h b/algo/ca/core/data/CaDataManager.h
index d210e52b93035516022b1b8411d1a1761b44dc5e..b8dac8630607d978e9152ce3ac392c7f765bf5c1 100644
--- a/algo/ca/core/data/CaDataManager.h
+++ b/algo/ca/core/data/CaDataManager.h
@@ -54,8 +54,9 @@ namespace cbm::algo::ca
/// \param nHits Number of hits to reserve
void ResetInputData(HitIndex_t nHits = 0) noexcept;
- /// \brief Pushes back a hit
- /// \param hit An ca::Hit object
+ /// \brief Pushes back a hit (with a data stream info)
+ /// \param hit An ca::Hit object
+ /// \param streamId Index of a data stream
void PushBackHit(const Hit& hit, int64_t streamId)
{
if (fInputData.fvStreamStartIndices.size() == 0 || fLastStreamId != streamId) { // new data stream
@@ -67,14 +68,14 @@ namespace cbm::algo::ca
fInputData.fvHits.push_back(hit);
}
+ /// \brief Pushes back a hit
+ /// \param hit An ca::Hit object
+ void PushBackHit(const Hit& hit) { fInputData.fvHits.push_back(hit); }
+
/// \brief Sets the number of hit keys
/// \param nKeys Number of hit keys
void SetNhitKeys(int nKeys) { fInputData.fNhitKeys = nKeys; }
- /// \brief Sends (moves) input data to an object (alternative method of data sending)
- /// \param destination Destination object of input data
- bool SendInputData(InputData& destination);
-
/// \brief Takes (moves) the instance of the input data object
InputData&& TakeInputData();
diff --git a/algo/ca/core/data/CaTrack.h b/algo/ca/core/data/CaTrack.h
index 29b7f3859db3513eaa525e474e5d2b3aa6f61d06..cf466361076dd8fcadcd22d80e53bf69832a8be8 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/CaEnumArray.h b/algo/ca/core/utils/CaEnumArray.h
index 5acd385893cffa7941e92182c268f44adbe0aed8..84e90a0eda1557717935c876b51dd6514b69235a 100644
--- a/algo/ca/core/utils/CaEnumArray.h
+++ b/algo/ca/core/utils/CaEnumArray.h
@@ -38,7 +38,7 @@ namespace cbm::algo::ca
public:
/// \brief Mutable access operator, indexed by enum entry
- T& operator[](const E& entry)
+ T& operator[](const E entry)
{
return std::array<T, static_cast<std::size_t>(E::END)>::operator[](static_cast<U>(entry));
}
diff --git a/algo/ca/core/utils/CaMonitor.h b/algo/ca/core/utils/CaMonitor.h
index 389b196fef4b51be0b5c8b054c0967cb6216fe7b..6d26469cfa6ceec81d71a9742786360776629f24 100644
--- a/algo/ca/core/utils/CaMonitor.h
+++ b/algo/ca/core/utils/CaMonitor.h
@@ -171,12 +171,13 @@ namespace cbm::algo::ca
using std::setfill;
using std::setw;
std::stringstream msg;
- constexpr size_t widthKey = 30;
constexpr size_t width = 24;
+ auto Cmp = [](const std::string& l, const std::string& r) { return l.size() < r.size(); };
msg << "\n===== Monitor: " << fsName << "\n";
if constexpr (!std::is_same_v<ETimerKey, EDummy>) {
+ size_t widthKeyTimer{std::max_element(faTimerNames.begin(), faTimerNames.end(), Cmp)->size() + 1};
msg << "\n----- Timers:\n";
- msg << setw(widthKey) << left << "Key" << ' ';
+ msg << setw(widthKeyTimer) << left << "Key" << ' ';
msg << setw(width) << left << "N Calls" << ' ';
msg << setw(width) << left << "Average [s]" << ' ';
msg << setw(width) << left << "Min [s]" << ' ';
@@ -185,7 +186,7 @@ namespace cbm::algo::ca
msg << right;
for (int iKey = 0; iKey < fMonitorData.GetNofTimers(); ++iKey) {
const Timer& timer = fMonitorData.GetTimer(static_cast<ETimerKey>(iKey));
- msg << setw(widthKey) << faTimerNames[iKey] << ' ';
+ msg << setw(widthKeyTimer) << faTimerNames[iKey] << ' ';
msg << setw(width) << timer.GetNofCalls() << ' ';
msg << setw(width) << timer.GetAverage() << ' ';
msg << setw(width) << timer.GetMin() << ' ';
@@ -195,7 +196,8 @@ namespace cbm::algo::ca
}
msg << "\n----- Counters:\n";
- msg << setw(widthKey) << left << "Key" << ' ';
+ size_t widthKeyCounter{std::max_element(faCounterNames.begin(), faCounterNames.end(), Cmp)->size() + 1};
+ msg << setw(widthKeyCounter) << left << "Key" << ' ';
msg << setw(width) << left << "Total" << ' ';
for (auto key : fvCounterRatioKeys) {
msg << setw(width) << left << std::string("per ") + faCounterNames[key] << ' ';
@@ -203,7 +205,7 @@ namespace cbm::algo::ca
msg << '\n';
for (int iKey = 0; iKey < fMonitorData.GetNofCounters(); ++iKey) {
auto counterValue = fMonitorData.GetCounterValue(static_cast<ECounterKey>(iKey));
- msg << setw(widthKey) << left << faCounterNames[iKey] << ' ';
+ msg << setw(widthKeyCounter) << left << faCounterNames[iKey] << ' ';
msg << setw(width) << right << counterValue << ' ';
for (auto keyDen : fvCounterRatioKeys) {
auto ratio = static_cast<double>(counterValue) / fMonitorData.GetCounterValue(keyDen);
diff --git a/algo/ca/core/utils/CaMonitorData.h b/algo/ca/core/utils/CaMonitorData.h
index 09b3d0e1174cc385357d17eaff108538cc7d80ec..5d93816c1fa0d8998957e98ff6f1b7d0651a9ff8 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/ca/core/utils/CaTrackingMonitor.h b/algo/ca/core/utils/CaTrackingMonitor.h
index f91de1d75ac356d14bbe01c330800d188b33788e..ffcee0c35371b8b30b40d84b82c347b497c58cd9 100644
--- a/algo/ca/core/utils/CaTrackingMonitor.h
+++ b/algo/ca/core/utils/CaTrackingMonitor.h
@@ -52,6 +52,11 @@ namespace cbm::algo::ca
{
TrackingChain,
PrepareInputData,
+ PrepareStsHits,
+ PrepareTrdHits,
+ PrepareTofHits,
+ InputDataTransmission,
+ CaHitCreation,
Tracking,
PrepareTimeslice,
TrackingThread,
@@ -103,6 +108,11 @@ namespace cbm::algo::ca
SetTimerName(ETimer::TrackingChain, "tracking chain");
SetTimerName(ETimer::PrepareInputData, "input data preparation");
+ SetTimerName(ETimer::PrepareStsHits, "STS hits preparation");
+ SetTimerName(ETimer::PrepareTrdHits, "TRD hits preparation");
+ SetTimerName(ETimer::PrepareTofHits, "TOF hits preparation");
+ SetTimerName(ETimer::InputDataTransmission, "input data transmission");
+ SetTimerName(ETimer::CaHitCreation, "CA hit creation");
SetTimerName(ETimer::Tracking, "algorithm execution");
SetTimerName(ETimer::PrepareTimeslice, "timeslice preparation");
SetTimerName(ETimer::TrackingThread, "tracking on one thread");
diff --git a/algo/ca/qa/CaQa.cxx b/algo/ca/qa/CaQa.cxx
index 05815198dd6e6a228f9db995496470a830d18b2c..8d081862c78d399aca6ec01d0dda38e432517185 100644
--- a/algo/ca/qa/CaQa.cxx
+++ b/algo/ca/qa/CaQa.cxx
@@ -187,22 +187,22 @@ void Qa::Init()
continue;
}
{
- auto sName = format("track_{}_theta", vsPointName[i]);
+ auto sName = format("{}_track_{}_theta", GetTaskName(), vsPointName[i]);
auto sTitl = format("#theta at {} hit; #theta", vsPointName[i]);
fvphTrkTheta[i] = MakeObj<H1D>(sName, sTitl, 62, 0., 90.);
}
{
- auto sName = format("track_{}_phi", vsPointName[i]);
+ auto sName = format("{}_track_{}_phi", GetTaskName(), vsPointName[i]);
auto sTitl = format("#phi at {} hit; #phi", vsPointName[i]);
fvphTrkPhi[i] = MakeObj<H1D>(sName, sTitl, 62, -180., 180.);
}
{
- auto sName = format("track_{}_thata_phi", vsPointName[i]);
+ auto sName = format("{}_track_{}_thata_phi", GetTaskName(), vsPointName[i]);
auto sTitl = format("#theta vs #phi at {} hit; #phi; #theta", vsPointName[i]);
fvphTrkPhiTheta[i] = MakeObj<H2D>(sName, sTitl, 62, -180., 180., 62, 0., 90.);
}
{
- auto sName = format("track_{}_chi2_ndf", vsPointName[i]);
+ auto sName = format("{}_track_{}_chi2_ndf", GetTaskName(), vsPointName[i]);
auto sTitl = format("#chi^{{2}}/NDF at {} hit; #chi^{{2}}/NDF", vsPointName[i]);
fvphTrkChi2Ndf[i] = MakeObj<H1D>(sName, sTitl, 100, 0., 20.);
}
@@ -212,7 +212,7 @@ void Qa::Init()
double xMin = -0.5;
double xMax = double(knStaMax) - 0.5;
{
- auto sName = "track_fst_lst_sta";
+ auto sName = format("{}_track_fst_lst_sta", GetTaskName());
auto sTitl = "First vs. last station index;ID^{last}_{station};ID^{first}_{station}";
fphTrkFstLstSta = MakeObj<H2D>(sName, sTitl, nBins, xMin, xMax, nBins, xMin, xMax);
}
@@ -227,7 +227,7 @@ void Qa::Init()
auto setNm = EHitSet::Input == hitSet ? "input" : "used";
auto setTl = EHitSet::Input == hitSet ? "Input" : "Used";
{ // XY
- auto name = format("ca_hit_{}_occupancy_xy", setNm);
+ auto name = format("{}_ca_hit_{}_occupancy_xy", GetTaskName(), setNm);
auto titl = format("{} hit occupancy in different stations in XY plane", setNm);
auto canv = CanvasConfig(name, titl);
for (int iSt = 0; iSt < nSt; ++iSt) {
@@ -238,7 +238,7 @@ void Qa::Init()
AddCanvasConfig(canv);
}
{ // ZX and ZY
- auto name = format("ca_hit_{}_occupancy_zx_zy", setNm);
+ auto name = format("{}_ca_hit_{}_occupancy_zx_zy", GetTaskName(), setNm);
auto titl = format("{} hit occupancy in different stations in ZX and ZY planes", setTl);
auto canv = CanvasConfig(name, titl);
{ // ZX
@@ -259,7 +259,7 @@ void Qa::Init()
}
}
if (kDebug) {
- auto name = format("ca_hit_usage_xy");
+ auto name = format("{}_ca_hit_usage_xy", GetTaskName());
auto titl = format("Hit usage in different stations in XY plane");
auto canv = CanvasConfig(name, titl);
for (int iSt = 0; iSt < nSt; ++iSt) {
@@ -273,7 +273,7 @@ void Qa::Init()
// Tracks canvas
{
- auto canv = CanvasConfig("ca_tracks", "Tracking output QA");
+ auto canv = CanvasConfig(format("{}_ca_tracks", GetTaskName()), "Tracking output QA");
{
auto pad = PadConfig(true, true, false, false, true);
pad.RegisterHistogram(fvphTrkPhiTheta[0], "colz");
diff --git a/algo/data/sts/Hit.h b/algo/data/sts/Hit.h
index ab1a733f2ca8372e1381433e61d9e94cdc8abb24..2df0d8370f9f0eb1edead87114bbaae959a8021e 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/bmon/Calibrate.cxx b/algo/detectors/bmon/Calibrate.cxx
index d5519854902605e3c84e77f584324d3811c88d04..fc2a17f29cfcf6ca0cef2b159f0bf6da3eb795b8 100644
--- a/algo/detectors/bmon/Calibrate.cxx
+++ b/algo/detectors/bmon/Calibrate.cxx
@@ -57,13 +57,6 @@ Calibrate::Calibrate(CalibrateSetup setup) : fSetup(setup), fSelectionBitsOffset
}
fChannelDeadTime = std::vector<double>(fChannelOffset.back(), std::numeric_limits<double>::quiet_NaN());
- // **** DEBUG: BEGIN
- for (size_t iO = 0; iO < fChannelOffset.size(); ++iO) {
- L_(info) << "Channel offset: diamond=" << iO << ", offset=" << fChannelOffset[iO];
- }
- L_(info) << "Selection mask: 0b" << std::bitset<32>(fSetup.selectionMask);
- L_(info) << "Bits offset: " << static_cast<uint32_t>(fSelectionBitsOffset);
- L_(info) << "Size of the channel dead time vector: " << fChannelDeadTime.size();
// **** DEBUG: END
}
diff --git a/algo/detectors/bmon/Hitfind.h b/algo/detectors/bmon/Hitfind.h
index 68d7d3a9b57e80099a6a4ce3854445622cca5384..0257a9b7483a588e0f06b89f54de477d08a1445a 100644
--- a/algo/detectors/bmon/Hitfind.h
+++ b/algo/detectors/bmon/Hitfind.h
@@ -45,11 +45,14 @@ namespace cbm::algo::bmon
/// \param iThread Index of thread
Output_t operator()(gsl::span<CbmBmonDigi> digisIn, uint32_t iThread = 0);
- private: // members
/// \brief Returns an index of the diamond by the address
/// \param address A hardware address of the digi
size_t GetDiamondIndex(uint32_t address) const { return ((fSelectionBitmask & address) >> fSelectionBitsOffset); }
+ /// \brief Gets diamond addresses vector
+ const PODVector<uint32_t>& GetDiamondAddresses() const { return fDiamondAddress; }
+
+ private: // members
uint32_t fNofThreads; ///< Number of threads
uint32_t fSelectionBitsOffset; ///< Number of bits to ther right from the first bit in the selection mask
uint32_t fSelectionBitmask; ///< Selection bitmask
diff --git a/algo/detectors/tof/Clusterizer.cxx b/algo/detectors/tof/Clusterizer.cxx
index caf6d34f758dcc976f9c70358f1c133605f1bbf6..fb36f7d47f3e924ac238574166128c4cf338b4dd 100644
--- a/algo/detectors/tof/Clusterizer.cxx
+++ b/algo/detectors/tof/Clusterizer.cxx
@@ -87,7 +87,6 @@ namespace cbm::algo::tof
auto digiIt = storDigi.begin();
while (1 < std::distance(digiIt, storDigi.end())) {
-
while (digiIt->first->GetSide() == std::next(digiIt)->first->GetSide()) { // Not one Digi of each end!
digiIt++;
if (2 > std::distance(digiIt, storDigi.end())) {
diff --git a/algo/detectors/tof/Hit.h b/algo/detectors/tof/Hit.h
index a7608fc72e94383cfc286e9fa7f005de276e2cb5..8e3a1a068531b531625f5118b2300643f0f3fc6d 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/detectors/tof/Hitfind.cxx b/algo/detectors/tof/Hitfind.cxx
index 119db83b1a43a26b37db80abe1e919ec39953f00..5abf282f3aecabec76ed6b541cd8f2877962360a 100644
--- a/algo/detectors/tof/Hitfind.cxx
+++ b/algo/detectors/tof/Hitfind.cxx
@@ -28,7 +28,6 @@ namespace cbm::algo::tof
for (int32_t Sm = 0; Sm < NbSm; Sm++) {
for (int32_t Rpc = 0; Rpc < NbRpc; Rpc++) {
-
auto par = std::make_unique<cbm::algo::tof::ClusterizerRpcPar>();
HitfindSetup::Rpc rpcPar = setup.rpcs[SmType][Sm * NbRpc + Rpc];
diff --git a/algo/evselector/RecoEventSelectorMonitor.h b/algo/evselector/RecoEventSelectorMonitor.h
new file mode 100644
index 0000000000000000000000000000000000000000..ae6524f804edd794bb3ec7bc0d3a377c64ba2946
--- /dev/null
+++ b/algo/evselector/RecoEventSelectorMonitor.h
@@ -0,0 +1,87 @@
+/* Copyright (C) 2025 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file RecoEventSelectorMonitor.h
+/// \date 28.01.2025
+/// \brief A monitor for reco event selector
+/// \author Sergei Zharko <s.zharko@gsi.de>
+
+#pragma once
+
+#include "CaMonitor.h"
+
+#include <boost/serialization/base_object.hpp>
+
+namespace cbm::algo::evselect
+{
+ /// \enum ECounter
+ /// \brief Counter keys for the event selector monitor
+ enum class ECounter
+ {
+ Timeslices, ///< number of processed timeslices
+ EventsTotal, ///< Total number of events processed
+ EventsNeStsHits, ///< Events with not enough STS hits
+ EventsNeTofHits, ///< Events with enough STS hits, but not enough TOF hits
+ EventsNeBmonHits, ///< Events with not enough BMon hits
+ EventsNeTracks, ///< Events with enough hits, but not enough tracks
+ EventsSelected, ///< Number of selected events
+ LambdaCandidates, ///< Number of lambda-candidates, returned by KFParticleFinder
+ END
+ };
+
+ /// \enum ETimer
+ /// \brief Timer keys for the event selector monitor
+ /* clang-format off */
+ enum class ETimer {
+ EventReconstruction,
+ BmonHitFinder,
+ StsHitFinder,
+ TofHitFinder,
+ TrdHitFinder,
+ TrackFinder,
+ V0Finder,
+ END
+ };
+ /* clang-format on */
+
+ /// \brief Specification of ca::MonitorData for the event selector
+ using MonitorData_t = ca::MonitorData<ECounter, ETimer>;
+
+ /// \class Monitor
+ /// \brief A monitor for the event selector
+ class Monitor : public ca::Monitor<ECounter, ETimer> {
+ public:
+ /// \brief Default constructor
+ Monitor() : ca::Monitor<ECounter, ETimer>("Event-selector Monitor")
+ {
+ SetCounterName(ECounter::Timeslices, "processed timeslices");
+ SetCounterName(ECounter::EventsTotal, "total events");
+ SetCounterName(ECounter::EventsNeStsHits, "events discarded by N STS hits");
+ SetCounterName(ECounter::EventsNeTofHits, "events discarded by N TOF hits");
+ SetCounterName(ECounter::EventsNeBmonHits, "events discarded by N BMon hits");
+ SetCounterName(ECounter::EventsNeTracks, "events discarded by N tracks");
+ SetCounterName(ECounter::LambdaCandidates, "potential lambda candidates");
+ SetCounterName(ECounter::EventsSelected, "selected events");
+
+ SetTimerName(ETimer::EventReconstruction, "event reconstruction");
+ SetTimerName(ETimer::BmonHitFinder, "hit finding in Bmon");
+ SetTimerName(ETimer::StsHitFinder, "hit finding in STS");
+ SetTimerName(ETimer::TofHitFinder, "hit finding in TOF");
+ SetTimerName(ETimer::TrdHitFinder, "hit finding in TRD");
+ SetTimerName(ETimer::TrackFinder, "track finding");
+ SetTimerName(ETimer::V0Finder, "V0 finding");
+
+ SetRatioKeys({ECounter::Timeslices});
+ }
+
+ private:
+ friend class boost::serialization::access;
+ template<typename Archive>
+ void serialize(Archive& ar, const unsigned int /*version*/)
+ {
+ ar& boost::serialization::base_object<ca::Monitor<ECounter, ETimer>>(*this);
+ }
+ };
+
+} // namespace cbm::algo::evselect
diff --git a/algo/global/ParFiles.cxx b/algo/global/ParFiles.cxx
index ac2c3b410accd1415bc4de89811785442eb9aca1..7b73fe77991c059b37d962a3694a552e36e278fd 100644
--- a/algo/global/ParFiles.cxx
+++ b/algo/global/ParFiles.cxx
@@ -45,6 +45,8 @@ ParFiles::ParFiles(uint32_t runId)
trd.hitfinder2d = "TrdHitfinder2DPar_mcbm2022.yaml";
ca.mainConfig = "TrackingChainConfig_mcbm2022.yaml";
+
+ kfp.V0FinderConfig = "kfp_lambda_v22a.yaml";
break;
case Setup::mCBM2024_03:
@@ -68,6 +70,8 @@ ParFiles::ParFiles(uint32_t runId)
trd.hitfinder2d = "TrdHitfinder2DPar_mcbm2024.yaml";
ca.mainConfig = "TrackingChainConfig_mcbm2024.yaml";
+
+ kfp.V0FinderConfig = "kfp_lambda_v24a.yaml";
break;
case Setup::mCBM2024_05:
@@ -91,6 +95,8 @@ ParFiles::ParFiles(uint32_t runId)
trd.hitfinder2d = "mcbm2024_05/TrdHitfinder2DPar.yaml";
ca.mainConfig = "mcbm2024_05/TrackingChainConfig.yaml";
+
+ kfp.V0FinderConfig = "kfp_lambda_v24b.yaml";
break;
case Setup::mCBM2025_02:
diff --git a/algo/global/ParFiles.h b/algo/global/ParFiles.h
index 691e27b39bf8af9b64b89cc209422f786048e1e7..0fe7b3094b1e25eaa0d7f3ccb207a6977f5d2836 100644
--- a/algo/global/ParFiles.h
+++ b/algo/global/ParFiles.h
@@ -54,6 +54,10 @@ namespace cbm::algo
struct {
fs::path mainConfig;
} ca;
+
+ struct {
+ fs::path V0FinderConfig;
+ } kfp;
};
} // namespace cbm::algo
diff --git a/algo/global/Reco.cxx b/algo/global/Reco.cxx
index 8460ea28bf34fa71071601370839602433199058..0f192df1b6c68934269a1eea7e5561b73dd6aa5d 100644
--- a/algo/global/Reco.cxx
+++ b/algo/global/Reco.cxx
@@ -22,6 +22,7 @@
#include "ca/core/data/CaTrack.h"
#include "compat/OpenMP.h"
#include "evbuild/Config.h"
+#include "kfp/KfpV0FinderChain.h"
#include "much/Unpack.h"
#include "qa/QaManager.h"
#include "qa/hitfind/BmonHitfindQa.h"
@@ -44,10 +45,6 @@
#include <xpu/host.h>
-// DEBUG: BEGIN
-#include <set>
-// DEBUG: END
-
using namespace cbm::algo;
using fles::Subsystem;
@@ -211,21 +208,6 @@ void Reco::Init(const Options& opts)
fStsHitFinder->SetParameters(hitFinderPars);
}
- // BMON Hitfinder
- if (Opts().Has(fles::Subsystem::BMON) && Opts().Has(Step::LocalReco)) {
- auto calibSetup = yaml::ReadFromFile<bmon::CalibrateSetup>(opts.ParamsDir() / parFiles.bmon.calibrate);
- fBmonCalibrator = std::make_unique<bmon::Calibrate>(calibSetup);
-
- auto hitfindSetup = yaml::ReadFromFile<bmon::HitfindSetup>(opts.ParamsDir() / parFiles.bmon.hitfinder);
- fBmonHitFinder = std::make_unique<bmon::Hitfind>(hitfindSetup);
-
- if (fQaManager != nullptr && Opts().Has(QaStep::RecoBmon)) {
- fBmonHitFinderQa = std::make_unique<bmon::HitfindQa>(fQaManager, "BmonHitfindEvent");
- fBmonHitFinderQa->InitParameters(calibSetup, hitfindSetup);
- fBmonHitFinderQa->Init();
- }
- }
-
// TOF Hitfinder
if (Opts().Has(fles::Subsystem::TOF) && Opts().Has(Step::LocalReco)) {
auto calibSetup = yaml::ReadFromFile<tof::CalibrateSetup>(opts.ParamsDir() / parFiles.tof.calibrate);
@@ -241,19 +223,65 @@ void Reco::Init(const Options& opts)
fTrdHitfind = std::make_unique<trd::Hitfind>(setup, setup2d);
}
+ // Digi event reconstruction:
+ {
+ fbReconstructDigiEvents = Opts().ReconstructDigiEvents();
+ // It makes no sence to reconstruct an event, if there is no STS, TRD or TOF
+ fbReconstructDigiEvents &= Opts().Has(fles::Subsystem::STS);
+ fbReconstructDigiEvents &= Opts().Has(fles::Subsystem::TRD);
+ fbReconstructDigiEvents &= Opts().Has(fles::Subsystem::TOF);
+ fbReconstructDigiEvents &= Opts().Has(fles::Subsystem::BMON);
+ }
+
// Tracking
if (Opts().Has(Step::Tracking)) {
if (fQaManager != nullptr && Opts().Has(QaStep::Tracking)) {
- fTracking = std::make_unique<TrackingChain>(fQaManager, "CaTimeslice");
+ fTracking = std::make_unique<TrackingChain>(ECbmRecoMode::Timeslice, fQaManager, "CaTimeslice");
}
else {
- fTracking = std::make_unique<TrackingChain>();
+ fTracking = std::make_unique<TrackingChain>(ECbmRecoMode::Timeslice);
}
fTracking->RegisterSetup(pTrackingSetup);
fTracking->SetContext(&fContext);
fTracking->Init();
}
+ if (fbReconstructDigiEvents) {
+ fEvSelectingMonitor.Reset();
+
+ // BMON hit finding in event reconstruction
+ auto bmonCalSetup = yaml::ReadFromFile<bmon::CalibrateSetup>(opts.ParamsDir() / parFiles.bmon.calibrate);
+ auto bmonHitSetup = yaml::ReadFromFile<bmon::HitfindSetup>(opts.ParamsDir() / parFiles.bmon.hitfinder);
+ fBmonCalibrator = std::make_unique<bmon::Calibrate>(bmonCalSetup);
+ fBmonHitFinder = std::make_unique<bmon::Hitfind>(bmonHitSetup);
+ if (fQaManager != nullptr && Opts().Has(QaStep::RecoBmon)) {
+ fBmonHitFinderQa = std::make_unique<bmon::HitfindQa>(fQaManager, "BmonHitfindEvent");
+ fBmonHitFinderQa->InitParameters(bmonCalSetup, bmonHitSetup);
+ fBmonHitFinderQa->Init();
+ }
+
+ // Tracking in event reconstruction
+ if (fQaManager != nullptr && Opts().Has(QaStep::Tracking)) {
+ fTrackingEvent = std::make_unique<TrackingChain>(ECbmRecoMode::EventByEvent, fQaManager, "CaEvent");
+ }
+ else {
+ fTrackingEvent = std::make_unique<TrackingChain>(ECbmRecoMode::EventByEvent);
+ }
+ fTrackingEvent->RegisterSetup(pTrackingSetup);
+ fTrackingEvent->SetContext(&fContext);
+ fTrackingEvent->Init();
+
+ if (fQaManager != nullptr && Opts().Has(QaStep::V0Finder)) {
+ fV0Finder = std::make_unique<V0FinderChain>(fQaManager);
+ }
+ else {
+ fV0Finder = std::make_unique<V0FinderChain>();
+ }
+ fV0Finder->SetContext(&fContext);
+ fV0Finder->SetBmonDefinedAddresses(fBmonHitFinder->GetDiamondAddresses());
+ fV0Finder->Init();
+ }
+
// Initialize the QA manager
if (fQaManager != nullptr) {
fQaManager->Init();
@@ -393,23 +421,19 @@ RecoResults Reco::Run(const fles::Timeslice& ts)
QueueEvbuildMetrics(evbuildMonitor);
}
- // ***** DEBUG: BEGIN
- if constexpr (0) {
- int nEvents = events.size();
- for (int iE = 0; iE < nEvents; ++iE) {
- const auto& event = events[iE];
- // Calibrate TOF digis:
- auto [bmonDigis, bmonCalMonitor] = (*fBmonCalibrator)(event.fBmon);
- auto [bmonHits, hitmonitor, digiIndices] = (*fBmonHitFinder)(bmonDigis);
- if (fBmonHitFinderQa != nullptr) {
- fBmonHitFinderQa->RegisterDigis(&bmonDigis);
- fBmonHitFinderQa->RegisterHits(&bmonHits);
- fBmonHitFinderQa->RegisterDigiIndices(&digiIndices);
- fBmonHitFinderQa->Exec();
- }
+ // --- Reconstruct and select digi events
+ if (fbReconstructDigiEvents) {
+ fEvSelectingMonitor.IncrementCounter(evselect::ECounter::Timeslices);
+ fEvSelectingMonitor.IncrementCounter(evselect::ECounter::EventsTotal, events.size());
+ for (auto& event : events) {
+ fEvSelectingMonitor.StartTimer(evselect::ETimer::EventReconstruction);
+ event.fSelectionTriggers = ReconstructEvent(event);
+ fEvSelectingMonitor.StopTimer(evselect::ETimer::EventReconstruction);
}
+ auto v0FinderMonitor = fV0Finder->GetMonitor();
+ fEvSelectingMonitor.IncrementCounter(evselect::ECounter::LambdaCandidates,
+ v0FinderMonitor.GetCounterValue(kfp::ECounter::KfpLambdaCandidates));
}
- // ***** DEBUG: END
// --- Filter data for output
if (Opts().HasOutput(RecoData::DigiTimeslice)) {
@@ -459,8 +483,18 @@ void Reco::Finalize()
fStsHitFinder->Finalize();
}
if (fTracking) {
+ L_(info) << "Track finding in a timeslice:";
fTracking->Finalize();
}
+ if (fTrackingEvent) {
+ L_(info) << "Track finding in digi events:";
+ fTrackingEvent->Finalize();
+ }
+ if (fbReconstructDigiEvents) {
+ fV0Finder->Finalize();
+ L_(info) << fEvSelectingMonitor.ToString();
+ }
+
if (Opts().Profiling() >= ProfilingSummary) {
L_(info) << MakeReportSubtimers("Run Summary", fTimesliceTimesAcc) << "\n"
@@ -487,6 +521,95 @@ void Reco::PrintTimings(xpu::timings& timings)
}
}
+CbmEventTriggers Reco::ReconstructEvent(const DigiEvent& digiEvent)
+{
+ CbmEventTriggers triggers(0);
+ RecoResults recoEvent;
+ //* BMON hit reconstruction
+ {
+ fEvSelectingMonitor.StartTimer(evselect::ETimer::BmonHitFinder);
+ auto [calDigis, calMonitor] = (*fBmonCalibrator)(digiEvent.fBmon);
+ auto [hits, hitMonitor, digiIndices] = (*fBmonHitFinder)(calDigis);
+ fEvSelectingMonitor.StopTimer(evselect::ETimer::BmonHitFinder);
+ if (fBmonHitFinderQa != nullptr) {
+ fBmonHitFinderQa->RegisterDigis(&calDigis);
+ fBmonHitFinderQa->RegisterHits(&hits);
+ fBmonHitFinderQa->RegisterDigiIndices(&digiIndices);
+ fBmonHitFinderQa->Exec();
+ }
+ recoEvent.bmonHits = std::move(hits);
+ }
+
+ //* STS hit reconstruction
+ {
+ fEvSelectingMonitor.StartTimer(evselect::ETimer::StsHitFinder);
+ auto stsResults = (*fStsHitFinder)(digiEvent.fSts);
+ fEvSelectingMonitor.StopTimer(evselect::ETimer::StsHitFinder);
+ if (stsResults.hits.NElements() < 4) { // TODO: Provide a config for cuts (testing mode for now)
+ fEvSelectingMonitor.IncrementCounter(evselect::ECounter::EventsNeStsHits);
+ return triggers;
+ }
+ recoEvent.stsHits = stsResults.hits;
+ }
+
+ //* TOF hit reconstruction
+ {
+ fEvSelectingMonitor.StartTimer(evselect::ETimer::TofHitFinder);
+ auto [caldigis, calmonitor] = (*fTofCalibrator)(digiEvent.fTof);
+ auto [hits, hitmonitor, digiindices] = (*fTofHitFinder)(caldigis);
+ fEvSelectingMonitor.StopTimer(evselect::ETimer::TofHitFinder);
+ if (hits.NElements() < 2) { // TODO: Provide a config for cuts (testing mode for now)
+ fEvSelectingMonitor.IncrementCounter(evselect::ECounter::EventsNeTofHits);
+ return triggers;
+ }
+ recoEvent.tofHits = std::move(hits);
+ }
+
+ //* TRD hit reconstruction
+ {
+ // FIXME: additional copy of digis, figure out how to pass 1d + 2d digis at once to hitfinder
+ fEvSelectingMonitor.StartTimer(evselect::ETimer::TrdHitFinder);
+ const auto& digis1d = digiEvent.fTrd;
+ const auto& digis2d = digiEvent.fTrd2d;
+ PODVector<CbmTrdDigi> allDigis{};
+ allDigis.reserve(digis1d.size() + digis2d.size());
+ std::copy(digis1d.begin(), digis1d.end(), std::back_inserter(allDigis));
+ std::copy(digis2d.begin(), digis2d.end(), std::back_inserter(allDigis));
+ auto trdResults = (*fTrdHitfind)(allDigis);
+ fEvSelectingMonitor.StopTimer(evselect::ETimer::TrdHitFinder);
+ recoEvent.trdHits = std::move(std::get<0>(trdResults));
+ }
+
+ //* Tracking
+ {
+ fEvSelectingMonitor.StartTimer(evselect::ETimer::TrackFinder);
+ TrackingChain::Input_t input{.stsHits = recoEvent.stsHits,
+ .tofHits = recoEvent.tofHits,
+ .trdHits = recoEvent.trdHits};
+ TrackingChain::Output_t output = fTrackingEvent->Run(input);
+ recoEvent.tracks = std::move(output.tracks);
+ recoEvent.trackStsHitIndices = std::move(output.stsHitIndices);
+ recoEvent.trackTofHitIndices = std::move(output.tofHitIndices);
+ recoEvent.trackTrdHitIndices = std::move(output.trdHitIndices);
+ fEvSelectingMonitor.StopTimer(evselect::ETimer::TrackFinder);
+ if (recoEvent.tracks.size() < 2) { // Reject all events with less then two tracks
+ fEvSelectingMonitor.IncrementCounter(evselect::ECounter::EventsNeTracks);
+ return triggers;
+ }
+ }
+
+ //* V0-selector
+ fEvSelectingMonitor.StartTimer(evselect::ETimer::V0Finder);
+ triggers = fV0Finder->ProcessEvent(recoEvent);
+ if (triggers.Test(CbmEventTriggers::ETrigger::Lambda)) {
+ L_(info) << "!!! Found event with potential lambda candidates";
+ fEvSelectingMonitor.IncrementCounter(evselect::ECounter::EventsSelected);
+ }
+ fEvSelectingMonitor.StopTimer(evselect::ETimer::V0Finder);
+ return triggers;
+}
+
+
template<class Unpacker>
auto Reco::RunUnpacker(const std::unique_ptr<Unpacker>& unpacker, const fles::Timeslice& ts) -> UnpackResult_t<Unpacker>
{
diff --git a/algo/global/Reco.h b/algo/global/Reco.h
index 6182decdad9d289e7a654aabbb6b5db8650efe29..35cb5c34626c8dc84ba994478e83a305522cfd98 100644
--- a/algo/global/Reco.h
+++ b/algo/global/Reco.h
@@ -5,11 +5,14 @@
#include "AlgoTraits.h"
#include "SubChain.h"
+#include "evselector/RecoEventSelectorMonitor.h"
#include "global/RecoResults.h"
#include <xpu/host.h>
// fwd declarations
+class CbmEventTriggers;
+
namespace fles
{
class Timeslice;
@@ -20,6 +23,7 @@ namespace cbm::algo
class HistogramSender;
class Options;
class TrackingChain;
+ class V0FinderChain;
template<class M>
struct UnpackMonitor;
@@ -130,13 +134,16 @@ namespace cbm::algo
void Init(const Options&);
RecoResults Run(const fles::Timeslice&);
+
+ CbmEventTriggers ReconstructEvent(const DigiEvent& event);
void Finalize();
void PrintTimings(xpu::timings&);
void QueueProcessingExtraMetrics(const ProcessingExtraMonitor&);
private:
- bool fInitialized = false;
+ bool fInitialized = false;
+ bool fbReconstructDigiEvents = false;
ChainContext fContext;
xpu::timings fTimesliceTimesAcc;
std::shared_ptr<HistogramSender> fSender;
@@ -178,7 +185,14 @@ namespace cbm::algo
std::unique_ptr<evbuild::EventbuildChain> fEventBuild;
// Tracking
- std::unique_ptr<TrackingChain> fTracking;
+ std::unique_ptr<TrackingChain> fTracking; ///< Tracking in timeslice
+ std::unique_ptr<TrackingChain> fTrackingEvent; ///< Tracking in event
+
+ // V0-finding
+ std::unique_ptr<V0FinderChain> fV0Finder; ///< V0-finding chain (in event or a bunch of events)
+
+ // Event selection
+ evselect::Monitor fEvSelectingMonitor; ///< Monitor for event selecting
// QA
std::unique_ptr<qa::Manager> fQaManager;
diff --git a/algo/global/RecoResults.h b/algo/global/RecoResults.h
index eb64fa2368569a398c21da0137a4678b13e46490..3bcee8be6acd66388ab8e0043313888c0f35c615 100644
--- a/algo/global/RecoResults.h
+++ b/algo/global/RecoResults.h
@@ -12,6 +12,7 @@
#include "DigiData.h"
#include "PartitionedSpan.h"
#include "PartitionedVector.h"
+#include "bmon/Hit.h"
#include "ca/core/data/CaTrack.h"
#include "ca/core/utils/CaVector.h"
#include "sts/Cluster.h"
@@ -26,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;
@@ -41,9 +44,11 @@ namespace cbm::algo
PartitionedSpan<sts::Hit> stsHits;
PartitionedVector<tof::Hit> tofHits;
PartitionedVector<trd::Hit> trdHits;
+ PartitionedVector<bmon::Hit> bmonHits;
ca::Vector<ca::Track> tracks;
- ca::Vector<std::vector<std::pair<uint32_t, uint32_t>>> trackStsHitIndices;
- ca::Vector<std::vector<std::pair<uint32_t, uint32_t>>> trackTofHitIndices;
+ 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/CMakeLists.txt b/algo/kfp/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..cc41f095bbc59bdb98f1ebcafb9ac8c09ec67927
--- /dev/null
+++ b/algo/kfp/CMakeLists.txt
@@ -0,0 +1 @@
+add_subdirectory(interface)
diff --git a/algo/kfp/KfpV0Finder.cxx b/algo/kfp/KfpV0Finder.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..c95965bc847347caccc52b1bcf23a0eaeda8fd92
--- /dev/null
+++ b/algo/kfp/KfpV0Finder.cxx
@@ -0,0 +1,441 @@
+/* Copyright (C) 2025 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file KfpV0Finder.cxx
+/// \date 01.02.2025
+/// \brief A V0 finding algorithm (implementation)
+/// \author Sergei Zharko <s.zharko@gsi.de>
+
+#include "kfp/KfpV0Finder.h"
+
+#include "global/RecoResults.h"
+
+#include <algorithm>
+#include <limits>
+#include <sstream>
+
+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.fCharge / (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);
+ }
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void V0Finder::CollectDca(const RecoResults& recoEvent)
+{
+ const auto& stsHitIndices = recoEvent.trackStsHitIndices;
+ for (size_t iTrk = 0; iTrk < stsHitIndices.size(); ++iTrk) {
+ const auto& stsHitIndicesInTrack = stsHitIndices[iTrk];
+ if (stsHitIndicesInTrack.size() < 2) { // DCA cannot be estimated
+ fEventMonitor.IncrementCounter(ECounter::TracksWoStsHits);
+ continue;
+ }
+ auto& particleInfo = fvParticleInfo[iTrk];
+ auto [iPtFst, iHitFst] = stsHitIndicesInTrack[0];
+ auto [iPtSnd, iHitSnd] = stsHitIndicesInTrack[1];
+ particleInfo.fDca = EstimateDca(recoEvent.stsHits[iPtFst][iHitFst], recoEvent.stsHits[iPtSnd][iHitSnd]);
+ AssignPid(particleInfo);
+ }
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void V0Finder::CollectT0(gsl::span<const bmon::Hit> bmonHits)
+{
+ fvT0s.clear();
+ fvT0s.reserve(bmonHits.size());
+ std::transform(bmonHits.begin(), bmonHits.end(), std::back_inserter(fvT0s),
+ [&](const auto& h) { return h.GetTime(); });
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+double V0Finder::EstimateDca(const sts::Hit& fst, const sts::Hit& snd) const
+{
+ double factor{(fst.Z() - fOrigin[2]) / (snd.Z() - fst.Z())};
+ double dcaX{fst.X() - fOrigin[0] - factor * (snd.X() - fst.X())};
+ double dcaY{fst.Y() - fOrigin[1] - factor * (snd.Y() - fst.Y())};
+ return std::sqrt(dcaX * dcaX + dcaY * dcaY);
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+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.Z() - 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);
+ 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->Init(kfpTracksFst, kfpTracksLst);
+ 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::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
+ fEventMonitor.StartTimer(ETimer::CollectT0);
+ CollectT0(recoEvent.bmonHits[fBmonPartitionIndex]);
+ if (fvT0s.empty()) {
+ fEventMonitor.IncrementCounter(ECounter::EventsWoTzero);
+ return res;
+ }
+ fEventMonitor.StopTimer(ETimer::CollectT0);
+
+ // ----- Estimate DCA of tracks and assign PID
+ // If a track has less then two STS hits, and undefined DCA value is stored
+ fEventMonitor.StartTimer(ETimer::CollectDca);
+ CollectDca(recoEvent);
+ fEventMonitor.StopTimer(ETimer::CollectDca);
+
+ // ----- 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::KfpEventsLambdaCand);
+ 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;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+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
new file mode 100644
index 0000000000000000000000000000000000000000..9605b581b2b980b0825cb68de1689fc68da2ad5e
--- /dev/null
+++ b/algo/kfp/KfpV0Finder.h
@@ -0,0 +1,276 @@
+/* Copyright (C) 2025 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file KfpV0Finder.h
+/// \date 01.02.2025
+/// \brief A V0 finding algorithm
+/// \author Sergei Zharko <s.zharko@gsi.de>
+
+#pragma once
+
+#include "CbmEventTriggers.h"
+#include "KFParticleTopoReconstructor.h"
+#include "global/RecoResults.h"
+#include "kfp/KfpV0FinderMonitor.h"
+
+#include <limits>
+#include <memory>
+
+
+namespace cbm::algo::kfp
+{
+ /// \class V0Finder
+ /// \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
+ };
+
+ //* Framework and physical constants (public)
+ 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
+
+ /// \brief Default constructor
+ V0Finder() = default;
+
+ /// \brief Copy constructor
+ V0Finder(const V0Finder&) = delete;
+
+ /// \brief Move constructor
+ V0Finder(V0Finder&&) = delete;
+
+ /// \brief Destructor
+ ~V0Finder() = default;
+
+ /// \brief Copy assignment operator
+ V0Finder& operator=(const V0Finder&) = delete;
+
+ /// \brief Move assignment operator
+ V0Finder& operator=(V0Finder&&) = delete;
+
+ /// \brief Adds particle to reconstruction list
+ /// \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 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 Accessor to topology reconstructor
+ const std::unique_ptr<KFParticleTopoReconstructor>& GetTopoReconstructor() const { return fpTopoReconstructor; }
+
+ /// \brief Gets track parameters
+ const auto& GetTrackAssignedParams() const { return fvTrackParam; }
+
+ /// \brief Initializes the instance (called in the beginning of the run)
+ void Init();
+
+ /// \brief Processes a reconstructed data sample, returns a collection of fired triggers
+ CbmEventTriggers Process(const RecoResults& recoEvent);
+
+ /// \brief Sets an address of a reference BMON diamond
+ /// \param iPartition An index of the BMON hit partition: used to selected the hits from the particular config
+ void SetBmonPartitionIndex(int iPartition) { fBmonPartitionIndex = iPartition; }
+
+ /// \param dca DCA [cm]
+ void SetMinPionDca(double dca) { fMinPionDca = dca; }
+
+ /// \brief Sets minimal proton DCA to primary vertex
+ /// \param dca DCA [cm]
+ void SetMinProtonDca(double dca) { fMinProtonDca = dca; }
+
+ /// \brief Sets origin
+ /// \param x X-coordinate of the origin [cm]
+ /// \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 = {float(x), float(y), float(z)}; }
+
+ /// \brief Sets minimal pion DCA to primary vertex
+ /// \brief Sets pion velocity range
+ /// \param vMin Minimal velocity [cm/ns]
+ /// \param vMax Maximal velocity [cm/ns]
+ void SetPionVelocityRange(double vMin, double vMax)
+ {
+ fMinBetaPion = vMin / kSpeedOfLight;
+ fMaxBetaPion = vMax / kSpeedOfLight;
+ }
+
+ /// \brief Sets proton velocity range
+ /// \param vMin Minimal velocity [cm/ns]
+ /// \param vMax Maximal velocity [cm/ns]
+ void SetProtonVelocityRange(double vMin, double vMax)
+ {
+ fMinBetaProton = vMin / kSpeedOfLight;
+ fMaxBetaProton = vMax / kSpeedOfLight;
+ }
+
+ /// \brief Sets the assigned PDG for primary particles
+ /// \param pdg PDG code of the particle
+ void SetPrimaryAssignedPdg(int pdg) { fPrimaryAssignedPdg = pdg; }
+
+ /// \brief Assignes an uncertainty to the momentum measurement
+ /// \param uncertainty Relative uncertainty ( = sqrt(var(q/p)) / (q/p))
+ void SetQpAssignedUncertainty(double uncertainty) { fQpAssignedUncertainty = uncertainty; }
+
+ /// \brief Sets an offset to t0
+ /// \param offset An offset [ns]
+ void SetTzeroOffset(double offset) { fTzeroOffset = offset; }
+
+ /// \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); }
+
+ /// \brief Sets cut on \f$\chi^2_{prim}\f$ of each track for 2-daughter decays
+ /// \param cut Cut value
+ void SetChiPrimaryCut2D(float cut) { GetKFParticleFinder()->SetChiPrimaryCut2D(cut); }
+
+ /// \brief Sets cut on \f$\chi^2_{geo}\f$ for 2-daughter decays
+ /// \param cut Cut value
+ void SetChi2Cut2D(float cut) { GetKFParticleFinder()->SetChi2Cut2D(cut); }
+
+ /// \brief Sets cut on \f$l/\Delta l\f$ for 2-daughter decays
+ /// \param cut Cut value
+ 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
+ ///
+ /// If multiple T0-s are found, the routine will run multiple times, until V0-candidates are found
+ void CollectT0(gsl::span<const bmon::Hit> bmonHits);
+
+ /// \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
+ /// \param snd second STS hit
+ /// \return dca [cm]
+ double EstimateDca(const sts::Hit& fst, const sts::Hit& snd) const;
+
+ /// \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 (private)
+ static constexpr bool kUseAverageSpeed{false}; ///< If an average speed of tof hits is used
+
+ 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
+ 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
+
+ //* 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
+
+ /// \brief An instance of the topology reconstructor
+ std::unique_ptr<KFParticleTopoReconstructor> fpTopoReconstructor{std::make_unique<KFParticleTopoReconstructor>()};
+ };
+} // namespace cbm::algo::kfp
diff --git a/algo/kfp/KfpV0FinderChain.cxx b/algo/kfp/KfpV0FinderChain.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..96dc60a593ee94df51cfde7a736e4f1c15c541a8
--- /dev/null
+++ b/algo/kfp/KfpV0FinderChain.cxx
@@ -0,0 +1,160 @@
+/* Copyright (C) 2025 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file KfpV0FinderChain.cxx
+/// \date 01.02.2025
+/// \brief A chain for V0 finding (implementation)
+/// \author Sergei Zharko <s.zharko@gsi.de>
+
+#include "kfp/KfpV0FinderChain.h"
+
+#include "global/ParFiles.h"
+#include "kfp/KfpV0FinderConfig.h"
+#include "log/AlgoFairloggerCompat.h"
+#include "yaml/Yaml.h"
+
+#include <sstream>
+
+using cbm::algo::V0FinderChain;
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+V0FinderChain::V0FinderChain(const std::unique_ptr<qa::Manager>& qaManager)
+ : fpFinderQa(qaManager != nullptr ? std::make_unique<kfp::V0FinderQa>(qaManager, "V0Finder") : nullptr)
+{
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void V0FinderChain::Finalize()
+{
+ GetMonitor(); // A hack to update the run monitor
+ L_(info) << fMonitorRun.ToString();
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+cbm::algo::kfp::V0FinderMonitorData_t V0FinderChain::GetMonitor()
+{
+ cbm::algo::kfp::V0FinderMonitorData_t monitorData = fMonitorTimeslice;
+ fMonitorTimeslice.Reset();
+ fMonitorRun.AddMonitorData(monitorData);
+ return monitorData;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void V0FinderChain::Init()
+try {
+ L_(info) << "kfp::V0FinderChain: initializing the V0-finder chain ...";
+
+ // ----- Read configuration file
+ ParFiles parFiles(Opts().RunId());
+ auto config = yaml::ReadFromFile<kfp::V0FinderConfig>(Opts().ParamsDir() / parFiles.kfp.V0FinderConfig);
+ L_(info) << config.ToString();
+
+ //* Check the V0 type
+ if (config.reconstructPdg != 3122) { // At the moment only Lambda analysis is possible
+ std::stringstream msg;
+ msg << "kfp::V0FinderChain: at the moment only lambda finding is possible. Provided PDG: " << config.reconstructPdg;
+ throw std::runtime_error(msg.str());
+ }
+
+ //* Define daughter particles
+ auto& particles{config.cuts.particles};
+ int iPion{-1};
+ int iProton{-1};
+ for (int iParticle = 0; iParticle < int(particles.size()); ++iParticle) {
+ const auto& particle = particles[iParticle];
+ auto CheckOutParticle = [&](int pdg, int& iParticleToFind) {
+ if (particle.pdg == pdg) {
+ if (iParticleToFind == -1) {
+ iParticleToFind = iParticle;
+ }
+ else {
+ std::stringstream msg;
+ msg << "kfp::V0FinderChain: entry for pdg= " << pdg << " is defined more then one time in the ";
+ msg << "config.cuts.particles";
+ throw std::runtime_error(msg.str());
+ }
+ }
+ };
+ CheckOutParticle(-211, iPion);
+ CheckOutParticle(2212, iProton);
+ }
+ if (iProton == -1 || iPion == -1) {
+ std::stringstream msg;
+ msg << "kfp::V0FinderChain: config cuts/particles: either pion or proton settings are not found";
+ throw std::runtime_error(msg.str());
+ }
+ const auto& pion{particles[iPion]};
+ const auto& proton{particles[iProton]};
+
+ // ----- Define a BMON diamond
+ if (fBmonDefinedAddresses.empty()) {
+ throw std::runtime_error("kfp::V0FinderChain: BMON available addresses were not set");
+ }
+ int iBmonPartitionSelect = -1;
+ for (int iPart = 0; iPart < static_cast<int>(fBmonDefinedAddresses.size()); ++iPart) {
+ if (config.bmonAddress == fBmonDefinedAddresses[iPart]) {
+ iBmonPartitionSelect = iPart;
+ break;
+ }
+ }
+ if (iBmonPartitionSelect < 0) {
+ std::stringstream msg;
+ msg << "kfp::V0FinderChain: a reference BMON address ( " << std::hex << config.bmonAddress << std::dec
+ << " ) differs from ones, provided by hitfinder. Please check your configuration";
+ throw std::runtime_error(msg.str());
+ }
+
+ // ----- Set the V0-finder properties
+ // TODO: In future, there are will be a several instances of the V0Finder, each for a particular thread
+ {
+ //* Set particle PID properties
+ fFinder.SetBmonPartitionIndex(iBmonPartitionSelect);
+ fFinder.SetMinPionDca(pion.minDca);
+ fFinder.SetMinProtonDca(proton.minDca);
+ fFinder.SetPionVelocityRange(pion.minVelocity, pion.maxVelocity);
+ fFinder.SetProtonVelocityRange(proton.minVelocity, proton.maxVelocity);
+
+ //* Set KFParticleFinder properties
+ auto& kfpCuts{config.cuts.kfp};
+ fFinder.SetLdLCut2D(kfpCuts.minDecayLDL);
+ fFinder.SetLCut(kfpCuts.minDecayLength);
+ fFinder.SetChi2Cut2D(kfpCuts.maxChi2NdfGeo);
+ fFinder.SetChiPrimaryCut2D(kfpCuts.maxChi2NdfPrim);
+
+ //* Set other properties
+ fFinder.SetTzeroOffset(config.tZeroOffset);
+ fFinder.SetQpAssignedUncertainty(config.qpAssignedUncertainty);
+ fFinder.AddDecayToReconstructionList(config.reconstructPdg);
+ fFinder.SetPrimaryAssignedPdg(config.primaryAssignedPdg);
+
+ //* Init the V0 finder
+ fFinder.Init();
+ }
+
+ if (fpFinderQa != nullptr) {
+ fpFinderQa->Init();
+ }
+
+ L_(info) << "kfp::V0FinderChain: initializing the V0-finder chain ... done";
+}
+catch (const std::exception& err) {
+ L_(info) << "kfp::V0FinderChain: initializing the V0-finder chain ... failed. Reason: " << err.what();
+ throw std::runtime_error("initialization of V0-finder chain failed");
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+V0FinderChain::EventOutput V0FinderChain::ProcessEvent(const RecoResults& recoEvent)
+{
+ EventOutput res = fFinder.Process(recoEvent);
+ if (fpFinderQa != nullptr) {
+ fpFinderQa->Exec(recoEvent, fFinder);
+ }
+ fMonitorTimeslice.AddMonitorData(fFinder.GetEventMonitor());
+ return res;
+}
diff --git a/algo/kfp/KfpV0FinderChain.h b/algo/kfp/KfpV0FinderChain.h
new file mode 100644
index 0000000000000000000000000000000000000000..95f10ebbab5e8c6c6086e8ce00d5a1c3ad4fa2ab
--- /dev/null
+++ b/algo/kfp/KfpV0FinderChain.h
@@ -0,0 +1,81 @@
+/* Copyright (C) 2025 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file KfpV0FinderChain.h
+/// \date 01.02.2025
+/// \brief A chain for V0 finding
+/// \author Sergei Zharko <s.zharko@gsi.de>
+
+#pragma once
+
+#include "CbmEventTriggers.h"
+#include "PODVector.h"
+#include "base/SubChain.h"
+#include "global/RecoResults.h"
+#include "kfp/KfpV0Finder.h"
+#include "kfp/KfpV0FinderMonitor.h"
+#include "kfp/KfpV0FinderQa.h"
+
+namespace cbm::algo
+{
+ namespace qa
+ {
+ class Manager;
+ }
+} // namespace cbm::algo
+
+namespace cbm::algo
+{
+ /// \class V0FinderChain
+ /// \brief A chain for the V0 finder
+ class V0FinderChain : public SubChain {
+ public:
+ using EventOutput = CbmEventTriggers;
+
+ /// \brief Default constructor
+ V0FinderChain() = default;
+
+ /// \brief Constructor from parameters
+ /// \param pQaManager A QA-manager
+ V0FinderChain(const std::unique_ptr<qa::Manager>& qaManager);
+
+ /// \brief Copy constructor
+ V0FinderChain(const V0FinderChain&) = delete;
+
+ /// \brief Move constructor
+ V0FinderChain(V0FinderChain&&) = delete;
+
+ /// \brief Destructor
+ ~V0FinderChain() = default;
+
+ /// \brief Copy assignment operator
+ V0FinderChain& operator=(const V0FinderChain&) = delete;
+
+ /// \brief Move assignment operator
+ V0FinderChain& operator=(V0FinderChain&&) = delete;
+
+ /// \brief Finalizes the instance (called in the end of the run)
+ void Finalize();
+
+ /// \brief Gets a monitor
+ kfp::V0FinderMonitorData_t GetMonitor();
+
+ /// \brief Sets BMON diamond addresses array
+ /// \note The addresses must be taken from bmon::Hitfind::GetDiamondAddresses()
+ void SetBmonDefinedAddresses(const PODVector<uint32_t>& addresses) { fBmonDefinedAddresses = addresses; }
+
+ /// \brief Initializes the instance (called in the beginning of the run)
+ void Init();
+
+ /// \brief Processes an event, returns a collection of fired triggers
+ EventOutput ProcessEvent(const RecoResults& recoEvent);
+
+ private:
+ kfp::V0Finder fFinder; ///< Instance of the V0-finding algorithm
+ kfp::V0FinderMonitor fMonitorRun; ///< Monitor per run
+ kfp::V0FinderMonitorData_t fMonitorTimeslice; ///< Monitor per timeslice
+ std::unique_ptr<kfp::V0FinderQa> fpFinderQa{nullptr}; ///< QA module
+ PODVector<uint32_t> fBmonDefinedAddresses; ///< Available addresses of BMON
+ };
+} // namespace cbm::algo
diff --git a/algo/kfp/KfpV0FinderConfig.cxx b/algo/kfp/KfpV0FinderConfig.cxx
index 54e08ec3234d4fd4126cfc612e07155bc7dcebf2..d74e4d8b2ca7f53755e941dcde73030ca27fbb4c 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/KfpV0FinderConfig.h b/algo/kfp/KfpV0FinderConfig.h
index dacc3705d843874c33d196909dd85b2ac655294e..4bd992e0cb406145ad29495f89acfc00827b6622 100644
--- a/algo/kfp/KfpV0FinderConfig.h
+++ b/algo/kfp/KfpV0FinderConfig.h
@@ -69,6 +69,7 @@ namespace cbm::algo::kfp
/// \brief Configuration for the V0 finder
struct V0FinderConfig {
Cuts cuts; ///< Different selection cuts
+ uint32_t bmonAddress; ///< Address of BMON diamond (if multiple alternative are present, only one must be selected)
double tZeroOffset; ///< Offset for T0 [ns]
double qpAssignedUncertainty; ///< Assigned relative uncertainty for q/p estimation
int primaryAssignedPdg; ///< Assigned PDG hypothesis for primary particles
@@ -76,6 +77,7 @@ namespace cbm::algo::kfp
CBM_YAML_PROPERTIES(
yaml::Property(&V0FinderConfig::cuts, "cuts", "Different selection cuts"),
+ yaml::Property(&V0FinderConfig::bmonAddress, "bmon_address", "Address of reference BMON diamond"),
yaml::Property(&V0FinderConfig::tZeroOffset, "t0_offset", "The t0 offset [ns]"),
yaml::Property(&V0FinderConfig::qpAssignedUncertainty, "qa_uncertainty", "Assigned relative uncertainty for q/p"),
yaml::Property(&V0FinderConfig::primaryAssignedPdg, "primary_pdg", "Assigned PDG code for primary tracks"),
diff --git a/algo/kfp/KfpV0FinderMonitor.h b/algo/kfp/KfpV0FinderMonitor.h
new file mode 100644
index 0000000000000000000000000000000000000000..d60d5b135478e2e6d379a26b69e8755d9b2bd96b
--- /dev/null
+++ b/algo/kfp/KfpV0FinderMonitor.h
@@ -0,0 +1,108 @@
+/* Copyright (C) 2025 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file KfpV0FinderMonitor.h
+/// \date 03.02.2025
+/// \brief A monitor for the V0Finder
+/// \author Sergei Zharko <s.zharko@gsi.de>
+
+#pragma once
+
+#include "CaMonitor.h"
+
+namespace cbm::algo::kfp
+{
+ /// \enum ECounter
+ /// \brief Counter keys for the V0FinderMonitor
+ enum class ECounter
+ {
+ TracksTotal, ///< Total number of tracks
+ TracksSelected, ///< Tracks, which satisfy topology PID applicability
+ TracksInfiniteParam, ///< Tracks, which have infinite parameters
+ TracksWoTofHits, ///< Tracks, which have no TOF hits
+ TracksWNegativeTofHitTime, ///< Tracks, the last TOF hit of which has a negative time (it's time is less then the t0)
+ TracksWoStsHits, ///< Tracks, which have no STS hits
+ 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
+ EventsWoTzero, ///< Number of events with undefined t-zero
+ EventsLambdaCand, ///< Events with at least one pion and one proton candidate
+ KfpEventsLambdaCand, ///< Events with lambda-candidates in KF-particle
+ KfpLambdaCandidates, ///< Number of lambda-candidates
+ END
+ };
+
+ /// \enum ETimer
+ /// \brief Timer keys for the V0FinderMonitor
+ /* clang-format off */
+ enum class ETimer
+ {
+ 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>;
+
+ /// \class Monitor
+ /// \brief A monitor for the V0Finder
+ class V0FinderMonitor : public ca::Monitor<ECounter, ETimer> {
+ public:
+ /// \brief Default constructor
+ V0FinderMonitor() : ca::Monitor<ECounter, ETimer>("V0 finder monitor")
+ {
+ SetCounterName(ECounter::TracksTotal, "all tracks");
+ SetCounterName(ECounter::TracksSelected, "pre-selected tracks");
+ SetCounterName(ECounter::TracksInfiniteParam, "tracks satisfying PID selection");
+ SetCounterName(ECounter::TracksWoTofHits, "tracks w/o TOF hits");
+ SetCounterName(ECounter::TracksWNegativeTofHitTime, "tracks w/ negative TOF time");
+ SetCounterName(ECounter::TracksWoStsHits, "tracks w/o STS hits");
+ 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");
+ SetCounterName(ECounter::EventsWoTzero, "events w/o t0");
+ SetCounterName(ECounter::EventsLambdaCand, "events passed to KFP");
+ SetCounterName(ECounter::KfpEventsLambdaCand, "events w/ lambda candidates");
+ SetCounterName(ECounter::KfpLambdaCandidates, "lambda candidates");
+
+ 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");
+
+ SetRatioKeys({ECounter::EventsTotal, ECounter::TracksTotal, ECounter::TracksSelected});
+ }
+
+ private:
+ friend class boost::serialization::access;
+ template<typename Archive>
+ void serialize(Archive& ar, const unsigned int /*version*/)
+ {
+ ar& boost::serialization::base_object<ca::Monitor<ECounter, ETimer>>(*this);
+ }
+ };
+
+} // namespace cbm::algo::kfp
diff --git a/algo/kfp/KfpV0FinderQa.cxx b/algo/kfp/KfpV0FinderQa.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..ffc903e2a2dcc2c12efa06b4bda11db6db29ffb9
--- /dev/null
+++ b/algo/kfp/KfpV0FinderQa.cxx
@@ -0,0 +1,100 @@
+/* Copyright (C) 2025 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file KfpV0FinderQa.cxx
+/// \date 13.02.2025
+/// \brief A V0 finding algorithm QA (implementation)
+/// \author Sergei Zharko <s.zharko@gsi.de>
+
+#include "kfp/KfpV0FinderQa.h"
+
+#include "global/RecoResults.h"
+#include "kfp/KfpV0Finder.h"
+
+
+using cbm::algo::kfp::V0FinderQa;
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void V0FinderQa::Init()
+{
+ using qa::CanvasConfig;
+ using qa::H1D;
+ using qa::PadConfig;
+
+ //* Histogram initialisation
+ fvphMassLambdaCand =
+ MakeObj<H1D>("kfp_mass_lambda", "Mass of #Lambda-candidates;m [GeV/c^{2}];Counts", kMassB, kMassL, kMassU);
+ fvphMassAll = MakeObj<H1D>("kfp_mass_all", "Mass of particles;m [GeV/c^{2}];Counts", 300, 0., 1.5);
+ fvphDcaAll = MakeObj<H1D>("kfp_dca_all", "DCA of tracks to origin;DCA [cm];Counts", kDcaB, kDcaL, kDcaU);
+ fvphBetaAll = MakeObj<H1D>("kfp_beta_all", "Speed of tracks;#beta;Counts", kBetaB, kBetaL, kBetaU);
+ fvphBetaPion = MakeObj<H1D>("kfp_beta_pion", "Speed of #pi-candidates;#beta;Counts", kBetaB, kBetaL, kBetaU);
+ fvphBetaProton = MakeObj<H1D>("kfp_beta_proton", "Speed of proton-candidates;#beta;Counts", kBetaB, kBetaL, kBetaU);
+ fvphMomAll = MakeObj<H1D>("kfp_mom_all", "Momentum of tracks;p [GeV/c];Counts", kBetaB, kBetaL, kBetaU);
+ fvphMomPion = MakeObj<H1D>("kfp_mom_pion", "Momentum of #pi-candidates;p [GeV/c];Counts", kBetaB, kBetaL, kBetaU);
+ fvphMomProton =
+ MakeObj<H1D>("kfp_mom_proton", "Momentum of proton-candidates;p [GeV/c];Counts", kBetaB, kBetaL, kBetaU);
+
+ //* Canvas initialisation
+ auto canv = CanvasConfig("kfp_lambda", "Lambda-trigger summary QA", 3, 2);
+ canv.AddPadConfig(PadConfig(fvphMassLambdaCand, "hist")); // (0, 0)
+ canv.AddPadConfig(PadConfig(fvphDcaAll, "hist")); // (0, 1)
+ canv.AddPadConfig(PadConfig(fvphBetaAll, "hist"));
+ canv.AddPadConfig(PadConfig(fvphBetaPion, "hist"));
+ canv.AddPadConfig(PadConfig(fvphBetaProton, "hist"));
+ AddCanvasConfig(canv);
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void V0FinderQa::Exec(const RecoResults& recoEvent, const V0Finder& v0Finder)
+{
+ //* Fill track distributions
+ const auto& tracks{recoEvent.tracks};
+ if (v0Finder.GetT0s().size() == 1) {
+ for (uint32_t iTrk = 0; iTrk < tracks.size(); ++iTrk) {
+ const auto& particleInfo{v0Finder.GetParticleInfo()[iTrk]};
+ const auto& trkParFst{(v0Finder.GetTrackAssignedParams()[iTrk]).first};
+ bool bPdgDefined = (particleInfo.fPdg != V0Finder::kUndefPdg);
+ fvphDcaAll->Fill(bPdgDefined ? particleInfo.fDca : -999);
+ fvphBetaAll->Fill(particleInfo.fBeta);
+ if (bPdgDefined) {
+ if (particleInfo.fPdg == -211) {
+ fvphBetaPion->Fill(particleInfo.fBeta);
+ fvphMomPion->Fill(trkParFst.GetP());
+ }
+ else if (particleInfo.fPdg == 2212) {
+ fvphBetaProton->Fill(particleInfo.fBeta);
+ fvphMomProton->Fill(trkParFst.GetP());
+ }
+ }
+ }
+ }
+ else {
+ for (uint32_t iTrk = 0; iTrk < tracks.size(); ++iTrk) {
+ const auto& particleInfo{v0Finder.GetParticleInfo()[iTrk]};
+ bool bPdgDefined = (particleInfo.fPdg != V0Finder::kUndefPdg);
+ fvphDcaAll->Fill(bPdgDefined ? particleInfo.fDca : -999);
+ fvphBetaAll->Fill(-999);
+ if (bPdgDefined) {
+ if (particleInfo.fPdg == -211) {
+ fvphBetaPion->Fill(-999);
+ fvphMomPion->Fill(-999);
+ }
+ else if (particleInfo.fPdg == 2212) {
+ fvphBetaProton->Fill(-999);
+ }
+ }
+ }
+ }
+
+ //* Fill particle distributions
+ const auto& particles = v0Finder.GetTopoReconstructor()->GetParticles();
+ for (const auto& particle : particles) {
+ fvphMassAll->Fill(particle.GetMass());
+ if (particle.GetPDG() == 3122) {
+ fvphMassLambdaCand->Fill(particle.GetMass());
+ }
+ }
+}
diff --git a/algo/kfp/KfpV0FinderQa.h b/algo/kfp/KfpV0FinderQa.h
new file mode 100644
index 0000000000000000000000000000000000000000..35b914e182bf5e1890a3ce2bd0f10b6598f40e63
--- /dev/null
+++ b/algo/kfp/KfpV0FinderQa.h
@@ -0,0 +1,84 @@
+/* Copyright (C) 2025 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file KfpV0FinderQa.h
+/// \date 13.02.2025
+/// \brief A V0 finding algorithm QA
+/// \author Sergei Zharko <s.zharko@gsi.de>
+
+#pragma once
+
+#include "qa/QaTaskHeader.h"
+
+namespace cbm::algo
+{
+ class RecoResults;
+
+ namespace qa
+ {
+ class H1D;
+ class H2D;
+ } // namespace qa
+
+ namespace kfp
+ {
+ class V0Finder;
+ }
+} // namespace cbm::algo
+
+namespace cbm::algo::kfp
+{
+ /// \class V0FinderQa
+ /// \brief A QA-task for the V0-finding algorithm
+ class V0FinderQa : public qa::TaskHeader {
+ public:
+ /// \brief Constructor
+ /// \param pManager Pointer to the QA manager
+ /// \param name Name of the QA
+ V0FinderQa(const std::unique_ptr<qa::Manager>& pManager, std::string_view name) : qa::TaskHeader(pManager, name) {}
+
+ /// \brief Copy constructor
+ V0FinderQa(const V0FinderQa&) = delete;
+
+ /// \brief Move constructor
+ V0FinderQa(V0FinderQa&&) = delete;
+
+ /// \brief Copy assignment operator
+ V0FinderQa& operator=(const V0FinderQa&) = delete;
+
+ /// \brief Move assignment operator
+ V0FinderQa& operator=(V0FinderQa&&) = delete;
+
+ /// \brief Executes the task, fills the histograms
+ /// \param recoEvent A reconstructed event instance
+ /// \param v0Finder A V0-finder instance
+ void Exec(const RecoResults& recoEvent, const V0Finder& v0Finder);
+
+ /// \brief Initialized the task
+ void Init();
+
+ private:
+ //* Constants
+ static constexpr int kMassB = 200; ///< Lambda-candidate mass: number of bins
+ static constexpr double kMassL = 1.08; ///< Lambda-candidate mass: lower bound [GeV/c2]
+ static constexpr double kMassU = 1.18; ///< Lambda-candidate mass: upper bound [GeV/c2]
+ static constexpr int kDcaB = 240; ///< DCA to origin: number of bins
+ static constexpr double kDcaL = 0.; ///< DCA to origin: lower bound [cm]
+ static constexpr double kDcaU = 12.; ///< DCA to origin: upper bound [cm]
+ static constexpr int kBetaB = 240; ///< Speed of particle: number of bins
+ static constexpr double kBetaL = 0.; ///< Speed of particle: lower bound [c]
+ static constexpr double kBetaU = 1.2; ///< Speed of particle: upper bound [c]
+
+ //* Histograms
+ qa::H1D* fvphMassLambdaCand{nullptr}; ///< Mass of Lambda-candidates
+ qa::H1D* fvphMassAll{nullptr}; ///< Mass of all particles in the topology
+ qa::H1D* fvphDcaAll{nullptr}; ///< DCA of particles to origin
+ qa::H1D* fvphBetaAll{nullptr}; ///< Speed of all particles
+ qa::H1D* fvphBetaPion{nullptr}; ///< Speed of pion-candidates
+ qa::H1D* fvphBetaProton{nullptr}; ///< Speed of proton-candidates
+ qa::H1D* fvphMomAll{nullptr}; ///< Speed of all particles
+ qa::H1D* fvphMomPion{nullptr}; ///< Speed of pion-candidates
+ qa::H1D* fvphMomProton{nullptr}; ///< Speed of proton-candidates
+ };
+} // namespace cbm::algo::kfp
diff --git a/algo/kfp/interface/CMakeLists.txt b/algo/kfp/interface/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..4beacabb176c8fdd9f21d310a7ffde164c2c6bd5
--- /dev/null
+++ b/algo/kfp/interface/CMakeLists.txt
@@ -0,0 +1,38 @@
+if(NOT CBM_ONLINE_STANDALONE)
+ ### CbmKFParticleOnlineInterface
+ # Just reusing the KFParticle library, compiled in the external
+ add_library(CbmKFParticleOnlineInterface INTERFACE)
+ target_include_directories(CbmKFParticleOnlineInterface INTERFACE)
+ target_compile_definitions(CbmKFParticleOnlineInterface
+ INTERFACE DO_TPCCATRACKER_EFF_PERFORMANCE NonhomogeneousField CBM USE_TIMERS)
+ target_link_libraries(CbmKFParticleOnlineInterface INTERFACE KFParticle)
+ install(TARGETS CbmKFParticleOnlineInterface DESTINATION lib)
+else()
+ # Creating a replacement of the KFParticle library for a standalone mode
+ set(KFP_SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/../../../external/KFParticle/KFParticle)
+
+ set(SRCS
+ ${KFP_SOURCE_DIR}/KFParticle.cxx
+ ${KFP_SOURCE_DIR}/KFParticleBase.cxx
+ ${KFP_SOURCE_DIR}/KFParticleBaseSIMD.cxx
+ ${KFP_SOURCE_DIR}/KFParticleDatabase.cxx
+ ${KFP_SOURCE_DIR}/KFParticleFinder.cxx
+ ${KFP_SOURCE_DIR}/KFParticlePVReconstructor.cxx
+ ${KFP_SOURCE_DIR}/KFParticleSIMD.cxx
+ ${KFP_SOURCE_DIR}/KFParticleTopoReconstructor.cxx
+ ${KFP_SOURCE_DIR}/KFPEmcCluster.cxx
+ ${KFP_SOURCE_DIR}/KFPTrack.cxx
+ ${KFP_SOURCE_DIR}/KFPTrackVector.cxx
+ ${KFP_SOURCE_DIR}/KFPVertex.cxx
+ ${KFP_SOURCE_DIR}/KFVertex.cxx
+ )
+
+ add_library(CbmKFParticleOnlineInterface SHARED ${SRCS})
+ target_include_directories(CbmKFParticleOnlineInterface PUBLIC ${KFP_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR})
+ target_compile_definitions(CbmKFParticleOnlineInterface PUBLIC NonhomogeneousField CBM CBM_ONLINE)
+ target_link_libraries(CbmKFParticleOnlineInterface
+ PUBLIC Vc::Vc
+ ROOT::Core
+ )
+ install(TARGETS CbmKFParticleOnlineInterface DESTINATION lib)
+endif()
diff --git a/algo/kfp/interface/RootTypesDef.h b/algo/kfp/interface/RootTypesDef.h
new file mode 100644
index 0000000000000000000000000000000000000000..e6eaa5500553524aaf4192d2a17108cd81b30da8
--- /dev/null
+++ b/algo/kfp/interface/RootTypesDef.h
@@ -0,0 +1,18 @@
+/* Copyright (C) 2025 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file RootTypesDef.h
+/// \date 11.02.2025
+/// \brief A compatibility header for the KFParticle code
+/// \author Sergei Zharko <s.zharko@gsi.de>
+
+#pragma once
+
+#if __has_include(<RtypesCore.h>)
+#include <RtypesCore.h>
+#else
+using Bool_t = bool;
+using Int_t = int;
+using Float_t = float;
+#endif
diff --git a/algo/qa/QaData.cxx b/algo/qa/QaData.cxx
index 520ed4e522e022e742179e2ecb4c5f564b95da4d..6c42b3a5237df1cad00eacd1d9b9a579cf661f6d 100644
--- a/algo/qa/QaData.cxx
+++ b/algo/qa/QaData.cxx
@@ -57,6 +57,7 @@ try {
<< " with inconsistent flags (see HistogramMetadata::CheckFlags for detailes)";
throw std::runtime_error(msg.str());
}
+ L_(info) << " - task: " << task.fsName << ", histogram: " << h.GetName();
vHistCfgs.emplace_back(h.GetName() + "!" + h.GetMetadataString(), task.fsName);
};
fsTaskNames += fmt::format("{} ", task.fsName);
diff --git a/algo/qa/hitfind/BmonHitfindQa.h b/algo/qa/hitfind/BmonHitfindQa.h
index 020f62fa821a2c8ec298c8aa758684aee023cc21..b2a8b3e30e709d099ef1059aaa5a619864d02c54 100644
--- a/algo/qa/hitfind/BmonHitfindQa.h
+++ b/algo/qa/hitfind/BmonHitfindQa.h
@@ -34,8 +34,6 @@ namespace cbm::algo::bmon
{
/// \class HitfindQa
/// \brief A QA module for the BMON hit-finder
- /// \param pManager Pointer to the QA manager
- /// \param name Name of the QA (directory)
class HitfindQa : public qa::TaskHeader {
public:
/// \brief Constructor
@@ -44,7 +42,6 @@ namespace cbm::algo::bmon
HitfindQa(const std::unique_ptr<qa::Manager>& pManager, std::string_view name) : qa::TaskHeader(pManager, name) {}
/// \brief Constructor from the configuration object
- /// \param config QA configuration object
HitfindQa() = default;
/// \brief Copy constructor
diff --git a/core/data/CbmEventTriggers.h b/core/data/CbmEventTriggers.h
index ce1b8de535825650cfb3e6cbf77661bbf44c36c5..a2c252df56bab630e79c7e776830b70e0fbbcd99 100644
--- a/core/data/CbmEventTriggers.h
+++ b/core/data/CbmEventTriggers.h
@@ -9,28 +9,37 @@
#pragma once
+#include <cstdint>
+#include <string>
+
#if !defined(NO_ROOT) && !XPU_IS_HIP_CUDA
#include <Rtypes.h> // for ClassDef
#endif
+#include <boost/serialization/access.hpp>
+
/// \class CbmEventTriggers
/// \brief Class to store different triggers for a given event
class CbmEventTriggers {
public:
/// \enum ETrigger
- /// \brief Defines a trigger address
- enum class ETrigger : uint8_t
+ /// \brief Defines a trigger bitmask
+ enum class ETrigger : uint32_t
{
- Lambda = 0b00000001, ///< Lambda-trigger
- Ks = 0b00000010 ///< Ks-trigger
+ Lambda = 0x00000001, ///< Lambda-trigger
+ Ks = 0x00000002 ///< Ks-trigger
};
//using Trigger_t = std::underlying_type_t<ETrigger>;
- using Trigger_t = uint8_t;
+ using Trigger_t = uint32_t;
/// \brief Default constructor
CbmEventTriggers() = default;
+ /// \brief A constructor from integer
+ /// \param bitmap A bitmap of the triggers
+ explicit CbmEventTriggers(uint32_t bitmap) { fBitmap = bitmap; }
+
/// \brief Copy constructor
CbmEventTriggers(const CbmEventTriggers&) = default;
@@ -46,33 +55,47 @@ class CbmEventTriggers {
/// \brief Move assignment operator
CbmEventTriggers& operator=(CbmEventTriggers&&) = default;
+ /// \brief Gets a bitmap
+ /// \return bitmap (integer)
+ Trigger_t GetBitmap() const { return fBitmap; }
+
/// \brief Sets a trigger
/// \param key Trigger key
- void Set(ETrigger key) { fTriggers |= static_cast<Trigger_t>(key); }
+ void Set(ETrigger key) { fBitmap |= static_cast<Trigger_t>(key); }
/// \brief Resets a trigger
/// \param key Trigger key
- void Reset(ETrigger key) { fTriggers &= ~static_cast<Trigger_t>(key); }
+ void Reset(ETrigger key) { fBitmap &= ~static_cast<Trigger_t>(key); }
+
+ /// \brief Resets all the triggers
+ void ResetAll() { fBitmap = 0; }
/// \brief Tests a particular single trigger
/// \param key Trigger key
- bool Test(ETrigger key) const { return static_cast<bool>(fTriggers & static_cast<Trigger_t>(key)); }
+ bool Test(ETrigger key) const { return static_cast<bool>(fBitmap & static_cast<Trigger_t>(key)); }
- /// \brief Tests, if ALL the triggers in the mask are on
- /// \param mask Trigger mask
- bool TestAll(Trigger_t mask) const { return mask == (fTriggers & mask); }
+ /// \brief Tests, if ALL the triggers in the bitmask are on
+ /// \param bitmask Trigger bitmask
+ bool TestAll(Trigger_t bitmask) const { return bitmask == (fBitmap & bitmask); }
- /// \brief Tests, if ANY of the triggers in the mask are on
- /// \param mask Trigger mask
- bool TestAny(Trigger_t mask) const { return static_cast<bool>(fTriggers | mask); }
+ /// \brief Tests, if ANY of the triggers in the bitmask are on
+ /// \param bitmask Trigger bitmask
+ bool TestAny(Trigger_t bitmask) const { return static_cast<bool>(fBitmap | bitmask); }
/// \brief String representation of the class content
std::string ToString() const;
private:
- Trigger_t fTriggers;
+ Trigger_t fBitmap{0}; ///< bitmap storing the triggers according to ETrigger
+
+ friend class boost::serialization::access;
+ template<class Archive>
+ void serialize(Archive& ar, const unsigned int /*version*/)
+ {
+ ar& fBitmap;
+ }
#if !defined(NO_ROOT) && !XPU_IS_HIP_CUDA
- ClassDefNV(CbmEventTriggers, 1);
+ ClassDefNV(CbmEventTriggers, 2);
#endif
};
diff --git a/core/data/global/CbmDigiEvent.h b/core/data/global/CbmDigiEvent.h
index 242740e57133ed58695947df108315761048f528..ec90985bbd4561a9c2a9ccfdef496caaa32b7ee7 100644
--- a/core/data/global/CbmDigiEvent.h
+++ b/core/data/global/CbmDigiEvent.h
@@ -1,4 +1,4 @@
-/* Copyright (C) 2021-22 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+/* Copyright (C) 2021-25 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
Authors: Volker Friese [committer] */
@@ -6,6 +6,7 @@
#define CBMDIGIEVENT_H 1
#include "CbmDigiData.h"
+#include "CbmEventTriggers.h"
#include <boost/serialization/access.hpp>
@@ -14,29 +15,32 @@
** @brief Collection of digis from all detector systems within one event
** @author Volker Friese <v.friese@gsi.de>
** @since 7.12.2022
- ** @version 2.0
+ ** @version 3.0
**/
class CbmDigiEvent {
public:
- CbmDigiData fData; ///< Event data
- uint64_t fNumber; ///< Event identifier
- double fTime; ///< Event trigger time [ns]
-
-
- friend class boost::serialization::access;
- /** @brief BOOST serializer**/
- template<class Archive>
- void serialize(Archive& ar, const unsigned int /*version*/)
- {
- ar& fData;
- ar& fNumber;
- ar& fTime;
- }
+ CbmDigiData fData; ///< Event data
+ uint64_t fNumber; ///< Event identifier
+ double fTime; ///< Event trigger time [ns]
+ CbmEventTriggers fSelectionTriggers; ///< Event selection triggers
+
+ friend class boost::serialization::access;
+ /** @brief BOOST serializer**/
+ template<class Archive>
+ void serialize(Archive& ar, const unsigned int version)
+ {
+ ar& fData;
+ ar& fNumber;
+ ar& fTime;
+ if (version >= 3) {
+ ar& fSelectionTriggers;
+ }
+ }
// --- ROOT serializer
#ifndef NO_ROOT
- ClassDefNV(CbmDigiEvent, 2);
+ ClassDefNV(CbmDigiEvent, 3);
#endif
/** @brief Clear content **/
@@ -45,6 +49,7 @@ public:
fData.Clear();
fNumber = 0;
fTime = 0.;
+ fSelectionTriggers.ResetAll();
}
};
diff --git a/external/CMakeLists.txt b/external/CMakeLists.txt
index d0843683605e1ffeec46e52a02158a289405d0ab..9328a585c1d7b5ed3640fea9f1126a30a34820bd 100644
--- a/external/CMakeLists.txt
+++ b/external/CMakeLists.txt
@@ -70,6 +70,7 @@ if(DOWNLOAD_EXTERNALS)
Include(InstallYamlCpp.cmake)
Include(InstallPoolSTL.cmake)
+
if (NOT CBM_ONLINE_STANDALONE) # Not required for online standalone
Include(InstallKFParticle.cmake)
@@ -86,8 +87,22 @@ if(DOWNLOAD_EXTERNALS)
Include(InstallQa.cmake)
endif()
+ else()
+
+ # Download KFParticle source from the repository, but do not install it
+ download_project_if_needed(PROJECT kfparticle_source
+ GIT_REPOSITORY "https://github.com/szharko/KFParticle.git"
+ GIT_TAG "d088e3e019ac588a0528aa560ceaca247580d881" # CBM_ONLINE directive to turn off dependencies on ROOT
+ SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/KFParticle
+ CONFIGURE_COMMAND ""
+ BUILD_COMMAND ""
+ INSTALL_COMMAND ""
+ )
+
endif()
+
+
else()
# Define targets which are needed by CbmRoot but are not available
# whithout the external packages
diff --git a/external/InstallParameter.cmake b/external/InstallParameter.cmake
index 77e3f424dc53079799a63e8be1dd5d1b94720f0d..7d02e04add0d89dfae88aea1d0ef2465fa3d1ee3 100644
--- a/external/InstallParameter.cmake
+++ b/external/InstallParameter.cmake
@@ -1,5 +1,7 @@
set(PARAMETER_VERSION cf6a880a56a8695b485efe358231270e9179cb82) # 2025/02/15
set(PARAMETER_SRC_URL "https://git.cbm.gsi.de/CbmSoft/cbmroot_parameter.git")
+#set(PARAMETER_VERSION 96b2189ac86b3d49c32f63da76e68f66544bba55) # 2025/05/02, BMON hitfinder parameters for online
+#set(PARAMETER_SRC_URL "https://git.cbm.gsi.de/s.zharko/cbmroot_parameter.git")
download_project_if_needed(PROJECT Parameter_source
GIT_REPOSITORY ${PARAMETER_SRC_URL}
diff --git a/macro/KF/configs/mcbm_kfpf_lambda.yaml b/macro/KF/configs/mcbm_kfpf_lambda.yaml
index 01e77b1855b37c5f3f763f6f1175faf7ff89af12..fc5599d3723479b6fc379a60f56e3d0d5983131c 100644
--- a/macro/KF/configs/mcbm_kfpf_lambda.yaml
+++ b/macro/KF/configs/mcbm_kfpf_lambda.yaml
@@ -7,11 +7,11 @@
# \brief Standard configuration for the lambda finding in mCBM
# \author Sergei Zharko <s.zharko@gsi.de>
-reconstruct_pdg: 3122 # Lambda -> pi- p
-
-t0_offset: 0.12 # Time offset for T0 [ns]
-qa_uncertainty: 0.1 # Assigned relative uncertainty for the q/p estimation
-primary_pdg: 321 # Assigned PDG code for primary tracks
+reconstruct_pdg: 3122 # Lambda -> pi- p
+bmon_address: 0x00202806 # Address of reference BMON diamond
+t0_offset: 0.12 # Time offset for T0 [ns]
+qa_uncertainty: 0.1 # Assigned relative uncertainty for the q/p estimation
+primary_pdg: 321 # Assigned PDG code for primary tracks
cuts: # Different cuts
kfp: # KFParticleFinder specific cuts
diff --git a/services/histserv/app/Application.cxx b/services/histserv/app/Application.cxx
index 5eb5362f051fc3f34fc1108e50bd17e1adf6e4f2..a4978c303efe8ae9925b2a65f017b8303c796a37 100644
--- a/services/histserv/app/Application.cxx
+++ b/services/histserv/app/Application.cxx
@@ -800,17 +800,21 @@ bool Application::PrepareCanvas(uint32_t uCanvIdx)
if ("nullptr" != sName) {
TObject* pObj = fArrayHisto[FindHistogram(sName)];
- if (nullptr != dynamic_cast<TProfile2D*>(pObj)) {
- dynamic_cast<TProfile2D*>(pObj)->Draw(conf.GetOption(uPadIdx, uObjIdx).data());
+ if (auto* pHistP2 = dynamic_cast<TProfile2D*>(pObj)) {
+ pHistP2->Draw(conf.GetOption(uPadIdx, uObjIdx).data());
} // if( nullptr != dynamic_cast< TProfile *>( pObj ) )
- else if (nullptr != dynamic_cast<TProfile*>(pObj)) {
- dynamic_cast<TProfile*>(pObj)->Draw(conf.GetOption(uPadIdx, uObjIdx).data());
+ else if (auto* pHistP1 = dynamic_cast<TProfile*>(pObj)) {
+ pHistP1->SetLineColor(uObjIdx + 1);
+ pHistP1->SetMarkerColor(uObjIdx + 1);
+ pHistP1->Draw(conf.GetOption(uPadIdx, uObjIdx).data());
} // if( nullptr != dynamic_cast< TProfile *>( pObj ) )
- else if (nullptr != dynamic_cast<TH2*>(pObj)) {
- dynamic_cast<TH2*>(pObj)->Draw(conf.GetOption(uPadIdx, uObjIdx).data());
+ else if (auto* pHistH2 = dynamic_cast<TH2*>(pObj)) {
+ pHistH2->Draw(conf.GetOption(uPadIdx, uObjIdx).data());
} // if( nullptr != dynamic_cast< TH2 *>( pObj ) )
- else if (nullptr != dynamic_cast<TH1*>(pObj)) {
- dynamic_cast<TH1*>(pObj)->Draw(conf.GetOption(uPadIdx, uObjIdx).data());
+ else if (auto* pHistH1 = dynamic_cast<TH1*>(pObj)) {
+ pHistH1->SetLineColor(uObjIdx + 1);
+ pHistH1->SetMarkerColor(uObjIdx + 1);
+ pHistH1->Draw(conf.GetOption(uPadIdx, uObjIdx).data());
} // if( nullptr != dynamic_cast< TH1 *>( pObj ) )
else
LOG(warning) << " Unsupported object type for " << sName << " when preparing canvas " << conf.GetName();