Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
LineSurface.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file LineSurface.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2016-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 
10 
22 
23 #include <algorithm>
24 #include <cmath>
25 #include <limits>
26 #include <stdexcept>
27 #include <utility>
28 
29 namespace Acts {
30 class DetectorElementBase;
31 } // namespace Acts
32 
34  double halez)
35  : GeometryObject(),
36  Surface(transform),
37  m_bounds(std::make_shared<const LineBounds>(radius, halez)) {}
38 
40  std::shared_ptr<const LineBounds> lbounds)
41  : GeometryObject(), Surface(transform), m_bounds(std::move(lbounds)) {}
42 
43 Acts::LineSurface::LineSurface(std::shared_ptr<const LineBounds> lbounds,
44  const DetectorElementBase& detelement)
45  : GeometryObject(), Surface(detelement), m_bounds(std::move(lbounds)) {
46  throw_assert(m_bounds, "LineBounds must not be nullptr");
47 }
48 
50  : GeometryObject(), Surface(other), m_bounds(other.m_bounds) {}
51 
53  const LineSurface& other,
54  const Transform3& shift)
55  : GeometryObject(), Surface(gctx, other, shift), m_bounds(other.m_bounds) {}
56 
58  if (this != &other) {
59  Surface::operator=(other);
60  m_bounds = other.m_bounds;
61  }
62  return *this;
63 }
64 
66  const Vector2& lposition,
67  const Vector3& direction) const {
68  Vector3 unitZ0 = lineDirection(gctx);
69 
70  // get the vector perpendicular to the momentum direction and the straw axis
71  Vector3 radiusAxisGlobal = unitZ0.cross(direction);
72  Vector3 locZinGlobal =
73  transform(gctx) * Vector3(0., 0., lposition[eBoundLoc1]);
74  // add eBoundLoc0 * radiusAxis
75  return Vector3(locZinGlobal +
76  lposition[eBoundLoc0] * radiusAxisGlobal.normalized());
77 }
78 
80  const GeometryContext& gctx, const Vector3& position,
81  const Vector3& direction, double tolerance) const {
82  using VectorHelpers::perp;
83 
84  // Bring the global position into the local frame. First remove the
85  // translation then the rotation.
86  Vector3 localPosition = referenceFrame(gctx, position, direction).inverse() *
87  (position - transform(gctx).translation());
88 
89  // `localPosition.z()` is not the distance to the PCA but the smallest
90  // distance between `position` and the imaginary plane surface defined by the
91  // local x,y axes in the global frame and the position of the line surface.
92  //
93  // This check is also done for the `PlaneSurface` so I aligned the
94  // `LineSurface` to do the same thing.
95  if (std::abs(localPosition.z()) > std::abs(tolerance)) {
96  return Result<Vector2>::failure(SurfaceError::GlobalPositionNotOnSurface);
97  }
98 
99  // Construct result from local x,y
100  Vector2 localXY = localPosition.head<2>();
101 
102  return Result<Vector2>::success(localXY);
103 }
104 
106  return "Acts::LineSurface";
107 }
108 
110  const GeometryContext& gctx, const Vector3& /*position*/,
111  const Vector3& direction) const {
112  Vector3 unitZ0 = lineDirection(gctx);
113  Vector3 unitD0 = unitZ0.cross(direction).normalized();
114  Vector3 unitDistance = unitD0.cross(unitZ0);
115 
116  RotationMatrix3 mFrame;
117  mFrame.col(0) = unitD0;
118  mFrame.col(1) = unitZ0;
119  mFrame.col(2) = unitDistance;
120 
121  return mFrame;
122 }
123 
125  const Vector3& /*pos*/,
126  const Vector3& /*mom*/) const {
127  return 1.;
128 }
129 
131  const GeometryContext& gctx, BinningValue /*bValue*/) const {
132  return center(gctx);
133 }
134 
136  const Vector2& /*lpos*/) const {
137  throw std::runtime_error(
138  "LineSurface: normal is undefined without known direction");
139 }
140 
142  if (m_bounds) {
143  return (*m_bounds.get());
144  }
145  return s_noBounds;
146 }
147 
149  const GeometryContext& gctx, const Vector3& position,
150  const Vector3& direction, const BoundaryCheck& bcheck,
151  ActsScalar tolerance) const {
152  // The nomenclature is following the header file and doxygen documentation
153 
154  const Vector3& ma = position;
155  const Vector3& ea = direction;
156 
157  // Origin of the line surface
158  Vector3 mb = transform(gctx).translation();
159  // Line surface axis
160  Vector3 eb = lineDirection(gctx);
161 
162  // Now go ahead and solve for the closest approach
163  Vector3 mab = mb - ma;
164  double eaTeb = ea.dot(eb);
165  double denom = 1 - eaTeb * eaTeb;
166 
167  // `tolerance` does not really have a meaning here it is just a sufficiently
168  // small number so `u` does not explode
169  if (std::abs(denom) < std::abs(tolerance)) {
170  // return a false intersection
171  return {{Intersection3D::invalid(), Intersection3D::invalid()}, this};
172  }
173 
174  double u = (mab.dot(ea) - mab.dot(eb) * eaTeb) / denom;
175  // Check if we are on the surface already
176  Intersection3D::Status status = std::abs(u) > std::abs(tolerance)
177  ? Intersection3D::Status::reachable
178  : Intersection3D::Status::onSurface;
179  Vector3 result = ma + u * ea;
180  // Evaluate the boundary check if requested
181  // m_bounds == nullptr prevents unnecessary calculations for PerigeeSurface
182  if (bcheck and m_bounds) {
183  // At closest approach: check inside R or and inside Z
184  Vector3 vecLocal = result - mb;
185  double cZ = vecLocal.dot(eb);
186  double hZ = m_bounds->get(LineBounds::eHalfLengthZ) + tolerance;
187  if ((std::abs(cZ) > std::abs(hZ)) or
188  ((vecLocal - cZ * eb).norm() >
189  m_bounds->get(LineBounds::eR) + tolerance)) {
190  status = Intersection3D::Status::missed;
191  }
192  }
193 
194  return {{Intersection3D(result, u, status), Intersection3D::invalid()}, this};
195 }
196 
198  const GeometryContext& gctx, const BoundVector& boundParams) const {
199  // Transform from bound to free parameters
200  FreeVector freeParams =
201  detail::transformBoundToFreeParameters(*this, gctx, boundParams);
202  // The global position
203  Vector3 position = freeParams.segment<3>(eFreePos0);
204  // The direction
205  Vector3 direction = freeParams.segment<3>(eFreeDir0);
206  // Get the sines and cosines directly
207  double cosTheta = std::cos(boundParams[eBoundTheta]);
208  double sinTheta = std::sin(boundParams[eBoundTheta]);
209  double cosPhi = std::cos(boundParams[eBoundPhi]);
210  double sinPhi = std::sin(boundParams[eBoundPhi]);
211  // retrieve the reference frame
212  auto rframe = referenceFrame(gctx, position, direction);
213 
214  // Initialize the jacobian from local to global
215  BoundToFreeMatrix jacToGlobal = BoundToFreeMatrix::Zero();
216 
217  // the local error components - given by the reference frame
218  jacToGlobal.topLeftCorner<3, 2>() = rframe.topLeftCorner<3, 2>();
219  // the time component
220  jacToGlobal(eFreeTime, eBoundTime) = 1;
221  // the momentum components
222  jacToGlobal(eFreeDir0, eBoundPhi) = -sinTheta * sinPhi;
223  jacToGlobal(eFreeDir0, eBoundTheta) = cosTheta * cosPhi;
224  jacToGlobal(eFreeDir1, eBoundPhi) = sinTheta * cosPhi;
225  jacToGlobal(eFreeDir1, eBoundTheta) = cosTheta * sinPhi;
226  jacToGlobal(eFreeDir2, eBoundTheta) = -sinTheta;
227  jacToGlobal(eFreeQOverP, eBoundQOverP) = 1;
228 
229  // the projection of direction onto ref frame normal
230  double ipdn = 1. / direction.dot(rframe.col(2));
231  // build the cross product of d(D)/d(eBoundPhi) components with y axis
232  Vector3 dDPhiY = rframe.block<3, 1>(0, 1).cross(
233  jacToGlobal.block<3, 1>(eFreeDir0, eBoundPhi));
234  // and the same for the d(D)/d(eTheta) components
235  Vector3 dDThetaY = rframe.block<3, 1>(0, 1).cross(
236  jacToGlobal.block<3, 1>(eFreeDir0, eBoundTheta));
237  // and correct for the x axis components
238  dDPhiY -= rframe.block<3, 1>(0, 0) * (rframe.block<3, 1>(0, 0).dot(dDPhiY));
239  dDThetaY -=
240  rframe.block<3, 1>(0, 0) * (rframe.block<3, 1>(0, 0).dot(dDThetaY));
241  // set the jacobian components for global d/ phi/Theta
242  jacToGlobal.block<3, 1>(eFreePos0, eBoundPhi) =
243  dDPhiY * boundParams[eBoundLoc0] * ipdn;
244  jacToGlobal.block<3, 1>(eFreePos0, eBoundTheta) =
245  dDThetaY * boundParams[eBoundLoc0] * ipdn;
246 
247  return jacToGlobal;
248 }
249 
251  const GeometryContext& gctx, const FreeVector& parameters) const {
252  // The global posiiton
253  Vector3 position = parameters.segment<3>(eFreePos0);
254  // The direction
255  Vector3 direction = parameters.segment<3>(eFreeDir0);
256  // The vector between position and center
257  Vector3 pcRowVec = position - center(gctx);
258  // The local frame z axis
259  Vector3 localZAxis = lineDirection(gctx);
260  // The local z coordinate
261  double pz = pcRowVec.dot(localZAxis);
262  // Cosine of angle between momentum direction and local frame z axis
263  double dz = localZAxis.dot(direction);
264  double norm = 1 / (1 - dz * dz);
265 
266  // Initialize the derivative of propagation path w.r.t. free parameter
267  FreeToPathMatrix freeToPath = FreeToPathMatrix::Zero();
268 
269  // The derivative of path w.r.t. position
270  freeToPath.segment<3>(eFreePos0) =
271  norm * (dz * localZAxis.transpose() - direction.transpose());
272 
273  // The derivative of path w.r.t. direction
274  freeToPath.segment<3>(eFreeDir0) =
275  norm * (pz * localZAxis.transpose() - pcRowVec.transpose());
276 
277  return freeToPath;
278 }
279 
281  const GeometryContext& gctx, const FreeVector& parameters) const {
282  // The global posiiton
283  Vector3 position = parameters.segment<3>(eFreePos0);
284  // The direction
285  Vector3 direction = parameters.segment<3>(eFreeDir0);
286  // The vector between position and center
287  Vector3 pcRowVec = position - center(gctx);
288  // The local frame z axis
289  Vector3 localZAxis = lineDirection(gctx);
290  // The local z coordinate
291  double pz = pcRowVec.dot(localZAxis);
292  // Cosine of angle between momentum direction and local frame z axis
293  double dz = localZAxis.dot(direction);
294  double norm = 1 / (1 - dz * dz);
295  // Calculate the derivative of local frame axes w.r.t its rotation
296  auto [rotToLocalXAxis, rotToLocalYAxis, rotToLocalZAxis] =
298 
299  // Initialize the derivative of propagation path w.r.t. local frame
300  // translation (origin) and rotation
301  AlignmentToPathMatrix alignToPath = AlignmentToPathMatrix::Zero();
302  alignToPath.segment<3>(eAlignmentCenter0) =
303  norm * (direction.transpose() - dz * localZAxis.transpose());
304  alignToPath.segment<3>(eAlignmentRotation0) =
305  norm * (dz * pcRowVec.transpose() + pz * direction.transpose()) *
306  rotToLocalZAxis;
307 
308  return alignToPath;
309 }
310 
312  const GeometryContext& gctx, const Vector3& position) const {
313  // calculate the transformation to local coordinates
314  Vector3 localPosition = transform(gctx).inverse() * position;
315  double localPhi = VectorHelpers::phi(localPosition);
316 
317  ActsMatrix<2, 3> loc3DToLocBound = ActsMatrix<2, 3>::Zero();
318  loc3DToLocBound << std::cos(localPhi), std::sin(localPhi), 0, 0, 0, 1;
319 
320  return loc3DToLocBound;
321 }
322 
324  const GeometryContext& gctx) const {
325  return transform(gctx).linear().col(2);
326 }