Skip to content

Commit 8889fda

Browse files
davidrohrdavidrohr
committed
GPU TPC: Refactor cluster merging
1 parent a63f272 commit 8889fda

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;
@@ -572,55 +559,86 @@ GPUd() void GPUTPCGMTrackParam::MirrorTo(GPUTPCGMPropagator& GPUrestrict() prop,
572559

573560
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)
574561
{
562+
const int32_t ihitFirst = ihit;
563+
{
564+
const ClusterNative& GPUrestrict() cl = merger.GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num];
565+
merger.GetConstantMem()->calibObjects.fastTransformHelper->Transform(clusters[ihit].sector, clusters[ihit].row, cl.getPad(), cl.getTime(), xx, yy, zz, mTOffset);
566+
}
575567
if (ihit + wayDirection >= 0 && ihit + wayDirection < maxN && clusters[ihit].row == clusters[ihit + wayDirection].row && clusters[ihit].sector == clusters[ihit + wayDirection].sector) {
576-
float maxDistY, maxDistZ;
577-
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
578-
maxDistY = (maxDistY + mC[0]) * 20.f;
579-
maxDistZ = (maxDistZ + mC[2]) * 20.f;
580-
int32_t noReject = 0; // Cannot reject if simple estimation of y/z fails (extremely unlike case)
568+
float maxDistY2, maxDistZ2;
569+
bool noReject = false; // Cannot reject if simple estimation of y/z fails (extremely unlike case)
581570
if (CAMath::Abs(clAlpha - prop.GetAlpha()) > 1.e-4f) {
582571
noReject = prop.RotateToAlpha(clAlpha);
583572
}
584573
float projY = 0, projZ = 0;
585-
if (noReject == 0) {
574+
if (!noReject) {
586575
noReject |= prop.GetPropagatedYZ(xx, projY, projZ);
587576
}
588-
float count = 0.f;
589-
xx = yy = zz = 0.f;
590-
clusterState = 0;
591-
while (true) {
592-
const ClusterNative& GPUrestrict() cl = merger.GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num];
593-
float clamp = cl.qTot;
594-
float clx, cly, clz;
595-
merger.GetConstantMem()->calibObjects.fastTransformHelper->Transform(clusters[ihit].sector, clusters[ihit].row, cl.getPad(), cl.getTime(), clx, cly, clz, mTOffset);
596-
float dy = cly - projY;
597-
float dz = clz - projZ;
598-
if (noReject == 0 && (dy * dy > maxDistY || dz * dz > maxDistZ)) {
599-
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));
577+
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
578+
const float kFactor = merger.GetConstantMem()->tpcTrackers[0].GetChiSeedFactor() * 4.f;
579+
maxDistY2 = (maxDistY2 + mC[0]) * kFactor;
580+
maxDistZ2 = (maxDistZ2 + mC[2]) * kFactor;
581+
auto chkFunction = [clusters, rejectChi2, maxDistY2, maxDistZ2, projY, projZ, noReject](int32_t ih, float y, float z) {
582+
float dy = y - projY;
583+
float dz = z - projZ;
584+
if (!noReject && (dy * dy > maxDistY2 || dz * dz > maxDistZ2)) {
585+
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));
600586
if (rejectChi2) {
601-
clusters[ihit].state |= GPUTPCGMMergedTrackHit::flagRejectDistance;
587+
clusters[ih].state |= GPUTPCGMMergedTrackHit::flagRejectDistance;
602588
}
589+
return false;
603590
} else {
604-
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)));
605-
xx += clx * clamp; // TODO: Weight in pad/time instead of XYZ
591+
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)));
592+
return true;
593+
}
594+
};
595+
const float tmpX = xx;
596+
float count;
597+
if (chkFunction(ihit, yy, zz)) {
598+
const ClusterNative& GPUrestrict() cl = merger.GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num];
599+
const float clamp = cl.qTot;
600+
xx *= clamp;
601+
yy *= clamp;
602+
zz *= clamp;
603+
clusterState = clusters[ihit].state;
604+
count = clamp;
605+
} else {
606+
xx = yy = zz = count = 0.f;
607+
clusterState = 0;
608+
}
609+
do {
610+
ihit += wayDirection;
611+
const ClusterNative& GPUrestrict() cl = merger.GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num];
612+
const float clamp = cl.qTot;
613+
float clx, cly, clz;
614+
merger.GetConstantMem()->calibObjects.fastTransformHelper->Transform(clusters[ihit].sector, clusters[ihit].row, cl.getPad(), cl.getTime(), clx, cly, clz, mTOffset);
615+
if (chkFunction(ihit, cly, clz)) {
616+
xx += clx * clamp;
606617
yy += cly * clamp;
607618
zz += clz * clamp;
608619
clusterState |= clusters[ihit].state;
609620
count += clamp;
610621
}
611-
if (!(ihit + wayDirection >= 0 && ihit + wayDirection < maxN && clusters[ihit].row == clusters[ihit + wayDirection].row && clusters[ihit].sector == clusters[ihit + wayDirection].sector)) {
612-
break;
613-
}
614-
ihit += wayDirection;
615-
}
622+
} while (ihit + wayDirection >= 0 && ihit + wayDirection < maxN && clusters[ihit].row == clusters[ihit + wayDirection].row && clusters[ihit].sector == clusters[ihit + wayDirection].sector);
616623
if (count < 0.1f) {
617624
CADEBUG(printf("\t\tNo matching cluster in double-row, skipping\n"));
625+
xx = tmpX;
618626
return -1;
619627
}
620628
xx /= count;
621629
yy /= count;
622630
zz /= count;
623631
}
632+
if (merger.Param().rec.tpc.rejectIFCLowRadiusCluster) {
633+
const float r2 = xx * xx + yy * yy;
634+
const float rmax2 = CAMath::Square(83.5f + merger.Param().rec.tpc.sysClusErrorMinDist);
635+
if (r2 < rmax2) {
636+
if (rejectChi2) {
637+
MarkClusters(clusters, ihitFirst, ihit, wayDirection, GPUTPCGMMergedTrackHit::flagRejectErr);
638+
}
639+
return -1;
640+
}
641+
}
624642
return 0;
625643
}
626644

0 commit comments

Comments
 (0)