ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CsvPlanarClusterReader.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CsvPlanarClusterReader.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2019 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 
11 #include <dfe/dfe_io_dsv.hpp>
12 
23 #include "Acts/Utilities/Units.hpp"
24 #include "TrackMlData.hpp"
25 
28  : m_cfg(cfg)
29  // TODO check that all files (hits,cells,truth) exists
30  ,
31  m_eventsRange(determineEventFilesRange(cfg.inputDir, "hits.csv")),
32  m_logger(Acts::getDefaultLogger("CsvPlanarClusterReader", lvl)) {
33  if (m_cfg.outputClusters.empty()) {
34  throw std::invalid_argument("Missing cluster output collection");
35  }
36  if (m_cfg.outputHitIds.empty()) {
37  throw std::invalid_argument("Missing hit id output collection");
38  }
39  if (m_cfg.outputHitParticlesMap.empty()) {
40  throw std::invalid_argument("Missing hit-particles map output collection");
41  }
42  if (m_cfg.outputSimulatedHits.empty()) {
43  throw std::invalid_argument("Missing simulated hits output collection");
44  }
45  if (not m_cfg.trackingGeometry) {
46  throw std::invalid_argument("Missing tracking geometry");
47  }
48  // fill the geo id to surface map once to speed up lookups later on
49  m_cfg.trackingGeometry->visitSurfaces([this](const Acts::Surface* surface) {
50  this->m_surfaces[surface->geoID()] = surface;
51  });
52 }
53 
55  return "CsvPlanarClusterReader";
56 }
57 
58 std::pair<size_t, size_t> FW::CsvPlanarClusterReader::availableEvents() const {
59  return m_eventsRange;
60 }
61 
62 namespace {
63 struct CompareHitId {
64  // support transparent comparision between identifiers and full objects
65  using is_transparent = void;
66  template <typename T>
67  constexpr bool operator()(const T& left, const T& right) const {
68  return left.hit_id < right.hit_id;
69  }
70  template <typename T>
71  constexpr bool operator()(uint64_t left_id, const T& right) const {
72  return left_id < right.hit_id;
73  }
74  template <typename T>
75  constexpr bool operator()(const T& left, uint64_t right_id) const {
76  return left.hit_id < right_id;
77  }
78 };
79 
81 inline Acts::GeometryID extractGeometryId(const FW::HitData& data) {
82  // if available, use the encoded geometry directly
83  if (data.geometry_id != 0u) {
84  return data.geometry_id;
85  }
86  // otherwise, reconstruct it from the available components
87  Acts::GeometryID geoId;
88  geoId.setVolume(data.volume_id);
89  geoId.setLayer(data.layer_id);
90  geoId.setSensitive(data.module_id);
91  return geoId;
92 }
93 
94 struct CompareGeometryId {
95  bool operator()(const FW::HitData& left, const FW::HitData& right) const {
96  auto leftId = extractGeometryId(left).value();
97  auto rightId = extractGeometryId(right).value();
98  return leftId < rightId;
99  }
100 };
101 
102 template <typename Data>
103 inline std::vector<Data> readEverything(
104  const std::string& inputDir, const std::string& filename,
105  const std::vector<std::string>& optionalColumns, size_t event) {
106  std::string path = FW::perEventFilepath(inputDir, filename, event);
107  dfe::NamedTupleCsvReader<Data> reader(path, optionalColumns);
108 
109  std::vector<Data> everything;
110  Data one;
111  while (reader.read(one)) {
112  everything.push_back(one);
113  }
114 
115  return everything;
116 }
117 
118 std::vector<FW::HitData> readHitsByGeoId(const std::string& inputDir,
119  size_t event) {
120  // geometry_id and t are optional columns
121  auto hits = readEverything<FW::HitData>(inputDir, "hits.csv",
122  {"geometry_id", "t"}, event);
123  // sort same way they will be sorted in the output container
124  std::sort(hits.begin(), hits.end(), CompareGeometryId{});
125  return hits;
126 }
127 
128 std::vector<FW::CellData> readCellsByHitId(const std::string& inputDir,
129  size_t event) {
130  // timestamp is an optional element
131  auto cells =
132  readEverything<FW::CellData>(inputDir, "cells.csv", {"timestamp"}, event);
133  // sort for fast hit id look up
134  std::sort(cells.begin(), cells.end(), CompareHitId{});
135  return cells;
136 }
137 
138 std::vector<FW::TruthHitData> readTruthHitsByHitId(const std::string& inputDir,
139  size_t event) {
140  // define all optional columns
141  std::vector<std::string> optionalColumns = {
142  "geometry_id", "tt", "te", "deltapx",
143  "deltapy", "deltapz", "deltae", "index",
144  };
145  auto truths = readEverything<FW::TruthHitData>(inputDir, "truth.csv",
146  optionalColumns, event);
147  // sort for fast hit id look up
148  std::sort(truths.begin(), truths.end(), CompareHitId{});
149  return truths;
150 }
151 
152 } // namespace
153 
155  const FW::AlgorithmContext& ctx) {
156  // hit_id in the files is not required to be neither continuous nor
157  // monotonic. internally, we want continous indices within [0,#hits)
158  // to simplify data handling. to be able to perform this mapping we first
159  // read all data into memory before converting to the internal event data
160  // types.
161  auto hits = readHitsByGeoId(m_cfg.inputDir, ctx.eventNumber);
162  auto cells = readCellsByHitId(m_cfg.inputDir, ctx.eventNumber);
163  auto truths = readTruthHitsByHitId(m_cfg.inputDir, ctx.eventNumber);
164 
165  // prepare containers for the hit data using the framework event data types
167  std::vector<uint64_t> hitIds;
168  IndexMultimap<ActsFatras::Barcode> hitParticlesMap;
169  SimHitContainer simHits;
170  clusters.reserve(hits.size());
171  hitIds.reserve(hits.size());
172  hitParticlesMap.reserve(truths.size());
173  simHits.reserve(truths.size());
174 
175  for (const HitData& hit : hits) {
176  Acts::GeometryID geoId = extractGeometryId(hit);
177 
178  // find associated truth/ simulation hits
179  std::vector<std::size_t> simHitIndices;
180  {
181  auto range = makeRange(std::equal_range(truths.begin(), truths.end(),
182  hit.hit_id, CompareHitId{}));
183  simHitIndices.reserve(range.size());
184  for (const auto& truth : range) {
185  const auto simGeometryId = Acts::GeometryID(truth.geometry_id);
186  // TODO validate geo id consistency
187  const auto simParticleId = ActsFatras::Barcode(truth.particle_id);
188  const auto simIndex = truth.index;
189  ActsFatras::Hit::Vector4 simPos4{
190  truth.tx * Acts::UnitConstants::mm,
191  truth.ty * Acts::UnitConstants::mm,
192  truth.tz * Acts::UnitConstants::mm,
193  truth.tt * Acts::UnitConstants::ns,
194  };
195  ActsFatras::Hit::Vector4 simMom4{
196  truth.tpx * Acts::UnitConstants::GeV,
197  truth.tpy * Acts::UnitConstants::GeV,
198  truth.tpz * Acts::UnitConstants::GeV,
199  truth.te * Acts::UnitConstants::GeV,
200  };
201  ActsFatras::Hit::Vector4 simDelta4{
202  truth.deltapx * Acts::UnitConstants::GeV,
203  truth.deltapy * Acts::UnitConstants::GeV,
204  truth.deltapz * Acts::UnitConstants::GeV,
205  truth.deltae * Acts::UnitConstants::GeV,
206  };
207 
208  // the cluster stores indices to the underlying simulation hits. thus
209  // their position in the container must be stable. the preordering of
210  // hits by geometry id should ensure that new sim hits are always added
211  // at the end and previously created ones rest at their existing
212  // locations.
213  auto inserted = simHits.emplace_hint(simHits.end(), simGeometryId,
214  simParticleId, simPos4, simMom4,
215  simMom4 + simDelta4, simIndex);
216  if (std::next(inserted) != simHits.end()) {
217  ACTS_FATAL("Truth hit sorting broke for input hit id " << hit.hit_id);
218  return ProcessCode::ABORT;
219  }
220  simHitIndices.push_back(simHits.index_of(inserted));
221  }
222  }
223 
224  // find matching pixel cell information
225  std::vector<Acts::DigitizationCell> digitizationCells;
226  {
227  auto range = makeRange(std::equal_range(cells.begin(), cells.end(),
228  hit.hit_id, CompareHitId{}));
229  for (const auto& c : range) {
230  digitizationCells.emplace_back(c.ch0, c.ch1, c.value);
231  }
232  }
233 
234  // identify hit surface
235  auto it = m_surfaces.find(geoId);
236  if (it == m_surfaces.end() or not it->second) {
237  ACTS_FATAL("Could not retrieve the surface for hit " << hit);
238  return ProcessCode::ABORT;
239  }
240  const Acts::Surface& surface = *(it->second);
241 
242  // transform global hit coordinates into local coordinates on the surface
244  hit.y * Acts::UnitConstants::mm,
245  hit.z * Acts::UnitConstants::mm);
246  double time = hit.t * Acts::UnitConstants::ns;
247  Acts::Vector3D mom(1, 1, 1); // fake momentum
248  Acts::Vector2D local(0, 0);
249  surface.globalToLocal(ctx.geoContext, pos, mom, local);
250  // TODO what to use as cluster uncertainty?
252  // create the planar cluster
254  surface.getSharedPtr(),
255  Identifier(identifier_type(geoId.value()), std::move(simHitIndices)),
256  std::move(cov), local[0], local[1], time, std::move(digitizationCells));
257 
258  // due to the previous sorting of the raw hit data by geometry id, new
259  // clusters should always end up at the end of the container. previous
260  // elements were not touched; cluster indices remain stable and can
261  // be used to identify the hit.
262  auto inserted =
263  clusters.emplace_hint(clusters.end(), geoId, std::move(cluster));
264  if (std::next(inserted) != clusters.end()) {
265  ACTS_FATAL("Something went horribly wrong with the hit sorting");
266  return ProcessCode::ABORT;
267  }
268  auto hitIndex = clusters.index_of(inserted);
269  auto truthRange = makeRange(std::equal_range(truths.begin(), truths.end(),
270  hit.hit_id, CompareHitId{}));
271  for (const auto& truth : truthRange) {
272  hitParticlesMap.emplace_hint(hitParticlesMap.end(), hitIndex,
273  truth.particle_id);
274  }
275 
276  // map internal hit/cluster index back to original, non-monotonic hit id
277  hitIds.push_back(hit.hit_id);
278  }
279 
280  // write the data to the EventStore
281  ctx.eventStore.add(m_cfg.outputClusters, std::move(clusters));
282  ctx.eventStore.add(m_cfg.outputHitIds, std::move(hitIds));
283  ctx.eventStore.add(m_cfg.outputHitParticlesMap, std::move(hitParticlesMap));
284  ctx.eventStore.add(m_cfg.outputSimulatedHits, std::move(simHits));
285 
287 }