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