@@ -1887,13 +1887,12 @@ GPUd() void GPUTPCGMMerger::Finalize2(int32_t nBlocks, int32_t nThreads, int32_t
18871887
18881888GPUd () void GPUTPCGMMerger::MergeLoopersInit(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread)
18891889{
1890- return ; // FIXME: !!!!
1891- const float lowPtThresh = Param ().rec .tpc .rejectQPtB5 * 1 .1f ; // Might need to merge tracks above the threshold with parts below the threshold
1890+ const float lowPtThresh = Param ().rec .tpc .rejectQPtB5 * 1 .1f ; // Might need to merge tracks above the threshold with parts below the rejection threshold
18921891 for (uint32_t i = get_global_id (0 ); i < mMemory ->nMergedTracks ; i += get_global_size (0 )) {
18931892 const auto & trk = mMergedTracks [i];
18941893 const auto & p = trk.GetParam ();
18951894 const float qptabs = CAMath::Abs (p.GetQPt ());
1896- if (trk.NClusters () && qptabs * Param ().qptB5Scaler > 5 .f && qptabs * Param ().qptB5Scaler <= lowPtThresh) {
1895+ if (trk.OK () && trk. NClusters () && trk. Leg () == 0 && qptabs * Param ().qptB5Scaler > 5 .f && qptabs * Param ().qptB5Scaler <= lowPtThresh) {
18971896 const int32_t sector = mClusters [trk.FirstClusterRef () + trk.NClusters () - 1 ].sector ;
18981897 const float refz = p.GetZ () + GetConstantMem ()->calibObjects .fastTransformHelper ->getCorrMap ()->convVertexTimeToZOffset (sector, p.GetTOffset (), Param ().continuousMaxTimeBin ) + (trk.CSide () ? -100 : 100 );
18991898 float sinA, cosA;
@@ -1942,12 +1941,12 @@ GPUd() void GPUTPCGMMerger::MergeLoopersSort(int32_t nBlocks, int32_t nThreads,
19421941
19431942GPUd () void GPUTPCGMMerger::MergeLoopersMain(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread)
19441943{
1945- const MergeLooperParam* params = mLooperCandidates ;
1944+ const MergeLooperParam* candidates = mLooperCandidates ;
19461945
19471946#if GPUCA_MERGE_LOOPER_MC && !defined(GPUCA_GPUCODE)
19481947 std::vector<int64_t > paramLabels (mMemory ->nLooperMatchCandidates );
19491948 for (uint32_t i = 0 ; i < mMemory ->nLooperMatchCandidates ; i++) {
1950- paramLabels[i] = GetTrackLabel (mMergedTracks [params [i].id ]);
1949+ paramLabels[i] = GetTrackLabel (mMergedTracks [candidates [i].id ]);
19511950 }
19521951 /* std::vector<bool> dropped(mMemory->nLooperMatchCandidates);
19531952 std::vector<bool> droppedMC(mMemory->nLooperMatchCandidates);
@@ -1961,26 +1960,47 @@ GPUd() void GPUTPCGMMerger::MergeLoopersMain(int32_t nBlocks, int32_t nThreads,
19611960 for (uint32_t i = get_global_id (0 ); i < mMemory ->nLooperMatchCandidates ; i += get_global_size (0 )) {
19621961 for (uint32_t j = i + 1 ; j < mMemory ->nLooperMatchCandidates ; j++) {
19631962 // int32_t bs = 0;
1964- if (CAMath::Abs (params[j].refz ) > CAMath::Abs (params[i].refz ) + 100 .f ) {
1963+ assert (CAMath::Abs (candidates[i].refz ) <= CAMath::Abs (candidates[j].refz ));
1964+ if (CAMath::Abs (candidates[j].refz ) > CAMath::Abs (candidates[i].refz ) + 100 .f ) {
19651965 break ;
19661966 }
1967- const float d2xy = CAMath::Sum2 (params [i].x - params [j].x , params [i].y - params [j].y );
1967+ const float d2xy = CAMath::Sum2 (candidates [i].x - candidates [j].x , candidates [i].y - candidates [j].y );
19681968 if (d2xy > 15 .f ) {
19691969 // bs |= 1;
19701970 continue ;
19711971 }
1972- const auto & trk1 = mMergedTracks [params[i].id ];
1973- const auto & trk2 = mMergedTracks [params[j].id ];
1972+
1973+ const GPUTPCGMMergedTrack* trkI = &mMergedTracks [candidates[i].id ];
1974+ float refZI = candidates[i].refz ;
1975+ {
1976+ const auto * tmp = trkI;
1977+ while (tmp->PrevSegment () >= 0 ) {
1978+ const auto * next = &mMergedTracks [tmp->PrevSegment ()];
1979+ if (next == trkI) {
1980+ break ;
1981+ }
1982+ tmp = next;
1983+ }
1984+ if (tmp != trkI && tmp->CSide () == trkI->CSide () && CAMath::Abs (tmp->GetParam ().GetZ ()) > CAMath::Abs (trkI->GetParam ().GetZ ())) {
1985+ float tmpRefZ = refZI + tmp->GetParam ().GetZ () - trkI->GetParam ().GetZ ();
1986+ if (CAMath::Abs (tmpRefZ) < CAMath::Abs (candidates[j].refz ) && CAMath::Abs (tmpRefZ) > CAMath::Abs (refZI)) {
1987+ trkI = tmp;
1988+ refZI = tmpRefZ;
1989+ }
1990+ }
1991+ };
1992+ const auto & trk1 = *trkI;
1993+ const auto & trk2 = mMergedTracks [candidates[j].id ];
19741994 const auto & param1 = trk1.GetParam ();
19751995 const auto & param2 = trk2.GetParam ();
19761996 if (CAMath::Abs (param1.GetDzDs ()) > 0 .03f && CAMath::Abs (param2.GetDzDs ()) > 0 .03f && param1.GetDzDs () * param2.GetDzDs () * param1.GetQPt () * param2.GetQPt () < 0 ) {
19771997 // bs |= 2;
19781998 continue ;
19791999 }
19802000
1981- 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 ())));
2001+ 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 ())));
19822002 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 ;
1983- const float phasecorrdirection = (params [j].refz * param1.GetQPt () * param1.GetDzDs ()) > 0 ? 1 : -1 ;
2003+ const float phasecorrdirection = (candidates [j].refz * param1.GetQPt () * param1.GetDzDs ()) > 0 ? 1 : -1 ;
19842004 const float dzcorr = dznormalized + phasecorr * phasecorrdirection;
19852005 const bool sameside = !(trk1.CSide () ^ trk2.CSide ());
19862006 const float dzcorrlimit[4 ] = {sameside ? 0 .018f : 0 .012f , sameside ? 0 .12f : 0 .025f , 0 .14f , 0 .15f };
@@ -2009,11 +2029,11 @@ GPUd() void GPUTPCGMMerger::MergeLoopersMain(int32_t nBlocks, int32_t nThreads,
20092029 const int64_t label2 = paramLabels[j];
20102030 bool labelEQ = label1 != -1 && label1 == label2;
20112031 if (1 || EQ || labelEQ) {
2012- // 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);
2032+ // 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);
20132033 static auto & tup = GPUROOTDump<TNtuple>::get (" mergeloopers" , " labeleq:sides:d2xy:tgl1:tgl2:qpt1:qpt2:dz:dzcorr:dtgl:dqpt:dznorm:bs" );
2014- 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);
2034+ 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);
20152035 static auto tup2 = GPUROOTDump<TNtuple>::getNew (" mergeloopers2" , " labeleq:refz1:refz2:tgl1:tgl2:qpt1:qpt2:snp1:snp2:a1:a2:dzn:phasecor:phasedir:dzcorr" );
2016- 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);
2036+ 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);
20172037 }
20182038 /* if (EQ) {
20192039 dropped[j] = true;
@@ -2027,9 +2047,9 @@ GPUd() void GPUTPCGMMerger::MergeLoopersMain(int32_t nBlocks, int32_t nThreads,
20272047 }*/
20282048#endif
20292049 if (EQ) {
2030- mMergedTracks [params [j].id ].SetMergedLooperUnconnected (true );
2050+ mMergedTracks [candidates [j].id ].SetMergedLooperUnconnected (true );
20312051 if (CAMath::Abs (param2.GetQPt () * Param ().qptB5Scaler ) >= Param ().rec .tpc .rejectQPtB5 ) {
2032- mMergedTracks [params [i].id ].SetMergedLooperUnconnected (true );
2052+ mMergedTracks [candidates [i].id ].SetMergedLooperUnconnected (true );
20332053 }
20342054 }
20352055 }
0 commit comments