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