Core/src/MagneticField.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "Acts/Definitions/Units.hpp"
0010 #include "Acts/MagneticField/BFieldMapUtils.hpp"
0011 #include "Acts/MagneticField/ConstantBField.hpp"
0012 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0013 #include "Acts/MagneticField/MagneticFieldProvider.hpp"
0014 #include "Acts/MagneticField/MultiRangeBField.hpp"
0015 #include "Acts/MagneticField/NullBField.hpp"
0016 #include "Acts/MagneticField/SolenoidBField.hpp"
0017 #include "Acts/MagneticField/TextMagneticFieldIo.hpp"
0018 #include "ActsPython/Utilities/Helpers.hpp"
0019 #include "ActsPython/Utilities/Macros.hpp"
0020
0021 #include <array>
0022 #include <cstddef>
0023 #include <filesystem>
0024 #include <memory>
0025 #include <stdexcept>
0026 #include <string>
0027 #include <tuple>
0028 #include <utility>
0029
0030 #include <pybind11/pybind11.h>
0031 #include <pybind11/stl.h>
0032
0033 namespace py = pybind11;
0034 using namespace pybind11::literals;
0035
0036 using namespace Acts;
0037 using namespace ActsExamples;
0038
0039 namespace ActsPython {
0040
0041 /// @brief Get the value of a field, throwing an exception if the result is
0042 /// invalid.
0043 Vector3 getField(MagneticFieldProvider& self, const Vector3& position,
0044                  MagneticFieldProvider::Cache& cache) {
0045   if (Result<Vector3> res = self.getField(position, cache); !res.ok()) {
0046     std::stringstream ss;
0047
0048     ss << "Field lookup failure with error: \"" << res.error() << "\"";
0049
0050     throw std::runtime_error{ss.str()};
0051   } else {
0052     return *res;
0053   }
0054 }
0055
0056 /// @brief Add the magnetic field bindings to a module.
0057 /// @param m the module to add the bindings to
0058 void addMagneticField(py::module_& m) {
0059   py::class_<Acts::MagneticFieldContext>(m, "MagneticFieldContext")
0060       .def(py::init<>());
0061
0062   py::class_<MagneticFieldProvider, std::shared_ptr<MagneticFieldProvider>>(
0063       m, "MagneticFieldProvider")
0064       .def("getField", &getField)
0065       .def("makeCache", &MagneticFieldProvider::makeCache);
0066
0067   py::class_<InterpolatedMagneticField,
0068              std::shared_ptr<InterpolatedMagneticField>>(
0069       m, "InterpolatedMagneticField");
0070
0071   m.def("solenoidFieldMap", &solenoidFieldMap, py::arg("rlim"), py::arg("zlim"),
0072         py::arg("nbins"), py::arg("field"));
0073
0074   py::class_<ConstantBField, MagneticFieldProvider,
0075              std::shared_ptr<ConstantBField>>(m, "ConstantBField")
0076       .def(py::init<Vector3>());
0077
0078   using InterpolatedMagneticField2 = InterpolatedBFieldMap<
0079       Grid<Vector2, Axis<AxisType::Equidistant>, Axis<AxisType::Equidistant>>>;
0080
0081   using InterpolatedMagneticField3 = InterpolatedBFieldMap<
0082       Grid<Vector3, Axis<AxisType::Equidistant>, Axis<AxisType::Equidistant>,
0083            Axis<AxisType::Equidistant>>>;
0084
0085   py::class_<InterpolatedMagneticField2, InterpolatedMagneticField,
0086              MagneticFieldProvider,
0087              std::shared_ptr<InterpolatedMagneticField2>>(
0088       m, "InterpolatedMagneticField2");
0089
0090   py::class_<InterpolatedMagneticField3, InterpolatedMagneticField,
0091              MagneticFieldProvider,
0092              std::shared_ptr<InterpolatedMagneticField3>>(
0093       m, "InterpolatedMagneticField3");
0094
0095   py::class_<NullBField, MagneticFieldProvider, std::shared_ptr<NullBField>>(
0096       m, "NullBField")
0097       .def(py::init<>());
0098
0099   py::class_<MultiRangeBField, MagneticFieldProvider,
0100              std::shared_ptr<MultiRangeBField>>(m, "MultiRangeBField")
0101       .def(py::init<std::vector<std::pair<RangeXD<3, double>, Vector3>>>());
0102
0103   {
0104     using Config = SolenoidBField::Config;
0105
0106     auto sol = py::class_<SolenoidBField, MagneticFieldProvider,
0107                           std::shared_ptr<SolenoidBField>>(m, "SolenoidBField")
0108                    .def(py::init<Config>())
0109                    .def(py::init([](double radius, double length,
0110                                     std::size_t nCoils, double bMagCenter) {
0111                           return SolenoidBField{
0112                               Config{radius, length, nCoils, bMagCenter}};
0113                         }),
0114                         py::arg("radius"), py::arg("length"), py::arg("nCoils"),
0115                         py::arg("bMagCenter"));
0116
0117     py::class_<Config>(sol, "Config")
0118         .def(py::init<>())
0119         .def_readwrite("radius", &Config::radius)
0120         .def_readwrite("length", &Config::length)
0121         .def_readwrite("nCoils", &Config::nCoils)
0122         .def_readwrite("bMagCenter", &Config::bMagCenter);
0123   }
0124
0125   m.def(
0126       "MagneticFieldMapXyz",
0127       [](const std::string& filename, double lengthUnit, double BFieldUnit,
0128          bool firstOctant) {
0129         const std::filesystem::path file = filename;
0130
0131         auto mapBins = [](std::array<std::size_t, 3> bins,
0132                           std::array<std::size_t, 3> sizes) {
0133           return (bins[0] * (sizes[1] * sizes[2]) + bins[1] * sizes[2] +
0134                   bins[2]);
0135         };
0136
0137         if (file.extension() == ".txt") {
0138           auto map = makeMagneticFieldMapXyzFromText(std::move(mapBins),
0139                                                      file.native(), lengthUnit,
0140                                                      BFieldUnit, firstOctant);
0141           return std::make_shared<decltype(map)>(std::move(map));
0142         } else {
0143           throw std::runtime_error("Unsupported magnetic field map file type");
0144         }
0145       },
0146       py::arg("file"), py::arg("lengthUnit") = UnitConstants::mm,
0147       py::arg("BFieldUnit") = UnitConstants::T, py::arg("firstOctant") = false);
0148
0149   m.def(
0150       "MagneticFieldMapRz",
0151       [](const std::string& filename, double lengthUnit, double BFieldUnit,
0152          bool firstQuadrant) {
0153         const std::filesystem::path file = filename;
0154
0155         auto mapBins = [](std::array<std::size_t, 2> bins,
0156                           std::array<std::size_t, 2> sizes) {
0157           return (bins[1] * sizes[0] + bins[0]);
0158         };
0159
0160         if (file.extension() == ".txt") {
0161           auto map = makeMagneticFieldMapRzFromText(std::move(mapBins),
0162                                                     file.native(), lengthUnit,
0163                                                     BFieldUnit, firstQuadrant);
0164           return std::make_shared<decltype(map)>(std::move(map));
0165         } else {
0166           throw std::runtime_error("Unsupported magnetic field map file type");
0167         }
0168       },
0169       py::arg("file"), py::arg("lengthUnit") = UnitConstants::mm,
0170       py::arg("BFieldUnit") = UnitConstants::T,
0171       py::arg("firstQuadrant") = false);
0172 }
0173
0174 }  // namespace ActsPython