Add in-memory DD4hep to Acts Geo Gen3 conversion

k4ActsTracking/CMakeLists.txt

@@ -41,6 +41,8 @@ set(_plugin_sources
   src/components/GeometryIdSelector.cxx
   src/components/Helpers.cxx
   src/components/TrackTruthAlg.cxx
+  src/components/ActsTestPropagator.cpp
+  src/components/DD4hepBlueprintConstruction.cpp
 )
 
 gaudi_add_module(k4ActsTrackingPlugins

k4ActsTracking/examples/plot_steps.py

@@ -0,0 +1,67 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2014-2024 Key4hep-Project.
+#
+# This file is part of Key4hep.
+# See https://key4hep.github.io/key4hep-doc/ for further info.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+
+import ROOT
+import argparse
+import sys
+
+
+def main():
+    parser = argparse.ArgumentParser(description="Plot step_r vs step_z from ROOT file")
+    parser.add_argument("input_file", help="Input ROOT file")
+    parser.add_argument(
+        "-o", "--output", default="step_plot.png", help="Output plot file"
+    )
+    args = parser.parse_args()
+
+    # Enable multi-threading for RDataFrame
+    ROOT.EnableImplicitMT()
+
+    # Open the ROOT file and create RDataFrame
+    try:
+        df = ROOT.RDataFrame("events", args.input_file)
+    except Exception as e:
+        print(f"Error opening file {args.input_file}: {e}")
+        sys.exit(1)
+
+    # Define step_r as sqrt(step_x^2 + step_y^2)
+    df = df.Define("step_r", "sqrt(step_x*step_x + step_y*step_y)")
+
+    # Create the scatter plot
+    hist = df.Graph("step_z", "step_r")
+
+    # Create canvas and draw
+    canvas = ROOT.TCanvas("canvas", "Step R vs Z", 1200, 800)
+    # hist.SetMarkerStyle(20)
+    # hist.SetMarkerSize(0.5)
+    hist.SetTitle("Step R vs Z;step_z;step_r")
+    hist.Draw("AP")
+
+    # Set axis limits
+    hist.GetXaxis().SetRangeUser(-2000, 2000)
+    hist.GetYaxis().SetRangeUser(0, 1500)
+
+    # Save the plot
+    canvas.SaveAs(args.output)
+    print(f"Plot saved as {args.output}")
+
+
+if __name__ == "__main__":
+    main()

k4ActsTracking/examples/test_visualize_acts_geo.py

@@ -0,0 +1,64 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2014-2024 Key4hep-Project.
+#
+# This file is part of Key4hep.
+# See https://key4hep.github.io/key4hep-doc/ for further info.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+
+import pathlib
+
+from Gaudi.Configuration import DEBUG, VERBOSE
+
+from Configurables import ActsGeoSvc, GeoSvc, ActsTestPropagator, EventDataSvc
+from k4FWCore import ApplicationMgr, IOSvc
+from k4FWCore.parseArgs import parser
+
+parser.add_argument("--compactFile", help="Compact file")
+parser.add_argument(
+    "--test-propagation",
+    help="Test propagation through the geometry using an ACTS particle gun",
+    action="store_true",
+    default=False,
+)
+
+args = parser.parse_known_args()[0]
+
+iosvc = IOSvc()
+if args.test_propagation:
+    iosvc.Output = f"{pathlib.Path(args.compactFile).stem}-steps.root"
+
+geoSvc = GeoSvc()
+geoSvc.detectors = [args.compactFile]
+
+actsGeoSvc = ActsGeoSvc("ActsGeoSvc")
+actsGeoSvc.DumpVisualization = True
+actsGeoSvc.ObjVisFileName = f"{pathlib.Path(args.compactFile).stem}-acts-geo.obj"
+actsGeoSvc.OutputLevel = VERBOSE
+
+alg_list = []
+
+if args.test_propagation:
+    propTest = ActsTestPropagator("TestPropagator")
+    propTest.OutputLevel = DEBUG
+    propTest.NumTracks = 100
+    alg_list.append(propTest)
+
+ApplicationMgr(
+    TopAlg=alg_list,
+    ExtSvc=[geoSvc, actsGeoSvc, EventDataSvc()],
+    EvtMax=1,
+    EvtSel="NONE",
+)

k4ActsTracking/include/k4ActsTracking/IActsGeoSvc.h

@@ -21,6 +21,9 @@
 #define IACTSGEOSVC_H
 
 #include <GaudiKernel/IService.h>
+
+#include <cstdint>
+#include <memory>
 #include <unordered_map>
 
 namespace dd4hep {
@@ -32,18 +35,21 @@ namespace dd4hep {
 namespace Acts {
   class TrackingGeometry;
   class Surface;
+  class MagneticFieldProvider;
 }  // namespace Acts
 
 class GAUDI_API IActsGeoSvc : virtual public IService {
 public:
-  using VolumeSurfaceMap = std::unordered_map<uint64_t, const Acts::Surface*>;
+  using CellIDSurfaceMap = std::unordered_map<uint64_t, const Acts::Surface*>;
 
 public:
   DeclareInterfaceID(IActsGeoSvc, 1, 0);
 
-  virtual const Acts::TrackingGeometry& trackingGeometry() const = 0;
+  virtual std::shared_ptr<const Acts::TrackingGeometry>      trackingGeometry() const   = 0;
+  virtual std::shared_ptr<const Acts::MagneticFieldProvider> magneticField() const      = 0;
+  virtual const CellIDSurfaceMap&                            cellIdToSurfaceMap() const = 0;
 
-  virtual ~IActsGeoSvc() {}
+  virtual ~IActsGeoSvc() = default;
 };
 
 #endif  // IACTSGEOSVC_H

k4ActsTracking/src/components/ActsGeoSvc.cpp

@@ -16,90 +16,145 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
+
 #include "ActsGeoSvc.h"
-#include <Acts/Geometry/GeometryContext.hpp>
+#include "DD4hepBlueprintConstruction.h"
 
 #include "k4ActsTracking/ActsGaudiLogger.h"
 
-#include "k4Interface/IGeoSvc.h"
+#include <k4FWCore/GaudiChecks.h>
+
+#include <Acts/Definitions/Algebra.hpp>
+#include <Acts/Definitions/Units.hpp>
+#include <Acts/Geometry/Blueprint.hpp>
+#include <Acts/Geometry/BlueprintNode.hpp>
+#include <Acts/Geometry/BlueprintOptions.hpp>
+#include <Acts/Geometry/ContainerBlueprintNode.hpp>
+#include <Acts/Geometry/CylinderVolumeBounds.hpp>
+#include <Acts/Geometry/Extent.hpp>
+#include <Acts/Geometry/GeometryContext.hpp>
+#include <Acts/Geometry/GeometryIdentifier.hpp>
+#include <Acts/Geometry/TrackingGeometry.hpp>
+#include <Acts/Geometry/VolumeAttachmentStrategy.hpp>
+#include <Acts/MagneticField/ConstantBField.hpp>
+#include <Acts/Surfaces/Surface.hpp>
+#include <Acts/Utilities/AxisDefinitions.hpp>
+#include <Acts/Visualization/ObjVisualization3D.hpp>
+#include <ActsPlugins/DD4hep/BlueprintBuilder.hpp>
+#include <ActsPlugins/DD4hep/DD4hepDetectorElement.hpp>
+
+#include <DD4hep/DD4hepUnits.h>
+#include <DD4hep/DetElement.h>
+#include <DD4hep/Detector.h>
+
+#include <GaudiKernel/StatusCode.h>
 
-#include "Acts/Geometry/TrackingGeometry.hpp"
-#include "Acts/Visualization/GeometryView3D.hpp"
-#include "Acts/Visualization/ObjVisualization3D.hpp"
-#if __has_include("ActsPlugins/DD4hep/ConvertDD4hepDetector.hpp")
-#include "ActsPlugins/DD4hep/ConvertDD4hepDetector.hpp"
-#else
-#include "Acts/Plugins/DD4hep/ConvertDD4hepDetector.hpp"
-namespace ActsPlugins {
-  using Acts::convertDD4hepDetector;
-  using Acts::sortDetElementsByID;
+#include <fmt/format.h>
+#include <fmt/ostream.h>
+#include <fmt/ranges.h>
 
-}  // namespace ActsPlugins
-#endif
+#include <array>
 
-using namespace Gaudi;
+template <> struct fmt::formatter<Acts::GeometryIdentifier> : fmt::ostream_formatter {};
 
 DECLARE_COMPONENT(ActsGeoSvc)
 
-ActsGeoSvc::ActsGeoSvc(const std::string& name, ISvcLocator* svc)
-    : base_class(name, svc), m_trackingGeoCtx(Acts::GeometryContext::dangerouslyDefaultConstruct()) {}
+ActsGeoSvc::ActsGeoSvc(const std::string& name, ISvcLocator* svcLoc) : base_class(name, svcLoc) {
+  m_bluePrintPopulationFuncs = {{"MAIA_v0", MuColl::MAIA_v0::populateBlueprint},
+                                {"MuSIC_v2", MuColl::MAIA_v0::populateBlueprint},
+                                {"ILD_FCCee_v01", FCCee::ILD_FCCee_v01::populateBlueprint},
+                                {"ILD_FCCee_v02", FCCee::ILD_FCCee_v02::populateBlueprint},
+                                {"CLD_o2_v07", FCCee::CLD_o2_v07::populateBlueprint},
+                                {"CLD_o2_v08", FCCee::CLD_o2_v07::populateBlueprint},
+                                {"LUXE_v0", LUXE::LUXE_v0::populateBlueprint}};
+}
 
 StatusCode ActsGeoSvc::initialize() {
-  m_dd4hepGeo = svcLocator()->service<IGeoSvc>(m_geoSvcName)->getDetector();
-  // necessary?
-  // m_dd4hepGeo->addExtension<IActsGeoSvc>(this);
-
-  Acts::BinningType bTypePhi              = Acts::equidistant;
-  Acts::BinningType bTypeR                = Acts::equidistant;
-  Acts::BinningType bTypeZ                = Acts::equidistant;
-  double            layerEnvelopeR        = Acts::UnitConstants::mm;
-  double            layerEnvelopeZ        = Acts::UnitConstants::mm;
-  double            defaultLayerThickness = Acts::UnitConstants::fm;
-
-  auto logger   = makeActsGaudiLogger(this);
-  m_trackingGeo = ActsPlugins::convertDD4hepDetector(
-      m_dd4hepGeo->world(), *logger, bTypePhi, bTypeR, bTypeZ, layerEnvelopeR, layerEnvelopeZ, defaultLayerThickness,
-      ActsPlugins::sortDetElementsByID, m_trackingGeoCtx, m_materialDeco);
-
-  /// Setting geometry debug option
-  if (m_debugGeometry == true) {
-    info() << "Geometry debugging is ON." << endmsg;
-
-    if (createGeoObj().isFailure()) {
-      error() << "Could not create geometry OBJ" << endmsg;
-      return StatusCode::FAILURE;
-    } else {
-      info() << "Geometry OBJ SUCCESSFULLY created" << endmsg;
-    }
+  m_geoSvc = Gaudi::svcLocator()->service<IGeoSvc>("GeoSvc");
+  K4_GAUDI_CHECK(m_geoSvc);
+
+  std::array<double, 3> magneticFieldVector = {0, 0, 0};
+  std::array<double, 3> position            = {0, 0, 0};
+  m_geoSvc->getDetector()->field().magneticField(position.data(), magneticFieldVector.data());
+  debug() << fmt::format("Retrieved magnetic field at position {}: {}", position, magneticFieldVector) << endmsg;
+  m_magneticField = std::make_shared<Acts::ConstantBField>(
+      Acts::Vector3(magneticFieldVector[0] / dd4hep::tesla * Acts::UnitConstants::T,
+                    magneticFieldVector[1] / dd4hep::tesla * Acts::UnitConstants::T,
+                    magneticFieldVector[2] / dd4hep::tesla * Acts::UnitConstants::T));
+
+  auto gaudiLogger = makeActsGaudiLogger(this);
+
+  info() << fmt::format("Acts::cm: {}, dd4hep::cm: {}", Acts::UnitConstants::cm, dd4hep::cm) << endmsg;
+
+  auto gctxt = Acts::GeometryContext::dangerouslyDefaultConstruct();
+
+  const auto* dd4hepDet = m_geoSvc->getDetector();
+  const auto  detName   = dd4hepDet->header().name();
+  info() << fmt::format("Constructing detector with name: {}", dd4hepDet->header().name()) << endmsg;
+
+  BlueprintBuilder builder{
+      {.dd4hepDetector = dd4hepDet, .lengthScale = Acts::UnitConstants::cm / dd4hep::cm, .gctx = gctxt},
+      gaudiLogger->cloneWithSuffix("|BlpBld")};
+
+  using Acts::Experimental::Blueprint;
+  using Acts::Experimental::BlueprintOptions;
+  using namespace Acts::UnitLiterals;
+  using enum Acts::AxisDirection;
+
+  Blueprint::Config cfg;
+  // Padding around subvolumes of the world volume
+  cfg.envelope[AxisZ] = {20_mm, 20_mm};
+  cfg.envelope[AxisR] = {0_mm, 20_mm};
+  Blueprint root{cfg};
+
+  debug() << fmt::format("Getting Blueprint construction function for detector: {}", detName) << endmsg;
+  if (const auto it = m_bluePrintPopulationFuncs.find(detName); it != m_bluePrintPopulationFuncs.end()) {
+    auto bluePrintFunc = it->second;
+    bluePrintFunc(detName, root, builder);
   } else {
-    info() << "Geometry converted without checking if GeoObj can be created" << endmsg;
+    error() << fmt::format("Cannot find a Blueprint construction function for detector: {}", detName) << endmsg;
+    return StatusCode::FAILURE;
   }
 
-  return StatusCode::SUCCESS;
-}
+  BlueprintOptions options;
 
-StatusCode ActsGeoSvc::execute() { return StatusCode::SUCCESS; }
+  debug() << "Constructing tracking geometry" << endmsg;
+  m_trackingGeo = root.construct(options, gctxt, *gaudiLogger->cloneWithSuffix("|Construct"));
 
-StatusCode ActsGeoSvc::finalize() { return StatusCode::SUCCESS; }
-
-/// Create a geometry OBJ file
-StatusCode ActsGeoSvc::createGeoObj() {
-  // Convert DD4Hep geometry to acts
-
-  Acts::ObjVisualization3D m_obj;
-
-  if (!m_trackingGeo) {
-    return StatusCode::FAILURE;
-  }
+  std::size_t nSurfaces = 0;
   m_trackingGeo->visitSurfaces([&](const Acts::Surface* surface) {
-    if (surface == nullptr) {
-      info() << "no surface??? " << endmsg;
-      return;
+    nSurfaces++;
+    const auto& actsDetElem = dynamic_cast<const ActsPlugins::DD4hepDetectorElement&>(*surface->surfacePlacement());
+    const auto& detElem     = actsDetElem.sourceElement();
+    verbose() << fmt::format("Adding Acts surface {} pointing to dd4hep DetElement {}", surface->geometryId(),
+                             detElem.volumeID())
+              << endmsg;
+    const auto& [existing, inserted] = m_cellIDToSurface.emplace(detElem.volumeID(), surface);
+    if (!inserted) {
+      error() << fmt::format(
+                     "The Acts surface {} pointing to dd4hep DetElement with cellID {} is already registered in the "
+                     "map for Acts surface {}",
+                     surface->geometryId(), detElem.volumeID(), existing->second->geometryId())
+              << endmsg;
     }
-    Acts::GeometryView3D::drawSurface(m_obj, *surface, m_trackingGeoCtx);
   });
-  m_obj.write(m_outputFileName.value());
-  info() << m_outputFileName << " SUCCESSFULLY written." << endmsg;
+
+  info() << fmt::format("Visited {} Surfaces and inserted {} pairs of CellID -> Acts::Surface* into the map.",
+                        nSurfaces, m_cellIDToSurface.size())
+         << endmsg;
+  if (nSurfaces != m_cellIDToSurface.size()) {
+    error() << fmt::format("{} Surfaces in the Tracking geometry but only {} distinct CellIDs found.", nSurfaces,
+                           m_cellIDToSurface.size());
+    return StatusCode::FAILURE;
+  }
+
+  if (m_dumpVisualization.value()) {
+    info() << "Creating visualiztion" << endmsg;
+    // Adjust the scale here to make it easier to import in blender
+    Acts::ObjVisualization3D vis{4, 0.001};
+    m_trackingGeo->visualize(vis, gctxt);
+    vis.write(m_objDumpFileName.value());
+  }
 
   return StatusCode::SUCCESS;
 }