Skip to content

Commit 8af91cc

Browse files
davidrohrdavidrohr
committed
GPU TPC: Refactor cluster merging
1 parent 0a2e17a commit 8af91cc

2 files changed

Lines changed: 63 additions & 45 deletions

File tree

‎GPU/GPUTracking/Merger/GPUTPCGMTrackParam.cxx‎

Lines changed: 62 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -121,26 +121,13 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger& GPUrestrict() merger, int32_
121121
uint8_t clusterState = clusters[ihit].state;
122122
const float clAlpha = param.Alpha(clusters[ihit].sector);
123123
float xx, yy, zz;
124-
{
125-
const ClusterNative& GPUrestrict() cl = merger.GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num];
126-
merger.GetConstantMem()->calibObjects.fastTransform->Transform(clusters[ihit].sector, clusters[ihit].row, cl.getPad(), cl.getTime(), xx, yy, zz, mTOffset);
127-
}
128124
CADEBUG(printf("\tHit %3d/%3d Row %3d: Cluster Alpha %8.3f %3d, X %8.3f - Y %8.3f, Z %8.3f (Missed %d)\n", ihit, maxN, (int32_t)clusters[ihit].row, clAlpha, (int32_t)clusters[ihit].sector, xx, yy, zz, nMissed));
129-
if (MergeDoubleRowClusters(ihit, wayDirection, clusters, merger, prop, xx, yy, zz, maxN, clAlpha, clusterState, !param.rec.tpc.rebuildTrackInFit && allowChangeClusters) == -1) {
125+
if (MergeDoubleRowClusters(ihit, wayDirection, clusters, merger, prop, xx, yy, zz, maxN, clAlpha, clusterState, param.rec.tpc.rebuildTrackInFit ? rebuilt : allowChangeClusters) == -1) {
130126
nMissed++;
131127
nMissed2++;
132128
continue;
133129
}
134130

135-
if (param.rec.tpc.rejectIFCLowRadiusCluster) {
136-
const float r2 = xx * xx + yy * yy;
137-
const float rmax = (83.5f + param.rec.tpc.sysClusErrorMinDist);
138-
if (r2 < rmax * rmax) {
139-
MarkClusters(clusters, ihitMergeFirst, ihit, wayDirection, GPUTPCGMMergedTrackHit::flagRejectErr);
140-
continue;
141-
}
142-
}
143-
144131
const auto& cluster = clusters[ihit];
145132
if (lastRow != 255 && CAMath::Abs(cluster.row - lastRow) > 1) {
146133
bool dodEdx = param.dodEdxEnabled && param.rec.tpc.adddEdxSubThresholdClusters && finalFit && CAMath::Abs(cluster.row - lastRow) == 2;
@@ -156,14 +143,14 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger& GPUrestrict() merger, int32_
156143
if (prop.GetPropagatedYZ(mX - GPUTPCGeometry::Row2X(iRow - step) + GPUTPCGeometry::Row2X(iRow), tmpY, tmpZ)) {
157144
break;
158145
}
159-
merger.GetConstantMem()->calibObjects.fastTransformHelper->InverseTransformYZtoX(cluster.sector, iRow, tmpY, tmpZ, tmpX);
146+
merger.GetConstantMem()->calibObjects.fastTransform->InverseTransformYZtoX(cluster.sector, iRow, tmpY, tmpZ, tmpX);
160147
if (prop.PropagateToXAlpha(tmpX, prop.GetAlpha(), inFlyDirection)) {
161148
break;
162149
}
163150
if GPUCA_RTC_CONSTEXPR (GPUCA_GET_CONSTEXPR(param.par, dodEdx)) {
164151
if (dodEdx) {
165152
float yUncorrected, zUncorrected;
166-
merger.GetConstantMem()->calibObjects.fastTransformHelper->InverseTransformYZtoNominalYZ(cluster.sector, iRow, mP[0], mP[1], yUncorrected, zUncorrected);
153+
merger.GetConstantMem()->calibObjects.fastTransform->InverseTransformYZtoNominalYZ(cluster.sector, iRow, mP[0], mP[1], yUncorrected, zUncorrected);
167154
uint32_t pad = CAMath::Float2UIntRn(GPUTPCGeometry::LinearY2Pad(cluster.sector, iRow, yUncorrected));
168155
if (!(pad >= GPUTPCGeometry::NPads(iRow) || (merger.GetConstantMem()->calibObjects.dEdxCalibContainer && merger.GetConstantMem()->calibObjects.dEdxCalibContainer->isDead(cluster.sector, iRow, pad)))) {
169156
dEdx.fillSubThreshold(iRow);
@@ -407,7 +394,7 @@ GPUdii() float GPUTPCGMTrackParam::FindBestInterpolatedHit(GPUTPCGMMerger& GPUre
407394
merger.GetConstantMem()->calibObjects.fastTransform->InverseTransformYZtoNominalYZ(sector, row, ImP0, ImP1, uncorrectedY, uncorrectedZ);
408395

409396
int32_t nCandidates = 0;
410-
while (nCandidates < param.rec.tpc.rebuildTrackInFitClusterCandidates && merger.ClusterCandidates()[(iTrk * GPUCA_ROW_COUNT + row) * param.rec.tpc.rebuildTrackInFitClusterCandidates + nCandidates].id > 1) {
397+
while (nCandidates < param.rec.tpc.rebuildTrackInFitClusterCandidates && merger.ClusterCandidates()[(iTrk * GPUTPCGeometry::NROWS + row) * param.rec.tpc.rebuildTrackInFitClusterCandidates + nCandidates].id > 1) {
411398
nCandidates++;
412399
}
413400
if (CAMath::Abs(uncorrectedY) <= rowData.getTPCMaxY()) {
@@ -569,55 +556,86 @@ GPUd() void GPUTPCGMTrackParam::MirrorTo(GPUTPCGMPropagator& GPUrestrict() prop,
569556

570557
GPUd() int32_t GPUTPCGMTrackParam::MergeDoubleRowClusters(int32_t& ihit, int32_t wayDirection, GPUTPCGMMergedTrackHit* GPUrestrict() clusters, const GPUTPCGMMerger& GPUrestrict() merger, GPUTPCGMPropagator& GPUrestrict() prop, float& GPUrestrict() xx, float& GPUrestrict() yy, float& GPUrestrict() zz, int32_t maxN, float clAlpha, uint8_t& GPUrestrict() clusterState, bool rejectChi2)
571558
{
559+
const int32_t ihitFirst = ihit;
560+
{
561+
const ClusterNative& GPUrestrict() cl = merger.GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num];
562+
merger.GetConstantMem()->calibObjects.fastTransform->Transform(clusters[ihit].sector, clusters[ihit].row, cl.getPad(), cl.getTime(), xx, yy, zz, mTOffset);
563+
}
572564
if (ihit + wayDirection >= 0 && ihit + wayDirection < maxN && clusters[ihit].row == clusters[ihit + wayDirection].row && clusters[ihit].sector == clusters[ihit + wayDirection].sector) {
573-
float maxDistY, maxDistZ;
574-
prop.GetErr2(maxDistY, maxDistZ, merger.Param(), zz, clusters[ihit].row, 0, clusters[ihit].sector, -1.f, 0.f, 0.f); // TODO: Use correct time, avgCharge
575-
maxDistY = (maxDistY + mC[0]) * 20.f;
576-
maxDistZ = (maxDistZ + mC[2]) * 20.f;
577-
int32_t noReject = 0; // Cannot reject if simple estimation of y/z fails (extremely unlike case)
565+
float maxDistY2, maxDistZ2;
566+
bool noReject = false; // Cannot reject if simple estimation of y/z fails (extremely unlike case)
578567
if (CAMath::Abs(clAlpha - prop.GetAlpha()) > 1.e-4f) {
579568
noReject = prop.RotateToAlpha(clAlpha);
580569
}
581570
float projY = 0, projZ = 0;
582-
if (noReject == 0) {
571+
if (!noReject) {
583572
noReject |= prop.GetPropagatedYZ(xx, projY, projZ);
584573
}
585-
float count = 0.f;
586-
xx = yy = zz = 0.f;
587-
clusterState = 0;
588-
while (true) {
589-
const ClusterNative& GPUrestrict() cl = merger.GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num];
590-
float clamp = cl.qTot;
591-
float clx, cly, clz;
592-
merger.GetConstantMem()->calibObjects.fastTransform->Transform(clusters[ihit].sector, clusters[ihit].row, cl.getPad(), cl.getTime(), clx, cly, clz, mTOffset);
593-
float dy = cly - projY;
594-
float dz = clz - projZ;
595-
if (noReject == 0 && (dy * dy > maxDistY || dz * dz > maxDistZ)) {
596-
CADEBUG(printf("\t\tRejecting double-row cluster: dy %f, dz %f, chiY %f, chiZ %f (Y: trk %f prj %f cl %f - Z: trk %f prj %f cl %f)\n", dy, dz, sqrtf(maxDistY), sqrtf(maxDistZ), mP[0], projY, cly, mP[1], projZ, clz));
574+
prop.GetErr2(maxDistY2, maxDistZ2, merger.Param(), zz, clusters[ihit].row, 0, clusters[ihit].sector, -1.f, 0.f, 0.f); // TODO: Use correct time, avgCharge
575+
const float kFactor = merger.GetConstantMem()->tpcTrackers[0].GetChiSeedFactor() * 4.f;
576+
maxDistY2 = (maxDistY2 + mC[0]) * kFactor;
577+
maxDistZ2 = (maxDistZ2 + mC[2]) * kFactor;
578+
auto chkFunction = [clusters, rejectChi2, maxDistY2, maxDistZ2, projY, projZ, noReject](int32_t ih, float y, float z) {
579+
float dy = y - projY;
580+
float dz = z - projZ;
581+
if (!noReject && (dy * dy > maxDistY2 || dz * dz > maxDistZ2)) {
582+
CADEBUG(printf("\t\tRejecting double-row cluster: dy %f, dz %f, chiY %f, chiZ %f (Y: trk %f prj %f cl %f - Z: trk %f prj %f cl %f)\n", dy, dz, sqrtf(maxDistY2), sqrtf(maxDistZ2), mP[0], projY, y, mP[1], projZ, z));
597583
if (rejectChi2) {
598-
clusters[ihit].state |= GPUTPCGMMergedTrackHit::flagRejectDistance;
584+
clusters[ih].state |= GPUTPCGMMergedTrackHit::flagRejectDistance;
599585
}
586+
return false;
600587
} else {
601-
CADEBUG(printf("\t\tMerging hit row %d X %f Y %f Z %f (dy %f, dz %f, chiY %f, chiZ %f)\n", clusters[ihit].row, clx, cly, clz, dy, dz, sqrtf(maxDistY), sqrtf(maxDistZ)));
602-
xx += clx * clamp; // TODO: Weight in pad/time instead of XYZ
588+
CADEBUG(printf("\t\tMerging hit row %d Y %f Z %f (dy %f, dz %f, chiY %f, chiZ %f)\n", clusters[ih].row, y, z, dy, dz, sqrtf(maxDistY2), sqrtf(maxDistZ2)));
589+
return true;
590+
}
591+
};
592+
const float tmpX = xx;
593+
float count;
594+
if (chkFunction(ihit, yy, zz)) {
595+
const ClusterNative& GPUrestrict() cl = merger.GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num];
596+
const float clamp = cl.qTot;
597+
xx *= clamp;
598+
yy *= clamp;
599+
zz *= clamp;
600+
clusterState = clusters[ihit].state;
601+
count = clamp;
602+
} else {
603+
xx = yy = zz = count = 0.f;
604+
clusterState = 0;
605+
}
606+
do {
607+
ihit += wayDirection;
608+
const ClusterNative& GPUrestrict() cl = merger.GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num];
609+
const float clamp = cl.qTot;
610+
float clx, cly, clz;
611+
merger.GetConstantMem()->calibObjects.fastTransform->Transform(clusters[ihit].sector, clusters[ihit].row, cl.getPad(), cl.getTime(), clx, cly, clz, mTOffset);
612+
if (chkFunction(ihit, cly, clz)) {
613+
xx += clx * clamp;
603614
yy += cly * clamp;
604615
zz += clz * clamp;
605616
clusterState |= clusters[ihit].state;
606617
count += clamp;
607618
}
608-
if (!(ihit + wayDirection >= 0 && ihit + wayDirection < maxN && clusters[ihit].row == clusters[ihit + wayDirection].row && clusters[ihit].sector == clusters[ihit + wayDirection].sector)) {
609-
break;
610-
}
611-
ihit += wayDirection;
612-
}
619+
} while (ihit + wayDirection >= 0 && ihit + wayDirection < maxN && clusters[ihit].row == clusters[ihit + wayDirection].row && clusters[ihit].sector == clusters[ihit + wayDirection].sector);
613620
if (count < 0.1f) {
614621
CADEBUG(printf("\t\tNo matching cluster in double-row, skipping\n"));
622+
xx = tmpX;
615623
return -1;
616624
}
617625
xx /= count;
618626
yy /= count;
619627
zz /= count;
620628
}
629+
if (merger.Param().rec.tpc.rejectIFCLowRadiusCluster) {
630+
const float r2 = xx * xx + yy * yy;
631+
const float rmax2 = CAMath::Square(83.5f + merger.Param().rec.tpc.sysClusErrorMinDist);
632+
if (r2 < rmax2) {
633+
if (rejectChi2) {
634+
MarkClusters(clusters, ihitFirst, ihit, wayDirection, GPUTPCGMMergedTrackHit::flagRejectErr);
635+
}
636+
return -1;
637+
}
638+
}
621639
return 0;
622640
}
623641

@@ -783,7 +801,7 @@ GPUdii() void GPUTPCGMTrackParam::PropagateLooper(const GPUTPCGMMerger& GPUrestr
783801
}
784802
}
785803

786-
GPUdi() void GPUTPCGMTrackParam::AttachClustersLooperFollow(const GPUTPCGMMerger& GPUrestrict() merger, GPUTPCGMPropagator& GPUrestrict() prop, int32_t sector, int32_t iRow, const int32_t iTrack, const bool up)
804+
GPUdi() void GPUTPCGMTrackParam::AttachClustersLooperFollow(const GPUTPCGMMerger& GPUrestrict() merger, GPUTPCGMPropagator& GPUrestrict() prop, int32_t sector, const int32_t iTrack, const bool up)
787805
{
788806
float toX = mX;
789807
bool inFlyDirection = (merger.MergedTracks()[iTrack].Leg() & 1) ^ up;

‎GPU/GPUTracking/Merger/GPUTPCGMTrackParam.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -166,7 +166,7 @@ class GPUTPCGMTrackParam
166166
GPUd() float AttachClusters(const GPUTPCGMMerger& GPUrestrict() merger, int32_t sector, int32_t iRow, int32_t iTrack, bool goodLeg, GPUTPCGMPropagator& prop); // Returns uncorrectedY for later use
167167
GPUd() float AttachClusters(const GPUTPCGMMerger& GPUrestrict() merger, int32_t sector, int32_t iRow, int32_t iTrack, bool goodLeg, float Y, float Z);
168168
GPUd() void AttachClustersLooper(const GPUTPCGMMerger& GPUrestrict() merger, int32_t sector, int32_t iRow, const int32_t iTrack, const bool up, const GPUTPCGMPropagator& prop);
169-
GPUd() void AttachClustersLooperFollow(const GPUTPCGMMerger& GPUrestrict() merger, GPUTPCGMPropagator& prop, int32_t sector, int32_t iRow, const int32_t iTrack, const bool up);
169+
GPUd() void AttachClustersLooperFollow(const GPUTPCGMMerger& GPUrestrict() merger, GPUTPCGMPropagator& prop, int32_t sector, const int32_t iTrack, const bool up);
170170
GPUd() void StoreLoopPropagation(const GPUTPCGMMerger& GPUrestrict() merger, int32_t sector, int32_t iRow, int32_t iTrack, bool outwards, float alpha);
171171
GPUd() void StoreOuter(gputpcgmmergertypes::GPUTPCOuterParam* outerParam, float alpha);
172172
GPUd() static void PropagateLooper(const GPUTPCGMMerger& merger, int32_t loopIdx);

0 commit comments

Comments
 (0)