algorithms/tracking/DD4hepBField.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 - 2025 Whitney Armstrong, Wouter Deconinck, Dmitry Kalinkin
0003
0004 #include "DD4hepBField.h"
0005
0006 #include <Acts/Definitions/Algebra.hpp>
0007 #include <Acts/Definitions/Units.hpp>
0008 #include <Acts/Utilities/Result.hpp>
0009 #include <DD4hep/Fields.h>
0010 #include <DD4hep/Objects.h>
0011 #include <Evaluator/DD4hepUnits.h>
0012 #include <Math/GenVector/Cartesian3D.h>
0013 #include <Math/GenVector/DisplacementVector3D.h>
0014 #include <Eigen/Core>
0015 #include <cmath>
0016 #include <limits>
0017
0018 namespace eicrecon::BField {
0019
0020 Acts::Result<Acts::Vector3>
0021 DD4hepBField::getField(const Acts::Vector3& position,
0022                        Acts::MagneticFieldProvider::Cache& /*cache*/) const {
0023   dd4hep::Position pos(position[0] * (dd4hep::mm / Acts::UnitConstants::mm),
0024                        position[1] * (dd4hep::mm / Acts::UnitConstants::mm),
0025                        position[2] * (dd4hep::mm / Acts::UnitConstants::mm));
0026
0027   // Avoid crash during lookup.
0028   // Infinite values are possible due to https://github.com/acts-project/acts/issues/4166.
0029   if ((!std::isfinite(position[0])) && (!std::isfinite(position[1])) &&
0030       (!std::isfinite(position[2]))) {
0031     return Acts::Result<Acts::Vector3>::success({0., 0., 0.});
0032   }
0033
0034   auto fieldObj = m_det->field();
0035   auto field    = fieldObj.magneticField(pos) * (Acts::UnitConstants::T / dd4hep::tesla);
0036
0037   // FIXME Acts doesn't seem to like exact zero components
0038   if (field.x() * field.y() * field.z() == 0) {
0039     static dd4hep::Direction epsilon{std::numeric_limits<double>::epsilon(),
0040                                      std::numeric_limits<double>::epsilon(),
0041                                      std::numeric_limits<double>::epsilon()};
0042     field += epsilon;
0043   }
0044
0045   return Acts::Result<Acts::Vector3>::success({field.x(), field.y(), field.z()});
0046 }
0047
0048 #if Acts_VERSION_MAJOR < 39
0049 Acts::Result<Acts::Vector3>
0050 DD4hepBField::getFieldGradient(const Acts::Vector3& position,
0051                                Acts::ActsMatrix<3, 3>& /*derivative*/,
0052                                Acts::MagneticFieldProvider::Cache& cache) const {
0053   return this->getField(position, cache);
0054 }
0055 #endif
0056
0057 } // namespace eicrecon::BField