Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:10:01 

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2017-2018 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 http://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/BinUtility.hpp>
0023 #include <Acts/Utilities/BinningType.hpp>
0024 
0025 #include <algorithm>
0026 #include <cstdio>
0027 #include <functional>
0028 #include <iostream>
0029 #include <stdexcept>
0030 #include <string>
0031 #include <tuple>
0032 #include <vector>
0033 
0034 #include <TFile.h>
0035 #include <TH1.h>
0036 #include <TH2.h>
0037 #include <TIterator.h>
0038 #include <TKey.h>
0039 #include <TList.h>
0040 #include <TObject.h>
0041 #include <boost/algorithm/string.hpp>
0042 #include <boost/algorithm/string/finder.hpp>
0043 #include <boost/algorithm/string/iter_find.hpp>
0044 
0045 namespace Acts {
0046 class ISurfaceMaterial;
0047 class IVolumeMaterial;
0048 }  // namespace Acts
0049 
0050 ActsExamples::RootMaterialDecorator::RootMaterialDecorator(
0051     const ActsExamples::RootMaterialDecorator::Config& config,
0052     Acts::Logging::Level level)
0053     : m_cfg(config),
0054       m_logger{Acts::getDefaultLogger("RootMaterialDecorator", level)} {
0055   // Validate the configuration
0056   if (m_cfg.folderSurfaceNameBase.empty()) {
0057     throw std::invalid_argument("Missing surface folder name base");
0058   } else if (m_cfg.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   // Setup 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   // Get the list of keys from the file
0071   TList* tlist = m_inputFile->GetListOfKeys();
0072   auto tIter = tlist->MakeIterator();
0073   tIter->Reset();
0074 
0075   // Iterate over the keys in the file
0076   while (TKey* key = (TKey*)(tIter->Next())) {
0077     // Remember the directory
0078     std::string tdName(key->GetName());
0079 
0080     ACTS_VERBOSE("Processing directory: " << tdName);
0081 
0082     // volume
0083     std::vector<std::string> splitNames;
0084     iter_split(splitNames, tdName,
0085                boost::algorithm::first_finder(m_cfg.voltag));
0086     // Surface Material
0087     if (splitNames[0] == m_cfg.folderSurfaceNameBase) {
0088       // The surface material to be read in for this
0089       std::shared_ptr<const Acts::ISurfaceMaterial> sMaterial = nullptr;
0090 
0091       boost::split(splitNames, splitNames[1], boost::is_any_of("_"));
0092       Acts::GeometryIdentifier::Value volID = std::stoi(splitNames[0]);
0093       // boundary
0094       iter_split(splitNames, tdName,
0095                  boost::algorithm::first_finder(m_cfg.boutag));
0096       boost::split(splitNames, splitNames[1], boost::is_any_of("_"));
0097       Acts::GeometryIdentifier::Value bouID = std::stoi(splitNames[0]);
0098       // layer
0099       iter_split(splitNames, tdName,
0100                  boost::algorithm::first_finder(m_cfg.laytag));
0101       boost::split(splitNames, splitNames[1], boost::is_any_of("_"));
0102       Acts::GeometryIdentifier::Value layID = std::stoi(splitNames[0]);
0103       // approach
0104       iter_split(splitNames, tdName,
0105                  boost::algorithm::first_finder(m_cfg.apptag));
0106       boost::split(splitNames, splitNames[1], boost::is_any_of("_"));
0107       Acts::GeometryIdentifier::Value appID = std::stoi(splitNames[0]);
0108       // sensitive
0109       iter_split(splitNames, tdName,
0110                  boost::algorithm::first_finder(m_cfg.sentag));
0111       Acts::GeometryIdentifier::Value senID = std::stoi(splitNames[1]);
0112 
0113       // Reconstruct the geometry ID
0114       Acts::GeometryIdentifier geoID;
0115       geoID.setVolume(volID);
0116       geoID.setBoundary(bouID);
0117       geoID.setLayer(layID);
0118       geoID.setApproach(appID);
0119       geoID.setSensitive(senID);
0120       ACTS_VERBOSE("GeometryIdentifier re-constructed as " << geoID);
0121 
0122       // Construct the names
0123       std::string nName = tdName + "/" + m_cfg.ntag;
0124       std::string vName = tdName + "/" + m_cfg.vtag;
0125       std::string oName = tdName + "/" + m_cfg.otag;
0126       std::string minName = tdName + "/" + m_cfg.mintag;
0127       std::string maxName = tdName + "/" + m_cfg.maxtag;
0128       std::string tName = tdName + "/" + m_cfg.ttag;
0129       std::string x0Name = tdName + "/" + m_cfg.x0tag;
0130       std::string l0Name = tdName + "/" + m_cfg.l0tag;
0131       std::string aName = tdName + "/" + m_cfg.atag;
0132       std::string zName = tdName + "/" + m_cfg.ztag;
0133       std::string rhoName = tdName + "/" + m_cfg.rhotag;
0134 
0135       // Get the histograms
0136       TH1F* n = dynamic_cast<TH1F*>(m_inputFile->Get(nName.c_str()));
0137       TH1F* v = dynamic_cast<TH1F*>(m_inputFile->Get(vName.c_str()));
0138       TH1F* o = dynamic_cast<TH1F*>(m_inputFile->Get(oName.c_str()));
0139       TH1F* min = dynamic_cast<TH1F*>(m_inputFile->Get(minName.c_str()));
0140       TH1F* max = dynamic_cast<TH1F*>(m_inputFile->Get(maxName.c_str()));
0141       TH2F* t = dynamic_cast<TH2F*>(m_inputFile->Get(tName.c_str()));
0142       TH2F* x0 = dynamic_cast<TH2F*>(m_inputFile->Get(x0Name.c_str()));
0143       TH2F* l0 = dynamic_cast<TH2F*>(m_inputFile->Get(l0Name.c_str()));
0144       TH2F* A = dynamic_cast<TH2F*>(m_inputFile->Get(aName.c_str()));
0145       TH2F* Z = dynamic_cast<TH2F*>(m_inputFile->Get(zName.c_str()));
0146       TH2F* rho = dynamic_cast<TH2F*>(m_inputFile->Get(rhoName.c_str()));
0147 
0148       std::vector<const TH1*> hists{n, v, o, min, max, t, x0, l0, A, Z, rho};
0149 
0150       // Only go on when you have all histograms
0151       if (std::all_of(hists.begin(), hists.end(),
0152                       [](const auto* hist) { return hist != nullptr; })) {
0153         // Get the number of bins
0154         int nbins0 = t->GetNbinsX();
0155         int nbins1 = t->GetNbinsY();
0156 
0157         // The material matrix
0158         Acts::MaterialSlabMatrix materialMatrix(
0159             nbins1, Acts::MaterialSlabVector(nbins0, Acts::MaterialSlab()));
0160 
0161         // We need binned material properties
0162         if (nbins0 * nbins1 > 1) {
0163           // Fill the matrix first
0164           for (int ib0 = 1; ib0 <= nbins0; ++ib0) {
0165             for (int ib1 = 1; ib1 <= nbins1; ++ib1) {
0166               double dt = t->GetBinContent(ib0, ib1);
0167               if (dt > 0.) {
0168                 double dx0 = x0->GetBinContent(ib0, ib1);
0169                 double dl0 = l0->GetBinContent(ib0, ib1);
0170                 double da = A->GetBinContent(ib0, ib1);
0171                 double dz = Z->GetBinContent(ib0, ib1);
0172                 double drho = rho->GetBinContent(ib0, ib1);
0173                 // Create material properties
0174                 const auto material =
0175                     Acts::Material::fromMassDensity(dx0, dl0, da, dz, drho);
0176                 materialMatrix[ib1 - 1][ib0 - 1] =
0177                     Acts::MaterialSlab(material, dt);
0178               }
0179             }
0180           }
0181 
0182           // Now reconstruct the bin untilities
0183           Acts::BinUtility bUtility;
0184           for (int ib = 1; ib < n->GetNbinsX() + 1; ++ib) {
0185             std::size_t nbins = std::size_t(n->GetBinContent(ib));
0186             Acts::BinningValue val = Acts::BinningValue(v->GetBinContent(ib));
0187             Acts::BinningOption opt = Acts::BinningOption(o->GetBinContent(ib));
0188             float rmin = min->GetBinContent(ib);
0189             float rmax = max->GetBinContent(ib);
0190             bUtility += Acts::BinUtility(nbins, rmin, rmax, opt, val);
0191           }
0192           ACTS_VERBOSE("Created " << bUtility);
0193 
0194           // Construct the binned material with the right bin utility
0195           sMaterial = std::make_shared<const Acts::BinnedSurfaceMaterial>(
0196               bUtility, std::move(materialMatrix));
0197 
0198         } else {
0199           // Only homogeneous material present
0200           double dt = t->GetBinContent(1, 1);
0201           double dx0 = x0->GetBinContent(1, 1);
0202           double dl0 = l0->GetBinContent(1, 1);
0203           double da = A->GetBinContent(1, 1);
0204           double dz = Z->GetBinContent(1, 1);
0205           double drho = rho->GetBinContent(1, 1);
0206           // Create and set the homogeneous surface material
0207           const auto material =
0208               Acts::Material::fromMassDensity(dx0, dl0, da, dz, drho);
0209           sMaterial = std::make_shared<const Acts::HomogeneousSurfaceMaterial>(
0210               Acts::MaterialSlab(material, dt));
0211         }
0212       }
0213       ACTS_VERBOSE("Successfully read Material for : " << geoID);
0214 
0215       // Insert into the new collection
0216       m_surfaceMaterialMap.insert({geoID, std::move(sMaterial)});
0217 
0218     } else if (splitNames[0] == m_cfg.folderVolumeNameBase) {
0219       // The volume material to be read in for this
0220       std::shared_ptr<const Acts::IVolumeMaterial> vMaterial = nullptr;
0221       // Volume key
0222       boost::split(splitNames, splitNames[1], boost::is_any_of("_"));
0223       Acts::GeometryIdentifier::Value volID = std::stoi(splitNames[0]);
0224 
0225       // Reconstruct the geometry ID
0226       Acts::GeometryIdentifier geoID;
0227       geoID.setVolume(volID);
0228       ACTS_VERBOSE("GeometryIdentifier re-constructed as " << geoID);
0229 
0230       // Construct the names
0231       std::string nName = tdName + "/" + m_cfg.ntag;
0232       std::string vName = tdName + "/" + m_cfg.vtag;
0233       std::string oName = tdName + "/" + m_cfg.otag;
0234       std::string minName = tdName + "/" + m_cfg.mintag;
0235       std::string maxName = tdName + "/" + m_cfg.maxtag;
0236       std::string x0Name = tdName + "/" + m_cfg.x0tag;
0237       std::string l0Name = tdName + "/" + m_cfg.l0tag;
0238       std::string aName = tdName + "/" + m_cfg.atag;
0239       std::string zName = tdName + "/" + m_cfg.ztag;
0240       std::string rhoName = tdName + "/" + m_cfg.rhotag;
0241 
0242       // Get the histograms
0243       TH1F* n = dynamic_cast<TH1F*>(m_inputFile->Get(nName.c_str()));
0244       TH1F* v = dynamic_cast<TH1F*>(m_inputFile->Get(vName.c_str()));
0245       TH1F* o = dynamic_cast<TH1F*>(m_inputFile->Get(oName.c_str()));
0246       TH1F* min = dynamic_cast<TH1F*>(m_inputFile->Get(minName.c_str()));
0247       TH1F* max = dynamic_cast<TH1F*>(m_inputFile->Get(maxName.c_str()));
0248       TH1F* x0 = dynamic_cast<TH1F*>(m_inputFile->Get(x0Name.c_str()));
0249       TH1F* l0 = dynamic_cast<TH1F*>(m_inputFile->Get(l0Name.c_str()));
0250       TH1F* A = dynamic_cast<TH1F*>(m_inputFile->Get(aName.c_str()));
0251       TH1F* Z = dynamic_cast<TH1F*>(m_inputFile->Get(zName.c_str()));
0252       TH1F* rho = dynamic_cast<TH1F*>(m_inputFile->Get(rhoName.c_str()));
0253 
0254       // Only go on when you have all the material histograms
0255       if ((x0 != nullptr) && (l0 != nullptr) && (A != nullptr) &&
0256           (Z != nullptr) && (rho != nullptr)) {
0257         // Get the number of grid points
0258         int points = x0->GetNbinsX();
0259         // If the bin information histograms are present the material is
0260         // either a 2D or a 3D grid
0261         if ((n != nullptr) && (v != nullptr) && (o != nullptr) &&
0262             (min != nullptr) && (max != nullptr)) {
0263           // Dimension of the grid
0264           int dim = n->GetNbinsX();
0265           // Now reconstruct the bin untilities
0266           Acts::BinUtility bUtility;
0267           for (int ib = 1; ib < dim + 1; ++ib) {
0268             std::size_t nbins = std::size_t(n->GetBinContent(ib));
0269             Acts::BinningValue val = Acts::BinningValue(v->GetBinContent(ib));
0270             Acts::BinningOption opt = Acts::BinningOption(o->GetBinContent(ib));
0271             float rmin = min->GetBinContent(ib);
0272             float rmax = max->GetBinContent(ib);
0273             bUtility += Acts::BinUtility(nbins, rmin, rmax, opt, val);
0274           }
0275           ACTS_VERBOSE("Created " << bUtility);
0276 
0277           if (dim == 2) {
0278             // 2D Grid material
0279             std::function<Acts::Vector2(Acts::Vector3)> transfoGlobalToLocal;
0280             Acts::Grid2D grid = createGrid2D(bUtility, transfoGlobalToLocal);
0281 
0282             Acts::Grid2D::point_t gMin = grid.minPosition();
0283             Acts::Grid2D::point_t gMax = grid.maxPosition();
0284             Acts::Grid2D::index_t gNBins = grid.numLocalBins();
0285 
0286             Acts::EAxis axis1(gMin[0], gMax[0], gNBins[0]);
0287             Acts::EAxis axis2(gMin[1], gMax[1], gNBins[1]);
0288 
0289             // Build the grid and fill it with data
0290             Acts::MaterialGrid2D mGrid(std::make_tuple(axis1, axis2));
0291 
0292             for (int p = 1; p <= points; p++) {
0293               double dx0 = x0->GetBinContent(p);
0294               double dl0 = l0->GetBinContent(p);
0295               double da = A->GetBinContent(p);
0296               double dz = Z->GetBinContent(p);
0297               double drho = rho->GetBinContent(p);
0298               // Create material properties
0299               const auto material =
0300                   Acts::Material::fromMassDensity(dx0, dl0, da, dz, drho);
0301               mGrid.at(p - 1) = material.parameters();
0302             }
0303             Acts::MaterialMapper<Acts::MaterialGrid2D> matMap(
0304                 transfoGlobalToLocal, mGrid);
0305             vMaterial = std::make_shared<Acts::InterpolatedMaterialMap<
0306                 Acts::MaterialMapper<Acts::MaterialGrid2D>>>(std::move(matMap),
0307                                                              bUtility);
0308           } else if (dim == 3) {
0309             // 3D Grid material
0310             std::function<Acts::Vector3(Acts::Vector3)> transfoGlobalToLocal;
0311             Acts::Grid3D grid = createGrid3D(bUtility, transfoGlobalToLocal);
0312 
0313             Acts::Grid3D::point_t gMin = grid.minPosition();
0314             Acts::Grid3D::point_t gMax = grid.maxPosition();
0315             Acts::Grid3D::index_t gNBins = grid.numLocalBins();
0316 
0317             Acts::EAxis axis1(gMin[0], gMax[0], gNBins[0]);
0318             Acts::EAxis axis2(gMin[1], gMax[1], gNBins[1]);
0319             Acts::EAxis axis3(gMin[2], gMax[2], gNBins[2]);
0320 
0321             // Build the grid and fill it with data
0322             Acts::MaterialGrid3D mGrid(std::make_tuple(axis1, axis2, axis3));
0323 
0324             for (int p = 1; p <= points; p++) {
0325               double dx0 = x0->GetBinContent(p);
0326               double dl0 = l0->GetBinContent(p);
0327               double da = A->GetBinContent(p);
0328               double dz = Z->GetBinContent(p);
0329               double drho = rho->GetBinContent(p);
0330               // Create material properties
0331               const auto material =
0332                   Acts::Material::fromMassDensity(dx0, dl0, da, dz, drho);
0333               mGrid.at(p - 1) = material.parameters();
0334             }
0335             Acts::MaterialMapper<Acts::MaterialGrid3D> matMap(
0336                 transfoGlobalToLocal, mGrid);
0337             vMaterial = std::make_shared<Acts::InterpolatedMaterialMap<
0338                 Acts::MaterialMapper<Acts::MaterialGrid3D>>>(std::move(matMap),
0339                                                              bUtility);
0340           }
0341         } else {
0342           // Homogeneous material
0343           double dx0 = x0->GetBinContent(1);
0344           double dl0 = l0->GetBinContent(1);
0345           double da = A->GetBinContent(1);
0346           double dz = Z->GetBinContent(1);
0347           double drho = rho->GetBinContent(1);
0348           // Create material properties
0349           const auto material =
0350               Acts::Material::fromMassDensity(dx0, dl0, da, dz, drho);
0351           vMaterial =
0352               std::make_shared<Acts::HomogeneousVolumeMaterial>(material);
0353         }
0354       }
0355       ACTS_VERBOSE("Successfully read Material for : " << geoID);
0356 
0357       // Insert into the new collection
0358       m_volumeMaterialMap.insert({geoID, std::move(vMaterial)});
0359 
0360       // Incorrect FolderName value
0361     } else {
0362       ACTS_ERROR(
0363           "Invalid FolderName, does not match any entry in the root file");
0364     }
0365   }
0366 }
0367 
0368 ActsExamples::RootMaterialDecorator::~RootMaterialDecorator() {
0369   if (m_inputFile != nullptr) {
0370     m_inputFile->Close();
0371   }
0372 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!