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