refactor edm4hep cov conversion

andiwand · Dec 12, 2023 · 28eeec6 · 28eeec6
1 parent 739adf0
commit 28eeec6
Showing 1 changed file with 17 additions and 11 deletions.
diff --git a/Plugins/EDM4hep/src/EDM4hepUtil.cpp b/Plugins/EDM4hep/src/EDM4hepUtil.cpp
@@ -9,6 +9,7 @@
 #include "Acts/Plugins/EDM4hep/EDM4hepUtil.hpp"
 
 #include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Definitions/TrackParametrization.hpp"
 #include "Acts/Definitions/Units.hpp"
 #include "Acts/EventData/MultiTrajectory.hpp"
 #include "Acts/EventData/MultiTrajectoryHelpers.hpp"
@@ -116,14 +117,8 @@ Parameters convertTrackParametersToEdm4hep(const Acts::GeometryContext& gctx,
 
   std::shared_ptr<const Acts::Surface> refSurface =
       params.referenceSurface().getSharedPtr();
-
   Acts::BoundVector targetPars = params.parameters();
-  std::optional<Acts::FreeVector> freePars;
-
-  auto makeFreePars = [&]() {
-    return Acts::detail::transformBoundToFreeParameters(
-        params.referenceSurface(), gctx, params.parameters());
-  };
+  std::optional<Acts::BoundSquareMatrix> targetCov = params.covariance();
 
   // If the reference surface is a perigee surface, we use that. Otherwise
   // we create a new perigee surface at the global position of the track
@@ -134,20 +129,31 @@ Parameters convertTrackParametersToEdm4hep(const Acts::GeometryContext& gctx,
     // We need to convert to the target parameters
     // Keep the free parameters around we might need them for the covariance
     // conversion
-    freePars = makeFreePars();
-    targetPars = Acts::detail::transformFreeToBoundParameters(freePars.value(),
+    Acts::FreeVector freePars = Acts::detail::transformBoundToFreeParameters(
+        params.referenceSurface(), gctx, params.parameters());
+    targetPars = Acts::detail::transformFreeToBoundParameters(freePars,
                                                               *refSurface, gctx)
                      .value();
+
+    if (targetCov) {
+      // We need to convert the covariance as well
+      Acts::BoundToFreeMatrix boundToFree =
+          params.referenceSurface().boundToFreeJacobian(gctx, targetPars);
+      Acts::FreeToBoundMatrix freeToBound =
+          refSurface->freeToBoundJacobian(gctx, freePars);
+      Acts::BoundSquareMatrix boundToBound = freeToBound * boundToFree;
+      targetCov = boundToBound * targetCov.value() * boundToBound.transpose();
+    }
   }
 
   Parameters result;
   result.surface = refSurface;
 
   // Only run covariance conversion if we have a covariance input
-  if (params.covariance()) {
+  if (targetCov) {
     Acts::ActsSquareMatrix<6> J = jacobianToEdm4hep(
         targetPars[eBoundTheta], targetPars[eBoundQOverP], Bz);
-    result.covariance = J * params.covariance().value() * J.transpose();
+    result.covariance = J * targetCov.value() * J.transpose();
   }
 
   result.values[0] = targetPars[Acts::eBoundLoc0];