64 for (
const auto&
vertex : vertices) {
65 mPrimaryVertices.emplace_back(
vertex);
66 mTotVertPerIteration[iteration]++;
67 if (!isBeamPositionOverridden) {
68 const float w =
vertex.getNContributors();
69 mBeamPos[0] = (mBeamPos[0] * mBeamPosWeight +
vertex.getX() *
w) / (mBeamPosWeight +
w);
70 mBeamPos[1] = (mBeamPos[1] * mBeamPosWeight +
vertex.getY() *
w) / (mBeamPosWeight +
w);
74 mROFramesPV.push_back(mPrimaryVertices.size());
80 mVerticesMCRecInfo.insert(mVerticesMCRecInfo.end(), labels.begin(), labels.end());
86 mVerticesContributorLabels.insert(mVerticesContributorLabels.end(), labels.begin(), labels.end());
92 mPrimaryVertices.insert(mPrimaryVertices.begin() + mROFramesPV[rofId], vertices.begin(), vertices.end());
93 for (
int i = rofId + 1;
i < mROFramesPV.size(); ++
i) {
94 mROFramesPV[
i] += vertices.size();
96 mTotVertPerIteration[iteration] += vertices.size();
102 mVerticesMCRecInfo.insert(mVerticesMCRecInfo.begin() + mROFramesPV[rofId], labels.begin(), labels.end());
105template <
int nLayers>
110 const auto& pvs = getPrimaryVertices(0, rofId);
111 for (
const auto& pv : pvs) {
112 n += pv.getNContributors();
114 mVerticesContributorLabels.insert(mVerticesContributorLabels.begin() +
n, labels.begin(), labels.end());
117template <
int nLayers>
119 gsl::span<const itsmft::CompClusterExt>
clusters,
120 gsl::span<const unsigned char>::iterator& pattIt,
131 for (
auto& rof : rofs) {
132 for (
int clusterId{rof.getFirstEntry()}; clusterId < rof.getFirstEntry() + rof.getNEntries(); ++clusterId) {
137 auto pattID =
c.getPatternID();
157 mClusterSize[clusterId] = std::clamp(
clusterSize, 0u, 255u);
158 auto sensorID =
c.getSensorID();
164 addTrackingFrameInfoToLayer(
layer, gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), trkXYZ.x(), geom->
getSensorRefAlpha(sensorID),
165 std::array<float, 2>{trkXYZ.y(), trkXYZ.z()},
166 std::array<float, 3>{sigmaY2, sigmaYZ, sigmaZ2});
169 addClusterToLayer(
layer, gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), mUnsortedClusters[
layer].size());
170 addClusterExternalIndexToLayer(
layer, clusterId);
172 for (
unsigned int iL{0}; iL < mUnsortedClusters.size(); ++iL) {
173 mROFramesClusters[iL].push_back(mUnsortedClusters[iL].
size());
178 for (
auto i = 0;
i < mNTrackletsPerCluster.size(); ++
i) {
179 mNTrackletsPerCluster[
i].resize(mUnsortedClusters[1].
size());
180 mNTrackletsPerClusterSum[
i].resize(mUnsortedClusters[1].
size() + 1);
183 if (mcLabels !=
nullptr) {
184 mClusterLabels = mcLabels;
190template <
int nLayers>
193 for (
int iLayer{0}; iLayer < nLayers; ++iLayer) {
207template <
int nLayers>
208void TimeFrame<nLayers>::prepareClusters(
const TrackingParameters& trkParam,
const int maxLayers)
210 const int numBins{trkParam.PhiBins * trkParam.ZBins};
211 const int stride{numBins + 1};
212 bounded_vector<ClusterHelper> cHelper(mMemoryPool.get());
213 bounded_vector<int> clsPerBin(numBins, 0, mMemoryPool.get());
214 bounded_vector<int> lutPerBin(numBins, 0, mMemoryPool.get());
215 for (
int rof{0}; rof < mNrof; ++rof) {
216 if ((
int)mMultiplicityCutMask.size() == mNrof && !mMultiplicityCutMask[rof]) {
219 for (
int iLayer{0}, stopLayer = std::min(trkParam.NLayers, maxLayers); iLayer < stopLayer; ++iLayer) {
220 const auto& unsortedClusters{getUnsortedClustersOnLayer(rof, iLayer)};
221 const int clustersNum{
static_cast<int>(unsortedClusters.size())};
222 auto* tableBase = mIndexTables[iLayer].data() + rof *
stride;
224 cHelper.resize(clustersNum);
226 for (
int iCluster{0}; iCluster < clustersNum; ++iCluster) {
227 const Cluster&
c = unsortedClusters[iCluster];
228 ClusterHelper&
h = cHelper[iCluster];
230 const float x =
c.xCoordinate - mBeamPos[0];
231 const float y =
c.yCoordinate - mBeamPos[1];
232 const float z =
c.zCoordinate;
234 float phi = math_utils::computePhi(
x,
y);
235 int zBin{mIndexTableUtils.getZBinIndex(iLayer,
z)};
236 if (zBin < 0 || zBin >= trkParam.ZBins) {
237 zBin = std::clamp(zBin, 0, trkParam.ZBins - 1);
238 mBogusClusters[iLayer]++;
240 int bin = mIndexTableUtils.getBinIndex(zBin, mIndexTableUtils.getPhiBinIndex(phi));
242 h.r = math_utils::hypot(
x,
y);
243 mMinR[iLayer] = o2::gpu::GPUCommonMath::Min(
h.r, mMinR[iLayer]);
244 mMaxR[iLayer] = o2::gpu::GPUCommonMath::Max(
h.r, mMaxR[iLayer]);
246 h.ind = clsPerBin[bin]++;
248 std::exclusive_scan(clsPerBin.begin(), clsPerBin.end(), lutPerBin.begin(), 0);
250 auto clusters2beSorted{getClustersOnLayer(rof, iLayer)};
251 for (
int iCluster{0}; iCluster < clustersNum; ++iCluster) {
252 const ClusterHelper&
h = cHelper[iCluster];
253 Cluster&
c = clusters2beSorted[lutPerBin[
h.bin] +
h.ind];
255 c = unsortedClusters[iCluster];
258 c.indexTableBinIndex =
h.bin;
260 std::copy_n(lutPerBin.data(), clsPerBin.size(), tableBase);
261 std::fill_n(tableBase + clsPerBin.size(),
stride - clsPerBin.size(), clustersNum);
263 std::fill(clsPerBin.begin(), clsPerBin.end(), 0);
269template <
int nLayers>
272 if (iteration == 0) {
273 if (maxLayers < trkParam.
NLayers && resetVertices) {
296 mIndexTableUtils.setTrackingParameters(trkParam);
300 for (
unsigned int iLayer{0}; iLayer < std::min((
int)mClusters.size(), maxLayers); ++iLayer) {
309 for (
int iLayer{0}; iLayer < trkParam.
NLayers; ++iLayer) {
311 for (
auto& tfInfo : mTrackingFrameInfo[iLayer]) {
313 tfInfo.covarianceTrackingFrame[0] += trkParam.
SystErrorY2[iLayer];
314 tfInfo.covarianceTrackingFrame[2] += trkParam.
SystErrorZ2[iLayer];
319 mNTrackletsPerROF.resize(2);
320 for (
auto&
v : mNTrackletsPerROF) {
321 v = bounded_vector<int>(mNrof + 1, 0, mMemoryPool.get());
323 if (iteration == 0 || iteration == 3) {
324 prepareClusters(trkParam, maxLayers);
326 mTotalTracklets = {0, 0};
327 if (maxLayers < trkParam.
NLayers) {
328 for (
size_t iLayer{0}; iLayer < maxLayers; ++iLayer) {
334 mTotVertPerIteration.resize(1 + iteration);
339 mMSangles.resize(trkParam.
NLayers);
340 mPhiCuts.resize(mClusters.size() - 1, 0.f);
342 float oneOverR{0.001f * 0.3f * std::abs(mBz) / trkParam.
TrackletMinPt};
343 for (
unsigned int iLayer{0}; iLayer <
nLayers; ++iLayer) {
346 if (iLayer < mClusters.size() - 1) {
347 const float& r1 = trkParam.
LayerRadii[iLayer];
348 const float& r2 = trkParam.
LayerRadii[iLayer + 1];
349 const float res1 = o2::gpu::CAMath::Hypot(trkParam.
PVres, mPositionResolution[iLayer]);
350 const float res2 = o2::gpu::CAMath::Hypot(trkParam.
PVres, mPositionResolution[iLayer + 1]);
351 const float cosTheta1half = o2::gpu::CAMath::Sqrt(1.f - math_utils::Sq(0.5f * r1 * oneOverR));
352 const float cosTheta2half = o2::gpu::CAMath::Sqrt(1.f - math_utils::Sq(0.5f * r2 * oneOverR));
353 float x = r2 * cosTheta1half - r1 * cosTheta2half;
354 float delta = o2::gpu::CAMath::Sqrt(1.f / (1.f - 0.25f * math_utils::Sq(
x * oneOverR)) * (math_utils::Sq(0.25f * r1 * r2 * math_utils::Sq(oneOverR) / cosTheta2half + cosTheta1half) * math_utils::Sq(res1) + math_utils::Sq(0.25f * r1 * r2 * math_utils::Sq(oneOverR) / cosTheta1half + cosTheta2half) * math_utils::Sq(res2)));
355 mPhiCuts[iLayer] = std::min(o2::gpu::CAMath::ASin(0.5f *
x * oneOverR) + 2.f * mMSangles[iLayer] + delta,
o2::constants::math::PI * 0.5f);
359 for (
int iLayer{0}; iLayer < std::min((
int)mTracklets.size(), maxLayers); ++iLayer) {
362 if (iLayer < (
int)mCells.size()) {
365 mTrackletsLookupTable[iLayer].resize(mClusters[iLayer + 1].
size() + 1, 0);
369 if (iLayer < (
int)mCells.size() - 1) {
377template <
int nLayers>
380 unsigned long size{0};
381 for (
const auto& trkl : mTracklets) {
384 for (
const auto&
cells : mCells) {
387 for (
const auto& cellsN : mCellsNeighbours) {
388 size +=
sizeof(
int) * cellsN.size();
390 return size +
sizeof(
Road<nLayers - 2>) * mRoads.size();
393template <
int nLayers>
396 LOGP(info,
"TimeFrame: Artefacts occupy {:.2f} MB", getArtefactsMemory() /
constants::MB);
399template <
int nLayers>
403 mPValphaX.reserve(mPrimaryVertices.size());
404 for (
auto& pv : mPrimaryVertices) {
405 mPValphaX.emplace_back(std::array<float, 2>{o2::gpu::CAMath::Hypot(pv.getX(), pv.getY()), math_utils::computePhi(pv.getX(), pv.getY())});
409template <
int nLayers>
412 for (ushort iLayer = 0; iLayer < 2; ++iLayer) {
413 for (
unsigned int iRof{0}; iRof < mNrof; ++iRof) {
414 if (mMultiplicityCutMask[iRof]) {
415 mTotalTracklets[iLayer] += mNTrackletsPerROF[iLayer][iRof];
418 std::exclusive_scan(mNTrackletsPerROF[iLayer].begin(), mNTrackletsPerROF[iLayer].
end(), mNTrackletsPerROF[iLayer].begin(), 0);
419 std::exclusive_scan(mNTrackletsPerCluster[iLayer].begin(), mNTrackletsPerCluster[iLayer].
end(), mNTrackletsPerClusterSum[iLayer].begin(), 0);
423template <
int nLayers>
426 for (uint32_t iLayer{0}; iLayer < getTracklets().size(); ++iLayer) {
429 for (uint32_t iTracklet{0}; iTracklet < getTracklets()[iLayer].size(); ++iTracklet) {
430 auto& trk = getTracklets()[iLayer][iTracklet];
431 int currentId{trk.firstClusterIndex};
432 if (currentId < prev) {
433 LOG(info) <<
"First Cluster Index not increasing monotonically on L:T:ID:Prev " << iLayer <<
"\t" << iTracklet <<
"\t" << currentId <<
"\t" << prev;
434 }
else if (currentId == prev) {
438 auto& lut{getTrackletsLookupTable()[iLayer - 1]};
439 if (
count != lut[prev + 1] - lut[prev]) {
440 LOG(info) <<
"LUT count broken " << iLayer - 1 <<
"\t" << prev <<
"\t" <<
count <<
"\t" << lut[prev + 1] <<
"\t" << lut[prev];
447 auto& lut{getTrackletsLookupTable()[iLayer - 1]};
448 if (iTracklet >= (uint32_t)(lut[currentId + 1]) || iTracklet < (uint32_t)(lut[currentId])) {
449 LOG(info) <<
"LUT broken: " << iLayer - 1 <<
"\t" << currentId <<
"\t" << iTracklet;
456template <
int nLayers>
461 for (
int iLayers{nLayers}; iLayers--;) {
462 mClusters[iLayers].clear();
463 mUnsortedClusters[iLayers].clear();
464 mTrackingFrameInfo[iLayers].clear();
465 mClusterExternalIndices[iLayers].clear();
466 mUsedClusters[iLayers].clear();
467 mROFramesClusters[iLayers].clear();
468 mNClustersPerROF[iLayers].clear();
470 for (
int i{2};
i--;) {
471 mTrackletsIndexROF[
i].clear();
475template <
int nLayers>
478 for (
auto& trkl : mTracklets) {
484template <
int nLayers>
487 LOG(info) <<
"-------- Tracklet LUT " <<
i;
489 for (
int j : mTrackletsLookupTable[
i]) {
492 LOG(info) << s.str();
493 LOG(info) <<
"--------";
496template <
int nLayers>
499 LOG(info) <<
"-------- Cell LUT " <<
i;
501 for (
int j : mCellsLookupTable[
i]) {
504 LOG(info) << s.str();
505 LOG(info) <<
"--------";
508template <
int nLayers>
511 for (
unsigned int i{0};
i < mTrackletsLookupTable.size(); ++
i) {
512 printTrackletLUTonLayer(
i);
516template <
int nLayers>
519 for (
unsigned int i{0};
i < mCellsLookupTable.size(); ++
i) {
520 printCellLUTonLayer(
i);
524template <
int nLayers>
527 LOG(info) <<
"Vertices in ROF (nROF = " << mNrof <<
", lut size = " << mROFramesPV.size() <<
")";
528 for (
unsigned int iR{0}; iR < mROFramesPV.size(); ++iR) {
529 LOG(info) << mROFramesPV[iR] <<
"\t";
531 LOG(info) <<
"\n\n Vertices:";
532 for (
unsigned int iV{0}; iV < mPrimaryVertices.size(); ++iV) {
533 LOG(info) << mPrimaryVertices[iV].getX() <<
"\t" << mPrimaryVertices[iV].getY() <<
"\t" << mPrimaryVertices[iV].getZ();
535 LOG(info) <<
"--------";
538template <
int nLayers>
541 LOG(info) <<
"--------";
542 for (
unsigned int iLayer{0}; iLayer < mROFramesClusters.size(); ++iLayer) {
543 LOG(info) <<
"Layer " << iLayer;
545 for (
auto value : mROFramesClusters[iLayer]) {
548 LOG(info) << s.str();
552template <
int nLayers>
555 LOG(info) <<
"--------";
556 for (
unsigned int iLayer{0}; iLayer < mNClustersPerROF.size(); ++iLayer) {
557 LOG(info) <<
"Layer " << iLayer;
559 for (
auto&
value : mNClustersPerROF[iLayer]) {
562 LOG(info) << s.str();
566template <
int nLayers>
569 LOG(info) <<
"Dumping slice of " << sliceSize <<
" rofs:";
570 for (
int iROF{startROF}; iROF < startROF + sliceSize; ++iROF) {
571 LOG(info) <<
"ROF " << iROF <<
" dump:";
572 for (
unsigned int iLayer{0}; iLayer < mClusters.size(); ++iLayer) {
573 LOG(info) <<
"Layer " << iLayer <<
" has: " << getClustersOnLayer(iROF, iLayer).size() <<
" clusters.";
575 LOG(info) <<
"Number of seeding vertices: " << getPrimaryVertices(iROF).size();
577 for (
auto&
v : getPrimaryVertices(iROF)) {
578 LOG(info) <<
"\t vertex " << iVertex++ <<
": x=" <<
v.getX() <<
" " <<
" y=" <<
v.getY() <<
" z=" <<
v.getZ() <<
" has " <<
v.getNContributors() <<
" contributors.";
583template <
int nLayers>
589 auto alloc =
vec.get_allocator().resource();
590 if (alloc != mMemoryPool.get()) {
594 auto initArrays = [&]<
typename T,
size_t S>(std::array<bounded_vector<T>,
S>& arr) {
595 for (
size_t i{0};
i <
S; ++
i) {
596 auto alloc = arr[
i].get_allocator().resource();
597 if (alloc != mMemoryPool.get()) {
602 auto initVectors = [&]<
typename T>(std::vector<bounded_vector<T>>&
vec) {
603 for (
size_t i{0};
i <
vec.size(); ++
i) {
604 auto alloc =
vec[
i].get_allocator().resource();
605 if (alloc != mMemoryPool.get()) {
611 initVector(mTotVertPerIteration);
612 initVector(mPrimaryVertices);
613 initVector(mROFramesPV);
614 initArrays(mClusters);
615 initArrays(mTrackingFrameInfo);
616 initArrays(mClusterExternalIndices);
617 initArrays(mROFramesClusters);
618 initArrays(mNTrackletsPerCluster);
619 initArrays(mNTrackletsPerClusterSum);
620 initArrays(mNClustersPerROF);
621 initArrays(mIndexTables);
622 initArrays(mUsedClusters);
623 initArrays(mUnsortedClusters);
624 initVector(mROFramesPV);
625 initVector(mPrimaryVertices);
627 initVector(mRoadLabels);
628 initVector(mMSangles);
629 initVector(mPhiCuts);
630 initVector(mPositionResolution);
631 initVector(mClusterSize);
632 initVector(mPValphaX);
633 initVector(mBogusClusters);
634 initVector(mVerticesContributorLabels);
635 initArrays(mTrackletsIndexROF);
636 initVectors(mTracks);
637 initVectors(mTracklets);
639 initVectors(mCellsNeighbours);
640 initVectors(mCellsLookupTable);
643template <
int nLayers>
Definition of the ITSMFT compact cluster.
Definition of the GeometryTGeo class.
Definition of the SegmentationAlpide class.
Class for time synchronization of RawReader instances.
const Mat3D & getMatrixL2G(int sensID) const
HMPID cluster implementation.
const Mat3D & getMatrixT2L(int lay, int hba, int sta, int det) const
int getLayer(int index) const
Get chip layer, from 0.
float getSensorRefAlpha(int isn) const
static GeometryTGeo * Instance()
void fillMatrixCache(int mask) override
int getNPixels() const
Returns the number of fired pixels.
static constexpr unsigned short InvalidPatternID
static constexpr float PitchCol
static constexpr float PitchRow
math_utils::Point3D< T > getClusterCoordinates(const CompCluster &cl) const
float getErr2X(int n) const
Returns the error^2 on the x position of the COG for the n_th element.
int getNpixels(int n) const
Returns the number of fired pixels of the n_th element.
bool isGroup(int n) const
Returns true if the element corresponds to a group of rare topologies.
float getErr2Z(int n) const
Returns the error^2 on the z position of the COG for the n_th element.
GLsizei const GLfloat * value
GLint GLenum GLboolean GLsizei stride
GLenum GLuint GLint GLint layer
GLubyte GLubyte GLubyte GLubyte w
GLdouble GLdouble GLdouble z
constexpr unsigned int nLayers
void deepVectorClear(std::vector< T > &vec)
std::pmr::vector< T > bounded_vector
void clearResizeBoundedArray(std::array< bounded_vector< T >, S > &arr, size_t size, BoundedMemoryResource *bmr, T def=T())
constexpr float DefClusError2Col
constexpr float DefClusError2Row
void clearResizeBoundedVector(bounded_vector< T > &vec, size_t size, BoundedMemoryResource *bmr, T def=T())
constexpr float DefClusErrorRow
constexpr float DefClusErrorCol
void checkTrackletLUTs()
Debug and printing.
void printTrackletLUTonLayer(int i)
void initialise(const int iteration, const TrackingParameters &trkParam, const int maxLayers=7, bool resetVertices=true)
void addPrimaryVertices(const bounded_vector< Vertex > &vertices, const int iteration)
void addPrimaryVerticesInROF(const bounded_vector< Vertex > &vertices, const int rofId, const int iteration)
void printSliceInfo(const int, const int)
void computeTrackletsPerROFScans()
void addPrimaryVerticesLabelsInROF(const bounded_vector< std::pair< MCCompLabel, float > > &labels, const int rofId)
void addPrimaryVerticesContributorLabelsInROF(const bounded_vector< MCCompLabel > &labels, const int rofId)
void addPrimaryVerticesLabels(bounded_vector< std::pair< MCCompLabel, float > > &labels)
int loadROFrameData(const o2::itsmft::ROFRecord &rof, gsl::span< const itsmft::Cluster > clusters, const dataformats::MCTruthContainer< MCCompLabel > *mcLabels=nullptr)
void fillPrimaryVerticesXandAlpha()
void setMemoryPool(std::shared_ptr< BoundedMemoryResource > &pool)
memory management
unsigned long getArtefactsMemory() const
void printCellLUTonLayer(int i)
void printArtefactsMemory() const
void addPrimaryVerticesContributorLabels(bounded_vector< MCCompLabel > &labels)
std::vector< float > LayerRadii
std::vector< float > SystErrorY2
std::vector< float > SystErrorZ2
float TrackletMinPt
Trackleting cuts.
std::vector< float > LayerResolution
std::vector< float > LayerxX0
int CellsPerRoad() const noexcept
int TrackletsPerRoad() const noexcept
std::vector< o2::ctf::BufferType > vec
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"
std::vector< Cluster > clusters
std::vector< Cell > cells