Project
Loading...
Searching...
No Matches
Clusterer.cxx
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
11
14
15#include <algorithm>
16
21
22#include "TTree.h"
23
24#ifdef WITH_OPENMP
25#include <omp.h>
26#endif
27using namespace o2::itsmft;
28
29namespace o2::its3
30{
31void Clusterer::process(int nThreads, PixelReader& reader, CompClusCont* compClus,
32 PatternCont* patterns, ROFRecCont* vecROFRec, MCTruth* labelsCl)
33{
34#ifdef _PERFORM_TIMING_
35 mTimer.Start(false);
36#endif
37 if (nThreads < 1) {
38 nThreads = 1;
39 }
40 auto autoDecode = reader.getDecodeNextAuto();
41 do {
42 if (autoDecode) {
43 reader.setDecodeNextAuto(false); // internally do not autodecode
44 if (!reader.decodeNextTrigger()) {
45 break; // on the fly decoding was requested, but there were no data left
46 }
47 }
48 if (reader.getInteractionRecord().isDummy()) {
49 continue; // No IR info was found
50 }
51 // pre-fetch all non-empty chips of current ROF
52 ChipPixelData* curChipData = nullptr;
53 mFiredChipsPtr.clear();
54 size_t nPix = 0;
55 while ((curChipData = reader.getNextChipData(mChips))) {
56 mFiredChipsPtr.push_back(curChipData);
57 nPix += curChipData->getData().size();
58 }
59
60 auto& rof = vecROFRec->emplace_back(reader.getInteractionRecord(), vecROFRec->size(), compClus->size(), 0); // create new ROF
61
62 uint16_t nFired = mFiredChipsPtr.size();
63 if (!nFired) {
64 if (autoDecode) {
65 continue;
66 }
67 break; // just 1 ROF was asked to be processed
68 }
69 if (nFired < nThreads) {
70 nThreads = nFired;
71 }
72#ifndef WITH_OPENMP
73 nThreads = 1;
74#endif
75 uint16_t chipStep = nThreads > 1 ? (nThreads == 2 ? 20 : 10) : nFired;
76 int dynGrp = std::min(4, std::max(1, nThreads / 2));
77 if (nThreads > mThreads.size()) {
78 int oldSz = mThreads.size();
79 mThreads.resize(nThreads);
80 for (size_t i = oldSz; i < nThreads; i++) {
81 mThreads[i] = std::make_unique<ClustererThread>(this, i);
82 }
83 }
84#ifdef WITH_OPENMP
85#pragma omp parallel for schedule(dynamic, dynGrp) num_threads(nThreads)
86 //>> start of MT region
87 for (uint16_t ic = 0; ic < nFired; ic += chipStep) {
88 auto ith = omp_get_thread_num();
89 if (nThreads > 1) {
90 mThreads[ith]->process(ic, std::min(chipStep, uint16_t(nFired - ic)),
91 &mThreads[ith]->compClusters,
92 patterns ? &mThreads[ith]->patterns : nullptr,
93 labelsCl ? reader.getDigitsMCTruth() : nullptr,
94 labelsCl ? &mThreads[ith]->labels : nullptr, rof);
95 } else { // put directly to the destination
96 mThreads[0]->process(0, nFired, compClus, patterns, labelsCl ? reader.getDigitsMCTruth() : nullptr, labelsCl, rof);
97 }
98 }
99 //<< end of MT region
100#else
101 mThreads[0]->process(0, nFired, compClus, patterns, labelsCl ? reader.getDigitsMCTruth() : nullptr, labelsCl, rof);
102#endif
103 // copy data of all threads but the 1st one to final destination
104 if (nThreads > 1) {
105#ifdef _PERFORM_TIMING_
106 mTimerMerge.Start(false);
107#endif
108 size_t nClTot = 0, nPattTot = 0;
109 int chid = 0;
110 std::vector<size_t> thrStatIdx(nThreads, 0);
111 for (int ith = 0; ith < nThreads; ith++) {
112 std::sort(mThreads[ith]->stats.begin(), mThreads[ith]->stats.end(), [](const ThreadStat& a, const ThreadStat& b) { return a.firstChip < b.firstChip; });
113 nClTot += mThreads[ith]->compClusters.size();
114 nPattTot += mThreads[ith]->patterns.size();
115 }
116 compClus->reserve(nClTot);
117 if (patterns) {
118 patterns->reserve(nPattTot);
119 }
120 while (chid < nFired) {
121 for (int ith = 0; ith < nThreads; ith++) {
122 if (thrStatIdx[ith] >= mThreads[ith]->stats.size()) {
123 continue;
124 }
125 const auto& stat = mThreads[ith]->stats[thrStatIdx[ith]];
126 if (stat.firstChip == chid) {
127 thrStatIdx[ith]++;
128 chid += stat.nChips; // next chip to look
129 const auto clbeg = mThreads[ith]->compClusters.begin() + stat.firstClus;
130 compClus->insert(compClus->end(), clbeg, clbeg + stat.nClus);
131 if (patterns) {
132 const auto ptbeg = mThreads[ith]->patterns.begin() + stat.firstPatt;
133 patterns->insert(patterns->end(), ptbeg, ptbeg + stat.nPatt);
134 }
135 if (labelsCl) {
136 labelsCl->mergeAtBack(mThreads[ith]->labels, stat.firstClus, stat.nClus);
137 }
138 }
139 }
140 }
141 for (int ith = 0; ith < nThreads; ith++) {
142 mThreads[ith]->patterns.clear();
143 mThreads[ith]->compClusters.clear();
144 mThreads[ith]->labels.clear();
145 mThreads[ith]->stats.clear();
146 }
147#ifdef _PERFORM_TIMING_
148 mTimerMerge.Stop();
149#endif
150 } else {
151 mThreads[0]->stats.clear();
152 }
153 rof.setNEntries(compClus->size() - rof.getFirstEntry()); // update
154 } while (autoDecode);
155 reader.setDecodeNextAuto(autoDecode); // restore setting
156#ifdef _PERFORM_TIMING_
157 mTimer.Stop();
158#endif
159}
160
161//__________________________________________________
162void Clusterer::ClustererThread::process(uint16_t chip, uint16_t nChips, CompClusCont* compClusPtr, PatternCont* patternsPtr,
163 const ConstMCTruth* labelsDigPtr, MCTruth* labelsClPtr, const ROFRecord& rofPtr)
164{
165 if (stats.empty() || stats.back().firstChip + stats.back().nChips != chip) { // there is a jump, register new block
166 stats.emplace_back(ThreadStat{chip, 0, uint32_t(compClusPtr->size()), patternsPtr ? uint32_t(patternsPtr->size()) : 0, 0, 0});
167 }
168 for (int ic = 0; ic < nChips; ic++) {
169 auto* curChipData = parent->mFiredChipsPtr[chip + ic];
170 auto chipID = curChipData->getChipID();
171 if (parent->mMaxBCSeparationToMask > 0) { // mask pixels fired from the previous ROF
172 const auto& chipInPrevROF = parent->mChipsOld[chipID];
173 if (std::abs(rofPtr.getBCData().differenceInBC(chipInPrevROF.getInteractionRecord())) < parent->mMaxBCSeparationToMask) {
174 parent->mMaxRowColDiffToMask != 0 ? curChipData->maskFiredInSample(parent->mChipsOld[chipID], parent->mMaxRowColDiffToMask) : curChipData->maskFiredInSample(parent->mChipsOld[chipID]);
175 }
176 }
177 auto validPixID = curChipData->getFirstUnmasked();
178 auto npix = curChipData->getData().size();
179 if (validPixID < npix) { // chip data may have all of its pixels masked!
180 auto valp = validPixID++;
181 if (validPixID == npix) { // special case of a single pixel fired on the chip
182 finishChipSingleHitFast(valp, curChipData, compClusPtr, patternsPtr, labelsDigPtr, labelsClPtr);
183 } else {
184 initChip(curChipData, valp);
185 for (; validPixID < npix; validPixID++) {
186 if (!curChipData->getData()[validPixID].isMasked()) {
187 updateChip(curChipData, validPixID);
188 }
189 }
190 finishChip(curChipData, compClusPtr, patternsPtr, labelsDigPtr, labelsClPtr);
191 }
192 }
193 if (parent->mMaxBCSeparationToMask > 0) { // current chip data will be used in the next ROF to mask overflow pixels
194 parent->mChipsOld[chipID].swap(*curChipData);
195 }
196 }
197 auto& currStat = stats.back();
198 currStat.nChips += nChips;
199 currStat.nClus = compClusPtr->size() - currStat.firstClus;
200 currStat.nPatt = patternsPtr ? (patternsPtr->size() - currStat.firstPatt) : 0;
201}
202
203//__________________________________________________
205 PatternCont* patternsPtr, const ConstMCTruth* labelsDigPtr, MCTruth* labelsClusPtr)
206{
207 const auto& pixData = curChipData->getData();
208 for (size_t i1 = 0; i1 < preClusterHeads.size(); ++i1) {
209 auto ci = preClusterIndices[i1];
210 if (ci < 0) {
211 continue;
212 }
213 BBox bbox(curChipData->getChipID());
214 int nlab = 0;
215 int next = preClusterHeads[i1];
216 pixArrBuff.clear();
217 while (next >= 0) {
218 const auto& pixEntry = pixels[next];
219 const auto pix = pixData[pixEntry.second];
220 pixArrBuff.push_back(pix); // needed for cluster topology
221 bbox.adjust(pix.getRowDirect(), pix.getCol());
222 if (labelsClusPtr) {
223 if (parent->mSquashingDepth) { // the MCtruth for this pixel is stored in chip data: due to squashing we lose contiguity
224 fetchMCLabels(curChipData->getOrderedPixId(pixEntry.second), labelsDigPtr, nlab);
225 } else { // the MCtruth for this pixel is at curChipData->startID+pixEntry.second
226 fetchMCLabels(pixEntry.second + curChipData->getStartID(), labelsDigPtr, nlab);
227 }
228 }
229 next = pixEntry.first;
230 }
231 preClusterIndices[i1] = -1;
232 for (size_t i2 = i1 + 1; i2 < preClusterHeads.size(); ++i2) {
233 if (preClusterIndices[i2] != ci) {
234 continue;
235 }
236 next = preClusterHeads[i2];
237 while (next >= 0) {
238 const auto& pixEntry = pixels[next];
239 const auto pix = pixData[pixEntry.second]; // PixelData
240 pixArrBuff.push_back(pix); // needed for cluster topology
241 bbox.adjust(pix.getRowDirect(), pix.getCol());
242 if (labelsClusPtr != nullptr) {
243 if (parent->mSquashingDepth) { // the MCtruth for this pixel is stored in chip data: due to squashing we lose contiguity
244 fetchMCLabels(curChipData->getOrderedPixId(pixEntry.second), labelsDigPtr, nlab);
245 } else { // the MCtruth for this pixel is at curChipData->startID+pixEntry.second
246 fetchMCLabels(pixEntry.second + curChipData->getStartID(), labelsDigPtr, nlab);
247 }
248 }
249 next = pixEntry.first;
250 }
251 preClusterIndices[i2] = -1;
252 }
253 if (bbox.isAcceptableSize()) {
254 parent->streamCluster(pixArrBuff, &labelsBuff, bbox, parent->mPattIdConverter, compClusPtr, patternsPtr, labelsClusPtr, nlab, constants::detID::isDetITS3(curChipData->getChipID()));
255 } else {
256 auto warnLeft = MaxHugeClusWarn - parent->mNHugeClus;
257 if (warnLeft > 0) {
258 LOGP(warn, "Splitting a huge cluster: chipID {}, rows {}:{} cols {}:{}{}", bbox.chipID, bbox.rowMin, bbox.rowMax, bbox.colMin, bbox.colMax,
259 warnLeft == 1 ? " (Further warnings will be muted)" : "");
260#ifdef WITH_OPENMP
261#pragma omp critical
262#endif
263 {
264 parent->mNHugeClus++;
265 }
266 }
267 BBox bboxT(bbox); // truncated box
268 std::vector<PixelData> pixbuf;
269 do {
270 bboxT.rowMin = bbox.rowMin;
271 bboxT.colMax = std::min(bbox.colMax, uint16_t(bboxT.colMin + o2::itsmft::ClusterPattern::MaxColSpan - 1));
272 do { // Select a subset of pixels fitting the reduced bounding box
273 bboxT.rowMax = std::min(bbox.rowMax, uint16_t(bboxT.rowMin + o2::itsmft::ClusterPattern::MaxRowSpan - 1));
274 for (const auto& pix : pixArrBuff) {
275 if (bboxT.isInside(pix.getRowDirect(), pix.getCol())) {
276 pixbuf.push_back(pix);
277 }
278 }
279 if (!pixbuf.empty()) { // Stream a piece of cluster only if the reduced bounding box is not empty
280 parent->streamCluster(pixbuf, &labelsBuff, bboxT, parent->mPattIdConverter, compClusPtr, patternsPtr, labelsClusPtr, nlab, constants::detID::isDetITS3(curChipData->getChipID()), true);
281 pixbuf.clear();
282 }
283 bboxT.rowMin = bboxT.rowMax + 1;
284 } while (bboxT.rowMin < bbox.rowMax);
285 bboxT.colMin = bboxT.colMax + 1;
286 } while (bboxT.colMin < bbox.colMax);
287 }
288 }
289}
290
291//__________________________________________________
293 PatternCont* patternsPtr, const ConstMCTruth* labelsDigPtr, MCTruth* labelsClusPtr)
294{
295 auto pix = curChipData->getData()[hit];
296 uint16_t row = pix.getRowDirect(), col = pix.getCol();
297
298 if (labelsClusPtr) { // MC labels were requested
299 int nlab = 0;
300 fetchMCLabels(curChipData->getStartID() + hit, labelsDigPtr, nlab);
301 auto cnt = compClusPtr->size();
302 for (int i = nlab; i--;) {
303 labelsClusPtr->addElement(cnt, labelsBuff[i]);
304 }
305 }
306
307 auto ib = constants::detID::isDetITS3(curChipData->getChipID());
308
309 // add to compact clusters, which must be always filled
310 unsigned char patt[ClusterPattern::MaxPatternBytes]{0x1 << (7 - (0 % 8))}; // unrolled 1 hit version of full loop in finishChip
311 uint16_t pattID = (parent->mPattIdConverter.size(ib) == 0) ? CompCluster::InvalidPatternID : parent->mPattIdConverter.findGroupID(1, 1, ib, patt);
312 if ((pattID == CompCluster::InvalidPatternID || parent->mPattIdConverter.isGroup(pattID, ib)) && patternsPtr) {
313 patternsPtr->emplace_back(1); // rowspan
314 patternsPtr->emplace_back(1); // colspan
315 patternsPtr->insert(patternsPtr->end(), std::begin(patt), std::begin(patt) + 1);
316 }
317 compClusPtr->emplace_back(row, col, pattID, curChipData->getChipID());
318}
319
320//__________________________________________________
322{
323#ifdef _PERFORM_TIMING_
324 mTimer.Stop();
325 mTimer.Reset();
326 mTimerMerge.Stop();
327 mTimerMerge.Reset();
328#endif
329}
330
331//__________________________________________________
333{
334 // init chip with the 1st unmasked pixel (entry "from" in the mChipData)
336 int chipId = curChipData->getChipID();
339 }
340
341 delete[] column1;
342 delete[] column2;
343 column1 = new int[size];
344 column2 = new int[size];
345 column1[0] = column1[size - 1] = -1;
346 column2[0] = column2[size - 1] = -1;
347 // init chip with the 1st unmasked pixel (entry "from" in the mChipData)
348 prev = column1 + 1;
349 curr = column2 + 1;
350 resetColumn(curr);
351
352 pixels.clear();
353 preClusterHeads.clear();
354 preClusterIndices.clear();
355 auto pix = curChipData->getData()[first];
356 currCol = pix.getCol();
357 curr[pix.getRowDirect()] = 0; // can use getRowDirect since the pixel is not masked
358 // start the first pre-cluster
359 preClusterHeads.push_back(0);
360 preClusterIndices.push_back(0);
361 pixels.emplace_back(-1, first); // id of current pixel
362 noLeftCol = true; // flag that there is no column on the left to check yet
363}
364
365//__________________________________________________
366void Clusterer::ClustererThread::updateChip(const ChipPixelData* curChipData, uint32_t ip)
367{
368 const auto pix = curChipData->getData()[ip];
369 uint16_t row = pix.getRowDirect(); // can use getRowDirect since the pixel is not masked
370 if (currCol != pix.getCol()) { // switch the buffers
371 swapColumnBuffers();
372 resetColumn(curr);
373 noLeftCol = false;
374 if (pix.getCol() > currCol + 1) {
375 // no connection with previous column, this pixel cannot belong to any of the
376 // existing preclusters, create a new precluster and flag to check only the row above for next pixels of this column
377 currCol = pix.getCol();
378 addNewPrecluster(ip, row);
379 noLeftCol = true;
380 return;
381 }
382 currCol = pix.getCol();
383 }
384
385 bool orphan = true;
386
387 if (noLeftCol) { // check only the row above
388 if (curr[row - 1] >= 0) {
389 expandPreCluster(ip, row, curr[row - 1]); // attach to the precluster of the previous row
390 return;
391 }
392 } else {
393#ifdef _ALLOW_DIAGONAL_ALPIDE_CLUSTERS_
394 int neighbours[]{curr[row - 1], prev[row], prev[row + 1], prev[row - 1]};
395#else
396 int neighbours[]{curr[row - 1], prev[row]};
397#endif
398 for (auto pci : neighbours) {
399 if (pci < 0) {
400 continue;
401 }
402 if (orphan) {
403 expandPreCluster(ip, row, pci); // attach to the adjascent precluster
404 orphan = false;
405 continue;
406 }
407 // reassign precluster index to smallest one
408 if (preClusterIndices[pci] < preClusterIndices[curr[row]]) {
409 preClusterIndices[curr[row]] = preClusterIndices[pci];
410 } else {
411 preClusterIndices[pci] = preClusterIndices[curr[row]];
412 }
413 }
414 }
415 if (orphan) {
416 addNewPrecluster(ip, row); // start new precluster
417 }
418}
419
420//__________________________________________________
421void Clusterer::ClustererThread::fetchMCLabels(int digID, const ConstMCTruth* labelsDig, int& nfilled)
422{
423 // transfer MC labels to cluster
424 if (nfilled >= MaxLabels) {
425 return;
426 }
427 const auto& lbls = labelsDig->getLabels(digID);
428 for (int i = lbls.size(); i--;) {
429 int ic = nfilled;
430 for (; ic--;) { // check if the label is already present
431 if (labelsBuff[ic] == lbls[i]) {
432 return; // label is found, do nothing
433 }
434 }
435 labelsBuff[nfilled++] = lbls[i];
436 if (nfilled >= MaxLabels) {
437 break;
438 }
439 }
440 //
441}
442
443//__________________________________________________
445{
446 // reset
447#ifdef _PERFORM_TIMING_
448 mTimer.Stop();
449 mTimer.Reset();
450 mTimerMerge.Stop();
451 mTimerMerge.Reset();
452#endif
453}
454
455//__________________________________________________
457{
458 LOGP(info, "Clusterizer has {} chips registered", mChips.size());
459 LOGP(info, "Clusterizer squashes overflow pixels separated by {} BC and <= {} in row/col seeking down to {} neighbour ROFs", mMaxBCSeparationToSquash, mMaxRowColDiffToMask, mSquashingDepth);
460 LOGP(info, "Clusterizer masks overflow pixels separated by < {} BC and <= {} in row/col", mMaxBCSeparationToMask, mMaxRowColDiffToMask);
461
462#ifdef _PERFORM_TIMING_
463 auto& tmr = const_cast<TStopwatch&>(mTimer); // ugly but this is what root does internally
464 auto& tmrm = const_cast<TStopwatch&>(mTimerMerge);
465 LOG(info) << "Inclusive clusterization timing (w/o disk IO): Cpu: " << tmr.CpuTime()
466 << " Real: " << tmr.RealTime() << " s in " << tmr.Counter() << " slots";
467 LOG(info) << "Threads output merging timing : Cpu: " << tmrm.CpuTime()
468 << " Real: " << tmrm.RealTime() << " s in " << tmrm.Counter() << " slots";
469
470#endif
471}
472} // namespace o2::its3
int32_t i
Definition of a container to keep Monte Carlo truth external to simulation objects.
uint32_t col
Definition RawData.h:4
Definition of the SegmentationMosaix class.
Definition of the ITS cluster finder.
gsl::span< const TruthElement > getLabels(uint32_t dataindex) const
void mergeAtBack(MCTruthContainer< TruthElement > const &other)
void addElement(uint32_t dataindex, TruthElement const &element, bool noElement=false)
static constexpr int MaxLabels
Definition Clusterer.h:76
void process(int nThreads, PixelReader &r, CompClusCont *compClus, PatternCont *patterns, ROFRecCont *vecROFRec, MCTruth *labelsCl=nullptr)
Definition Clusterer.cxx:31
void clear()
load the dictionary of cluster topologies
void print() const
static constexpr int MaxHugeClusWarn
Definition Clusterer.h:77
uint32_t getOrderedPixId(int pos) const
Definition PixelData.h:275
uint32_t getStartID() const
Definition PixelData.h:110
uint16_t getChipID() const
Definition PixelData.h:108
const std::vector< PixelData > & getData() const
Definition PixelData.h:115
static constexpr int MaxPatternBytes
static constexpr uint8_t MaxRowSpan
static constexpr uint8_t MaxColSpan
static constexpr unsigned short InvalidPatternID
Definition CompCluster.h:46
PixelReader class for the ITSMFT.
Definition PixelReader.h:34
virtual int decodeNextTrigger()=0
bool getDecodeNextAuto() const
Definition PixelReader.h:67
void setDecodeNextAuto(bool v)
Definition PixelReader.h:68
const o2::InteractionRecord & getInteractionRecord() const
Definition PixelReader.h:58
virtual const o2::dataformats::ConstMCTruthContainerView< o2::MCCompLabel > * getDigitsMCTruth() const
Definition PixelReader.h:49
virtual bool getNextChipData(ChipPixelData &chipData)=0
const BCData & getBCData() const
Definition ROFRecord.h:58
GLsizeiptr size
Definition glcorearb.h:659
GLint first
Definition glcorearb.h:399
GLboolean GLboolean GLboolean b
Definition glcorearb.h:1233
GLint GLint GLsizei GLint GLenum GLenum const void * pixels
Definition glcorearb.h:275
GLboolean GLboolean GLboolean GLboolean a
Definition glcorearb.h:1233
bool isDetITS3(T detID)
Definition SpecsV2.h:211
std::vector< unsigned char > PatternCont
Definition Clusterer.h:59
std::vector< itsmft::CompClusterExt > CompClusCont
Definition Clusterer.h:58
std::vector< itsmft::ROFRecord > ROFRecCont
Definition Clusterer.h:60
int64_t differenceInBC(const InteractionRecord &other) const
bool isInside(uint16_t row, uint16_t col) const
Definition Clusterer.h:86
std::vector< ThreadStat > stats
Definition Clusterer.h:147
void finishChipSingleHitFast(uint32_t hit, ChipPixelData *curChipData, CompClusCont *compClusPtr, PatternCont *patternsPtr, const ConstMCTruth *labelsDigPtr, MCTruth *labelsClusPTr)
void process(uint16_t chip, uint16_t nChips, CompClusCont *compClusPtr, PatternCont *patternsPtr, const ConstMCTruth *labelsDigPtr, MCTruth *labelsClPtr, const itsmft::ROFRecord &rofPtr)
void updateChip(const ChipPixelData *curChipData, uint32_t ip)
void fetchMCLabels(int digID, const ConstMCTruth *labelsDig, int &nfilled)
void initChip(const ChipPixelData *curChipData, uint32_t first)
void finishChip(ChipPixelData *curChipData, CompClusCont *compClus, PatternCont *patterns, const ConstMCTruth *labelsDig, MCTruth *labelsClus)
methods and transient data used within a thread
Definition Clusterer.h:114
std::vector< o2::mch::DsChannelId > chid
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"
std::vector< int > row