ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
VerticesHelper.hpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file VerticesHelper.hpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2020 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 #pragma once
10 
11 #include <utility>
12 #include <vector>
14 
15 namespace Acts {
16 
17 namespace detail {
18 
21 namespace VerticesHelper {
22 
32 static std::vector<double> phiSegments(double phiMin = -M_PI,
33  double phiMax = M_PI,
34  std::vector<double> phiRefs = {},
35  double phiTolerance = 1e-6) {
36  // This is to ensure that the extrema are built regardless of number
37  // of segments
38  std::vector<double> phiSegments;
39  std::vector<double> quarters = {-M_PI, -0.5 * M_PI, 0., 0.5 * M_PI, M_PI};
40  // It does not cover the full azimuth
41  if (phiMin != -M_PI or phiMax != M_PI) {
42  phiSegments.push_back(phiMin);
43  for (unsigned int iq = 1; iq < 4; ++iq) {
44  if (phiMin < quarters[iq] and phiMax > quarters[iq]) {
45  phiSegments.push_back(quarters[iq]);
46  }
47  }
48  phiSegments.push_back(phiMax);
49  } else {
50  phiSegments = quarters;
51  }
52  // Insert the reference phis if
53  if (not phiRefs.empty()) {
54  for (const auto& phiRef : phiRefs) {
55  // Trying to find the right patch
56  auto match = std::find_if(
57  phiSegments.begin(), phiSegments.end(), [&](double phiSeg) {
58  return std::abs(phiSeg - phiRef) < phiTolerance;
59  });
60  if (match == phiSegments.end()) {
61  phiSegments.push_back(phiRef);
62  }
63  }
64  std::sort(phiSegments.begin(), phiSegments.end());
65  }
66  return phiSegments;
67 }
68 
83 template <typename vertex_t, typename transform_t>
84 void createSegment(std::vector<vertex_t>& vertices,
85  std::pair<double, double> rxy, double phi1, double phi2,
86  unsigned int lseg, int addon = 0,
87  const vertex_t& offset = vertex_t::Zero(),
88  const transform_t& transform = transform_t::Identity()) {
89  // Calculate the number of segments - 1 is the minimum
90  unsigned int segs = std::abs(phi2 - phi1) / (2 * M_PI) * lseg;
91  segs = segs > 0 ? segs : 1;
92  double phistep = (phi2 - phi1) / segs;
93  // Create the segments
94  for (unsigned int iphi = 0; iphi < segs + addon; ++iphi) {
95  double phi = phi1 + iphi * phistep;
96  vertex_t vertex = vertex_t::Zero();
97  vertex(0) = rxy.first * std::cos(phi);
98  vertex(1) = rxy.second * std::sin(phi);
99  vertex = vertex + offset;
100  vertices.push_back(transform * vertex);
101  }
102 }
103 
115 static std::vector<Vector2D> ellispoidVertices(double innerRx, double innerRy,
116  double outerRx, double outerRy,
117  double avgPhi = 0.,
118  double halfPhi = M_PI,
119  unsigned int lseg = 1) {
120  // List of vertices counter-clockwise starting at smallest phi w.r.t center,
121  // for both inner/outer ring/segment
122  std::vector<Acts::Vector2D> rvertices; // return vertices
123  std::vector<Acts::Vector2D> ivertices; // inner vertices
124  std::vector<Acts::Vector2D> overtices; // outer verices
125 
126  bool innerExists = (innerRx > 0. and innerRy > 0.);
127  bool closed = std::abs(halfPhi - M_PI) < s_onSurfaceTolerance;
128 
129  // Get the phi segments from the helper method
131  avgPhi - halfPhi, avgPhi + halfPhi, {avgPhi});
132 
133  // The inner (if exists) and outer bow
134  for (unsigned int iseg = 0; iseg < phiSegs.size() - 1; ++iseg) {
135  int addon = (iseg == phiSegs.size() - 2 and not closed) ? 1 : 0;
136  if (innerExists) {
137  detail::VerticesHelper::createSegment<Vector2D, Eigen::Affine2d>(
138  ivertices, {innerRx, innerRy}, phiSegs[iseg], phiSegs[iseg + 1], lseg,
139  addon);
140  }
141  detail::VerticesHelper::createSegment<Vector2D, Eigen::Affine2d>(
142  overtices, {outerRx, outerRy}, phiSegs[iseg], phiSegs[iseg + 1], lseg,
143  addon);
144  }
145 
146  // We want to keep the same counter-clockwise orientation for displaying
147  if (not innerExists) {
148  if (not closed) {
149  // Add the center case we have a sector
150  rvertices.push_back(Vector2D(0., 0.));
151  }
152  rvertices.insert(rvertices.end(), overtices.begin(), overtices.end());
153  } else if (not closed) {
154  rvertices.insert(rvertices.end(), overtices.begin(), overtices.end());
155  rvertices.insert(rvertices.end(), ivertices.rbegin(), ivertices.rend());
156  } else {
157  rvertices.insert(rvertices.end(), overtices.begin(), overtices.end());
158  rvertices.insert(rvertices.end(), ivertices.begin(), ivertices.end());
159  }
160  return rvertices;
161 }
162 
172 static std::vector<Vector2D> circularVertices(double innerR, double outerR,
173  double avgPhi = 0.,
174  double halfPhi = M_PI,
175  unsigned int lseg = 1) {
176  return ellispoidVertices(innerR, innerR, outerR, outerR, avgPhi, halfPhi,
177  lseg);
178 }
179 
190 template <typename vertex_t, typename vertex_container_t>
191 bool isInsidePolygon(const vertex_t& point,
192  const vertex_container_t& vertices) {
193  // when we move along the edges of a convex polygon, a point on the inside of
194  // the polygon will always appear on the same side of each edge.
195  // a point on the outside will switch sides at least once.
196 
197  // returns which side of the connecting line between `ll0` and `ll1` the point
198  // `p` is on. computes the sign of the z-component of the cross-product
199  // between the line normal vector and the vector from `ll0` to `p`.
200  auto lineSide = [&](auto&& ll0, auto&& ll1) {
201  auto normal = ll1 - ll0;
202  auto delta = point - ll0;
203  return std::signbit((normal[0] * delta[1]) - (normal[1] * delta[0]));
204  };
205 
206  auto iv = std::begin(vertices);
207  auto l0 = *iv;
208  auto l1 = *(++iv);
209  // use vertex0 to vertex1 to define reference sign and compare w/ all edges
210  auto reference = lineSide(l0, l1);
211  for (++iv; iv != std::end(vertices); ++iv) {
212  l0 = l1;
213  l1 = *iv;
214  if (lineSide(l0, l1) != reference) {
215  return false;
216  }
217  }
218  // manual check for last edge from last vertex back to the first vertex
219  if (lineSide(l1, *std::begin(vertices)) != reference) {
220  return false;
221  }
222  // point was always on the same side. point must be inside.
223  return true;
224 }
225 
236 template <typename vertex_t>
237 bool isInsideRectangle(const vertex_t& point, const vertex_t& lowerLeft,
238  const vertex_t& upperRight) {
239  return (lowerLeft[0] <= point[0]) && (point[0] < upperRight[0]) &&
240  (lowerLeft[1] <= point[1]) && (point[1] < upperRight[1]);
241 }
242 
243 }; // namespace VerticesHelper
244 
245 } // namespace detail
246 
247 } // namespace Acts