ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ImpactPointEstimator.ipp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ImpactPointEstimator.ipp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2019 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 
14 
15 template <typename input_track_t, typename propagator_t,
16  typename propagator_options_t>
18  input_track_t, propagator_t,
19  propagator_options_t>::calculate3dDistance(const GeometryContext& gctx,
20  const BoundParameters& trkParams,
21  const Vector3D& vtxPos) const {
23  Vector3D momDir;
24 
25  auto res = getDistanceAndMomentum(gctx, trkParams, vtxPos, deltaR, momDir);
26 
27  if (!res.ok()) {
28  return res.error();
29  }
30 
31  // Return distance
32  return deltaR.norm();
33 }
34 
35 template <typename input_track_t, typename propagator_t,
36  typename propagator_options_t>
40  const Acts::MagneticFieldContext& mctx,
41  const BoundParameters& trkParams,
42  const Vector3D& vtxPos) const {
44  Vector3D momDir;
45 
46  auto res = getDistanceAndMomentum(gctx, trkParams, vtxPos, deltaR, momDir);
47 
48  if (!res.ok()) {
49  return res.error();
50  }
51 
52  // Normalize deltaR
53  deltaR.normalize();
54 
55  // corrected deltaR for small deviations from orthogonality
56  Vector3D corrDeltaR = deltaR - (deltaR.dot(momDir)) * momDir;
57 
58  // get perpendicular direction vector
59  Vector3D perpDir = momDir.cross(corrDeltaR);
60 
61  Transform3D thePlane;
62  // rotation matrix
63  thePlane.matrix().block(0, 0, 3, 1) = corrDeltaR;
64  thePlane.matrix().block(0, 1, 3, 1) = perpDir;
65  thePlane.matrix().block(0, 2, 3, 1) = momDir;
66  // translation
67  thePlane.matrix().block(0, 3, 3, 1) = vtxPos;
68 
69  std::shared_ptr<PlaneSurface> planeSurface =
70  Surface::makeShared<PlaneSurface>(
71  std::make_shared<Transform3D>(thePlane));
72 
73  // Create propagator options
74  propagator_options_t pOptions(gctx, mctx);
75  pOptions.direction = backward;
76 
77  // Do the propagation to linPointPos
78  auto result = m_cfg.propagator->propagate(trkParams, *planeSurface, pOptions);
79  if (result.ok()) {
80  return std::move((*result).endParameters);
81  } else {
82  return result.error();
83  }
84 }
85 
86 template <typename input_track_t, typename propagator_t,
87  typename propagator_options_t>
91  const BoundParameters* trkParams,
92  const Vector3D& vertexPos) const {
93  if (trkParams == nullptr) {
94  return VertexingError::EmptyInput;
95  }
96 
97  // surface rotation
98  RotationMatrix3D myRotation =
99  trkParams->referenceSurface().transform(gctx).rotation();
100  // Surface translation
101  Vector3D myTranslation =
102  trkParams->referenceSurface().transform(gctx).translation();
103 
104  // x and y direction of plane
105  Vector3D xDirPlane = myRotation.col(0);
106  Vector3D yDirPlane = myRotation.col(1);
107 
108  // transform vertex position in local plane reference frame
109  Vector3D vertexLocPlane = vertexPos - myTranslation;
110 
111  // local x/y vertex position
112  Vector2D vertexLocXY{vertexLocPlane.dot(xDirPlane),
113  vertexLocPlane.dot(yDirPlane)};
114 
115  // track covariance
116  auto cov = trkParams->covariance();
117  ActsSymMatrixD<2> myWeightXY = cov->block<2, 2>(0, 0).inverse();
118 
119  // 2-dim residual
120  Vector2D myXYpos =
121  Vector2D(trkParams->parameters()[eX], trkParams->parameters()[eY]) -
122  vertexLocXY;
123 
124  // return chi2
125  return myXYpos.dot(myWeightXY * myXYpos);
126 }
127 
128 template <typename input_track_t, typename propagator_t,
129  typename propagator_options_t>
131  input_track_t, propagator_t,
132  propagator_options_t>::performNewtonApproximation(const Vector3D& trkPos,
133  const Vector3D& vtxPos,
134  double phi, double theta,
135  double r) const {
136  double sinNewPhi = -std::sin(phi);
137  double cosNewPhi = std::cos(phi);
138 
139  int nIter = 0;
140  bool hasConverged = false;
141 
142  double cotTheta = 1. / std::tan(theta);
143 
144  // start iteration until convergence or max iteration reached
145  while (!hasConverged && nIter < m_cfg.maxIterations) {
146  double x0 = trkPos.x();
147  double y0 = trkPos.y();
148  double z0 = trkPos.z();
149 
150  double xc = vtxPos.x();
151  double yc = vtxPos.y();
152  double zc = vtxPos.z();
153 
154  double derivative = (x0 - xc) * (-r * cosNewPhi) +
155  (y0 - yc) * r * sinNewPhi +
156  (z0 - zc - r * phi * cotTheta) * (-r * cotTheta);
157  double secDerivative = r * (-(x0 - xc) * sinNewPhi - (y0 - yc) * cosNewPhi +
158  r * cotTheta * cotTheta);
159 
160  if (secDerivative < 0.) {
161  return VertexingError::NumericFailure;
162  }
163 
164  double deltaPhi = -derivative / secDerivative;
165 
166  phi += deltaPhi;
167  sinNewPhi = -std::sin(phi);
168  cosNewPhi = std::cos(phi);
169 
170  nIter += 1;
171 
172  if (std::abs(deltaPhi) < m_cfg.precision) {
173  hasConverged = true;
174  }
175  } // end while loop
176 
177  if (!hasConverged) {
178  // max iterations reached but did not converge
179  return VertexingError::NotConverged;
180  }
181  return phi;
182 }
183 
184 template <typename input_track_t, typename propagator_t,
185  typename propagator_options_t>
189  const BoundParameters& trkParams,
190  const Vector3D& vtxPos, Vector3D& deltaR,
191  Vector3D& momDir) const {
192  Vector3D trkSurfaceCenter = trkParams.referenceSurface().center(gctx);
193 
194  double d0 = trkParams.parameters()[ParID_t::eLOC_D0];
195  double z0 = trkParams.parameters()[ParID_t::eLOC_Z0];
196  double phi = trkParams.parameters()[ParID_t::ePHI];
197  double theta = trkParams.parameters()[ParID_t::eTHETA];
198  double qOvP = trkParams.parameters()[ParID_t::eQOP];
199 
200  double cotTheta = 1. / std::tan(theta);
201 
202  // get B-field z-component at current position
203  double bZ = m_cfg.bField.getField(trkSurfaceCenter)[eZ];
204 
205  // The radius
206  double r;
207  // Curvature is infinite w/o b field
208  if (bZ == 0. || std::abs(qOvP) < m_cfg.minQoP) {
209  r = m_cfg.maxRho;
210  } else {
211  // signed(!) r
212  r = std::sin(theta) * (1. / qOvP) / bZ;
213  }
214 
215  Vector3D vec0 = trkSurfaceCenter + Vector3D(-(d0 - r) * std::sin(phi),
216  (d0 - r) * std::cos(phi),
217  z0 + r * phi * cotTheta);
218 
219  // Perform newton approximation method
220  // this will change the value of phi
221  auto res = performNewtonApproximation(vec0, vtxPos, phi, theta, r);
222  if (!res.ok()) {
223  return res.error();
224  }
225  // Set new phi value
226  phi = *res;
227 
228  // Set momentum direction
229  momDir = Vector3D(std::sin(theta) * std::cos(phi),
230  std::sin(phi) * std::sin(theta), std::cos(theta));
231 
232  // point of closest approach in 3D
233  Vector3D pointCA3d =
234  vec0 + r * Vector3D(-std::sin(phi), std::cos(phi), -cotTheta * phi);
235 
236  // Set deltaR
237  deltaR = pointCA3d - vtxPos;
238 
239  return {};
240 }
241 
242 template <typename input_track_t, typename propagator_t,
243  typename propagator_options_t>
247  const Vertex<input_track_t>& vtx,
248  const GeometryContext& gctx,
249  const Acts::MagneticFieldContext& mctx) const {
250  // estimating the d0 and its significance by propagating the trajectory state
251  // towards
252  // the vertex position. By this time the vertex should NOT contain this
253  // trajectory anymore
254  const std::shared_ptr<PerigeeSurface> perigeeSurface =
255  Surface::makeShared<PerigeeSurface>(vtx.position());
256 
257  // Create propagator options
258  propagator_options_t pOptions(gctx, mctx);
259  pOptions.direction = backward;
260 
261  // Do the propagation to linPoint
262  auto result = m_cfg.propagator->propagate(track, *perigeeSurface, pOptions);
263 
264  if (!result.ok()) {
265  return result.error();
266  }
267 
268  const auto& propRes = *result;
269  const auto& params = propRes.endParameters->parameters();
270  const double d0 = params[ParID_t::eLOC_D0];
271  const double z0 = params[ParID_t::eLOC_Z0];
272  const double phi = params[ParID_t::ePHI];
273  const double theta = params[ParID_t::eTHETA];
274 
275  ActsSymMatrixD<2> vrtXYCov = vtx.covariance().template block<2, 2>(0, 0);
276 
277  // Covariance of perigee parameters after propagation to perigee surface
278  const auto& perigeeCov = *(propRes.endParameters->covariance());
279 
280  ActsVectorD<2> d0JacXY(-std::sin(phi), std::cos(phi));
281 
282  ImpactParametersAndSigma newIPandSigma;
283 
284  newIPandSigma.IPd0 = d0;
285  double d0_PVcontrib = d0JacXY.transpose() * (vrtXYCov * d0JacXY);
286  if (d0_PVcontrib >= 0) {
287  newIPandSigma.sigmad0 = std::sqrt(
288  d0_PVcontrib + perigeeCov(ParID_t::eLOC_D0, ParID_t::eLOC_D0));
289  newIPandSigma.PVsigmad0 = std::sqrt(d0_PVcontrib);
290  } else {
291  newIPandSigma.sigmad0 =
292  std::sqrt(perigeeCov(ParID_t::eLOC_D0, ParID_t::eLOC_D0));
293  newIPandSigma.PVsigmad0 = 0;
294  }
295 
296  ActsSymMatrixD<2> covPerigeeZ0Theta;
297  covPerigeeZ0Theta(0, 0) = perigeeCov(ParID_t::eLOC_Z0, ParID_t::eLOC_Z0);
298  covPerigeeZ0Theta(0, 1) = perigeeCov(ParID_t::eLOC_Z0, ParID_t::eTHETA);
299  covPerigeeZ0Theta(1, 0) = perigeeCov(ParID_t::eTHETA, ParID_t::eLOC_Z0);
300  covPerigeeZ0Theta(1, 1) = perigeeCov(ParID_t::eTHETA, ParID_t::eTHETA);
301 
302  double vtxZZCov = vtx.covariance()(eZ, eZ);
303 
304  ActsVectorD<2> z0JacZ0Theta(std::sin(theta), z0 * std::cos(theta));
305 
306  if (vtxZZCov >= 0) {
307  newIPandSigma.IPz0SinTheta = z0 * std::sin(theta);
308  newIPandSigma.sigmaz0SinTheta = std::sqrt(
309  z0JacZ0Theta.transpose() * (covPerigeeZ0Theta * z0JacZ0Theta) +
310  std::sin(theta) * vtxZZCov * std::sin(theta));
311 
312  newIPandSigma.PVsigmaz0SinTheta =
313  std::sqrt(std::sin(theta) * vtxZZCov * std::sin(theta));
314  newIPandSigma.IPz0 = z0;
315  newIPandSigma.sigmaz0 =
316  std::sqrt(vtxZZCov + perigeeCov(ParID_t::eLOC_Z0, ParID_t::eLOC_Z0));
317  newIPandSigma.PVsigmaz0 = std::sqrt(vtxZZCov);
318  } else {
319  // Remove contribution from PV
320  newIPandSigma.IPz0SinTheta = z0 * std::sin(theta);
321  double sigma2z0sinTheta =
322  (z0JacZ0Theta.transpose() * (covPerigeeZ0Theta * z0JacZ0Theta));
323  newIPandSigma.sigmaz0SinTheta = std::sqrt(sigma2z0sinTheta);
324  newIPandSigma.PVsigmaz0SinTheta = 0;
325 
326  newIPandSigma.IPz0 = z0;
327  newIPandSigma.sigmaz0 =
328  std::sqrt(perigeeCov(ParID_t::eLOC_Z0, ParID_t::eLOC_Z0));
329  newIPandSigma.PVsigmaz0 = 0;
330  }
331 
332  return newIPandSigma;
333 }