Plugins/DD4hep/DD4hepCylindricalDetectorGen3Tests.cpp

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 <boost/test/unit_test.hpp>
0010
0011 #include "Acts/Definitions/Units.hpp"
0012 #include "Acts/Detector/Detector.hpp"
0013 #include "Acts/Geometry/Blueprint.hpp"
0014 #include "Acts/Geometry/BlueprintNode.hpp"
0015 #include "Acts/Geometry/ContainerBlueprintNode.hpp"
0016 #include "Acts/Geometry/GeometryContext.hpp"
0017 #include "Acts/Geometry/GeometryIdentifierBlueprintNode.hpp"
0018 #include "Acts/Geometry/LayerBlueprintNode.hpp"
0019 #include "Acts/Geometry/MaterialDesignatorBlueprintNode.hpp"
0020 #include "Acts/Geometry/VolumeAttachmentStrategy.hpp"
0021 #include "Acts/Navigation/SurfaceArrayNavigationPolicy.hpp"
0022 #include "Acts/Navigation/TryAllNavigationPolicy.hpp"
0023 #include "Acts/Surfaces/Surface.hpp"
0024 #include "Acts/Surfaces/SurfaceArray.hpp"
0025 #include "Acts/Surfaces/SurfaceBounds.hpp"
0026 #include "Acts/Utilities/Enumerate.hpp"
0027 #include "Acts/Utilities/Logger.hpp"
0028 #include "Acts/Utilities/TransformRange.hpp"
0029 #include "ActsPlugins/DD4hep/DD4hepConversionHelpers.hpp"
0030 #include "ActsPlugins/DD4hep/DD4hepDetectorElement.hpp"
0031 #include "ActsTests/CommonHelpers/CylindricalTrackingGeometry.hpp"
0032 #include "ActsTests/CommonHelpers/TemporaryDirectory.hpp"
0033
0034 #include <format>
0035 #include <fstream>
0036 #include <string>
0037
0038 #include <DD4hep/DetElement.h>
0039 #include <DD4hep/DetFactoryHelper.h>
0040 #include <DD4hep/Detector.h>
0041 #include <XML/Utilities.h>
0042 #include <XMLFragments.hpp>
0043
0044 #include "DD4hepTestsHelper.hpp"
0045
0046 using namespace Acts;
0047 using namespace ActsPlugins;
0048
0049 namespace ActsTests {
0050
0051 GeometryContext tContext;
0052 CylindricalTrackingGeometry cGeometry = CylindricalTrackingGeometry(tContext);
0053
0054 const char* beampipe_head_xml =
0055     R""""(
0056     <detectors>
0057         <detector id="0" name="BeamPipe" type="BarrelDetector">
0058             <type_flags type="DetType_TRACKER + DetType_BEAMPIPE"/>
0059             <layers>
0060                 <layer id="0"  name="BeamPipeLayer">
0061 )"""";
0062
0063 const char* nec_head_xml =
0064     R""""(
0065         <detector id="1" name="PixelNegativeEndcap" type="EndcapDetector" readout="PixelReadout">
0066             <type_flags type="DetType_TRACKER + DetType_BARREL"/>
0067 )"""";
0068
0069 const char* barrel_head_xml =
0070     R""""(
0071         <detector id="2" name="PixelBarrel" type="BarrelDetector" readout="PixelReadout">
0072             <type_flags type="DetType_TRACKER + DetType_BARREL"/>
0073 )"""";
0074
0075 const char* pec_head_xml =
0076     R""""(
0077         <detector id="3" name="PixelPositiveEndcap" type="EndcapDetector" readout="PixelReadout">
0078             <type_flags type="DetType_TRACKER + DetType_BARREL"/>
0079 )"""";
0080
0081 const char* plugin_xml =
0082     R"""(<plugins>
0083     <plugin name="DD4hep_ParametersPlugin">
0084       <argument value="world"/>
0085       <argument value="acts_world: bool = true"/>
0086       <argument value="acts_world_type: int = 3"/>
0087       <argument value="acts_world_bvalues_n: int = 3"/>
0088       <argument value="acts_world_bvalues_0: double = 0"/>
0089       <argument value="acts_world_bvalues_1: double = 150*mm"/>
0090       <argument value="acts_world_bvalues_2: double = 1000*mm"/>
0091       <argument value="acts_world_binning: str = r"/>
0092       <argument value="acts_world_geo_id: str = incremental"/>
0093       <argument value="acts_world_root_volume_finder: str = indexed"/>
0094     </plugin>
0095     <plugin name="DD4hep_ParametersPlugin">
0096       <argument value="BeamPipe"/>
0097       <argument value="acts_volume: bool = true"/>
0098       <argument value="acts_volume_type: int = 3"/>
0099       <argument value="acts_volume_bvalues_n: int = 3"/>
0100       <argument value="acts_volume_bvalues_0: double = 0"/>
0101       <argument value="acts_volume_bvalues_1: double = 20*mm"/>
0102       <argument value="acts_volume_bvalues_2: double = 1000*mm"/>
0103       <argument value="acts_volume_internals: bool = true"/>
0104       <argument value="acts_volume_internals_type: str = layer"/>
0105     </plugin>
0106     <plugin name="DD4hep_ParametersPlugin">
0107       <argument value="Pixel"/>
0108       <argument value="acts_container: bool = true"/>
0109       <argument value="acts_container_type: int = 3"/>
0110       <argument value="acts_container_bvalues_n: int = 3"/>
0111       <argument value="acts_container_bvalues_0: double = 20*mm"/>
0112       <argument value="acts_container_bvalues_1: double = 150*mm"/>
0113       <argument value="acts_container_bvalues_2: double = 1000*mm"/>
0114       <argument value="acts_container_binning: str = z"/>
0115     </plugin>
0116     <plugin name="DD4hep_ParametersPlugin">
0117       <argument value="PixelNegativeEndcap"/>
0118       <argument value="acts_container: bool = true"/>
0119       <argument value="acts_container_type: int = 3"/>
0120       <argument value="acts_container_bvalues_n: int = 3"/>
0121       <argument value="acts_container_bvalues_0: double = 20*mm"/>
0122       <argument value="acts_container_bvalues_1: double = 150*mm"/>
0123       <argument value="acts_container_bvalues_2: double = 200*mm"/>
0124       <argument value="acts_container_binning: str = z"/>
0125       <argument value="acts_container_z: double = -800*mm"/>
0126     </plugin>
0127     <plugin name="DD4hep_ParametersPlugin">
0128       <argument value="PixelBarrel"/>
0129       <argument value="acts_container: bool = true"/>
0130       <argument value="acts_container_type: int = 3"/>
0131       <argument value="acts_container_bvalues_n: int = 3"/>
0132       <argument value="acts_container_bvalues_0: double = 20*mm"/>
0133       <argument value="acts_container_bvalues_1: double = 150*mm"/>
0134       <argument value="acts_container_bvalues_2: double = 600*mm"/>
0135       <argument value="acts_container_binning: str = r"/>
0136     </plugin>
0137     <plugin name="DD4hep_ParametersPlugin">
0138       <argument value="PixelPositiveEndcap"/>
0139       <argument value="acts_container: bool = true"/>
0140       <argument value="acts_container_type: int = 3"/>
0141       <argument value="acts_container_bvalues_n: int = 3"/>
0142       <argument value="acts_container_bvalues_0: double = 20*mm"/>
0143       <argument value="acts_container_bvalues_1: double = 150*mm"/>
0144       <argument value="acts_container_bvalues_2: double = 200*mm"/>
0145       <argument value="acts_container_binning: str = z"/>
0146       <argument value="acts_container_z: double = 800*mm"/>
0147     </plugin>
0148   </plugins>
0149   )""";
0150
0151 const char* indent_4_xml = "    ";
0152 const char* indent_8_xml = "        ";
0153 const char* indent_12_xml = "            ";
0154
0155 void generateXML(const std::filesystem::path& xmlPath) {
0156   CylindricalTrackingGeometry::DetectorStore dStore;
0157
0158   // Nec surfaces
0159   double necZ = -800.;
0160   auto necR0Surfaces = cGeometry.surfacesRing(dStore, 6.4, 12.4, 18., 0.125, 0.,
0161                                               42., necZ, 2., 22u);
0162
0163   auto necR1Surfaces = cGeometry.surfacesRing(dStore, 12.4, 20.4, 30., 0.125,
0164                                               0., 80., necZ, 2., 22u);
0165
0166   std::vector<std::vector<Surface*>> necSurfaces = {necR0Surfaces,
0167                                                     necR1Surfaces};
0168
0169   // Barrel surfaces
0170   std::vector<std::array<double, 2u>> innerOuter = {
0171       {25., 35.}, {65., 75.}, {110., 120.}};
0172   auto b0Surfaces = cGeometry.surfacesCylinder(dStore, 8.4, 36., 0.15, 0.14,
0173                                                31., 3., 2., {16, 14});
0174
0175   auto b1Surfaces = cGeometry.surfacesCylinder(dStore, 8.4, 36., 0.15, 0.14,
0176                                                71., 3., 2., {32, 14});
0177
0178   auto b2Surfaces = cGeometry.surfacesCylinder(dStore, 8.4, 36., 0.15, 0.14,
0179                                                116., 3., 2., {52, 14});
0180
0181   std::vector<std::vector<Surface*>> barrelSurfaces = {b0Surfaces, b1Surfaces,
0182                                                        b2Surfaces};
0183
0184   // Nec surfaces
0185   double pecZ = 800.;
0186   auto pecR0Surfaces = cGeometry.surfacesRing(dStore, 6.4, 12.4, 18., 0.125, 0.,
0187                                               42., pecZ, 2., 22u);
0188
0189   auto pecR1Surfaces = cGeometry.surfacesRing(dStore, 12.4, 20.4, 30., 0.125,
0190                                               0., 80., pecZ, 2., 22u);
0191
0192   std::vector<std::vector<Surface*>> pecSurfaces = {pecR0Surfaces,
0193                                                     pecR1Surfaces};
0194
0195   // Create an XML from it
0196   std::ofstream cxml;
0197   cxml.open(xmlPath);
0198   cxml << head_xml;
0199   cxml << segmentation_xml;
0200
0201   // Add the beam pipe first
0202   cxml << beampipe_head_xml << '\n';
0203   cxml << indent_12_xml << "<acts_passive_surface>" << '\n';
0204   cxml << indent_12_xml
0205        << "<tubs rmin=\"19*mm\" rmax=\"20*mm\" dz=\"800*mm\" "
0206           "material=\"Air\"/>"
0207        << '\n';
0208   cxml << indent_12_xml << "</acts_passive_surface>" << '\n';
0209   cxml << indent_12_xml << "</layer> " << '\n';
0210   cxml << indent_8_xml << "</layers>" << '\n';
0211   cxml << indent_8_xml << "</detector>" << '\n';
0212
0213   // Declare the assembly
0214   const char* pixel_assembly_xml =
0215       R"""(<detector id="4" name="Pixel" type="DD4hep_SubdetectorAssembly" vis="invisible">
0216             <shape name="PixelVolume" type="Tube" rmin="20*mm" rmax="150*mm" dz="1000*mm" material="Air"/>
0217             <composite name="PixelNegativeEndcap"/>
0218             <composite name="PixelBarrel"/>
0219             <composite name="PixelPositiveEndcap"/>
0220         </detector>)""";
0221
0222   cxml << pixel_assembly_xml << '\n';
0223   cxml << "</detectors>" << '\n';
0224
0225   // The Pixel detector
0226
0227   // Nec - 1 Layer only
0228   cxml << "<detectors>" << '\n';
0229
0230   cxml << nec_head_xml << '\n';
0231   cxml << indent_8_xml << "<layers>" << '\n';
0232   cxml << indent_8_xml << "<acts_container/> " << '\n';
0233   cxml << indent_4_xml << "<layer name=\"NegEndcapLayer_0\" id=\"0\">" << '\n';
0234   cxml << indent_4_xml << "<acts_volume dz=\"20*mm\" cz=\"-800.*mm\"/>" << '\n';
0235   cxml << indent_8_xml << "<modules>" << '\n';
0236   for (const auto& ring : necSurfaces) {
0237     for (const auto& s : ring) {
0238       cxml << indent_12_xml
0239            << DD4hepTestsHelper::surfaceToXML(tContext, *s,
0240                                               Transform3::Identity())
0241            << "\n";
0242     }
0243   }
0244   cxml << indent_8_xml << "</modules>" << '\n';
0245   cxml << indent_8_xml << "</layer> " << '\n';
0246   cxml << indent_8_xml << "</layers>" << '\n';
0247   cxml << indent_8_xml << "</detector>" << '\n';
0248
0249   // Barrel
0250   cxml << barrel_head_xml << '\n';
0251   cxml << indent_8_xml << "<layers>" << '\n';
0252   cxml << indent_8_xml << "<acts_container/> " << '\n';
0253   for (const auto [ib, bs] : enumerate(barrelSurfaces)) {
0254     cxml << indent_4_xml << "<layer name=\"PixelBarrel_" << ib << "\" id=\""
0255          << ib << "\">" << '\n';
0256     cxml << indent_4_xml << "<acts_volume rmin=\"" << innerOuter[ib][0u]
0257          << "*mm\" rmax=\"" << innerOuter[ib][1u] << "*mm\"/>";
0258     cxml << indent_8_xml << "<modules>" << '\n';
0259     for (const auto& s : bs) {
0260       cxml << indent_12_xml
0261            << DD4hepTestsHelper::surfaceToXML(tContext, *s,
0262                                               Transform3::Identity())
0263            << "\n";
0264     }
0265     cxml << indent_8_xml << "</modules>" << '\n';
0266     cxml << indent_8_xml << "</layer>" << '\n';
0267   }
0268   cxml << indent_8_xml << "</layers>" << '\n';
0269   cxml << indent_8_xml << "</detector>" << '\n';
0270
0271   // Pec - 1 layer only
0272   cxml << pec_head_xml << '\n';
0273   cxml << indent_8_xml << "<layers>" << '\n';
0274   cxml << indent_8_xml << "<acts_container/> " << '\n';
0275   cxml << indent_4_xml << "<layer name=\"PosEndcapLayer_0\" id=\"0\">" << '\n';
0276   cxml << indent_4_xml << "<acts_volume dz=\"20*mm\" cz=\"800.*mm\"/>" << '\n';
0277   cxml << indent_8_xml << "<modules>" << '\n';
0278   for (const auto& ring : pecSurfaces) {
0279     for (const auto& s : ring) {
0280       cxml << indent_12_xml
0281            << DD4hepTestsHelper::surfaceToXML(tContext, *s,
0282                                               Transform3::Identity())
0283            << "\n";
0284     }
0285   }
0286   cxml << indent_8_xml << "</modules>" << '\n';
0287   cxml << indent_8_xml << "</layer> " << '\n';
0288   cxml << indent_8_xml << "</layers>" << '\n';
0289   cxml << indent_8_xml << "</detector>" << '\n';
0290   cxml << indent_8_xml << "</detectors>" << '\n';
0291
0292   cxml << plugin_xml << '\n';
0293   cxml << end_xml << '\n';
0294   cxml.close();
0295 }
0296
0297 using namespace dd4hep;
0298 using namespace UnitLiterals;
0299
0300 using enum AxisBoundaryType;
0301 using enum AxisDirection;
0302 using enum CylinderVolumeBounds::Face;
0303 using enum SurfaceArrayNavigationPolicy::LayerType;
0304 using AttachmentStrategy = VolumeAttachmentStrategy;
0305 using ResizeStrategy = VolumeResizeStrategy;
0306
0307 // --------- Helper functions --------------
0308 // Print Element Tree
0309 void print_elements(const DetElement& el, unsigned int level = 0) {
0310   for (const auto& [name, child] : el.children()) {
0311     for (unsigned int i = 0; i < level; ++i) {
0312       std::cout << "\t";
0313     }
0314     std::cout << "-> " << name << std::endl;
0315     print_elements(child, level + 1);
0316   }
0317 }
0318
0319 // Find the first element with a given name
0320 const DetElement* find_element(const DetElement& el, const std::string& el_name,
0321                                unsigned int search_depth = 0) {
0322   static unsigned int level = 0;
0323   static bool abort_search = false;
0324   for (const auto& [name, child] : el.children()) {
0325     if (el_name == name) {
0326       abort_search = true;
0327       return &child;
0328     }
0329     if (++level <= search_depth) {
0330       auto el_child = find_element(child, el_name, search_depth);
0331       if (abort_search) {
0332         return el_child;
0333       }
0334     }
0335   }
0336   return nullptr;
0337 }
0338
0339 // Helper function to convert shared_ptr vector to const ptr vector
0340 std::vector<const Surface*> makeConstPtrVector(
0341     const std::vector<std::shared_ptr<Surface>>& surfs) {
0342   std::vector<const Surface*> constPtrs;
0343   constPtrs.reserve(surfs.size());
0344   for (const auto& surf : surfs) {
0345     constPtrs.push_back(surf.get());
0346   }
0347   return constPtrs;
0348 }
0349
0350 // Helper struct to keep ProtoLayer and its associated surfaces together
0351 struct LayerData {
0352   ProtoLayer protoLayer;
0353   std::vector<std::shared_ptr<Surface>> surfaces;
0354
0355   LayerData(const GeometryContext& gctx,
0356             std::vector<std::shared_ptr<Surface>> surfs)
0357       : protoLayer(gctx, makeConstPtrVector(surfs)),
0358         surfaces(std::move(surfs)) {}
0359 };
0360
0361 // Helper function to merge layers that overlap in z
0362 std::vector<LayerData> mergeLayers(const GeometryContext& gctx,
0363                                    std::vector<LayerData> layers) {
0364   using enum AxisDirection;
0365
0366   std::vector<LayerData> mergedLayers;
0367   if (layers.empty()) {
0368     return mergedLayers;
0369   }
0370
0371   mergedLayers.push_back(std::move(layers.front()));
0372
0373   for (std::size_t i = 1; i < layers.size(); i++) {
0374     auto& current = layers[i];
0375     auto& prev = mergedLayers.back();
0376
0377     // Check if they overlap in z
0378     bool overlap =
0379         (current.protoLayer.min(AxisZ) <= prev.protoLayer.max(AxisZ) &&
0380          current.protoLayer.max(AxisZ) >= prev.protoLayer.min(AxisZ));
0381
0382     if (overlap) {
0383       // Merge surfaces
0384       std::vector<std::shared_ptr<Surface>> mergedSurfaces;
0385       mergedSurfaces.reserve(current.surfaces.size() + prev.surfaces.size());
0386       mergedSurfaces.insert(mergedSurfaces.end(), current.surfaces.begin(),
0387                             current.surfaces.end());
0388       mergedSurfaces.insert(mergedSurfaces.end(), prev.surfaces.begin(),
0389                             prev.surfaces.end());
0390
0391       mergedLayers.pop_back();
0392       mergedLayers.emplace_back(gctx, std::move(mergedSurfaces));
0393       auto& merged = mergedLayers.back();
0394       merged.protoLayer.envelope[AxisR] = current.protoLayer.envelope[AxisR];
0395       merged.protoLayer.envelope[AxisZ] = current.protoLayer.envelope[AxisZ];
0396     } else {
0397       mergedLayers.push_back(std::move(current));
0398     }
0399   }
0400
0401   return mergedLayers;
0402 }
0403
0404 // --------- Helper functions --------------
0405
0406 BOOST_AUTO_TEST_SUITE(DD4hepPlugin)
0407
0408 BOOST_AUTO_TEST_CASE(DD4hepCylidricalDetectorExplicit) {
0409   TemporaryDirectory tempDir{};
0410   auto xmlPath = tempDir.path() / "CylindricalDetector.xml";
0411   generateXML(xmlPath);
0412
0413   auto lcdd = &(dd4hep::Detector::getInstance());
0414   lcdd->fromCompact(xmlPath.string());
0415   lcdd->volumeManager();
0416   lcdd->apply("DD4hepVolumeManager", 0, nullptr);
0417
0418   constexpr std::size_t s_beamPipeVolumeId =
0419       1;  // where do volume IDs come from?
0420   constexpr std::size_t s_pixelVolumeId = 10;
0421
0422   DetElement world = lcdd->world();
0423
0424   // std::cout << "DD4Hep Elements: " << std::endl;
0425   // print_elements(world);
0426
0427   auto worldSolidDim = lcdd->worldVolume().solid().dimensions();  // better way?
0428
0429   Experimental::Blueprint::Config cfg;
0430   // Is the following correct?
0431   cfg.envelope[AxisX] = {worldSolidDim[0], worldSolidDim[0]};
0432   cfg.envelope[AxisY] = {worldSolidDim[1], worldSolidDim[1]};
0433   cfg.envelope[AxisZ] = {worldSolidDim[2], worldSolidDim[2]};
0434   // cfg.envelope[AxisR] = {0_mm, worldSolidDim[1]};
0435
0436   auto blueprint = std::make_unique<Experimental::Blueprint>(cfg);
0437   auto& cylinder = blueprint->addCylinderContainer("Detector", AxisR);
0438
0439   // ------- Add Beam Pipe to Blueprint -------
0440   auto level1_elements = world.children();
0441   DetElement bpipe_top = level1_elements["BeamPipe"];
0442   auto bpipe_top_position = bpipe_top.placement().position();
0443
0444   Transform3 beamPipeTransform;
0445   beamPipeTransform.setIdentity();
0446   beamPipeTransform = Translation3(
0447       bpipe_top_position.x(), bpipe_top_position.y(), bpipe_top_position.z());
0448
0449   auto beamPipeDim = bpipe_top.solid().dimensions();
0450   // Do I use the values retrieved from DD4Hep correctly?
0451   double beamPipeRMax = beamPipeDim[1] * 1_cm;
0452   double beamPipeHalfZ = beamPipeDim[2] * 1_cm;
0453
0454   std::vector<std::shared_ptr<DD4hepDetectorElement>> detectorElements;
0455
0456   cylinder.withGeometryIdentifier([&](auto& geoId) {
0457     geoId.setAllVolumeIdsTo(s_beamPipeVolumeId);
0458     geoId.addMaterial("BeamPipe_Material", [&](auto& mat) {
0459       mat.configureFace(
0460           OuterCylinder,
0461           {AxisRPhi, Bound,
0462            20},  // Where do these 20-s come from? (here and on the next line)
0463           {AxisZ, Bound, 20});
0464       mat.addStaticVolume(beamPipeTransform,
0465                           std::make_shared<CylinderVolumeBounds>(
0466                               0, beamPipeRMax, beamPipeHalfZ),
0467                           "BeamPipe");
0468     });
0469   });
0470   // ------- Add Beam Pipe to Blueprint -------
0471
0472   // ------- Add Pixel to Blueprint -------
0473   auto pixelElement = find_element(world, "Pixel");
0474   auto pixelBarrelElement = find_element(*pixelElement, "PixelBarrel");
0475
0476   cylinder.addMaterial("Pixel_Material", [&](auto& mat) {
0477     mat.configureFace(OuterCylinder, {AxisRPhi, Bound, 20}, {AxisZ, Bound, 20});
0478     auto& pixelContainer = mat.addCylinderContainer("Pixel", AxisZ);
0479
0480     // Add barrel container
0481     auto& barrelGeoId = pixelContainer.withGeometryIdentifier();
0482     barrelGeoId.setAllVolumeIdsTo(s_pixelVolumeId)
0483         .incrementLayerIds(1)
0484         .sortBy([](auto& a, auto& b) {
0485           auto& boundsA =
0486               dynamic_cast<const CylinderVolumeBounds&>(a.volumeBounds());
0487           auto& boundsB =
0488               dynamic_cast<const CylinderVolumeBounds&>(b.volumeBounds());
0489           using enum CylinderVolumeBounds::BoundValues;
0490           double aMidR = (boundsA.get(eMinR) + boundsA.get(eMaxR)) / 2.0;
0491           double bMidR = (boundsB.get(eMinR) + boundsB.get(eMaxR)) / 2.0;
0492           return aMidR < bMidR;
0493         });
0494
0495     auto& barrel = barrelGeoId.addCylinderContainer("Pixel_Barrel", AxisR);
0496     barrel.setAttachmentStrategy(AttachmentStrategy::Gap);
0497     barrel.setResizeStrategy(ResizeStrategy::Gap);
0498
0499     std::map<int, std::vector<std::shared_ptr<Surface>>> layers{};
0500     int layerId{0};
0501     for (const auto& [nameLayer, layer] : pixelBarrelElement->children()) {
0502       for (const auto& [nameModule, module] : layer.children()) {
0503         std::string detAxis =
0504             getParamOr<std::string>("axis_definitions", module, "XYZ");
0505         auto dd4hepDetEl = std::make_shared<DD4hepDetectorElement>(
0506             module, detAxis, 1_cm, false, nullptr);
0507         detectorElements.push_back(dd4hepDetEl);
0508         layers[layerId].push_back(dd4hepDetEl->surface().getSharedPtr());
0509       }
0510       layerId++;
0511     }
0512
0513     for (const auto& [ilayer, surfaces] : layers) {
0514       Experimental::BlueprintNode* lparent = nullptr;
0515
0516       // Outermost layer can't have material, because it will get merged
0517       // with the outer cylinder shell of the endcap cylinders
0518       //
0519       //                          Material here is fine   |
0520       //                                                  |
0521       //                                                  |
0522       //                Will get merged. Therefore none   |
0523       //                   of them can have material      |
0524       //                                                  |
0525       //                               |                  |
0526       //                               |                  |
0527       //                               |                  |
0528       //       +-----------------------+-+----------------+---------+
0529       //       |                         |                |         |
0530       //       |                         |                |         |
0531       //       v                         v                |         v
0532       // +----------+-------------------------------------+---+----------+
0533       // |          |                Barrel L2            |   |          |
0534       // |          +-------------------------------------v---+          |
0535       // |          |                   Gap                   |          |
0536       // |          +-----------------------------------------+          |
0537       // |  Neg EC  |                Barrel L1                |  Pos EC  |
0538       // |          +-----------------------------------------+          |
0539       // |          |                   Gap                   |          |
0540       // |          +-----------------------------------------+          |
0541       // |          |                Barrel L0                |          |
0542       // +----------+-----------------------------------------+----------+
0543       if (ilayer < static_cast<int>(layers.size() - 1)) {
0544         auto& lmat = barrel.addMaterial(
0545             std::format("Pixel_Barrel_L{}_Material", ilayer));
0546         lmat.configureFace(OuterCylinder, {AxisRPhi, Bound, 40},
0547                            {AxisZ, Bound, 20});
0548         lparent = &lmat;
0549       } else {
0550         lparent = &barrel;
0551       }
0552
0553       // Add layer with surfaces
0554       auto& layer = lparent->addLayer(std::format("Pixel_Barrel_L{}", ilayer));
0555
0556       // Set navigation policy for efficient surface lookup
0557       layer.setNavigationPolicyFactory(
0558           NavigationPolicyFactory{}
0559               .add<SurfaceArrayNavigationPolicy>(
0560                   SurfaceArrayNavigationPolicy::Config{.layerType = Cylinder,
0561                                                        .bins = {30, 10}})
0562               .add<TryAllNavigationPolicy>(
0563                   TryAllNavigationPolicy::Config{.sensitives = false})
0564               .asUniquePtr());
0565
0566       layer.setSurfaces(surfaces);
0567       layer.setEnvelope(ExtentEnvelope{{
0568           .z = {5_mm, 5_mm},  // ???
0569           .r = {2_mm, 2_mm},  // ???
0570       }});
0571     }
0572
0573     // Add endcap containers
0574     for (int ecid : {-1, 1}) {
0575       const std::string_view s = ecid == 1 ? "p" : "n";
0576       auto& ecGeoId = pixelContainer.withGeometryIdentifier();
0577       ecGeoId.setAllVolumeIdsTo(s_pixelVolumeId + ecid).incrementLayerIds(1);
0578       auto& ec =
0579           ecGeoId.addCylinderContainer(std::format("Pixel_{}EC", s), AxisZ);
0580       ec.setAttachmentStrategy(AttachmentStrategy::Gap);
0581       ec.setResizeStrategy(ResizeStrategy::Expand);
0582
0583       std::map<int, std::vector<std::shared_ptr<Surface>>> initialLayers{};
0584       const DetElement* pixelEndcapElement =
0585           ecid == 1 ? find_element(*pixelElement, "PixelPositiveEndcap")
0586                     : find_element(*pixelElement, "PixelNegativeEndcap");
0587       layerId = 0;
0588       for (const auto& [nameLayer, layer] : pixelEndcapElement->children()) {
0589         for (const auto& [nameModule, module] : layer.children()) {
0590           std::string detAxis =
0591               getParamOr<std::string>("axis_definitions", module, "XYZ");
0592           auto dd4hepDetEl = std::make_shared<DD4hepDetectorElement>(
0593               module, detAxis, 1_cm, false, nullptr);
0594           detectorElements.push_back(dd4hepDetEl);
0595           initialLayers[layerId].push_back(
0596               dd4hepDetEl->surface().getSharedPtr());
0597         }
0598         layerId++;
0599       }
0600
0601       // Create proto layers from surfaces
0602       std::vector<LayerData> protoLayers;
0603       protoLayers.reserve(initialLayers.size());
0604       for (const auto& [key, surfaces] : initialLayers) {
0605         auto& layer = protoLayers.emplace_back(GeometryContext(), surfaces);
0606         layer.protoLayer.envelope[AxisR] = {2_mm, 2_mm};
0607         layer.protoLayer.envelope[AxisZ] = {1_mm, 1_mm};
0608       }
0609       // Sort by z position
0610       std::ranges::sort(protoLayers,
0611                         [](const LayerData& a, const LayerData& b) {
0612                           return std::abs(a.protoLayer.medium(AxisZ)) <
0613                                  std::abs(b.protoLayer.medium(AxisZ));
0614                         });
0615
0616       std::vector<LayerData> mergedLayers =
0617           mergeLayers(GeometryContext(), protoLayers);
0618
0619       // Create layers from proto layers
0620       for (const auto& [key, pl] : enumerate(mergedLayers)) {
0621         pl.protoLayer.medium(AxisZ);
0622         auto layerName = std::format("Pixel_{}EC_L{}", s, key);
0623         auto addLayer = [&layerName, &pl](auto& parent) {
0624           // Add layer with surfaces
0625           auto& layer = parent.addLayer(layerName);
0626
0627           layer.setNavigationPolicyFactory(
0628               NavigationPolicyFactory{}
0629                   .add<SurfaceArrayNavigationPolicy>(
0630                       SurfaceArrayNavigationPolicy::Config{.layerType = Disc,
0631                                                            .bins = {30, 30}})
0632                   .add<TryAllNavigationPolicy>(
0633                       TryAllNavigationPolicy::Config{.sensitives = false})
0634                   .asUniquePtr());
0635
0636           layer.setSurfaces(pl.surfaces);
0637           layer.setEnvelope(ExtentEnvelope{{
0638               .z = {1_mm, 1_mm},
0639               .r = {2_mm, 2_mm},
0640           }});
0641         };
0642
0643         if (key < mergedLayers.size() - 1) {
0644           ec.addMaterial(layerName + "_Material", [&](auto& lmat) {
0645             lmat.configureFace(ecid < 0 ? NegativeDisc : PositiveDisc,
0646                                {AxisR, Bound, 40}, {AxisPhi, Bound, 40});
0647             addLayer(lmat);
0648           });
0649         } else {
0650           addLayer(ec);
0651         }
0652       }
0653     }
0654   });
0655   // ------- Add Pixel to Blueprint -------
0656
0657   std::cout << tempDir.path() << std::endl;
0658   std::ofstream dotOut{tempDir.path() / "CylindricalDetector.dot"};
0659   blueprint->graphviz(dotOut);
0660
0661   // Final step
0662   GeometryContext gctx;
0663   auto logger = getDefaultLogger("Geo", Logging::VERBOSE);
0664   auto trackingGeometry = blueprint->construct({}, gctx, *logger);
0665
0666   BOOST_REQUIRE_NE(&world, nullptr);
0667   BOOST_REQUIRE_NE(trackingGeometry.get(), nullptr);
0668
0669   // Kill that instance before going into the next test
0670   lcdd->destroyInstance();
0671 }
0672
0673 BOOST_AUTO_TEST_SUITE_END()
0674
0675 }  // namespace ActsTests