ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RootBFieldWriter.hpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RootBFieldWriter.hpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2017 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
9 #pragma once
10 
15 #include <Acts/Utilities/Units.hpp>
16 #include <TFile.h>
17 #include <TTree.h>
18 #include <array>
19 #include <boost/optional.hpp>
20 #include <ios>
21 #include <mutex>
22 #include <sstream>
23 #include <stdexcept>
24 
28 
29 namespace FW {
30 
35 template <typename bfield_t>
37  public:
39  enum class GridType { rz = 0, xyz = 1 };
40 
41  struct Config {
43  std::string treeName = "TTree";
45  std::string fileName = "TFile.root";
47  std::string fileMode = "recreate";
49  std::shared_ptr<const bfield_t> bField = nullptr;
51  GridType gridType = GridType::xyz;
56  boost::optional<std::array<double, 2>> rBounds;
60  boost::optional<std::array<double, 2>> zBounds;
64  size_t rBins = 200;
66  // @note setting this parameter is optional, in case no bin numbers are
68  size_t zBins = 300;
70  // @note setting this parameter is optional, in case no bin numbers are
72  size_t phiBins = 100;
73  };
74 
76  static void run(const Config& cfg,
77  std::unique_ptr<const Acts::Logger> p_logger =
78  Acts::getDefaultLogger("RootBFieldWriter",
80  // Set up (local) logging
81  // @todo Remove dangerous using declaration once the logger macro
82  // tolerates it
83  using namespace Acts;
84  ACTS_LOCAL_LOGGER(std::move(p_logger))
85 
86  // Check basic configuration
87  if (cfg.treeName.empty()) {
88  throw std::invalid_argument("Missing tree name");
89  } else if (cfg.fileName.empty()) {
90  throw std::invalid_argument("Missing file name");
91  } else if (!cfg.bField) {
92  throw std::invalid_argument("Missing interpolated magnetic field");
93  }
94 
95  // Setup ROOT I/O
96  ACTS_INFO("Registering new ROOT output File : " << cfg.fileName);
97  TFile* outputFile = TFile::Open(cfg.fileName.c_str(), cfg.fileMode.c_str());
98  if (!outputFile) {
99  throw std::ios_base::failure("Could not open '" + cfg.fileName);
100  }
101  outputFile->cd();
102  TTree* outputTree = new TTree(cfg.treeName.c_str(), cfg.treeName.c_str());
103  if (!outputTree)
104  throw std::bad_alloc();
105 
106  // The position values
107  double z;
108  outputTree->Branch("z", &z);
109 
110  // The BField values
111  double Bz;
112  outputTree->Branch("Bz", &Bz);
113 
114  // Get the underlying mapper of the InterpolatedBFieldMap
115  auto mapper = cfg.bField->getMapper();
116 
117  // Access the minima and maxima of all axes
118  auto minima = mapper.getMin();
119  auto maxima = mapper.getMax();
120  auto nBins = mapper.getNBins();
121 
122  if (cfg.gridType == GridType::xyz) {
123  ACTS_INFO("Map will be written out in cartesian coordinates (x,y,z).");
124 
125  // Write out the interpolated magnetic field map
126  double stepX = 0., stepY = 0., stepZ = 0.;
127  double minX = 0., minY = 0., minZ = 0.;
128  double maxX = 0., maxY = 0., maxZ = 0.;
129  size_t nBinsX = 0, nBinsY = 0, nBinsZ = 0;
130 
131  // The position values in xy
132  double x;
133  outputTree->Branch("x", &x);
134  double y;
135  outputTree->Branch("y", &y);
136  // The BField values in xy
137  double Bx;
138  outputTree->Branch("Bx", &Bx);
139  double By;
140  outputTree->Branch("By", &By);
141 
142  // check if range is user defined
143  if (cfg.rBounds && cfg.zBounds) {
144  ACTS_INFO("User defined ranges handed over.");
145 
146  // print out map in user defined range
147  minX = cfg.rBounds->at(0);
148  minY = cfg.rBounds->at(0);
149  minZ = cfg.zBounds->at(0);
150 
151  maxX = cfg.rBounds->at(1);
152  maxY = cfg.rBounds->at(1);
153  maxZ = cfg.zBounds->at(1);
154 
155  nBinsX = cfg.rBins;
156  nBinsY = cfg.rBins;
157  nBinsZ = cfg.zBins;
158 
159  } else {
160  ACTS_INFO(
161  "No user defined ranges handed over - printing out whole map.");
162  // print out whole map
163  // check dimension of Bfieldmap
164  if (minima.size() == 3 && maxima.size() == 3) {
165  minX = minima.at(0);
166  minY = minima.at(1);
167  minZ = minima.at(2);
168 
169  maxX = maxima.at(0);
170  maxY = maxima.at(1);
171  maxZ = maxima.at(2);
172 
173  nBinsX = nBins.at(0);
174  nBinsY = nBins.at(1);
175  nBinsZ = nBins.at(2);
176 
177  } else if (minima.size() == 2 && maxima.size() == 2) {
178  minX = -maxima.at(0);
179  minY = -maxima.at(0);
180  minZ = minima.at(1);
181 
182  maxX = maxima.at(0);
183  maxY = maxima.at(0);
184  maxZ = maxima.at(1);
185 
186  nBinsX = nBins.at(0);
187  nBinsY = nBins.at(0);
188  nBinsZ = nBins.at(1);
189  } else {
190  std::ostringstream errorMsg;
191  errorMsg
192  << "BField has wrong dimension. The dimension needs to be "
193  "either 2 (r,z,Br,Bz) or 3(x,y,z,Bx,By,Bz) in order to be "
194  "written out by this writer.";
195  throw std::invalid_argument(errorMsg.str());
196  }
197  }
198 
199  stepX = fabs(minX - maxX) / nBinsX;
200  stepY = fabs(minY - maxY) / nBinsY;
201  stepZ = fabs(minZ - maxZ) / nBinsZ;
202 
203  for (size_t i = 0; i <= nBinsX; i++) {
204  double raw_x = minX + i * stepX;
205  for (size_t j = 0; j <= nBinsY; j++) {
206  double raw_y = minY + j * stepY;
207  for (size_t k = 0; k <= nBinsZ; k++) {
208  double raw_z = minZ + k * stepZ;
209  Acts::Vector3D position(raw_x, raw_y, raw_z);
210  if (cfg.bField->isInside(position)) {
211  auto bField = cfg.bField->getField(position);
212 
213  x = raw_x / Acts::units::_mm;
214  y = raw_y / Acts::units::_mm;
215  z = raw_z / Acts::units::_mm;
216  Bx = bField.x() / Acts::units::_T;
217  By = bField.y() / Acts::units::_T;
218  Bz = bField.z() / Acts::units::_T;
219  outputTree->Fill();
220  }
221  } // for z
222  } // for y
223  } // for x
224  } else {
225  ACTS_INFO("Map will be written out in cylinder coordinates (r,z).");
226  // The position value in r
227  double r;
228  outputTree->Branch("r", &r);
229  // The BField value in r
230  double Br;
231  outputTree->Branch("Br", &Br);
232 
233  double minR = 0, maxR = 0;
234  double minZ = 0, maxZ = 0;
235  size_t nBinsR = 0, nBinsZ = 0, nBinsPhi = 0;
236  double stepR = 0, stepZ = 0;
237 
238  if (cfg.rBounds && cfg.zBounds) {
239  ACTS_INFO("User defined ranges handed over.");
240 
241  minR = cfg.rBounds->at(0);
242  minZ = cfg.zBounds->at(0);
243 
244  maxR = cfg.rBounds->at(1);
245  maxZ = cfg.zBounds->at(1);
246 
247  nBinsR = cfg.rBins;
248  nBinsZ = cfg.zBins;
249  nBinsPhi = cfg.phiBins;
250  } else {
251  ACTS_INFO(
252  "No user defined ranges handed over - printing out whole map.");
253 
254  if (minima.size() == 3 && maxima.size() == 3) {
255  minR = 0.;
256  minZ = minima.at(2);
257 
258  maxR = maxima.at(0);
259  maxZ = maxima.at(2);
260 
261  nBinsR = nBins.at(0);
262  nBinsZ = nBins.at(2);
263  nBinsPhi = 100.;
264 
265  } else if (minima.size() == 2 || maxima.size() == 2) {
266  minR = minima.at(0);
267  minZ = minima.at(1);
268 
269  maxR = maxima.at(0);
270  maxZ = maxima.at(1);
271 
272  nBinsR = nBins.at(0);
273  nBinsZ = nBins.at(1);
274  nBinsPhi = 100.;
275 
276  } else {
277  std::ostringstream errorMsg;
278  errorMsg
279  << "BField has wrong dimension. The dimension needs to be "
280  "either 2 (r,z,Br,Bz) or 3(x,y,z,Bx,By,Bz) in order to be "
281  "written out by this writer.";
282  throw std::invalid_argument(errorMsg.str());
283  }
284  }
285  double minPhi = -M_PI;
286  stepR = fabs(minR - maxR) / nBinsR;
287  stepZ = fabs(minZ - maxZ) / nBinsZ;
288  double stepPhi = (2 * M_PI) / nBinsPhi;
289 
290  for (size_t i = 0; i < nBinsPhi; i++) {
291  double phi = minPhi + i * stepPhi;
292  for (size_t k = 0; k < nBinsZ; k++) {
293  double raw_z = minZ + k * stepZ;
294  for (size_t j = 0; j < nBinsR; j++) {
295  double raw_r = minR + j * stepR;
296  Acts::Vector3D position(raw_r * cos(phi), raw_r * sin(phi), raw_z);
297  if (cfg.bField->isInside(position)) {
298  auto bField = cfg.bField->getField(position);
299  z = raw_z / Acts::units::_mm;
300  r = raw_r / Acts::units::_mm;
301  Bz = bField.z() / Acts::units::_T;
303  outputTree->Fill();
304  }
305  }
306  }
307  } // for
308  }
309 
310  // Tear down ROOT I/O
311  ACTS_INFO("Closing and Writing ROOT output File : " << cfg.fileName);
312  outputFile->cd();
313  outputTree->Write();
314  outputFile->Close();
315  }
316 };
317 
318 } // namespace FW