Improvements to Acts / Trk convertors and tests

Need to have surface to calculate the global position for ACTS parameters

Print warning if no surface found

Cache position, and improve debug output

Event/DumpEventDataToJSON/CMakeLists.txt

@@ -5,6 +5,7 @@ atlas_subdir( DumpEventDataToJSON )
 
 # External dependencies:
 find_package( nlohmann_json )
+find_package( Acts COMPONENTS Core )
 
 # Component(s) in the package:
 atlas_add_component( DumpEventDataToJSON

Event/DumpEventDataToJSON/src/DumpEventDataToJsonAlg.cxx

@@ -7,8 +7,10 @@
 #include <cmath>
 #include <fstream>
 #include <string>
+#include <algorithm>    // std::reverse
 
 #include "Acts/EventData/TrackContainer.hpp"
+#include "Acts/EventData/TrackParameters.hpp"
 #include "ActsEvent/MultiTrajectory.h"
 #include "Gaudi/Property.h"
 #include "GaudiKernel/Algorithm.h"
@@ -68,6 +70,9 @@ StatusCode DumpEventDataToJsonAlg::initialize() {
   } else {
     m_extrapolator.disable();
   }
+  
+  ATH_CHECK(m_trackingGeometryTool.retrieve());
+
   return StatusCode::SUCCESS;
 }
 
@@ -75,23 +80,56 @@ StatusCode DumpEventDataToJsonAlg::initialize() {
 template <>
 nlohmann::json DumpEventDataToJsonAlg::getData(const Acts::TrackProxy<Acts::ConstVectorTrackContainer, ActsTrk::ConstMultiTrajectory, Acts::detail::ConstRefHolder, true> &track) {
   nlohmann::json data;
-  
+
+  Acts::GeometryContext gctx = m_trackingGeometryTool->getGeometryContext(getContext()).context();
+
   // ACTS units are GeV, whilst ATLAS is MeV. So we need to convert.
   data["dparams"] = {track.loc0(), track.loc1(), track.phi(), track.theta(), track.qOverP() * 0.001};
+  ATH_MSG_VERBOSE(" Track has dparams"<<data["dparams"][0]<<", "<<data["dparams"][1]<<", "<<data["dparams"][2]<<", "<<data["dparams"][3]<<", "<<data["dparams"][4]);
+  ATH_MSG_VERBOSE(track.referenceSurface().toString(gctx));
+  
+  // Add dparams positions to the output
+  const Acts::BoundTrackParameters trackparams(track.referenceSurface().getSharedPtr(),
+                                                         track.parameters());
+  auto trackPosition = trackparams.position(gctx);;
+  data["pos"].push_back(trackPosition.x());
+  data["pos"].push_back(trackPosition.y());
+  data["pos"].push_back(trackPosition.z());
+  
+  unsigned int nTrackStates = track.nTrackStates();
+  ATH_MSG_VERBOSE("Track has " << nTrackStates << " states.");
+  // Unfortunately actsTracks are stored in reverse order, so we need to do some gymnastics
+  // (There is certainly a more elegant way to do this, but since this will all be changed soon I don't think it matters)
+  std::vector<ActsTrk::ConstMultiTrajectory::ConstTrackStateProxy> trackStates;
+  trackStates.reserve(nTrackStates);
+  for (auto trackstate : track.trackStates()) {
+    trackStates.push_back(trackstate);
+  }
+
+  std::reverse(trackStates.begin(), trackStates.end());
 
-  unsigned int count=0;
-  for ( auto tsos : track.trackStates() ) {
+  unsigned int count = 0;
+  for (auto trackstate : trackStates) {
     // Currently only converting smoothed states, but we will extend this later.
-    if (tsos.hasSmoothed()) {
-      data["pos"].push_back(tsos.smoothed().x());
-      data["pos"].push_back(tsos.smoothed().y());
-      data["pos"].push_back(tsos.smoothed().z());
+    if (trackstate.hasSmoothed() && trackstate.hasReferenceSurface()) {
+      const Acts::BoundTrackParameters params(trackstate.referenceSurface().getSharedPtr(),
+                                                         trackstate.smoothed(),
+                                                         trackstate.smoothedCovariance());
+      ATH_MSG_VERBOSE("Track parameters: "<<params.parameters());
+
+      auto pos = params.position(gctx);
+      data["pos"].push_back(pos.x());
+      data["pos"].push_back(pos.y());
+      data["pos"].push_back(pos.z());
+      ATH_MSG_VERBOSE("TrackState "<<count<<" has smoothed state and reference surface. Position is "<<pos.x()<<", "<<pos.y()<<", "<<pos.z());  
+      ATH_MSG_VERBOSE(params.referenceSurface().toString(gctx));  
+      ATH_MSG_VERBOSE("GeometryId "<<params.referenceSurface().geometryId().value());  
+    } else {
+      ATH_MSG_WARNING("TrackState "<<count<<" does not have smoothed state ["<<trackstate.hasSmoothed()<<"] or reference surface ["<<trackstate.hasReferenceSurface()<<"]. Skipping.");
     }
-    count++;
-
     // TODO: Add measurements etc
+    count++;
   }
-  ATH_MSG_VERBOSE("Track has " << count << " states.");
 
   return data;
 }
@@ -145,7 +183,7 @@ StatusCode DumpEventDataToJsonAlg::execute() {
     ATH_MSG_VERBOSE("Trying to load " << vtcHandle.key() << " with " << vtcHandle->size_impl() << " tracks");
 
     auto multiTraj = std::make_unique<ActsTrk::ConstMultiTrajectory>(&(*tsHandle), &(*pHandle), &(*jHandle), &(*mHandle));
-
+    multiTraj->fillSurfaces(m_trackingGeometryTool->trackingGeometry().get());
     Acts::TrackContainer<Acts::ConstVectorTrackContainer,
                          ActsTrk::ConstMultiTrajectory, Acts::detail::ConstRefHolder>
         tc{*vtcHandle, *multiTraj};
@@ -362,8 +400,6 @@ nlohmann::json DumpEventDataToJsonAlg::getData(const Trk::Track &track) {
                        peri->parameters()[Trk::theta],
                        peri->parameters()[Trk::qOverP]};
 
-    data["pos"] = {peri->position().x(), peri->position().y(),
-                   peri->position().z()};
   } else {
     data["pos"] = {};
   }

Event/DumpEventDataToJSON/src/DumpEventDataToJsonAlg.h

@@ -5,8 +5,11 @@
 #ifndef DUMPEVENTDATATOJSONALG_H
 #define DUMPEVENTDATATOJSONALG_H
 
+// Core
 #include "AthenaBaseComps/AthAlgorithm.h"
 #include "GaudiKernel/ToolHandle.h"
+#include "StoreGate/ReadHandleKey.h"
+//EDM
 #include "xAODEventInfo/EventInfo.h"
 #include "xAODTracking/TrackParticleContainer.h"
 #include "xAODJet/JetContainer.h"
@@ -14,8 +17,10 @@
 #include "xAODCaloEvent/CaloClusterContainer.h"
 #include "CaloEvent/CaloCellContainer.h"
 #include "TrkTrack/TrackCollection.h"
-#include "StoreGate/ReadHandleKey.h"
+//Tools
 #include "TrkExInterfaces/IExtrapolationEngine.h"
+#include "ActsGeometryInterfaces/IActsTrackingGeometryTool.h"
+// PRDs
 #include "MuonPrepRawData/CscPrepDataContainer.h"
 #include "MuonPrepRawData/MdtPrepDataContainer.h"
 #include "MuonPrepRawData/RpcPrepDataContainer.h"
@@ -25,8 +30,10 @@
 #include "InDetPrepRawData/PixelClusterContainer.h"
 #include "InDetPrepRawData/SCT_ClusterContainer.h"
 #include "InDetPrepRawData/TRT_DriftCircleContainer.h"
+// ACTS
 #include "ActsEvent/MultiTrajectory.h"
 #include "Acts/EventData/VectorTrackContainer.hpp"
+// Misc
 #include <nlohmann/json.hpp>
 #include <string>
 #include <vector>
@@ -99,6 +106,7 @@ protected:
   Gaudi::Property<bool> m_extrapolateTrackParticless{this, "ExtrapolateTrackParticles", false, "If true, attempt to extrapolate tracks and add additional positions."};
 
   ToolHandle<Trk::IExtrapolationEngine> m_extrapolator{this, "Extrapolator", "Trk::ExtrapolationEngine/AtlasExtrapolation"};
+  ToolHandle<IActsTrackingGeometryTool> m_trackingGeometryTool{this, "TrackingGeometryTool", "ActsTrackingGeometryTool"};
 
   Gaudi::Property<std::string> m_outputJSON_Name{this, "OutputLocation", "EventData.json", "Default filename for "};
 

Tracking/Acts/ActsEventCnv/src/ActsToTrkConverterTool.cxx

@@ -435,7 +435,7 @@ void ActsTrk::ActsToTrkConverterTool::trkTrackCollectionToActsTrackContainer(
     const Acts::GeometryContext & gctx) const {
   ATH_MSG_VERBOSE("Calling trkTrackCollectionToActsTrackContainer with "
                   << trackColl.size() << " tracks.");
-
+  unsigned int trkCount = 0;
   for (auto trk : trackColl) {
     // Do conversions!
     const Trk::TrackStates *trackStates =
@@ -448,7 +448,7 @@ void ActsTrk::ActsToTrkConverterTool::trkTrackCollectionToActsTrackContainer(
     auto track = tc.getTrack(tc.addTrack());
     auto trackStateContainer = tc.trackStateContainer();
 
-    ATH_MSG_VERBOSE("Track has " << trackStates->size()
+    ATH_MSG_VERBOSE("Track "<<trkCount++<<" has " << trackStates->size()
                                  << " track states on surfaces.");
 
     // loop over track states on surfaces, convert and add them to the ACTS
@@ -477,14 +477,12 @@ void ActsTrk::ActsToTrkConverterTool::trkTrackCollectionToActsTrackContainer(
                                    << actsTSOS.index());
       track.tipIndex() = actsTSOS.index();
 
-      ATH_MSG_VERBOSE("TrackProxy has " << track.nTrackStates()
-                                 << " track states on surfaces.");
       if (tsos->trackParameters()) {
         ATH_MSG_VERBOSE("Converting track parameters.");
         // TODO - work out whether we should set predicted, filtered, smoothed
         const Acts::BoundTrackParameters parameters =
             trkTrackParametersToActsParameters(*(tsos->trackParameters()), gctx);
-
+        ATH_MSG_VERBOSE("Track parameters: "<<parameters.parameters());
         // Sanity check on positions
         actsTrackParameterPositionCheck(parameters, *(tsos->trackParameters()), gctx);
 
@@ -503,6 +501,12 @@ void ActsTrk::ActsToTrkConverterTool::trkTrackCollectionToActsTrackContainer(
           actsTSOS.smoothedCovariance() = *parameters.covariance();
           // Not yet implemented in MultiTrajectory.icc
           // actsTSOS.typeFlags() |= Acts::TrackStateFlag::ParameterFlag;
+          if (!(actsTSOS.hasSmoothed() && actsTSOS.hasReferenceSurface())) {
+            ATH_MSG_WARNING("TrackState does not have smoothed state ["<<actsTSOS.hasSmoothed()<<"] or reference surface ["<<actsTSOS.hasReferenceSurface()<<"].");
+          } else {
+            ATH_MSG_VERBOSE("TrackState has smoothed state and reference surface.");  
+          }
+
         }
       }
       if (tsos->measurementOnTrack()) {
@@ -539,14 +543,21 @@ void ActsTrk::ActsToTrkConverterTool::actsTrackParameterPositionCheck(
     const Acts::BoundTrackParameters &parameters,
     const Trk::TrackParameters &trkparameters,
     const Acts::GeometryContext &gctx) const {
-  if (std::fabs(parameters.position(gctx).x() - trkparameters.position().x()) >
+  auto actsPos = parameters.position(gctx);
+  ATH_MSG_VERBOSE("Acts position: "
+                    << actsPos << " vs "
+                    << trkparameters.position());
+  ATH_MSG_VERBOSE(parameters.referenceSurface().toString(gctx));
+  ATH_MSG_VERBOSE("GeometryId "<<parameters.referenceSurface().geometryId().value());  
+
+  if (std::fabs(actsPos.x() - trkparameters.position().x()) >
           0.01 ||
-      std::fabs(parameters.position(gctx).y() - trkparameters.position().y()) >
+      std::fabs(actsPos.y() - trkparameters.position().y()) >
           0.01 ||
-      std::fabs(parameters.position(gctx).z() - trkparameters.position().z()) >
+      std::fabs(actsPos.z() - trkparameters.position().z()) >
           0.01) {
     ATH_MSG_WARNING("Parameter position mismatch: "
-                    << parameters.position(gctx) << " vs "
+                    << actsPos << " vs "
                     << trkparameters.position());
     ATH_MSG_WARNING("Acts surface:");
     ATH_MSG_WARNING(parameters.referenceSurface().toString(gctx));

Tracking/Acts/ActsEventCnv/test/testConversion.py

@@ -8,42 +8,53 @@ import math
 if "__main__" == __name__:
     RunConversion()
 
-    def _dparamsEqual(acts, trk):
-        for acts_track_params, trk_track_params in zip(acts, trk):
-            if (not math.isclose(acts_track_params, trk_track_params)):
+    def _valuesEqual(acts, trk):
+        for acts_values, trk_values in zip(acts, trk):
+            if (not math.isclose(acts_values, trk_values)):
                 return False
         return True
-
+   
     # Now compare outputs
     import json
     success = False
     with open('dump.json') as f:
         data = json.load(f)
         for event in data:
+            found_ni_differences = 0
             print('Processing', event)
             acts = data[event]['TrackContainers']['ConvertedVectorTrackContainer']
             trk = data[event]['Tracks']['CombinedITkTracks']
             if (acts != trk):
-                # We might need to make this comparison more sophisticated
-                # Since we do not necessarily expect binary identical results
-                # ... but we are a ways off of worrying about this right now
+                found_difference = False
+                # Okay, so simple comparison fails... let's try to find where
                 print('WARNING: Acts and Trk tracks differ')
                 print('We have ', len(acts), '[',
                       len(trk), '] Acts [ Trk ] tracks')
                 for i, (acts_track, trk_track) in enumerate(zip(acts, trk)):
-                    if (not _dparamsEqual(acts_track['dparams'], trk_track['dparams'])):
+                    if (not _valuesEqual(acts_track['dparams'], trk_track['dparams'])):
                         print('ERROR: Acts and Trk dparams differ for track', i)
                         print('Acts dparams:', acts_track['dparams'])
                         print('Trk dparams:', trk_track['dparams'])
-                    if (acts_track['pos'] != trk_track['pos']):
+                        found_difference = True
+                    if (not _valuesEqual(acts_track['pos'], trk_track['pos'])):
                         print('ERROR: Acts and Trk pos differ for track', i)
                         print('Acts pos:', acts_track['pos'])
                         print('Trk pos:', trk_track['pos'])
+                        found_difference = True
+                    if not found_difference:
+                        # Simple comparison failed, but no numerically significant difference found (in what we compare, at least)
+                        # (Of course, this does not exclude that there are important differences in quantities we are not checking!)
+                        found_ni_differences += 1
+                        success=True
             else:
                 print('SUCCESS: Acts and Trk tracks agree')
                 success = True
+    if found_ni_differences > 0:
+        print('WARNING: Found', found_ni_differences, 'tracks which have differences.')
     if not success:
         print('ERROR: the output of the conversion is not correct')
-        # Will uncomment this once the conversion is fixed
-        # import sys
-        # sys.exit(1)
+        import sys
+        sys.exit(1)
+    else:
+        print ('SUCCESS: the output of the conversion is correct (at least, to the precision we check)')
+