Root/src/RootMaterialWriter.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2017-2018 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 http://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, int(binningData[b - 1].bins()));
0134         v->SetBinContent(b, int(binningData[b - 1].binvalue));
0135         o->SetBinContent(b, 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
0196     // get the geometry ID
0197     Acts::GeometryIdentifier geoID = key;
0198     // decode the geometryID
0199     const auto gvolID = geoID.volume();
0200
0201     // create the directory
0202     std::string tdName = m_cfg.folderVolumeNameBase.c_str();
0203     tdName += m_cfg.voltag + std::to_string(gvolID);
0204
0205     // create a new directory
0206     m_outputFile->mkdir(tdName.c_str());
0207     m_outputFile->cd(tdName.c_str());
0208
0209     ACTS_VERBOSE("Writing out map at " << tdName);
0210
0211     // understand what sort of material you have in mind
0212     auto bvMaterial3D = dynamic_cast<const Acts::InterpolatedMaterialMap<
0213         Acts::MaterialMapper<Acts::MaterialGrid3D>>*>(vMaterial);
0214     auto bvMaterial2D = dynamic_cast<const Acts::InterpolatedMaterialMap<
0215         Acts::MaterialMapper<Acts::MaterialGrid2D>>*>(vMaterial);
0216
0217     std::size_t points = 1;
0218     if (bvMaterial3D != nullptr || bvMaterial2D != nullptr) {
0219       // Get the binning data
0220       std::vector<Acts::BinningData> binningData;
0221       if (bvMaterial3D != nullptr) {
0222         binningData = bvMaterial3D->binUtility().binningData();
0223         Acts::MaterialGrid3D grid = bvMaterial3D->getMapper().getGrid();
0224         points = grid.size();
0225       } else {
0226         binningData = bvMaterial2D->binUtility().binningData();
0227         Acts::MaterialGrid2D grid = bvMaterial2D->getMapper().getGrid();
0228         points = grid.size();
0229       }
0230
0231       // 2-D or 3-D maps
0232       std::size_t binningBins = binningData.size();
0233
0234       // The bin number information
0235       TH1F* n = new TH1F(m_cfg.ntag.c_str(), "bins; bin", binningBins, -0.5,
0236                          binningBins - 0.5);
0237
0238       // The binning value information
0239       TH1F* v = new TH1F(m_cfg.vtag.c_str(), "binning values; bin", binningBins,
0240                          -0.5, binningBins - 0.5);
0241
0242       // The binning option information
0243       TH1F* o = new TH1F(m_cfg.otag.c_str(), "binning options; bin",
0244                          binningBins, -0.5, binningBins - 0.5);
0245
0246       // The binning option information
0247       TH1F* min = new TH1F(m_cfg.mintag.c_str(), "min; bin", binningBins, -0.5,
0248                            binningBins - 0.5);
0249
0250       // The binning option information
0251       TH1F* max = new TH1F(m_cfg.maxtag.c_str(), "max; bin", binningBins, -0.5,
0252                            binningBins - 0.5);
0253
0254       // Now fill the histogram content
0255       std::size_t b = 1;
0256       for (auto bData : binningData) {
0257         // Fill: nbins, value, option, min, max
0258         n->SetBinContent(b, int(binningData[b - 1].bins()));
0259         v->SetBinContent(b, int(binningData[b - 1].binvalue));
0260         o->SetBinContent(b, int(binningData[b - 1].option));
0261         min->SetBinContent(b, binningData[b - 1].min);
0262         max->SetBinContent(b, binningData[b - 1].max);
0263         ++b;
0264       }
0265       n->Write();
0266       v->Write();
0267       o->Write();
0268       min->Write();
0269       max->Write();
0270     }
0271
0272     TH1F* x0 = new TH1F(m_cfg.x0tag.c_str(), "X_{0} [mm] ;gridPoint", points,
0273                         -0.5, points - 0.5);
0274     TH1F* l0 = new TH1F(m_cfg.l0tag.c_str(), "#Lambda_{0} [mm] ;gridPoint",
0275                         points, -0.5, points - 0.5);
0276     TH1F* A = new TH1F(m_cfg.atag.c_str(), "X_{0} [mm] ;gridPoint", points,
0277                        -0.5, points - 0.5);
0278     TH1F* Z = new TH1F(m_cfg.ztag.c_str(), "#Lambda_{0} [mm] ;gridPoint",
0279                        points, -0.5, points - 0.5);
0280     TH1F* rho = new TH1F(m_cfg.rhotag.c_str(), "#rho [g/mm^3] ;gridPoint",
0281                          points, -0.5, points - 0.5);
0282     // homogeneous volume
0283     if (points == 1) {
0284       auto mat = vMaterial->material({0, 0, 0});
0285       x0->SetBinContent(1, mat.X0());
0286       l0->SetBinContent(1, mat.L0());
0287       A->SetBinContent(1, mat.Ar());
0288       Z->SetBinContent(1, mat.Z());
0289       rho->SetBinContent(1, mat.massDensity());
0290     } else {
0291       // 3d grid volume
0292       if (bvMaterial3D != nullptr) {
0293         Acts::MaterialGrid3D grid = bvMaterial3D->getMapper().getGrid();
0294         for (std::size_t point = 0; point < points; point++) {
0295           auto mat = Acts::Material(grid.at(point));
0296           if (mat) {
0297             x0->SetBinContent(point + 1, mat.X0());
0298             l0->SetBinContent(point + 1, mat.L0());
0299             A->SetBinContent(point + 1, mat.Ar());
0300             Z->SetBinContent(point + 1, mat.Z());
0301             rho->SetBinContent(point + 1, mat.massDensity());
0302           }
0303         }
0304       }
0305       // 2d grid volume
0306       else if (bvMaterial2D != nullptr) {
0307         Acts::MaterialGrid2D grid = bvMaterial2D->getMapper().getGrid();
0308         for (std::size_t point = 0; point < points; point++) {
0309           auto mat = Acts::Material(grid.at(point));
0310           if (mat) {
0311             x0->SetBinContent(point + 1, mat.X0());
0312             l0->SetBinContent(point + 1, mat.L0());
0313             A->SetBinContent(point + 1, mat.Ar());
0314             Z->SetBinContent(point + 1, mat.Z());
0315             rho->SetBinContent(point + 1, mat.massDensity());
0316           }
0317         }
0318       }
0319     }
0320     x0->Write();
0321     l0->Write();
0322     A->Write();
0323     Z->Write();
0324     rho->Write();
0325   }
0326 }
0327
0328 ActsExamples::RootMaterialWriter::~RootMaterialWriter() {
0329   if (m_outputFile != nullptr) {
0330     m_outputFile->Close();
0331   }
0332 }
0333
0334 void ActsExamples::RootMaterialWriter::write(
0335     const Acts::TrackingGeometry& tGeometry) {
0336   // Create a detector material map and loop recursively through it
0337   Acts::DetectorMaterialMaps detMatMap;
0338   auto hVolume = tGeometry.highestTrackingVolume();
0339   if (hVolume != nullptr) {
0340     collectMaterial(*hVolume, detMatMap);
0341   }
0342   // Write the resulting map to the file
0343   writeMaterial(detMatMap);
0344 }
0345
0346 void ActsExamples::RootMaterialWriter::collectMaterial(
0347     const Acts::TrackingVolume& tVolume,
0348     Acts::DetectorMaterialMaps& detMatMap) {
0349   // If the volume has volume material, write that
0350   if (tVolume.volumeMaterialSharedPtr() != nullptr && m_cfg.processVolumes) {
0351     detMatMap.second[tVolume.geometryId()] = tVolume.volumeMaterialSharedPtr();
0352   }
0353
0354   // If confined layers exist, loop over them and collect the layer material
0355   if (tVolume.confinedLayers() != nullptr) {
0356     ACTS_VERBOSE("Collecting material for " << tVolume.volumeName()
0357                                             << " layers");
0358     for (auto& lay : tVolume.confinedLayers()->arrayObjects()) {
0359       collectMaterial(*lay, detMatMap);
0360     }
0361   }
0362
0363   // If any of the boundary surfaces has material collect that
0364   if (m_cfg.processBoundaries) {
0365     for (auto& bou : tVolume.boundarySurfaces()) {
0366       const auto& bSurface = bou->surfaceRepresentation();
0367       if (bSurface.surfaceMaterialSharedPtr() != nullptr) {
0368         detMatMap.first[bSurface.geometryId()] =
0369             bSurface.surfaceMaterialSharedPtr();
0370       }
0371     }
0372   }
0373
0374   // If the volume has sub volumes, step down
0375   if (tVolume.confinedVolumes() != nullptr) {
0376     for (auto& tvol : tVolume.confinedVolumes()->arrayObjects()) {
0377       collectMaterial(*tvol, detMatMap);
0378     }
0379   }
0380 }
0381
0382 void ActsExamples::RootMaterialWriter::collectMaterial(
0383     const Acts::Layer& tLayer, Acts::DetectorMaterialMaps& detMatMap) {
0384   // If the representing surface has material, collect it
0385   const auto& rSurface = tLayer.surfaceRepresentation();
0386   if (rSurface.surfaceMaterialSharedPtr() != nullptr &&
0387       m_cfg.processRepresenting) {
0388     detMatMap.first[rSurface.geometryId()] =
0389         rSurface.surfaceMaterialSharedPtr();
0390   }
0391
0392   // Check the approach surfaces
0393   if (tLayer.approachDescriptor() != nullptr && m_cfg.processApproaches) {
0394     for (auto& aSurface : tLayer.approachDescriptor()->containedSurfaces()) {
0395       if (aSurface->surfaceMaterialSharedPtr() != nullptr) {
0396         detMatMap.first[aSurface->geometryId()] =
0397             aSurface->surfaceMaterialSharedPtr();
0398       }
0399     }
0400   }
0401
0402   // Check the sensitive surfaces
0403   if (tLayer.surfaceArray() != nullptr && m_cfg.processSensitives) {
0404     // sensitive surface loop
0405     for (auto& sSurface : tLayer.surfaceArray()->surfaces()) {
0406       if (sSurface->surfaceMaterialSharedPtr() != nullptr) {
0407         detMatMap.first[sSurface->geometryId()] =
0408             sSurface->surfaceMaterialSharedPtr();
0409       }
0410     }
0411   }
0412 }