Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DigitizationConfigurator.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DigitizationConfigurator.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2023 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/Digitization/DigitizationConfigurator.hpp"
10 
11 #include "Acts/Definitions/Algebra.hpp"
12 #include "Acts/Definitions/TrackParametrization.hpp"
13 #include "Acts/Surfaces/AnnulusBounds.hpp"
14 #include "Acts/Surfaces/DiscTrapezoidBounds.hpp"
15 #include "Acts/Surfaces/RadialBounds.hpp"
16 #include "Acts/Surfaces/RectangleBounds.hpp"
17 #include "Acts/Surfaces/Surface.hpp"
18 #include "Acts/Surfaces/SurfaceBounds.hpp"
19 #include "Acts/Surfaces/TrapezoidBounds.hpp"
20 #include "Acts/Utilities/BinUtility.hpp"
21 #include "Acts/Utilities/BinningType.hpp"
22 #include "Acts/Utilities/Zip.hpp"
23 #include "ActsExamples/Digitization/SmearingConfig.hpp"
24 #include "ActsFatras/Digitization/UncorrelatedHitSmearer.hpp"
25 
26 #include <algorithm>
27 #include <cmath>
28 #include <memory>
29 
30 namespace {
34 bool digiConfigMaybeEqual(ActsExamples::DigiComponentsConfig &a,
35  ActsExamples::DigiComponentsConfig &b) {
36  // Check smearing config
37  for (const auto &[as, bs] :
38  Acts::zip(a.smearingDigiConfig, b.smearingDigiConfig)) {
39  if (as.index != bs.index) {
40  return false;
41  }
42  }
43  // Check geometric config
44  const auto &ag = a.geometricDigiConfig;
45  const auto &bg = b.geometricDigiConfig;
46  return (ag.indices == bg.indices and ag.segmentation == bg.segmentation and
47  ag.thickness == bg.thickness and ag.threshold == bg.threshold and
48  ag.digital == bg.digital);
49 }
50 } // namespace
51 
52 void ActsExamples::DigitizationConfigurator::operator()(
53  const Acts::Surface *surface) {
54  if (surface->associatedDetectorElement() != nullptr) {
55  Acts::GeometryIdentifier geoId = surface->geometryId();
56  auto dInputConfig = inputDigiComponents.find((geoId));
57  if (dInputConfig != inputDigiComponents.end()) {
58  // The output config, copy over the smearing part
59  DigiComponentsConfig dOutputConfig;
60  dOutputConfig.smearingDigiConfig = dInputConfig->smearingDigiConfig;
61 
62  if (not dInputConfig->geometricDigiConfig.indices.empty()) {
63  // Copy over what can be done
64  dOutputConfig.geometricDigiConfig.indices =
65  dInputConfig->geometricDigiConfig.indices;
66  dOutputConfig.geometricDigiConfig.thickness =
67  dInputConfig->geometricDigiConfig.thickness;
68  dOutputConfig.geometricDigiConfig.chargeSmearer =
69  dInputConfig->geometricDigiConfig.chargeSmearer;
70  dOutputConfig.geometricDigiConfig.digital =
71  dInputConfig->geometricDigiConfig.digital;
72 
73  const Acts::SurfaceBounds &sBounds = surface->bounds();
74  auto boundValues = sBounds.values();
75 
76  const auto &inputSegmentation =
77  dInputConfig->geometricDigiConfig.segmentation;
78  Acts::BinUtility outputSegmentation;
79 
80  switch (sBounds.type()) {
81  // The module is a rectangle module
82  case Acts::SurfaceBounds::eRectangle: {
83  if (inputSegmentation.binningData()[0].binvalue == Acts::binX) {
84  Acts::ActsScalar minX = boundValues[Acts::RectangleBounds::eMinX];
85  Acts::ActsScalar maxX = boundValues[Acts::RectangleBounds::eMaxX];
86  unsigned int nBins = static_cast<unsigned int>(std::round(
87  (maxX - minX) / inputSegmentation.binningData()[0].step));
88  outputSegmentation +=
89  Acts::BinUtility(nBins, minX, maxX, Acts::open, Acts::binX);
90  }
91  if (inputSegmentation.binningData()[0].binvalue == Acts::binY or
92  inputSegmentation.dimensions() == 2) {
93  unsigned int accessBin =
94  inputSegmentation.dimensions() == 2 ? 1 : 0;
95  Acts::ActsScalar minY = boundValues[Acts::RectangleBounds::eMinY];
96  Acts::ActsScalar maxY = boundValues[Acts::RectangleBounds::eMaxY];
97  unsigned int nBins = static_cast<unsigned int>(
98  std::round((maxY - minY) /
99  inputSegmentation.binningData()[accessBin].step));
100  outputSegmentation +=
101  Acts::BinUtility(nBins, minY, maxY, Acts::open, Acts::binY);
102  }
103  } break;
104 
105  // The module is a trapezoid module
106  case Acts::SurfaceBounds::eTrapezoid: {
107  if (inputSegmentation.binningData()[0].binvalue == Acts::binX) {
108  Acts::ActsScalar maxX = std::max(
109  boundValues[Acts::TrapezoidBounds::eHalfLengthXnegY],
110  boundValues[Acts::TrapezoidBounds::eHalfLengthXposY]);
111  unsigned int nBins = static_cast<unsigned int>(std::round(
112  2 * maxX / inputSegmentation.binningData()[0].step));
113  outputSegmentation +=
114  Acts::BinUtility(nBins, -maxX, maxX, Acts::open, Acts::binX);
115  }
116  if (inputSegmentation.binningData()[0].binvalue == Acts::binY or
117  inputSegmentation.dimensions() == 2) {
118  unsigned int accessBin =
119  inputSegmentation.dimensions() == 2 ? 1 : 0;
120  Acts::ActsScalar maxY =
121  boundValues[Acts::TrapezoidBounds::eHalfLengthY];
122  unsigned int nBins = static_cast<unsigned int>(
123  std::round((2 * maxY) /
124  inputSegmentation.binningData()[accessBin].step));
125  outputSegmentation +=
126  Acts::BinUtility(nBins, -maxY, maxY, Acts::open, Acts::binY);
127  }
128  } break;
129 
130  // The module is an annulus module
131  case Acts::SurfaceBounds::eAnnulus: {
132  if (inputSegmentation.binningData()[0].binvalue == Acts::binR) {
133  Acts::ActsScalar minR = boundValues[Acts::AnnulusBounds::eMinR];
134  Acts::ActsScalar maxR = boundValues[Acts::AnnulusBounds::eMaxR];
135  unsigned int nBins = static_cast<unsigned int>(std::round(
136  (maxR - minR) / inputSegmentation.binningData()[0].step));
137  outputSegmentation +=
138  Acts::BinUtility(nBins, minR, maxR, Acts::open, Acts::binR);
139  }
140  if (inputSegmentation.binningData()[0].binvalue == Acts::binPhi or
141  inputSegmentation.dimensions() == 2) {
142  unsigned int accessBin =
143  inputSegmentation.dimensions() == 2 ? 1 : 0;
144  Acts::ActsScalar averagePhi =
145  boundValues[Acts::AnnulusBounds::eAveragePhi];
146  Acts::ActsScalar minPhi =
147  averagePhi - boundValues[Acts::AnnulusBounds::eMinPhiRel];
148  Acts::ActsScalar maxPhi =
149  averagePhi + boundValues[Acts::AnnulusBounds::eMaxPhiRel];
150  unsigned int nBins = static_cast<unsigned int>(
151  std::round((maxPhi - minPhi) /
152  inputSegmentation.binningData()[accessBin].step));
153  outputSegmentation += Acts::BinUtility(nBins, minPhi, maxPhi,
154  Acts::open, Acts::binPhi);
155  }
156 
157  } break;
158 
159  // The module is a Disc Trapezoid
160  case Acts::SurfaceBounds::eDiscTrapezoid: {
161  Acts::ActsScalar minR =
162  boundValues[Acts::DiscTrapezoidBounds::eMinR];
163  Acts::ActsScalar maxR =
164  boundValues[Acts::DiscTrapezoidBounds::eMaxR];
165 
166  if (inputSegmentation.binningData()[0].binvalue == Acts::binR) {
167  unsigned int nBins = static_cast<unsigned int>(std::round(
168  (maxR - minR) / inputSegmentation.binningData()[0].step));
169  outputSegmentation +=
170  Acts::BinUtility(nBins, minR, maxR, Acts::open, Acts::binR);
171  }
172  if (inputSegmentation.binningData()[0].binvalue == Acts::binPhi or
173  inputSegmentation.dimensions() == 2) {
174  unsigned int accessBin =
175  inputSegmentation.dimensions() == 2 ? 1 : 0;
176  Acts::ActsScalar hxMinR =
177  boundValues[Acts::DiscTrapezoidBounds::eHalfLengthXminR];
178  Acts::ActsScalar hxMaxR =
179  boundValues[Acts::DiscTrapezoidBounds::eHalfLengthXmaxR];
180 
181  Acts::ActsScalar averagePhi =
182  boundValues[Acts::DiscTrapezoidBounds::eAveragePhi];
183  Acts::ActsScalar alphaMinR = std::atan2(minR, hxMinR);
184  Acts::ActsScalar alphaMaxR = std::atan2(maxR, hxMaxR);
185  Acts::ActsScalar alpha = std::max(alphaMinR, alphaMaxR);
186  unsigned int nBins = static_cast<unsigned int>(std::round(
187  2 * alpha / inputSegmentation.binningData()[accessBin].step));
188  outputSegmentation += Acts::BinUtility(nBins, averagePhi - alpha,
189  averagePhi + alpha,
190  Acts::open, Acts::binPhi);
191  }
192 
193  } break;
194 
195  case Acts::SurfaceBounds::eDisc: {
196  if (inputSegmentation.binningData()[0].binvalue == Acts::binR) {
197  Acts::ActsScalar minR = boundValues[Acts::RadialBounds::eMinR];
198  Acts::ActsScalar maxR = boundValues[Acts::RadialBounds::eMaxR];
199  unsigned int nBins = static_cast<unsigned int>(std::round(
200  (maxR - minR) / inputSegmentation.binningData()[0].step));
201  outputSegmentation +=
202  Acts::BinUtility(nBins, minR, maxR, Acts::open, Acts::binR);
203  }
204  if (inputSegmentation.binningData()[0].binvalue == Acts::binPhi or
205  inputSegmentation.dimensions() == 2) {
206  unsigned int accessBin =
207  inputSegmentation.dimensions() == 2 ? 1 : 0;
208 
209  Acts::ActsScalar averagePhi =
210  boundValues[Acts::RadialBounds::eAveragePhi];
211  Acts::ActsScalar halfPhiSector =
212  boundValues[Acts::RadialBounds::eHalfPhiSector];
213  Acts::ActsScalar minPhi = averagePhi - halfPhiSector;
214  Acts::ActsScalar maxPhi = averagePhi + halfPhiSector;
215 
216  unsigned int nBins = static_cast<unsigned int>(
217  std::round((maxPhi - minPhi) /
218  inputSegmentation.binningData()[accessBin].step));
219  outputSegmentation += Acts::BinUtility(nBins, minPhi, maxPhi,
220  Acts::open, Acts::binPhi);
221  }
222 
223  } break;
224 
225  default:
226  break;
227  }
228  // Set the adapted segmentation class
229  dOutputConfig.geometricDigiConfig.segmentation = outputSegmentation;
230  }
231 
232  // Compactify the output map where possible
233  if (compactify) {
234  // Check for a representing volume configuration, insert if not
235  // present
236  Acts::GeometryIdentifier volGeoId =
237  Acts::GeometryIdentifier().setVolume(geoId.volume());
238 
239  auto volRep = volumeLayerComponents.find(volGeoId);
240  if (volRep != volumeLayerComponents.end() and
241  digiConfigMaybeEqual(dOutputConfig, volRep->second)) {
242  // return if the volume representation already covers this one
243  return;
244  } else {
245  volumeLayerComponents[volGeoId] = dOutputConfig;
246  outputDigiComponents.push_back({volGeoId, dOutputConfig});
247  }
248 
249  // Check for a representing layer configuration, insert if not present
250  Acts::GeometryIdentifier volLayGeoId =
251  Acts::GeometryIdentifier(volGeoId).setLayer(geoId.layer());
252  auto volLayRep = volumeLayerComponents.find(volLayGeoId);
253 
254  if (volLayRep != volumeLayerComponents.end() and
255  digiConfigMaybeEqual(dOutputConfig, volLayRep->second)) {
256  return;
257  } else {
258  volumeLayerComponents[volLayGeoId] = dOutputConfig;
259  outputDigiComponents.push_back({volLayGeoId, dOutputConfig});
260  }
261  }
262 
263  // Insert into the output list
264  outputDigiComponents.push_back({geoId, dOutputConfig});
265  }
266  }
267 }