Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:11:56 

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