ECCE @ EIC 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 
10 
13 #include <TFile.h>
14 #include <TH2F.h>
15 #include <ios>
16 #include <iostream>
17 #include <stdexcept>
18 
21  : m_cfg(cfg) {
22  // Validate the configuration
23  if (m_cfg.folderNameBase.empty()) {
24  throw std::invalid_argument("Missing folder name base");
25  } else if (m_cfg.fileName.empty()) {
26  throw std::invalid_argument("Missing file name");
27  } else if (!m_cfg.logger) {
28  throw std::invalid_argument("Missing logger");
29  } else if (m_cfg.name.empty()) {
30  throw std::invalid_argument("Missing service name");
31  }
32 }
33 
35  const Acts::DetectorMaterialMaps& detMaterial) {
36  // Setup ROOT I/O
37  TFile* outputFile = TFile::Open(m_cfg.fileName.c_str(), "recreate");
38  if (!outputFile) {
39  throw std::ios_base::failure("Could not open '" + m_cfg.fileName);
40  }
41 
42  // Change to the output file
43  outputFile->cd();
44 
45  auto& surfaceMaps = detMaterial.first;
46  for (auto& [key, value] : surfaceMaps) {
47  // Get the Surface material
48  const Acts::ISurfaceMaterial* sMaterial = value.get();
49 
50  // get the geometry ID
51  Acts::GeometryID geoID = key;
52  // decode the geometryID
53  const auto gvolID = geoID.volume();
54  const auto gbouID = geoID.boundary();
55  const auto glayID = geoID.layer();
56  const auto gappID = geoID.approach();
57  const auto gsenID = geoID.sensitive();
58  // create the directory
59  std::string tdName = m_cfg.folderNameBase.c_str();
60  tdName += m_cfg.voltag + std::to_string(gvolID);
61  tdName += m_cfg.boutag + std::to_string(gbouID);
62  tdName += m_cfg.laytag + std::to_string(glayID);
63  tdName += m_cfg.apptag + std::to_string(gappID);
64  tdName += m_cfg.sentag + std::to_string(gsenID);
65  // create a new directory
66  outputFile->mkdir(tdName.c_str());
67  outputFile->cd(tdName.c_str());
68 
69  ACTS_VERBOSE("Writing out map at " << tdName);
70 
71  size_t bins0 = 1, bins1 = 1;
72  // understand what sort of material you have in mind
73  const Acts::BinnedSurfaceMaterial* bsm =
74  dynamic_cast<const Acts::BinnedSurfaceMaterial*>(sMaterial);
75  if (bsm) {
76  // overwrite the bin numbers
77  bins0 = bsm->binUtility().bins(0);
78  bins1 = bsm->binUtility().bins(1);
79 
80  // Get the binning data
81  auto& binningData = bsm->binUtility().binningData();
82  // 1-D or 2-D maps
83  size_t binningBins = binningData.size();
84 
85  // The bin number information
86  TH1F* n = new TH1F(m_cfg.ntag.c_str(), "bins; bin", binningBins, -0.5,
87  binningBins - 0.5);
88 
89  // The binning value information
90  TH1F* v = new TH1F(m_cfg.vtag.c_str(), "binning values; bin", binningBins,
91  -0.5, binningBins - 0.5);
92 
93  // The binning option information
94  TH1F* o = new TH1F(m_cfg.otag.c_str(), "binning options; bin",
95  binningBins, -0.5, binningBins - 0.5);
96 
97  // The binning option information
98  TH1F* min = new TH1F(m_cfg.mintag.c_str(), "min; bin", binningBins, -0.5,
99  binningBins - 0.5);
100 
101  // The binning option information
102  TH1F* max = new TH1F(m_cfg.maxtag.c_str(), "max; bin", binningBins, -0.5,
103  binningBins - 0.5);
104 
105  // Now fill the histogram content
106  size_t b = 1;
107  for (auto bData : binningData) {
108  // Fill: nbins, value, option, min, max
109  n->SetBinContent(b, int(binningData[b - 1].bins()));
110  v->SetBinContent(b, int(binningData[b - 1].binvalue));
111  o->SetBinContent(b, int(binningData[b - 1].option));
112  min->SetBinContent(b, binningData[b - 1].min);
113  max->SetBinContent(b, binningData[b - 1].max);
114  ++b;
115  }
116  n->Write();
117  v->Write();
118  o->Write();
119  min->Write();
120  max->Write();
121  }
122 
123  TH2F* t = new TH2F(m_cfg.ttag.c_str(), "thickness [mm] ;b0 ;b1", bins0,
124  -0.5, bins0 - 0.5, bins1, -0.5, bins1 - 0.5);
125  TH2F* x0 = new TH2F(m_cfg.x0tag.c_str(), "X_{0} [mm] ;b0 ;b1", bins0, -0.5,
126  bins0 - 0.5, bins1, -0.5, bins1 - 0.5);
127  TH2F* l0 = new TH2F(m_cfg.l0tag.c_str(), "#Lambda_{0} [mm] ;b0 ;b1", bins0,
128  -0.5, bins0 - 0.5, bins1, -0.5, bins1 - 0.5);
129  TH2F* A = new TH2F(m_cfg.atag.c_str(), "X_{0} [mm] ;b0 ;b1", bins0, -0.5,
130  bins0 - 0.5, bins1, -0.5, bins1 - 0.5);
131  TH2F* Z = new TH2F(m_cfg.ztag.c_str(), "#Lambda_{0} [mm] ;b0 ;b1", bins0,
132  -0.5, bins0 - 0.5, bins1, -0.5, bins1 - 0.5);
133  TH2F* rho = new TH2F(m_cfg.rhotag.c_str(), "#rho [g/mm^3] ;b0 ;b1", bins0,
134  -0.5, bins0 - 0.5, bins1, -0.5, bins1 - 0.5);
135 
136  // loop over the material and fill
137  for (size_t b0 = 0; b0 < bins0; ++b0) {
138  for (size_t b1 = 0; b1 < bins1; ++b1) {
139  // get the material for the bin
140  auto& mat = sMaterial->materialProperties(b0, b1);
141  if (mat) {
142  t->SetBinContent(b0 + 1, b1 + 1, mat.thickness());
143  x0->SetBinContent(b0 + 1, b1 + 1, mat.material().X0());
144  l0->SetBinContent(b0 + 1, b1 + 1, mat.material().L0());
145  A->SetBinContent(b0 + 1, b1 + 1, mat.material().Ar());
146  Z->SetBinContent(b0 + 1, b1 + 1, mat.material().Z());
147  rho->SetBinContent(b0 + 1, b1 + 1, mat.material().massDensity());
148  }
149  }
150  }
151  t->Write();
152  x0->Write();
153  l0->Write();
154  A->Write();
155  Z->Write();
156  rho->Write();
157  }
158 
159  outputFile->Close();
160 }
161 
163  // Create a detector material map and loop recursively through it
164  Acts::DetectorMaterialMaps detMatMap;
165  auto hVolume = tGeometry.highestTrackingVolume();
166  if (hVolume != nullptr) {
167  collectMaterial(*hVolume, detMatMap);
168  }
169  // Write the resulting map to the file
170  write(detMatMap);
171 }
172 
174  const Acts::TrackingVolume& tVolume,
175  Acts::DetectorMaterialMaps& detMatMap) {
176  // If the volume has volume material, write that
177  if (tVolume.volumeMaterialSharedPtr() != nullptr and m_cfg.processVolumes) {
178  detMatMap.second[tVolume.geoID()] = tVolume.volumeMaterialSharedPtr();
179  }
180 
181  // If confined layers exist, loop over them and collect the layer material
182  if (tVolume.confinedLayers() != nullptr) {
183  for (auto& lay : tVolume.confinedLayers()->arrayObjects()) {
184  collectMaterial(*lay, detMatMap);
185  }
186  }
187 
188  // If any of the boundary surfaces has material collect that
189  if (m_cfg.processBoundaries) {
190  for (auto& bou : tVolume.boundarySurfaces()) {
191  const auto& bSurface = bou->surfaceRepresentation();
192  if (bSurface.surfaceMaterialSharedPtr() != nullptr) {
193  detMatMap.first[bSurface.geoID()] = bSurface.surfaceMaterialSharedPtr();
194  }
195  }
196  }
197 
198  // If the volume has sub volumes, step down
199  if (tVolume.confinedVolumes() != nullptr) {
200  for (auto& tvol : tVolume.confinedVolumes()->arrayObjects()) {
201  collectMaterial(*tvol, detMatMap);
202  }
203  }
204 }
205 
207  const Acts::Layer& tLayer, Acts::DetectorMaterialMaps& detMatMap) {
208  // If the representing surface has material, collect it
209  const auto& rSurface = tLayer.surfaceRepresentation();
210  if (rSurface.surfaceMaterialSharedPtr() != nullptr and
211  m_cfg.processRepresenting) {
212  detMatMap.first[rSurface.geoID()] = rSurface.surfaceMaterialSharedPtr();
213  }
214 
215  // Check the approach surfaces
216  if (tLayer.approachDescriptor() != nullptr and m_cfg.processApproaches) {
217  for (auto& aSurface : tLayer.approachDescriptor()->containedSurfaces()) {
218  if (aSurface->surfaceMaterialSharedPtr() != nullptr) {
219  detMatMap.first[aSurface->geoID()] =
220  aSurface->surfaceMaterialSharedPtr();
221  }
222  }
223  }
224 
225  // Check the sensitive surfaces
226  if (tLayer.surfaceArray() != nullptr and m_cfg.processSensitives) {
227  // sensitive surface loop
228  for (auto& sSurface : tLayer.surfaceArray()->surfaces()) {
229  if (sSurface->surfaceMaterialSharedPtr() != nullptr) {
230  detMatMap.first[sSurface->geoID()] =
231  sSurface->surfaceMaterialSharedPtr();
232  }
233  }
234  }
235 }