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);
569 value_type crv = track.getCurvature(bZ);
570 value_type tgfv = -(crv *
x - snp) / (crv *
y + csp);
571 sn = tgfv / math_utils::detail::sqrt<value_type>(1.f + tgfv * tgfv);
572 cs = math_utils::detail::sqrt<value_type>((1. - sn) * (1. + sn));
575 x = xv * cs + yv * sn;
576 yv = -xv * sn + yv * cs;
580 alp += math_utils::detail::asin<value_type>(sn);
581 if (!tmpT.rotate(alp) || !propagateToX(tmpT, xv, bZ, 0.85, maxStep, matCorr, tofInfo, signCorr)) {
582#ifndef GPUCA_ALIGPUCODE
583 LOG(
debug) <<
"failed to propagate to alpha=" << alp <<
" X=" << xv << vtx <<
" | Track is: " << tmpT.asString();
584#elif !defined(GPUCA_NO_FMT)
585 LOG(
debug) <<
"failed to propagate to alpha=" << alp <<
" X=" << xv << vtx;
591 math_utils::detail::sincos<value_type>(alp, sn, cs);
592 auto s2ylocvtx = vtx.getSigmaX2() * sn * sn + vtx.getSigmaY2() * cs * cs - 2. * vtx.getSigmaXY() * cs * sn;
593 dca->set(track.getY() - yv, track.getZ() - zv,
594 track.getSigmaY2() + s2ylocvtx, track.getSigmaZY(), track.getSigmaZ2() + vtx.getSigmaZ2());
600template <
typename value_T>
601GPUd() bool
PropagatorImpl<value_T>::propagateToDCABxByBz(const
o2::dataformats::VertexBase& vtx, TrackParCov_t& track,
603 o2::dataformats::DCA* dca, track::TrackLTIntegral* tofInfo,
604 int signCorr, value_type maxD)
const
607 value_type sn, cs, alp = track.getAlpha();
608 math_utils::detail::sincos<value_type>(alp, sn, cs);
609 value_type
x = track.getX(),
y = track.getY(), snp = track.getSnp(), csp = math_utils::detail::sqrt<value_type>((1.f - snp) * (1.f + snp));
610 value_type xv = vtx.getX() * cs + vtx.getY() * sn, yv = -vtx.getX() * sn + vtx.getY() * cs, zv = vtx.getZ();
614 value_type d = math_utils::detail::abs<value_type>(
x * snp -
y * csp);
618 value_type crv = track.getCurvature(mNominalBz);
619 value_type tgfv = -(crv *
x - snp) / (crv *
y + csp);
620 sn = tgfv / math_utils::detail::sqrt<value_type>(1.f + tgfv * tgfv);
621 cs = math_utils::detail::sqrt<value_type>((1. - sn) * (1. + sn));
624 x = xv * cs + yv * sn;
625 yv = -xv * sn + yv * cs;
629 alp += math_utils::detail::asin<value_type>(sn);
630 if (!tmpT.rotate(alp) || !PropagateToXBxByBz(tmpT, xv, 0.85, maxStep, matCorr, tofInfo, signCorr)) {
631#ifndef GPUCA_ALIGPUCODE
632 LOG(
debug) <<
"failed to propagate to alpha=" << alp <<
" X=" << xv << vtx <<
" | Track is: " << tmpT.asString();
633#elif !defined(GPUCA_NO_FMT)
634 LOG(
debug) <<
"failed to propagate to alpha=" << alp <<
" X=" << xv << vtx;
640 math_utils::detail::sincos<value_type>(alp, sn, cs);
641 auto s2ylocvtx = vtx.getSigmaX2() * sn * sn + vtx.getSigmaY2() * cs * cs - 2. * vtx.getSigmaXY() * cs * sn;
642 dca->set(track.getY() - yv, track.getZ() - zv,
643 track.getSigmaY2() + s2ylocvtx, track.getSigmaZY(), track.getSigmaZ2() + vtx.getSigmaZ2());
649template <
typename value_T>
650GPUd() bool
PropagatorImpl<value_T>::propagateToDCA(const math_utils::
Point3D<value_type>& vtx, TrackPar_t& track, value_type bZ,
652 std::
array<value_type, 2>* dca, track::TrackLTIntegral* tofInfo,
653 int signCorr, value_type maxD)
const
656 value_type sn, cs, alp = track.getAlpha();
657 math_utils::detail::sincos<value_type>(alp, sn, cs);
658 value_type
x = track.getX(),
y = track.getY(), snp = track.getSnp(), csp = math_utils::detail::sqrt<value_type>((1.f - snp) * (1.f + snp));
659 value_type xv = vtx.X() * cs + vtx.Y() * sn, yv = -vtx.X() * sn + vtx.Y() * cs, zv = vtx.Z();
663 value_type d = math_utils::detail::abs<value_type>(
x * snp -
y * csp);
667 value_type crv = track.getCurvature(bZ);
668 value_type tgfv = -(crv *
x - snp) / (crv *
y + csp);
669 sn = tgfv / math_utils::detail::sqrt<value_type>(1.f + tgfv * tgfv);
670 cs = math_utils::detail::sqrt<value_type>((1. - sn) * (1. + sn));
673 x = xv * cs + yv * sn;
674 yv = -xv * sn + yv * cs;
678 alp += math_utils::detail::asin<value_type>(sn);
679 if (!tmpT.rotateParam(alp) || !propagateToX(tmpT, xv, bZ, 0.85, maxStep, matCorr, tofInfo, signCorr)) {
680#ifndef GPUCA_ALIGPUCODE
681 LOG(
debug) <<
"failed to propagate to alpha=" << alp <<
" X=" << xv <<
" for vertex "
682 << vtx.X() <<
' ' << vtx.Y() <<
' ' << vtx.Z() <<
" | Track is: " << tmpT.asString();
684 LOG(
debug) <<
"failed to propagate to alpha=" << alp <<
" X=" << xv <<
" for vertex " << vtx.X() <<
' ' << vtx.Y() <<
' ' << vtx.Z();
690 (*dca)[0] = track.getY() - yv;
691 (*dca)[1] = track.getZ() - zv;
697template <
typename value_T>
700 std::
array<value_type, 2>* dca, track::TrackLTIntegral* tofInfo,
701 int signCorr, value_type maxD)
const
704 value_type sn, cs, alp = track.getAlpha();
705 math_utils::detail::sincos<value_type>(alp, sn, cs);
706 value_type
x = track.getX(),
y = track.getY(), snp = track.getSnp(), csp = math_utils::detail::sqrt<value_type>((1.f - snp) * (1.f + snp));
707 value_type xv = vtx.X() * cs + vtx.Y() * sn, yv = -vtx.X() * sn + vtx.Y() * cs, zv = vtx.Z();
711 value_type d = math_utils::detail::abs<value_type>(
x * snp -
y * csp);
715 value_type crv = track.getCurvature(mNominalBz);
716 value_type tgfv = -(crv *
x - snp) / (crv *
y + csp);
717 sn = tgfv / math_utils::detail::sqrt<value_type>(1.f + tgfv * tgfv);
718 cs = math_utils::detail::sqrt<value_type>((1. - sn) * (1. + sn));
721 x = xv * cs + yv * sn;
722 yv = -xv * sn + yv * cs;
726 alp += math_utils::detail::asin<value_type>(sn);
727 if (!tmpT.rotateParam(alp) || !PropagateToXBxByBz(tmpT, xv, 0.85, maxStep, matCorr, tofInfo, signCorr)) {
728#ifndef GPUCA_ALIGPUCODE
729 LOG(
debug) <<
"failed to propagate to alpha=" << alp <<
" X=" << xv <<
" for vertex "
730 << vtx.X() <<
' ' << vtx.Y() <<
' ' << vtx.Z() <<
" | Track is: " << tmpT.asString();
732 LOG(
debug) <<
"failed to propagate to alpha=" << alp <<
" X=" << xv <<
" for vertex " << vtx.X() <<
' ' << vtx.Y() <<
' ' << vtx.Z();
738 (*dca)[0] = track.getY() - yv;
739 (*dca)[1] = track.getZ() - zv;
745template <
typename value_T>
746GPUd() float
PropagatorImpl<value_T>::estimateLTIncrement(const
o2::track::TrackParametrization<value_type>& trc,
747 const
o2::math_utils::
Point3D<value_type>& posStart,
748 const
o2::math_utils::
Point3D<value_type>& posEnd)
const
751 float dX = posEnd.X() - posStart.X(), dY = posEnd.Y() - posStart.Y(), dZ = posEnd.Z() - posStart.Z(), d2XY = dX * dX + dY * dY;
752 if (getNominalBz() != 0) {
755 getFieldXYZ(posAv,
b);
756 float curvH = math_utils::detail::abs<value_type>(0.5f * trc.getCurvature(
b[2])), asArg = curvH * math_utils::detail::sqrt<value_type>(d2XY);
762 return math_utils::detail::sqrt<value_type>(d2XY + dZ * dZ);
766template <
typename value_T>
767GPUd() value_T
PropagatorImpl<value_T>::estimateLTFast(
o2::track::TrackLTIntegral& lt, const
o2::track::TrackParametrization<value_type>& trc)
const
769 value_T xdca = 0., ydca = 0.,
length = 0.;
771 constexpr float TinyF = 1e-9;
772 auto straigh_line_approx = [&]() {
773 auto csp2 = (1.f - trc.getSnp()) * (1.f + trc.getSnp());
775 auto csp = math_utils::detail::sqrt<value_type>(csp2);
776 auto tgp = trc.getSnp() / csp,
f = trc.getX() * tgp - trc.getY();
777 xdca = tgp *
f * csp2;
779 auto dx = xdca - trc.getX(), dy = ydca - trc.getY(), dz = dx * trc.getTgl() / csp;
780 return math_utils::detail::sqrt<value_type>(dx * dx + dy * dy + dz * dz);
783 return math_utils::detail::abs<value_type>(trc.getY() * math_utils::detail::sqrt<value_type>(1. + trc.getTgl() * trc.getTgl()));
786 trc.getCircleParamsLoc(mNominalBz,
c);
788 auto distC = math_utils::detail::sqrt<value_type>(
c.getCenterD2());
790 auto nrm = (distC -
c.rC) / distC;
793 auto v0x = trc.getX() -
c.xC, v0y = trc.getY() -
c.yC, v1x = xdca -
c.xC, v1y = ydca -
c.yC;
794 auto angcos = (v0x * v1x + v0y * v1y) / (
c.rC *
c.rC);
795 if (math_utils::detail::abs<value_type>(angcos) < 1.f) {
796 auto ang = math_utils::detail::acos<value_type>(angcos);
797 if ((trc.getSign() > 0.f) == (mNominalBz > 0.f)) {
802 length = math_utils::detail::abs<value_type>(
c.rC * ang * math_utils::detail::sqrt<value_type>(1. + trc.getTgl() * trc.getTgl()));
804 length = straigh_line_approx();
811 length = straigh_line_approx();
814 value_T dcaT = math_utils::detail::sqrt<value_type>(xdca * xdca + ydca * ydca);
816 lt.addStep(
length, trc.getQ2P2());
821template <
typename value_T>
824#if !defined(GPUCA_STANDALONE) && !defined(GPUCA_GPUCODE)
825 if (corrType == MatCorrType::USEMatCorrTGeo) {
829 if (mTGeoFallBackAllowed) {
832 throw std::runtime_error(
"requested MatLUT is absent and fall-back to TGeo is disabled");
836 return mMatLUT->getMatBudget(
p0.X(),
p0.Y(),
p0.Z(),
p1.X(),
p1.Y(),
p1.Z());
839template <
typename value_T>
844#if defined(GPUCA_GPUCODE_DEVICE) && defined(GPUCA_HAS_GLOBAL_SYMBOL_CONSTANT_MEM)
847 const auto*
f = mGPUField;
850 f->GetField(xyz.X(), xyz.Y(), xyz.Z(), bxyzF);
852 constexpr value_type kCLight1 = 1. / o2::gpu::gpu_common_constants::kCLight;
853 for (uint
i = 0;
i < 3; ++
i) {
854 bxyz[
i] =
static_cast<value_type
>(bxyzF[
i]) * kCLight1;
859 mFieldFast->Field(xyz, bxyz);
861#ifdef GPUCA_STANDALONE
862 LOG(fatal) <<
"Normal field cannot be used in standalone benchmark";
864 mField->field(xyz, bxyz);
871template <
typename value_T>
877#if defined(GPUCA_GPUCODE_DEVICE) && defined(GPUCA_HAS_GLOBAL_SYMBOL_CONSTANT_MEM)
880 const auto*
f = mGPUField;
882 constexpr value_type kCLight1 = 1. / o2::gpu::gpu_common_constants::kCLight;
883 bz =
f->GetFieldBz(xyz.X(), xyz.Y(), xyz.Z()) * kCLight1;
887 mFieldFast->GetBz(xyz, bz);
889#ifdef GPUCA_STANDALONE
890 LOG(fatal) <<
"Normal field cannot be used in standalone benchmark";
892 bz = mField->GetBz(xyz.X(), xyz.Y(), xyz.Z());
900template <
typename value_T>
903 getFieldXYZImpl<float>(xyz, bxyz);
906template <
typename value_T>
909 getFieldXYZImpl<double>(xyz, bxyz);
912template <
typename value_T>
915 return getBzImpl<float>(xyz);
918template <
typename value_T>
921 return getBzImpl<double>(xyz);
926#if !defined(GPUCA_GPUCODE) || defined(GPUCA_GPUCODE_DEVICE)