@@ -1860,13 +1860,12 @@ GPUd() void GPUTPCGMMerger::Finalize2(int32_t nBlocks, int32_t nThreads, int32_t
18601860
18611861GPUd () void GPUTPCGMMerger::MergeLoopersInit(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread)
18621862{
1863- return ; // FIXME: !!!!
1864- const float lowPtThresh = Param ().rec .tpc .rejectQPtB5 * 1 .1f ; // Might need to merge tracks above the threshold with parts below the threshold
1863+ const float lowPtThresh = Param ().rec .tpc .rejectQPtB5 * 1 .1f ; // Might need to merge tracks above the threshold with parts below the rejection threshold
18651864 for (uint32_t i = get_global_id (0 ); i < mMemory ->nMergedTracks ; i += get_global_size (0 )) {
18661865 const auto & trk = mMergedTracks [i];
18671866 const auto & p = trk.GetParam ();
18681867 const float qptabs = CAMath::Abs (p.GetQPt ());
1869- if (trk.NClusters () && qptabs * Param ().qptB5Scaler > 5 .f && qptabs * Param ().qptB5Scaler <= lowPtThresh) {
1868+ if (trk.OK () && trk. NClusters () && trk. Leg () == 0 && qptabs * Param ().qptB5Scaler > 5 .f && qptabs * Param ().qptB5Scaler <= lowPtThresh) {
18701869 const int32_t sector = mClusters [trk.FirstClusterRef () + trk.NClusters () - 1 ].sector ;
18711870 const float refz = p.GetZ () + GetConstantMem ()->calibObjects .fastTransformHelper ->getCorrMap ()->convVertexTimeToZOffset (sector, p.GetTOffset (), Param ().continuousMaxTimeBin ) + (trk.CSide () ? -100 : 100 );
18721871 float sinA, cosA;
@@ -1915,12 +1914,12 @@ GPUd() void GPUTPCGMMerger::MergeLoopersSort(int32_t nBlocks, int32_t nThreads,
19151914
19161915GPUd () void GPUTPCGMMerger::MergeLoopersMain(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread)
19171916{
1918- const MergeLooperParam* params = mLooperCandidates ;
1917+ const MergeLooperParam* candidates = mLooperCandidates ;
19191918
19201919#if GPUCA_MERGE_LOOPER_MC && !defined(GPUCA_GPUCODE)
19211920 std::vector<int64_t > paramLabels (mMemory ->nLooperMatchCandidates );
19221921 for (uint32_t i = 0 ; i < mMemory ->nLooperMatchCandidates ; i++) {
1923- paramLabels[i] = GetTrackLabel (mMergedTracks [params [i].id ]);
1922+ paramLabels[i] = GetTrackLabel (mMergedTracks [candidates [i].id ]);
19241923 }
19251924 /* std::vector<bool> dropped(mMemory->nLooperMatchCandidates);
19261925 std::vector<bool> droppedMC(mMemory->nLooperMatchCandidates);
@@ -1934,26 +1933,47 @@ GPUd() void GPUTPCGMMerger::MergeLoopersMain(int32_t nBlocks, int32_t nThreads,
19341933 for (uint32_t i = get_global_id (0 ); i < mMemory ->nLooperMatchCandidates ; i += get_global_size (0 )) {
19351934 for (uint32_t j = i + 1 ; j < mMemory ->nLooperMatchCandidates ; j++) {
19361935 // int32_t bs = 0;
1937- if (CAMath::Abs (params[j].refz ) > CAMath::Abs (params[i].refz ) + 100 .f ) {
1936+ assert (CAMath::Abs (candidates[i].refz ) <= CAMath::Abs (candidates[j].refz ));
1937+ if (CAMath::Abs (candidates[j].refz ) > CAMath::Abs (candidates[i].refz ) + 100 .f ) {
19381938 break ;
19391939 }
1940- const float d2xy = CAMath::Sum2 (params [i].x - params [j].x , params [i].y - params [j].y );
1940+ const float d2xy = CAMath::Sum2 (candidates [i].x - candidates [j].x , candidates [i].y - candidates [j].y );
19411941 if (d2xy > 15 .f ) {
19421942 // bs |= 1;
19431943 continue ;
19441944 }
1945- const auto & trk1 = mMergedTracks [params[i].id ];
1946- const auto & trk2 = mMergedTracks [params[j].id ];
1945+
1946+ const GPUTPCGMMergedTrack* trkI = &mMergedTracks [candidates[i].id ];
1947+ float refZI = candidates[i].refz ;
1948+ {
1949+ const auto * tmp = trkI;
1950+ while (tmp->PrevSegment () >= 0 ) {
1951+ const auto * next = &mMergedTracks [tmp->PrevSegment ()];
1952+ if (next == trkI) {
1953+ break ;
1954+ }
1955+ tmp = next;
1956+ }
1957+ if (tmp != trkI && tmp->CSide () == trkI->CSide () && CAMath::Abs (tmp->GetParam ().GetZ ()) > CAMath::Abs (trkI->GetParam ().GetZ ())) {
1958+ float tmpRefZ = refZI + tmp->GetParam ().GetZ () - trkI->GetParam ().GetZ ();
1959+ if (CAMath::Abs (tmpRefZ) < CAMath::Abs (candidates[j].refz ) && CAMath::Abs (tmpRefZ) > CAMath::Abs (refZI)) {
1960+ trkI = tmp;
1961+ refZI = tmpRefZ;
1962+ }
1963+ }
1964+ };
1965+ const auto & trk1 = *trkI;
1966+ const auto & trk2 = mMergedTracks [candidates[j].id ];
19471967 const auto & param1 = trk1.GetParam ();
19481968 const auto & param2 = trk2.GetParam ();
19491969 if (CAMath::Abs (param1.GetDzDs ()) > 0 .03f && CAMath::Abs (param2.GetDzDs ()) > 0 .03f && param1.GetDzDs () * param2.GetDzDs () * param1.GetQPt () * param2.GetQPt () < 0 ) {
19501970 // bs |= 2;
19511971 continue ;
19521972 }
19531973
1954- const float dznormalized = (CAMath::Abs (params [j].refz ) - CAMath::Abs (params[i]. refz )) / (CAMath::TwoPi () * 0 .5f * (CAMath::Abs (param1.GetDzDs ()) + CAMath::Abs (param2.GetDzDs ())) * 1 .f / (0 .5f * (CAMath::Abs (param1.GetQPt ()) + CAMath::Abs (param2.GetQPt ())) * CAMath::Abs (Param ().polynomialField .GetNominalBz ())));
1974+ const float dznormalized = (CAMath::Abs (candidates [j].refz ) - CAMath::Abs (refZI )) / (CAMath::TwoPi () * 0 .5f * (CAMath::Abs (param1.GetDzDs ()) + CAMath::Abs (param2.GetDzDs ())) * 1 .f / (0 .5f * (CAMath::Abs (param1.GetQPt ()) + CAMath::Abs (param2.GetQPt ())) * CAMath::Abs (Param ().polynomialField .GetNominalBz ())));
19551975 const float phasecorr = CAMath::Modf ((CAMath::ASin (param1.GetSinPhi ()) + trk1.GetAlpha () - CAMath::ASin (param2.GetSinPhi ()) - trk2.GetAlpha ()) / CAMath::TwoPi () + 5 .5f , 1 .f ) - 0 .5f ;
1956- const float phasecorrdirection = (params [j].refz * param1.GetQPt () * param1.GetDzDs ()) > 0 ? 1 : -1 ;
1976+ const float phasecorrdirection = (candidates [j].refz * param1.GetQPt () * param1.GetDzDs ()) > 0 ? 1 : -1 ;
19571977 const float dzcorr = dznormalized + phasecorr * phasecorrdirection;
19581978 const bool sameside = !(trk1.CSide () ^ trk2.CSide ());
19591979 const float dzcorrlimit[4 ] = {sameside ? 0 .018f : 0 .012f , sameside ? 0 .12f : 0 .025f , 0 .14f , 0 .15f };
@@ -1982,11 +2002,11 @@ GPUd() void GPUTPCGMMerger::MergeLoopersMain(int32_t nBlocks, int32_t nThreads,
19822002 const int64_t label2 = paramLabels[j];
19832003 bool labelEQ = label1 != -1 && label1 == label2;
19842004 if (1 || EQ || labelEQ) {
1985- // printf("Matching track %d/%d %u-%u (%ld/%ld): dist %f side %d %d, tgl %f %f, qpt %f %f, x %f %f, y %f %f\n", (int32_t)EQ, (int32_t)labelEQ, i, j, label1, label2, d, (int32_t)mMergedTracks[params [i].id].CSide(), (int32_t)mMergedTracks[params [j].id].CSide(), params [i].tgl, params [j].tgl, params [i].qpt, params [j].qpt, params [i].x, params [j].x, params [i].y, params [j].y);
2005+ // printf("Matching track %d/%d %u-%u (%ld/%ld): dist %f side %d %d, tgl %f %f, qpt %f %f, x %f %f, y %f %f\n", (int32_t)EQ, (int32_t)labelEQ, i, j, label1, label2, d, (int32_t)mMergedTracks[candidates [i].id].CSide(), (int32_t)mMergedTracks[candidates [j].id].CSide(), candidates [i].tgl, candidates [j].tgl, candidates [i].qpt, candidates [j].qpt, candidates [i].x, candidates [j].x, candidates [i].y, candidates [j].y);
19862006 static auto & tup = GPUROOTDump<TNtuple>::get (" mergeloopers" , " labeleq:sides:d2xy:tgl1:tgl2:qpt1:qpt2:dz:dzcorr:dtgl:dqpt:dznorm:bs" );
1987- tup.Fill ((float )labelEQ, (trk1.CSide () ? 1 : 0 ) | (trk2.CSide () ? 2 : 0 ), d2xy, param1.GetDzDs (), param2.GetDzDs (), param1.GetQPt (), param2.GetQPt (), CAMath::Abs (params [j].refz ) - CAMath::Abs (params[i]. refz ), dzcorr, dtgl, dqpt, dznorm, bs);
2007+ tup.Fill ((float )labelEQ, (trk1.CSide () ? 1 : 0 ) | (trk2.CSide () ? 2 : 0 ), d2xy, param1.GetDzDs (), param2.GetDzDs (), param1.GetQPt (), param2.GetQPt (), CAMath::Abs (candidates [j].refz ) - CAMath::Abs (refZI ), dzcorr, dtgl, dqpt, dznorm, bs);
19882008 static auto tup2 = GPUROOTDump<TNtuple>::getNew (" mergeloopers2" , " labeleq:refz1:refz2:tgl1:tgl2:qpt1:qpt2:snp1:snp2:a1:a2:dzn:phasecor:phasedir:dzcorr" );
1989- tup2.Fill ((float )labelEQ, params[i]. refz , params [j].refz , param1.GetDzDs (), param2.GetDzDs (), param1.GetQPt (), param2.GetQPt (), param1.GetSinPhi (), param2.GetSinPhi (), trk1.GetAlpha (), trk2.GetAlpha (), dznormalized, phasecorr, phasecorrdirection, dzcorr);
2009+ tup2.Fill ((float )labelEQ, refZI, candidates [j].refz , param1.GetDzDs (), param2.GetDzDs (), param1.GetQPt (), param2.GetQPt (), param1.GetSinPhi (), param2.GetSinPhi (), trk1.GetAlpha (), trk2.GetAlpha (), dznormalized, phasecorr, phasecorrdirection, dzcorr);
19902010 }
19912011 /* if (EQ) {
19922012 dropped[j] = true;
@@ -2000,9 +2020,9 @@ GPUd() void GPUTPCGMMerger::MergeLoopersMain(int32_t nBlocks, int32_t nThreads,
20002020 }*/
20012021#endif
20022022 if (EQ) {
2003- mMergedTracks [params [j].id ].SetMergedLooperUnconnected (true );
2023+ mMergedTracks [candidates [j].id ].SetMergedLooperUnconnected (true );
20042024 if (CAMath::Abs (param2.GetQPt () * Param ().qptB5Scaler ) >= Param ().rec .tpc .rejectQPtB5 ) {
2005- mMergedTracks [params [i].id ].SetMergedLooperUnconnected (true );
2025+ mMergedTracks [candidates [i].id ].SetMergedLooperUnconnected (true );
20062026 }
20072027 }
20082028 }
0 commit comments