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