ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DiscSurface.ipp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DiscSurface.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 // DiscSurface.ipp, Acts project
12 
13 inline const Vector2D DiscSurface::localPolarToCartesian(
14  const Vector2D& lpolar) const {
15  return Vector2D(lpolar[eLOC_R] * cos(lpolar[eLOC_PHI]),
16  lpolar[eLOC_R] * sin(lpolar[eLOC_PHI]));
17 }
18 
19 inline const Vector2D DiscSurface::localCartesianToPolar(
20  const Vector2D& lcart) const {
21  return Vector2D(
22  sqrt(lcart[eLOC_X] * lcart[eLOC_X] + lcart[eLOC_Y] * lcart[eLOC_Y]),
23  atan2(lcart[eLOC_Y], lcart[eLOC_X]));
24 }
25 
26 inline void DiscSurface::initJacobianToGlobal(const GeometryContext& gctx,
27  BoundToFreeMatrix& jacobian,
28  const Vector3D& position,
29  const Vector3D& direction,
30  const BoundVector& pars) const {
31  // The trigonometry required to convert the direction to spherical
32  // coordinates and then compute the sines and cosines again can be
33  // surprisingly expensive from a performance point of view.
34  //
35  // Here, we can avoid it because the direction is by definition a unit
36  // vector, with the following coordinate conversions...
37  const double x = direction(0); // == cos(phi) * sin(theta)
38  const double y = direction(1); // == sin(phi) * sin(theta)
39  const double z = direction(2); // == cos(theta)
40 
41  // ...which we can invert to directly get the sines and cosines:
42  const double cos_theta = z;
43  const double sin_theta = sqrt(x * x + y * y);
44  const double inv_sin_theta = 1. / sin_theta;
45  const double cos_phi = x * inv_sin_theta;
46  const double sin_phi = y * inv_sin_theta;
47  // retrieve the reference frame
48  const auto rframe = referenceFrame(gctx, position, direction);
49 
50  // special polar coordinates for the Disc
51  double lrad = pars[eLOC_0];
52  double lphi = pars[eLOC_1];
53  double lcos_phi = cos(lphi);
54  double lsin_phi = sin(lphi);
55  // the local error components - rotated from reference frame
56  jacobian.block<3, 1>(0, eLOC_0) =
57  lcos_phi * rframe.block<3, 1>(0, 0) + lsin_phi * rframe.block<3, 1>(0, 1);
58  jacobian.block<3, 1>(0, eLOC_1) =
59  lrad * (lcos_phi * rframe.block<3, 1>(0, 1) -
60  lsin_phi * rframe.block<3, 1>(0, 0));
61  // the time component
62  jacobian(3, eT) = 1;
63  // the momentum components
64  jacobian(4, ePHI) = (-sin_theta) * sin_phi;
65  jacobian(4, eTHETA) = cos_theta * cos_phi;
66  jacobian(5, ePHI) = sin_theta * cos_phi;
67  jacobian(5, eTHETA) = cos_theta * sin_phi;
68  jacobian(6, eTHETA) = (-sin_theta);
69  jacobian(7, eQOP) = 1;
70 }
71 
72 inline const RotationMatrix3D DiscSurface::initJacobianToLocal(
73  const GeometryContext& gctx, FreeToBoundMatrix& jacobian,
74  const Vector3D& position, const Vector3D& direction) const {
75  using VectorHelpers::perp;
76  using VectorHelpers::phi;
77  // Optimized trigonometry on the propagation direction
78  const double x = direction(0); // == cos(phi) * sin(theta)
79  const double y = direction(1); // == sin(phi) * sin(theta)
80  const double z = direction(2); // == cos(theta)
81  // can be turned into cosine/sine
82  const double cosTheta = z;
83  const double sinTheta = sqrt(x * x + y * y);
84  const double invSinTheta = 1. / sinTheta;
85  const double cosPhi = x * invSinTheta;
86  const double sinPhi = y * invSinTheta;
87  // The measurement frame of the surface
88  RotationMatrix3D rframeT =
89  referenceFrame(gctx, position, direction).transpose();
90  // calculate the transformation to local coorinates
91  const Vector3D pos_loc = transform(gctx).inverse() * position;
92  const double lr = perp(pos_loc);
93  const double lphi = phi(pos_loc);
94  const double lcphi = cos(lphi);
95  const double lsphi = sin(lphi);
96  // rotate into the polar coorindates
97  auto lx = rframeT.block<1, 3>(0, 0);
98  auto ly = rframeT.block<1, 3>(1, 0);
99  jacobian.block<1, 3>(0, 0) = lcphi * lx + lsphi * ly;
100  jacobian.block<1, 3>(1, 0) = (lcphi * ly - lsphi * lx) / lr;
101  // Time element
102  jacobian(eT, 3) = 1;
103  // Directional and momentum elements for reference frame surface
104  jacobian(ePHI, 4) = -sinPhi * invSinTheta;
105  jacobian(ePHI, 5) = cosPhi * invSinTheta;
106  jacobian(eTHETA, 4) = cosPhi * cosTheta;
107  jacobian(eTHETA, 5) = sinPhi * cosTheta;
108  jacobian(eTHETA, 6) = -sinTheta;
109  jacobian(eQOP, 7) = 1;
110  // return the transposed reference frame
111  return rframeT;
112 }
113 
115  const GeometryContext& gctx, const Vector3D& position,
116  const Vector3D& direction, const BoundaryCheck& bcheck) const {
117  // Get the contextual transform
118  auto gctxTransform = transform(gctx);
119  // Use the intersection helper for planar surfaces
120  auto intersection =
121  PlanarHelper::intersectionEstimate(gctxTransform, position, direction);
122  // Evaluate boundary check if requested (and reachable)
123  if (intersection.status != Intersection::Status::unreachable and bcheck and
124  m_bounds != nullptr) {
125  // Built-in local to global for speed reasons
126  const auto& tMatrix = gctxTransform.matrix();
127  const Vector3D vecLocal(intersection.position - tMatrix.block<3, 1>(0, 3));
128  const Vector2D lcartesian =
129  tMatrix.block<3, 2>(0, 0).transpose() * vecLocal;
130  if (bcheck.type() == BoundaryCheck::Type::eAbsolute and
131  m_bounds->coversFullAzimuth()) {
132  double tolerance = s_onSurfaceTolerance + bcheck.tolerance()[eLOC_R];
133  if (not m_bounds->insideRadialBounds(VectorHelpers::perp(lcartesian),
134  tolerance)) {
135  intersection.status = Intersection::Status::missed;
136  }
137  } else if (not insideBounds(localCartesianToPolar(lcartesian), bcheck)) {
138  intersection.status = Intersection::Status::missed;
139  }
140  }
141  return intersection;
142 }
143 
144 inline const Vector3D DiscSurface::normal(const GeometryContext& gctx,
145  const Vector2D& /*unused*/) const {
146  // fast access via tranform matrix (and not rotation())
147  const auto& tMatrix = transform(gctx).matrix();
148  return Vector3D(tMatrix(0, 2), tMatrix(1, 2), tMatrix(2, 2));
149 }
150 
151 inline const Vector3D DiscSurface::binningPosition(const GeometryContext& gctx,
152  BinningValue bValue) const {
153  if (bValue == binR) {
154  double r = m_bounds->binningValueR();
155  double phi = m_bounds->binningValuePhi();
156  return Vector3D(r * cos(phi), r * sin(phi), center(gctx).z());
157  }
158  return center(gctx);
159 }
160 
161 inline double DiscSurface::binningPositionValue(const GeometryContext& gctx,
162  BinningValue bValue) const {
163  // only modify binR
164  if (bValue == binR) {
165  return VectorHelpers::perp(center(gctx)) + m_bounds->binningValueR();
166  }
167  return GeometryObject::binningPositionValue(gctx, bValue);
168 }
169 
170 inline double DiscSurface::pathCorrection(const GeometryContext& gctx,
171  const Vector3D& position,
172  const Vector3D& direction) const {
174  return 1. / std::abs(Surface::normal(gctx, position).dot(direction));
175 }