168 auto dx = xToGo - track.getX();
169 int dir = dx > 0.f ? 1 : -1;
174 std::array<value_type, 3>
b{};
175 while (math_utils::detail::abs<value_type>(dx) > Epsilon) {
176 auto step = math_utils::detail::min<value_type>(math_utils::detail::abs<value_type>(dx), maxStep);
180 auto x = track.getX() + step;
181 auto xyz0 = track.getXYZGlo();
182 getFieldXYZ(xyz0, &
b[0]);
184 auto correct = [&track, &xyz0, tofInfo, matCorr, signCorr,
this]() {
186 if (matCorr != MatCorrType::USEMatCorrNONE) {
187 auto xyz1 = track.getXYZGlo();
188 auto mb = this->getMatBudget(matCorr, xyz0, xyz1);
189 if (!track.correctForMaterial(mb.meanX2X0, mb.getXRho(signCorr))) {
193 tofInfo->addStep(mb.length, track.getQ2P2());
194 tofInfo->addX2X0(mb.meanX2X0);
195 tofInfo->addXRho(mb.getXRho(signCorr));
197 }
else if (tofInfo) {
198 auto xyz1 = track.getXYZGlo();
200 tofInfo->addStep(stepV.R(), track.getQ2P2());
205 if (!track.propagateTo(
x,
b)) {
208 if (maxSnp > 0 && math_utils::detail::abs<value_type>(track.getSnp()) >= maxSnp) {
215 dx = xToGo - track.getX();
222template <
typename value_T>
223GPUd() bool
PropagatorImpl<value_T>::PropagateToXBxByBz(TrackPar_t& track, value_type xToGo, value_type maxSnp, value_type maxStep,
237 auto dx = xToGo - track.getX();
238 int dir = dx > 0.f ? 1 : -1;
243 std::array<value_type, 3>
b{};
244 while (math_utils::detail::abs<value_type>(dx) > Epsilon) {
245 auto step = math_utils::detail::min<value_type>(math_utils::detail::abs<value_type>(dx), maxStep);
249 auto x = track.getX() +
step;
250 auto xyz0 = track.getXYZGlo();
251 getFieldXYZ(xyz0, &
b[0]);
253 auto correct = [&track, &xyz0, tofInfo, matCorr, signCorr,
this]() {
255 if (matCorr != MatCorrType::USEMatCorrNONE) {
256 auto xyz1 = track.getXYZGlo();
257 auto mb = this->getMatBudget(matCorr, xyz0, xyz1);
258 if (!track.correctForELoss(((signCorr < 0) ? -mb.length : mb.
length) * mb.meanRho)) {
262 tofInfo->addStep(mb.length, track.getQ2P2());
263 tofInfo->addX2X0(mb.meanX2X0);
264 tofInfo->addXRho(mb.getXRho(signCorr));
266 }
else if (tofInfo) {
267 auto xyz1 = track.getXYZGlo();
269 tofInfo->addStep(stepV.R(), track.getQ2P2());
274 if (!track.propagateParamTo(
x,
b)) {
277 if (maxSnp > 0 && math_utils::detail::abs<value_type>(track.getSnp()) >= maxSnp) {
284 dx = xToGo - track.getX();
291template <
typename value_T>
292GPUd() bool
PropagatorImpl<value_T>::propagateToX(TrackParCov_t& track, value_type xToGo, value_type bZ, value_type maxSnp, value_type maxStep,
306 auto dx = xToGo - track.getX();
307 int dir = dx > 0.f ? 1 : -1;
312 while (math_utils::detail::abs<value_type>(dx) > Epsilon) {
313 auto step = math_utils::detail::min<value_type>(math_utils::detail::abs<value_type>(dx), maxStep);
317 auto x = track.getX() +
step;
318 auto xyz0 = track.getXYZGlo();
319 auto correct = [&track, &xyz0, tofInfo, matCorr, signCorr,
this]() {
321 if (matCorr != MatCorrType::USEMatCorrNONE) {
322 auto xyz1 = track.getXYZGlo();
323 auto mb = this->getMatBudget(matCorr, xyz0, xyz1);
324 if (!track.correctForMaterial(mb.meanX2X0, mb.getXRho(signCorr))) {
328 tofInfo->addStep(mb.length, track.getQ2P2());
329 tofInfo->addX2X0(mb.meanX2X0);
330 tofInfo->addXRho(mb.getXRho(signCorr));
332 }
else if (tofInfo) {
333 auto xyz1 = track.getXYZGlo();
335 tofInfo->addStep(stepV.R(), track.getQ2P2());
339 if (!track.propagateTo(
x, bZ)) {
342 if (maxSnp > 0 && math_utils::detail::abs<value_type>(track.getSnp()) >= maxSnp) {
349 dx = xToGo - track.getX();
356template <
typename value_T>
357GPUd() bool
PropagatorImpl<value_T>::propagateToX(TrackPar_t& track, value_type xToGo, value_type bZ, value_type maxSnp, value_type maxStep,
371 auto dx = xToGo - track.getX();
372 int dir = dx > 0.f ? 1 : -1;
377 while (math_utils::detail::abs<value_type>(dx) > Epsilon) {
378 auto step = math_utils::detail::min<value_type>(math_utils::detail::abs<value_type>(dx), maxStep);
382 auto x = track.getX() +
step;
383 auto xyz0 = track.getXYZGlo();
385 auto correct = [&track, &xyz0, tofInfo, matCorr, signCorr,
this]() {
387 if (matCorr != MatCorrType::USEMatCorrNONE) {
388 auto xyz1 = track.getXYZGlo();
389 auto mb = this->getMatBudget(matCorr, xyz0, xyz1);
390 if (!track.correctForELoss(mb.getXRho(signCorr))) {
394 tofInfo->addStep(mb.length, track.getQ2P2());
395 tofInfo->addX2X0(mb.meanX2X0);
396 tofInfo->addXRho(mb.getXRho(signCorr));
398 }
else if (tofInfo) {
399 auto xyz1 = track.getXYZGlo();
401 tofInfo->addStep(stepV.R(), track.getQ2P2());
406 if (!track.propagateParamTo(
x, bZ)) {
409 if (maxSnp > 0 && math_utils::detail::abs<value_type>(track.getSnp()) >= maxSnp) {
416 dx = xToGo - track.getX();
423template <
typename value_T>
424template <
typename track_T>
425GPUd() bool
PropagatorImpl<value_T>::propagateToR(track_T& track, value_type
r,
bool bzOnly, value_type maxSnp, value_type maxStep,
426 MatCorrType matCorr, track::TrackLTIntegral* tofInfo,
int signCorr)
const
428 const value_T MaxPhiLoc = math_utils::detail::asin<value_T>(maxSnp), MaxPhiLocSafe = 0.95 * MaxPhiLoc;
429 auto bz = getNominalBz();
430 if (math_utils::detail::abs(bz) > constants::math::Almost0) {
433 value_type r0 = math_utils::detail::sqrt<value_T>(track.getX() * track.getX() + track.getY() * track.getY());
434 value_type dr = (
r - r0);
435 value_type rTmp =
r - (math_utils::detail::abs<value_T>(dr) > 1. ? (dr > 0 ? 0.5 : -0.5) : 0.5 * dr);
437 crad.
c = crad.
cc = 1.f;
438 crad.
s = crad.
ss = crad.
cs = 0.f;
440 cross.circlesCrossInfo(crad, traux, 0.);
441 if (cross.
nDCA < 1) {
444 double phiCross[2] = {}, dphi[2] = {};
445 auto curv = track.getCurvature(bz);
446 bool clockwise = curv < 0;
447 auto phiLoc = math_utils::detail::asin<double>(track.getSnp());
448 auto phi0 = phiLoc + track.getAlpha();
450 for (
int i = 0;
i < cross.
nDCA;
i++) {
455 auto normX = double(cross.
yDCA[
i]) - double(traux.yC), normY = -(double(cross.
xDCA[
i]) - double(traux.xC));
460 phiCross[
i] = math_utils::detail::atan2<double>(normY, normX);
462 dphi[
i] = phiCross[
i] - phi0;
469 int sel = cross.
nDCA == 1 ? 0 : (clockwise ? (dphi[0] < dphi[1] ? 0 : 1) : (dphi[1] < dphi[0] ? 0 : 1));
470 auto deltaPhi = dphi[sel];
473 auto phiLocFin = phiLoc + deltaPhi;
475 if (math_utils::detail::abs<value_type>(phiLocFin) < MaxPhiLocSafe) {
476 auto deltaX = (math_utils::detail::sin<double>(phiLocFin) - track.getSnp()) / track.getCurvature(bz);
477 if (!track.propagateTo(track.getX() + deltaX, bz)) {
482 if (math_utils::detail::abs<value_type>(deltaPhi) < (2 * MaxPhiLocSafe)) {
483 auto rot = phiLoc + 0.5 * deltaPhi;
484 if (!track.rotate(track.getAlpha() + rot)) {
491 auto rot = phiLoc + (deltaPhi > 0 ? MaxPhiLocSafe : -MaxPhiLocSafe);
492 if (!track.rotate(track.getAlpha() + rot)) {
498 auto tgtPhiLoc = deltaPhi > 0 ? MaxPhiLocSafe : -MaxPhiLocSafe;
499 auto deltaX = (math_utils::detail::sin<double>(tgtPhiLoc) - track.getSnp()) / track.getCurvature(bz);
500 if (!track.propagateTo(track.getX() + deltaX, bz)) {
503 deltaPhi -= tgtPhiLoc - phiLoc;
504 phiLoc = deltaPhi > 0 ? MaxPhiLocSafe : -MaxPhiLocSafe;
511 if (!track.getXatLabR(
r, xfin, bz)) {
514 return propagateToX(track, xfin, bzOnly, maxSnp, maxStep, matCorr, tofInfo, signCorr);
518template <
typename value_T>
519template <
typename track_T>
520GPUd() bool
PropagatorImpl<value_T>::propagateToAlphaX(track_T& track, value_type
alpha, value_type
x,
bool bzOnly, value_type maxSnp, value_type maxStep,
int minSteps,
521 MatCorrType matCorr, track::TrackLTIntegral* tofInfo,
int signCorr)
const
524 auto snp = track.getSnpAt(
alpha,
x, getNominalBz());
526 if (math_utils::detail::abs<value_type>(snp) < maxSnp * 0.9 && track.rotate(
alpha)) {
527 auto dx = math_utils::detail::abs<value_type>(
x - track.getX());
531 return propagateTo(track,
x, bzOnly, maxSnp, math_utils::detail::min<value_type>(dx / minSteps, maxStep), matCorr, tofInfo, signCorr);
551template <
typename value_T>
552GPUd() bool
PropagatorImpl<value_T>::propagateToDCA(const
o2::dataformats::VertexBase& vtx, TrackParCov_t& track, value_type bZ,
554 o2::dataformats::DCA* dca, track::TrackLTIntegral* tofInfo,
555 int signCorr, value_type maxD)
const
558 value_type sn, cs, alp = track.getAlpha();
559 math_utils::detail::sincos<value_type>(alp, sn, cs);
560 value_type
x = track.getX(),
y = track.getY(), snp = track.getSnp(), csp = math_utils::detail::sqrt<value_type>((1.f - snp) * (1.f + snp));
561 value_type xv = vtx.getX() * cs + vtx.getY() * sn, yv = -vtx.getX() * sn + vtx.getY() * cs, zv = vtx.getZ();
565 value_type d = math_utils::detail::abs<value_type>(
x * snp -
y * csp);
573 value_type crv = track.getCurvature(bZ);
574 value_type tgfv = -(crv *
x - snp) / (crv *
y + csp);
575 sn = tgfv / math_utils::detail::sqrt<value_type>(1.f + tgfv * tgfv);
576 cs = math_utils::detail::sqrt<value_type>((1. - sn) * (1. + sn));
579 x = xv * cs + yv * sn;
580 yv = -xv * sn + yv * cs;
584 alp += math_utils::detail::asin<value_type>(sn);
585 if (!tmpT.rotate(alp) || !propagateToX(tmpT, xv, bZ, 0.85, maxStep, matCorr, tofInfo, signCorr)) {
586#ifndef GPUCA_ALIGPUCODE
587 LOG(
debug) <<
"failed to propagate to alpha=" << alp <<
" X=" << xv << vtx <<
" | Track is: " << tmpT.asString();
588#elif !defined(GPUCA_NO_FMT)
589 LOG(
debug) <<
"failed to propagate to alpha=" << alp <<
" X=" << xv << vtx;
599 math_utils::detail::sincos<value_type>(alp, sn, cs);
600 auto s2ylocvtx = vtx.getSigmaX2() * sn * sn + vtx.getSigmaY2() * cs * cs - 2. * vtx.getSigmaXY() * cs * sn;
601 dca->set(track.getY() - yv, track.getZ() - zv,
602 track.getSigmaY2() + s2ylocvtx, track.getSigmaZY(), track.getSigmaZ2() + vtx.getSigmaZ2());
608template <
typename value_T>
609GPUd() bool
PropagatorImpl<value_T>::propagateToDCABxByBz(const
o2::dataformats::VertexBase& vtx, TrackParCov_t& track,
611 o2::dataformats::DCA* dca, track::TrackLTIntegral* tofInfo,
612 int signCorr, value_type maxD)
const
615 value_type sn, cs, alp = track.getAlpha();
616 math_utils::detail::sincos<value_type>(alp, sn, cs);
617 value_type
x = track.getX(),
y = track.getY(), snp = track.getSnp(), csp = math_utils::detail::sqrt<value_type>((1.f - snp) * (1.f + snp));
618 value_type xv = vtx.getX() * cs + vtx.getY() * sn, yv = -vtx.getX() * sn + vtx.getY() * cs, zv = vtx.getZ();
622 value_type d = math_utils::detail::abs<value_type>(
x * snp -
y * csp);
630 value_type crv = track.getCurvature(mNominalBz);
631 value_type tgfv = -(crv *
x - snp) / (crv *
y + csp);
632 sn = tgfv / math_utils::detail::sqrt<value_type>(1.f + tgfv * tgfv);
633 cs = math_utils::detail::sqrt<value_type>((1. - sn) * (1. + sn));
636 x = xv * cs + yv * sn;
637 yv = -xv * sn + yv * cs;
641 alp += math_utils::detail::asin<value_type>(sn);
642 if (!tmpT.rotate(alp) || !PropagateToXBxByBz(tmpT, xv, 0.85, maxStep, matCorr, tofInfo, signCorr)) {
643#ifndef GPUCA_ALIGPUCODE
644 LOG(
debug) <<
"failed to propagate to alpha=" << alp <<
" X=" << xv << vtx <<
" | Track is: " << tmpT.asString();
645#elif !defined(GPUCA_NO_FMT)
646 LOG(
debug) <<
"failed to propagate to alpha=" << alp <<
" X=" << xv << vtx;
656 math_utils::detail::sincos<value_type>(alp, sn, cs);
657 auto s2ylocvtx = vtx.getSigmaX2() * sn * sn + vtx.getSigmaY2() * cs * cs - 2. * vtx.getSigmaXY() * cs * sn;
658 dca->set(track.getY() - yv, track.getZ() - zv,
659 track.getSigmaY2() + s2ylocvtx, track.getSigmaZY(), track.getSigmaZ2() + vtx.getSigmaZ2());
665template <
typename value_T>
666GPUd() bool
PropagatorImpl<value_T>::propagateToDCA(const math_utils::
Point3D<value_type>& vtx, TrackPar_t& track, value_type bZ,
668 std::
array<value_type, 2>* dca, track::TrackLTIntegral* tofInfo,
669 int signCorr, value_type maxD)
const
672 value_type sn, cs, alp = track.getAlpha();
673 math_utils::detail::sincos<value_type>(alp, sn, cs);
674 value_type
x = track.getX(),
y = track.getY(), snp = track.getSnp(), csp = math_utils::detail::sqrt<value_type>((1.f - snp) * (1.f + snp));
675 value_type xv = vtx.X() * cs + vtx.Y() * sn, yv = -vtx.X() * sn + vtx.Y() * cs, zv = vtx.Z();
679 value_type d = math_utils::detail::abs<value_type>(
x * snp -
y * csp);
687 value_type crv = track.getCurvature(bZ);
688 value_type tgfv = -(crv *
x - snp) / (crv *
y + csp);
689 sn = tgfv / math_utils::detail::sqrt<value_type>(1.f + tgfv * tgfv);
690 cs = math_utils::detail::sqrt<value_type>((1. - sn) * (1. + sn));
693 x = xv * cs + yv * sn;
694 yv = -xv * sn + yv * cs;
698 alp += math_utils::detail::asin<value_type>(sn);
699 if (!tmpT.rotateParam(alp) || !propagateToX(tmpT, xv, bZ, 0.85, maxStep, matCorr, tofInfo, signCorr)) {
700#ifndef GPUCA_ALIGPUCODE
701 LOG(
debug) <<
"failed to propagate to alpha=" << alp <<
" X=" << xv <<
" for vertex "
702 << vtx.X() <<
' ' << vtx.Y() <<
' ' << vtx.Z() <<
" | Track is: " << tmpT.asString();
704 LOG(
debug) <<
"failed to propagate to alpha=" << alp <<
" X=" << xv <<
" for vertex " << vtx.X() <<
' ' << vtx.Y() <<
' ' << vtx.Z();
714 (*dca)[0] = track.getY() - yv;
715 (*dca)[1] = track.getZ() - zv;
721template <
typename value_T>
724 std::
array<value_type, 2>* dca, track::TrackLTIntegral* tofInfo,
725 int signCorr, value_type maxD)
const
728 value_type sn, cs, alp = track.getAlpha();
729 math_utils::detail::sincos<value_type>(alp, sn, cs);
730 value_type
x = track.getX(),
y = track.getY(), snp = track.getSnp(), csp = math_utils::detail::sqrt<value_type>((1.f - snp) * (1.f + snp));
731 value_type xv = vtx.X() * cs + vtx.Y() * sn, yv = -vtx.X() * sn + vtx.Y() * cs, zv = vtx.Z();
735 value_type d = math_utils::detail::abs<value_type>(
x * snp -
y * csp);
743 value_type crv = track.getCurvature(mNominalBz);
744 value_type tgfv = -(crv *
x - snp) / (crv *
y + csp);
745 sn = tgfv / math_utils::detail::sqrt<value_type>(1.f + tgfv * tgfv);
746 cs = math_utils::detail::sqrt<value_type>((1. - sn) * (1. + sn));
749 x = xv * cs + yv * sn;
750 yv = -xv * sn + yv * cs;
754 alp += math_utils::detail::asin<value_type>(sn);
755 if (!tmpT.rotateParam(alp) || !PropagateToXBxByBz(tmpT, xv, 0.85, maxStep, matCorr, tofInfo, signCorr)) {
756#ifndef GPUCA_ALIGPUCODE
757 LOG(
debug) <<
"failed to propagate to alpha=" << alp <<
" X=" << xv <<
" for vertex "
758 << vtx.X() <<
' ' << vtx.Y() <<
' ' << vtx.Z() <<
" | Track is: " << tmpT.asString();
760 LOG(
debug) <<
"failed to propagate to alpha=" << alp <<
" X=" << xv <<
" for vertex " << vtx.X() <<
' ' << vtx.Y() <<
' ' << vtx.Z();
770 (*dca)[0] = track.getY() - yv;
771 (*dca)[1] = track.getZ() - zv;
777template <
typename value_T>
778GPUd() float
PropagatorImpl<value_T>::estimateLTIncrement(const
o2::track::TrackParametrization<value_type>& trc,
779 const
o2::math_utils::
Point3D<value_type>& posStart,
780 const
o2::math_utils::
Point3D<value_type>& posEnd)
const
783 float dX = posEnd.X() - posStart.X(), dY = posEnd.Y() - posStart.Y(), dZ = posEnd.Z() - posStart.Z(), d2XY = dX * dX + dY * dY;
784 if (getNominalBz() != 0) {
787 getFieldXYZ(posAv,
b);
788 float curvH = math_utils::detail::abs<value_type>(0.5f * trc.getCurvature(
b[2])), asArg = curvH * math_utils::detail::sqrt<value_type>(d2XY);
794 return math_utils::detail::sqrt<value_type>(d2XY + dZ * dZ);
798template <
typename value_T>
799GPUd() value_T
PropagatorImpl<value_T>::estimateLTFast(
o2::track::TrackLTIntegral& lt, const
o2::track::TrackParametrization<value_type>& trc)
const
801 value_T xdca = 0., ydca = 0.,
length = 0.;
803 constexpr float TinyF = 1e-9;
804 auto straigh_line_approx = [&]() {
805 auto csp2 = (1.f - trc.getSnp()) * (1.f + trc.getSnp());
807 auto csp = math_utils::detail::sqrt<value_type>(csp2);
808 auto tgp = trc.getSnp() / csp,
f = trc.getX() * tgp - trc.getY();
809 xdca = tgp *
f * csp2;
811 auto dx = xdca - trc.getX(), dy = ydca - trc.getY(), dz = dx * trc.getTgl() / csp;
812 return math_utils::detail::sqrt<value_type>(dx * dx + dy * dy + dz * dz);
815 return math_utils::detail::abs<value_type>(trc.getY() * math_utils::detail::sqrt<value_type>(1. + trc.getTgl() * trc.getTgl()));
818 trc.getCircleParamsLoc(mNominalBz,
c);
820 auto distC = math_utils::detail::sqrt<value_type>(
c.getCenterD2());
822 auto nrm = (distC -
c.rC) / distC;
825 auto v0x = trc.getX() -
c.xC, v0y = trc.getY() -
c.yC, v1x = xdca -
c.xC, v1y = ydca -
c.yC;
826 auto angcos = (v0x * v1x + v0y * v1y) / (
c.rC *
c.rC);
827 if (math_utils::detail::abs<value_type>(angcos) < 1.f) {
828 auto ang = math_utils::detail::acos<value_type>(angcos);
829 if ((trc.getSign() > 0.f) == (mNominalBz > 0.f)) {
834 length = math_utils::detail::abs<value_type>(
c.rC * ang * math_utils::detail::sqrt<value_type>(1. + trc.getTgl() * trc.getTgl()));
836 length = straigh_line_approx();
843 length = straigh_line_approx();
846 value_T dcaT = math_utils::detail::sqrt<value_type>(xdca * xdca + ydca * ydca);
848 lt.addStep(
length, trc.getQ2P2());
853template <
typename value_T>
856#if !defined(GPUCA_STANDALONE) && !defined(GPUCA_GPUCODE)
857 if (corrType == MatCorrType::USEMatCorrTGeo) {
861 if (mTGeoFallBackAllowed) {
864 throw std::runtime_error(
"requested MatLUT is absent and fall-back to TGeo is disabled");
868 return mMatLUT->getMatBudget(
p0.X(),
p0.Y(),
p0.Z(),
p1.X(),
p1.Y(),
p1.Z());
871template <
typename value_T>
876#if defined(GPUCA_GPUCODE_DEVICE) && defined(GPUCA_HAS_GLOBAL_SYMBOL_CONSTANT_MEM)
879 const auto*
f = mGPUField;
882 f->GetField(xyz.X(), xyz.Y(), xyz.Z(), bxyzF);
884 constexpr value_type kCLight1 = 1. / o2::gpu::gpu_common_constants::kCLight;
885 for (uint
i = 0;
i < 3; ++
i) {
886 bxyz[
i] =
static_cast<value_type
>(bxyzF[
i]) * kCLight1;
891 mFieldFast->Field(xyz, bxyz);
893#ifdef GPUCA_STANDALONE
894 LOG(fatal) <<
"Normal field cannot be used in standalone benchmark";
896 mField->field(xyz, bxyz);
903template <
typename value_T>
909#if defined(GPUCA_GPUCODE_DEVICE) && defined(GPUCA_HAS_GLOBAL_SYMBOL_CONSTANT_MEM)
912 const auto*
f = mGPUField;
914 constexpr value_type kCLight1 = 1. / o2::gpu::gpu_common_constants::kCLight;
915 bz =
f->GetFieldBz(xyz.X(), xyz.Y(), xyz.Z()) * kCLight1;
919 mFieldFast->GetBz(xyz, bz);
921#ifdef GPUCA_STANDALONE
922 LOG(fatal) <<
"Normal field cannot be used in standalone benchmark";
924 bz = mField->GetBz(xyz.X(), xyz.Y(), xyz.Z());
932template <
typename value_T>
935 getFieldXYZImpl<float>(xyz, bxyz);
938template <
typename value_T>
941 getFieldXYZImpl<double>(xyz, bxyz);
944template <
typename value_T>
947 return getBzImpl<float>(xyz);
950template <
typename value_T>
953 return getBzImpl<double>(xyz);
958#if !defined(GPUCA_GPUCODE) || defined(GPUCA_GPUCODE_DEVICE)