Root/src/RootMaterialDecorator.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/RootMaterialDecorator.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Geometry/GeometryIdentifier.hpp"
0013 #include "Acts/Material/BinnedSurfaceMaterial.hpp"
0014 #include "Acts/Material/HomogeneousSurfaceMaterial.hpp"
0015 #include "Acts/Material/HomogeneousVolumeMaterial.hpp"
0016 #include "Acts/Material/InterpolatedMaterialMap.hpp"
0017 #include "Acts/Material/Material.hpp"
0018 #include "Acts/Material/MaterialGridHelper.hpp"
0019 #include "Acts/Material/MaterialSlab.hpp"
0020 #include "Acts/Utilities/AxisDefinitions.hpp"
0021 #include "Acts/Utilities/BinUtility.hpp"
0022 #include "Acts/Utilities/BinningType.hpp"
0023 #include "Acts/Utilities/Grid.hpp"
0024 #include "Acts/Utilities/Logger.hpp"
0025
0026 #include <iostream>
0027 #include <stdexcept>
0028 #include <string>
0029 #include <tuple>
0030 #include <vector>
0031
0032 #include <TFile.h>
0033 #include <TH1.h>
0034 #include <TH2.h>
0035 #include <TIterator.h>
0036 #include <TKey.h>
0037 #include <TList.h>
0038 #include <TObject.h>
0039 #include <boost/algorithm/string.hpp>
0040 #include <boost/algorithm/string/finder.hpp>
0041 #include <boost/algorithm/string/iter_find.hpp>
0042
0043 namespace Acts {
0044 class ISurfaceMaterial;
0045 class IVolumeMaterial;
0046 }  // namespace Acts
0047
0048 using namespace Acts;
0049
0050 ActsPlugins::RootMaterialDecorator::RootMaterialDecorator(
0051     const ActsPlugins::RootMaterialDecorator::Config& config,
0052     Logging::Level level)
0053     : m_cfg(config),
0054       m_logger{getDefaultLogger("RootMaterialDecorator", level)} {
0055   // Validate the configuration
0056   if (m_cfg.accessorOptions.folderSurfaceNameBase.empty()) {
0057     throw std::invalid_argument("Missing surface folder name base");
0058   } else if (m_cfg.accessorOptions.folderVolumeNameBase.empty()) {
0059     throw std::invalid_argument("Missing volume folder name base");
0060   } else if (m_cfg.fileName.empty()) {
0061     throw std::invalid_argument("Missing file name");
0062   }
0063
0064   // Set up ROOT I/O
0065   m_inputFile = TFile::Open(m_cfg.fileName.c_str());
0066   if (m_inputFile == nullptr) {
0067     throw std::ios_base::failure("Could not open '" + m_cfg.fileName + "'");
0068   }
0069
0070   ActsPlugins::RootMaterialMapIo accessor(m_cfg.accessorConfig,
0071                                           m_logger->clone("RootMaterialMapIo"));
0072   auto [surfaceMaps, volumeMaps] =
0073       accessor.read(*m_inputFile, m_cfg.accessorOptions);
0074
0075   m_surfaceMaterialMap = std::move(surfaceMaps);
0076
0077   // Get the list of keys from the file
0078   TList* tlist = m_inputFile->GetListOfKeys();
0079   auto tIter = tlist->MakeIterator();
0080   tIter->Reset();
0081
0082   // Iterate over the keys in the file
0083   while (auto* key = static_cast<TKey*>(tIter->Next())) {
0084     // Remember the directory
0085     std::string tdName(key->GetName());
0086
0087     std::vector<std::string> splitNames;
0088     iter_split(
0089         splitNames, tdName,
0090         boost::algorithm::first_finder(m_cfg.accessorConfig.volumePrefix));
0091
0092     ACTS_VERBOSE("Processing directory: " << tdName);
0093     if (splitNames[0] == m_cfg.accessorOptions.folderVolumeNameBase) {
0094       // The volume material to be read in for this
0095       std::shared_ptr<const IVolumeMaterial> vMaterial = nullptr;
0096       // Volume key
0097       boost::split(splitNames, splitNames[1], boost::is_any_of("_"));
0098       GeometryIdentifier::Value volID = std::stoi(splitNames[0]);
0099
0100       // Reconstruct the geometry ID
0101       auto geoID = GeometryIdentifier().withVolume(volID);
0102       ACTS_VERBOSE("GeometryIdentifier re-constructed as " << geoID);
0103
0104       // Construct the names
0105       std::string nName = tdName + "/" + m_cfg.accessorConfig.nBinsHistName;
0106       std::string vName = tdName + "/" + m_cfg.accessorConfig.axisDirHistName;
0107       std::string oName =
0108           tdName + "/" + m_cfg.accessorConfig.axisBoundaryTypeHistName;
0109       std::string minName =
0110           tdName + "/" + m_cfg.accessorConfig.minRangeHistName;
0111       std::string maxName =
0112           tdName + "/" + m_cfg.accessorConfig.maxRangeHistName;
0113       std::string x0Name = tdName + "/" + m_cfg.accessorConfig.x0HistName;
0114       std::string l0Name = tdName + "/" + m_cfg.accessorConfig.l0HistName;
0115       std::string aName = tdName + "/" + m_cfg.accessorConfig.aHistName;
0116       std::string zName = tdName + "/" + m_cfg.accessorConfig.zHistName;
0117       std::string rhoName = tdName + "/" + m_cfg.accessorConfig.rhoHistName;
0118
0119       // Get the histograms
0120       TH1F* n = dynamic_cast<TH1F*>(m_inputFile->Get(nName.c_str()));
0121       TH1F* v = dynamic_cast<TH1F*>(m_inputFile->Get(vName.c_str()));
0122       TH1F* o = dynamic_cast<TH1F*>(m_inputFile->Get(oName.c_str()));
0123       TH1F* min = dynamic_cast<TH1F*>(m_inputFile->Get(minName.c_str()));
0124       TH1F* max = dynamic_cast<TH1F*>(m_inputFile->Get(maxName.c_str()));
0125       TH1F* x0 = dynamic_cast<TH1F*>(m_inputFile->Get(x0Name.c_str()));
0126       TH1F* l0 = dynamic_cast<TH1F*>(m_inputFile->Get(l0Name.c_str()));
0127       TH1F* A = dynamic_cast<TH1F*>(m_inputFile->Get(aName.c_str()));
0128       TH1F* Z = dynamic_cast<TH1F*>(m_inputFile->Get(zName.c_str()));
0129       TH1F* rho = dynamic_cast<TH1F*>(m_inputFile->Get(rhoName.c_str()));
0130
0131       // Only go on when you have all the material histograms
0132       if ((x0 != nullptr) && (l0 != nullptr) && (A != nullptr) &&
0133           (Z != nullptr) && (rho != nullptr)) {
0134         // Get the number of grid points
0135         int points = x0->GetNbinsX();
0136         // If the bin information histograms are present the material is
0137         // either a 2D or a 3D grid
0138         if ((n != nullptr) && (v != nullptr) && (o != nullptr) &&
0139             (min != nullptr) && (max != nullptr)) {
0140           // Dimension of the grid
0141           int dim = n->GetNbinsX();
0142           // Now reconstruct the bin utilities
0143           BinUtility bUtility;
0144           for (int ib = 1; ib < dim + 1; ++ib) {
0145             auto nbins = static_cast<std::size_t>(n->GetBinContent(ib));
0146             auto val = static_cast<AxisDirection>(v->GetBinContent(ib));
0147             auto opt = static_cast<BinningOption>(o->GetBinContent(ib));
0148             float rmin = min->GetBinContent(ib);
0149             float rmax = max->GetBinContent(ib);
0150             bUtility += BinUtility(nbins, rmin, rmax, opt, val);
0151           }
0152           ACTS_VERBOSE("Created " << bUtility);
0153
0154           if (dim == 2) {
0155             // 2D Grid material
0156             std::function<Vector2(Vector3)> transfoGlobalToLocal;
0157             Grid2D grid = createGrid2D(bUtility, transfoGlobalToLocal);
0158
0159             Grid2D::point_t gMin = grid.minPosition();
0160             Grid2D::point_t gMax = grid.maxPosition();
0161             Grid2D::index_t gNBins = grid.numLocalBins();
0162
0163             EAxis axis1(gMin[0], gMax[0], gNBins[0]);
0164             EAxis axis2(gMin[1], gMax[1], gNBins[1]);
0165
0166             // Build the grid and fill it with data
0167             MaterialGrid2D mGrid(std::make_tuple(axis1, axis2));
0168
0169             for (int p = 1; p <= points; p++) {
0170               float dx0 = x0->GetBinContent(p);
0171               float dl0 = l0->GetBinContent(p);
0172               float da = A->GetBinContent(p);
0173               float dz = Z->GetBinContent(p);
0174               float drho = rho->GetBinContent(p);
0175               // Create material properties
0176               const auto material =
0177                   Material::fromMassDensity(dx0, dl0, da, dz, drho);
0178               mGrid.at(p - 1) = material.parameters();
0179             }
0180             MaterialMapLookup<MaterialGrid2D> matMap(transfoGlobalToLocal,
0181                                                      mGrid);
0182             vMaterial = std::make_shared<
0183                 InterpolatedMaterialMap<MaterialMapLookup<MaterialGrid2D>>>(
0184                 std::move(matMap), bUtility);
0185           } else if (dim == 3) {
0186             // 3D Grid material
0187             std::function<Vector3(Vector3)> transfoGlobalToLocal;
0188             Grid3D grid = createGrid3D(bUtility, transfoGlobalToLocal);
0189
0190             Grid3D::point_t gMin = grid.minPosition();
0191             Grid3D::point_t gMax = grid.maxPosition();
0192             Grid3D::index_t gNBins = grid.numLocalBins();
0193
0194             EAxis axis1(gMin[0], gMax[0], gNBins[0]);
0195             EAxis axis2(gMin[1], gMax[1], gNBins[1]);
0196             EAxis axis3(gMin[2], gMax[2], gNBins[2]);
0197
0198             // Build the grid and fill it with data
0199             MaterialGrid3D mGrid(std::make_tuple(axis1, axis2, axis3));
0200
0201             for (int p = 1; p <= points; p++) {
0202               float dx0 = x0->GetBinContent(p);
0203               float dl0 = l0->GetBinContent(p);
0204               float da = A->GetBinContent(p);
0205               float dz = Z->GetBinContent(p);
0206               float drho = rho->GetBinContent(p);
0207               // Create material properties
0208               const auto material =
0209                   Material::fromMassDensity(dx0, dl0, da, dz, drho);
0210               mGrid.at(p - 1) = material.parameters();
0211             }
0212             MaterialMapLookup<MaterialGrid3D> matMap(transfoGlobalToLocal,
0213                                                      mGrid);
0214             vMaterial = std::make_shared<
0215                 InterpolatedMaterialMap<MaterialMapLookup<MaterialGrid3D>>>(
0216                 std::move(matMap), bUtility);
0217           }
0218         } else {
0219           // Homogeneous material
0220           double dx0 = x0->GetBinContent(1);
0221           double dl0 = l0->GetBinContent(1);
0222           double da = A->GetBinContent(1);
0223           double dz = Z->GetBinContent(1);
0224           double drho = rho->GetBinContent(1);
0225           // Create material properties
0226           const auto material =
0227               Material::fromMassDensity(dx0, dl0, da, dz, drho);
0228           vMaterial = std::make_shared<HomogeneousVolumeMaterial>(material);
0229         }
0230       }
0231       ACTS_VERBOSE("Successfully read Material for : " << geoID);
0232
0233       // Insert into the new collection
0234       m_volumeMaterialMap.try_emplace(geoID, std::move(vMaterial));
0235     }
0236   }
0237 }
0238
0239 ActsPlugins::RootMaterialDecorator::~RootMaterialDecorator() {
0240   if (m_inputFile != nullptr) {
0241     m_inputFile->Close();
0242   }
0243 }