MagneticField/src/FieldMapRootIo.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2017-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/MagneticField/FieldMapRootIo.hpp"
0010
0011 #include "Acts/MagneticField/BFieldMapUtils.hpp"
0012
0013 #include <stdexcept>
0014 #include <vector>
0015
0016 #include <RtypesCore.h>
0017 #include <TFile.h>
0018 #include <TTree.h>
0019
0020 ActsExamples::detail::InterpolatedMagneticField2
0021 ActsExamples::makeMagneticFieldMapRzFromRoot(
0022     const std::function<std::size_t(std::array<std::size_t, 2> binsRZ,
0023                                     std::array<std::size_t, 2> nBinsRZ)>&
0024         localToGlobalBin,
0025     const std::string& fieldMapFile, const std::string& treeName,
0026     Acts::ActsScalar lengthUnit, Acts::ActsScalar BFieldUnit,
0027     bool firstQuadrant) {
0028   /// [1] Read in field map file
0029   // Grid position points in r and z
0030   std::vector<double> rPos;
0031   std::vector<double> zPos;
0032   // components of magnetic field on grid points
0033   std::vector<Acts::Vector2> bField;
0034   // [1] Read in file and fill values
0035   TFile* inputFile = TFile::Open(fieldMapFile.c_str());
0036   if (inputFile == nullptr) {
0037     throw std::runtime_error("file does not exist");
0038   }
0039   TTree* tree = inputFile->Get<TTree>(treeName.c_str());
0040   if (tree == nullptr) {
0041     throw std::runtime_error("object not found in file");
0042   }
0043   Int_t entries = tree->GetEntries();
0044
0045   double r = 0, z = 0;
0046   double Br = 0, Bz = 0;
0047
0048   tree->SetBranchAddress("r", &r);
0049   tree->SetBranchAddress("z", &z);
0050
0051   tree->SetBranchAddress("Br", &Br);
0052   tree->SetBranchAddress("Bz", &Bz);
0053
0054   // reserve size
0055   rPos.reserve(entries);
0056   zPos.reserve(entries);
0057   bField.reserve(entries);
0058
0059   for (int i = 0; i < entries; i++) {
0060     tree->GetEvent(i);
0061     rPos.push_back(r);
0062     zPos.push_back(z);
0063     bField.push_back(Acts::Vector2(Br, Bz));
0064   }
0065   delete inputFile;
0066   /// [2] use helper function in core
0067   return Acts::fieldMapRZ(localToGlobalBin, rPos, zPos, bField, lengthUnit,
0068                           BFieldUnit, firstQuadrant);
0069 }
0070
0071 ActsExamples::detail::InterpolatedMagneticField3
0072 ActsExamples::makeMagneticFieldMapXyzFromRoot(
0073     const std::function<std::size_t(std::array<std::size_t, 3> binsXYZ,
0074                                     std::array<std::size_t, 3> nBinsXYZ)>&
0075         localToGlobalBin,
0076     const std::string& fieldMapFile, const std::string& treeName,
0077     Acts::ActsScalar lengthUnit, Acts::ActsScalar BFieldUnit,
0078     bool firstOctant) {
0079   /// [1] Read in field map file
0080   // Grid position points in x, y and z
0081   std::vector<double> xPos;
0082   std::vector<double> yPos;
0083   std::vector<double> zPos;
0084   // components of magnetic field on grid points
0085   std::vector<Acts::Vector3> bField;
0086   // [1] Read in file and fill values
0087   TFile* inputFile = TFile::Open(fieldMapFile.c_str());
0088   if (inputFile == nullptr) {
0089     throw std::runtime_error("file does not exist");
0090   }
0091   TTree* tree = inputFile->Get<TTree>(treeName.c_str());
0092   if (tree == nullptr) {
0093     throw std::runtime_error("object not found in file");
0094   }
0095   Int_t entries = tree->GetEntries();
0096
0097   float x = 0, y = 0, z = 0;
0098   float Bx = 0, By = 0, Bz = 0;
0099
0100   tree->SetBranchAddress("x", &x);
0101   tree->SetBranchAddress("y", &y);
0102   tree->SetBranchAddress("z", &z);
0103
0104   tree->SetBranchAddress("bx", &Bx);
0105   tree->SetBranchAddress("by", &By);
0106   tree->SetBranchAddress("bz", &Bz);
0107
0108   // reserve size
0109   xPos.reserve(entries);
0110   yPos.reserve(entries);
0111   zPos.reserve(entries);
0112   bField.reserve(entries);
0113
0114   for (int i = 0; i < entries; i++) {
0115     tree->GetEvent(i);
0116     xPos.push_back(x);
0117     yPos.push_back(y);
0118     zPos.push_back(z);
0119     bField.push_back(Acts::Vector3(Bx, By, Bz));
0120   }
0121   delete inputFile;
0122
0123   return Acts::fieldMapXYZ(localToGlobalBin, xPos, yPos, zPos, bField,
0124                            lengthUnit, BFieldUnit, firstOctant);
0125 }