Root/src/RootMaterialWriter.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/RootMaterialWriter.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Geometry/ApproachDescriptor.hpp"
0013 #include "Acts/Geometry/BoundarySurfaceT.hpp"
0014 #include "Acts/Geometry/GeometryIdentifier.hpp"
0015 #include "Acts/Geometry/Layer.hpp"
0016 #include "Acts/Geometry/TrackingGeometry.hpp"
0017 #include "Acts/Geometry/TrackingVolume.hpp"
0018 #include "Acts/Material/BinnedSurfaceMaterial.hpp"
0019 #include "Acts/Material/IVolumeMaterial.hpp"
0020 #include "Acts/Material/InterpolatedMaterialMap.hpp"
0021 #include "Acts/Material/Material.hpp"
0022 #include "Acts/Material/MaterialGridHelper.hpp"
0023 #include "Acts/Material/TrackingGeometryMaterial.hpp"
0024 #include "Acts/Surfaces/Surface.hpp"
0025 #include "Acts/Surfaces/SurfaceArray.hpp"
0026 #include "Acts/Utilities/BinUtility.hpp"
0027 #include "Acts/Utilities/BinningData.hpp"
0028 #include "Acts/Utilities/Enumerate.hpp"
0029 #include "Acts/Utilities/Logger.hpp"
0030
0031 #include <cstddef>
0032 #include <ios>
0033 #include <stdexcept>
0034 #include <vector>
0035
0036 #include <TFile.h>
0037 #include <TH1.h>
0038 #include <TH2.h>
0039
0040 namespace ActsExamples {
0041
0042 RootMaterialWriter::RootMaterialWriter(const RootMaterialWriter::Config& config,
0043                                        Acts::Logging::Level level)
0044     : m_cfg(config),
0045       m_logger{Acts::getDefaultLogger("RootMaterialWriter", level)} {
0046   // Validate the configuration
0047   if (m_cfg.accessorOptions.folderSurfaceNameBase.empty()) {
0048     throw std::invalid_argument("Missing surface folder name base");
0049   } else if (m_cfg.accessorOptions.folderVolumeNameBase.empty()) {
0050     throw std::invalid_argument("Missing volume folder name base");
0051   } else if (m_cfg.filePath.empty()) {
0052     throw std::invalid_argument("Missing file name");
0053   }
0054
0055   // Setup ROOT I/O
0056   m_outputFile = TFile::Open(m_cfg.filePath.c_str(), m_cfg.fileMode.c_str());
0057   if (m_outputFile == nullptr) {
0058     throw std::ios_base::failure("Could not open '" + m_cfg.filePath + "'");
0059   }
0060 }
0061
0062 void RootMaterialWriter::writeMaterial(
0063     const Acts::TrackingGeometryMaterial& detMaterial) {
0064   // Change to the output file
0065   m_outputFile->cd();
0066
0067   const auto& volumeMaps = detMaterial.volumeMaterials;
0068   if (!detMaterial.keyedSurfaces.empty()) {
0069     throw std::invalid_argument(
0070         "ROOT material output does not support stable keys; use JSON");
0071   }
0072
0073   // Write the surface material maps
0074   ActsPlugins::RootMaterialMapIo accessor(m_cfg.accessorConfig,
0075                                           m_logger->clone("RootMaterialMapIo"));
0076
0077   accessor.write(*m_outputFile, detMaterial, m_cfg.accessorOptions);
0078
0079   // Write the volume material maps
0080   for (auto& [key, value] : volumeMaps) {
0081     // Get the Volume material
0082     const Acts::IVolumeMaterial* vMaterial = value.get();
0083     if (vMaterial == nullptr) {
0084       ACTS_WARNING("No material for volume " << key << " skipping");
0085       continue;
0086     }
0087
0088     // get the geometry ID
0089     Acts::GeometryIdentifier geoID = key;
0090     // decode the geometryID
0091     const auto gvolID = geoID.volume();
0092
0093     // create the directory
0094     std::string tdName = m_cfg.accessorOptions.folderVolumeNameBase.c_str();
0095     tdName += m_cfg.accessorConfig.volumePrefix + std::to_string(gvolID);
0096
0097     // create a new directory
0098     m_outputFile->mkdir(tdName.c_str());
0099     m_outputFile->cd(tdName.c_str());
0100
0101     ACTS_VERBOSE("Writing out map at " << tdName);
0102
0103     // understand what sort of material you have in mind
0104     auto bvMaterial3D = dynamic_cast<const Acts::InterpolatedMaterialMap<
0105         Acts::MaterialMapLookup<Acts::MaterialGrid3D>>*>(vMaterial);
0106     auto bvMaterial2D = dynamic_cast<const Acts::InterpolatedMaterialMap<
0107         Acts::MaterialMapLookup<Acts::MaterialGrid2D>>*>(vMaterial);
0108
0109     int points = 1;
0110     if (bvMaterial3D != nullptr || bvMaterial2D != nullptr) {
0111       // Get the binning data
0112       std::vector<Acts::BinningData> binningData;
0113       if (bvMaterial3D != nullptr) {
0114         binningData = bvMaterial3D->binUtility().binningData();
0115         Acts::MaterialGrid3D grid = bvMaterial3D->getMapper().getGrid();
0116         points = static_cast<int>(grid.size());
0117       } else {
0118         binningData = bvMaterial2D->binUtility().binningData();
0119         Acts::MaterialGrid2D grid = bvMaterial2D->getMapper().getGrid();
0120         points = static_cast<int>(grid.size());
0121       }
0122
0123       // 2-D or 3-D maps
0124       auto bins = static_cast<int>(binningData.size());
0125       auto fBins = static_cast<float>(bins);
0126
0127       // The bin number information
0128       TH1F n(m_cfg.accessorConfig.nBinsHistName.c_str(), "bins; bin", bins,
0129              -0.5, fBins - 0.5);
0130
0131       // The binning value information
0132       TH1F v(m_cfg.accessorConfig.axisDirHistName.c_str(),
0133              "binning values; bin", bins, -0.5, fBins - 0.5);
0134
0135       // The binning option information
0136       TH1F o(m_cfg.accessorConfig.axisBoundaryTypeHistName.c_str(),
0137              "binning options; bin", bins, -0.5, fBins - 0.5);
0138
0139       // The binning option information
0140       TH1F rmin(m_cfg.accessorConfig.minRangeHistName.c_str(), "min; bin", bins,
0141                 -0.5, fBins - 0.5);
0142
0143       // The binning option information
0144       TH1F rmax(m_cfg.accessorConfig.maxRangeHistName.c_str(), "max; bin", bins,
0145                 -0.5, fBins - 0.5);
0146
0147       // Now fill the histogram content
0148       for (const auto& [b, bData] : enumerate(binningData)) {
0149         // Fill: nbins, value, option, min, max
0150         n.SetBinContent(static_cast<int>(b),
0151                         static_cast<int>(binningData[b - 1].bins()));
0152         v.SetBinContent(static_cast<int>(b),
0153                         static_cast<int>(binningData[b - 1].binvalue));
0154         o.SetBinContent(static_cast<int>(b),
0155                         static_cast<int>(binningData[b - 1].option));
0156         rmin.SetBinContent(static_cast<int>(b), binningData[b - 1].min);
0157         rmax.SetBinContent(static_cast<int>(b), binningData[b - 1].max);
0158       }
0159       n.Write();
0160       v.Write();
0161       o.Write();
0162       rmin.Write();
0163       rmax.Write();
0164     }
0165
0166     auto fPoints = static_cast<float>(points);
0167     TH1F x0(m_cfg.accessorConfig.x0HistName.c_str(), "X_{0} [mm] ;gridPoint",
0168             points, -0.5, fPoints - 0.5);
0169     TH1F l0(m_cfg.accessorConfig.l0HistName.c_str(),
0170             "#Lambda_{0} [mm] ;gridPoint", points, -0.5, fPoints - 0.5);
0171     TH1F A(m_cfg.accessorConfig.aHistName.c_str(), "X_{0} [mm] ;gridPoint",
0172            points, -0.5, fPoints - 0.5);
0173     TH1F Z(m_cfg.accessorConfig.zHistName.c_str(),
0174            "#Lambda_{0} [mm] ;gridPoint", points, -0.5, fPoints - 0.5);
0175     TH1F rho(m_cfg.accessorConfig.rhoHistName.c_str(),
0176              "#rho [g/mm^3] ;gridPoint", points, -0.5, fPoints - 0.5);
0177     // homogeneous volume
0178     if (points == 1) {
0179       auto mat = vMaterial->material({0, 0, 0});
0180       x0.SetBinContent(1, mat.X0());
0181       l0.SetBinContent(1, mat.L0());
0182       A.SetBinContent(1, mat.Ar());
0183       Z.SetBinContent(1, mat.Z());
0184       rho.SetBinContent(1, mat.massDensity());
0185     } else {
0186       // 3d grid volume
0187       if (bvMaterial3D != nullptr) {
0188         Acts::MaterialGrid3D grid = bvMaterial3D->getMapper().getGrid();
0189         for (int point = 0; point < points; point++) {
0190           auto mat = Acts::Material(grid.at(point));
0191           if (!mat.isVacuum()) {
0192             x0.SetBinContent(point + 1, mat.X0());
0193             l0.SetBinContent(point + 1, mat.L0());
0194             A.SetBinContent(point + 1, mat.Ar());
0195             Z.SetBinContent(point + 1, mat.Z());
0196             rho.SetBinContent(point + 1, mat.massDensity());
0197           }
0198         }
0199       }
0200       // 2d grid volume
0201       else if (bvMaterial2D != nullptr) {
0202         Acts::MaterialGrid2D grid = bvMaterial2D->getMapper().getGrid();
0203         for (int point = 0; point < points; point++) {
0204           auto mat = Acts::Material(grid.at(point));
0205           if (!mat.isVacuum()) {
0206             x0.SetBinContent(point + 1, mat.X0());
0207             l0.SetBinContent(point + 1, mat.L0());
0208             A.SetBinContent(point + 1, mat.Ar());
0209             Z.SetBinContent(point + 1, mat.Z());
0210             rho.SetBinContent(point + 1, mat.massDensity());
0211           }
0212         }
0213       }
0214     }
0215     x0.Write();
0216     l0.Write();
0217     A.Write();
0218     Z.Write();
0219     rho.Write();
0220   }
0221 }
0222
0223 RootMaterialWriter::~RootMaterialWriter() {
0224   if (m_outputFile != nullptr) {
0225     m_outputFile->Close();
0226   }
0227 }
0228
0229 void RootMaterialWriter::write(const Acts::TrackingGeometry& tGeometry) {
0230   // Create a detector material map and loop recursively through it
0231   Acts::TrackingGeometryMaterial detMatMap;
0232   auto hVolume = tGeometry.highestTrackingVolume();
0233   if (hVolume != nullptr) {
0234     collectMaterial(*hVolume, detMatMap);
0235   }
0236   // Write the resulting map to the file
0237   writeMaterial(detMatMap);
0238 }
0239
0240 void RootMaterialWriter::collectMaterial(
0241     const Acts::TrackingVolume& tVolume,
0242     Acts::TrackingGeometryMaterial& detMatMap) {
0243   // If the volume has volume material, write that
0244   if (tVolume.volumeMaterialPtr() != nullptr && m_cfg.processVolumes) {
0245     detMatMap.volumeMaterials[tVolume.geometryId()] =
0246         tVolume.volumeMaterialPtr();
0247   }
0248
0249   // If confined layers exist, loop over them and collect the layer material
0250   if (tVolume.confinedLayers() != nullptr) {
0251     ACTS_VERBOSE("Collecting material for " << tVolume.volumeName()
0252                                             << " layers");
0253     for (auto& lay : tVolume.confinedLayers()->arrayObjects()) {
0254       collectMaterial(*lay, detMatMap);
0255     }
0256   }
0257
0258   // If any of the boundary surfaces has material collect that
0259   if (m_cfg.processBoundaries) {
0260     for (auto& bou : tVolume.boundarySurfaces()) {
0261       const auto& bSurface = bou->surfaceRepresentation();
0262       if (bSurface.surfaceMaterialSharedPtr() != nullptr) {
0263         detMatMap.surfaceMaterials[bSurface.geometryId()] =
0264             bSurface.surfaceMaterialSharedPtr();
0265       }
0266     }
0267   }
0268
0269   // If the volume has sub volumes, step down
0270   if (tVolume.confinedVolumes() != nullptr) {
0271     for (auto& tvol : tVolume.confinedVolumes()->arrayObjects()) {
0272       collectMaterial(*tvol, detMatMap);
0273     }
0274   }
0275 }
0276
0277 void RootMaterialWriter::collectMaterial(
0278     const Acts::Layer& tLayer, Acts::TrackingGeometryMaterial& detMatMap) {
0279   // If the representing surface has material, collect it
0280   const auto& rSurface = tLayer.surfaceRepresentation();
0281   if (rSurface.surfaceMaterialSharedPtr() != nullptr &&
0282       m_cfg.processRepresenting) {
0283     detMatMap.surfaceMaterials[rSurface.geometryId()] =
0284         rSurface.surfaceMaterialSharedPtr();
0285   }
0286
0287   // Check the approach surfaces
0288   if (tLayer.approachDescriptor() != nullptr && m_cfg.processApproaches) {
0289     for (auto& aSurface : tLayer.approachDescriptor()->containedSurfaces()) {
0290       if (aSurface->surfaceMaterialSharedPtr() != nullptr) {
0291         detMatMap.surfaceMaterials[aSurface->geometryId()] =
0292             aSurface->surfaceMaterialSharedPtr();
0293       }
0294     }
0295   }
0296
0297   // Check the sensitive surfaces
0298   if (tLayer.surfaceArray() != nullptr && m_cfg.processSensitives) {
0299     // sensitive surface loop
0300     for (auto& sSurface : tLayer.surfaceArray()->surfaces()) {
0301       if (sSurface->surfaceMaterialSharedPtr() != nullptr) {
0302         detMatMap.surfaceMaterials[sSurface->geometryId()] =
0303             sSurface->surfaceMaterialSharedPtr();
0304       }
0305     }
0306   }
0307 }
0308
0309 }  // namespace ActsExamples