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