packages/trackbase/MagneticFieldOptions.cc

0001
0002
0003 #include "MagneticFieldOptions.h"
0004
0005 #include <Acts/Definitions/Units.hpp>
0006 #include <Acts/MagneticField/BFieldMapUtils.hpp>
0007 #include <Acts/MagneticField/MagneticFieldProvider.hpp>
0008 #include <Acts/MagneticField/SolenoidBField.hpp>
0009 #include <Acts/Utilities/Logger.hpp>
0010 #include <ActsExamples/MagneticField/FieldMapRootIo.hpp>
0011 #include <ActsExamples/MagneticField/FieldMapTextIo.hpp>
0012 #include <ActsExamples/MagneticField/ScalableBFieldService.hpp>
0013 #include <ActsExamples/Utilities/Options.hpp>
0014
0015 #include <filesystem>
0016 #include <memory>
0017 #include <stdexcept>
0018 #include <string>
0019
0020 #include <boost/program_options.hpp>
0021
0022 void ActsExamples::Options::addMagneticFieldOptions(Description& desc) {
0023   using boost::program_options::bool_switch;
0024   using boost::program_options::value;
0025
0026   // avoid adding the options twice
0027   if (desc.find_nothrow("bf-constant-tesla", true) != nullptr) {
0028     return;
0029   }
0030
0031   auto opt = desc.add_options();
0032   opt("bf-constant-tesla", value<Reals<3>>(),
0033       "Set a constant magnetic field vector in Tesla. If given, this takes "
0034       "preference over all other options.");
0035   opt("bf-scalable", bool_switch(),
0036       "If given, the constant field strength will be scaled differently in "
0037       "every event. This is for testing only.");
0038   opt("bf-scalable-scalor", value<double>()->default_value(1.25),
0039       "Scaling factor for the event-dependent field strength scaling. A unit "
0040       "value means that the field strength stays the same for every event.");
0041   opt("bf-map-file", value<std::string>(),
0042       "Read a magnetic field map from the given file. ROOT and text file "
0043       "formats are supported. Only used if no constant field is given.");
0044   opt("bf-map-tree", value<std::string>()->default_value("bField"),
0045       "Name of the TTree in the ROOT file. Only used if the field map is read "
0046       "from a ROOT file.");
0047   opt("bf-map-type", value<std::string>()->default_value("xyz"),
0048       "Either 'xyz' or 'rz' to define the type of the field map.");
0049   opt("bf-map-octantonly", bool_switch(),
0050       "If given, the field map is assumed to describe only the first "
0051       "octant/quadrant and the field is symmetrically extended to the full "
0052       "space.");
0053   opt("bf-map-lengthscale-mm", value<double>()->default_value(1.),
0054       "Optional length scale modifier for the field map grid. This options "
0055       "only needs to be set if the length unit in the field map file is not "
0056       "`mm`. The value must scale from the stored unit to the equivalent value "
0057       "in `mm`.");
0058   opt("bf-map-fieldscale-tesla", value<double>()->default_value(1.),
0059       "Optional field value scale modifier for the field map value. This "
0060       "option only needs to be set if the field value unit in the field map "
0061       "file is not `Tesla`. The value must scale from the stored unit to the "
0062       "equivalent value in `Tesla`.");
0063   opt("bf-solenoid-mag-tesla", value<double>()->default_value(0.),
0064       "The magnitude of a solenoid magnetic field in the center in `Tesla`. "
0065       "Only used "
0066       "if neither constant field nor a magnetic field map is given.");
0067   opt("bf-solenoid-length", value<double>()->default_value(6000),
0068       "The length of the solenoid magnetic field in `mm`.");
0069   opt("bf-solenoid-radius", value<double>()->default_value(1200),
0070       "The radius of the solenoid magnetic field in `mm`.");
0071   opt("bf-solenoid-ncoils", value<size_t>()->default_value(1194),
0072       "Number of coils for the solenoid magnetic field.");
0073   opt("bf-solenoid-map-rlim",
0074       value<Interval>()->value_name("MIN:MAX")->default_value({0, 1200}),
0075       "The length bounds of the grid created from the analytical solenoid "
0076       "field in `mm`.");
0077   opt("bf-solenoid-map-zlim",
0078       value<Interval>()->value_name("MIN:MAX")->default_value({-3000, 3000}),
0079       "The radius bounds of the grid created from the analytical solenoid "
0080       "field in `mm`.");
0081   opt("bf-solenoid-map-nbins", value<Reals<2>>()->default_value({{150, 200}}),
0082       "The number of bins in r-z directions for the grid created from the "
0083       "analytical solenoid field.");
0084 }
0085
0086
0087 std::shared_ptr<Acts::MagneticFieldProvider>
0088 ActsExamples::Options::readMagneticField(const Variables& vars) {
0089   using namespace ActsExamples::detail;
0090   using std::filesystem::path;
0091
0092   // first option: create a constant field
0093   if (vars.count("bf-constant-tesla") != 0u) {
0094     const auto values = vars["bf-constant-tesla"].as<Reals<3>>();
0095     Acts::Vector3 field(values[0] * Acts::UnitConstants::T,
0096                         values[1] * Acts::UnitConstants::T,
0097                         values[2] * Acts::UnitConstants::T);
0098     if (vars["bf-scalable"].as<bool>()) {
0099       return std::make_shared<ScalableBField>(field);
0100     } else {
0101       return std::make_shared<Acts::ConstantBField>(field);
0102     }
0103   }
0104
0105   // second option: read a field map from a file
0106   if (vars.count("bf-map-file") != 0u) {
0107     const path file = vars["bf-map-file"].as<std::string>();
0108     const auto tree = vars["bf-map-tree"].as<std::string>();
0109     const auto type = vars["bf-map-type"].as<std::string>();
0110     const auto useOctantOnly = vars["bf-map-octantonly"].as<bool>();
0111     const auto lengthUnit =
0112         vars["bf-map-lengthscale-mm"].as<double>() * Acts::UnitConstants::mm;
0113     const auto fieldUnit =
0114         vars["bf-map-fieldscale-tesla"].as<double>() * Acts::UnitConstants::T;
0115
0116     bool readRoot = false;
0117     if (file.extension() == ".root") {
0118       readRoot = true;
0119     } else if (file.extension() == ".txt") {
0120       readRoot = false;
0121     } else {
0122       throw std::runtime_error("Unsupported magnetic field map file type");
0123     }
0124
0125     if (type == "xyz") {
0126       auto mapBins = [](const std::array<size_t, 3>& bins,
0127                         const std::array<size_t, 3>& sizes) {
0128         return (bins[0] * (sizes[1] * sizes[2]) + bins[1] * sizes[2] + bins[2]);
0129       };
0130
0131       if (readRoot) {
0132         auto map = makeMagneticFieldMapXyzFromRoot(
0133             std::move(mapBins), file.native(), tree, lengthUnit, fieldUnit,
0134             useOctantOnly);
0135         return std::make_shared<InterpolatedMagneticField3>(std::move(map));
0136
0137       } else {
0138         auto map = makeMagneticFieldMapXyzFromText(std::move(mapBins),
0139                                                    file.native(), lengthUnit,
0140                                                    fieldUnit, useOctantOnly);
0141         return std::make_shared<InterpolatedMagneticField3>(std::move(map));
0142       }
0143
0144     } else if (type == "rz") {
0145       auto mapBins = [](std::array<size_t, 2> bins,
0146                         std::array<size_t, 2> sizes) {
0147         return (bins[1] * sizes[0] + bins[0]);
0148       };
0149
0150       if (readRoot) {
0151         auto map = makeMagneticFieldMapRzFromRoot(
0152             std::move(mapBins), file.native(), tree, lengthUnit, fieldUnit,
0153             useOctantOnly);
0154         return std::make_shared<InterpolatedMagneticField2>(std::move(map));
0155
0156       } else {
0157         auto map = makeMagneticFieldMapRzFromText(std::move(mapBins),
0158                                                   file.native(), lengthUnit,
0159                                                   fieldUnit, useOctantOnly);
0160         return std::make_shared<InterpolatedMagneticField2>(std::move(map));
0161       }
0162
0163     } else {
0164       throw std::runtime_error("Unknown magnetic field map type");
0165     }
0166   }
0167
0168   // third option: create a solenoid field
0169   if (vars["bf-solenoid-mag-tesla"].as<double>() > 0) {
0170     // Construct a solenoid field
0171     Acts::SolenoidBField::Config solenoidConfig{};
0172     solenoidConfig.length =
0173         vars["bf-solenoid-length"].as<double>() * Acts::UnitConstants::mm;
0174     solenoidConfig.radius =
0175         vars["bf-solenoid-radius"].as<double>() * Acts::UnitConstants::mm;
0176     solenoidConfig.nCoils = vars["bf-solenoid-ncoils"].as<size_t>();
0177     solenoidConfig.bMagCenter =
0178         vars["bf-solenoid-mag-tesla"].as<double>() * Acts::UnitConstants::T;
0179
0180     const auto solenoidField = Acts::SolenoidBField(solenoidConfig);
0181     // The parameters for creating a field map
0182     auto getRange = [&](const char* name, auto unit, auto& lower, auto& upper) {
0183       auto interval = vars[name].as<Options::Interval>();
0184       lower = interval.lower.value() * unit;
0185       upper = interval.upper.value() * unit;
0186     };
0187     std::pair<double, double> rlim, zlim;
0188     getRange("bf-solenoid-map-rlim", Acts::UnitConstants::mm, rlim.first,
0189              rlim.second);
0190     getRange("bf-solenoid-map-zlim", Acts::UnitConstants::mm, zlim.first,
0191              zlim.second);
0192     const auto nbins = vars["bf-solenoid-map-nbins"].as<Reals<2>>();
0193     auto map =
0194         Acts::solenoidFieldMap(rlim, zlim, {nbins[0], nbins[1]}, solenoidField);
0195     return std::make_shared<InterpolatedMagneticField2>(std::move(map));
0196   }
0197
0198   // default option: no field
0199   return std::make_shared<Acts::NullBField>();
0200 }