Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RootMaterialWriter.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RootMaterialWriter.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2017-2018 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
9 #include "ActsExamples/Io/Root/RootMaterialWriter.hpp"
10 
11 #include "Acts/Geometry/ApproachDescriptor.hpp"
12 #include "Acts/Geometry/BoundarySurfaceT.hpp"
13 #include "Acts/Geometry/Layer.hpp"
14 #include "Acts/Geometry/TrackingGeometry.hpp"
15 #include "Acts/Geometry/TrackingVolume.hpp"
16 #include "Acts/Material/ISurfaceMaterial.hpp"
17 #include "Acts/Material/IVolumeMaterial.hpp"
18 #include "Acts/Material/InterpolatedMaterialMap.hpp"
19 #include "Acts/Material/Material.hpp"
20 #include "Acts/Material/MaterialGridHelper.hpp"
21 #include "Acts/Material/MaterialSlab.hpp"
22 #include "Acts/Surfaces/Surface.hpp"
23 #include "Acts/Surfaces/SurfaceArray.hpp"
24 #include "Acts/Utilities/BinUtility.hpp"
25 #include "Acts/Utilities/BinnedArray.hpp"
26 #include "Acts/Utilities/BinningData.hpp"
27 #include "Acts/Utilities/Logger.hpp"
28 #include <Acts/Geometry/GeometryIdentifier.hpp>
29 #include <Acts/Material/BinnedSurfaceMaterial.hpp>
30 
31 #include <cstddef>
32 #include <ios>
33 #include <stdexcept>
34 #include <type_traits>
35 #include <vector>
36 
37 #include <TFile.h>
38 #include <TH1.h>
39 #include <TH2.h>
40 
41 ActsExamples::RootMaterialWriter::RootMaterialWriter(
42  const ActsExamples::RootMaterialWriter::Config& config,
43  Acts::Logging::Level level)
44  : m_cfg(config),
45  m_logger{Acts::getDefaultLogger("RootMaterialWriter", level)} {
46  // Validate the configuration
47  if (m_cfg.folderSurfaceNameBase.empty()) {
48  throw std::invalid_argument("Missing surface folder name base");
49  } else if (m_cfg.folderVolumeNameBase.empty()) {
50  throw std::invalid_argument("Missing volume folder name base");
51  } else if (m_cfg.filePath.empty()) {
52  throw std::invalid_argument("Missing file name");
53  }
54 
55  // Setup ROOT I/O
56  m_outputFile = TFile::Open(m_cfg.filePath.c_str(), m_cfg.fileMode.c_str());
57  if (m_outputFile == nullptr) {
58  throw std::ios_base::failure("Could not open '" + m_cfg.filePath);
59  }
60 }
61 
62 void ActsExamples::RootMaterialWriter::writeMaterial(
63  const Acts::DetectorMaterialMaps& detMaterial) {
64  // Change to the output file
65  m_outputFile->cd();
66 
67  auto& surfaceMaps = detMaterial.first;
68  for (auto& [key, value] : surfaceMaps) {
69  // Get the Surface material
70  const Acts::ISurfaceMaterial* sMaterial = value.get();
71 
72  // get the geometry ID
73  Acts::GeometryIdentifier geoID = key;
74  // decode the geometryID
75  const auto gvolID = geoID.volume();
76  const auto gbouID = geoID.boundary();
77  const auto glayID = geoID.layer();
78  const auto gappID = geoID.approach();
79  const auto gsenID = geoID.sensitive();
80  // create the directory
81  std::string tdName = m_cfg.folderSurfaceNameBase.c_str();
82  tdName += m_cfg.voltag + std::to_string(gvolID);
83  tdName += m_cfg.boutag + std::to_string(gbouID);
84  tdName += m_cfg.laytag + std::to_string(glayID);
85  tdName += m_cfg.apptag + std::to_string(gappID);
86  tdName += m_cfg.sentag + std::to_string(gsenID);
87  // create a new directory
88  m_outputFile->mkdir(tdName.c_str());
89  m_outputFile->cd(tdName.c_str());
90 
91  ACTS_VERBOSE("Writing out map at " << tdName);
92 
93  size_t bins0 = 1, bins1 = 1;
94  // understand what sort of material you have in mind
95  const Acts::BinnedSurfaceMaterial* bsm =
96  dynamic_cast<const Acts::BinnedSurfaceMaterial*>(sMaterial);
97  if (bsm != nullptr) {
98  // overwrite the bin numbers
99  bins0 = bsm->binUtility().bins(0);
100  bins1 = bsm->binUtility().bins(1);
101 
102  // Get the binning data
103  auto& binningData = bsm->binUtility().binningData();
104  // 1-D or 2-D maps
105  size_t binningBins = binningData.size();
106 
107  // The bin number information
108  TH1F* n = new TH1F(m_cfg.ntag.c_str(), "bins; bin", binningBins, -0.5,
109  binningBins - 0.5);
110 
111  // The binning value information
112  TH1F* v = new TH1F(m_cfg.vtag.c_str(), "binning values; bin", binningBins,
113  -0.5, binningBins - 0.5);
114 
115  // The binning option information
116  TH1F* o = new TH1F(m_cfg.otag.c_str(), "binning options; bin",
117  binningBins, -0.5, binningBins - 0.5);
118 
119  // The binning option information
120  TH1F* min = new TH1F(m_cfg.mintag.c_str(), "min; bin", binningBins, -0.5,
121  binningBins - 0.5);
122 
123  // The binning option information
124  TH1F* max = new TH1F(m_cfg.maxtag.c_str(), "max; bin", binningBins, -0.5,
125  binningBins - 0.5);
126 
127  // Now fill the histogram content
128  size_t b = 1;
129  for (auto bData : binningData) {
130  // Fill: nbins, value, option, min, max
131  n->SetBinContent(b, int(binningData[b - 1].bins()));
132  v->SetBinContent(b, int(binningData[b - 1].binvalue));
133  o->SetBinContent(b, int(binningData[b - 1].option));
134  min->SetBinContent(b, binningData[b - 1].min);
135  max->SetBinContent(b, binningData[b - 1].max);
136  ++b;
137  }
138  n->Write();
139  v->Write();
140  o->Write();
141  min->Write();
142  max->Write();
143  }
144 
145  TH2F* t = new TH2F(m_cfg.ttag.c_str(), "thickness [mm] ;b0 ;b1", bins0,
146  -0.5, bins0 - 0.5, bins1, -0.5, bins1 - 0.5);
147  TH2F* x0 = new TH2F(m_cfg.x0tag.c_str(), "X_{0} [mm] ;b0 ;b1", bins0, -0.5,
148  bins0 - 0.5, bins1, -0.5, bins1 - 0.5);
149  TH2F* l0 = new TH2F(m_cfg.l0tag.c_str(), "#Lambda_{0} [mm] ;b0 ;b1", bins0,
150  -0.5, bins0 - 0.5, bins1, -0.5, bins1 - 0.5);
151  TH2F* A = new TH2F(m_cfg.atag.c_str(), "X_{0} [mm] ;b0 ;b1", bins0, -0.5,
152  bins0 - 0.5, bins1, -0.5, bins1 - 0.5);
153  TH2F* Z = new TH2F(m_cfg.ztag.c_str(), "#Lambda_{0} [mm] ;b0 ;b1", bins0,
154  -0.5, bins0 - 0.5, bins1, -0.5, bins1 - 0.5);
155  TH2F* rho = new TH2F(m_cfg.rhotag.c_str(), "#rho [g/mm^3] ;b0 ;b1", bins0,
156  -0.5, bins0 - 0.5, bins1, -0.5, bins1 - 0.5);
157 
158  // loop over the material and fill
159  for (size_t b0 = 0; b0 < bins0; ++b0) {
160  for (size_t b1 = 0; b1 < bins1; ++b1) {
161  // get the material for the bin
162  auto& mat = sMaterial->materialSlab(b0, b1);
163  if (mat) {
164  t->SetBinContent(b0 + 1, b1 + 1, mat.thickness());
165  x0->SetBinContent(b0 + 1, b1 + 1, mat.material().X0());
166  l0->SetBinContent(b0 + 1, b1 + 1, mat.material().L0());
167  A->SetBinContent(b0 + 1, b1 + 1, mat.material().Ar());
168  Z->SetBinContent(b0 + 1, b1 + 1, mat.material().Z());
169  rho->SetBinContent(b0 + 1, b1 + 1, mat.material().massDensity());
170  }
171  }
172  }
173  t->Write();
174  x0->Write();
175  l0->Write();
176  A->Write();
177  Z->Write();
178  rho->Write();
179  }
180 
181  auto& volumeMaps = detMaterial.second;
182  for (auto& [key, value] : volumeMaps) {
183  // Get the Volume material
184  const Acts::IVolumeMaterial* vMaterial = value.get();
185 
186  // get the geometry ID
187  Acts::GeometryIdentifier geoID = key;
188  // decode the geometryID
189  const auto gvolID = geoID.volume();
190 
191  // create the directory
192  std::string tdName = m_cfg.folderVolumeNameBase.c_str();
193  tdName += m_cfg.voltag + std::to_string(gvolID);
194 
195  // create a new directory
196  m_outputFile->mkdir(tdName.c_str());
197  m_outputFile->cd(tdName.c_str());
198 
199  ACTS_VERBOSE("Writing out map at " << tdName);
200 
201  // understand what sort of material you have in mind
202  auto bvMaterial3D = dynamic_cast<const Acts::InterpolatedMaterialMap<
203  Acts::MaterialMapper<Acts::MaterialGrid3D>>*>(vMaterial);
204  auto bvMaterial2D = dynamic_cast<const Acts::InterpolatedMaterialMap<
205  Acts::MaterialMapper<Acts::MaterialGrid2D>>*>(vMaterial);
206 
207  size_t points = 1;
208  if (bvMaterial3D != nullptr || bvMaterial2D != nullptr) {
209  // Get the binning data
210  std::vector<Acts::BinningData> binningData;
211  if (bvMaterial3D != nullptr) {
212  binningData = bvMaterial3D->binUtility().binningData();
213  Acts::MaterialGrid3D grid = bvMaterial3D->getMapper().getGrid();
214  points = grid.size();
215  } else {
216  binningData = bvMaterial2D->binUtility().binningData();
217  Acts::MaterialGrid2D grid = bvMaterial2D->getMapper().getGrid();
218  points = grid.size();
219  }
220 
221  // 2-D or 3-D maps
222  size_t binningBins = binningData.size();
223 
224  // The bin number information
225  TH1F* n = new TH1F(m_cfg.ntag.c_str(), "bins; bin", binningBins, -0.5,
226  binningBins - 0.5);
227 
228  // The binning value information
229  TH1F* v = new TH1F(m_cfg.vtag.c_str(), "binning values; bin", binningBins,
230  -0.5, binningBins - 0.5);
231 
232  // The binning option information
233  TH1F* o = new TH1F(m_cfg.otag.c_str(), "binning options; bin",
234  binningBins, -0.5, binningBins - 0.5);
235 
236  // The binning option information
237  TH1F* min = new TH1F(m_cfg.mintag.c_str(), "min; bin", binningBins, -0.5,
238  binningBins - 0.5);
239 
240  // The binning option information
241  TH1F* max = new TH1F(m_cfg.maxtag.c_str(), "max; bin", binningBins, -0.5,
242  binningBins - 0.5);
243 
244  // Now fill the histogram content
245  size_t b = 1;
246  for (auto bData : binningData) {
247  // Fill: nbins, value, option, min, max
248  n->SetBinContent(b, int(binningData[b - 1].bins()));
249  v->SetBinContent(b, int(binningData[b - 1].binvalue));
250  o->SetBinContent(b, int(binningData[b - 1].option));
251  min->SetBinContent(b, binningData[b - 1].min);
252  max->SetBinContent(b, binningData[b - 1].max);
253  ++b;
254  }
255  n->Write();
256  v->Write();
257  o->Write();
258  min->Write();
259  max->Write();
260  }
261 
262  TH1F* x0 = new TH1F(m_cfg.x0tag.c_str(), "X_{0} [mm] ;gridPoint", points,
263  -0.5, points - 0.5);
264  TH1F* l0 = new TH1F(m_cfg.l0tag.c_str(), "#Lambda_{0} [mm] ;gridPoint",
265  points, -0.5, points - 0.5);
266  TH1F* A = new TH1F(m_cfg.atag.c_str(), "X_{0} [mm] ;gridPoint", points,
267  -0.5, points - 0.5);
268  TH1F* Z = new TH1F(m_cfg.ztag.c_str(), "#Lambda_{0} [mm] ;gridPoint",
269  points, -0.5, points - 0.5);
270  TH1F* rho = new TH1F(m_cfg.rhotag.c_str(), "#rho [g/mm^3] ;gridPoint",
271  points, -0.5, points - 0.5);
272  // homogeneous volume
273  if (points == 1) {
274  auto mat = vMaterial->material({0, 0, 0});
275  x0->SetBinContent(1, mat.X0());
276  l0->SetBinContent(1, mat.L0());
277  A->SetBinContent(1, mat.Ar());
278  Z->SetBinContent(1, mat.Z());
279  rho->SetBinContent(1, mat.massDensity());
280  } else {
281  // 3d grid volume
282  if (bvMaterial3D != nullptr) {
283  Acts::MaterialGrid3D grid = bvMaterial3D->getMapper().getGrid();
284  for (size_t point = 0; point < points; point++) {
285  auto mat = Acts::Material(grid.at(point));
286  if (mat) {
287  x0->SetBinContent(point + 1, mat.X0());
288  l0->SetBinContent(point + 1, mat.L0());
289  A->SetBinContent(point + 1, mat.Ar());
290  Z->SetBinContent(point + 1, mat.Z());
291  rho->SetBinContent(point + 1, mat.massDensity());
292  }
293  }
294  }
295  // 2d grid volume
296  else if (bvMaterial2D != nullptr) {
297  Acts::MaterialGrid2D grid = bvMaterial2D->getMapper().getGrid();
298  for (size_t point = 0; point < points; point++) {
299  auto mat = Acts::Material(grid.at(point));
300  if (mat) {
301  x0->SetBinContent(point + 1, mat.X0());
302  l0->SetBinContent(point + 1, mat.L0());
303  A->SetBinContent(point + 1, mat.Ar());
304  Z->SetBinContent(point + 1, mat.Z());
305  rho->SetBinContent(point + 1, mat.massDensity());
306  }
307  }
308  }
309  }
310  x0->Write();
311  l0->Write();
312  A->Write();
313  Z->Write();
314  rho->Write();
315  }
316 }
317 
318 ActsExamples::RootMaterialWriter::~RootMaterialWriter() {
319  if (m_outputFile != nullptr) {
320  m_outputFile->Close();
321  }
322 }
323 
324 void ActsExamples::RootMaterialWriter::write(
325  const Acts::TrackingGeometry& tGeometry) {
326  // Create a detector material map and loop recursively through it
327  Acts::DetectorMaterialMaps detMatMap;
328  auto hVolume = tGeometry.highestTrackingVolume();
329  if (hVolume != nullptr) {
330  collectMaterial(*hVolume, detMatMap);
331  }
332  // Write the resulting map to the file
333  writeMaterial(detMatMap);
334 }
335 
336 void ActsExamples::RootMaterialWriter::collectMaterial(
337  const Acts::TrackingVolume& tVolume,
338  Acts::DetectorMaterialMaps& detMatMap) {
339  // If the volume has volume material, write that
340  if (tVolume.volumeMaterialSharedPtr() != nullptr and m_cfg.processVolumes) {
341  detMatMap.second[tVolume.geometryId()] = tVolume.volumeMaterialSharedPtr();
342  }
343 
344  // If confined layers exist, loop over them and collect the layer material
345  if (tVolume.confinedLayers() != nullptr) {
346  for (auto& lay : tVolume.confinedLayers()->arrayObjects()) {
347  collectMaterial(*lay, detMatMap);
348  }
349  }
350 
351  // If any of the boundary surfaces has material collect that
352  if (m_cfg.processBoundaries) {
353  for (auto& bou : tVolume.boundarySurfaces()) {
354  const auto& bSurface = bou->surfaceRepresentation();
355  if (bSurface.surfaceMaterialSharedPtr() != nullptr) {
356  detMatMap.first[bSurface.geometryId()] =
357  bSurface.surfaceMaterialSharedPtr();
358  }
359  }
360  }
361 
362  // If the volume has sub volumes, step down
363  if (tVolume.confinedVolumes() != nullptr) {
364  for (auto& tvol : tVolume.confinedVolumes()->arrayObjects()) {
365  collectMaterial(*tvol, detMatMap);
366  }
367  }
368 }
369 
370 void ActsExamples::RootMaterialWriter::collectMaterial(
371  const Acts::Layer& tLayer, Acts::DetectorMaterialMaps& detMatMap) {
372  // If the representing surface has material, collect it
373  const auto& rSurface = tLayer.surfaceRepresentation();
374  if (rSurface.surfaceMaterialSharedPtr() != nullptr and
375  m_cfg.processRepresenting) {
376  detMatMap.first[rSurface.geometryId()] =
377  rSurface.surfaceMaterialSharedPtr();
378  }
379 
380  // Check the approach surfaces
381  if (tLayer.approachDescriptor() != nullptr and m_cfg.processApproaches) {
382  for (auto& aSurface : tLayer.approachDescriptor()->containedSurfaces()) {
383  if (aSurface->surfaceMaterialSharedPtr() != nullptr) {
384  detMatMap.first[aSurface->geometryId()] =
385  aSurface->surfaceMaterialSharedPtr();
386  }
387  }
388  }
389 
390  // Check the sensitive surfaces
391  if (tLayer.surfaceArray() != nullptr and m_cfg.processSensitives) {
392  // sensitive surface loop
393  for (auto& sSurface : tLayer.surfaceArray()->surfaces()) {
394  if (sSurface->surfaceMaterialSharedPtr() != nullptr) {
395  detMatMap.first[sSurface->geometryId()] =
396  sSurface->surfaceMaterialSharedPtr();
397  }
398  }
399  }
400 }