Don't use radiusOfInnermostHit for EDM4hep tracks (#72)

* Don't use radiusOfInnermostHit for EDM4hep tracks and compute it from hits when going from EDM4hep to LCIO. See key4hep/EDM4hep#326 * Introdue utility function to determine radius of innermost hit * Add documentation about radius calculation --------- Co-authored-by: jmcarcell <[email protected]> Co-authored-by: tmadlener <[email protected]>
key4hep · Jul 3, 2024 · e102e37 · e102e37
1 parent f6343b2
commit e102e37
Show file tree

Hide file tree

Showing 7 changed files with 33 additions and 6 deletions.
diff --git a/k4EDM4hep2LcioConv/include/k4EDM4hep2LcioConv/k4EDM4hep2LcioConv.h b/k4EDM4hep2LcioConv/include/k4EDM4hep2LcioConv/k4EDM4hep2LcioConv.h
@@ -134,6 +134,10 @@ std::optional<int32_t> attachParticleIDMetaData(IMPL::LCEventImpl* lcEvent, cons
 /**
  * Convert EDM4hep Tracks to LCIO. Simultaneously populate the mapping from
  * EDM4hep to LCIO objects for relation resolving in a second step.
+ *
+ * NOTE: Since the edm4hep::Track does not have a radiusOfInnermostHit field,
+ * this quantity is calculated on the fly from the attached TrackState using the
+ * getRadiusOfStateAtFirstHit function with the default 2D version.
  */
 template <typename TrackMapT>
 std::unique_ptr<lcio::LCCollectionVec> convertTracks(const edm4hep::TrackCollection* const edmCollection,
@@ -430,6 +434,12 @@ convEvent(const podio::Frame& edmEvent, const podio::Frame& metadata = podio::Fr
   return convertEvent(edmEvent, metadata);
 }
 
+/**
+ * Get the radius of the TrackState at the first from the given track. This is
+ * used to set the radiusOfInnermostHit in the LCIO track during the conversion
+ */
+std::optional<double> getRadiusOfStateAtFirstHit(const edm4hep::Track& track, bool use3D = false);
+
 } // namespace EDM4hep2LCIOConv
 
 #include "k4EDM4hep2LcioConv/k4EDM4hep2LcioConv.ipp"

diff --git a/k4EDM4hep2LcioConv/include/k4EDM4hep2LcioConv/k4EDM4hep2LcioConv.ipp b/k4EDM4hep2LcioConv/include/k4EDM4hep2LcioConv/k4EDM4hep2LcioConv.ipp
@@ -1,6 +1,7 @@
 #include "k4EDM4hep2LcioConv/MappingUtils.h"
 
 #include <cassert>
+#include <cmath>
 
 #include "TMath.h"
 
@@ -25,7 +26,7 @@ std::unique_ptr<lcio::LCCollectionVec> convertTracks(const edm4hep::TrackCollect
       lcio_tr->setNdf(edm_tr.getNdf());
       lcio_tr->setdEdx(edm_tr.getDEdx());
       lcio_tr->setdEdxError(edm_tr.getDEdxError());
-      lcio_tr->setRadiusOfInnermostHit(edm_tr.getRadiusOfInnermostHit());
+      lcio_tr->setRadiusOfInnermostHit(getRadiusOfStateAtFirstHit(edm_tr).value_or(-1.0));
 
       // Loop over the hit Numbers in the track
       lcio_tr->subdetectorHitNumbers().resize(edm_tr.subdetectorHitNumbers_size());

diff --git a/k4EDM4hep2LcioConv/include/k4EDM4hep2LcioConv/k4Lcio2EDM4hepConv.ipp b/k4EDM4hep2LcioConv/include/k4EDM4hep2LcioConv/k4Lcio2EDM4hepConv.ipp
@@ -348,7 +348,6 @@ std::unique_ptr<edm4hep::TrackCollection> convertTracks(const std::string& name,
     lval.setNdf(rval->getNdf());
     lval.setDEdx(rval->getdEdx());
     lval.setDEdxError(rval->getdEdxError());
-    lval.setRadiusOfInnermostHit(rval->getRadiusOfInnermostHit());
 
     auto subdetectorHitNum = rval->getSubdetectorHitNumbers();
     for (auto hitNum : subdetectorHitNum) {

diff --git a/k4EDM4hep2LcioConv/src/k4EDM4hep2LcioConv.cpp b/k4EDM4hep2LcioConv/src/k4EDM4hep2LcioConv.cpp
@@ -147,4 +147,14 @@ std::unique_ptr<lcio::LCEventImpl> convertEvent(const podio::Frame& edmEvent, co
   return lcioEvent;
 }
 
+std::optional<double> getRadiusOfStateAtFirstHit(const edm4hep::Track& track, bool use3D) {
+  for (const auto& state : track.getTrackStates()) {
+    if (state.location == edm4hep::TrackState::AtFirstHit) {
+      const auto refP = state.referencePoint;
+      return std::sqrt(refP.x * refP.x + refP.y * refP.y + use3D * refP.z * refP.z);
+    }
+  }
+  return std::nullopt;
+}
+
 } // namespace EDM4hep2LCIOConv
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
@@ -10,7 +10,7 @@ target_link_libraries(TestUtils PUBLIC EDM4HEP::edm4hep LCIO::lcio)
 target_include_directories(TestUtils PUBLIC ${LCIO_INCLUDE_DIRS})
 
 add_executable(compare-contents compare_contents.cpp)
-target_link_libraries(compare-contents PRIVATE edmCompare podio::podioRootIO)
+target_link_libraries(compare-contents PRIVATE edmCompare podio::podioRootIO k4EDM4hep2LcioConv)
 target_include_directories(compare-contents PRIVATE
   $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}/src>)
 

diff --git a/tests/src/CompareEDM4hepLCIO.cc b/tests/src/CompareEDM4hepLCIO.cc
@@ -1,8 +1,10 @@
 #include "CompareEDM4hepLCIO.h"
 #include "ComparisonUtils.h"
+#include "k4EDM4hep2LcioConv/k4EDM4hep2LcioConv.h"
 
 #include "IMPL/TrackerHitImpl.h"
 
+#include <cmath>
 #include <cstdint>
 
 #include "TMath.h"
@@ -314,7 +316,13 @@ bool compare(const EVENT::Track* lcioElem, const edm4hep::Track& edm4hepElem, co
   ASSERT_COMPARE_VALS(lcioElem->getdEdx(), dxQuantities[0].value, "dEdx in DxQuantities in Track");
   ASSERT_COMPARE_VALS(lcioElem->getdEdxError(), dxQuantities[0].error, "dEdxError in DxQuantities in Track");
 
-  ASSERT_COMPARE(lcioElem, edm4hepElem, getRadiusOfInnermostHit, "radiusOfInnermostHit in Track");
+  double radius = EDM4hep2LCIOConv::getRadiusOfStateAtFirstHit(edm4hepElem).value_or(-1.0);
+  double radius3D = EDM4hep2LCIOConv::getRadiusOfStateAtFirstHit(edm4hepElem, true).value_or(-1.0);
+  const double radiusLCIO = lcioElem->getRadiusOfInnermostHit();
+  if (std::abs(radius - radiusLCIO) > radiusLCIO / 1e6 && std::abs(radius3D - radiusLCIO) > radiusLCIO / 1e6) {
+    std::cerr << "radiusOfInnermostHit in Track (LCIO: " << radiusLCIO << "), EDM4hep: 2d: " << radius
+              << ", 3d: " << radius3D << ")" << std::endl;
+  }
 
   ASSERT_COMPARE_RELATION(lcioElem, edm4hepElem, getTracks, objectMaps.tracks, "Tracks in Track");
 
@@ -323,7 +331,7 @@ bool compare(const EVENT::Track* lcioElem, const edm4hep::Track& edm4hepElem, co
   ASSERT_COMPARE_VALS(lcioTrackStates.size(), edm4hepTrackStates.size(), "number of TrackStates in Track");
   for (size_t i = 0; i < lcioTrackStates.size(); ++i) {
     if (!compare(lcioTrackStates[i], edm4hepTrackStates[i])) {
-      std::cerr << " " << i << " in Track" << std::endl;
+      std::cerr << "TrackState " << i << " in Track" << std::endl;
       return false;
     }
   }

diff --git a/tests/src/EDM4hep2LCIOUtilities.cc b/tests/src/EDM4hep2LCIOUtilities.cc
@@ -199,7 +199,6 @@ edm4hep::TrackCollection createTracks(const int num_elements, const int subdetec
     elem.setType(2); // TODO specific type
     elem.setChi2(i * 10.f);
     elem.setNdf(i * 12);
-    elem.setRadiusOfInnermostHit(i * 5.f);
 
     elem.setDEdx(i);
     elem.setDEdxError(i / std::sqrt(i + 1));