Skip to content

Commit 5d625cc

Browse files
davidrohrdavidrohr
committed
GPU TPC: Improve / Fix looper following over sector boundaries
1 parent 96079f8 commit 5d625cc

File tree

2 files changed

+61
-95
lines changed

2 files changed

+61
-95
lines changed

‎GPU/GPUTracking/Merger/GPUTPCGMTrackParam.cxx‎

Lines changed: 57 additions & 90 deletions
Original file line numberDiff line numberDiff line change
@@ -190,7 +190,7 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int32_
190190
lastSector = cluster.sector;
191191
}
192192
// clang-format off
193-
CADEBUG(printf("\t%21sPropaga Alpha %8.3f , X %8.3f - Y %8.3f, Z %8.3f - QPt %7.2f (%7.2f), SP %5.2f (%5.2f) --- Res %8.3f %8.3f --- Cov sY %8.3f sZ %8.3f sSP %8.3f sPt %8.3f - YPt %8.3f - Err %d", "", prop.GetAlpha(), mX, mP[0], mP[1], mP[4], prop.GetQPt0(), mP[2], prop.GetSinPhi0(), mP[0] - yy, mP[1] - zz, sqrtf(mC[0]), sqrtf(mC[2]), sqrtf(mC[5]), sqrtf(mC[14]), mC[10], err));
193+
CADEBUG(printf("\t%21sPropaga Alpha %8.3f , X %8.3f - Y %8.3f, Z %8.3f - QPt %7.2f (%7.2f), SP %5.2f (%5.2f) --- Res %8.3f %8.3f --- Cov sY %8.3f sZ %8.3f sSP %8.3f sPt %8.3f - YPt %8.3f - Err %d", "", prop.GetAlpha(), mX, mP[0], mP[1], mP[4], prop.GetQPt0(), mP[2], prop.GetSinPhi0(), mP[0] - yy, mP[1] - zz, sqrtf(mC[0]), sqrtf(mC[2]), sqrtf(mC[5]), sqrtf(mC[14]), mC[10], retValProp));
194194
// clang-format on
195195

196196
if (crossCE) {
@@ -324,7 +324,8 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int32_
324324
}
325325
}
326326
if (nWays - iWay <= 2 && !(merger->Param().rec.tpc.disableRefitAttachment & 4)) {
327-
StoreLoopPropagation(merger, lastSector, lastRow, iTrk, lastRow > clusters[0].row, prop.GetAlpha());
327+
StoreLoopPropagation(merger, lastSector, lastRow, iTrk, lastRow > clusters[(iWay & 1) ? (maxN - 1) : 0].row, prop.GetAlpha());
328+
CADEBUG(printf("STORING %d lastRow %d row %d out %d\n", iTrk, (int)lastRow, (int)clusters[(iWay & 1) ? (maxN - 1) : 0].row, lastRow > clusters[(iWay & 1) ? (maxN - 1) : 0].row));
328329
}
329330
if (((nWays - iWay) & 1) && (iWay != nWays - 1) && !track.CCE() && !track.Looper()) {
330331
ShiftZ(clusters, merger, maxN);
@@ -587,16 +588,9 @@ GPUd() bool GPUTPCGMTrackParam::AttachClustersPropagate(const GPUTPCGMMerger* GP
587588
return dodEdx;
588589
}
589590

590-
GPUd() bool GPUTPCGMTrackParam::FollowCircleChk(float lrFactor, float toY, float toX, bool up, bool right)
591-
{
592-
return CAMath::Abs(mX * lrFactor - toY) > 1.f && // transport further in Y
593-
CAMath::Abs(mP[2]) < 0.7f && // rotate back
594-
(up ? (-mP[0] * lrFactor > toX || (right ^ (mP[2] > 0))) : (-mP[0] * lrFactor < toX || (right ^ (mP[2] < 0)))); // don't overshoot in X
595-
}
596-
597591
GPUdii() void GPUTPCGMTrackParam::StoreOuter(gputpcgmmergertypes::GPUTPCOuterParam* outerParam, float alpha)
598592
{
599-
CADEBUG(printf("\t%21sStorO Alpha %8.3f , X %8.3f - Y %8.3f, Z %8.3f - QPt %7.2f (%7.2f), SP %5.2f (%5.2f) --- Cov sY %8.3f sZ %8.3f sSP %8.3f sPt %8.3f\n", "", prop.GetAlpha(), mX, mP[0], mP[1], mP[4], prop.GetQPt0(), mP[2], prop.GetSinPhi0(), sqrtf(mC[0]), sqrtf(mC[2]), sqrtf(mC[5]), sqrtf(mC[14])));
593+
CADEBUG(printf("\t%21sStorO Alpha %8.3f , X %8.3f - Y %8.3f, Z %8.3f - QPt %7.2f, SP %5.2f --- Cov sY %8.3f sZ %8.3f sSP %8.3f sPt %8.3f\n", "", alpha, mX, mP[0], mP[1], mP[4], mP[2], sqrtf(mC[0]), sqrtf(mC[2]), sqrtf(mC[5]), sqrtf(mC[14])));
600594
for (int32_t i = 0; i < 5; i++) {
601595
outerParam->P[i] = mP[i];
602596
}
@@ -607,18 +601,18 @@ GPUdii() void GPUTPCGMTrackParam::StoreOuter(gputpcgmmergertypes::GPUTPCOuterPar
607601
outerParam->alpha = alpha;
608602
}
609603

610-
GPUdic(0, 1) void GPUTPCGMTrackParam::StoreLoopPropagation(const GPUTPCGMMerger* GPUrestrict() Merger, int32_t sector, int32_t iRow, int32_t iTrack, bool outerParam, float alpha)
604+
GPUdic(0, 1) void GPUTPCGMTrackParam::StoreLoopPropagation(const GPUTPCGMMerger* GPUrestrict() Merger, int32_t sector, int32_t iRow, int32_t iTrack, bool outwards, float alpha)
611605
{
612606
if (iRow == 0 || iRow == GPUCA_ROW_COUNT - 1) {
613607
return;
614608
}
615-
if (CAMath::Abs(mP[2]) >= GPUCA_MAX_SIN_PHI_LOW) {
609+
if (CAMath::Abs(mP[2]) >= GPUCA_MAX_SIN_PHI) { // TODO: How can we avoid this?
616610
return;
617611
}
618612
if (CAMath::Abs(mP[2]) < 0.75) {
619613
return;
620614
}
621-
if ((mP[2] * mP[4] < 0) ^ outerParam) {
615+
if ((mP[2] * mP[4] < 0) ^ outwards) {
622616
return;
623617
}
624618

@@ -634,7 +628,7 @@ GPUdic(0, 1) void GPUTPCGMTrackParam::StoreLoopPropagation(const GPUTPCGMMerger*
634628
data.alpha = alpha;
635629
data.sector = sector;
636630
data.row = iRow;
637-
data.outerParam = outerParam;
631+
data.outwards = outwards;
638632
Merger->LoopData()[nLoopData] = data;
639633
}
640634

@@ -651,107 +645,80 @@ GPUdii() void GPUTPCGMTrackParam::PropagateLooper(const GPUTPCGMMerger* GPUrestr
651645

652646
GPUTPCGMLoopData& data = Merger->LoopData()[loopIdx];
653647
prop.SetTrack(&data.param, data.alpha);
654-
data.param.AttachClustersLooper(Merger, data.sector, data.row, data.track, data.outerParam, prop);
655-
// data.param.FollowCircle(Merger, prop, data.sector, data.row, data.track, data.toAlpha, data.toX, data.toY, data.toSector, data.toRow, data.inFlyDirection);
648+
if (false) {
649+
data.param.AttachClustersLooper(Merger, data.sector, data.row, data.track, data.outwards, prop);
650+
} else {
651+
data.param.AttachClustersLooperFollow(Merger, prop, data.sector, data.row, data.track, data.outwards);
652+
}
656653
}
657654

658-
GPUdi() int32_t GPUTPCGMTrackParam::FollowCircle(const GPUTPCGMMerger* GPUrestrict() Merger, GPUTPCGMPropagator& GPUrestrict() prop, int32_t sector, int32_t iRow, int32_t iTrack, float toAlpha, float toX, float toY, int32_t toSector, int32_t toRow, bool inFlyDirection)
655+
GPUdi() void GPUTPCGMTrackParam::AttachClustersLooperFollow(const GPUTPCGMMerger* GPUrestrict() Merger, GPUTPCGMPropagator& GPUrestrict() prop, int32_t sector, int32_t iRow, int32_t iTrack, bool up)
659656
{
657+
float toX = mX;
658+
bool inFlyDirection = (Merger->MergedTracks()[iTrack].Leg() & 1) ^ up;
659+
660660
static constexpr float kSectAngle = 2 * M_PI / 18.f;
661661
const GPUParam& GPUrestrict() param = Merger->Param();
662-
bool right;
663-
float dAlpha = toAlpha - prop.GetAlpha();
664-
int32_t sectorSide = sector >= (GPUCA_NSECTORS / 2) ? (GPUCA_NSECTORS / 2) : 0;
665-
if (CAMath::Abs(dAlpha) > 0.001f) {
666-
right = CAMath::Abs(dAlpha) < CAMath::Pi() ? (dAlpha > 0) : (dAlpha < 0);
667-
} else {
668-
right = toY > mP[0];
669-
}
670-
bool up = (mP[2] < 0) ^ right;
671-
int32_t targetRow = up ? (GPUCA_ROW_COUNT - 1) : 0;
672-
float lrFactor = mP[2] < 0 ? -1.f : 1.f; // !(right ^ down) // TODO: shouldn't it be "right ? 1.f : -1.f", but that gives worse results...
662+
bool right = (mP[2] < 0) ^ up;
663+
const int32_t sectorSide = sector >= (GPUCA_NSECTORS / 2) ? (GPUCA_NSECTORS / 2) : 0;
664+
float lrFactor = right ^ !up ? 1.f : -1.f;
673665
// clang-format off
674-
CADEBUG(printf("CIRCLE Track %d: Sector %d Alpha %f X %f Y %f Z %f SinPhi %f DzDs %f - Next hit: Sector %d Alpha %f X %f Y %f - Right %d Up %d dAlpha %f lrFactor %f\n", iTrack, sector, prop.GetAlpha(), mX, mP[0], mP[1], mP[2], mP[3], toSector, toAlpha, toX, toY, (int32_t)right, (int32_t)up, dAlpha, lrFactor));
666+
CADEBUG(printf("\nCIRCLE Track %d: Sector %d Alpha %f X %f Y %f Z %f SinPhi %f DzDs %f QPt %f - Right %d Up %d lrFactor %f\n", iTrack, sector, prop.GetAlpha(), mX, mP[0], mP[1], mP[2], mP[3], mP[4], (int32_t)right, (int32_t)up, lrFactor));
675667
// clang-format on
676668

677-
AttachClustersPropagate(Merger, sector, iRow, targetRow, iTrack, false, prop, inFlyDirection, 0.7f);
678669
if (prop.RotateToAlpha(prop.GetAlpha() + (CAMath::Pi() / 2.f) * lrFactor)) {
679-
return 1;
670+
return;
680671
}
681672
CADEBUG(printf("\tRotated: X %f Y %f Z %f SinPhi %f (Alpha %f / %f)\n", mP[0], mX, mP[1], mP[2], prop.GetAlpha(), prop.GetAlpha() + CAMath::Pi() / 2.f));
682-
while (sector != toSector || FollowCircleChk(lrFactor, toY, toX, up, right)) {
683-
while ((sector != toSector) ? (CAMath::Abs(mX) <= CAMath::Abs(mP[0]) * CAMath::Tan(kSectAngle / 2.f)) : FollowCircleChk(lrFactor, toY, toX, up, right)) {
684-
int32_t err = prop.PropagateToXAlpha(mX + 1.f, prop.GetAlpha(), inFlyDirection);
673+
uint32_t maxTries = 100;
674+
while (true) {
675+
while (CAMath::Abs(mX) <= CAMath::Abs(mP[0]) * CAMath::Tan(kSectAngle / 2.f)) {
676+
if (maxTries-- == 0) {
677+
return;
678+
}
679+
if (CAMath::Abs(mP[2]) > 0.7f) {
680+
return;
681+
}
682+
if (up ? (-mP[0] * lrFactor > GPUTPCGeometry::Row2X(GPUCA_ROW_COUNT - 1)) : (-mP[0] * lrFactor < GPUTPCGeometry::Row2X(0))) {
683+
return;
684+
}
685+
if (!((up ? (-mP[0] * lrFactor >= toX) : (-mP[0] * lrFactor <= toX)) || (right ^ (mP[2] > 0)))) {
686+
return;
687+
}
688+
int32_t err = prop.PropagateToXAlpha(mX + (up ? 1.f : -1.f), prop.GetAlpha(), inFlyDirection);
685689
if (err) {
686690
CADEBUG(printf("\t\tpropagation error (%d)\n", err));
687-
prop.RotateToAlpha(prop.GetAlpha() - (CAMath::Pi() / 2.f) * lrFactor);
688-
return 1;
691+
return;
689692
}
690693
CADEBUG(printf("\tPropagated to y = %f: X %f Z %f SinPhi %f\n", mX, mP[0], mP[1], mP[2]));
691-
for (int32_t j = 0; j < GPUCA_ROW_COUNT; j++) {
694+
for (int32_t j = 0; j < GPUCA_ROW_COUNT; j++) { // TODO: Avoid iterating over all rows
692695
float rowX = GPUTPCGeometry::Row2X(j);
693696
if (CAMath::Abs(rowX - (-mP[0] * lrFactor)) < 1.5f) {
694-
CADEBUG(printf("\t\tAttempt row %d (Y %f Z %f)\n", j, mX * lrFactor, mP[1]));
697+
CADEBUG(printf("\t\tAttempt row %d (X %f Y %f Z %f)\n", j, rowX, mX * lrFactor, mP[1]));
695698
AttachClusters(Merger, sector, j, iTrack, false, mX * lrFactor, mP[1]);
696699
}
697700
}
698701
}
699-
if (sector != toSector) {
700-
if (right) {
701-
if (++sector >= sectorSide + 18) {
702-
sector -= 18;
703-
}
704-
} else {
705-
if (--sector < sectorSide) {
706-
sector += 18;
707-
}
708-
}
709-
CADEBUG(printf("\tRotating to sector %d\n", sector));
710-
if (prop.RotateToAlpha(param.Alpha(sector) + (CAMath::Pi() / 2.f) * lrFactor)) {
711-
CADEBUG(printf("\t\trotation error\n"));
712-
prop.RotateToAlpha(prop.GetAlpha() - (CAMath::Pi() / 2.f) * lrFactor);
713-
return 1;
714-
}
715-
CADEBUG(printf("\tAfter Rotatin Alpha %f Position X %f Y %f Z %f SinPhi %f\n", prop.GetAlpha(), mP[0], mX, mP[1], mP[2]));
716-
}
717-
}
718-
CADEBUG(printf("\tRotating back\n"));
719-
for (int32_t i = 0; i < 2; i++) {
720-
if (prop.RotateToAlpha(prop.GetAlpha() + (CAMath::Pi() / 2.f) * lrFactor) == 0) {
721-
break;
722-
}
723-
if (i) {
724-
CADEBUG(printf("Final rotation failed\n"));
725-
return 1;
726-
}
727-
CADEBUG(printf("\tresetting physical model\n"));
728-
prop.SetTrack(this, prop.GetAlpha());
729-
}
730-
prop.Rotate180();
731-
CADEBUG(printf("\tMirrored position: Alpha %f X %f Y %f Z %f SinPhi %f DzDs %f\n", prop.GetAlpha(), mX, mP[0], mP[1], mP[2], mP[3]));
732-
iRow = toRow;
733-
float dx = toX - GPUTPCGeometry::Row2X(toRow);
734-
if (up ^ (toX > mX)) {
735-
if (up) {
736-
while (iRow < GPUCA_ROW_COUNT - 2 && GPUTPCGeometry::Row2X(iRow + 1) + dx <= mX) {
737-
iRow++;
702+
if (right) {
703+
if (++sector >= sectorSide + 18) {
704+
sector -= 18;
738705
}
739706
} else {
740-
while (iRow > 1 && GPUTPCGeometry::Row2X(iRow - 1) + dx >= mX) {
741-
iRow--;
707+
if (--sector < sectorSide) {
708+
sector += 18;
742709
}
743710
}
744-
prop.PropagateToXAlpha(GPUTPCGeometry::Row2X(iRow) + dx, prop.GetAlpha(), inFlyDirection);
745-
AttachClustersPropagate(Merger, sector, iRow, toRow, iTrack, false, prop, inFlyDirection);
746-
}
747-
if (prop.PropagateToXAlpha(toX, prop.GetAlpha(), inFlyDirection)) {
748-
mX = toX;
711+
CADEBUG(printf("\tRotating to sector %d: %f --> %f\n", sector, prop.GetAlpha(), param.Alpha(sector) + (CAMath::Pi() / 2.f) * lrFactor));
712+
int32_t err = prop.RotateToAlpha(param.Alpha(sector) + (CAMath::Pi() / 2.f) * lrFactor);
713+
if (err) {
714+
CADEBUG(printf("Rotation Error %d\n", err));
715+
return;
716+
}
717+
CADEBUG(printf("\tAfter Rotating Alpha %f Position X %f Y %f Z %f SinPhi %f\n", prop.GetAlpha(), mP[0], mX, mP[1], mP[2]));
749718
}
750-
CADEBUG(printf("Final position: Alpha %f X %f Y %f Z %f SinPhi %f DzDs %f\n", prop.GetAlpha(), mX, mP[0], mP[1], mP[2], mP[3]));
751-
return (0);
752719
}
753720

754-
GPUdi() void GPUTPCGMTrackParam::AttachClustersLooper(const GPUTPCGMMerger* GPUrestrict() Merger, int32_t sector, int32_t iRow, int32_t iTrack, bool outer, GPUTPCGMPropagator& GPUrestrict() prop)
721+
GPUdi() void GPUTPCGMTrackParam::AttachClustersLooper(const GPUTPCGMMerger* GPUrestrict() Merger, int32_t sector, int32_t iRow, int32_t iTrack, bool outwards, GPUTPCGMPropagator& GPUrestrict() prop)
755722
{
756723
static constexpr float kSectAngle = 2 * M_PI / 18.f;
757724
// Note that the coordinate system is rotated by 90 degree swapping X and Y!
@@ -761,9 +728,9 @@ GPUdi() void GPUTPCGMTrackParam::AttachClustersLooper(const GPUTPCGMMerger* GPUr
761728
float SinPhi = CAMath::Sqrt(1 - mP[2] * mP[2]) * (mP[2] > 0 ? -1 : 1);
762729
float b = prop.GetBz(prop.GetAlpha(), mX, mP[0], mP[1]);
763730

764-
float dx = outer ? 1.f : -1.f;
731+
float dx = outwards ? 1.f : -1.f;
765732
const float myRowX = GPUTPCGeometry::Row2X(iRow);
766-
// printf("\nAttachMirror sector %d row %d outer %d\n", (int)sector, (int)iRow, (int)outer);
733+
// printf("\nAttachMirror sector %d row %d outwards %d\n", (int)sector, (int)iRow, (int)outwards);
767734
// printf("X %f Y %f Z %f SinPhi %f -->\n", mX, mP[0], mP[1], mP[2]);
768735
// printf("X %f Y %f Z %f SinPhi %f, dx %f\n", X, Y, Z, SinPhi, dx);
769736
uint32_t maxTries = 100;
@@ -795,7 +762,7 @@ GPUdi() void GPUTPCGMTrackParam::AttachClustersLooper(const GPUTPCGMMerger* GPUr
795762

796763
// printf("count %d: At X %f Y %f Z %f SinPhi %f\n", maxTries, mP[2] > 0 ? -Y : Y, mP[2] > 0 ? X : -X, Z, SinPhi);
797764
float paramX = mP[2] > 0 ? -Y : Y;
798-
int32_t step = outer ? 1 : -1;
765+
int32_t step = outwards ? 1 : -1;
799766
int32_t found = 0;
800767
for (int32_t j = iRow; j >= 0 && j < GPUCA_ROW_COUNT && found < 3; j += step) {
801768
float rowX = mX + GPUTPCGeometry::Row2X(j) - myRowX;

‎GPU/GPUTracking/Merger/GPUTPCGMTrackParam.h‎

Lines changed: 4 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -150,9 +150,9 @@ class GPUTPCGMTrackParam
150150
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
151151
GPUd() float AttachClusters(const GPUTPCGMMerger* GPUrestrict() Merger, int32_t sector, int32_t iRow, int32_t iTrack, bool goodLeg, float Y, float Z);
152152
// We force to compile these twice, for PropagateLooper and for Fit, for better optimization
153-
GPUd() void AttachClustersLooper(const GPUTPCGMMerger* GPUrestrict() Merger, int32_t sector, int32_t iRow, int32_t iTrack, bool outer, GPUTPCGMPropagator& prop);
154-
GPUd() int32_t FollowCircle(const GPUTPCGMMerger* GPUrestrict() Merger, GPUTPCGMPropagator& prop, int32_t sector, int32_t iRow, int32_t iTrack, float toAlpha, float toX, float toY, int32_t toSector, int32_t toRow, bool inFlyDirection);
155-
GPUd() void StoreLoopPropagation(const GPUTPCGMMerger* GPUrestrict() Merger, int32_t sector, int32_t iRow, int32_t iTrack, bool outerParam, float alpha);
153+
GPUd() void AttachClustersLooper(const GPUTPCGMMerger* GPUrestrict() Merger, int32_t sector, int32_t iRow, int32_t iTrack, bool outwards, GPUTPCGMPropagator& prop);
154+
GPUd() void AttachClustersLooperFollow(const GPUTPCGMMerger* GPUrestrict() Merger, GPUTPCGMPropagator& prop, int32_t sector, int32_t iRow, int32_t iTrack, bool outwards);
155+
GPUd() void StoreLoopPropagation(const GPUTPCGMMerger* GPUrestrict() Merger, int32_t sector, int32_t iRow, int32_t iTrack, bool outwards, float alpha);
156156
GPUd() void StoreOuter(gputpcgmmergertypes::GPUTPCOuterParam* outerParam, float alpha);
157157
GPUd() static void PropagateLooper(const GPUTPCGMMerger* GPUrestrict() Merger, int32_t loopIdx);
158158

@@ -213,7 +213,6 @@ class GPUTPCGMTrackParam
213213
}
214214

215215
private:
216-
GPUd() bool FollowCircleChk(float lrFactor, float toY, float toX, bool up, bool right);
217216
GPUd() int32_t initResetT0();
218217

219218
float mX; // x position
@@ -230,7 +229,7 @@ struct GPUTPCGMLoopData {
230229
float alpha;
231230
uint8_t sector;
232231
uint8_t row;
233-
uint8_t outerParam;
232+
uint8_t outwards;
234233
};
235234

236235
GPUdi() int32_t GPUTPCGMTrackParam::initResetT0()

0 commit comments

Comments
 (0)