EIC 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/Geometry/GeometryIdentifier.hpp>
12 #include <Acts/Material/BinnedSurfaceMaterial.hpp>
13 #include <Acts/Material/HomogeneousSurfaceMaterial.hpp>
14 #include <Acts/Utilities/BinUtility.hpp>
15 #include <Acts/Utilities/BinningType.hpp>
16 
17 #include <cstdio>
18 #include <iostream>
19 #include <sstream>
20 #include <stdexcept>
21 #include <string>
22 
23 #include <TFile.h>
24 #include <TH2F.h>
25 #include <TIterator.h>
26 #include <TKey.h>
27 #include <TList.h>
28 #include <boost/algorithm/string.hpp>
29 #include <boost/algorithm/string/finder.hpp>
30 #include <boost/algorithm/string/iter_find.hpp>
31 
32 ActsExamples::RootMaterialDecorator::RootMaterialDecorator(
33  const ActsExamples::RootMaterialDecorator::Config& cfg)
34  : m_cfg(cfg), m_inputFile(nullptr) {
35  // Validate the configuration
36  if (m_cfg.folderNameBase.empty()) {
37  throw std::invalid_argument("Missing ROOT folder name");
38  } else if (m_cfg.fileName.empty()) {
39  throw std::invalid_argument("Missing file name");
40  } else if (!m_cfg.logger) {
41  throw std::invalid_argument("Missing logger");
42  } else if (m_cfg.name.empty()) {
43  throw std::invalid_argument("Missing service name");
44  }
45 
46  // Setup ROOT I/O
47  m_inputFile = TFile::Open(m_cfg.fileName.c_str());
48  if (!m_inputFile) {
49  throw std::ios_base::failure("Could not open '" + m_cfg.fileName);
50  }
51 
52  // Get the list of keys from the file
53  TList* tlist = m_inputFile->GetListOfKeys();
54  auto tIter = tlist->MakeIterator();
55  tIter->Reset();
56 
57  // Iterate over the keys in the file
58  while (TKey* key = (TKey*)(tIter->Next())) {
59  // The surface material to be read in for this
60  std::shared_ptr<const Acts::ISurfaceMaterial> sMaterial = nullptr;
61 
62  // Remember the directory
63  std::string tdName(key->GetName());
64 
65  ACTS_VERBOSE("Processing directory: " << tdName);
66 
67  // volume
68  std::vector<std::string> splitNames;
69  iter_split(splitNames, tdName,
70  boost::algorithm::first_finder(m_cfg.voltag));
71  boost::split(splitNames, splitNames[1], boost::is_any_of("_"));
72  Acts::GeometryIdentifier::Value volID = std::stoi(splitNames[0]);
73  // boundary
74  iter_split(splitNames, tdName,
75  boost::algorithm::first_finder(m_cfg.boutag));
76  boost::split(splitNames, splitNames[1], boost::is_any_of("_"));
77  Acts::GeometryIdentifier::Value bouID = std::stoi(splitNames[0]);
78  // layer
79  iter_split(splitNames, tdName,
80  boost::algorithm::first_finder(m_cfg.laytag));
81  boost::split(splitNames, splitNames[1], boost::is_any_of("_"));
82  Acts::GeometryIdentifier::Value layID = std::stoi(splitNames[0]);
83  // approach
84  iter_split(splitNames, tdName,
85  boost::algorithm::first_finder(m_cfg.apptag));
86  boost::split(splitNames, splitNames[1], boost::is_any_of("_"));
87  Acts::GeometryIdentifier::Value appID = std::stoi(splitNames[0]);
88  // sensitive
89  iter_split(splitNames, tdName,
90  boost::algorithm::first_finder(m_cfg.sentag));
91  Acts::GeometryIdentifier::Value senID = std::stoi(splitNames[1]);
92 
93  // Reconstruct the geometry ID
94  Acts::GeometryIdentifier geoID;
95  geoID.setVolume(volID);
96  geoID.setBoundary(bouID);
97  geoID.setLayer(layID);
98  geoID.setApproach(appID);
99  geoID.setSensitive(senID);
100  ACTS_VERBOSE("GeometryIdentifier re-constructed as " << geoID);
101 
102  // Construct the names
103  std::string nName = tdName + "/" + m_cfg.ntag;
104  std::string vName = tdName + "/" + m_cfg.vtag;
105  std::string oName = tdName + "/" + m_cfg.otag;
106  std::string minName = tdName + "/" + m_cfg.mintag;
107  std::string maxName = tdName + "/" + m_cfg.maxtag;
108  std::string tName = tdName + "/" + m_cfg.ttag;
109  std::string x0Name = tdName + "/" + m_cfg.x0tag;
110  std::string l0Name = tdName + "/" + m_cfg.l0tag;
111  std::string aName = tdName + "/" + m_cfg.atag;
112  std::string zName = tdName + "/" + m_cfg.ztag;
113  std::string rhoName = tdName + "/" + m_cfg.rhotag;
114 
115  // Get the histograms
116  TH1F* n = dynamic_cast<TH1F*>(m_inputFile->Get(nName.c_str()));
117  TH1F* v = dynamic_cast<TH1F*>(m_inputFile->Get(vName.c_str()));
118  TH1F* o = dynamic_cast<TH1F*>(m_inputFile->Get(oName.c_str()));
119  TH1F* min = dynamic_cast<TH1F*>(m_inputFile->Get(minName.c_str()));
120  TH1F* max = dynamic_cast<TH1F*>(m_inputFile->Get(maxName.c_str()));
121  TH2F* t = dynamic_cast<TH2F*>(m_inputFile->Get(tName.c_str()));
122  TH2F* x0 = dynamic_cast<TH2F*>(m_inputFile->Get(x0Name.c_str()));
123  TH2F* l0 = dynamic_cast<TH2F*>(m_inputFile->Get(l0Name.c_str()));
124  TH2F* A = dynamic_cast<TH2F*>(m_inputFile->Get(aName.c_str()));
125  TH2F* Z = dynamic_cast<TH2F*>(m_inputFile->Get(zName.c_str()));
126  TH2F* rho = dynamic_cast<TH2F*>(m_inputFile->Get(rhoName.c_str()));
127 
128  // Only go on when you have all histograms
129  if (n and v and o and min and max and t and x0 and l0 and A and Z and rho) {
130  // Get the number of bins
131  int nbins0 = t->GetNbinsX();
132  int nbins1 = t->GetNbinsY();
133 
134  // The material matrix
135  Acts::MaterialSlabMatrix materialMatrix(
136  nbins1, Acts::MaterialSlabVector(nbins0, Acts::MaterialSlab()));
137 
138  // We need binned material properties
139  if (nbins0 * nbins1 > 1) {
140  // Fill the matrix first
141  for (int ib0 = 1; ib0 <= nbins0; ++ib0) {
142  for (int ib1 = 1; ib1 <= nbins1; ++ib1) {
143  double dt = t->GetBinContent(ib0, ib1);
144  if (dt > 0.) {
145  double dx0 = x0->GetBinContent(ib0, ib1);
146  double dl0 = l0->GetBinContent(ib0, ib1);
147  double da = A->GetBinContent(ib0, ib1);
148  double dz = Z->GetBinContent(ib0, ib1);
149  double drho = rho->GetBinContent(ib0, ib1);
150  // Create material properties
151  const auto material =
152  Acts::Material::fromMassDensity(dx0, dl0, da, dz, drho);
153  materialMatrix[ib1 - 1][ib0 - 1] =
154  Acts::MaterialSlab(material, dt);
155  }
156  }
157  }
158 
159  // Now reconstruct the bin untilities
160  Acts::BinUtility bUtility;
161  for (int ib = 1; ib < n->GetNbinsX() + 1; ++ib) {
162  size_t nbins = size_t(n->GetBinContent(ib));
163  Acts::BinningValue val = Acts::BinningValue(v->GetBinContent(ib));
164  Acts::BinningOption opt = Acts::BinningOption(o->GetBinContent(ib));
165  float rmin = min->GetBinContent(ib);
166  float rmax = max->GetBinContent(ib);
167  bUtility += Acts::BinUtility(nbins, rmin, rmax, opt, val);
168  }
169  ACTS_VERBOSE("Created " << bUtility);
170 
171  // Construct the binned material with the right bin utility
172  sMaterial = std::make_shared<const Acts::BinnedSurfaceMaterial>(
173  bUtility, std::move(materialMatrix));
174 
175  } else {
176  // Only homogeneous material present
177  double dt = t->GetBinContent(1, 1);
178  double dx0 = x0->GetBinContent(1, 1);
179  double dl0 = l0->GetBinContent(1, 1);
180  double da = A->GetBinContent(1, 1);
181  double dz = Z->GetBinContent(1, 1);
182  double drho = rho->GetBinContent(1, 1);
183  // Create and set the homogenous surface material
184  const auto material =
185  Acts::Material::fromMassDensity(dx0, dl0, da, dz, drho);
186  sMaterial = std::make_shared<const Acts::HomogeneousSurfaceMaterial>(
187  Acts::MaterialSlab(material, dt));
188  }
189  }
190  ACTS_VERBOSE("Successfully read Material for : " << geoID);
191 
192  // Insert into the new collection
193  m_surfaceMaterialMap.insert({geoID, std::move(sMaterial)});
194  }
195 }
196 
197 ActsExamples::RootMaterialDecorator::~RootMaterialDecorator() {
198  m_inputFile->Close();
199 }