File indexing completed on 2025-11-03 08:57:58
0001 
0002 
0003 
0004 
0005 
0006 
0007 
0008 
0009 #include "ActsExamples/Io/Csv/CsvBFieldWriter.hpp"
0010 
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Units.hpp"
0013 #include "Acts/MagneticField/InterpolatedBFieldMap.hpp"
0014 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0015 #include "Acts/Utilities/Logger.hpp"
0016 #include "Acts/Utilities/Result.hpp"
0017 #include "Acts/Utilities/VectorHelpers.hpp"
0018 #include "ActsExamples/Io/Csv/CsvInputOutput.hpp"
0019 
0020 #include <iomanip>
0021 #include <ostream>
0022 #include <stdexcept>
0023 #include <vector>
0024 
0025 namespace ActsExamples {
0026 template <CsvBFieldWriter::CoordinateType Coord, bool Grid>
0027 void CsvBFieldWriter::run(const Config<Coord, Grid>& config,
0028                           Acts::Logging::Level level) {
0029   ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("CsvBFieldWriter", level));
0030 
0031   
0032   using ConfigType = std::decay_t<decltype(config)>;
0033   using FieldType = typename decltype(ConfigType::bField)::element_type;
0034   using Vector = Acts::ActsVector<ConfigType::NDims>;
0035 
0036   FieldType& field = *config.bField;
0037 
0038   
0039   
0040   
0041   std::vector<std::string> fields;
0042 
0043   if constexpr (Coord == CoordinateType::XYZ) {
0044     fields = {"x", "y", "z", "Bx", "By", "Bz"};
0045   } else {
0046     fields = {"r", "z", "Br", "Bz"};
0047   }
0048 
0049   
0050   
0051   CsvWriter writer(fields, config.fileName);
0052 
0053   
0054   
0055   
0056   std::array<std::size_t, ConfigType::NDims> bins{};
0057   Vector min, max;
0058 
0059   if constexpr (Grid) {
0060     
0061     
0062     if (bins.size() != ConfigType::NDims ||
0063         field.getMin().size() != ConfigType::NDims ||
0064         field.getMax().size() != ConfigType::NDims) {
0065       throw std::invalid_argument("Incorrect input B-field dimensionality.");
0066     }
0067 
0068     
0069     
0070     
0071     for (std::size_t i = 0; i < ConfigType::NDims; ++i) {
0072       
0073       
0074       bins[i] = field.getNBins()[i];
0075       min[i] = field.getMin()[i];
0076       max[i] = field.getMax()[i];
0077 
0078       
0079       if (config.bins[i]) {
0080         bins[i] = *config.bins[i];
0081       }
0082 
0083       if (config.range[i][0]) {
0084         min[i] = *config.range[i][0];
0085       }
0086 
0087       if (config.range[i][1]) {
0088         max[i] = *config.range[i][1];
0089       }
0090     }
0091   } else {
0092     
0093     
0094     
0095     for (std::size_t i = 0; i < ConfigType::NDims; ++i) {
0096       bins[i] = config.bins[i];
0097       min[i] = config.range[i][0];
0098       max[i] = config.range[i][1];
0099     }
0100   }
0101 
0102   
0103   Vector delta;
0104 
0105   for (std::size_t i = 0; i < ConfigType::NDims; ++i) {
0106     delta[i] = (max[i] - min[i]) / (bins[i] - 1);
0107   }
0108 
0109   
0110   
0111   Acts::MagneticFieldContext mctx{};
0112   typename FieldType::Cache cache = field.makeCache(mctx);
0113 
0114   
0115   
0116   
0117   if constexpr (Coord == CoordinateType::XYZ) {
0118     ACTS_INFO("Writing XYZ field of size " << bins[0] << " x " << bins[1]
0119                                            << " x " << bins[2] << " to file "
0120                                            << config.fileName << "...");
0121 
0122     std::size_t total_items = bins[0] * bins[1] * bins[2];
0123 
0124     
0125     
0126     
0127     for (std::size_t x = 0; x < bins[0]; ++x) {
0128       for (std::size_t y = 0; y < bins[1]; ++y) {
0129         for (std::size_t z = 0; z < bins[2]; ++z) {
0130           
0131           
0132           Acts::Vector3 pos = {x * delta[0] + min[0], y * delta[1] + min[1],
0133                                z * delta[2] + min[2]};
0134 
0135           Acts::Vector3 bField;
0136           if (auto fieldMap =
0137                   dynamic_cast<const Acts::InterpolatedMagneticField*>(
0138                       &field)) {
0139             
0140             
0141             
0142             bField = fieldMap->getFieldUnchecked(pos);
0143           } else {
0144             Acts::Result<Acts::Vector3> flx = field.getField(pos, cache);
0145 
0146             
0147             
0148             
0149             if (flx.ok()) {
0150               bField = *flx;
0151             } else {
0152               throw std::runtime_error("B-field returned a non-extant value!");
0153             }
0154           }
0155 
0156           writer.append(pos[0] / Acts::UnitConstants::mm,
0157                         pos[1] / Acts::UnitConstants::mm,
0158                         pos[2] / Acts::UnitConstants::mm,
0159                         bField[0] / Acts::UnitConstants::T,
0160                         bField[1] / Acts::UnitConstants::T,
0161                         bField[2] / Acts::UnitConstants::T);
0162 
0163           
0164           
0165           
0166           std::size_t idx = (x * bins[1] * bins[2]) + (y * bins[2]) + z + 1;
0167 
0168           if (idx % 10000 == 0 || idx == total_items) {
0169             ACTS_VERBOSE("Wrote " << idx << " out of " << total_items
0170                                   << " items (" << std::setprecision(3)
0171                                   << ((100.f * idx) / total_items) << "%).");
0172           }
0173         }
0174       }
0175     }
0176   } else {
0177     ACTS_INFO("Writing RZ field of size " << bins[0] << " x " << bins[1]
0178                                           << " to file " << config.fileName
0179                                           << "...");
0180 
0181     std::size_t total_items = bins[0] * bins[1];
0182 
0183     
0184     
0185     
0186     for (std::size_t r = 0; r < bins[0]; ++r) {
0187       for (std::size_t z = 0; z < bins[1]; ++z) {
0188         
0189         
0190         Acts::Vector3 pos(min[0] + r * delta[0], 0.f, min[1] + z * delta[1]);
0191 
0192         Acts::Vector3 bField;
0193         if (auto fieldMap =
0194                 dynamic_cast<const Acts::InterpolatedMagneticField*>(&field)) {
0195           
0196           
0197           
0198           bField = fieldMap->getFieldUnchecked(pos);
0199         } else {
0200           Acts::Result<Acts::Vector3> flx = field.getField(pos, cache);
0201 
0202           
0203           
0204           
0205           
0206           if (flx.ok()) {
0207             bField = *flx;
0208           } else {
0209             throw std::runtime_error("B-field returned a non-extant value!");
0210           }
0211         }
0212 
0213         writer.append(
0214             pos[0] / Acts::UnitConstants::mm, pos[2] / Acts::UnitConstants::mm,
0215             Acts::VectorHelpers::perp(bField) / Acts::UnitConstants::T,
0216             bField[2] / Acts::UnitConstants::T);
0217 
0218         
0219         
0220         std::size_t idx = (r * bins[1]) + z + 1;
0221 
0222         if (idx % 10000 == 0 || idx == total_items) {
0223           ACTS_VERBOSE("Wrote " << idx << " out of " << total_items
0224                                 << " items (" << std::setprecision(3)
0225                                 << ((100.f * idx) / total_items) << "%).");
0226         }
0227       }
0228     }
0229   }
0230 }
0231 
0232 
0233 
0234 
0235 
0236 
0237 template void CsvBFieldWriter::run<CsvBFieldWriter::CoordinateType::XYZ, true>(
0238     const Config<CoordinateType::XYZ, true>&, Acts::Logging::Level);
0239 template void CsvBFieldWriter::run<CsvBFieldWriter::CoordinateType::RZ, true>(
0240     const Config<CoordinateType::RZ, true>&, Acts::Logging::Level);
0241 template void CsvBFieldWriter::run<CsvBFieldWriter::CoordinateType::XYZ, false>(
0242     const Config<CoordinateType::XYZ, false>&, Acts::Logging::Level);
0243 template void CsvBFieldWriter::run<CsvBFieldWriter::CoordinateType::RZ, false>(
0244     const Config<CoordinateType::RZ, false>&, Acts::Logging::Level);
0245 }