diff --git a/algo/CMakeLists.txt b/algo/CMakeLists.txt
index f7e03a8eede7ebc43055a16f1219ca2a113096fc..ee9040be2fb508c78b0fece2e4c865d0a098a47f 100644
--- a/algo/CMakeLists.txt
+++ b/algo/CMakeLists.txt
@@ -93,6 +93,7 @@ set(SRCS
detectors/tof/ReadoutConfig.cxx
detectors/tof/Unpack.cxx
detectors/tof/Hitfind.cxx
+ detectors/tof/HitfinderChain.cxx
detectors/bmon/ReadoutConfig.cxx
detectors/bmon/Unpack.cxx
detectors/trd/ReadoutConfig.cxx
diff --git a/algo/base/PartitionedSpan.h b/algo/base/PartitionedSpan.h
index 988cc6f2daba2dedbb37730f7951cffdd23af62f..e004f92b485f91e09092e3b9aebc439d43413005 100644
--- a/algo/base/PartitionedSpan.h
+++ b/algo/base/PartitionedSpan.h
@@ -101,6 +101,14 @@ namespace cbm::algo
EnsureDimensions();
}
+ template<typename U, typename = detail::EnableOnConst<U, T>>
+ PartitionedSpan(PartitionedSpan<U> other)
+ : fData(other.Data())
+ , fOffsets(other.Offsets())
+ , fAdresses(other.Addresses())
+ {
+ }
+
gsl::span<T> operator[](size_t i) const
{
EnsureBounds(i);
diff --git a/algo/base/config/Yaml.h b/algo/base/config/Yaml.h
index 50e6684eb05b54ad3be7f2a48708a88f05747fc6..dbb60801885e1110366278a343cd1ab9711042a2 100644
--- a/algo/base/config/Yaml.h
+++ b/algo/base/config/Yaml.h
@@ -13,11 +13,15 @@
#include "BaseTypes.h"
#include "Definitions.h"
#include "Property.h"
+#include "compat/Filesystem.h"
#include "util/EnumDict.h"
namespace cbm::algo::config
{
+ template<typename T>
+ T Read(const YAML::Node& node);
+
template<typename T, T... Values, typename Func>
constexpr void ForEach(std::integer_sequence<T, Values...>, Func&& func)
{
@@ -34,6 +38,13 @@ namespace cbm::algo::config
static constexpr std::optional<YAML::EMITTER_MANIP> value = T::FormatAs;
};
+ template<typename T>
+ T ReadFromFile(fs::path path)
+ {
+ YAML::Node node = YAML::LoadFile(path.string());
+ return Read<T>(node);
+ }
+
template<typename T>
T Read(const YAML::Node& node)
{
diff --git a/algo/ca/TrackingChain.cxx b/algo/ca/TrackingChain.cxx
index 05d1f1d0cd909eb948ad3d8157913964c03d2bcd..39f3b7a7af2a235970eb3307d08e7eeac7488d80 100644
--- a/algo/ca/TrackingChain.cxx
+++ b/algo/ca/TrackingChain.cxx
@@ -16,7 +16,8 @@
#include "CaInitManager.h"
#include "CaParameters.h"
-using cbm::algo::PartitionedPODVector;
+using namespace cbm::algo;
+
using cbm::algo::TrackingChain;
using cbm::algo::ca::EDetectorID;
using cbm::algo::ca::Framework;
@@ -32,7 +33,7 @@ using cbm::algo::ca::constants::clrs::GNb; // grin bald text
void TrackingChain::Init()
{
// ------ Read parameters from binary
- std::string paramFileBase = "mcbm_beam_2022_05_23_nickel_w_time.ca.par"; // TODO: Get the setup name from Opts()
+ std::string paramFileBase = "mcbm_beam_2022_05_23_nickel.ca.par"; // TODO: Get the setup name from Opts()
auto paramFile = Opts().ParamsDir();
paramFile /= paramFileBase;
L_(info) << "Tracking Chain: reading CA parameters file " << GNb << paramFile.string() << CL << '\n';
@@ -48,7 +49,7 @@ void TrackingChain::Init()
// ---------------------------------------------------------------------------------------------------------------------
//
-TrackingChain::Return_t TrackingChain::Run(const RecoResults& recoResults)
+TrackingChain::Return_t TrackingChain::Run(Input_t recoResults)
{
// ----- Init input data ---------------------------------------------------------------------------------------------
this->PrepareInput(recoResults);
@@ -71,13 +72,14 @@ void TrackingChain::Finalize() {}
// ---------------------------------------------------------------------------------------------------------------------
//
-void TrackingChain::PrepareInput(const RecoResults& recoResults)
+void TrackingChain::PrepareInput(Input_t recoResults)
{
fNofHitKeys = 0;
- int nHitsTot = recoResults.stsHits.NElements();
+ int nHitsTot = recoResults.stsHits.NElements() + recoResults.tofHits.NElements();
L_(info) << "Tracking chain: input has " << nHitsTot << " hits";
fCaDataManager.ResetInputData(nHitsTot);
ReadHits<EDetectorID::Sts>(recoResults.stsHits);
+ ReadHits<EDetectorID::Tof>(recoResults.tofHits);
fCaDataManager.SetNhitKeys(fNofHitKeys);
L_(info) << "Tracking chain:" << fCaDataManager.GetNofHits() << " will be passed to the ca::Framework";
fCaFramework.ReceiveInputData(fCaDataManager.TakeInputData());
@@ -86,8 +88,10 @@ void TrackingChain::PrepareInput(const RecoResults& recoResults)
// ---------------------------------------------------------------------------------------------------------------------
//
template<EDetectorID DetID>
-void TrackingChain::ReadHits(const PartitionedPODVector<ca::HitTypes_t::at<DetID>>& hits)
+void TrackingChain::ReadHits(PartitionedSpan<const ca::HitTypes_t::at<DetID>> hits)
{
+ constexpr bool IsSTS = (DetID == EDetectorID::Sts);
+
ca::HitKeyIndex_t firstHitKey = fNofHitKeys;
int64_t dataStreamDet = static_cast<int64_t>(DetID) << 60; // detector part of the data stream
for (size_t iPartition = 0; iPartition < hits.NPartitions(); ++iPartition) {
@@ -97,7 +101,7 @@ void TrackingChain::ReadHits(const PartitionedPODVector<ca::HitTypes_t::at<DetID
int iStLocal = -1;
// FIXME: This definition of the station index works only for STS, and there is no any guaranty, that it will
// work for other mCBM setups.
- if constexpr (DetID == EDetectorID::Sts) { iStLocal = (extHitAddress >> 4) & 0xF; }
+ if constexpr (IsSTS) { iStLocal = (extHitAddress >> 4) & 0xF; }
int iStActive = (iStLocal != -1) ? fCaFramework.GetParameters().GetStationIndexActive(iStLocal, DetID) : -1;
size_t iOffset = hits.Offsets()[iPartition];
if (iStActive < 0) { continue; }
@@ -107,7 +111,7 @@ void TrackingChain::ReadHits(const PartitionedPODVector<ca::HitTypes_t::at<DetID
int iHitExt = iOffset + iPartHit;
// ---- Fill ca::Hit
ca::Hit caHit;
- if constexpr (DetID == EDetectorID::Sts) {
+ if constexpr (IsSTS) {
caHit.SetFrontKey(firstHitKey + hit.fFrontClusterId);
caHit.SetBackKey(firstHitKey + hit.fBackClusterId);
}
@@ -115,18 +119,18 @@ void TrackingChain::ReadHits(const PartitionedPODVector<ca::HitTypes_t::at<DetID
caHit.SetFrontKey(firstHitKey + iHitExt);
caHit.SetBackKey(caHit.FrontKey());
}
- caHit.SetX(hit.fX);
- caHit.SetY(hit.fY);
- caHit.SetZ(hit.fZ);
- caHit.SetT(hit.fTime);
- caHit.SetDx2(hit.fDx * hit.fDx);
- caHit.SetDy2(hit.fDy * hit.fDy);
- caHit.SetDxy(hit.fDxy);
- caHit.SetDt2(hit.fTimeError * hit.fTimeError);
+ caHit.SetX(hit.X());
+ caHit.SetY(hit.Y());
+ caHit.SetZ(hit.Z());
+ caHit.SetT(hit.Time());
+ caHit.SetDx2(hit.Dx() * hit.Dx());
+ caHit.SetDy2(hit.Dy() * hit.Dy());
+ if constexpr (IsSTS) caHit.SetDxy(hit.fDxy);
+ caHit.SetDt2(hit.TimeError() * hit.TimeError());
/// FIXME: Define ranges from the hit, when will be available
- caHit.SetRangeX(3.5 * hit.fDx);
- caHit.SetRangeY(3.5 * hit.fDy);
- caHit.SetRangeT(3.5 * hit.fTimeError);
+ caHit.SetRangeX(3.5 * hit.Dx());
+ caHit.SetRangeY(3.5 * hit.Dy());
+ caHit.SetRangeT(3.5 * hit.TimeError());
caHit.SetStation(iStActive);
caHit.SetId(fCaDataManager.GetNofHits());
fCaDataManager.PushBackHit(caHit, dataStream);
@@ -137,4 +141,4 @@ void TrackingChain::ReadHits(const PartitionedPODVector<ca::HitTypes_t::at<DetID
} // iPartition
}
-template void TrackingChain::ReadHits<EDetectorID::Sts>(const PartitionedPODVector<HitTypes_t::at<EDetectorID::Sts>>&);
+// template void TrackingChain::ReadHits<EDetectorID::Sts>(const PartitionedPODVector<HitTypes_t::at<EDetectorID::Sts>>&);
diff --git a/algo/ca/TrackingChain.h b/algo/ca/TrackingChain.h
index c59f81f400304af2d52f2a896dc51c53ac2dbb29..e5dcaccec3154cbba13c89b2a8c760e5bb41e340 100644
--- a/algo/ca/TrackingChain.h
+++ b/algo/ca/TrackingChain.h
@@ -10,6 +10,7 @@
#pragma once
#include "TrackingDefs.h"
+#include "tof/Hit.h"
#include <vector>
@@ -18,10 +19,10 @@
#include "CaTrack.h"
#include "CaTrackingMonitor.h"
#include "CaVector.h"
-#include "PartitionedVector.h"
+#include "PartitionedSpan.h"
#include "RecoResults.h"
#include "SubChain.h"
-#include "data/sts/Hit.h"
+#include "sts/Hit.h"
namespace cbm::algo
{
@@ -31,6 +32,11 @@ namespace cbm::algo
/// The class executes a tracking algorithm in the online data reconstruction chain.
class TrackingChain : public SubChain {
public:
+ struct Input_t {
+ PartitionedSpan<sts::Hit> stsHits;
+ PartitionedSpan<tof::Hit> tofHits;
+ };
+
using Return_t = std::pair<ca::Vector<ca::Track>, ca::TrackingMonitorData>;
/// \brief Provides action in the initialization of the run
@@ -39,7 +45,7 @@ namespace cbm::algo
/// \brief Provides action for a given time-slice
/// \param recoResults Structure of reconstruction results
/// \return A pair (vector of tracks, tracking monitor)
- Return_t Run(const RecoResults& recoResults);
+ Return_t Run(Input_t recoResults);
/// \brief Provides action in the end of the run
void Finalize();
@@ -47,13 +53,13 @@ namespace cbm::algo
private:
/// \brief Prepares input data
/// \param recoResults Structure of reconstruction results
- void PrepareInput(const RecoResults& recoResults);
+ void PrepareInput(Input_t recoResults);
/// \brief Reads from different detector subsystems
/// \tparam DetID Detector ID
/// \param hits Hits vector
template<ca::EDetectorID DetID>
- void ReadHits(const PartitionedPODVector<ca::HitTypes_t::at<DetID>>& hits);
+ void ReadHits(PartitionedSpan<const ca::HitTypes_t::at<DetID>> hits);
ca::Framework fCaFramework {}; ///< CA framework instance
ca::DataManager fCaDataManager {}; ///< CA data manager
diff --git a/algo/data/sts/Hit.h b/algo/data/sts/Hit.h
index 0a33355535b615cca69a1f6bee867468295b3bd5..091ba8db3c1688a28b6f76462a18c1a6cd329f4c 100644
--- a/algo/data/sts/Hit.h
+++ b/algo/data/sts/Hit.h
@@ -26,6 +26,16 @@ namespace cbm::algo::sts
u32 fFrontClusterId; ///< Index of front-side cluster, used by tracking to reduce combinatorics
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; }
+
private: // serialization
friend class boost::serialization::access;
diff --git a/algo/detectors/tof/Clusterizer.cxx b/algo/detectors/tof/Clusterizer.cxx
index 14a9aa1e35e331aa5226352f5cb2cbaa395b94ef..ee370d980c6dc0f72d5bf1efe81ccf29209d84a3 100644
--- a/algo/detectors/tof/Clusterizer.cxx
+++ b/algo/detectors/tof/Clusterizer.cxx
@@ -120,8 +120,8 @@ namespace cbm::algo::tof
Clusterizer::resultType Clusterizer::buildClusters(std::vector<inputType>& input)
{
// Hit variables
- TofCluster cluster;
- std::vector<TofCluster> clustersOut;
+ Hit cluster;
+ std::vector<Hit> clustersOut;
// Reference cell of a cluster
TofCell* cell = nullptr;
@@ -231,8 +231,8 @@ namespace cbm::algo::tof
return clustersOut;
}
- bool Clusterizer::AddNextChan(std::vector<inputType>& input, int32_t lastChan, TofCluster& cluster,
- std::vector<TofCluster>& clustersOut, std::vector<inputType::iterator>* lastChanPos)
+ bool Clusterizer::AddNextChan(std::vector<inputType>& input, int32_t lastChan, Hit& cluster,
+ std::vector<Hit>& clustersOut, std::vector<inputType::iterator>* lastChanPos)
{
//D.Smith 25.8.23: Why are "C" digis (position "2") not considered here?
diff --git a/algo/detectors/tof/Clusterizer.h b/algo/detectors/tof/Clusterizer.h
index 611a34c1d44ba53b40a558a4722d0a48cfc1cde3..6771ba8a410bdb0b829349a2387b8819af8f3100 100644
--- a/algo/detectors/tof/Clusterizer.h
+++ b/algo/detectors/tof/Clusterizer.h
@@ -12,9 +12,8 @@
// TOF Classes and includes
class CbmTofDigi;
-// ROOT Classes and includes
-#include "Math/Rotation3D.h"
-#include "Math/Vector3Dfwd.h"
+#include "ClusterizerPars.h"
+#include "Hit.h"
// C++ Classes and includes
#include <memory>
@@ -24,117 +23,9 @@ class CbmTofDigi;
namespace cbm::algo::tof
{
- struct TofCell {
- double sizeX, sizeY;
- ROOT::Math::XYZVector pos;
- ROOT::Math::Rotation3D rotation;
- };
-
- struct ClusterizerChanPar {
- std::vector<double> fvCPTOff; //[nbSide]
- std::vector<double> fvCPTotGain; //[nbSide]
- std::vector<double> fvCPTotOff; //[nbSide]
- std::vector<std::vector<double>> fvCPWalk; //[nbSide][nbWalkBins]
- int32_t address; //unique address
- TofCell cell;
- };
-
- struct ClusterizerRpcPar {
- uint32_t fDeadStrips;
- double fPosYMaxScal;
- double fdMaxTimeDist;
- double fdMaxSpaceDist;
- double fSigVel;
- double fTimeRes;
- int32_t numClWalkBinX;
- double TOTMax;
- double TOTMin;
- bool swapChannelSides;
- double channelDeadtime;
- double fCPTOffYBinWidth;
- double fCPTOffYRange;
- std::vector<double> fCPTOffY; //[nBin]
- std::vector<ClusterizerChanPar> fChanPar = {};
- };
-
- struct TofCluster {
- ROOT::Math::XYZVector hitPos = ROOT::Math::XYZVector(0.0, 0.0, 0.0);
- ROOT::Math::XYZVector hitPosErr = ROOT::Math::XYZVector(0.0, 0.0, 0.0);
- double hitTime = 0.0;
- double hitTimeErr = 0.0;
- int32_t address = 0;
-
- std::vector<int32_t> vDigiIndRef;
- double weightsSum = 0.0;
-
- int32_t numChan() { return vDigiIndRef.size() / 2; }
-
- void reset()
- {
- vDigiIndRef.clear();
- hitPos = ROOT::Math::XYZVector(0.0, 0.0, 0.0);
- hitPosErr = ROOT::Math::XYZVector(0.5, 0.5, 0.5); //D.Smith 15.8.23: Set this back to zero eventually.
- hitTime = 0.0;
- hitTimeErr = 0.0;
- weightsSum = 0.0;
- address = 0;
- }
-
- void add(ROOT::Math::XYZVector pos, double time, double timeErr, double weight, int32_t digiIndA, int32_t digiIndB)
- {
- vDigiIndRef.push_back(digiIndA);
- vDigiIndRef.push_back(digiIndB);
- hitPos += pos * weight;
- hitTime += time * weight;
- hitTimeErr += timeErr * weight;
- weightsSum += weight;
- }
-
- void normalize(double timeErr)
- {
- // a/=b is translated to a := a*(1/b) in the ROOT::Math::XYZVector class, which has a different
- // rounding behavior than a := (a/b). In rare cases this leads to 1.000... becoming 0.999... inside
- // the floor() operation in the finalize() function of this class, and results in a different
- // channel being associated with the cluster. To reproduce the output of the old hit finder, we
- // divide element-wise instead. Perhaps floor() should be replaced by round().
- //////// hitPos /= weightsSum;
-
- hitPos.SetXYZ(hitPos.X() / weightsSum, hitPos.Y() / weightsSum, hitPos.Z() / weightsSum);
- hitTime /= weightsSum;
- hitTimeErr = timeErr;
- }
-
- void finalize(const TofCell& trafoCell, const ClusterizerRpcPar& par)
- {
- //get hit channel
- const int32_t channel = par.fChanPar.size() / 2 + floor(hitPos.X() / trafoCell.sizeX);
-
- // D.Smith 15.8.23: This channel correction doesn't seem to be needed
- // if (channel < 0) channel = 0;
- // if (channel > par.fChanPar.size() - 1) channel = par.fChanPar.size() - 1;
- address = par.fChanPar[channel].address;
-
- //Calibrate hit time if applicable
- if (par.fCPTOffYRange != 0) {
- const double dBin = (hitPos.Y() + par.fCPTOffYRange) / par.fCPTOffYBinWidth;
- const int i1 = floor(dBin);
- const int i2 = ceil(dBin);
- hitTime -= par.fCPTOffY[i1] + (dBin - i1) * (par.fCPTOffY[i2] - par.fCPTOffY[i1]);
- }
-
- //get offset by rotation to master frame. get hit position by adding offset to cell coordinates
- hitPos = trafoCell.pos + trafoCell.rotation(hitPos);
-
- // Simple errors, not properly done at all for now
- // Right way of doing it should take into account the weight distribution
- // and real system time resolution
- hitPosErr = ROOT::Math::XYZVector(0.5, 0.5, 0.5);
- }
- };
-
class Clusterizer {
public:
- typedef std::vector<TofCluster> resultType;
+ typedef std::vector<Hit> resultType;
typedef std::vector<std::pair<CbmTofDigi*, int32_t>> inputType;
/**
@@ -165,8 +56,8 @@ namespace cbm::algo::tof
resultType buildClusters(std::vector<inputType>& input);
- bool AddNextChan(std::vector<inputType>& input, int32_t iLastChan, TofCluster& cluster,
- std::vector<TofCluster>& clustersOut, std::vector<inputType::iterator>* lastChanPos = nullptr);
+ bool AddNextChan(std::vector<inputType>& input, int32_t iLastChan, Hit& cluster, std::vector<Hit>& clustersOut,
+ std::vector<inputType::iterator>* lastChanPos = nullptr);
};
} // namespace cbm::algo::tof
diff --git a/algo/detectors/tof/ClusterizerPars.h b/algo/detectors/tof/ClusterizerPars.h
new file mode 100644
index 0000000000000000000000000000000000000000..fea7d8036cd7e808bc55005de9cef1e696689455
--- /dev/null
+++ b/algo/detectors/tof/ClusterizerPars.h
@@ -0,0 +1,47 @@
+/* Copyright (C) 2022 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer], Pierre-Alain Loizeau, Felix Weiglhofer */
+#pragma once
+
+#include <vector>
+
+#include "Math/Rotation3D.h"
+#include "Math/Vector3Dfwd.h"
+
+
+namespace cbm::algo::tof
+{
+ struct TofCell {
+ double sizeX, sizeY;
+ ROOT::Math::XYZVector pos;
+ ROOT::Math::Rotation3D rotation;
+ };
+
+ struct ClusterizerChanPar {
+ std::vector<double> fvCPTOff; //[nbSide]
+ std::vector<double> fvCPTotGain; //[nbSide]
+ std::vector<double> fvCPTotOff; //[nbSide]
+ std::vector<std::vector<double>> fvCPWalk; //[nbSide][nbWalkBins]
+ int32_t address; //unique address
+ TofCell cell;
+ };
+
+ struct ClusterizerRpcPar {
+ uint32_t fDeadStrips;
+ double fPosYMaxScal;
+ double fdMaxTimeDist;
+ double fdMaxSpaceDist;
+ double fSigVel;
+ double fTimeRes;
+ int32_t numClWalkBinX;
+ double TOTMax;
+ double TOTMin;
+ bool swapChannelSides;
+ double channelDeadtime;
+ double fCPTOffYBinWidth;
+ double fCPTOffYRange;
+ std::vector<double> fCPTOffY; //[nBin]
+ std::vector<ClusterizerChanPar> fChanPar = {};
+ };
+
+} // namespace cbm::algo::tof
diff --git a/algo/detectors/tof/Hit.h b/algo/detectors/tof/Hit.h
new file mode 100644
index 0000000000000000000000000000000000000000..23899e9a5d08fee05d13b35d0d949b2573a06914
--- /dev/null
+++ b/algo/detectors/tof/Hit.h
@@ -0,0 +1,105 @@
+/* Copyright (C) 2022 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer], Pierre-Alain Loizeau, Felix Weiglhofer */
+#pragma once
+
+#include <cstdint>
+#include <vector>
+
+#include "ClusterizerPars.h"
+#include "Definitions.h"
+#include "Math/Rotation3D.h"
+#include "Math/Vector3Dfwd.h"
+
+namespace cbm::algo::tof
+{
+
+ struct Hit {
+ ROOT::Math::XYZVector hitPos = ROOT::Math::XYZVector(0.0, 0.0, 0.0);
+ ROOT::Math::XYZVector hitPosErr = ROOT::Math::XYZVector(0.0, 0.0, 0.0);
+ double hitTime = 0.0;
+ double hitTimeErr = 0.0;
+ int32_t address = 0;
+
+ std::vector<int32_t> vDigiIndRef;
+ double weightsSum = 0.0;
+
+ // 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; }
+
+ // Interface end
+
+ int32_t numChan() { return vDigiIndRef.size() / 2; }
+
+ void reset()
+ {
+ vDigiIndRef.clear();
+ hitPos = ROOT::Math::XYZVector(0.0, 0.0, 0.0);
+ hitPosErr = ROOT::Math::XYZVector(0.5, 0.5, 0.5); //D.Smith 15.8.23: Set this back to zero eventually.
+ hitTime = 0.0;
+ hitTimeErr = 0.0;
+ weightsSum = 0.0;
+ address = 0;
+ }
+
+ void add(ROOT::Math::XYZVector pos, double time, double timeErr, double weight, int32_t digiIndA, int32_t digiIndB)
+ {
+ vDigiIndRef.push_back(digiIndA);
+ vDigiIndRef.push_back(digiIndB);
+ hitPos += pos * weight;
+ hitTime += time * weight;
+ hitTimeErr += timeErr * weight;
+ weightsSum += weight;
+ }
+
+ void normalize(double timeErr)
+ {
+ // a/=b is translated to a := a*(1/b) in the ROOT::Math::XYZVector class, which has a different
+ // rounding behavior than a := (a/b). In rare cases this leads to 1.000... becoming 0.999... inside
+ // the floor() operation in the finalize() function of this class, and results in a different
+ // channel being associated with the cluster. To reproduce the output of the old hit finder, we
+ // divide element-wise instead. Perhaps floor() should be replaced by round().
+ //////// hitPos /= weightsSum;
+
+ hitPos.SetXYZ(hitPos.X() / weightsSum, hitPos.Y() / weightsSum, hitPos.Z() / weightsSum);
+ hitTime /= weightsSum;
+ hitTimeErr = timeErr;
+ }
+
+ void finalize(const TofCell& trafoCell, const ClusterizerRpcPar& par)
+ {
+ //get hit channel
+ const int32_t channel = par.fChanPar.size() / 2 + floor(hitPos.X() / trafoCell.sizeX);
+
+ // D.Smith 15.8.23: This channel correction doesn't seem to be needed
+ // if (channel < 0) channel = 0;
+ // if (channel > par.fChanPar.size() - 1) channel = par.fChanPar.size() - 1;
+ address = par.fChanPar[channel].address;
+
+ //Calibrate hit time if applicable
+ if (par.fCPTOffYRange != 0) {
+ const double dBin = (hitPos.Y() + par.fCPTOffYRange) / par.fCPTOffYBinWidth;
+ const int i1 = floor(dBin);
+ const int i2 = ceil(dBin);
+ hitTime -= par.fCPTOffY[i1] + (dBin - i1) * (par.fCPTOffY[i2] - par.fCPTOffY[i1]);
+ }
+
+ //get offset by rotation to master frame. get hit position by adding offset to cell coordinates
+ hitPos = trafoCell.pos + trafoCell.rotation(hitPos);
+
+ // Simple errors, not properly done at all for now
+ // Right way of doing it should take into account the weight distribution
+ // and real system time resolution
+ hitPosErr = ROOT::Math::XYZVector(0.5, 0.5, 0.5);
+ }
+ };
+
+} // namespace cbm::algo::tof
diff --git a/algo/detectors/tof/Hitfind.cxx b/algo/detectors/tof/Hitfind.cxx
index ba1f4bdae6c39cce97f27654c0cd2710b8047a59..e1c1379171370714b3bf52ed95774df405683d0a 100644
--- a/algo/detectors/tof/Hitfind.cxx
+++ b/algo/detectors/tof/Hitfind.cxx
@@ -16,31 +16,39 @@ using fles::Subsystem;
namespace cbm::algo::tof
{
// ----- Execution -------------------------------------------------------
- Hitfind::resultType Hitfind::operator()(std::vector<CbmTofDigi>* digiIn)
+ Hitfind::resultType Hitfind::operator()(gsl::span<CbmTofDigi> digiIn)
{
xpu::push_timer("TofHitfind");
// --- Output data
resultType result = {};
- std::vector<TofCluster>& clusterTs = result.first;
- HitfindMonitorData& monitor = result.second;
+ auto& clusterTs = result.first;
+ auto& monitor = result.second;
// Do calibration globally (optional, should not be used together with RPC-wise calibration)
- *digiIn = CalibRawDigis(*digiIn);
+ CalibRawDigis(digiIn, monitor);
// Loop over the digis array and store the Digis in separate vectors for
// each RPC modules
- for (int32_t idigi = 0; idigi < digiIn->size(); idigi++) {
- CbmTofDigi* pDigi = &(digiIn->at(idigi));
+ for (int32_t idigi = 0; idigi < digiIn.size(); idigi++) {
+ CbmTofDigi* pDigi = &(digiIn[idigi]);
// These are doubles in the digi class
const double SmType = pDigi->GetType();
const double Sm = pDigi->GetSm();
const double Rpc = pDigi->GetRpc();
const int NbRpc = fNbRpc[SmType];
+ if (SmType >= fNbSm.size() || Sm * NbRpc + Rpc >= fNbSm[SmType] * NbRpc) {
+ // Error already counted for monitoring during Digis calibration, so just discard here
+ continue;
+ }
fStorDigi[SmType][Sm * NbRpc + Rpc].push_back(std::make_pair(pDigi, idigi));
}
+ std::vector<Hit> clustersFlat;
+ std::vector<size_t> rpcSizes; // nClusters per RPC
+ std::vector<u32> rpcAddresses; // RPC addresses
+
// --- RPC loop
for (uint32_t SmType = 0; SmType < fNbSm.size(); SmType++) {
const int NbRpc = fNbRpc[SmType];
@@ -52,10 +60,18 @@ namespace cbm::algo::tof
std::vector<std::pair<CbmTofDigi*, int32_t>>& digiExp = fStorDigi[SmType][Sm * NbRpc + Rpc];
// Build clusters
- std::vector<cbm::algo::tof::TofCluster> clusters = fAlgo[SmType][Sm * NbRpc + Rpc](digiExp);
+ std::vector<cbm::algo::tof::Hit> clusters = fAlgo[SmType][Sm * NbRpc + Rpc](digiExp);
+
+ // Store hw address of partition
+ // TODO: is this correct? Does this address make sense?
+ u32 rpcAddress = CbmTofAddress::GetUniqueAddress(Sm, Rpc, 0, 0, SmType);
+ rpcAddresses.push_back(rpcAddress);
+
+ // store partition size
+ rpcSizes.push_back(clusters.size());
// Append clusters to output
- clusterTs.insert(clusterTs.end(), clusters.begin(), clusters.end());
+ clustersFlat.insert(clustersFlat.end(), clusters.begin(), clusters.end());
// Clear digi storage
digiExp.clear();
@@ -63,10 +79,15 @@ namespace cbm::algo::tof
}
}
+ // Monitoring
xpu::timings timings = xpu::pop_timer();
monitor.fTimeHitfind = timings.wall();
- monitor.fNumDigis = digiIn->size();
- monitor.fNumHits = clusterTs.size();
+ monitor.fNumDigis = digiIn.size();
+ monitor.fNumHits = clustersFlat.size();
+
+ // Create ouput vector
+ clusterTs = PartitionedVector(std::move(clustersFlat), rpcSizes, rpcAddresses);
+
return result;
}
// ----------------------------------------------------------------------------
@@ -141,22 +162,28 @@ namespace cbm::algo::tof
// -------------------------- Calibration -------------------------------------
- std::vector<CbmTofDigi> Hitfind::CalibRawDigis(const std::vector<CbmTofDigi>& digiVec)
+ void Hitfind::CalibRawDigis(gsl::span<CbmTofDigi> digiVec, HitfindMonitorData& monitor)
{
// channel deadtime map
std::map<int32_t, double> mChannelDeadTime;
- std::vector<CbmTofDigi> calDigiVecOut;
- for (int32_t iDigi = 0; iDigi < digiVec.size(); iDigi++) {
+ for (size_t iDigi = 0; iDigi < digiVec.size(); iDigi++) {
- CbmTofDigi pDigi = digiVec.at(iDigi);
+ CbmTofDigi pDigi = digiVec[iDigi];
const double SmType = pDigi.GetType();
const double Sm = pDigi.GetSm();
const double Rpc = pDigi.GetRpc();
const double Chan = pDigi.GetChannel();
const double Side = pDigi.GetSide();
const int NbRpc = fNbRpc[SmType];
- HitfindSetup::Rpc rpcPar = fTofConfig.rpcs[SmType][Sm * NbRpc + Rpc];
+
+ auto& rpcs = fTofConfig.rpcs;
+ if (SmType >= rpcs.size() || Sm * NbRpc + Rpc >= rpcs.at(SmType).size()) {
+ monitor.fDigiCalibUnknownRPC++;
+ continue;
+ }
+
+ HitfindSetup::Rpc rpcPar = fTofConfig.rpcs.at(SmType).at(Sm * NbRpc + Rpc);
HitfindSetup::Channel chanPar = rpcPar.chanPar[Chan];
if (rpcPar.swapChannelSides && 5 != SmType && 8 != SmType) {
@@ -203,14 +230,12 @@ namespace cbm::algo::tof
if (0 < iWx) { dWT -= dDTot * (walk[iWx - 1] - walk[iWx]); }
}
pCalDigi.SetTime(pCalDigi.GetTime() - dWT); // calibrate Digi Time
- calDigiVecOut.push_back(pCalDigi);
+ digiVec[iDigi] = pCalDigi;
}
/// Sort the buffers of hits due to the time offsets applied
- std::sort(calDigiVecOut.begin(), calDigiVecOut.end(),
+ std::sort(digiVec.begin(), digiVec.end(),
[](const CbmTofDigi& a, const CbmTofDigi& b) -> bool { return a.GetTime() < b.GetTime(); });
-
- return calDigiVecOut;
}
} // namespace cbm::algo::tof
diff --git a/algo/detectors/tof/Hitfind.h b/algo/detectors/tof/Hitfind.h
index 09f761f4d524f7637124f71efba53504a0cfb35b..897ce36b8d5daefef64947ac24407219a6bcb37f 100644
--- a/algo/detectors/tof/Hitfind.h
+++ b/algo/detectors/tof/Hitfind.h
@@ -5,7 +5,7 @@
#ifndef TOFHITFIND_H
#define TOFHITFIND_H 1
-#include "Timeslice.hpp"
+#include "CbmTofDigi.h"
#include "tof/Clusterizer.h"
#include "tof/HitfindSetup.h"
@@ -15,8 +15,7 @@
#include <sstream>
#include <vector>
-#include "DigiData.h"
-#include "PODVector.h"
+#include "PartitionedVector.h"
namespace cbm::algo::tof
{
@@ -31,6 +30,7 @@ namespace cbm::algo::tof
double fTimeHitfind = 0;
size_t fNumDigis = 0;
size_t fNumHits = 0;
+ size_t fDigiCalibUnknownRPC = 0;
std::string print() const
{
@@ -42,7 +42,7 @@ namespace cbm::algo::tof
};
/** @class Hitfind
- ** @brief Algo class for hitfinding
+ ** @brief Algo class for hitfinding
** @author Dominik Smith <d.smith@gsi.de>
** @since 16.10.2023
**
@@ -50,13 +50,15 @@ namespace cbm::algo::tof
class Hitfind {
public:
- typedef std::pair<std::vector<TofCluster>, HitfindMonitorData> resultType;
+ typedef std::pair<PartitionedVector<Hit>, HitfindMonitorData> resultType;
/** @brief Algorithm execution
** @param fles timeslice to hitfind
** @return pair: digi timeslice, monitoring data
+ **
+ ** @note Modifies input digis for time calibration
**/
- resultType operator()(std::vector<CbmTofDigi>* digiIn);
+ resultType operator()(gsl::span<CbmTofDigi> digiIn);
/** @brief Default constructor **/
Hitfind() {};
@@ -85,7 +87,7 @@ namespace cbm::algo::tof
std::vector<int32_t> fNbRpc;
/** @brief Applies calibration to input digis **/
- std::vector<CbmTofDigi> CalibRawDigis(const std::vector<CbmTofDigi>& digiVec);
+ void CalibRawDigis(gsl::span<CbmTofDigi> digiVec, HitfindMonitorData& monitor);
};
} // namespace cbm::algo::tof
diff --git a/algo/detectors/tof/HitfinderChain.cxx b/algo/detectors/tof/HitfinderChain.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..6d78c43b7cf88751373602ac87280a587edce73b
--- /dev/null
+++ b/algo/detectors/tof/HitfinderChain.cxx
@@ -0,0 +1,26 @@
+/* Copyright (C) 2023 FIAS Frankfurt Institute for Advanced Studies, Frankfurt / Main
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Felix Weiglhofer [committer] */
+#include "HitfinderChain.h"
+
+#include "config/Yaml.h"
+
+using namespace cbm::algo;
+using namespace cbm::algo::tof;
+
+void HitfinderChain::Init()
+{
+ auto setup = config::ReadFromFile<HitfindSetup>(Opts().ParamsDir() / "TofHitfinderPar.yaml");
+ fHitfind = std::make_unique<Hitfind>();
+ fHitfind->fTofConfig = setup;
+ fHitfind->Init();
+}
+
+HitfinderChain::ReturnType HitfinderChain::Run(gsl::span<CbmTofDigi> digis)
+{
+ auto ret = (*fHitfind)(digis);
+
+ L_(info) << "TS contains " << ret.first.NElements() << " TOF Hits";
+ L_(error) << "TOF Digis with unknown RPCs: " << ret.second.fDigiCalibUnknownRPC;
+ return ret;
+}
diff --git a/algo/detectors/tof/HitfinderChain.h b/algo/detectors/tof/HitfinderChain.h
new file mode 100644
index 0000000000000000000000000000000000000000..20d650aafcfdea44d587ed8d9652501a10e6a139
--- /dev/null
+++ b/algo/detectors/tof/HitfinderChain.h
@@ -0,0 +1,37 @@
+/* Copyright (C) 2023 FIAS Frankfurt Institute for Advanced Studies, Frankfurt / Main
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Felix Weiglhofer [committer] */
+#pragma once
+
+#include "CbmTofDigi.h"
+
+#include <gsl/span>
+#include <memory>
+
+#include "Hitfind.h"
+#include "HitfindSetup.h"
+#include "SubChain.h"
+
+namespace cbm::algo::tof
+{
+
+ class HitfinderChain : public SubChain {
+
+ public:
+ using ReturnType = Hitfind::resultType;
+
+ HitfinderChain() = default;
+ ~HitfinderChain() = default;
+
+ void Init();
+
+ /**
+ * @brief Run the hitfinder and apply time calibration to digis
+ */
+ ReturnType Run(gsl::span<CbmTofDigi> digis);
+
+ private:
+ std::unique_ptr<Hitfind> fHitfind;
+ };
+
+} // namespace cbm::algo::tof
diff --git a/algo/global/Reco.cxx b/algo/global/Reco.cxx
index c6ab97da629207f37c61e9369783b32c068de35c..54f169c521395cd154706668d8e5cdbc59551fa5 100644
--- a/algo/global/Reco.cxx
+++ b/algo/global/Reco.cxx
@@ -45,6 +45,7 @@ void Reco::Init(const Options& opts)
SetContext(&fContext);
fUnpack.SetContext(&fContext);
fStsHitFinder.SetContext(&fContext);
+ fTofHitFinder.SetContext(&fContext);
fTracking.SetContext(&fContext);
if (Opts().HistogramUri() != "") {
@@ -84,6 +85,9 @@ void Reco::Init(const Options& opts)
hitFinderPars.memory = Params().sts.memory;
fStsHitFinder.SetParameters(hitFinderPars);
+ // TOF Hitfinder
+ fTofHitFinder.Init();
+
// Tracking
fTracking.Init();
@@ -124,6 +128,7 @@ RecoResults Reco::Run(const fles::Timeslice& ts)
case RecoParams::UnpackMode::CPU:
unpackResult = fUnpack.Run(ts);
unpackMonitor = unpackResult.second;
+ QueueUnpackerMetrics(ts, unpackMonitor, unpackResult.first);
break;
}
}
@@ -135,27 +140,36 @@ RecoResults Reco::Run(const fles::Timeslice& ts)
auto [ev, mon] = fEventBuild->Run(unpackResult.first);
events = std::move(ev);
evbuildMonitor = mon;
+ QueueEvbuildMetrics(evbuildMonitor);
}
sts::HitfinderMonitor stsHitfinderMonitor;
PartitionedSpan<sts::Hit> stsHits;
if (Opts().HasStep(Step::LocalReco) && Opts().HasDetector(fles::Subsystem::STS)) {
auto stsResults = fStsHitFinder(unpackResult.first.fSts);
- stsHits = std::move(stsResults.hits);
+ stsHits = stsResults.hits;
stsHitfinderMonitor = stsResults.monitor;
+ QueueStsRecoMetrics(stsHitfinderMonitor);
+ }
+
+ PartitionedVector<tof::Hit> tofHits;
+ if (Opts().HasStep(Step::LocalReco) && Opts().HasDetector(fles::Subsystem::TOF)) {
+ auto [hits, monitor] = fTofHitFinder.Run(unpackResult.first.fTof);
+ tofHits = std::move(hits);
+ QueueTofRecoMetrics(monitor);
}
// --- Tracking
if (Opts().HasStep(Step::Tracking)) {
- TrackingChain::Return_t trackingOutput = fTracking.Run(results);
+ TrackingChain::Input_t input {
+ .stsHits = stsHits,
+ .tofHits = tofHits,
+ };
+ TrackingChain::Return_t trackingOutput = fTracking.Run(input);
if (Opts().HasOutput(RecoData::Track)) { results.tracks = std::move(trackingOutput.first); }
QueueTrackingMetrics(trackingOutput.second);
}
- QueueUnpackerMetrics(ts, unpackMonitor, unpackResult.first);
- QueueStsRecoMetrics(stsHitfinderMonitor);
- QueueEvbuildMetrics(evbuildMonitor);
-
if (Opts().HasOutput(RecoData::DigiEvent)) results.events = std::move(events);
if (Opts().HasOutput(RecoData::Hit)) results.stsHits = std::move(stsHits);
}
@@ -261,6 +275,19 @@ void Reco::QueueStsRecoMetrics(const sts::HitfinderMonitor& monitor)
});
}
+void Reco::QueueTofRecoMetrics(const tof::HitfindMonitorData& mon)
+{
+ if (!HasMonitor()) return;
+
+ GetMonitor().QueueMetric("cbmreco", {{"hostname", fles::system::current_hostname()}, {"child", Opts().ChildId()}},
+ {
+ {"tofRecoTimeTotal", mon.fTimeHitfind},
+ {"tofRecoNumDigisIn", mon.fNumDigis},
+ {"tofRecoNumHits", mon.fNumHits},
+ {"tofRecoUnknownRPC", mon.fDigiCalibUnknownRPC},
+ });
+}
+
void Reco::QueueEvbuildMetrics(const evbuild::EventbuildChainMonitorData& mon)
{
if (!HasMonitor()) return;
diff --git a/algo/global/Reco.h b/algo/global/Reco.h
index 7932985b6bb93d755d58416755a708690fbe74ba..3fe138a22bb0e646b7d4254c2aa009fcc489dab0 100644
--- a/algo/global/Reco.h
+++ b/algo/global/Reco.h
@@ -4,6 +4,8 @@
#ifndef CBM_ALGO_GLOBAL_RECO_H
#define CBM_ALGO_GLOBAL_RECO_H
+#include "tof/HitfinderChain.h"
+
#include <xpu/host.h>
#include "EventbuildChain.h"
@@ -50,6 +52,10 @@ namespace cbm::algo
UnpackChain fUnpack;
sts::HitfinderChain fStsHitFinder;
+ // TOF
+ tof::HitfinderChain fTofHitFinder;
+
+ // Eventbuilding
std::unique_ptr<evbuild::EventbuildChain> fEventBuild;
// Tracking
@@ -59,6 +65,7 @@ namespace cbm::algo
void QueueUnpackerMetrics(const fles::Timeslice&, const UnpackMonitorData&, const DigiData&);
void QueueStsRecoMetrics(const sts::HitfinderMonitor&);
+ void QueueTofRecoMetrics(const tof::HitfindMonitorData&);
void QueueEvbuildMetrics(const evbuild::EventbuildChainMonitorData&);
void QueueTrackingMetrics(const ca::TrackingMonitorData&);
};
diff --git a/algo/test/_GTestPartitionedSpan.cxx b/algo/test/_GTestPartitionedSpan.cxx
index 7afc70f50d6eb6a1728e606dab79701f62c0e032..4a730e27672644fb2b52a96e742713930325ef5f 100644
--- a/algo/test/_GTestPartitionedSpan.cxx
+++ b/algo/test/_GTestPartitionedSpan.cxx
@@ -9,7 +9,7 @@
using namespace cbm::algo;
template<typename T>
-void EXPECT_CONTAINER_EQ(gsl::span<T> a, std::vector<T> b)
+void EXPECT_CONTAINER_EQ(gsl::span<T> a, std::vector<i32> b)
{
EXPECT_EQ(a.size(), b.size());
for (size_t i = 0; i < a.size(); ++i) {
@@ -36,7 +36,8 @@ protected:
fAddressesInvalid = {0x0, 0x100};
}
- void Ensure(const PartitionedSpan<i32>& vec)
+ template<typename T> // T is either i32 or 'const i32'
+ void Ensure(const PartitionedSpan<T>& vec)
{
EXPECT_EQ(vec.NElements(), 9);
EXPECT_EQ(vec.NPartitions(), 3);
@@ -101,6 +102,14 @@ TEST_F(PartitionedSpanTest, IsCopyConstructable)
Ensure(vecCopy);
}
+TEST_F(PartitionedSpanTest, IsConstCopyConstructable)
+{
+ PartitionedSpan vec(fData, fOffsets, fAddresses);
+
+ PartitionedSpan<const i32> vecCopy(vec);
+ Ensure(vecCopy);
+}
+
TEST_F(PartitionedSpanTest, ThrowsOnNumAddressesMismatchesNumPartions)
{
EXPECT_THROW(PartitionedSpan vec(fData, fOffsets, fAddressesInvalid), std::runtime_error);
diff --git a/reco/tasks/CbmTaskTofClusterizer.cxx b/reco/tasks/CbmTaskTofClusterizer.cxx
index 9868752d191c8fa3a9774d365f17de8323469413..a6eaf631d2644c8cfaaff78700f337973b8a9a66 100644
--- a/reco/tasks/CbmTaskTofClusterizer.cxx
+++ b/reco/tasks/CbmTaskTofClusterizer.cxx
@@ -352,10 +352,10 @@ bool CbmTaskTofClusterizer::BuildClusters()
*fTofCalDigiVec = fTofDigiVec;
//call cluster finder
- std::vector<cbm::algo::tof::TofCluster> clusters = fAlgo(fTofCalDigiVec).first;
+ auto clusters = fAlgo(*fTofCalDigiVec).first;
//Store hits and match
- for (auto const& cluster : clusters) {
+ for (auto const& cluster : clusters.Data()) {
const int32_t hitIndex = fTofHitsColl->GetEntriesFast();
TVector3 hitpos = TVector3(cluster.hitPos.X(), cluster.hitPos.Y(), cluster.hitPos.Z());
TVector3 hiterr = TVector3(cluster.hitPosErr.X(), cluster.hitPosErr.Y(), cluster.hitPosErr.Z());