EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DoubleHitSpacePointBuilder.ipp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DoubleHitSpacePointBuilder.ipp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 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 
11 #include <cmath>
12 #include <limits>
13 
14 namespace Acts {
15 namespace detail {
34  double qmag = 0.;
36  double m = 0.;
38  double n = 0.;
41  double limit = 1.;
44  double limitExtended = 0.;
45 };
46 
56 double differenceOfClustersChecked(const Vector3D& pos1, const Vector3D& pos2,
57  const Vector3D& posVertex,
58  const double maxDistance,
59  const double maxAngleTheta2,
60  const double maxAnglePhi2) {
61  // Check if measurements are close enough to each other
62  if ((pos1 - pos2).norm() > maxDistance) {
63  return -1.;
64  }
65 
66  // Calculate the angles of the vectors
67  double phi1, theta1, phi2, theta2;
68  phi1 = VectorHelpers::phi(pos1 - posVertex);
69  theta1 = VectorHelpers::theta(pos1 - posVertex);
70  phi2 = VectorHelpers::phi(pos2 - posVertex);
71  theta2 = VectorHelpers::theta(pos2 - posVertex);
72 
73  // Calculate the squared difference between the theta angles
74  double diffTheta2 = (theta1 - theta2) * (theta1 - theta2);
75  if (diffTheta2 > maxAngleTheta2) {
76  return -1.;
77  }
78 
79  // Calculate the squared difference between the phi angles
80  double diffPhi2 = (phi1 - phi2) * (phi1 - phi2);
81  if (diffPhi2 > maxAnglePhi2) {
82  return -1.;
83  }
84 
85  // Return the squared distance between both vector
86  return diffTheta2 + diffPhi2;
87 }
88 
96 std::pair<Vector2D, Vector2D> findLocalTopAndBottomEnd(
97  const Vector2D& local, const CartesianSegmentation* segment) {
98  auto& binData = segment->binUtility().binningData();
99  auto& boundariesX = binData[0].boundaries();
100  auto& boundariesY = binData[1].boundaries();
101 
102  // Search the x-/y-bin of the cluster
103  size_t binX = binData[0].searchLocal(local);
104  size_t binY = binData[1].searchLocal(local);
105 
106  // Storage of the local top (first) and bottom (second) end
107  std::pair<Vector2D, Vector2D> topBottomLocal;
108 
109  if (boundariesX[binX + 1] - boundariesX[binX] <
110  boundariesY[binY + 1] - boundariesY[binY]) {
111  // Set the top and bottom end of the strip in local coordinates
112  topBottomLocal.first = {(boundariesX[binX] + boundariesX[binX + 1]) / 2,
113  boundariesY[binY + 1]};
114  topBottomLocal.second = {(boundariesX[binX] + boundariesX[binX + 1]) / 2,
115  boundariesY[binY]};
116  } else {
117  // Set the top and bottom end of the strip in local coordinates
118  topBottomLocal.first = {boundariesX[binX],
119  (boundariesY[binY] + boundariesY[binY + 1]) / 2};
120  topBottomLocal.second = {boundariesX[binX + 1],
121  (boundariesY[binY] + boundariesY[binY + 1]) / 2};
122  }
123  return topBottomLocal;
124 }
125 
137  const Vector3D& q, const Vector3D& r) {
148  Vector3D ac = c - a;
149  double qr = q.dot(r);
150  double denom = q.dot(q) - qr * qr;
151 
152  // Check for numerical stability
153  if (fabs(denom) > 10e-7) {
154  // Return lambda0
155  return (ac.dot(r) * qr - ac.dot(q) * r.dot(r)) / denom;
156  }
157  // lambda0 is in the interval [-1,0]. This return serves as error check.
158  return 1.;
159 }
160 
172  double stripLengthGapTolerance) {
174  // Check if the limits are allowed to be increased
175  if (stripLengthGapTolerance <= 0.) {
176  return false;
177  }
178  spaPoPa.qmag = spaPoPa.q.norm();
179  // Increase the limits. This allows a check if the point is just slightly
180  // outside the SDE
181  spaPoPa.limitExtended =
182  spaPoPa.limit + stripLengthGapTolerance / spaPoPa.qmag;
183  // Check if m is just slightly outside
184  if (fabs(spaPoPa.m) > spaPoPa.limitExtended) {
185  return false;
186  }
187  // Calculate n if not performed previously
188  if (spaPoPa.n == 0.) {
189  spaPoPa.n = -spaPoPa.t.dot(spaPoPa.qs) / spaPoPa.r.dot(spaPoPa.qs);
190  }
191  // Check if n is just slightly outside
192  if (fabs(spaPoPa.n) > spaPoPa.limitExtended) {
193  return false;
194  }
195 
213 
214  // Calculate the scaling factor to project lengths of the second SDE on the
215  // first SDE
216  double secOnFirstScale =
217  spaPoPa.q.dot(spaPoPa.r) / (spaPoPa.qmag * spaPoPa.qmag);
218  // Check if both overshoots are in the same direction
219  if (spaPoPa.m > 1. && spaPoPa.n > 1.) {
220  // Calculate the overshoots
221  double mOvershoot = spaPoPa.m - 1.;
222  double nOvershoot =
223  (spaPoPa.n - 1.) * secOnFirstScale; // Perform projection
224  // Resolve worse overshoot
225  double biggerOvershoot = std::max(mOvershoot, nOvershoot);
226  // Move m and n towards 0
227  spaPoPa.m -= biggerOvershoot;
228  spaPoPa.n -= (biggerOvershoot / secOnFirstScale);
229  // Check if this recovered the space point
230  return fabs(spaPoPa.m) < spaPoPa.limit && fabs(spaPoPa.n) < spaPoPa.limit;
231  }
232  // Check if both overshoots are in the same direction
233  if (spaPoPa.m < -1. && spaPoPa.n < -1.) {
234  // Calculate the overshoots
235  double mOvershoot = -(spaPoPa.m + 1.);
236  double nOvershoot =
237  -(spaPoPa.n + 1.) * secOnFirstScale; // Perform projection
238  // Resolve worse overshoot
239  double biggerOvershoot = std::max(mOvershoot, nOvershoot);
240  // Move m and n towards 0
241  spaPoPa.m += biggerOvershoot;
242  spaPoPa.n += (biggerOvershoot / secOnFirstScale);
243  // Check if this recovered the space point
244  return fabs(spaPoPa.m) < spaPoPa.limit && fabs(spaPoPa.n) < spaPoPa.limit;
245  }
246  // No solution could be found
247  return false;
248 }
249 
263 bool calculateSpacePoint(const std::pair<Vector3D, Vector3D>& stripEnds1,
264  const std::pair<Vector3D, Vector3D>& stripEnds2,
265  const Vector3D& posVertex,
266  SpacePointParameters& spaPoPa,
267  const double stripLengthTolerance) {
288 
289  spaPoPa.s = stripEnds1.first + stripEnds1.second - 2 * posVertex;
290  spaPoPa.t = stripEnds2.first + stripEnds2.second - 2 * posVertex;
291  spaPoPa.qs = spaPoPa.q.cross(spaPoPa.s);
292  spaPoPa.rt = spaPoPa.r.cross(spaPoPa.t);
293  spaPoPa.m = -spaPoPa.s.dot(spaPoPa.rt) / spaPoPa.q.dot(spaPoPa.rt);
294 
295  // Set the limit for the parameter
296  if (spaPoPa.limit == 1. && stripLengthTolerance != 0.) {
297  spaPoPa.limit = 1. + stripLengthTolerance;
298  }
299 
300  // Check if m and n can be resolved in the interval (-1, 1)
301  return (fabs(spaPoPa.m) <= spaPoPa.limit &&
302  fabs(spaPoPa.n = -spaPoPa.t.dot(spaPoPa.qs) /
303  spaPoPa.r.dot(spaPoPa.qs)) <= spaPoPa.limit);
304 }
305 } // namespace detail
306 } // namespace Acts
307 
311 template <typename Cluster>
313  DoubleHitSpacePointConfig cfg)
314  : m_cfg(std::move(cfg)) {}
315 
316 template <typename Cluster>
318  const Cluster& cluster) const {
319  // Local position information
320  auto par = cluster.parameters();
323  return local;
324 }
325 
326 template <typename Cluster>
328  const GeometryContext& gctx, const Cluster& cluster) const {
329  // Receive corresponding surface
330  auto& clusterSurface = cluster.referenceObject();
331 
332  // Transform local into global position information
333  Vector3D mom(1., 1., 1.);
334  return clusterSurface.localToGlobal(gctx, localCoords(cluster), mom);
335 }
336 
337 template <typename Cluster>
339  const GeometryContext& gctx,
340  const std::vector<const Cluster*>& clustersFront,
341  const std::vector<const Cluster*>& clustersBack,
342  std::vector<std::pair<const Cluster*, const Cluster*>>& clusterPairs)
343  const {
344  // Return if no clusters are given in a vector
345  if (clustersFront.empty() || clustersBack.empty()) {
346  return;
347  }
348 
349  // Declare helper variables
350  double currentDiff;
351  double diffMin;
352  unsigned int clusterMinDist;
353 
354  // Walk through all clusters on both surfaces
355  for (unsigned int iClustersFront = 0; iClustersFront < clustersFront.size();
356  iClustersFront++) {
357  // Set the closest distance to the maximum of double
359  // Set the corresponding index to an element not in the list of clusters
360  clusterMinDist = clustersBack.size();
361  for (unsigned int iClustersBack = 0; iClustersBack < clustersBack.size();
362  iClustersBack++) {
363  // Calculate the distances between the hits
365  globalCoords(gctx, *(clustersFront[iClustersFront])),
366  globalCoords(gctx, *(clustersBack[iClustersBack])), m_cfg.vertex,
367  m_cfg.diffDist, m_cfg.diffPhi2, m_cfg.diffTheta2);
368  // Store the closest clusters (distance and index) calculated so far
369  if (currentDiff < diffMin && currentDiff >= 0.) {
370  diffMin = currentDiff;
371  clusterMinDist = iClustersBack;
372  }
373  }
374 
375  // Store the best (=closest) result
376  if (clusterMinDist < clustersBack.size()) {
377  std::pair<const Cluster*, const Cluster*> clusterPair;
378  clusterPair = std::make_pair(clustersFront[iClustersFront],
379  clustersBack[clusterMinDist]);
380  clusterPairs.push_back(clusterPair);
381  }
382  }
383 }
384 
385 template <typename Cluster>
386 std::pair<Acts::Vector3D, Acts::Vector3D>
388  const GeometryContext& gctx, const Cluster& cluster) const {
389  // Calculate the local coordinates of the cluster
390  const Acts::Vector2D local = localCoords(cluster);
391 
392  // Receive the binning
393  auto segment = dynamic_cast<const Acts::CartesianSegmentation*>(
394  &(cluster.digitizationModule()->segmentation()));
395 
396  std::pair<Vector2D, Vector2D> topBottomLocal =
397  detail::findLocalTopAndBottomEnd(local, segment);
398 
399  // Calculate the global coordinates of the top and bottom end of the strip
400  Acts::Vector3D mom(1., 1., 1); // mom is a dummy variable
401  const auto* sur = &cluster.referenceObject();
402  Acts::Vector3D topGlobal =
403  sur->localToGlobal(gctx, topBottomLocal.first, mom);
404  Acts::Vector3D bottomGlobal =
405  sur->localToGlobal(gctx, topBottomLocal.second, mom);
406 
407  // Return the top and bottom end of the strip in global coordinates
408  return std::make_pair(topGlobal, bottomGlobal);
409 }
410 
411 template <typename Cluster>
412 void Acts::SpacePointBuilder<Acts::SpacePoint<Cluster>>::calculateSpacePoints(
413  const GeometryContext& gctx,
414  const std::vector<std::pair<const Cluster*, const Cluster*>>& clusterPairs,
415  std::vector<Acts::SpacePoint<Cluster>>& spacePoints) const {
417 
418  detail::SpacePointParameters spaPoPa;
419 
420  // Walk over every found candidate pair
421  for (const auto& cp : clusterPairs) {
422  // Calculate the ends of the SDEs
423  const auto& ends1 = endsOfStrip(gctx, *(cp.first));
424  const auto& ends2 = endsOfStrip(gctx, *(cp.second));
425 
426  spaPoPa.q = ends1.first - ends1.second;
427  spaPoPa.r = ends2.first - ends2.second;
428 
429  // Fast skipping if a perpendicular projection should be used
430  double resultPerpProj;
431  if (m_cfg.usePerpProj) {
432  resultPerpProj = detail::calcPerpendicularProjection(
433  ends1.first, ends2.first, spaPoPa.q, spaPoPa.r);
434  if (resultPerpProj <= 0.) {
436  sp.clusterModule.push_back(cp.first);
437  sp.clusterModule.push_back(cp.second);
438  Vector3D pos = ends1.first + resultPerpProj * spaPoPa.q;
439  sp.vector = pos;
440  spacePoints.push_back(std::move(sp));
441  continue;
442  }
443  }
444 
445  if (calculateSpacePoint(ends1, ends2, m_cfg.vertex, spaPoPa,
446  m_cfg.stripLengthTolerance)) {
447  // Store the space point
449  sp.clusterModule.push_back(cp.first);
450  sp.clusterModule.push_back(cp.second);
451  Vector3D pos = 0.5 * (ends1.first + ends1.second + spaPoPa.m * spaPoPa.q);
452  // TODO: Clusters should deliver timestamp
453  sp.vector = pos;
454  spacePoints.push_back(std::move(sp));
455  } else {
462  // Check if a recovery the point(s) and store them if successful
463  if (detail::recoverSpacePoint(spaPoPa, m_cfg.stripLengthGapTolerance)) {
465  sp.clusterModule.push_back(cp.first);
466  sp.clusterModule.push_back(cp.second);
467  Vector3D pos =
468  0.5 * (ends1.first + ends1.second + spaPoPa.m * spaPoPa.q);
469  sp.vector = pos;
470  spacePoints.push_back(std::move(sp));
471  }
472  }
473  }
474 }