examples/tools/convert_csv_bfield.cpp

0001 /** TRACCC library, part of the ACTS project (R&D line)
0002  *
0003  * (c) 2025 CERN for the benefit of the ACTS project
0004  *
0005  * Mozilla Public License Version 2.0
0006  */
0007
0008 #include <fstream>
0009 #include <iostream>
0010 #include <sstream>
0011
0012 #include <boost/program_options.hpp>
0013
0014 #include "bfield_type.hpp"
0015 #include "traccc/definitions/common.hpp"
0016
0017 // TODO: Remove this when we have logging
0018 #define TRACCC_INFO(x)                        \
0019   do {                                        \
0020     std::cout << "INFO : " << x << std::endl; \
0021   } while (0)
0022 #define TRACCC_FATAL(x)                        \
0023   do {                                         \
0024     std::cout << "FATAL : " << x << std::endl; \
0025   } while (0)
0026
0027 #define CHECK_IOSTATE(x)                                               \
0028   do {                                                                 \
0029     if (!x.good()) {                                                   \
0030       std::stringstream _exception_ss;                                 \
0031       _exception_ss << "Stringstream is not good on line " << __LINE__ \
0032                     << " of file " << __FILE__;                        \
0033       throw std::runtime_error(_exception_ss.str());                   \
0034     }                                                                  \
0035   } while (0)
0036
0037 void parse_opts(int argc, char* argv[],
0038                 boost::program_options::variables_map& vm) {
0039   boost::program_options::options_description opts("general options");
0040
0041   opts.add_options()("help", "produce help message")(
0042       "input,i", boost::program_options::value<std::string>()->required(),
0043       "input magnetic field to read")(
0044       "output,o", boost::program_options::value<std::string>()->required(),
0045       "output magnetic field to write")(
0046       "delimiter,d",
0047       boost::program_options::value<char>()->default_value(' ')->required(),
0048       "delimiter character");
0049
0050   boost::program_options::parsed_options parsed =
0051       boost::program_options::command_line_parser(argc, argv)
0052           .options(opts)
0053           .run();
0054
0055   boost::program_options::store(parsed, vm);
0056
0057   if (vm.count("help")) {
0058     std::cout << opts << std::endl;
0059     std::exit(0);
0060   }
0061
0062   try {
0063     boost::program_options::notify(vm);
0064   } catch (boost::program_options::required_option& e) {
0065     TRACCC_FATAL(e.what());
0066     std::exit(1);
0067   }
0068 }
0069
0070 field_t read_bfield(const std::string& fn, char delimiter) {
0071   std::ifstream f;
0072
0073   float minx = std::numeric_limits<float>::max();
0074   float maxx = std::numeric_limits<float>::lowest();
0075   float miny = std::numeric_limits<float>::max();
0076   float maxy = std::numeric_limits<float>::lowest();
0077   float minz = std::numeric_limits<float>::max();
0078   float maxz = std::numeric_limits<float>::lowest();
0079
0080   {
0081     TRACCC_INFO("Opening magnetic field to compute field limits");
0082
0083     f.open(fn);
0084
0085     if (!f.good()) {
0086       TRACCC_FATAL("Failed to open input file " << fn << "!");
0087       std::exit(1);
0088     }
0089
0090     std::string line;
0091
0092     TRACCC_INFO("Skipping the first four lines (headers)");
0093
0094     for (std::size_t i = 0; i < 4; ++i) {
0095       std::getline(f, line);
0096     }
0097
0098     float xp, yp, zp;
0099     float Bx, By, Bz;
0100
0101     (void)Bx, (void)By, (void)Bz;
0102
0103     std::size_t n_lines = 0;
0104
0105     TRACCC_INFO("Iterating over lines in the magnetic field file");
0106
0107     /*
0108      * Read every line, and update our current minima and maxima
0109      * appropriately.
0110      */
0111     while (std::getline(f, line)) {
0112       CHECK_IOSTATE(f);
0113
0114       std::string word;
0115       std::stringstream ss(line);
0116       CHECK_IOSTATE(ss);
0117
0118       std::getline(ss, word, delimiter);
0119       CHECK_IOSTATE(ss);
0120       xp = static_cast<float>(std::atof(word.c_str()));
0121       std::getline(ss, word, delimiter);
0122       CHECK_IOSTATE(ss);
0123       yp = static_cast<float>(std::atof(word.c_str()));
0124       std::getline(ss, word, delimiter);
0125       CHECK_IOSTATE(ss);
0126       zp = static_cast<float>(std::atof(word.c_str()));
0127       std::getline(ss, word, delimiter);
0128       CHECK_IOSTATE(ss);
0129       Bx = static_cast<float>(std::atof(word.c_str()));
0130       std::getline(ss, word, delimiter);
0131       CHECK_IOSTATE(ss);
0132       By = static_cast<float>(std::atof(word.c_str()));
0133       std::getline(ss, word);
0134       Bz = static_cast<float>(std::atof(word.c_str()));
0135
0136       minx = std::min(minx, xp);
0137       maxx = std::max(maxx, xp);
0138
0139       miny = std::min(miny, yp);
0140       maxy = std::max(maxy, yp);
0141
0142       minz = std::min(minz, zp);
0143       maxz = std::max(maxz, zp);
0144
0145       ++n_lines;
0146     }
0147
0148     TRACCC_INFO("Read " << n_lines << " lines of magnetic field data");
0149
0150     TRACCC_INFO("Closing magnetic field file");
0151
0152     f.close();
0153   }
0154
0155   TRACCC_INFO("Field dimensions in x = [" << minx << ", " << maxx << "]");
0156   TRACCC_INFO("Field dimensions in y = [" << miny << ", " << maxy << "]");
0157   TRACCC_INFO("Field dimensions in z = [" << minz << ", " << maxz << "]");
0158
0159   TRACCC_INFO("Assuming sample spacing of 100.0 in each dimension");
0160
0161   /*
0162    * Now that we have the limits of our field, compute the size in each
0163    * dimension.
0164    */
0165   std::size_t sx =
0166       static_cast<std::size_t>(std::lround((maxx - minx) / 100.0)) + 1;
0167   std::size_t sy =
0168       static_cast<std::size_t>(std::lround((maxy - miny) / 100.0)) + 1;
0169   std::size_t sz =
0170       static_cast<std::size_t>(std::lround((maxz - minz) / 100.0)) + 1;
0171
0172   TRACCC_INFO("Magnetic field size is " << sx << "x" << sy << "x" << sz);
0173
0174   TRACCC_INFO("Constructing matching vector field...");
0175
0176   covfie::algebra::affine<3> translation =
0177       covfie::algebra::affine<3>::translation(-minx, -miny, -minz);
0178   covfie::algebra::affine<3> scaling = covfie::algebra::affine<3>::scaling(
0179       static_cast<float>(sx - 1) / (maxx - minx),
0180       static_cast<float>(sy - 1) / (maxy - miny),
0181       static_cast<float>(sz - 1) / (maxz - minz));
0182
0183   field_t field(covfie::make_parameter_pack(
0184       field_t::backend_t::configuration_t(scaling * translation),
0185       field_t::backend_t::backend_t::configuration_t{},
0186       field_t::backend_t::backend_t::backend_t::configuration_t{sx, sy, sz}));
0187   field_t::view_t fv(field);
0188
0189   {
0190     TRACCC_INFO("Re-opening magnetic field to gather data");
0191
0192     f.open(fn);
0193
0194     if (!f.good()) {
0195       TRACCC_FATAL("Failed to open input file " << fn << "!");
0196       std::exit(1);
0197     }
0198
0199     std::string line;
0200
0201     TRACCC_INFO("Skipping the first line (header)");
0202
0203     std::getline(f, line);
0204
0205     float xp, yp, zp;
0206     float Bx, By, Bz;
0207
0208     std::size_t n_lines = 0;
0209
0210     TRACCC_INFO("Iterating over lines in the magnetic field file");
0211
0212     /*
0213      * Read every line, and update our current minima and maxima
0214      * appropriately.
0215      */
0216     while (std::getline(f, line)) {
0217       CHECK_IOSTATE(f);
0218
0219       std::string word;
0220       std::stringstream ss(line);
0221       CHECK_IOSTATE(ss);
0222
0223       std::getline(ss, word, delimiter);
0224       CHECK_IOSTATE(ss);
0225       xp = static_cast<float>(std::atof(word.c_str()));
0226       std::getline(ss, word, delimiter);
0227       CHECK_IOSTATE(ss);
0228       yp = static_cast<float>(std::atof(word.c_str()));
0229       std::getline(ss, word, delimiter);
0230       CHECK_IOSTATE(ss);
0231       zp = static_cast<float>(std::atof(word.c_str()));
0232       std::getline(ss, word, delimiter);
0233       CHECK_IOSTATE(ss);
0234       Bx = static_cast<float>(std::atof(word.c_str()));
0235       std::getline(ss, word, delimiter);
0236       CHECK_IOSTATE(ss);
0237       By = static_cast<float>(std::atof(word.c_str()));
0238       std::getline(ss, word);
0239       Bz = static_cast<float>(std::atof(word.c_str()));
0240
0241       field_t::view_t::output_t& p = fv.at(xp, yp, zp);
0242
0243       p[0] = Bx * traccc::unit<float>::T;
0244       p[1] = By * traccc::unit<float>::T;
0245       p[2] = Bz * traccc::unit<float>::T;
0246
0247       n_lines++;
0248     }
0249
0250     TRACCC_INFO("Read " << n_lines << " lines of magnetic field data");
0251
0252     TRACCC_INFO("Closing magnetic field file");
0253
0254     f.close();
0255   }
0256
0257   return field;
0258 }
0259
0260 int main(int argc, char** argv) {
0261   boost::program_options::variables_map vm;
0262   parse_opts(argc, argv, vm);
0263
0264   TRACCC_INFO("Welcome to the traccc magnetic field converter!");
0265   TRACCC_INFO("Using magnetic field file \"" << vm["input"].as<std::string>()
0266                                              << "\"");
0267   TRACCC_INFO("Starting read of input file...");
0268
0269   field_t fb =
0270       read_bfield(vm["input"].as<std::string>(), vm["delimiter"].as<char>());
0271
0272   TRACCC_INFO("Writing magnetic field to file \""
0273               << vm["output"].as<std::string>() << "\"...");
0274
0275   std::ofstream fs(vm["output"].as<std::string>(), std::ofstream::binary);
0276
0277   if (!fs.good()) {
0278     TRACCC_FATAL("Failed to open output file " << vm["output"].as<std::string>()
0279                                                << "!");
0280     std::exit(1);
0281   }
0282
0283   fb.dump(fs);
0284
0285   fs.close();
0286
0287   TRACCC_INFO("Conversion complete, goodbye!");
0288
0289   return 0;
0290 }