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& [surfaceMaps, volumeMaps] = detMaterial;
0068
0069   // Write the surface material maps
0070   ActsPlugins::RootMaterialMapIo accessor(m_cfg.accessorConfig,
0071                                           m_logger->clone("RootMaterialMapIo"));
0072
0073   for (const auto& [geoId, sMap] : surfaceMaps) {
0074     // Get the Surface material
0075     accessor.write(*m_outputFile, geoId, *sMap, m_cfg.accessorOptions);
0076   }
0077
0078   // Write the volume material maps
0079   for (auto& [key, value] : volumeMaps) {
0080     // Get the Volume material
0081     const Acts::IVolumeMaterial* vMaterial = value.get();
0082     if (vMaterial == nullptr) {
0083       ACTS_WARNING("No material for volume " << key << " skipping");
0084       continue;
0085     }
0086
0087     // get the geometry ID
0088     Acts::GeometryIdentifier geoID = key;
0089     // decode the geometryID
0090     const auto gvolID = geoID.volume();
0091
0092     // create the directory
0093     std::string tdName = m_cfg.accessorOptions.folderVolumeNameBase.c_str();
0094     tdName += m_cfg.accessorConfig.volumePrefix + std::to_string(gvolID);
0095
0096     // create a new directory
0097     m_outputFile->mkdir(tdName.c_str());
0098     m_outputFile->cd(tdName.c_str());
0099
0100     ACTS_VERBOSE("Writing out map at " << tdName);
0101
0102     // understand what sort of material you have in mind
0103     auto bvMaterial3D = dynamic_cast<const Acts::InterpolatedMaterialMap<
0104         Acts::MaterialMapLookup<Acts::MaterialGrid3D>>*>(vMaterial);
0105     auto bvMaterial2D = dynamic_cast<const Acts::InterpolatedMaterialMap<
0106         Acts::MaterialMapLookup<Acts::MaterialGrid2D>>*>(vMaterial);
0107
0108     int points = 1;
0109     if (bvMaterial3D != nullptr || bvMaterial2D != nullptr) {
0110       // Get the binning data
0111       std::vector<Acts::BinningData> binningData;
0112       if (bvMaterial3D != nullptr) {
0113         binningData = bvMaterial3D->binUtility().binningData();
0114         Acts::MaterialGrid3D grid = bvMaterial3D->getMapper().getGrid();
0115         points = static_cast<int>(grid.size());
0116       } else {
0117         binningData = bvMaterial2D->binUtility().binningData();
0118         Acts::MaterialGrid2D grid = bvMaterial2D->getMapper().getGrid();
0119         points = static_cast<int>(grid.size());
0120       }
0121
0122       // 2-D or 3-D maps
0123       auto bins = static_cast<int>(binningData.size());
0124       auto fBins = static_cast<float>(bins);
0125
0126       // The bin number information
0127       TH1F n(m_cfg.accessorConfig.nBinsHistName.c_str(), "bins; bin", bins,
0128              -0.5, fBins - 0.5);
0129
0130       // The binning value information
0131       TH1F v(m_cfg.accessorConfig.axisDirHistName.c_str(),
0132              "binning values; bin", bins, -0.5, fBins - 0.5);
0133
0134       // The binning option information
0135       TH1F o(m_cfg.accessorConfig.axisBoundaryTypeHistName.c_str(),
0136              "binning options; bin", bins, -0.5, fBins - 0.5);
0137
0138       // The binning option information
0139       TH1F rmin(m_cfg.accessorConfig.minRangeHistName.c_str(), "min; bin", bins,
0140                 -0.5, fBins - 0.5);
0141
0142       // The binning option information
0143       TH1F rmax(m_cfg.accessorConfig.maxRangeHistName.c_str(), "max; bin", bins,
0144                 -0.5, fBins - 0.5);
0145
0146       // Now fill the histogram content
0147       for (const auto& [b, bData] : enumerate(binningData)) {
0148         // Fill: nbins, value, option, min, max
0149         n.SetBinContent(static_cast<int>(b),
0150                         static_cast<int>(binningData[b - 1].bins()));
0151         v.SetBinContent(static_cast<int>(b),
0152                         static_cast<int>(binningData[b - 1].binvalue));
0153         o.SetBinContent(static_cast<int>(b),
0154                         static_cast<int>(binningData[b - 1].option));
0155         rmin.SetBinContent(static_cast<int>(b), binningData[b - 1].min);
0156         rmax.SetBinContent(static_cast<int>(b), binningData[b - 1].max);
0157       }
0158       n.Write();
0159       v.Write();
0160       o.Write();
0161       rmin.Write();
0162       rmax.Write();
0163     }
0164
0165     auto fPoints = static_cast<float>(points);
0166     TH1F x0(m_cfg.accessorConfig.x0HistName.c_str(), "X_{0} [mm] ;gridPoint",
0167             points, -0.5, fPoints - 0.5);
0168     TH1F l0(m_cfg.accessorConfig.l0HistName.c_str(),
0169             "#Lambda_{0} [mm] ;gridPoint", points, -0.5, fPoints - 0.5);
0170     TH1F A(m_cfg.accessorConfig.aHistName.c_str(), "X_{0} [mm] ;gridPoint",
0171            points, -0.5, fPoints - 0.5);
0172     TH1F Z(m_cfg.accessorConfig.zHistName.c_str(),
0173            "#Lambda_{0} [mm] ;gridPoint", points, -0.5, fPoints - 0.5);
0174     TH1F rho(m_cfg.accessorConfig.rhoHistName.c_str(),
0175              "#rho [g/mm^3] ;gridPoint", points, -0.5, fPoints - 0.5);
0176     // homogeneous volume
0177     if (points == 1) {
0178       auto mat = vMaterial->material({0, 0, 0});
0179       x0.SetBinContent(1, mat.X0());
0180       l0.SetBinContent(1, mat.L0());
0181       A.SetBinContent(1, mat.Ar());
0182       Z.SetBinContent(1, mat.Z());
0183       rho.SetBinContent(1, mat.massDensity());
0184     } else {
0185       // 3d grid volume
0186       if (bvMaterial3D != nullptr) {
0187         Acts::MaterialGrid3D grid = bvMaterial3D->getMapper().getGrid();
0188         for (int point = 0; point < points; point++) {
0189           auto mat = Acts::Material(grid.at(point));
0190           if (!mat.isVacuum()) {
0191             x0.SetBinContent(point + 1, mat.X0());
0192             l0.SetBinContent(point + 1, mat.L0());
0193             A.SetBinContent(point + 1, mat.Ar());
0194             Z.SetBinContent(point + 1, mat.Z());
0195             rho.SetBinContent(point + 1, mat.massDensity());
0196           }
0197         }
0198       }
0199       // 2d grid volume
0200       else if (bvMaterial2D != nullptr) {
0201         Acts::MaterialGrid2D grid = bvMaterial2D->getMapper().getGrid();
0202         for (int point = 0; point < points; point++) {
0203           auto mat = Acts::Material(grid.at(point));
0204           if (!mat.isVacuum()) {
0205             x0.SetBinContent(point + 1, mat.X0());
0206             l0.SetBinContent(point + 1, mat.L0());
0207             A.SetBinContent(point + 1, mat.Ar());
0208             Z.SetBinContent(point + 1, mat.Z());
0209             rho.SetBinContent(point + 1, mat.massDensity());
0210           }
0211         }
0212       }
0213     }
0214     x0.Write();
0215     l0.Write();
0216     A.Write();
0217     Z.Write();
0218     rho.Write();
0219   }
0220 }
0221
0222 RootMaterialWriter::~RootMaterialWriter() {
0223   if (m_outputFile != nullptr) {
0224     m_outputFile->Close();
0225   }
0226 }
0227
0228 void RootMaterialWriter::write(const Acts::TrackingGeometry& tGeometry) {
0229   // Create a detector material map and loop recursively through it
0230   Acts::TrackingGeometryMaterial detMatMap;
0231   auto hVolume = tGeometry.highestTrackingVolume();
0232   if (hVolume != nullptr) {
0233     collectMaterial(*hVolume, detMatMap);
0234   }
0235   // Write the resulting map to the file
0236   writeMaterial(detMatMap);
0237 }
0238
0239 void RootMaterialWriter::collectMaterial(
0240     const Acts::TrackingVolume& tVolume,
0241     Acts::TrackingGeometryMaterial& detMatMap) {
0242   // If the volume has volume material, write that
0243   if (tVolume.volumeMaterialPtr() != nullptr && m_cfg.processVolumes) {
0244     detMatMap.second[tVolume.geometryId()] = tVolume.volumeMaterialPtr();
0245   }
0246
0247   // If confined layers exist, loop over them and collect the layer material
0248   if (tVolume.confinedLayers() != nullptr) {
0249     ACTS_VERBOSE("Collecting material for " << tVolume.volumeName()
0250                                             << " layers");
0251     for (auto& lay : tVolume.confinedLayers()->arrayObjects()) {
0252       collectMaterial(*lay, detMatMap);
0253     }
0254   }
0255
0256   // If any of the boundary surfaces has material collect that
0257   if (m_cfg.processBoundaries) {
0258     for (auto& bou : tVolume.boundarySurfaces()) {
0259       const auto& bSurface = bou->surfaceRepresentation();
0260       if (bSurface.surfaceMaterialSharedPtr() != nullptr) {
0261         detMatMap.first[bSurface.geometryId()] =
0262             bSurface.surfaceMaterialSharedPtr();
0263       }
0264     }
0265   }
0266
0267   // If the volume has sub volumes, step down
0268   if (tVolume.confinedVolumes() != nullptr) {
0269     for (auto& tvol : tVolume.confinedVolumes()->arrayObjects()) {
0270       collectMaterial(*tvol, detMatMap);
0271     }
0272   }
0273 }
0274
0275 void RootMaterialWriter::collectMaterial(
0276     const Acts::Layer& tLayer, Acts::TrackingGeometryMaterial& detMatMap) {
0277   // If the representing surface has material, collect it
0278   const auto& rSurface = tLayer.surfaceRepresentation();
0279   if (rSurface.surfaceMaterialSharedPtr() != nullptr &&
0280       m_cfg.processRepresenting) {
0281     detMatMap.first[rSurface.geometryId()] =
0282         rSurface.surfaceMaterialSharedPtr();
0283   }
0284
0285   // Check the approach surfaces
0286   if (tLayer.approachDescriptor() != nullptr && m_cfg.processApproaches) {
0287     for (auto& aSurface : tLayer.approachDescriptor()->containedSurfaces()) {
0288       if (aSurface->surfaceMaterialSharedPtr() != nullptr) {
0289         detMatMap.first[aSurface->geometryId()] =
0290             aSurface->surfaceMaterialSharedPtr();
0291       }
0292     }
0293   }
0294
0295   // Check the sensitive surfaces
0296   if (tLayer.surfaceArray() != nullptr && m_cfg.processSensitives) {
0297     // sensitive surface loop
0298     for (auto& sSurface : tLayer.surfaceArray()->surfaces()) {
0299       if (sSurface->surfaceMaterialSharedPtr() != nullptr) {
0300         detMatMap.first[sSurface->geometryId()] =
0301             sSurface->surfaceMaterialSharedPtr();
0302       }
0303     }
0304   }
0305 }
0306
0307 }  // namespace ActsExamples