Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ModuleClusters.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ModuleClusters.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2021 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 
14 
15 #include <array>
16 #include <cmath>
17 #include <cstdint>
18 #include <cstdlib>
19 #include <memory>
20 #include <stdexcept>
21 #include <type_traits>
22 
23 namespace ActsExamples {
24 
26  ModuleValue mval;
27  mval.paramIndices = std::move(params.indices);
28  mval.paramValues = std::move(params.values);
29  mval.paramVariances = std::move(params.variances);
30  mval.sources = {simhit};
31 
32  if (m_merge and not params.cluster.channels.empty()) {
33  // Break-up the cluster
34  for (auto cell : params.cluster.channels) {
35  ModuleValue mval_cell = mval;
36  mval_cell.value = cell;
37  m_moduleValues.push_back(std::move(mval_cell));
38  }
39  } else {
40  // pass-through mode or smeared indices only
41  mval.value = std::move(params.cluster);
42  m_moduleValues.push_back(std::move(mval));
43  }
44 }
45 
46 std::vector<std::pair<DigitizedParameters, std::set<ModuleClusters::simhit_t>>>
48  if (m_merge) { // (re-)build the clusters
49  merge();
50  }
51  std::vector<std::pair<DigitizedParameters, std::set<simhit_t>>> retv;
52  for (ModuleValue& mval : m_moduleValues) {
53  if (std::holds_alternative<Cluster::Cell>(mval.value)) {
54  // Should never happen! Either the cluster should have
55  // passed-through or the cells merged into clusters in the
56  // merge() step. Treat this as a bug.
57  throw std::runtime_error("Invalid cluster!");
58  }
59  DigitizedParameters dpars;
60  dpars.indices = mval.paramIndices;
61  dpars.values = mval.paramValues;
62  dpars.variances = mval.paramVariances;
63  dpars.cluster = std::get<Cluster>(mval.value);
64  retv.emplace_back(std::move(dpars), mval.sources);
65  }
66  return retv;
67 }
68 
69 // Needed for clusterization
70 int getCellRow(const ModuleValue& mval) {
71  if (std::holds_alternative<ActsExamples::Cluster::Cell>(mval.value)) {
72  return std::get<ActsExamples::Cluster::Cell>(mval.value).bin[0];
73  }
74  throw std::domain_error("ModuleValue does not contain cell!");
75 }
76 
78  if (std::holds_alternative<ActsExamples::Cluster::Cell>(mval.value)) {
79  return std::get<ActsExamples::Cluster::Cell>(mval.value).bin[1];
80  }
81  throw std::domain_error("ModuleValue does not contain cell!");
82 }
83 
85  return mval.label;
86 }
87 
88 void clusterAddCell(std::vector<ModuleValue>& cl, const ModuleValue& ce) {
89  cl.push_back(ce);
90 }
91 
92 std::vector<ModuleValue> ModuleClusters::createCellCollection() {
93  std::vector<ModuleValue> cells;
94  for (ModuleValue& mval : m_moduleValues) {
95  if (std::holds_alternative<Cluster::Cell>(mval.value)) {
96  cells.push_back(mval);
97  }
98  }
99  return cells;
100 }
101 
103  std::vector<ModuleValue> cells = createCellCollection();
104 
105  std::vector<ModuleValue> newVals;
106 
107  if (not cells.empty()) {
108  // Case where we actually have geometric clusters
109  std::vector<std::vector<ModuleValue>> merged =
110  Acts::Ccl::createClusters<std::vector<ModuleValue>,
111  std::vector<std::vector<ModuleValue>>>(
113 
114  for (std::vector<ModuleValue>& cellv : merged) {
115  // At this stage, the cellv vector contains cells that form a
116  // consistent cluster based on a connected component analysis
117  // only. Still have to check if they match based on the other
118  // indices (a good example of this would a for a timing
119  // detector).
120 
121  for (std::vector<ModuleValue>& remerged : mergeParameters(cellv)) {
122  newVals.push_back(squash(remerged));
123  }
124  }
125  m_moduleValues = std::move(newVals);
126  } else {
127  // no geo clusters
128  for (std::vector<ModuleValue>& merged : mergeParameters(m_moduleValues)) {
129  newVals.push_back(squash(merged));
130  }
131  m_moduleValues = std::move(newVals);
132  }
133 }
134 
135 // ATTN: returns vector of index into `indices'
136 std::vector<size_t> ModuleClusters::nonGeoEntries(
137  std::vector<Acts::BoundIndices>& indices) {
138  std::vector<size_t> retv;
139  for (size_t i = 0; i < indices.size(); i++) {
140  auto idx = indices.at(i);
141  if (std::find(m_geoIndices.begin(), m_geoIndices.end(), idx) ==
142  m_geoIndices.end()) {
143  retv.push_back(i);
144  }
145  }
146  return retv;
147 }
148 
149 // Merging based on parameters
150 std::vector<std::vector<ModuleValue>> ModuleClusters::mergeParameters(
151  std::vector<ModuleValue> values) {
152  std::vector<std::vector<ModuleValue>> retv;
153 
154  std::vector<bool> used(values.size(), false);
155  for (size_t i = 0; i < values.size(); i++) {
156  if (used.at(i)) {
157  continue;
158  }
159 
160  retv.emplace_back();
161  std::vector<ModuleValue>& thisvec = retv.back();
162 
163  // Value has not yet been claimed, so claim it
164  thisvec.push_back(std::move(values.at(i)));
165  used.at(i) = true;
166 
167  // Values previously visited by index `i' have already been added
168  // to a cluster or used to seed a new cluster, so start at the
169  // next unseen one
170  for (size_t j = i + 1; j < values.size(); j++) {
171  // Still may have already been used, so check it
172  if (used.at(j)) {
173  continue;
174  }
175 
176  // Now look for a match between current cluster and value `j' Consider
177  // them matched until we find evidence to the contrary. This
178  // way, merging still works when digitization is done by
179  // clusters only
180  bool matched = true;
181 
182  // The cluster to be merged into can have more than one
183  // associated value at this point, so we have to consider them
184  // all
185  for (ModuleValue& thisval : thisvec) {
186  // Loop over non-geometric dimensions
187  for (auto k : nonGeoEntries(thisval.paramIndices)) {
188  Acts::ActsScalar p_i = thisval.paramValues.at(k);
189  Acts::ActsScalar p_j = values.at(j).paramValues.at(k);
190  Acts::ActsScalar v_i = thisval.paramVariances.at(k);
191  Acts::ActsScalar v_j = values.at(j).paramVariances.at(k);
192 
193  Acts::ActsScalar left = 0, right = 0;
194  if (p_i < p_j) {
195  left = p_i + m_nsigma * std::sqrt(v_i);
196  right = p_j - m_nsigma * std::sqrt(v_j);
197  } else {
198  left = p_j + m_nsigma * std::sqrt(v_j);
199  right = p_i - m_nsigma * std::sqrt(v_i);
200  }
201  if (left < right) {
202  // We know these two don't match, so break out of the
203  // dimension loop
204  matched = false;
205  break;
206  }
207  } // Loop over `k' (non-geo dimensions)
208  if (matched) {
209  // The value under consideration matched at least one
210  // associated to the current cluster so no need to keep
211  // checking others in current cluster
212  break;
213  }
214  } // Loop on current cluster
215  if (matched) {
216  // Claim value `j'
217  used.at(j) = true;
218  thisvec.push_back(std::move(values.at(j)));
219  }
220  } // Loop on `j'
221  } // Loop on `i'
222  return retv;
223 }
224 
225 ModuleValue ModuleClusters::squash(std::vector<ModuleValue>& values) {
226  ModuleValue mval;
227  Acts::ActsScalar tot = 0;
228  Acts::ActsScalar tot2 = 0;
229  std::vector<Acts::ActsScalar> weights;
230 
231  // First, start by computing cell weights
232  for (ModuleValue& other : values) {
233  if (std::holds_alternative<Cluster::Cell>(other.value)) {
234  weights.push_back(std::get<Cluster::Cell>(other.value).activation);
235  } else {
236  weights.push_back(1);
237  }
238  tot += weights.back();
239  tot2 += weights.back() * weights.back();
240  }
241 
242  // Now, go over the non-geometric indices
243  for (size_t i = 0; i < values.size(); i++) {
244  ModuleValue& other = values.at(i);
245  for (size_t j = 0; j < other.paramIndices.size(); j++) {
246  auto idx = other.paramIndices.at(j);
247  if (std::find(m_geoIndices.begin(), m_geoIndices.end(), idx) ==
248  m_geoIndices.end()) {
249  if (std::find(mval.paramIndices.begin(), mval.paramIndices.end(),
250  idx) == mval.paramIndices.end()) {
251  mval.paramIndices.push_back(idx);
252  }
253  if (mval.paramValues.size() < (j + 1)) {
254  mval.paramValues.push_back(0);
255  mval.paramVariances.push_back(0);
256  }
257  Acts::ActsScalar f = weights.at(i) / (tot > 0 ? tot : 1);
259  weights.at(i) * weights.at(i) / (tot2 > 0 ? tot2 : 1);
260  mval.paramValues.at(j) += f * other.paramValues.at(j);
261  mval.paramVariances.at(j) += f2 * other.paramVariances.at(j);
262  }
263  }
264  }
265 
266  // Now do the geometric indices
267  Cluster clus;
268 
269  const auto& binningData = m_segmentation.binningData();
270  Acts::Vector2 pos(0., 0.);
271  Acts::Vector2 var(0., 0.);
272 
273  size_t b0min = SIZE_MAX;
274  size_t b0max = 0;
275  size_t b1min = SIZE_MAX;
276  size_t b1max = 0;
277 
278  for (size_t i = 0; i < values.size(); i++) {
279  ModuleValue& other = values.at(i);
280  if (not std::holds_alternative<Cluster::Cell>(other.value)) {
281  continue;
282  }
283 
284  Cluster::Cell ch = std::get<Cluster::Cell>(other.value);
285  auto bin = ch.bin;
286 
287  size_t b0 = bin[0];
288  size_t b1 = bin[1];
289 
290  b0min = std::min(b0min, b0);
291  b0max = std::max(b0max, b0);
292  b1min = std::min(b1min, b1);
293  b1max = std::max(b1max, b1);
294 
295  float p0 = binningData[0].center(b0);
296  float w0 = binningData[0].width(b0);
297  float p1 = binningData[1].center(b1);
298  float w1 = binningData[1].width(b1);
299 
300  pos += Acts::Vector2(weights.at(i) * p0, weights.at(i) * p1);
301  // Assume uniform distribution to compute error
302  // N.B. This will overestimate the variance
303  // but it's better than nothing for now
304  var += Acts::Vector2(weights.at(i) * weights.at(i) * w0 * w0 / 12,
305  weights.at(i) * weights.at(i) * w1 * w1 / 12);
306 
307  clus.channels.push_back(std::move(ch));
308 
309  // Will have the right value at last iteration Do it here to
310  // avoid having bogus values when there are no clusters
311  clus.sizeLoc0 = b0max - b0min + 1;
312  clus.sizeLoc1 = b1max - b1min + 1;
313  }
314 
315  if (tot > 0) {
316  pos /= tot;
317  var /= (tot * tot);
318  }
319 
320  for (auto idx : m_geoIndices) {
321  mval.paramIndices.push_back(idx);
322  mval.paramValues.push_back(pos[idx]);
323  mval.paramVariances.push_back(var[idx]);
324  }
325 
326  mval.value = std::move(clus);
327 
328  // Finally do the hit association
329  for (ModuleValue& other : values) {
330  mval.sources.merge(other.sources);
331  }
332 
333  return mval;
334 }
335 
336 } // namespace ActsExamples