212 GPUd()
bool recalculatePCAWithErrors(
int cand = 0);
218 auto m = calcPCACovMatrix(
cand);
219 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))};
228 GPUdi()
void setPropagateToPCA(
bool v = true) { mPropagateToPCA =
v; }
229 GPUdi()
void setMaxIter(
int n = 20) { mMaxIter =
n > 2 ?
n : 2; }
230 GPUdi()
void setMaxR(
float r = 200.) { mMaxR2 =
r *
r; }
231 GPUdi()
void setMaxDZIni(
float d = 4.) { mMaxDZIni = d; }
232 GPUdi()
void setMaxDXYIni(
float d = 4.) { mMaxDXYIni = d > 0 ? d : 1e9; }
233 GPUdi()
void setMaxChi2(
float chi2 = 999.) { mMaxChi2 = chi2; }
235 GPUdi()
void setMinParamChange(
float x = 1e-3) { mMinParamChange =
x > 1e-4 ?
x : 1.e-4; }
236 GPUdi()
void setMinRelChi2Change(
float r = 0.9) { mMinRelChi2Change =
r > 0.1 ?
r : 999.; }
237 GPUdi()
void setUseAbsDCA(
bool v) { mUseAbsDCA =
v; }
238 GPUdi()
void setWeightedFinalPCA(
bool v) { mWeightedFinalPCA =
v; }
239 GPUdi()
void setMaxDistance2ToMerge(
float v) { mMaxDist2ToMergeSeeds =
v; }
241 GPUdi()
void setUsePropagator(
bool v) { mUsePropagator =
v; }
242 GPUdi()
void setRefitWithMatCorr(
bool v) { mRefitWithMatCorr =
v; }
243 GPUdi()
void setMaxSnp(
float s) { mMaxSnp =
s; }
244 GPUdi()
void setMaxStep(
float s) { mMaxStep =
s; }
245 GPUdi()
void setMinXSeed(
float x) { mMinXSeed =
x; }
246 GPUdi()
void setCollinear(
bool isCollinear) { mIsCollinear = isCollinear; }
248 GPUdi()
int getNCandidates()
const {
return mCurHyp; }
249 GPUdi()
int getMaxIter()
const {
return mMaxIter; }
250 GPUdi() float getMaxR()
const {
return o2::gpu::GPUCommonMath::Sqrt(mMaxR2); }
251 GPUdi() float getMaxDZIni()
const {
return mMaxDZIni; }
252 GPUdi() float getMaxDXYIni()
const {
return mMaxDXYIni; }
253 GPUdi() float getMaxChi2()
const {
return mMaxChi2; }
254 GPUdi() float getMinParamChange()
const {
return mMinParamChange; }
255 GPUdi() float getBz()
const {
return mBz; }
256 GPUdi() float getMaxDistance2ToMerge()
const {
return mMaxDist2ToMergeSeeds; }
257 GPUdi() bool getUseAbsDCA()
const {
return mUseAbsDCA; }
258 GPUdi() bool getWeightedFinalPCA()
const {
return mWeightedFinalPCA; }
259 GPUdi() bool getPropagateToPCA()
const {
return mPropagateToPCA; }
261 GPUdi() bool getUsePropagator()
const {
return mUsePropagator; }
262 GPUdi() bool getRefitWithMatCorr()
const {
return mRefitWithMatCorr; }
263 GPUdi() float getMaxSnp()
const {
return mMaxSnp; }
264 GPUdi() float getMasStep()
const {
return mMaxStep; }
265 GPUdi() float getMinXSeed()
const {
return mMinXSeed; }
267 template <
class... Tr>
271 GPUdi()
int getFitterID()
const {
return mFitterID; }
273 GPUdi() size_t getCallID()
const {
return mCallID; }
277 GPUd()
bool calcInverseWeight();
278 GPUd()
void calcResidDerivatives();
279 GPUd()
void calcResidDerivativesNoErr();
281 GPUd()
void calcChi2Derivatives();
282 GPUd()
void calcChi2DerivativesNoErr();
285 GPUd()
void calcTrackResiduals();
286 GPUd()
void calcTrackDerivatives();
287 GPUd()
double calcChi2() const;
288 GPUd()
double calcChi2NoErr() const;
291 GPUd()
bool minimizeChi2NoErr();
292 GPUd()
bool roughDZCut() const;
293 GPUd()
bool closerToAlternative() const;
295 GPUd()
bool propagateParamToX(
o2::track::TrackPar&
t,
float x);
302 GPUd() float getTrackX(
int i,
int cand = 0)
const {
return getTrackPos(
i,
cand)[0]; }
304 GPUd() MatStd3D getTrackRotMatrix(
int i) const
308 mat(0, 0) = mat(1, 1) = mTrAux[
i].c;
309 mat(0, 1) = -mTrAux[
i].s;
310 mat(1, 0) = mTrAux[
i].s;
314 GPUd() MatSym3D getTrackCovMatrix(
int i,
int cand = 0) const
316 const auto& trc = mCandTr[mOrder[
cand]][
i];
318 mat(0, 0) = trc.getSigmaY2() * XerrFactor;
319 mat(1, 1) = trc.getSigmaY2();
320 mat(2, 2) = trc.getSigmaZ2();
321 mat(2, 1) = trc.getSigmaZY();
326 template <
class T,
class... Tr>
329#ifndef GPUCA_GPUCODE_DEVICE
330 static_assert(std::is_convertible<T, Track>(),
"Wrong track type");
339 mAllowAltPreference =
true;
341 mPropFailed.fill(
false);
342 mTrPropDone.fill(
false);
357 mLoggerBadCov.clear();
358 mLoggerBadInv.clear();
359 mLoggerBadProp.clear();
377 CrossInfo mCrossings;
390 LogLogThrottler mLoggerBadCov{};
391 LogLogThrottler mLoggerBadInv{};
392 LogLogThrottler mLoggerBadProp{};
397 int mCrossIDAlt = -1;
400 bool mAllowAltPreference =
true;
401 bool mUseAbsDCA =
false;
402 bool mWeightedFinalPCA =
false;
403 bool mPropagateToPCA =
true;
404 bool mUsePropagator =
false;
405 bool mRefitWithMatCorr =
false;
406 bool mIsCollinear =
false;
410 float mMaxR2 = 200. * 200.;
411 float mMinXSeed = -50.;
412 float mMaxDZIni = 4.;
413 float mMaxDXYIni = 4.;
414 float mMinParamChange = 1e-3;
415 float mMinRelChi2Change = 0.9;
416 float mMaxChi2 = 100;
417 float mMaxDist2ToMergeSeeds = 1.;
418 float mMaxSnp = 0.95;
419 float mMaxStep = 2.0;
426template <
int N,
typename... Args>
427template <
class... Tr>
432 static_assert(
sizeof...(args) == N,
"incorrect number of input tracks");
435 for (
int i = 0;
i < N;
i++) {
436 mTrAux[
i].set(*mOrigTrPtr[
i], mBz);
438 if (!mCrossings.set(mTrAux[0], *mOrigTrPtr[0], mTrAux[1], *mOrigTrPtr[1], mMaxDXYIni, mIsCollinear)) {
439 mFitStatus[mCurHyp] = FitStatus::NoCrossing;
445 if (mCrossings.nDCA == MAXHYP) {
446 auto dst2 = (mCrossings.xDCA[0] - mCrossings.xDCA[1]) * (mCrossings.xDCA[0] - mCrossings.xDCA[1]) +
447 (mCrossings.yDCA[0] - mCrossings.yDCA[1]) * (mCrossings.yDCA[0] - mCrossings.yDCA[1]);
448 if (dst2 < mMaxDist2ToMergeSeeds) {
450 mCrossings.xDCA[0] = 0.5 * (mCrossings.xDCA[0] + mCrossings.xDCA[1]);
451 mCrossings.yDCA[0] = 0.5 * (mCrossings.yDCA[0] + mCrossings.yDCA[1]);
455 for (
int ic = 0; ic < mCrossings.nDCA; ic++) {
457 if (mCrossings.xDCA[ic] * mCrossings.xDCA[ic] + mCrossings.yDCA[ic] * mCrossings.yDCA[ic] > mMaxR2) {
458 mFitStatus[mCurHyp] = FitStatus::RejRadius;
462 mCrossIDAlt = (mCrossings.nDCA == 2 && mAllowAltPreference) ? 1 - ic : -1;
463 mPCA[mCurHyp][0] = mCrossings.xDCA[ic];
464 mPCA[mCurHyp][1] = mCrossings.yDCA[ic];
466 if (mUseAbsDCA ? minimizeChi2NoErr() : minimizeChi2()) {
467 mOrder[mCurHyp] = mCurHyp;
468 if (mPropagateToPCA && !propagateTracksToVertex(mCurHyp)) {
475 for (
int i = mCurHyp;
i--;) {
476 for (
int j =
i;
j--;) {
477 if (mChi2[mOrder[
i]] < mChi2[mOrder[
j]]) {
478 o2::gpu::GPUCommonMath::Swap(mOrder[
i], mOrder[
j]);
482 if (mUseAbsDCA && mWeightedFinalPCA) {
483 for (
int i = mCurHyp;
i--;) {
484 recalculatePCAWithErrors(
i);
491template <
int N,
typename... Args>
492GPUd() bool DCAFitterN<N, Args...>::calcPCACoefs()
495 if (!calcInverseWeight()) {
496 mFitStatus[mCurHyp] = FitStatus::FailInvWeight;
499 for (
int i = N;
i--;) {
500 const auto& taux = mTrAux[
i];
501 const auto& tcov = mTrcEInv[mCurHyp][
i];
503 miei[0][0] = taux.c * tcov.sxx;
504 miei[0][1] = -taux.s * tcov.syy;
505 miei[0][2] = -taux.s * tcov.syz;
506 miei[1][0] = taux.s * tcov.sxx;
507 miei[1][1] = taux.c * tcov.syy;
508 miei[1][2] = taux.c * tcov.syz;
510 miei[2][1] = tcov.syz;
511 miei[2][2] = tcov.szz;
512 mTrCFVT[mCurHyp][
i] = mWeightInv * miei;
518template <
int N,
typename... Args>
519GPUd() bool DCAFitterN<N, Args...>::calcInverseWeight()
522 auto* arrmat = mWeightInv.Array();
523 memset(arrmat, 0,
sizeof(mWeightInv));
530 for (
int i = N;
i--;) {
531 const auto& taux = mTrAux[
i];
532 const auto& tcov = mTrcEInv[mCurHyp][
i];
533 arrmat[XX] += taux.cc * tcov.sxx + taux.ss * tcov.syy;
534 arrmat[XY] += taux.cs * (tcov.sxx - tcov.syy);
535 arrmat[XZ] += -taux.s * tcov.syz;
536 arrmat[YY] += taux.cc * tcov.syy + taux.ss * tcov.sxx;
537 arrmat[YZ] += taux.c * tcov.syz;
538 arrmat[ZZ] += tcov.szz;
541 return mWeightInv.Invert();
545template <
int N,
typename... Args>
546GPUd()
void DCAFitterN<N, Args...>::calcResidDerivatives()
550 for (
int i = N;
i--;) {
551 const auto& taux = mTrAux[
i];
552 for (
int j = N;
j--;) {
553 const auto& matT = mTrCFVT[mCurHyp][
j];
554 const auto& trDx = mTrDer[mCurHyp][
j];
555 auto& dr1 = mDResidDx[
i][
j];
556 auto& dr2 = mD2ResidDx2[
i][
j];
558 matMT[0][0] = taux.c * matT[0][0] + taux.s * matT[1][0];
559 matMT[0][1] = taux.c * matT[0][1] + taux.s * matT[1][1];
560 matMT[0][2] = taux.c * matT[0][2] + taux.s * matT[1][2];
561 matMT[1][0] = -taux.s * matT[0][0] + taux.c * matT[1][0];
562 matMT[1][1] = -taux.s * matT[0][1] + taux.c * matT[1][1];
563 matMT[1][2] = -taux.s * matT[0][2] + taux.c * matT[1][2];
564 matMT[2][0] = matT[2][0];
565 matMT[2][1] = matT[2][1];
566 matMT[2][2] = matT[2][2];
569 dr1[0] = -(matMT[0][0] + matMT[0][1] * trDx.dydx + matMT[0][2] * trDx.dzdx);
570 dr1[1] = -(matMT[1][0] + matMT[1][1] * trDx.dydx + matMT[1][2] * trDx.dzdx);
571 dr1[2] = -(matMT[2][0] + matMT[2][1] * trDx.dydx + matMT[2][2] * trDx.dzdx);
574 dr2[0] = -(matMT[0][1] * trDx.d2ydx2 + matMT[0][2] * trDx.d2zdx2);
575 dr2[1] = -(matMT[1][1] * trDx.d2ydx2 + matMT[1][2] * trDx.d2zdx2);
576 dr2[2] = -(matMT[2][1] * trDx.d2ydx2 + matMT[2][2] * trDx.d2zdx2);
583 dr2[1] += trDx.d2ydx2;
584 dr2[2] += trDx.d2zdx2;
591template <
int N,
typename... Args>
592GPUd()
void DCAFitterN<N, Args...>::calcResidDerivativesNoErr()
595 constexpr double NInv1 = 1. - NInv;
596 for (
int i = N;
i--;) {
597 const auto& trDxi = mTrDer[mCurHyp][
i];
598 auto& dr1ii = mDResidDx[
i][
i];
599 auto& dr2ii = mD2ResidDx2[
i][
i];
601 dr1ii[1] = NInv1 * trDxi.dydx;
602 dr1ii[2] = NInv1 * trDxi.dzdx;
605 dr2ii[1] = NInv1 * trDxi.d2ydx2;
606 dr2ii[2] = NInv1 * trDxi.d2zdx2;
608 for (
int j =
i;
j--;) {
609 auto& dr1ij = mDResidDx[
i][
j];
610 auto& dr1ji = mDResidDx[
j][
i];
611 const auto& trDxj = mTrDer[mCurHyp][
j];
612 auto cij = mCosDif[
i][
j], sij = mSinDif[
i][
j];
615 dr1ij[0] = -(cij + sij * trDxj.dydx);
616 dr1ij[1] = -(-sij + cij * trDxj.dydx);
617 dr1ij[2] = -trDxj.dzdx * NInv;
620 dr1ji[0] = -(cij - sij * trDxi.dydx);
621 dr1ji[1] = -(sij + cij * trDxi.dydx);
622 dr1ji[2] = -trDxi.dzdx * NInv;
624 auto& dr2ij = mD2ResidDx2[
i][
j];
625 auto& dr2ji = mD2ResidDx2[
j][
i];
627 dr2ij[0] = -sij * trDxj.d2ydx2;
628 dr2ij[1] = -cij * trDxj.d2ydx2;
629 dr2ij[2] = -trDxj.d2zdx2 * NInv;
632 dr2ji[0] = sij * trDxi.d2ydx2;
633 dr2ji[1] = -cij * trDxi.d2ydx2;
634 dr2ji[2] = -trDxi.d2zdx2 * NInv;
641template <
int N,
typename... Args>
642GPUd()
void DCAFitterN<N, Args...>::calcRMatrices()
645 for (
int i = N;
i--;) {
646 const auto& mi = mTrAux[
i];
647 for (
int j =
i;
j--;) {
648 const auto& mj = mTrAux[
j];
649 mCosDif[
i][
j] = (mi.c * mj.c + mi.s * mj.s) * NInv;
650 mSinDif[
i][
j] = (mi.s * mj.c - mi.c * mj.s) * NInv;
656template <
int N,
typename... Args>
657GPUd()
void DCAFitterN<N, Args...>::calcChi2Derivatives()
663 for (
int i = N;
i--;) {
664 auto& dchi1 = mDChi2Dx[
i];
666 for (
int j = N;
j--;) {
667 const auto&
res = mTrRes[mCurHyp][
j];
668 const auto& covI = mTrcEInv[mCurHyp][
j];
669 const auto& dr1 = mDResidDx[
j][
i];
670 auto& cidr = covIDrDx[
i][
j];
671 cidr[0] = covI.sxx * dr1[0];
672 cidr[1] = covI.syy * dr1[1] + covI.syz * dr1[2];
673 cidr[2] = covI.syz * dr1[1] + covI.szz * dr1[2];
679 for (
int i = N;
i--;) {
680 for (
int j =
i + 1;
j--;) {
681 auto& dchi2 = mD2Chi2Dx2[
i][
j];
683 for (
int k = N; k--;) {
684 const auto& dr1j = mDResidDx[k][
j];
685 const auto& cidrkj = covIDrDx[
i][k];
688 const auto&
res = mTrRes[mCurHyp][k];
689 const auto& covI = mTrcEInv[mCurHyp][k];
690 const auto& dr2ij = mD2ResidDx2[k][
j];
691 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]);
699template <
int N,
typename... Args>
700GPUd()
void DCAFitterN<N, Args...>::calcChi2DerivativesNoErr()
703 for (
int i = N;
i--;) {
704 auto& dchi1 = mDChi2Dx[
i];
706 for (
int j = N;
j--;) {
707 const auto&
res = mTrRes[mCurHyp][
j];
708 const auto& dr1 = mDResidDx[
j][
i];
712 auto& dchi2 = mD2Chi2Dx2[
i][
j];
714 for (
int k = N; k--;) {
723template <
int N,
typename... Args>
724GPUd()
void DCAFitterN<N, Args...>::calcPCA()
727 mPCA[mCurHyp] = mTrCFVT[mCurHyp][N - 1] * mTrPos[mCurHyp][N - 1];
728 for (
int i = N - 1;
i--;) {
729 mPCA[mCurHyp] += mTrCFVT[mCurHyp][
i] * mTrPos[mCurHyp][
i];
734template <
int N,
typename... Args>
735GPUd() bool DCAFitterN<N, Args...>::recalculatePCAWithErrors(
int cand)
738 if (isPropagateTracksToVertexDone(cand) && !propagateTracksToVertex(cand)) {
741 int saveCurHyp = mCurHyp;
742 mCurHyp = mOrder[cand];
744 for (
int i = N;
i--;) {
745 if (!mTrcEInv[mCurHyp][
i].
set(mCandTr[mCurHyp][
i], XerrFactor)) {
746 if (mLoggerBadCov.needToLog()) {
748 printf(
"fitter %d: error (%ld muted): overrode invalid track covariance from %s\n",
749 mFitterID, mLoggerBadCov.evCount, mCandTr[mCurHyp][
i].asString().c_str());
751 printf(
"fitter %d: error (%ld muted): overrode invalid track covariance cyy:%e czz:%e cyz:%e\n",
752 mFitterID, mLoggerBadCov.evCount, mCandTr[mCurHyp][
i].getSigmaY2(), mCandTr[mCurHyp][
i].getSigmaZ2(), mCandTr[mCurHyp][
i].getSigmaZY());
755 mFitStatus[mCurHyp] = FitStatus::FailInvCov;
756 if (mBadCovPolicy == Discard) {
758 }
else if (mBadCovPolicy == OverrideAndFlag) {
759 mPropFailed[mCurHyp] =
true;
763 if (!calcPCACoefs()) {
764 mCurHyp = saveCurHyp;
768 auto oldPCA = mPCA[mOrder[cand]];
770 mCurHyp = saveCurHyp;
775template <
int N,
typename... Args>
776GPUd()
void DCAFitterN<N, Args...>::calcPCANoErr()
779 auto& pca = mPCA[mCurHyp];
780 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);
782 pca[2] = mTrPos[mCurHyp][N - 1][2];
783 for (
int i = N - 1;
i--;) {
789 pca[2] += mTrPos[mCurHyp][
i][2];
797template <
int N,
typename... Args>
798GPUd()
o2::math_utils::SMatrix<
double, 3, 3,
o2::math_utils::MatRepSym<
double, 3>>
DCAFitterN<N, Args...>::calcPCACovMatrix(
int cand)
const
803 for (
int i = N;
i--;) {
807 if (covTr.Invert()) {
812 if (nAdded && covm.Invert()) {
817 for (
int i = N;
i--;) {
824template <
int N,
typename... Args>
825GPUd()
void DCAFitterN<N, Args...>::calcTrackResiduals()
829 for (
int i = N;
i--;) {
830 mTrRes[mCurHyp][
i] = mTrPos[mCurHyp][
i];
831 vtxLoc = mPCA[mCurHyp];
833 mTrRes[mCurHyp][
i] -= vtxLoc;
838template <
int N,
typename... Args>
839GPUdi()
void DCAFitterN<N, Args...>::calcTrackDerivatives()
842 for (
int i = N;
i--;) {
843 mTrDer[mCurHyp][
i].set(mCandTr[mCurHyp][
i], mBz);
848template <
int N,
typename... Args>
849GPUdi() double DCAFitterN<N, Args...>::calcChi2()
const
853 for (
int i = N;
i--;) {
854 const auto&
res = mTrRes[mCurHyp][
i];
855 const auto& covI = mTrcEInv[mCurHyp][
i];
856 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;
862template <
int N,
typename... Args>
863GPUdi() double DCAFitterN<N, Args...>::calcChi2NoErr()
const
867 for (
int i = N;
i--;) {
868 const auto&
res = mTrRes[mCurHyp][
i];
875template <
int N,
typename... Args>
876GPUd() bool DCAFitterN<N, Args...>::correctTracks(const VecND& corrX)
879 for (
int i = N;
i--;) {
880 const auto& trDer = mTrDer[mCurHyp][
i];
881 auto dx2h = 0.5 * corrX[
i] * corrX[
i];
882 mTrPos[mCurHyp][
i][0] -= corrX[
i];
883 mTrPos[mCurHyp][
i][1] -= trDer.dydx * corrX[
i] - dx2h * trDer.d2ydx2;
884 mTrPos[mCurHyp][
i][2] -= trDer.dzdx * corrX[
i] - dx2h * trDer.d2zdx2;
890template <
int N,
typename... Args>
891GPUd() bool DCAFitterN<N, Args...>::propagateTracksToVertex(
int icand)
894 int ord = mOrder[icand];
895 if (mTrPropDone[ord]) {
900 if (mRefitWithMatCorr) {
901 int curHypSav = mCurHyp, curCrosIDAlt = mCrossIDAlt;
904 auto restore = [
this, curHypSav, curCrosIDAlt]() { this->mCurHyp = curHypSav; this->mCrossIDAlt = curCrosIDAlt; };
905 if (!(mUseAbsDCA ? minimizeChi2NoErr() : minimizeChi2())) {
912 for (
int i = N;
i--;) {
913 if (mUseAbsDCA || mUsePropagator || mMatCorr != o2::base::Propagator::MatCorrType::USEMatCorrNONE) {
914 mCandTr[ord][
i] = *mOrigTrPtr[
i];
916 auto x = mTrAux[
i].c * mPCA[ord][0] + mTrAux[
i].s * mPCA[ord][1];
917 if (!propagateToX(mCandTr[ord][
i],
x)) {
922 mTrPropDone[ord] =
true;
927template <
int N,
typename... Args>
931 int ord = mOrder[icand];
933 if (!mTrPropDone[ord]) {
934 auto x = mTrAux[
i].c * mPCA[ord][0] + mTrAux[
i].s * mPCA[ord][1];
935 if (!propagateParamToX(trc,
x)) {
943template <
int N,
typename... Args>
944GPUdi() double DCAFitterN<N, Args...>::getAbsMax(const VecND&
v)
947 for (
int i = N;
i--;) {
948 auto vai = o2::gpu::GPUCommonMath::Abs(
v[
i]);
957template <
int N,
typename... Args>
958GPUd() bool DCAFitterN<N, Args...>::minimizeChi2()
961 for (
int i = N;
i--;) {
962 mCandTr[mCurHyp][
i] = *mOrigTrPtr[
i];
963 auto x = mTrAux[
i].c * mPCA[mCurHyp][0] + mTrAux[
i].s * mPCA[mCurHyp][1];
965 mFitStatus[mCurHyp] = FitStatus::RejTrackX;
968 if (!propagateToX(mCandTr[mCurHyp][
i],
x)) {
971 setTrackPos(mTrPos[mCurHyp][
i], mCandTr[mCurHyp][
i]);
972 if (!mTrcEInv[mCurHyp][
i].
set(mCandTr[mCurHyp][
i], XerrFactor)) {
973 if (mLoggerBadCov.needToLog()) {
975 printf(
"fitter %d: error (%ld muted): overrode invalid track covariance from %s\n",
976 mFitterID, mLoggerBadCov.evCount, mCandTr[mCurHyp][
i].asString().c_str());
978 printf(
"fitter %d: error (%ld muted): overrode invalid track covariance cyy:%e czz:%e cyz:%e\n",
979 mFitterID, mLoggerBadCov.evCount, mCandTr[mCurHyp][
i].getSigmaY2(), mCandTr[mCurHyp][
i].getSigmaZ2(), mCandTr[mCurHyp][
i].getSigmaZY());
982 mFitStatus[mCurHyp] = FitStatus::FailInvCov;
983 if (mBadCovPolicy == Discard) {
985 }
else if (mBadCovPolicy == OverrideAndFlag) {
986 mPropFailed[mCurHyp] =
true;
991 if (mMaxDZIni > 0 && !roughDZCut()) {
992 mFitStatus[mCurHyp] = FitStatus::RejTrackX;
996 if (!calcPCACoefs()) {
1000 calcTrackResiduals();
1001 float chi2Upd, chi2 = calcChi2();
1003 calcTrackDerivatives();
1004 calcResidDerivatives();
1005 calcChi2Derivatives();
1008 if (!mD2Chi2Dx2.Invert()) {
1009 if (mLoggerBadInv.needToLog()) {
1010 printf(
"fitter %d: error (%ld muted): Inversion failed\n", mFitterID, mLoggerBadCov.evCount);
1012 mFitStatus[mCurHyp] = FitStatus::FailInv2ndDeriv;
1015 VecND dx = mD2Chi2Dx2 * mDChi2Dx;
1016 if (!correctTracks(dx)) {
1017 mFitStatus[mCurHyp] = FitStatus::FailCorrTracks;
1021 if (mCrossIDAlt >= 0 && closerToAlternative()) {
1022 mFitStatus[mCurHyp] = FitStatus::FailCloserAlt;
1023 mAllowAltPreference =
false;
1026 calcTrackResiduals();
1027 chi2Upd = calcChi2();
1028 if (getAbsMax(dx) < mMinParamChange || chi2Upd > chi2 * mMinRelChi2Change) {
1030 mFitStatus[mCurHyp] = FitStatus::Converged;
1034 }
while (++mNIters[mCurHyp] < mMaxIter);
1035 if (mNIters[mCurHyp] == mMaxIter) {
1036 mFitStatus[mCurHyp] = FitStatus::MaxIter;
1039 mChi2[mCurHyp] = chi2 * NInv;
1040 if (mChi2[mCurHyp] >= mMaxChi2) {
1041 mFitStatus[mCurHyp] = FitStatus::RejChi2Max;
1048template <
int N,
typename... Args>
1049GPUd() bool DCAFitterN<N, Args...>::minimizeChi2NoErr()
1053 for (
int i = N;
i--;) {
1054 mCandTr[mCurHyp][
i] = *mOrigTrPtr[
i];
1055 auto x = mTrAux[
i].c * mPCA[mCurHyp][0] + mTrAux[
i].s * mPCA[mCurHyp][1];
1056 if (
x < mMinXSeed) {
1057 mFitStatus[mCurHyp] = FitStatus::RejTrackX;
1060 if (!propagateParamToX(mCandTr[mCurHyp][
i],
x)) {
1063 setTrackPos(mTrPos[mCurHyp][
i], mCandTr[mCurHyp][
i]);
1065 if (mMaxDZIni > 0 && !roughDZCut()) {
1066 mFitStatus[mCurHyp] = FitStatus::RejTrackX;
1071 calcTrackResiduals();
1072 float chi2Upd, chi2 = calcChi2NoErr();
1074 calcTrackDerivatives();
1075 calcResidDerivativesNoErr();
1076 calcChi2DerivativesNoErr();
1079 if (!mD2Chi2Dx2.Invert()) {
1080 if (mLoggerBadInv.needToLog()) {
1081 printf(
"fitter %d: error (%ld muted): Inversion failed\n", mFitterID, mLoggerBadCov.evCount);
1083 mFitStatus[mCurHyp] = FitStatus::FailInv2ndDeriv;
1086 VecND dx = mD2Chi2Dx2 * mDChi2Dx;
1087 if (!correctTracks(dx)) {
1088 mFitStatus[mCurHyp] = FitStatus::FailCorrTracks;
1092 if (mCrossIDAlt >= 0 && closerToAlternative()) {
1093 mFitStatus[mCurHyp] = FitStatus::FailCloserAlt;
1094 mAllowAltPreference =
false;
1097 calcTrackResiduals();
1098 chi2Upd = calcChi2NoErr();
1099 if (getAbsMax(dx) < mMinParamChange || chi2Upd > chi2 * mMinRelChi2Change) {
1101 mFitStatus[mCurHyp] = FitStatus::Converged;
1105 }
while (++mNIters[mCurHyp] < mMaxIter);
1106 if (mNIters[mCurHyp] == mMaxIter) {
1107 mFitStatus[mCurHyp] = FitStatus::MaxIter;
1110 mChi2[mCurHyp] = chi2 * NInv;
1111 if (mChi2[mCurHyp] >= mMaxChi2) {
1112 mFitStatus[mCurHyp] = FitStatus::RejChi2Max;
1119template <
int N,
typename... Args>
1120GPUd() bool DCAFitterN<N, Args...>::roughDZCut()
const
1124 for (
int i = N; accept &&
i--;) {
1125 for (
int j =
i;
j--;) {
1126 if (o2::gpu::GPUCommonMath::Abs(mCandTr[mCurHyp][
i].getZ() - mCandTr[mCurHyp][
j].getZ()) > mMaxDZIni) {
1136template <
int N,
typename... Args>
1137GPUd() bool DCAFitterN<N, Args...>::closerToAlternative()
const
1140 auto dxCur = mPCA[mCurHyp][0] - mCrossings.xDCA[mCrossIDCur], dyCur = mPCA[mCurHyp][1] - mCrossings.yDCA[mCrossIDCur];
1141 auto dxAlt = mPCA[mCurHyp][0] - mCrossings.xDCA[mCrossIDAlt], dyAlt = mPCA[mCurHyp][1] - mCrossings.yDCA[mCrossIDAlt];
1142 return dxCur * dxCur + dyCur * dyCur > dxAlt * dxAlt + dyAlt * dyAlt;
1146template <
int N,
typename... Args>
1149#ifndef GPUCA_GPUCODE_DEVICE
1150 LOG(info) << N <<
"-prong vertex fitter in " << (mUseAbsDCA ?
"abs." :
"weighted") <<
" distance minimization mode, collinear tracks mode: " << (mIsCollinear ?
"ON" :
"OFF");
1151 LOG(info) <<
"Bz: " << mBz <<
" MaxIter: " << mMaxIter <<
" MaxChi2: " << mMaxChi2 <<
" MatCorrType: " <<
int(mMatCorr);
1152 LOG(info) <<
"Stopping condition: Max.param change < " << mMinParamChange <<
" Rel.Chi2 change > " << mMinRelChi2Change;
1153 LOG(info) <<
"Discard candidates for : Rvtx > " << getMaxR() <<
" DZ between tracks > " << mMaxDZIni;
1154 LOG(info) <<
"PropagateToPCA:" << mPropagateToPCA <<
" WeightedFinalPCA:" << mWeightedFinalPCA <<
" UsePropagator:" << mUsePropagator <<
" RefitWithMatCorr:" << mRefitWithMatCorr;
1156 for (
int i = 0;
i < mCrossings.nDCA;
i++) {
1157 rep += fmt::format(
"seed{}:{}/{} ",
i, mTrPropDone[
i], mPropFailed[
i]);
1159 LOG(info) <<
"Last call: NCand:" << mCurHyp <<
" from " << mCrossings.nDCA <<
" seeds, prop.done/failed: " << rep;
1162 printf(
"%d-prong vertex fitter in abs. distance minimization mode\n", N);
1164 printf(
"%d-prong vertex fitter in weighted distance minimization mode\n", N);
1166 printf(
"Bz: %1.f MaxIter: %3.d MaxChi2: %2.3f\n", mBz, mMaxIter, mMaxChi2);
1167 printf(
"Stopping condition: Max.param change < %2.3f Rel.Chi2 change > %2.3f\n", mMinParamChange, mMinRelChi2Change);
1168 printf(
"Discard candidates for : Rvtx > %2.3f DZ between tracks > %2.3f\n", getMaxR(), mMaxDZIni);
1173template <
int N,
typename... Args>
1174GPUd()
o2::track::
TrackParCov DCAFitterN<N, Args...>::createParentTrackParCov(
int cand,
bool sectorAlpha)
const
1176 const auto& trP = getTrack(0, cand);
1177 const auto& trN = getTrack(1, cand);
1181 for (
int it = 0; it < N; it++) {
1182 const auto& trc = getTrack(it, cand);
1185 trc.getPxPyPzGlo(pvecT);
1186 trc.getCovXYZPxPyPzGlo(covT);
1187 constexpr int MomInd[6] = {9, 13, 14, 18, 19, 20};
1188 for (
int i = 0;
i < 6;
i++) {
1189 covV[MomInd[
i]] += covT[MomInd[
i]];
1191 for (
int i = 0;
i < 3;
i++) {
1192 pvecV[
i] += pvecT[
i];
1194 q += trc.getCharge();
1196 auto covVtxV = calcPCACovMatrix(cand);
1197 covV[0] = covVtxV(0, 0);
1198 covV[1] = covVtxV(1, 0);
1199 covV[2] = covVtxV(1, 1);
1200 covV[3] = covVtxV(2, 0);
1201 covV[4] = covVtxV(2, 1);
1202 covV[5] = covVtxV(2, 2);
1207template <
int N,
typename... Args>
1208GPUd()
o2::track::TrackPar DCAFitterN<N, Args...>::createParentTrackPar(
int cand,
bool sectorAlpha)
const
1210 const auto& trP = getTrack(0, cand);
1211 const auto& trN = getTrack(1, cand);
1212 const auto& wvtx = getPCACandidate(cand);
1215 for (
int it = 0; it < N; it++) {
1216 const auto& trc = getTrack(it, cand);
1218 trc.getPxPyPzGlo(pvecT);
1219 for (
int i = 0;
i < 3;
i++) {
1220 pvecV[
i] += pvecT[
i];
1222 q += trc.getCharge();
1229template <
int N,
typename... Args>
1230GPUdi() bool DCAFitterN<N, Args...>::propagateParamToX(
o2::track::TrackPar& t,
float x)
1233 if (mUsePropagator || mMatCorr != o2::base::Propagator::MatCorrType::USEMatCorrNONE) {
1234#ifndef GPUCA_GPUCODE
1238 res = t.propagateParamTo(
x, mBz);
1241 mFitStatus[mCurHyp] = FitStatus::FailProp;
1242 mPropFailed[mCurHyp] =
true;
1243 if (mLoggerBadProp.needToLog()) {
1244#ifndef GPUCA_GPUCODE
1245 printf(
"fitter %d: error (%ld muted): propagation failed for %s\n", mFitterID, mLoggerBadProp.evCount, t.asString().c_str());
1247 printf(
"fitter %d: error (%ld muted): propagation failed\n", mFitterID, mLoggerBadProp.evCount);
1255template <
int N,
typename... Args>
1259 if (mUsePropagator || mMatCorr != o2::base::Propagator::MatCorrType::USEMatCorrNONE) {
1260#ifndef GPUCA_GPUCODE
1264 res = t.propagateTo(
x, mBz);
1267 mFitStatus[mCurHyp] = FitStatus::FailProp;
1268 mPropFailed[mCurHyp] =
true;
1269 if (mLoggerBadProp.needToLog()) {
1270#ifndef GPUCA_GPUCODE
1271 printf(
"fitter %d: error (%ld muted): propagation failed for %s\n", mFitterID, mLoggerBadProp.evCount, t.asString().c_str());
1273 printf(
"fitter %d: error (%ld muted): propagation failed\n", mFitterID, mLoggerBadProp.evCount);
1285template <
typename Fitter>
1286void print(
const int nBlocks,
const int nThreads, Fitter& ft);
1288template <
typename Fitter,
class... Tr>
1289int process(
const int nBlocks,
const int nThreads, Fitter&, Tr&... args);
1291template <
class Fitter,
class... Tr>
1292void processBulk(
const int nBlocks,
const int nThreads,
const int nBatches, std::vector<Fitter>& fitters, std::vector<int>& results, std::vector<Tr>&... args);