Root/src/RootBFieldWriter.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2021 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 http://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsExamples/Io/Root/RootBFieldWriter.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/VectorHelpers.hpp"
0016
0017 #include <cassert>
0018 #include <ios>
0019 #include <sstream>
0020 #include <stdexcept>
0021 #include <utility>
0022 #include <vector>
0023
0024 #include <TFile.h>
0025 #include <TTree.h>
0026
0027 namespace ActsExamples {
0028
0029 /// Write down an interpolated magnetic field map
0030 void RootBFieldWriter::run(const Config& config,
0031                            std::unique_ptr<const Acts::Logger> p_logger) {
0032   // Set up (local) logging
0033   // @todo Remove dangerous using declaration once the logger macro
0034   // tolerates it
0035   using namespace Acts;
0036   ACTS_LOCAL_LOGGER(std::move(p_logger))
0037
0038   Acts::MagneticFieldContext bFieldContext;
0039
0040   // Check basic configuration
0041   if (config.treeName.empty()) {
0042     throw std::invalid_argument("Missing tree name");
0043   } else if (config.fileName.empty()) {
0044     throw std::invalid_argument("Missing file name");
0045   } else if (config.bField == nullptr) {
0046     throw std::invalid_argument("Missing interpolated magnetic field");
0047   }
0048
0049   // Setup ROOT I/O
0050   ACTS_INFO("Registering new ROOT output File : " << config.fileName);
0051   TFile* outputFile =
0052       TFile::Open(config.fileName.c_str(), config.fileMode.c_str());
0053   if (outputFile == nullptr) {
0054     throw std::ios_base::failure("Could not open '" + config.fileName + "'");
0055   }
0056   TTree* outputTree = new TTree(config.treeName.c_str(),
0057                                 config.treeName.c_str(), 99, outputFile);
0058   if (outputTree == nullptr) {
0059     throw std::bad_alloc();
0060   }
0061
0062   // Access the minima and maxima of all axes
0063   auto minima = config.bField->getMin();
0064   auto maxima = config.bField->getMax();
0065   auto nBins = config.bField->getNBins();
0066
0067   if (logger().doPrint(Acts::Logging::VERBOSE)) {
0068     std::stringstream ss;
0069     ss << "Minimum:";
0070     for (auto m : minima) {
0071       ss << " " << m;
0072     }
0073     ACTS_VERBOSE(ss.str());
0074     ss.str("");
0075     ss << "Maximum:";
0076     for (auto m : maxima) {
0077       ss << " " << m;
0078     }
0079     ACTS_VERBOSE(ss.str());
0080     ss.str("");
0081     ss << "nBins:";
0082     for (auto m : nBins) {
0083       ss << " " << m;
0084     }
0085     ACTS_VERBOSE(ss.str());
0086   }
0087
0088   if (config.gridType == GridType::xyz) {
0089     ACTS_INFO("Map will be written out in cartesian coordinates (x,y,z).");
0090
0091     // Write out the interpolated magnetic field map
0092     double minX = 0., minY = 0., minZ = 0.;
0093     double maxX = 0., maxY = 0., maxZ = 0.;
0094     std::size_t nBinsX = 0, nBinsY = 0, nBinsZ = 0;
0095
0096     // The position values in xyz
0097     double x = 0;
0098     outputTree->Branch("x", &x);
0099     double y = 0;
0100     outputTree->Branch("y", &y);
0101     double z = 0;
0102     outputTree->Branch("z", &z);
0103
0104     // The BField values in xyz
0105     double Bx = 0;
0106     outputTree->Branch("Bx", &Bx);
0107     double By = 0;
0108     outputTree->Branch("By", &By);
0109     double Bz = 0;
0110     outputTree->Branch("Bz", &Bz);
0111
0112     // check if range is user defined
0113     if (config.rBounds && config.zBounds) {
0114       ACTS_INFO("User defined ranges handed over.");
0115
0116       // print out map in user defined range
0117       minX = config.rBounds->at(0);
0118       minY = config.rBounds->at(0);
0119       minZ = config.zBounds->at(0);
0120
0121       maxX = config.rBounds->at(1);
0122       maxY = config.rBounds->at(1);
0123       maxZ = config.zBounds->at(1);
0124
0125       nBinsX = config.rBins;
0126       nBinsY = config.rBins;
0127       nBinsZ = config.zBins;
0128
0129     } else {
0130       ACTS_INFO("No user defined ranges handed over - write out whole map.");
0131       // write out whole map
0132
0133       // check dimension of Bfieldmap
0134       if (minima.size() == 3 && maxima.size() == 3) {
0135         minX = minima.at(0);
0136         minY = minima.at(1);
0137         minZ = minima.at(2);
0138
0139         maxX = maxima.at(0);
0140         maxY = maxima.at(1);
0141         maxZ = maxima.at(2);
0142
0143         nBinsX = nBins.at(0);
0144         nBinsY = nBins.at(1);
0145         nBinsZ = nBins.at(2);
0146
0147       } else if (minima.size() == 2 && maxima.size() == 2) {
0148         minX = -maxima.at(0);
0149         minY = -maxima.at(0);
0150         minZ = minima.at(1);
0151
0152         maxX = maxima.at(0);
0153         maxY = maxima.at(0);
0154         maxZ = maxima.at(1);
0155
0156         nBinsX = nBins.at(0);
0157         nBinsY = nBins.at(0);
0158         nBinsZ = nBins.at(1);
0159       } else {
0160         throw std::invalid_argument(
0161             "BField has wrong dimension. The dimension needs to be "
0162             "either 2 (r,z,Br,Bz) or 3(x,y,z,Bx,By,Bz) in order to be "
0163             "written out by this writer.");
0164       }
0165     }
0166
0167     assert(maxX > minX);
0168     assert(maxY > minY);
0169     assert(maxZ > minZ);
0170
0171     double stepX = (maxX - minX) / (nBinsX - 1);
0172     double stepY = (maxY - minY) / (nBinsY - 1);
0173     double stepZ = (maxZ - minZ) / (nBinsZ - 1);
0174
0175     for (std::size_t i = 0; i < nBinsX; i++) {
0176       double raw_x = minX + i * stepX;
0177       for (std::size_t j = 0; j < nBinsY; j++) {
0178         double raw_y = minY + j * stepY;
0179         for (std::size_t k = 0; k < nBinsZ; k++) {
0180           double raw_z = minZ + k * stepZ;
0181           Acts::Vector3 position(raw_x, raw_y, raw_z);
0182           Vector3 bField = config.bField->getFieldUnchecked(position);
0183
0184           x = raw_x / Acts::UnitConstants::mm;
0185           y = raw_y / Acts::UnitConstants::mm;
0186           z = raw_z / Acts::UnitConstants::mm;
0187           Bx = bField.x() / Acts::UnitConstants::T;
0188           By = bField.y() / Acts::UnitConstants::T;
0189           Bz = bField.z() / Acts::UnitConstants::T;
0190           outputTree->Fill();
0191         }  // for z
0192       }    // for y
0193     }      // for x
0194
0195   } else {
0196     ACTS_INFO("Map will be written out in cylinder coordinates (r,z).");
0197
0198     // The position value in rz
0199     double r = 0;
0200     outputTree->Branch("r", &r);
0201     double z = 0;
0202     outputTree->Branch("z", &z);
0203     // The BField value in rz
0204     double Br = 0;
0205     outputTree->Branch("Br", &Br);
0206     double Bz = 0;
0207     outputTree->Branch("Bz", &Bz);
0208
0209     double minR = 0, maxR = 0;
0210     double minZ = 0, maxZ = 0;
0211     std::size_t nBinsR = 0, nBinsZ = 0;
0212
0213     if (config.rBounds && config.zBounds) {
0214       ACTS_INFO("User defined ranges handed over.");
0215
0216       minR = config.rBounds->at(0);
0217       minZ = config.zBounds->at(0);
0218
0219       maxR = config.rBounds->at(1);
0220       maxZ = config.zBounds->at(1);
0221
0222       nBinsR = config.rBins;
0223       nBinsZ = config.zBins;
0224     } else {
0225       ACTS_INFO("No user defined ranges handed over - printing out whole map.");
0226
0227       if (minima.size() == 3 && maxima.size() == 3) {
0228         minR = 0.;
0229         minZ = minima.at(2);
0230
0231         maxR = maxima.at(0);
0232         maxZ = maxima.at(2);
0233
0234         nBinsR = nBins.at(0);
0235         nBinsZ = nBins.at(2);
0236
0237       } else if (minima.size() == 2 || maxima.size() == 2) {
0238         minR = minima.at(0);
0239         minZ = minima.at(1);
0240
0241         maxR = maxima.at(0);
0242         maxZ = maxima.at(1);
0243
0244         nBinsR = nBins.at(0);
0245         nBinsZ = nBins.at(1);
0246
0247       } else {
0248         throw std::invalid_argument(
0249             "BField has wrong dimension. The dimension needs to be "
0250             "either 2 (r,z,Br,Bz) or 3(x,y,z,Bx,By,Bz) in order to be "
0251             "written out by this writer.");
0252       }
0253     }
0254
0255     assert(maxR > minR);
0256     assert(maxZ > minZ);
0257
0258     double stepR = (maxR - minR) / (nBinsR - 1);
0259     double stepZ = (maxZ - minZ) / (nBinsZ - 1);
0260
0261     for (std::size_t k = 0; k < nBinsZ; k++) {
0262       double raw_z = minZ + k * stepZ;
0263       for (std::size_t j = 0; j < nBinsR; j++) {
0264         double raw_r = minR + j * stepR;
0265         Acts::Vector3 position(raw_r, 0.0, raw_z);  // position at phi=0
0266         ACTS_VERBOSE("Requesting position: " << position.transpose());
0267         auto bField = config.bField->getFieldUnchecked(position);
0268         z = raw_z / Acts::UnitConstants::mm;
0269         r = raw_r / Acts::UnitConstants::mm;
0270         Bz = bField.z() / Acts::UnitConstants::T;
0271         Br = VectorHelpers::perp(bField) / Acts::UnitConstants::T;
0272         outputTree->Fill();
0273       }  // for R
0274     }    // for z
0275   }
0276
0277   // Tear down ROOT I/O
0278   ACTS_INFO("Closing and Writing ROOT output File : " << config.fileName);
0279   outputTree->Write();
0280   delete outputFile;
0281 }
0282 }  // namespace ActsExamples