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/ForwardRICH_geo.cpp

-//==========================================================================
-//  Forward Ring Imaging Cherenkov Detector
-//--------------------------------------------------------------------------
-//
-// Author: C. Peng (ANL)
-// Date: 09/30/2020
-//
-//==========================================================================
-
-#include <XML/Helper.h>
-#include "TMath.h"
-#include "TString.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 std;
-using namespace dd4hep;
-using namespace dd4hep::rec;
-
-typedef ROOT::Math::XYPoint Point;
-
-
-// check if a square in a ring
-inline bool in_ring(const Point &pt, double side, double rmin, double rmax, double phmin, double phmax)
-{
-    if (pt.r() > rmax || pt.r() < rmin) {
-        return false;
-    }
-
-    // check four corners
-    std::vector<Point> pts {
-        Point(pt.x() - side/2., pt.y() - side/2.),
-        Point(pt.x() - side/2., pt.y() + side/2.),
-        Point(pt.x() + side/2., pt.y() - side/2.),
-        Point(pt.x() + side/2., pt.y() + side/2.),
-    };
-    for (auto &p : pts) {
-        if (p.r() > rmax || p.r() < rmin || p.phi() > phmax || p.phi() < phmin) {
-            return false;
-        }
-    }
-    return true;
-}
-
-// check if a square is overlapped with the others
-inline bool overlap(const Point &pt, double side, const std::vector<Point> &pts)
-{
-    for (auto &p : pts) {
-        auto pn = (p - pt)/side;
-        if ((std::abs(pn.x()) < 1. - 1e-6) && (std::abs(pn.y()) < 1. - 1e-6)) {
-            return true;
-        }
-    }
-    return false;
-}
-
-// a helper function to recursively fill square in a ring
-void add_square(Point p, std::vector<Point> &res, double lside, double rmin, double rmax,
-                double phmin, double phmax)
-{
-    // outside of the ring or overlapping
-    if (!in_ring(p, lside, rmin, rmax, phmin, phmax) || overlap(p, lside, res)) {
-        return;
-    }
-
-    res.emplace_back(p);
-
-    // check adjacent squares
-    add_square(Point(p.x() + lside, p.y()), res, lside, rmin, rmax, phmin, phmax);
-    add_square(Point(p.x() - lside, p.y()), res, lside, rmin, rmax, phmin, phmax);
-    add_square(Point(p.x(), p.y() + lside), res, lside, rmin, rmax, phmin, phmax);
-    add_square(Point(p.x(), p.y() - lside), res, lside, rmin, rmax, phmin, phmax);
-}
-
-// fill squares
-std::vector<Point> fill_squares(Point ref, double lside, double rmin, double rmax,
-                                double phmin = 0., double phmax = 2.*M_PI)
-{
-    // start with a seed square and find one in the ring
-    // move to center
-    ref = ref - Point(int(ref.x()/lside)*lside, int(ref.y()/lside)*lside);
-
-    auto find_seed = [] (const Point &ref, int n, double side, double rmin, double rmax, double phmin, double phmax) {
-        for (int ix = -n; ix < n; ++ix) {
-            for (int iy = -n; iy < n; ++iy) {
-                Point pt(ref.x() + ix*side, ref.y() + iy*side);
-                if (in_ring(pt, side, rmin, rmax, phmin, phmax)) {
-                    return pt;
-                }
-            }
-        }
-        return ref;
-    };
-
-    std::vector<Point> res;
-    ref = find_seed(ref, int(rmax/lside) + 2, lside, rmin, rmax, phmin, phmax);
-    add_square(ref, res, lside, rmin, rmax, phmin, phmax);
-    return res;
-}
-
-// 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();
-    xml::Component rads = detElem.child(_Unicode(radiator));
-    xml::Component mir = detElem.child(_Unicode(mirror));
-    xml::Component mcp = detElem.child(_Unicode(mcppmt));
-
-    // dimensions
-    double z0 = dims.z0();
-    double length = dims.length();
-    double rmin = dims.rmin();
-    double rmax1 = dims.attr<double>(_Unicode(rmax1));
-    double rmax2 = dims.attr<double>(_Unicode(rmax2));
-
-    // mirror setting
-    auto mThick = mir.thickness();
-    auto mirZ = mir.attr<double>(_Unicode(zdiff));
-
-    // mcppmt setting
-    auto pRmin = mcp.rmin();
-    auto pRmax = mcp.rmax();
-    auto pThick = mcp.thickness();
-    auto pSize = mcp.attr<double>(_Unicode(module_size));
-    auto pGap = mcp.attr<double>(_Unicode(module_gap));
-    auto pTol = mcp.attr<double>(_Unicode(rtol));
-    auto pZ = mcp.attr<double>(_Unicode(zdiff));
-
-    // materials
-    auto mirMat = desc.material(mir.materialStr());
-    auto gasMat = desc.material(rads.materialStr());
-    auto mcpMat = desc.material(mcp.materialStr());
-
-    // constants
-    auto richCenterAngle = std::atan((rmin + (rmax2 - rmin)/2.)/mirZ);
-    //std::cout << richCenterAngle*180./M_PI << std::endl;
-
-    // an envelope for the detector
-    // use a complicated shape to avoid conflict with the other parts
-    // cone for radiator and the first set of mirrors
-    double halfLength = length/2.;
-    Cone env1(halfLength, rmin, rmax1, rmin, rmax2);
-    // envelope for detection plane
-    // Cone env2(halfLength - pZ/2., rmin, pRmax, rmin, rmax2);
-    Tube env2(rmin, pRmax + pTol + pGap + 1.0*cm, (length - pZ)/2., 0., 2*M_PI);
-
-    UnionSolid envShape(env1, env2, Position(0., 0., pZ));
-
-    Volume envVol(detName + "_envelope", envShape, gasMat);
-    envVol.setVisAttributes(desc.visAttributes(detElem.visStr()));
-
-    // ---------------
-    // spherical mirrors inside it
-    int ilayer = 1;
-
-    // optical surface
-    OpticalSurfaceManager surfMgr = desc.surfaceManager();
-    OpticalSurface mirSurf  = surfMgr.opticalSurface("MirrorOpticalSurface");
-    // mirror slices
-    int imod = 1;
-    for (xml::Collection_t sl(mir, _Unicode(slice)); sl; ++sl, ++imod) {
-        auto focus = sl.attr<double>(_Unicode(focus));
-        auto wphi = sl.attr<double>(_Unicode(phiw));
-        auto rotZ = sl.attr<double>(_Unicode(rotz));
-        auto mRmin = sl.attr<double>(_Unicode(rmin));
-        auto mRmax = sl.attr<double>(_Unicode(rmax));
-        double curve = 0.;
-        if (sl.hasAttr(_Unicode(curve))) {
-            curve = sl.attr<double>(_Unicode(curve));
-        }
-        // geometry of mirror slice
-        PlacedVolume mirPV;
-        Volume mirVol(Form("mirror_v_dummy%d", imod));
-        mirVol.setMaterial(mirMat);
-        mirVol.setVisAttributes(desc.visAttributes(mir.visStr()));
-        // spherical mirror
-        if (curve > 0.) {
-            // somehow geant4 does not support -wphi/2. to wphi/2., so additonal rotation in Z
-            double mTheta1 = std::asin(mRmin/curve);
-            double mTheta2 = std::asin(mRmax/curve);
-            double rotY = -std::asin(focus/curve);
-            mirVol.setSolid(Sphere(curve, curve + mThick, mTheta1, mTheta2, 0., wphi));
-            // action is in a reverse order
-            Transform3D tr = Translation3D(0., 0., mirZ - halfLength)   // move for z position
-                           * RotationZ(rotZ)                            // rotate phi angle
-                           * RotationY(rotY)                            // rotate for focus point
-                           * RotationX(180*degree)
-                           * Translation3D(0., 0., -curve)              // move spherical shell to origin
-                           * RotationZ(-wphi/2.);                       // center phi angle to 0. (-wphi/2., wphi/2.)
-            mirPV = envVol.placeVolume(mirVol, tr);
-        // plane mirror
-        } else {
-            mirVol.setSolid(Tube(mRmin, mRmax, mThick/2.0, 0., wphi));
-            Transform3D tr = Translation3D(0., 0., mirZ - halfLength)   // move for z position
-                           * RotationZ(rotZ)                            // rotate phi angle
-                           * RotationZ(-wphi/2.);                       // center phi angle to 0. (-wphi/2., wphi/2.)
-            mirPV = envVol.placeVolume(mirVol, tr);
-        }
-        mirPV.addPhysVolID("layer", ilayer).addPhysVolID("module", imod);
-        DetElement mirDE(det, Form("Mirror_DE%d", imod), imod);
-        mirDE.setPlacement(mirPV);
-        SkinSurface mirSurfBorder(desc, mirDE, Form("RICHmirror%d", imod), mirSurf, mirVol);
-        mirSurfBorder.isValid();
-    }
-    ilayer++;
-
-    // ---------------
-    // photo-detector unit
-    // Fill the photo-detection plane with square shape MCP-PMTs
-    Box mcpShape1(pSize/2.0, pSize/2.0, pThick/2.0);
-    Volume mcpVol1("mcppmt_v_material", mcpShape1, mcpMat);
-
-    // a thin layer of cherenkov gas for accepting optical photons
-    Box mcpShape(pSize/2.0, pSize/2.0, pThick/2.0 + 0.1*mm);
-    Volume mcpVol("mcppmt_v", mcpShape, gasMat);
-    mcpVol.placeVolume(mcpVol1, Position(0., 0., -0.1*mm));
-
-    mcpVol.setVisAttributes(desc.visAttributes(mcp.visStr()));
-    sens.setType("photoncounter");
-    mcpVol.setSensitiveDetector(sens);
-
-    // photo-detector plane envelope
-    for (size_t ipd = 0; ipd < 6; ++ipd) {
-        double phmin = -M_PI/6.;
-        double phmax = M_PI/6.;
-        Tube pdEnvShape(pRmin - pTol - pGap, pRmax + pTol + pGap, pThick/2.0 + 0.1*cm, phmin, phmax);
-        Volume pdVol("pd_envelope", pdEnvShape, desc.material("AirOptical"));
-        auto points = fill_squares(Point(0., 0.), pSize + pGap, pRmin - pTol - pGap, pRmax + pTol + pGap, phmin, phmax);
-        for (size_t i = 0; i < points.size(); ++i) {
-            auto pt = points[i];
-            auto mcpPV = pdVol.placeVolume(mcpVol, Position(pt.x(), pt.y(), 0.));
-            mcpPV.addPhysVolID("layer", ilayer).addPhysVolID("module", i + 1);
-            DetElement mcpDE(det, Form("MCPPMT_DE%d_%d", ipd + 1, i + 1), i + 1);
-            mcpDE.setPlacement(mcpPV);
-        }
-        Transform3D tr = Translation3D(0., 0., -halfLength + pZ + pThick/2.0)   // move for z position
-                        * RotationZ(ipd*M_PI/3.)        // rotate phi angle
-                        * RotationY(-richCenterAngle);  // rotate to perpendicular position
-        auto pdPV = envVol.placeVolume(pdVol, tr);
-        pdPV.addPhysVolID("layer", ilayer).addPhysVolID("piece", ipd + 1);
-    }
-    Volume motherVol = desc.pickMotherVolume(det);
-    PlacedVolume envPV = motherVol.placeVolume(envVol, Position(0, 0, z0 + halfLength));
-    envPV.addPhysVolID("system", detID);
-    det.setPlacement(envPV);
-
-    return det;
-}
-//@}
-
-// clang-format off
-DECLARE_DETELEMENT(refdet_ForwardRICH, createDetector)
-

src/GaseousRICH_geo.cpp

+//==========================================================================
+//  Gaseous Ring Imaging Cherenkov Detector
+//--------------------------------------------------------------------------
+//
+// Author: C. Peng (ANL)
+// Date: 09/30/2020
+//
+//==========================================================================
+
+#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 std;
+using namespace dd4hep;
+using namespace dd4hep::rec;
+
+// headers
+void build_radiator(Detector &desc, Volume &env, xml::Component plm, const Position &offset);
+void build_mirrors(Detector &desc, DetElement &sdet, Volume &env, xml::Component plm, const Position &offset);
+void build_sensors(Detector &desc, Volume &env, xml::Component plm, const Position &offset, SensitiveDetector &sens);
+
+// helper function to get x, y, z if defined in a xml component
+template<class XmlComp>
+Position get_xml_xyz(XmlComp &comp, dd4hep::xml::Strng_t name)
+{
+    Position pos(0., 0., 0.);
+    if (comp.hasChild(name)) {
+        auto child = comp.child(name);
+        pos.SetX(dd4hep::getAttrOrDefault<double>(child, _Unicode(x), 0.));
+        pos.SetY(dd4hep::getAttrOrDefault<double>(child, _Unicode(y), 0.));
+        pos.SetZ(dd4hep::getAttrOrDefault<double>(child, _Unicode(z), 0.));
+    }
+    return pos;
+}
+
+// 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();
+
+    // build a big envelope for all the components, filled with optical material to ensure transport of optical photons
+    double z0 = dims.z0();
+    double length = dims.length();
+    double rmin = dims.rmin();
+    double rmax0 = dims.attr<double>(_Unicode(rmax0));
+    double rmax1 = dims.attr<double>(_Unicode(rmax1));
+    double rmax2 = dims.attr<double>(_Unicode(rmax2));
+    double snout_length = dims.attr<double>(_Unicode(snout_length));
+    // fill envelope with radiator materials (need optical property for optical photons)
+    auto gasMat = desc.material(dd4hep::getAttrOrDefault(detElem, _Unicode(gas), "AirOptical"));
+
+    Cone snout(snout_length/2.0, rmin, rmax0, rmin, rmax1);
+    Tube tank(rmin, rmax2, (length - snout_length)/2., 0., 2*M_PI);
+    // shift the snout to the left side of tank
+    UnionSolid envShape(tank, snout, Position(0., 0., -length/2.));
+    // some reference points
+    auto snout_front = Position(0., 0., -(length + snout_length)/2.);
+    // auto snout_end = Position(0., 0., -(length - snout_length)/2.);  // tank_front
+    // auto tank_end = Position(0., 0., (length - snout_length)/2.);
+
+    Volume envVol(detName + "_envelope", envShape, gasMat);
+    envVol.setVisAttributes(desc.visAttributes(detElem.visStr()));
+
+    // sensitive detector type
+    sens.setType("tracker");
+
+    // @TODO: place a radiator
+    // build_radiator(desc, envVol, detElement.child(_Unicode(radiator)), snout_front);
+
+    // place mirrors
+    build_mirrors(desc, det, envVol, detElem.child(_Unicode(mirrors)), snout_front);
+
+    // place photo-sensors
+    build_sensors(desc, envVol, detElem.child(_Unicode(sensors)), snout_front, sens);
+
+    Volume motherVol = desc.pickMotherVolume(det);
+    PlacedVolume envPV = motherVol.placeVolume(envVol, Position(0, 0, z0) - snout_front);
+    envPV.addPhysVolID("system", detID);
+    det.setPlacement(envPV);
+
+    return det;
+}
+
+// @TODO: implement a radiator, now the envelope serves as the radiator
+void build_radiator(Detector &desc, Volume &env, xml::Component plm, const Position &offset)
+{
+    // place holder
+}
+
+// place mirrors
+void build_mirrors(Detector &desc, DetElement &sdet, Volume &env, xml::Component plm, const Position &offset)
+{
+    double thickness = dd4hep::getAttrOrDefault<double>(plm, _Unicode(dz), 1.*dd4hep::mm);
+    auto mat = desc.material(plm.attr<std::string>(_Unicode(material)));
+    auto vis = desc.visAttributes(plm.attr<std::string>(_Unicode(vis)));
+
+    // optical surface
+    OpticalSurfaceManager surfMgr = desc.surfaceManager();
+    auto surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault(plm, _Unicode(surface), "MirrorOpticalSurface"));
+
+    // placements
+    auto gpos = get_xml_xyz(plm, _Unicode(position)) + offset;
+    auto grot = get_xml_xyz(plm, _Unicode(position));
+    int imod = 1;
+    for (xml::Collection_t mir(plm, _Unicode(mirror)); mir; ++mir, ++imod) {
+        double rmin = mir.attr<double>(_Unicode(rmin));
+        double rmax = mir.attr<double>(_Unicode(rmax));
+        double phiw = dd4hep::getAttrOrDefault<double>(mir, _Unicode(phiw), 2.*M_PI);
+        double rotz = dd4hep::getAttrOrDefault<double>(mir, _Unicode(rotz), 0.);
+        double roty = dd4hep::getAttrOrDefault<double>(mir, _Unicode(roty), 0.);
+        double rotx = dd4hep::getAttrOrDefault<double>(mir, _Unicode(rotx), 0.);
+
+        Volume vol(Form("mirror_v_dummy%d", imod));
+        vol.setMaterial(mat);
+        vol.setVisAttributes(vis);
+        // mirror curvature
+        double curve = dd4hep::getAttrOrDefault<double>(mir, _Unicode(curve), 0.);
+        // spherical mirror
+        if (curve > 0.) {
+            double th1 = std::asin(rmin/curve);
+            double th2 = std::asin(rmax/curve);
+            vol.setSolid(Sphere(curve, curve + thickness, th1, th2, 0., phiw));
+        // plane mirror
+        } else {
+            vol.setSolid(Tube(rmin, rmax, thickness/2., 0., phiw));
+        }
+        // transforms are in a reverse order
+        Transform3D tr = Translation3D(gpos.x(), gpos.y(), gpos.z())
+                       * RotationZYX(grot.z(), grot.y(), grot.x())
+                       * RotationZ(rotz)                // rotation of the piece
+                       * RotationY(roty)                // rotation of the piece
+                       * RotationX(rotx)                // rotation of the piece
+                       * Translation3D(0., 0., -curve)  // move spherical shell to origin (no move for planes)
+                       * RotationZ(-phiw/2.);           // center phi angle to 0. (-phiw/2., phiw/2.)
+        auto pv = env.placeVolume(vol, tr);
+        DetElement de(sdet, Form("mirror_de%d", imod), imod);
+        de.setPlacement(pv);
+        SkinSurface skin(desc, de, Form("mirror_optical_surface%d", imod), surf, vol);
+        skin.isValid();
+    }
+}
+
+// place photo-sensors
+void build_sensors(Detector &desc, Volume &env, xml::Component plm, const Position &offset, SensitiveDetector &sens)
+{
+    // build sensor unit geometry
+    auto mod = plm.child(_Unicode(module));
+    double sx = mod.attr<double>(_Unicode(sx));
+    double sy = mod.attr<double>(_Unicode(sy));
+    double sz = mod.attr<double>(_Unicode(sz));
+    double gap = mod.attr<double>(_Unicode(gap));
+    auto mat = desc.material(mod.attr<std::string>(_Unicode(material)));
+    auto vis = desc.visAttributes(mod.attr<std::string>(_Unicode(vis)));
+
+    Box sensor(sx/2., sy/2., sz/2.);
+    Volume svol("sensor_v", sensor, mat);
+    svol.setVisAttributes(vis);
+
+    // a thin layer of cherenkov gas for accepting optical photons
+    auto opt = plm.child(_Unicode(optical));
+    double opthick = opt.attr<double>(_Unicode(thickness));
+    auto opmat = desc.material(opt.attr<std::string>(_Unicode(material)));
+
+    Box opshape(sx/2., sy/2., sz/2. + opthick/2.);
+    Volume opvol("sensor_v_optical", opshape, opmat);
+    opvol.placeVolume(svol, Position(0., 0., 0.));
+    opvol.setSensitiveDetector(sens);
+
+    // photo-detector plane envelope
+    auto gpos = get_xml_xyz(plm, _Unicode(position)) + offset;
+    auto grot = get_xml_xyz(plm, _Unicode(position));
+    int isec = 1;
+    for (xml::Collection_t sec(plm, _Unicode(sector)); sec; ++sec, ++isec) {
+        double rmin = sec.attr<double>(_Unicode(rmin));
+        double rmax = sec.attr<double>(_Unicode(rmax));
+        double phiw = dd4hep::getAttrOrDefault<double>(sec, _Unicode(phiw), 2.*M_PI);
+        double rotz = dd4hep::getAttrOrDefault<double>(sec, _Unicode(rotz), 0.);
+        double roty = dd4hep::getAttrOrDefault<double>(sec, _Unicode(roty), 0.);
+        double rotx = dd4hep::getAttrOrDefault<double>(sec, _Unicode(rotx), 0.);
+
+        // fill sensors to the piece
+        auto points = athena::geo::fillRectangles({0., 0.}, sx + gap, sy + gap, rmin - gap, rmax + gap, -phiw/2., phiw/2.);
+        int imod = 1;
+        for (auto &p : points) {
+            // transofrms are in a reversed order
+            Transform3D tr = Translation3D(gpos.x(), gpos.y(), gpos.z())
+                           * RotationZYX(grot.z(), grot.y(), grot.x())
+                           * RotationZ(rotz)                    // rotation of the sector
+                           * RotationY(roty)                    // rotation of the sector
+                           * RotationX(rotx)                    // rotation of the sector
+                           * Translation3D(p.x(), p.y(), 0.);   // place modules in each sector
+            auto pv = env.placeVolume(opvol, tr);
+            pv.addPhysVolID("sector", isec).addPhysVolID("module", imod++);
+        }
+    }
+}
+
+//@}
+
+// clang-format off
+DECLARE_DETELEMENT(athena_GaseousRICH, createDetector)
+

src/GeometryHelpers.cpp

+#include "GeometryHelpers.h"
+
+// some utility functions that can be shared
+namespace athena::geo {
+
+typedef ROOT::Math::XYPoint Point;
+
+// check if a 2d point is already in the container
+bool already_placed(const Point &p, const std::vector<Point> &vec, double xs = 1.0, double ys = 1.0, double tol = 1e-6)
+{
+    for (auto &pt : vec) {
+        if ((std::abs(pt.x() - p.x())/xs < tol) && std::abs(pt.y() - p.y())/ys < tol) {
+            return true;
+        }
+    }
+    return false;
+}
+
+// check if a square in a ring
+inline bool rec_in_ring(const Point &pt, double sx, double sy, double rmin, double rmax, double phmin, double phmax)
+{
+    if (pt.r() > rmax || pt.r() < rmin) {
+        return false;
+    }
+
+    // check four corners
+    std::vector<Point> pts {
+        Point(pt.x() - sx/2., pt.y() - sy/2.),
+        Point(pt.x() - sx/2., pt.y() + sy/2.),
+        Point(pt.x() + sx/2., pt.y() - sy/2.),
+        Point(pt.x() + sx/2., pt.y() + sy/2.),
+    };
+    for (auto &p : pts) {
+        if (p.r() > rmax || p.r() < rmin || p.phi() > phmax || p.phi() < phmin) {
+            return false;
+        }
+    }
+    return true;
+}
+
+// a helper function to recursively fill square in a ring
+void add_rectangle(Point p, std::vector<Point> &res, double sx, double sy,
+                   double rmin, double rmax, double phmin, double phmax, int max_depth = 20, int depth = 0)
+{
+    // std::cout << depth << "/" << max_depth << std::endl;
+    // exceeds the maximum depth in searching or already placed
+    if ((depth > max_depth) || (already_placed(p, res, sx, sy))) {
+        return;
+    }
+
+    bool in_ring = rec_in_ring(p, sx, sy, rmin, rmax, phmin, phmax);
+    if (in_ring) {
+        res.emplace_back(p);
+    }
+
+    // continue search for a good placement or if no placement found yet
+    if (in_ring || res.empty()) {
+        // check adjacent squares
+        add_rectangle(Point(p.x() + sx, p.y()), res, sx, sy, rmin, rmax, phmin, phmax, max_depth, depth + 1);
+        add_rectangle(Point(p.x() - sx, p.y()), res, sx, sy, rmin, rmax, phmin, phmax, max_depth, depth + 1);
+        add_rectangle(Point(p.x(), p.y() + sy), res, sx, sy, rmin, rmax, phmin, phmax, max_depth, depth + 1);
+        add_rectangle(Point(p.x(), p.y() - sy), res, sx, sy, rmin, rmax, phmin, phmax, max_depth, depth + 1);
+    }
+}
+
+// fill squares
+std::vector<Point> fillRectangles(Point ref, double sx, double sy, double rmin, double rmax, double phmin, double phmax)
+{
+    // convert (0, 2pi) to (-pi, pi)
+    if (phmax > M_PI) {
+        phmin -= M_PI;
+        phmax -= M_PI;
+    }
+    // start with a seed square and find one in the ring
+    // move to center
+    ref = ref - Point(int(ref.x()/sx)*sx, int(ref.y()/sy)*sy);
+
+    std::vector<Point> res;
+    add_rectangle(ref, res, sx, sy, rmin, rmax, phmin, phmax, (int(rmax/sx) + 1)*(int(rmax/sy) + 1)*2);
+    return res;
+}
+
+// check if a regular polygon is inside a ring
+bool poly_in_ring(const Point &p, int nsides, double lside, double rmin, double rmax, double phmin, double phmax)
+{
+    // outer radius is contained
+    if ((p.r() + lside <= rmax) && (p.r() - lside >= rmin)) {
+        return true;
+    }
+
+    // inner radius is not contained
+    double rin = std::cos(M_PI/nsides)*lside;
+    if ((p.r() + rin > rmax) || (p.r() - rin < rmin)) {
+        return false;
+    }
+
+    // in between, check every corner
+    for (int i = 0; i < nsides; ++i) {
+        double phi = (i + 0.5)*2.*M_PI/static_cast<double>(nsides);
+        Point p2(p.x() + 2.*lside*std::sin(phi), p.y() + 2.*lside*std::cos(phi));
+        if ((p2.r() > rmax) || (p2.r() < rmin)) {
+            return false;
+        }
+    }
+    return true;
+}
+
+// recursively fill square (nside=4) or hexagon (nside=6) in a ring, other polygons won't work
+void add_poly(Point p, std::vector<Point> &res, int nsides, double lside,
+              double rmin, double rmax, double phmin, double phmax, int max_depth = 20, int depth = 0)
+{
+    // std::cout << depth << "/" << max_depth << std::endl;
+    // exceeds the maximum depth in searching or already placed
+    if ((depth > max_depth) || (already_placed(p, res, lside, lside))) {
+        return;
+    }
+
+    bool in_ring = poly_in_ring(p, nsides, lside, rmin, rmax, phmin, phmax);
+    if (in_ring) {
+        res.emplace_back(p);
+    }
+
+    // recursively add neigbors, continue if it was a good placement or no placement found yet
+    if (in_ring || res.empty()) {
+        double d = 2.*std::cos(M_PI/static_cast<double>(nsides))*lside;
+        for (int i = 0; i < nsides; ++i) {
+            double phi = i*2.*M_PI/static_cast<double>(nsides);
+            add_poly(Point(p.x() + 2.*lside*std::sin(phi), p.y() + 2.*lside*std::cos(phi)), res, nsides, lside,
+                     rmin, rmax, phmin, phmax, max_depth, depth + 1);
+        }
+    }
+}
+
+std::vector<Point> fillHexagons(Point ref, double lside, double rmin, double rmax, double phmin, double phmax)
+{
+    // convert (0, 2pi) to (-pi, pi)
+    if (phmax > M_PI) {
+        phmin -= M_PI;
+        phmax -= M_PI;
+    }
+    // start with a seed and find one in the ring
+    // move to center
+    ref = ref - Point(int(ref.x()/lside)*lside, int(ref.y()/lside)*lside);
+
+    std::vector<Point> res;
+    add_poly(ref, res, 6, lside, rmin, rmax, phmin, phmax, std::pow(int(rmax/lside) + 1, 2)*2);
+    return res;
+}
+
+} // athena::geo

src/GeometryHelpers.h

+#pragma once
+#include <vector>
+#include "Math/Point2D.h"
+
+// some utility functions that can be shared
+namespace athena::geo {
+
+using Point =  ROOT::Math::XYPoint;
+
+/** 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
+inline std::vector<Point> fillSquares(Point ref, double size, double rmin, double rmax,
+                                      double phmin = -M_PI, double phmax = M_PI)
+{
+    return fillRectangles(ref, size, size, rmin, rmax, phmin, phmax);
+}
+
+std::vector<Point> fillHexagons(Point ref, double lside, double rmin, double rmax,
+                                double phmin = -M_PI, double phmax = M_PI);
+
+} // athena::geo

src/HomogeneousCalorimeter_geo.cpp

+//==========================================================================
+//  A general implementation for homogeneous calorimeter
+//  it supports three types of placements
+//  1. Individual module placement with module dimensions and positions
+//  2. Array placement with module dimensions and numbers of row and columns
+//  3. Disk placement with module dimensions and (Rmin, Rmax), and (Phimin, Phimax)
+//  4. Lines placement with module dimensions and (mirrorx, mirrory)
+//     (NOTE: anchor point is the 0th block of the line instead of line center)
+//--------------------------------------------------------------------------
+//  Author: Chao Peng (ANL)
+//  Date: 06/09/2021
+//==========================================================================
+
+#include "GeometryHelpers.h"
+#include "DD4hep/DetFactoryHelper.h"
+#include <XML/Helper.h>
+#include <iostream>
+#include <algorithm>
+#include <tuple>
+#include <math.h>
+
+using namespace dd4hep;
+
+/** \addtogroup calorimeters Calorimeters
+ */
+
+/** \addtogroup Homogeneous Calorimeter
+ * \brief Type: **HomogeneousCalorimeter**.
+ * \author C. Peng
+ * \ingroup calorimeters
+ *
+ *
+ * \code
+ *   <detector id="1" name="HyCal" type="HomogeneousCalorimeter" readout="EcalHits">
+ *     <position x="0" y="0" z="0"/>
+ *     <rotation x="0" y="0" z="0"/>
+ *     <placements>
+ *       <array nrow="34" ncol="34" sector="1">
+ *         <position x="0" y="0" z="-9.73*cm"/>
+ *         <module sizex="2.05*cm" sizey="2.05*cm" sizez="18*cm" vis="GreenVis" material="PbWO4"/>
+ *         <wrapper thickness="0.015*cm" material="Epoxy" vis="WhiteVis"/>
+ *         <removal row="16" col="16"/>
+ *         <removal row="16" col="17"/>
+ *         <removal row="17" col="16"/>
+ *         <removal row="17" col="17"/>
+ *       </array>
+ *       <array nrow="6" ncol="24" sector="2">
+ *         <position x="-17*(2.05+0.015)*cm+12*(3.8+0.015)*cm" y="17*(2.05+0.015)*cm+3*(3.8+0.015)*cm" z="0"/>
+ *         <module sizex="3.8*cm" sizey="3.8*cm" sizez="45*cm" vis="BlueVis" material="PbGlass"/>
+ *         <wrapper thickness="0.015*cm" material="Epoxy" vis="WhiteVis"/>
+ *       </array>
+ *       <array nrow="24" ncol="6" sector="3">
+ *         <position x="17*(2.05+0.015)*cm+3*(3.8+0.015)*cm" y="17*(2.05+0.015)*cm-12*(3.8+0.015)*cm" z="0"/>
+ *         <module sizex="3.8*cm" sizey="3.8*cm" sizez="45*cm" vis="BlueVis" material="PbGlass"/>
+ *         <wrapper thickness="0.015*cm" material="Epoxy" vis="WhiteVis"/>
+ *       </array>
+ *       <array nrow="6" ncol="24" sector="4">
+ *         <position x="17*(2.05+0.015)*cm-12*(3.8+0.015)*cm" y="-17*(2.05+0.015)*cm-3*(3.8+0.015)*cm" z="0"/>
+ *         <module sizex="3.8*cm" sizey="3.8*cm" sizez="45*cm" vis="BlueVis" material="PbGlass"/>
+ *         <wrapper thickness="0.015*cm" material="Epoxy" vis="WhiteVis"/>
+ *       </array>
+ *       <array nrow="24" ncol="6" sector="3">
+ *         <position x="-17*(2.05+0.015)*cm-3*(3.8+0.015)*cm" y="-17*(2.05+0.015)*cm+12*(3.8+0.015)*cm" z="0"/>
+ *         <module sizex="3.8*cm" sizey="3.8*cm" sizez="45*cm" vis="BlueVis" material="PbGlass"/>
+ *         <wrapper thickness="0.015*cm" material="Epoxy" vis="WhiteVis"/>
+ *       </array>
+ *     </placements>
+ *   </detector>
+ *
+ *   <detector id="2" name="SomeBlocks" type="HomogeneousCalorimeter" readout="EcalHits">
+ *     <position x="0" y="0" z="30*cm"/>
+ *     <rotation x="0" y="0" z="0"/>
+ *     <placements>
+ *       <individuals sector="1"/>
+ *         <module sizex="2.05*cm" sizey="2.05*cm" sizez="20*cm" vis="GreenVis" material="PbWO4"/>
+ *         <wrapper thickness="0.015*cm" material="Epoxy" vis="WhiteVis"/>
+ *         <placement x="1*cm" y="1*cm" z="0" id="1"/>
+ *         <placement x="-1*cm" y="1*cm" z="0" id="2"/>
+ *         <placement x="1*cm" y="-1*cm" z="0" id="3"/>
+ *         <placement x="-1*cm" y="-1*cm" z="0" id="4"/>
+ *       </individuals>
+ *     </placements>
+ *   </detector>
+ *
+ *   <detector id="2" name="DiskShapeCalorimeter" type="HomogeneousCalorimeter" readout="EcalHits">
+ *     <position x="0" y="0" z="-30*cm"/>
+ *     <rotation x="0" y="0" z="0"/>
+ *     <placements>
+ *       <disk rmin="25*cm" rmax="125*cm" length="20.5*cm" phimin="0" phimax="360*degree" sector="1"/>
+ *         <module sizex="2.05*cm" sizey="2.05*cm" sizez="20*cm" vis="GreenVis" material="PbWO4"/>
+ *         <wrapper thickness="0.015*cm" material="Epoxy" vis="WhiteVis"/>
+ *     </placements>
+ *   </detector>
+ *
+ *   <detector id="3" name="SomeLines" type="HomogeneousCalorimeter" readout="EcalHits">
+ *     <position x="0" y="0" z="60*cm"/>
+ *     <rotation x="0" y="0" z="0"/>
+ *     <placements>
+ *       <lines sector="1" mirrorx="true" mirrory="true"/>
+ *         <module sizex="2.05*cm" sizey="2.05*cm" sizez="20*cm" vis="GreenVis" material="PbWO4"/>
+ *         <wrapper thickness="0.015*cm" material="Epoxy" vis="WhiteVis"/>
+ *         <line x="10.25*mm" y="10.25*mm" axis="x" begin="8" nmods="16"/>
+ *         <line x="10.25*mm" y="30.75*mm" axis="y" begin="8" nmods="16"/>
+ *         <line x="10.25*mm" y="51.25*mm" axis="z" begin="8" nmods="16"/>
+ *       </individuals>
+ *     </placements>
+ *   </detector>
+ * \endcode
+ *
+ * @{
+ */
+
+// headers
+static std::tuple<int, int> add_individuals(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int id);
+static std::tuple<int, int> add_array(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int id);
+static std::tuple<int, int> add_disk(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int id);
+static std::tuple<int, int> add_lines(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int id);
+
+// helper function to get x, y, z if defined in a xml component
+template<class XmlComp>
+Position get_xml_xyz(XmlComp &comp, dd4hep::xml::Strng_t name)
+{
+    Position pos(0., 0., 0.);
+    if (comp.hasChild(name)) {
+        auto child = comp.child(name);
+        pos.SetX(dd4hep::getAttrOrDefault<double>(child, _Unicode(x), 0.));
+        pos.SetY(dd4hep::getAttrOrDefault<double>(child, _Unicode(y), 0.));
+        pos.SetZ(dd4hep::getAttrOrDefault<double>(child, _Unicode(z), 0.));
+    }
+    return pos;
+}
+
+// main
+static Ref_t create_detector(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);
+    sens.setType("calorimeter");
+    // envelope
+    Assembly assembly(detName);
+
+    // module placement
+    xml::Component plm = detElem.child(_Unicode(placements));
+    std::map<int, int> sectorModuleNumbers;
+    auto addModuleNumbers = [&sectorModuleNumbers] (int sector, int nmod) {
+        auto it = sectorModuleNumbers.find(sector);
+        if (it != sectorModuleNumbers.end()) {
+            it->second += nmod;
+        } else {
+            sectorModuleNumbers[sector] = nmod;
+        }
+    };
+    int sector_id = 1;
+    for (xml::Collection_t mod(plm, _Unicode(individuals)); mod; ++mod) {
+        auto [sector, nmod] = add_individuals(desc, assembly, mod, sens, sector_id++);
+        addModuleNumbers(sector, nmod);
+    }
+    for (xml::Collection_t arr(plm, _Unicode(array)); arr; ++arr) {
+        auto [sector, nmod] = add_array(desc, assembly, arr, sens, sector_id++);
+        addModuleNumbers(sector, nmod);
+    }
+    for (xml::Collection_t disk(plm, _Unicode(disk)); disk; ++disk) {
+        auto [sector, nmod] = add_disk(desc, assembly, disk, sens, sector_id++);
+        addModuleNumbers(sector, nmod);
+    }
+    for (xml::Collection_t lines(plm, _Unicode(lines)); lines; ++lines) {
+        auto [sector, nmod] = add_lines(desc, assembly, lines, sens, sector_id++);
+        addModuleNumbers(sector, nmod);
+    }
+
+    for (auto [sector, nmods] : sectorModuleNumbers) {
+        desc.add(Constant(Form((detName + "_NModules_Sector%d").c_str(), sector), std::to_string(nmods)));
+    }
+
+    // detector position and rotation
+    auto pos = get_xml_xyz(detElem, _Unicode(position));
+    auto rot = get_xml_xyz(detElem, _Unicode(rotation));
+    Volume motherVol = desc.pickMotherVolume(det);
+    Transform3D tr = Translation3D(pos.x(), pos.y(), pos.z()) * RotationZYX(rot.z(), rot.y(), rot.x());
+    PlacedVolume envPV = motherVol.placeVolume(assembly, tr);
+    envPV.addPhysVolID("system", detID);
+    det.setPlacement(envPV);
+    return det;
+}
+
+// helper function to build module with or w/o wrapper
+std::tuple<Volume, Position> build_module(Detector &desc, xml::Collection_t &plm, SensitiveDetector &sens)
+{
+    auto mod = plm.child(_Unicode(module));
+    auto sx = mod.attr<double>(_Unicode(sizex));
+    auto sy = mod.attr<double>(_Unicode(sizey));
+    auto sz = mod.attr<double>(_Unicode(sizez));
+    Box modShape(sx/2., sy/2., sz/2.);
+    auto modMat = desc.material(mod.attr<std::string>(_Unicode(material)));
+    Volume modVol("module_vol", modShape, modMat);
+    modVol.setSensitiveDetector(sens);
+    modVol.setVisAttributes(desc.visAttributes(mod.attr<std::string>(_Unicode(vis))));
+
+    // no wrapper
+    if (!plm.hasChild(_Unicode(wrapper))) {
+        return std::make_tuple(modVol, Position{sx, sy, sz});
+    // build wrapper
+    } else {
+        auto wrp = plm.child(_Unicode(wrapper));
+        auto thickness = wrp.attr<double>(_Unicode(thickness));
+        if (thickness < 1e-12*mm) {
+            return std::make_tuple(modVol, Position{sx, sy, sz});
+        }
+        auto wrpMat = desc.material(wrp.attr<std::string>(_Unicode(material)));
+        Box wrpShape((sx + thickness)/2., (sy + thickness)/2., sz/2.);
+        Volume wrpVol("wrapper_vol", wrpShape, wrpMat);
+        wrpVol.placeVolume(modVol, Position(0., 0., 0.));
+        wrpVol.setVisAttributes(desc.visAttributes(wrp.attr<std::string>(_Unicode(vis))));
+        return std::make_tuple(wrpVol, Position{sx + thickness, sy + thickness, sz});
+    }
+}
+
+// place modules, id must be provided
+static std::tuple<int, int> add_individuals(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
+{
+    auto [modVol, modSize] = build_module(desc, plm, sens);
+    int sector_id = dd4hep::getAttrOrDefault<int>(plm, _Unicode(sector), sid);
+    int nmodules = 0;
+    for (xml::Collection_t pl(plm, _Unicode(placement)); pl; ++pl) {
+        Position pos(dd4hep::getAttrOrDefault<double>(pl, _Unicode(x), 0.),
+                     dd4hep::getAttrOrDefault<double>(pl, _Unicode(y), 0.),
+                     dd4hep::getAttrOrDefault<double>(pl, _Unicode(z), 0.));
+        Position rot(dd4hep::getAttrOrDefault<double>(pl, _Unicode(rotx), 0.),
+                     dd4hep::getAttrOrDefault<double>(pl, _Unicode(roty), 0.),
+                     dd4hep::getAttrOrDefault<double>(pl, _Unicode(rotz), 0.));
+        auto mid = pl.attr<int>(_Unicode(id));
+        Transform3D tr = Translation3D(pos.x(), pos.y(), pos.z())
+                       * RotationZYX(rot.z(), rot.y(), rot.x());
+        auto modPV = env.placeVolume(modVol, tr);
+        modPV.addPhysVolID("sector", sector_id).addPhysVolID("module", mid);
+        nmodules ++;
+    }
+
+    return {sector_id, nmodules};
+}
+
+// place array of modules
+static std::tuple<int, int> add_array(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
+{
+    auto [modVol, modSize] = build_module(desc, plm, sens);
+    int sector_id = dd4hep::getAttrOrDefault<int>(plm, _Unicode(sector), sid);
+    int id_begin = dd4hep::getAttrOrDefault<int>(plm, _Unicode(id_begin), 1);
+    int nrow = plm.attr<int>(_Unicode(nrow));
+    int ncol = plm.attr<int>(_Unicode(ncol));
+
+    // compute array position
+    double begx = -modSize.x()*ncol/2. + modSize.x()/2.;
+    double begy = modSize.y()*nrow/2. - modSize.y()/2.;
+
+    std::vector<std::pair<int, int>> removals;
+    // get the removal list
+    for (xml::Collection_t rm(plm, _Unicode(removal)); rm; ++rm) {
+        removals.push_back({rm.attr<int>(_Unicode(row)), rm.attr<int>(_Unicode(col))});
+    }
+
+    // placement to mother
+    auto pos = get_xml_xyz(plm, _Unicode(position));
+    auto rot = get_xml_xyz(plm, _Unicode(rotation));
+    int nmodules = 0;
+    for (int i = 0; i < nrow; ++i) {
+        for (int j = 0; j < ncol; ++j) {
+            if (std::find(removals.begin(), removals.end(), std::pair<int, int>(i, j)) != removals.end()) {
+                continue;
+            }
+            double px = begx + modSize.x()*j;
+            double py = begy - modSize.y()*i;
+            Transform3D tr = RotationZYX(rot.z(), rot.y(), rot.x())
+                           * Translation3D(pos.x() + px, pos.y() + py, pos.z());
+            auto modPV = env.placeVolume(modVol, tr);
+            modPV.addPhysVolID("sector", sector_id).addPhysVolID("module", i*ncol + j + id_begin);
+            nmodules ++;
+        }
+    }
+    return {sector_id, nmodules};
+}
+
+// place disk of modules
+static std::tuple<int, int> add_disk(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
+{
+    auto [modVol, modSize] = build_module(desc, plm, sens);
+    int sector_id = dd4hep::getAttrOrDefault<int>(plm, _Unicode(sector), sid);
+    int id_begin = dd4hep::getAttrOrDefault<int>(plm, _Unicode(id_begin), 1);
+    double rmin = plm.attr<double>(_Unicode(rmin));
+    double rmax = plm.attr<double>(_Unicode(rmax));
+    double phimin = dd4hep::getAttrOrDefault<double>(plm, _Unicode(phimin), 0.);
+    double phimax = dd4hep::getAttrOrDefault<double>(plm, _Unicode(phimax), 2.*M_PI);
+
+    auto points = athena::geo::fillRectangles({0., 0.}, modSize.x(), modSize.y(), rmin, rmax, phimin, phimax);
+    // placement to mother
+    auto pos = get_xml_xyz(plm, _Unicode(position));
+    auto rot = get_xml_xyz(plm, _Unicode(rotation));
+    int mid = 0;
+    for (auto &p : points) {
+        Transform3D tr = RotationZYX(rot.z(), rot.y(), rot.x())
+                       * Translation3D(pos.x() + p.x(), pos.y() + p.y(), pos.z());
+        auto modPV = env.placeVolume(modVol, tr);
+        modPV.addPhysVolID("sector", sector_id).addPhysVolID("module", id_begin + mid++);
+    }
+    return {sector_id, mid};
+}
+
+// place lines of modules (anchor point is the 0th module of this line)
+static std::tuple<int, int> add_lines(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
+{
+    auto [modVol, modSize] = build_module(desc, plm, sens);
+    int sector_id = dd4hep::getAttrOrDefault<int>(plm, _Unicode(sector), sid);
+    int id_begin = dd4hep::getAttrOrDefault<int>(plm, _Unicode(id_begin), 1);
+    bool mirrorx = dd4hep::getAttrOrDefault<bool>(plm, _Unicode(mirrorx), false);
+    bool mirrory = dd4hep::getAttrOrDefault<bool>(plm, _Unicode(mirrory), false);
+    bool mirrorz = dd4hep::getAttrOrDefault<bool>(plm, _Unicode(mirrorz), false);
+
+    // line placement
+    int mid = 0;
+    for (xml::Collection_t pl(plm, _Unicode(line)); pl; ++pl) {
+        Position pos(dd4hep::getAttrOrDefault<double>(pl, _Unicode(x), 0.),
+                     dd4hep::getAttrOrDefault<double>(pl, _Unicode(y), 0.),
+                     dd4hep::getAttrOrDefault<double>(pl, _Unicode(z), 0.));
+        Position rot(dd4hep::getAttrOrDefault<double>(pl, _Unicode(rotx), 0.),
+                     dd4hep::getAttrOrDefault<double>(pl, _Unicode(roty), 0.),
+                     dd4hep::getAttrOrDefault<double>(pl, _Unicode(rotz), 0.));
+        // determine axis
+        std::string axis = dd4hep::getAttrOrDefault(pl, _Unicode(axis), "x");
+        std::transform(axis.begin(), axis.end(), axis.begin(), [](char c) { return std::tolower(c); });
+        if ((axis != "x") && (axis != "y") && (axis != "z")) {
+            std::cerr << "HomogeneousCalorimeter Error: cannot determine axis of line " << axis
+                      << ", abort placement of this line." << std::endl;
+            continue;
+        }
+
+        int begin = dd4hep::getAttrOrDefault<int>(pl, _Unicode(begin), 0);
+        int nmods = pl.attr<int>(_Unicode(nmods));
+
+        std::vector<Position> trans;
+        for (int i = 0; i < nmods; ++i) {
+            Position tran{ (axis == "x") ? pos.x() + (begin + i)*modSize.x() : pos.x(),
+                           (axis == "y") ? pos.y() + (begin + i)*modSize.y() : pos.y(),
+                           (axis == "z") ? pos.z() + (begin + i)*modSize.z() : pos.z() };
+            trans.push_back(tran);
+        }
+
+        // mirror placement
+        if (mirrorx) {
+            int curr_size = trans.size();
+            for (size_t i = 0; i < curr_size; ++i) {
+                trans.push_back(Position{-trans[i].x(), trans[i].y(), trans[i].z()});
+            }
+        }
+        if (mirrory) {
+            int curr_size = trans.size();
+            for (size_t i = 0; i < curr_size; ++i) {
+                trans.push_back(Position{trans[i].x(), -trans[i].y(), trans[i].z()});
+            }
+        }
+        if (mirrorz) {
+            int curr_size = trans.size();
+            for (size_t i = 0; i < curr_size; ++i) {
+                trans.push_back(Position{trans[i].x(), trans[i].y(), -trans[i].z()});
+            }
+        }
+
+        // place volume
+        for (auto &p : trans) {
+            Transform3D tr = RotationZYX(rot.z(), rot.y(), rot.x())
+                           * Translation3D(p.x(), p.y(), p.z());
+            auto modPV = env.placeVolume(modVol, tr);
+            modPV.addPhysVolID("sector", sector_id).addPhysVolID("module", id_begin + mid++);
+        }
+    }
+    return {sector_id, mid};
+}
+//@}
+DECLARE_DETELEMENT(HomogeneousCalorimeter, create_detector)
+

src/HybridCalorimeter_geo.cpp

+//==========================================================================
+//  A general implementation for homogeneous calorimeter
+//  it supports three types of placements
+//  1. Individual module placement with module dimensions and positions
+//  2. Array placement with module dimensions and numbers of row and columns
+//  3. Disk placement with module dimensions and (Rmin, Rmax), and (Phimin, Phimax)
+//  4. Lines placement with module dimensions and (mirrorx, mirrory)
+//     (NOTE: anchor point is the 0th block of the line instead of line center)
+//--------------------------------------------------------------------------
+//  Author: Chao Peng (ANL)
+//  Date: 06/09/2021
+//==========================================================================
+
+#include "GeometryHelpers.h"
+#include "DD4hep/DetFactoryHelper.h"
+#include <XML/Helper.h>
+#include <iostream>
+#include <algorithm>
+#include <tuple>
+#include <math.h>
+#include <fmt/core.h>
+
+using namespace dd4hep;
+
+// main
+static Ref_t create_detector(Detector& desc, xml::Handle_t handle, SensitiveDetector sens)
+{
+
+  using namespace std;
+  using namespace fmt;
+
+    xml::DetElement detElem = handle;
+    std::string detName = detElem.nameStr();
+    int detID = detElem.id();
+    DetElement det(detName, detID);
+    sens.setType("calorimeter");
+
+    auto glass_material = desc.material("SciGlass");
+    auto crystal_material = desc.material("PbWO4");
+    auto air_material = desc.material("Air");
+
+    double ROut = desc.constantAsDouble("EcalEndcapN_rmax");
+    double RIn_el = desc.constantAsDouble("EcalEndcapN_rmin1");
+    double ionCutout_dphi = desc.constantAsDouble("EcalEndcapNIonCutout_dphi");
+    double RIn = desc.constantAsDouble("EcalEndcapN_rmin2");
+    double SizeZ = desc.constantAsDouble("EcalEndcapN_thickness");
+    double thickness = desc.constantAsDouble("EcalEndcapN_thickness");
+    double trans_radius = desc.constantAsDouble("EcalEndcapNCrystal_rmax");
+    double Glass_z0 = desc.constantAsDouble("GlassModule_z0");
+    double Glass_Width = desc.constantAsDouble("GlassModule_width");
+    double Glass_thickness = desc.constantAsDouble("GlassModule_length");
+    double Glass_Gap = desc.constantAsDouble("GlassModule_wrap");
+    double glass_distance = desc.constantAsDouble("GlassModule_distance");
+
+    double Crystal_Width = desc.constantAsDouble("CrystalModule_width");
+    double Crystal_thickness = desc.constantAsDouble("CrystalModule_length");
+    double Crystal_Gap = desc.constantAsDouble("CrystalModule_wrap");
+    double crystal_distance = desc.constantAsDouble("CrystalModule_distance");
+    double Crystal_z0 = desc.constantAsDouble("CrystalModule_z0");
+
+    // RIn and ROut will define outer tube embedding the calorimeter
+    // centers_rmin/out define the maximum radius of module centers
+    // so that modules are not overlapping with mother tube volume
+    const double glassHypotenuse = std::hypot(glass_distance, glass_distance)/2;
+    const double crystalHypotenuse = glassHypotenuse/2;
+    // Offset these values a bit so we don't miss edge-blocks
+    const double glassCenters_rmax  = ROut - glassHypotenuse + 1 * mm;
+    const double crystalCenters_rmin = RIn + crystalHypotenuse - 1 * mm;
+    // Also limits of the inner crystal blocks fill
+    const double cutout_tan = tan(ionCutout_dphi/2);
+    const double cutout_rmin = RIn_el + crystalHypotenuse - 1 * mm;
+
+    // Offset to align the modules at the zmin of the endcap,
+    const double Crystal_offset = -0.5 * (Crystal_thickness - thickness);
+    const double Glass_offset = -0.5*(Glass_thickness - thickness);
+
+    // envelope
+    // consists of an glass tube of the full thickness, and a crystal inner tube
+    // for the crystal that allows us to get closet to the beampipe without
+    // overlaps, and a partial electron tube that allows us to get closer to the
+    // electron beampipe in the region where there is no ion beampipe
+    Tube glass_tube(min(RIn + glassHypotenuse*2, trans_radius), ROut, SizeZ / 2.0, 0., 360.0 * deg);
+    Tube crystal_tube(RIn, min(RIn + glassHypotenuse*2, trans_radius), Crystal_thickness/ 2.0, 0., 360.0 * deg);
+    Tube electron_tube(RIn_el, RIn, Crystal_thickness / 2., ionCutout_dphi / 2., 360.0 * deg - ionCutout_dphi / 2);
+    UnionSolid outer_envelope(glass_tube, crystal_tube, Position(0,0,Crystal_offset));
+    UnionSolid envelope(outer_envelope, electron_tube, Position(0,0,Crystal_offset));
+    Volume ecal_vol("negative_ecal", envelope, air_material);
+    ecal_vol.setVisAttributes(desc.visAttributes("HybridEcalOuterVis"));
+
+    // TODO why 1cm and not something else?
+    double Glass_OuterR = ROut - 1 * cm ;
+    double Glass_InnerR = trans_radius;
+
+    // Geometry of modules
+    Box glass_box("glass_box", Glass_Width * 0.5, Glass_Width * 0.5, Glass_thickness * 0.5);
+    Volume glass_module("glass_module", glass_box, glass_material);
+    glass_module.setVisAttributes(desc.visAttributes("EcalEndcapNModuleVis"));
+    glass_module.setSensitiveDetector(sens);
+    
+    Box crystal_box("crystal_box",  Crystal_Width* 0.5, Crystal_Width * 0.5, Crystal_thickness * 0.5);
+    Volume crystal_module("crystal_module", crystal_box, crystal_material);
+    crystal_module.setVisAttributes(desc.visAttributes("EcalEndcapNModuleVis"));
+    crystal_module.setSensitiveDetector(sens);
+
+    // GLASS
+    double diameter = 2 * Glass_OuterR;
+
+    // Can we fit an even or odd amount of glass blocks within our rmax?
+    // This determines the transition points between crystal and glass as we need the
+    // outer crystal arangement to be in multiples of 2 (aligned with glass)
+    const int towersInRow = std::ceil((diameter + Glass_Gap) /  (Glass_Width + Glass_Gap));
+    const bool align_even = (towersInRow % 2);
+
+    // Is it odd or even number of towersInRow
+    double leftTowerPos, topTowerPos;
+    if(towersInRow%2) {
+      //             |
+      //      [ ][ ][ ][ ][ ]
+      //       ^     |
+      int towersInHalfRow = std::ceil(towersInRow/2.0);
+      topTowerPos = leftTowerPos = -towersInHalfRow * (Glass_Width + Glass_Gap);
+
+    } else {
+      //               |
+      //      [ ][ ][ ][ ][ ][ ]
+      //       ^      |
+      int towersInHalfRow = towersInRow/2;
+      topTowerPos = leftTowerPos = -(towersInHalfRow - 0.5) * (Glass_Width + Glass_Gap);
+    }
+
+    int moduleIndex = 0;
+
+    int glass_module_index = 0;
+    int cryst_module_index = 0;
+    for(int rowIndex=0; rowIndex < towersInRow; rowIndex++) {
+      for(int colIndex=0; colIndex < towersInRow; colIndex++) {
+        const double glass_x = leftTowerPos + colIndex * (Glass_Width + Glass_Gap);
+        const double glass_y = topTowerPos + rowIndex * (Glass_Width + Glass_Gap);
+        const double r = std::hypot(glass_x, glass_y);
+        // crystal if within the transition radius (as defined by the equivalent glass
+        // block)
+        if (r < trans_radius) {
+          for (const auto dx : {-1, 1}) {
+            for (const auto& dy : {-1, 1}) {
+              const double crystal_x = glass_x + dx * crystal_distance / 2;
+              const double crystal_y = glass_y + dy * crystal_distance / 2;
+              const double crystal_r = std::hypot(crystal_x, crystal_y);
+              // check if crystal in the main crystal ring?
+              const bool mainRing = (crystal_r > crystalCenters_rmin);
+              const bool innerRing = !mainRing && crystal_r > cutout_rmin;
+              const bool ionCutout = crystal_x > 0 && fabs(crystal_y / crystal_x) < cutout_tan;
+              if (mainRing || (innerRing && !ionCutout)) {
+                auto placement =
+                    ecal_vol.placeVolume(crystal_module, Position(crystal_x, crystal_y, Crystal_z0 + Crystal_offset));
+                placement.addPhysVolID("sector", 1);
+                placement.addPhysVolID("module", cryst_module_index++);
+              }
+            }
+          }
+        // Glass block if within the rmax
+        } else if (r < glassCenters_rmax) {
+          // glass module
+          auto placement = ecal_vol.placeVolume(glass_module, Position(glass_x, glass_y, Glass_z0 + Glass_offset));
+          placement.addPhysVolID("sector", 2);
+          placement.addPhysVolID("module", glass_module_index++);
+        }
+      }
+    }
+
+    desc.add(Constant("EcalEndcapN_NModules_Sector1", std::to_string(cryst_module_index)));
+    desc.add(Constant("EcalEndcapN_NModules_Sector2", std::to_string(glass_module_index)));
+//    fmt::print("Total Glass modules: {}\n", towerIndex);
+//    fmt::print("CE EMCAL GLASS END\n\n");
+
+     // detector position and rotation
+     auto pos = detElem.position();
+     auto rot = detElem.rotation();
+     Volume motherVol = desc.pickMotherVolume(det);
+     Transform3D tr(RotationZYX(rot.z(), rot.y(), rot.x()), Position(pos.x(), pos.y(), pos.z()));
+     PlacedVolume envPV = motherVol.placeVolume(ecal_vol, tr);
+     envPV.addPhysVolID("system", detID);
+     det.setPlacement(envPV);
+     return det;
+
+}
+
+//@}
+DECLARE_DETELEMENT(HybridCalorimeter, create_detector)
+

src/MRich_geo.cpp

+//
+// Author     : Whit Armstrong (warmstrong@anl.gov)
+//
+#include <XML/Helper.h>
+#include "TMath.h"
+#include "TString.h"
+#include "DDRec/Surface.h"
+#include "DDRec/DetectorData.h"
+#include "DD4hep/OpticalSurfaces.h"
+#include "DD4hep/DetFactoryHelper.h"
+#include "DD4hep/Printout.h"
+#include "GeometryHelpers.h"
+#include "Math/Vector3D.h"
+#include "Math/AxisAngle.h"
+#include "Math/VectorUtil.h"
+
+using namespace std;
+using namespace dd4hep;
+using namespace dd4hep::rec;
+
+using Placements = vector<PlacedVolume>;
+
+
+
+static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDetector sens){
+  xml_det_t      x_det    = e;
+  Material       air      = description.material("AirOptical");
+  Material       vacuum   = description.vacuum();
+  string         det_name = x_det.nameStr();
+  DetElement     sdet(det_name, x_det.id());
+  Assembly       assembly(det_name);
+  sens.setType("tracker");
+  OpticalSurfaceManager surfMgr = description.surfaceManager();
+
+  bool projective = getAttrOrDefault(x_det, _Unicode(projective), false);
+  bool reflect    = x_det.reflect(true);
+
+  PlacedVolume pv;
+
+  map<string, Volume>     modules;
+  map<string, Placements> sensitives;
+  map<string, Volume>     module_assemblies;
+  std::map<std::string,DetElement> module_assembly_delements;
+
+  int                     n_sensor = 1;
+
+  // dimensions
+  xml::Component dims   = x_det.dimensions();
+  auto           rmin   = dims.rmin();
+  auto           rmax   = dims.rmax();
+  auto           length = dims.length();
+  auto           zmin   = dims.zmin();
+  auto           zpos   = zmin + length / 2;
+
+  // envelope
+  Tube   envShape(rmin, rmax, length / 2., 0., 2 * M_PI);
+  Volume envVol("MRICH_Envelope", envShape, air);
+  envVol.setVisAttributes(description.visAttributes(x_det.visStr()));
+  if (x_det.hasChild(_Unicode(envelope))) {
+    xml_comp_t x_envelope = x_det.child(_Unicode(envelope));
+    double thickness  = x_envelope.thickness();
+    Material material = description.material(x_envelope.materialStr());
+    Tube   envInsideShape(rmin + thickness, rmax - thickness, length / 2. - thickness);
+    SubtractionSolid envShellShape(envShape, envInsideShape);
+    Volume envShell("MRICH_Envelope_Inside", envShellShape, material);
+    envVol.placeVolume(envShell);
+  }
+
+  // expect only one module (for now)
+  xml_comp_t x_mod = x_det.child(_U(module));
+  string     mod_name            = x_mod.nameStr();
+  double     mod_width           = getAttrOrDefault(x_mod, _U(width), 130.0 * mm);
+  double     mod_height          = getAttrOrDefault(x_mod, _U(height), 130.0 * mm);
+  double     mod_length          = getAttrOrDefault(x_mod, _U(length), 130.0 * mm);
+
+  // module
+  Box    m_solid(mod_width / 2.0, mod_height / 2.0, mod_length / 2.0);
+  Volume m_volume(mod_name, m_solid, air);
+  m_volume.setVisAttributes(description.visAttributes(x_mod.visStr()));
+  DetElement mod_de( mod_name + std::string("_mod_") + std::to_string(1), 1);
+  double z_placement = - mod_length / 2.0;
+
+  // todo module frame
+  if (x_mod.hasChild(_Unicode(frame))) {
+    xml_comp_t x_frame    = x_mod.child(_Unicode(frame));
+    double     frame_thickness = getAttrOrDefault(x_frame, _U(thickness), 2.0 * mm);
+    Box        frame_inside(mod_width / 2.0 - frame_thickness, mod_height / 2.0 - frame_thickness, mod_length / 2.0 - frame_thickness);
+    SubtractionSolid frame_solid(m_solid, frame_inside);
+    Material   frame_mat       = description.material(x_frame.materialStr());
+    Volume     frame_vol(mod_name+"_frame", frame_solid, frame_mat);
+    auto       frame_vis       = getAttrOrDefault<std::string>(x_frame, _U(vis), std::string("GrayVis"));
+    frame_vol.setVisAttributes(description.visAttributes(frame_vis));
+    // update position
+    z_placement += frame_thickness / 2.0;
+    // place volume
+    m_volume.placeVolume(frame_vol);
+    // update position
+    z_placement += frame_thickness / 2.0;
+  }
+
+  // aerogel box
+  if (x_mod.hasChild(_Unicode(aerogel))) {
+    xml_comp_t x_aerogel  = x_mod.child(_Unicode(aerogel));
+    double     aerogel_width   = getAttrOrDefault(x_aerogel, _U(width), 130.0 * mm);
+    double     aerogel_length  = getAttrOrDefault(x_aerogel, _U(length), 130.0 * mm);
+    Material   aerogel_mat     = description.material(x_aerogel.materialStr());
+    auto       aerogel_vis     = getAttrOrDefault<std::string>(x_aerogel, _U(vis), std::string("InvisibleWithDaughters"));
+
+    xml_comp_t x_aerogel_frame = x_aerogel.child(_Unicode(frame));
+    double     foam_thickness  = getAttrOrDefault(x_aerogel_frame, _U(thickness), 2.0 * mm);
+    Material   foam_mat        = description.material(x_aerogel_frame.materialStr());
+    auto       foam_vis        = getAttrOrDefault<std::string>(x_aerogel_frame, _U(vis), std::string("RedVis"));
+
+    // foam frame
+    Box foam_box(aerogel_width / 2.0 + foam_thickness, aerogel_width / 2.0 + foam_thickness, (aerogel_length + foam_thickness) / 2.0);
+    Box foam_sub_box(aerogel_width / 2.0, aerogel_width / 2.0, (aerogel_length + foam_thickness) / 2.0);
+    SubtractionSolid foam_frame_solid(foam_box, foam_sub_box, Position(0, 0, foam_thickness));
+    Volume           foam_vol(mod_name+"_aerogel_frame", foam_frame_solid, foam_mat);
+    foam_vol.setVisAttributes(description.visAttributes(foam_vis));
+
+    // aerogel
+    Box              aerogel_box(aerogel_width / 2.0, aerogel_width / 2.0, (aerogel_length) / 2.0);
+    Volume           aerogel_vol(mod_name+"_aerogel", aerogel_box, aerogel_mat);
+    aerogel_vol.setVisAttributes(description.visAttributes(aerogel_vis));
+
+    // update position
+    z_placement += (aerogel_length + foam_thickness) / 2.0;
+    // place foam frame
+    pv = m_volume.placeVolume(foam_vol,Position(0,0,z_placement));
+    // place aerogel
+    z_placement += foam_thickness / 2.0;
+    pv = m_volume.placeVolume(aerogel_vol,Position(0,0,z_placement));
+    DetElement aerogel_de(mod_de, mod_name + std::string("_aerogel_de") + std::to_string(1), 1);
+    aerogel_de.setPlacement(pv);
+    // update position
+    z_placement += aerogel_length / 2.0;
+
+    // optical surfaces
+    auto aerogel_surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault(x_aerogel, _Unicode(surface), "MRICH_AerogelOpticalSurface"));
+    SkinSurface skin_surf(description, aerogel_de, Form("MRICH_aerogel_skin_surface_%d", 1), aerogel_surf, aerogel_vol);
+    skin_surf.isValid();
+  }
+
+  // Fresnel Lens
+  if (x_mod.hasChild(_Unicode(lens))) {
+    xml_comp_t x_lens     = x_mod.child(_Unicode(lens));
+
+    //  - The lens has a constant groove pitch (delta r) as opposed to fixing the groove height.
+    //  - The lens area outside of the effective diamtere is flat.
+    //  - The grooves are not curved, rather they are polycone shaped, ie a flat approximating the curvature.
+    auto   lens_vis       = getAttrOrDefault<std::string>(x_lens, _U(vis), std::string("AnlBlue"));
+    double groove_pitch   = getAttrOrDefault(x_lens, _Unicode(pitch), 0.2 * mm);// 0.5 * mm);
+    double lens_f         = getAttrOrDefault(x_lens, _Unicode(focal_length), 6.0*2.54*cm);
+    double eff_diameter   = getAttrOrDefault(x_lens, _Unicode(effective_diameter), 152.4 * mm);
+    double lens_width     = getAttrOrDefault(x_lens, _Unicode(width), 6.7*2.54*cm);
+    double center_thickness = getAttrOrDefault(x_lens, _U(thickness), 0.068 * 2.54 * cm);//2.0 * mm);
+
+    double n_acrylic        = 1.49;
+    double lens_curvature   = 1.0 / (lens_f*(n_acrylic - 1.0)); //confirmed
+    double full_ring_rmax   = std::min(eff_diameter / 2.0, lens_width/2.0);
+
+    double N_grooves        = std::ceil((full_ring_rmax) / groove_pitch);
+    double groove_last_rmin = (N_grooves - 1) * groove_pitch;
+    double groove_last_rmax = N_grooves * groove_pitch;
+
+    auto   groove_sagitta = [&](double r) { return lens_curvature * std::pow(r, 2) / (1.0 + 1.0); };
+    double lens_thickness = groove_sagitta(groove_last_rmax) - groove_sagitta(groove_last_rmin) + center_thickness;
+
+    Material         lens_mat = description.material(x_lens.materialStr());
+    Box              lens_box(lens_width / 2.0, lens_width / 2.0, (center_thickness) / 2.0);
+    SubtractionSolid flat_lens(lens_box, Tube(0.0, full_ring_rmax, 2 * center_thickness));
+
+    Assembly lens_vol(mod_name + "_lens");
+    Volume   flatpart_lens_vol( "flatpart_lens", flat_lens, lens_mat);
+    lens_vol.placeVolume(flatpart_lens_vol);
+
+    int    i_groove           = 0;
+    double groove_rmax        = groove_pitch;
+    double groove_rmin        = 0;
+
+    while ( groove_rmax <= full_ring_rmax ) {
+      double   dZ = groove_sagitta(groove_rmax) - groove_sagitta(groove_rmin);
+      Polycone groove_solid(0, 2.0 * M_PI,
+                            {groove_rmin, groove_rmin, groove_rmin},
+                            {groove_rmax, groove_rmax, groove_rmin},
+                            {-lens_thickness/2.0, lens_thickness/2.0-dZ, lens_thickness/2.0});
+      Volume   lens_groove_vol("lens_groove_" + std::to_string(i_groove), groove_solid, lens_mat);
+      lens_vol.placeVolume(lens_groove_vol);
+
+      i_groove++;
+      groove_rmin = (i_groove  )*groove_pitch;
+      groove_rmax = (i_groove+1)*groove_pitch;
+    }
+
+    lens_vol.setVisAttributes(description.visAttributes(lens_vis));
+
+    // update position
+    z_placement += lens_thickness/2.0;
+    // place volume
+    pv = m_volume.placeVolume(lens_vol,Position(0,0,z_placement));
+    DetElement lens_de(mod_de, mod_name + std::string("_lens_de") + std::to_string(1), 1);
+    lens_de.setPlacement(pv);
+    // update position
+    z_placement += lens_thickness/2.0;
+
+    // optical surfaces
+    auto lens_surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault(x_lens, _Unicode(surface), "MRICH_LensOpticalSurface"));
+    SkinSurface skin_surf(description, lens_de, Form("MRichFresnelLens_skin_surface_%d", 1), lens_surf, lens_vol);
+    skin_surf.isValid();
+  }
+
+  // mirror
+  if (x_mod.hasChild(_Unicode(space))) {
+    xml_comp_t x_space = x_mod.child(_Unicode(space));
+    z_placement += getAttrOrDefault(x_space, _U(thickness), 0.0 * mm);
+  }
+
+  // mirror
+  if (x_mod.hasChild(_Unicode(mirror))) {
+    xml_comp_t x_mirror   = x_mod.child(_Unicode(mirror));
+    auto   mirror_vis     = getAttrOrDefault<std::string>(x_mirror, _U(vis), std::string("AnlGray"));
+    double mirror_x1      = getAttrOrDefault(x_mirror, _U(x1), 100.0 * mm);
+    double mirror_x2      = getAttrOrDefault(x_mirror, _U(x2), 80.0 * mm);
+    double mirror_length  = getAttrOrDefault(x_mirror, _U(length), 130.0 * mm);
+    double mirror_thickness  = getAttrOrDefault(x_mirror, _U(thickness), 2.0 * mm);
+    double outer_x1 = (mirror_x1+mirror_thickness)/2.0;
+    double outer_x2 = (mirror_x2+mirror_thickness)/2.0;
+    Trd2   outer_mirror_trd(outer_x1, outer_x2, outer_x1,  outer_x2, mirror_length/2.0);
+    Trd2   inner_mirror_trd(mirror_x1 / 2.0,  mirror_x2 / 2.0, mirror_x1 / 2.0,mirror_x2 / 2.0, mirror_length/2.0+0.1*mm);
+    SubtractionSolid mirror_solid(outer_mirror_trd, inner_mirror_trd);
+    Material mirror_mat        = description.material(x_mirror.materialStr());
+    Volume   mirror_vol(mod_name+"_mirror", mirror_solid, mirror_mat);
+
+    // update position
+    z_placement += mirror_length/2.0;
+    // place volume
+    pv = m_volume.placeVolume(mirror_vol,Position(0,0,z_placement));
+    DetElement mirror_de(mod_de, mod_name + std::string("_mirror_de") + std::to_string(1), 1);
+    mirror_de.setPlacement(pv);
+    // update position
+    z_placement += mirror_length/2.0;
+
+    // optical surfaces
+    auto mirror_surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault(x_mirror, _Unicode(surface), "MRICH_MirrorOpticalSurface"));
+    SkinSurface skin_surf(description, mirror_de, Form("MRICH_mirror_skin_surface_%d", 1), mirror_surf, mirror_vol);
+    skin_surf.isValid();
+  }
+
+  // photon detector
+  if (x_mod.hasChild(_Unicode(photodet))) {
+    xml_comp_t x_photodet = x_mod.child(_Unicode(photodet));
+    auto       photodet_vis       = getAttrOrDefault<std::string>(x_photodet, _U(vis), std::string("AnlRed"));
+    double     photodet_width     = getAttrOrDefault(x_photodet, _U(width), 130.0 * mm);
+    double     photodet_thickness = getAttrOrDefault(x_photodet, _U(thickness), 2.0 * mm);
+    Material   photodet_mat       = description.material(x_photodet.materialStr());
+    Box        window_box(photodet_width/2.0,photodet_width/2.0,photodet_thickness/2.0);
+    Volume     window_vol(mod_name+"_window", window_box, photodet_mat);
+
+    // update position
+    z_placement += photodet_thickness/2.0;
+    // place volume
+    pv = m_volume.placeVolume(window_vol,Position(0,0,z_placement));
+    DetElement   comp_de(mod_de, mod_name + std::string("_sensor_de_") + std::to_string(1) ,  1);
+    comp_de.setPlacement(pv);
+    // update position
+    z_placement += photodet_thickness/2.0;
+
+    // sensitive
+    pv.addPhysVolID("sensor", n_sensor);
+    window_vol.setSensitiveDetector(sens);
+    sensitives[mod_name].push_back(pv);
+    ++n_sensor;
+
+    // sensor
+    xml_comp_t x_sensor  = x_photodet.child(_Unicode(sensor));
+    double     sensor_thickness   = getAttrOrDefault(x_sensor, _U(thickness), 2.0 * mm);
+    Material   sensor_mat         = description.material(x_sensor.materialStr());
+    int        sensor_nx          = getAttrOrDefault(x_sensor, _Unicode(nx), 2);
+    int        sensor_ny          = getAttrOrDefault(x_sensor, _Unicode(ny), 2);
+
+    // layers
+    int i_layer = 1;
+    for (xml_coll_t li(x_photodet, _Unicode(layer)); li; ++li) {
+      xml_comp_t x_layer = li;
+      Material layer_mat = description.material(x_layer.materialStr());
+      double   layer_thickness = x_layer.thickness();
+      Box      layer_box(photodet_width/2.0,photodet_width/2.0,layer_thickness/2.0);
+      Volume   layer_vol(mod_name + "_layer_" + std::to_string(i_layer), layer_box, layer_mat);
+
+      // update position
+      z_placement += layer_thickness / 2.0;
+      // place volume
+      pv = m_volume.placeVolume(layer_vol,Position(0,0,z_placement));
+      DetElement layer_de(mod_de, mod_name + std::string("_layer_de_") + std::to_string(i_layer),  1);
+      layer_de.setPlacement(pv);
+      // update position
+      z_placement += layer_thickness / 2.0;
+
+      i_layer++;
+    }
+  }
+
+  //for (size_t ic = 0; ic < sensVols.size(); ++ic) {
+  //  PlacedVolume sens_pv = sensVols[ic];
+  //  DetElement   comp_de(mod_de, std::string("de_") + sens_pv.volume().name(), ic + 1);
+  //  comp_de.setPlacement(sens_pv);
+  //  // 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());
+  //}
+  //DetElement window_de(sdet, mod_name + std::string("_window_de") + std::to_string(1), 1);
+  //window_de.setPlacement(pv);
+
+  modules[mod_name]                   = m_volume;
+  module_assembly_delements[mod_name] = mod_de;
+  // end module
+
+  // place modules in the sectors (disk)
+  auto points = athena::geo::fillSquares({0., 0.}, mod_width, rmin, rmax);
+
+  // mod_name = ...
+  Placements& sensVols = sensitives[mod_name];
+  auto        mod_v    = modules[mod_name];
+  // determine module direction, always facing z = 0
+  double roty = dims.zmin() < 0. ? -M_PI : 0 ;
+
+  // read module positions
+  std::vector<std::tuple<double,double,double>> positions;
+  for (xml_coll_t x_positions_i(x_det, _Unicode(positions)); x_positions_i; ++x_positions_i) {
+    xml_comp_t x_positions = x_positions_i;
+    for (xml_coll_t x_position_i(x_positions, _U(position)); x_position_i; ++x_position_i) {
+      xml_comp_t x_position = x_position_i;
+      positions.push_back(
+        std::make_tuple(x_positions.scale() * x_position.x() * mm,
+                        x_positions.scale() * x_position.y() * mm,
+                        -x_positions.z0()));
+    }
+  }
+  // if no positions, then autoplacement
+  if (positions.empty()) {
+    for (double x = mod_width / 2.0; x < rmax - mod_width / 2.0; x += mod_width) {
+      for (double y = mod_width / 2.0; y < rmax - mod_width / 2.0; y += mod_width) {
+        if (pow(x + mod_width / 2.0,2) + pow(y + mod_width / 2.0,2) > rmax*rmax) continue;
+        if (pow(x - mod_width / 2.0,2) + pow(y - mod_width / 2.0,2) < rmin*rmin) continue;
+        positions.push_back(std::make_tuple(x, y, 0));
+      }
+    }
+  }
+
+  // place modules
+  int i_mod = 1; // starts at 1
+  for (auto& p: positions) {
+
+    // get positions in one quadrant
+    double x = std::get<0>(p);
+    double y = std::get<1>(p);
+    double z0 = std::get<2>(p);
+
+    // and place in all quadrants (intentional shadowing)
+    for (auto& p: decltype(positions){{x,y,z0}, {y,-x,z0}, {-x,-y,z0}, {-y,x,z0}}) {
+
+      // get positions (intentional shadowing)
+      double x = std::get<0>(p);
+      double y = std::get<1>(p);
+      double z0 = std::get<2>(p);
+
+      Transform3D tr;
+      if(projective) {
+        double rotAngX = atan(y/z0);
+        double rotAngY = -1.*atan(x/z0);
+        tr = Translation3D(x, y, 0) * RotationX(rotAngX) * RotationY(rotAngY);
+      } else {
+        tr = Translation3D(x, y, 0) * RotationX(0);
+      }
+
+      // mod placement
+      pv = envVol.placeVolume(mod_v, tr);
+      pv.addPhysVolID("module", i_mod);
+
+      auto mod_det_element  = module_assembly_delements[mod_name].clone(mod_name + "__" + std::to_string(i_mod));
+      mod_det_element.setPlacement(pv);
+      sdet.add(mod_det_element);
+
+      i_mod++;
+    }
+  }
+
+  // additional layers
+  if (x_det.hasChild(_Unicode(layer))) {
+    xml_comp_t x_layer = x_det.child(_Unicode(layer));
+    double   layer_thickness = x_layer.thickness();
+    Material layer_mat = description.material(x_layer.materialStr());
+    Tube   frameShape(rmin, rmax, layer_thickness / 2., 0., 2 * M_PI);
+    Volume frameVol("MRICH_Frame", frameShape, layer_mat);
+    pv = envVol.placeVolume(frameVol, Position(0, 0, (length - layer_thickness) / 2.0));
+  }
+
+  // place envelope
+  Volume motherVol = description.pickMotherVolume(sdet);
+  if (reflect) {
+    pv = motherVol.placeVolume(envVol, Transform3D(RotationZYX(0, M_PI, 0), Position(0, 0, -zpos)));
+  } else {
+    pv = motherVol.placeVolume(envVol, Transform3D(RotationZYX(0,    0, 0), Position(0, 0, +zpos)));
+  }
+  pv.addPhysVolID("system", x_det.id());
+  sdet.setPlacement(pv);
+  return sdet;
+}
+
+
+//void addModules(Volume &mother, xml::DetElement &detElem, Detector &description, SensitiveDetector &sens)
+//{
+//    xml::Component dims = detElem.dimensions();
+//    xml::Component mods = detElem.child(_Unicode(modules));
+//
+//    auto rmin = dims.rmin();
+//    auto rmax = dims.rmax();
+//
+//    auto mThick = mods.attr<double>(_Unicode(thickness));
+//    auto mWidth = mods.attr<double>(_Unicode(width));
+//    auto mGap = mods.attr<double>(_Unicode(gap));
+//
+//    auto modMat = description.material(mods.materialStr());
+//    auto gasMat = description.material("AirOptical");
+//
+//    // single module
+//    Box mShape(mWidth/2., mWidth/2., mThick/2. - 0.1*mm);
+//    Volume mVol("ce_MRICH_mod_Solid", mShape, modMat);
+//
+//    // a thin gas layer to detect optical photons
+//    Box modShape(mWidth/2., mWidth/2., mThick/2.);
+//    Volume modVol("ce_MRICH_mod_Solid_v", modShape, gasMat);
+//    // thin gas layer is on top (+z) of the material
+//    modVol.placeVolume(mVol, Position(0., 0., -0.1*mm));
+//
+//    modVol.setVisAttributes(description.visAttributes(mods.visStr()));
+//    sens.setType("tracker");
+//    modVol.setSensitiveDetector(sens);
+//
+//    // place modules in the sectors (disk)
+//    auto points = ref::utils::fillSquares({0., 0.}, mWidth + mGap, rmin - mGap, rmax + mGap);
+//
+//    // determine module direction, always facing z = 0
+//    double roty = dims.z() > 0. ? M_PI/2. : -M_PI/2.;
+//    int imod = 1;
+//    for (auto &p : points) {
+//        // operations are inversely ordered
+//        Transform3D tr = Translation3D(p.x(), p.y(), 0.)        // move to position
+//                       * RotationY(roty);                       // facing z = 0.
+//        auto modPV = mother.placeVolume(modVol, tr);
+//        modPV.addPhysVolID("sector", 1).addPhysVolID("module", imod ++);
+//    }
+//}
+
+// clang-format off
+DECLARE_DETELEMENT(athena_MRICH, createDetector)
+

src/PFRICH_geo.cpp

+//----------------------------------
+//  pfRICH: Proximity Focusing RICH
+//  Author: C. Dilks
+//----------------------------------
+
+#include "DD4hep/DetFactoryHelper.h"
+#include "DD4hep/OpticalSurfaces.h"
+#include "DD4hep/Printout.h"
+#include "DDRec/DetectorData.h"
+#include "DDRec/Surface.h"
+#include "GeometryHelpers.h"
+#include "Math/Point2D.h"
+#include "TMath.h"
+#include "TString.h"
+#include <XML/Helper.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 vesselLength    = dims.attr<double>(_Unicode(length));
+  double vesselZmin      = dims.attr<double>(_Unicode(zmin));
+  double vesselZmax      = dims.attr<double>(_Unicode(zmax));
+  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 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));
+  // - 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 plane
+  auto   sensorPlaneElem = detElem.child(_Unicode(sensors)).child(_Unicode(plane));
+  double sensorPlaneDist = sensorPlaneElem.attr<double>(_Unicode(sensordist));
+  double sensorPlaneRmin = sensorPlaneElem.attr<double>(_Unicode(rmin));
+  double sensorPlaneRmax = sensorPlaneElem.attr<double>(_Unicode(rmax));
+  // - debugging switches
+  int debug_optics_mode = detElem.attr<int>(_Unicode(debug_optics));
+
+  // if debugging optics, override some settings
+  bool debug_optics = debug_optics_mode > 0;
+  if (debug_optics) {
+    printout(WARNING, "PFRICH_geo", "DEBUGGING PFRICH 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, "PFRICH_geo", "UNKNOWN debug_optics_mode");
+      return det;
+    };
+    aerogelVis = sensorVis;
+    gasvolVis = vesselVis = desc.invisible();
+  };
+
+  // BUILD VESSEL //////////////////////////////////////
+  /* - `vessel`: aluminum enclosure, the mother volume of the pfRICH
+   * - `gasvol`: gas volume, which fills `vessel`; all other volumes defined below
+   *   are children of `gasvol`
+   */
+
+  // tank solids
+  double boreDelta = vesselRmin1 - vesselRmin0;
+  Cone vesselTank(vesselLength / 2.0, vesselRmin1, vesselRmax1, vesselRmin0, vesselRmax0);
+  Cone gasvolTank(vesselLength / 2.0 - windowThickness, vesselRmin1 + wallThickness, vesselRmax1 - wallThickness,
+                  vesselRmin0 + wallThickness, vesselRmax0 - wallThickness);
+
+  //  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 = vesselTank;
+    gasvolSolid = gasvolTank;
+    break; // `!debug_optics`
+  case 1:
+    vesselSolid = vesselBox;
+    gasvolSolid = gasvolBox;
+    break;
+  case 2:
+    vesselSolid = vesselBox;
+    gasvolSolid = gasvolTank;
+    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., vesselLength / 2.0);
+  auto originBack  = Position(0., 0., -vesselLength / 2.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 + boreDelta * aerogelThickness / vesselLength, /* at backplane */
+      radiatorRmax,
+      radiatorRmin, /* at frontplane */
+      radiatorRmax
+      );
+  Cone filterSolid(
+      filterThickness/2,
+      radiatorRmin + boreDelta * (aerogelThickness + airGap + filterThickness) / vesselLength, /* at backplane */
+      radiatorRmax,
+      radiatorRmin + boreDelta * (aerogelThickness + airGap) / vesselLength, /* at frontplane */
+      radiatorRmax
+      );
+  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 - 0.5 * aerogelThickness) + 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, currently not in use)
+  );
+  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 airgap (FIXME: actually a gas 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();
+  };
+
+
+  // BUILD SENSORS ///////////////////////
+
+  // solid and volume: single sensor module
+  Box    sensorSolid(sensorSide / 2., sensorSide / 2., sensorThickness / 2.);
+  Volume sensorVol(detName + "_sensor", sensorSolid, sensorMat);
+  sensorVol.setVisAttributes(sensorVis);
+
+  // sensitivity
+  if (!debug_optics)
+    sensorVol.setSensitiveDetector(sens);
+
+  // sensor plane positioning: we want `sensorPlaneDist` to be the distance between the
+  // aerogel backplane (i.e., aerogel/filter boundary) and the sensor active surface (e.g, photocathode)
+  double sensorZpos     = radiatorFrontplane - aerogelThickness - sensorPlaneDist - 0.5 * sensorThickness;
+  auto   sensorPlanePos = Position(0., 0., sensorZpos) + originFront; // reference position
+  // miscellaneous
+  int    imod    = 0;           // module number
+  double tBoxMax = vesselRmax1; // sensors will be tiled in tBox, within annular limits
+
+  // SENSOR MODULE LOOP ------------------------
+  /* cartesian tiling loop
+   * - start at (x=0,y=0), to center the grid
+   * - loop over positive-x positions; for each, place the corresponding negative-x sensor too
+   * - nested similar loop over y positions
+   */
+  double sx, sy;
+  for (double usx = 0; usx <= tBoxMax; usx += sensorSide + sensorGap) {
+    for (int sgnx = 1; sgnx >= (usx > 0 ? -1 : 1); sgnx -= 2) {
+      for (double usy = 0; usy <= tBoxMax; usy += sensorSide + sensorGap) {
+        for (int sgny = 1; sgny >= (usy > 0 ? -1 : 1); sgny -= 2) {
+
+          // sensor (x,y) center
+          sx = sgnx * usx;
+          sy = sgny * usy;
+
+          // annular cut
+          if (std::hypot(sx, sy) < sensorPlaneRmin || std::hypot(sx, sy) > sensorPlaneRmax)
+            continue;
+
+          // placement (note: transformations are in reverse order)
+          auto sensorPV = gasvolVol.placeVolume(
+              sensorVol, Transform3D(Translation3D(sensorPlanePos.x(), sensorPlanePos.y(),
+                                                   sensorPlanePos.z()) // move to reference position
+                                     * Translation3D(sx, sy, 0.)       // move to grid position
+                                     ));
+
+          // generate LUT for module number -> sensor position, for readout mapping tests
+          // printf("%d %f %f\n",imod,sensorPV.position().x(),sensorPV.position().y());
+
+          // properties
+          sensorPV.addPhysVolID("module", imod);
+          DetElement sensorDE(det, Form("sensor_de_%d", imod), imod);
+          sensorDE.setPlacement(sensorPV);
+          if (!debug_optics) {
+            SkinSurface sensorSkin(desc, sensorDE, "sensor_optical_surface", sensorSurf,
+                                   sensorVol); // TODO: 3rd arg needs `imod`?
+            sensorSkin.isValid();
+          };
+
+          // increment sensor module number
+          imod++;
+        };
+      };
+    };
+  };
+  // END SENSOR MODULE LOOP ------------------------
+  //
+  // Add service material if desired
+  if (detElem.child("sensors").hasChild("services")) {
+    xml_comp_t x_service = detElem.child("sensors").child(_Unicode(services));
+    Assembly   service_vol("services");
+    service_vol.setVisAttributes(desc, x_service.visStr());
+
+    // Compute service total thickness from components
+    double     total_thickness = 0;
+    xml_coll_t ci(x_service, _Unicode(component));
+    for (ci.reset(), total_thickness = 0.0; ci; ++ci) {
+      total_thickness += xml_comp_t(ci).thickness();
+    }
+
+    int    ncomponents   = 0;
+    double thickness_sum = -total_thickness / 2.0;
+    for (xml_coll_t ci(x_service, _Unicode(component)); ci; ++ci, ncomponents++) {
+      xml_comp_t x_comp    = ci;
+      double     thickness = x_comp.thickness();
+      Tube       c_tube{sensorPlaneRmin, sensorPlaneRmax, thickness/2};
+      Volume     c_vol{_toString(ncomponents, "component%d"), c_tube, desc.material(x_comp.materialStr())};
+      c_vol.setVisAttributes(desc, x_comp.visStr());
+      service_vol.placeVolume(c_vol, Position(0, 0, thickness_sum + thickness / 2.0));
+      thickness_sum += thickness;
+    }
+    gasvolVol.placeVolume(service_vol,
+                          Transform3D(Translation3D(sensorPlanePos.x(), sensorPlanePos.y(),
+                                                    sensorPlanePos.z() - sensorThickness - total_thickness)));
+  }
+
+  // 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_PFRICH, createDetector)

src/PolyhedraEndcapCalorimeter2_geo.cpp

@@ -11,7 +11,7 @@
 //
 //==========================================================================
 //
-// Modified for TOPSiDE detector
+// Modified for ATHENA detector
 //
 //==========================================================================
 #include "DD4hep/DetFactoryHelper.h"
@@ -21,102 +21,112 @@ using namespace std;
 using namespace dd4hep;
 using namespace dd4hep::detail;
 
-static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector sens) {
-    xml_det_t x_det = e;
-    xml_dim_t dim = x_det.dimensions();
-    int det_id = x_det.id();
-    bool reflect = x_det.reflect(true);
-    string det_name = x_det.nameStr();
-    Material air = description.air();
-    int numsides = dim.numsides();
-    double rmin = dim.rmin();
-    double rmax = dim.rmax() * std::cos(M_PI / numsides);
-    double zmin = dim.zmin();
-    Layering layering(x_det);
-    double totalThickness = layering.totalThickness();
-    Volume endcapVol("endcap", PolyhedraRegular(numsides, rmin, rmax, totalThickness), air);
-    DetElement endcap("endcap", det_id);
+static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector sens)
+{
+  xml_det_t      x_det    = e;
+  xml_dim_t      dim      = x_det.dimensions();
+  int            det_id   = x_det.id();
+  bool           reflect  = x_det.reflect(true);
+  string         det_name = x_det.nameStr();
+  Material       air      = description.air();
+  int            numsides = dim.numsides();
+  xml::Component pos      = x_det.position();
+  double         rmin     = dim.rmin();
+  double         rmax     = dim.rmax();
+  double         zmin     = dim.zmin();
+  Layering       layering(x_det);
+  double         totalThickness = layering.totalThickness();
+  Volume         endcapVol("endcap", PolyhedraRegular(numsides, rmin, rmax, totalThickness), air);
+  DetElement     endcap("endcap", det_id);
 
-    int l_num = 1;
-    int layerType = 0;
-    double layerZ = -totalThickness / 2;
+  //std::cout << "totalThickness = " << totalThickness << "\n";
+  //std::cout << "zmin = " << zmin << "\n";
+  //std::cout << "rmin = " << rmin << "\n";
+  //std::cout << "rmax = " << rmax << "\n";
+  //std::cout << "nlayers = " << std::size(layering.layers()) << "\n";
+  int    l_num     = 1;
+  int    layerType = 0;
+  double layerZ    = -totalThickness / 2;
 
-    endcapVol.setAttributes(description, x_det.regionStr(), x_det.limitsStr(), x_det.visStr());
+  endcapVol.setAttributes(description, x_det.regionStr(), x_det.limitsStr(), x_det.visStr());
 
-    for (xml_coll_t xc(x_det, _U(layer)); xc; ++xc) {
-        xml_comp_t x_layer = xc;
-        double l_thick = layering.layer(l_num - 1)->thickness();
-        string l_name = _toString(layerType, "layer%d");
-        int l_repeat = x_layer.repeat();
-        Volume l_vol(l_name, PolyhedraRegular(numsides, rmin, rmax, l_thick), air);
-        vector<PlacedVolume> sensitives;
+  for (xml_coll_t xc(x_det, _U(layer)); xc; ++xc) {
+    //std::cout << "l_num = " << l_num << "\n";
+    //std::cout << "xc = " << xc << "\n";
+    xml_comp_t x_layer = xc;
+    double     l_thick = layering.layer(l_num - 1)->thickness();
+    //std::cout << "xc = " << xc << "\n";
+    string               l_name   = _toString(layerType, "layer%d");
+    int                  l_repeat = x_layer.repeat();
+    Volume               l_vol(l_name, PolyhedraRegular(numsides, rmin, rmax, l_thick), air);
+    vector<PlacedVolume> sensitives;
 
-        int s_num = 1;
-        double sliceZ = -l_thick / 2;
-        for (xml_coll_t xs(x_layer, _U(slice)); xs; ++xs) {
-            xml_comp_t x_slice = xs;
-            string s_name = _toString(s_num, "slice%d");
-            double s_thick = x_slice.thickness();
-            Material s_mat = description.material(x_slice.materialStr());
-            Volume s_vol(s_name, PolyhedraRegular(numsides, rmin, rmax, s_thick), s_mat);
+    int    s_num  = 1;
+    double sliceZ = -l_thick / 2;
+    for (xml_coll_t xs(x_layer, _U(slice)); xs; ++xs) {
+      xml_comp_t x_slice = xs;
+      string     s_name  = _toString(s_num, "slice%d");
+      double     s_thick = x_slice.thickness();
+      Material   s_mat   = description.material(x_slice.materialStr());
+      Volume     s_vol(s_name, PolyhedraRegular(numsides, rmin, rmax, s_thick), s_mat);
 
-            s_vol.setVisAttributes(description.visAttributes(x_slice.visStr()));
-            sliceZ += s_thick / 2;
-            PlacedVolume s_phv = l_vol.placeVolume(s_vol, Position(0, 0, sliceZ));
-            s_phv.addPhysVolID("slice", s_num);
-            if (x_slice.isSensitive()) {
-                sens.setType("calorimeter");
-                s_vol.setSensitiveDetector(sens);
-                sensitives.push_back(s_phv);
-            }
-            sliceZ += s_thick / 2;
-            s_num++;
-        }
-        l_vol.setVisAttributes(description.visAttributes(x_layer.visStr()));
-        if (l_repeat <= 0) throw std::runtime_error(x_det.nameStr() + "> Invalid repeat value");
-        for (int j = 0; j < l_repeat; ++j) {
-            string phys_lay = _toString(l_num, "layer%d");
-            layerZ += l_thick / 2;
-            DetElement layer_elt(endcap, phys_lay, l_num);
-            PlacedVolume pv = endcapVol.placeVolume(l_vol, Position(0, 0, layerZ));
-            pv.addPhysVolID("layer", l_num);
-            layer_elt.setPlacement(pv);
-            for (size_t ic = 0; ic < sensitives.size(); ++ic) {
-                PlacedVolume sens_pv = sensitives[ic];
-                DetElement comp_elt(layer_elt, sens_pv.volume().name(), l_num);
-                comp_elt.setPlacement(sens_pv);
-            }
-            layerZ += l_thick / 2;
-            ++l_num;
-        }
-        ++layerType;
+      s_vol.setVisAttributes(description.visAttributes(x_slice.visStr()));
+      sliceZ += s_thick / 2;
+      PlacedVolume s_phv = l_vol.placeVolume(s_vol, Position(0, 0, sliceZ));
+      s_phv.addPhysVolID("slice", s_num);
+      if (x_slice.isSensitive()) {
+        sens.setType("calorimeter");
+        s_vol.setSensitiveDetector(sens);
+        sensitives.push_back(s_phv);
+      }
+      sliceZ += s_thick / 2;
+      s_num++;
     }
-
-    double z_pos = zmin + totalThickness / 2;
-    PlacedVolume pv;
-    // Reflect it.
-    Assembly assembly(det_name);
-    DetElement endcapAssyDE(det_name, det_id);
-    Volume motherVol = description.pickMotherVolume(endcapAssyDE);
-    if (reflect) {
-        pv = assembly.placeVolume(endcapVol,
-                                  Transform3D(RotationZYX(M_PI / numsides, M_PI, 0), Position(0, 0, -z_pos)));
-        pv.addPhysVolID("barrel", 2);
-        Ref_t(endcap)->SetName((det_name + "_backward").c_str());
-        endcap.setPlacement(pv);
-    } else {
-        pv = assembly.placeVolume(endcapVol,
-                                  Transform3D(RotationZYX(M_PI / numsides, 0, 0), Position(0, 0, z_pos)));
-        pv.addPhysVolID("barrel", 1);
-        Ref_t(endcap)->SetName((det_name + "_forward").c_str());
-        endcap.setPlacement(pv);
+    l_vol.setVisAttributes(description.visAttributes(x_layer.visStr()));
+    if (l_repeat <= 0)
+      throw std::runtime_error(x_det.nameStr() + "> Invalid repeat value");
+    for (int j = 0; j < l_repeat; ++j) {
+      string phys_lay = _toString(l_num, "layer%d");
+      layerZ += l_thick / 2;
+      DetElement   layer_elt(endcap, phys_lay, l_num);
+      PlacedVolume pv = endcapVol.placeVolume(l_vol, Position(0, 0, layerZ));
+      pv.addPhysVolID("layer", l_num);
+      layer_elt.setPlacement(pv);
+      for (size_t ic = 0; ic < sensitives.size(); ++ic) {
+        PlacedVolume sens_pv = sensitives[ic];
+        DetElement   comp_elt(layer_elt, sens_pv.volume().name(), l_num);
+        comp_elt.setPlacement(sens_pv);
+      }
+      layerZ += l_thick / 2;
+      ++l_num;
     }
-    endcapAssyDE.add(endcap);
-    pv = motherVol.placeVolume(assembly);
-    pv.addPhysVolID("system", det_id);
-    endcapAssyDE.setPlacement(pv);
-    return endcapAssyDE;
+    ++layerType;
+  }
+
+  double       z_pos = zmin + totalThickness / 2;
+  PlacedVolume pv;
+  // Reflect it.
+  Assembly   assembly(det_name);
+  DetElement endcapAssyDE(det_name, det_id);
+  Volume     motherVol = description.pickMotherVolume(endcapAssyDE);
+  if (reflect) {
+    pv = assembly.placeVolume(endcapVol, Transform3D(RotationZYX(M_PI / numsides, M_PI, 0), Position(0, 0, -z_pos)));
+    pv.addPhysVolID("barrel", 2);
+    Ref_t(endcap)->SetName((det_name + "_backward").c_str());
+    endcap.setPlacement(pv);
+  } else {
+    pv = assembly.placeVolume(endcapVol, Transform3D(RotationZYX(M_PI / numsides, 0, 0), Position(0, 0, z_pos)));
+    pv.addPhysVolID("barrel", 1);
+    Ref_t(endcap)->SetName((det_name + "_forward").c_str());
+    endcap.setPlacement(pv);
+  }
+  endcapAssyDE.add(endcap);
+  pv = motherVol.placeVolume(assembly,Position(pos.x(),pos.y(),pos.z()));
+  pv.addPhysVolID("system", det_id);
+  endcapAssyDE.setPlacement(pv);
+  return endcapAssyDE;
 }
 
 // clang-format off
-DECLARE_DETELEMENT(refdet_PolyhedraEndcapCalorimeter2, create_detector)
+DECLARE_DETELEMENT(athena_PolyhedraEndcapCalorimeter2, create_detector)
+DECLARE_DETELEMENT(athena_PolyhedraEndcapCalorimeter, create_detector)

src/ScFiCalorimeter_geo.cpp

+//==========================================================================
+//  Scintillating fiber calorimeter with tower shape blocks
+//  reference: https://github.com/adamjaro/lmon/blob/master/calo/src/WScFiZXv3.cxx
+//  Support disk placement
+//--------------------------------------------------------------------------
+//  Author: Chao Peng (ANL)
+//  Date: 07/19/2021
+//==========================================================================
+
+#include "GeometryHelpers.h"
+#include "DD4hep/DetFactoryHelper.h"
+#include <XML/Helper.h>
+#include <iostream>
+#include <algorithm>
+#include <tuple>
+#include <math.h>
+
+using namespace dd4hep;
+using Point =  ROOT::Math::XYPoint;
+
+std::tuple<Volume, Position> build_module(const Detector &desc, const xml::Component &mod_x, SensitiveDetector &sens);
+
+// helper function to get x, y, z if defined in a xml component
+template<class XmlComp>
+Position get_xml_xyz(const XmlComp &comp, dd4hep::xml::Strng_t name)
+{
+    Position pos(0., 0., 0.);
+    if (comp.hasChild(name)) {
+        auto child = comp.child(name);
+        pos.SetX(dd4hep::getAttrOrDefault<double>(child, _Unicode(x), 0.));
+        pos.SetY(dd4hep::getAttrOrDefault<double>(child, _Unicode(y), 0.));
+        pos.SetZ(dd4hep::getAttrOrDefault<double>(child, _Unicode(z), 0.));
+    }
+    return pos;
+}
+
+// main
+static Ref_t create_detector(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);
+    sens.setType("calorimeter");
+    auto dim = detElem.dimensions();
+    auto rmin = dim.rmin();
+    auto rmax = dim.rmax();
+    auto length = dim.length();
+    auto phimin = dd4hep::getAttrOrDefault<double>(dim, _Unicode(phimin), 0.);
+    auto phimax = dd4hep::getAttrOrDefault<double>(dim, _Unicode(phimax), 2.*M_PI);
+    // envelope
+    Tube envShape(rmin, rmax, length/2., phimin, phimax);
+    Volume env(detName + "_envelope", envShape, desc.material("Air"));
+    env.setVisAttributes(desc.visAttributes(detElem.visStr()));
+
+    // build module
+    auto [modVol, modSize] = build_module(desc, detElem.child(_Unicode(module)), sens);
+    double modSizeR = std::sqrt(modSize.x() * modSize.x() + modSize.y() * modSize.y());
+    double assembly_rwidth = modSizeR*2.;
+    int nas = int((rmax - rmin) / assembly_rwidth) + 1;
+    std::vector<Assembly> assemblies;
+    // calorimeter block z-offsets (as blocks are shorter than the volume length)
+    const double block_offset = -0.5*(length - modSize.z());
+    for (int i = 0; i < nas; ++i) {
+        Assembly assembly(detName + Form("_ring%d", i + 1));
+        auto assemblyPV = env.placeVolume(assembly, Position{0., 0., block_offset});
+        assemblyPV.addPhysVolID("ring", i + 1);
+        assemblies.emplace_back(std::move(assembly));
+    }
+    // std::cout << assemblies.size() << std::endl;
+
+    int modid = 1;
+    for (int ix = 0; ix < int(2.*rmax / modSize.x()) + 1; ++ix) {
+        double mx = modSize.x() * ix - rmax;
+        for (int iy = 0; iy < int(2.*rmax / modSize.y()) + 1; ++iy) {
+            double my = modSize.y() * iy - rmax;
+            double mr = std::sqrt(mx*mx + my*my);
+            if (mr - modSizeR >= rmin && mr + modSizeR <= rmax) {
+                int ias = int((mr - rmin) / assembly_rwidth);
+                auto &assembly = assemblies[ias];
+                auto modPV = assembly.placeVolume(modVol, Position(mx, my, 0.));
+                modPV.addPhysVolID("module", modid++);
+            }
+        }
+    }
+
+    desc.add(Constant(detName + "_NModules", std::to_string(modid - 1)));
+
+    for (auto &assembly : assemblies) {
+        assembly.ptr()->Voxelize("");
+    }
+
+    // detector position and rotation
+    auto pos = get_xml_xyz(detElem, _Unicode(position));
+    auto rot = get_xml_xyz(detElem, _Unicode(rotation));
+    Volume motherVol = desc.pickMotherVolume(det);
+    Transform3D tr = Translation3D(pos.x(), pos.y(), pos.z()) * RotationZYX(rot.z(), rot.y(), rot.x());
+    PlacedVolume envPV = motherVol.placeVolume(env, tr);
+    envPV.addPhysVolID("system", detID);
+    det.setPlacement(envPV);
+    return det;
+}
+
+// helper function to build module with scintillating fibers
+std::tuple<Volume, Position> build_module(const Detector &desc, const xml::Component &mod_x, SensitiveDetector &sens)
+{
+    auto sx = mod_x.attr<double>(_Unicode(sizex));
+    auto sy = mod_x.attr<double>(_Unicode(sizey));
+    auto sz = mod_x.attr<double>(_Unicode(sizez));
+
+    Box modShape(sx/2., sy/2., sz/2.);
+    auto modMat = desc.material(mod_x.attr<std::string>(_Unicode(material)));
+    Volume modVol("module_vol", modShape, modMat);
+    if (mod_x.hasAttr(_Unicode(vis))) {
+        modVol.setVisAttributes(desc.visAttributes(mod_x.attr<std::string>(_Unicode(vis))));
+    }
+
+    if (mod_x.hasChild("fiber")) {
+      auto fiber_x  = mod_x.child(_Unicode(fiber));
+      auto fr       = fiber_x.attr<double>(_Unicode(radius));
+      auto fsx      = fiber_x.attr<double>(_Unicode(spacex));
+      auto fsy      = fiber_x.attr<double>(_Unicode(spacey));
+      auto foff     = dd4hep::getAttrOrDefault<double>(fiber_x, _Unicode(offset), 0.5*mm);
+      auto fiberMat = desc.material(fiber_x.attr<std::string>(_Unicode(material)));
+      Tube fiberShape(0., fr, sz/2.);
+      Volume fiberVol("fiber_vol", fiberShape, fiberMat);
+      fiberVol.setSensitiveDetector(sens);
+
+      // Fibers are placed in a honeycomb with the radius = sqrt(3)/2. * hexagon side length
+      // So each fiber is fully contained in a regular hexagon, which are placed as
+      //           ______________________________________
+      //           |          ____        ____          |
+      // even:     |         /    \      /    \         |
+      //           |    ____/      \____/      \____    |
+      //           |   /    \      /    \      /    \   |
+      // odd:      |  /      \____/      \____/      \  |
+      //           |  \      /    \      /    \      /  |
+      //           |   \____/      \____/      \____/   |
+      // even:     |        \      /    \      /        |
+      //           |         \____/      \____/      ___|___
+      //           |____________________________________|___offset
+      //                                              | |
+      //                                              |offset
+      // the parameters space x and space y are used to add additional spaces between the hexagons
+      double fside  = 2. / std::sqrt(3.) * fr;
+      double fdistx = 2. * fside + fsx;
+      double fdisty = 2. * fr + fsy;
+
+      // maximum numbers of the fibers, help narrow the loop range
+      int nx = int(sx / (2.*fr)) + 1;
+      int ny = int(sy / (2.*fr)) + 1;
+
+      // std::cout << sx << ", " << sy << ", " << fr << ", " << nx << ", " << ny << std::endl;
+
+      // place the fibers
+      double y0 = (foff + fside);
+      int nfibers = 0;
+      for (int iy = 0; iy < ny; ++iy) {
+          double y = y0 + fdisty * iy;
+          // about to touch the boundary
+          if ((sy - y) < y0) { break; }
+          double x0 = (iy % 2) ? (foff + fside) : (foff + fside + fdistx / 2.);
+          for (int ix = 0; ix < nx; ++ix) {
+              double x = x0 + fdistx * ix;
+              // about to touch the boundary
+              if ((sx - x) < x0) { break; }
+              auto fiberPV = modVol.placeVolume(fiberVol, nfibers++, Position{x - sx/2., y - sy/2., 0});
+              //std::cout << "(" << ix << ", " << iy << ", " << x - sx/2. << ", " << y - sy/2. << ", " << fr << "),\n";
+              fiberPV.addPhysVolID("fiber_x", ix + 1).addPhysVolID("fiber_y", iy + 1);
+          }
+      }
+    // if no fibers we make the module itself sensitive
+    } else {
+      modVol.setSensitiveDetector(sens);
+    }
+
+
+    return std::make_tuple(modVol, Position{sx, sy, sz});
+}
+
+DECLARE_DETELEMENT(ScFiCalorimeter, create_detector)
+

src/ShashlikCalorimeter_geo.cpp

+//==========================================================================
+//  Implementation for shashlik calorimeter modules
+//  it supports disk placements with (rmin, rmax), and (phimin, phimax)
+//--------------------------------------------------------------------------
+//  Author: Chao Peng (ANL)
+//  Date: 06/22/2021
+//==========================================================================
+
+#include "GeometryHelpers.h"
+#include "DD4hep/DetFactoryHelper.h"
+#include <XML/Layering.h>
+#include <XML/Helper.h>
+#include <iostream>
+#include <algorithm>
+#include <tuple>
+#include <math.h>
+
+using namespace dd4hep;
+
+static void add_disk_shashlik(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int id);
+
+// helper function to get x, y, z if defined in a xml component
+template<class XmlComp>
+Position get_xml_xyz(XmlComp &comp, dd4hep::xml::Strng_t name)
+{
+    Position pos(0., 0., 0.);
+    if (comp.hasChild(name)) {
+        auto child = comp.child(name);
+        pos.SetX(dd4hep::getAttrOrDefault<double>(child, _Unicode(x), 0.));
+        pos.SetY(dd4hep::getAttrOrDefault<double>(child, _Unicode(y), 0.));
+        pos.SetZ(dd4hep::getAttrOrDefault<double>(child, _Unicode(z), 0.));
+    }
+    return pos;
+}
+
+static Ref_t create_detector(Detector& desc, xml::Handle_t handle, SensitiveDetector sens)
+{
+    static const std::string func = "ShashlikCalorimeter";
+    xml::DetElement detElem = handle;
+    std::string detName = detElem.nameStr();
+    int detID = detElem.id();
+    DetElement det(detName, detID);
+    sens.setType("calorimeter");
+    // envelope
+    Assembly assembly(detName);
+
+    // module placement
+    xml::Component plm = detElem.child(_Unicode(placements));
+    int sector = 1;
+    for (xml::Collection_t mod(plm, _Unicode(disk)); mod; ++mod) {
+        add_disk_shashlik(desc, assembly, mod, sens, sector++);
+    }
+
+    // detector position and rotation
+    auto pos = get_xml_xyz(detElem, _Unicode(position));
+    auto rot = get_xml_xyz(detElem, _Unicode(rotation));
+    Volume motherVol = desc.pickMotherVolume(det);
+    Transform3D tr = Translation3D(pos.x(), pos.y(), pos.z()) * RotationZYX(rot.z(), rot.y(), rot.x());
+    PlacedVolume envPV = motherVol.placeVolume(assembly, tr);
+    envPV.addPhysVolID("system", detID);
+    det.setPlacement(envPV);
+    return det;
+}
+
+// helper function to build module with or w/o wrapper
+std::tuple<Volume, int, double, double> build_shashlik(Detector &desc, xml::Collection_t &plm, SensitiveDetector &sens)
+{
+    auto mod = plm.child(_Unicode(module));
+    // a modular volume
+    std::string shape = dd4hep::getAttrOrDefault(mod, _Unicode(shape), "square");
+    std::transform(shape.begin(), shape.end(), shape.begin(), [](char c) { return std::tolower(c); });
+    int nsides = 4;
+    if (shape == "hexagon") {
+        nsides = 6;
+    } else if (shape != "square") {
+        std::cerr << "ShashlikCalorimeter Error: Unsupported shape of module " << shape
+                  << ". Please choose from (square, hexagon). Proceed with square shape." << std::endl;
+    }
+    double slen = mod.attr<double>(_Unicode(side_length));
+    double rmax = slen/2./std::sin(M_PI/nsides);
+    Layering layering(mod);
+    auto len = layering.totalThickness();
+
+    // wrapper info
+    PolyhedraRegular mpoly(nsides, 0., rmax, len);
+    Volume mvol("shashlik_module_vol", mpoly, desc.air());
+    mvol.setVisAttributes(desc.visAttributes(dd4hep::getAttrOrDefault(mod, _Unicode(vis), "GreenVis")));
+
+    double gap = 0.;
+    Volume wvol("shashlik_wrapper_vol");
+    if (plm.hasChild(_Unicode(wrapper))) {
+        auto wrap = plm.child(_Unicode(wrapper));
+        gap = wrap.attr<double>(_Unicode(thickness));
+        if (gap > 1e-6*mm) {
+            wvol.setSolid(PolyhedraRegular(nsides, 0., rmax + gap, len));
+            wvol.setMaterial(desc.material(wrap.attr<std::string>(_Unicode(material))));
+            wvol.setVisAttributes(desc.visAttributes(dd4hep::getAttrOrDefault(wrap, _Unicode(vis), "WhiteVis")));
+            wvol.placeVolume(mvol, Position{0., 0., 0.});
+        }
+    }
+
+    // layer start point
+    double lz  = -len/2.;
+    int lnum = 1;
+    // Loop over the sets of layer elements in the detector.
+    for (xml_coll_t li(mod, _U(layer)); li; ++li) {
+        int repeat = li.attr<int>(_Unicode(repeat));
+        // Loop over number of repeats for this layer.
+        for (int j = 0; j < repeat; j++) {
+            std::string lname = Form("layer%d", lnum);
+            double lthick = layering.layer(lnum - 1)->thickness();  // Layer's thickness.
+            PolyhedraRegular lpoly(nsides, 0., rmax, lthick);
+            Volume lvol(lname, lpoly, desc.air());
+
+            // Loop over the sublayers or slices for this layer.
+            int snum = 1;
+            double sz = -lthick/2.;
+            for (xml_coll_t si(li, _U(slice)); si; ++si) {
+                std::string sname = Form("slice%d", snum);
+                double sthick = si.attr<double>(_Unicode(thickness));
+                PolyhedraRegular spoly(nsides, 0., rmax, sthick);
+                Volume svol(sname, spoly, desc.material(si.attr<std::string>(_Unicode(material))));
+
+                std::string issens = dd4hep::getAttrOrDefault(si, _Unicode(sensitive), "no");
+                std::transform(issens.begin(), issens.end(), issens.begin(), [](char c) { return std::tolower(c); });
+                if ((issens == "yes") || (issens == "y") || (issens == "true")) {
+                    svol.setSensitiveDetector(sens);
+                }
+                svol.setAttributes(desc,
+                        dd4hep::getAttrOrDefault(si, _Unicode(region), ""),
+                        dd4hep::getAttrOrDefault(si, _Unicode(limits), ""),
+                        dd4hep::getAttrOrDefault(si, _Unicode(vis), "InvisibleNoDaughters"));
+
+                // Slice placement.
+                auto slicePV = lvol.placeVolume(svol, Position(0, 0, sz + sthick/2.));
+                slicePV.addPhysVolID("slice", snum++);
+                // Increment Z position of slice.
+                sz += sthick;
+            }
+
+            // Set region, limitset, and vis of layer.
+            lvol.setAttributes(desc,
+                    dd4hep::getAttrOrDefault(li, _Unicode(region), ""),
+                    dd4hep::getAttrOrDefault(li, _Unicode(limits), ""),
+                    dd4hep::getAttrOrDefault(li, _Unicode(vis), "InvisibleNoDaughters"));
+
+            auto layerPV = mvol.placeVolume(lvol, Position(0, 0, lz + lthick/2));
+            layerPV.addPhysVolID("layer", lnum++);
+            // Increment to next layer Z position.
+            lz += lthick;
+        }
+    }
+
+    if (gap > 1e-6*mm) {
+        return std::make_tuple(wvol, nsides, 2.*std::sin(M_PI/nsides)*(rmax + gap), len);
+    } else {
+        return std::make_tuple(mvol, nsides, slen, len);
+    }
+}
+
+// place disk of modules
+static void add_disk_shashlik(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
+{
+    auto [mvol, nsides, sidelen, len]  = build_shashlik(desc, plm, sens);
+    int sector_id = dd4hep::getAttrOrDefault<int>(plm, _Unicode(sector), sid);
+    int id_begin = dd4hep::getAttrOrDefault<int>(plm, _Unicode(id_begin), 1);
+    double rmin = plm.attr<double>(_Unicode(rmin));
+    double rmax = plm.attr<double>(_Unicode(rmax));
+    double phimin = dd4hep::getAttrOrDefault<double>(plm, _Unicode(phimin), 0.);
+    double phimax = dd4hep::getAttrOrDefault<double>(plm, _Unicode(phimax), 2.*M_PI);
+
+    auto points = (nsides == 6) ? athena::geo::fillHexagons({0., 0.}, sidelen, rmin, rmax, phimin, phimax)
+                                : athena::geo::fillSquares({0., 0.}, sidelen*1.414, rmin, rmax, phimin, phimax);
+    // placement to mother
+    auto pos = get_xml_xyz(plm, _Unicode(position));
+    auto rot = get_xml_xyz(plm, _Unicode(rotation));
+    int mid = 0;
+    for (auto &p : points) {
+        Transform3D tr = RotationZYX(rot.z(), rot.y(), rot.x())
+                       * Translation3D(pos.x() + p.x(), pos.y() + p.y(), pos.z() + len/2.)
+                       * RotationZ((nsides == 4) ? 45*degree : 0);
+        auto modPV = env.placeVolume(mvol, tr);
+        modPV.addPhysVolID("sector", sector_id).addPhysVolID("module", id_begin + mid++);
+    }
+}
+
+
+DECLARE_DETELEMENT(ShashlikCalorimeter, create_detector)
+

src/SimpleDiskDetector_geo.cpp

+//==========================================================================
+//  AIDA Detector description implementation
+//--------------------------------------------------------------------------
+// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
+// All rights reserved.
+//
+// For the licensing terms see $DD4hepINSTALL/LICENSE.
+// For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
+//
+// Author     : M.Frank
+//
+//==========================================================================
+//
+// Specialized generic detector constructor
+//
+//==========================================================================
+#include "DD4hep/DetFactoryHelper.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;
+
+static Ref_t SimpleDiskDetector_create_detector(Detector& description, xml_h e, SensitiveDetector sens)
+{
+  xml_det_t      x_det    = e;
+  Material       air      = description.air();
+  string         det_name = x_det.nameStr();
+  bool           reflect  = x_det.reflect();
+  DetElement     sdet(det_name, x_det.id());
+  Assembly       assembly(det_name);
+  PlacedVolume   pv;
+  int            l_num = 0;
+  xml::Component pos   = x_det.position();
+
+  Acts::ActsExtension* detWorldExt = new Acts::ActsExtension();
+  detWorldExt->addType("endcap", "detector");
+  sdet.addExtension<Acts::ActsExtension>(detWorldExt);
+
+  for (xml_coll_t i(x_det, _U(layer)); i; ++i, ++l_num) {
+    xml_comp_t x_layer    = i;
+    string     l_nam      = det_name + _toString(l_num, "_layer%d");
+    double     zmin       = x_layer.inner_z();
+    double     rmin       = x_layer.inner_r();
+    double     rmax       = x_layer.outer_r();
+    double     z          = zmin;
+    double     layerWidth = 0.;
+    int        s_num      = 0;
+
+    for (xml_coll_t j(x_layer, _U(slice)); j; ++j) {
+      double thickness = xml_comp_t(j).thickness();
+      layerWidth += thickness;
+    }
+    Tube   l_tub(rmin, rmax, layerWidth/2.0, 2 * M_PI);
+    Volume l_vol(l_nam, l_tub, air);
+    l_vol.setVisAttributes(description, x_layer.visStr());
+    DetElement layer;
+    PlacedVolume layer_pv;
+    if (!reflect) {
+      layer = DetElement(sdet, l_nam + "_pos", l_num);
+      layer_pv = assembly.placeVolume(l_vol, Position(0, 0, zmin + layerWidth / 2.));
+      layer_pv.addPhysVolID("barrel", 3).addPhysVolID("layer", l_num);
+      layer.setPlacement(layer_pv);
+      Acts::ActsExtension* layerExtension = new Acts::ActsExtension();
+      layerExtension->addType("sensitive disk", "layer");
+      //layerExtension->addType("axes", "definitions", "XzY");
+      // need all four of these or else it is ignored.
+      //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");
+
+      layer.addExtension<Acts::ActsExtension>(layerExtension);
+    } else {
+      layer = DetElement(sdet, l_nam + "_neg", l_num);
+      layer_pv = assembly.placeVolume(l_vol, Transform3D(RotationY(M_PI), Position(0, 0, -zmin - layerWidth / 2)));
+      layer_pv.addPhysVolID("barrel", 2).addPhysVolID("layer", l_num);
+      layer.setPlacement(layer_pv);
+      // DetElement layerR = layer.clone(l_nam+"_neg");
+      // sdet.add(layerR.setPlacement(pv));
+      Acts::ActsExtension* layerExtension = new Acts::ActsExtension();
+      layerExtension->addType("sensitive disk", "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");
+      layer.addExtension<Acts::ActsExtension>(layerExtension);
+    }
+
+    double tot_thickness = -layerWidth / 2.0;
+    for (xml_coll_t j(x_layer, _U(slice)); j; ++j, ++s_num) {
+      xml_comp_t x_slice = j;
+      double     thick   = x_slice.thickness();
+      Material   mat     = description.material(x_slice.materialStr());
+      string     s_nam   = l_nam + _toString(s_num, "_slice%d");
+      Volume     s_vol(s_nam, Tube(rmin, rmax, thick/2.0), mat);
+      if(!reflect){
+        s_nam += "_pos";
+      } else {
+        s_nam += "_neg";
+      }
+      DetElement slice_de(layer, s_nam , s_num);
+      if (x_slice.isSensitive()) {
+        sens.setType("tracker");
+        s_vol.setSensitiveDetector(sens);
+        Acts::ActsExtension* sensorExtension = new Acts::ActsExtension();
+        //sensorExtension->addType("sensor", "detector");
+        slice_de.addExtension<Acts::ActsExtension>(sensorExtension);
+      }
+      s_vol.setAttributes(description, x_slice.regionStr(), x_slice.limitsStr(), x_slice.visStr());
+      pv = l_vol.placeVolume(s_vol, Position(0, 0, tot_thickness + thick / 2));
+      pv.addPhysVolID("slice", s_num);
+      slice_de.setPlacement(pv);
+      tot_thickness = tot_thickness + thick;
+    }
+
+  }
+  if (x_det.hasAttr(_U(combineHits))) {
+    sdet.setCombineHits(x_det.attr<bool>(_U(combineHits)), sens);
+  }
+  pv = description.pickMotherVolume(sdet).placeVolume(assembly, Position(pos.x(), pos.y(), pos.z()));
+  pv.addPhysVolID("system", x_det.id()); // Set the subdetector system ID.
+  sdet.setPlacement(pv);
+  return sdet;
+}
+
+DECLARE_DETELEMENT(ref_SolenoidEndcap, SimpleDiskDetector_create_detector)
+DECLARE_DETELEMENT(athena_SolenoidEndcap, SimpleDiskDetector_create_detector)

src/SolenoidCoil_geo.cpp

+//==========================================================================
+//  AIDA Detector description implementation 
+//--------------------------------------------------------------------------
+// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
+// All rights reserved.
+//
+// For the licensing terms see $DD4hepINSTALL/LICENSE.
+// For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
+//
+// Author     : M.Frank
+//
+//==========================================================================
+//
+// Specialized generic detector constructor
+// 
+//==========================================================================
+#include "DD4hep/DetFactoryHelper.h"
+
+using namespace std;
+using namespace dd4hep;
+using namespace dd4hep::detail;
+
+static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector sens)  {
+  xml_det_t  x_det     = e;
+  string     det_name  = x_det.nameStr();
+  Material   air       = description.air();
+  DetElement sdet        (det_name,x_det.id());
+  Assembly   assembly    (det_name+"_assembly");
+  PlacedVolume pv;
+  int n = 0;
+  xml::Component  pos  = x_det.position();
+
+  for(xml_coll_t i(x_det,_U(layer)); i; ++i, ++n)  {
+    xml_comp_t x_layer = i;
+    string  l_name = det_name+_toString(n,"_layer%d");
+    double  outer_z = x_layer.outer_z();
+    double  rmin = x_layer.inner_r();
+    double  r    = rmin;
+    DetElement layer(sdet,_toString(n,"layer%d"),x_layer.id());
+    Tube    l_tub (rmin,2*rmin,outer_z);
+    Volume  l_vol(l_name,l_tub,air);
+    int im = 0;
+
+    for(xml_coll_t j(x_layer,_U(slice)); j; ++j, ++im)  {
+      // If slices are only given a thickness attribute, they are radially concentric slices
+      // If slices are given an inner_z attribute, they are longitudinal slices with equal rmin
+      xml_comp_t x_slice = j;
+      Material mat = description.material(x_slice.materialStr());
+      string s_name= l_name+_toString(im,"_slice%d");
+      double thickness = x_slice.thickness();
+      double s_outer_z = dd4hep::getAttrOrDefault(x_slice, _Unicode(outer_z), outer_z);
+      double s_inner_z = dd4hep::getAttrOrDefault(x_slice, _Unicode(inner_z), 0.0*dd4hep::cm);
+      Tube   s_tub(r,r+thickness,(s_inner_z > 0? 0.5*(s_outer_z-s_inner_z): s_outer_z),2*M_PI);
+      Volume s_vol(s_name, s_tub, mat);
+
+      if ( x_slice.isSensitive() ) {
+        sens.setType("tracker");
+        s_vol.setSensitiveDetector(sens);
+      }
+      // Set Attributes
+      s_vol.setAttributes(description,x_slice.regionStr(),x_slice.limitsStr(),x_slice.visStr());
+      if (s_inner_z > 0) {
+        // Place off-center volumes twice
+        Position s_pos(0, 0, 0.5*(s_outer_z+s_inner_z));
+        pv = l_vol.placeVolume(s_vol, -s_pos);
+        pv = l_vol.placeVolume(s_vol, +s_pos);
+      } else {
+        r += thickness;
+        pv = l_vol.placeVolume(s_vol);
+      }
+      // Slices have no extra id. Take the ID of the layer!
+      pv.addPhysVolID("slice",im);
+    }
+    l_tub.setDimensions(rmin,r,outer_z);
+    //cout << l_name << " " << rmin << " " << r << " " << z << endl;
+    l_vol.setVisAttributes(description,x_layer.visStr());
+      
+    pv = assembly.placeVolume(l_vol);
+    pv.addPhysVolID("layer",n);
+    layer.setPlacement(pv);
+  }
+  if ( x_det.hasAttr(_U(combineHits)) ) {
+    sdet.setCombineHits(x_det.combineHits(),sens);
+  }
+
+  pv = description.pickMotherVolume(sdet).placeVolume(assembly,Position(pos.x(),pos.y(),pos.z()));
+  pv.addPhysVolID("system",sdet.id()).addPhysVolID("barrel",0);
+  sdet.setPlacement(pv);
+  return sdet;
+}
+
+DECLARE_DETELEMENT(athena_SolenoidCoil,create_detector)

src/TestDetector.cpp

-#include <XML/Helper.h>
-
-using namespace dd4hep;
-
-static Ref_t createDetector(Detector& desc, xml_h handle, SensitiveDetector sens) {
-  xml::DetElement detElem = handle;
-  std::string     detName = detElem.nameStr();
-  int             detID   = detElem.id();
-
-  xml::Component dims      = detElem.dimensions();
-  double         rInner    = dims.inner_radius();
-  double         rMin      = dims.rmin();
-  double         thickness = dims.thickness();
-  double         innerZ    = dims.inner_z();
-  double         angle     = dims.angle();
-
-  Material mat = desc.material(detElem.materialStr());
-
-  Tube             outerTubeShape(rMin, rInner + thickness, innerZ + thickness);
-  Tube             innerTubeShape(0, rInner, innerZ);
-  SubtractionSolid unchamferedShape(outerTubeShape, innerTubeShape);
-  Cone             chamferShape(thickness, 0, rMin, 0, rMin + 2 * tan(angle) * thickness);
-  SubtractionSolid detShape(unchamferedShape, chamferShape, Position(0, 0, innerZ + thickness));
-  Volume           detVol(detName, detShape, mat);
-
-  detVol.setVisAttributes(desc.visAttributes(detElem.visStr()));
-
-  DetElement   det(detName, detID);
-  Volume       motherVol = desc.pickMotherVolume(det);
-  PlacedVolume detPV     = motherVol.placeVolume(detVol);
-  det.setPlacement(detPV);
-  return det;
-}
-
-// clang-format off
-DECLARE_DETELEMENT(TestDetector, createDetector)

src/TrapEndcapTracker_geo.cpp

+/** \addtogroup Trackers Trackers
+ * \brief Type: **BarrelTrackerWithFrame**.
+ * \author W. Armstrong
+ *
+ * \ingroup trackers
+ *
+ * @{
+ */
+#include <array>
+#include <map>
+#include "DD4hep/DetFactoryHelper.h"
+#include "DD4hep/Printout.h"
+#include "DD4hep/Shapes.h"
+#include "DDRec/Surface.h"
+#include "DDRec/DetectorData.h"
+#include "XML/Utilities.h"
+#include "XML/Layering.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::rec;
+using namespace dd4hep::detail;
+
+/** Endcap Trapezoidal Tracker.
+ *
+ * @author Whitney Armstrong
+ *
+ */
+static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector sens)
+{
+  typedef vector<PlacedVolume> Placements;
+  xml_det_t                    x_det    = e;
+  Material                     vacuum   = description.vacuum();
+  int                          det_id   = x_det.id();
+  string                       det_name = x_det.nameStr();
+  bool                         reflect  = x_det.reflect(false);
+  DetElement                   sdet(det_name, det_id);
+  Assembly                     assembly(det_name);
+
+  Material  air  = description.material("Air");
+  // Volume      assembly    (det_name,Box(10000,10000,10000),vacuum);
+  Volume                  motherVol = description.pickMotherVolume(sdet);
+  int                     m_id = 0, c_id = 0, n_sensor = 0;
+  map<string, Volume>     modules;
+  map<string, Placements> sensitives;
+  map<string, std::vector<VolPlane>>        volplane_surfaces;
+  map<string, std::array<double, 2>> module_thicknesses;
+  PlacedVolume            pv;
+
+
+  // ACTS extension
+  {
+    Acts::ActsExtension* detWorldExt = new Acts::ActsExtension();
+    detWorldExt->addType("endcap", "detector");
+    // SJJ probably need to set the envelope here, as ACTS can't figure
+    // that out for Assembly volumes. May also need binning to properly pick up
+    // on the support material @TODO
+    //
+    // 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);
+  }
+
+  assembly.setVisAttributes(description.invisible());
+  sens.setType("tracker");
+
+  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,(reflect?-1.0:1.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));
+  }
+
+  for (xml_coll_t mi(x_det, _U(module)); mi; ++mi, ++m_id) {
+    xml_comp_t x_mod = mi;
+    string     m_nam = x_mod.nameStr();
+    xml_comp_t trd   = x_mod.trd();
+
+    double     posY;
+    double     x1 = trd.x1();
+    double     x2 = trd.x2();
+    double     z  = trd.z();
+    double     total_thickness = 0.;
+    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();
+
+    double     thickness_so_far = 0.0;
+    double     thickness_sum    = -total_thickness / 2.0;
+    double     y1               = total_thickness / 2;
+    double     y2               = total_thickness / 2;
+    Trapezoid m_solid(x1, x2, y1, y2, z);
+    Volume m_volume(m_nam, m_solid, vacuum);
+    m_volume.setVisAttributes(description.visAttributes(x_mod.visStr()));
+
+    Solid  frame_s;
+    if(x_mod.hasChild("frame")){
+      // build frame from trd (assumed to be smaller)
+      xml_comp_t m_frame         = x_mod.child(_U(frame));
+      xml_comp_t f_pos           = m_frame.child(_U(position));
+      xml_comp_t frame_trd       = m_frame.trd();
+      double     frame_thickness = getAttrOrDefault(m_frame, _U(thickness), total_thickness);
+      double     frame_x1        = frame_trd.x1();
+      double     frame_x2        = frame_trd.x2();
+      double     frame_z         = frame_trd.z();
+      // make the frame match the total thickness if thickness attribute is not given
+      Trapezoid        f_solid1(x1, x2,frame_thickness / 2.0, frame_thickness / 2.0, z);
+      Trapezoid        f_solid(frame_x1, frame_x2, frame_thickness / 2.0, frame_thickness / 2.0, frame_z) ;
+      SubtractionSolid frame_shape(f_solid1, f_solid);
+      frame_s = frame_shape;
+
+      Material f_mat  = description.material(m_frame.materialStr());
+      Volume f_vol(m_nam + "_frame", frame_shape, f_mat);
+      f_vol.setVisAttributes(description.visAttributes(m_frame.visStr()));
+
+      // figure out how to best place
+      pv = m_volume.placeVolume(f_vol, Position(f_pos.x(), f_pos.y(),  f_pos.z()));
+    }
+
+    for (ci.reset(), n_sensor = 1, c_id = 0, posY = -y1; ci; ++ci, ++c_id) {
+      xml_comp_t c           = ci;
+      double     c_thick     = c.thickness();
+      auto       comp_x1     = getAttrOrDefault(c, _Unicode(x1), x1);
+      auto       comp_x2     = getAttrOrDefault(c, _Unicode(x2), x2);
+      auto       comp_height = getAttrOrDefault(c, _Unicode(height), z);
+
+      Material c_mat  = description.material(c.materialStr());
+      string   c_name = _toString(c_id, "component%d");
+
+      Trapezoid comp_s1(comp_x1, comp_x2, c_thick / 2e0, c_thick / 2e0, comp_height);
+      Solid  comp_shape = comp_s1;
+      if(frame_s.isValid()) {
+        comp_shape = SubtractionSolid( comp_s1, frame_s); 
+      }
+      Volume   c_vol(c_name, comp_shape, c_mat);
+
+      c_vol.setVisAttributes(description.visAttributes(c.visStr()));
+      pv = m_volume.placeVolume(c_vol, Position(0, posY + c_thick / 2, 0));
+      if (c.isSensitive()) {
+        module_thicknesses[m_nam] = {thickness_so_far + c_thick/2.0, total_thickness-thickness_so_far - c_thick/2.0};
+        //std::cout << " adding sensitive volume" << c_name << "\n";
+        sdet.check(n_sensor > 2, "SiTrackerEndcap2::fromCompact: " + c_name + " Max of 2 modules allowed!");
+        pv.addPhysVolID("sensor", n_sensor);
+        sens.setType("tracker");
+        c_vol.setSensitiveDetector(sens);
+        sensitives[m_nam].push_back(pv);
+        ++n_sensor;
+        // -------- create a measurement plane for the tracking surface attched to the sensitive volume -----
+        Vector3D u(0., 0., -1.);
+        Vector3D v(-1., 0., 0.);
+        Vector3D n(0., 1., 0.);
+        // 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);
+
+        //--------------------------------------------
+      }
+      posY += c_thick;
+      thickness_sum += c_thick;
+      thickness_so_far += c_thick;
+    }
+    modules[m_nam] = m_volume;
+  }
+
+  for (xml_coll_t li(x_det, _U(layer)); li; ++li) {
+    xml_comp_t x_layer(li);
+    int        l_id    = x_layer.id();
+    int        mod_num = 1;
+
+    xml_comp_t l_env      = x_layer.child(_U(envelope));
+    string     layer_name = det_name + std::string("_layer") + std::to_string(l_id);
+
+    std::string layer_vis    = l_env.attr<std::string>(_Unicode(vis));
+    double      layer_rmin   = l_env.attr<double>(_Unicode(rmin));
+    double      layer_rmax   = l_env.attr<double>(_Unicode(rmax));
+    double      layer_length = l_env.attr<double>(_Unicode(length));
+    double      layer_zstart = l_env.attr<double>(_Unicode(zstart));
+    double      layer_center_z =  layer_zstart + layer_length/2.0;
+    //printout(INFO,"ROOTGDMLParse","+++ Read geometry from GDML file file:%s",input.c_str());
+    //std::cout << "SiTracker Endcap layer " << l_id << " zstart = " << layer_zstart/dd4hep::mm << "mm ( " << layer_length/dd4hep::mm << " mm thick )\n";
+
+    //Assembly    layer_assembly(layer_name);
+    //assembly.placeVolume(layer_assembly);
+    Tube       layer_tub(layer_rmin, layer_rmax, layer_length / 2);
+    Volume     layer_vol(layer_name, layer_tub, air); // Create the layer envelope volume.
+    layer_vol.setVisAttributes(description.visAttributes(layer_vis));
+
+    PlacedVolume layer_pv;
+    if (reflect) {
+      layer_pv =
+          assembly.placeVolume(layer_vol, Transform3D(RotationZYX(0.0, -M_PI, 0.0), Position(0, 0, -layer_center_z)));
+      layer_pv.addPhysVolID("layer", l_id);
+      layer_name += "_N";
+    } else {
+      layer_pv = assembly.placeVolume(layer_vol, Position(0, 0, layer_center_z));
+      layer_pv.addPhysVolID("layer", l_id);
+      layer_name += "_P";
+    }
+    DetElement layer_element(sdet, layer_name, l_id);
+    layer_element.setPlacement(layer_pv);
+
+    Acts::ActsExtension* layerExtension = new Acts::ActsExtension();
+    layerExtension->addType("sensitive disk", "layer");
+    //layerExtension->addType("axes", "definitions", "XZY");
+    //layerExtension->addType("sensitive disk", "layer");
+    //layerExtension->addType("axes", "definitions", "XZY");
+    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");
+    }
+    layer_element.addExtension<Acts::ActsExtension>(layerExtension);
+
+    for (xml_coll_t ri(x_layer, _U(ring)); ri; ++ri) {
+      xml_comp_t  x_ring   = ri;
+      double      r        = x_ring.r();
+      double      phi0     = x_ring.phi0(0);
+      double      zstart   = x_ring.zstart();
+      double      dz       = x_ring.dz(0);
+      int         nmodules = x_ring.nmodules();
+      string      m_nam    = x_ring.moduleStr();
+      Volume      m_vol    = modules[m_nam];
+      double      iphi     = 2 * M_PI / nmodules;
+      double      phi      = phi0;
+      Placements& sensVols = sensitives[m_nam];
+
+      for (int k = 0; k < nmodules; ++k) {
+        string     m_base = _toString(l_id, "layer%d") + _toString(mod_num, "_module%d");
+        double     x      = -r * std::cos(phi);
+        double     y      = -r * std::sin(phi);
+
+        if (!reflect) {
+          DetElement module(layer_element, m_base + "_pos", det_id);
+          pv = layer_vol.placeVolume(
+              m_vol, Transform3D(RotationZYX(0, -M_PI / 2 - phi, -M_PI / 2), Position(x, y, zstart + dz)));
+          pv.addPhysVolID("module", mod_num);
+          module.setPlacement(pv);
+          for (size_t ic = 0; ic < sensVols.size(); ++ic) {
+            PlacedVolume sens_pv = sensVols[ic];
+            DetElement   comp_elt(module, sens_pv.volume().name(), mod_num);
+            comp_elt.setPlacement(sens_pv);
+        //std::cout << " adding ACTS extension" << "\n";
+            Acts::ActsExtension* moduleExtension = new Acts::ActsExtension("XZY");
+            comp_elt.addExtension<Acts::ActsExtension>(moduleExtension);
+            volSurfaceList(comp_elt)->push_back(volplane_surfaces[m_nam][ic]);
+          }
+        } else {
+          pv = layer_vol.placeVolume(
+              m_vol, Transform3D(RotationZYX(0, -M_PI / 2 - phi, -M_PI / 2), Position(x, y, -zstart - dz)));
+          pv.addPhysVolID("module", mod_num);
+          DetElement r_module(layer_element, m_base + "_neg", det_id);
+          r_module.setPlacement(pv);
+          for (size_t ic = 0; ic < sensVols.size(); ++ic) {
+            PlacedVolume sens_pv = sensVols[ic];
+            DetElement   comp_elt(r_module, sens_pv.volume().name(), mod_num);
+            comp_elt.setPlacement(sens_pv);
+        //std::cout << " adding ACTS extension" << "\n";
+            Acts::ActsExtension* moduleExtension = new Acts::ActsExtension("XZY");
+            comp_elt.addExtension<Acts::ActsExtension>(moduleExtension);
+            volSurfaceList(comp_elt)->push_back(volplane_surfaces[m_nam][ic]);
+          }
+        }
+        dz = -dz;
+        phi += iphi;
+        ++mod_num;
+      }
+    }
+  }
+  pv = motherVol.placeVolume(assembly,Position(0,0,(reflect?-1.0e-9:1.0e-9)) );
+  pv.addPhysVolID("system", det_id);
+  sdet.setPlacement(pv);
+  return sdet;
+}
+
+//@}
+// clang-format off
+DECLARE_DETELEMENT(athena_TrapEndcapTracker, create_detector)
+DECLARE_DETELEMENT(athena_GEMTrackerEndcap, create_detector)

src/cb_Solenoid.cpp

-#include <XML/Helper.h>
-///////////////////////////
-// Central Barrel Solenoid
-///////////////////////////
-
-using namespace dd4hep;
-
-static Ref_t createDetector(Detector& desc, xml_h handle, SensitiveDetector sens)
-{
-  xml::DetElement detElem = handle;
-  std::string     detName = detElem.nameStr();
-  int             detID   = detElem.id();
-
-  xml::Component dims = detElem.dimensions();
-  double SizeZ  = dims.z();     // Size in Z direction
-  double ROut   = dims.rmax();  // Outer diameter
-  double RIn    = dims.rmin();  // Inner diameter
-  double ShiftZ = dims.delta();
-  Material mat = desc.material(detElem.materialStr());
-
-  Tube   cb_Solenoid_GVol_Solid(RIn, ROut, SizeZ / 2., 0., 360 * deg);
-  Volume detVol("cb_Solenoid_GVol_Logic", cb_Solenoid_GVol_Solid, mat);
-
-  detVol.setVisAttributes(desc.visAttributes(detElem.visStr()));
-
-  DetElement   det(detName, detID);
-  Volume       motherVol = desc.pickMotherVolume(det);
-  Transform3D  tr(RotationZYX(0., 0., 0.), Position(0., 0., ShiftZ));
-  PlacedVolume detPV = motherVol.placeVolume(detVol, tr);
-  det.setPlacement(detPV);
-  return det;
-}
-
-// clang-format off
-DECLARE_DETELEMENT(cb_Solenoid, createDetector)

src/deprecated/BackwardGEMTracker_geo.cpp

+#include "DDRec/Surface.h"
+#include "DDRec/DetectorData.h"
+#include "DD4hep/OpticalSurfaces.h"
+#include "DD4hep/DetFactoryHelper.h"
+#include "DD4hep/Printout.h"
+#include <XML/Helper.h>
+//////////////////////////////////
+// Electron Endcap GEM Tracking
+//////////////////////////////////
+
+using namespace std;
+using namespace dd4hep;
+
+static Ref_t createDetector(Detector& desc, xml_h e, SensitiveDetector sens)
+{
+  xml_det_t  x_det      = e;
+  string     detName    = x_det.nameStr();
+  int        detID      = x_det.id();
+
+  xml_dim_t  dim        = x_det.dimensions();
+  double     RIn        = dim.rmin();
+  double     ROut       = dim.rmax();
+  double     SizeZ      = dim.length();
+  xml_dim_t  pos        = x_det.position();
+
+  Material   Vacuum     = desc.material("Vacuum");
+
+  // Create Global Volume 
+  Tube ce_GEM_GVol_Solid(RIn, ROut, SizeZ / 2.0, 0., 360.0 * deg);
+  Volume detVol("ce_GEM_GVol_Logic", ce_GEM_GVol_Solid, Vacuum);
+  detVol.setVisAttributes(desc.visAttributes(x_det.visStr()));
+
+  // Construct Layers
+  xml_comp_t x_layer = x_det.child(_U(layer));
+  const int repeat   = x_layer.repeat();
+  
+  xml_comp_t x_slice = x_layer.child(_U(slice));
+  Material slice_mat = desc.material(x_slice.materialStr());
+  double layerSizeZ = x_slice.thickness();
+  double layerRIn;
+  double layerROut;
+  double layerPosZ;
+
+  // Loop over layers
+  for(int i = 0; i < repeat; i++) {
+    layerRIn  = RIn + 1.0 * cm + ((double)i * 0.5) * cm;
+    layerROut = ROut;//RIn + ((double)i * 0.5) * cm;
+    layerPosZ = SizeZ / 2.0 - 5.0 * cm - ((double)i * 3.0) * cm;
+    layerSizeZ = 1.0 * cm;
+
+    string logic_layer_name = detName + _toString(i, "_Logic_lay_%d");
+    Volume layerVol(logic_layer_name,Tube(layerRIn, layerROut, layerSizeZ / 2.0, 0.0, 360.0 * deg), slice_mat);
+    layerVol.setVisAttributes(desc,x_layer.visStr());
+    sens.setType("tracker");
+    layerVol.setSensitiveDetector(sens);
+
+    Position     layer_pos = Position(0.0, 0.0, layerPosZ);
+    PlacedVolume layerPV = detVol.placeVolume(layerVol, layer_pos);
+    layerPV.addPhysVolID("layer", i+1);
+  }
+
+  DetElement   det(detName, detID);
+  Volume       motherVol = desc.pickMotherVolume(det);
+  Transform3D  tr(RotationZYX(0.0, 0.0, 0.0), Position(pos.x(), pos.y(), pos.z()));
+  PlacedVolume detPV = motherVol.placeVolume(detVol, tr);
+  detPV.addPhysVolID("system", detID);
+  det.setPlacement(detPV);
+  return det;
+}
+
+DECLARE_DETELEMENT(ce_GEM, createDetector)

src/deprecated/BarrelCalorimeterHybrid_geo.cpp

+//==========================================================================
+//  AIDA Detector description implementation
+//--------------------------------------------------------------------------
+// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
+// All rights reserved.
+//
+// For the licensing terms see $DD4hepINSTALL/LICENSE.
+// For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
+//
+// Author     : M.Frank
+//
+//==========================================================================
+//
+// Specialized generic detector constructor
+//
+//==========================================================================
+//
+// Implementation of the Sci Fiber geometry: M. Żurek 07/19/2021
+#include "DD4hep/DetFactoryHelper.h"
+#include "XML/Layering.h"
+#include "Math/Point2D.h"
+#include "TGeoPolygon.h"
+#include "TMath.h"
+
+using namespace std;
+using namespace dd4hep;
+using namespace dd4hep::detail;
+
+typedef ROOT::Math::XYPoint Point;
+
+// 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 = 1e-2) {
+  // 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;
+
+}
+
+// Create detector
+static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector sens)  {
+  static double tolerance = 0e0;
+  Layering      layering (e);
+  xml_det_t     x_det     = e;
+  Material      air       = description.air();
+  int           det_id    = x_det.id();
+  string        det_name  = x_det.nameStr();
+  xml_comp_t    x_staves  = x_det.staves();
+  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;
+  double        support_thickness = 0.0;
+  if(x_staves.hasChild("support")){
+    support_thickness = getAttrOrDefault(x_staves.child(_U(support)), _U(thickness), 5.0 * cm);
+  }
+  double        mod_z     = layering.totalThickness() + support_thickness;
+  double        outer_r   = inner_r + mod_z;
+  double        totThick  = mod_z;
+  double        offset    = x_det.attr<double>(_Unicode(offset));
+  DetElement    sdet      (det_name,det_id);
+  Volume        motherVol = description.pickMotherVolume(sdet);
+  PolyhedraRegular hedra  (nsides,inner_r,inner_r+totThick+tolerance*2e0,x_dim.z());
+  Volume        envelope  (det_name,hedra,air);
+  PlacedVolume  env_phv   = motherVol.placeVolume(envelope,Transform3D(Translation3D(0,0,offset)*RotationZ(M_PI/nsides)));
+
+  env_phv.addPhysVolID("system",det_id);
+  sdet.setPlacement(env_phv);
+
+  DetElement    stave_det("stave0",det_id);
+  double dx = 0.0; //mod_z / std::sin(dphi); // dx per layer
+
+  // Compute the top and bottom face measurements.
+  double trd_x2 = (2 * std::tan(hphi) * outer_r - dx)/2 - tolerance;
+  double trd_x1 = (2 * std::tan(hphi) * inner_r + dx)/2 - tolerance;
+  double trd_y1 = x_dim.z()/2 - tolerance;
+  double trd_y2 = trd_y1;
+  double trd_z  = mod_z/2 - tolerance;
+
+  // Create the trapezoid for the stave.
+  Trapezoid trd(trd_x1, // Outer side, i.e. the "long" X side.
+                trd_x2, // Inner side, i.e. the "short"  X side.
+                trd_y1, // Corresponds to subdetector (or module) Z.
+                trd_y2, //
+                trd_z); // Thickness, in Y for top stave, when rotated.
+
+  Volume mod_vol("stave",trd,air);
+  double l_pos_z = -(layering.totalThickness() / 2) - support_thickness/2.0;
+
+  //double trd_x2_support = trd_x1;
+  double trd_x1_support = (2 * std::tan(hphi) * outer_r - dx- support_thickness)/2 - tolerance;
+
+  Solid  support_frame_s;
+  // optional stave support
+  if(x_staves.hasChild("support")){
+    xml_comp_t x_support         = x_staves.child(_U(support));
+    // is the support on the inside surface?
+    bool       is_inside_support = getAttrOrDefault<bool>(x_support, _Unicode(inside), true);
+    // number of "beams" running the length of the stave.
+    int    n_beams        = getAttrOrDefault<int>(x_support, _Unicode(n_beams), 3);
+    double beam_thickness = support_thickness / 4.0; // maybe a parameter later...
+    trd_x1_support        = (2 * std::tan(hphi) * (outer_r - support_thickness + beam_thickness)) / 2 - tolerance;
+    double grid_size      = getAttrOrDefault(x_support, _Unicode(grid_size), 25.0 * cm);
+    double beam_width     = 2.0 * trd_x1_support / (n_beams + 1); // quick hack to make some gap between T beams
+
+    double cross_beam_thickness    = support_thickness/4.0;
+    //double trd_x1_support    = (2 * std::tan(hphi) * (inner_r + beam_thickness)) / 2 - tolerance;
+    double trd_x2_support = trd_x2;
+
+    int n_cross_supports = std::floor((trd_y1-cross_beam_thickness)/grid_size);
+
+    Box        beam_vert_s(beam_thickness / 2.0 - tolerance, trd_y1, support_thickness / 2.0 - tolerance);
+    Box        beam_hori_s(beam_width / 2.0 - tolerance, trd_y1, beam_thickness / 2.0 - tolerance);
+    UnionSolid T_beam_s(beam_vert_s, beam_hori_s, Position(0, 0, -support_thickness / 2.0 + beam_thickness / 2.0));
+
+    // cross supports
+    Trapezoid  trd_support(trd_x1_support,trd_x2_support,
+                           beam_thickness / 2.0 - tolerance, beam_thickness / 2.0 - tolerance,
+                          support_thickness / 2.0 - tolerance - cross_beam_thickness/2.0);
+    UnionSolid support_array_start_s(T_beam_s,trd_support,Position(0,0,cross_beam_thickness/2.0));
+    for (int isup = 0; isup < n_cross_supports; isup++) {
+      support_array_start_s = UnionSolid(support_array_start_s, trd_support, Position(0, -1.0 * isup * grid_size, cross_beam_thickness/2.0));
+      support_array_start_s = UnionSolid(support_array_start_s, trd_support, Position(0, 1.0 * isup * grid_size, cross_beam_thickness/2.0));
+    }
+    support_array_start_s =
+        UnionSolid(support_array_start_s, beam_hori_s,
+                   Position(-1.8 * 0.5*(trd_x1+trd_x2_support) / n_beams, 0, -support_thickness / 2.0 + beam_thickness / 2.0));
+    support_array_start_s =
+        UnionSolid(support_array_start_s, beam_hori_s,
+                   Position(1.8 * 0.5*(trd_x1+trd_x2_support) / n_beams, 0, -support_thickness / 2.0 + beam_thickness / 2.0));
+    support_array_start_s =
+        UnionSolid(support_array_start_s, beam_vert_s, Position(-1.8 * 0.5*(trd_x1+trd_x2_support) / n_beams, 0, 0));
+    support_array_start_s =
+        UnionSolid(support_array_start_s, beam_vert_s, Position(1.8 * 0.5*(trd_x1+trd_x2_support) / n_beams, 0, 0));
+
+    support_frame_s = support_array_start_s;
+
+    Material support_mat = description.material(x_support.materialStr());
+    Volume   support_vol("support_frame_v", support_frame_s, support_mat);
+    support_vol.setVisAttributes(description,x_support.visStr());
+
+    // figure out how to best place
+    //auto pv = mod_vol.placeVolume(support_vol, Position(0.0, 0.0, l_pos_z + support_thickness / 2.0));
+    auto pv = mod_vol.placeVolume(support_vol, Position(0.0, 0.0, -l_pos_z - support_thickness / 2.0));
+  }
+  //l_pos_z += support_thickness;
+
+  sens.setType("calorimeter");
+  { // =====  buildBarrelStave(description, sens, module_volume) =====
+    // Parameters for computing the layer X dimension:
+    double stave_z  = trd_y1;
+    double tan_hphi = std::tan(hphi);
+    double l_dim_x  = trd_x1; // Starting X dimension for the layer.
+
+    // 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();
+      // Loop over number of repeats for this layer.
+      for (int j=0; j<repeat; j++)    {
+        string l_name = _toString(l_num,"layer%d");
+        double l_thickness = layering.layer(l_num-1)->thickness();  // Layer's thickness.
+
+        Position   l_pos(0,0,l_pos_z+l_thickness/2);      // Position of the layer.
+	      double l_trd_x1 = l_dim_x - tolerance;
+	      double l_trd_x2 = l_dim_x + l_thickness*tan_hphi - tolerance;
+	      double l_trd_y1 = stave_z-tolerance;
+	      double l_trd_y2 = l_trd_y1;
+	      double l_trd_z  = l_thickness/2-tolerance;
+
+        Trapezoid  l_trd(l_trd_x1,l_trd_x2,l_trd_y1,l_trd_y2,l_trd_z);
+        Volume     l_vol(l_name,l_trd,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  = _toString(s_num,"slice%d");
+          double     s_thick = x_slice.thickness();
+          Volume     s_vol(s_name);
+          DetElement slice(layer,s_name,det_id);
+
+	          double s_trd_x1 = l_dim_x + (s_pos_z+l_thickness/2)*tan_hphi - tolerance;
+	          double s_trd_x2 = l_dim_x + (s_pos_z+l_thickness/2+s_thick)*tan_hphi - tolerance;
+	          double s_trd_y1 = stave_z-tolerance;
+	          double s_trd_y2 = s_trd_y1;
+	          double s_trd_z  = s_thick/2-tolerance;
+
+
+            Trapezoid  s_trd(s_trd_x1, s_trd_x2, s_trd_y1, s_trd_y2, s_trd_z);
+	          s_vol.setSolid(s_trd);
+	          s_vol.setMaterial(description.material(x_slice.materialStr()));
+
+
+          if (x_slice.hasChild("fiber")) {
+            xml_comp_t x_fiber = x_slice.child(_Unicode(fiber));
+            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");
+
+            // Calculate fiber positions inside the slice
+            Tube f_tube(0, f_radius, stave_z-tolerance);
+
+            // 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
+            pair<int, int> grid_div = getNdivisions(s_trd_x1, s_thick-tolerance, 2.0*cm, 2.0*cm);
+            // Calculate polygonal grid coordinates (vertices)
+            vector<tuple<int, Point, Point, Point, Point>> grid_vtx = gridPoints(grid_div.first, grid_div.second, s_trd_x1, s_thick-tolerance, hphi);
+
+            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-tolerance, hphi);
+            for (auto &line : f_pos) {
+              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) {
+                  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]++;
+                  }
+                }
+
+                string f_name = "fiber" + to_string(f_grid_id) + "_" + to_string(f_id);
+                Volume f_vol(f_name, f_tube, description.material(x_fiber.materialStr()));
+                DetElement fiber(slice, f_name, det_id);
+                if ( x_fiber.isSensitive() ) {
+                  f_vol.setSensitiveDetector(sens);
+                }
+                fiber.setAttributes(description,f_vol,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, f_tr);
+                fiber_phv.addPhysVolID(f_id_grid, f_grid_id + 1).addPhysVolID(f_id_fiber, f_id + 1);
+                fiber.setPlacement(fiber_phv);
+              }
+	        }
+          }
+
+          if ( x_slice.isSensitive() ) {
+            s_vol.setSensitiveDetector(sens);
+          }
+
+          slice.setAttributes(description,s_vol,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;
+
+          // Increment slice number.
+          ++s_num;
+        }
+
+        // Set region, limitset, and vis of layer.
+        layer.setAttributes(description,l_vol,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.
+        double xcut = l_thickness * tan_hphi;
+        l_dim_x += xcut;
+        l_pos_z += l_thickness;
+        ++l_num;
+      }
+    }
+  }
+
+  // Set stave visualization.
+  if ( x_staves )   {
+    mod_vol.setVisAttributes(description.visAttributes(x_staves.visStr()));
+  }
+  // Phi start for a stave.
+  double phi = M_PI / nsides;
+  double mod_x_off = dx / 2;             // Stave X offset, derived from the dx.
+  double mod_y_off = inner_r + mod_z/2;  // Stave Y offset
+
+  // Create nsides staves.
+  for (int i = 0; i < nsides; i++, phi -= dphi)      { // i is module number
+    // Compute the stave position
+    double m_pos_x = mod_x_off * std::cos(phi) - mod_y_off * std::sin(phi);
+    double m_pos_y = mod_x_off * std::sin(phi) + mod_y_off * std::cos(phi);
+    Transform3D tr(RotationZYX(0,phi,M_PI*0.5),Translation3D(-m_pos_x,-m_pos_y,0));
+    PlacedVolume pv = envelope.placeVolume(mod_vol,tr);
+    pv.addPhysVolID("system",det_id);
+    pv.addPhysVolID("module",i+1);
+    DetElement sd = i==0 ? stave_det : stave_det.clone(_toString(i,"stave%d"));
+    sd.setPlacement(pv);
+    sdet.add(sd);
+  }
+
+  // Set envelope volume attributes.
+  envelope.setAttributes(description,x_det.regionStr(),x_det.limitsStr(),x_det.visStr());
+  return sdet;
+}
+
+DECLARE_DETELEMENT(athena_EcalBarrelHybrid,create_detector)