EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DigitizationAlgorithm.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DigitizationAlgorithm.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2017 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 
22 #include "Acts/Utilities/Units.hpp"
27 
28 #include <iostream>
29 #include <stdexcept>
30 
33  : ActsExamples::BareAlgorithm("DigitizationAlgorithm", lvl),
34  m_cfg(std::move(cfg)) {
35  if (m_cfg.inputSimulatedHits.empty()) {
36  throw std::invalid_argument("Missing input hits collection");
37  }
38  if (m_cfg.outputClusters.empty()) {
39  throw std::invalid_argument("Missing output clusters collection");
40  }
41  if (not m_cfg.trackingGeometry) {
42  throw std::invalid_argument("Missing tracking geometry");
43  }
45  throw std::invalid_argument("Missing planar module stepper");
46  }
47  if (!m_cfg.randomNumbers) {
48  throw std::invalid_argument("Missing random numbers tool");
49  }
50  // fill the digitizables map to allow lookup by geometry id only
51  m_cfg.trackingGeometry->visitSurfaces([this](const Acts::Surface* surface) {
52  Digitizable dg;
53  // require a valid surface
54  dg.surface = surface;
55  if (not dg.surface) {
56  return;
57  }
58  // require an associated detector element
59  dg.detectorElement = dynamic_cast<const Acts::IdentifiedDetectorElement*>(
60  dg.surface->associatedDetectorElement());
61  if (not dg.detectorElement) {
62  return;
63  }
64  // require an associated digitization module
65  dg.digitizer = dg.detectorElement->digitizationModule().get();
66  if (not dg.digitizer) {
67  return;
68  }
69  // record all valid surfaces
70  this->m_digitizables.insert_or_assign(surface->geometryId(), dg);
71  });
72 }
73 
75  const AlgorithmContext& ctx) const {
76  // Prepare the input and output collections
77  const auto& hits =
78  ctx.eventStore.get<SimHitContainer>(m_cfg.inputSimulatedHits);
80 
81  for (auto&& [moduleGeoId, moduleHits] : groupByModule(hits)) {
82  // can only digitize hits on digitizable surfaces
83  const auto it = m_digitizables.find(moduleGeoId);
84  if (it == m_digitizables.end()) {
85  continue;
86  }
87 
88  const auto& dg = it->second;
89  // local intersection / direction
90  const auto invTransfrom = dg.surface->transform(ctx.geoContext).inverse();
91 
92  // use iterators manually so we can retrieve the hit index in the container
93  for (auto ih = moduleHits.begin(); ih != moduleHits.end(); ++ih) {
94  const auto& hit = *ih;
95  const auto idx = hits.index_of(ih);
96 
97  Acts::Vector2D localIntersect = (invTransfrom * hit.position()).head<2>();
98  Acts::Vector3D localDirection =
99  invTransfrom.linear() * hit.unitDirection();
100 
101  // compute digitization steps
102  const auto thickness = dg.detectorElement->thickness();
103  const auto lorentzAngle = dg.digitizer->lorentzAngle();
104  auto lorentzShift = thickness * std::tan(lorentzAngle);
105  lorentzShift *= -(dg.digitizer->readoutDirection());
106  // now calculate the steps through the silicon
107  std::vector<Acts::DigitizationStep> dSteps =
108  m_cfg.planarModuleStepper->cellSteps(ctx.geoContext, *dg.digitizer,
109  localIntersect, localDirection);
110  // everything under threshold or edge effects
111  if (!dSteps.size()) {
112  ACTS_VERBOSE("No steps returned from stepper.");
113  continue;
114  }
115 
116  // lets create a cluster - centroid method
117  double localX = 0.;
118  double localY = 0.;
119  double totalPath = 0.;
120  // the cells to be used
121  std::vector<Acts::DigitizationCell> usedCells;
122  usedCells.reserve(dSteps.size());
123  // loop over the steps
124  for (auto dStep : dSteps) {
125  // @todo implement smearing
126  localX += dStep.stepLength * dStep.stepCellCenter.x();
127  localY += dStep.stepLength * dStep.stepCellCenter.y();
128  totalPath += dStep.stepLength;
129  usedCells.push_back(Acts::DigitizationCell(dStep.stepCell.channel0,
130  dStep.stepCell.channel1,
131  dStep.stepLength));
132  }
133  // divide by the total path
134  localX /= totalPath;
135  localX += lorentzShift;
136  localY /= totalPath;
137 
138  // get the segmentation & find the corresponding cell id
139  const Acts::Segmentation& segmentation = dg.digitizer->segmentation();
140  auto binUtility = segmentation.binUtility();
141  Acts::Vector2D localPosition(localX, localY);
142  // @todo remove unneccesary conversion
143  // size_t bin0 = binUtility.bin(localPosition, 0);
144  // size_t bin1 = binUtility.bin(localPosition, 1);
145  // size_t binSerialized = binUtility.serialize({{bin0, bin1, 0}});
146 
147  // the covariance is currently set to 0.
149  cov << 0.05, 0., 0., 0., 0.05, 0., 0., 0.,
151 
152  // create the planar cluster
153  Acts::PlanarModuleCluster pCluster(
154  dg.surface->getSharedPtr(), Identifier(identifier_type(idx), {idx}),
155  std::move(cov), localX, localY, hit.time(), std::move(usedCells));
156 
157  // insert into the cluster container. since the input data is already
158  // sorted by geoId, we should always be able to add at the end.
159  clusters.emplace_hint(clusters.end(), hit.geometryId(),
160  std::move(pCluster));
161  }
162  }
163 
164  ACTS_DEBUG("digitized " << hits.size() << " hits into " << clusters.size()
165  << " clusters");
166 
167  // write the clusters to the EventStore
168  ctx.eventStore.add(m_cfg.outputClusters, std::move(clusters));
170 }