EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
LineSurface.ipp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file LineSurface.ipp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2018-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 inline Vector3D LineSurface::localToGlobal(const GeometryContext& gctx,
10  const Vector2D& lposition,
11  const Vector3D& momentum) const {
12  const auto& sTransform = transform(gctx);
13  const auto& tMatrix = sTransform.matrix();
14  Vector3D lineDirection(tMatrix(0, 2), tMatrix(1, 2), tMatrix(2, 2));
15 
16  // get the vector perpendicular to the momentum and the straw axis
17  Vector3D radiusAxisGlobal(lineDirection.cross(momentum));
18  Vector3D locZinGlobal = sTransform * Vector3D(0., 0., lposition[eBoundLoc1]);
19  // add eBoundLoc0 * radiusAxis
20  return Vector3D(locZinGlobal +
21  lposition[eBoundLoc0] * radiusAxisGlobal.normalized());
22 }
23 
24 inline Result<Vector2D> LineSurface::globalToLocal(
25  const GeometryContext& gctx, const Vector3D& position,
26  const Vector3D& momentum) const {
27  using VectorHelpers::perp;
28 
29  const auto& sTransform = transform(gctx);
30  const auto& tMatrix = sTransform.matrix();
31  Vector3D lineDirection(tMatrix(0, 2), tMatrix(1, 2), tMatrix(2, 2));
32  // Bring the global position into the local frame
33  Vector3D loc3Dframe = sTransform.inverse() * position;
34  // construct localPosition with sign*perp(candidate) and z.()
35  Vector2D lposition(perp(loc3Dframe), loc3Dframe.z());
36  Vector3D sCenter(tMatrix(0, 3), tMatrix(1, 3), tMatrix(2, 3));
37  Vector3D decVec(position - sCenter);
38  // assign the right sign
39  double sign = ((lineDirection.cross(momentum)).dot(decVec) < 0.) ? -1. : 1.;
40  lposition[eBoundLoc0] *= sign;
41  return Result<Vector2D>::success(lposition);
42 }
43 
44 inline std::string LineSurface::name() const {
45  return "Acts::LineSurface";
46 }
47 
48 inline RotationMatrix3D LineSurface::referenceFrame(
49  const GeometryContext& gctx, const Vector3D& /*unused*/,
50  const Vector3D& momentum) const {
51  RotationMatrix3D mFrame;
52  const auto& tMatrix = transform(gctx).matrix();
53  Vector3D measY(tMatrix(0, 2), tMatrix(1, 2), tMatrix(2, 2));
54  Vector3D measX(measY.cross(momentum).normalized());
55  Vector3D measDepth(measX.cross(measY));
56  // assign the columnes
57  mFrame.col(0) = measX;
58  mFrame.col(1) = measY;
59  mFrame.col(2) = measDepth;
60  // return the rotation matrix
61  return mFrame;
62 }
63 
64 inline double LineSurface::pathCorrection(const GeometryContext& /*unused*/,
65  const Vector3D& /*pos*/,
66  const Vector3D& /*mom*/) const {
67  return 1.;
68 }
69 
70 inline Vector3D LineSurface::binningPosition(const GeometryContext& gctx,
71  BinningValue /*bValue*/) const {
72  return center(gctx);
73 }
74 
75 inline Vector3D LineSurface::normal(const GeometryContext& gctx,
76  const Vector2D& /*lpos*/) const {
77  const auto& tMatrix = transform(gctx).matrix();
78  return Vector3D(tMatrix(0, 2), tMatrix(1, 2), tMatrix(2, 2));
79 }
80 
81 inline const SurfaceBounds& LineSurface::bounds() const {
82  if (m_bounds) {
83  return (*m_bounds.get());
84  }
85  return s_noBounds;
86 }
87 
89  const GeometryContext& gctx, const Vector3D& position,
90  const Vector3D& direction, const BoundaryCheck& bcheck) const {
91  // following nominclature found in header file and doxygen documentation
92  // line one is the straight track
93  const Vector3D& ma = position;
94  const Vector3D& ea = direction;
95  // line two is the line surface
96  const auto& tMatrix = transform(gctx).matrix();
97  Vector3D mb = tMatrix.block<3, 1>(0, 3).transpose();
98  Vector3D eb = tMatrix.block<3, 1>(0, 2).transpose();
99  // now go ahead and solve for the closest approach
100  Vector3D mab(mb - ma);
101  double eaTeb = ea.dot(eb);
102  double denom = 1 - eaTeb * eaTeb;
103  // validity parameter
104  Intersection3D::Status status = Intersection3D::Status::unreachable;
105  if (denom * denom > s_onSurfaceTolerance * s_onSurfaceTolerance) {
106  double u = (mab.dot(ea) - mab.dot(eb) * eaTeb) / denom;
107  // Check if we are on the surface already
108  status = (u * u < s_onSurfaceTolerance * s_onSurfaceTolerance)
109  ? Intersection3D::Status::onSurface
110  : Intersection3D::Status::reachable;
111  Vector3D result = (ma + u * ea);
112  // Evaluate the boundary check if requested
113  // m_bounds == nullptr prevents unecessary calulations for PerigeeSurface
114  if (bcheck and m_bounds) {
115  // At closest approach: check inside R or and inside Z
116  const Vector3D vecLocal(result - mb);
117  double cZ = vecLocal.dot(eb);
118  double hZ =
119  m_bounds->get(LineBounds::eHalfLengthZ) + s_onSurfaceTolerance;
120  if ((cZ * cZ > hZ * hZ) or
121  ((vecLocal - cZ * eb).norm() >
122  m_bounds->get(LineBounds::eR) + s_onSurfaceTolerance)) {
123  status = Intersection3D::Status::missed;
124  }
125  }
126  return {Intersection3D(result, u, status), this};
127  }
128  // return a false intersection
129  return {Intersection3D(position, std::numeric_limits<double>::max(), status),
130  this};
131 }
132 
133 inline void LineSurface::initJacobianToGlobal(const GeometryContext& gctx,
134  BoundToFreeMatrix& jacobian,
135  const Vector3D& position,
136  const Vector3D& direction,
137  const BoundVector& pars) const {
138  // The trigonometry required to convert the direction to spherical
139  // coordinates and then compute the sines and cosines again can be
140  // surprisingly expensive from a performance point of view.
141  //
142  // Here, we can avoid it because the direction is by definition a unit
143  // vector, with the following coordinate conversions...
144  const double x = direction(0); // == cos(phi) * sin(theta)
145  const double y = direction(1); // == sin(phi) * sin(theta)
146  const double z = direction(2); // == cos(theta)
147 
148  // ...which we can invert to directly get the sines and cosines:
149  const double cos_theta = z;
150  const double sin_theta = sqrt(x * x + y * y);
151  const double inv_sin_theta = 1. / sin_theta;
152  const double cos_phi = x * inv_sin_theta;
153  const double sin_phi = y * inv_sin_theta;
154  // retrieve the reference frame
155  const auto rframe = referenceFrame(gctx, position, direction);
156  // the local error components - given by the reference frame
157  jacobian.topLeftCorner<3, 2>() = rframe.topLeftCorner<3, 2>();
158  // the time component
159  jacobian(3, eBoundTime) = 1;
160  // the momentum components
161  jacobian(4, eBoundPhi) = (-sin_theta) * sin_phi;
162  jacobian(4, eBoundTheta) = cos_theta * cos_phi;
163  jacobian(5, eBoundPhi) = sin_theta * cos_phi;
164  jacobian(5, eBoundTheta) = cos_theta * sin_phi;
165  jacobian(6, eBoundTheta) = (-sin_theta);
166  jacobian(7, eBoundQOverP) = 1;
167 
168  // the projection of direction onto ref frame normal
169  double ipdn = 1. / direction.dot(rframe.col(2));
170  // build the cross product of d(D)/d(eBoundPhi) components with y axis
171  auto dDPhiY =
172  rframe.block<3, 1>(0, 1).cross(jacobian.block<3, 1>(4, eBoundPhi));
173  // and the same for the d(D)/d(eTheta) components
174  auto dDThetaY =
175  rframe.block<3, 1>(0, 1).cross(jacobian.block<3, 1>(4, eBoundTheta));
176  // and correct for the x axis components
177  dDPhiY -= rframe.block<3, 1>(0, 0) * (rframe.block<3, 1>(0, 0).dot(dDPhiY));
178  dDThetaY -=
179  rframe.block<3, 1>(0, 0) * (rframe.block<3, 1>(0, 0).dot(dDThetaY));
180  // set the jacobian components for global d/ phi/Theta
181  jacobian.block<3, 1>(0, eBoundPhi) = dDPhiY * pars[eBoundLoc0] * ipdn;
182  jacobian.block<3, 1>(0, eBoundTheta) = dDThetaY * pars[eBoundLoc0] * ipdn;
183 }
184 
185 inline BoundRowVector LineSurface::derivativeFactors(
186  const GeometryContext& gctx, const Vector3D& position,
187  const Vector3D& direction, const RotationMatrix3D& rft,
188  const BoundToFreeMatrix& jacobian) const {
189  // the vector between position and center
190  ActsRowVectorD<3> pc = (position - center(gctx)).transpose();
191  // the longitudinal component vector (alogn local z)
192  ActsRowVectorD<3> locz = rft.block<1, 3>(1, 0);
193  // build the norm vector comonent by subtracting the longitudinal one
194  double long_c = locz * direction;
195  ActsRowVectorD<3> norm_vec = direction.transpose() - long_c * locz;
196  // calculate the s factors for the dependency on X
197  const BoundRowVector s_vec =
198  norm_vec * jacobian.topLeftCorner<3, eBoundSize>();
199  // calculate the d factors for the dependency on Tx
200  const BoundRowVector d_vec = locz * jacobian.block<3, eBoundSize>(4, 0);
201  // normalisation of normal & longitudinal components
202  double norm = 1. / (1. - long_c * long_c);
203  // create a matrix representation
205  long_mat.colwise() += locz.transpose();
206  // build the combined normal & longitudinal components
207  return (norm * (s_vec - pc * (long_mat * d_vec.asDiagonal() -
208  jacobian.block<3, eBoundSize>(4, 0))));
209 }
210 
211 inline AlignmentRowVector LineSurface::alignmentToPathDerivative(
212  const GeometryContext& gctx, const RotationMatrix3D& rotToLocalZAxis,
213  const Vector3D& position, const Vector3D& direction) const {
214  // The vector between position and center
215  const ActsRowVector<double, 3> pcRowVec =
216  (position - center(gctx)).transpose();
217  // The local frame transform
218  const auto& sTransform = transform(gctx);
219  const auto& rotation = sTransform.rotation();
220  // The local frame z axis
221  const Vector3D localZAxis = rotation.col(2);
222  // The local z coordinate
223  const double localZ = pcRowVec * localZAxis;
224 
225  // Cosine of angle between momentum direction and local frame z axis
226  const double dirZ = localZAxis.dot(direction);
227  const double norm = 1. / (1. - dirZ * dirZ);
228  // Initialize the derivative of propagation path w.r.t. local frame
229  // translation (origin) and rotation
230  AlignmentRowVector alignToPath = AlignmentRowVector::Zero();
231  alignToPath.segment<3>(eAlignmentCenter0) =
232  norm * (direction.transpose() - dirZ * localZAxis.transpose());
233  alignToPath.segment<3>(eAlignmentRotation0) =
234  norm * (dirZ * pcRowVec + localZ * direction.transpose()) *
235  rotToLocalZAxis;
236 
237  return alignToPath;
238 }
239 
241 LineSurface::localCartesianToBoundLocalDerivative(
242  const GeometryContext& gctx, const Vector3D& position) const {
243  using VectorHelpers::phi;
244  // The local frame transform
245  const auto& sTransform = transform(gctx);
246  // calculate the transformation to local coorinates
247  const Vector3D localPos = sTransform.inverse() * position;
248  const double lphi = phi(localPos);
249  const double lcphi = std::cos(lphi);
250  const double lsphi = std::sin(lphi);
251  LocalCartesianToBoundLocalMatrix loc3DToLocBound =
252  LocalCartesianToBoundLocalMatrix::Zero();
253  loc3DToLocBound << lcphi, lsphi, 0, 0, 0, 1;
254 
255  return loc3DToLocBound;
256 }