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