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 "ActsPlugins/Root/RootMagneticFieldIo.hpp"
0019 #include "ActsPython/Utilities/Helpers.hpp"
0020 #include "ActsPython/Utilities/Macros.hpp"
0021
0022 #include <array>
0023 #include <cstddef>
0024 #include <filesystem>
0025 #include <memory>
0026 #include <stdexcept>
0027 #include <string>
0028 #include <tuple>
0029 #include <utility>
0030
0031 #include <pybind11/pybind11.h>
0032 #include <pybind11/stl.h>
0033
0034 namespace py = pybind11;
0035 using namespace pybind11::literals;
0036
0037 using namespace Acts;
0038 using namespace ActsExamples;
0039
0040 namespace ActsPython {
0041
0042 /// @brief Get the value of a field, throwing an exception if the result is
0043 /// invalid.
0044 Vector3 getField(MagneticFieldProvider& self, const Vector3& position,
0045                  MagneticFieldProvider::Cache& cache) {
0046   if (Result<Vector3> res = self.getField(position, cache); !res.ok()) {
0047     std::stringstream ss;
0048
0049     ss << "Field lookup failure with error: \"" << res.error() << "\"";
0050
0051     throw std::runtime_error{ss.str()};
0052   } else {
0053     return *res;
0054   }
0055 }
0056
0057 /// @brief Add the magnetic field bindings to a module.
0058 /// @param m the module to add the bindings to
0059 void addMagneticField(py::module_& m) {
0060   py::class_<Acts::MagneticFieldContext>(m, "MagneticFieldContext")
0061       .def(py::init<>());
0062
0063   py::class_<MagneticFieldProvider, std::shared_ptr<MagneticFieldProvider>>(
0064       m, "MagneticFieldProvider")
0065       .def("getField", &getField)
0066       .def("makeCache", &MagneticFieldProvider::makeCache);
0067
0068   py::class_<InterpolatedMagneticField,
0069              std::shared_ptr<InterpolatedMagneticField>>(
0070       m, "InterpolatedMagneticField");
0071
0072   m.def("solenoidFieldMap", &solenoidFieldMap, py::arg("rlim"), py::arg("zlim"),
0073         py::arg("nbins"), py::arg("field"));
0074
0075   py::class_<ConstantBField, MagneticFieldProvider,
0076              std::shared_ptr<ConstantBField>>(m, "ConstantBField")
0077       .def(py::init<Vector3>());
0078
0079   using InterpolatedMagneticField2 = InterpolatedBFieldMap<
0080       Grid<Vector2, Axis<AxisType::Equidistant>, Axis<AxisType::Equidistant>>>;
0081
0082   using InterpolatedMagneticField3 = InterpolatedBFieldMap<
0083       Grid<Vector3, Axis<AxisType::Equidistant>, Axis<AxisType::Equidistant>,
0084            Axis<AxisType::Equidistant>>>;
0085
0086   py::class_<InterpolatedMagneticField2, InterpolatedMagneticField,
0087              MagneticFieldProvider,
0088              std::shared_ptr<InterpolatedMagneticField2>>(
0089       m, "InterpolatedMagneticField2");
0090
0091   py::class_<InterpolatedMagneticField3, InterpolatedMagneticField,
0092              MagneticFieldProvider,
0093              std::shared_ptr<InterpolatedMagneticField3>>(
0094       m, "InterpolatedMagneticField3");
0095
0096   py::class_<NullBField, MagneticFieldProvider, std::shared_ptr<NullBField>>(
0097       m, "NullBField")
0098       .def(py::init<>());
0099
0100   py::class_<MultiRangeBField, MagneticFieldProvider,
0101              std::shared_ptr<MultiRangeBField>>(m, "MultiRangeBField")
0102       .def(py::init<std::vector<std::pair<RangeXD<3, double>, Vector3>>>());
0103
0104   {
0105     using Config = SolenoidBField::Config;
0106
0107     auto sol = py::class_<SolenoidBField, MagneticFieldProvider,
0108                           std::shared_ptr<SolenoidBField>>(m, "SolenoidBField")
0109                    .def(py::init<Config>())
0110                    .def(py::init([](double radius, double length,
0111                                     std::size_t nCoils, double bMagCenter) {
0112                           return SolenoidBField{
0113                               Config{radius, length, nCoils, bMagCenter}};
0114                         }),
0115                         py::arg("radius"), py::arg("length"), py::arg("nCoils"),
0116                         py::arg("bMagCenter"));
0117
0118     py::class_<Config>(sol, "Config")
0119         .def(py::init<>())
0120         .def_readwrite("radius", &Config::radius)
0121         .def_readwrite("length", &Config::length)
0122         .def_readwrite("nCoils", &Config::nCoils)
0123         .def_readwrite("bMagCenter", &Config::bMagCenter);
0124   }
0125
0126   m.def(
0127       "MagneticFieldMapXyz",
0128       [](const std::string& filename, const std::string& tree,
0129          double lengthUnit, double BFieldUnit, bool firstOctant) {
0130         const std::filesystem::path file = filename;
0131
0132         auto mapBins = [](std::array<std::size_t, 3> bins,
0133                           std::array<std::size_t, 3> sizes) {
0134           return (bins[0] * (sizes[1] * sizes[2]) + bins[1] * sizes[2] +
0135                   bins[2]);
0136         };
0137
0138         if (file.extension() == ".root") {
0139           auto map = ActsPlugins::makeMagneticFieldMapXyzFromRoot(
0140               std::move(mapBins), file.native(), tree, lengthUnit, BFieldUnit,
0141               firstOctant);
0142           return std::make_shared<decltype(map)>(std::move(map));
0143         } else if (file.extension() == ".txt") {
0144           auto map = makeMagneticFieldMapXyzFromText(std::move(mapBins),
0145                                                      file.native(), lengthUnit,
0146                                                      BFieldUnit, firstOctant);
0147           return std::make_shared<decltype(map)>(std::move(map));
0148         } else {
0149           throw std::runtime_error("Unsupported magnetic field map file type");
0150         }
0151       },
0152       py::arg("file"), py::arg("tree") = "bField",
0153       py::arg("lengthUnit") = UnitConstants::mm,
0154       py::arg("BFieldUnit") = UnitConstants::T, py::arg("firstOctant") = false);
0155
0156   m.def(
0157       "MagneticFieldMapRz",
0158       [](const std::string& filename, const std::string& tree,
0159          double lengthUnit, double BFieldUnit, bool firstQuadrant) {
0160         const std::filesystem::path file = filename;
0161
0162         auto mapBins = [](std::array<std::size_t, 2> bins,
0163                           std::array<std::size_t, 2> sizes) {
0164           return (bins[1] * sizes[0] + bins[0]);
0165         };
0166
0167         if (file.extension() == ".root") {
0168           auto map = ActsPlugins::makeMagneticFieldMapRzFromRoot(
0169               std::move(mapBins), file.native(), tree, lengthUnit, BFieldUnit,
0170               firstQuadrant);
0171           return std::make_shared<decltype(map)>(std::move(map));
0172         } else if (file.extension() == ".txt") {
0173           auto map = makeMagneticFieldMapRzFromText(std::move(mapBins),
0174                                                     file.native(), lengthUnit,
0175                                                     BFieldUnit, firstQuadrant);
0176           return std::make_shared<decltype(map)>(std::move(map));
0177         } else {
0178           throw std::runtime_error("Unsupported magnetic field map file type");
0179         }
0180       },
0181       py::arg("file"), py::arg("tree") = "bField",
0182       py::arg("lengthUnit") = UnitConstants::mm,
0183       py::arg("BFieldUnit") = UnitConstants::T,
0184       py::arg("firstQuadrant") = false);
0185 }
0186
0187 }  // namespace ActsPython