Skip to content

Commit cfac307

Browse files
davidrohrdavidrohr
committed
GPU TPC: Refactor cluster merging
1 parent e459444 commit cfac307

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
@@ -122,26 +122,13 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger& GPUrestrict() merger, int32_
122122
uint8_t clusterState = clusters[ihit].state;
123123
const float clAlpha = param.Alpha(clusters[ihit].sector);
124124
float xx, yy, zz;
125-
{
126-
const ClusterNative& GPUrestrict() cl = merger.GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num];
127-
merger.GetConstantMem()->calibObjects.fastTransformHelper->Transform(clusters[ihit].sector, clusters[ihit].row, cl.getPad(), cl.getTime(), xx, yy, zz, mTOffset);
128-
}
129125
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));
130-
if (MergeDoubleRowClusters(ihit, wayDirection, clusters, merger, prop, xx, yy, zz, maxN, clAlpha, clusterState, !param.rec.tpc.rebuildTrackInFit && allowChangeClusters) == -1) {
126+
if (MergeDoubleRowClusters(ihit, wayDirection, clusters, merger, prop, xx, yy, zz, maxN, clAlpha, clusterState, param.rec.tpc.rebuildTrackInFit ? rebuilt : allowChangeClusters) == -1) {
131127
nMissed++;
132128
nMissed2++;
133129
continue;
134130
}
135131

136-
if (param.rec.tpc.rejectIFCLowRadiusCluster) {
137-
const float r2 = xx * xx + yy * yy;
138-
const float rmax = (83.5f + param.rec.tpc.sysClusErrorMinDist);
139-
if (r2 < rmax * rmax) {
140-
MarkClusters(clusters, ihitMergeFirst, ihit, wayDirection, GPUTPCGMMergedTrackHit::flagRejectErr);
141-
continue;
142-
}
143-
}
144-
145132
const auto& cluster = clusters[ihit];
146133
if (lastRow != 255 && CAMath::Abs(cluster.row - lastRow) > 1) {
147134
bool dodEdx = param.dodEdxEnabled && param.rec.tpc.adddEdxSubThresholdClusters && finalFit && CAMath::Abs(cluster.row - lastRow) == 2;
@@ -570,55 +557,86 @@ GPUd() void GPUTPCGMTrackParam::MirrorTo(GPUTPCGMPropagator& GPUrestrict() prop,
570557

571558
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)
572559
{
560+
const int32_t ihitFirst = ihit;
561+
{
562+
const ClusterNative& GPUrestrict() cl = merger.GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num];
563+
merger.GetConstantMem()->calibObjects.fastTransformHelper->Transform(clusters[ihit].sector, clusters[ihit].row, cl.getPad(), cl.getTime(), xx, yy, zz, mTOffset);
564+
}
573565
if (ihit + wayDirection >= 0 && ihit + wayDirection < maxN && clusters[ihit].row == clusters[ihit + wayDirection].row && clusters[ihit].sector == clusters[ihit + wayDirection].sector) {
574-
float maxDistY, maxDistZ;
575-
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
576-
maxDistY = (maxDistY + mC[0]) * 20.f;
577-
maxDistZ = (maxDistZ + mC[2]) * 20.f;
578-
int32_t noReject = 0; // Cannot reject if simple estimation of y/z fails (extremely unlike case)
566+
float maxDistY2, maxDistZ2;
567+
bool noReject = false; // Cannot reject if simple estimation of y/z fails (extremely unlike case)
579568
if (CAMath::Abs(clAlpha - prop.GetAlpha()) > 1.e-4f) {
580569
noReject = prop.RotateToAlpha(clAlpha);
581570
}
582571
float projY = 0, projZ = 0;
583-
if (noReject == 0) {
572+
if (!noReject) {
584573
noReject |= prop.GetPropagatedYZ(xx, projY, projZ);
585574
}
586-
float count = 0.f;
587-
xx = yy = zz = 0.f;
588-
clusterState = 0;
589-
while (true) {
590-
const ClusterNative& GPUrestrict() cl = merger.GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num];
591-
float clamp = cl.qTot;
592-
float clx, cly, clz;
593-
merger.GetConstantMem()->calibObjects.fastTransformHelper->Transform(clusters[ihit].sector, clusters[ihit].row, cl.getPad(), cl.getTime(), clx, cly, clz, mTOffset);
594-
float dy = cly - projY;
595-
float dz = clz - projZ;
596-
if (noReject == 0 && (dy * dy > maxDistY || dz * dz > maxDistZ)) {
597-
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));
575+
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
576+
const float kFactor = merger.GetConstantMem()->tpcTrackers[0].GetChiSeedFactor() * 4.f;
577+
maxDistY2 = (maxDistY2 + mC[0]) * kFactor;
578+
maxDistZ2 = (maxDistZ2 + mC[2]) * kFactor;
579+
auto chkFunction = [clusters, rejectChi2, maxDistY2, maxDistZ2, projY, projZ, noReject](int32_t ih, float y, float z) {
580+
float dy = y - projY;
581+
float dz = z - projZ;
582+
if (!noReject && (dy * dy > maxDistY2 || dz * dz > maxDistZ2)) {
583+
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));
598584
if (rejectChi2) {
599-
clusters[ihit].state |= GPUTPCGMMergedTrackHit::flagRejectDistance;
585+
clusters[ih].state |= GPUTPCGMMergedTrackHit::flagRejectDistance;
600586
}
587+
return false;
601588
} else {
602-
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)));
603-
xx += clx * clamp; // TODO: Weight in pad/time instead of XYZ
589+
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)));
590+
return true;
591+
}
592+
};
593+
const float tmpX = xx;
594+
float count;
595+
if (chkFunction(ihit, yy, zz)) {
596+
const ClusterNative& GPUrestrict() cl = merger.GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num];
597+
const float clamp = cl.qTot;
598+
xx *= clamp;
599+
yy *= clamp;
600+
zz *= clamp;
601+
clusterState = clusters[ihit].state;
602+
count = clamp;
603+
} else {
604+
xx = yy = zz = count = 0.f;
605+
clusterState = 0;
606+
}
607+
do {
608+
ihit += wayDirection;
609+
const ClusterNative& GPUrestrict() cl = merger.GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num];
610+
const float clamp = cl.qTot;
611+
float clx, cly, clz;
612+
merger.GetConstantMem()->calibObjects.fastTransformHelper->Transform(clusters[ihit].sector, clusters[ihit].row, cl.getPad(), cl.getTime(), clx, cly, clz, mTOffset);
613+
if (chkFunction(ihit, cly, clz)) {
614+
xx += clx * clamp;
604615
yy += cly * clamp;
605616
zz += clz * clamp;
606617
clusterState |= clusters[ihit].state;
607618
count += clamp;
608619
}
609-
if (!(ihit + wayDirection >= 0 && ihit + wayDirection < maxN && clusters[ihit].row == clusters[ihit + wayDirection].row && clusters[ihit].sector == clusters[ihit + wayDirection].sector)) {
610-
break;
611-
}
612-
ihit += wayDirection;
613-
}
620+
} while (ihit + wayDirection >= 0 && ihit + wayDirection < maxN && clusters[ihit].row == clusters[ihit + wayDirection].row && clusters[ihit].sector == clusters[ihit + wayDirection].sector);
614621
if (count < 0.1f) {
615622
CADEBUG(printf("\t\tNo matching cluster in double-row, skipping\n"));
623+
xx = tmpX;
616624
return -1;
617625
}
618626
xx /= count;
619627
yy /= count;
620628
zz /= count;
621629
}
630+
if (merger.Param().rec.tpc.rejectIFCLowRadiusCluster) {
631+
const float r2 = xx * xx + yy * yy;
632+
const float rmax2 = CAMath::Square(83.5f + merger.Param().rec.tpc.sysClusErrorMinDist);
633+
if (r2 < rmax2) {
634+
if (rejectChi2) {
635+
MarkClusters(clusters, ihitFirst, ihit, wayDirection, GPUTPCGMMergedTrackHit::flagRejectErr);
636+
}
637+
return -1;
638+
}
639+
}
622640
return 0;
623641
}
624642

0 commit comments

Comments
 (0)