EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AnnulusBounds.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file AnnulusBounds.cpp
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 #include "Acts/Surfaces/AnnulusBounds.hpp"
10 
11 #include "Acts/Surfaces/detail/VerticesHelper.hpp"
12 #include "Acts/Utilities/Helpers.hpp"
13 #include "Acts/Utilities/detail/periodic.hpp"
14 
15 #include <cmath>
16 #include <iomanip>
17 #include <iostream>
18 
19 Acts::AnnulusBounds::AnnulusBounds(
20  const std::array<double, eSize>& values) noexcept(false)
21  : m_values(values), m_moduleOrigin({values[eOriginX], values[eOriginY]}) {
22  checkConsistency();
23 
24  m_rotationStripPC =
25  Eigen::Translation<double, 2>(Vector2D(0, -get(eAveragePhi)));
26  m_translation = Eigen::Translation<double, 2>(m_moduleOrigin);
27 
28  m_shiftXY = m_moduleOrigin * -1;
29  m_shiftPC =
30  Vector2D(VectorHelpers::perp(m_shiftXY), VectorHelpers::phi(m_shiftXY));
31 
32  // we need the corner points of the module to do the inside
33  // checking, calculate them here once, they don't change
34 
35  // find inner outer radius at edges in STRIP PC
36  auto circIx = [](double O_x, double O_y, double r, double phi) -> Vector2D {
37  // _____________________________________________
38  // / 2 2 2 2 2 2
39  // O_x + O_y*m - \/ - O_x *m + 2*O_x*O_y*m - O_y + m *r + r
40  // x = --------------------------------------------------------------
41  // 2
42  // m + 1
43  //
44  // y = m*x
45  //
46  double m = std::tan(phi);
47  Vector2D dir(std::cos(phi), std::sin(phi));
48  double x1 = (O_x + O_y * m -
49  std::sqrt(-std::pow(O_x, 2) * std::pow(m, 2) +
50  2 * O_x * O_y * m - std::pow(O_y, 2) +
51  std::pow(m, 2) * std::pow(r, 2) + std::pow(r, 2))) /
52  (std::pow(m, 2) + 1);
53  double x2 = (O_x + O_y * m +
54  std::sqrt(-std::pow(O_x, 2) * std::pow(m, 2) +
55  2 * O_x * O_y * m - std::pow(O_y, 2) +
56  std::pow(m, 2) * std::pow(r, 2) + std::pow(r, 2))) /
57  (std::pow(m, 2) + 1);
58 
59  Vector2D v1(x1, m * x1);
60  if (v1.dot(dir) > 0)
61  return v1;
62  return {x2, m * x2};
63  };
64 
65  // calculate corners in STRIP XY, keep them we need them for minDistance()
66  m_outLeftStripXY =
67  circIx(m_moduleOrigin[eBoundLoc0], m_moduleOrigin[eBoundLoc1], get(eMaxR),
68  get(eMaxPhiRel));
69  m_inLeftStripXY =
70  circIx(m_moduleOrigin[eBoundLoc0], m_moduleOrigin[eBoundLoc1], get(eMinR),
71  get(eMaxPhiRel));
72  m_outRightStripXY =
73  circIx(m_moduleOrigin[eBoundLoc0], m_moduleOrigin[eBoundLoc1], get(eMaxR),
74  get(eMinPhiRel));
75  m_inRightStripXY =
76  circIx(m_moduleOrigin[eBoundLoc0], m_moduleOrigin[eBoundLoc1], get(eMinR),
77  get(eMinPhiRel));
78 
79  m_outLeftStripPC = {m_outLeftStripXY.norm(),
80  VectorHelpers::phi(m_outLeftStripXY)};
81  m_inLeftStripPC = {m_inLeftStripXY.norm(),
82  VectorHelpers::phi(m_inLeftStripXY)};
83  m_outRightStripPC = {m_outRightStripXY.norm(),
84  VectorHelpers::phi(m_outRightStripXY)};
85  m_inRightStripPC = {m_inRightStripXY.norm(),
86  VectorHelpers::phi(m_inRightStripXY)};
87 
88  m_outLeftModulePC = stripXYToModulePC(m_outLeftStripXY);
89  m_inLeftModulePC = stripXYToModulePC(m_inLeftStripXY);
90  m_outRightModulePC = stripXYToModulePC(m_outRightStripXY);
91  m_inRightModulePC = stripXYToModulePC(m_inRightStripXY);
92 }
93 
94 std::vector<Acts::Vector2D> Acts::AnnulusBounds::corners() const {
95  auto rot = m_rotationStripPC.inverse();
96 
97  return {rot * m_outRightStripPC, rot * m_outLeftStripPC,
98  rot * m_inLeftStripPC, rot * m_inRightStripPC};
99 }
100 
101 std::vector<Acts::Vector2D> Acts::AnnulusBounds::vertices(
102  unsigned int lseg) const {
103  // List of vertices counter-clockwise starting with left inner
104  std::vector<Acts::Vector2D> rvertices;
105 
106  double phiMinInner = VectorHelpers::phi(m_inLeftStripXY);
107  double phiMaxInner = VectorHelpers::phi(m_inRightStripXY);
108  double phiMinOuter = VectorHelpers::phi(m_outRightStripXY);
109  double phiMaxOuter = VectorHelpers::phi(m_outLeftStripXY);
110 
111  std::vector<double> phisInner =
112  detail::VerticesHelper::phiSegments(phiMinInner, phiMaxInner);
113  std::vector<double> phisOuter =
114  detail::VerticesHelper::phiSegments(phiMinOuter, phiMaxOuter);
115 
116  // Inner bow from phi_min -> phi_max
117  for (unsigned int iseg = 0; iseg < phisInner.size() - 1; ++iseg) {
118  int addon = (iseg == phisInner.size() - 2) ? 1 : 0;
119  detail::VerticesHelper::createSegment<Vector2D, Transform2D>(
120  rvertices, {get(eMinR), get(eMinR)}, phisInner[iseg],
121  phisInner[iseg + 1], lseg, addon);
122  }
123  // Upper bow from phi_min -> phi_max
124  for (unsigned int iseg = 0; iseg < phisOuter.size() - 1; ++iseg) {
125  int addon = (iseg == phisOuter.size() - 2) ? 1 : 0;
126  detail::VerticesHelper::createSegment<Vector2D, Transform2D>(
127  rvertices, {get(eMaxR), get(eMaxR)}, phisOuter[iseg],
128  phisOuter[iseg + 1], lseg, addon);
129  }
130 
131  return rvertices;
132 }
133 
134 bool Acts::AnnulusBounds::inside(const Vector2D& lposition, double tolR,
135  double tolPhi) const {
136  // locpo is PC in STRIP SYSTEM
137  // need to perform internal rotation induced by m_phiAvg
138  Vector2D locpo_rotated = m_rotationStripPC * lposition;
139  double phiLoc = locpo_rotated[eBoundLoc1];
140  double rLoc = locpo_rotated[eBoundLoc0];
141 
142  if (phiLoc < (get(eMinPhiRel) - tolPhi) ||
143  phiLoc > (get(eMaxPhiRel) + tolPhi)) {
144  return false;
145  }
146 
147  // calculate R in MODULE SYSTEM to evaluate R-bounds
148  if (tolR == 0.) {
149  // don't need R, can use R^2
150  double r_mod2 =
151  m_shiftPC[eBoundLoc0] * m_shiftPC[eBoundLoc0] + rLoc * rLoc +
152  2 * m_shiftPC[eBoundLoc0] * rLoc * cos(phiLoc - m_shiftPC[eBoundLoc1]);
153 
154  if (r_mod2 < get(eMinR) * get(eMinR) || r_mod2 > get(eMaxR) * get(eMaxR)) {
155  return false;
156  }
157  } else {
158  // use R
159  double r_mod = sqrt(
160  m_shiftPC[eBoundLoc0] * m_shiftPC[eBoundLoc0] + rLoc * rLoc +
161  2 * m_shiftPC[eBoundLoc0] * rLoc * cos(phiLoc - m_shiftPC[eBoundLoc1]));
162 
163  if (r_mod < (get(eMinR) - tolR) || r_mod > (get(eMaxR) + tolR)) {
164  return false;
165  }
166  }
167  return true;
168 }
169 
170 bool Acts::AnnulusBounds::inside(const Vector2D& lposition,
171  const BoundaryCheck& bcheck) const {
172  // locpo is PC in STRIP SYSTEM
173  if (bcheck.type() == BoundaryCheck::Type::eAbsolute) {
174  return inside(lposition, bcheck.tolerance()[eBoundLoc0],
175  bcheck.tolerance()[eBoundLoc1]);
176  } else {
177  // first check if inside. We don't need to look into the covariance if
178  // inside
179  if (inside(lposition, 0., 0.)) {
180  return true;
181  }
182 
183  // we need to rotated the locpo
184  Vector2D locpo_rotated = m_rotationStripPC * lposition;
185 
186  // covariance is given in STRIP SYSTEM in PC
187  // we need to convert the covariance to the MODULE SYSTEM in PC
188  // via jacobian.
189  // The following transforms into STRIP XY, does the shift into MODULE XY,
190  // and then transforms into MODULE PC
191  double dphi = m_phiAvg;
192  double phi_strip = locpo_rotated[eBoundLoc1];
193  double r_strip = locpo_rotated[eBoundLoc0];
194  double O_x = m_shiftXY[eBoundLoc0];
195  double O_y = m_shiftXY[eBoundLoc1];
196 
197  // For a transformation from cartesian into polar coordinates
198  //
199  // [ _________ ]
200  // [ / 2 2 ]
201  // [ \/ x + y ]
202  // [ r' ] [ ]
203  // v = [ ] = [ / y \]
204  // [phi'] [2*atan|----------------|]
205  // [ | _________|]
206  // [ | / 2 2 |]
207  // [ \x + \/ x + y /]
208  //
209  // Where x, y are polar coordinates that can be rotated by dPhi
210  //
211  // [x] [O_x + r*cos(dPhi - phi)]
212  // [ ] = [ ]
213  // [y] [O_y - r*sin(dPhi - phi)]
214  //
215  // The general jacobian is:
216  //
217  // [d d ]
218  // [--(f_x) --(f_x)]
219  // [dx dy ]
220  // Jgen = [ ]
221  // [d d ]
222  // [--(f_y) --(f_y)]
223  // [dx dy ]
224  //
225  // which means in this case:
226  //
227  // [ d d ]
228  // [ ----------(rMod) ---------(rMod) ]
229  // [ dr_{strip} dphiStrip ]
230  // J = [ ]
231  // [ d d ]
232  // [----------(phiMod) ---------(phiMod)]
233  // [dr_{strip} dphiStrip ]
234  //
235  // Performing the derivative one gets:
236  //
237  // [B*O_x + C*O_y + rStrip rStrip*(B*O_y + O_x*sin(dPhi - phiStrip))]
238  // [---------------------- -----------------------------------------]
239  // [ ___ ___ ]
240  // [ \/ A \/ A ]
241  // J = [ ]
242  // [ -(B*O_y - C*O_x) rStrip*(B*O_x + C*O_y + rStrip) ]
243  // [ ----------------- ------------------------------- ]
244  // [ A A ]
245  //
246  // where
247  // 2 2 2
248  // A = O_x + 2*O_x*rStrip*cos(dPhi - phiStrip) + O_y -
249  // 2*O_y*rStrip*sin(dPhi - phiStrip) + rStrip B = cos(dPhi - phiStrip) C =
250  // -sin(dPhi - phiStrip)
251 
252  double cosDPhiPhiStrip = std::cos(dphi - phi_strip);
253  double sinDPhiPhiStrip = std::sin(dphi - phi_strip);
254 
255  double A = O_x * O_x + 2 * O_x * r_strip * cosDPhiPhiStrip + O_y * O_y -
256  2 * O_y * r_strip * sinDPhiPhiStrip + r_strip * r_strip;
257  double sqrtA = std::sqrt(A);
258 
259  double B = cosDPhiPhiStrip;
260  double C = -sinDPhiPhiStrip;
261  Eigen::Matrix<double, 2, 2> jacobianStripPCToModulePC;
262  jacobianStripPCToModulePC(0, 0) = (B * O_x + C * O_y + r_strip) / sqrtA;
263  jacobianStripPCToModulePC(0, 1) =
264  r_strip * (B * O_y + O_x * sinDPhiPhiStrip) / sqrtA;
265  jacobianStripPCToModulePC(1, 0) = -(B * O_y - C * O_x) / A;
266  jacobianStripPCToModulePC(1, 1) =
267  r_strip * (B * O_x + C * O_y + r_strip) / A;
268 
269  // covariance is given in STRIP PC
270  auto covStripPC = bcheck.covariance();
271  // calculate covariance in MODULE PC using jacobian from above
272  auto covModulePC = jacobianStripPCToModulePC * covStripPC *
273  jacobianStripPCToModulePC.transpose();
274 
275  // Mahalanobis distance uses inverse covariance as weights
276  auto weightStripPC = covStripPC.inverse();
277  auto weightModulePC = covModulePC.inverse();
278 
279  double minDist = std::numeric_limits<double>::max();
280 
281  Vector2D currentClosest;
282  double currentDist;
283 
284  // do projection in STRIP PC
285 
286  // first: STRIP system. locpo is in STRIP PC already
287  currentClosest = closestOnSegment(m_inLeftStripPC, m_outLeftStripPC,
288  locpo_rotated, weightStripPC);
289  currentDist = squaredNorm(locpo_rotated - currentClosest, weightStripPC);
290  minDist = currentDist;
291 
292  currentClosest = closestOnSegment(m_inRightStripPC, m_outRightStripPC,
293  locpo_rotated, weightStripPC);
294  currentDist = squaredNorm(locpo_rotated - currentClosest, weightStripPC);
295  if (currentDist < minDist) {
296  minDist = currentDist;
297  }
298 
299  // now: MODULE system. Need to transform locpo to MODULE PC
300  // transform is STRIP PC -> STRIP XY -> MODULE XY -> MODULE PC
301  Vector2D locpoStripXY(
302  locpo_rotated[eBoundLoc0] * std::cos(locpo_rotated[eBoundLoc1]),
303  locpo_rotated[eBoundLoc0] * std::sin(locpo_rotated[eBoundLoc1]));
304  Vector2D locpoModulePC = stripXYToModulePC(locpoStripXY);
305 
306  // now check edges in MODULE PC (inner and outer circle)
307  // assuming Mahalanobis distances are of same unit if covariance
308  // is correctly transformed
309  currentClosest = closestOnSegment(m_inLeftModulePC, m_inRightModulePC,
310  locpoModulePC, weightModulePC);
311  currentDist = squaredNorm(locpoModulePC - currentClosest, weightModulePC);
312  if (currentDist < minDist) {
313  minDist = currentDist;
314  }
315 
316  currentClosest = closestOnSegment(m_outLeftModulePC, m_outRightModulePC,
317  locpoModulePC, weightModulePC);
318  currentDist = squaredNorm(locpoModulePC - currentClosest, weightModulePC);
319  if (currentDist < minDist) {
320  minDist = currentDist;
321  }
322 
323  // compare resulting Mahalanobis distance to configured
324  // "number of sigmas"
325  // we square it b/c we never took the square root of the distance
326  return minDist < bcheck.tolerance()[0] * bcheck.tolerance()[0];
327  }
328 }
329 
330 Acts::Vector2D Acts::AnnulusBounds::stripXYToModulePC(
331  const Vector2D& vStripXY) const {
332  Vector2D vecModuleXY = vStripXY + m_shiftXY;
333  return {vecModuleXY.norm(), VectorHelpers::phi(vecModuleXY)};
334 }
335 
336 Acts::Vector2D Acts::AnnulusBounds::closestOnSegment(
337  const Vector2D& a, const Vector2D& b, const Vector2D& p,
338  const Eigen::Matrix<double, 2, 2>& weight) const {
339  // connecting vector
340  auto n = b - a;
341  // squared norm of line
342  auto f = (n.transpose() * weight * n).value();
343  // weighted scalar product of line to point and segment line
344  auto u = ((p - a).transpose() * weight * n).value() / f;
345  // clamp to [0, 1], convert to point
346  return std::min(std::max(u, 0.0), 1.0) * n + a;
347 }
348 
349 double Acts::AnnulusBounds::squaredNorm(
350  const Vector2D& v, const Eigen::Matrix<double, 2, 2>& weight) const {
351  return (v.transpose() * weight * v).value();
352 }
353 
354 Acts::Vector2D Acts::AnnulusBounds::moduleOrigin() const {
355  return Eigen::Rotation2D<double>(m_phiAvg) * m_moduleOrigin;
356 }
357 
358 // Ostream operator overload
359 std::ostream& Acts::AnnulusBounds::toStream(std::ostream& sl) const {
360  sl << std::setiosflags(std::ios::fixed);
361  sl << std::setprecision(7);
362  sl << "Acts::AnnulusBounds: (innerRadius, outerRadius, minPhi, maxPhi) = ";
363  sl << "(" << get(eMinR) << ", " << get(eMaxR) << ", " << phiMin() << ", "
364  << phiMax() << ")" << '\n';
365  sl << " - shift xy = " << m_shiftXY.x() << ", " << m_shiftXY.y() << '\n';
366  ;
367  sl << " - shift pc = " << m_shiftPC.x() << ", " << m_shiftPC.y() << '\n';
368  ;
369  sl << std::setprecision(-1);
370  return sl;
371 }