64 for (
const auto&
vertex : vertices) {
65 mPrimaryVertices.emplace_back(
vertex);
66 if (!isBeamPositionOverridden) {
67 const float w =
vertex.getNContributors();
68 mBeamPos[0] = (mBeamPos[0] * mBeamPosWeight +
vertex.getX() *
w) / (mBeamPosWeight +
w);
69 mBeamPos[1] = (mBeamPos[1] * mBeamPosWeight +
vertex.getY() *
w) / (mBeamPosWeight +
w);
73 mROFramesPV.push_back(mPrimaryVertices.size());
85 mVerticesMCRecInfo.insert(mVerticesMCRecInfo.end(), labels.begin(), labels.end());
91 mPrimaryVertices.insert(mPrimaryVertices.begin() + mROFramesPV[rofId], vertices.begin(), vertices.end());
92 for (
int i = rofId + 1;
i < mROFramesPV.size(); ++
i) {
93 mROFramesPV[
i] += vertices.size();
95 mTotVertPerIteration[iteration] += vertices.size();
101 mVerticesMCRecInfo.insert(mVerticesMCRecInfo.begin() + mROFramesPV[rofId], labels.begin(), labels.end());
104template <
int nLayers>
108 for (
const auto&
vertex : vertices) {
109 if (
vertex.getTimeStamp().getTimeStamp() < rofId) {
110 insertPastVertex(
vertex, iteration);
112 if (
vertex.getTimeStamp().getTimeStamp() > rofId) {
113 futureVertices.emplace_back(
vertex);
116 mPrimaryVertices.emplace_back(
vertex);
117 mTotVertPerIteration[iteration]++;
118 if (!isBeamPositionOverridden) {
119 const float w =
vertex.getNContributors();
120 mBeamPos[0] = (mBeamPos[0] * mBeamPosWeight +
vertex.getX() *
w) / (mBeamPosWeight +
w);
121 mBeamPos[1] = (mBeamPos[1] * mBeamPosWeight +
vertex.getY() *
w) / (mBeamPosWeight +
w);
125 mROFramesPV.push_back(mPrimaryVertices.size());
126 for (
auto&
vertex : futureVertices) {
127 mPrimaryVertices.emplace_back(
vertex);
128 mTotVertPerIteration[iteration]++;
132template <
int nLayers>
134 gsl::span<const itsmft::CompClusterExt>
clusters,
135 gsl::span<const unsigned char>::iterator& pattIt,
146 for (
auto& rof : rofs) {
147 for (
int clusterId{rof.getFirstEntry()}; clusterId < rof.getFirstEntry() + rof.getNEntries(); ++clusterId) {
152 auto pattID =
c.getPatternID();
172 mClusterSize[clusterId] = std::clamp(
clusterSize, 0u, 255u);
173 auto sensorID =
c.getSensorID();
179 addTrackingFrameInfoToLayer(
layer, gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), trkXYZ.x(), geom->
getSensorRefAlpha(sensorID),
180 std::array<float, 2>{trkXYZ.y(), trkXYZ.z()},
181 std::array<float, 3>{sigmaY2, sigmaYZ, sigmaZ2});
184 addClusterToLayer(
layer, gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), mUnsortedClusters[
layer].size());
185 addClusterExternalIndexToLayer(
layer, clusterId);
187 for (
unsigned int iL{0}; iL < mUnsortedClusters.size(); ++iL) {
188 mROFramesClusters[iL].push_back(mUnsortedClusters[iL].
size());
193 for (
auto i = 0;
i < mNTrackletsPerCluster.size(); ++
i) {
194 mNTrackletsPerCluster[
i].resize(mUnsortedClusters[1].
size());
195 mNTrackletsPerClusterSum[
i].resize(mUnsortedClusters[1].
size() + 1);
198 if (mcLabels !=
nullptr) {
199 mClusterLabels = mcLabels;
205template <
int nLayers>
208 for (
int iLayer{0}; iLayer < nLayers; ++iLayer) {
222template <
int nLayers>
226 const int stride{numBins + 1};
227 bounded_vector<ClusterHelper> cHelper(mMemoryPool.get());
228 bounded_vector<int> clsPerBin(numBins, 0, mMemoryPool.get());
229 bounded_vector<int> lutPerBin(numBins, 0, mMemoryPool.get());
230 for (
int rof{0}; rof < mNrof; ++rof) {
231 if ((
int)mMultiplicityCutMask.size() == mNrof && !mMultiplicityCutMask[rof]) {
234 for (
int iLayer{0}, stopLayer = std::min(trkParam.
NLayers, maxLayers); iLayer < stopLayer; ++iLayer) {
235 const auto& unsortedClusters{getUnsortedClustersOnLayer(rof, iLayer)};
236 const int clustersNum{
static_cast<int>(unsortedClusters.size())};
237 auto* tableBase = mIndexTables[iLayer].data() + rof *
stride;
239 cHelper.resize(clustersNum);
241 for (
int iCluster{0}; iCluster < clustersNum; ++iCluster) {
242 const Cluster&
c = unsortedClusters[iCluster];
243 ClusterHelper&
h = cHelper[iCluster];
245 const float x =
c.xCoordinate - mBeamPos[0];
246 const float y =
c.yCoordinate - mBeamPos[1];
247 const float z =
c.zCoordinate;
249 float phi = math_utils::computePhi(
x,
y);
250 int zBin{mIndexTableUtils.getZBinIndex(iLayer,
z)};
251 if (zBin < 0 || zBin >= trkParam.
ZBins) {
252 zBin = std::clamp(zBin, 0, trkParam.
ZBins - 1);
253 mBogusClusters[iLayer]++;
255 int bin = mIndexTableUtils.getBinIndex(zBin, mIndexTableUtils.getPhiBinIndex(phi));
257 h.r = math_utils::hypot(
x,
y);
258 mMinR[iLayer] = o2::gpu::GPUCommonMath::Min(
h.r, mMinR[iLayer]);
259 mMaxR[iLayer] = o2::gpu::GPUCommonMath::Max(
h.r, mMaxR[iLayer]);
261 h.ind = clsPerBin[bin]++;
263 std::exclusive_scan(clsPerBin.begin(), clsPerBin.end(), lutPerBin.begin(), 0);
265 auto clusters2beSorted{getClustersOnLayer(rof, iLayer)};
266 for (
int iCluster{0}; iCluster < clustersNum; ++iCluster) {
267 const ClusterHelper&
h = cHelper[iCluster];
270 c = unsortedClusters[iCluster];
273 c.indexTableBinIndex =
h.bin;
275 std::copy_n(lutPerBin.data(), clsPerBin.size(), tableBase);
276 std::fill_n(tableBase + clsPerBin.size(),
stride - clsPerBin.size(), clustersNum);
278 std::fill(clsPerBin.begin(), clsPerBin.end(), 0);
284template <
int nLayers>
287 if (iteration == 0) {
288 if (maxLayers < trkParam.
NLayers && resetVertices) {
310 mIndexTableUtils.setTrackingParameters(trkParam);
314 for (
unsigned int iLayer{0}; iLayer < std::min((
int)mClusters.size(), maxLayers); ++iLayer) {
323 for (
int iLayer{0}; iLayer < trkParam.
NLayers; ++iLayer) {
325 for (
auto& tfInfo : mTrackingFrameInfo[iLayer]) {
327 tfInfo.covarianceTrackingFrame[0] += trkParam.
SystErrorY2[iLayer];
328 tfInfo.covarianceTrackingFrame[2] += trkParam.
SystErrorZ2[iLayer];
333 mNTrackletsPerROF.resize(2);
334 for (
auto&
v : mNTrackletsPerROF) {
335 v = bounded_vector<int>(mNrof + 1, 0, mMemoryPool.get());
337 if (iteration == 0 || iteration == 3) {
338 prepareClusters(trkParam, maxLayers);
340 mTotalTracklets = {0, 0};
341 if (maxLayers < trkParam.
NLayers) {
342 for (
size_t iLayer{0}; iLayer < maxLayers; ++iLayer) {
348 mTotVertPerIteration.resize(1 + iteration);
353 mMSangles.resize(trkParam.
NLayers);
354 mPhiCuts.resize(mClusters.size() - 1, 0.f);
356 float oneOverR{0.001f * 0.3f * std::abs(mBz) / trkParam.
TrackletMinPt};
357 for (
unsigned int iLayer{0}; iLayer <
nLayers; ++iLayer) {
360 if (iLayer < mClusters.size() - 1) {
361 const float& r1 = trkParam.
LayerRadii[iLayer];
362 const float& r2 = trkParam.
LayerRadii[iLayer + 1];
363 const float res1 = o2::gpu::CAMath::Hypot(trkParam.
PVres, mPositionResolution[iLayer]);
364 const float res2 = o2::gpu::CAMath::Hypot(trkParam.
PVres, mPositionResolution[iLayer + 1]);
365 const float cosTheta1half = o2::gpu::CAMath::Sqrt(1.f - math_utils::Sq(0.5f * r1 * oneOverR));
366 const float cosTheta2half = o2::gpu::CAMath::Sqrt(1.f - math_utils::Sq(0.5f * r2 * oneOverR));
367 float x = r2 * cosTheta1half - r1 * cosTheta2half;
368 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)));
369 mPhiCuts[iLayer] = std::min(o2::gpu::CAMath::ASin(0.5f *
x * oneOverR) + 2.f * mMSangles[iLayer] + delta,
o2::constants::math::PI * 0.5f);
373 for (
int iLayer{0}; iLayer < std::min((
int)mTracklets.size(), maxLayers); ++iLayer) {
376 if (iLayer < (
int)mCells.size()) {
379 mTrackletsLookupTable[iLayer].resize(mClusters[iLayer + 1].
size(), 0);
383 if (iLayer < (
int)mCells.size() - 1) {
391template <
int nLayers>
394 unsigned long size{0};
395 for (
const auto& trkl : mTracklets) {
398 for (
const auto&
cells : mCells) {
401 for (
const auto& cellsN : mCellsNeighbours) {
402 size +=
sizeof(
int) * cellsN.size();
404 return size +
sizeof(
Road<nLayers - 2>) * mRoads.size();
407template <
int nLayers>
410 LOGP(info,
"TimeFrame: Artefacts occupy {:.2f} MB", getArtefactsMemory() /
constants::MB);
413template <
int nLayers>
417 mPValphaX.reserve(mPrimaryVertices.size());
418 for (
auto& pv : mPrimaryVertices) {
419 mPValphaX.emplace_back(std::array<float, 2>{o2::gpu::CAMath::Hypot(pv.getX(), pv.getY()), math_utils::computePhi(pv.getX(), pv.getY())});
423template <
int nLayers>
426 for (ushort iLayer = 0; iLayer < 2; ++iLayer) {
427 for (
unsigned int iRof{0}; iRof < mNrof; ++iRof) {
428 if (mMultiplicityCutMask[iRof]) {
429 mTotalTracklets[iLayer] += mNTrackletsPerROF[iLayer][iRof];
432 std::exclusive_scan(mNTrackletsPerROF[iLayer].begin(), mNTrackletsPerROF[iLayer].
end(), mNTrackletsPerROF[iLayer].begin(), 0);
433 std::exclusive_scan(mNTrackletsPerCluster[iLayer].begin(), mNTrackletsPerCluster[iLayer].
end(), mNTrackletsPerClusterSum[iLayer].begin(), 0);
437template <
int nLayers>
440 for (uint32_t iLayer{0}; iLayer < getTracklets().size(); ++iLayer) {
443 for (uint32_t iTracklet{0}; iTracklet < getTracklets()[iLayer].size(); ++iTracklet) {
444 auto& trk = getTracklets()[iLayer][iTracklet];
445 int currentId{trk.firstClusterIndex};
446 if (currentId < prev) {
447 LOG(info) <<
"First Cluster Index not increasing monotonically on L:T:ID:Prev " << iLayer <<
"\t" << iTracklet <<
"\t" << currentId <<
"\t" << prev;
448 }
else if (currentId == prev) {
452 auto& lut{getTrackletsLookupTable()[iLayer - 1]};
453 if (
count != lut[prev + 1] - lut[prev]) {
454 LOG(info) <<
"LUT count broken " << iLayer - 1 <<
"\t" << prev <<
"\t" <<
count <<
"\t" << lut[prev + 1] <<
"\t" << lut[prev];
461 auto& lut{getTrackletsLookupTable()[iLayer - 1]};
462 if (iTracklet >= (uint32_t)(lut[currentId + 1]) || iTracklet < (uint32_t)(lut[currentId])) {
463 LOG(info) <<
"LUT broken: " << iLayer - 1 <<
"\t" << currentId <<
"\t" << iTracklet;
470template <
int nLayers>
475 for (
int iLayers{nLayers}; iLayers--;) {
476 mClusters[iLayers].clear();
477 mUnsortedClusters[iLayers].clear();
478 mTrackingFrameInfo[iLayers].clear();
479 mClusterExternalIndices[iLayers].clear();
480 mUsedClusters[iLayers].clear();
481 mROFramesClusters[iLayers].clear();
482 mNClustersPerROF[iLayers].clear();
484 for (
int i{2};
i--;) {
485 mTrackletsIndexROF[
i].clear();
489template <
int nLayers>
492 for (
auto& trkl : mTracklets) {
498template <
int nLayers>
501 LOG(info) <<
"-------- Tracklet LUT " <<
i;
503 for (
int j : mTrackletsLookupTable[
i]) {
506 LOG(info) << s.str();
507 LOG(info) <<
"--------";
510template <
int nLayers>
513 LOG(info) <<
"-------- Cell LUT " <<
i;
515 for (
int j : mCellsLookupTable[
i]) {
518 LOG(info) << s.str();
519 LOG(info) <<
"--------";
522template <
int nLayers>
525 for (
unsigned int i{0};
i < mTrackletsLookupTable.size(); ++
i) {
526 printTrackletLUTonLayer(
i);
530template <
int nLayers>
533 for (
unsigned int i{0};
i < mCellsLookupTable.size(); ++
i) {
534 printCellLUTonLayer(
i);
538template <
int nLayers>
541 LOG(info) <<
"Vertices in ROF (nROF = " << mNrof <<
", lut size = " << mROFramesPV.size() <<
")";
542 for (
unsigned int iR{0}; iR < mROFramesPV.size(); ++iR) {
543 LOG(info) << mROFramesPV[iR] <<
"\t";
545 LOG(info) <<
"\n\n Vertices:";
546 for (
unsigned int iV{0}; iV < mPrimaryVertices.size(); ++iV) {
547 LOG(info) << mPrimaryVertices[iV].getX() <<
"\t" << mPrimaryVertices[iV].getY() <<
"\t" << mPrimaryVertices[iV].getZ();
549 LOG(info) <<
"--------";
552template <
int nLayers>
555 LOG(info) <<
"--------";
556 for (
unsigned int iLayer{0}; iLayer < mROFramesClusters.size(); ++iLayer) {
557 LOG(info) <<
"Layer " << iLayer;
559 for (
auto value : mROFramesClusters[iLayer]) {
562 LOG(info) << s.str();
566template <
int nLayers>
569 LOG(info) <<
"--------";
570 for (
unsigned int iLayer{0}; iLayer < mNClustersPerROF.size(); ++iLayer) {
571 LOG(info) <<
"Layer " << iLayer;
573 for (
auto&
value : mNClustersPerROF[iLayer]) {
576 LOG(info) << s.str();
580template <
int nLayers>
583 LOG(info) <<
"Dumping slice of " << sliceSize <<
" rofs:";
584 for (
int iROF{startROF}; iROF < startROF + sliceSize; ++iROF) {
585 LOG(info) <<
"ROF " << iROF <<
" dump:";
586 for (
unsigned int iLayer{0}; iLayer < mClusters.size(); ++iLayer) {
587 LOG(info) <<
"Layer " << iLayer <<
" has: " << getClustersOnLayer(iROF, iLayer).size() <<
" clusters.";
589 LOG(info) <<
"Number of seeding vertices: " << getPrimaryVertices(iROF).size();
591 for (
auto&
v : getPrimaryVertices(iROF)) {
592 LOG(info) <<
"\t vertex " << iVertex++ <<
": x=" <<
v.getX() <<
" " <<
" y=" <<
v.getY() <<
" z=" <<
v.getZ() <<
" has " <<
v.getNContributors() <<
" contributors.";
597template <
int nLayers>
604 auto alloc =
vec.get_allocator().resource();
605 if (alloc != mMemoryPool.get()) {
609 auto initArrays = [&]<
typename T,
size_t S>(std::array<bounded_vector<T>,
S>& arr) {
610 for (
size_t i{0};
i <
S; ++
i) {
611 auto alloc = arr[
i].get_allocator().resource();
612 if (alloc != mMemoryPool.get()) {
617 auto initVectors = [&]<
typename T>(std::vector<bounded_vector<T>>&
vec) {
618 for (
size_t i{0};
i <
vec.size(); ++
i) {
619 auto alloc =
vec[
i].get_allocator().resource();
620 if (alloc != mMemoryPool.get()) {
626 initVector(mTotVertPerIteration);
627 initVector(mPrimaryVertices);
628 initVector(mROFramesPV);
629 initArrays(mClusters);
630 initArrays(mTrackingFrameInfo);
631 initArrays(mClusterExternalIndices);
632 initArrays(mROFramesClusters);
633 initArrays(mNTrackletsPerCluster);
634 initArrays(mNTrackletsPerClusterSum);
635 initArrays(mNClustersPerROF);
636 initArrays(mIndexTables);
637 initArrays(mUsedClusters);
638 initArrays(mUnsortedClusters);
639 initVector(mROFramesPV);
640 initVector(mPrimaryVertices);
642 initVector(mRoadLabels);
643 initVector(mMSangles);
644 initVector(mPhiCuts);
645 initVector(mPositionResolution);
646 initVector(mClusterSize);
647 initVector(mPValphaX);
648 initVector(mBogusClusters);
649 initArrays(mTrackletsIndexROF);
650 initVectors(mTracks);
651 initVectors(mTracklets);
653 initVectors(mCellsNeighbours);
654 initVectors(mCellsLookupTable);
657template <
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
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 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 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 addPrimaryVertices(const bounded_vector< Vertex > &vertices)
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