File indexing completed on 2025-07-06 07:52:16
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include "Acts/Vertexing/FullBilloirVertexFitter.hpp"
0010
0011 #include "Acts/Definitions/TrackParametrization.hpp"
0012 #include "Acts/MagneticField/ConstantBField.hpp"
0013 #include "Acts/Surfaces/PerigeeSurface.hpp"
0014 #include "Acts/Utilities/detail/periodic.hpp"
0015 #include "Acts/Vertexing/TrackAtVertex.hpp"
0016 #include "Acts/Vertexing/VertexingError.hpp"
0017
0018 namespace {
0019
0020
0021
0022
0023 struct BilloirTrack {
0024 explicit BilloirTrack(const Acts::InputTrack& params)
0025 : originalTrack(params) {}
0026
0027 BilloirTrack(const BilloirTrack& arg) = default;
0028
0029 Acts::InputTrack originalTrack;
0030 double chi2 = 0;
0031
0032
0033 Acts::ActsMatrix<Acts::eBoundSize, Acts::eBoundSize> W;
0034 Acts::ActsMatrix<Acts::eBoundSize, 4> D;
0035 Acts::ActsMatrix<Acts::eBoundSize, 3> E;
0036 Acts::ActsSquareMatrix<3> C;
0037 Acts::ActsMatrix<4, 3> B;
0038 Acts::ActsSquareMatrix<3> Cinv;
0039 Acts::Vector3 U;
0040 Acts::ActsMatrix<4, 3> BCinv;
0041 Acts::BoundVector deltaQ;
0042 };
0043
0044
0045
0046
0047 struct BilloirVertex {
0048
0049 Acts::SquareMatrix4 A = Acts::SquareMatrix4::Zero();
0050
0051 Acts::Vector4 T = Acts::Vector4::Zero();
0052
0053 Acts::SquareMatrix4 sumBCinvBt = Acts::SquareMatrix4::Zero();
0054
0055 Acts::Vector4 sumBCinvU = Acts::Vector4::Zero();
0056 };
0057
0058 }
0059
0060 Acts::Result<Acts::Vertex> Acts::FullBilloirVertexFitter::fit(
0061 const std::vector<InputTrack>& paramVector,
0062 const VertexingOptions& vertexingOptions,
0063 MagneticFieldProvider::Cache& fieldCache) const {
0064 unsigned int nTracks = paramVector.size();
0065 double chi2 = std::numeric_limits<double>::max();
0066
0067 if (nTracks == 0) {
0068 return Vertex(Vector3(0., 0., 0.));
0069 }
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079 int ndf = 3 * nTracks - 3;
0080
0081
0082 if (nTracks < 2) {
0083 ndf = 1;
0084 }
0085
0086
0087
0088 if (vertexingOptions.useConstraintInFit) {
0089 ndf += 3;
0090 }
0091
0092 std::vector<BilloirTrack> billoirTracks;
0093 std::vector<Vector3> trackMomenta;
0094
0095 Vector4 linPoint = vertexingOptions.constraint.fullPosition();
0096 Vertex fittedVertex;
0097
0098 for (int nIter = 0; nIter < m_cfg.maxIterations; ++nIter) {
0099 billoirTracks.clear();
0100 double newChi2 = 0;
0101 BilloirVertex billoirVertex;
0102
0103 Vector3 linPointPos = VectorHelpers::position(linPoint);
0104
0105
0106 const std::shared_ptr<PerigeeSurface> perigeeSurface =
0107 Surface::makeShared<PerigeeSurface>(linPointPos);
0108
0109
0110 for (std::size_t iTrack = 0; iTrack < nTracks; ++iTrack) {
0111 const InputTrack& trackContainer = paramVector[iTrack];
0112
0113 const auto& trackParams = m_cfg.extractParameters(trackContainer);
0114
0115 auto result =
0116 m_cfg.trackLinearizer(trackParams, linPoint[3], *perigeeSurface,
0117 vertexingOptions.geoContext,
0118 vertexingOptions.magFieldContext, fieldCache);
0119 if (!result.ok()) {
0120 return result.error();
0121 }
0122
0123 const auto& linTrack = *result;
0124 const auto& parametersAtPCA = linTrack.parametersAtPCA;
0125 double d0 = parametersAtPCA[BoundIndices::eBoundLoc0];
0126 double z0 = parametersAtPCA[BoundIndices::eBoundLoc1];
0127 double phi = parametersAtPCA[BoundIndices::eBoundPhi];
0128 double theta = parametersAtPCA[BoundIndices::eBoundTheta];
0129 double qOverP = parametersAtPCA[BoundIndices::eBoundQOverP];
0130 double t0 = parametersAtPCA[BoundIndices::eBoundTime];
0131
0132
0133
0134 if (nIter == 0) {
0135 trackMomenta.push_back(Vector3(phi, theta, qOverP));
0136 }
0137
0138
0139
0140
0141 double fPhi = trackMomenta[iTrack][0];
0142 double fTheta = trackMomenta[iTrack][1];
0143 double fQOvP = trackMomenta[iTrack][2];
0144 double fTime = linPoint[FreeIndices::eFreeTime];
0145 BilloirTrack billoirTrack(trackContainer);
0146
0147 billoirTrack.deltaQ << d0, z0, phi - fPhi, theta - fTheta, qOverP - fQOvP,
0148 t0 - fTime;
0149
0150
0151 ActsMatrix<eBoundSize, 4> D = linTrack.positionJacobian;
0152
0153
0154 ActsMatrix<eBoundSize, 3> E = linTrack.momentumJacobian;
0155
0156
0157 BoundSquareMatrix W = linTrack.weightAtPCA;
0158 ActsMatrix<4, eBoundSize> DtW = D.transpose() * W;
0159 ActsMatrix<3, eBoundSize> EtW = E.transpose() * W;
0160
0161
0162 billoirTrack.D = D;
0163 billoirTrack.E = E;
0164 billoirTrack.W = W;
0165 billoirTrack.C = EtW * E;
0166 billoirTrack.B = DtW * E;
0167 billoirTrack.U = EtW * billoirTrack.deltaQ;
0168 billoirTrack.Cinv = (billoirTrack.C).inverse();
0169 billoirTrack.BCinv =
0170 billoirTrack.B * billoirTrack.Cinv;
0171
0172
0173 billoirVertex.T += DtW * billoirTrack.deltaQ;
0174 billoirVertex.A += DtW * D;
0175 billoirVertex.sumBCinvU +=
0176 billoirTrack.BCinv * billoirTrack.U;
0177 billoirVertex.sumBCinvBt +=
0178 billoirTrack.BCinv *
0179 billoirTrack.B.transpose();
0180
0181 billoirTracks.push_back(billoirTrack);
0182 }
0183
0184
0185
0186
0187 Vector4 deltaVFac = billoirVertex.T - billoirVertex.sumBCinvU;
0188
0189
0190
0191
0192 SquareMatrix4 invCovV = billoirVertex.A - billoirVertex.sumBCinvBt;
0193 if (vertexingOptions.useConstraintInFit) {
0194
0195
0196
0197
0198 Vector4 posInBilloirFrame =
0199 vertexingOptions.constraint.fullPosition() - linPoint;
0200
0201
0202
0203
0204 deltaVFac += vertexingOptions.constraint.fullCovariance().inverse() *
0205 posInBilloirFrame;
0206
0207 invCovV += vertexingOptions.constraint.fullCovariance().inverse();
0208 }
0209
0210
0211 SquareMatrix4 covV = invCovV.inverse();
0212
0213 Vector4 deltaV = covV * deltaVFac;
0214
0215
0216
0217 std::vector<std::optional<BoundSquareMatrix>> covDeltaP(nTracks);
0218
0219
0220
0221 for (std::size_t iTrack = 0; iTrack < billoirTracks.size(); ++iTrack) {
0222 auto& billoirTrack = billoirTracks[iTrack];
0223
0224 Vector3 deltaP = (billoirTrack.Cinv) *
0225 (billoirTrack.U - billoirTrack.B.transpose() * deltaV);
0226
0227
0228 trackMomenta[iTrack] += deltaP;
0229
0230
0231 const auto correctedPhiTheta = detail::normalizePhiTheta(
0232 trackMomenta[iTrack][0], trackMomenta[iTrack][1]);
0233 trackMomenta[iTrack][0] = correctedPhiTheta.first;
0234 trackMomenta[iTrack][1] = correctedPhiTheta.second;
0235
0236
0237 Acts::BoundVector diffQ = billoirTrack.deltaQ - billoirTrack.D * deltaV -
0238 billoirTrack.E * deltaP;
0239 billoirTrack.chi2 = diffQ.transpose().dot(billoirTrack.W * diffQ);
0240
0241 newChi2 += billoirTrack.chi2;
0242
0243
0244
0245
0246
0247
0248
0249
0250
0251
0252
0253 ActsMatrix<eBoundSize, 7> transMat;
0254 transMat.setZero();
0255 transMat.block<2, 2>(0, 0) = billoirTrack.D.template block<2, 2>(0, 0);
0256 transMat(1, 2) = 1.;
0257 transMat(2, 4) = 1.;
0258 transMat(3, 5) = 1.;
0259 transMat(4, 6) = 1.;
0260 transMat(5, 3) = 1.;
0261
0262
0263 ActsMatrix<4, 3> covVP = -covV * billoirTrack.BCinv;
0264
0265
0266 ActsSquareMatrix<3> covP =
0267 billoirTrack.Cinv +
0268 billoirTrack.BCinv.transpose() * covV * billoirTrack.BCinv;
0269
0270 ActsSquareMatrix<7> cov;
0271 cov.setZero();
0272 cov.block<4, 4>(0, 0) = covV;
0273 cov.block<4, 3>(0, 4) = covVP;
0274 cov.block<3, 4>(4, 0) = covVP.transpose();
0275 cov.block<3, 3>(4, 4) = covP;
0276
0277 covDeltaP[iTrack] = transMat * cov * transMat.transpose();
0278 }
0279
0280
0281 linPoint += deltaV;
0282
0283
0284
0285 if (vertexingOptions.useConstraintInFit) {
0286
0287
0288 Vector4 posInBilloirFrame =
0289 vertexingOptions.constraint.fullPosition() - linPoint;
0290
0291 newChi2 +=
0292 (posInBilloirFrame.transpose())
0293 .dot(vertexingOptions.constraint.fullCovariance().inverse() *
0294 posInBilloirFrame);
0295 }
0296
0297 if (!std::isnormal(newChi2)) {
0298 ACTS_ERROR("Encountered non-normal chi2 value during the fit.");
0299 return VertexingError::NumericFailure;
0300 }
0301
0302 if (newChi2 < chi2) {
0303 chi2 = newChi2;
0304
0305 fittedVertex.setFullPosition(linPoint);
0306 fittedVertex.setFullCovariance(covV);
0307 fittedVertex.setFitQuality(chi2, ndf);
0308
0309 std::vector<TrackAtVertex> tracksAtVertex;
0310
0311 std::shared_ptr<PerigeeSurface> perigee =
0312 Surface::makeShared<PerigeeSurface>(
0313 VectorHelpers::position(linPoint));
0314
0315 for (std::size_t iTrack = 0; iTrack < billoirTracks.size(); ++iTrack) {
0316 const auto& billoirTrack = billoirTracks[iTrack];
0317
0318
0319
0320
0321
0322
0323 BoundVector paramVec = BoundVector::Zero();
0324 paramVec[eBoundPhi] = trackMomenta[iTrack](0);
0325 paramVec[eBoundTheta] = trackMomenta[iTrack](1);
0326 paramVec[eBoundQOverP] = trackMomenta[iTrack](2);
0327 paramVec[eBoundTime] = linPoint[FreeIndices::eFreeTime];
0328 BoundTrackParameters refittedParams(
0329 perigee, paramVec, covDeltaP[iTrack],
0330 m_cfg.extractParameters(billoirTrack.originalTrack)
0331 .particleHypothesis());
0332 TrackAtVertex trackAtVertex(billoirTrack.chi2, refittedParams,
0333 billoirTrack.originalTrack);
0334 tracksAtVertex.push_back(std::move(trackAtVertex));
0335 }
0336
0337 fittedVertex.setTracksAtVertex(std::move(tracksAtVertex));
0338 }
0339 }
0340
0341 return fittedVertex;
0342 }