scripts/view15/.DAWN_1.history

+0.0
+72.5
+189.5
+0
+0
+0
+1500
+50
+3
+0.001
+0
+1
+1
+1
+0.5
+0.5
+0.5
+25.5
+71
+0.001
+0.001
+0.001
+1
+70
+0.001
+1
+0
+1
+evince
+0
+0

scripts/view15/generate_eps

@@ -23,7 +23,7 @@
 
 
 function print_the_help {
-  echo "USAGE: $0 -i <PRIM_FILE>  "
+  echo "USAGE: $0 -i <PRIM_FILE> <slices ...> "
   echo "  OPTIONS: "
   echo "            -t,--tag           filename tag (default: view1)"
   exit 
@@ -53,12 +53,16 @@ do
       shift # past argument
       shift # past value
       ;;
-    *)    # unknown option
-      #POSITIONAL+=("$1") # save it in an array for later
+    -[a-zA-Z]*) # unknown option
+      POSITIONAL+=("$1") # save it in an array for later
       echo "unknown option $1"
       print_the_help
       shift # past argument
       ;;
+    *)     # positional options
+      POSITIONAL+=("$1") # save it in an array for later
+      shift # past argument
+      ;;
   esac
 done
 set -- "${POSITIONAL[@]}" # restore positional parameters
@@ -104,7 +108,7 @@ make_slice(){
   rm "${FILE_TAG}.prim"
 }
 
-for zzz in $(seq 150 50 2000) ;
+for zzz in $@ ;
 do
   make_slice ${zzz} 
 done

scripts/view2/.DAWN_1.history

@@ -6,7 +6,7 @@
 0
 0
 8
-3
+1
 0.001
 0
 1

scripts/view20/.DAWN_1.history

+1.34392e+07
+90
+180
+0
+0
+0
+0
+1
+1
+0.001
+0
+1
+1
+1
+0.5
+0.5
+0.5
+19
+71
+0.001
+0.001
+0.001
+3
+71
+0.001
+0
+0
+1
+evince
+0
+0

scripts/view20/generate_eps

+#!/bin/bash
+
+#export DAWN_PS_PREVIEWER="derp"
+
+function print_the_help {
+  echo "USAGE: $0 <PRIM_FILE>  "
+  echo "  OPTIONS: "
+  echo "            -t,--tag           filename tag (default: view1)"
+  exit 
+}
+
+FILE_TAG="view20"
+INPUT_FILE="../../g4_0000.prim"
+
+
+POSITIONAL=()
+while [[ $# -gt 0 ]]
+do
+  key="$1"
+
+  case $key in
+    -h|--help)
+      shift # past argument
+      print_the_help
+      ;;
+    -t|--tag)
+      FILE_TAG="$2"
+      shift # past argument
+      shift # past value
+      ;;
+    -i|--input)
+      INPUT_FILE="$2"
+      shift # past argument
+      shift # past value
+      ;;
+    *)    # unknown option
+      #POSITIONAL+=("$1") # save it in an array for later
+      echo "unknown option $1"
+      print_the_help
+      shift # past argument
+      ;;
+  esac
+done
+set -- "${POSITIONAL[@]}" # restore positional parameters
+
+
+# Top view with a thin slice down the middle
+dawncut 0 1 0 1 ${INPUT_FILE} ${FILE_TAG}_top_temp0.prim   
+dawncut 0 -1 0 1 ${FILE_TAG}_top_temp0.prim  ${FILE_TAG}_top.prim
+../../bin/dawn_tweak --mag 2 --draw 1 --theta 90 --phi 90
+dawn -d ${FILE_TAG}_top.prim 
+ps2pdf ${FILE_TAG}_top.eps ${FILE_TAG}_top_full.pdf
+gs -o ${FILE_TAG}_top.pdf -sDEVICE=pdfwrite \
+  -c "[/CropBox [51 250 550 590] /PAGES pdfmark" \
+  -f ${FILE_TAG}_top_full.pdf
+pdftoppm ${FILE_TAG}_top.pdf ${FILE_TAG}_top -png -singlefile -cropbox -thinlinemode solid -aaVector yes
+
+# Side view (lines)
+dawncut 1 0 0 1 ${INPUT_FILE} ${FILE_TAG}_temp0.prim 
+dawncut -1 0 0 1 ${FILE_TAG}_temp0.prim  ${FILE_TAG}.prim
+../../bin/dawn_tweak --mag 1 --draw 1 --theta 180 --phi 90 
+dawn -d ${FILE_TAG}.prim
+ps2pdf ${FILE_TAG}.eps ${FILE_TAG}_full.pdf
+gs -o ${FILE_TAG}.pdf -sDEVICE=pdfwrite \
+  -c "[/CropBox [51 250 550 590] /PAGES pdfmark" \
+  -f ${FILE_TAG}_full.pdf
+pdftoppm ${FILE_TAG}.pdf ${FILE_TAG} -png -singlefile -cropbox -thinlinemode solid -aaVector yes
+
+#npdet_info print EcalEndcapN_z0 --value-only ../../athena.xml
+#180.5 cm
+zcut=$(npdet_info print EcalEndcapN_z0 --value-only ${DETECTOR_PATH}/athena.xml )
+NMOD1=$(npdet_info print EcalEndcapN_NModules_Sector1 --value-only ${DETECTOR_PATH}/calorimeters.xml )
+NMOD2=$(npdet_info print EcalEndcapN_NModules_Sector2 --value-only ${DETECTOR_PATH}/calorimeters.xml )
+
+echo "NMOD1 = ${NMOD1}"
+echo "NMOD2 = ${NMOD2}"
+echo "zcut  = ${zcut}"
+
+# Top view with a thin slice down the middle
+dawncut 0 0  1 -1800 ${INPUT_FILE} ${FILE_TAG}_endcapN_temp0.prim   
+dawncut 0 0 -1 2200 ${FILE_TAG}_endcapN_temp0.prim  ${FILE_TAG}_endcapN.prim
+../../bin/dawn_tweak --mag 5 --draw 3 --theta 180 --phi 0
+dawn -d ${FILE_TAG}_endcapN.prim 
+ps2pdf ${FILE_TAG}_endcapN.eps ${FILE_TAG}_endcapN_full.pdf
+gs -o ${FILE_TAG}_endcapN.pdf -sDEVICE=pdfwrite \
+  -c "[/CropBox [50 170 550 670] /PAGES pdfmark" \
+  -f ${FILE_TAG}_endcapN_full.pdf
+pdftoppm ${FILE_TAG}_endcapN.pdf ${FILE_TAG}_endcapN -png -singlefile -cropbox -thinlinemode solid \
+         -aaVector yes -r 1200
+convert -pointsize 180 -fill black -draw "text 200,200 \"$NMOD1 Crystals\"" \
+        ${FILE_TAG}_endcapN.png ${FILE_TAG}_endcapN.png
+convert -pointsize 180 -fill black -draw "text 200,400 \"$NMOD2 Glasses\"" \
+        ${FILE_TAG}_endcapN.png ${FILE_TAG}_endcapN.png
+

scripts/view5/.DAWN_1.history

+1.34392e+07
+0
+0
+1
+0
+0
+491.1
+1.2
+5
+0.001
+0
+1
+1
+1
+0.5
+0.5
+0.5
+19
+71
+0.01
+0.01
+0.01
+3
+70
+0.01
+1
+1
+1
+evince
+0
+0

scripts/view6/.DAWN_1.history

+1.34392e+07
+0
+180
+0
+0
+0
+0
+8.6
+1
+0.001
+0
+1
+1
+1
+0.5
+0.5
+0.5
+19
+71
+0.01
+0.01
+0.01
+3
+70
+0.01
+1
+1
+1
+evince
+0
+0

scripts/view6/generate_eps

@@ -59,6 +59,7 @@ INPUT_FILE=${FILE_TAG}_input.prim
 # units are mm
 dawncut 0 0 -1 1 ${INPUT_FILE} ${FILE_TAG}a_temp0.prim 
 dawncut 0 0 1 1  ${FILE_TAG}a_temp0.prim  ${FILE_TAG}a.prim
+../../bin/dawn_tweak --mag 14 
 dawn -d ${FILE_TAG}a.prim 
 ps2pdf ${FILE_TAG}a.eps ${FILE_TAG}a_full.pdf
 gs -o ${FILE_TAG}a.pdf -sDEVICE=pdfwrite \
@@ -69,6 +70,7 @@ pdftoppm ${FILE_TAG}a.pdf ${FILE_TAG}a -png -singlefile -cropbox
 #SiTracker Endcap layer 5 zstart = 860mm ( 90 mm thick )
 dawncut 0 0 1   945   ${INPUT_FILE} ${FILE_TAG}b_temp0.prim 
 dawncut 0 0 -1 -865 ${FILE_TAG}b_temp0.prim  ${FILE_TAG}b.prim
+../../bin/dawn_tweak --mag 14 
 dawn -d ${FILE_TAG}b.prim 
 ps2pdf ${FILE_TAG}b.eps ${FILE_TAG}b_full.pdf
 gs -o ${FILE_TAG}b.pdf -sDEVICE=pdfwrite \
@@ -79,6 +81,7 @@ pdftoppm ${FILE_TAG}b.pdf ${FILE_TAG}b -png -singlefile -cropbox
 #SiTracker Endcap layer 4 zstart = 695mm ( 90 mm thick )
 dawncut 0 0 1 780   ${INPUT_FILE} ${FILE_TAG}c_temp0.prim 
 dawncut 0 0 -1 -700 ${FILE_TAG}c_temp0.prim  ${FILE_TAG}c.prim
+../../bin/dawn_tweak --mag 14 
 dawn -d ${FILE_TAG}c.prim 
 ps2pdf ${FILE_TAG}c.eps ${FILE_TAG}c_full.pdf
 gs -o ${FILE_TAG}c.pdf -sDEVICE=pdfwrite \
@@ -101,6 +104,7 @@ pdftoppm ${FILE_TAG}d.pdf ${FILE_TAG}d -png -singlefile -cropbox
 # slice at z = -2m
 dawncut 0 0 1 430 ${INPUT_FILE} ${FILE_TAG}e_temp0.prim 
 dawncut 0 0 -1 -370 ${FILE_TAG}e_temp0.prim  ${FILE_TAG}e.prim
+../../bin/dawn_tweak --mag 14 
 dawn -d ${FILE_TAG}e.prim 
 ps2pdf ${FILE_TAG}e.eps ${FILE_TAG}e_full.pdf
 gs -o ${FILE_TAG}e.pdf -sDEVICE=pdfwrite \
@@ -111,6 +115,7 @@ pdftoppm ${FILE_TAG}e.pdf ${FILE_TAG}e -png -singlefile -cropbox
 #SiTracker Endcap layer 1 zstart = 200mm ( 30 mm thick )
 dawncut 0 0 1 225 ${INPUT_FILE} ${FILE_TAG}f_temp0.prim 
 dawncut 0 0 -1 205 ${FILE_TAG}f_temp0.prim  ${FILE_TAG}f.prim
+../../bin/dawn_tweak --mag 14 
 dawn -d ${FILE_TAG}f.prim 
 ps2pdf ${FILE_TAG}f.eps ${FILE_TAG}f_full.pdf
 gs -o ${FILE_TAG}f.pdf -sDEVICE=pdfwrite \

src/BarrelBarDetectorWithSideFrame_geo.cpp

+/** \addtogroup PID
+ * \brief Type: **Barrel bar detector (think DIRC) with longitudinal frame**.
+ * \author S. Joosten
+ *
+ * \ingroup PID
+ *
+ *
+ * \code
+ * \endcode
+ *
+ * @{
+ */
+#include "DD4hep/DetFactoryHelper.h"
+#include "DD4hep/Printout.h"
+#include "DD4hep/Shapes.h"
+#include "DDRec/DetectorData.h"
+#include "DDRec/Surface.h"
+#include "XML/Layering.h"
+#include <array>
+
+using namespace std;
+using namespace dd4hep;
+
+/** Barrel Bar detector with optional framek
+ *
+ * - Optional "frame" tag within the module element (frame eats part of the module width
+ *   and is longitudinal at both sides)
+ * - Detector is setup as a "tracker" so we can use the hits to perform fast MC
+ *   for e.g. the DIRC
+ *
+ */
+static Ref_t create_BarrelBarDetectorWithSideFrame(Detector& description, xml_h e, SensitiveDetector sens)
+{
+  typedef vector<PlacedVolume>            Placements;
+  xml_det_t                               x_det    = e;
+  Material                                air      = description.air();
+  int                                     det_id   = x_det.id();
+  string                                  det_name = x_det.nameStr();
+  DetElement                              sdet(det_name, det_id);
+  map<string, Volume>                     volumes;
+  map<string, Placements>                 sensitives;
+  map<string, std::vector<rec::VolPlane>> volplane_surfaces;
+  PlacedVolume                            pv;
+  dd4hep::xml::Dimension                  dimensions(x_det.dimensions());
+  xml_dim_t                               dirc_pos = x_det.position();
+
+  map<string, std::array<double, 2>> module_thicknesses;
+
+  Tube   topVolumeShape(dimensions.rmin(), dimensions.rmax(), dimensions.length() * 0.5);
+  Volume assembly(det_name, topVolumeShape, air);
+
+  sens.setType("tracker");
+
+  // loop over the modules
+  for (xml_coll_t mi(x_det, _U(module)); mi; ++mi) {
+    xml_comp_t x_mod = mi;
+    string     m_nam = x_mod.nameStr();
+
+    if (volumes.find(m_nam) != volumes.end()) {
+      printout(ERROR, "BarrelBarDetectorWithSideFrame",
+               string((string("Module with named ") + m_nam + string(" already exists."))).c_str());
+      throw runtime_error("Logics error in building modules.");
+    }
+
+    int    ncomponents     = 0;
+    int    sensor_number   = 1;
+    double total_thickness = 0;
+
+    // Compute module total thickness from components
+    xml_coll_t ci(x_mod, _U(module_component));
+    for (ci.reset(), total_thickness = 0.0; ci; ++ci) {
+      total_thickness += xml_comp_t(ci).thickness();
+    }
+    // the module assembly volume
+    Assembly m_vol(m_nam);
+    volumes[m_nam] = m_vol;
+    m_vol.setVisAttributes(description, x_mod.visStr());
+
+    // Optional module frame.
+    // frame is 2 longitudinal bars at the side of the module
+    //
+    // |___| <-- example module cross section (x-y plane), frame is flush with the
+    //           bottom of the module and protrudes on the top if needed
+    //
+
+    // Get frame width, as it impacts the main module for being built. We
+    // construct the actual frame structure later (once we know the module width)
+    double frame_width = 0;
+    if (x_mod.hasChild("frame")) {
+      xml_comp_t m_frame = x_mod.child(_U(frame));
+      frame_width        = m_frame.width();
+    }
+
+    double thickness_so_far    = 0.0;
+    double thickness_sum       = -total_thickness / 2.0;
+    double max_component_width = 0;
+    for (xml_coll_t ci(x_mod, _U(module_component)); ci; ++ci, ++ncomponents) {
+      xml_comp_t x_comp = ci;
+      xml_comp_t x_pos  = x_comp.position(false);
+      string     c_nam  = _toString(ncomponents, "component%d");
+
+      double box_width    = x_comp.width() - 2 * frame_width;
+      max_component_width = fmax(max_component_width, box_width);
+
+      Box    c_box{box_width / 2, x_comp.length() / 2, x_comp.thickness() / 2};
+      Volume c_vol{c_nam, c_box, description.material(x_comp.materialStr())};
+
+      pv = m_vol.placeVolume(c_vol, Position(0, 0, thickness_sum + x_comp.thickness() / 2.0));
+
+      c_vol.setRegion(description, x_comp.regionStr());
+      c_vol.setLimitSet(description, x_comp.limitsStr());
+      c_vol.setVisAttributes(description, x_comp.visStr());
+      if (x_comp.isSensitive()) {
+        c_vol.setSensitiveDetector(sens);
+        sensitives[m_nam].push_back(pv);
+        module_thicknesses[m_nam] = {thickness_so_far + x_comp.thickness() / 2.0,
+                                     total_thickness - thickness_so_far - x_comp.thickness() / 2.0};
+        // -------- create a measurement plane for the tracking surface attched to the sensitive volume -----
+        rec::Vector3D u(0., 1., 0.);
+        rec::Vector3D v(0., 0., 1.);
+        rec::Vector3D n(1., 0., 0.);
+
+        // compute the inner and outer thicknesses that need to be assigned to the tracking surface
+        // depending on wether the support is above or below the sensor
+        double inner_thickness = module_thicknesses[m_nam][0];
+        double outer_thickness = module_thicknesses[m_nam][1];
+
+        rec::SurfaceType type(rec::SurfaceType::Sensitive);
+
+        rec::VolPlane surf(c_vol, type, inner_thickness, outer_thickness, u, v, n); //,o ) ;
+        volplane_surfaces[m_nam].push_back(surf);
+      }
+      thickness_sum += x_comp.thickness();
+      thickness_so_far += x_comp.thickness();
+    }
+    // Now add-on the frame
+    if (x_mod.hasChild("frame")) {
+      xml_comp_t m_frame         = x_mod.child(_U(frame));
+      double     frame_thickness = getAttrOrDefault<double>(m_frame, _U(thickness), total_thickness);
+
+      Box lframe_box{m_frame.width() / 2., m_frame.length() / 2., frame_thickness / 2.};
+      Box rframe_box{m_frame.width() / 2., m_frame.length() / 2., frame_thickness / 2.};
+
+      // Keep track of frame with so we can adjust the module bars appropriately
+
+      Volume lframe_vol{"left_frame", lframe_box, description.material(m_frame.materialStr())};
+      Volume rframe_vol{"right_frame", rframe_box, description.material(m_frame.materialStr())};
+      lframe_vol.setVisAttributes(description, m_frame.visStr());
+      rframe_vol.setVisAttributes(description, m_frame.visStr());
+
+      m_vol.placeVolume(lframe_vol, Position(frame_width / 2. + max_component_width / 2, 0.,
+                                             frame_thickness / 2. - total_thickness / 2.0));
+      m_vol.placeVolume(rframe_vol, Position(-frame_width / 2. - max_component_width / 2, 0.,
+                                             frame_thickness / 2. - total_thickness / 2.0));
+    }
+  }
+
+  // now build the layers
+  for (xml_coll_t li(x_det, _U(layer)); li; ++li) {
+    xml_comp_t x_layer  = li;
+    xml_comp_t x_barrel = x_layer.child(_U(barrel_envelope));
+    xml_comp_t x_layout = x_layer.child(_U(rphi_layout));
+    xml_comp_t z_layout = x_layer.child(_U(z_layout)); // Get the <z_layout> element.
+    int        lay_id   = x_layer.id();
+    string     m_nam    = x_layer.moduleStr();
+    string     lay_nam  = _toString(x_layer.id(), "layer%d");
+    Tube       lay_tub(x_barrel.inner_r(), x_barrel.outer_r(), x_barrel.z_length() / 2.0);
+    Volume     lay_vol(lay_nam, lay_tub, air); // Create the layer envelope volume.
+    lay_vol.setVisAttributes(description, x_layer.visStr());
+
+    double phi0     = x_layout.phi0();     // Starting phi of first module.
+    double phi_tilt = x_layout.phi_tilt(); // Phi tilt of a module.
+    double rc       = x_layout.rc();       // Radius of the module center.
+    int    nphi     = x_layout.nphi();     // Number of modules in phi.
+    double rphi_dr  = x_layout.dr();       // The delta radius of every other module.
+    double phi_incr = (M_PI * 2) / nphi;   // Phi increment for one module.
+    double phic     = phi0;                // Phi of the module center.
+    double z0       = z_layout.z0();       // Z position of first module in phi.
+    double nz       = z_layout.nz();       // Number of modules to place in z.
+    double z_dr     = z_layout.dr();       // Radial displacement parameter, of every other module.
+
+    Volume      module_env = volumes[m_nam];
+    DetElement  lay_elt(sdet, _toString(x_layer.id(), "layer%d"), lay_id);
+    Placements& sensVols = sensitives[m_nam];
+
+    // Z increment for module placement along Z axis.
+    // Adjust for z0 at center of module rather than
+    // the end of cylindrical envelope.
+    double z_incr = nz > 1 ? (2.0 * z0) / (nz - 1) : 0.0;
+    // Starting z for module placement along Z axis.
+    double module_z = -z0;
+    int    module   = 1;
+
+    // Loop over the number of modules in phi.
+    for (int ii = 0; ii < nphi; ii++) {
+      double dx = z_dr * std::cos(phic + phi_tilt); // Delta x of module position.
+      double dy = z_dr * std::sin(phic + phi_tilt); // Delta y of module position.
+      double x  = rc * std::cos(phic);              // Basic x module position.
+      double y  = rc * std::sin(phic);              // Basic y module position.
+
+      // Loop over the number of modules in z.
+      for (int j = 0; j < nz; j++) {
+        string     module_name = _toString(module, "module%d");
+        DetElement mod_elt(lay_elt, module_name, module);
+
+        Transform3D tr(RotationZYX(0, ((M_PI / 2) - phic - phi_tilt), -M_PI / 2), Position(x, y, module_z));
+
+        pv = lay_vol.placeVolume(module_env, tr);
+        pv.addPhysVolID("module", module);
+        pv.addPhysVolID("section", j);
+        mod_elt.setPlacement(pv);
+        for (size_t ic = 0; ic < sensVols.size(); ++ic) {
+          PlacedVolume sens_pv = sensVols[ic];
+          DetElement   comp_de(mod_elt, std::string("de_") + sens_pv.volume().name(), module);
+          comp_de.setPlacement(sens_pv);
+          rec::volSurfaceList(comp_de)->push_back(volplane_surfaces[m_nam][ic]);
+        }
+
+        /// Increase counters etc.
+        module++;
+        // Adjust the x and y coordinates of the module.
+        x += dx;
+        y += dy;
+        // Flip sign of x and y adjustments.
+        dx *= -1;
+        dy *= -1;
+        // Add z increment to get next z placement pos.
+        module_z += z_incr;
+      }
+      phic += phi_incr; // Increment the phi placement of module.
+      rc += rphi_dr;    // Increment the center radius according to dr parameter.
+      rphi_dr *= -1;    // Flip sign of dr parameter.
+      module_z = -z0;   // Reset the Z placement parameter for module.
+    }
+    // Create the PhysicalVolume for the layer.
+    pv = assembly.placeVolume(lay_vol); // Place layer in mother
+    pv.addPhysVolID("layer", lay_id);   // Set the layer ID.
+    lay_elt.setAttributes(description, lay_vol, x_layer.regionStr(), x_layer.limitsStr(), x_layer.visStr());
+    lay_elt.setPlacement(pv);
+  }
+  sdet.setAttributes(description, assembly, x_det.regionStr(), x_det.limitsStr(), x_det.visStr());
+  assembly.setVisAttributes(description.invisible());
+  pv = description.pickMotherVolume(sdet).placeVolume(assembly, Position(0, 0, dirc_pos.z()));
+  pv.addPhysVolID("system", det_id); // Set the subdetector system ID.
+  sdet.setPlacement(pv);
+  return sdet;
+}
+
+//@}
+// clang-format off
+DECLARE_DETELEMENT(BarrelBarDetectorWithSideFrame, create_BarrelBarDetectorWithSideFrame)
+DECLARE_DETELEMENT(athena_FakeDIRC, create_BarrelBarDetectorWithSideFrame)

src/BarrelCalorimeterInterlayers_geo.cpp

+// Detector plugin to support a hybrid central barrel calorimeter
+// The detector consists of interlayers of Pb/ScFi (segmentation in global r, phi) and W/Si (segmentation in local x, y)
+// Assembly is used as the envelope so two different detectors can be interlayered with each other
+//
+//
+// 06/19/2021: Implementation of the Sci Fiber geometry. M. Żurek
+// 07/09/2021: Support interlayers between multiple detectors. C. Peng
+// 07/23/2021: Add assemblies as mother volumes of fibers to reduce the number of daughter volumes. C. Peng, M. Żurek
+//     Reference: TGeo performance issue with large number of daughter volumes
+//     https://indico.cern.ch/event/967418/contributions/4075358/attachments/2128099/3583278/201009_shKo_dd4hep.pdf
+// 07/24/2021: Changed support implementation to avoid too many uses of boolean geometries. DAWN view seems to have
+//     issue dealing with it. C. Peng
+
+#include "DD4hep/DetFactoryHelper.h"
+#include "XML/Layering.h"
+#include "Math/Point2D.h"
+#include "TGeoPolygon.h"
+
+using namespace std;
+using namespace dd4hep;
+using namespace dd4hep::detail;
+
+typedef ROOT::Math::XYPoint Point;
+// fiber placement helpers, defined below
+vector<vector<Point>> fiberPositions(double radius, double x_spacing, double z_spacing,
+                                     double x, double z, double phi, double spacing_tol = 1e-2);
+std::pair<int, int> getNdivisions(double x, double z, double dx, double dz);
+vector<tuple<int, Point, Point, Point, Point>> gridPoints(int div_x, int div_z, double x, double z, double phi);
+
+// geometry helpers
+void buildFibers(Detector& desc, SensitiveDetector &sens, Volume &mother, xml_comp_t x_fiber,
+                 const std::tuple<double, double, double, double> &dimensions);
+void buildSupport(Detector& desc, Volume &mother, xml_comp_t x_support,
+                  const std::tuple<double, double, double, double> &dimensions);
+
+
+// barrel ecal layers contained in an assembly
+static Ref_t create_detector(Detector& desc, xml_h e, SensitiveDetector sens)  {
+  Layering      layering (e);
+  xml_det_t     x_det     = e;
+  Material      air       = desc.air();
+  int           det_id    = x_det.id();
+  string        det_name  = x_det.nameStr();
+  double        offset    = x_det.attr<double>(_Unicode(offset));
+  xml_comp_t    x_dim     = x_det.dimensions();
+  int           nsides    = x_dim.numsides();
+  double        inner_r   = x_dim.rmin();
+  double        dphi      = (2*M_PI/nsides);
+  double        hphi      = dphi/2;
+
+  DetElement    sdet      (det_name, det_id);
+  Volume        motherVol = desc.pickMotherVolume(sdet);
+
+  Assembly      envelope  (det_name);
+  Transform3D   tr        = Translation3D(0, 0, offset) * RotationZ(hphi);
+  PlacedVolume  env_phv   = motherVol.placeVolume(envelope, tr);
+  sens.setType("calorimeter");
+
+  env_phv.addPhysVolID("system",det_id);
+  sdet.setPlacement(env_phv);
+
+  // build a single stave
+  DetElement    stave_det("stave0", det_id);
+  Assembly      mod_vol("stave");
+
+  // keep tracking of the total thickness
+  double l_pos_z = inner_r;
+  { // =====  buildBarrelStave(desc, sens, module_volume) =====
+    // Parameters for computing the layer X dimension:
+    double tan_hphi = std::tan(hphi);
+    double l_dim_y  = x_dim.z()/2.;
+
+    // Loop over the sets of layer elements in the detector.
+    int l_num = 1;
+    for(xml_coll_t li(x_det, _U(layer)); li; ++li)  {
+      xml_comp_t x_layer = li;
+      int repeat = x_layer.repeat();
+      double l_space_between = dd4hep::getAttrOrDefault(x_layer, _Unicode(space_between), 0.);
+      double l_space_before = dd4hep::getAttrOrDefault(x_layer, _Unicode(space_before), 0.);
+      l_pos_z += l_space_before;
+      // Loop over number of repeats for this layer.
+      for (int j = 0; j < repeat; j++)    {
+        string l_name = Form("layer%d", l_num);
+        double l_thickness = layering.layer(l_num - 1)->thickness();  // Layer's thickness.
+        double l_dim_x = tan_hphi* l_pos_z;
+        l_pos_z += l_thickness;
+
+        Position   l_pos(0, 0, l_pos_z - l_thickness/2.);      // Position of the layer.
+	    double l_trd_x1 = l_dim_x;
+	    double l_trd_x2 = l_dim_x + l_thickness*tan_hphi;
+	    double l_trd_y1 = l_dim_y;
+	    double l_trd_y2 = l_trd_y1;
+	    double l_trd_z  = l_thickness/2;
+        Trapezoid  l_shape(l_trd_x1, l_trd_x2, l_trd_y1, l_trd_y2, l_trd_z);
+        Volume     l_vol(l_name, l_shape, air);
+        DetElement layer(stave_det, l_name, det_id);
+
+        // Loop over the sublayers or slices for this layer.
+        int s_num = 1;
+        double s_pos_z = -(l_thickness / 2.);
+        for(xml_coll_t si(x_layer,_U(slice)); si; ++si)  {
+          xml_comp_t x_slice = si;
+          string     s_name  = Form("slice%d", s_num);
+          double     s_thick = x_slice.thickness();
+	      double s_trd_x1 = l_dim_x + (s_pos_z + l_thickness/2)*tan_hphi;
+	      double s_trd_x2 = l_dim_x + (s_pos_z + l_thickness/2 + s_thick)*tan_hphi;
+	      double s_trd_y1 = l_trd_y1;
+	      double s_trd_y2 = s_trd_y1;
+	      double s_trd_z  = s_thick/2.;
+          Trapezoid  s_shape(s_trd_x1, s_trd_x2, s_trd_y1, s_trd_y2, s_trd_z);
+          Volume     s_vol(s_name, s_shape, desc.material(x_slice.materialStr()));
+          DetElement slice(layer, s_name, det_id);
+
+          // build fibers
+          if (x_slice.hasChild(_Unicode(fiber))) {
+              buildFibers(desc, sens, s_vol, x_slice.child(_Unicode(fiber)), {s_trd_x1, s_thick, l_dim_y, hphi});
+          }
+
+          if ( x_slice.isSensitive() ) {
+            s_vol.setSensitiveDetector(sens);
+          }
+          s_vol.setAttributes(desc, x_slice.regionStr(), x_slice.limitsStr(), x_slice.visStr());
+
+          // Slice placement.
+          PlacedVolume slice_phv = l_vol.placeVolume(s_vol, Position(0, 0, s_pos_z + s_thick/2));
+          slice_phv.addPhysVolID("slice", s_num);
+          slice.setPlacement(slice_phv);
+          // Increment Z position of slice.
+          s_pos_z += s_thick;
+          ++s_num;
+        }
+
+        // Set region, limitset, and vis of layer.
+        l_vol.setAttributes(desc, x_layer.regionStr(), x_layer.limitsStr(), x_layer.visStr());
+
+        PlacedVolume layer_phv = mod_vol.placeVolume(l_vol, l_pos);
+        layer_phv.addPhysVolID("layer", l_num);
+        layer.setPlacement(layer_phv);
+        // Increment to next layer Z position. Do not add space_between for the last layer
+        if (j < repeat - 1) {
+          l_pos_z += l_space_between;
+        }
+        ++l_num;
+      }
+    }
+  }
+  // Phi start for a stave.
+  double phi = M_PI / nsides;
+  // Create nsides staves.
+  for (int i = 0; i < nsides; i++, phi -= dphi)      { // i is module number
+    // Compute the stave position
+    Transform3D tr(RotationZYX(0, phi, M_PI*0.5), Translation3D(0, 0, 0));
+    PlacedVolume pv = envelope.placeVolume(mod_vol, tr);
+    pv.addPhysVolID("module", i + 1);
+    DetElement sd = (i == 0) ? stave_det : stave_det.clone(Form("stave%d", i));
+    sd.setPlacement(pv);
+    sdet.add(sd);
+  }
+
+  // optional stave support
+  if (x_det.hasChild(_U(staves))) {
+    xml_comp_t x_staves = x_det.staves();
+    mod_vol.setVisAttributes(desc.visAttributes(x_staves.visStr()));
+    if (x_staves.hasChild(_U(support))) {
+      buildSupport(desc, mod_vol, x_staves.child(_U(support)), {inner_r, l_pos_z, x_dim.z(), hphi});
+    }
+  }
+
+  // Set envelope volume attributes.
+  envelope.setAttributes(desc, x_det.regionStr(), x_det.limitsStr(), x_det.visStr());
+  return sdet;
+}
+
+void buildFibers(Detector& desc, SensitiveDetector &sens, Volume &s_vol, xml_comp_t x_fiber,
+                 const std::tuple<double, double, double, double> &dimensions)
+{
+  auto [s_trd_x1, s_thick, s_length, hphi] = dimensions;
+  double f_radius = getAttrOrDefault(x_fiber, _U(radius), 0.1 * cm);
+  double f_spacing_x = getAttrOrDefault(x_fiber, _Unicode(spacing_x), 0.122 * cm);
+  double f_spacing_z = getAttrOrDefault(x_fiber, _Unicode(spacing_z), 0.134 * cm);
+  std::string f_id_grid = getAttrOrDefault(x_fiber, _Unicode(identifier_grid), "grid");
+  std::string f_id_fiber = getAttrOrDefault(x_fiber, _Unicode(identifier_fiber), "fiber");
+
+  // Set up the readout grid for the fiber layers
+  // Trapezoid is divided into segments with equal dz and equal number of divisions in x
+  // Every segment is a polygon that can be attached later to the lightguide
+  // The grid size is assumed to be ~2x2 cm (starting values). This is to be larger than
+  // SiPM chip (for GlueX 13mmx13mm: 4x4 grid 3mmx3mm with 3600 50×50 μm pixels each)
+  // See, e.g., https://arxiv.org/abs/1801.03088 Fig. 2d
+
+  // Calculate number of divisions
+  auto grid_div = getNdivisions(s_trd_x1, s_thick, 2.0*cm, 2.0*cm);
+  // Calculate polygonal grid coordinates (vertices)
+  auto grid_vtx = gridPoints(grid_div.first, grid_div.second, s_trd_x1, s_thick, hphi);
+  Tube f_tube(0, f_radius, s_length);
+  Volume f_vol("fiber_vol", f_tube, desc.material(x_fiber.materialStr()));
+
+  vector<int> f_id_count(grid_div.first*grid_div.second, 0);
+  auto f_pos = fiberPositions(f_radius, f_spacing_x, f_spacing_z, s_trd_x1, s_thick, hphi);
+  // std::cout << f_pos.size() << " lines, ~" << f_pos.front().size() << " fibers each line" << std::endl;
+  for (size_t il = 0; il < f_pos.size(); ++il) {
+    auto &line = f_pos[il];
+    if (line.empty()) {
+      continue;
+    }
+    double l_pos_y = line.front().y();
+    // use assembly as intermediate volume container to reduce number of daughter volumes
+    Assembly lfibers(Form("fiber_array_line_%lu", il));
+    for (auto &p : line) {
+      int f_grid_id = -1;
+      int f_id = -1;
+      // Check to which grid fiber belongs to
+      for (auto &poly_vtx : grid_vtx) {
+        if (p.y() != l_pos_y) {
+          std::cerr << Form("Expected the same y position from a same line: %.2f, but got %.2f", l_pos_y, p.y())
+                    << std::endl;
+          continue;
+        }
+        auto [grid_id, vtx_a, vtx_b, vtx_c, vtx_d] = poly_vtx;
+        double poly_x[4] = {vtx_a.x(), vtx_b.x(), vtx_c.x(), vtx_d.x()};
+        double poly_y[4] = {vtx_a.y(), vtx_b.y(), vtx_c.y(), vtx_d.y()};
+        double f_xy[2] = {p.x(), p.y()};
+
+        TGeoPolygon poly(4);
+        poly.SetXY(poly_x, poly_y);
+        poly.FinishPolygon();
+
+        if(poly.Contains(f_xy)) {
+          f_grid_id = grid_id;
+          f_id = f_id_count[grid_id];
+          f_id_count[grid_id]++;
+        }
+      }
+
+      if ( x_fiber.isSensitive() ) {
+        f_vol.setSensitiveDetector(sens);
+      }
+      f_vol.setAttributes(desc, x_fiber.regionStr(), x_fiber.limitsStr(), x_fiber.visStr());
+
+      // Fiber placement
+      // Transform3D f_tr(RotationZYX(0,0,M_PI*0.5),Position(p.x(), 0, p.y()));
+      // PlacedVolume fiber_phv = s_vol.placeVolume(f_vol, Position(p.x(), 0., p.y()));
+      PlacedVolume fiber_phv = lfibers.placeVolume(f_vol, Position(p.x(), 0., 0.));
+      fiber_phv.addPhysVolID(f_id_grid, f_grid_id + 1).addPhysVolID(f_id_fiber, f_id + 1);
+    }
+    lfibers.ptr()->Voxelize("");
+    Transform3D l_tr(RotationZYX(0,0,M_PI*0.5),Position(0., 0, l_pos_y));
+    s_vol.placeVolume(lfibers, l_tr);
+  }
+}
+
+// DAWN view seems to have some issue with overlapping solids even if they were unions
+// The support is now built without overlapping
+void buildSupport(Detector& desc, Volume &mod_vol, xml_comp_t x_support,
+                  const std::tuple<double, double, double, double> &dimensions)
+{
+  auto [inner_r, l_pos_z, stave_length, hphi] = dimensions;
+
+  double support_thickness = getAttrOrDefault(x_support, _Unicode(thickness), 5. * cm);
+  double beam_thickness    = getAttrOrDefault(x_support, _Unicode(beam_thickness), support_thickness/4.);
+  // sanity check
+  if (beam_thickness > support_thickness/3.) {
+    std::cerr << Form("beam_thickness (%.2f) cannot be greater than support_thickness/3 (%.2f), shrink it to fit",
+                      beam_thickness, support_thickness/3.) << std::endl;
+    beam_thickness = support_thickness/3.;
+  }
+  double   trd_x1_support = std::tan(hphi) * l_pos_z;
+  double   trd_x2_support = std::tan(hphi) * (l_pos_z + support_thickness);
+  double   trd_y          = stave_length / 2.;
+  Assembly env_vol        ("support_envelope");
+
+  double grid_size        = getAttrOrDefault(x_support, _Unicode(grid_size), 25. * cm);
+  int    n_cross_supports = std::floor(trd_y - beam_thickness)/grid_size;
+  // number of "beams" running the length of the stave.
+  // @TODO make it configurable
+  int    n_beams          = getAttrOrDefault(x_support, _Unicode(n_beams), 3);;
+  double beam_width       = 2. * trd_x1_support / (n_beams + 1); // quick hack to make some gap between T beams
+  double beam_gap         = getAttrOrDefault(x_support, _Unicode(beam_gap), 3.*cm);
+
+  // build T-shape beam
+  double beam_space_x    = beam_width + beam_gap;
+  double beam_space_z    = support_thickness - beam_thickness;
+  double cross_thickness = support_thickness - beam_thickness;
+  double beam_pos_z      = beam_thickness / 2.;
+  double beam_center_z   = support_thickness / 2. - beam_pos_z;
+
+  Box        beam_vert_s(beam_thickness / 2., trd_y, cross_thickness / 2.);
+  Box        beam_hori_s(beam_width / 2., trd_y, beam_thickness / 2.);
+  UnionSolid T_beam_s(beam_hori_s, beam_vert_s, Position(0., 0., support_thickness / 2.));
+  Volume H_beam_vol("H_beam", T_beam_s, desc.material(x_support.materialStr()));
+  H_beam_vol.setVisAttributes(desc, x_support.visStr());
+  // place H beams first
+  double beam_start_x = - (n_beams - 1) * (beam_width + beam_gap) / 2.;
+  for (int i = 0; i < n_beams; ++i) {
+    Position beam_pos(beam_start_x + i * (beam_width + beam_gap), 0., - support_thickness / 2. + beam_pos_z);
+    env_vol.placeVolume(H_beam_vol, beam_pos);
+  }
+
+  // place central crossing beams that connects the H beams
+  double cross_x = beam_space_x - beam_thickness;
+  Box cross_s(cross_x / 2., beam_thickness / 2., cross_thickness / 2.);
+  Volume cross_vol("cross_center_beam", cross_s, desc.material(x_support.materialStr()));
+  cross_vol.setVisAttributes(desc, x_support.visStr());
+  for (int i = 0; i < n_beams - 1; ++i) {
+    env_vol.placeVolume(cross_vol, Position(beam_start_x + beam_space_x * (i + 0.5), 0., beam_pos_z));
+    for (int j = 1; j < n_cross_supports; j++) {
+      env_vol.placeVolume(cross_vol, Position(beam_start_x + beam_space_x * (i + 0.5), -j * grid_size, beam_pos_z));
+      env_vol.placeVolume(cross_vol, Position(beam_start_x + beam_space_x * (i + 0.5), j * grid_size, beam_pos_z));
+    }
+  }
+
+  // place edge crossing beams that connects the neighbour support
+  // @TODO: connection part is still using boolean volumes, maybe problematic to DAWN
+  double    cross_edge_x = trd_x1_support + beam_start_x - beam_thickness / 2.;
+  double    cross_trd_x1 = cross_edge_x + std::tan(hphi) * beam_thickness;
+  double    cross_trd_x2 = cross_trd_x1 + 2.* std::tan(hphi) * cross_thickness;
+  double    edge_pos_x   = beam_start_x - cross_trd_x1 / 2. - beam_thickness / 2;
+  Trapezoid cross_s2_trd (cross_trd_x1 / 2., cross_trd_x2 / 2.,
+                          beam_thickness / 2., beam_thickness / 2., cross_thickness / 2.);
+  Box       cross_s2_box ((cross_trd_x2 - cross_trd_x1)/4., beam_thickness / 2., cross_thickness / 2.);
+  SubtractionSolid cross_s2(cross_s2_trd, cross_s2_box, Position((cross_trd_x2 + cross_trd_x1)/4., 0., 0.));
+  Volume cross_vol2("cross_edge_beam", cross_s2, desc.material(x_support.materialStr()));
+  cross_vol2.setVisAttributes(desc, x_support.visStr());
+  env_vol.placeVolume(cross_vol2, Position(edge_pos_x, 0., beam_pos_z));
+  env_vol.placeVolume(cross_vol2, Transform3D(Translation3D(-edge_pos_x, 0., beam_pos_z) * RotationZ(M_PI)));
+  for (int j = 1; j < n_cross_supports; j++) {
+    env_vol.placeVolume(cross_vol2, Position(edge_pos_x, -j * grid_size, beam_pos_z));
+    env_vol.placeVolume(cross_vol2, Position(edge_pos_x, j * grid_size, beam_pos_z));
+    env_vol.placeVolume(cross_vol2, Transform3D(Translation3D(-edge_pos_x, -j * grid_size, beam_pos_z) * RotationZ(M_PI)));
+    env_vol.placeVolume(cross_vol2, Transform3D(Translation3D(-edge_pos_x, j * grid_size, beam_pos_z) * RotationZ(M_PI)));
+  }
+
+  mod_vol.placeVolume(env_vol, Position(0.0, 0.0, l_pos_z + support_thickness/2.));
+}
+
+// Fill fiber lattice into trapezoid starting from position (0,0) in x-z coordinate system
+vector<vector<Point>> fiberPositions(double radius, double x_spacing, double z_spacing, double x, double z, double phi,
+                                     double spacing_tol)
+{
+  // z_spacing - distance between fiber layers in z
+  // x_spacing - distance between fiber centers in x
+  // x - half-length of the shorter (bottom) base of the trapezoid
+  // z - height of the trapezoid
+  // phi - angle between z and trapezoid arm
+
+  vector<vector<Point>> positions;
+  int z_layers = floor((z / 2 - radius - spacing_tol) / z_spacing); // number of layers that fit in z/2
+
+  double z_pos = 0.;
+  double x_pos = 0.;
+
+  for (int l = -z_layers; l < z_layers + 1; l++) {
+    vector<Point> xline;
+    z_pos                = l * z_spacing;
+    double x_max         = x + (z / 2. + z_pos) * tan(phi) - spacing_tol; // calculate max x at particular z_pos
+    (l % 2 == 0) ? x_pos = 0. : x_pos = x_spacing / 2;                    // account for spacing/2 shift
+
+    while (x_pos < (x_max - radius)) {
+      xline.push_back(Point(x_pos, z_pos));
+      if (x_pos != 0.)
+        xline.push_back(Point(-x_pos, z_pos)); // using symmetry around x=0
+      x_pos += x_spacing;
+    }
+    // Sort fiber IDs for a better organization
+    sort(xline.begin(), xline.end(), [](const Point& p1, const Point& p2) { return p1.x() < p2.x(); });
+    positions.emplace_back(std::move(xline));
+  }
+  return positions;
+}
+
+// Calculate number of divisions for the readout grid for the fiber layers
+std::pair<int, int> getNdivisions(double x, double z, double dx, double dz)
+{
+  // x and z defined as in vector<Point> fiberPositions
+  // dx, dz - size of the grid in x and z we want to get close to with the polygons
+  // See also descripltion when the function is called
+
+  double SiPMsize = 13.0 * mm;
+  double grid_min = SiPMsize + 3.0 * mm;
+
+  if (dz < grid_min) {
+    dz = grid_min;
+  }
+
+  if (dx < grid_min) {
+    dx = grid_min;
+  }
+
+  int nfit_cells_z = floor(z / dz);
+  int n_cells_z    = nfit_cells_z;
+
+  if (nfit_cells_z == 0)
+    n_cells_z++;
+
+  int nfit_cells_x = floor((2 * x) / dx);
+  int n_cells_x    = nfit_cells_x;
+
+  if (nfit_cells_x == 0)
+    n_cells_x++;
+
+  return std::make_pair(n_cells_x, n_cells_z);
+}
+
+// Calculate dimensions of the polygonal grid in the cartesian coordinate system x-z
+vector<tuple<int, Point, Point, Point, Point>> gridPoints(int div_x, int div_z, double x, double z, double phi)
+{
+  // x, z and phi defined as in vector<Point> fiberPositions
+  // div_x, div_z - number of divisions in x and z
+  double dz = z / div_z;
+
+  std::vector<std::tuple<int, Point, Point, Point, Point>> points;
+
+  for (int iz = 0; iz < div_z + 1; iz++) {
+    for (int ix = 0; ix < div_x + 1; ix++) {
+      double A_z = -z / 2 + iz * dz;
+      double B_z = -z / 2 + (iz + 1) * dz;
+
+      double len_x_for_z        = 2 * (x + iz * dz * tan(phi));
+      double len_x_for_z_plus_1 = 2 * (x + (iz + 1) * dz * tan(phi));
+
+      double dx_for_z        = len_x_for_z / div_x;
+      double dx_for_z_plus_1 = len_x_for_z_plus_1 / div_x;
+
+      double A_x = -len_x_for_z / 2. + ix * dx_for_z;
+      double B_x = -len_x_for_z_plus_1 / 2. + ix * dx_for_z_plus_1;
+
+      double C_z = B_z;
+      double D_z = A_z;
+      double C_x = B_x + dx_for_z_plus_1;
+      double D_x = A_x + dx_for_z;
+
+      int id = ix + div_x * iz;
+
+      auto A = Point(A_x, A_z);
+      auto B = Point(B_x, B_z);
+      auto C = Point(C_x, C_z);
+      auto D = Point(D_x, D_z);
+
+      // vertex points filled in the clock-wise direction
+      points.push_back(make_tuple(id, A, B, C, D));
+    }
+  }
+
+  return points;
+}
+
+DECLARE_DETELEMENT(athena_EcalBarrelInterlayers, create_detector)
+

src/BarrelTrackerWithFrame_geo.cpp

-/** \addtogroup VertexTracker Vertex Trackers
- * \brief Type: **SiVertexBarrel**.
+/** \addtogroup Trackers Trackers
+ * \brief Type: **BarrelTrackerWithFrame**.
  * \author W. Armstrong
- * \ingroup trackers
- *
  *
- * \code
- * \endcode
+ * \ingroup trackers
  *
  * @{
  */
@@ -15,131 +12,185 @@
 #include "DDRec/Surface.h"
 #include "DDRec/DetectorData.h"
 #include "XML/Layering.h"
-#include "Acts/Plugins/DD4hep/ActsExtension.hpp"
-#include "Acts/Surfaces/PlanarBounds.hpp"
-#include "Acts/Surfaces/RectangleBounds.hpp"
-#include "Acts/Surfaces/TrapezoidBounds.hpp"
-#include "Acts/Definitions/Units.hpp"
+#include "XML/Utilities.h"
+#include <array>
 
+#if defined(USE_ACTSDD4HEP)
+#include "ActsDD4hep/ActsExtension.hpp"
+#include "ActsDD4hep/ConvertMaterial.hpp"
+#else
+#include "Acts/Plugins/DD4hep/ActsExtension.hpp"
+#include "Acts/Plugins/DD4hep/ConvertDD4hepMaterial.hpp"
+#endif
 
 using namespace std;
 using namespace dd4hep;
 using namespace dd4hep::rec;
 using namespace dd4hep::detail;
 
-static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector sens) {
+
+/** Barrel Tracker with space frame.
+ *
+ * - Optional "support" tag within the detector element.
+ * 
+ * The shapes are created using createShape which can be one of many basic geomtries. 
+ * See the examples Check_shape_*.xml in 
+ * [dd4hep's examples/ClientTests/compact](https://github.com/AIDASoft/DD4hep/tree/master/examples/ClientTests/compact)
+ * directory.
+ *
+ *
+ * - Optional "frame" tag within the module element.
+ *
+ * \ingroup trackers
+ *
+ * \code
+ * \endcode
+ * 
+ *
+ * @author Whitney Armstrong
+ */
+static Ref_t create_BarrelTrackerWithFrame(Detector& description, xml_h e, SensitiveDetector sens) {
   typedef vector<PlacedVolume> Placements;
   xml_det_t                    x_det    = e;
   Material                     air      = description.air();
   int                          det_id   = x_det.id();
   string                       det_name = x_det.nameStr();
   DetElement                   sdet(det_name, det_id);
-  //Assembly                     assembly(det_name);
-  map<string, Volume>          volumes;
-  map<string, Placements>      sensitives;
-  map<string, xml_h>      xmleles;
-  PlacedVolume                 pv;
-  dd4hep::xml::Dimension dimensions(x_det.dimensions());
-
-  Acts::ActsExtension* detWorldExt = new Acts::ActsExtension();
-  detWorldExt->addType("barrel", "detector");
-  sdet.addExtension<Acts::ActsExtension>(detWorldExt);
-  Tube topVolumeShape(dimensions.rmin(), dimensions.rmax(), dimensions.length() * 0.5);
-  Volume assembly(det_name,topVolumeShape,air);
 
+  map<string, Volume>                volumes;
+  map<string, Placements>            sensitives;
+  map<string, std::vector<VolPlane>> volplane_surfaces;
+  map<string, std::array<double, 2>> module_thicknesses;
 
-//  The Cold Plate is approximately 30 mm wide and is based on the same carbon-ply layup
-//as for the IB Stave. Two pipes with an inner diameter of 2.67 mm and a wall thickness
-//of 64 μm have been used. The two pipes are interconnected at one end of the Cold Plate
-//providing a loop, whose inlet and outlet are on the same side and correspond to the
-//
-//requirements have led to an equilateral section of the frame with a 42 mm wide side, that
-//provides almost the same rigidity for all the possible Stave positions.
-//
-//                     module mat                um
-//
-//Aluminium                                      50
-//Polyimide                                      100
-//Carbon fibre                                   120
-//Silicon                                        50
-//Eccobond-45                                    100
-//
-//Metal layers                 Aluminium        200
-//Insulating layers            Polyimide        200
-//Glue Cooling tube wall       Eccobond-45      100
-//
-//Carbon fleece            40
-//Carbon paper             30
-//Polyimide                64
-//Water                    
-//Carbon fibre             120
-//Eccobond-45              100
+  PlacedVolume           pv;
+  dd4hep::xml::Dimension dimensions(x_det.dimensions());
 
+  // ACTS extension
+  {
+    Acts::ActsExtension* detWorldExt = new Acts::ActsExtension();
+    detWorldExt->addType("barrel", "detector");
+    // Add the volume boundary material if configured
+    for (xml_coll_t bmat(x_det, _Unicode(boundary_material)); bmat; ++bmat) {
+      xml_comp_t x_boundary_material = bmat;
+      Acts::xmlToProtoSurfaceMaterial(x_boundary_material, *detWorldExt, "boundary_material");
+    }
+    sdet.addExtension<Acts::ActsExtension>(detWorldExt);
+  }
 
+  Tube topVolumeShape(dimensions.rmin(), dimensions.rmax(), dimensions.length() * 0.5);
+  Volume assembly(det_name,topVolumeShape,air);
 
   sens.setType("tracker");
 
+  // Loop over the suports
+  for (xml_coll_t su(x_det, _U(support)); su; ++su) {
+    xml_comp_t x_support = su;
+    double      support_thickness = getAttrOrDefault(x_support, _U(thickness), 2.0 * mm);
+    double      support_length    = getAttrOrDefault(x_support, _U(length), 2.0 * mm);
+    double      support_rmin      = getAttrOrDefault(x_support, _U(rmin), 2.0 * mm);
+    double      support_zstart    = getAttrOrDefault(x_support, _U(zstart), 2.0 * mm);
+    std::string support_name      = getAttrOrDefault<std::string>(x_support, _Unicode(name), "support_tube");
+    std::string support_vis       = getAttrOrDefault<std::string>(x_support, _Unicode(vis), "AnlRed");
+    xml_dim_t  pos        (x_support.child(_U(position), false));
+    xml_dim_t  rot        (x_support.child(_U(rotation), false));
+    Solid support_solid;
+    if(x_support.hasChild("shape")){
+      xml_comp_t shape(x_support.child(_U(shape)));
+      string     shape_type = shape.typeStr();
+      support_solid  = xml::createShape(description, shape_type, shape);
+    } else {
+      support_solid = Tube(support_rmin, support_rmin + support_thickness, support_length / 2);
+    }
+    Transform3D tr = Transform3D(Rotation3D(),Position(0,0,(support_zstart + support_length / 2)));
+    if ( pos.ptr() && rot.ptr() )  {
+      Rotation3D  rot3D(RotationZYX(rot.z(0),rot.y(0),rot.x(0)));
+      Position    pos3D(pos.x(0),pos.y(0),pos.z(0));
+      tr = Transform3D(rot3D, pos3D);
+    }
+    else if ( pos.ptr() )  {
+      tr = Transform3D(Rotation3D(),Position(pos.x(0),pos.y(0),pos.z(0)));
+    }
+    else if ( rot.ptr() )  {
+      Rotation3D rot3D(RotationZYX(rot.z(0),rot.y(0),rot.x(0)));
+      tr = Transform3D(rot3D,Position());
+    }
+    Material    support_mat       = description.material(x_support.materialStr());
+    Volume      support_vol(support_name, support_solid, support_mat);
+    support_vol.setVisAttributes(description.visAttributes(support_vis));
+    pv = assembly.placeVolume(support_vol, tr);
+    // pv = assembly.placeVolume(support_vol, Position(0, 0, support_zstart + support_length / 2));
+  }
+
   // loop over the modules
   for (xml_coll_t mi(x_det, _U(module)); mi; ++mi) {
     xml_comp_t x_mod = mi;
-    xml_comp_t m_env = x_mod.child(_U(frame));
     string     m_nam = x_mod.nameStr();
-    xmleles[m_nam]  = x_mod;
-
-    // triangular volume envelope
-    double frame_thickness = m_env.thickness();
-    double frame_width     = m_env.width();
-    double frame_height    = getAttrOrDefault<double>(m_env, _U(height), 5.0 * mm);
-    double tanth           = frame_height/(frame_width/2.0);
-    double frame_height2   = frame_height-frame_thickness-frame_thickness/tanth;
-    double frame_width2    = 2.0*frame_height2/tanth;
-
-    Trd1 moduleframe_part1(frame_width / 2, 0.001 * mm, m_env.length() / 2,
-                           frame_height / 2);
-    Trd1 moduleframe_part2(frame_width2/2, 0.001 * mm,
-                           m_env.length() / 2 + 0.01 * mm, frame_height2/2);
-    SubtractionSolid moduleframe(moduleframe_part1, moduleframe_part2,Position(0.0,frame_thickness,0.0));
-    Volume v_module(m_nam+"_vol", moduleframe, description.material(m_env.materialStr()));
-    v_module.setVisAttributes(description, m_env.visStr());
-
-    int ncomponents   = 0;
-    int sensor_number = 1;
 
     if (volumes.find(m_nam) != volumes.end()) {
-      printout(ERROR, "SiTrackerBarrel", "Logics error in building modules.");
+      printout(ERROR, "BarrelTrackerWithFrame", string((string("Module with named ") + m_nam + string(" already exists."))).c_str() );
       throw runtime_error("Logics error in building modules.");
     }
+
+    int ncomponents   = 0;
+    int sensor_number = 1;
     double total_thickness = 0;
+
+    // Compute module total thickness from components
     xml_coll_t ci(x_mod, _U(module_component));
     for (ci.reset(), total_thickness = 0.0; ci; ++ci) {
       total_thickness += xml_comp_t(ci).thickness();
     }
-    // module assembly
+    // the module assembly volume
     Assembly m_vol( m_nam );
-    m_vol.placeVolume(v_module, Position(0.0,0.0,frame_height/2+total_thickness/2.0));
     volumes[m_nam] = m_vol;
     m_vol.setVisAttributes(description.visAttributes(x_mod.visStr()));
 
+    // Optional module frame.
+    if(x_mod.hasChild("frame")){
+      xml_comp_t m_frame = x_mod.child(_U(frame));
+      //xmleles[m_nam]  = x_mod;
+      double frame_thickness = m_frame.thickness();
+      double frame_width     = m_frame.width();
+      double frame_height    = getAttrOrDefault<double>(m_frame, _U(height), 5.0 * mm);
+      double tanth           = frame_height/(frame_width/2.0);
+      double frame_height2   = frame_height-frame_thickness-frame_thickness/tanth;
+      double frame_width2    = 2.0*frame_height2/tanth;
+
+      Trd1 moduleframe_part1(frame_width / 2, 0.001 * mm, m_frame.length() / 2,
+                             frame_height / 2);
+      Trd1 moduleframe_part2(frame_width2/2, 0.001 * mm,
+                             m_frame.length() / 2 + 0.01 * mm, frame_height2/2);
+
+      SubtractionSolid moduleframe(moduleframe_part1, moduleframe_part2,Position(0.0,frame_thickness,0.0));
+      Volume v_moduleframe(m_nam+"_vol", moduleframe, description.material(m_frame.materialStr()));
+      v_moduleframe.setVisAttributes(description, m_frame.visStr());
+      m_vol.placeVolume(v_moduleframe, Position(0.0, 0.0, frame_height / 2 + total_thickness / 2.0));
+    }
+
+
+    double thickness_so_far = 0.0;
     double thickness_sum = -total_thickness/2.0;
     for (xml_coll_t ci(x_mod, _U(module_component)); ci; ++ci, ++ncomponents) {
-      xml_comp_t x_comp = ci;
-      xml_comp_t x_pos  = x_comp.position(false);
-      xml_comp_t x_rot  = x_comp.rotation(false);
-      string     c_nam  = _toString(ncomponents, "component%d");
-      Box        c_box(x_comp.width() / 2, x_comp.length() / 2, x_comp.thickness() / 2);
-      Volume     c_vol(c_nam, c_box, description.material(x_comp.materialStr()));
+      xml_comp_t   x_comp = ci;
+      xml_comp_t   x_pos  = x_comp.position(false);
+      xml_comp_t   x_rot  = x_comp.rotation(false);
+      const string c_nam  = _toString(ncomponents, "component%d");
+      Box          c_box(x_comp.width() / 2, x_comp.length() / 2, x_comp.thickness() / 2);
+      Volume       c_vol(c_nam, c_box, description.material(x_comp.materialStr()));
 
+      // Utility variable for the relative z-offset based off the previous components
+      const double zoff = thickness_sum+x_comp.thickness() / 2.0;
       if (x_pos && x_rot) {
-        Position c_pos(x_pos.x(0), x_pos.y(0), x_pos.z(0));
+        Position c_pos(x_pos.x(0), x_pos.y(0), x_pos.z(0) + zoff);
         RotationZYX c_rot(x_rot.z(0), x_rot.y(0), x_rot.x(0));
         pv = m_vol.placeVolume(c_vol, Transform3D(c_rot, c_pos));
       } else if (x_rot) {
-        Position c_pos(0, 0, thickness_sum + x_comp.thickness() / 2.0);
-        pv = m_vol.placeVolume(c_vol, Transform3D(RotationZYX(x_rot.z(0), x_rot.y(0), x_rot.x(0)),c_pos));
+        Position c_pos(0, 0, zoff);
+        pv = m_vol.placeVolume(c_vol, Transform3D(RotationZYX(x_rot.z(0), x_rot.y(0), x_rot.x(0)), c_pos));
       } else if (x_pos) {
-        pv = m_vol.placeVolume(c_vol, Position(x_pos.x(0), x_pos.y(0), x_pos.z(0)));
+        pv = m_vol.placeVolume(c_vol, Position(x_pos.x(0), x_pos.y(0), x_pos.z(0) + zoff));
       } else {
-        pv = m_vol.placeVolume(c_vol, Position(0,0,thickness_sum+x_comp.thickness()/2.0));
+        pv = m_vol.placeVolume(c_vol, Position(0, 0, zoff));
       }
       c_vol.setRegion(description, x_comp.regionStr());
       c_vol.setLimitSet(description, x_comp.limitsStr());
@@ -148,8 +199,36 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
         pv.addPhysVolID("sensor", sensor_number++);
         c_vol.setSensitiveDetector(sens);
         sensitives[m_nam].push_back(pv);
+        module_thicknesses[m_nam] = {thickness_so_far + x_comp.thickness()/2.0, total_thickness-thickness_so_far - x_comp.thickness()/2.0};
+
+        // -------- create a measurement plane for the tracking surface attched to the sensitive volume -----
+        Vector3D u(-1., 0., 0.);
+        Vector3D v(0., -1., 0.);
+        Vector3D n(0., 0., 1.);
+        //    Vector3D o( 0. , 0. , 0. ) ;
+
+        // compute the inner and outer thicknesses that need to be assigned to the tracking surface
+        // depending on wether the support is above or below the sensor
+        double inner_thickness = module_thicknesses[m_nam][0];
+        double outer_thickness = module_thicknesses[m_nam][1];
+
+        SurfaceType type(SurfaceType::Sensitive);
+
+        // if( isStripDetector )
+        //  type.setProperty( SurfaceType::Measurement1D , true ) ;
+
+        VolPlane surf(c_vol, type, inner_thickness, outer_thickness, u, v, n); //,o ) ;
+        volplane_surfaces[m_nam].push_back(surf);
+
+        //--------------------------------------------
       }
       thickness_sum += x_comp.thickness();
+      thickness_so_far += x_comp.thickness();
+      // apply relative offsets in z-position used to stack components side-by-side
+      if (x_pos) {
+        thickness_sum += x_pos.z(0);
+        thickness_so_far += x_pos.z(0);
+      }
     }
   }
 
@@ -178,19 +257,22 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
     double      z_dr       = z_layout.dr();       // Radial displacement parameter, of every other module.
 
     Volume      module_env = volumes[m_nam];
-    DetElement  lay_elt(sdet, _toString(x_layer.id(), "layer%d"), lay_id);
+    DetElement  lay_elt(sdet, lay_nam, lay_id);
     Placements& sensVols = sensitives[m_nam];
 
     // the local coordinate systems of modules in dd4hep and acts differ
     // see http://acts.web.cern.ch/ACTS/latest/doc/group__DD4hepPlugins.html
-    Acts::ActsExtension* layerExtension = new Acts::ActsExtension();
-    layerExtension->addType("sensitive cylinder", "layer");
-    //layerExtension->addValue(0, "r_min", "envelope");
-    //layerExtension->addValue(0, "r_max", "envelope");
-    //layerExtension->addValue(0, "z_min", "envelope");
-    //layerExtension->addValue(0, "z_max", "envelope");
-    //layerExtension->addType("axes", "definitions", "XzY");
-    lay_elt.addExtension<Acts::ActsExtension>(layerExtension);
+    {
+      Acts::ActsExtension* layerExtension = new Acts::ActsExtension();
+      // layer is simple tube so no need to set envelope
+      layerExtension->addType("sensitive cylinder", "layer");
+      // Add the proto layer material
+      for(xml_coll_t lmat(x_layer, _Unicode(layer_material)); lmat; ++lmat) {
+        xml_comp_t x_layer_material = lmat;
+        xmlToProtoSurfaceMaterial(x_layer_material, *layerExtension, "layer_material");
+      }
+      lay_elt.addExtension<Acts::ActsExtension>(layerExtension);
+    }
 
     // Z increment for module placement along Z axis.
     // Adjust for z0 at center of module rather than
@@ -222,11 +304,17 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
           PlacedVolume sens_pv = sensVols[ic];
           DetElement comp_de(mod_elt, std::string("de_") + sens_pv.volume().name(), module);
           comp_de.setPlacement(sens_pv);
-          Acts::ActsExtension* sensorExtension = new Acts::ActsExtension();
-          //sensorExtension->addType("sensor", "detector");
-          comp_de.addExtension<Acts::ActsExtension>(sensorExtension);
+          // ACTS extension
+          {
+            Acts::ActsExtension* sensorExtension = new Acts::ActsExtension();
+            //sensorExtension->addType("sensor", "detector");
+            comp_de.addExtension<Acts::ActsExtension>(sensorExtension);
+          }
           //comp_de.setAttributes(description, sens_pv.volume(), x_layer.regionStr(), x_layer.limitsStr(),
           //                      xml_det_t(xmleles[m_nam]).visStr());
+          //
+
+          volSurfaceList(comp_de)->push_back(volplane_surfaces[m_nam][ic]);
         }
 
         /// Increase counters etc.
@@ -256,14 +344,13 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
   assembly.setVisAttributes(description.invisible());
   pv = description.pickMotherVolume(sdet).placeVolume(assembly);
   pv.addPhysVolID("system", det_id);  // Set the subdetector system ID.
-  pv.addPhysVolID("barrel", 1);       // Flag this as a barrel subdetector.
   sdet.setPlacement(pv);
   return sdet;
 }
 
 //@}
 // clang-format off
-DECLARE_DETELEMENT(BarrelTrackerWithFrame, create_detector)
-DECLARE_DETELEMENT(athena_VertexBarrel, create_detector)
-DECLARE_DETELEMENT(athena_TrackerBarrel, create_detector)
-DECLARE_DETELEMENT(refdet_SiVertexBarrel, create_detector)
+DECLARE_DETELEMENT(BarrelTrackerWithFrame, create_BarrelTrackerWithFrame)
+DECLARE_DETELEMENT(athena_TrackerBarrel,   create_BarrelTrackerWithFrame)
+DECLARE_DETELEMENT(athena_VertexBarrel,    create_BarrelTrackerWithFrame)
+DECLARE_DETELEMENT(athena_TOFBarrel,       create_BarrelTrackerWithFrame)

src/CompositeTracker_geo.cpp

+//==========================================================================
+//  Based off DD4hep_SubDetectorAssembly
+//
+//  This is a simple plugin to allow compositing different detectors
+//  into a single barrel or endcap for ACTS
+//  Note: positive/negative position strings differentiate between
+//        positive/negative endcaps
+//--------------------------------------------------------------------------
+
+#include "DD4hep/DetFactoryHelper.h"
+#include "DD4hep/Printout.h"
+#include "XML/Utilities.h"
+
+#if defined(USE_ACTSDD4HEP)
+#include "ActsDD4hep/ActsExtension.hpp"
+#else
+#include "Acts/Plugins/DD4hep/ActsExtension.hpp"
+#endif
+
+using namespace dd4hep;
+using namespace dd4hep::detail;
+
+static Ref_t create_element(Detector& description, xml_h e, Ref_t)
+{
+  xml_det_t         x_det(e);
+  const std::string det_name = x_det.nameStr();
+  DetElement        sdet(det_name, x_det.id());
+  Volume            vol;
+  Position          pos;
+
+  const bool usePos = x_det.hasChild(_U(position));
+
+  sdet.setType("compound");
+  xml::setDetectorTypeFlag(e, sdet);
+
+  const std::string actsType = getAttrOrDefault(x_det, _Unicode(actsType), "endcap");
+  printout(DEBUG, det_name, "+++ Creating composite tracking detector (type: " + actsType + ")");
+  assert(actsType == "barrel" || actsType == "endcap");
+
+  // ACTS extension
+  {
+    Acts::ActsExtension* detWorldExt = new Acts::ActsExtension();
+    detWorldExt->addType(actsType, "detector");
+    sdet.addExtension<Acts::ActsExtension>(detWorldExt);
+  }
+
+  if (usePos) {
+    pos = Position(x_det.position().x(), x_det.position().y(), x_det.position().z());
+  }
+  vol = Assembly(det_name);
+  vol.setAttributes(description, x_det.regionStr(), x_det.limitsStr(), x_det.visStr());
+
+  Volume       mother = description.pickMotherVolume(sdet);
+  PlacedVolume pv;
+  if (usePos) {
+    pv = mother.placeVolume(vol, pos);
+  } else {
+    pv = mother.placeVolume(vol);
+  }
+  sdet.setPlacement(pv);
+  for (xml_coll_t c(x_det, _U(composite)); c; ++c) {
+    xml_dim_t         component = c;
+    const std::string nam       = component.nameStr();
+    description.declareParent(nam, sdet);
+  }
+  return sdet;
+}
+
+DECLARE_DETELEMENT(athena_CompositeTracker, create_element)

src/CylinderTrackerBarrel_geo.cpp

@@ -2,20 +2,26 @@
 // Specialized generic detector constructor
 //==========================================================================
 
-#include "Acts/Plugins/DD4hep/ActsExtension.hpp"
-#include "Acts/Plugins/DD4hep/ConvertDD4hepMaterial.hpp"
 #include "DD4hep/DetFactoryHelper.h"
 #include "DD4hep/Printout.h"
 
+#if defined(USE_ACTSDD4HEP)
+#include "ActsDD4hep/ActsExtension.hpp"
+#include "ActsDD4hep/ConvertMaterial.hpp"
+#else
+#include "Acts/Plugins/DD4hep/ActsExtension.hpp"
+#include "Acts/Plugins/DD4hep/ConvertDD4hepMaterial.hpp"
+#endif
+
 using namespace std;
 using namespace dd4hep;
 using namespace dd4hep::detail;
 
 /** A barrel tracker with a module that is curved (not flat).
  *
- *
+ * \ingroup tracking
  */
-static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector sens)
+static Ref_t CylinderTrackerBarrel_create_detector(Detector& description, xml_h e, SensitiveDetector sens)
 {
   typedef vector<PlacedVolume> Placements;
   xml_det_t                    x_det = e;
@@ -48,7 +54,9 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
     auto     module_length    = m_env.length();
     auto     module_phi       = getAttrOrDefault(m_env, _Unicode(phi), 90.0);
 
-    Volume m_vol(m_nam, Tube(module_rmin, module_rmin + module_thickness, module_length / 2), air);
+    Volume m_vol(
+        m_nam,
+        Tube(module_rmin, module_rmin + module_thickness, module_length / 2, -module_phi / 2.0, module_phi / 2.0), air);
     int    ncomponents = 0, sensor_number = 1;
     module_assembly.placeVolume(m_vol, Position(-module_rmin, 0, 0));
     mod_volumes[m_nam] = module_assembly;
@@ -118,10 +126,10 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
     Volume     m_env    = mod_volumes[m_nam];
     DetElement lay_elt(sdet, _toString(x_layer.id(), "layer%d"), lay_id);
 
-    //Acts::ActsExtension* layerExtension = new Acts::ActsExtension();
-    //layerExtension->addType("sensitive cylinder", "layer");
-    //// layerExtension->addValue(10. * Acts::UnitConstants::mm, "r", "envelope");
-    //lay_elt.addExtension<Acts::ActsExtension>(layerExtension);
+    ///Acts::ActsExtension* layerExtension = new Acts::ActsExtension();
+    ///layerExtension->addType("sensitive cylinder", "layer");
+    /////// layerExtension->addValue(10. * Acts::UnitConstants::mm, "r", "envelope");
+    ///lay_elt.addExtension<Acts::ActsExtension>(layerExtension);
 
     Placements& sensVols = sensitives[m_nam];
 
@@ -158,7 +166,7 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
           PlacedVolume sens_pv = sensVols[ic];
           DetElement   comp_elt(mod_elt, sens_pv.volume().name(), module);
           comp_elt.setPlacement(sens_pv);
-          //Acts::ActsExtension* moduleExtension = new Acts::ActsExtension("YZX");
+          //Acts::ActsExtension* moduleExtension = new Acts::ActsExtension();
           //comp_elt.addExtension<Acts::ActsExtension>(moduleExtension);
         }
 
@@ -194,6 +202,8 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
 }
 
 // clang-format off
-DECLARE_DETELEMENT(refdet_CylinderTrackerBarrel, create_detector)
-DECLARE_DETELEMENT(refdet_MMTrackerBarrel, create_detector)
-DECLARE_DETELEMENT(refdet_RWellTrackerBarrel, create_detector)
+DECLARE_DETELEMENT(athena_CylinderTrackerBarrel, CylinderTrackerBarrel_create_detector)
+DECLARE_DETELEMENT(athena_MMTrackerBarrel,       CylinderTrackerBarrel_create_detector)
+DECLARE_DETELEMENT(athena_RWellTrackerBarrel,    CylinderTrackerBarrel_create_detector)
+DECLARE_DETELEMENT(athena_CylinderVertexBarrel,  CylinderTrackerBarrel_create_detector)
+

src/DIRC_geo.cpp

+#include "DD4hep/DetFactoryHelper.h"
+#include "DD4hep/OpticalSurfaces.h"
+#include "DD4hep/Printout.h"
+#include "DDRec/DetectorData.h"
+#include "DDRec/Surface.h"
+#include <XML/Helper.h>
+
+//////////////////////////////////
+// Central Barrel DIRC
+//////////////////////////////////
+
+using namespace std;
+using namespace dd4hep;
+
+// Fixed Trap constructor. This function is a workaround of this bug:
+// https://github.com/AIDASoft/DD4hep/issues/850
+// Should be used instead of dd4hep::Trap(pName, pZ, pY, pX, pLTX) constructor
+dd4hep::Trap MakeTrap( const std::string& pName, double pZ, double pY, double pX, double pLTX );
+
+static Ref_t createDetector(Detector& desc, xml_h e, SensitiveDetector sens)
+{
+  xml_det_t xml_det   = e;
+  string    det_name = xml_det.nameStr();
+  int       det_id   = xml_det.id();
+
+  // Main detector xml element
+  xml_dim_t dirc_dim   = xml_det.dimensions();
+  xml_dim_t dirc_pos   = xml_det.position();
+  xml_dim_t dirc_rot   = xml_det.rotation();
+  double    det_rin   = dirc_dim.rmin();
+  double    det_rout  = dirc_dim.rmax();
+  double    SizeZ = dirc_dim.length();
+
+  // DEBUG
+  // double mirror_r1 = x_det.attr<double>(_Unicode(r1));
+
+  // DIRC box:
+  xml_comp_t xml_box_module = xml_det.child(_U(module));
+
+  Material Vacuum = desc.material("Vacuum");
+  Material air = desc.material("AirOptical");
+  Material quartz = desc.material("Quartz");
+  Material epotek = desc.material("Epotek");
+  Material nlak33a = desc.material("Nlak33a");
+  auto& bar_material = quartz;
+  auto mirror_material = desc.material("Aluminum");  // mirror material
+
+  Tube     det_geo(det_rin, det_rout, SizeZ / 2., 0., 360.0 * deg);
+  //Volume   det_volume("DIRC", det_geo, Vacuum);
+  Assembly   det_volume("DIRC");
+  det_volume.setVisAttributes(desc.visAttributes(xml_det.visStr()));
+
+  DetElement   det(det_name, det_id);
+  Volume       mother_vol = desc.pickMotherVolume(det);
+
+  Transform3D  tr(RotationZYX(0, dirc_rot.theta(), 0.0), Position(0.0, 0.0, dirc_pos.z()));
+  PlacedVolume det_plvol = mother_vol.placeVolume(det_volume, tr);
+
+  det_plvol.addPhysVolID("system", det_id);
+  det.setPlacement(det_plvol);
+
+  // Parts Dimentions
+
+  int fLensId = 6; // focusing system
+                   // 0    no lens
+                   // 1    spherical lens
+                   // 3    3-layer spherical lens
+                   // 6    3-layer cylindrical lens
+                   // 10   ideal lens (thickness = 0, ideal focusing)
+  
+  int fGeomType = 0;  // Full DIRC - 0, 1 only one plate
+  int fRunType = 0;   // 0, 10 - simulation, 1, 5 - lookup table, 2,3,4 - reconstruction
+
+
+  double fPrizm[4];
+  fPrizm[0] = 360 * mm;
+  fPrizm[1] = 300 * mm;
+  fPrizm[3] = 50 * mm;
+  fPrizm[2] = fPrizm[3] + 300 * tan(32 * deg) * mm;
+  double fBarsGap = 0.15 * mm;
+
+  std::cout << "DIRC: fPrizm[2] " << fPrizm[2] << std::endl;
+
+  double fdTilt = 80 * deg;
+  double fPrizmT[6];
+  fPrizmT[0] = 390 * mm;
+  fPrizmT[1] = (400 - 290 * cos(fdTilt)) * mm; //
+  fPrizmT[2] = 290 * sin(fdTilt) * mm;       // hight
+  fPrizmT[3] = 50 * mm;                      // face
+  fPrizmT[4] = 290 * mm;
+  fPrizmT[5] = 290 * cos(fdTilt)* mm;
+
+  double fMirror[3];
+  fMirror[0] = 20 * mm;
+  fMirror[1] = fPrizm[0];
+  fMirror[2] = 1 * mm;
+  //  fPrizm[0] = 170; fPrizm[1] = 300; fPrizm[2] = 50+300*tan(45*deg); fPrizm[3] = 50;
+
+//  double fBar[3];
+//  fBar[0] = 17 * mm;
+//  fBar[1] = (fPrizm[0] - (fNBar - 1) * fBarsGap) / fNBar;
+//  fBar[2] = 1050 * mm; // 4200; //4200
+
+  double fMcpTotal[3];
+  double fMcpActive[3];
+  fMcpTotal[0] = fMcpTotal[1] = 53 + 4;
+  fMcpTotal[2] = 1*mm;
+  fMcpActive[0] = fMcpActive[1] = 53;
+  fMcpActive[2] = 1*mm;
+
+  double fLens[4];
+  fLens[0] = fLens[1] = 40 * mm;
+  fLens[2] = 10 * mm;
+  double fRadius = (det_rin + det_rout)/2;
+
+  double fBoxWidth = fPrizm[0];
+
+  double fFd[3];
+  fFd[0] = fBoxWidth;
+  fFd[1] = fPrizm[2];
+  fFd[2] = 1 * mm;
+
+  fLens[0] = fPrizm[3];
+  fLens[1] = fPrizm[0];
+  fLens[2] = 12 * mm;
+
+  // Getting box XML
+  const int fNBoxes = xml_box_module.repeat();
+  const double box_width = xml_box_module.width();
+  const double box_height = xml_box_module.height();
+  const double box_length = xml_box_module.length() + 550*mm;
+
+
+  // The DIRC
+  Assembly dirc_module("DIRCModule");
+
+  //Volume lDirc("lDirc", gDirc, air);
+  dirc_module.setVisAttributes(desc.visAttributes(xml_box_module.visStr()));
+
+
+  // FD... whatever F and D is
+  xml_comp_t xml_fd = xml_box_module.child(_Unicode(fd));
+  Box gFd("gFd", xml_fd.height()/2, xml_fd.width()/2, xml_fd.thickness()/2);
+  Volume lFd ("lFd", gFd, desc.material(xml_fd.materialStr()));
+  lFd.setVisAttributes(desc.visAttributes(xml_fd.visStr()));
+  //lFd.setSensitiveDetector(sens);
+
+
+  // The Bar
+  xml_comp_t xml_bar = xml_box_module.child(_Unicode(bar));
+  double bar_height = xml_bar.height();
+  double bar_width = xml_bar.width();
+  double bar_length = xml_bar.length();
+  Box gBar("gBar", bar_height/2, bar_width/2, bar_length/2);
+  Volume lBar("lBar", gBar, desc.material(xml_bar.materialStr()));
+  lBar.setVisAttributes(desc.visAttributes(xml_bar.visStr()));
+
+  // Glue
+  xml_comp_t xml_glue = xml_box_module.child(_Unicode(glue));
+  double glue_thickness = xml_glue.thickness();  // 0.05 * mm;
+  Box gGlue("gGlue", bar_height/2, bar_width/2, glue_thickness/2);
+  Volume lGlue("lGlue", gGlue, desc.material(xml_glue.materialStr()));
+  lGlue.setVisAttributes(desc.visAttributes(xml_glue.visStr()));
+
+  sens.setType("tracker");
+  lBar.setSensitiveDetector(sens);
+
+  int bars_repeat_z = 4;    // TODO parametrize!
+  double bar_assm_length = (bar_length + glue_thickness) * bars_repeat_z;
+  int fNBar = xml_bar.repeat();
+  double bar_gap = xml_bar.gap();
+  for (int y_index = 0; y_index < fNBar; y_index++) {
+    double shift_y = y_index * (bar_width + bar_gap) - 0.5 * box_width + 0.5 * bar_width;
+    for (int z_index = 0; z_index < bars_repeat_z; z_index++) {
+      double z = -0.5 * bar_assm_length + 0.5 * bar_length + (bar_length + glue_thickness) * z_index;
+      auto placed_bar = dirc_module.placeVolume(lBar, Position(0, shift_y, z));
+      dirc_module.placeVolume(lGlue, Position(0, shift_y, z + 0.5 * (bar_length + glue_thickness)));
+      placed_bar.addPhysVolID("section", z_index);
+      placed_bar.addPhysVolID("bar", y_index);
+    }
+  }
+
+
+  // The Mirror
+  xml_comp_t xml_mirror = xml_box_module.child(_Unicode(mirror));
+  Box gMirror("gMirror", xml_mirror.height()/2, xml_mirror.width()/2, xml_mirror.thickness()/2);
+  Volume lMirror("lMirror", gMirror, desc.material(xml_mirror.materialStr()));
+  dirc_module.placeVolume(lMirror, Position(0, 0, -0.5 * (bar_assm_length - xml_mirror.thickness())));
+  lMirror.setVisAttributes(desc.visAttributes(xml_mirror.visStr()));
+
+  // The mirror optical surface
+  OpticalSurfaceManager surfMgr = desc.surfaceManager();
+  auto surf = surfMgr.opticalSurface("MirrorOpticalSurface");
+  SkinSurface skin(desc, det, Form("dirc_mirror_optical_surface"), surf, lMirror);
+  skin.isValid();
+
+  // LENS
+  // Lens volumes
+  Volume lLens1;
+  Volume lLens2;
+  Volume lLens3;
+
+  double lensMinThikness = 2.0 * mm;
+  double layer12 = lensMinThikness * 2;
+
+  // r1 = (r1==0)? 27.45: r1;
+  // r2 = (r2==0)? 20.02: r2;
+
+  double r1 = 33 * mm;
+  double r2 = 24 * mm;
+  double shight = 25 * mm;
+
+  Position zTrans1(0, 0, -r1 - fLens[2] / 2. + r1 - sqrt(r1 * r1 - shight / 2. * shight / 2.) + lensMinThikness);
+  Position zTrans2(0, 0, -r2 - fLens[2] / 2. + r2 - sqrt(r2 * r2 - shight / 2. * shight / 2.) + layer12);
+
+  Box gfbox("fbox", 0.5 * fLens[0], 0.5 * fLens[1], 0.5 * fLens[2]);
+  Box gcbox("cbox", 0.5 * fLens[0], 0.5 * fLens[1] + 1*mm, 0.5 * fLens[2]);
+
+//  Volume gfbox_volume("gfbox_volume", gfbox, bar_material);
+//  lDirc.placeVolume(gfbox_volume, Position(0, 0, 0));
+//
+//  Volume gcbox_volume("gcbox_volume", gcbox, bar_material);
+//  lDirc.placeVolume(gcbox_volume, Position(0, 0, 50));
+
+
+  Position tTrans1(0.5 * (fLens[0] + shight), 0, -fLens[2] + layer12);
+  Position tTrans0(-0.5 * (fLens[0] + shight), 0, -fLens[2] + layer12);
+  SubtractionSolid tubox("tubox", gfbox, gcbox, tTrans1);
+  SubtractionSolid gubox("gubox", tubox, gcbox, tTrans0);
+
+//  Volume tubox_volume("tubox_volume", tubox, bar_material);
+//  lDirc.placeVolume(tubox_volume, Position(0, 0, 100));
+//
+//  Volume gubox_volume("gubox_volume", gubox, bar_material);
+//  lDirc.placeVolume(gubox_volume, Position(0, 0, 150));
+
+  Tube gcylinder1("Cylinder1", 0, r1, 0.5 * fLens[1], 0 * deg, 360 * deg);
+  Tube gcylinder2("Cylinder2", 0, r2, 0.5 * fLens[1] - 0.5*mm, 0 * deg, 360 * deg);
+  Tube gcylinder1c("Cylinder1c", 0, r1, 0.5 * fLens[1] + 0.5*mm, 0 * deg, 360 * deg);
+  Tube gcylinder2c("Cylinder2c", 0, r2, 0.5 * fLens[1] + 0.5*mm, 0 * deg, 360 * deg);
+  RotationX xRot(-M_PI / 2.);
+
+  IntersectionSolid gLens1("Lens1", gubox, gcylinder1, Transform3D(xRot, zTrans1));
+  SubtractionSolid gLenst("temp", gubox, gcylinder1c, Transform3D(xRot, zTrans1));
+
+//  Volume gLens1_volume("gLens1_volume", gLens1, bar_material);
+//  lDirc.placeVolume(gLens1_volume, Position(0, 0, 200));
+//
+//  Volume gLenst_volume("gLenst_volume", gLenst, bar_material);
+//  lDirc.placeVolume(gLenst_volume, Position(0, 0, 250));
+
+  IntersectionSolid gLens2("Lens2", gLenst, gcylinder2, Transform3D(xRot, zTrans2));
+  SubtractionSolid gLens3("Lens3", gLenst, gcylinder2c, Transform3D(xRot, zTrans2));
+
+  lLens1 = Volume("lLens1", gLens1, bar_material);
+  lLens2 = Volume("lLens2", gLens2, nlak33a);
+  lLens3 = Volume("lLens3", gLens3, bar_material);
+
+  lLens1.setVisAttributes(desc.visAttributes("DIRCLens1"));
+  lLens2.setVisAttributes(desc.visAttributes("DIRCLens2"));
+  lLens3.setVisAttributes(desc.visAttributes("DIRCLens3"));
+
+
+  double shifth = 0.5 * (bar_assm_length + fLens[2]);
+  // fmt::print("LENS HERE shifth={}\n", shifth);
+
+  lLens1.setVisAttributes(desc.visAttributes("AnlTeal"));
+  dirc_module.placeVolume(lLens1, Position(0, 0, shifth));
+  dirc_module.placeVolume(lLens2, Position(0, 0, shifth));
+  dirc_module.placeVolume(lLens3, Position(0, 0, shifth));
+  
+  // The Prizm
+  Trap gPrizm = MakeTrap("gPrizm", fPrizm[0], fPrizm[1], fPrizm[2], fPrizm[3]);
+  Volume lPrizm("lPrizm", gPrizm, bar_material);
+  lPrizm.setVisAttributes(desc.visAttributes("DIRCPrism"));
+
+  //G4RotationMatrix *fdRot = new G4RotationMatrix();
+  //G4RotationMatrix *fdrot = new G4RotationMatrix();
+  double evshiftz = 0.5 * bar_assm_length + fPrizm[1] + fMcpActive[2] / 2. + fLens[2];
+  double evshiftx = -3*mm;
+
+  double prism_shift_x = (fPrizm[2] + fPrizm[3]) / 4. - 0.5 * fPrizm[3] + 1.5*mm;
+  double prism_shift_z = 0.5 * (bar_assm_length + fPrizm[1]) + fLens[2];
+
+  Position  fPrismShift(prism_shift_x, 0, prism_shift_z);
+  dirc_module.placeVolume(lPrizm, Transform3D(xRot, fPrismShift));
+  dirc_module.placeVolume(lFd, Position(0.5 * fFd[1] - 0.5 * fPrizm[3] - evshiftx, 0, evshiftz));
+
+  double dphi = 2 * M_PI / (double)fNBoxes;
+  for (int i = 0; i < fNBoxes; i++) {
+    double phi = dphi * i;
+    double dx = -fRadius * cos(phi);
+    double dy = -fRadius * sin(phi);
+
+    //G4RotationMatrix *tRot = new G4RotationMatrix();
+
+    Transform3D  tr(RotationZ(phi+M_PI), Position(dx, dy, 0));
+    PlacedVolume box_placement = det_volume.placeVolume(dirc_module, tr);
+    box_placement.addPhysVolID("module", i);
+
+
+    // fmt::print("placing dircbox # {} -tphi={:.0f} dx={:.0f}, dy={:.0f}\n", i, phi/deg, dx/cm, dy/cm);
+
+    //new G4PVPlacement(tRot, G4ThreeVector(dx, dy, 0), lDirc, "wDirc", lExpHall, false, i);
+  }
+
+    
+  //////////////////
+  // DIRC Bars
+  //////////////////
+
+  // double bar_radius   = 83.65 * cm;
+  // double bar_length   = SizeZ;
+  // double bar_width    = 42. * cm;
+  // double bar_thicknes = 1.7 * cm;
+  // int    bar_count    = 2 * M_PI * bar_radius / bar_width;
+  // double bar_dphi     = 2 * 3.1415926 / bar_count;
+  // Material bar_material = desc.material("Quartz");
+
+  // Box    bar_geo(bar_thicknes / 2., bar_width / 2., bar_length / 2.);
+  // Volume bar_volume("cb_DIRC_bars_Logix", bar_geo, bar_material);
+  // bar_volume.setVisAttributes(desc.visAttributes(xml_det.visStr()));
+  // sens.setType("tracker");
+  // bar_volume.setSensitiveDetector(sens);
+
+  // for (int ia = 0; ia < bar_count; ia++) {
+  //   double phi = (ia * (bar_dphi));
+  //   double x   = -bar_radius * cos(phi);
+  //   double y   = -bar_radius * sin(phi);
+
+  //   Transform3D  tr(RotationZ(bar_dphi * ia), Position(x, y, 0));
+  //   PlacedVolume barPV = det_volume.placeVolume(bar_volume, tr);
+  //   barPV.addPhysVolID("module", ia);
+  // }
+  return det;
+}
+
+DECLARE_DETELEMENT(cb_DIRC, createDetector)
+
+dd4hep::Trap MakeTrap( const std::string& pName, double pZ, double pY, double pX, double pLTX )
+{
+    // Fixed Trap constructor. This function is a workaround of this bug:
+    // https://github.com/AIDASoft/DD4hep/issues/850
+    // Should be used instead of dd4hep::Trap(pName, pZ, pY, pX, pLTX) constructor
+
+    double fDz  = 0.5*pZ;
+    double fTthetaCphi = 0;
+    double fTthetaSphi = 0;
+    double fDy1 = 0.5*pY;
+    double fDx1 = 0.5*pX;
+    double fDx2 = 0.5*pLTX;
+    double fTalpha1 = 0.5*(pLTX - pX)/pY;
+    double fDy2 = fDy1;
+    double fDx3 = fDx1;
+    double fDx4 = fDx2;
+    double fTalpha2 = fTalpha1;
+
+    return Trap(pName, fDz,  fTthetaCphi,  fTthetaSphi,
+                fDy1,  fDx1,  fDx2,  fTalpha1,
+                fDy2,  fDx3,  fDx4,  fTalpha2);
+}
+

src/DRICH_geo.cpp

+//==========================================================================
+//  dRICh: Dual Ring Imaging Cherenkov Detector
+//--------------------------------------------------------------------------
+//
+// Author: Christopher Dilks (Duke University)
+//
+// - Design Adapted from Standalone Fun4all and GEMC implementations
+//   [ Evaristo Cisbani, Cristiano Fanelli, Alessio Del Dotto, et al. ]
+//
+//==========================================================================
+
+#include <XML/Helper.h>
+#include "TMath.h"
+#include "TString.h"
+#include "GeometryHelpers.h"
+#include "Math/Point2D.h"
+#include "DDRec/Surface.h"
+#include "DDRec/DetectorData.h"
+#include "DD4hep/OpticalSurfaces.h"
+#include "DD4hep/DetFactoryHelper.h"
+#include "DD4hep/Printout.h"
+
+using namespace dd4hep;
+using namespace dd4hep::rec;
+
+// create the detector
+static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetector sens) {
+  xml::DetElement detElem = handle;
+  std::string detName = detElem.nameStr();
+  int detID = detElem.id();
+
+  DetElement det(detName, detID);
+  xml::Component dims = detElem.dimensions();
+  OpticalSurfaceManager surfMgr = desc.surfaceManager();
+
+  // attributes -----------------------------------------------------------
+  // - vessel
+  double  vesselZmin       =  dims.attr<double>(_Unicode(zmin));
+  double  vesselLength     =  dims.attr<double>(_Unicode(length));
+  double  vesselRmin0      =  dims.attr<double>(_Unicode(rmin0));
+  double  vesselRmin1      =  dims.attr<double>(_Unicode(rmin1));
+  double  vesselRmax0      =  dims.attr<double>(_Unicode(rmax0));
+  double  vesselRmax1      =  dims.attr<double>(_Unicode(rmax1));
+  double  vesselRmax2      =  dims.attr<double>(_Unicode(rmax2));
+  double  snoutLength      =  dims.attr<double>(_Unicode(snout_length));
+  int     nSectors         =  dims.attr<int>(_Unicode(nsectors));
+  double  wallThickness    =  dims.attr<double>(_Unicode(wall_thickness));
+  double  windowThickness  =  dims.attr<double>(_Unicode(window_thickness));
+  auto    vesselMat        =  desc.material(detElem.attr<std::string>(_Unicode(material)));
+  auto    gasvolMat        =  desc.material(detElem.attr<std::string>(_Unicode(gas)));
+  auto    vesselVis        =  desc.visAttributes(detElem.attr<std::string>(_Unicode(vis_vessel)));
+  auto    gasvolVis        =  desc.visAttributes(detElem.attr<std::string>(_Unicode(vis_gas)));
+  // - radiator (applies to aerogel and filter)
+  auto    radiatorElem        =  detElem.child(_Unicode(radiator));
+  double  radiatorRmin        =  radiatorElem.attr<double>(_Unicode(rmin));
+  double  radiatorRmax        =  radiatorElem.attr<double>(_Unicode(rmax));
+  double  radiatorPhiw        =  radiatorElem.attr<double>(_Unicode(phiw));
+  double  radiatorPitch       =  radiatorElem.attr<double>(_Unicode(pitch));
+  double  radiatorFrontplane  =  radiatorElem.attr<double>(_Unicode(frontplane));
+  // - aerogel
+  auto    aerogelElem       =  radiatorElem.child(_Unicode(aerogel));
+  auto    aerogelMat        =  desc.material(aerogelElem.attr<std::string>(_Unicode(material)));
+  auto    aerogelVis        =  desc.visAttributes(aerogelElem.attr<std::string>(_Unicode(vis)));
+  double  aerogelThickness  =  aerogelElem.attr<double>(_Unicode(thickness));
+  // - filter
+  auto    filterElem       =  radiatorElem.child(_Unicode(filter));
+  auto    filterMat        =  desc.material(filterElem.attr<std::string>(_Unicode(material)));
+  auto    filterVis        =  desc.visAttributes(filterElem.attr<std::string>(_Unicode(vis)));
+  double  filterThickness  =  filterElem.attr<double>(_Unicode(thickness));
+  // - mirror
+  auto    mirrorElem       =  detElem.child(_Unicode(mirror));
+  auto    mirrorMat        =  desc.material(mirrorElem.attr<std::string>(_Unicode(material)));
+  auto    mirrorVis        =  desc.visAttributes(mirrorElem.attr<std::string>(_Unicode(vis)));
+  auto    mirrorSurf       =  surfMgr.opticalSurface(mirrorElem.attr<std::string>(_Unicode(surface)));
+  double  mirrorBackplane  =  mirrorElem.attr<double>(_Unicode(backplane));
+  double  mirrorThickness  =  mirrorElem.attr<double>(_Unicode(thickness));
+  double  mirrorRmin       =  mirrorElem.attr<double>(_Unicode(rmin));
+  double  mirrorRmax       =  mirrorElem.attr<double>(_Unicode(rmax));
+  double  mirrorPhiw       =  mirrorElem.attr<double>(_Unicode(phiw));
+  double  focusTuneZ       =  mirrorElem.attr<double>(_Unicode(focus_tune_z));
+  double  focusTuneX       =  mirrorElem.attr<double>(_Unicode(focus_tune_x));
+  // - sensor module
+  auto    sensorElem       =  detElem.child(_Unicode(sensors)).child(_Unicode(module));
+  auto    sensorMat        =  desc.material(sensorElem.attr<std::string>(_Unicode(material)));
+  auto    sensorVis        =  desc.visAttributes(sensorElem.attr<std::string>(_Unicode(vis)));
+  auto    sensorSurf       =  surfMgr.opticalSurface(sensorElem.attr<std::string>(_Unicode(surface)));
+  double  sensorSide       =  sensorElem.attr<double>(_Unicode(side));
+  double  sensorGap        =  sensorElem.attr<double>(_Unicode(gap));
+  double  sensorThickness  =  sensorElem.attr<double>(_Unicode(thickness));
+  // - sensor sphere
+  auto    sensorSphElem     =  detElem.child(_Unicode(sensors)).child(_Unicode(sphere));
+  double  sensorSphRadius   =  sensorSphElem.attr<double>(_Unicode(radius));
+  double  sensorSphCenterX  =  sensorSphElem.attr<double>(_Unicode(centerx));
+  double  sensorSphCenterZ  =  sensorSphElem.attr<double>(_Unicode(centerz));
+  // - sensor sphere patch cuts
+  auto    sensorSphPatchElem  =  detElem.child(_Unicode(sensors)).child(_Unicode(sphericalpatch));
+  double  sensorSphPatchPhiw  =  sensorSphPatchElem.attr<double>(_Unicode(phiw));
+  double  sensorSphPatchRmin  =  sensorSphPatchElem.attr<double>(_Unicode(rmin));
+  double  sensorSphPatchRmax  =  sensorSphPatchElem.attr<double>(_Unicode(rmax));
+  double  sensorSphPatchZmin  =  sensorSphPatchElem.attr<double>(_Unicode(zmin));
+  // - debugging switches
+  int   debug_optics_mode  =  detElem.attr<int>(_Unicode(debug_optics));
+  bool  debug_mirror       =  mirrorElem.attr<bool>(_Unicode(debug));
+  bool  debug_sensors      =  sensorSphElem.attr<bool>(_Unicode(debug));
+
+  // if debugging optics, override some settings
+  bool debug_optics = debug_optics_mode > 0;
+  if(debug_optics) {
+    printout(WARNING,"DRich_geo","DEBUGGING DRICH OPTICS");
+    switch(debug_optics_mode) {
+      case 1: vesselMat = aerogelMat = filterMat = sensorMat = gasvolMat = desc.material("VacuumOptical"); break;
+      case 2: vesselMat = aerogelMat = filterMat = sensorMat = desc.material("VacuumOptical"); break;
+      default: printout(FATAL,"DRich_geo","UNKNOWN debug_optics_mode"); return det;
+    };
+    aerogelVis = sensorVis = mirrorVis;
+    gasvolVis = vesselVis = desc.invisible();
+  };
+
+
+  // BUILD VESSEL ====================================================================
+  /* - `vessel`: aluminum enclosure, the mother volume of the dRICh
+   * - `gasvol`: gas volume, which fills `vessel`; all other volumes defined below
+   *   are children of `gasvol`
+   * - the dRICh vessel geometry has two regions: the snout refers to the conic region
+   *   in the front, housing the aerogel, while the tank refers to the cylindrical
+   *   region, housing the rest of the detector components
+   */
+
+  // derived attributes
+  double tankLength = vesselLength - snoutLength;
+  double vesselZmax = vesselZmin + vesselLength;
+
+  // snout solids
+  double boreDelta = vesselRmin1 - vesselRmin0;
+  double snoutDelta = vesselRmax1 - vesselRmax0;
+  Cone vesselSnout(
+      snoutLength/2.0,
+      vesselRmin0,
+      vesselRmax0,
+      vesselRmin0 + boreDelta * snoutLength / vesselLength,
+      vesselRmax1
+      );
+  Cone gasvolSnout(
+      /* note: `gasvolSnout` extends a bit into the tank, so it touches `gasvolTank`
+       * - the extension distance is equal to the tank `windowThickness`, so the
+       *   length of `gasvolSnout` == length of `vesselSnout`
+       * - the extension backplane radius is calculated using similar triangles
+       */
+      snoutLength/2.0,
+      vesselRmin0 + wallThickness,
+      vesselRmax0 - wallThickness,
+      vesselRmin0 + boreDelta * (snoutLength-windowThickness) / vesselLength + wallThickness,
+      vesselRmax1 - wallThickness + windowThickness * (vesselRmax1 - vesselRmax0) / snoutLength
+      );
+
+  // tank solids
+  Cone vesselTank(
+      tankLength/2.0,
+      vesselSnout.rMin2(),
+      vesselRmax2,
+      vesselRmin1,
+      vesselRmax2
+      );
+  Cone gasvolTank(
+      tankLength/2.0 - windowThickness,
+      gasvolSnout.rMin2(),
+      vesselRmax2 - wallThickness,
+      vesselRmin1 + wallThickness,
+      vesselRmax2 - wallThickness
+      );
+
+  // snout + tank solids
+  UnionSolid vesselUnion(
+      vesselTank,
+      vesselSnout,
+      Position(0., 0., -vesselLength/2.)
+      );
+  UnionSolid gasvolUnion(
+      gasvolTank,
+      gasvolSnout,
+      Position(0., 0., -vesselLength/2. + windowThickness)
+      );
+
+  //  extra solids for `debug_optics` only
+  Box vesselBox(1001,1001,1001);
+  Box gasvolBox(1000,1000,1000);
+
+  // choose vessel and gasvol solids (depending on `debug_optics_mode` (0=disabled))
+  Solid vesselSolid, gasvolSolid;
+  switch(debug_optics_mode) {
+    case 0:  vesselSolid=vesselUnion;  gasvolSolid=gasvolUnion;  break; // `!debug_optics`
+    case 1:  vesselSolid=vesselBox;    gasvolSolid=gasvolBox;    break;
+    case 2:  vesselSolid=vesselBox;    gasvolSolid=gasvolUnion;  break;
+  };
+
+  // volumes
+  Volume vesselVol(detName, vesselSolid, vesselMat);
+  Volume gasvolVol(detName+"_gas", gasvolSolid, gasvolMat);
+  vesselVol.setVisAttributes(vesselVis);
+  gasvolVol.setVisAttributes(gasvolVis);
+
+  // reference positions
+  // - the vessel is created such that the center of the cylindrical tank volume
+  //   coincides with the origin; this is called the "origin position" of the vessel
+  // - when the vessel (and its children volumes) is placed, it is translated in
+  //   the z-direction to be in the proper ATHENA-integration location
+  // - these reference positions are for the frontplane and backplane of the vessel,
+  //   with respect to the vessel origin position
+  auto originFront = Position(0., 0., -tankLength/2.0 - snoutLength );
+  auto originBack =  Position(0., 0., tankLength/2.0 );
+
+  // initialize sensor centroids (used for mirror parameterization below); this is
+  // the average (x,y,z) of the placed sensors, w.r.t. originFront
+  double sensorCentroidX = 0;
+  double sensorCentroidZ = 0;
+  int sensorCount = 0;
+
+
+  // sensitive detector type
+  sens.setType("tracker");
+
+
+  // BUILD RADIATOR ====================================================================
+
+  // attributes
+  double airGap = 0.01*mm; // air gap between aerogel and filter (FIXME? actually it's currently a gas gap)
+
+  // solid and volume: create aerogel and filter
+  Cone aerogelSolid(
+      aerogelThickness/2,
+      radiatorRmin,
+      radiatorRmax,
+      radiatorRmin + boreDelta  * aerogelThickness / vesselLength,
+      radiatorRmax + snoutDelta * aerogelThickness / snoutLength
+      );
+  Cone filterSolid(
+      filterThickness/2,
+      radiatorRmin + boreDelta  * (aerogelThickness + airGap) / vesselLength,
+      radiatorRmax + snoutDelta * (aerogelThickness + airGap) / snoutLength,
+      radiatorRmin + boreDelta  * (aerogelThickness + airGap + filterThickness) / vesselLength,
+      radiatorRmax + snoutDelta * (aerogelThickness + airGap + filterThickness) / snoutLength
+      );
+
+  Volume aerogelVol( detName+"_aerogel", aerogelSolid, aerogelMat );
+  Volume filterVol(  detName+"_filter",  filterSolid,  filterMat );
+  aerogelVol.setVisAttributes(aerogelVis);
+  filterVol.setVisAttributes(filterVis);
+
+  // aerogel placement and surface properties
+  // TODO [low-priority]: define skin properties for aerogel and filter
+  auto radiatorPos = Position(0., 0., radiatorFrontplane) + originFront;
+  auto aerogelPV = gasvolVol.placeVolume(aerogelVol,
+        Translation3D(radiatorPos.x(), radiatorPos.y(), radiatorPos.z()) // re-center to originFront
+      * RotationY(radiatorPitch) // change polar angle to specified pitch // FIXME: probably broken (not yet in use anyway)
+      );
+  DetElement aerogelDE(det, "aerogel_de", 0);
+  aerogelDE.setPlacement(aerogelPV);
+  //SkinSurface aerogelSkin(desc, aerogelDE, "mirror_optical_surface", aerogelSurf, aerogelVol);
+  //aerogelSkin.isValid();
+
+  // filter placement and surface properties
+  if(!debug_optics) {
+    auto filterPV = gasvolVol.placeVolume(filterVol,
+          Translation3D(0., 0., airGap) // add an air gap
+        * Translation3D(radiatorPos.x(), radiatorPos.y(), radiatorPos.z()) // re-center to originFront
+        * RotationY(radiatorPitch) // change polar angle
+        * Translation3D(0., 0., (aerogelThickness+filterThickness)/2.) // move to aerogel backplane
+        );
+    DetElement filterDE(det, "filter_de", 0);
+    filterDE.setPlacement(filterPV);
+    //SkinSurface filterSkin(desc, filterDE, "mirror_optical_surface", filterSurf, filterVol);
+    //filterSkin.isValid();
+  };
+
+
+  // SECTOR LOOP //////////////////////////////////
+  for(int isec=0; isec<nSectors; isec++) {
+
+    // debugging filters, limiting the number of sectors
+    if( (debug_mirror||debug_sensors||debug_optics) && isec!=0) continue;
+
+    // sector rotation about z axis
+    double sectorRotation = isec * 360/nSectors * degree;
+    std::string secName = "sec" + std::to_string(isec);
+
+
+
+    // BUILD SENSORS ====================================================================
+
+    // if debugging sphere properties, restrict number of sensors drawn
+    if(debug_sensors) { sensorSide = 2*M_PI*sensorSphRadius / 64; };
+
+    // solid and volume: single sensor module
+    Box sensorSolid(sensorSide/2., sensorSide/2., sensorThickness/2.);
+    Volume sensorVol(detName+"_sensor_"+secName, sensorSolid, sensorMat);
+    sensorVol.setVisAttributes(sensorVis);
+
+    auto sensorSphPos = Position(sensorSphCenterX, 0., sensorSphCenterZ) + originFront;
+
+    // sensitivity
+    if(!debug_optics) sensorVol.setSensitiveDetector(sens);
+
+    // SENSOR MODULE LOOP ------------------------
+    /* ALGORITHM: generate sphere of positions
+     * - NOTE: there are two coordinate systems here:
+     *   - "global" the main ATHENA coordinate system
+     *   - "generator" (vars end in `Gen`) is a local coordinate system for
+     *     generating points on a sphere; it is related to the global system by
+     *     a rotation; we do this so the "patch" (subset of generated
+     *     positions) of sensors we choose to build is near the equator, where
+     *     point distribution is more uniform
+     * - PROCEDURE: loop over `thetaGen`, with subloop over `phiGen`, each divided evenly
+     *   - the number of points to generate depends how many sensors (+`sensorGap`)
+     *     can fit within each ring of constant `thetaGen` or `phiGen`
+     *   - we divide the relevant circumference by the sensor
+     *     size(+`sensorGap`), and this number is allowed to be a fraction,
+     *     because likely we don't care about generating a full sphere and
+     *     don't mind a "seam" at the overlap point
+     *   - if we pick a patch of the sphere near the equator, and not near
+     *     the poles or seam, the sensor distribution will appear uniform
+     */
+
+    // initialize module number for this sector
+    int imod=0;
+
+    // thetaGen loop: iterate less than "0.5 circumference / sensor size" times
+    double nTheta = M_PI*sensorSphRadius / (sensorSide+sensorGap);
+    for(int t=0; t<(int)(nTheta+0.5); t++) {
+      double thetaGen = t/((double)nTheta) * M_PI;
+
+      // phiGen loop: iterate less than "circumference at this latitude / sensor size" times
+      double nPhi = 2*M_PI * sensorSphRadius * std::sin(thetaGen) / (sensorSide+sensorGap);
+      for(int p=0; p<(int)(nPhi+0.5); p++) {
+        double phiGen = p/((double)nPhi) * 2*M_PI - M_PI; // shift to [-pi,pi]
+
+        // determine global phi and theta
+        // - convert {radius,thetaGen,phiGen} -> {xGen,yGen,zGen}
+        double xGen = sensorSphRadius * std::sin(thetaGen) * std::cos(phiGen);
+        double yGen = sensorSphRadius * std::sin(thetaGen) * std::sin(phiGen);
+        double zGen = sensorSphRadius * std::cos(thetaGen);
+        // - convert {xGen,yGen,zGen} -> global {x,y,z} via rotation
+        double x = zGen;
+        double y = xGen;
+        double z = yGen;
+        // - convert global {x,y,z} -> global {phi,theta}
+        double phi = std::atan2(y,x);
+        double theta = std::acos(z/sensorSphRadius);
+
+        // shift global coordinates so we can apply spherical patch cuts
+        double zCheck = z + sensorSphCenterZ;
+        double xCheck = x + sensorSphCenterX;
+        double yCheck = y;
+        double rCheck = std::hypot(xCheck,yCheck);
+        double phiCheck = std::atan2(yCheck,xCheck);
+
+        // patch cut
+        bool patchCut =
+          std::fabs(phiCheck) < sensorSphPatchPhiw
+          && zCheck > sensorSphPatchZmin
+          && rCheck > sensorSphPatchRmin
+          && rCheck < sensorSphPatchRmax;
+        if(debug_sensors) patchCut = std::fabs(phiCheck) < sensorSphPatchPhiw;
+        if(patchCut) {
+
+          // append sensor position to centroid calculation
+          if(isec==0) {
+            sensorCentroidX += xCheck;
+            sensorCentroidZ += zCheck;
+            sensorCount++;
+          };
+
+          // placement (note: transformations are in reverse order)
+          // - transformations operate on global coordinates; the corresponding
+          //   generator coordinates are provided in the comments
+          auto sensorPV = gasvolVol.placeVolume(sensorVol,
+                RotationZ(sectorRotation) // rotate about beam axis to sector
+              * Translation3D(sensorSphPos.x(), sensorSphPos.y(), sensorSphPos.z()) // move sphere to reference position
+              * RotationX(phiGen) // rotate about `zGen`
+              * RotationZ(thetaGen) // rotate about `yGen`
+              * Translation3D(sensorSphRadius, 0., 0.) // push radially to spherical surface
+              * RotationY(M_PI/2) // rotate sensor to be compatible with generator coords
+              * RotationZ(-M_PI/2) // correction for readout segmentation mapping
+              );
+
+          // generate LUT for module number -> sensor position, for readout mapping tests
+          //if(isec==0) printf("%d %f %f\n",imod,sensorPV.position().x(),sensorPV.position().y());
+
+          // properties
+          sensorPV.addPhysVolID("sector", isec).addPhysVolID("module", imod);
+          DetElement sensorDE(det, Form("sensor_de%d_%d", isec, imod), (imod<<3)|isec ); // id must match IRTAlgorithm usage
+          sensorDE.setPlacement(sensorPV);
+          if(!debug_optics) {
+            SkinSurface sensorSkin(desc, sensorDE, Form("sensor_optical_surface%d", isec), sensorSurf, sensorVol);
+            sensorSkin.isValid();
+          };
+
+          // increment sensor module number
+          imod++;
+
+        }; // end patch cuts
+      }; // end phiGen loop
+    }; // end thetaGen loop
+
+    // calculate centroid sensor position
+    if(isec==0) {
+      sensorCentroidX /= sensorCount;
+      sensorCentroidZ /= sensorCount;
+    };
+
+    // END SENSOR MODULE LOOP ------------------------
+
+
+    // BUILD MIRRORS ====================================================================
+
+    // derive spherical mirror parameters `(zM,xM,rM)`, for given image point
+    // coordinates `(zI,xI)` and `dO`, defined as the z-distance between the
+    // object and the mirror surface
+    // - all coordinates are specified w.r.t. the object point coordinates
+    // - this is point-to-point focusing, but it can be used to effectively steer
+    //   parallel-to-point focusing
+    double zM,xM,rM;
+    auto FocusMirror = [&zM,&xM,&rM](double zI, double xI, double dO) {
+      zM = dO*zI / (2*dO-zI);
+      xM = dO*xI / (2*dO-zI);
+      rM = dO - zM;
+    };
+
+    // attributes, re-defined w.r.t. IP, needed for mirror positioning
+    double zS = sensorSphCenterZ + vesselZmin; // sensor sphere attributes
+    double xS = sensorSphCenterX;
+    double rS = sensorSphRadius;
+    double B = vesselZmax - mirrorBackplane; // distance between IP and mirror back plane
+
+    // focus 1: set mirror to focus IP on center of sensor sphere `(zS,xS)`
+    /*double zF = zS;
+    double xF = xS;
+    FocusMirror(zF,xF,B);*/
+
+    // focus 2: move focal region along sensor sphere radius, according to `focusTuneLong`
+    // - specifically, along the radial line which passes through the approximate centroid
+    //   of the sensor region `(sensorCentroidZ,sensorCentroidX)`
+    // - `focusTuneLong` is the distance to move, given as a fraction of `sensorSphRadius`
+    // - `focusTuneLong==0` means `(zF,xF)==(zS,xS)`
+    // - `focusTuneLong==1` means `(zF,xF)` will be on the sensor sphere, near the centroid
+    /*
+    double zC = sensorCentroidZ + vesselZmin;
+    double xC = sensorCentroidX;
+    double slopeF = (xC-xS) / (zC-zS);
+    double thetaF = std::atan(std::fabs(slopeF));
+    double zF = zS + focusTuneLong * sensorSphRadius * std::cos(thetaF);
+    double xF = xS - focusTuneLong * sensorSphRadius * std::sin(thetaF);
+    //FocusMirror(zF,xF,B);
+
+    // focus 3: move along line perpendicular to focus 2's radial line,
+    // according to `focusTunePerp`, with the same numerical scale as `focusTuneLong`
+    zF += focusTunePerp * sensorSphRadius * std::cos(M_PI/2-thetaF);
+    xF += focusTunePerp * sensorSphRadius * std::sin(M_PI/2-thetaF);
+    FocusMirror(zF,xF,B);
+    */
+
+    // focus 4: use (z,x) coordinates for tune parameters
+    double zF = zS + focusTuneZ;
+    double xF = xS + focusTuneX;
+    FocusMirror(zF,xF,B);
+
+    // re-define mirror attributes to be w.r.t vessel front plane
+    double mirrorCenterZ = zM - vesselZmin;
+    double mirrorCenterX = xM;
+    double mirrorRadius = rM;
+
+    // spherical mirror patch cuts and rotation
+    double mirrorThetaRot = std::asin(mirrorCenterX/mirrorRadius);
+    double mirrorTheta1 = mirrorThetaRot - std::asin((mirrorCenterX-mirrorRmin)/mirrorRadius);
+    double mirrorTheta2 = mirrorThetaRot + std::asin((mirrorRmax-mirrorCenterX)/mirrorRadius);
+
+    // if debugging, draw full sphere
+    if(debug_mirror) { mirrorTheta1=0; mirrorTheta2=M_PI; /*mirrorPhiw=2*M_PI;*/ };
+
+    // solid : create sphere at origin, with specified angular limits;
+    // phi limits are increased to fill gaps (overlaps are cut away later)
+    Sphere mirrorSolid1(
+        mirrorRadius,
+        mirrorRadius + mirrorThickness,
+        mirrorTheta1,
+        mirrorTheta2,
+        -40*degree,
+        40*degree
+        );
+
+    /* CAUTION: if any of the relative placements or boolean operations below
+     * are changed, you MUST make sure this does not break access to the sphere
+     * primitive and positioning in Juggler `IRTAlgorithm`; cross check the
+     * mirror sphere attributes carefully!
+     */
+    /*
+    // PRINT MIRROR ATTRIBUTES (before any sector z-rotation)
+    printf("zM = %f\n",zM); // sphere centerZ, w.r.t. IP
+    printf("xM = %f\n",xM); // sphere centerX, w.r.t. IP
+    printf("rM = %f\n",rM); // sphere radius
+    */
+
+    // mirror placement transformation (note: transformations are in reverse order)
+    auto mirrorPos = Position(mirrorCenterX, 0., mirrorCenterZ) + originFront;
+    Transform3D mirrorPlacement(
+          Translation3D(mirrorPos.x(), mirrorPos.y(), mirrorPos.z()) // re-center to specified position
+        * RotationY(-mirrorThetaRot) // rotate about vertical axis, to be within vessel radial walls
+        );
+
+    // cut overlaps with other sectors using "pie slice" wedges, to the extent specified
+    // by `mirrorPhiw`
+    Tube pieSlice( 0.01*cm, vesselRmax2, tankLength/2.0, -mirrorPhiw/2.0, mirrorPhiw/2.0);
+    IntersectionSolid mirrorSolid2( pieSlice, mirrorSolid1, mirrorPlacement );
+
+    // mirror volume, attributes, and placement
+    Volume mirrorVol(detName+"_mirror_"+secName, mirrorSolid2, mirrorMat);
+    mirrorVol.setVisAttributes(mirrorVis);
+    auto mirrorPV2 = gasvolVol.placeVolume(mirrorVol,
+          RotationZ(sectorRotation) // rotate about beam axis to sector
+        * Translation3D(0,0,0)
+        );
+
+    // properties
+    DetElement mirrorDE(det, Form("mirror_de%d", isec), isec);
+    mirrorDE.setPlacement(mirrorPV2);
+    SkinSurface mirrorSkin(desc, mirrorDE, Form("mirror_optical_surface%d", isec), mirrorSurf, mirrorVol);
+    mirrorSkin.isValid();
+
+
+  }; // END SECTOR LOOP //////////////////////////
+
+
+  // place gas volume
+  PlacedVolume gasvolPV = vesselVol.placeVolume(gasvolVol,Position(0, 0, 0));
+  DetElement gasvolDE(det, "gasvol_de", 0);
+  gasvolDE.setPlacement(gasvolPV);
+
+  // place mother volume (vessel)
+  Volume motherVol = desc.pickMotherVolume(det);
+  PlacedVolume vesselPV = motherVol.placeVolume(vesselVol,
+      Position(0, 0, vesselZmin) - originFront
+      );
+  vesselPV.addPhysVolID("system", detID);
+  det.setPlacement(vesselPV);
+
+  return det;
+};
+
+// clang-format off
+DECLARE_DETELEMENT(athena_DRICH, createDetector)

src/FieldMapBrBz.cpp

+#include <DD4hep/DetFactoryHelper.h>
+#include <DD4hep/FieldTypes.h>
+#include <DD4hep/Printout.h>
+#include <XML/Utilities.h>
+
+#include <cstdlib>
+#include <filesystem>
+#include <fstream>
+#include <iostream>
+#include <sstream>
+#include <stdexcept>
+#include <string>
+#include <tuple>
+namespace fs = std::filesystem;
+
+#include "FileLoaderHelper.h"
+
+using namespace dd4hep;
+
+
+// implementation of the field map
+class FieldMapBrBz : public dd4hep::CartesianField::Object
+{
+public:
+    FieldMapBrBz(const std::string &field_type = "magnetic");
+    void Configure(double rmin, double rmax, double rstep, double zmin, double zmax, double zstep);
+    void LoadMap(const std::string &map_file, double scale);
+    void GetIndices(double r, double z, int &ir, int &iz, double &dr, double &dz);
+    void SetTransform(const Transform3D &tr) { trans = tr; trans_inv = tr.Inverse(); }
+
+    virtual void fieldComponents(const double *pos, double *field);
+
+
+private:
+    Transform3D trans, trans_inv;
+    double rmin, rmax, rstep, zmin, zmax, zstep;
+    std::vector<std::vector<std::vector<double>>> Bvals;
+};
+
+// constructor
+FieldMapBrBz::FieldMapBrBz(const std::string &field_type)
+{
+    std::string ftype = field_type;
+    for (auto &c : ftype) { c = tolower(c); }
+
+    // set type
+    if (ftype == "magnetic") {
+        type = CartesianField::MAGNETIC;
+    } else if (ftype == "electric") {
+        type = CartesianField::ELECTRIC;
+    } else {
+        type = CartesianField::UNKNOWN;
+        std::cout << "FieldMapBrBz Warning: Unknown field type " << field_type << "!" << std::endl;
+    }
+}
+
+void FieldMapBrBz::Configure(double r1, double r2, double rs, double z1, double z2, double zs)
+{
+    rmin = r1;
+    rmax = r2;
+    rstep = rs;
+    zmin = z1;
+    zmax = z2;
+    zstep = zs;
+
+    int nr = int((r2 - r1)/rs) + 2;
+    int nz = int((z2 - z1)/zs) + 2;
+
+    Bvals.resize(nr);
+    for (auto &B2 : Bvals) {
+        B2.resize(nz);
+        for (auto &B : B2) {
+            B.resize(2, 0.);
+        }
+    }
+}
+
+void FieldMapBrBz::GetIndices(double r, double z, int &ir, int &iz, double &dr, double &dz)
+{
+    // boundary check
+    if (r > rmax || r < rmin || z > zmax || z < zmin) {
+        ir = -1;
+        iz = -1;
+        return;
+    }
+    // get indices
+    double idr, idz;
+    dr = std::modf((r - rmin)/rstep, &idr);
+    dz = std::modf((z - zmin)/zstep, &idz);
+
+    ir = static_cast<int>(idr);
+    iz = static_cast<int>(idz);
+}
+
+// load data
+void FieldMapBrBz::LoadMap(const std::string &map_file, double scale)
+{
+    std::string line;
+    std::ifstream input(map_file);
+    if (!input) {
+        std::cout << "FieldMapBrBz Error: file \"" << map_file << "\" cannot be read." << std::endl;
+    }
+
+    double r, z, br, bz;
+    int ir, iz;
+    double dr, dz;
+    while (std::getline(input, line).good()) {
+        std::istringstream iss(line);
+        iss >> r >> z >> br >> bz;
+        GetIndices(r, z, ir, iz, dr, dz);
+        if (ir < 0 || iz < 0) {
+            std::cout << "FieldMapBrBz Warning: coordinates out of range ("
+                      << r << ", " << z << "), skipped it." << std::endl;
+        } else {
+            Bvals[ir][iz] = {br*scale, bz*scale};
+            // ROOT::Math::XYZPoint p(r, 0, z);
+            // std::cout << p << " -> " << trans*p << std::endl;
+            // std::cout << ir << ", " << iz << ", " << br << ", " << bz << std::endl;
+        }
+    }
+}
+
+// get field components
+void FieldMapBrBz::fieldComponents(const double *pos, double *field)
+{
+    // coordinate conversion
+    auto p = trans_inv*ROOT::Math::XYZPoint(pos[0], pos[1], pos[2]);
+
+    // coordinates conversion
+    const double r = sqrt(p.x()*p.x() + p.y()*p.y());
+    const double z = p.z();
+    const double phi = atan2(p.y(), p.x());
+
+    int ir, iz;
+    double dr, dz;
+    GetIndices(r, z, ir, iz, dr, dz);
+
+    // out of the range
+    if (ir < 0 || iz < 0) { return; }
+
+    // p1    p3
+    //    p
+    // p0    p2
+    auto &p0 = Bvals[ir][iz];
+    auto &p1 = Bvals[ir][iz + 1];
+    auto &p2 = Bvals[ir + 1][iz];
+    auto &p3 = Bvals[ir + 1][iz + 1];
+
+    // linear interpolation
+    double Br = p0[0] * (1-dr) * (1-dz)
+              + p1[0] * (1-dr) *    dz
+              + p2[0] *    dr  * (1-dz)
+              + p3[0] *    dr  *    dz;
+
+    double Bz = p0[1] * (1-dr) * (1-dz)
+              + p1[1] * (1-dr) *    dz
+              + p2[1] *    dr  * (1-dz)
+              + p3[1] *    dr  *    dz;
+
+    // convert Br Bz to Bx By Bz
+    auto B = trans*ROOT::Math::XYZPoint(Br*sin(phi), Br*cos(phi), Bz);
+    field[0] += B.x()*tesla;
+    field[1] += B.y()*tesla;
+    field[2] += B.z()*tesla;
+    return;
+}
+
+
+// assign the field map to CartesianField
+static Ref_t create_field_map_brbz(Detector & /*lcdd*/, xml::Handle_t handle)
+{
+    xml_comp_t x_par(handle);
+
+    if (!x_par.hasAttr(_Unicode(field_map))) {
+        throw std::runtime_error("FieldMapBrBz Error: must have an xml attribute \"field_map\" for the field map.");
+    }
+
+    CartesianField field;
+    std::string field_type = x_par.attr<std::string>(_Unicode(field_type));
+
+    // dimensions
+    xml_comp_t x_dim = x_par.dimensions();
+
+    // min, max, step
+    xml_comp_t r_dim = x_dim.child(_Unicode(transverse));
+    xml_comp_t z_dim = x_dim.child(_Unicode(longitudinal));
+
+    std::string field_map_file = x_par.attr<std::string>(_Unicode(field_map));
+    std::string field_map_url = x_par.attr<std::string>(_Unicode(url));
+
+    EnsureFileFromURLExists(field_map_url,field_map_file);
+
+    double field_map_scale = x_par.attr<double>(_Unicode(scale));
+
+    if( !fs::exists(fs::path(field_map_file))  ) {
+        printout(ERROR, "FieldMapBrBz", "file " + field_map_file + " does not exist");
+        printout(ERROR, "FieldMapBrBz", "use a FileLoader plugin before the field element");
+        std::_Exit(EXIT_FAILURE);
+    }
+
+    auto map = new FieldMapBrBz(field_type);
+    map->Configure(r_dim.rmin(), r_dim.rmax(), r_dim.step(), z_dim.zmin(), z_dim.zmax(), z_dim.step());
+
+    // translation, rotation
+    static double deg2r = ROOT::Math::Pi()/180.;
+    RotationZYX rot(0., 0., 0.);
+    if (x_dim.hasChild(_Unicode(rotation))) {
+        xml_comp_t rot_dim = x_dim.child(_Unicode(rotation));
+        rot = RotationZYX(rot_dim.z()*deg2r, rot_dim.y()*deg2r, rot_dim.x()*deg2r);
+    }
+
+    Translation3D trans(0., 0., 0.);
+    if (x_dim.hasChild(_Unicode(translation))) {
+        xml_comp_t trans_dim = x_dim.child(_Unicode(translation));
+        trans = Translation3D(trans_dim.x(), trans_dim.y(), trans_dim.z());
+    }
+    map->SetTransform(trans*rot);
+
+    map->LoadMap(field_map_file, field_map_scale);
+    field.assign(map, x_par.nameStr(), "FieldMapBrBz");
+
+    return field;
+}
+
+DECLARE_XMLELEMENT(FieldMapBrBz, create_field_map_brbz)
+

src/FileLoader.cpp

+#include <DD4hep/DetFactoryHelper.h>
+#include <DD4hep/Primitives.h>
+#include <DD4hep/Factories.h>
+#include <DD4hep/Printout.h>
+#include <XML/Utilities.h>
+
+#include <fmt/core.h>
+
+#include <filesystem>
+#include <iostream>
+#include <cstdlib>
+#include <string>
+
+#include "FileLoaderHelper.h"
+
+using namespace dd4hep;
+
+void usage(int argc, char** argv) {
+  std::cout <<
+    "Usage: -plugin <name> -arg [-arg]                                                  \n"
+    "     name:   factory name     FileLoader                                           \n"
+    "     cache:<string>           cache location (may be read-only)                    \n"
+    "     file:<string>            file location                                        \n"
+    "     url:<string>             url location                                         \n"
+    "     cmd:<string>             download command with {0} for url, {1} for output    \n"
+    "\tArguments given: " << arguments(argc,argv) << std::endl;
+  std::exit(EINVAL);
+}
+
+// Plugin to download files
+long load_file(
+    Detector& /* desc */,
+    int argc,
+    char** argv
+) {
+  // argument parsing
+  std::string cache, file, url;
+  std::string cmd("curl --retry 5 -f {0} -o {1}");
+  for (int i = 0; i < argc && argv[i]; ++i) {
+    if      (0 == std::strncmp("cache:", argv[i], 6)) cache = (argv[i] + 6);
+    else if (0 == std::strncmp("file:", argv[i], 5)) file = (argv[i] + 5);
+    else if (0 == std::strncmp("url:", argv[i], 4)) url = (argv[i] + 4);
+    else if (0 == std::strncmp("cmd:", argv[i], 4)) cmd = (argv[i] + 4);
+    else usage(argc, argv);
+  }
+  printout(DEBUG, "FileLoader", "arg cache: " + cache);
+  printout(DEBUG, "FileLoader", "arg file: " + file);
+  printout(DEBUG, "FileLoader", "arg url: " + url);
+  printout(DEBUG, "FileLoader", "arg cmd: " + cmd);
+
+  // if file or url is empty, do nothing
+  if (file.empty()) {
+    printout(WARNING, "FileLoader", "no file specified");
+    return 0;
+  }
+  if (url.empty()) {
+    printout(WARNING, "FileLoader", "no url specified");
+    return 0;
+  }
+
+  EnsureFileFromURLExists(url, file, cache, cmd);
+
+  return 0;
+}
+
+DECLARE_APPLY(FileLoader, load_file)

src/FileLoaderHelper.h

+#pragma once
+
+#include <DD4hep/DetFactoryHelper.h>
+#include <DD4hep/Primitives.h>
+#include <DD4hep/Factories.h>
+#include <DD4hep/Printout.h>
+
+#include <fmt/core.h>
+
+#include <filesystem>
+#include <iostream>
+#include <cstdlib>
+#include <string>
+
+namespace fs = std::filesystem;
+
+using namespace dd4hep;
+
+// Function to download files
+inline void
+EnsureFileFromURLExists(
+  std::string url,
+  std::string file,
+  std::string cache = "",
+  std::string cmd = "curl --retry 5 -f {0} -o {1}"
+) {
+  // parse cache for environment variables
+  auto pos = std::string::npos;
+  while ((pos = cache.find('$')) != std::string::npos) {
+    auto after = cache.find_first_not_of(
+      "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
+      "abcdefghijklmnopqrstuvwxyz"
+      "0123456789"
+      "_",
+      pos + 1);
+    if (after == std::string::npos) after = cache.size(); // cache ends on env var
+    const std::string env_name(cache.substr(pos + 1, after - pos - 1));
+    auto env_ptr = std::getenv(env_name.c_str());
+    const std::string env_value(env_ptr != nullptr ? env_ptr : "");
+    cache.erase(pos, after - pos);
+    cache.insert(pos, env_value);
+    printout(INFO, "FileLoader", "$" + env_name + " -> " + env_value);
+  }
+
+  // create file path
+  fs::path file_path(file);
+
+  // create hash from url, hex of unsigned long long
+  std::string hash = fmt::format("{:016x}", dd4hep::detail::hash64(url)); // TODO: Use c++20 std::fmt
+
+  // create file parent path, if not exists
+  fs::path parent_path = file_path.parent_path();
+  if (!fs::exists(parent_path)) {
+    if (fs::create_directories(parent_path) == false) {
+      printout(ERROR, "FileLoader", "parent path " + parent_path.string() + " cannot be created");
+      printout(ERROR, "FileLoader", "check permissions and retry");
+      std::_Exit(EXIT_FAILURE);
+    }
+  }
+
+  // if file exists and is symlink to correct hash
+  fs::path hash_path(parent_path / hash);
+  if (fs::exists(file_path)
+   && fs::equivalent(file_path, hash_path)) {
+    printout(INFO, "FileLoader", "Link " + file + " -> hash " + hash + " already exists");
+    return;
+  }
+
+  // if hash does not exist, we try to retrieve file from cache
+  if (!fs::exists(hash_path)) {
+    // recursive loop into cache directory
+    fs::path cache_path(cache);
+    printout(INFO, "FileLoader", "Cache " + cache_path.string());
+    if (fs::exists(cache_path)) {
+      for (auto const& dir_entry: fs::recursive_directory_iterator(cache_path)) {
+        if (!dir_entry.is_directory()) continue;
+        fs::path cache_dir_path = cache_path / dir_entry;
+        printout(INFO, "FileLoader", "Checking " + cache_dir_path.string());
+        fs::path cache_hash_path = cache_dir_path / hash;
+        if (fs::exists(cache_hash_path)) {
+          // symlink hash to cache/.../hash
+          printout(INFO, "FileLoader", "File " + file + " with hash " + hash + " found in " + cache_hash_path.string());
+          try {
+            fs::create_symlink(cache_hash_path, hash_path);
+          } catch (const fs::filesystem_error&) {
+            printout(ERROR, "FileLoader", "unable to link from " + hash_path.string() + " to " + cache_hash_path.string());
+            printout(ERROR, "FileLoader", "check permissions and retry");
+            std::_Exit(EXIT_FAILURE);
+          }
+          break;
+        }
+      }
+    }
+  }
+
+  // if hash does not exist, we try to retrieve file from url
+  if (!fs::exists(hash_path)) {
+    cmd = fmt::format(cmd, url, hash_path.c_str()); // TODO: Use c++20 std::fmt
+    printout(INFO, "FileLoader", "Downloading " + file + " as hash " + hash + " with " + cmd);
+    // run cmd
+    auto ret = std::system(cmd.c_str());
+    if (!fs::exists(hash_path)) {
+      printout(ERROR, "FileLoader", "unable to run cmd " + cmd);
+      printout(ERROR, "FileLoader", "check command and retry");
+      printout(ERROR, "FileLoader", "hint: allow insecure connections with -k");
+      std::_Exit(EXIT_FAILURE);
+    }
+  }
+
+  // check if file already exists
+  if (fs::exists(file_path)) {
+    // file already exists
+    if (fs::is_symlink(file_path)) {
+      // file is symlink
+      if (fs::equivalent(hash_path, fs::read_symlink(file_path))) {
+        // link points to correct path
+        return;
+      } else {
+        // link points to incorrect path 
+        if (fs::remove(file_path) == false) {
+          printout(ERROR, "FileLoader", "unable to remove symlink " + file_path.string());
+          printout(ERROR, "FileLoader", "check permissions or remove manually");
+          std::_Exit(EXIT_FAILURE);
+        }
+      }
+    } else {
+      // file exists but not symlink
+      printout(ERROR, "FileLoader", "will not remove actual file " + file_path.string());
+      printout(ERROR, "FileLoader", "check content, remove manually, and retry");
+      std::_Exit(EXIT_FAILURE);
+    }
+  }
+  // file_path now does not exist
+
+  // symlink file_path to hash_path
+  try {
+    // use new path from hash so file link is local
+    fs::create_symlink(fs::path(hash), file_path);
+  } catch (const fs::filesystem_error&) {
+    printout(ERROR, "FileLoader", "unable to link from " + file_path.string() + " to " + hash_path.string());
+    printout(ERROR, "FileLoader", "check permissions and retry");
+    std::_Exit(EXIT_FAILURE);
+  }
+}

src/GaseousRICH_geo.cpp

@@ -75,7 +75,7 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
     envVol.setVisAttributes(desc.visAttributes(detElem.visStr()));
 
     // sensitive detector type
-    sens.setType("photoncounter");
+    sens.setType("tracker");
 
     // @TODO: place a radiator
     // build_radiator(desc, envVol, detElement.child(_Unicode(radiator)), snout_front);

src/GeometryHelpers.h

@@ -5,9 +5,18 @@
 // some utility functions that can be shared
 namespace athena::geo {
 
-typedef ROOT::Math::XYPoint Point;
+using Point =  ROOT::Math::XYPoint;
 
-// fill rectangles in a ring
+/** Fill rectangles in a ring (disk).
+ *
+ * @param ref  2D reference point.
+ * @param sx   x side length
+ * @param sy   y side length
+ * @param rmin inner radius of disk
+ * @param rmax  outer radius of disk to fill
+ * @param phmin  phi min
+ * @param phmax phi max
+ */
 std::vector<Point> fillRectangles(Point ref, double sx, double sy, double rmin, double rmax,
                                   double phmin = -M_PI, double phmax = M_PI);
 // fill squares in a ring