Root/src/RootMaterialMapIo.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/RootMaterialMapIo.hpp"
0010
0011 #include "Acts/Geometry/GeometryIdentifier.hpp"
0012 #include "Acts/Material/BinnedSurfaceMaterial.hpp"
0013 #include "Acts/Material/GridSurfaceMaterial.hpp"
0014 #include "Acts/Material/HomogeneousSurfaceMaterial.hpp"
0015 #include "Acts/Material/Material.hpp"
0016 #include "Acts/Material/MaterialSlab.hpp"
0017 #include "Acts/Surfaces/Surface.hpp"
0018 #include "Acts/Utilities/Enumerate.hpp"
0019
0020 #include <TFile.h>
0021 #include <TH1I.h>
0022 #include <TH2F.h>
0023 #include <TKey.h>
0024 #include <TList.h>
0025 #include <TObject.h>
0026 #include <TTree.h>
0027 #include <boost/algorithm/string.hpp>
0028 #include <boost/algorithm/string/finder.hpp>
0029 #include <boost/algorithm/string/iter_find.hpp>
0030
0031 using namespace Acts;
0032
0033 void ActsPlugins::RootMaterialMapIo::write(
0034     TFile& rFile, const GeometryIdentifier& geoID,
0035     const ISurfaceMaterial& surfaceMaterial, const Options& options) {
0036   /// Change to the file
0037   rFile.cd();
0038
0039   // Homogeneous surface material writing into one tree
0040   auto homogeneousMaterial =
0041       dynamic_cast<const HomogeneousSurfaceMaterial*>(&surfaceMaterial);
0042   if (homogeneousMaterial != nullptr) {
0043     if (m_hTree == nullptr) {
0044       m_hTree = new TTree(options.homogeneousMaterialTreeName.c_str(),
0045                           "Homogeneous Material Tree");
0046       connectForWrite(*m_hTree, m_homogenousMaterialTreePayload);
0047     }
0048     fillMaterialSlab(m_homogenousMaterialTreePayload,
0049                      homogeneousMaterial->materialSlab());
0050     m_homogenousMaterialTreePayload.hGeoId = geoID.value();
0051     m_hTree->Fill();
0052     return;
0053   }
0054
0055   // Binned surface material writing
0056   auto bsMaterial =
0057       dynamic_cast<const BinnedSurfaceMaterial*>(&surfaceMaterial);
0058   if (bsMaterial != nullptr) {
0059     // decode the geometryID
0060     const auto gvolID = geoID.volume();
0061     const auto gbouID = geoID.boundary();
0062     const auto glayID = geoID.layer();
0063     const auto gappID = geoID.approach();
0064     const auto gsenID = geoID.sensitive();
0065     // create the directory
0066     std::string tdName = options.folderSurfaceNameBase.c_str();
0067     tdName += m_cfg.volumePrefix + std::to_string(gvolID);
0068     tdName += m_cfg.portalPrefix + std::to_string(gbouID);
0069     tdName += m_cfg.layerPrefix + std::to_string(glayID);
0070     tdName += m_cfg.passivePrefix + std::to_string(gappID);
0071     tdName += m_cfg.sensitivePrefix + std::to_string(gsenID);
0072     // create a new directory
0073     rFile.mkdir(tdName.c_str());
0074     rFile.cd(tdName.c_str());
0075
0076     // Boundary condistions
0077     // Get the binning data
0078     auto& binningData = bsMaterial->binUtility().binningData();
0079     // 1-D or 2-D maps
0080     auto bins = static_cast<int>(binningData.size());
0081     auto fBins = static_cast<float>(bins);
0082
0083     // The bin number information
0084     TH1F n(m_cfg.nBinsHistName.c_str(), "bins; bin", bins, -0.5, fBins - 0.5);
0085
0086     // The binning value information
0087     TH1F v(m_cfg.axisDirHistName.c_str(), "binning values; bin", bins, -0.5,
0088            fBins - 0.5);
0089
0090     // The binning option information
0091     TH1F o(m_cfg.axisBoundaryTypeHistName.c_str(), "binning options; bin", bins,
0092            -0.5, fBins - 0.5);
0093
0094     // The binning option information - range min
0095     TH1F rmin(m_cfg.minRangeHistName.c_str(), "min; bin", bins, -0.5,
0096               fBins - 0.5);
0097
0098     // The binning option information - range max
0099     TH1F rmax(m_cfg.maxRangeHistName.c_str(), "max; bin", bins, -0.5,
0100               fBins - 0.5);
0101
0102     // Now fill the histogram content
0103     for (auto [b, bData] : enumerate(binningData)) {
0104       // Fill: nbins, value, option, min, max
0105       n.SetBinContent(static_cast<int>(b) + 1, static_cast<int>(bData.bins()));
0106       v.SetBinContent(static_cast<int>(b) + 1,
0107                       static_cast<int>(bData.binvalue));
0108       o.SetBinContent(static_cast<int>(b) + 1, static_cast<int>(bData.option));
0109       rmin.SetBinContent(static_cast<int>(b) + 1, bData.min);
0110       rmax.SetBinContent(static_cast<int>(b) + 1, bData.max);
0111     }
0112     n.Write();
0113     v.Write();
0114     o.Write();
0115     rmin.Write();
0116     rmax.Write();
0117
0118     // If compressed writing is not enabled, write the binned surface material
0119     // as histograms
0120     if (!options.indexedMaterial) {
0121       fillBinnedSurfaceMaterial(*bsMaterial);
0122       return;
0123     }
0124
0125     // Otherwise, write the binned surface material into the TTree
0126     if (m_gTree == nullptr) {
0127       // Back to file level
0128       rFile.cd();
0129       m_gTree = new TTree(options.indexedMaterialTreeName.c_str(),
0130                           "Indexed Material Tree");
0131       connectForWrite(*m_gTree, m_indexedMaterialTreePayload);
0132       // Back to the directory
0133       rFile.cd(tdName.c_str());
0134     }
0135     fillBinnedSurfaceMaterial(m_indexedMaterialTreePayload, *bsMaterial);
0136     return;
0137   }
0138 }
0139
0140 void ActsPlugins::RootMaterialMapIo::write(
0141     TFile& rFile, const TrackingGeometryMaterial& detectorMaterial,
0142     const Options& options) {
0143   const auto& [surfaceMaterials, volumeMaterials] = detectorMaterial;
0144   for (const auto& [geoID, sMaterial] : surfaceMaterials) {
0145     write(rFile, geoID, *sMaterial, options);
0146   }
0147   if (m_hTree != nullptr) {
0148     m_hTree->Write();
0149   }
0150   if (m_gTree != nullptr) {
0151     m_gTree->Write();
0152   }
0153 }
0154
0155 void ActsPlugins::RootMaterialMapIo::connectForWrite(
0156     TTree& rTree, MaterialTreePayload& treePayload) {
0157   if (&treePayload == &m_homogenousMaterialTreePayload) {
0158     rTree.Branch("hGeoId", &treePayload.hGeoId);
0159   }
0160   rTree.Branch(m_cfg.thicknessHistName.c_str(), &treePayload.ht);
0161   rTree.Branch(m_cfg.x0HistName.c_str(), &treePayload.hX0);
0162   rTree.Branch(m_cfg.l0HistName.c_str(), &treePayload.hL0);
0163   rTree.Branch(m_cfg.aHistName.c_str(), &treePayload.hA);
0164   rTree.Branch(m_cfg.zHistName.c_str(), &treePayload.hZ);
0165   rTree.Branch(m_cfg.rhoHistName.c_str(), &treePayload.hRho);
0166 }
0167
0168 void ActsPlugins::RootMaterialMapIo::connectForRead(
0169     TTree& rTree, MaterialTreePayload& treePayload) {
0170   if (&treePayload == &m_homogenousMaterialTreePayload) {
0171     rTree.SetBranchAddress("hGeoId", &treePayload.hGeoId);
0172   }
0173   rTree.SetBranchAddress(m_cfg.thicknessHistName.c_str(), &treePayload.ht);
0174   rTree.SetBranchAddress(m_cfg.x0HistName.c_str(), &treePayload.hX0);
0175   rTree.SetBranchAddress(m_cfg.l0HistName.c_str(), &treePayload.hL0);
0176   rTree.SetBranchAddress(m_cfg.aHistName.c_str(), &treePayload.hA);
0177   rTree.SetBranchAddress(m_cfg.zHistName.c_str(), &treePayload.hZ);
0178   rTree.SetBranchAddress(m_cfg.rhoHistName.c_str(), &treePayload.hRho);
0179 }
0180
0181 void ActsPlugins::RootMaterialMapIo::fillMaterialSlab(
0182     MaterialTreePayload& payload, const MaterialSlab& materialSlab) {
0183   payload.ht = materialSlab.thickness();
0184   payload.hX0 = materialSlab.material().X0();
0185   payload.hL0 = materialSlab.material().L0();
0186   payload.hA = materialSlab.material().Ar();
0187   payload.hZ = materialSlab.material().Z();
0188   payload.hRho = materialSlab.material().massDensity();
0189 }
0190
0191 void ActsPlugins::RootMaterialMapIo::fillBinnedSurfaceMaterial(
0192     const BinnedSurfaceMaterial& bsMaterial) {
0193   auto bins0 = static_cast<int>(bsMaterial.binUtility().bins(0));
0194   auto bins1 = static_cast<int>(bsMaterial.binUtility().bins(1));
0195   auto fBins0 = static_cast<float>(bins0);
0196   auto fBins1 = static_cast<float>(bins1);
0197
0198   TH2F t(m_cfg.thicknessHistName.c_str(), "thickness [mm] ;b0 ;b1", bins0, -0.5,
0199          fBins0 - 0.5, bins1, -0.5, fBins1 - 0.5);
0200   TH2F x0(m_cfg.x0HistName.c_str(), "X_{0} [mm] ;b0 ;b1", bins0, -0.5,
0201           fBins0 - 0.5, bins1, -0.5, fBins1 - 0.5);
0202   TH2F l0(m_cfg.l0HistName.c_str(), "#Lambda_{0} [mm] ;b0 ;b1", bins0, -0.5,
0203           fBins0 - 0.5, bins1, -0.5, fBins1 - 0.5);
0204   TH2F A(m_cfg.aHistName.c_str(), "X_{0} [mm] ;b0 ;b1", bins0, -0.5,
0205          fBins0 - 0.5, bins1, -0.5, fBins1 - 0.5);
0206   TH2F Z(m_cfg.zHistName.c_str(), "#Lambda_{0} [mm] ;b0 ;b1", bins0, -0.5,
0207          fBins0 - 0.5, bins1, -0.5, fBins1 - 0.5);
0208   TH2F rho(m_cfg.rhoHistName.c_str(), "#rho [g/mm^3] ;b0 ;b1", bins0, -0.5,
0209            fBins0 - 0.5, bins1, -0.5, fBins1 - 0.5);
0210
0211   // Loop over the material matrix and fill the histograms
0212   const auto& materialMatrix = bsMaterial.fullMaterial();
0213   for (auto [b1, materialVector] : enumerate(materialMatrix)) {
0214     for (auto [b0, mat] : enumerate(materialVector)) {
0215       t.SetBinContent(static_cast<int>(b0) + 1, static_cast<int>(b1) + 1,
0216                       mat.thickness());
0217       x0.SetBinContent(static_cast<int>(b0) + 1, static_cast<int>(b1) + 1,
0218                        mat.material().X0());
0219       l0.SetBinContent(static_cast<int>(b0) + 1, static_cast<int>(b1) + 1,
0220                        mat.material().L0());
0221       A.SetBinContent(static_cast<int>(b0) + 1, static_cast<int>(b1) + 1,
0222                       mat.material().Ar());
0223       Z.SetBinContent(static_cast<int>(b0) + 1, static_cast<int>(b1) + 1,
0224                       mat.material().Z());
0225       rho.SetBinContent(static_cast<int>(b0) + 1, static_cast<int>(b1) + 1,
0226                         mat.material().massDensity());
0227     }
0228   }
0229   t.Write();
0230   x0.Write();
0231   l0.Write();
0232   A.Write();
0233   Z.Write();
0234   rho.Write();
0235 }
0236
0237 void ActsPlugins::RootMaterialMapIo::fillBinnedSurfaceMaterial(
0238     MaterialTreePayload& payload, const BinnedSurfaceMaterial& bsMaterial) {
0239   std::size_t bins0 = bsMaterial.binUtility().bins(0);
0240   std::size_t bins1 = bsMaterial.binUtility().bins(1);
0241
0242   TH2I idx(m_cfg.indexHistName.c_str(), "indices; bin0; bin1",
0243            static_cast<int>(bins0), -0.5, static_cast<float>(bins0) - 0.5,
0244            static_cast<int>(bins1), -0.5, static_cast<float>(bins1) - 0.5);
0245   // lLop over the material matrix, record the index and fill the indexed tree
0246   const auto& materialMatrix = bsMaterial.fullMaterial();
0247   for (auto [b1, materialVector] : enumerate(materialMatrix)) {
0248     for (auto [b0, mat] : enumerate(materialVector)) {
0249       idx.SetBinContent(static_cast<int>(b0) + 1, static_cast<int>(b1) + 1,
0250                         static_cast<float>(payload.index));
0251       payload.index++;
0252       fillMaterialSlab(payload, mat);
0253       m_gTree->Fill();
0254     }
0255   }
0256   idx.Write();
0257 }
0258
0259 TrackingGeometryMaterial ActsPlugins::RootMaterialMapIo::read(
0260     TFile& rFile, const Options& options) {
0261   TrackingGeometryMaterial detectorMaterial;
0262
0263   auto& [surfaceMaterials, volumeMaterials] = detectorMaterial;
0264
0265   auto homogeneousMaterialTree = dynamic_cast<TTree*>(
0266       rFile.Get(options.homogeneousMaterialTreeName.c_str()));
0267
0268   // Read homogeneous material tree
0269   if (homogeneousMaterialTree != nullptr) {
0270     connectForRead(*homogeneousMaterialTree, m_homogenousMaterialTreePayload);
0271     for (int i = 0; i < homogeneousMaterialTree->GetEntries(); ++i) {
0272       homogeneousMaterialTree->GetEntry(i);
0273       GeometryIdentifier geoID(m_homogenousMaterialTreePayload.hGeoId);
0274       MaterialSlab materialSlab(
0275           Material::fromMassDensity(m_homogenousMaterialTreePayload.hX0,
0276                                     m_homogenousMaterialTreePayload.hL0,
0277                                     m_homogenousMaterialTreePayload.hA,
0278                                     m_homogenousMaterialTreePayload.hZ,
0279                                     m_homogenousMaterialTreePayload.hRho),
0280           m_homogenousMaterialTreePayload.ht);
0281       auto homogeneousMaterial =
0282           std::make_shared<HomogeneousSurfaceMaterial>(materialSlab);
0283       surfaceMaterials.try_emplace(geoID, homogeneousMaterial);
0284     }
0285   }
0286
0287   // Read the binned surface material, if there - connect it to the payload
0288   auto indexedMaterialTree =
0289       dynamic_cast<TTree*>(rFile.Get(options.indexedMaterialTreeName.c_str()));
0290   if (indexedMaterialTree != nullptr) {
0291     connectForRead(*indexedMaterialTree, m_indexedMaterialTreePayload);
0292   }
0293
0294   // Get the list of keys from the file
0295   TList* tlist = rFile.GetListOfKeys();
0296   auto tIter = tlist->MakeIterator();
0297   tIter->Reset();
0298
0299   // Iterate over the keys in the file
0300   while (auto key = static_cast<TKey*>(tIter->Next())) {
0301     // Remember the directory
0302     std::string tdName(key->GetName());
0303
0304     ACTS_VERBOSE("Processing directory: " << tdName);
0305
0306     // volume
0307     std::vector<std::string> splitNames;
0308     iter_split(splitNames, tdName,
0309                boost::algorithm::first_finder(m_cfg.volumePrefix));
0310     // Surface Material
0311     if (splitNames[0] == options.folderSurfaceNameBase) {
0312       // The surface material to be read in for this
0313       std::shared_ptr<const ISurfaceMaterial> sMaterial = nullptr;
0314
0315       boost::split(splitNames, splitNames[1], boost::is_any_of("_"));
0316       GeometryIdentifier::Value volID = std::stoi(splitNames[0]);
0317       // boundary
0318       iter_split(splitNames, tdName,
0319                  boost::algorithm::first_finder(m_cfg.portalPrefix));
0320       boost::split(splitNames, splitNames[1], boost::is_any_of("_"));
0321       GeometryIdentifier::Value bouID = std::stoi(splitNames[0]);
0322       // layer
0323       iter_split(splitNames, tdName,
0324                  boost::algorithm::first_finder(m_cfg.layerPrefix));
0325       boost::split(splitNames, splitNames[1], boost::is_any_of("_"));
0326       GeometryIdentifier::Value layID = std::stoi(splitNames[0]);
0327       // approach
0328       iter_split(splitNames, tdName,
0329                  boost::algorithm::first_finder(m_cfg.passivePrefix));
0330       boost::split(splitNames, splitNames[1], boost::is_any_of("_"));
0331       GeometryIdentifier::Value appID = std::stoi(splitNames[0]);
0332       // sensitive
0333       iter_split(splitNames, tdName,
0334                  boost::algorithm::first_finder(m_cfg.sensitivePrefix));
0335       GeometryIdentifier::Value senID = std::stoi(splitNames[1]);
0336
0337       // Reconstruct the geometry ID
0338       auto geoID = GeometryIdentifier()
0339                        .withVolume(volID)
0340                        .withBoundary(bouID)
0341                        .withLayer(layID)
0342                        .withApproach(appID)
0343                        .withSensitive(senID);
0344
0345       ACTS_VERBOSE("GeometryIdentifier re-constructed as " << geoID);
0346
0347       auto texturedSurfaceMaterial =
0348           readTextureSurfaceMaterial(rFile, tdName, indexedMaterialTree);
0349       surfaceMaterials.try_emplace(geoID, texturedSurfaceMaterial);
0350     }
0351   }
0352   return detectorMaterial;
0353 }
0354
0355 std::shared_ptr<const ISurfaceMaterial>
0356 ActsPlugins::RootMaterialMapIo::readTextureSurfaceMaterial(
0357     TFile& rFile, const std::string& tdName, TTree* indexedMaterialTree) {
0358   std::shared_ptr<const ISurfaceMaterial> texturedSurfaceMaterial = nullptr;
0359
0360   // Construct the common names & get the common histograms
0361   std::string nName = tdName + "/" + m_cfg.nBinsHistName;
0362   std::string vName = tdName + "/" + m_cfg.axisDirHistName;
0363   std::string oName = tdName + "/" + m_cfg.axisBoundaryTypeHistName;
0364   std::string minName = tdName + "/" + m_cfg.minRangeHistName;
0365   std::string maxName = tdName + "/" + m_cfg.maxRangeHistName;
0366   // Get the histograms
0367   auto n = dynamic_cast<TH1F*>(rFile.Get(nName.c_str()));
0368   auto v = dynamic_cast<TH1F*>(rFile.Get(vName.c_str()));
0369   auto o = dynamic_cast<TH1F*>(rFile.Get(oName.c_str()));
0370   auto minh = dynamic_cast<TH1F*>(rFile.Get(minName.c_str()));
0371   auto maxh = dynamic_cast<TH1F*>(rFile.Get(maxName.c_str()));
0372
0373   std::vector<const TH1*> hists{n, v, o, minh, maxh};
0374   if (std::ranges::any_of(hists,
0375                           [](const auto* hist) { return hist == nullptr; })) {
0376     ACTS_ERROR(
0377         "Failed to read all required histograms for textured surface "
0378         "material from file: "
0379         << rFile.GetName());
0380     return nullptr;
0381   }
0382
0383   // Now reconstruct the bin utilities
0384   BinUtility bUtility;
0385   for (int ib = 1; ib < n->GetNbinsX() + 1; ++ib) {
0386     auto nbins = static_cast<std::size_t>(n->GetBinContent(ib));
0387     auto val = static_cast<AxisDirection>(v->GetBinContent(ib));
0388     auto opt = static_cast<BinningOption>(o->GetBinContent(ib));
0389     auto rmin = static_cast<float>(minh->GetBinContent(ib));
0390     auto rmax = static_cast<float>(maxh->GetBinContent(ib));
0391     bUtility += BinUtility(nbins, rmin, rmax, opt, val);
0392   }
0393   ACTS_VERBOSE("Created " << bUtility);
0394
0395   /// Draw from histogram only source
0396   if (indexedMaterialTree == nullptr) {
0397     // Construct the names for histogram type storage
0398     std::string tName = tdName + "/" + m_cfg.thicknessHistName;
0399     std::string x0Name = tdName + "/" + m_cfg.x0HistName;
0400     std::string l0Name = tdName + "/" + m_cfg.l0HistName;
0401     std::string aName = tdName + "/" + m_cfg.aHistName;
0402     std::string zName = tdName + "/" + m_cfg.zHistName;
0403     std::string rhoName = tdName + "/" + m_cfg.rhoHistName;
0404
0405     // Get the histograms
0406     auto t = dynamic_cast<TH2F*>(rFile.Get(tName.c_str()));
0407     auto x0 = dynamic_cast<TH2F*>(rFile.Get(x0Name.c_str()));
0408     auto l0 = dynamic_cast<TH2F*>(rFile.Get(l0Name.c_str()));
0409     auto A = dynamic_cast<TH2F*>(rFile.Get(aName.c_str()));
0410     auto Z = dynamic_cast<TH2F*>(rFile.Get(zName.c_str()));
0411     auto rho = dynamic_cast<TH2F*>(rFile.Get(rhoName.c_str()));
0412
0413     hists = {t, x0, l0, A, Z, rho};
0414
0415     // Only go on when you have all histograms
0416     if (std::ranges::all_of(hists,
0417                             [](const auto* hist) { return hist != nullptr; })) {
0418       // Get the number of bins
0419       int nbins0 = t->GetNbinsX();
0420       int nbins1 = t->GetNbinsY();
0421       // The material matrix
0422       MaterialSlabMatrix materialMatrix(
0423           nbins1, MaterialSlabVector(nbins0, MaterialSlab::Nothing()));
0424       // Fill the matrix from the histogram content
0425       for (int ib0 = 1; ib0 <= nbins0; ++ib0) {
0426         for (int ib1 = 1; ib1 <= nbins1; ++ib1) {
0427           auto dt = static_cast<float>(t->GetBinContent(ib0, ib1));
0428           if (dt > 0.) {
0429             auto dx0 = static_cast<float>(x0->GetBinContent(ib0, ib1));
0430             auto dl0 = static_cast<float>(l0->GetBinContent(ib0, ib1));
0431             auto da = static_cast<float>(A->GetBinContent(ib0, ib1));
0432             auto dz = static_cast<float>(Z->GetBinContent(ib0, ib1));
0433             auto drho = static_cast<float>(rho->GetBinContent(ib0, ib1));
0434             // Create material properties
0435             const auto material =
0436                 Material::fromMassDensity(dx0, dl0, da, dz, drho);
0437             materialMatrix[ib1 - 1][ib0 - 1] = MaterialSlab(material, dt);
0438           }
0439         }
0440       }  // Construct the binned material with the right bin utility
0441       texturedSurfaceMaterial = std::make_shared<const BinnedSurfaceMaterial>(
0442           bUtility, std::move(materialMatrix));
0443     }
0444   } else {
0445     // Construct the names for histogram type storage
0446     std::string indexName = tdName + "/" + m_cfg.indexHistName;
0447     // Get the histograms
0448     auto ih = dynamic_cast<TH2I*>(rFile.Get(indexName.c_str()));
0449     if (ih != nullptr) {
0450       // Get the number of bins
0451       int nbins0 = ih->GetNbinsX();
0452       int nbins1 = ih->GetNbinsY();
0453       // The material matrix
0454       MaterialSlabMatrix materialMatrix(
0455           nbins1, MaterialSlabVector(nbins0, MaterialSlab::Nothing()));
0456       // Fill the matrix from the tree entries
0457       for (int ib0 = 1; ib0 <= nbins0; ++ib0) {
0458         for (int ib1 = 1; ib1 <= nbins1; ++ib1) {
0459           auto idx = static_cast<int>(ih->GetBinContent(ib0, ib1));
0460           indexedMaterialTree->GetEntry(idx);
0461           const auto material = Material::fromMassDensity(
0462               m_indexedMaterialTreePayload.hX0,
0463               m_indexedMaterialTreePayload.hL0, m_indexedMaterialTreePayload.hA,
0464               m_indexedMaterialTreePayload.hZ,
0465               m_indexedMaterialTreePayload.hRho);
0466           materialMatrix[ib1 - 1][ib0 - 1] =
0467               MaterialSlab(material, m_indexedMaterialTreePayload.ht);
0468         }
0469       }  // Construct the binned material with the right bin utility
0470       texturedSurfaceMaterial = std::make_shared<const BinnedSurfaceMaterial>(
0471           bUtility, std::move(materialMatrix));
0472     }
0473   }
0474
0475   return texturedSurfaceMaterial;
0476 }