Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RootMaterialDecorator.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RootMaterialDecorator.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2017-2018 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 #include "ActsExamples/Io/Root/RootMaterialDecorator.hpp"
10 
11 #include "Acts/Definitions/Algebra.hpp"
12 #include "Acts/Material/InterpolatedMaterialMap.hpp"
13 #include "Acts/Material/Material.hpp"
14 #include "Acts/Material/MaterialGridHelper.hpp"
15 #include "Acts/Material/MaterialSlab.hpp"
16 #include "Acts/Utilities/Logger.hpp"
17 #include "Acts/Utilities/detail/Grid.hpp"
18 #include <Acts/Geometry/GeometryIdentifier.hpp>
19 #include <Acts/Material/BinnedSurfaceMaterial.hpp>
20 #include <Acts/Material/HomogeneousSurfaceMaterial.hpp>
21 #include <Acts/Material/HomogeneousVolumeMaterial.hpp>
22 #include <Acts/Utilities/BinUtility.hpp>
23 #include <Acts/Utilities/BinningType.hpp>
24 
25 #include <algorithm>
26 #include <cstdio>
27 #include <functional>
28 #include <iostream>
29 #include <stdexcept>
30 #include <string>
31 #include <tuple>
32 #include <vector>
33 
34 #include <TFile.h>
35 #include <TH1.h>
36 #include <TH2.h>
37 #include <TIterator.h>
38 #include <TKey.h>
39 #include <TList.h>
40 #include <TObject.h>
41 #include <boost/algorithm/string.hpp>
42 #include <boost/algorithm/string/finder.hpp>
43 #include <boost/algorithm/string/iter_find.hpp>
44 
45 namespace Acts {
46 class ISurfaceMaterial;
47 class IVolumeMaterial;
48 } // namespace Acts
49 
50 ActsExamples::RootMaterialDecorator::RootMaterialDecorator(
51  const ActsExamples::RootMaterialDecorator::Config& config,
52  Acts::Logging::Level level)
53  : m_cfg(config),
54  m_logger{Acts::getDefaultLogger("RootMaterialDecorator", level)} {
55  // Validate the configuration
56  if (m_cfg.folderSurfaceNameBase.empty()) {
57  throw std::invalid_argument("Missing surface folder name base");
58  } else if (m_cfg.folderVolumeNameBase.empty()) {
59  throw std::invalid_argument("Missing volume folder name base");
60  } else if (m_cfg.fileName.empty()) {
61  throw std::invalid_argument("Missing file name");
62  }
63 
64  // Setup ROOT I/O
65  m_inputFile = TFile::Open(m_cfg.fileName.c_str());
66  if (m_inputFile == nullptr) {
67  throw std::ios_base::failure("Could not open '" + m_cfg.fileName);
68  }
69 
70  // Get the list of keys from the file
71  TList* tlist = m_inputFile->GetListOfKeys();
72  auto tIter = tlist->MakeIterator();
73  tIter->Reset();
74 
75  // Iterate over the keys in the file
76  while (TKey* key = (TKey*)(tIter->Next())) {
77  // Remember the directory
78  std::string tdName(key->GetName());
79 
80  ACTS_VERBOSE("Processing directory: " << tdName);
81 
82  // volume
83  std::vector<std::string> splitNames;
84  iter_split(splitNames, tdName,
85  boost::algorithm::first_finder(m_cfg.voltag));
86  // Surface Material
87  if (splitNames[0] == m_cfg.folderSurfaceNameBase) {
88  // The surface material to be read in for this
89  std::shared_ptr<const Acts::ISurfaceMaterial> sMaterial = nullptr;
90 
91  boost::split(splitNames, splitNames[1], boost::is_any_of("_"));
92  Acts::GeometryIdentifier::Value volID = std::stoi(splitNames[0]);
93  // boundary
94  iter_split(splitNames, tdName,
95  boost::algorithm::first_finder(m_cfg.boutag));
96  boost::split(splitNames, splitNames[1], boost::is_any_of("_"));
97  Acts::GeometryIdentifier::Value bouID = std::stoi(splitNames[0]);
98  // layer
99  iter_split(splitNames, tdName,
100  boost::algorithm::first_finder(m_cfg.laytag));
101  boost::split(splitNames, splitNames[1], boost::is_any_of("_"));
102  Acts::GeometryIdentifier::Value layID = std::stoi(splitNames[0]);
103  // approach
104  iter_split(splitNames, tdName,
105  boost::algorithm::first_finder(m_cfg.apptag));
106  boost::split(splitNames, splitNames[1], boost::is_any_of("_"));
107  Acts::GeometryIdentifier::Value appID = std::stoi(splitNames[0]);
108  // sensitive
109  iter_split(splitNames, tdName,
110  boost::algorithm::first_finder(m_cfg.sentag));
111  Acts::GeometryIdentifier::Value senID = std::stoi(splitNames[1]);
112 
113  // Reconstruct the geometry ID
114  Acts::GeometryIdentifier geoID;
115  geoID.setVolume(volID);
116  geoID.setBoundary(bouID);
117  geoID.setLayer(layID);
118  geoID.setApproach(appID);
119  geoID.setSensitive(senID);
120  ACTS_VERBOSE("GeometryIdentifier re-constructed as " << geoID);
121 
122  // Construct the names
123  std::string nName = tdName + "/" + m_cfg.ntag;
124  std::string vName = tdName + "/" + m_cfg.vtag;
125  std::string oName = tdName + "/" + m_cfg.otag;
126  std::string minName = tdName + "/" + m_cfg.mintag;
127  std::string maxName = tdName + "/" + m_cfg.maxtag;
128  std::string tName = tdName + "/" + m_cfg.ttag;
129  std::string x0Name = tdName + "/" + m_cfg.x0tag;
130  std::string l0Name = tdName + "/" + m_cfg.l0tag;
131  std::string aName = tdName + "/" + m_cfg.atag;
132  std::string zName = tdName + "/" + m_cfg.ztag;
133  std::string rhoName = tdName + "/" + m_cfg.rhotag;
134 
135  // Get the histograms
136  TH1F* n = dynamic_cast<TH1F*>(m_inputFile->Get(nName.c_str()));
137  TH1F* v = dynamic_cast<TH1F*>(m_inputFile->Get(vName.c_str()));
138  TH1F* o = dynamic_cast<TH1F*>(m_inputFile->Get(oName.c_str()));
139  TH1F* min = dynamic_cast<TH1F*>(m_inputFile->Get(minName.c_str()));
140  TH1F* max = dynamic_cast<TH1F*>(m_inputFile->Get(maxName.c_str()));
141  TH2F* t = dynamic_cast<TH2F*>(m_inputFile->Get(tName.c_str()));
142  TH2F* x0 = dynamic_cast<TH2F*>(m_inputFile->Get(x0Name.c_str()));
143  TH2F* l0 = dynamic_cast<TH2F*>(m_inputFile->Get(l0Name.c_str()));
144  TH2F* A = dynamic_cast<TH2F*>(m_inputFile->Get(aName.c_str()));
145  TH2F* Z = dynamic_cast<TH2F*>(m_inputFile->Get(zName.c_str()));
146  TH2F* rho = dynamic_cast<TH2F*>(m_inputFile->Get(rhoName.c_str()));
147 
148  std::vector<const TH1*> hists{n, v, o, min, max, t, x0, l0, A, Z, rho};
149 
150  // Only go on when you have all histograms
151  if (std::all_of(hists.begin(), hists.end(),
152  [](const auto* hist) { return hist != nullptr; })) {
153  // Get the number of bins
154  int nbins0 = t->GetNbinsX();
155  int nbins1 = t->GetNbinsY();
156 
157  // The material matrix
158  Acts::MaterialSlabMatrix materialMatrix(
159  nbins1, Acts::MaterialSlabVector(nbins0, Acts::MaterialSlab()));
160 
161  // We need binned material properties
162  if (nbins0 * nbins1 > 1) {
163  // Fill the matrix first
164  for (int ib0 = 1; ib0 <= nbins0; ++ib0) {
165  for (int ib1 = 1; ib1 <= nbins1; ++ib1) {
166  double dt = t->GetBinContent(ib0, ib1);
167  if (dt > 0.) {
168  double dx0 = x0->GetBinContent(ib0, ib1);
169  double dl0 = l0->GetBinContent(ib0, ib1);
170  double da = A->GetBinContent(ib0, ib1);
171  double dz = Z->GetBinContent(ib0, ib1);
172  double drho = rho->GetBinContent(ib0, ib1);
173  // Create material properties
174  const auto material =
175  Acts::Material::fromMassDensity(dx0, dl0, da, dz, drho);
176  materialMatrix[ib1 - 1][ib0 - 1] =
177  Acts::MaterialSlab(material, dt);
178  }
179  }
180  }
181 
182  // Now reconstruct the bin untilities
183  Acts::BinUtility bUtility;
184  for (int ib = 1; ib < n->GetNbinsX() + 1; ++ib) {
185  size_t nbins = size_t(n->GetBinContent(ib));
186  Acts::BinningValue val = Acts::BinningValue(v->GetBinContent(ib));
187  Acts::BinningOption opt = Acts::BinningOption(o->GetBinContent(ib));
188  float rmin = min->GetBinContent(ib);
189  float rmax = max->GetBinContent(ib);
190  bUtility += Acts::BinUtility(nbins, rmin, rmax, opt, val);
191  }
192  ACTS_VERBOSE("Created " << bUtility);
193 
194  // Construct the binned material with the right bin utility
195  sMaterial = std::make_shared<const Acts::BinnedSurfaceMaterial>(
196  bUtility, std::move(materialMatrix));
197 
198  } else {
199  // Only homogeneous material present
200  double dt = t->GetBinContent(1, 1);
201  double dx0 = x0->GetBinContent(1, 1);
202  double dl0 = l0->GetBinContent(1, 1);
203  double da = A->GetBinContent(1, 1);
204  double dz = Z->GetBinContent(1, 1);
205  double drho = rho->GetBinContent(1, 1);
206  // Create and set the homogeneous surface material
207  const auto material =
208  Acts::Material::fromMassDensity(dx0, dl0, da, dz, drho);
209  sMaterial = std::make_shared<const Acts::HomogeneousSurfaceMaterial>(
210  Acts::MaterialSlab(material, dt));
211  }
212  }
213  ACTS_VERBOSE("Successfully read Material for : " << geoID);
214 
215  // Insert into the new collection
216  m_surfaceMaterialMap.insert({geoID, std::move(sMaterial)});
217 
218  } else if (splitNames[0] == m_cfg.folderVolumeNameBase) {
219  // The volume material to be read in for this
220  std::shared_ptr<const Acts::IVolumeMaterial> vMaterial = nullptr;
221  // Volume key
222  boost::split(splitNames, splitNames[1], boost::is_any_of("_"));
223  Acts::GeometryIdentifier::Value volID = std::stoi(splitNames[0]);
224 
225  // Reconstruct the geometry ID
226  Acts::GeometryIdentifier geoID;
227  geoID.setVolume(volID);
228  ACTS_VERBOSE("GeometryIdentifier re-constructed as " << geoID);
229 
230  // Construct the names
231  std::string nName = tdName + "/" + m_cfg.ntag;
232  std::string vName = tdName + "/" + m_cfg.vtag;
233  std::string oName = tdName + "/" + m_cfg.otag;
234  std::string minName = tdName + "/" + m_cfg.mintag;
235  std::string maxName = tdName + "/" + m_cfg.maxtag;
236  std::string x0Name = tdName + "/" + m_cfg.x0tag;
237  std::string l0Name = tdName + "/" + m_cfg.l0tag;
238  std::string aName = tdName + "/" + m_cfg.atag;
239  std::string zName = tdName + "/" + m_cfg.ztag;
240  std::string rhoName = tdName + "/" + m_cfg.rhotag;
241 
242  // Get the histograms
243  TH1F* n = dynamic_cast<TH1F*>(m_inputFile->Get(nName.c_str()));
244  TH1F* v = dynamic_cast<TH1F*>(m_inputFile->Get(vName.c_str()));
245  TH1F* o = dynamic_cast<TH1F*>(m_inputFile->Get(oName.c_str()));
246  TH1F* min = dynamic_cast<TH1F*>(m_inputFile->Get(minName.c_str()));
247  TH1F* max = dynamic_cast<TH1F*>(m_inputFile->Get(maxName.c_str()));
248  TH1F* x0 = dynamic_cast<TH1F*>(m_inputFile->Get(x0Name.c_str()));
249  TH1F* l0 = dynamic_cast<TH1F*>(m_inputFile->Get(l0Name.c_str()));
250  TH1F* A = dynamic_cast<TH1F*>(m_inputFile->Get(aName.c_str()));
251  TH1F* Z = dynamic_cast<TH1F*>(m_inputFile->Get(zName.c_str()));
252  TH1F* rho = dynamic_cast<TH1F*>(m_inputFile->Get(rhoName.c_str()));
253 
254  // Only go on when you have all the material histograms
255  if ((x0 != nullptr) and (l0 != nullptr) and (A != nullptr) and
256  (Z != nullptr) and (rho != nullptr)) {
257  // Get the number of grid points
258  int points = x0->GetNbinsX();
259  // If the bin information histograms are present the material is
260  // either a 2D or a 3D grid
261  if ((n != nullptr) and (v != nullptr) and (o != nullptr) and
262  (min != nullptr) and (max != nullptr)) {
263  // Dimension of the grid
264  int dim = n->GetNbinsX();
265  // Now reconstruct the bin untilities
266  Acts::BinUtility bUtility;
267  for (int ib = 1; ib < dim + 1; ++ib) {
268  size_t nbins = size_t(n->GetBinContent(ib));
269  Acts::BinningValue val = Acts::BinningValue(v->GetBinContent(ib));
270  Acts::BinningOption opt = Acts::BinningOption(o->GetBinContent(ib));
271  float rmin = min->GetBinContent(ib);
272  float rmax = max->GetBinContent(ib);
273  bUtility += Acts::BinUtility(nbins, rmin, rmax, opt, val);
274  }
275  ACTS_VERBOSE("Created " << bUtility);
276 
277  if (dim == 2) {
278  // 2D Grid material
279  std::function<Acts::Vector2(Acts::Vector3)> transfoGlobalToLocal;
280  Acts::Grid2D grid = createGrid2D(bUtility, transfoGlobalToLocal);
281 
282  Acts::Grid2D::point_t gMin = grid.minPosition();
283  Acts::Grid2D::point_t gMax = grid.maxPosition();
284  Acts::Grid2D::index_t gNBins = grid.numLocalBins();
285 
286  Acts::EAxis axis1(gMin[0], gMax[0], gNBins[0]);
287  Acts::EAxis axis2(gMin[1], gMax[1], gNBins[1]);
288 
289  // Build the grid and fill it with data
290  Acts::MaterialGrid2D mGrid(std::make_tuple(axis1, axis2));
291 
292  for (int p = 1; p <= points; p++) {
293  double dx0 = x0->GetBinContent(p);
294  double dl0 = l0->GetBinContent(p);
295  double da = A->GetBinContent(p);
296  double dz = Z->GetBinContent(p);
297  double drho = rho->GetBinContent(p);
298  // Create material properties
299  const auto material =
300  Acts::Material::fromMassDensity(dx0, dl0, da, dz, drho);
301  mGrid.at(p - 1) = material.parameters();
302  }
303  Acts::MaterialMapper<Acts::MaterialGrid2D> matMap(
304  transfoGlobalToLocal, mGrid);
305  vMaterial = std::make_shared<Acts::InterpolatedMaterialMap<
306  Acts::MaterialMapper<Acts::MaterialGrid2D>>>(std::move(matMap),
307  bUtility);
308  } else if (dim == 3) {
309  // 3D Grid material
310  std::function<Acts::Vector3(Acts::Vector3)> transfoGlobalToLocal;
311  Acts::Grid3D grid = createGrid3D(bUtility, transfoGlobalToLocal);
312 
313  Acts::Grid3D::point_t gMin = grid.minPosition();
314  Acts::Grid3D::point_t gMax = grid.maxPosition();
315  Acts::Grid3D::index_t gNBins = grid.numLocalBins();
316 
317  Acts::EAxis axis1(gMin[0], gMax[0], gNBins[0]);
318  Acts::EAxis axis2(gMin[1], gMax[1], gNBins[1]);
319  Acts::EAxis axis3(gMin[2], gMax[2], gNBins[2]);
320 
321  // Build the grid and fill it with data
322  Acts::MaterialGrid3D mGrid(std::make_tuple(axis1, axis2, axis3));
323 
324  for (int p = 1; p <= points; p++) {
325  double dx0 = x0->GetBinContent(p);
326  double dl0 = l0->GetBinContent(p);
327  double da = A->GetBinContent(p);
328  double dz = Z->GetBinContent(p);
329  double drho = rho->GetBinContent(p);
330  // Create material properties
331  const auto material =
332  Acts::Material::fromMassDensity(dx0, dl0, da, dz, drho);
333  mGrid.at(p - 1) = material.parameters();
334  }
335  Acts::MaterialMapper<Acts::MaterialGrid3D> matMap(
336  transfoGlobalToLocal, mGrid);
337  vMaterial = std::make_shared<Acts::InterpolatedMaterialMap<
338  Acts::MaterialMapper<Acts::MaterialGrid3D>>>(std::move(matMap),
339  bUtility);
340  }
341  } else {
342  // Homogeneous material
343  double dx0 = x0->GetBinContent(1);
344  double dl0 = l0->GetBinContent(1);
345  double da = A->GetBinContent(1);
346  double dz = Z->GetBinContent(1);
347  double drho = rho->GetBinContent(1);
348  // Create material properties
349  const auto material =
350  Acts::Material::fromMassDensity(dx0, dl0, da, dz, drho);
351  vMaterial =
352  std::make_shared<Acts::HomogeneousVolumeMaterial>(material);
353  }
354  }
355  ACTS_VERBOSE("Successfully read Material for : " << geoID);
356 
357  // Insert into the new collection
358  m_volumeMaterialMap.insert({geoID, std::move(vMaterial)});
359 
360  // Incorrect FolderName value
361  } else {
362  ACTS_ERROR(
363  "Invalid FolderName, does not match any entry in the root file");
364  }
365  }
366 }
367 
368 ActsExamples::RootMaterialDecorator::~RootMaterialDecorator() {
369  if (m_inputFile != nullptr) {
370  m_inputFile->Close();
371  }
372 }