Skip to content

Commit 58e662d

Browse files
davidrohrdavidrohr
committed
GPU TPC: Refactor cluster merging
1 parent 7d543d9 commit 58e662d

File tree

1 file changed

+58
-40
lines changed

1 file changed

+58
-40
lines changed

‎GPU/GPUTracking/Merger/GPUTPCGMTrackParam.cxx‎

Lines changed: 58 additions & 40 deletions
Original file line numberDiff line numberDiff line change
@@ -127,26 +127,13 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger& GPUrestrict() merger, int32_
127127
uint8_t clusterState = clusters[ihit].state;
128128
const float clAlpha = param.Alpha(clusters[ihit].sector);
129129
float xx, yy, zz;
130-
{
131-
const ClusterNative& GPUrestrict() cl = merger.GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num];
132-
merger.GetConstantMem()->calibObjects.fastTransformHelper->Transform(clusters[ihit].sector, clusters[ihit].row, cl.getPad(), cl.getTime(), xx, yy, zz, mTOffset);
133-
}
134130
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));
135-
if (MergeDoubleRowClusters(ihit, wayDirection, clusters, merger, prop, xx, yy, zz, maxN, clAlpha, clusterState, !param.rec.tpc.rebuildTrackInFit && allowChangeClusters) == -1) {
131+
if (MergeDoubleRowClusters(ihit, wayDirection, clusters, merger, prop, xx, yy, zz, maxN, clAlpha, clusterState, param.rec.tpc.rebuildTrackInFit ? rebuilt : allowChangeClusters) == -1) {
136132
nMissed++;
137133
nMissed2++;
138134
continue;
139135
}
140136

141-
if (param.rec.tpc.rejectIFCLowRadiusCluster) {
142-
const float r2 = xx * xx + yy * yy;
143-
const float rmax = (83.5f + param.rec.tpc.sysClusErrorMinDist);
144-
if (r2 < rmax * rmax) {
145-
MarkClusters(clusters, ihitMergeFirst, ihit, wayDirection, GPUTPCGMMergedTrackHit::flagRejectErr);
146-
continue;
147-
}
148-
}
149-
150137
const auto& cluster = clusters[ihit];
151138
if (lastRow != 255 && CAMath::Abs(cluster.row - lastRow) > 1) {
152139
bool dodEdx = param.dodEdxEnabled && param.rec.tpc.adddEdxSubThresholdClusters && finalFit && CAMath::Abs(cluster.row - lastRow) == 2;
@@ -575,55 +562,86 @@ GPUd() void GPUTPCGMTrackParam::MirrorTo(GPUTPCGMPropagator& GPUrestrict() prop,
575562

576563
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)
577564
{
565+
const int32_t ihitFirst = ihit;
566+
{
567+
const ClusterNative& GPUrestrict() cl = merger.GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num];
568+
merger.GetConstantMem()->calibObjects.fastTransformHelper->Transform(clusters[ihit].sector, clusters[ihit].row, cl.getPad(), cl.getTime(), xx, yy, zz, mTOffset);
569+
}
578570
if (ihit + wayDirection >= 0 && ihit + wayDirection < maxN && clusters[ihit].row == clusters[ihit + wayDirection].row && clusters[ihit].sector == clusters[ihit + wayDirection].sector) {
579-
float maxDistY, maxDistZ;
580-
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
581-
maxDistY = (maxDistY + mC[0]) * 20.f;
582-
maxDistZ = (maxDistZ + mC[2]) * 20.f;
583-
int32_t noReject = 0; // Cannot reject if simple estimation of y/z fails (extremely unlike case)
571+
float maxDistY2, maxDistZ2;
572+
bool noReject = false; // Cannot reject if simple estimation of y/z fails (extremely unlike case)
584573
if (CAMath::Abs(clAlpha - prop.GetAlpha()) > 1.e-4f) {
585574
noReject = prop.RotateToAlpha(clAlpha);
586575
}
587576
float projY = 0, projZ = 0;
588-
if (noReject == 0) {
577+
if (!noReject) {
589578
noReject |= prop.GetPropagatedYZ(xx, projY, projZ);
590579
}
591-
float count = 0.f;
592-
xx = yy = zz = 0.f;
593-
clusterState = 0;
594-
while (true) {
595-
const ClusterNative& GPUrestrict() cl = merger.GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num];
596-
float clamp = cl.qTot;
597-
float clx, cly, clz;
598-
merger.GetConstantMem()->calibObjects.fastTransformHelper->Transform(clusters[ihit].sector, clusters[ihit].row, cl.getPad(), cl.getTime(), clx, cly, clz, mTOffset);
599-
float dy = cly - projY;
600-
float dz = clz - projZ;
601-
if (noReject == 0 && (dy * dy > maxDistY || dz * dz > maxDistZ)) {
602-
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));
580+
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
581+
const float kFactor = merger.GetConstantMem()->tpcTrackers[0].GetChiSeedFactor() * 4.f;
582+
maxDistY2 = (maxDistY2 + mC[0]) * kFactor;
583+
maxDistZ2 = (maxDistZ2 + mC[2]) * kFactor;
584+
auto chkFunction = [clusters, rejectChi2, maxDistY2, maxDistZ2, projY, projZ, noReject](int32_t ih, float y, float z) {
585+
float dy = y - projY;
586+
float dz = z - projZ;
587+
if (!noReject && (dy * dy > maxDistY2 || dz * dz > maxDistZ2)) {
588+
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));
603589
if (rejectChi2) {
604-
clusters[ihit].state |= GPUTPCGMMergedTrackHit::flagRejectDistance;
590+
clusters[ih].state |= GPUTPCGMMergedTrackHit::flagRejectDistance;
605591
}
592+
return false;
606593
} else {
607-
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)));
608-
xx += clx * clamp; // TODO: Weight in pad/time instead of XYZ
594+
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)));
595+
return true;
596+
}
597+
};
598+
const float tmpX = xx;
599+
float count;
600+
if (chkFunction(ihit, yy, zz)) {
601+
const ClusterNative& GPUrestrict() cl = merger.GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num];
602+
const float clamp = cl.qTot;
603+
xx *= clamp;
604+
yy *= clamp;
605+
zz *= clamp;
606+
clusterState = clusters[ihit].state;
607+
count = clamp;
608+
} else {
609+
xx = yy = zz = count = 0.f;
610+
clusterState = 0;
611+
}
612+
do {
613+
ihit += wayDirection;
614+
const ClusterNative& GPUrestrict() cl = merger.GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num];
615+
const float clamp = cl.qTot;
616+
float clx, cly, clz;
617+
merger.GetConstantMem()->calibObjects.fastTransformHelper->Transform(clusters[ihit].sector, clusters[ihit].row, cl.getPad(), cl.getTime(), clx, cly, clz, mTOffset);
618+
if (chkFunction(ihit, cly, clz)) {
619+
xx += clx * clamp;
609620
yy += cly * clamp;
610621
zz += clz * clamp;
611622
clusterState |= clusters[ihit].state;
612623
count += clamp;
613624
}
614-
if (!(ihit + wayDirection >= 0 && ihit + wayDirection < maxN && clusters[ihit].row == clusters[ihit + wayDirection].row && clusters[ihit].sector == clusters[ihit + wayDirection].sector)) {
615-
break;
616-
}
617-
ihit += wayDirection;
618-
}
625+
} while (ihit + wayDirection >= 0 && ihit + wayDirection < maxN && clusters[ihit].row == clusters[ihit + wayDirection].row && clusters[ihit].sector == clusters[ihit + wayDirection].sector);
619626
if (count < 0.1f) {
620627
CADEBUG(printf("\t\tNo matching cluster in double-row, skipping\n"));
628+
xx = tmpX;
621629
return -1;
622630
}
623631
xx /= count;
624632
yy /= count;
625633
zz /= count;
626634
}
635+
if (merger.Param().rec.tpc.rejectIFCLowRadiusCluster) {
636+
const float r2 = xx * xx + yy * yy;
637+
const float rmax2 = CAMath::Square(83.5f + merger.Param().rec.tpc.sysClusErrorMinDist);
638+
if (r2 < rmax2) {
639+
if (rejectChi2) {
640+
MarkClusters(clusters, ihitFirst, ihit, wayDirection, GPUTPCGMMergedTrackHit::flagRejectErr);
641+
}
642+
return -1;
643+
}
644+
}
627645
return 0;
628646
}
629647

0 commit comments

Comments
 (0)