Estimate the full track parameters from three space points.
This method is based on the conformal map transformation. It estimates the full bound track parameters, i.e. (loc0, loc1, phi, theta, q/p, t) at the bottom space point. The bottom space is assumed to be the first element in the range defined by the iterators. The magnetic field (which might be along any direction) is also necessary for the momentum estimation.
96 {
97
98 size_t num_sp = std::distance(spBegin, spEnd);
99 if (num_sp != 3) {
100 std::cout << "ERROR::less than 3 point provided" << std::endl;
101 return std::nullopt;
102 }
103
104
105 Acts::ActsScalar b_field_in_tesla = bField.norm() / Acts::UnitConstants::T;
106 Acts::ActsScalar b_field_min_in_tesla = bFieldMin / Acts::UnitConstants::T;
107
108 if (b_field_in_tesla < b_field_min_in_tesla) {
109
110
111 std::cout << "The magnetic field at the bottom space point: B = "
112 << b_field_in_tesla
113 << " T is smaller than |B|_min = " << b_field_min_in_tesla
114 << " T. Estimation is not performed." << std::endl;
115 return std::nullopt;
116 }
117
118
119 std::array<Acts::Vector3, 3> sp_global_positions = {
120 Acts::Vector3::Zero(), Acts::Vector3::Zero(), Acts::Vector3::Zero()};
121
122
123 for (size_t isp = 0; isp < 3; ++isp) {
124 spacepoint_iterator_t it = std::next(spBegin, isp);
125 if (*it == nullptr) {
126 std::cout << "Empty space point found. This should not happen."
127 << std::endl;
128 return std::nullopt;
129 }
130 const auto& sp = *it;
131 sp_global_positions[isp] = Acts::Vector3(sp->x(), sp->y(), sp->z());
132 }
133
134
135
136
137
138
139 Acts::Vector3 rel_vec = sp_global_positions[1] - sp_global_positions[0];
140 Acts::Vector3 new_z_axis = bField.normalized();
141 Acts::Vector3 new_y_axis = new_z_axis.cross(rel_vec).normalized();
142 Acts::Vector3 new_x_axis = new_y_axis.cross(new_z_axis);
143 Acts::RotationMatrix3 rotation;
144 rotation.col(0) = new_x_axis;
145 rotation.col(1) = new_y_axis;
146 rotation.col(2) = new_z_axis;
147
148 Acts::Translation3 trans(sp_global_positions[0]);
149
150 Acts::Transform3 transform(trans * rotation);
151
152
153 Acts::Vector3 local1 = transform.inverse() * sp_global_positions[1];
154 Acts::Vector3 local2 = transform.inverse() * sp_global_positions[2];
155
156
157 auto uv_transform = [](const Acts::Vector3& local) -> Acts::Vector2 {
158 Acts::Vector2 uv;
159 Acts::ActsScalar denominator =
160 local.x() * local.x() + local.y() * local.y();
161 uv.x() = local.x() / denominator;
162 uv.y() = local.y() / denominator;
163 return uv;
164 };
165
166 Acts::Vector2 uv1 = uv_transform(local1);
167 Acts::Vector2 uv2 = uv_transform(local2);
168
169
170
171 Acts::ActsScalar a = (uv2.y() - uv1.y()) / (uv2.x() - uv1.x());
172 Acts::ActsScalar b = uv2.y() - a * uv2.x();
173
174 Acts::ActsScalar rho = -2.0 * b / std::hypot(1., a);
175
176
177 Acts::ActsScalar rn = local2.x() * local2.x() + local2.y() * local2.y();
178
179 Acts::ActsScalar inv_tan_theta =
180 local2.z() * std::sqrt(1. / rn) / (1. + rho * rho * rn);
181
182
183 Acts::Vector3 trans_direction(1., a, std::hypot(1, a) * inv_tan_theta);
184
185 [[maybe_unused]] Acts::Vector3 direction =
186 rotation * trans_direction.normalized();
187
188
189 Acts::BoundVector params = Acts::BoundVector::Zero();
190
191
192
193
194
195 Acts::Vector3 bottom_local_pos = Tp.inverse() * sp_global_positions[0];
196
197
198 params[Acts::eBoundLoc0] = bottom_local_pos.x();
199 params[Acts::eBoundLoc1] = bottom_local_pos.y();
200
201
202
203 Acts::ActsScalar q_over_pt =
204 rho * (Acts::UnitConstants::m) / (0.3 * b_field_in_tesla);
205
206 params[Acts::eBoundQOverP] = q_over_pt / std::hypot(1., inv_tan_theta);
207
208
209
210 Acts::ActsScalar p_in_ge_v = std::abs(1.0 / params[Acts::eBoundQOverP]);
211 Acts::ActsScalar pz_in_ge_v = 1.0 / std::abs(q_over_pt) * inv_tan_theta;
212 Acts::ActsScalar mass_in_ge_v = mass / Acts::UnitConstants::GeV;
213
214
215 Acts::ActsScalar v = p_in_ge_v / std::hypot(p_in_ge_v, mass_in_ge_v);
216 Acts::ActsScalar vz = pz_in_ge_v / std::hypot(p_in_ge_v, mass_in_ge_v);
217
218
219 Acts::ActsScalar pathz = sp_global_positions[0].dot(bField) / bField.norm();
220
221
222 if (pathz != 0) {
223 params[Acts::eBoundTime] = pathz / vz;
224 } else {
225 params[Acts::eBoundTime] = sp_global_positions[0].norm() / v;
226 }
227
228 return params;
229 }