Skip to content

Commit bfa2ec2

Browse files
f3schmconcas
f3sch
authored and
mconcas
committed
ITS: stream comp computeLayerCells
Signed-off-by: Felix Schlepper <felix.schlepper@cern.ch>
1 parent 188af4a commit bfa2ec2

File tree

1 file changed

+156
-65
lines changed

1 file changed

+156
-65
lines changed

‎Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx‎

Lines changed: 156 additions & 65 deletions
Original file line numberDiff line numberDiff line change
@@ -271,11 +271,16 @@ void TrackerTraits<nLayers>::computeLayerCells(const int iteration)
271271
std::ofstream off(std::format("cells{}.txt", iter++));
272272
#endif
273273

274+
constexpr float radl = 9.36f; // Radiation length of Si [cm]
275+
constexpr float rho = 2.33f; // Density of Si [g/cm^3]
276+
274277
for (int iLayer = 0; iLayer < mTrkParams[iteration].CellsPerRoad(); ++iLayer) {
275-
mTimeFrame->getCells()[iLayer].clear();
276-
mTimeFrame->getCellsLabel(iLayer).clear();
278+
deepVectorClear(mTimeFrame->getCells()[iLayer]);
277279
if (iLayer > 0) {
278-
mTimeFrame->getCellsLookupTable()[iLayer - 1].clear();
280+
deepVectorClear(mTimeFrame->getCellsLookupTable()[iLayer - 1]);
281+
}
282+
if (mTimeFrame->hasMCinformation()) {
283+
deepVectorClear(mTimeFrame->getCellsLabel(iLayer));
279284
}
280285
}
281286

@@ -295,87 +300,173 @@ void TrackerTraits<nLayers>::computeLayerCells(const int iteration)
295300
resolution = resolution > 1.e-12 ? resolution : 1.f;
296301
#endif
297302

303+
// count number of cells found
298304
const int currentLayerTrackletsNum{static_cast<int>(mTimeFrame->getTracklets()[iLayer].size())};
299-
for (int iTracklet{0}; iTracklet < currentLayerTrackletsNum; ++iTracklet) {
300-
301-
const Tracklet& currentTracklet{mTimeFrame->getTracklets()[iLayer][iTracklet]};
302-
const int nextLayerClusterIndex{currentTracklet.secondClusterIndex};
303-
const int nextLayerFirstTrackletIndex{
304-
mTimeFrame->getTrackletsLookupTable()[iLayer][nextLayerClusterIndex]};
305-
const int nextLayerLastTrackletIndex{
306-
mTimeFrame->getTrackletsLookupTable()[iLayer][nextLayerClusterIndex + 1]};
307-
308-
if (nextLayerFirstTrackletIndex == nextLayerLastTrackletIndex) {
309-
continue;
310-
}
305+
bounded_vector<int> perTrackletCount(currentLayerTrackletsNum + 1, 0, mMemoryPool.get());
306+
tbb::parallel_for(
307+
tbb::blocked_range<int>(0, currentLayerTrackletsNum),
308+
[&](const tbb::blocked_range<int>& Tracklets) {
309+
for (int iTracklet = Tracklets.begin(); iTracklet < Tracklets.end(); ++iTracklet) {
310+
const Tracklet& currentTracklet{mTimeFrame->getTracklets()[iLayer][iTracklet]};
311+
const int nextLayerClusterIndex{currentTracklet.secondClusterIndex};
312+
const int nextLayerFirstTrackletIndex{
313+
mTimeFrame->getTrackletsLookupTable()[iLayer][nextLayerClusterIndex]};
314+
const int nextLayerLastTrackletIndex{
315+
mTimeFrame->getTrackletsLookupTable()[iLayer][nextLayerClusterIndex + 1]};
316+
317+
if (nextLayerFirstTrackletIndex == nextLayerLastTrackletIndex) {
318+
continue;
319+
}
311320

312-
for (int iNextTracklet{nextLayerFirstTrackletIndex}; iNextTracklet < nextLayerLastTrackletIndex; ++iNextTracklet) {
313-
if (mTimeFrame->getTracklets()[iLayer + 1][iNextTracklet].firstClusterIndex != nextLayerClusterIndex) {
314-
break;
315-
}
316-
const Tracklet& nextTracklet{mTimeFrame->getTracklets()[iLayer + 1][iNextTracklet]};
317-
const float deltaTanLambda{std::abs(currentTracklet.tanLambda - nextTracklet.tanLambda)};
321+
int foundCells{0};
322+
for (int iNextTracklet{nextLayerFirstTrackletIndex}; iNextTracklet < nextLayerLastTrackletIndex; ++iNextTracklet) {
323+
if (mTimeFrame->getTracklets()[iLayer + 1][iNextTracklet].firstClusterIndex != nextLayerClusterIndex) {
324+
break;
325+
}
326+
const Tracklet& nextTracklet{mTimeFrame->getTracklets()[iLayer + 1][iNextTracklet]};
327+
const float deltaTanLambda{std::abs(currentTracklet.tanLambda - nextTracklet.tanLambda)};
318328

319329
#ifdef OPTIMISATION_OUTPUT
320-
bool good{mTimeFrame->getTrackletsLabel(iLayer)[iTracklet] == mTimeFrame->getTrackletsLabel(iLayer + 1)[iNextTracklet]};
321-
float signedDelta{currentTracklet.tanLambda - nextTracklet.tanLambda};
322-
off << std::format("{}\t{:d}\t{}\t{}\t{}\t{}", iLayer, good, signedDelta, signedDelta / (mTrkParams[iteration].CellDeltaTanLambdaSigma), tanLambda, resolution) << std::endl;
330+
bool good{mTimeFrame->getTrackletsLabel(iLayer)[iTracklet] == mTimeFrame->getTrackletsLabel(iLayer + 1)[iNextTracklet]};
331+
float signedDelta{currentTracklet.tanLambda - nextTracklet.tanLambda};
332+
off << std::format("{}\t{:d}\t{}\t{}\t{}\t{}", iLayer, good, signedDelta, signedDelta / (mTrkParams[iteration].CellDeltaTanLambdaSigma), tanLambda, resolution) << std::endl;
323333
#endif
324334

325-
if (deltaTanLambda / mTrkParams[iteration].CellDeltaTanLambdaSigma < mTrkParams[iteration].NSigmaCut) {
335+
if (deltaTanLambda / mTrkParams[iteration].CellDeltaTanLambdaSigma < mTrkParams[iteration].NSigmaCut) {
326336

327-
/// Track seed preparation. Clusters are numbered progressively from the innermost going outward.
328-
const int clusId[3]{
329-
mTimeFrame->getClusters()[iLayer][currentTracklet.firstClusterIndex].clusterId,
330-
mTimeFrame->getClusters()[iLayer + 1][nextTracklet.firstClusterIndex].clusterId,
331-
mTimeFrame->getClusters()[iLayer + 2][nextTracklet.secondClusterIndex].clusterId};
332-
const auto& cluster1_glo = mTimeFrame->getUnsortedClusters()[iLayer].at(clusId[0]);
333-
const auto& cluster2_glo = mTimeFrame->getUnsortedClusters()[iLayer + 1].at(clusId[1]);
334-
const auto& cluster3_tf = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer + 2).at(clusId[2]);
335-
auto track{buildTrackSeed(cluster1_glo, cluster2_glo, cluster3_tf)};
337+
/// Track seed preparation. Clusters are numbered progressively from the innermost going outward.
338+
const int clusId[3]{
339+
mTimeFrame->getClusters()[iLayer][currentTracklet.firstClusterIndex].clusterId,
340+
mTimeFrame->getClusters()[iLayer + 1][nextTracklet.firstClusterIndex].clusterId,
341+
mTimeFrame->getClusters()[iLayer + 2][nextTracklet.secondClusterIndex].clusterId};
342+
const auto& cluster1_glo = mTimeFrame->getUnsortedClusters()[iLayer].at(clusId[0]);
343+
const auto& cluster2_glo = mTimeFrame->getUnsortedClusters()[iLayer + 1].at(clusId[1]);
344+
const auto& cluster3_tf = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer + 2).at(clusId[2]);
345+
auto track{buildTrackSeed(cluster1_glo, cluster2_glo, cluster3_tf)};
336346

337-
float chi2{0.f};
338-
bool good{false};
339-
for (int iC{2}; iC--;) {
340-
const TrackingFrameInfo& trackingHit = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer + iC).at(clusId[iC]);
347+
float chi2{0.f};
348+
bool good{false};
349+
for (int iC{2}; iC--;) {
350+
const TrackingFrameInfo& trackingHit = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer + iC).at(clusId[iC]);
341351

342-
if (!track.rotate(trackingHit.alphaTrackingFrame)) {
343-
break;
344-
}
352+
if (!track.rotate(trackingHit.alphaTrackingFrame)) {
353+
break;
354+
}
345355

346-
if (!track.propagateTo(trackingHit.xTrackingFrame, getBz())) {
347-
break;
348-
}
356+
if (!track.propagateTo(trackingHit.xTrackingFrame, getBz())) {
357+
break;
358+
}
349359

350-
constexpr float radl = 9.36f; // Radiation length of Si [cm]
351-
constexpr float rho = 2.33f; // Density of Si [g/cm^3]
352-
if (!track.correctForMaterial(mTrkParams[0].LayerxX0[iLayer + iC], mTrkParams[0].LayerxX0[iLayer] * radl * rho, true)) {
353-
break;
354-
}
360+
if (!track.correctForMaterial(mTrkParams[0].LayerxX0[iLayer + iC], mTrkParams[0].LayerxX0[iLayer] * radl * rho, true)) {
361+
break;
362+
}
355363

356-
auto predChi2{track.getPredictedChi2Quiet(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)};
357-
if (!track.o2::track::TrackParCov::update(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)) {
358-
break;
359-
}
360-
if (!iC && predChi2 > mTrkParams[iteration].MaxChi2ClusterAttachment) {
361-
break;
364+
const auto predChi2{track.getPredictedChi2Quiet(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)};
365+
if (!iC && predChi2 > mTrkParams[iteration].MaxChi2ClusterAttachment) {
366+
break;
367+
}
368+
369+
if (!track.o2::track::TrackParCov::update(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)) {
370+
break;
371+
}
372+
373+
good = !iC;
374+
chi2 += predChi2;
375+
}
376+
if (good) {
377+
++foundCells;
378+
}
362379
}
363-
good = !iC;
364-
chi2 += predChi2;
365380
}
366-
if (!good) {
381+
perTrackletCount[iTracklet] = foundCells;
382+
}
383+
});
384+
385+
// calculate offset table and check if any cells where found
386+
std::exclusive_scan(perTrackletCount.begin(), perTrackletCount.end(), perTrackletCount.begin(), 0);
387+
auto totalCells{perTrackletCount.back()};
388+
if (totalCells == 0) {
389+
continue;
390+
}
391+
auto& layerCells = mTimeFrame->getCells()[iLayer];
392+
layerCells.resize(totalCells);
393+
394+
tbb::parallel_for(
395+
tbb::blocked_range<int>(0, currentLayerTrackletsNum),
396+
[&](const tbb::blocked_range<int>& Tracklets) {
397+
for (int iTracklet = Tracklets.begin(); iTracklet < Tracklets.end(); ++iTracklet) {
398+
if (perTrackletCount[iTracklet] == perTrackletCount[iTracklet + 1]) {
367399
continue;
368400
}
369-
if (iLayer > 0 && (int)mTimeFrame->getCellsLookupTable()[iLayer - 1].size() <= iTracklet) {
370-
mTimeFrame->getCellsLookupTable()[iLayer - 1].resize(iTracklet + 1, mTimeFrame->getCells()[iLayer].size());
401+
402+
const Tracklet& currentTracklet{mTimeFrame->getTracklets()[iLayer][iTracklet]};
403+
const int nextLayerClusterIndex{currentTracklet.secondClusterIndex};
404+
const int nextLayerFirstTrackletIndex{
405+
mTimeFrame->getTrackletsLookupTable()[iLayer][nextLayerClusterIndex]};
406+
const int nextLayerLastTrackletIndex{
407+
mTimeFrame->getTrackletsLookupTable()[iLayer][nextLayerClusterIndex + 1]};
408+
409+
int position = perTrackletCount[iTracklet];
410+
for (int iNextTracklet{nextLayerFirstTrackletIndex}; iNextTracklet < nextLayerLastTrackletIndex; ++iNextTracklet) {
411+
if (mTimeFrame->getTracklets()[iLayer + 1][iNextTracklet].firstClusterIndex != nextLayerClusterIndex) {
412+
break;
413+
}
414+
const Tracklet& nextTracklet{mTimeFrame->getTracklets()[iLayer + 1][iNextTracklet]};
415+
const float deltaTanLambda{std::abs(currentTracklet.tanLambda - nextTracklet.tanLambda)};
416+
417+
if (deltaTanLambda / mTrkParams[iteration].CellDeltaTanLambdaSigma < mTrkParams[iteration].NSigmaCut) {
418+
419+
/// Track seed preparation. Clusters are numbered progressively from the innermost going outward.
420+
const int clusId[3]{
421+
mTimeFrame->getClusters()[iLayer][currentTracklet.firstClusterIndex].clusterId,
422+
mTimeFrame->getClusters()[iLayer + 1][nextTracklet.firstClusterIndex].clusterId,
423+
mTimeFrame->getClusters()[iLayer + 2][nextTracklet.secondClusterIndex].clusterId};
424+
const auto& cluster1_glo = mTimeFrame->getUnsortedClusters()[iLayer].at(clusId[0]);
425+
const auto& cluster2_glo = mTimeFrame->getUnsortedClusters()[iLayer + 1].at(clusId[1]);
426+
const auto& cluster3_tf = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer + 2).at(clusId[2]);
427+
auto track{buildTrackSeed(cluster1_glo, cluster2_glo, cluster3_tf)};
428+
429+
float chi2{0.f};
430+
bool good{false};
431+
for (int iC{2}; iC--;) {
432+
const TrackingFrameInfo& trackingHit = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer + iC).at(clusId[iC]);
433+
434+
if (!track.rotate(trackingHit.alphaTrackingFrame)) {
435+
break;
436+
}
437+
438+
if (!track.propagateTo(trackingHit.xTrackingFrame, getBz())) {
439+
break;
440+
}
441+
442+
if (!track.correctForMaterial(mTrkParams[0].LayerxX0[iLayer + iC], mTrkParams[0].LayerxX0[iLayer] * radl * rho, true)) {
443+
break;
444+
}
445+
446+
const auto predChi2{track.getPredictedChi2Quiet(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)};
447+
if (!iC && predChi2 > mTrkParams[iteration].MaxChi2ClusterAttachment) {
448+
break;
449+
}
450+
451+
if (!track.o2::track::TrackParCov::update(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)) {
452+
break;
453+
}
454+
455+
good = !iC;
456+
chi2 += predChi2;
457+
}
458+
if (good) {
459+
layerCells[position++] = CellSeed(iLayer, clusId[0], clusId[1], clusId[2], iTracklet, iNextTracklet, track, chi2);
460+
}
461+
}
371462
}
372-
mTimeFrame->getCells()[iLayer].emplace_back(iLayer, clusId[0], clusId[1], clusId[2],
373-
iTracklet, iNextTracklet, track, chi2);
374463
}
375-
}
376-
}
464+
});
465+
377466
if (iLayer > 0) {
378-
mTimeFrame->getCellsLookupTable()[iLayer - 1].resize(currentLayerTrackletsNum + 1, mTimeFrame->getCells()[iLayer].size());
467+
auto& lut = mTimeFrame->getCellsLookupTable()[iLayer - 1];
468+
lut.resize(currentLayerTrackletsNum + 1);
469+
std::copy_n(perTrackletCount.begin(), currentLayerTrackletsNum + 1, lut.begin());
379470
}
380471
}
381472
});

0 commit comments

Comments
 (0)