214 GPUd()
bool recalculatePCAWithErrors(
int cand = 0);
220 auto m = calcPCACovMatrix(
cand);
221 return {
static_cast<float>(
m(0, 0)),
static_cast<float>(
m(1, 0)),
static_cast<float>(
m(1, 1)),
static_cast<float>(
m(2, 0)),
static_cast<float>(
m(2, 1)),
static_cast<float>(
m(2, 2))};
230 GPUdi()
void setPropagateToPCA(
bool v = true) { mPropagateToPCA =
v; }
231 GPUdi()
void setMaxIter(
int n = 20) { mMaxIter =
n > 2 ?
n : 2; }
232 GPUdi()
void setMaxR(
float r = 200.) { mMaxR2 =
r *
r; }
233 GPUdi()
void setMaxDZIni(
float d = 4.) { mMaxDZIni = d; }
234 GPUdi()
void setMaxDXYIni(
float d = 4.) { mMaxDXYIni = d > 0 ? d : 1e9; }
235 GPUdi()
void setMaxChi2(
float chi2 = 999.) { mMaxChi2 = chi2; }
237 GPUdi()
void setMinParamChange(
float x = 1e-3) { mMinParamChange =
x > 1e-4 ?
x : 1.e-4; }
238 GPUdi()
void setMinRelChi2Change(
float r = 0.9) { mMinRelChi2Change =
r > 0.1 ?
r : 999.; }
239 GPUdi()
void setUseAbsDCA(
bool v) { mUseAbsDCA =
v; }
240 GPUdi()
void setWeightedFinalPCA(
bool v) { mWeightedFinalPCA =
v; }
241 GPUdi()
void setMaxDistance2ToMerge(
float v) { mMaxDist2ToMergeSeeds =
v; }
243 GPUdi()
void setUsePropagator(
bool v) { mUsePropagator =
v; }
244 GPUdi()
void setRefitWithMatCorr(
bool v) { mRefitWithMatCorr =
v; }
245 GPUdi()
void setMaxSnp(
float s) { mMaxSnp =
s; }
246 GPUdi()
void setMaxStep(
float s) { mMaxStep =
s; }
247 GPUdi()
void setMinXSeed(
float x) { mMinXSeed =
x; }
248 GPUdi()
void setCollinear(
bool isCollinear) { mIsCollinear = isCollinear; }
250 GPUdi()
int getNCandidates()
const {
return mCurHyp; }
251 GPUdi()
int getMaxIter()
const {
return mMaxIter; }
252 GPUdi() float getMaxR()
const {
return o2::gpu::GPUCommonMath::Sqrt(mMaxR2); }
253 GPUdi() float getMaxDZIni()
const {
return mMaxDZIni; }
254 GPUdi() float getMaxDXYIni()
const {
return mMaxDXYIni; }
255 GPUdi() float getMaxChi2()
const {
return mMaxChi2; }
256 GPUdi() float getMinParamChange()
const {
return mMinParamChange; }
257 GPUdi() float getBz()
const {
return mBz; }
258 GPUdi() float getMaxDistance2ToMerge()
const {
return mMaxDist2ToMergeSeeds; }
259 GPUdi() bool getUseAbsDCA()
const {
return mUseAbsDCA; }
260 GPUdi() bool getWeightedFinalPCA()
const {
return mWeightedFinalPCA; }
261 GPUdi() bool getPropagateToPCA()
const {
return mPropagateToPCA; }
263 GPUdi() bool getUsePropagator()
const {
return mUsePropagator; }
264 GPUdi() bool getRefitWithMatCorr()
const {
return mRefitWithMatCorr; }
265 GPUdi() float getMaxSnp()
const {
return mMaxSnp; }
266 GPUdi() float getMasStep()
const {
return mMaxStep; }
267 GPUdi() float getMinXSeed()
const {
return mMinXSeed; }
269 template <
class... Tr>
273 GPUdi()
int getFitterID()
const {
return mFitterID; }
275 GPUdi() size_t getCallID()
const {
return mCallID; }
279 GPUd()
bool calcInverseWeight();
280 GPUd()
void calcResidDerivatives();
281 GPUd()
void calcResidDerivativesNoErr();
283 GPUd()
void calcChi2Derivatives();
284 GPUd()
void calcChi2DerivativesNoErr();
287 GPUd()
void calcTrackResiduals();
288 GPUd()
void calcTrackDerivatives();
289 GPUd()
double calcChi2() const;
290 GPUd()
double calcChi2NoErr() const;
293 GPUd()
bool minimizeChi2NoErr();
294 GPUd()
bool roughDZCut() const;
295 GPUd()
bool closerToAlternative() const;
297 GPUd()
bool propagateParamToX(
o2::track::TrackPar&
t,
float x);
304 GPUd() float getTrackX(
int i,
int cand = 0)
const {
return getTrackPos(
i,
cand)[0]; }
306 GPUd() MatStd3D getTrackRotMatrix(
int i) const
310 mat(0, 0) = mat(1, 1) = mTrAux[
i].c;
311 mat(0, 1) = -mTrAux[
i].s;
312 mat(1, 0) = mTrAux[
i].s;
316 GPUd() MatSym3D getTrackCovMatrix(
int i,
int cand = 0) const
318 const auto& trc = mCandTr[mOrder[
cand]][
i];
320 mat(0, 0) = trc.getSigmaY2() * XerrFactor;
321 mat(1, 1) = trc.getSigmaY2();
322 mat(2, 2) = trc.getSigmaZ2();
323 mat(2, 1) = trc.getSigmaZY();
328 template <
class T,
class... Tr>
331#ifndef GPUCA_GPUCODE_DEVICE
332 static_assert(std::is_convertible<T, Track>(),
"Wrong track type");
341 mAllowAltPreference =
true;
343 mPropFailed.fill(
false);
344 mTrPropDone.fill(
false);
359 mLoggerBadCov.clear();
360 mLoggerBadInv.clear();
361 mLoggerBadProp.clear();
369 std::array<std::array<Vec3D, N>, N> mDResidDx;
372 std::array<std::array<Vec3D, N>, N> mD2ResidDx2;
378 std::array<TrackAuxPar, N> mTrAux;
379 CrossInfo mCrossings;
381 std::array<ArrTrackCovI, MAXHYP> mTrcEInv;
382 std::array<ArrTrack, MAXHYP> mCandTr;
383 std::array<ArrTrCoef, MAXHYP> mTrCFVT;
384 std::array<ArrTrDer, MAXHYP> mTrDer;
385 std::array<ArrTrPos, MAXHYP> mTrPos;
386 std::array<ArrTrPos, MAXHYP> mTrRes;
387 std::array<Vec3D, MAXHYP> mPCA;
388 std::array<float, MAXHYP> mChi2 = {0};
389 std::array<int, MAXHYP> mNIters;
390 std::array<bool, MAXHYP> mTrPropDone{};
391 std::array<bool, MAXHYP> mPropFailed{};
392 LogLogThrottler mLoggerBadCov{};
393 LogLogThrottler mLoggerBadInv{};
394 LogLogThrottler mLoggerBadProp{};
396 std::array<int, MAXHYP> mOrder{0};
399 int mCrossIDAlt = -1;
401 std::array<FitStatus, MAXHYP> mFitStatus{};
402 bool mAllowAltPreference =
true;
403 bool mUseAbsDCA =
false;
404 bool mWeightedFinalPCA =
false;
405 bool mPropagateToPCA =
true;
406 bool mUsePropagator =
false;
407 bool mRefitWithMatCorr =
false;
408 bool mIsCollinear =
false;
412 float mMaxR2 = 200. * 200.;
413 float mMinXSeed = -50.;
414 float mMaxDZIni = 4.;
415 float mMaxDXYIni = 4.;
416 float mMinParamChange = 1e-3;
417 float mMinRelChi2Change = 0.9;
418 float mMaxChi2 = 100;
419 float mMaxDist2ToMergeSeeds = 1.;
420 float mMaxSnp = 0.95;
421 float mMaxStep = 2.0;
428template <
int N,
typename... Args>
429template <
class... Tr>
434 static_assert(
sizeof...(args) == N,
"incorrect number of input tracks");
437 for (
int i = 0;
i < N;
i++) {
438 mTrAux[
i].set(*mOrigTrPtr[
i], mBz);
440 if (!mCrossings.set(mTrAux[0], *mOrigTrPtr[0], mTrAux[1], *mOrigTrPtr[1], mMaxDXYIni, mIsCollinear)) {
441 mFitStatus[mCurHyp] = FitStatus::NoCrossing;
447 if (mCrossings.nDCA == MAXHYP) {
448 auto dst2 = (mCrossings.xDCA[0] - mCrossings.xDCA[1]) * (mCrossings.xDCA[0] - mCrossings.xDCA[1]) +
449 (mCrossings.yDCA[0] - mCrossings.yDCA[1]) * (mCrossings.yDCA[0] - mCrossings.yDCA[1]);
450 if (dst2 < mMaxDist2ToMergeSeeds) {
452 mCrossings.xDCA[0] = 0.5 * (mCrossings.xDCA[0] + mCrossings.xDCA[1]);
453 mCrossings.yDCA[0] = 0.5 * (mCrossings.yDCA[0] + mCrossings.yDCA[1]);
457 for (
int ic = 0; ic < mCrossings.nDCA; ic++) {
459 if (mCrossings.xDCA[ic] * mCrossings.xDCA[ic] + mCrossings.yDCA[ic] * mCrossings.yDCA[ic] > mMaxR2) {
460 mFitStatus[mCurHyp] = FitStatus::RejRadius;
464 mCrossIDAlt = (mCrossings.nDCA == 2 && mAllowAltPreference) ? 1 - ic : -1;
465 mPCA[mCurHyp][0] = mCrossings.xDCA[ic];
466 mPCA[mCurHyp][1] = mCrossings.yDCA[ic];
468 if (mUseAbsDCA ? minimizeChi2NoErr() : minimizeChi2()) {
469 mOrder[mCurHyp] = mCurHyp;
470 if (mPropagateToPCA && !propagateTracksToVertex(mCurHyp)) {
477 for (
int i = mCurHyp;
i--;) {
478 for (
int j =
i;
j--;) {
479 if (mChi2[mOrder[
i]] < mChi2[mOrder[
j]]) {
480 o2::gpu::GPUCommonMath::Swap(mOrder[
i], mOrder[
j]);
484 if (mUseAbsDCA && mWeightedFinalPCA) {
485 for (
int i = mCurHyp;
i--;) {
486 recalculatePCAWithErrors(
i);
493template <
int N,
typename... Args>
494GPUd() bool DCAFitterN<N, Args...>::calcPCACoefs()
497 if (!calcInverseWeight()) {
498 mFitStatus[mCurHyp] = FitStatus::FailInvWeight;
501 for (
int i = N;
i--;) {
502 const auto& taux = mTrAux[
i];
503 const auto& tcov = mTrcEInv[mCurHyp][
i];
505 miei[0][0] = taux.c * tcov.sxx;
506 miei[0][1] = -taux.s * tcov.syy;
507 miei[0][2] = -taux.s * tcov.syz;
508 miei[1][0] = taux.s * tcov.sxx;
509 miei[1][1] = taux.c * tcov.syy;
510 miei[1][2] = taux.c * tcov.syz;
512 miei[2][1] = tcov.syz;
513 miei[2][2] = tcov.szz;
514 mTrCFVT[mCurHyp][
i] = mWeightInv * miei;
520template <
int N,
typename... Args>
521GPUd() bool DCAFitterN<N, Args...>::calcInverseWeight()
524 auto* arrmat = mWeightInv.Array();
525 memset(arrmat, 0,
sizeof(mWeightInv));
532 for (
int i = N;
i--;) {
533 const auto& taux = mTrAux[
i];
534 const auto& tcov = mTrcEInv[mCurHyp][
i];
535 arrmat[XX] += taux.cc * tcov.sxx + taux.ss * tcov.syy;
536 arrmat[XY] += taux.cs * (tcov.sxx - tcov.syy);
537 arrmat[XZ] += -taux.s * tcov.syz;
538 arrmat[YY] += taux.cc * tcov.syy + taux.ss * tcov.sxx;
539 arrmat[YZ] += taux.c * tcov.syz;
540 arrmat[ZZ] += tcov.szz;
543 return mWeightInv.Invert();
547template <
int N,
typename... Args>
548GPUd()
void DCAFitterN<N, Args...>::calcResidDerivatives()
552 for (
int i = N;
i--;) {
553 const auto& taux = mTrAux[
i];
554 for (
int j = N;
j--;) {
555 const auto& matT = mTrCFVT[mCurHyp][
j];
556 const auto& trDx = mTrDer[mCurHyp][
j];
557 auto& dr1 = mDResidDx[
i][
j];
558 auto& dr2 = mD2ResidDx2[
i][
j];
560 matMT[0][0] = taux.c * matT[0][0] + taux.s * matT[1][0];
561 matMT[0][1] = taux.c * matT[0][1] + taux.s * matT[1][1];
562 matMT[0][2] = taux.c * matT[0][2] + taux.s * matT[1][2];
563 matMT[1][0] = -taux.s * matT[0][0] + taux.c * matT[1][0];
564 matMT[1][1] = -taux.s * matT[0][1] + taux.c * matT[1][1];
565 matMT[1][2] = -taux.s * matT[0][2] + taux.c * matT[1][2];
566 matMT[2][0] = matT[2][0];
567 matMT[2][1] = matT[2][1];
568 matMT[2][2] = matT[2][2];
571 dr1[0] = -(matMT[0][0] + matMT[0][1] * trDx.dydx + matMT[0][2] * trDx.dzdx);
572 dr1[1] = -(matMT[1][0] + matMT[1][1] * trDx.dydx + matMT[1][2] * trDx.dzdx);
573 dr1[2] = -(matMT[2][0] + matMT[2][1] * trDx.dydx + matMT[2][2] * trDx.dzdx);
576 dr2[0] = -(matMT[0][1] * trDx.d2ydx2 + matMT[0][2] * trDx.d2zdx2);
577 dr2[1] = -(matMT[1][1] * trDx.d2ydx2 + matMT[1][2] * trDx.d2zdx2);
578 dr2[2] = -(matMT[2][1] * trDx.d2ydx2 + matMT[2][2] * trDx.d2zdx2);
585 dr2[1] += trDx.d2ydx2;
586 dr2[2] += trDx.d2zdx2;
593template <
int N,
typename... Args>
594GPUd()
void DCAFitterN<N, Args...>::calcResidDerivativesNoErr()
597 constexpr double NInv1 = 1. - NInv;
598 for (
int i = N;
i--;) {
599 const auto& trDxi = mTrDer[mCurHyp][
i];
600 auto& dr1ii = mDResidDx[
i][
i];
601 auto& dr2ii = mD2ResidDx2[
i][
i];
603 dr1ii[1] = NInv1 * trDxi.dydx;
604 dr1ii[2] = NInv1 * trDxi.dzdx;
607 dr2ii[1] = NInv1 * trDxi.d2ydx2;
608 dr2ii[2] = NInv1 * trDxi.d2zdx2;
610 for (
int j =
i;
j--;) {
611 auto& dr1ij = mDResidDx[
i][
j];
612 auto& dr1ji = mDResidDx[
j][
i];
613 const auto& trDxj = mTrDer[mCurHyp][
j];
614 auto cij = mCosDif[
i][
j], sij = mSinDif[
i][
j];
617 dr1ij[0] = -(cij + sij * trDxj.dydx);
618 dr1ij[1] = -(-sij + cij * trDxj.dydx);
619 dr1ij[2] = -trDxj.dzdx * NInv;
622 dr1ji[0] = -(cij - sij * trDxi.dydx);
623 dr1ji[1] = -(sij + cij * trDxi.dydx);
624 dr1ji[2] = -trDxi.dzdx * NInv;
626 auto& dr2ij = mD2ResidDx2[
i][
j];
627 auto& dr2ji = mD2ResidDx2[
j][
i];
629 dr2ij[0] = -sij * trDxj.d2ydx2;
630 dr2ij[1] = -cij * trDxj.d2ydx2;
631 dr2ij[2] = -trDxj.d2zdx2 * NInv;
634 dr2ji[0] = sij * trDxi.d2ydx2;
635 dr2ji[1] = -cij * trDxi.d2ydx2;
636 dr2ji[2] = -trDxi.d2zdx2 * NInv;
643template <
int N,
typename... Args>
644GPUd()
void DCAFitterN<N, Args...>::calcRMatrices()
647 for (
int i = N;
i--;) {
648 const auto& mi = mTrAux[
i];
649 for (
int j =
i;
j--;) {
650 const auto& mj = mTrAux[
j];
651 mCosDif[
i][
j] = (mi.c * mj.c + mi.s * mj.s) * NInv;
652 mSinDif[
i][
j] = (mi.s * mj.c - mi.c * mj.s) * NInv;
658template <
int N,
typename... Args>
659GPUd()
void DCAFitterN<N, Args...>::calcChi2Derivatives()
662 std::array<std::array<Vec3D, N>, N> covIDrDx;
665 for (
int i = N;
i--;) {
666 auto& dchi1 = mDChi2Dx[
i];
668 for (
int j = N;
j--;) {
669 const auto&
res = mTrRes[mCurHyp][
j];
670 const auto& covI = mTrcEInv[mCurHyp][
j];
671 const auto& dr1 = mDResidDx[
j][
i];
672 auto& cidr = covIDrDx[
i][
j];
673 cidr[0] = covI.sxx * dr1[0];
674 cidr[1] = covI.syy * dr1[1] + covI.syz * dr1[2];
675 cidr[2] = covI.syz * dr1[1] + covI.szz * dr1[2];
681 for (
int i = N;
i--;) {
682 for (
int j =
i + 1;
j--;) {
683 auto& dchi2 = mD2Chi2Dx2[
i][
j];
685 for (
int k = N; k--;) {
686 const auto& dr1j = mDResidDx[k][
j];
687 const auto& cidrkj = covIDrDx[
i][k];
690 const auto&
res = mTrRes[mCurHyp][k];
691 const auto& covI = mTrcEInv[mCurHyp][k];
692 const auto& dr2ij = mD2ResidDx2[k][
j];
693 dchi2 +=
res[0] * covI.sxx * dr2ij[0] +
res[1] * (covI.syy * dr2ij[1] + covI.syz * dr2ij[2]) +
res[2] * (covI.syz * dr2ij[1] + covI.szz * dr2ij[2]);
701template <
int N,
typename... Args>
702GPUd()
void DCAFitterN<N, Args...>::calcChi2DerivativesNoErr()
705 for (
int i = N;
i--;) {
706 auto& dchi1 = mDChi2Dx[
i];
708 for (
int j = N;
j--;) {
709 const auto&
res = mTrRes[mCurHyp][
j];
710 const auto& dr1 = mDResidDx[
j][
i];
714 auto& dchi2 = mD2Chi2Dx2[
i][
j];
716 for (
int k = N; k--;) {
725template <
int N,
typename... Args>
726GPUd()
void DCAFitterN<N, Args...>::calcPCA()
729 mPCA[mCurHyp] = mTrCFVT[mCurHyp][N - 1] * mTrPos[mCurHyp][N - 1];
730 for (
int i = N - 1;
i--;) {
731 mPCA[mCurHyp] += mTrCFVT[mCurHyp][
i] * mTrPos[mCurHyp][
i];
736template <
int N,
typename... Args>
737GPUd() bool DCAFitterN<N, Args...>::recalculatePCAWithErrors(
int cand)
740 if (isPropagateTracksToVertexDone(cand) && !propagateTracksToVertex(cand)) {
743 int saveCurHyp = mCurHyp;
744 mCurHyp = mOrder[cand];
746 for (
int i = N;
i--;) {
747 if (!mTrcEInv[mCurHyp][
i].
set(mCandTr[mCurHyp][
i], XerrFactor)) {
748 if (mLoggerBadCov.needToLog()) {
750 printf(
"fitter %d: error (%ld muted): overrode invalid track covariance from %s\n",
751 mFitterID, mLoggerBadCov.evCount, mCandTr[mCurHyp][
i].asString().c_str());
753 printf(
"fitter %d: error (%ld muted): overrode invalid track covariance cyy:%e czz:%e cyz:%e\n",
754 mFitterID, mLoggerBadCov.evCount, mCandTr[mCurHyp][
i].getSigmaY2(), mCandTr[mCurHyp][
i].getSigmaZ2(), mCandTr[mCurHyp][
i].getSigmaZY());
757 mFitStatus[mCurHyp] = FitStatus::FailInvCov;
758 if (mBadCovPolicy == Discard) {
760 }
else if (mBadCovPolicy == OverrideAndFlag) {
761 mPropFailed[mCurHyp] =
true;
765 if (!calcPCACoefs()) {
766 mCurHyp = saveCurHyp;
770 auto oldPCA = mPCA[mOrder[cand]];
772 mCurHyp = saveCurHyp;
777template <
int N,
typename... Args>
778GPUd()
void DCAFitterN<N, Args...>::calcPCANoErr()
781 auto& pca = mPCA[mCurHyp];
782 o2::math_utils::rotateZd(mTrPos[mCurHyp][N - 1][0], mTrPos[mCurHyp][N - 1][1], pca[0], pca[1], mTrAux[N - 1].s, mTrAux[N - 1].
c);
784 pca[2] = mTrPos[mCurHyp][N - 1][2];
785 for (
int i = N - 1;
i--;) {
791 pca[2] += mTrPos[mCurHyp][
i][2];
799template <
int N,
typename... Args>
800GPUd()
o2::math_utils::SMatrix<
double, 3, 3,
o2::math_utils::MatRepSym<
double, 3>>
DCAFitterN<N, Args...>::calcPCACovMatrix(
int cand)
const
805 for (
int i = N;
i--;) {
809 if (covTr.Invert()) {
814 if (nAdded && covm.Invert()) {
819 for (
int i = N;
i--;) {
826template <
int N,
typename... Args>
827GPUd()
void DCAFitterN<N, Args...>::calcTrackResiduals()
831 for (
int i = N;
i--;) {
832 mTrRes[mCurHyp][
i] = mTrPos[mCurHyp][
i];
833 vtxLoc = mPCA[mCurHyp];
835 mTrRes[mCurHyp][
i] -= vtxLoc;
840template <
int N,
typename... Args>
841GPUdi()
void DCAFitterN<N, Args...>::calcTrackDerivatives()
844 for (
int i = N;
i--;) {
845 mTrDer[mCurHyp][
i].set(mCandTr[mCurHyp][
i], mBz);
850template <
int N,
typename... Args>
851GPUdi() double DCAFitterN<N, Args...>::calcChi2()
const
855 for (
int i = N;
i--;) {
856 const auto&
res = mTrRes[mCurHyp][
i];
857 const auto& covI = mTrcEInv[mCurHyp][
i];
858 chi2 +=
res[0] *
res[0] * covI.sxx +
res[1] *
res[1] * covI.syy +
res[2] *
res[2] * covI.szz + 2. *
res[1] *
res[2] * covI.syz;
864template <
int N,
typename... Args>
865GPUdi() double DCAFitterN<N, Args...>::calcChi2NoErr()
const
869 for (
int i = N;
i--;) {
870 const auto&
res = mTrRes[mCurHyp][
i];
877template <
int N,
typename... Args>
878GPUd() bool DCAFitterN<N, Args...>::correctTracks(const VecND& corrX)
881 for (
int i = N;
i--;) {
882 const auto& trDer = mTrDer[mCurHyp][
i];
883 auto dx2h = 0.5 * corrX[
i] * corrX[
i];
884 mTrPos[mCurHyp][
i][0] -= corrX[
i];
885 mTrPos[mCurHyp][
i][1] -= trDer.dydx * corrX[
i] - dx2h * trDer.d2ydx2;
886 mTrPos[mCurHyp][
i][2] -= trDer.dzdx * corrX[
i] - dx2h * trDer.d2zdx2;
892template <
int N,
typename... Args>
893GPUd() bool DCAFitterN<N, Args...>::propagateTracksToVertex(
int icand)
896 int ord = mOrder[icand];
897 if (mTrPropDone[ord]) {
902 if (mRefitWithMatCorr) {
903 int curHypSav = mCurHyp, curCrosIDAlt = mCrossIDAlt;
906 auto restore = [
this, curHypSav, curCrosIDAlt]() { this->mCurHyp = curHypSav; this->mCrossIDAlt = curCrosIDAlt; };
907 if (!(mUseAbsDCA ? minimizeChi2NoErr() : minimizeChi2())) {
914 for (
int i = N;
i--;) {
915 if (mUseAbsDCA || mUsePropagator || mMatCorr != o2::base::Propagator::MatCorrType::USEMatCorrNONE) {
916 mCandTr[ord][
i] = *mOrigTrPtr[
i];
918 auto x = mTrAux[
i].c * mPCA[ord][0] + mTrAux[
i].s * mPCA[ord][1];
919 if (!propagateToX(mCandTr[ord][
i],
x)) {
924 mTrPropDone[ord] =
true;
929template <
int N,
typename... Args>
933 int ord = mOrder[icand];
935 if (!mTrPropDone[ord]) {
936 auto x = mTrAux[
i].c * mPCA[ord][0] + mTrAux[
i].s * mPCA[ord][1];
937 if (!propagateParamToX(trc,
x)) {
945template <
int N,
typename... Args>
946GPUdi() double DCAFitterN<N, Args...>::getAbsMax(const VecND&
v)
949 for (
int i = N;
i--;) {
950 auto vai = o2::gpu::GPUCommonMath::Abs(
v[
i]);
959template <
int N,
typename... Args>
960GPUd() bool DCAFitterN<N, Args...>::minimizeChi2()
963 for (
int i = N;
i--;) {
964 mCandTr[mCurHyp][
i] = *mOrigTrPtr[
i];
965 auto x = mTrAux[
i].c * mPCA[mCurHyp][0] + mTrAux[
i].s * mPCA[mCurHyp][1];
967 mFitStatus[mCurHyp] = FitStatus::RejTrackX;
970 if (!propagateToX(mCandTr[mCurHyp][
i],
x)) {
973 setTrackPos(mTrPos[mCurHyp][
i], mCandTr[mCurHyp][
i]);
974 if (!mTrcEInv[mCurHyp][
i].
set(mCandTr[mCurHyp][
i], XerrFactor)) {
975 if (mLoggerBadCov.needToLog()) {
977 printf(
"fitter %d: error (%ld muted): overrode invalid track covariance from %s\n",
978 mFitterID, mLoggerBadCov.evCount, mCandTr[mCurHyp][
i].asString().c_str());
980 printf(
"fitter %d: error (%ld muted): overrode invalid track covariance cyy:%e czz:%e cyz:%e\n",
981 mFitterID, mLoggerBadCov.evCount, mCandTr[mCurHyp][
i].getSigmaY2(), mCandTr[mCurHyp][
i].getSigmaZ2(), mCandTr[mCurHyp][
i].getSigmaZY());
984 mFitStatus[mCurHyp] = FitStatus::FailInvCov;
985 if (mBadCovPolicy == Discard) {
987 }
else if (mBadCovPolicy == OverrideAndFlag) {
988 mPropFailed[mCurHyp] =
true;
993 if (mMaxDZIni > 0 && !roughDZCut()) {
994 mFitStatus[mCurHyp] = FitStatus::RejTrackRoughZ;
998 if (!calcPCACoefs()) {
1002 calcTrackResiduals();
1003 float chi2Upd, chi2 = calcChi2();
1005 calcTrackDerivatives();
1006 calcResidDerivatives();
1007 calcChi2Derivatives();
1010 if (!mD2Chi2Dx2.Invert()) {
1011 if (mLoggerBadInv.needToLog()) {
1012 printf(
"fitter %d: error (%ld muted): Inversion failed\n", mFitterID, mLoggerBadCov.evCount);
1014 mFitStatus[mCurHyp] = FitStatus::FailInv2ndDeriv;
1017 VecND dx = mD2Chi2Dx2 * mDChi2Dx;
1018 if (!correctTracks(dx)) {
1019 mFitStatus[mCurHyp] = FitStatus::FailCorrTracks;
1023 if (mCrossIDAlt >= 0 && closerToAlternative()) {
1024 mFitStatus[mCurHyp] = FitStatus::FailCloserAlt;
1025 mAllowAltPreference =
false;
1028 calcTrackResiduals();
1029 chi2Upd = calcChi2();
1030 if (getAbsMax(dx) < mMinParamChange || chi2Upd > chi2 * mMinRelChi2Change) {
1032 mFitStatus[mCurHyp] = FitStatus::Converged;
1036 }
while (++mNIters[mCurHyp] < mMaxIter);
1037 if (mNIters[mCurHyp] == mMaxIter) {
1038 mFitStatus[mCurHyp] = FitStatus::MaxIter;
1041 mChi2[mCurHyp] = chi2 * NInv;
1042 if (mChi2[mCurHyp] >= mMaxChi2) {
1043 mFitStatus[mCurHyp] = FitStatus::RejChi2Max;
1050template <
int N,
typename... Args>
1051GPUd() bool DCAFitterN<N, Args...>::minimizeChi2NoErr()
1055 for (
int i = N;
i--;) {
1056 mCandTr[mCurHyp][
i] = *mOrigTrPtr[
i];
1057 auto x = mTrAux[
i].c * mPCA[mCurHyp][0] + mTrAux[
i].s * mPCA[mCurHyp][1];
1058 if (
x < mMinXSeed) {
1059 mFitStatus[mCurHyp] = FitStatus::RejTrackX;
1062 if (!propagateParamToX(mCandTr[mCurHyp][
i],
x)) {
1065 setTrackPos(mTrPos[mCurHyp][
i], mCandTr[mCurHyp][
i]);
1067 if (mMaxDZIni > 0 && !roughDZCut()) {
1068 mFitStatus[mCurHyp] = FitStatus::RejTrackRoughZ;
1073 calcTrackResiduals();
1074 float chi2Upd, chi2 = calcChi2NoErr();
1076 calcTrackDerivatives();
1077 calcResidDerivativesNoErr();
1078 calcChi2DerivativesNoErr();
1081 if (!mD2Chi2Dx2.Invert()) {
1082 if (mLoggerBadInv.needToLog()) {
1083 printf(
"fitter %d: error (%ld muted): Inversion failed\n", mFitterID, mLoggerBadCov.evCount);
1085 mFitStatus[mCurHyp] = FitStatus::FailInv2ndDeriv;
1088 VecND dx = mD2Chi2Dx2 * mDChi2Dx;
1089 if (!correctTracks(dx)) {
1090 mFitStatus[mCurHyp] = FitStatus::FailCorrTracks;
1094 if (mCrossIDAlt >= 0 && closerToAlternative()) {
1095 mFitStatus[mCurHyp] = FitStatus::FailCloserAlt;
1096 mAllowAltPreference =
false;
1099 calcTrackResiduals();
1100 chi2Upd = calcChi2NoErr();
1101 if (getAbsMax(dx) < mMinParamChange || chi2Upd > chi2 * mMinRelChi2Change) {
1103 mFitStatus[mCurHyp] = FitStatus::Converged;
1107 }
while (++mNIters[mCurHyp] < mMaxIter);
1108 if (mNIters[mCurHyp] == mMaxIter) {
1109 mFitStatus[mCurHyp] = FitStatus::MaxIter;
1112 mChi2[mCurHyp] = chi2 * NInv;
1113 if (mChi2[mCurHyp] >= mMaxChi2) {
1114 mFitStatus[mCurHyp] = FitStatus::RejChi2Max;
1121template <
int N,
typename... Args>
1122GPUd() bool DCAFitterN<N, Args...>::roughDZCut()
const
1126 for (
int i = N; accept &&
i--;) {
1127 for (
int j =
i;
j--;) {
1128 if (o2::gpu::GPUCommonMath::Abs(mCandTr[mCurHyp][
i].getZ() - mCandTr[mCurHyp][
j].getZ()) > mMaxDZIni) {
1138template <
int N,
typename... Args>
1139GPUd() bool DCAFitterN<N, Args...>::closerToAlternative()
const
1142 auto dxCur = mPCA[mCurHyp][0] - mCrossings.xDCA[mCrossIDCur], dyCur = mPCA[mCurHyp][1] - mCrossings.yDCA[mCrossIDCur];
1143 auto dxAlt = mPCA[mCurHyp][0] - mCrossings.xDCA[mCrossIDAlt], dyAlt = mPCA[mCurHyp][1] - mCrossings.yDCA[mCrossIDAlt];
1144 return dxCur * dxCur + dyCur * dyCur > dxAlt * dxAlt + dyAlt * dyAlt;
1148template <
int N,
typename... Args>
1151#ifndef GPUCA_GPUCODE_DEVICE
1152 LOG(info) << N <<
"-prong vertex fitter in " << (mUseAbsDCA ?
"abs." :
"weighted") <<
" distance minimization mode, collinear tracks mode: " << (mIsCollinear ?
"ON" :
"OFF");
1153 LOG(info) <<
"Bz: " << mBz <<
" MaxIter: " << mMaxIter <<
" MaxChi2: " << mMaxChi2 <<
" MatCorrType: " <<
int(mMatCorr);
1154 LOG(info) <<
"Stopping condition: Max.param change < " << mMinParamChange <<
" Rel.Chi2 change > " << mMinRelChi2Change;
1155 LOG(info) <<
"Discard candidates for : Rvtx > " << getMaxR() <<
" DZ between tracks > " << mMaxDZIni;
1156 LOG(info) <<
"PropagateToPCA:" << mPropagateToPCA <<
" WeightedFinalPCA:" << mWeightedFinalPCA <<
" UsePropagator:" << mUsePropagator <<
" RefitWithMatCorr:" << mRefitWithMatCorr;
1158 for (
int i = 0;
i < mCrossings.nDCA;
i++) {
1159 rep += fmt::format(
"seed{}:{}/{} ",
i, mTrPropDone[
i], mPropFailed[
i]);
1161 LOG(info) <<
"Last call: NCand:" << mCurHyp <<
" from " << mCrossings.nDCA <<
" seeds, prop.done/failed: " << rep;
1164 printf(
"%d-prong vertex fitter in abs. distance minimization mode\n", N);
1166 printf(
"%d-prong vertex fitter in weighted distance minimization mode\n", N);
1168 printf(
"Bz: %1.f MaxIter: %3.d MaxChi2: %2.3f\n", mBz, mMaxIter, mMaxChi2);
1169 printf(
"Stopping condition: Max.param change < %2.3f Rel.Chi2 change > %2.3f\n", mMinParamChange, mMinRelChi2Change);
1170 printf(
"Discard candidates for : Rvtx > %2.3f DZ between tracks > %2.3f\n", getMaxR(), mMaxDZIni);
1175template <
int N,
typename... Args>
1176GPUd()
o2::track::
TrackParCov DCAFitterN<N, Args...>::createParentTrackParCov(
int cand,
bool sectorAlpha)
const
1178 const auto& trP = getTrack(0, cand);
1179 const auto& trN = getTrack(1, cand);
1180 std::array<float, 21> covV = {0.};
1181 std::array<float, 3> pvecV = {0.};
1183 for (
int it = 0; it < N; it++) {
1184 const auto& trc = getTrack(it, cand);
1185 std::array<float, 3> pvecT = {0.};
1186 std::array<float, 21> covT = {0.};
1187 trc.getPxPyPzGlo(pvecT);
1188 trc.getCovXYZPxPyPzGlo(covT);
1189 constexpr int MomInd[6] = {9, 13, 14, 18, 19, 20};
1190 for (
int i = 0;
i < 6;
i++) {
1191 covV[MomInd[
i]] += covT[MomInd[
i]];
1193 for (
int i = 0;
i < 3;
i++) {
1194 pvecV[
i] += pvecT[
i];
1196 q += trc.getCharge();
1198 auto covVtxV = calcPCACovMatrix(cand);
1199 covV[0] = covVtxV(0, 0);
1200 covV[1] = covVtxV(1, 0);
1201 covV[2] = covVtxV(1, 1);
1202 covV[3] = covVtxV(2, 0);
1203 covV[4] = covVtxV(2, 1);
1204 covV[5] = covVtxV(2, 2);
1209template <
int N,
typename... Args>
1210GPUd()
o2::track::TrackPar DCAFitterN<N, Args...>::createParentTrackPar(
int cand,
bool sectorAlpha)
const
1212 const auto& trP = getTrack(0, cand);
1213 const auto& trN = getTrack(1, cand);
1214 const auto& wvtx = getPCACandidate(cand);
1215 std::array<float, 3> pvecV = {0.};
1217 for (
int it = 0; it < N; it++) {
1218 const auto& trc = getTrack(it, cand);
1219 std::array<float, 3> pvecT = {0.};
1220 trc.getPxPyPzGlo(pvecT);
1221 for (
int i = 0;
i < 3;
i++) {
1222 pvecV[
i] += pvecT[
i];
1224 q += trc.getCharge();
1226 const std::array<float, 3>
vertex = {(
float)wvtx[0], (
float)wvtx[1], (
float)wvtx[2]};
1231template <
int N,
typename... Args>
1232GPUdi() bool DCAFitterN<N, Args...>::propagateParamToX(
o2::track::TrackPar& t,
float x)
1235 if (mUsePropagator || mMatCorr != o2::base::Propagator::MatCorrType::USEMatCorrNONE) {
1236#ifndef GPUCA_GPUCODE
1240 res = t.propagateParamTo(
x, mBz);
1243 mFitStatus[mCurHyp] = FitStatus::FailProp;
1244 mPropFailed[mCurHyp] =
true;
1245 if (mLoggerBadProp.needToLog()) {
1246#ifndef GPUCA_GPUCODE
1247 printf(
"fitter %d: error (%ld muted): propagation failed for %s\n", mFitterID, mLoggerBadProp.evCount, t.asString().c_str());
1249 printf(
"fitter %d: error (%ld muted): propagation failed\n", mFitterID, mLoggerBadProp.evCount);
1257template <
int N,
typename... Args>
1261 if (mUsePropagator || mMatCorr != o2::base::Propagator::MatCorrType::USEMatCorrNONE) {
1262#ifndef GPUCA_GPUCODE
1266 res = t.propagateTo(
x, mBz);
1269 mFitStatus[mCurHyp] = FitStatus::FailProp;
1270 mPropFailed[mCurHyp] =
true;
1271 if (mLoggerBadProp.needToLog()) {
1272#ifndef GPUCA_GPUCODE
1273 printf(
"fitter %d: error (%ld muted): propagation failed for %s\n", mFitterID, mLoggerBadProp.evCount, t.asString().c_str());
1275 printf(
"fitter %d: error (%ld muted): propagation failed\n", mFitterID, mLoggerBadProp.evCount);
1287template <
typename Fitter>
1288void print(
const int nBlocks,
const int nThreads, Fitter& ft);
1290template <
typename Fitter,
class... Tr>
1291int process(
const int nBlocks,
const int nThreads, Fitter&, Tr&... args);
1293template <
class Fitter,
class... Tr>
1294void processBulk(
const int nBlocks,
const int nThreads,
const int nBatches, std::vector<Fitter>& fitters, std::vector<int>& results, std::vector<Tr>&... args);