Project
Loading...
Searching...
No Matches
GPUTPCGMMerger.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#define GPUCA_CADEBUG 0
16#define GPUCA_MERGE_LOOPER_MC 0
17// #define GPUCA_CADEBUG_ENABLED
18
19#include "GPUCommonDef.h"
20
21#if !defined(GPUCA_GPUCODE) && (defined(GPUCA_MERGER_BY_MC_LABEL) || defined(GPUCA_CADEBUG_ENABLED) || GPUCA_MERGE_LOOPER_MC)
23#include "GPUROOTDump.h"
24#endif
25
26#ifndef GPUCA_GPUCODE_DEVICE
27#include <cstdio>
28#include <cstring>
29#include <cmath>
30#include "GPUReconstruction.h"
31#endif
32
33#include "GPUTPCTracker.h"
34#include "GPUTPCTrackParam.h"
35#include "GPUTPCGMMerger.h"
36#include "GPUO2DataTypes.h"
37#include "TPCFastTransform.h"
38#include "GPUTPCConvertImpl.h"
39#include "GPUTPCGeometry.h"
40#include "GPUDefParametersRuntime.h"
41
42#include "GPUCommonMath.h"
43#include "GPUCommonAlgorithm.h"
44#include "GPUCommonConstants.h"
45
46#include "GPUTPCTrackParam.h"
47#include "GPUTPCGMMergedTrack.h"
48#include "GPUParam.h"
50
51#include "GPUTPCGMTrackParam.h"
52#include "GPUTPCGMSectorTrack.h"
53#include "GPUTPCGMBorderTrack.h"
54
57#ifndef GPUCA_GPUCODE
60#endif
61
62using namespace o2::gpu;
63using namespace o2::tpc;
64using namespace gputpcgmmergertypes;
65
66namespace o2::gpu::internal
67{
69 float refz;
70 float x;
71 float y;
72 uint32_t id;
73};
74
76 GPUd() bool operator()(const GPUTPCGMBorderRange& a, const GPUTPCGMBorderRange& b)
77 {
78 return GPUCA_DETERMINISTIC_CODE((a.fMax != b.fMax) ? (a.fMax < b.fMax) : (a.fId < b.fId), a.fMax < b.fMax);
79 }
80};
82 GPUd() bool operator()(const GPUTPCGMBorderRange& a, const GPUTPCGMBorderRange& b)
83 {
84 return GPUCA_DETERMINISTIC_CODE((a.fMin != b.fMin) ? (a.fMin < b.fMin) : (a.fId < b.fId), a.fMin < b.fMin);
85 }
86};
87
91 GPUd() bool operator()(const int32_t aa, const int32_t bb)
92 {
95 if (a.CCE() != b.CCE()) {
96 return a.CCE() > b.CCE();
97 }
98 GPUCA_DETERMINISTIC_CODE( // clang-format off
99 if (a.NClusters() != b.NClusters()) {
100 return a.NClusters() > b.NClusters();
101 } if (CAMath::Abs(a.GetParam().GetQPt()) != CAMath::Abs(b.GetParam().GetQPt())) {
102 return CAMath::Abs(a.GetParam().GetQPt()) > CAMath::Abs(b.GetParam().GetQPt());
103 } if (a.GetParam().GetY() != b.GetParam().GetY()) {
104 return a.GetParam().GetY() > b.GetParam().GetY();
105 }
106 return aa > bb;
107 , // !GPUCA_DETERMINISTIC_CODE
108 return a.NClusters() > b.NClusters();
109 ) // clang-format on
110 }
111};
112
116 GPUd() bool operator()(const int32_t aa, const int32_t bb)
117 {
120 GPUCA_DETERMINISTIC_CODE( // clang-format off
121 if (CAMath::Abs(a.GetParam().GetQPt()) != CAMath::Abs(b.GetParam().GetQPt())) {
122 return CAMath::Abs(a.GetParam().GetQPt()) > CAMath::Abs(b.GetParam().GetQPt());
123 } if (a.GetParam().GetY() != b.GetParam().GetY()) {
124 return a.GetParam().GetY() > b.GetParam().GetY();
125 }
126 return a.GetParam().GetZ() > b.GetParam().GetZ();
127 , // !GPUCA_DETERMINISTIC_CODE
128 return CAMath::Abs(a.GetParam().GetQPt()) > CAMath::Abs(b.GetParam().GetQPt());
129 ) // clang-format on
130 }
131};
132
134 GPUd() bool operator()(const MergeLooperParam& a, const MergeLooperParam& b)
135 {
136 return GPUCA_DETERMINISTIC_CODE(CAMath::Abs(a.refz) != CAMath::Abs(b.refz) ? CAMath::Abs(a.refz) < CAMath::Abs(b.refz) : a.id < b.id, CAMath::Abs(a.refz) < CAMath::Abs(b.refz));
137 }
138};
139
140} // namespace o2::gpu::internal
141
142using namespace o2::gpu::internal;
143
144#ifndef GPUCA_GPUCODE
145
146#include "GPUQA.h"
147#include "GPUMemorySizeScalers.h"
148
150{
151 for (int32_t iSector = 0; iSector < NSECTORS; iSector++) {
152 mNextSectorInd[iSector] = iSector + 1;
153 mPrevSectorInd[iSector] = iSector - 1;
154 }
155 int32_t mid = NSECTORS / 2 - 1;
156 int32_t last = NSECTORS - 1;
157 mNextSectorInd[mid] = 0;
158 mPrevSectorInd[0] = mid;
159 mNextSectorInd[last] = NSECTORS / 2;
160 mPrevSectorInd[NSECTORS / 2] = last;
161}
162
163// DEBUG CODE
164#if !defined(GPUCA_GPUCODE) && (defined(GPUCA_MERGER_BY_MC_LABEL) || defined(GPUCA_CADEBUG_ENABLED) || GPUCA_MERGE_LOOPER_MC)
165#include "GPUQAHelper.h"
166
167template <class T>
168inline const auto* resolveMCLabels(const o2::dataformats::ConstMCTruthContainerView<o2::MCCompLabel>* a, const AliHLTTPCClusterMCLabel* b)
169{
170 return a;
171}
172template <>
173inline const auto* resolveMCLabels<AliHLTTPCClusterMCLabel>(const o2::dataformats::ConstMCTruthContainerView<o2::MCCompLabel>* a, const AliHLTTPCClusterMCLabel* b)
174{
175 return b;
176}
177
178template <class T, class S>
179int64_t GPUTPCGMMerger::GetTrackLabelA(const S& trk) const
180{
181 GPUTPCGMSectorTrack* sectorTrack = nullptr;
182 int32_t nClusters = 0;
183 if constexpr (std::is_same_v<S, GPUTPCGMBorderTrack&>) {
184 sectorTrack = &mSectorTrackInfos[trk.TrackID()];
185 nClusters = sectorTrack->OrigTrack()->NHits();
186 } else {
187 nClusters = trk.NClusters();
188 }
189 auto acc = GPUTPCTrkLbl<false, GPUTPCTrkLbl_ret>(resolveMCLabels<T>(GetConstantMem()->ioPtrs.clustersNative ? GetConstantMem()->ioPtrs.clustersNative->clustersMCTruth : nullptr, GetConstantMem()->ioPtrs.mcLabelsTPC), 0.5f);
190 for (int32_t i = 0; i < nClusters; i++) {
191 int32_t id;
192 if constexpr (std::is_same_v<S, GPUTPCGMBorderTrack&>) {
193 const GPUTPCTracker& tracker = GetConstantMem()->tpcTrackers[sectorTrack->Sector()];
194 const GPUTPCHitId& ic = tracker.TrackHits()[sectorTrack->OrigTrack()->FirstHitID() + i];
195 id = tracker.Data().ClusterDataIndex(tracker.Data().Row(ic.RowIndex()), ic.HitIndex()) + GetConstantMem()->ioPtrs.clustersNative->clusterOffset[sectorTrack->Sector()][0];
196 } else {
197 id = mClusters[trk.FirstClusterRef() + i].num;
198 }
199 acc.addLabel(id);
200 }
201 return acc.computeLabel().id;
202}
203
204template <class S>
205int64_t GPUTPCGMMerger::GetTrackLabel(const S& trk) const
206{
207#ifdef GPUCA_TPC_GEOMETRY_O2
208 if (GetConstantMem()->ioPtrs.clustersNative->clustersMCTruth) {
209 return GetTrackLabelA<o2::dataformats::ConstMCTruthContainerView<o2::MCCompLabel>, S>(trk);
210 } else
211#endif
212 {
213 return GetTrackLabelA<AliHLTTPCClusterMCLabel, S>(trk);
214 }
215}
216
217#endif
218// END DEBUG CODE
219
221{
222 uint32_t nErr = 0;
223 for (uint32_t i = 0; i < mMemory->nMergedTracks; i++) {
224 const GPUTPCGMMergedTrack& trk = mMergedTracks[i];
225 if (trk.OK()) {
226 if (trk.NClusters() == 0) {
227 GPUError("FAILURE: Track marked ok but has 0 clusters");
228 nErr++;
229 }
230 if (!trk.CCE() && !trk.MergedLooper()) {
231 const GPUTPCGMMergedTrack* updTrk = &trk;
232 while (updTrk->PrevSegment() >= 0) {
233 auto next = &mMergedTracks[updTrk->PrevSegment()];
234 if (!next->MergedLooper()) {
235 GPUError("FAILURE: prev segment not marked as merged looper\n");
236 nErr++;
237 }
238 if (next == &trk) {
239 GPUError("FAILURE: segment cycle found\n");
240 break;
241 }
242 updTrk = next;
243 }
244 if (updTrk->NClusters() == 0) {
245 printf("FAILURE: segment leg has 0 clusters");
246 }
247 }
248 }
249 }
250
251 if (nErr == 0) {
252 GPUInfo("Merged Tracks OK");
253 } else {
254 throw std::runtime_error("Error during track merging");
255 }
256}
257
259{
260 uint32_t nErr = 0;
261 std::vector<bool> trkUsed(SectorTrackInfoLocalTotal());
262 for (int32_t i = 0; i < SectorTrackInfoLocalTotal(); i++) {
263 trkUsed[i] = false;
264 }
265
266 for (int32_t itr = 0; itr < SectorTrackInfoLocalTotal(); itr++) {
267 GPUTPCGMSectorTrack& track = mSectorTrackInfos[itr];
268 if (track.PrevSegmentNeighbour() >= 0 && mSectorTrackInfos[track.PrevSegmentNeighbour()].NextSegmentNeighbour() != itr) {
269 GPUError("FAILURE: Invalid reciprocal segment link: %d PrevSegmentNeighbour %d NextSegmentNeighbour %d", itr, track.PrevSegmentNeighbour(), mSectorTrackInfos[track.PrevSegmentNeighbour()].NextSegmentNeighbour());
270 nErr++;
271 }
272 if (track.NextSegmentNeighbour() >= 0 && mSectorTrackInfos[track.NextSegmentNeighbour()].PrevSegmentNeighbour() != itr) {
273 GPUError("FAILURE: Invalid reciprocal segment link: %d NextSegmentNeighbour %d PrevSegmentNeighbour %d", itr, track.NextSegmentNeighbour(), mSectorTrackInfos[track.NextSegmentNeighbour()].PrevSegmentNeighbour());
274 nErr++;
275 }
276 if (track.PrevNeighbour() >= 0 && mSectorTrackInfos[track.PrevNeighbour()].NextNeighbour() != itr) {
277 GPUError("FAILURE: Invalid reciprocal link: %d PrevNeighbour %d NextNeighbour %d", itr, track.PrevNeighbour(), mSectorTrackInfos[track.PrevNeighbour()].NextNeighbour());
278 nErr++;
279 }
280 if (track.NextNeighbour() >= 0 && mSectorTrackInfos[track.NextNeighbour()].PrevNeighbour() != itr) {
281 GPUError("FAILURE: Invalid reciprocal link: %d NextNeighbour %d PrevNeighbour %d", itr, track.NextNeighbour(), mSectorTrackInfos[track.NextNeighbour()].PrevNeighbour());
282 nErr++;
283 }
284 if (track.PrevSegmentNeighbour() >= 0) {
285 continue;
286 }
287 if (track.PrevNeighbour() >= 0) {
288 continue;
289 }
290 GPUTPCGMSectorTrack *trbase = &track, *tr = &track;
291 while (true) {
292 int32_t iTrk = tr - mSectorTrackInfos;
293 if (trkUsed[iTrk]) {
294 GPUError("FAILURE: double use");
295 nErr++;
296 break;
297 }
298 trkUsed[iTrk] = true;
299
300 int32_t jtr = tr->NextSegmentNeighbour();
301 if (jtr >= 0) {
302 tr = &(mSectorTrackInfos[jtr]);
303 if (tr->PrevNeighbour() >= 0) {
304 GPUError("FAILURE: Non-base segment has previous leg");
305 nErr++;
306 }
307 continue;
308 }
309 jtr = trbase->NextNeighbour();
310 if (jtr >= 0) {
311 trbase = &(mSectorTrackInfos[jtr]);
312 tr = trbase;
313 if (tr->PrevSegmentNeighbour() >= 0) {
314 GPUError("FAILURE: Neibhbour leg has previous segment neightbout");
315 nErr++;
316 break;
317 }
318 continue;
319 }
320 break;
321 }
322 }
323 for (int32_t i = 0; i < SectorTrackInfoLocalTotal(); i++) {
324 if (trkUsed[i] == false) {
325 GPUError("FAILURE: trk missed");
326 nErr++;
327 }
328 }
329 if (nErr == 0) {
330 GPUInfo("Merged Track Graph OK");
331 } else {
332 throw std::runtime_error("Invalid merge graph");
333 }
334}
335
336void GPUTPCGMMerger::PrintMergeGraph(const GPUTPCGMSectorTrack* trk, std::ostream& out) const
337{
338 const GPUTPCGMSectorTrack* orgTrack = trk;
339 while (trk->PrevSegmentNeighbour() >= 0) {
340 trk = &mSectorTrackInfos[trk->PrevSegmentNeighbour()];
341 }
342 const GPUTPCGMSectorTrack* orgTower = trk;
343 while (trk->PrevNeighbour() >= 0) {
344 trk = &mSectorTrackInfos[trk->PrevNeighbour()];
345 }
346
347 int32_t nextId = trk - mSectorTrackInfos;
348 out << "Graph of track " << (orgTrack - mSectorTrackInfos) << "\n";
349 while (nextId >= 0) {
350 trk = &mSectorTrackInfos[nextId];
351 if (trk->PrevSegmentNeighbour() >= 0) {
352 out << "TRACK TREE INVALID!!! " << trk->PrevSegmentNeighbour() << " --> " << nextId << "\n";
353 }
354 out << (trk == orgTower ? "--" : " ");
355 while (nextId >= 0) {
356 GPUTPCGMSectorTrack* trk2 = &mSectorTrackInfos[nextId];
357 if (trk != trk2 && (trk2->PrevNeighbour() >= 0 || trk2->NextNeighbour() >= 0)) {
358 out << " (TRACK TREE INVALID!!! " << trk2->PrevNeighbour() << " <-- " << nextId << " --> " << trk2->NextNeighbour() << ") ";
359 }
360 char tmp[128];
361 snprintf(tmp, 128, " %s%5d(%5.2f)", trk2 == orgTrack ? "!" : " ", nextId, trk2->QPt());
362 out << tmp;
363 nextId = trk2->NextSegmentNeighbour();
364 }
365 out << "\n";
366 nextId = trk->NextNeighbour();
367 }
368}
369
371
373{
374 computePointerWithAlignment(mem, mSectorTrackInfos, mNTotalSectorTracks);
375 computePointerWithAlignment(mem, mSectorTrackInfoIndex, NSECTORS * 2 + 1);
376 if (mRec->GetProcessingSettings().deterministicGPUReconstruction) {
377 computePointerWithAlignment(mem, mTmpSortMemory, std::max(mNTotalSectorTracks, mNMaxTracks * 2));
378 }
379
380 void* memBase = mem;
381 computePointerWithAlignment(mem, mBorderMemory, 2 * mNTotalSectorTracks); // MergeBorders & Resolve
382 computePointerWithAlignment(mem, mBorderRangeMemory, 2 * mNTotalSectorTracks);
383 int32_t nTracks = 0;
384 for (int32_t iSector = 0; iSector < NSECTORS; iSector++) {
385 const int32_t n = *mRec->GetConstantMem().tpcTrackers[iSector].NTracks();
386 mBorder[iSector] = mBorderMemory + 2 * nTracks;
387 mBorder[NSECTORS + iSector] = mBorderMemory + 2 * nTracks + n;
388 mBorderRange[iSector] = mBorderRangeMemory + 2 * nTracks;
389 nTracks += n;
390 }
391 computePointerWithAlignment(mem, mTrackLinks, mNTotalSectorTracks);
392 computePointerWithAlignment(mem, mTrackCCRoots, mNTotalSectorTracks);
393 void* memMax = mem;
394 mem = memBase;
395 computePointerWithAlignment(mem, mTrackIDs, GPUCA_NSECTORS * mNMaxSingleSectorTracks); // UnpackResetIds - RefitSectorTracks - UnpackSectorGlobal
396 memMax = (void*)std::max((size_t)mem, (size_t)memMax);
397 mem = memBase;
398 computePointerWithAlignment(mem, mTrackSort, mNMaxTracks); // PrepareForFit0 - SortTracksQPt - PrepareForFit1 - PrepareForFit1 / Finalize0 - Finalize2
399 computePointerWithAlignment(mem, mSharedCount, mNMaxClusters);
400 memMax = (void*)std::max((size_t)mem, (size_t)memMax);
401 mem = memBase;
402 computePointerWithAlignment(mem, mLoopData, mNMaxTracks); // GPUTPCGMMergerTrackFit - GPUTPCGMMergerFollowLoopers
403 computePointerWithAlignment(mem, mRetryRefitIds, mNMaxTracks); // Reducing mNMaxTracks for mLoopData / mRetryRefitIds does not save memory, since the other parts are larger anyway
404 memMax = (void*)std::max((size_t)mem, (size_t)memMax);
405 mem = memBase;
406 computePointerWithAlignment(mem, mLooperCandidates, mNMaxLooperMatches); // MergeLoopers 1-3
407 memMax = (void*)std::max((size_t)mem, (size_t)memMax);
408 return memMax;
409}
410
412{
413 computePointerWithAlignment(mem, mMemory);
414 return mem;
415}
416
418{
419 computePointerWithAlignment(mem, mTrackOrderAttach, mNMaxTracks);
420 const bool mergerSortTracks = mRec->GetProcessingSettings().mergerSortTracks == -1 ? mRec->getGPUParameters(mRec->GetRecoStepsGPU() & GPUDataTypes::RecoStep::TPCMerging).par_SORT_BEFORE_FIT : mRec->GetProcessingSettings().mergerSortTracks;
421 if (mergerSortTracks) {
422 computePointerWithAlignment(mem, mTrackOrderProcess, mNMaxTracks);
423 }
424 return mem;
425}
426
428{
429 computePointerWithAlignment(mem, mMergedTracks, mNMaxTracks);
430 if (mRec->GetParam().dodEdxEnabled) {
431 computePointerWithAlignment(mem, mMergedTracksdEdx, mNMaxTracks);
432 if (mRec->GetParam().rec.tpc.dEdxClusterRejectionFlagMask != mRec->GetParam().rec.tpc.dEdxClusterRejectionFlagMaskAlt) {
433 computePointerWithAlignment(mem, mMergedTracksdEdxAlt, mNMaxTracks);
434 }
435 }
436 computePointerWithAlignment(mem, mClusters, mNMaxMergedTrackClusters);
437 computePointerWithAlignment(mem, mClusterAttachment, mNMaxClusters);
438 return mem;
439}
440
442{
443 if ((mRec->GetRecoSteps() & GPUDataTypes::RecoStep::Refit) || mRec->GetProcessingSettings().outputSharedClusterMap) {
444 computePointerWithAlignment(mem, mClusterStateExt, mNMaxClusters);
445 } else {
446 mClusterStateExt = nullptr;
447 }
448 return mem;
449}
450
452{
453 computePointerWithAlignment(mem, mOutputTracksTPCO2, HostProcessor(this).NOutputTracksTPCO2());
454 return mem;
455}
456
458{
459 computePointerWithAlignment(mem, mOutputClusRefsTPCO2, HostProcessor(this).NOutputClusRefsTPCO2());
460 return mem;
461}
462
464{
465 computePointerWithAlignment(mem, mOutputTracksTPCO2MC, HostProcessor(this).NOutputTracksTPCO2());
466 return mem;
467}
468
470{
471 computePointerWithAlignment(mem, mTrackSortO2, mNMaxTracks);
472 computePointerWithAlignment(mem, mClusRefTmp, mNMaxTracks);
473 return mem;
474}
475
477{
483 if (mRec->GetProcessingSettings().createO2Output) {
487 if (mRec->GetProcessingSettings().runMC) {
489 }
490 }
492}
493
495{
496 mNTotalSectorTracks = 0;
497 mNClusters = 0;
498 mNMaxSingleSectorTracks = 0;
499 for (int32_t iSector = 0; iSector < NSECTORS; iSector++) {
500 uint32_t ntrk = *mRec->GetConstantMem().tpcTrackers[iSector].NTracks();
501 mNTotalSectorTracks += ntrk;
502 mNClusters += *mRec->GetConstantMem().tpcTrackers[iSector].NTrackHits();
503 if (mNMaxSingleSectorTracks < ntrk) {
504 mNMaxSingleSectorTracks = ntrk;
505 }
506 }
507 mNMaxMergedTrackClusters = mRec->MemoryScalers()->NTPCMergedTrackHits(mNClusters);
508 if (CAMath::Abs(Param().polynomialField.GetNominalBz()) < (gpu_common_constants::kZeroFieldCut * gpu_common_constants::kCLight)) {
509 mNMaxTracks = mRec->MemoryScalers()->getValue(mNTotalSectorTracks, mNTotalSectorTracks); // 0 magnetic field
510 } else {
511 mNMaxTracks = mRec->MemoryScalers()->NTPCMergedTracks(mNTotalSectorTracks);
512 }
513 if (io.clustersNative) {
514 mNMaxClusters = io.clustersNative->nClustersTotal;
516 mNMaxClusters = 0;
517 for (int32_t i = 0; i < NSECTORS; i++) {
518 mNMaxClusters += mRec->GetConstantMem().tpcTrackers[i].NHitsTotal();
519 }
520 } else {
521 mNMaxClusters = mNClusters;
522 }
523 mNMaxLooperMatches = mNMaxClusters / 4; // We have that much scratch memory anyway
524}
525
527{
528 for (int32_t i = 0; i < NSECTORS; i++) {
529 if (mRec->GetConstantMem().tpcTrackers[i].CommonMemory()->nLocalTracks > (int32_t)mNMaxSingleSectorTracks) {
530 throw std::runtime_error("mNMaxSingleSectorTracks too small");
531 }
532 }
534 throw std::runtime_error("Must run also sector tracking");
535 }
536 return 0;
537}
538
539#endif // GPUCA_GPUCODE
540
541GPUd() void GPUTPCGMMerger::ClearTrackLinks(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, bool output)
542{
543 const int32_t n = output ? mMemory->nMergedTracks : SectorTrackInfoLocalTotal();
544 for (int32_t i = iBlock * nThreads + iThread; i < n; i += nThreads * nBlocks) {
545 mTrackLinks[i] = -1;
546 }
547}
548
549GPUd() int32_t GPUTPCGMMerger::RefitSectorTrack(GPUTPCGMSectorTrack& sectorTrack, const GPUTPCTrack* inTrack, float alpha, int32_t sector)
550{
552 prop.SetMaterialTPC();
553 prop.SetMaxSinPhi(GPUCA_MAX_SIN_PHI);
554 prop.SetSeedingErrors(true); // Larger errors for seeds, better since we don't start with good hypothesis
555 prop.SetFitInProjections(false);
556 prop.SetPolynomialField(&Param().polynomialField);
558 trk.X() = inTrack->Param().GetX();
559 trk.Y() = inTrack->Param().GetY();
560 trk.Z() = inTrack->Param().GetZ();
561 trk.SinPhi() = inTrack->Param().GetSinPhi();
562 trk.DzDs() = inTrack->Param().GetDzDs();
563 trk.QPt() = inTrack->Param().GetQPt();
564 trk.TOffset() = GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->convZOffsetToVertexTime(sector, inTrack->Param().GetZOffset(), Param().continuousMaxTimeBin);
565 const auto tmp = sectorTrack.ClusterTN() > sectorTrack.ClusterT0() ? std::array<float, 2>{sectorTrack.ClusterTN(), sectorTrack.ClusterT0()} : std::array<float, 2>{sectorTrack.ClusterT0(), sectorTrack.ClusterTN()};
566 trk.ShiftZ(this, sector, tmp[0], tmp[1], inTrack->Param().GetX()); // We do not store the inner / outer cluster X, so we just use the track X instead
567 sectorTrack.SetX2(0.f);
568 for (int32_t way = 0; way < 2; way++) {
569 if (way) {
570 prop.SetFitInProjections(true);
571 prop.SetPropagateBzOnly(true);
572 }
573 trk.ResetCovariance();
574 prop.SetTrack(&trk, alpha);
575 int32_t start = way ? inTrack->NHits() - 1 : 0;
576 int32_t end = way ? 0 : (inTrack->NHits() - 1);
577 int32_t incr = way ? -1 : 1;
578 for (int32_t i = start; i != end; i += incr) {
579 float x, y, z;
580 int32_t row, flags;
581 const GPUTPCTracker& tracker = GetConstantMem()->tpcTrackers[sector];
582 const GPUTPCHitId& ic = tracker.TrackHits()[inTrack->FirstHitID() + i];
583 int32_t clusterIndex = tracker.Data().ClusterDataIndex(tracker.Data().Row(ic.RowIndex()), ic.HitIndex());
584 row = ic.RowIndex();
585 const ClusterNative& cl = GetConstantMem()->ioPtrs.clustersNative->clustersLinear[GetConstantMem()->ioPtrs.clustersNative->clusterOffset[sector][0] + clusterIndex];
586 flags = cl.getFlags();
587 GetConstantMem()->calibObjects.fastTransformHelper->Transform(sector, row, cl.getPad(), cl.getTime(), x, y, z, trk.TOffset());
588 if (prop.PropagateToXAlpha(x, alpha, way == 0)) {
589 return way == 0;
590 }
591 trk.ConstrainSinPhi();
592 if (prop.Update(y, z, row, Param(), flags & GPUTPCGMMergedTrackHit::clustererAndSharedFlags, 0, false, sector, -1.f, 0.f, 0.f)) { // TODO: Use correct time / avgCharge
593 return way == 0;
594 }
595 trk.ConstrainSinPhi();
596 }
597 if (way) {
598 sectorTrack.SetParam2(trk);
599 } else {
600 sectorTrack.Set(trk, inTrack, alpha, sector);
601 }
602 }
603 return 0;
604}
605
606GPUd() void GPUTPCGMMerger::SetTrackClusterT(GPUTPCGMSectorTrack& track, int32_t iSector, const GPUTPCTrack* sectorTr)
607{
608 const GPUTPCTracker& trk = GetConstantMem()->tpcTrackers[iSector];
609 const GPUTPCHitId& ic1 = trk.TrackHits()[sectorTr->FirstHitID()];
610 const GPUTPCHitId& ic2 = trk.TrackHits()[sectorTr->FirstHitID() + sectorTr->NHits() - 1];
611 int32_t clusterIndex1 = trk.Data().ClusterDataIndex(trk.Data().Row(ic1.RowIndex()), ic1.HitIndex());
612 int32_t clusterIndex2 = trk.Data().ClusterDataIndex(trk.Data().Row(ic2.RowIndex()), ic2.HitIndex());
613 const ClusterNative* cl = GetConstantMem()->ioPtrs.clustersNative->clustersLinear + GetConstantMem()->ioPtrs.clustersNative->clusterOffset[iSector][0];
614 track.SetClusterT(cl[clusterIndex1].getTime(), cl[clusterIndex2].getTime());
615}
616
617GPUd() void GPUTPCGMMerger::UnpackSaveNumber(int32_t id)
618{
619 mSectorTrackInfoIndex[id] = mMemory->nUnpackedTracks;
620}
621
622GPUd() void GPUTPCGMMerger::UnpackSectorGlobal(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, int32_t iSector)
623{
624 const GPUTPCTracker& trk = GetConstantMem()->tpcTrackers[iSector];
625 float alpha = Param().Alpha(iSector);
626 const GPUTPCTrack* sectorTr = mMemory->firstExtrapolatedTracks[iSector];
627 uint32_t nLocalTracks = trk.CommonMemory()->nLocalTracks;
628 uint32_t nTracks = *trk.NTracks();
629 for (uint32_t itr = nLocalTracks + iBlock * nThreads + iThread; itr < nTracks; itr += nBlocks * nThreads) {
630 sectorTr = &trk.Tracks()[itr];
631 int32_t localId = mTrackIDs[((sectorTr->LocalTrackId() >> 24) & 0x3F) * mNMaxSingleSectorTracks + (sectorTr->LocalTrackId() & 0xFFFFFF)];
632 if (localId == -1) {
633 continue;
634 }
635 uint32_t myTrack = CAMath::AtomicAdd(&mMemory->nUnpackedTracks, 1u);
636 GPUTPCGMSectorTrack& track = mSectorTrackInfos[myTrack];
637 SetTrackClusterT(track, iSector, sectorTr);
638 track.Set(this, sectorTr, alpha, iSector);
639 track.SetGlobalSectorTrackCov();
640 track.SetPrevNeighbour(-1);
641 track.SetNextNeighbour(-1);
642 track.SetNextSegmentNeighbour(-1);
643 track.SetPrevSegmentNeighbour(-1);
644 track.SetLocalTrackId(localId | (sectorTr->LocalTrackId() & 0x40000000));
645 }
646}
647
648GPUd() void GPUTPCGMMerger::UnpackResetIds(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, int32_t iSector)
649{
650 const GPUTPCTracker& trk = GetConstantMem()->tpcTrackers[iSector];
651 uint32_t nLocalTracks = trk.CommonMemory()->nLocalTracks;
652 for (uint32_t i = iBlock * nThreads + iThread; i < nLocalTracks; i += nBlocks * nThreads) {
653 mTrackIDs[iSector * mNMaxSingleSectorTracks + i] = -1;
654 }
655}
656
657GPUd() void GPUTPCGMMerger::RefitSectorTracks(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, int32_t iSector)
658{
659 const GPUTPCTracker& trk = GetConstantMem()->tpcTrackers[iSector];
660 uint32_t nLocalTracks = trk.CommonMemory()->nLocalTracks;
661
662 float alpha = Param().Alpha(iSector);
663 const GPUTPCTrack* sectorTr = nullptr;
664
665 for (uint32_t itr = iBlock * nThreads + iThread; itr < nLocalTracks; itr += nBlocks * nThreads) {
666 sectorTr = &trk.Tracks()[itr];
668 SetTrackClusterT(track, iSector, sectorTr);
669 if (RefitSectorTrack(track, sectorTr, alpha, iSector)) {
670 track.Set(this, sectorTr, alpha, iSector); // TODO: Why does the refit fail, it shouldn't, this workaround should be removed
671 if (!track.FilterErrors(this, iSector, GPUCA_MAX_SIN_PHI, 0.1f)) {
672 continue;
673 }
674 }
675
676 CADEBUG(GPUInfo("INPUT Sector %d, Track %u, QPt %f DzDs %f", iSector, itr, track.QPt(), track.DzDs()));
677 track.SetPrevNeighbour(-1);
678 track.SetNextNeighbour(-1);
679 track.SetNextSegmentNeighbour(-1);
680 track.SetPrevSegmentNeighbour(-1);
681 track.SetExtrapolatedTrackId(0, -1);
682 track.SetExtrapolatedTrackId(1, -1);
683 uint32_t myTrack = CAMath::AtomicAdd(&mMemory->nUnpackedTracks, 1u);
684 mTrackIDs[iSector * mNMaxSingleSectorTracks + sectorTr->LocalTrackId()] = myTrack;
685 mSectorTrackInfos[myTrack] = track;
686 }
687}
688
689GPUd() void GPUTPCGMMerger::LinkExtrapolatedTracks(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread)
690{
691 for (int32_t itr = SectorTrackInfoGlobalFirst(0) + iBlock * nThreads + iThread; itr < SectorTrackInfoGlobalLast(NSECTORS - 1); itr += nThreads * nBlocks) {
692 GPUTPCGMSectorTrack& extrapolatedTrack = mSectorTrackInfos[itr];
693 GPUTPCGMSectorTrack& localTrack = mSectorTrackInfos[extrapolatedTrack.LocalTrackId() & 0xFFFFFF];
694 int up = (extrapolatedTrack.LocalTrackId() & 0x40000000) ? 1 : 0;
695 localTrack.SetExtrapolatedTrackId(up, itr);
696 }
697}
698
699GPUd() void GPUTPCGMMerger::MergeSectorsPrepareStep2(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, int32_t iBorder, GPUTPCGMBorderTrack** B, GPUAtomic(uint32_t) * nB, bool useOrigTrackParam)
700{
701 //* prepare sector tracks for merging with next/previous/same sector
702 //* each track transported to the border line
703
704 float fieldBz = Param().bzCLight;
705
706 float dAlpha = Param().dAlpha / 2;
707 float x0 = 0;
708
709 if (iBorder == 0) { // transport to the left edge of the sector and rotate horizontally
710 dAlpha = dAlpha - CAMath::Pi() / 2;
711 } else if (iBorder == 1) { // transport to the right edge of the sector and rotate horizontally
712 dAlpha = -dAlpha - CAMath::Pi() / 2;
713 } else if (iBorder == 2) { // transport to the middle of the sector and rotate vertically to the border on the left
714 x0 = GPUTPCGeometry::Row2X(63);
715 } else if (iBorder == 3) { // transport to the middle of the sector and rotate vertically to the border on the right
716 dAlpha = -dAlpha;
717 x0 = GPUTPCGeometry::Row2X(63);
718 } else if (iBorder == 4) { // transport to the middle of the sßector, w/o rotation
719 dAlpha = 0;
720 x0 = GPUTPCGeometry::Row2X(63);
721 }
722
723 const float maxSin = CAMath::Sin(60.f / 180.f * CAMath::Pi());
724 float cosAlpha = CAMath::Cos(dAlpha);
725 float sinAlpha = CAMath::Sin(dAlpha);
726
727 GPUTPCGMSectorTrack trackTmp;
728 for (int32_t itr = iBlock * nThreads + iThread; itr < SectorTrackInfoLocalTotal(); itr += nThreads * nBlocks) {
729 const GPUTPCGMSectorTrack* track = &mSectorTrackInfos[itr];
730 int32_t iSector = track->Sector();
731
732 if (track->PrevSegmentNeighbour() >= 0 && track->Sector() == mSectorTrackInfos[track->PrevSegmentNeighbour()].Sector()) {
733 continue;
734 }
735 if (useOrigTrackParam) { // TODO: Check how far this makes sense with sector track refit
736 if (CAMath::Abs(track->QPt()) * Param().qptB5Scaler < Param().rec.tpc.mergerLooperQPtB5Limit) {
737 continue;
738 }
739 const GPUTPCGMSectorTrack* trackMin = track;
740 while (track->NextSegmentNeighbour() >= 0 && track->Sector() == mSectorTrackInfos[track->NextSegmentNeighbour()].Sector()) {
741 track = &mSectorTrackInfos[track->NextSegmentNeighbour()];
742 if (track->OrigTrack()->Param().X() < trackMin->OrigTrack()->Param().X()) {
743 trackMin = track;
744 }
745 }
746 trackTmp = *trackMin;
747 track = &trackTmp;
748 if (trackTmp.X2() != 0.f) {
749 trackTmp.UseParam2();
750 } else {
751 trackTmp.Set(this, trackMin->OrigTrack(), trackMin->Alpha(), trackMin->Sector());
752 }
753 } else {
754 if (CAMath::Abs(track->QPt()) * Param().qptB5Scaler < Param().rec.tpc.mergerLooperSecondHorizontalQPtB5Limit) {
755 if (iBorder == 0 && track->NextNeighbour() >= 0) {
756 continue;
757 }
758 if (iBorder == 1 && track->PrevNeighbour() >= 0) {
759 continue;
760 }
761 }
762 }
764
765 if (track->TransportToXAlpha(this, x0, sinAlpha, cosAlpha, fieldBz, b, maxSin)) {
766 b.SetTrackID(itr);
767 b.SetNClusters(track->NClusters());
768 for (int32_t i = 0; i < 4; i++) {
769 if (CAMath::Abs(b.Cov()[i]) >= 5.0f) {
770 b.SetCov(i, 5.0f);
771 }
772 }
773 if (CAMath::Abs(b.Cov()[4]) >= 0.5f) {
774 b.SetCov(4, 0.5f);
775 }
776 uint32_t myTrack = CAMath::AtomicAdd(&nB[iSector], 1u);
777 B[iSector][myTrack] = b;
778 }
779 }
780}
781
782template <>
783GPUd() void GPUTPCGMMerger::MergeBorderTracks<0>(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, int32_t iSector1, GPUTPCGMBorderTrack* B1, int32_t N1, int32_t iSector2, GPUTPCGMBorderTrack* B2, int32_t N2, int32_t mergeMode)
784{
785 CADEBUG(GPUInfo("\nMERGING Sectors %d %d NTracks %d %d CROSS %d", iSector1, iSector2, N1, N2, mergeMode));
786 GPUTPCGMBorderRange* range1 = mBorderRange[iSector1];
787 GPUTPCGMBorderRange* range2 = mBorderRange[iSector2] + *GetConstantMem()->tpcTrackers[iSector2].NTracks();
788 bool sameSector = (iSector1 == iSector2);
789 for (int32_t itr = iBlock * nThreads + iThread; itr < N1; itr += nThreads * nBlocks) {
790 GPUTPCGMBorderTrack& b = B1[itr];
791 float d = CAMath::Max(0.5f, 3.5f * CAMath::Sqrt(b.Cov()[1]));
792 if (CAMath::Abs(b.Par()[4]) * Param().qptB5Scaler >= 20) {
793 d *= 2;
794 } else if (d > 3) {
795 d = 3;
796 }
797 CADEBUG(printf(" Input Sector 1 %d Track %d: ", iSector1, itr); for (int32_t i = 0; i < 5; i++) { printf("%8.3f ", b.Par()[i]); } printf(" - "); for (int32_t i = 0; i < 5; i++) { printf("%8.3f ", b.Cov()[i]); } printf(" - D %8.3f\n", d));
799 range.fId = itr;
800 range.fMin = b.Par()[1] + b.ZOffsetLinear() - d;
801 range.fMax = b.Par()[1] + b.ZOffsetLinear() + d;
802 range1[itr] = range;
803 if (sameSector) {
804 range2[itr] = range;
805 }
806 }
807 if (!sameSector) {
808 for (int32_t itr = iBlock * nThreads + iThread; itr < N2; itr += nThreads * nBlocks) {
809 GPUTPCGMBorderTrack& b = B2[itr];
810 float d = CAMath::Max(0.5f, 3.5f * CAMath::Sqrt(b.Cov()[1]));
811 if (CAMath::Abs(b.Par()[4]) * Param().qptB5Scaler >= 20) {
812 d *= 2;
813 } else if (d > 3) {
814 d = 3;
815 }
816 CADEBUG(printf(" Input Sector 2 %d Track %d: ", iSector2, itr); for (int32_t i = 0; i < 5; i++) { printf("%8.3f ", b.Par()[i]); } printf(" - "); for (int32_t i = 0; i < 5; i++) { printf("%8.3f ", b.Cov()[i]); } printf(" - D %8.3f\n", d));
818 range.fId = itr;
819 range.fMin = b.Par()[1] + b.ZOffsetLinear() - d;
820 range.fMax = b.Par()[1] + b.ZOffsetLinear() + d;
821 range2[itr] = range;
822 }
823 }
824}
825
826template <>
827GPUd() void GPUTPCGMMerger::MergeBorderTracks<1>(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, int32_t iSector1, GPUTPCGMBorderTrack* B1, int32_t N1, int32_t iSector2, GPUTPCGMBorderTrack* B2, int32_t N2, int32_t mergeMode)
828{
829#if !defined(GPUCA_GPUCODE_COMPILEKERNELS)
830 GPUTPCGMBorderRange* range1 = mBorderRange[iSector1];
831 GPUTPCGMBorderRange* range2 = mBorderRange[iSector2] + *GetConstantMem()->tpcTrackers[iSector2].NTracks();
832
833 if (iThread == 0) {
834 if (iBlock == 0) {
835 GPUCommonAlgorithm::sortDeviceDynamic(range1, range1 + N1, [](const GPUTPCGMBorderRange& a, const GPUTPCGMBorderRange& b) { return GPUCA_DETERMINISTIC_CODE((a.fMin != b.fMin) ? (a.fMin < b.fMin) : (a.fId < b.fId), a.fMin < b.fMin); });
836 } else if (iBlock == 1) {
837 GPUCommonAlgorithm::sortDeviceDynamic(range2, range2 + N2, [](const GPUTPCGMBorderRange& a, const GPUTPCGMBorderRange& b) { return GPUCA_DETERMINISTIC_CODE((a.fMax != b.fMax) ? (a.fMax < b.fMax) : (a.fId < b.fId), a.fMax < b.fMax); });
838 }
839 }
840#else
841 printf("This sorting variant is disabled for RTC");
842#endif
843}
844
845template <>
846GPUd() void GPUTPCGMMerger::MergeBorderTracks<3>(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, GPUTPCGMBorderRange* range, int32_t N, int32_t cmpMax)
847{
848#ifndef GPUCA_SPECIALIZE_THRUST_SORTS
849 if (iThread == 0 && iBlock == 0) {
850 if (cmpMax) {
851 GPUCommonAlgorithm::sortDeviceDynamic(range, range + N, MergeBorderTracks_compMax());
852 } else {
853 GPUCommonAlgorithm::sortDeviceDynamic(range, range + N, MergeBorderTracks_compMin());
854 }
855 }
856#endif
857}
858
859template <>
860GPUd() void GPUTPCGMMerger::MergeBorderTracks<2>(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, int32_t iSector1, GPUTPCGMBorderTrack* B1, int32_t N1, int32_t iSector2, GPUTPCGMBorderTrack* B2, int32_t N2, int32_t mergeMode)
861{
862 // int32_t statAll = 0, statMerged = 0;
863 float factor2ys = Param().rec.tpc.trackMergerFactor2YS;
864 float factor2zt = Param().rec.tpc.trackMergerFactor2ZT;
865 float factor2k = Param().rec.tpc.trackMergerFactor2K;
866 float factor2General = Param().rec.tpc.trackMergerFactor2General;
867
868 factor2k = factor2General * factor2k;
869 factor2ys = factor2General * factor2ys;
870 factor2zt = factor2General * factor2zt;
871
872 int32_t minNPartHits = Param().rec.tpc.trackMergerMinPartHits;
873 int32_t minNTotalHits = Param().rec.tpc.trackMergerMinTotalHits;
874
875 bool sameSector = (iSector1 == iSector2);
876
877 GPUTPCGMBorderRange* range1 = mBorderRange[iSector1];
878 GPUTPCGMBorderRange* range2 = mBorderRange[iSector2] + *GetConstantMem()->tpcTrackers[iSector2].NTracks();
879
880 int32_t i2 = 0;
881 for (int32_t i1 = iBlock * nThreads + iThread; i1 < N1; i1 += nThreads * nBlocks) {
882 GPUTPCGMBorderRange r1 = range1[i1];
883 while (i2 < N2 && range2[i2].fMax < r1.fMin) {
884 i2++;
885 }
886
887 GPUTPCGMBorderTrack& b1 = B1[r1.fId];
888 if (b1.NClusters() < minNPartHits) {
889 continue;
890 }
891 int32_t iBest2 = -1;
892 int32_t lBest2 = 0;
893 // statAll++;
894 for (int32_t k2 = i2; k2 < N2; k2++) {
895 GPUTPCGMBorderRange r2 = range2[k2];
896 if (r2.fMin > r1.fMax) {
897 break;
898 }
899 if (sameSector && (r1.fId >= r2.fId)) {
900 continue;
901 }
902 // do check
903
904 GPUTPCGMBorderTrack& b2 = B2[r2.fId];
905#if defined(GPUCA_MERGER_BY_MC_LABEL) && !defined(GPUCA_GPUCODE)
906 int64_t label1 = GetTrackLabel(b1);
907 int64_t label2 = GetTrackLabel(b2);
908 if (label1 != label2 && label1 != -1) // DEBUG CODE, match by MC label
909#endif
910 {
911 CADEBUG(if (GetConstantMem()->ioPtrs.mcLabelsTPC) {printf("Comparing track %3d to %3d: ", r1.fId, r2.fId); for (int32_t i = 0; i < 5; i++) { printf("%8.3f ", b1.Par()[i]); } printf(" - "); for (int32_t i = 0; i < 5; i++) { printf("%8.3f ", b1.Cov()[i]); } printf("\n%28s", ""); });
912 CADEBUG(if (GetConstantMem()->ioPtrs.mcLabelsTPC) {for (int32_t i = 0; i < 5; i++) { printf("%8.3f ", b2.Par()[i]); } printf(" - "); for (int32_t i = 0; i < 5; i++) { printf("%8.3f ", b2.Cov()[i]); } printf(" - %5s - ", GetTrackLabel(b1) == GetTrackLabel(b2) ? "CLONE" : "FAKE"); });
913 if (b2.NClusters() < lBest2) {
914 CADEBUG2(continue, printf("!NCl1\n"));
915 }
916 if (mergeMode > 0) {
917 // Merging CE tracks
918 int32_t maxRowDiff = mergeMode == 2 ? 1 : 3; // TODO: check cut
919 if (CAMath::Abs(b1.Row() - b2.Row()) > maxRowDiff) {
920 CADEBUG2(continue, printf("!ROW\n"));
921 }
922 if (CAMath::Abs(b1.Par()[2] - b2.Par()[2]) > 0.5f || CAMath::Abs(b1.Par()[3] - b2.Par()[3]) > 0.5f) {
923 CADEBUG2(continue, printf("!CE SinPhi/Tgl\n")); // Crude cut to avoid totally wrong matches, TODO: check cut
924 }
925 }
926
927 GPUCA_DEBUG_STREAMER_CHECK(float weight = b1.Par()[4] * b1.Par()[4]; if (o2::utils::DebugStreamer::checkStream(o2::utils::StreamFlags::streamMergeBorderTracksAll, b1.TrackID(), weight)) { MergedTrackStreamer(b1, b2, "merge_all_tracks", iSector1, iSector2, mergeMode, weight, o2::utils::DebugStreamer::getSamplingFrequency(o2::utils::StreamFlags::streamMergeBorderTracksAll)); });
928
929 if (!b1.CheckChi2Y(b2, factor2ys)) {
930 CADEBUG2(continue, printf("!Y\n"));
931 }
932 // if( !b1.CheckChi2Z(b2, factor2zt ) ) CADEBUG2(continue, printf("!NCl1\n"));
933 if (!b1.CheckChi2QPt(b2, factor2k)) {
934 CADEBUG2(continue, printf("!QPt\n"));
935 }
936 float fys = CAMath::Abs(b1.Par()[4]) * Param().qptB5Scaler < 20 ? factor2ys : (2.f * factor2ys);
937 float fzt = CAMath::Abs(b1.Par()[4]) * Param().qptB5Scaler < 20 ? factor2zt : (2.f * factor2zt);
938 if (!b1.CheckChi2YS(b2, fys)) {
939 CADEBUG2(continue, printf("!YS\n"));
940 }
941 if (!b1.CheckChi2ZT(b2, fzt)) {
942 CADEBUG2(continue, printf("!ZT\n"));
943 }
944 if (CAMath::Abs(b1.Par()[4]) * Param().qptB5Scaler < 20) {
945 if (b2.NClusters() < minNPartHits) {
946 CADEBUG2(continue, printf("!NCl2\n"));
947 }
948 if (b1.NClusters() + b2.NClusters() < minNTotalHits) {
949 CADEBUG2(continue, printf("!NCl3\n"));
950 }
951 }
952 CADEBUG(printf("OK: dZ %8.3f D1 %8.3f D2 %8.3f\n", CAMath::Abs(b1.Par()[1] - b2.Par()[1]), 3.5f * sqrt(b1.Cov()[1]), 3.5f * sqrt(b2.Cov()[1])));
953 } // DEBUG CODE, match by MC label
954 lBest2 = b2.NClusters();
955 iBest2 = b2.TrackID();
956 }
957
958 if (iBest2 < 0) {
959 continue;
960 }
961 GPUCA_DEBUG_STREAMER_CHECK(float weight = b1.Par()[4] * b1.Par()[4]; if (o2::utils::DebugStreamer::checkStream(o2::utils::StreamFlags::streamMergeBorderTracksBest, b1.TrackID(), weight)) { MergedTrackStreamer(b1, MergedTrackStreamerFindBorderTrack(B2, N2, iBest2), "merge_best_track", iSector1, iSector2, mergeMode, weight, o2::utils::DebugStreamer::getSamplingFrequency(o2::utils::StreamFlags::streamMergeBorderTracksBest)); });
962
963 // statMerged++;
964
965 CADEBUG(GPUInfo("Found match %d %d", b1.TrackID(), iBest2));
966
967 mTrackLinks[b1.TrackID()] = iBest2;
968 if (mergeMode > 0) {
969 GPUCA_DETERMINISTIC_CODE(CAMath::AtomicMax(&mTrackLinks[iBest2], b1.TrackID()), mTrackLinks[iBest2] = b1.TrackID());
970 }
971 }
972 // GPUInfo("STAT: sectors %d, %d: all %d merged %d", iSector1, iSector2, statAll, statMerged);
973}
974
975GPUdii() void GPUTPCGMMerger::MergeBorderTracksSetup(int32_t& n1, int32_t& n2, GPUTPCGMBorderTrack*& b1, GPUTPCGMBorderTrack*& b2, int32_t& jSector, int32_t iSector, int8_t withinSector, int8_t mergeMode) const
976{
977 if (withinSector == 1) { // Merge tracks within the same sector
978 jSector = iSector;
979 n1 = n2 = mMemory->tmpCounter[iSector];
980 b1 = b2 = mBorder[iSector];
981 } else if (withinSector == -1) { // Merge tracks accross the central electrode
982 jSector = (iSector + NSECTORS / 2);
983 const int32_t offset = mergeMode == 2 ? NSECTORS : 0;
984 n1 = mMemory->tmpCounter[iSector + offset];
985 n2 = mMemory->tmpCounter[jSector + offset];
986 b1 = mBorder[iSector + offset];
987 b2 = mBorder[jSector + offset];
988 } else { // Merge tracks of adjacent sectors
989 jSector = mNextSectorInd[iSector];
990 n1 = mMemory->tmpCounter[iSector];
991 n2 = mMemory->tmpCounter[NSECTORS + jSector];
992 b1 = mBorder[iSector];
993 b2 = mBorder[NSECTORS + jSector];
994 }
995}
996
997template <int32_t I>
998GPUd() void GPUTPCGMMerger::MergeBorderTracks(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, int32_t iSector, int8_t withinSector, int8_t mergeMode)
999{
1000 int32_t n1, n2;
1001 GPUTPCGMBorderTrack *b1, *b2;
1002 int32_t jSector;
1003 MergeBorderTracksSetup(n1, n2, b1, b2, jSector, iSector, withinSector, mergeMode);
1004 MergeBorderTracks<I>(nBlocks, nThreads, iBlock, iThread, iSector, b1, n1, jSector, b2, n2, mergeMode);
1005}
1006
1007#if !defined(GPUCA_GPUCODE) || defined(GPUCA_GPUCODE_DEVICE) // FIXME: DR: WORKAROUND to avoid CUDA bug creating host symbols for device code.
1008template GPUdni() void GPUTPCGMMerger::MergeBorderTracks<0>(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, int32_t iSector, int8_t withinSector, int8_t mergeMode);
1009template GPUdni() void GPUTPCGMMerger::MergeBorderTracks<1>(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, int32_t iSector, int8_t withinSector, int8_t mergeMode);
1010template GPUdni() void GPUTPCGMMerger::MergeBorderTracks<2>(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, int32_t iSector, int8_t withinSector, int8_t mergeMode);
1011#endif
1012
1013GPUd() void GPUTPCGMMerger::MergeWithinSectorsPrepare(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread)
1014{
1015 float x0 = GPUTPCGeometry::Row2X(63);
1016 const float maxSin = CAMath::Sin(60.f / 180.f * CAMath::Pi());
1017
1018 for (int32_t itr = iBlock * nThreads + iThread; itr < SectorTrackInfoLocalTotal(); itr += nThreads * nBlocks) {
1019 GPUTPCGMSectorTrack& track = mSectorTrackInfos[itr];
1020 int32_t iSector = track.Sector();
1022 if (track.TransportToX(this, x0, Param().bzCLight, b, maxSin)) {
1023 b.SetTrackID(itr);
1024 CADEBUG(printf("WITHIN SECTOR %d Track %d - ", iSector, itr); for (int32_t i = 0; i < 5; i++) { printf("%8.3f ", b.Par()[i]); } printf(" - "); for (int32_t i = 0; i < 5; i++) { printf("%8.3f ", b.Cov()[i]); } printf("\n"));
1025 b.SetNClusters(track.NClusters());
1026 uint32_t myTrack = CAMath::AtomicAdd(&mMemory->tmpCounter[iSector], 1u);
1027 mBorder[iSector][myTrack] = b;
1028 }
1029 }
1030}
1031
1032GPUd() void GPUTPCGMMerger::MergeSectorsPrepare(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, int32_t border0, int32_t border1, int8_t useOrigTrackParam)
1033{
1034 bool part2 = iBlock & 1;
1035 int32_t border = part2 ? border1 : border0;
1036 GPUAtomic(uint32_t)* n = mMemory->tmpCounter;
1037 GPUTPCGMBorderTrack** b = mBorder;
1038 if (part2) {
1039 n += NSECTORS;
1040 b += NSECTORS;
1041 }
1042 MergeSectorsPrepareStep2((nBlocks + !part2) >> 1, nThreads, iBlock >> 1, iThread, border, b, n, useOrigTrackParam);
1043}
1044
1045GPUdi() void GPUTPCGMMerger::setBlockRange(int32_t elems, int32_t nBlocks, int32_t iBlock, int32_t& start, int32_t& end)
1046{
1047 start = (elems + nBlocks - 1) / nBlocks * iBlock;
1048 end = (elems + nBlocks - 1) / nBlocks * (iBlock + 1);
1049 end = CAMath::Min(elems, end);
1050}
1051
1052GPUd() void GPUTPCGMMerger::hookEdge(int32_t u, int32_t v)
1053{
1054 if (v < 0) {
1055 return;
1056 }
1057 while (true) {
1058 u = mTrackCCRoots[u];
1059 v = mTrackCCRoots[v];
1060 if (u == v) {
1061 break;
1062 }
1063 int32_t h = CAMath::Max(u, v);
1064 int32_t l = CAMath::Min(u, v);
1065
1066 int32_t old = CAMath::AtomicCAS(&mTrackCCRoots[h], h, l);
1067 if (old == h) {
1068 break;
1069 }
1070
1071 u = mTrackCCRoots[h];
1072 v = l;
1073 }
1074}
1075
1076GPUd() void GPUTPCGMMerger::ResolveFindConnectedComponentsSetup(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread)
1077{
1078 int32_t start, end;
1079 setBlockRange(SectorTrackInfoLocalTotal(), nBlocks, iBlock, start, end);
1080 for (int32_t i = start + iThread; i < end; i += nThreads) {
1081 mTrackCCRoots[i] = i;
1082 }
1083}
1084
1085GPUd() void GPUTPCGMMerger::ResolveFindConnectedComponentsHookLinks(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread)
1086{
1087 // Compute connected components in parallel, step 1.
1088 // Source: Adaptive Work-Efficient Connected Components on the GPU, Sutton et al, 2016 (https://arxiv.org/pdf/1612.01178.pdf)
1089 int32_t start, end;
1090 setBlockRange(SectorTrackInfoLocalTotal(), nBlocks, iBlock, start, end);
1091 for (int32_t itr = start + iThread; itr < end; itr += nThreads) {
1092 hookEdge(itr, mTrackLinks[itr]);
1093 }
1094}
1095
1096GPUd() void GPUTPCGMMerger::ResolveFindConnectedComponentsHookNeighbors(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread)
1097{
1098 // Compute connected components in parallel, step 1 - Part 2.
1099 nBlocks = nBlocks / 4 * 4;
1100 if (iBlock >= nBlocks) {
1101 return;
1102 }
1103
1104 int32_t start, end;
1105 setBlockRange(SectorTrackInfoLocalTotal(), nBlocks / 4, iBlock / 4, start, end);
1106
1107 int32_t myNeighbor = iBlock % 4;
1108
1109 for (int32_t itr = start + iThread; itr < end; itr += nThreads) {
1110 int32_t v = mSectorTrackInfos[itr].AnyNeighbour(myNeighbor);
1111 hookEdge(itr, v);
1112 }
1113}
1114
1115GPUd() void GPUTPCGMMerger::ResolveFindConnectedComponentsMultiJump(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread)
1116{
1117 // Compute connected components in parallel, step 2.
1118 int32_t start, end;
1119 setBlockRange(SectorTrackInfoLocalTotal(), nBlocks, iBlock, start, end);
1120 for (int32_t itr = start + iThread; itr < end; itr += nThreads) {
1121 int32_t root = itr;
1122 int32_t next = mTrackCCRoots[root];
1123 if (root == next) {
1124 continue;
1125 }
1126 do {
1127 root = next;
1128 next = mTrackCCRoots[next];
1129 } while (root != next);
1130 mTrackCCRoots[itr] = root;
1131 }
1132}
1133
1134GPUd() void GPUTPCGMMerger::ResolveMergeSectors(GPUResolveSharedMemory& smem, int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, int8_t useOrigTrackParam, int8_t mergeAll)
1135{
1136 if (!mergeAll) {
1137 /*int32_t neighborType = useOrigTrackParam ? 1 : 0;
1138 int32_t old1 = newTrack2.PrevNeighbour(0);
1139 int32_t old2 = newTrack1.NextNeighbour(0);
1140 if (old1 < 0 && old2 < 0) neighborType = 0;
1141 if (old1 == itr) continue;
1142 if (neighborType) old1 = newTrack2.PrevNeighbour(1);
1143 if ( old1 >= 0 )
1144 {
1145 GPUTPCGMSectorTrack &oldTrack1 = mSectorTrackInfos[old1];
1146 if ( oldTrack1.NClusters() < newTrack1.NClusters() ) {
1147 newTrack2.SetPrevNeighbour( -1, neighborType );
1148 oldTrack1.SetNextNeighbour( -1, neighborType );
1149 } else continue;
1150 }
1151
1152 if (old2 == itr2) continue;
1153 if (neighborType) old2 = newTrack1.NextNeighbour(1);
1154 if ( old2 >= 0 )
1155 {
1156 GPUTPCGMSectorTrack &oldTrack2 = mSectorTrackInfos[old2];
1157 if ( oldTrack2.NClusters() < newTrack2.NClusters() )
1158 {
1159 oldTrack2.SetPrevNeighbour( -1, neighborType );
1160 } else continue;
1161 }
1162 newTrack1.SetNextNeighbour( itr2, neighborType );
1163 newTrack2.SetPrevNeighbour( itr, neighborType );*/
1164 }
1165
1166 int32_t start, end;
1167 setBlockRange(SectorTrackInfoLocalTotal(), nBlocks, iBlock, start, end);
1168
1169 for (int32_t baseIdx = 0; baseIdx < SectorTrackInfoLocalTotal(); baseIdx += nThreads) {
1170 int32_t itr = baseIdx + iThread;
1171 bool inRange = itr < SectorTrackInfoLocalTotal();
1172
1173 int32_t itr2 = -1;
1174 if (inRange) {
1175 itr2 = mTrackLinks[itr];
1176 }
1177
1178 bool resolveSector = (itr2 > -1);
1179 if (resolveSector) {
1180 int32_t root = mTrackCCRoots[itr];
1181 resolveSector &= (start <= root) && (root < end);
1182 }
1183
1184 int16_t smemIdx = work_group_scan_inclusive_add(int16_t(resolveSector));
1185
1186 if (resolveSector) {
1187 smem.iTrack1[smemIdx - 1] = itr;
1188 smem.iTrack2[smemIdx - 1] = itr2;
1189 }
1190 GPUbarrier();
1191
1192 if (iThread < nThreads - 1) {
1193 continue;
1194 }
1195
1196 const int32_t nSectors = smemIdx;
1197
1198 for (int32_t i = 0; i < nSectors; i++) {
1199 itr = smem.iTrack1[i];
1200 itr2 = smem.iTrack2[i];
1201
1202 GPUTPCGMSectorTrack* track1 = &mSectorTrackInfos[itr];
1203 GPUTPCGMSectorTrack* track2 = &mSectorTrackInfos[itr2];
1204 GPUTPCGMSectorTrack* track1Base = track1;
1205 GPUTPCGMSectorTrack* track2Base = track2;
1206
1207 bool sameSegment = CAMath::Abs(track1->NClusters() > track2->NClusters() ? track1->QPt() : track2->QPt()) * Param().qptB5Scaler < 2 || track1->QPt() * track2->QPt() > 0;
1208 // GPUInfo("\nMerge %d with %d - same segment %d", itr, itr2, (int32_t) sameSegment);
1209 // PrintMergeGraph(track1, std::cout);
1210 // PrintMergeGraph(track2, std::cout);
1211
1212 while (track2->PrevSegmentNeighbour() >= 0) {
1213 track2 = &mSectorTrackInfos[track2->PrevSegmentNeighbour()];
1214 }
1215 if (sameSegment) {
1216 if (track1 == track2) {
1217 continue;
1218 }
1219 while (track1->PrevSegmentNeighbour() >= 0) {
1220 track1 = &mSectorTrackInfos[track1->PrevSegmentNeighbour()];
1221 if (track1 == track2) {
1222 goto NextTrack;
1223 }
1224 }
1225 GPUCommonAlgorithm::swap(track1, track1Base);
1226 for (int32_t k = 0; k < 2; k++) {
1227 GPUTPCGMSectorTrack* tmp = track1Base;
1228 while (tmp->Neighbour(k) >= 0) {
1229 tmp = &mSectorTrackInfos[tmp->Neighbour(k)];
1230 if (tmp == track2) {
1231 goto NextTrack;
1232 }
1233 }
1234 }
1235
1236 while (track1->NextSegmentNeighbour() >= 0) {
1237 track1 = &mSectorTrackInfos[track1->NextSegmentNeighbour()];
1238 if (track1 == track2) {
1239 goto NextTrack;
1240 }
1241 }
1242 } else {
1243 while (track1->PrevSegmentNeighbour() >= 0) {
1244 track1 = &mSectorTrackInfos[track1->PrevSegmentNeighbour()];
1245 }
1246
1247 if (track1 == track2) {
1248 continue;
1249 }
1250 for (int32_t k = 0; k < 2; k++) {
1251 GPUTPCGMSectorTrack* tmp = track1;
1252 while (tmp->Neighbour(k) >= 0) {
1253 tmp = &mSectorTrackInfos[tmp->Neighbour(k)];
1254 if (tmp == track2) {
1255 goto NextTrack;
1256 }
1257 }
1258 }
1259
1260 float z1min, z1max, z2min, z2max;
1261 z1min = track1->MinClusterT();
1262 z1max = track1->MaxClusterT();
1263 z2min = track2->MinClusterT();
1264 z2max = track2->MaxClusterT();
1265 if (track1 != track1Base) {
1266 z1min = CAMath::Min(z1min, track1Base->MinClusterT());
1267 z1max = CAMath::Max(z1max, track1Base->MaxClusterT());
1268 }
1269 if (track2 != track2Base) {
1270 z2min = CAMath::Min(z2min, track2Base->MinClusterT());
1271 z2max = CAMath::Max(z2max, track2Base->MaxClusterT());
1272 }
1273 bool goUp = z2max - z1min > z1max - z2min;
1274
1275 if (track1->Neighbour(goUp) < 0 && track2->Neighbour(!goUp) < 0) {
1276 track1->SetNeighbor(track2 - mSectorTrackInfos, goUp);
1277 track2->SetNeighbor(track1 - mSectorTrackInfos, !goUp);
1278 // GPUInfo("Result (simple neighbor)");
1279 // PrintMergeGraph(track1, std::cout);
1280 continue;
1281 } else if (track1->Neighbour(goUp) < 0) {
1282 track2 = &mSectorTrackInfos[track2->Neighbour(!goUp)];
1283 GPUCommonAlgorithm::swap(track1, track2);
1284 } else if (track2->Neighbour(!goUp) < 0) {
1285 track1 = &mSectorTrackInfos[track1->Neighbour(goUp)];
1286 } else { // Both would work, but we use the simpler one
1287 track1 = &mSectorTrackInfos[track1->Neighbour(goUp)];
1288 }
1289 track1Base = track1;
1290 }
1291
1292 track2Base = track2;
1293 if (!sameSegment) {
1294 while (track1->NextSegmentNeighbour() >= 0) {
1295 track1 = &mSectorTrackInfos[track1->NextSegmentNeighbour()];
1296 }
1297 }
1298 track1->SetNextSegmentNeighbour(track2 - mSectorTrackInfos);
1299 track2->SetPrevSegmentNeighbour(track1 - mSectorTrackInfos);
1300 // k = 0: Merge right side
1301 // k = 1: Merge left side
1302 for (int32_t k = 0; k < 2; k++) {
1303 track1 = track1Base;
1304 track2 = track2Base;
1305 while (track2->Neighbour(k) >= 0) {
1306 if (track1->Neighbour(k) >= 0) {
1307 GPUTPCGMSectorTrack* track1new = &mSectorTrackInfos[track1->Neighbour(k)];
1308 GPUTPCGMSectorTrack* track2new = &mSectorTrackInfos[track2->Neighbour(k)];
1309 track2->SetNeighbor(-1, k);
1310 track2new->SetNeighbor(-1, k ^ 1);
1311 track1 = track1new;
1312 while (track1->NextSegmentNeighbour() >= 0) {
1313 track1 = &mSectorTrackInfos[track1->NextSegmentNeighbour()];
1314 }
1315 track1->SetNextSegmentNeighbour(track2new - mSectorTrackInfos);
1316 track2new->SetPrevSegmentNeighbour(track1 - mSectorTrackInfos);
1317 track1 = track1new;
1318 track2 = track2new;
1319 } else {
1320 GPUTPCGMSectorTrack* track2new = &mSectorTrackInfos[track2->Neighbour(k)];
1321 track1->SetNeighbor(track2->Neighbour(k), k);
1322 track2->SetNeighbor(-1, k);
1323 track2new->SetNeighbor(track1 - mSectorTrackInfos, k ^ 1);
1324 }
1325 }
1326 }
1327 // GPUInfo("Result");
1328 // PrintMergeGraph(track1, std::cout);
1329 NextTrack:;
1330 }
1331 }
1332}
1333
1334GPUd() void GPUTPCGMMerger::MergeCEFill(const GPUTPCGMSectorTrack* track, const GPUTPCGMMergedTrackHit& cls, int32_t itr)
1335{
1336 if (Param().rec.tpc.mergerCERowLimit > 0 && CAMath::Abs(track->QPt()) * Param().qptB5Scaler < 0.3f && (cls.row < Param().rec.tpc.mergerCERowLimit || cls.row >= GPUCA_ROW_COUNT - Param().rec.tpc.mergerCERowLimit)) {
1337 return;
1338 }
1339
1340 float z = 0;
1341 {
1342 float x, y;
1343 auto& cln = mConstantMem->ioPtrs.clustersNative->clustersLinear[cls.num];
1344 GPUTPCConvertImpl::convert(*mConstantMem, cls.sector, cls.row, cln.getPad(), cln.getTime(), x, y, z);
1345 }
1346
1347 if (!Param().par.continuousTracking && CAMath::Abs(z) > 10) {
1348 return;
1349 }
1350 int32_t sector = track->Sector();
1351 for (int32_t attempt = 0; attempt < 2; attempt++) {
1353 const float x0 = GPUTPCGeometry::Row2X(attempt == 0 ? 63 : cls.row);
1354 if (track->TransportToX(this, x0, Param().bzCLight, b, GPUCA_MAX_SIN_PHI_LOW)) {
1355 b.SetTrackID(itr);
1356 b.SetNClusters(mMergedTracks[itr].NClusters());
1357 if (CAMath::Abs(b.Cov()[4]) >= 0.5f) {
1358 b.SetCov(4, 0.5f); // TODO: Is this needed and better than the cut in BorderTrack?
1359 }
1360 if (track->CSide()) {
1361 b.SetPar(1, b.Par()[1] - 2 * (z - b.ZOffsetLinear()));
1362 b.SetZOffsetLinear(-b.ZOffsetLinear());
1363 }
1364 b.SetRow(cls.row);
1365 uint32_t id = sector + attempt * NSECTORS;
1366 uint32_t myTrack = CAMath::AtomicAdd(&mMemory->tmpCounter[id], 1u);
1367 mBorder[id][myTrack] = b;
1368 break;
1369 }
1370 }
1371}
1372
1373GPUd() void GPUTPCGMMerger::MergeCE(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread)
1374{
1375 const ClusterNative* cls = mConstantMem->ioPtrs.clustersNative->clustersLinear;
1376 for (uint32_t i = iBlock * nThreads + iThread; i < mMemory->nMergedTracks; i += nThreads * nBlocks) {
1377 if (mMergedTracks[i].CSide() == 0 && mTrackLinks[i] >= 0) {
1378 if (mTrackLinks[mTrackLinks[i]] != (int32_t)i) {
1379 continue;
1380 }
1381 GPUTPCGMMergedTrack* trk[2] = {&mMergedTracks[i], &mMergedTracks[mTrackLinks[i]]};
1382
1383 if (!trk[1]->OK() || trk[1]->CCE()) {
1384 continue;
1385 }
1386 bool celooper = (trk[0]->GetParam().GetQPt() * Param().qptB5Scaler > 1 && trk[0]->GetParam().GetQPt() * trk[1]->GetParam().GetQPt() < 0);
1387 celooper |= trk[0]->PrevSegment() != -1 && trk[1]->PrevSegment() != -1;
1388 if (!celooper && trk[0]->GetParam().GetPar(3) * trk[1]->GetParam().GetPar(3) < 0) {
1389 continue;
1390 }
1391
1392 bool needswap = false;
1393 if (trk[0]->PrevSegment() == -1 && trk[1]->PrevSegment() >= 0) {
1394 needswap = true;
1395 } else if (celooper) {
1396 const float z0max = -CAMath::Min(cls[mClusters[trk[0]->FirstClusterRef()].num].getTime(), cls[mClusters[trk[0]->FirstClusterRef() + trk[0]->NClusters() - 1].num].getTime());
1397 const float z1max = -CAMath::Min(cls[mClusters[trk[1]->FirstClusterRef()].num].getTime(), cls[mClusters[trk[1]->FirstClusterRef() + trk[1]->NClusters() - 1].num].getTime());
1398 if (z1max < z0max) {
1399 needswap = true;
1400 }
1401 } else {
1402 if (mClusters[trk[0]->FirstClusterRef()].row > mClusters[trk[1]->FirstClusterRef()].row) {
1403 needswap = true;
1404 }
1405 }
1406 if (needswap) {
1407 GPUCommonAlgorithm::swap(trk[0], trk[1]);
1408 }
1409
1410 if (Param().par.continuousTracking) {
1411 GPUTPCGMMergedTrackHit* clsmax;
1412 const float tmax = CAMath::MaxWithRef(cls[mClusters[trk[0]->FirstClusterRef()].num].getTime(), cls[mClusters[trk[0]->FirstClusterRef() + trk[0]->NClusters() - 1].num].getTime(),
1413 cls[mClusters[trk[1]->FirstClusterRef()].num].getTime(), cls[mClusters[trk[1]->FirstClusterRef() + trk[1]->NClusters() - 1].num].getTime(),
1414 &mClusters[trk[0]->FirstClusterRef()], &mClusters[trk[0]->FirstClusterRef() + trk[0]->NClusters() - 1],
1415 &mClusters[trk[1]->FirstClusterRef()], &mClusters[trk[1]->FirstClusterRef() + trk[1]->NClusters() - 1], clsmax);
1416 const float offset = CAMath::Max(tmax - mConstantMem->calibObjects.fastTransformHelper->getCorrMap()->getMaxDriftTime(clsmax->sector, clsmax->row, cls[clsmax->num].getPad()), 0.f);
1417 trk[1]->Param().Z() += mConstantMem->calibObjects.fastTransformHelper->getCorrMap()->convDeltaTimeToDeltaZinTimeFrame(trk[1]->CSide() * NSECTORS / 2, trk[1]->Param().TOffset() - offset);
1418 trk[1]->Param().TOffset() = offset;
1419 if (celooper) {
1420 trk[0]->Param().Z() += mConstantMem->calibObjects.fastTransformHelper->getCorrMap()->convDeltaTimeToDeltaZinTimeFrame(trk[0]->CSide() * NSECTORS / 2, trk[0]->Param().TOffset() - offset);
1421 trk[0]->Param().TOffset() = offset;
1422 }
1423 }
1424
1425 if (celooper) { // TODO: Need propper handling, avoid falsely flagging the primary leg as looper
1426 trk[0]->SetMergedLooperConnected(true);
1427 trk[0]->SetCCE(true);
1428 trk[0]->SetLooper(true);
1429 trk[1]->SetMergedLooperConnected(true);
1430 trk[1]->SetCCE(true);
1431 trk[1]->SetLooper(true);
1432 continue;
1433 }
1434
1435 uint32_t newRef = CAMath::AtomicAdd(&mMemory->nMergedTrackClusters, trk[0]->NClusters() + trk[1]->NClusters());
1436 if (newRef + trk[0]->NClusters() + trk[1]->NClusters() >= mNMaxMergedTrackClusters) {
1437 raiseError(GPUErrors::ERROR_MERGER_CE_HIT_OVERFLOW, newRef + trk[0]->NClusters() + trk[1]->NClusters(), mNMaxMergedTrackClusters);
1438 for (uint32_t k = newRef; k < mNMaxMergedTrackClusters; k++) {
1439 mClusters[k].num = 0;
1440 mClusters[k].state = 0;
1441 }
1442 CAMath::AtomicExch(&mMemory->nMergedTrackClusters, mNMaxMergedTrackClusters);
1443 return;
1444 }
1445
1446 int32_t pos = newRef;
1447#pragma unroll
1448 for (int32_t k = 1; k >= 0; k--) {
1449 for (uint32_t j = 0; j != trk[k]->NClusters(); j++) {
1450 mClusters[pos++] = mClusters[trk[k]->FirstClusterRef() + j];
1451 }
1452 }
1453 trk[1]->SetFirstClusterRef(newRef);
1454 trk[1]->SetNClusters(trk[0]->NClusters() + trk[1]->NClusters());
1455 if (trk[1]->NClusters() > GPUCA_MERGER_MAX_TRACK_CLUSTERS) {
1456 trk[1]->SetFirstClusterRef(trk[1]->FirstClusterRef() + trk[1]->NClusters() - GPUCA_MERGER_MAX_TRACK_CLUSTERS);
1457 trk[1]->SetNClusters(GPUCA_MERGER_MAX_TRACK_CLUSTERS);
1458 }
1459 trk[1]->SetCCE(true);
1460 trk[0]->SetNClusters(0);
1461 trk[0]->SetOK(false);
1462 }
1463 }
1464
1465 // for (int32_t i = 0;i < mMemory->nMergedTracks;i++) {if (mMergedTracks[i].CCE() == false) {mMergedTracks[i].SetNClusters(0);mMergedTracks[i].SetOK(false);}} //Remove all non-CE tracks
1466}
1467
1468namespace o2::gpu::internal
1469{
1470namespace // anonymous
1471{
1472struct GPUTPCGMMerger_CompareClusterIds {
1474 const bool revert;
1475 GPUd() GPUTPCGMMerger_CompareClusterIds(const GPUTPCGMMerger::trackCluster* cmp, bool r) : mCmp(cmp), revert(r) {}
1476 GPUd() bool operator()(const int16_t aa, const int16_t bb)
1477 {
1480 if (a.row != b.row) {
1481 return (a.row > b.row) ^ revert;
1482 }
1483 return GPUCA_DETERMINISTIC_CODE((a.id != b.id) ? (a.id > b.id) : (aa > bb), a.id > b.id);
1484 }
1485};
1486} // anonymous namespace
1487} // namespace o2::gpu::internal
1488
1489GPUd() void GPUTPCGMMerger::CollectMergedTracks(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread)
1490{
1491 static constexpr int32_t kMaxParts = 16;
1492 static constexpr int32_t kMaxClusters = GPUCA_MERGER_MAX_TRACK_CLUSTERS;
1493
1494 GPUTPCGMSectorTrack* trackParts[kMaxParts];
1495
1496 int32_t itr = iBlock * nThreads + iThread;
1497 GPUTPCGMSectorTrack* trbase = nullptr;
1498 int32_t leg = 0;
1499 int32_t lastMergedSegment = -1;
1500 bool revertSegments = false;
1501 bool revertInSegment = false;
1502 while (true) {
1503 if (trbase && !Param().rec.tpc.dropLoopers) {
1504 int32_t jtr = trbase->NextNeighbour();
1505 if (jtr >= 0) {
1506 trbase = &(mSectorTrackInfos[jtr]);
1507 if (trbase->PrevSegmentNeighbour() >= 0) {
1508 trbase = nullptr;
1509 } else {
1510 if (Param().rec.enableCyclicGraphWorkarounds) {
1511 trbase->SetPrevSegmentNeighbour(1000000001);
1512 }
1513 leg += revertSegments ? 1 : -1;
1514 }
1515 } else {
1516 trbase = nullptr;
1517 }
1518 }
1519
1520 if (trbase == nullptr) {
1521 while (itr < SectorTrackInfoLocalTotal()) {
1522 trbase = &mSectorTrackInfos[itr];
1523 if (trbase->PrevSegmentNeighbour() >= 0 || trbase->PrevNeighbour() >= 0) {
1524 itr += nThreads * nBlocks;
1525 continue;
1526 }
1527 break;
1528 }
1529 if (itr >= SectorTrackInfoLocalTotal()) {
1530 break;
1531 }
1532 revertSegments = false;
1533 revertInSegment = false;
1534 if (Param().rec.enableCyclicGraphWorkarounds) {
1535 trbase->SetPrevSegmentNeighbour(1000000000);
1536 }
1537 int32_t jtr = trbase->NextNeighbour();
1538 leg = 0;
1539 if (jtr >= 0) {
1540 int32_t lasttr = itr;
1541 while (jtr >= 0) { // --------------- count segments ---------------
1542 if (Param().rec.enableCyclicGraphWorkarounds && &mSectorTrackInfos[jtr] == trbase) {
1543 break; // Break cyclic graph
1544 }
1545 lasttr = jtr;
1546 leg++;
1547 jtr = mSectorTrackInfos[jtr].NextNeighbour();
1548 }
1549
1550 float mainT = 1e9;
1551 revertSegments = true;
1552 for (uint32_t k = 0; k < 2; k++) { // --------------- check if first or last segment is primary ---------------
1553 int32_t ichk = k ? lasttr : itr;
1554 const GPUTPCGMSectorTrack* trchk = &mSectorTrackInfos[ichk];
1555 while (true) {
1556 float t = -trchk->MinClusterT();
1557 if (t < mainT) {
1558 if (k) {
1559 revertSegments = false;
1560 break;
1561 }
1562 mainT = t;
1563 }
1564 int32_t next = trchk->NextSegmentNeighbour();
1565 if (next < 0 || (Param().rec.enableCyclicGraphWorkarounds && next == ichk)) {
1566 break; // Breaks also cycles
1567 }
1568 trchk = &mSectorTrackInfos[next];
1569 }
1570 }
1571 if (revertSegments) {
1572 leg = 0;
1573 }
1574
1575 { // --------------- find longest sector track of main segment ---------------
1576 int32_t length = 0;
1577 int32_t ichk = revertSegments ? itr : lasttr;
1578 const GPUTPCGMSectorTrack* trchk = &mSectorTrackInfos[ichk];
1579 const GPUTPCGMSectorTrack* longest = trchk;
1580 while (true) {
1581 if (trchk->OrigTrack()->NHits() > length) {
1582 longest = trchk;
1583 length = trchk->OrigTrack()->NHits();
1584 }
1585 int32_t next = trchk->NextSegmentNeighbour();
1586 if (next < 0 || (Param().rec.enableCyclicGraphWorkarounds && next == ichk)) {
1587 break; // Breaks also cycles
1588 }
1589 trchk = &mSectorTrackInfos[next];
1590 }
1591 revertInSegment = longest->ClusterT0() < longest->ClusterTN();
1592 }
1593 }
1594 lastMergedSegment = -1;
1595 itr += nThreads * nBlocks;
1596 }
1597
1598 do {
1599 int32_t nParts = 0;
1600 int32_t nHits = 0;
1601
1602 GPUTPCGMSectorTrack* tr = trbase;
1603 while (true) {
1604 if (nParts >= kMaxParts) {
1605 break;
1606 }
1607 if (nHits + tr->NClusters() > kMaxClusters) {
1608 break;
1609 }
1610 nHits += tr->NClusters();
1611
1612 trackParts[nParts++] = tr;
1613 for (int32_t i = 0; i < 2; i++) {
1614 if (tr->ExtrapolatedTrackId(i) != -1) {
1615 if (nParts >= kMaxParts) {
1616 break;
1617 }
1618 if (nHits + mSectorTrackInfos[tr->ExtrapolatedTrackId(i)].NClusters() > kMaxClusters) {
1619 break;
1620 }
1621 trackParts[nParts++] = &mSectorTrackInfos[tr->ExtrapolatedTrackId(i)];
1622 nHits += mSectorTrackInfos[tr->ExtrapolatedTrackId(i)].NClusters();
1623 }
1624 }
1625 int32_t jtr = tr->NextSegmentNeighbour();
1626 if (jtr >= 0) {
1627 tr = &(mSectorTrackInfos[jtr]);
1628 if (Param().rec.enableCyclicGraphWorkarounds) {
1629 tr->SetPrevSegmentNeighbour(1000000002);
1630 }
1631 continue;
1632 }
1633 break;
1634 }
1635
1636 // unpack and sort clusters
1637 if (nParts > 1 && (!revertInSegment ^ (leg & 1))) {
1638 GPUCommonAlgorithm::sort(trackParts, trackParts + nParts, [](const GPUTPCGMSectorTrack* a, const GPUTPCGMSectorTrack* b) {
1639 GPUCA_DETERMINISTIC_CODE( // clang-format off
1640 if (a->X() != b->X()) {
1641 return (a->X() > b->X());
1642 }
1643 if (a->Y() != b->Y()) {
1644 return (a->Y() > b->Y());
1645 }
1646 if (a->Z() != b->Z()) {
1647 return (a->Z() > b->Z());
1648 }
1649 return a->QPt() > b->QPt();
1650 , // !GPUCA_DETERMINISTIC_CODE
1651 return (a->X() > b->X());
1652 ) // clang-format on
1653 });
1654 }
1655
1656 trackCluster trackClusters[kMaxClusters];
1657 nHits = 0;
1658 for (int32_t ipart = 0; ipart < nParts; ipart++) {
1659 const GPUTPCGMSectorTrack* t = trackParts[ipart];
1660 CADEBUG(printf("Collect Track %d Part %d QPt %f DzDs %f\n", mMemory->nMergedTracks, ipart, t->QPt(), t->DzDs()));
1661 int32_t nTrackHits = t->NClusters();
1662 trackCluster* c2 = trackClusters + nHits + nTrackHits - 1;
1663 for (int32_t i = 0; i < nTrackHits; i++, c2--) {
1664 const GPUTPCTracker& trk = GetConstantMem()->tpcTrackers[t->Sector()];
1665 const GPUTPCHitId& ic = trk.TrackHits()[t->OrigTrack()->FirstHitID() + i];
1666 uint32_t id = trk.Data().ClusterDataIndex(trk.Data().Row(ic.RowIndex()), ic.HitIndex()) + GetConstantMem()->ioPtrs.clustersNative->clusterOffset[t->Sector()][0];
1667 *c2 = trackCluster{id, (uint8_t)ic.RowIndex(), t->Sector()};
1668 }
1669 nHits += nTrackHits;
1670 }
1671 if (nHits < GPUCA_TRACKLET_SELECTOR_MIN_HITS_B5(trbase->QPt() * Param().qptB5Scaler)) {
1672 break;
1673 }
1674
1675 const bool mustReverse = revertInSegment ^ (leg & 1);
1676 bool ordered = !mustReverse;
1677 if (ordered) {
1678 for (int32_t i = 1; i < nHits; i++) {
1679 if ((trackClusters[i].row > trackClusters[i - 1].row) ^ mustReverse || trackClusters[i].id == trackClusters[i - 1].id) {
1680 ordered = false;
1681 break;
1682 }
1683 }
1684 }
1685 int32_t firstTrackIndex = 0;
1686 int32_t lastTrackIndex = nParts - 1;
1687 if (ordered == 0) {
1688 int32_t nTmpHits = 0;
1689 trackCluster trackClustersUnsorted[kMaxClusters];
1690 int16_t clusterIndices[kMaxClusters];
1691 for (int32_t i = 0; i < nHits; i++) {
1692 trackClustersUnsorted[i] = trackClusters[i];
1693 clusterIndices[i] = i;
1694 }
1695
1696 GPUCommonAlgorithm::sort(clusterIndices, clusterIndices + nHits, GPUTPCGMMerger_CompareClusterIds(trackClusters, mustReverse));
1697
1698 nTmpHits = 0;
1699 firstTrackIndex = lastTrackIndex = -1;
1700 for (int32_t i = 0; i < nParts; i++) {
1701 nTmpHits += trackParts[i]->NClusters();
1702 if (nTmpHits > clusterIndices[0] && firstTrackIndex == -1) {
1703 firstTrackIndex = i;
1704 }
1705 if (nTmpHits > clusterIndices[nHits - 1] && lastTrackIndex == -1) {
1706 lastTrackIndex = i;
1707 }
1708 }
1709
1710 int32_t nFilteredHits = 0;
1711 int32_t indPrev = -1;
1712 for (int32_t i = 0; i < nHits; i++) {
1713 int32_t ind = clusterIndices[i];
1714 if (indPrev >= 0 && trackClustersUnsorted[ind].id == trackClustersUnsorted[indPrev].id) {
1715 continue;
1716 }
1717 indPrev = ind;
1718 trackClusters[nFilteredHits] = trackClustersUnsorted[ind];
1719 nFilteredHits++;
1720 }
1721 nHits = nFilteredHits;
1722 }
1723
1724 const uint32_t iMergedTrackFirstCluster = CAMath::AtomicAdd(&mMemory->nMergedTrackClusters, (uint32_t)nHits);
1725 if (iMergedTrackFirstCluster + nHits > mNMaxMergedTrackClusters) {
1726 raiseError(GPUErrors::ERROR_MERGER_HIT_OVERFLOW, iMergedTrackFirstCluster, mNMaxMergedTrackClusters);
1727 CAMath::AtomicExch(&mMemory->nMergedTrackClusters, mNMaxMergedTrackClusters);
1728 break;
1729 }
1730
1731 GPUTPCGMMergedTrackHit* const cl = mClusters + iMergedTrackFirstCluster;
1732
1733 for (int32_t i = 0; i < nHits; i++) {
1734 uint8_t state;
1735 const ClusterNative& c = GetConstantMem()->ioPtrs.clustersNative->clustersLinear[trackClusters[i].id];
1736 state = c.getFlags();
1737 cl[i].state = state & GPUTPCGMMergedTrackHit::clustererAndSharedFlags; // Only allow edge, deconvoluted, and shared flags
1738 cl[i].row = trackClusters[i].row;
1739 cl[i].num = trackClusters[i].id;
1740 cl[i].sector = trackClusters[i].sector;
1741 }
1742
1743 uint32_t iOutputTrack = CAMath::AtomicAdd(&mMemory->nMergedTracks, 1u);
1744 if (iOutputTrack >= mNMaxTracks) {
1745 raiseError(GPUErrors::ERROR_MERGER_TRACK_OVERFLOW, iOutputTrack, mNMaxTracks);
1746 CAMath::AtomicExch(&mMemory->nMergedTracks, mNMaxTracks);
1747 break;
1748 }
1749
1750 GPUTPCGMMergedTrack& mergedTrack = mMergedTracks[iOutputTrack];
1751 GPUTPCGMTrackParam& p1 = mergedTrack.Param();
1752 const GPUTPCGMSectorTrack& p2 = *trackParts[firstTrackIndex];
1753 mergedTrack.SetFlags(0);
1754 mergedTrack.SetOK(true);
1755 mergedTrack.SetLeg(leg);
1756 mergedTrack.SetLooper(leg > 0);
1757 mergedTrack.SetNClusters(nHits);
1758 mergedTrack.SetFirstClusterRef(iMergedTrackFirstCluster);
1759 mergedTrack.SetCSide(p2.CSide());
1760 mergedTrack.SetMergedLooperConnected(leg > 0);
1761 if (revertSegments) {
1762 mergedTrack.SetPrevSegment(-1);
1763 if (lastMergedSegment >= 0) {
1764 mMergedTracks[lastMergedSegment].SetPrevSegment(iOutputTrack);
1765 }
1766 } else {
1767 mergedTrack.SetPrevSegment(lastMergedSegment);
1768 }
1769 lastMergedSegment = iOutputTrack;
1770
1772 const float toX = GPUTPCGeometry::Row2X(cl[0].row);
1773 if (p2.TransportToX<2>(this, toX, Param().bzCLight, b, GPUCA_MAX_SIN_PHI, false)) {
1774 p1.X() = toX;
1775 p1.Y() = b.Par()[0];
1776 p1.Z() = b.Par()[1];
1777 p1.SinPhi() = b.Par()[2];
1778 } else {
1779 p1.X() = p2.X();
1780 p1.Y() = p2.Y();
1781 p1.Z() = p2.Z();
1782 p1.SinPhi() = p2.SinPhi();
1783 }
1784 p1.TOffset() = p2.TOffset();
1785 p1.DzDs() = p2.DzDs();
1786 p1.QPt() = p2.QPt();
1787 mergedTrack.SetAlpha(p2.Alpha());
1788 if (CAMath::Abs(Param().polynomialField.GetNominalBz()) < (gpu_common_constants::kZeroFieldCut * gpu_common_constants::kCLight)) {
1789 p1.QPt() = 100.f / Param().rec.bz0Pt10MeV;
1790 }
1791
1792 // if (nParts > 1) printf("Merged %d: QPt %f %d parts %d hits\n", mMemory->nMergedTracks, p1.QPt(), nParts, nHits);
1793
1794 /*if (GPUQA::QAAvailable() && mRec->GetQA() && mRec->GetQA()->SuppressTrack(mMemory->nMergedTracks))
1795 {
1796 mergedTrack.SetOK(0);
1797 mergedTrack.SetNClusters(0);
1798 }
1799 if (mergedTrack.NClusters() && mergedTrack.OK()) */
1800 if (leg == 0 && Param().rec.tpc.mergeCE) {
1801 auto& cls = mConstantMem->ioPtrs.clustersNative->clustersLinear;
1802 bool CEside = cls[cl[0].num].getTime() < cls[cl[nHits - 1].num].getTime();
1803 MergeCEFill(trackParts[CEside ? lastTrackIndex : firstTrackIndex], cl[CEside ? (nHits - 1) : 0], iOutputTrack);
1804 }
1805 } while (false);
1806 }
1807}
1808
1809GPUd() void GPUTPCGMMerger::SortTracksPrepare(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread)
1810{
1811 for (uint32_t i = iBlock * nThreads + iThread; i < mMemory->nMergedTracks; i += nThreads * nBlocks) {
1812 mTrackOrderProcess[i] = i;
1813 }
1814}
1815
1816GPUd() void GPUTPCGMMerger::PrepareForFit0(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread)
1817{
1818 for (uint32_t i = iBlock * nThreads + iThread; i < mMemory->nMergedTracks; i += nBlocks * nThreads) {
1819 mTrackSort[i] = i;
1820 }
1821}
1822
1823GPUd() void GPUTPCGMMerger::SortTracks(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread)
1824{
1825#ifndef GPUCA_SPECIALIZE_THRUST_SORTS
1826 if (iThread == 0 && iBlock == 0) {
1827 GPUCommonAlgorithm::sortDeviceDynamic(mTrackOrderProcess, mTrackOrderProcess + mMemory->nMergedTracks, GPUTPCGMMergerSortTracks_comp(mMergedTracks));
1828 }
1829#endif
1830}
1831
1832GPUd() void GPUTPCGMMerger::SortTracksQPt(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread)
1833{
1834#ifndef GPUCA_SPECIALIZE_THRUST_SORTS
1835 if (iThread == 0 && iBlock == 0) {
1836 GPUCommonAlgorithm::sortDeviceDynamic(mTrackSort, mTrackSort + mMemory->nMergedTracks, GPUTPCGMMergerSortTracksQPt_comp(mMergedTracks));
1837 }
1838#endif
1839}
1840
1841GPUd() void GPUTPCGMMerger::PrepareForFit1(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread)
1842{
1843 for (uint32_t i = iBlock * nThreads + iThread; i < mMemory->nMergedTracks; i += nBlocks * nThreads) {
1844 mTrackOrderAttach[mTrackSort[i]] = i;
1845 GPUTPCGMMergedTrack& trk = mMergedTracks[i];
1846 if (trk.OK()) {
1847 for (uint32_t j = 0; j < trk.NClusters(); j++) {
1848 mClusterAttachment[mClusters[trk.FirstClusterRef() + j].num] = attachAttached | attachGood;
1849 CAMath::AtomicAdd(&mSharedCount[mClusters[trk.FirstClusterRef() + j].num], 1u);
1850 }
1851 if (!trk.CCE() && !trk.MergedLooper()) {
1852 GPUTPCGMMergedTrack* updTrk = trk.GetFirstSegment(mMergedTracks, Param().rec.enableCyclicGraphWorkarounds);
1853 const auto &cl0 = mClusters[trk.FirstClusterRef()], &cln = mClusters[updTrk->FirstClusterRef() + updTrk->NClusters() - 1];
1854 const auto& GPUrestrict() cls = GetConstantMem()->ioPtrs.clustersNative->clustersLinear;
1855 float z0 = cls[cl0.num].getTime(), zn = cls[cln.num].getTime();
1856 const auto tmp = zn > z0 ? std::array<float, 3>{zn, z0, GPUTPCGeometry::Row2X(cln.row)} : std::array<float, 3>{z0, zn, GPUTPCGeometry::Row2X(cl0.row)};
1857 trk.Param().ShiftZ(this, cl0.sector, tmp[0], tmp[1], tmp[2]);
1858 updTrk = &trk;
1859 while (updTrk->PrevSegment() >= 0) {
1860 auto next = &mMergedTracks[updTrk->PrevSegment()];
1861 if (Param().rec.enableCyclicGraphWorkarounds && next == &trk) {
1862 break;
1863 }
1864 updTrk = next;
1865 updTrk->Param().TOffset() = trk.Param().TOffset();
1866 }
1867 }
1868 }
1869 }
1870}
1871
1872GPUd() void GPUTPCGMMerger::PrepareForFit2(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread)
1873{
1874 for (uint32_t i = iBlock * nThreads + iThread; i < mMemory->nMergedTrackClusters; i += nBlocks * nThreads) {
1875 if (mSharedCount[mClusters[i].num] > 1) {
1876 mClusters[i].state |= GPUTPCGMMergedTrackHit::flagShared;
1877 }
1878 }
1879 if (mClusterStateExt) {
1880 for (uint32_t i = iBlock * nThreads + iThread; i < mNMaxClusters; i += nBlocks * nThreads) {
1881 uint8_t state = GetConstantMem()->ioPtrs.clustersNative->clustersLinear[i].getFlags();
1882 if (mSharedCount[i] > 1) {
1884 }
1885 mClusterStateExt[i] = state;
1886 }
1887 }
1888}
1889
1890GPUd() void GPUTPCGMMerger::Finalize0(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread)
1891{
1892 for (uint32_t i = iBlock * nThreads + iThread; i < mMemory->nMergedTracks; i += nThreads * nBlocks) {
1893 mTrackSort[mTrackOrderAttach[i]] = i;
1894 }
1895 for (uint32_t i = iBlock * nThreads + iThread; i < mMemory->nMergedTrackClusters; i += nThreads * nBlocks) {
1896 mClusterAttachment[mClusters[i].num] = 0; // Reset adjacent attachment for attached clusters, set correctly below
1897 }
1898}
1899
1900GPUd() void GPUTPCGMMerger::Finalize1(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread)
1901{
1902 for (uint32_t i = iBlock * nThreads + iThread; i < mMemory->nMergedTracks; i += nThreads * nBlocks) {
1903 const GPUTPCGMMergedTrack& trk = mMergedTracks[i];
1904 if (!trk.OK() || trk.NClusters() == 0) {
1905 continue;
1906 }
1907 for (uint32_t j = 0; j < trk.NClusters(); j++) {
1908 int32_t id = mClusters[trk.FirstClusterRef() + j].num;
1909 uint32_t weight = mTrackOrderAttach[i] | attachAttached;
1910 uint8_t clusterState = mClusters[trk.FirstClusterRef() + j].state;
1911 if (!(clusterState & GPUTPCGMMergedTrackHit::flagReject)) {
1912 weight |= attachGood;
1913 } else if (clusterState & GPUTPCGMMergedTrackHit::flagHighIncl) {
1915 }
1916 if (trk.Leg() == 0) {
1918 }
1919 CAMath::AtomicMax(&mClusterAttachment[id], weight);
1920 }
1921 }
1922}
1923
1924GPUd() void GPUTPCGMMerger::Finalize2(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread)
1925{
1926 for (uint32_t i = iBlock * nThreads + iThread; i < mNMaxClusters; i += nThreads * nBlocks) {
1927 if (mClusterAttachment[i] != 0) {
1928 mClusterAttachment[i] = (mClusterAttachment[i] & attachFlagMask) | mTrackSort[mClusterAttachment[i] & attachTrackMask];
1929 }
1930 }
1931}
1932
1933GPUd() void GPUTPCGMMerger::MergeLoopersInit(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread)
1934{
1935 const float lowPtThresh = Param().rec.tpc.rejectQPtB5 * 1.1f; // Might need to merge tracks above the threshold with parts below the rejection threshold
1936 for (uint32_t i = get_global_id(0); i < mMemory->nMergedTracks; i += get_global_size(0)) {
1937 const auto& trk = mMergedTracks[i];
1938 const auto& p = trk.GetParam();
1939 const float qptabs = CAMath::Abs(p.GetQPt());
1940 if (trk.OK() && trk.NClusters() && trk.Leg() == 0 && qptabs * Param().qptB5Scaler > 5.f && qptabs * Param().qptB5Scaler <= lowPtThresh) {
1941 const int32_t sector = mClusters[trk.FirstClusterRef() + trk.NClusters() - 1].sector;
1942 const float refz = p.GetZ() + GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->convVertexTimeToZOffset(sector, p.GetTOffset(), Param().continuousMaxTimeBin) + (trk.CSide() ? -100 : 100);
1943 float sinA, cosA;
1944 CAMath::SinCos(trk.GetAlpha(), sinA, cosA);
1945 float gx = cosA * p.GetX() - sinA * p.GetY();
1946 float gy = cosA * p.GetY() + sinA * p.GetX();
1947 float bz = Param().polynomialField.GetFieldBz(gx, gy, p.GetZ());
1948 const float r1 = p.GetQPt() * bz;
1949 const float r = CAMath::Abs(r1) > 0.0001f ? (1.f / r1) : 10000;
1950 const float mx = p.GetX() + r * p.GetSinPhi();
1951 const float my = p.GetY() - r * CAMath::Sqrt(1 - p.GetSinPhi() * p.GetSinPhi());
1952 const float gmx = cosA * mx - sinA * my;
1953 const float gmy = cosA * my + sinA * mx;
1954 uint32_t myId = CAMath::AtomicAdd(&mMemory->nLooperMatchCandidates, 1u);
1955 if (myId >= mNMaxLooperMatches) {
1956 raiseError(GPUErrors::ERROR_LOOPER_MATCH_OVERFLOW, myId, mNMaxLooperMatches);
1957 CAMath::AtomicExch(&mMemory->nLooperMatchCandidates, mNMaxLooperMatches);
1958 return;
1959 }
1960 mLooperCandidates[myId] = MergeLooperParam{refz, gmx, gmy, i};
1961
1962 /*printf("Track %u Sanity qpt %f snp %f bz %f\n", mMemory->nLooperMatchCandidates, p.GetQPt(), p.GetSinPhi(), bz);
1963 for (uint32_t k = 0;k < trk.NClusters();k++) {
1964 float xx, yy, zz;
1965 const ClusterNative& GPUrestrict() cl = GetConstantMem()->ioPtrs.clustersNative->clustersLinear[mClusters[trk.FirstClusterRef() + k].num];
1966 GetConstantMem()->calibObjects.fastTransformHelper->Transform(mClusters[trk.FirstClusterRef() + k].sector, mClusters[trk.FirstClusterRef() + k].row, cl.getPad(), cl.getTime(), xx, yy, zz, p.GetTOffset());
1967 float sa2, ca2;
1968 CAMath::SinCos(Param().Alpha(mClusters[trk.FirstClusterRef() + k].sector), sa2, ca2);
1969 float cx = ca2 * xx - sa2 * yy;
1970 float cy = ca2 * yy + sa2 * xx;
1971 float dist = CAMath::Sqrt((cx - gmx) * (cx - gmx) + (cy - gmy) * (cy - gmy));
1972 printf("Hit %3d/%3d sector %d xy %f %f R %f\n", k, trk.NClusters(), (int32_t)mClusters[trk.FirstClusterRef() + k].sector, cx, cy, dist);
1973 }*/
1974 }
1975 }
1976}
1977
1978GPUd() void GPUTPCGMMerger::MergeLoopersSort(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread)
1979{
1980#ifndef GPUCA_SPECIALIZE_THRUST_SORTS
1981 if (iThread == 0 && iBlock == 0) {
1982 GPUCommonAlgorithm::sortDeviceDynamic(mLooperCandidates, mLooperCandidates + mMemory->nLooperMatchCandidates, GPUTPCGMMergerMergeLoopers_comp());
1983 }
1984#endif
1985}
1986
1987GPUd() void GPUTPCGMMerger::MergeLoopersMain(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread)
1988{
1989 const MergeLooperParam* candidates = mLooperCandidates;
1990
1991#if GPUCA_MERGE_LOOPER_MC && !defined(GPUCA_GPUCODE)
1992 std::vector<int64_t> paramLabels(mMemory->nLooperMatchCandidates);
1993 for (uint32_t i = 0; i < mMemory->nLooperMatchCandidates; i++) {
1994 paramLabels[i] = GetTrackLabel(mMergedTracks[candidates[i].id]);
1995 }
1996 /*std::vector<bool> dropped(mMemory->nLooperMatchCandidates);
1997 std::vector<bool> droppedMC(mMemory->nLooperMatchCandidates);
1998 std::vector<int32_t> histMatch(101);
1999 std::vector<int32_t> histFail(101);*/
2000 if (!mRec->GetProcessingSettings().runQA) {
2001 throw std::runtime_error("Need QA enabled for the Merge Loopers MC QA");
2002 }
2003#endif
2004
2005 for (uint32_t i = get_global_id(0); i < mMemory->nLooperMatchCandidates; i += get_global_size(0)) {
2006 for (uint32_t j = i + 1; j < mMemory->nLooperMatchCandidates; j++) {
2007 // int32_t bs = 0;
2008 assert(CAMath::Abs(candidates[i].refz) <= CAMath::Abs(candidates[j].refz));
2009 if (CAMath::Abs(candidates[j].refz) > CAMath::Abs(candidates[i].refz) + 100.f) {
2010 break;
2011 }
2012 const float d2xy = CAMath::Sum2(candidates[i].x - candidates[j].x, candidates[i].y - candidates[j].y);
2013 if (d2xy > 15.f) {
2014 // bs |= 1;
2015 continue;
2016 }
2017
2018 const GPUTPCGMMergedTrack* trkI = &mMergedTracks[candidates[i].id];
2019 float refZI = candidates[i].refz;
2020 {
2021 const auto* tmp = trkI->GetFirstSegment(mMergedTracks, Param().rec.enableCyclicGraphWorkarounds);
2022 if (tmp != trkI && tmp->CSide() == trkI->CSide() && CAMath::Abs(tmp->GetParam().GetZ()) > CAMath::Abs(trkI->GetParam().GetZ())) {
2023 float tmpRefZ = refZI + tmp->GetParam().GetZ() - trkI->GetParam().GetZ();
2024 if (CAMath::Abs(tmpRefZ) < CAMath::Abs(candidates[j].refz) && CAMath::Abs(tmpRefZ) > CAMath::Abs(refZI)) {
2025 trkI = tmp;
2026 refZI = tmpRefZ;
2027 }
2028 }
2029 };
2030 const auto& trk1 = *trkI;
2031 const auto& trk2 = mMergedTracks[candidates[j].id];
2032 const auto& param1 = trk1.GetParam();
2033 const auto& param2 = trk2.GetParam();
2034 if (CAMath::Abs(param1.GetDzDs()) > 0.03f && CAMath::Abs(param2.GetDzDs()) > 0.03f && param1.GetDzDs() * param2.GetDzDs() * param1.GetQPt() * param2.GetQPt() < 0) {
2035 // bs |= 2;
2036 continue;
2037 }
2038
2039 const float dznormalized = (CAMath::Abs(candidates[j].refz) - CAMath::Abs(refZI)) / (CAMath::TwoPi() * 0.5f * (CAMath::Abs(param1.GetDzDs()) + CAMath::Abs(param2.GetDzDs())) * 1.f / (0.5f * (CAMath::Abs(param1.GetQPt()) + CAMath::Abs(param2.GetQPt())) * CAMath::Abs(Param().polynomialField.GetNominalBz())));
2040 const float phasecorr = CAMath::Modf((CAMath::ASin(param1.GetSinPhi()) + trk1.GetAlpha() - CAMath::ASin(param2.GetSinPhi()) - trk2.GetAlpha()) / CAMath::TwoPi() + 5.5f, 1.f) - 0.5f;
2041 const float phasecorrdirection = (candidates[j].refz * param1.GetQPt() * param1.GetDzDs()) > 0 ? 1 : -1;
2042 const float dzcorr = dznormalized + phasecorr * phasecorrdirection;
2043 const bool sameside = !(trk1.CSide() ^ trk2.CSide());
2044 const float dzcorrlimit[4] = {sameside ? 0.018f : 0.012f, sameside ? 0.12f : 0.025f, 0.14f, 0.15f};
2045 const int32_t dzcorrcount = sameside ? 4 : 2;
2046 bool dzcorrok = false;
2047 float dznorm = 0.f;
2048 for (int32_t k = 0; k < dzcorrcount; k++) {
2049 const float d = CAMath::Abs(dzcorr - 0.5f * k);
2050 if (d <= dzcorrlimit[k]) {
2051 dzcorrok = true;
2052 dznorm = d / dzcorrlimit[k];
2053 break;
2054 }
2055 }
2056 if (!dzcorrok) {
2057 // bs |= 4;
2058 continue;
2059 }
2060
2061 const float dtgl = param1.GetDzDs() - (param1.GetQPt() * param2.GetQPt() > 0 ? param2.GetDzDs() : -param2.GetDzDs());
2062 const float dqpt = (CAMath::Abs(param1.GetQPt()) - CAMath::Abs(param2.GetQPt())) / CAMath::Min(param1.GetQPt(), param2.GetQPt());
2063 float d = CAMath::Sum2(dtgl * (1.f / 0.03f), dqpt * (1.f / 0.04f)) + d2xy * (1.f / 4.f) + dznorm * (1.f / 0.3f);
2064 bool EQ = d < 6.f;
2065#if GPUCA_MERGE_LOOPER_MC && !defined(GPUCA_GPUCODE)
2066 const int64_t label1 = paramLabels[i];
2067 const int64_t label2 = paramLabels[j];
2068 bool labelEQ = label1 != -1 && label1 == label2;
2069 if (1 || EQ || labelEQ) {
2070 // printf("Matching track %d/%d %u-%u (%ld/%ld): dist %f side %d %d, tgl %f %f, qpt %f %f, x %f %f, y %f %f\n", (int32_t)EQ, (int32_t)labelEQ, i, j, label1, label2, d, (int32_t)mMergedTracks[candidates[i].id].CSide(), (int32_t)mMergedTracks[candidates[j].id].CSide(), candidates[i].tgl, candidates[j].tgl, candidates[i].qpt, candidates[j].qpt, candidates[i].x, candidates[j].x, candidates[i].y, candidates[j].y);
2071 static auto& tup = GPUROOTDump<TNtuple>::get("mergeloopers", "labeleq:sides:d2xy:tgl1:tgl2:qpt1:qpt2:dz:dzcorr:dtgl:dqpt:dznorm:bs");
2072 tup.Fill((float)labelEQ, (trk1.CSide() ? 1 : 0) | (trk2.CSide() ? 2 : 0), d2xy, param1.GetDzDs(), param2.GetDzDs(), param1.GetQPt(), param2.GetQPt(), CAMath::Abs(candidates[j].refz) - CAMath::Abs(refZI), dzcorr, dtgl, dqpt, dznorm, bs);
2073 static auto tup2 = GPUROOTDump<TNtuple>::getNew("mergeloopers2", "labeleq:refz1:refz2:tgl1:tgl2:qpt1:qpt2:snp1:snp2:a1:a2:dzn:phasecor:phasedir:dzcorr");
2074 tup2.Fill((float)labelEQ, refZI, candidates[j].refz, param1.GetDzDs(), param2.GetDzDs(), param1.GetQPt(), param2.GetQPt(), param1.GetSinPhi(), param2.GetSinPhi(), trk1.GetAlpha(), trk2.GetAlpha(), dznormalized, phasecorr, phasecorrdirection, dzcorr);
2075 }
2076 /*if (EQ) {
2077 dropped[j] = true;
2078 }
2079 if (labelEQ) {
2080 droppedMC[j] = true;
2081 histMatch[CAMath::Min<int32_t>(100, d * 10.f)]++;
2082 }
2083 if (d < 10.f && !labelEQ) {
2084 histFail[CAMath::Min<int32_t>(100, d * 10.f)]++;
2085 }*/
2086#endif
2087 if (EQ) {
2088 mMergedTracks[candidates[j].id].SetMergedLooperUnconnected(true);
2089 if (CAMath::Abs(param2.GetQPt() * Param().qptB5Scaler) >= Param().rec.tpc.rejectQPtB5) {
2090 mMergedTracks[candidates[i].id].SetMergedLooperUnconnected(true);
2091 }
2092 }
2093 }
2094 }
2095 /*#if GPUCA_MERGE_LOOPER_MC && !defined(GPUCA_GPUCODE)
2096 int32_t total = 0, totalmc = 0, good = 0, missed = 0, fake = 0;
2097 for (uint32_t i = 0; i < mMemory->nLooperMatchCandidates; i++) {
2098 total += dropped[i];
2099 totalmc += droppedMC[i];
2100 good += dropped[i] && droppedMC[i];
2101 missed += droppedMC[i] && !dropped[i];
2102 fake += dropped[i] && !droppedMC[i];
2103 }
2104 if (good) {
2105 printf("%20s: %8d\n", "Total", total);
2106 printf("%20s: %8d\n", "TotalMC", totalmc);
2107 printf("%20s: %8d (%8.3f%% %8.3f%%)\n", "Good", good, 100.f * good / total, 100.f * good / totalmc);
2108 printf("%20s: %8d (%8.3f%%)\n", "Missed", missed, 100.f * missed / totalmc);
2109 printf("%20s: %8d (%8.3f%%)\n", "Fake", fake, 100.f * fake / total);
2110 }
2111 printf("Match histo\n");
2112 for (uint32_t i = 0; i < histMatch.size(); i++) {
2113 printf("%8.3f: %3d\n", i / 10.f + 0.05f, histMatch[i]);
2114 }
2115 printf("Fake histo\n");
2116 for (uint32_t i = 0; i < histFail.size(); i++) {
2117 printf("%8.3f: %3d\n", i / 10.f + 0.05f, histFail[i]);
2118 }
2119#endif*/
2120}
benchmark::State & state
Class of a TPC cluster in TPC-native coordinates (row, time)
A const (ready only) version of MCTruthContainer.
int32_t i
#define GPUdii()
#define get_global_size(dim)
#define GPUAtomic(type)
#define GPUbarrier()
#define GPUdni()
#define GPUrestrict()
#define get_global_id(dim)
#define GPUCA_DEBUG_STREAMER_CHECK(...)
#define GPUCA_DETERMINISTIC_CODE(det, indet)
#define GPUCA_MERGER_MAX_TRACK_CLUSTERS
#define GPUCA_MAX_SIN_PHI_LOW
#define GPUCA_TRACKLET_SELECTOR_MIN_HITS_B5(QPTB5)
#define GPUCA_MAX_SIN_PHI
#define CADEBUG(...)
Definition GPUDef.h:72
#define CADEBUG2(cmd,...)
Definition GPUDef.h:73
const GPUTPCGMMerger::trackCluster *const mCmp
const bool revert
#define GPUCA_NSECTORS
#define GPUCA_ROW_COUNT
constexpr int p2()
constexpr int p1()
constexpr to accelerate the coordinates changing
void output(const std::map< std::string, ChannelStat > &channels)
Definition rawdump.cxx:197
uint16_t pos
Definition RawData.h:3
uint32_t j
Definition RawData.h:0
uint32_t c
Definition RawData.h:2
Definition of TPCFastTransform class.
int nClusters
double num
Definition B.h:16
Class for time synchronization of RawReader instances.
GPUReconstruction * mRec
static void computePointerWithAlignment(T *&basePtr, S *&objPtr, size_t nEntries=1)
static GPUROOTDump< T, Args... > getNew(const char *name1, Names... names)
Definition GPUROOTDump.h:64
static GPUROOTDump< T, Args... > & get(const char *name1, Names... names)
Definition GPUROOTDump.h:58
RecoStepField GetRecoStepsGPU() const
int16_t RegisterMemoryAllocation(T *proc, void *(T::*setPtr)(void *), int32_t type, const char *name="", const GPUMemoryReuse &re=GPUMemoryReuse())
RecoStepField GetRecoSteps() const
const GPUParam & GetParam() const
const GPUConstantMem & GetConstantMem() const
GPUMemorySizeScalers * MemoryScalers()
const GPUSettingsProcessing & GetProcessingSettings() const
virtual const GPUDefParameters & getGPUParameters(bool doGPU) const =0
void * SetPointersRefitScratch(void *mem)
void * SetPointersMemory(void *mem)
static constexpr const int32_t NSECTORS
void * SetPointersOutput(void *mem)
void * SetPointersOutputO2MC(void *mem)
void * SetPointersOutputO2Scratch(void *mem)
void * SetPointersMerger(void *mem)
void * SetPointersOutputO2Clus(void *mem)
void * SetPointersOutputState(void *mem)
void SetMaxData(const GPUTrackingInOutPointers &io)
void * SetPointersOutputO2(void *mem)
int32_t int32_t int32_t bool float float Z
int32_t int32_t int32_t bool float Y
GLint GLint GLsizei GLint border
Definition glcorearb.h:275
GLdouble n
Definition glcorearb.h:1982
GLfloat GLfloat GLfloat alpha
Definition glcorearb.h:279
GLint GLenum GLint x
Definition glcorearb.h:403
GLuint GLuint end
Definition glcorearb.h:469
const GLdouble * v
Definition glcorearb.h:832
GLenum array
Definition glcorearb.h:4274
GLuint GLuint GLfloat weight
Definition glcorearb.h:5477
GLboolean GLboolean GLboolean b
Definition glcorearb.h:1233
GLenum GLint * range
Definition glcorearb.h:1899
GLint y
Definition glcorearb.h:270
GLintptr offset
Definition glcorearb.h:660
GLuint GLsizei GLsizei * length
Definition glcorearb.h:790
typedef void(APIENTRYP PFNGLCULLFACEPROC)(GLenum mode)
GLuint GLfloat x0
Definition glcorearb.h:5034
GLbitfield flags
Definition glcorearb.h:1570
GLboolean r
Definition glcorearb.h:1233
GLuint start
Definition glcorearb.h:469
GLboolean GLboolean GLboolean GLboolean a
Definition glcorearb.h:1233
GLuint id
Definition glcorearb.h:650
GLdouble GLdouble GLdouble z
Definition glcorearb.h:843
uint8_t itsSharedClusterMap uint8_t
Node par(int index)
Parameters.
GPUdi() o2
Definition TrackTRD.h:38
GPUd() const expr uint32_t MultivariatePolynomialHelper< Dim
if(!okForPhiMin(phi0, phi1))
Global TPC definitions and constants.
Definition SimTraits.h:167
@ streamMergeBorderTracksBest
stream MergeBorderTracks best track
@ streamMergeBorderTracksAll
stream MergeBorderTracks all tracks
Defining DataPointCompositeObject explicitly as copiable.
std::string getTime(uint64_t ts)
GPUReconstruction * rec
GPUTPCTracker tpcTrackers[GPUCA_NSECTORS]
size_t getValue(size_t maxVal, size_t val)
size_t NTPCMergedTrackHits(size_t tpcSectorTrackHitss)
size_t NTPCMergedTracks(size_t tpcSectorTracks)
const o2::tpc::ClusterNativeAccess * clustersNative
const GPUTPCGMMergedTrack & GPUrestrict() b
GPUd() bool operator()(const int32_t aa
GPUhd() GPUTPCGMMergerSortTracksQPt_comp(GPUTPCGMMergedTrack *cmp)
const GPUTPCGMMergedTrack & GPUrestrict() b
GPUhd() GPUTPCGMMergerSortTracks_comp(GPUTPCGMMergedTrack *cmp)
GPUd() bool operator()(const int32_t aa
std::vector< int > row
char const *restrict const cmp
Definition x9.h:96