Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:12:19

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 "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   // Calculate jacobian from our internal parametrization (d0, z0, phi, theta,
0032   // q/p) to the LCIO / edm4hep (see:
0033   // https://bib-pubdb1.desy.de/record/81214/files/LC-DET-2006-004%5B1%5D.pdf)
0034   // one (d0, z0, phi, tan(lambda), omega). Top left 3x3 matrix in the
0035   // jacobian is 1. Bottom right 2x2 matrix is:
0036   //
0037   // [  d                                 ]
0038   // [------(cot(theta))         0        ]
0039   // [dtheta                              ]
0040   // [                                    ]
0041   // [  d   /  B*q/p   \   d  /  B*q/p   \]
0042   // [------|----------|  ----|----------|]
0043   // [dtheta\sin(theta)/  dq/p\sin(theta)/]
0044   //
0045   // =
0046   //
0047   // [     2                        ]
0048   // [- csc (theta)           0     ]
0049   // [                              ]
0050   // [-B*q/p*cos(theta)       B     ]
0051   // [------------------  ----------]
0052   // [      2             sin(theta)]
0053   // [   sin (theta)                ]
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;  // dOmega / d(qop)
0060   J(4, 3) = -Bz * qOverP * std::cos(theta) /
0061             (sinTheta * sinTheta);  // dOmega / dTheta
0062   return J;
0063 }
0064 
0065 ActsSquareMatrix<6> jacobianFromEdm4hep(double tanLambda, double omega,
0066                                         double Bz) {
0067   // [     d      /                     pi\                                  ]
0068   // [------------|-atan(\tan\lambda) + --|                 0                ]
0069   // [d\tan\lambda\                     2 /                                  ]
0070   // [                                                                       ]
0071   // [     d      /         \Omega        \     d   /         \Omega        \]
0072   // [------------|-----------------------|  -------|-----------------------|]
0073   // [d\tan\lambda|     __________________|  d\Omega|     __________________|]
0074   // [            |    /            2     |         |    /            2     |]
0075   // [            \B*\/  \tan\lambda  + 1 /         \B*\/  \tan\lambda  + 1 /]
0076   //
0077   // =
0078   //
0079   // [         -1                                     ]
0080   // [   ----------------                 0           ]
0081   // [              2                                 ]
0082   // [   \tan\lambda  + 1                             ]
0083   // [                                                ]
0084   // [  -\Omega*\tan\lambda               1           ]
0085   // [-----------------------  -----------------------]
0086   // [                    3/2       __________________]
0087   // [  /           2    \         /            2     ]
0088   // [B*\\tan\lambda  + 1/     B*\/  \tan\lambda  + 1 ]
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   // If the reference surface is a perigee surface, we use that. Otherwise
0130   // we create a new perigee surface at the global position of the track
0131   // parameters.
0132   if (dynamic_cast<const Acts::PerigeeSurface*>(refSurface.get()) == nullptr) {
0133     refSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(global);
0134 
0135     // We need to convert to the target parameters
0136     // Keep the free parameters around we might need them for the covariance
0137     // conversion
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   // Only run covariance conversion if we have a covariance input
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 }  // namespace detail
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 }  // namespace Acts::EDM4hepUtil