File indexing completed on 2025-01-18 09:12:19
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include "Acts/Plugins/EDM4hep/EDM4hepUtil.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/Definitions/Units.hpp"
0014 #include "Acts/EventData/MultiTrajectory.hpp"
0015 #include "Acts/EventData/MultiTrajectoryHelpers.hpp"
0016 #include "Acts/Propagator/detail/CovarianceEngine.hpp"
0017 #include "Acts/Propagator/detail/JacobianEngine.hpp"
0018
0019 #include <numbers>
0020
0021 #include <edm4hep/EDM4hepVersion.h>
0022 #include <edm4hep/MCParticle.h>
0023 #include <edm4hep/MutableSimTrackerHit.h>
0024 #include <edm4hep/SimTrackerHit.h>
0025 #include <edm4hep/TrackState.h>
0026
0027 namespace Acts::EDM4hepUtil {
0028 namespace detail {
0029
0030 ActsSquareMatrix<6> jacobianToEdm4hep(double theta, double qOverP, double Bz) {
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055 ActsSquareMatrix<6> J;
0056 J.setIdentity();
0057 double sinTheta = std::sin(theta);
0058 J(3, 3) = -1.0 / (sinTheta * sinTheta);
0059 J(4, 4) = Bz / sinTheta;
0060 J(4, 3) = -Bz * qOverP * std::cos(theta) /
0061 (sinTheta * sinTheta);
0062 return J;
0063 }
0064
0065 ActsSquareMatrix<6> jacobianFromEdm4hep(double tanLambda, double omega,
0066 double Bz) {
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090 ActsSquareMatrix<6> J;
0091 J.setIdentity();
0092 J(3, 3) = -1 / (tanLambda * tanLambda + 1);
0093 J(4, 3) = -1 * omega * tanLambda /
0094 (Bz * std::pow(tanLambda * tanLambda + 1, 3. / 2.));
0095 J(4, 4) = 1 / (Bz * std::hypot(tanLambda, 1));
0096
0097 return J;
0098 }
0099
0100 void packCovariance(const ActsSquareMatrix<6>& from, float* to) {
0101 for (int i = 0; i < from.rows(); i++) {
0102 for (int j = 0; j <= i; j++) {
0103 std::size_t k = (i + 1) * i / 2 + j;
0104 to[k] = from(i, j);
0105 }
0106 }
0107 }
0108
0109 void unpackCovariance(const float* from, ActsSquareMatrix<6>& to) {
0110 auto k = [](std::size_t i, std::size_t j) { return (i + 1) * i / 2 + j; };
0111 for (int i = 0; i < to.rows(); i++) {
0112 for (int j = 0; j < to.cols(); j++) {
0113 to(i, j) = from[j <= i ? k(i, j) : k(j, i)];
0114 }
0115 }
0116 }
0117
0118 Parameters convertTrackParametersToEdm4hep(const Acts::GeometryContext& gctx,
0119 double Bz,
0120 const BoundTrackParameters& params) {
0121 Acts::Vector3 global = params.referenceSurface().localToGlobal(
0122 gctx, params.parameters().template head<2>(), params.direction());
0123
0124 std::shared_ptr<const Acts::Surface> refSurface =
0125 params.referenceSurface().getSharedPtr();
0126 Acts::BoundVector targetPars = params.parameters();
0127 std::optional<Acts::BoundSquareMatrix> targetCov = params.covariance();
0128
0129
0130
0131
0132 if (dynamic_cast<const Acts::PerigeeSurface*>(refSurface.get()) == nullptr) {
0133 refSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(global);
0134
0135
0136
0137
0138
0139 auto perigeeParams = Acts::detail::boundToBoundConversion(
0140 gctx, params, *refSurface, Vector3{0, 0, Bz})
0141 .value();
0142 targetPars = perigeeParams.parameters();
0143 targetCov = perigeeParams.covariance();
0144 }
0145
0146 Parameters result;
0147 result.surface = refSurface;
0148
0149
0150 if (targetCov) {
0151 Acts::ActsSquareMatrix<6> J = jacobianToEdm4hep(
0152 targetPars[eBoundTheta], targetPars[eBoundQOverP], Bz);
0153 result.covariance = J * targetCov.value() * J.transpose();
0154 }
0155
0156 result.values[0] = targetPars[Acts::eBoundLoc0];
0157 result.values[1] = targetPars[Acts::eBoundLoc1];
0158 result.values[2] = targetPars[Acts::eBoundPhi];
0159 result.values[3] =
0160 std::tan(std::numbers::pi / 2. - targetPars[Acts::eBoundTheta]);
0161 result.values[4] = targetPars[Acts::eBoundQOverP] /
0162 std::sin(targetPars[Acts::eBoundTheta]) * Bz;
0163 result.values[5] = targetPars[Acts::eBoundTime];
0164
0165 result.particleHypothesis = params.particleHypothesis();
0166
0167 return result;
0168 }
0169
0170 BoundTrackParameters convertTrackParametersFromEdm4hep(
0171 double Bz, const Parameters& params) {
0172 BoundVector targetPars;
0173
0174 ActsSquareMatrix<6> J =
0175 jacobianFromEdm4hep(params.values[3], params.values[4], Bz);
0176
0177 std::optional<BoundMatrix> cov;
0178 if (params.covariance.has_value()) {
0179 cov = J * params.covariance.value() * J.transpose();
0180 }
0181
0182 targetPars[eBoundLoc0] = params.values[0];
0183 targetPars[eBoundLoc1] = params.values[1];
0184 targetPars[eBoundPhi] = params.values[2];
0185 targetPars[eBoundTheta] = std::numbers::pi / 2. - std::atan(params.values[3]);
0186 targetPars[eBoundQOverP] =
0187 params.values[4] * std::sin(targetPars[eBoundTheta]) / Bz;
0188 targetPars[eBoundTime] = params.values[5];
0189
0190 return {params.surface, targetPars, cov, params.particleHypothesis};
0191 }
0192
0193 }
0194
0195 #if EDM4HEP_VERSION_MAJOR >= 1 || \
0196 (EDM4HEP_VERSION_MAJOR == 0 && EDM4HEP_VERSION_MINOR == 99)
0197 edm4hep::MCParticle getParticle(const edm4hep::SimTrackerHit& hit) {
0198 return hit.getParticle();
0199 }
0200
0201 void setParticle(edm4hep::MutableSimTrackerHit& hit,
0202 const edm4hep::MCParticle& particle) {
0203 hit.setParticle(particle);
0204 }
0205 #else
0206 edm4hep::MCParticle getParticle(const edm4hep::SimTrackerHit& hit) {
0207 return hit.getMCParticle();
0208 }
0209
0210 void setParticle(edm4hep::MutableSimTrackerHit& hit,
0211 const edm4hep::MCParticle& particle) {
0212 hit.setMCParticle(particle);
0213 }
0214 #endif
0215
0216 }