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/Surfaces/Surface.hpp"
0024 #include "Acts/Surfaces/SurfaceArray.hpp"
0025 #include "Acts/Utilities/BinUtility.hpp"
0026 #include "Acts/Utilities/BinnedArray.hpp"
0027 #include "Acts/Utilities/BinningData.hpp"
0028 #include "Acts/Utilities/Enumerate.hpp"
0029 #include "Acts/Utilities/Logger.hpp"
0030 #include <Acts/Geometry/GeometryIdentifier.hpp>
0031 #include <Acts/Material/BinnedSurfaceMaterial.hpp>
0032
0033 #include <cstddef>
0034 #include <ios>
0035 #include <stdexcept>
0036 #include <type_traits>
0037 #include <vector>
0038
0039 #include <TFile.h>
0040 #include <TH1.h>
0041 #include <TH2.h>
0042
0043 ActsExamples::RootMaterialWriter::RootMaterialWriter(
0044     const ActsExamples::RootMaterialWriter::Config& config,
0045     Acts::Logging::Level level)
0046     : m_cfg(config),
0047       m_logger{Acts::getDefaultLogger("RootMaterialWriter", level)} {
0048   // Validate the configuration
0049   if (m_cfg.folderSurfaceNameBase.empty()) {
0050     throw std::invalid_argument("Missing surface folder name base");
0051   } else if (m_cfg.folderVolumeNameBase.empty()) {
0052     throw std::invalid_argument("Missing volume folder name base");
0053   } else if (m_cfg.filePath.empty()) {
0054     throw std::invalid_argument("Missing file name");
0055   }
0056
0057   // Setup ROOT I/O
0058   m_outputFile = TFile::Open(m_cfg.filePath.c_str(), m_cfg.fileMode.c_str());
0059   if (m_outputFile == nullptr) {
0060     throw std::ios_base::failure("Could not open '" + m_cfg.filePath + "'");
0061   }
0062 }
0063
0064 void ActsExamples::RootMaterialWriter::writeMaterial(
0065     const Acts::DetectorMaterialMaps& detMaterial) {
0066   // Change to the output file
0067   m_outputFile->cd();
0068
0069   auto& surfaceMaps = detMaterial.first;
0070   for (auto& [key, value] : surfaceMaps) {
0071     // Get the Surface material
0072     const Acts::ISurfaceMaterial* sMaterial = value.get();
0073
0074     // get the geometry ID
0075     Acts::GeometryIdentifier geoID = key;
0076     // decode the geometryID
0077     const auto gvolID = geoID.volume();
0078     const auto gbouID = geoID.boundary();
0079     const auto glayID = geoID.layer();
0080     const auto gappID = geoID.approach();
0081     const auto gsenID = geoID.sensitive();
0082     // create the directory
0083     std::string tdName = m_cfg.folderSurfaceNameBase.c_str();
0084     tdName += m_cfg.voltag + std::to_string(gvolID);
0085     tdName += m_cfg.boutag + std::to_string(gbouID);
0086     tdName += m_cfg.laytag + std::to_string(glayID);
0087     tdName += m_cfg.apptag + std::to_string(gappID);
0088     tdName += m_cfg.sentag + std::to_string(gsenID);
0089     // create a new directory
0090     m_outputFile->mkdir(tdName.c_str());
0091     m_outputFile->cd(tdName.c_str());
0092
0093     ACTS_VERBOSE("Writing out map at " << tdName);
0094
0095     std::size_t bins0 = 1, bins1 = 1;
0096     // understand what sort of material you have in mind
0097     const Acts::BinnedSurfaceMaterial* bsm =
0098         dynamic_cast<const Acts::BinnedSurfaceMaterial*>(sMaterial);
0099     if (bsm != nullptr) {
0100       // overwrite the bin numbers
0101       bins0 = bsm->binUtility().bins(0);
0102       bins1 = bsm->binUtility().bins(1);
0103
0104       // Get the binning data
0105       auto& binningData = bsm->binUtility().binningData();
0106       // 1-D or 2-D maps
0107       std::size_t binningBins = binningData.size();
0108
0109       // The bin number information
0110       TH1F* n = new TH1F(m_cfg.ntag.c_str(), "bins; bin", binningBins, -0.5,
0111                          binningBins - 0.5);
0112
0113       // The binning value information
0114       TH1F* v = new TH1F(m_cfg.vtag.c_str(), "binning values; bin", binningBins,
0115                          -0.5, binningBins - 0.5);
0116
0117       // The binning option information
0118       TH1F* o = new TH1F(m_cfg.otag.c_str(), "binning options; bin",
0119                          binningBins, -0.5, binningBins - 0.5);
0120
0121       // The binning option information
0122       TH1F* min = new TH1F(m_cfg.mintag.c_str(), "min; bin", binningBins, -0.5,
0123                            binningBins - 0.5);
0124
0125       // The binning option information
0126       TH1F* max = new TH1F(m_cfg.maxtag.c_str(), "max; bin", binningBins, -0.5,
0127                            binningBins - 0.5);
0128
0129       // Now fill the histogram content
0130       std::size_t b = 1;
0131       for (auto bData : binningData) {
0132         // Fill: nbins, value, option, min, max
0133         n->SetBinContent(b, static_cast<int>(binningData[b - 1].bins()));
0134         v->SetBinContent(b, static_cast<int>(binningData[b - 1].binvalue));
0135         o->SetBinContent(b, static_cast<int>(binningData[b - 1].option));
0136         min->SetBinContent(b, binningData[b - 1].min);
0137         max->SetBinContent(b, binningData[b - 1].max);
0138         ++b;
0139       }
0140       n->Write();
0141       v->Write();
0142       o->Write();
0143       min->Write();
0144       max->Write();
0145     }
0146
0147     TH2F* t = new TH2F(m_cfg.ttag.c_str(), "thickness [mm] ;b0 ;b1", bins0,
0148                        -0.5, bins0 - 0.5, bins1, -0.5, bins1 - 0.5);
0149     TH2F* x0 = new TH2F(m_cfg.x0tag.c_str(), "X_{0} [mm] ;b0 ;b1", bins0, -0.5,
0150                         bins0 - 0.5, bins1, -0.5, bins1 - 0.5);
0151     TH2F* l0 = new TH2F(m_cfg.l0tag.c_str(), "#Lambda_{0} [mm] ;b0 ;b1", bins0,
0152                         -0.5, bins0 - 0.5, bins1, -0.5, bins1 - 0.5);
0153     TH2F* A = new TH2F(m_cfg.atag.c_str(), "X_{0} [mm] ;b0 ;b1", bins0, -0.5,
0154                        bins0 - 0.5, bins1, -0.5, bins1 - 0.5);
0155     TH2F* Z = new TH2F(m_cfg.ztag.c_str(), "#Lambda_{0} [mm] ;b0 ;b1", bins0,
0156                        -0.5, bins0 - 0.5, bins1, -0.5, bins1 - 0.5);
0157     TH2F* rho = new TH2F(m_cfg.rhotag.c_str(), "#rho [g/mm^3] ;b0 ;b1", bins0,
0158                          -0.5, bins0 - 0.5, bins1, -0.5, bins1 - 0.5);
0159
0160     // loop over the material and fill
0161     if (bsm != nullptr) {
0162       const auto& materialMatrix = bsm->fullMaterial();
0163       for (auto [b1, materialVector] : Acts::enumerate(materialMatrix)) {
0164         for (auto [b0, mat] : Acts::enumerate(materialVector)) {
0165           t->SetBinContent(b0 + 1, b1 + 1, mat.thickness());
0166           x0->SetBinContent(b0 + 1, b1 + 1, mat.material().X0());
0167           l0->SetBinContent(b0 + 1, b1 + 1, mat.material().L0());
0168           A->SetBinContent(b0 + 1, b1 + 1, mat.material().Ar());
0169           Z->SetBinContent(b0 + 1, b1 + 1, mat.material().Z());
0170           rho->SetBinContent(b0 + 1, b1 + 1, mat.material().massDensity());
0171         }
0172       }
0173     } else if (bins1 == 1 && bins0 == 1) {
0174       // homogeneous surface
0175       auto mat = sMaterial->materialSlab(Acts::Vector3{0, 0, 0});
0176       t->SetBinContent(1, 1, mat.thickness());
0177       x0->SetBinContent(1, 1, mat.material().X0());
0178       l0->SetBinContent(1, 1, mat.material().L0());
0179       A->SetBinContent(1, 1, mat.material().Ar());
0180       Z->SetBinContent(1, 1, mat.material().Z());
0181       rho->SetBinContent(1, 1, mat.material().massDensity());
0182     }
0183     t->Write();
0184     x0->Write();
0185     l0->Write();
0186     A->Write();
0187     Z->Write();
0188     rho->Write();
0189   }
0190
0191   auto& volumeMaps = detMaterial.second;
0192   for (auto& [key, value] : volumeMaps) {
0193     // Get the Volume material
0194     const Acts::IVolumeMaterial* vMaterial = value.get();
0195     if (vMaterial == nullptr) {
0196       ACTS_WARNING("No material for volume " << key << " skipping");
0197       continue;
0198     }
0199
0200     // get the geometry ID
0201     Acts::GeometryIdentifier geoID = key;
0202     // decode the geometryID
0203     const auto gvolID = geoID.volume();
0204
0205     // create the directory
0206     std::string tdName = m_cfg.folderVolumeNameBase.c_str();
0207     tdName += m_cfg.voltag + std::to_string(gvolID);
0208
0209     // create a new directory
0210     m_outputFile->mkdir(tdName.c_str());
0211     m_outputFile->cd(tdName.c_str());
0212
0213     ACTS_VERBOSE("Writing out map at " << tdName);
0214
0215     // understand what sort of material you have in mind
0216     auto bvMaterial3D = dynamic_cast<const Acts::InterpolatedMaterialMap<
0217         Acts::MaterialMapper<Acts::MaterialGrid3D>>*>(vMaterial);
0218     auto bvMaterial2D = dynamic_cast<const Acts::InterpolatedMaterialMap<
0219         Acts::MaterialMapper<Acts::MaterialGrid2D>>*>(vMaterial);
0220
0221     std::size_t points = 1;
0222     if (bvMaterial3D != nullptr || bvMaterial2D != nullptr) {
0223       // Get the binning data
0224       std::vector<Acts::BinningData> binningData;
0225       if (bvMaterial3D != nullptr) {
0226         binningData = bvMaterial3D->binUtility().binningData();
0227         Acts::MaterialGrid3D grid = bvMaterial3D->getMapper().getGrid();
0228         points = grid.size();
0229       } else {
0230         binningData = bvMaterial2D->binUtility().binningData();
0231         Acts::MaterialGrid2D grid = bvMaterial2D->getMapper().getGrid();
0232         points = grid.size();
0233       }
0234
0235       // 2-D or 3-D maps
0236       std::size_t binningBins = binningData.size();
0237
0238       // The bin number information
0239       TH1F* n = new TH1F(m_cfg.ntag.c_str(), "bins; bin", binningBins, -0.5,
0240                          binningBins - 0.5);
0241
0242       // The binning value information
0243       TH1F* v = new TH1F(m_cfg.vtag.c_str(), "binning values; bin", binningBins,
0244                          -0.5, binningBins - 0.5);
0245
0246       // The binning option information
0247       TH1F* o = new TH1F(m_cfg.otag.c_str(), "binning options; bin",
0248                          binningBins, -0.5, binningBins - 0.5);
0249
0250       // The binning option information
0251       TH1F* min = new TH1F(m_cfg.mintag.c_str(), "min; bin", binningBins, -0.5,
0252                            binningBins - 0.5);
0253
0254       // The binning option information
0255       TH1F* max = new TH1F(m_cfg.maxtag.c_str(), "max; bin", binningBins, -0.5,
0256                            binningBins - 0.5);
0257
0258       // Now fill the histogram content
0259       std::size_t b = 1;
0260       for (auto bData : binningData) {
0261         // Fill: nbins, value, option, min, max
0262         n->SetBinContent(b, static_cast<int>(binningData[b - 1].bins()));
0263         v->SetBinContent(b, static_cast<int>(binningData[b - 1].binvalue));
0264         o->SetBinContent(b, static_cast<int>(binningData[b - 1].option));
0265         min->SetBinContent(b, binningData[b - 1].min);
0266         max->SetBinContent(b, binningData[b - 1].max);
0267         ++b;
0268       }
0269       n->Write();
0270       v->Write();
0271       o->Write();
0272       min->Write();
0273       max->Write();
0274     }
0275
0276     TH1F* x0 = new TH1F(m_cfg.x0tag.c_str(), "X_{0} [mm] ;gridPoint", points,
0277                         -0.5, points - 0.5);
0278     TH1F* l0 = new TH1F(m_cfg.l0tag.c_str(), "#Lambda_{0} [mm] ;gridPoint",
0279                         points, -0.5, points - 0.5);
0280     TH1F* A = new TH1F(m_cfg.atag.c_str(), "X_{0} [mm] ;gridPoint", points,
0281                        -0.5, points - 0.5);
0282     TH1F* Z = new TH1F(m_cfg.ztag.c_str(), "#Lambda_{0} [mm] ;gridPoint",
0283                        points, -0.5, points - 0.5);
0284     TH1F* rho = new TH1F(m_cfg.rhotag.c_str(), "#rho [g/mm^3] ;gridPoint",
0285                          points, -0.5, points - 0.5);
0286     // homogeneous volume
0287     if (points == 1) {
0288       auto mat = vMaterial->material({0, 0, 0});
0289       x0->SetBinContent(1, mat.X0());
0290       l0->SetBinContent(1, mat.L0());
0291       A->SetBinContent(1, mat.Ar());
0292       Z->SetBinContent(1, mat.Z());
0293       rho->SetBinContent(1, mat.massDensity());
0294     } else {
0295       // 3d grid volume
0296       if (bvMaterial3D != nullptr) {
0297         Acts::MaterialGrid3D grid = bvMaterial3D->getMapper().getGrid();
0298         for (std::size_t point = 0; point < points; point++) {
0299           auto mat = Acts::Material(grid.at(point));
0300           if (mat.isValid()) {
0301             x0->SetBinContent(point + 1, mat.X0());
0302             l0->SetBinContent(point + 1, mat.L0());
0303             A->SetBinContent(point + 1, mat.Ar());
0304             Z->SetBinContent(point + 1, mat.Z());
0305             rho->SetBinContent(point + 1, mat.massDensity());
0306           }
0307         }
0308       }
0309       // 2d grid volume
0310       else if (bvMaterial2D != nullptr) {
0311         Acts::MaterialGrid2D grid = bvMaterial2D->getMapper().getGrid();
0312         for (std::size_t point = 0; point < points; point++) {
0313           auto mat = Acts::Material(grid.at(point));
0314           if (mat.isValid()) {
0315             x0->SetBinContent(point + 1, mat.X0());
0316             l0->SetBinContent(point + 1, mat.L0());
0317             A->SetBinContent(point + 1, mat.Ar());
0318             Z->SetBinContent(point + 1, mat.Z());
0319             rho->SetBinContent(point + 1, mat.massDensity());
0320           }
0321         }
0322       }
0323     }
0324     x0->Write();
0325     l0->Write();
0326     A->Write();
0327     Z->Write();
0328     rho->Write();
0329   }
0330 }
0331
0332 ActsExamples::RootMaterialWriter::~RootMaterialWriter() {
0333   if (m_outputFile != nullptr) {
0334     m_outputFile->Close();
0335   }
0336 }
0337
0338 void ActsExamples::RootMaterialWriter::write(
0339     const Acts::TrackingGeometry& tGeometry) {
0340   // Create a detector material map and loop recursively through it
0341   Acts::DetectorMaterialMaps detMatMap;
0342   auto hVolume = tGeometry.highestTrackingVolume();
0343   if (hVolume != nullptr) {
0344     collectMaterial(*hVolume, detMatMap);
0345   }
0346   // Write the resulting map to the file
0347   writeMaterial(detMatMap);
0348 }
0349
0350 void ActsExamples::RootMaterialWriter::collectMaterial(
0351     const Acts::TrackingVolume& tVolume,
0352     Acts::DetectorMaterialMaps& detMatMap) {
0353   // If the volume has volume material, write that
0354   if (tVolume.volumeMaterialPtr() != nullptr && m_cfg.processVolumes) {
0355     detMatMap.second[tVolume.geometryId()] = tVolume.volumeMaterialPtr();
0356   }
0357
0358   // If confined layers exist, loop over them and collect the layer material
0359   if (tVolume.confinedLayers() != nullptr) {
0360     ACTS_VERBOSE("Collecting material for " << tVolume.volumeName()
0361                                             << " layers");
0362     for (auto& lay : tVolume.confinedLayers()->arrayObjects()) {
0363       collectMaterial(*lay, detMatMap);
0364     }
0365   }
0366
0367   // If any of the boundary surfaces has material collect that
0368   if (m_cfg.processBoundaries) {
0369     for (auto& bou : tVolume.boundarySurfaces()) {
0370       const auto& bSurface = bou->surfaceRepresentation();
0371       if (bSurface.surfaceMaterialSharedPtr() != nullptr) {
0372         detMatMap.first[bSurface.geometryId()] =
0373             bSurface.surfaceMaterialSharedPtr();
0374       }
0375     }
0376   }
0377
0378   // If the volume has sub volumes, step down
0379   if (tVolume.confinedVolumes() != nullptr) {
0380     for (auto& tvol : tVolume.confinedVolumes()->arrayObjects()) {
0381       collectMaterial(*tvol, detMatMap);
0382     }
0383   }
0384 }
0385
0386 void ActsExamples::RootMaterialWriter::collectMaterial(
0387     const Acts::Layer& tLayer, Acts::DetectorMaterialMaps& detMatMap) {
0388   // If the representing surface has material, collect it
0389   const auto& rSurface = tLayer.surfaceRepresentation();
0390   if (rSurface.surfaceMaterialSharedPtr() != nullptr &&
0391       m_cfg.processRepresenting) {
0392     detMatMap.first[rSurface.geometryId()] =
0393         rSurface.surfaceMaterialSharedPtr();
0394   }
0395
0396   // Check the approach surfaces
0397   if (tLayer.approachDescriptor() != nullptr && m_cfg.processApproaches) {
0398     for (auto& aSurface : tLayer.approachDescriptor()->containedSurfaces()) {
0399       if (aSurface->surfaceMaterialSharedPtr() != nullptr) {
0400         detMatMap.first[aSurface->geometryId()] =
0401             aSurface->surfaceMaterialSharedPtr();
0402       }
0403     }
0404   }
0405
0406   // Check the sensitive surfaces
0407   if (tLayer.surfaceArray() != nullptr && m_cfg.processSensitives) {
0408     // sensitive surface loop
0409     for (auto& sSurface : tLayer.surfaceArray()->surfaces()) {
0410       if (sSurface->surfaceMaterialSharedPtr() != nullptr) {
0411         detMatMap.first[sSurface->geometryId()] =
0412             sSurface->surfaceMaterialSharedPtr();
0413       }
0414     }
0415   }
0416 }