EDM4hep/src/EDM4hepUtil.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsPlugins/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 using namespace Acts;
0028 using namespace Acts::detail;
0029
0030 namespace ActsPlugins::EDM4hepUtil {
0031 namespace detail {
0032
0033 ActsSquareMatrix<6> jacobianToEdm4hep(double theta, double qOverP, double Bz) {
0034   // Calculate jacobian from our internal parametrization (d0, z0, phi, theta,
0035   // q/p) to the LCIO / edm4hep (see:
0036   // https://bib-pubdb1.desy.de/record/81214/files/LC-DET-2006-004%5B1%5D.pdf)
0037   // one (d0, z0, phi, tan(lambda), omega). Top left 3x3 matrix in the
0038   // jacobian is 1. Bottom right 2x2 matrix is:
0039   //
0040   // [  d                                 ]
0041   // [------(cot(theta))         0        ]
0042   // [dtheta                              ]
0043   // [                                    ]
0044   // [  d   /  B*q/p   \   d  /  B*q/p   \]
0045   // [------|----------|  ----|----------|]
0046   // [dtheta\sin(theta)/  dq/p\sin(theta)/]
0047   //
0048   // =
0049   //
0050   // [     2                        ]
0051   // [- csc (theta)           0     ]
0052   // [                              ]
0053   // [-B*q/p*cos(theta)       B     ]
0054   // [------------------  ----------]
0055   // [      2             sin(theta)]
0056   // [   sin (theta)                ]
0057
0058   ActsSquareMatrix<6> J;
0059   J.setIdentity();
0060   double sinTheta = std::sin(theta);
0061   J(3, 3) = -1.0 / (sinTheta * sinTheta);
0062   J(4, 4) = Bz / sinTheta;  // dOmega / d(qop)
0063   J(4, 3) = -Bz * qOverP * std::cos(theta) /
0064             (sinTheta * sinTheta);  // dOmega / dTheta
0065   return J;
0066 }
0067
0068 ActsSquareMatrix<6> jacobianFromEdm4hep(double tanLambda, double omega,
0069                                         double Bz) {
0070   // [     d      /                     pi\                                  ]
0071   // [------------|-atan(\tan\lambda) + --|                 0                ]
0072   // [d\tan\lambda\                     2 /                                  ]
0073   // [                                                                       ]
0074   // [     d      /         \Omega        \     d   /         \Omega        \]
0075   // [------------|-----------------------|  -------|-----------------------|]
0076   // [d\tan\lambda|     __________________|  d\Omega|     __________________|]
0077   // [            |    /            2     |         |    /            2     |]
0078   // [            \B*\/  \tan\lambda  + 1 /         \B*\/  \tan\lambda  + 1 /]
0079   //
0080   // =
0081   //
0082   // [         -1                                     ]
0083   // [   ----------------                 0           ]
0084   // [              2                                 ]
0085   // [   \tan\lambda  + 1                             ]
0086   // [                                                ]
0087   // [  -\Omega*\tan\lambda               1           ]
0088   // [-----------------------  -----------------------]
0089   // [                    3/2       __________________]
0090   // [  /           2    \         /            2     ]
0091   // [B*\\tan\lambda  + 1/     B*\/  \tan\lambda  + 1 ]
0092
0093   ActsSquareMatrix<6> J;
0094   J.setIdentity();
0095   J(3, 3) = -1 / (tanLambda * tanLambda + 1);
0096   J(4, 3) = -1 * omega * tanLambda /
0097             (Bz * std::pow(tanLambda * tanLambda + 1, 3. / 2.));
0098   J(4, 4) = 1 / (Bz * std::hypot(tanLambda, 1));
0099
0100   return J;
0101 }
0102
0103 void packCovariance(const ActsSquareMatrix<6>& from, float* to) {
0104   for (int i = 0; i < from.rows(); i++) {
0105     for (int j = 0; j <= i; j++) {
0106       std::size_t k = (i + 1) * i / 2 + j;
0107       to[k] = from(i, j);
0108     }
0109   }
0110 }
0111
0112 void unpackCovariance(const float* from, ActsSquareMatrix<6>& to) {
0113   auto k = [](std::size_t i, std::size_t j) { return (i + 1) * i / 2 + j; };
0114   for (int i = 0; i < to.rows(); i++) {
0115     for (int j = 0; j < to.cols(); j++) {
0116       to(i, j) = from[j <= i ? k(i, j) : k(j, i)];
0117     }
0118   }
0119 }
0120
0121 Parameters convertTrackParametersToEdm4hep(const GeometryContext& gctx,
0122                                            double Bz,
0123                                            const BoundTrackParameters& params) {
0124   Vector3 global = params.referenceSurface().localToGlobal(
0125       gctx, params.parameters().template head<2>(), params.direction());
0126
0127   std::shared_ptr<const Surface> refSurface =
0128       params.referenceSurface().getSharedPtr();
0129   BoundVector targetPars = params.parameters();
0130   std::optional<BoundSquareMatrix> targetCov = params.covariance();
0131
0132   // If the reference surface is a perigee surface, we use that. Otherwise
0133   // we create a new perigee surface at the global position of the track
0134   // parameters.
0135   if (dynamic_cast<const PerigeeSurface*>(refSurface.get()) == nullptr) {
0136     refSurface = Surface::makeShared<PerigeeSurface>(global);
0137
0138     // We need to convert to the target parameters
0139     // Keep the free parameters around we might need them for the covariance
0140     // conversion
0141
0142     auto perigeeParams =
0143         boundToBoundConversion(gctx, params, *refSurface, Vector3{0, 0, Bz})
0144             .value();
0145     targetPars = perigeeParams.parameters();
0146     targetCov = perigeeParams.covariance();
0147   }
0148
0149   Parameters result;
0150   result.surface = refSurface;
0151
0152   // Only run covariance conversion if we have a covariance input
0153   if (targetCov) {
0154     ActsSquareMatrix<6> J = jacobianToEdm4hep(targetPars[eBoundTheta],
0155                                               targetPars[eBoundQOverP], Bz);
0156     result.covariance = J * targetCov.value() * J.transpose();
0157   }
0158
0159   result.values[0] = targetPars[eBoundLoc0];
0160   result.values[1] = targetPars[eBoundLoc1];
0161   result.values[2] = targetPars[eBoundPhi];
0162   result.values[3] = std::tan(std::numbers::pi / 2. - targetPars[eBoundTheta]);
0163   result.values[4] =
0164       targetPars[eBoundQOverP] / std::sin(targetPars[eBoundTheta]) * Bz;
0165   result.values[5] = targetPars[eBoundTime];
0166
0167   result.particleHypothesis = params.particleHypothesis();
0168
0169   return result;
0170 }
0171
0172 BoundTrackParameters convertTrackParametersFromEdm4hep(
0173     double Bz, const Parameters& params) {
0174   BoundVector targetPars;
0175
0176   ActsSquareMatrix<6> J =
0177       jacobianFromEdm4hep(params.values[3], params.values[4], Bz);
0178
0179   std::optional<BoundMatrix> cov;
0180   if (params.covariance.has_value()) {
0181     cov = J * params.covariance.value() * J.transpose();
0182   }
0183
0184   targetPars[eBoundLoc0] = params.values[0];
0185   targetPars[eBoundLoc1] = params.values[1];
0186   targetPars[eBoundPhi] = params.values[2];
0187   targetPars[eBoundTheta] = std::numbers::pi / 2. - std::atan(params.values[3]);
0188   targetPars[eBoundQOverP] =
0189       params.values[4] * std::sin(targetPars[eBoundTheta]) / Bz;
0190   targetPars[eBoundTime] = params.values[5];
0191
0192   return {params.surface, targetPars, cov, params.particleHypothesis};
0193 }
0194
0195 }  // namespace detail
0196
0197 #if EDM4HEP_VERSION_MAJOR >= 1 || \
0198     (EDM4HEP_VERSION_MAJOR == 0 && EDM4HEP_VERSION_MINOR == 99)
0199 edm4hep::MCParticle getParticle(const edm4hep::SimTrackerHit& hit) {
0200   return hit.getParticle();
0201 }
0202
0203 void setParticle(edm4hep::MutableSimTrackerHit& hit,
0204                  const edm4hep::MCParticle& particle) {
0205   hit.setParticle(particle);
0206 }
0207 #else
0208 edm4hep::MCParticle getParticle(const edm4hep::SimTrackerHit& hit) {
0209   return hit.getMCParticle();
0210 }
0211
0212 void setParticle(edm4hep::MutableSimTrackerHit& hit,
0213                  const edm4hep::MCParticle& particle) {
0214   hit.setMCParticle(particle);
0215 }
0216 #endif
0217
0218 }  // namespace ActsPlugins::EDM4hepUtil