35 #include "CLHEP/Random/RandFlat.h"
36 #include "CLHEP/Random/RandGauss.h"
37 #include "GaudiKernel/PhysicalConstants.h"
57 if (cache.m_solenoid) {
60 cache.m_fieldCache.getField(
R,
H);
67 if (cache.m_solenoid) {
68 cache.m_fieldCache.getFieldZR(
R,
H, dH);
70 cache.m_fieldCache.getField(
R,
H, dH);
82 getMagneticField(Cache& cache,
87 constexpr
double magScale = 10000.;
88 const double*
R = position.data();
93 cache.m_includeBgradients) {
95 getFieldGradient(cache,
R,
H, dH);
96 BG[0] =
H[0] * magScale;
97 BG[1] =
H[1] * magScale;
98 BG[2] =
H[2] * magScale;
99 BG[3] = dH[0] * magScale;
100 BG[4] = dH[1] * magScale;
101 BG[5] = dH[2] * magScale;
102 BG[6] = dH[3] * magScale;
103 BG[7] = dH[4] * magScale;
104 BG[8] = dH[5] * magScale;
105 BG[9] = dH[6] * magScale;
106 BG[10] = dH[7] * magScale;
107 BG[11] = dH[8] * magScale;
111 getField(cache,
R,
H);
112 BG[0] =
H[0] * magScale;
113 BG[1] =
H[1] * magScale;
114 BG[2] =
H[2] * magScale;
116 for (
int i = 3;
i < 12; ++
i) {
126 std::unique_ptr<Trk::TrackParameters>
133 return std::make_unique<Trk::CurvilinearParameters>(
138 (inputTrackParameters->covariance() ? std::optional<
AmgSymMatrix(5)>(*inputTrackParameters->covariance())
147 (inputTrackParameters->covariance() ? std::optional<
AmgSymMatrix(5)>(*inputTrackParameters->covariance())
151 std::unique_ptr<Trk::TrackParameters>
153 std::vector<Trk::DestSurf>& targetSurfaces,
155 std::vector<unsigned int>& solutions,
164 std::vector<std::pair<unsigned int, double>> currentDist;
166 std::vector<std::pair<const Trk::Surface*, Trk::BoundaryCheck>>
::iterator sIter = targetSurfaces.begin();
167 std::vector<std::pair<unsigned int, double>>
::iterator oIter = currentDist.begin();
173 for (; sIter != targetSurfaces.end(); ++sIter) {
176 unsigned int numSf = sIter - sBeg;
177 if (distSol.
first() > tol) {
178 if (!currentDist.empty()) {
179 oIter = currentDist.end();
180 while (oIter != currentDist.begin() && distSol.
first() < (*(oIter - 1)).second)
182 oIter = currentDist.insert(oIter, std::pair<unsigned int, double>(numSf, distSol.
first()));
184 currentDist.emplace_back(numSf, distSol.
first());
188 if (!currentDist.empty()) {
189 oIter = currentDist.end();
190 while (oIter != currentDist.begin() && distSol.
second() < (*(oIter - 1)).second)
192 oIter = currentDist.insert(oIter, std::pair<unsigned int, double>(numSf, distSol.
second()));
194 currentDist.emplace_back(numSf, distSol.
second());
201 for (oIter = currentDist.begin(); oIter != currentDist.end(); ++oIter) {
202 Amg::Vector3D xsct = position + propDir * direction * ((*oIter).second);
203 if (targetSurfaces[(*oIter).first].first->isOnSurface(xsct, (*oIter).second, 0.001, 0.001)) {
204 if (usePathLimit &&
path > 0. &&
path < ((*oIter).second)) {
206 return std::make_unique<Trk::CurvilinearParameters>(endpoint, parm.
momentum(), parm.
charge());
208 path = (*oIter).second;
209 solutions.push_back((*oIter).first);
213 return std::make_unique<Trk::CurvilinearParameters>(xsct, parm.
momentum(), parm.
charge());
215 return sf->createUniqueTrackParameters(xsct, parm.
momentum(), parm.
charge(), std::nullopt);
235 dEds(Cache& cache,
double p)
237 cache.m_delIoni = 0.;
243 if (cache.m_material->x0() == 0 || cache.m_material->averageZ() == 0)
251 double eLoss = cache.m_MPV ? 0.9 * cache.m_delIoni + 0.15 * cache.m_delRad : cache.m_delIoni + cache.m_delRad;
261 dgdlambda(Cache& cache,
double l)
265 if (cache.m_material->x0() == 0 || cache.m_material->averageZ() == 0)
267 if (cache.m_material->zOverAtimesRho() == 0)
270 double p = std::abs(1. /
l);
273 double E = std::sqrt(
p *
p +
m *
m);
276 double I = 16.e-6 *
std::pow(cache.m_material->averageZ(), 0.9);
277 double kaz = 0.5 * 30.7075 * cache.m_material->zOverAtimesRho();
280 double lnCore = 4. * me * me / (
m *
m *
m *
m *
I *
I *
l *
l *
l *
l) / (1. + 2. *
gamma * me /
m + me * me /
m *
m);
281 double lnCore_deriv =
282 -4. * me * me / (
m *
m *
m *
m *
I *
I) *
285 4. *
l *
l *
l * me * me / (
m *
m));
286 double ln_deriv = 2. *
l *
m *
m *
std::log(lnCore) + lnCore_deriv / (lnCore *
beta *
beta);
287 double Bethe_Bloch_deriv = -kaz * ln_deriv;
291 double delta = 2. *
std::log(28.816
e-6 * std::sqrt(1000. * cache.m_material->zOverAtimesRho()) /
I) +
293 double delta_deriv = -2. / (
l *
beta *
beta) + 2. * delta *
l *
m *
m;
294 Bethe_Bloch_deriv += kaz * delta_deriv;
298 double Bethe_Heitler_deriv = me * me / (
m *
m * cache.m_material->x0() *
l *
l *
l *
E);
301 double radiative_deriv = 0.;
302 if ((cache.m_particle ==
Trk::muon) && (
E > 8000.)) {
304 radiative_deriv = 6.803e-5 / (cache.m_material->x0() *
l *
l *
l *
E) +
305 2. * 2.278
e-11 / (cache.m_material->x0() *
l *
l *
l) -
306 3. * 9.899
e-18 *
E / (cache.m_material->x0() *
l *
l *
l);
308 radiative_deriv = 9.253e-5 / (cache.m_material->x0() *
l *
l *
l *
E);
314 return 0.9 * Bethe_Bloch_deriv + 0.15 * Bethe_Heitler_deriv + 0.15 * radiative_deriv;
316 return Bethe_Bloch_deriv + Bethe_Heitler_deriv + radiative_deriv;
324 updateMaterialEffects(Cache& cache,
double mom,
double sinTheta,
double path)
329 if (cache.m_straggling) {
330 cache.m_covariance(4, 4) += cache.m_stragglingVariance;
331 cache.m_stragglingVariance = 0.;
334 if (!cache.m_multipleScattering)
337 double totalMomentumLoss =
mom - cache.m_matupd_lastmom;
338 double pathSinceLastUpdate =
path - cache.m_matupd_lastpath;
340 double pathAbs = std::abs(pathSinceLastUpdate);
342 if (pathAbs < 1.
e-03)
345 double matX0 = cache.m_material->x0();
348 int msLayers =
int(pathAbs / cache.m_layXmax / matX0) + 1;
351 double layerThickness = pathAbs / msLayers;
352 double radiationLengths = layerThickness / matX0;
354 double remainingPath = pathAbs;
359 double dLambdads = 0.;
360 double thetaVariance = 0.;
362 double average_dEds = totalMomentumLoss / pathAbs;
364 double cumulatedVariance = cache.m_inputThetaVariance + cache.m_combinedCovariance(3, 3) + cache.m_covariance(3, 3);
365 if (cumulatedVariance < 0) {
366 cumulatedVariance = 0;
368 double cumulatedX0 = 0.;
370 bool useCache = cache.m_extrapolationCache !=
nullptr;
372 double dX0 = std::abs(cache.m_combinedThickness) - pathAbs / matX0;
375 if (cache.m_extrapolationCache->x0tot() > 0)
376 cumulatedX0 = cache.m_extrapolationCache->x0tot() + dX0;
383 momentum = cache.m_matupd_lastmom + totalMomentumLoss * (
layer - 0.5) / msLayers;
386 E = std::sqrt(mom2 + massSquared);
389 double C0 = 13.6 * 13.6 / mom2 /
beta /
beta;
390 double C1 = 2 * 0.038 * C0;
391 double C2 = 0.038 * 0.038 * C0;
393 double MS2 = radiationLengths;
396 remainingPath -= layerThickness;
403 cumulatedX0 = cumulatedVariance / C0;
404 if (cumulatedX0 > 0.001) {
405 double lX0 =
log(cumulatedX0);
406 cumulatedX0 = cumulatedX0 / (1 + 2 * 0.038 * lX0 + 0.038 * 0.038 * lX0 * lX0);
410 double MS2s = MS2 + cumulatedX0;
411 thetaVariance = C0 * MS2 + C1 * MS2s *
log(MS2s) + C2 * MS2s *
log(MS2s) *
log(MS2s);
413 if (cumulatedX0 > 0.001) {
414 double lX0 =
log(cumulatedX0);
415 thetaVariance += -C1 * cumulatedX0 * lX0 - C2 * cumulatedX0 * lX0 * lX0;
419 double varScale = cache.m_scatteringScale * cache.m_scatteringScale;
420 double positionVariance = thetaVariance * (layerThickness * layerThickness / 3. + remainingPath * layerThickness +
421 remainingPath * remainingPath);
422 cache.m_covariance(0, 0) += varScale * positionVariance;
423 cache.m_covariance(1, 1) += varScale * positionVariance;
424 cache.m_covariance.fillSymmetric(
425 2, 0, cache.m_covariance(2, 0) + varScale * thetaVariance / (sinTheta) * (layerThickness / 2. + remainingPath));
426 cache.m_covariance(2, 2) += varScale * thetaVariance / (sinTheta * sinTheta);
427 cache.m_covariance.fillSymmetric(
428 3, 1, cache.m_covariance(3, 1) + varScale * thetaVariance * (-layerThickness / 2. - remainingPath));
429 cache.m_covariance(3, 3) += varScale * thetaVariance;
430 cache.m_covariance(4, 4) +=
431 varScale * 3. * thetaVariance * thetaVariance * layerThickness * layerThickness * dLambdads * dLambdads;
432 cumulatedVariance += varScale * thetaVariance;
436 cache.m_matupd_lastmom =
mom;
437 cache.m_matupd_lastpath =
path;
445 covarianceContribution(Cache& cache,
452 double finalMomentum = targetParms->
momentum().mag();
455 updateMaterialEffects(cache, finalMomentum,
sin(trackParameters->
momentum().theta()),
path);
462 Trk::
RungeKuttaUtils::newCovarianceMatrix(Jac, cache.m_combinedCovariance + cache.m_covariance);
464 measurementCovariance += localMSCov;
472 covarianceContribution(Cache& cache,
475 double finalMomentum,
479 updateMaterialEffects(cache, finalMomentum,
sin(trackParameters->momentum().theta()),
path);
482 measurementCovariance += cache.m_combinedCovariance + cache.m_covariance;
499 rungeKuttaStep(Cache& cache,
500 bool errorPropagation,
506 double& distanceStepped)
508 constexpr
double sol = 0.0299792458;
511 double lambda1 =
P[6];
512 double lambda2 =
P[6];
513 double lambda3 =
P[6];
514 double lambda4 =
P[6];
523 double initialMomentum = std::abs(1. /
P[6]);
533 double BG23[12] = { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. };
535 double BG4[12] = { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. };
544 if (cache.m_energyLoss && cache.m_matPropOK) {
548 if (errorPropagation) {
549 if (cache.m_includeGgradient) {
550 dgdl = dgdlambda(cache, lambda1);
553 -lambda1 * lambda1 *
g *
E * (3. - (
momentum *
momentum) / (
E *
E)) - lambda1 * lambda1 * lambda1 *
E * dgdl;
563 getMagneticField(cache,
570 dir.y() * BG1[2] -
dir.z() * BG1[1],
dir.z() * BG1[0] -
dir.x() * BG1[2],
dir.x() * BG1[1] -
dir.y() * BG1[0]);
571 k1 = sol * lambda1 * k1;
574 if (
steps++ > cache.m_maxSteps)
577 if (cache.m_energyLoss && cache.m_matPropOK) {
578 momentum = initialMomentum + (
h / 2.) * dP1;
579 if (
momentum <= cache.m_momentumCutOff)
584 if (errorPropagation) {
586 -lambda2 * lambda2 *
g *
E * (3. - (
momentum *
momentum) / (
E *
E)) - lambda2 * lambda2 * lambda2 *
E * dgdl;
589 pos = initialPos + (
h / 2.) * initialDir + (
h *
h / 8.) * k1;
590 dir = initialDir + (
h / 2.) * k1;
593 getMagneticField(cache,
pos, errorPropagation, BG23);
595 dir.z() * BG23[0] -
dir.x() * BG23[2],
596 dir.x() * BG23[1] -
dir.y() * BG23[0]);
597 k2 = sol * lambda2 * k2;
600 if (cache.m_energyLoss && cache.m_matPropOK) {
601 momentum = initialMomentum + (
h / 2.) * dP2;
602 if (
momentum <= cache.m_momentumCutOff)
607 if (errorPropagation) {
609 -lambda3 * lambda3 *
g *
E * (3. - (
momentum *
momentum) / (
E *
E)) - lambda3 * lambda3 * lambda3 *
E * dgdl;
612 dir = initialDir + (
h / 2.) * k2;
615 dir.z() * BG23[0] -
dir.x() * BG23[2],
616 dir.x() * BG23[1] -
dir.y() * BG23[0]);
617 k3 = sol * lambda3 * k3;
620 if (cache.m_energyLoss && cache.m_matPropOK) {
622 if (
momentum <= cache.m_momentumCutOff)
627 if (errorPropagation) {
629 -lambda4 * lambda4 *
g *
E * (3. - (
momentum *
momentum) / (
E *
E)) - lambda4 * lambda4 * lambda4 *
E * dgdl;
632 pos = initialPos +
h * initialDir + (
h *
h / 2.) * k3;
633 dir = initialDir +
h * k3;
637 getMagneticField(cache,
pos, errorPropagation, BG4);
639 dir.y() * BG4[2] -
dir.z() * BG4[1],
dir.z() * BG4[0] -
dir.x() * BG4[2],
dir.x() * BG4[1] -
dir.y() * BG4[0]);
640 k4 = sol * lambda4 * k4;
644 double errorEstimate =
std::max(errorPos.mag(), 1
e-20);
648 h =
h *
std::min(
std::max(0.25, std::sqrt(std::sqrt(cache.m_tolerance / errorEstimate))), 4.);
651 if (errorEstimate > 4. * cache.m_tolerance)
656 pos = initialPos + distanceStepped * initialDir + (distanceStepped * distanceStepped / 6.) * (k1 + k2 + k3);
662 dir = initialDir + (distanceStepped / 6.) * (k1 + 2. * k2 + 2. * k3 + k4);
668 double norm = 1. / std::sqrt(
P[3] *
P[3] +
P[4] *
P[4] +
P[5] *
P[5]);
674 if (cache.m_energyLoss && cache.m_matPropOK) {
675 momentum = initialMomentum + (distanceStepped / 6.) * (dP1 + 2. * dP2 + 2. * dP3 + dP4);
676 if (
momentum <= cache.m_momentumCutOff)
689 if (errorPropagation) {
690 double initialjLambda =
P[41];
691 double jLambda = initialjLambda;
697 for (
int i = 7;
i < 42;
i += 7) {
710 jDir.z() * BG1[0] - jDir.x() * BG1[2],
711 jDir.x() * BG1[1] - jDir.y() * BG1[0]);
714 if (cache.m_includeBgradients) {
716 (dir1.y() * BG1[9] - dir1.z() * BG1[6]) * jPos.x() + (dir1.y() * BG1[10] - dir1.z() * BG1[7]) * jPos.y() +
717 (dir1.y() * BG1[11] - dir1.z() * BG1[8]) * jPos.z(),
718 (dir1.z() * BG1[3] - dir1.x() * BG1[9]) * jPos.x() + (dir1.z() * BG1[4] - dir1.x() * BG1[10]) * jPos.y() +
719 (dir1.z() * BG1[5] - dir1.x() * BG1[11]) * jPos.z(),
720 (dir1.x() * BG1[6] - dir1.y() * BG1[3]) * jPos.x() + (dir1.x() * BG1[7] - dir1.y() * BG1[4]) * jPos.y() +
721 (dir1.x() * BG1[8] - dir1.y() * BG1[5]) * jPos.z());
724 jk1 = sol * lambda1 * jk1;
728 jdL1 = dL1 * jLambda;
729 jk1 = jk1 + k1 * jLambda / lambda1;
733 jPos = initialjPos + (distanceStepped / 2.) * initialjDir + (distanceStepped * distanceStepped / 8.) * jk1;
734 jDir = initialjDir + (distanceStepped / 2.) * jk1;
737 jDir.z() * BG23[0] - jDir.x() * BG23[2],
738 jDir.x() * BG23[1] - jDir.y() * BG23[0]);
740 if (cache.m_includeBgradients) {
741 Amg::Vector3D C((dir2.y() * BG23[9] - dir2.z() * BG23[6]) * jPos.x() +
742 (dir2.y() * BG23[10] - dir2.z() * BG23[7]) * jPos.y() +
743 (dir2.y() * BG23[11] - dir2.z() * BG23[8]) * jPos.z(),
744 (dir2.z() * BG23[3] - dir2.x() * BG23[9]) * jPos.x() +
745 (dir2.z() * BG23[4] - dir2.x() * BG23[10]) * jPos.y() +
746 (dir2.z() * BG23[5] - dir2.x() * BG23[11]) * jPos.z(),
747 (dir2.x() * BG23[6] - dir2.y() * BG23[3]) * jPos.x() +
748 (dir2.x() * BG23[7] - dir2.y() * BG23[4]) * jPos.y() +
749 (dir2.x() * BG23[8] - dir2.y() * BG23[5]) * jPos.z());
752 jk2 = sol * lambda2 * jk2;
755 jLambda = initialjLambda + (distanceStepped / 2.) * jdL1;
756 jdL2 = dL2 * jLambda;
757 jk2 = jk2 + k2 * jLambda / lambda2;
761 jDir = initialjDir + (distanceStepped / 2.) * jk2;
764 jDir.z() * BG23[0] - jDir.x() * BG23[2],
765 jDir.x() * BG23[1] - jDir.y() * BG23[0]);
767 if (cache.m_includeBgradients) {
768 Amg::Vector3D C((dir3.y() * BG23[9] - dir3.z() * BG23[6]) * jPos.x() +
769 (dir3.y() * BG23[10] - dir3.z() * BG23[7]) * jPos.y() +
770 (dir3.y() * BG23[11] - dir3.z() * BG23[8]) * jPos.z(),
771 (dir3.z() * BG23[3] - dir3.x() * BG23[9]) * jPos.x() +
772 (dir3.z() * BG23[4] - dir3.x() * BG23[10]) * jPos.y() +
773 (dir3.z() * BG23[5] - dir3.x() * BG23[11]) * jPos.z(),
774 (dir3.x() * BG23[6] - dir3.y() * BG23[3]) * jPos.x() +
775 (dir3.x() * BG23[7] - dir3.y() * BG23[4]) * jPos.y() +
776 (dir3.x() * BG23[8] - dir3.y() * BG23[5]) * jPos.z());
779 jk3 = sol * lambda3 * jk3;
782 jLambda = initialjLambda + (distanceStepped / 2.) * jdL2;
783 jdL3 = dL3 * jLambda;
784 jk3 = jk3 + k3 * jLambda / lambda3;
788 jPos = initialjPos + distanceStepped * initialjDir + (distanceStepped * distanceStepped / 2.) * jk3;
789 jDir = initialjDir + distanceStepped * jk3;
792 jDir.z() * BG4[0] - jDir.x() * BG4[2],
793 jDir.x() * BG4[1] - jDir.y() * BG4[0]);
795 if (cache.m_includeBgradients) {
797 (dir4.y() * BG4[9] - dir4.z() * BG4[6]) * jPos.x() + (dir4.y() * BG4[10] - dir4.z() * BG4[7]) * jPos.y() +
798 (dir4.y() * BG4[11] - dir4.z() * BG4[8]) * jPos.z(),
799 (dir4.z() * BG4[3] - dir4.x() * BG4[9]) * jPos.x() + (dir4.z() * BG4[4] - dir4.x() * BG4[10]) * jPos.y() +
800 (dir4.z() * BG4[5] - dir4.x() * BG4[11]) * jPos.z(),
801 (dir4.x() * BG4[6] - dir4.y() * BG4[3]) * jPos.x() + (dir4.x() * BG4[7] - dir4.y() * BG4[4]) * jPos.y() +
802 (dir4.x() * BG4[8] - dir4.y() * BG4[5]) * jPos.z());
805 jk4 = sol * lambda4 * jk4;
808 jLambda = initialjLambda + distanceStepped * jdL3;
809 jdL4 = dL4 * jLambda;
810 jk4 = jk4 + k4 * jLambda / lambda4;
815 initialjPos + distanceStepped * initialjDir + (distanceStepped * distanceStepped / 6.) * (jk1 + jk2 + jk3);
816 jDir = initialjDir + (distanceStepped / 6.) * (jk1 + 2. * jk2 + 2. * jk3 + jk4);
818 jLambda = initialjLambda + (distanceStepped / 6.) * (jdL1 + 2. * jdL2 + 2. * jdL3 + jdL4);
835 for (
int i = 0;
i < 12; ++
i) {
845 propagateWithJacobianImpl(Cache& cache,
846 bool errorPropagation,
852 double maxPath = cache.m_maxPath;
855 double dDir[3] = { 0., 0., 0. };
856 int targetPassed = 0;
857 double previousDistance = 0.;
858 double distanceStepped = 0.;
859 bool distanceEstimationSuccessful =
false;
860 double BG1[12] = { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. };
862 bool firstStep =
true;
863 double absolutePath = 0.;
870 surfaceType, targetSurface,
P, distanceEstimationSuccessful);
871 if (!distanceEstimationSuccessful) {
877 distanceToTarget > 100. ?
h = 100. : distanceToTarget < -100. ?
h = -100. :
h = distanceToTarget;
883 double distanceTolerance =
std::min(
std::max(std::abs(distanceToTarget) * cache.m_tolerance, 1
e-6), 1
e-2);
885 while (std::abs(distanceToTarget) > distanceTolerance) {
887 if (!rungeKuttaStep(cache, errorPropagation,
h,
P,
dDir, BG1, firstStep, distanceStepped))
889 path += distanceStepped;
890 absolutePath += std::abs(distanceStepped);
892 if (std::abs(distanceStepped) > 0.001) {
893 cache.m_sigmaIoni = cache.m_sigmaIoni - cache.m_kazL *
log(std::abs(distanceStepped));
896 if (errorPropagation && cache.m_straggling) {
897 double sigTot2 = cache.m_sigmaIoni * cache.m_sigmaIoni + cache.m_sigmaRad * cache.m_sigmaRad;
899 mom = std::abs(1. /
P[6]);
900 beta =
mom / std::sqrt(
mom *
mom + cache.m_particleMass * cache.m_particleMass);
902 cache.m_stragglingVariance += sigTot2 / (bp2 * bp2) * distanceStepped * distanceStepped;
904 if (cache.m_matstates || errorPropagation)
905 cache.m_combinedEloss.update(cache.m_delIoni * distanceStepped,
906 cache.m_sigmaIoni * std::abs(distanceStepped),
907 cache.m_delRad * distanceStepped,
908 cache.m_sigmaRad * std::abs(distanceStepped),
911 previousDistance = std::abs(distanceToTarget);
913 surfaceType, targetSurface,
P, distanceEstimationSuccessful);
914 if (!distanceEstimationSuccessful)
918 if (
h * distanceToTarget < 0.) {
923 if (std::abs(
h) > std::abs(distanceToTarget))
924 h = distanceToTarget;
927 if ((targetPassed > 3 && std::abs(distanceToTarget) >= previousDistance) || (absolutePath > maxPath))
930 if (
steps++ > cache.m_maxSteps)
934 if (cache.m_material && cache.m_material->x0() != 0.)
935 cache.m_combinedThickness += std::abs(
path) / cache.m_material->x0();
941 pos[0] =
pos[0] + distanceToTarget * (
dir[0] + 0.5 * distanceToTarget *
dDir[0]);
942 pos[1] =
pos[1] + distanceToTarget * (
dir[1] + 0.5 * distanceToTarget *
dDir[1]);
943 pos[2] =
pos[2] + distanceToTarget * (
dir[2] + 0.5 * distanceToTarget *
dDir[2]);
964 std::unique_ptr<Trk::TrackParameters>
965 propagateRungeKuttaImpl(Cache& cache,
966 bool errorPropagation,
978 cache.m_charge = inputTrackParameters.
charge();
980 std::unique_ptr<Trk::TrackParameters> trackParameters{};
986 if (cache.m_tolerance <= 0.)
988 if (cache.m_momentumCutOff < 0.)
992 if (inputTrackParameters.parameters()[
Trk::qOverP] == 0) {
993 trackParameters = createStraightLine(&inputTrackParameters);
994 if (!trackParameters) {
998 trackParameters.reset(inputTrackParameters.
clone());
1001 if (std::abs(1. / trackParameters->parameters()[
Trk::qOverP]) <= cache.m_momentumCutOff) {
1006 if (cache.m_matPropOK && ((cache.m_material->zOverAtimesRho() == 0.) || (cache.m_material->x0() == 0.) ||
1007 (cache.m_material->zOverAtimesRho() != cache.m_material->zOverAtimesRho()))) {
1008 cache.m_matPropOK =
false;
1010 if (cache.m_matPropOK && cache.m_straggling)
1011 cache.m_stragglingVariance = 0.;
1012 if (errorPropagation || cache.m_matstates) {
1013 cache.m_combinedCovariance.setZero();
1014 cache.m_covariance.setZero();
1016 cache.m_inputThetaVariance = trackParameters->covariance() ? (*trackParameters->covariance())(3, 3) : 0.;
1017 cache.m_combinedEloss.set(0., 0., 0., 0., 0., 0.);
1018 cache.m_combinedThickness = 0.;
1031 double d =
T(0, 3) *
T(0, 2) +
T(1, 3) *
T(1, 2) +
T(2, 3) *
T(2, 2);
1044 if (!propagateWithJacobianImpl(cache, errorPropagation, ty,
s, cache.m_P,
path)) {
1051 double s[6] = {
T(0, 3),
T(1, 3),
T(2, 3),
T(0, 2),
T(1, 2),
T(2, 2) };
1052 if (!propagateWithJacobianImpl(cache, errorPropagation, ty,
s, cache.m_P,
path)) {
1060 double s[9] = {
T(0, 3),
T(1, 3),
T(2, 3),
T(0, 2),
1063 if (!propagateWithJacobianImpl(cache, errorPropagation, ty,
s, cache.m_P,
path)) {
1070 double k =
static_cast<const Trk::ConeSurface*
>(&targetSurface)->bounds().tanAlpha();
1072 double s[9] = {
T(0, 3),
T(1, 3),
T(2, 3),
T(0, 2),
T(1, 2),
T(2, 2),
k, (
double)propagationDirection, 0. };
1073 if (!propagateWithJacobianImpl(cache, errorPropagation, ty,
s, cache.m_P,
path)) {
1080 double s[6] = {
T(0, 3),
T(1, 3),
T(2, 3), 0., 0., 1. };
1081 if (!propagateWithJacobianImpl(cache, errorPropagation, ty,
s, cache.m_P,
path)) {
1089 double d =
T(0, 3) *
T(0, 2) +
T(1, 3) *
T(1, 2) +
T(2, 3) *
T(2, 2);
1102 if (!propagateWithJacobianImpl(cache, errorPropagation, ty,
s, cache.m_P,
path)) {
1107 if (propagationDirection *
path < 0.) {
1118 if (boundaryCheck && !targetSurface.
isOnSurface(gp)) {
1122 if (!errorPropagation || !trackParameters->covariance()) {
1123 return std::make_unique<Trk::CurvilinearParameters>(gp, localp[2], localp[3], localp[4]);
1131 if (cache.m_matPropOK && (cache.m_multipleScattering || cache.m_straggling) && std::abs(
path) > 0.)
1132 covarianceContribution(cache, trackParameters.
get(),
path, std::abs(1. / cache.m_P[6]), measurementCovariance);
1135 gp, localp[2], localp[3], localp[4], std::move(measurementCovariance));
1141 if (boundaryCheck) {
1143 if (!targetSurface.insideBounds(localPosition, 0.)) {
1149 targetSurface.createUniqueTrackParameters(localp[0], localp[1], localp[2], localp[3], localp[4], std::nullopt);
1151 if (!errorPropagation || !trackParameters->covariance()) {
1152 return onTargetSurf;
1160 if (cache.m_matPropOK && (cache.m_multipleScattering || cache.m_straggling) && std::abs(
path) > 0.)
1161 covarianceContribution(cache, trackParameters.
get(),
path, onTargetSurf.
get(), measurementCovariance);
1163 return targetSurface.createUniqueTrackParameters(
1164 localp[0], localp[1], localp[2], localp[3], localp[4], std::move(measurementCovariance));
1176 declareInterface<Trk::IPropagator>(
this);
1201 ATH_MSG_WARNING(
"Simulation mode requested but material updator not found - no brem photon emission.");
1210 return StatusCode::SUCCESS;
1215 return StatusCode::SUCCESS;
1224 ATH_MSG_WARNING(
"[STEP_Propagator] STEP_Propagator does not handle neutral track parameters."
1225 <<
"Use the StraightLinePropagator instead.");
1240 double Jacobian[25];
1244 getFieldCacheObject(cache, ctx);
1245 setCacheFromProperties(cache);
1246 clearMaterialEffects(cache);
1263 return propagateRungeKuttaImpl(cache,
true, trackParameters, targetSurface, propagationDirection,
1264 magneticFieldProperties,
particle, boundaryCheck, Jacobian, returnCurv);
1275 std::vector<unsigned int>& solutions,
double&
path,
bool usePathLimit,
bool returnCurv,
1281 getFieldCacheObject(cache, ctx);
1282 setCacheFromProperties(cache);
1283 clearMaterialEffects(cache);
1311 return propagateNeutral(trackParameters, targetSurfaces, propagationDirection, solutions,
path,
1312 usePathLimit, returnCurv);
1314 return propagateRungeKutta(cache,
true, trackParameters, targetSurfaces, propagationDirection,
1315 magneticFieldProperties,
particle, solutions,
path, returnCurv);
1326 std::vector<unsigned int>& solutions,
PathLimit& pathLim,
TimeLimit& timeLim,
bool returnCurv,
1332 getFieldCacheObject(cache, ctx);
1333 setCacheFromProperties(cache);
1334 clearMaterialEffects(cache);
1362 double dTim = timeLim.
tMax - timeLim.
time;
1370 if (timMax > 0. && timMax <
path)
1372 bool usePathLimit = (
path > 0.);
1385 std::unique_ptr<Trk::TrackParameters> nextPar{};
1388 nextPar = propagateNeutral(trackParameters, targetSurfaces, propagationDirection, solutions,
path,
1389 usePathLimit, returnCurv);
1391 nextPar = propagateRungeKutta(cache,
true, trackParameters, targetSurfaces, propagationDirection,
1392 magneticFieldProperties,
particle, solutions,
path, returnCurv);
1412 std::vector<unsigned int>& solutions, std::vector<const Trk::TrackStateOnSurface*>* matstates,
1413 std::vector<std::pair<std::unique_ptr<Trk::TrackParameters>,
int>>* intersections,
double&
path,
1420 getFieldCacheObject(cache, ctx);
1421 setCacheFromProperties(cache);
1422 clearMaterialEffects(cache);
1454 return propagateNeutral(trackParameters, targetSurfaces, propagationDirection, solutions,
path,
1455 usePathLimit, returnCurv);
1457 return propagateRungeKutta(cache,
true, trackParameters, targetSurfaces, propagationDirection,
1458 magneticFieldProperties,
particle, solutions,
path, returnCurv);
1465 const EventContext& ctx,
1471 std::optional<Trk::TransportJacobian>& jacobian,
1477 double Jacobian[25];
1481 getFieldCacheObject(cache, ctx);
1482 setCacheFromProperties(cache);
1483 clearMaterialEffects(cache);
1500 std::unique_ptr<Trk::TrackParameters>
parameters =
1501 propagateRungeKuttaImpl(cache,
true, trackParameters, targetSurface, propagationDirection,
1502 magneticFieldProperties,
particle, boundaryCheck, Jacobian, returnCurv);
1505 Jacobian[24] = Jacobian[20];
1510 jacobian = std::make_optional<Trk::TransportJacobian>(Jacobian);
1528 double Jacobian[25];
1532 getFieldCacheObject(cache, ctx);
1533 setCacheFromProperties(cache);
1534 clearMaterialEffects(cache);
1546 return propagateRungeKuttaImpl(cache,
false, trackParameters, targetSurface, propagationDirection,
1547 magneticFieldProperties,
particle, boundaryCheck, Jacobian, returnCurv);
1554 const EventContext& ctx,
1560 std::optional<Trk::TransportJacobian>& jacobian,
1565 double Jacobian[25];
1570 getFieldCacheObject(cache, ctx);
1571 setCacheFromProperties(cache);
1572 clearMaterialEffects(cache);
1585 std::unique_ptr<Trk::TrackParameters>
parameters =
1586 propagateRungeKuttaImpl(cache,
true, trackParameters, targetSurface, propagationDirection,
1587 magneticFieldProperties,
particle, boundaryCheck, Jacobian, returnCurv);
1590 Jacobian[24] = Jacobian[20];
1595 jacobian = std::make_optional<Trk::TransportJacobian>(Jacobian);
1607 const EventContext& ctx,
1617 getFieldCacheObject(cache, ctx);
1618 setCacheFromProperties(cache);
1619 clearMaterialEffects(cache);
1639 return std::nullopt;
1642 return std::nullopt;
1645 return std::nullopt;
1656 return std::nullopt;
1665 double d = T(0, 3) * T(0, 2) + T(1, 3) * T(1, 2) + T(2, 3) * T(2, 2);
1678 if (!propagateWithJacobianImpl(cache,
false, ty,
s, cache.
m_P,
path)){
1679 return std::nullopt;
1685 double s[6] = {T(0, 3), T(1, 3), T(2, 3), T(0, 2), T(1, 2), T(2, 2)};
1686 if (!propagateWithJacobianImpl(cache,
false, ty,
s, cache.
m_P,
path)) {
1687 return std::nullopt;
1695 T(0, 3), T(1, 3), T(2, 3), T(0, 2), T(1, 2), T(2, 2), cyl->
bounds().
r(),
Trk::alongMomentum, 0.};
1696 if (!propagateWithJacobianImpl(cache,
false, ty,
s, cache.
m_P,
path)) {
1697 return std::nullopt;
1703 double k =
static_cast<const Trk::ConeSurface*
>(&targetSurface)->bounds().tanAlpha();
1705 double s[9] = {T(0, 3), T(1, 3), T(2, 3), T(0, 2), T(1, 2), T(2, 2),
k,
Trk::alongMomentum, 0.};
1706 if (!propagateWithJacobianImpl(cache,
false, ty,
s, cache.
m_P,
path)) {
1707 return std::nullopt;
1713 double s[6] = {T(0, 3), T(1, 3), T(2, 3), 0., 0., 1.};
1714 if (!propagateWithJacobianImpl(cache,
false, ty,
s, cache.
m_P,
path)) {
1715 return std::nullopt;
1722 double d = T(0, 3) * T(0, 2) + T(1, 3) * T(1, 2) + T(2, 3) * T(2, 2);
1735 if (!propagateWithJacobianImpl(cache,
false, ty,
s, cache.
m_P,
path)) {
1736 return std::nullopt;
1742 return std::make_optional<Trk::TrackSurfaceIntersection>(globalPosition, direction,
path);
1745 std::optional<Trk::TrackSurfaceIntersection>
1758 auto tmpTrackParameters =
1759 Trk::Perigee(0., 0., direction.phi(), direction.theta(),
qOverP, perigeeSurface, std::nullopt);
1761 std::optional<Trk::TrackSurfaceIntersection> solution =
1763 ?
intersect(ctx, tmpTrackParameters, surface,
1783 getFieldCacheObject(cache, ctx);
1784 setCacheFromProperties(cache);
1785 clearMaterialEffects(cache);
1811 double dDir[3] = {0., 0., 0.};
1812 double distanceStepped = 0.;
1813 double BG1[12] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
1815 bool firstStep =
true;
1817 double radius2Max = cylinderBounds.
r() * cylinderBounds.
r();
1819 double radius2 = PP[0] * PP[0] + PP[1] * PP[1];
1820 double direction = PP[0] * PP[3] + PP[1] * PP[4];
1821 double h = maxStepSize;
1824 if ((std::abs(PP[2]) > zMax) || (radius2 > radius2Max))
1829 positionsList.push_back(initialPosition);
1831 bool perigee =
false;
1832 if (std::abs(direction) < 0.00001) {
1836 for (
int i = 0;
i != 2; ++
i) {
1842 double p[7] = {PP[0], PP[1], PP[2], PP[3], PP[4], PP[5], PP[6]};
1844 while (std::abs(
path) < maxPath) {
1846 if (!rungeKuttaStep(cache,
false,
h,
p,
dDir, BG1, firstStep, distanceStepped))
1851 if (
h > maxStepSize) {
1853 }
else if (
h < -maxStepSize) {
1860 positionsList.push_back(globalPosition);
1862 positionsList.push_front(globalPosition);
1866 radius2 =
p[0] *
p[0] +
p[1] *
p[1];
1867 if ((std::abs(
p[2]) > zMax) || (radius2 > radius2Max))
1871 if ((
p[0] *
p[3] +
p[1] *
p[4]) * direction < 0.) {
1888 std::vector<unsigned int>& solutions,
double& totalPath,
bool returnCurv)
const {
1894 std::unique_ptr<Trk::TrackParameters> trackParameters{};
1906 if (inputTrackParameters.parameters()[
Trk::qOverP] == 0) {
1907 trackParameters = createStraightLine(&inputTrackParameters);
1908 if (!trackParameters) {
1912 trackParameters.reset(inputTrackParameters.
clone());
1925 if (errorPropagation && !trackParameters->covariance()) {
1926 errorPropagation =
false;
1931 cache.m_combinedCovariance.setZero();
1932 cache.m_covariance.setZero();
1936 cache.
m_inputThetaVariance = trackParameters->covariance() ? (*trackParameters->covariance())(3, 3) : 0.;
1952 bool validStep =
true;
1957 double Jacobian[21];
1960 validStep = propagateWithJacobian(cache, errorPropagation, targetSurfaces, cache.
m_P,
1961 propagationDirection, solutions,
path, totalPath);
1965 if (propagationDirection *
path <= 0.) {
1973 std::unique_ptr<Trk::CurvilinearParameters> cPar =
nullptr;
1975 if (!errorPropagation) {
1976 cPar = std::make_unique<Trk::CurvilinearParameters>(
1982 Trk::RungeKuttaUtils::newCovarianceMatrix(Jacobian, *trackParameters->covariance());
1985 std::abs(totalPath) > 0.) {
1986 covarianceContribution(cache, trackParameters.get(), totalPath, std::abs(1. / cache.
m_P[6]),
1987 measurementCovariance);
1989 cPar = std::make_unique<Trk::CurvilinearParameters>(
1991 std::move(measurementCovariance));
1997 dumpMaterialEffects(cache, cPar.get(), totalPath);
2006 bool solution =
false;
2007 std::vector<unsigned int> valid_solutions;
2008 valid_solutions.reserve(solutions.size());
2011 while (iSol != solutions.end()) {
2012 if (targetSurfaces[*iSol].
first->isOnSurface(gp, targetSurfaces[*iSol].second, 0.001, 0.001)) {
2019 cache.
m_P, localp, Jacobian);
2022 valid_solutions.push_back(*iSol);
2026 solutions = std::move(valid_solutions);
2031 if (solutions.empty()) {
2038 smear(cache, localp[2], localp[3], trackParameters.get(), radDist);
2041 std::unique_ptr<Trk::TrackParameters> onTargetSurf =
2044 : targetSurfaces[solutions[0]].
first->createUniqueTrackParameters(
2045 localp[0], localp[1], localp[2], localp[3], localp[4], std::nullopt);
2047 if (!errorPropagation) {
2050 return std::make_unique<Trk::CurvilinearParameters>(gp, localp[2], localp[3], localp[4]);
2052 return onTargetSurf;
2057 for (
double i : Jacobian) {
2064 Trk::RungeKuttaUtils::newCovarianceMatrix(Jacobian, *trackParameters->covariance());
2071 covarianceContribution(cache, trackParameters.get(), totalPath, std::abs(1. / cache.
m_P[6]),
2072 measurementCovariance);
2074 covarianceContribution(cache, trackParameters.get(), totalPath, onTargetSurf.get(),
2075 measurementCovariance);
2081 return std::make_unique<Trk::CurvilinearParameters>(gp, localp[2], localp[3], localp[4],
2082 std::move(measurementCovariance));
2086 return targetSurfaces[solutions[0]].first->createUniqueTrackParameters(
2087 localp[0], localp[1], localp[2], localp[3], localp[4], std::move(measurementCovariance));
2095 std::vector<DestSurf>& sfs,
double*
P,
2097 std::vector<unsigned int>& solutions,
double&
path,
2098 double sumPath)
const {
2100 double*
pos = &
P[0];
2101 double*
dir = &
P[3];
2102 double dDir[3] = {0., 0., 0.};
2104 double previousDistance = 0.;
2105 double distanceStepped = 0.;
2106 double BG1[12] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
2108 bool firstStep =
true;
2116 double helpSoft = 1.;
2119 int restartLimit = 10;
2139 double distanceToTarget = propDir * maxPath;
2142 int nextSf = sfs.size();
2143 int nextSfCand = nextSf;
2146 unsigned int numSf = 0;
2147 unsigned int iCurr = 0;
2149 for (; sIter != sfs.end(); ++sIter) {
2151 double distEst = -propDir * maxPath;
2152 double dist1Est = distEst;
2154 distEst = distSol.
first();
2155 dist1Est = distSol.
first();
2157 (std::abs(distEst) < tol || (propDir * distEst < -tol && propDir * distSol.
second() > tol)))
2158 distEst = distSol.
second();
2163 if (distEst * propDir > -tol) {
2165 cache.
m_currentDist[iCurr] = std::pair<int, std::pair<double, double>>(
2166 0, std::pair<double, double>(distEst, distSol.
currentDistance(
true)));
2168 cache.
m_currentDist[iCurr] = std::pair<int, std::pair<double, double>>(
2169 1, std::pair<double, double>(distEst, distSol.
currentDistance(
true)));
2171 if (tol < propDir * distEst && propDir * distEst < propDir * distanceToTarget) {
2172 distanceToTarget = distEst;
2178 cache.
m_currentDist[iCurr] = std::pair<int, std::pair<double, double>>(
2181 if (std::abs(dist1Est) < tol)
2182 startSf = (
int)iCurr;
2186 if (distanceToTarget == maxPath || numSf == 0)
2196 double absPath = 0.;
2197 distanceToTarget > 100. ?
h = 100. : distanceToTarget < -100. ?
h = -100. :
h = distanceToTarget;
2206 std::pair<size_t, float> dist2next = lbu->
distanceToNext(position, propDir * direction0);
2207 if (dist2next.first < lbu->
bins() && std::abs(dist2next.second) > 1. &&
2208 std::abs(dist2next.second) < std::abs(
h)) {
2209 h = dist2next.second * propDir;
2224 mom = std::abs(1. /
P[6]);
2225 sampleBrem(cache,
mom);
2228 while (numSf > 0 && (std::abs(distanceToTarget) > distanceTolerance ||
2229 std::abs(
path + distanceStepped) < tol)) {
2231 if (!rungeKuttaStep(cache, errorPropagation,
h,
P,
dDir, BG1, firstStep, distanceStepped)) {
2238 m_momentumCutOff, cache.
m_bremMom, position, direction,
2241 for (
int i = 0;
i < 3; ++
i)
2242 P[3 +
i] = direction[
i];
2249 mom = std::abs(1. /
P[6]);
2253 if (std::abs(distanceStepped) > 0.001) {
2276 absPath += std::abs(distanceStepped);
2279 mom = std::abs(1. /
P[6]);
2283 if (std::abs(distanceStepped) > 0.001)
2302 if (absPath > maxPath)
2311 bool restart =
false;
2313 if (propDir *
path < -tol || absPath - std::abs(
path) > 10.) {
2314 helpSoft = std::abs(
path) / absPath > 0.5 ? std::abs(
path) / absPath : 0.5;
2321 float distanceToNextBin =
h;
2327 std::pair<size_t, float> dist2next = lbu->
distanceToNext(position, propDir * direction);
2328 distanceToNextBin = dist2next.second;
2331 std::pair<size_t, float> dist2previous = lbu->
distanceToNext(position, -propDir * direction);
2332 float stepOver = dist2previous.second;
2335 auto cPar = std::make_unique<Trk::CurvilinearParameters>(
Amg::Vector3D(
P[0],
P[1],
P[2]), localp[2],
2336 localp[3], localp[4]);
2338 if (binIDMat && binIDMat->second > 0 && !iMat) {
2340 }
else if (binIDMat && binIDMat->second > 0 &&
2341 (iMat->second == 0 || iMat->second == binIDMat->second)) {
2343 }
else if (iMat && iMat->second > 0) {
2349 if (binIDMat && binIDMat->second > 0 && !iMat) {
2350 cache.
m_hitVector->emplace_back(cPar->uniqueClone(), hitTiming, -binIDMat->second, 0.);
2351 }
else if (binIDMat && binIDMat->second > 0 &&
2352 (iMat->second == 0 || iMat->second == binIDMat->second)) {
2353 cache.
m_hitVector->emplace_back(cPar->uniqueClone(), hitTiming, -binIDMat->second, 0.);
2354 }
else if (iMat && iMat->second > 0) {
2355 cache.
m_hitVector->emplace_back(cPar->uniqueClone(), hitTiming, iMat->second, 0.);
2368 updateMaterialEffects(cache,
mom,
sin(direction.theta()), sumPath +
path - stepOver);
2373 if (distanceToNextBin <
h) {
2374 Amg::Vector3D probe = position + (distanceToNextBin +
h) * propDir * direction;
2375 std::pair<size_t, float> d2n = lbu->
distanceToNext(probe, propDir * direction);
2376 distanceToNextBin += d2n.second +
h;
2378 }
else if (dist2next.first < lbu->
bins() && std::abs(distanceToNextBin) < 0.01 &&
2382 auto cPar = std::make_unique<Trk::CurvilinearParameters>(
Amg::Vector3D(
P[0],
P[1],
P[2]), localp[2],
2383 localp[3], localp[4]);
2387 Amg::Vector3D probe = position + (distanceToNextBin + 0.01) * propDir * direction.normalized();
2391 if (binIDMat && binIDMat->second > 0 && !nextMat) {
2393 }
else if (binIDMat && binIDMat->second > 0 &&
2394 (nextMat->second == 0 || nextMat->second == binIDMat->second)) {
2397 }
else if (nextMat && nextMat->second > 0) {
2403 if (binIDMat && binIDMat->second > 0 && !nextMat) {
2404 cache.
m_hitVector->emplace_back(cPar->uniqueClone(), hitTiming, -binIDMat->second, 0.);
2405 }
else if (binIDMat && binIDMat->second > 0 &&
2406 (nextMat->second == 0 ||
2407 nextMat->second == binIDMat->second)) {
2408 cache.
m_hitVector->emplace_back(cPar->uniqueClone(), hitTiming, -binIDMat->second, 0.);
2409 }
else if (nextMat && nextMat->second > 0) {
2410 cache.
m_hitVector->emplace_back(cPar->uniqueClone(), hitTiming, nextMat->second, 0.);
2415 if (binIDMat != nextMat) {
2419 updateMaterialEffects(cache,
mom,
sin(direction.theta()), sumPath +
path);
2424 std::pair<size_t, float> d2n = lbu->
distanceToNext(probe, propDir * direction.normalized());
2425 distanceToNextBin += d2n.second + 0.01;
2433 bool flipDirection =
false;
2435 nextSfCand = nextSf;
2436 double dev = direction0.dot(direction);
2438 std::vector<std::pair<int, std::pair<double, double>>>
::iterator vsIter = vsBeg;
2454 for (; vsIter != vsEnd; ++vsIter) {
2457 if (numRestart > restartLimit)
2463 distanceToTarget = propDir * maxPath;
2469 if ((*vsIter).first != -1 &&
2470 (
ic == nextSf || (*vsIter).first == 1 || nextSf < 0 || std::abs((*vsIter).second.first) < 500. ||
2471 std::abs(
path) > 0.5 * std::abs((*vsIter).second.second))) {
2472 previousDistance = (*vsIter).second.first;
2474 (*sIter).
first->straightLineDistanceEstimate(position, propDir * direction);
2475 double distanceEst = -propDir * maxPath;
2477 distanceEst = distSol.
first();
2479 std::abs(distSol.
first() * propDir + distanceStepped - previousDistance) >
2480 std::abs(distSol.
second() * propDir + distanceStepped - previousDistance)) {
2481 distanceEst = distSol.
second();
2488 if (
ic == startSf && distanceEst < 0 && distSol.
first() > 0)
2489 distanceEst = distSol.
first();
2492 if (
ic == nextSf && std::abs(distanceEst) < tol && std::abs(
path) < tol) {
2493 (*vsIter).first = -1;
2496 distanceToTarget = maxPath;
2504 if ((*vsIter).second.first * propDir * distanceEst < 0. &&
2505 std::abs(distanceEst) > distanceTolerance) {
2509 std::abs((*vsIter).second.second) < tol ||
2512 ((*vsIter).first)++;
2514 if ((*vsIter).first > 3) {
2515 helpSoft = fmax(0.05, 1. - 0.05 * (*vsIter).first);
2516 if ((*vsIter).first > 20)
2517 helpSoft = 1. / (*vsIter).first;
2520 if ((*vsIter).first > 50 &&
h * propDir > 0) {
2522 (*vsIter).first = -1;
2527 if ((*vsIter).first != -1)
2528 flipDirection =
true;
2529 }
else if (std::abs((*vsIter).second.second) > tol &&
2532 if (
ic > nextSf && nextSf != -1) {
2533 if (propDir * distanceEst < (cache.
m_currentDist.at(nextSf)).second.first - tol) {
2534 if ((*vsIter).first != -1) {
2535 ((*vsIter).first)++;
2536 flipDirection =
true;
2540 }
else if (distanceToTarget >
2542 if ((*vsIter).first != -1) {
2543 ((*vsIter).first)++;
2544 flipDirection =
true;
2549 }
else if (
ic == nextSf) {
2557 (*vsIter).second.first = propDir * distanceEst;
2563 if ((*vsIter).first != -1 && (distanceEst > 0. ||
ic == nextSf)) {
2565 if (distanceEst < std::abs(distanceToTarget)) {
2566 distanceToTarget = propDir * distanceEst;
2570 }
else if (std::abs(
path) > std::abs((*vsIter).second.second) || dev < 0.985 ||
2573 (*sIter).
first->straightLineDistanceEstimate(position, propDir * direction);
2574 double distanceEst = -propDir * maxPath;
2576 distanceEst = distSol.
first();
2579 (*vsIter).second.first = propDir * distanceEst;
2582 if (distanceEst > tol && distanceEst < maxPath) {
2583 (*vsIter).first = 0;
2587 if ((*vsIter).first != -1 && distanceEst > 0.) {
2589 if (distanceEst < std::abs(distanceToTarget)) {
2590 distanceToTarget = propDir * distanceEst;
2598 if (std::abs(distanceToTarget) <= distanceTolerance &&
path * propDir < distanceTolerance) {
2599 (*vsIter).first = -1;
2608 if (nextSf < 0 && nextSfCand < 0)
2611 if (flipDirection) {
2613 if (nextSf < 0 || nextSf >= num_vs_dist)
2615 distanceToTarget = (*(vsBeg + nextSf)).second.first;
2617 }
else if (nextSfCand != nextSf) {
2618 nextSf = nextSfCand;
2620 if (nextSf < 0 || nextSf >= num_vs_dist)
2627 if (std::abs(
h) > std::abs(distanceToTarget))
2628 h = distanceToTarget;
2631 if (cache.
m_binMat && std::abs(
h) > std::abs(distanceToNextBin) + 0.001) {
2632 if (distanceToNextBin > 0) {
2633 h = distanceToNextBin * propDir;
2645 if (std::abs(
path) > maxPath)
2659 mom = std::abs(1. /
P[6]);
2664 pos[0] =
pos[0] + distanceToTarget * (
dir[0] + 0.5 * distanceToTarget *
dDir[0]);
2665 pos[1] =
pos[1] + distanceToTarget * (
dir[1] + 0.5 * distanceToTarget *
dDir[1]);
2666 pos[2] =
pos[2] + distanceToTarget * (
dir[2] + 0.5 * distanceToTarget *
dDir[2]);
2669 dir[0] =
dir[0] + distanceToTarget *
dDir[0];
2670 dir[1] =
dir[1] + distanceToTarget *
dDir[1];
2671 dir[2] =
dir[2] + distanceToTarget *
dDir[2];
2683 std::vector<std::pair<int, std::pair<double, double>>>
::iterator vsIter = vsBeg;
2686 for (; vsIter != vsEnd; ++vsIter) {
2687 if ((*vsIter).first != -1 && propDir * (*vsIter).second.first >= propDir * distanceToTarget - tol &&
2688 propDir * (*vsIter).second.first < 0.01 &&
index != nextSf) {
2689 solutions.push_back(
index);
2691 if (
index == nextSf)
2692 solutions.push_back(
index);
2703 double path)
const {
2719 auto eloss = !m_detailedEloss
2722 : std::make_unique<Trk::EnergyLoss>(
2729 std::sqrt(cache.m_covariance(3, 3)));
2733 std::move(eloss), cvlTP->associatedSurface());
2741 cache.m_combinedCovariance += cache.m_covariance;
2742 cache.m_covariance.setZero();
2750 double radDist)
const {
2766 double th = std::sqrt(2.) * 15. * std::sqrt(radDist) / (
beta *
momentum) *
2780 theta = rotated.theta();
2781 phi = rotated.phi();
2803 if (fieldCondObj ==
nullptr) {
2804 ATH_MSG_ERROR(
"extrapolate: Failed to retrieve AtlasFieldCacheCondObj with key "
2805 << m_fieldCacheCondObjInputKey.key());
2811 CLHEP::HepRandomEngine*
2814 if (!m_simulation || !m_rngWrapper) {
2817 if (m_rngWrapper->evtSeeded(ctx) != ctx.evt()) {
2821 wrapper_nc->setSeed (this->
name(), ctx);
2823 return m_rngWrapper->getEngine (ctx);