EDM4hep/src/EDM4hepUtil.cpp

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