ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HelicalTrackLinearizer.ipp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file HelicalTrackLinearizer.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 
10 
11 template <typename propagator_t, typename propagator_options_t>
14  const BoundParameters& params, const SpacePointVector& linPoint,
16  const Acts::MagneticFieldContext& mctx) const {
17  Vector3D linPointPos = VectorHelpers::position(linPoint);
18 
19  const std::shared_ptr<PerigeeSurface> perigeeSurface =
20  Surface::makeShared<PerigeeSurface>(linPointPos);
21 
22  // Create propagator options
23  propagator_options_t pOptions(gctx, mctx);
24  pOptions.direction = backward;
25 
26  const BoundParameters* endParams = nullptr;
27  // Do the propagation to linPointPos
28  auto result = m_cfg.propagator->propagate(params, *perigeeSurface, pOptions);
29  if (result.ok()) {
30  endParams = (*result).endParameters.get();
31 
32  } else {
33  return result.error();
34  }
35 
36  BoundVector paramsAtPCA = endParams->parameters();
37  SpacePointVector positionAtPCA = SpacePointVector::Zero();
38  VectorHelpers::position(positionAtPCA) = endParams->position();
39  BoundSymMatrix parCovarianceAtPCA = *(endParams->covariance());
40 
41  if (endParams->covariance()->determinant() == 0) {
42  // Use the original parameters
43  paramsAtPCA = params.parameters();
44  VectorHelpers::position(positionAtPCA) = params.position();
45  parCovarianceAtPCA = *(params.covariance());
46  }
47 
48  // phiV and functions
49  double phiV = paramsAtPCA(ParID_t::ePHI);
50  double sinPhiV = std::sin(phiV);
51  double cosPhiV = std::cos(phiV);
52 
53  // theta and functions
54  double th = paramsAtPCA(ParID_t::eTHETA);
55  const double sinTh = std::sin(th);
56  const double tanTh = std::tan(th);
57 
58  // q over p
59  double qOvP = paramsAtPCA(ParID_t::eQOP);
60  double sgnH = (qOvP < 0.) ? -1 : 1;
61 
62  Vector3D momentumAtPCA(phiV, th, qOvP);
63 
64  // get B-field z-component at current position
65  double Bz = m_cfg.bField.getField(VectorHelpers::position(positionAtPCA))[eZ];
66 
67  double rho;
68  double rho2;
69  // Curvature is infinite w/o b field
70  if (Bz == 0. || std::abs(qOvP) < m_cfg.minQoP) {
71  rho = m_cfg.maxRho;
72  } else {
73  // signed(!) rho
74  rho2 = sinTh / (qOvP * Bz);
75  rho = sinTh * (1. / qOvP) / Bz;
76  }
77 
78  // Eq. 5.34 in Ref(1) (see .hpp)
79  double X = positionAtPCA(0) - linPointPos.x() + rho * sinPhiV;
80  double Y = positionAtPCA(1) - linPointPos.y() - rho * cosPhiV;
81  const double S2 = (X * X + Y * Y);
82  const double S = std::sqrt(S2);
83 
87  BoundVector predParamsAtPCA;
88 
89  int sgnX = (X < 0.) ? -1 : 1;
90  int sgnY = (Y < 0.) ? -1 : 1;
91 
92  double phiAtPCA;
93  if (std::abs(X) > std::abs(Y)) {
94  phiAtPCA = sgnH * sgnX * std::acos(-sgnH * Y / S);
95  } else {
96  phiAtPCA = std::asin(sgnH * X / S);
97  if ((sgnH * sgnY) > 0) {
98  phiAtPCA = sgnH * sgnX * M_PI - phiAtPCA;
99  }
100  }
101 
102  // Eq. 5.33 in Ref(1) (see .hpp)
103  predParamsAtPCA[0] = rho - sgnH * S;
104  predParamsAtPCA[1] =
105  positionAtPCA[eZ] - linPointPos.z() + rho * (phiV - phiAtPCA) / tanTh;
106  predParamsAtPCA[2] = phiAtPCA;
107  predParamsAtPCA[3] = th;
108  predParamsAtPCA[4] = qOvP;
109  predParamsAtPCA[5] = 0.;
110 
111  // Fill position jacobian (D_k matrix), Eq. 5.36 in Ref(1)
112  SpacePointToBoundMatrix positionJacobian;
113  positionJacobian.setZero();
114  // First row
115  positionJacobian(0, 0) = -sgnH * X / S;
116  positionJacobian(0, 1) = -sgnH * Y / S;
117 
118  const double S2tanTh = S2 * tanTh;
119 
120  // Second row
121  positionJacobian(1, 0) = rho * Y / S2tanTh;
122  positionJacobian(1, 1) = -rho * X / S2tanTh;
123  positionJacobian(1, 2) = 1.;
124 
125  // Third row
126  positionJacobian(2, 0) = -Y / S2;
127  positionJacobian(2, 1) = X / S2;
128 
129  // TODO: include timing in track linearization
130  // Last row
131  positionJacobian(5, 3) = 1;
132 
133  // Fill momentum jacobian (E_k matrix), Eq. 5.37 in Ref(1)
134  ActsMatrixD<eBoundParametersSize, 3> momentumJacobian;
135  momentumJacobian.setZero();
136 
137  double R = X * cosPhiV + Y * sinPhiV;
138  double Q = X * sinPhiV - Y * cosPhiV;
139  double dPhi = phiAtPCA - phiV;
140 
141  // First row
142  momentumJacobian(0, 0) = -sgnH * rho * R / S;
143 
144  double qOvSred = 1 - sgnH * Q / S;
145 
146  momentumJacobian(0, 1) = qOvSred * rho / tanTh;
147  momentumJacobian(0, 2) = -qOvSred * rho / qOvP;
148 
149  const double rhoOverS2 = rho / S2;
150 
151  // Second row
152  momentumJacobian(1, 0) = (1 - rhoOverS2 * Q) * rho / tanTh;
153  momentumJacobian(1, 1) = (dPhi + rho * R / (S2tanTh * tanTh)) * rho;
154  momentumJacobian(1, 2) = (dPhi - rhoOverS2 * R) * rho / (qOvP * tanTh);
155 
156  // Third row
157  momentumJacobian(2, 0) = rhoOverS2 * Q;
158  momentumJacobian(2, 1) = -rho * R / S2tanTh;
159  momentumJacobian(2, 2) = rhoOverS2 * R / qOvP;
160 
161  // Last two rows:
162  momentumJacobian(3, 1) = 1.;
163  momentumJacobian(4, 2) = 1.;
164 
165  // const term F(V_0, p_0) in Talyor expansion
166  BoundVector constTerm = predParamsAtPCA - positionJacobian * positionAtPCA -
167  momentumJacobian * momentumAtPCA;
168 
169  // The parameter weight
170  ActsSymMatrixD<5> parWeight =
171  (parCovarianceAtPCA.block<5, 5>(0, 0)).inverse();
172 
173  BoundSymMatrix weightAtPCA{BoundSymMatrix::Identity()};
174  weightAtPCA.block<5, 5>(0, 0) = parWeight;
175 
176  return LinearizedTrack(paramsAtPCA, parCovarianceAtPCA, weightAtPCA, linPoint,
177  positionJacobian, momentumJacobian, positionAtPCA,
178  momentumAtPCA, constTerm);
179 }