scripts/view14/generate_eps

@@ -69,22 +69,6 @@ set -- "${POSITIONAL[@]}" # restore positional parameters
 
 # units are mm
 
-echo "view14 produces a series of slightly rotated XY slices a different z locations. Along beamline"
-
-# slice at z = 2m
-# note the offset has to change with sign of the direction to cut in the opposite direction.
-dawncut 0  0  1  10005 ${INPUT_FILE} ${FILE_TAG}b_temp0.prim 
-dawncut 0  0 -1 -1000 ${FILE_TAG}b_temp0.prim  ${FILE_TAG}b.prim
-dawn -d ${FILE_TAG}b.prim 
-ps2pdf ${FILE_TAG}b.eps ${FILE_TAG}b_full.pdf
-gs -o ${FILE_TAG}b.pdf -sDEVICE=pdfwrite \
-  -c "[/CropBox [50 175 550 675] /PAGES pdfmark" \
-  -f ${FILE_TAG}b_full.pdf
-pdftoppm ${FILE_TAG}b.pdf ${FILE_TAG}b -png -singlefile -cropbox
-
-echo "done ..."
-
-
 original_file_tag="${FILE_TAG}"
 
 make_slice(){

scripts/view15/.DAWN_1.history

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

scripts/view2/.DAWN_1.history

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

scripts/view20/.DAWN_1.history

 1.34392e+07
-0
+90
 180
 0
 0
 0
 0
-5
-3
+1
+1
 0.001
 0
 1

scripts/view20/generate_eps

@@ -47,7 +47,7 @@ set -- "${POSITIONAL[@]}" # restore positional parameters
 # Top view with a thin slice down the middle
 dawncut 0 1 0 1 ${INPUT_FILE} ${FILE_TAG}_top_temp0.prim   
 dawncut 0 -1 0 1 ${FILE_TAG}_top_temp0.prim  ${FILE_TAG}_top.prim
-../../bin/dawn_tweak --mag 2 --draw 1 --theta 270 --phi 90
+../../bin/dawn_tweak --mag 2 --draw 1 --theta 90 --phi 90
 dawn -d ${FILE_TAG}_top.prim 
 ps2pdf ${FILE_TAG}_top.eps ${FILE_TAG}_top_full.pdf
 gs -o ${FILE_TAG}_top.pdf -sDEVICE=pdfwrite \
@@ -68,13 +68,18 @@ pdftoppm ${FILE_TAG}.pdf ${FILE_TAG} -png -singlefile -cropbox -thinlinemode sol
 
 #npdet_info print EcalEndcapN_z0 --value-only ../../athena.xml
 #180.5 cm
-NMOD1=`npdet_info print EcalEndcapN_NModules_Sector1 ${DETECTOR_PATH}/calorimeters.xml | tail -1`
-NMOD2=`npdet_info print EcalEndcapN_NModules_Sector2 ${DETECTOR_PATH}/calorimeters.xml | tail -1`
+zcut=$(npdet_info print EcalEndcapN_z0 --value-only ${DETECTOR_PATH}/athena.xml )
+NMOD1=$(npdet_info print EcalEndcapN_NModules_Sector1 --value-only ${DETECTOR_PATH}/calorimeters.xml )
+NMOD2=$(npdet_info print EcalEndcapN_NModules_Sector2 --value-only ${DETECTOR_PATH}/calorimeters.xml )
+
+echo "NMOD1 = ${NMOD1}"
+echo "NMOD2 = ${NMOD2}"
+echo "zcut  = ${zcut}"
 
 # Top view with a thin slice down the middle
 dawncut 0 0  1 -1800 ${INPUT_FILE} ${FILE_TAG}_endcapN_temp0.prim   
-dawncut 0 0 -1 1805 ${FILE_TAG}_endcapN_temp0.prim  ${FILE_TAG}_endcapN.prim
-../../bin/dawn_tweak --mag 5 --draw 1 --theta 180 --phi 0 --draw 3
+dawncut 0 0 -1 2200 ${FILE_TAG}_endcapN_temp0.prim  ${FILE_TAG}_endcapN.prim
+../../bin/dawn_tweak --mag 5 --draw 3 --theta 180 --phi 0
 dawn -d ${FILE_TAG}_endcapN.prim 
 ps2pdf ${FILE_TAG}_endcapN.eps ${FILE_TAG}_endcapN_full.pdf
 gs -o ${FILE_TAG}_endcapN.pdf -sDEVICE=pdfwrite \

scripts/view5/.DAWN_1.history

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

scripts/view6/.DAWN_1.history

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

src/BarrelBarDetectorWithSideFrame_geo.cpp

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

src/BarrelCalorimeterInterlayers_geo.cpp

@@ -21,7 +21,7 @@ using namespace dd4hep;
 using namespace dd4hep::detail;
 
 typedef ROOT::Math::XYPoint Point;
-// fiber placement helpers, defined in BarrelCalorimeterHybrid_geo
+// fiber placement helpers, defined below
 vector<vector<Point>> fiberPositions(double radius, double x_spacing, double z_spacing,
                                      double x, double z, double phi, double spacing_tol = 1e-2);
 std::pair<int, int> getNdivisions(double x, double z, double dx, double dz);
@@ -277,23 +277,22 @@ void buildSupport(Detector& desc, Volume &mod_vol, xml_comp_t x_support,
   double beam_width       = 2. * trd_x1_support / (n_beams + 1); // quick hack to make some gap between T beams
   double beam_gap         = getAttrOrDefault(x_support, _Unicode(beam_gap), 3.*cm);
 
-  // build H-shape beam
+  // build T-shape beam
   double beam_space_x    = beam_width + beam_gap;
   double beam_space_z    = support_thickness - beam_thickness;
-  double cross_thickness = support_thickness - 2.*beam_thickness;
-  double beam_pos_z      = -beam_thickness / 2.;
-  double beam_center_z   = support_thickness / 2. + beam_pos_z;
+  double cross_thickness = support_thickness - beam_thickness;
+  double beam_pos_z      = beam_thickness / 2.;
+  double beam_center_z   = support_thickness / 2. - beam_pos_z;
 
-  Box        beam_vert_s(beam_thickness / 2. , trd_y, cross_thickness / 2.);
+  Box        beam_vert_s(beam_thickness / 2., trd_y, cross_thickness / 2.);
   Box        beam_hori_s(beam_width / 2., trd_y, beam_thickness / 2.);
-  UnionSolid T_beam_s(beam_hori_s, beam_vert_s, Position(0., 0., beam_space_z / 2.));
-  UnionSolid H_beam_s(T_beam_s, beam_hori_s, Position(0., 0., support_thickness - beam_thickness));
-  Volume H_beam_vol("H_beam", H_beam_s, desc.material(x_support.materialStr()));
+  UnionSolid T_beam_s(beam_hori_s, beam_vert_s, Position(0., 0., support_thickness / 2.));
+  Volume H_beam_vol("H_beam", T_beam_s, desc.material(x_support.materialStr()));
   H_beam_vol.setVisAttributes(desc, x_support.visStr());
   // place H beams first
   double beam_start_x = - (n_beams - 1) * (beam_width + beam_gap) / 2.;
   for (int i = 0; i < n_beams; ++i) {
-    Position beam_pos(beam_start_x + i * (beam_width + beam_gap), 0., - support_thickness / 2. - beam_pos_z);
+    Position beam_pos(beam_start_x + i * (beam_width + beam_gap), 0., - support_thickness / 2. + beam_pos_z);
     env_vol.placeVolume(H_beam_vol, beam_pos);
   }
 
@@ -303,10 +302,10 @@ void buildSupport(Detector& desc, Volume &mod_vol, xml_comp_t x_support,
   Volume cross_vol("cross_center_beam", cross_s, desc.material(x_support.materialStr()));
   cross_vol.setVisAttributes(desc, x_support.visStr());
   for (int i = 0; i < n_beams - 1; ++i) {
-    env_vol.placeVolume(cross_vol, Position(beam_start_x + beam_space_x * (i + 0.5), 0., 0.));
+    env_vol.placeVolume(cross_vol, Position(beam_start_x + beam_space_x * (i + 0.5), 0., beam_pos_z));
     for (int j = 1; j < n_cross_supports; j++) {
-      env_vol.placeVolume(cross_vol, Position(beam_start_x + beam_space_x * (i + 0.5), -j * grid_size, 0.));
-      env_vol.placeVolume(cross_vol, Position(beam_start_x + beam_space_x * (i + 0.5), j * grid_size, 0.));
+      env_vol.placeVolume(cross_vol, Position(beam_start_x + beam_space_x * (i + 0.5), -j * grid_size, beam_pos_z));
+      env_vol.placeVolume(cross_vol, Position(beam_start_x + beam_space_x * (i + 0.5), j * grid_size, beam_pos_z));
     }
   }
 
@@ -315,24 +314,135 @@ void buildSupport(Detector& desc, Volume &mod_vol, xml_comp_t x_support,
   double    cross_edge_x = trd_x1_support + beam_start_x - beam_thickness / 2.;
   double    cross_trd_x1 = cross_edge_x + std::tan(hphi) * beam_thickness;
   double    cross_trd_x2 = cross_trd_x1 + 2.* std::tan(hphi) * cross_thickness;
-  double    edge_pos_x   = beam_start_x - beam_thickness / 2. - cross_trd_x1 / 2.;
+  double    edge_pos_x   = beam_start_x - cross_trd_x1 / 2. - beam_thickness / 2;
   Trapezoid cross_s2_trd (cross_trd_x1 / 2., cross_trd_x2 / 2.,
                           beam_thickness / 2., beam_thickness / 2., cross_thickness / 2.);
-  Box       cross_s2_box ((cross_trd_x2 - cross_trd_x1)/2., beam_thickness / 2., cross_thickness / 2.);
-  SubtractionSolid cross_s2(cross_s2_trd, cross_s2_box, Position((cross_trd_x1 + cross_trd_x2) / 4., 0., 0.));
+  Box       cross_s2_box ((cross_trd_x2 - cross_trd_x1)/4., beam_thickness / 2., cross_thickness / 2.);
+  SubtractionSolid cross_s2(cross_s2_trd, cross_s2_box, Position((cross_trd_x2 + cross_trd_x1)/4., 0., 0.));
   Volume cross_vol2("cross_edge_beam", cross_s2, desc.material(x_support.materialStr()));
   cross_vol2.setVisAttributes(desc, x_support.visStr());
-  env_vol.placeVolume(cross_vol2, Position(edge_pos_x, 0., 0.));
-  env_vol.placeVolume(cross_vol2, Transform3D(Translation3D(-edge_pos_x, 0., 0.) * RotationZ(M_PI)));
+  env_vol.placeVolume(cross_vol2, Position(edge_pos_x, 0., beam_pos_z));
+  env_vol.placeVolume(cross_vol2, Transform3D(Translation3D(-edge_pos_x, 0., beam_pos_z) * RotationZ(M_PI)));
   for (int j = 1; j < n_cross_supports; j++) {
-    env_vol.placeVolume(cross_vol2, Position(edge_pos_x, -j * grid_size, 0.));
-    env_vol.placeVolume(cross_vol2, Position(edge_pos_x, j * grid_size, 0.));
-    env_vol.placeVolume(cross_vol2, Transform3D(Translation3D(-edge_pos_x, -j * grid_size, 0.) * RotationZ(M_PI)));
-    env_vol.placeVolume(cross_vol2, Transform3D(Translation3D(-edge_pos_x, j * grid_size, 0.) * RotationZ(M_PI)));
+    env_vol.placeVolume(cross_vol2, Position(edge_pos_x, -j * grid_size, beam_pos_z));
+    env_vol.placeVolume(cross_vol2, Position(edge_pos_x, j * grid_size, beam_pos_z));
+    env_vol.placeVolume(cross_vol2, Transform3D(Translation3D(-edge_pos_x, -j * grid_size, beam_pos_z) * RotationZ(M_PI)));
+    env_vol.placeVolume(cross_vol2, Transform3D(Translation3D(-edge_pos_x, j * grid_size, beam_pos_z) * RotationZ(M_PI)));
   }
 
   mod_vol.placeVolume(env_vol, Position(0.0, 0.0, l_pos_z + support_thickness/2.));
 }
 
+// Fill fiber lattice into trapezoid starting from position (0,0) in x-z coordinate system
+vector<vector<Point>> fiberPositions(double radius, double x_spacing, double z_spacing, double x, double z, double phi,
+                                     double spacing_tol)
+{
+  // z_spacing - distance between fiber layers in z
+  // x_spacing - distance between fiber centers in x
+  // x - half-length of the shorter (bottom) base of the trapezoid
+  // z - height of the trapezoid
+  // phi - angle between z and trapezoid arm
+
+  vector<vector<Point>> positions;
+  int z_layers = floor((z / 2 - radius - spacing_tol) / z_spacing); // number of layers that fit in z/2
+
+  double z_pos = 0.;
+  double x_pos = 0.;
+
+  for (int l = -z_layers; l < z_layers + 1; l++) {
+    vector<Point> xline;
+    z_pos                = l * z_spacing;
+    double x_max         = x + (z / 2. + z_pos) * tan(phi) - spacing_tol; // calculate max x at particular z_pos
+    (l % 2 == 0) ? x_pos = 0. : x_pos = x_spacing / 2;                    // account for spacing/2 shift
+
+    while (x_pos < (x_max - radius)) {
+      xline.push_back(Point(x_pos, z_pos));
+      if (x_pos != 0.)
+        xline.push_back(Point(-x_pos, z_pos)); // using symmetry around x=0
+      x_pos += x_spacing;
+    }
+    // Sort fiber IDs for a better organization
+    sort(xline.begin(), xline.end(), [](const Point& p1, const Point& p2) { return p1.x() < p2.x(); });
+    positions.emplace_back(std::move(xline));
+  }
+  return positions;
+}
+
+// Calculate number of divisions for the readout grid for the fiber layers
+std::pair<int, int> getNdivisions(double x, double z, double dx, double dz)
+{
+  // x and z defined as in vector<Point> fiberPositions
+  // dx, dz - size of the grid in x and z we want to get close to with the polygons
+  // See also descripltion when the function is called
+
+  double SiPMsize = 13.0 * mm;
+  double grid_min = SiPMsize + 3.0 * mm;
+
+  if (dz < grid_min) {
+    dz = grid_min;
+  }
+
+  if (dx < grid_min) {
+    dx = grid_min;
+  }
+
+  int nfit_cells_z = floor(z / dz);
+  int n_cells_z    = nfit_cells_z;
+
+  if (nfit_cells_z == 0)
+    n_cells_z++;
+
+  int nfit_cells_x = floor((2 * x) / dx);
+  int n_cells_x    = nfit_cells_x;
+
+  if (nfit_cells_x == 0)
+    n_cells_x++;
+
+  return std::make_pair(n_cells_x, n_cells_z);
+}
+
+// Calculate dimensions of the polygonal grid in the cartesian coordinate system x-z
+vector<tuple<int, Point, Point, Point, Point>> gridPoints(int div_x, int div_z, double x, double z, double phi)
+{
+  // x, z and phi defined as in vector<Point> fiberPositions
+  // div_x, div_z - number of divisions in x and z
+  double dz = z / div_z;
+
+  std::vector<std::tuple<int, Point, Point, Point, Point>> points;
+
+  for (int iz = 0; iz < div_z + 1; iz++) {
+    for (int ix = 0; ix < div_x + 1; ix++) {
+      double A_z = -z / 2 + iz * dz;
+      double B_z = -z / 2 + (iz + 1) * dz;
+
+      double len_x_for_z        = 2 * (x + iz * dz * tan(phi));
+      double len_x_for_z_plus_1 = 2 * (x + (iz + 1) * dz * tan(phi));
+
+      double dx_for_z        = len_x_for_z / div_x;
+      double dx_for_z_plus_1 = len_x_for_z_plus_1 / div_x;
+
+      double A_x = -len_x_for_z / 2. + ix * dx_for_z;
+      double B_x = -len_x_for_z_plus_1 / 2. + ix * dx_for_z_plus_1;
+
+      double C_z = B_z;
+      double D_z = A_z;
+      double C_x = B_x + dx_for_z_plus_1;
+      double D_x = A_x + dx_for_z;
+
+      int id = ix + div_x * iz;
+
+      auto A = Point(A_x, A_z);
+      auto B = Point(B_x, B_z);
+      auto C = Point(C_x, C_z);
+      auto D = Point(D_x, D_z);
+
+      // vertex points filled in the clock-wise direction
+      points.push_back(make_tuple(id, A, B, C, D));
+    }
+  }
+
+  return points;
+}
+
 DECLARE_DETELEMENT(athena_EcalBarrelInterlayers, create_detector)
 

src/BarrelTrackerWithFrame_geo.cpp

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

src/CompositeTracker_geo.cpp

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

src/CylinderTrackerBarrel_geo.cpp

@@ -2,11 +2,17 @@
 // Specialized generic detector constructor
 //==========================================================================
 
-#include "Acts/Plugins/DD4hep/ActsExtension.hpp"
-#include "Acts/Plugins/DD4hep/ConvertDD4hepMaterial.hpp"
 #include "DD4hep/DetFactoryHelper.h"
 #include "DD4hep/Printout.h"
 
+#if defined(USE_ACTSDD4HEP)
+#include "ActsDD4hep/ActsExtension.hpp"
+#include "ActsDD4hep/ConvertMaterial.hpp"
+#else
+#include "Acts/Plugins/DD4hep/ActsExtension.hpp"
+#include "Acts/Plugins/DD4hep/ConvertDD4hepMaterial.hpp"
+#endif
+
 using namespace std;
 using namespace dd4hep;
 using namespace dd4hep::detail;

src/DIRC_geo.cpp

@@ -162,7 +162,7 @@ static Ref_t createDetector(Detector& desc, xml_h e, SensitiveDetector sens)
   Volume lGlue("lGlue", gGlue, desc.material(xml_glue.materialStr()));
   lGlue.setVisAttributes(desc.visAttributes(xml_glue.visStr()));
 
-  sens.setType("photoncounter");
+  sens.setType("tracker");
   lBar.setSensitiveDetector(sens);
 
   int bars_repeat_z = 4;    // TODO parametrize!
@@ -320,7 +320,7 @@ static Ref_t createDetector(Detector& desc, xml_h e, SensitiveDetector sens)
   // Box    bar_geo(bar_thicknes / 2., bar_width / 2., bar_length / 2.);
   // Volume bar_volume("cb_DIRC_bars_Logix", bar_geo, bar_material);
   // bar_volume.setVisAttributes(desc.visAttributes(xml_det.visStr()));
-  // sens.setType("photoncounter");
+  // sens.setType("tracker");
   // bar_volume.setSensitiveDetector(sens);
 
   // for (int ia = 0; ia < bar_count; ia++) {

src/DRich_geo.cpp → src/DRICH_geo.cpp

@@ -4,7 +4,7 @@
 //
 // Author: Christopher Dilks (Duke University)
 //
-// - Design Adapted from Standalone Fun4all and GEMC implementations 
+// - Design Adapted from Standalone Fun4all and GEMC implementations
 //   [ Evaristo Cisbani, Cristiano Fanelli, Alessio Del Dotto, et al. ]
 //
 //==========================================================================
@@ -20,7 +20,6 @@
 #include "DD4hep/DetFactoryHelper.h"
 #include "DD4hep/Printout.h"
 
-using namespace std;
 using namespace dd4hep;
 using namespace dd4hep::rec;
 
@@ -35,88 +34,105 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
   OpticalSurfaceManager surfMgr = desc.surfaceManager();
 
   // attributes -----------------------------------------------------------
-  // TODO [low priority]: some attributes have default values, some do not;
-  //                      make this more consistent
   // - vessel
-  double vesselZ0 = dims.attr<double>(_Unicode(z0));
-  double vesselLength = dims.attr<double>(_Unicode(length));
-  double vesselRmin0 = dims.attr<double>(_Unicode(rmin0));
-  double vesselRmin1 = dims.attr<double>(_Unicode(rmin1));
-  double vesselRmax0 = dims.attr<double>(_Unicode(rmax0));
-  double vesselRmax1 = dims.attr<double>(_Unicode(rmax1));
-  double vesselRmax2 = dims.attr<double>(_Unicode(rmax2));
-  double snoutLength = dims.attr<double>(_Unicode(snout_length));
-  int nSectors = getAttrOrDefault<int>(dims, _Unicode(nsectors), 6);
-  double wallThickness = getAttrOrDefault<double>(dims, _Unicode(wall_thickness), 0.5*cm);
-  double windowThickness = getAttrOrDefault<double>(dims, _Unicode(window_thickness), 0.1*cm);
-  auto vesselMat = desc.material(getAttrOrDefault(detElem, _Unicode(material), "Aluminum"));
-  auto gasvolMat = desc.material(getAttrOrDefault(detElem, _Unicode(gas), "AirOptical"));
-  auto vesselVis = desc.visAttributes(detElem.attr<std::string>(_Unicode(vis_vessel)));
-  auto gasvolVis = desc.visAttributes(detElem.attr<std::string>(_Unicode(vis_gas)));
+  double  vesselZmin       =  dims.attr<double>(_Unicode(zmin));
+  double  vesselLength     =  dims.attr<double>(_Unicode(length));
+  double  vesselRmin0      =  dims.attr<double>(_Unicode(rmin0));
+  double  vesselRmin1      =  dims.attr<double>(_Unicode(rmin1));
+  double  vesselRmax0      =  dims.attr<double>(_Unicode(rmax0));
+  double  vesselRmax1      =  dims.attr<double>(_Unicode(rmax1));
+  double  vesselRmax2      =  dims.attr<double>(_Unicode(rmax2));
+  double  snoutLength      =  dims.attr<double>(_Unicode(snout_length));
+  int     nSectors         =  dims.attr<int>(_Unicode(nsectors));
+  double  wallThickness    =  dims.attr<double>(_Unicode(wall_thickness));
+  double  windowThickness  =  dims.attr<double>(_Unicode(window_thickness));
+  auto    vesselMat        =  desc.material(detElem.attr<std::string>(_Unicode(material)));
+  auto    gasvolMat        =  desc.material(detElem.attr<std::string>(_Unicode(gas)));
+  auto    vesselVis        =  desc.visAttributes(detElem.attr<std::string>(_Unicode(vis_vessel)));
+  auto    gasvolVis        =  desc.visAttributes(detElem.attr<std::string>(_Unicode(vis_gas)));
   // - radiator (applies to aerogel and filter)
-  auto radiatorElem = detElem.child(_Unicode(radiator));
-  double radiatorRmin = getAttrOrDefault<double>(radiatorElem, _Unicode(rmin), 10.*cm);
-  double radiatorRmax = getAttrOrDefault<double>(radiatorElem, _Unicode(rmax), 80.*cm);
-  double radiatorPhiw = getAttrOrDefault<double>(radiatorElem, _Unicode(phiw), 2*M_PI/nSectors);
-  double radiatorPitch = getAttrOrDefault<double>(radiatorElem, _Unicode(pitch), 0.*degree);
-  double radiatorFrontplane = getAttrOrDefault<double>(radiatorElem, _Unicode(frontplane), 2.5*cm);
+  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 = getAttrOrDefault<double>(aerogelElem, _Unicode(thickness), 2.5*cm);
+  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 = getAttrOrDefault<double>(filterElem, _Unicode(thickness), 2.5*cm);
+  auto    filterElem       =  radiatorElem.child(_Unicode(filter));
+  auto    filterMat        =  desc.material(filterElem.attr<std::string>(_Unicode(material)));
+  auto    filterVis        =  desc.visAttributes(filterElem.attr<std::string>(_Unicode(vis)));
+  double  filterThickness  =  filterElem.attr<double>(_Unicode(thickness));
   // - mirror
-  auto mirrorElem = detElem.child(_Unicode(mirror));
-  auto mirrorMat = desc.material(mirrorElem.attr<std::string>(_Unicode(material)));
-  auto mirrorVis = desc.visAttributes(mirrorElem.attr<std::string>(_Unicode(vis)));
-  auto mirrorSurf = surfMgr.opticalSurface(getAttrOrDefault(mirrorElem, _Unicode(surface), "MirrorOpticalSurface"));
-  double mirrorBackplane = getAttrOrDefault<double>(mirrorElem, _Unicode(backplane), 240.*cm);
-  double mirrorThickness = getAttrOrDefault<double>(mirrorElem, _Unicode(thickness), 2.*mm);
-  double mirrorRadius = getAttrOrDefault<double>(mirrorElem, _Unicode(radius), 190*cm);
-  double mirrorCenterX = getAttrOrDefault<double>(mirrorElem, _Unicode(centerx), 95*cm);
-  double mirrorRmin = getAttrOrDefault<double>(mirrorElem, _Unicode(rmin), 10.*cm);
-  double mirrorRmax = getAttrOrDefault<double>(mirrorElem, _Unicode(rmax), 150.*cm);
-  double mirrorPhiw = getAttrOrDefault<double>(mirrorElem, _Unicode(phiw), 2*M_PI/nSectors);
-  int mirrorDebug = getAttrOrDefault<int>(mirrorElem, _Unicode(debug), 0);
+  auto    mirrorElem       =  detElem.child(_Unicode(mirror));
+  auto    mirrorMat        =  desc.material(mirrorElem.attr<std::string>(_Unicode(material)));
+  auto    mirrorVis        =  desc.visAttributes(mirrorElem.attr<std::string>(_Unicode(vis)));
+  auto    mirrorSurf       =  surfMgr.opticalSurface(mirrorElem.attr<std::string>(_Unicode(surface)));
+  double  mirrorBackplane  =  mirrorElem.attr<double>(_Unicode(backplane));
+  double  mirrorThickness  =  mirrorElem.attr<double>(_Unicode(thickness));
+  double  mirrorRmin       =  mirrorElem.attr<double>(_Unicode(rmin));
+  double  mirrorRmax       =  mirrorElem.attr<double>(_Unicode(rmax));
+  double  mirrorPhiw       =  mirrorElem.attr<double>(_Unicode(phiw));
+  double  focusTuneZ       =  mirrorElem.attr<double>(_Unicode(focus_tune_z));
+  double  focusTuneX       =  mirrorElem.attr<double>(_Unicode(focus_tune_x));
   // - sensor module
-  auto sensorElem = detElem.child(_Unicode(sensors)).child(_Unicode(module));
-  auto sensorMat = desc.material(sensorElem.attr<std::string>(_Unicode(material)));
-  auto sensorVis = desc.visAttributes(sensorElem.attr<std::string>(_Unicode(vis)));
-  auto sensorSurf = surfMgr.opticalSurface(getAttrOrDefault(sensorElem, _Unicode(surface), "MirrorOpticalSurface"));
-  double sensorSide = sensorElem.attr<double>(_Unicode(side));
-  double sensorGap = sensorElem.attr<double>(_Unicode(gap));
-  double sensorThickness = sensorElem.attr<double>(_Unicode(thickness));
+  auto    sensorElem       =  detElem.child(_Unicode(sensors)).child(_Unicode(module));
+  auto    sensorMat        =  desc.material(sensorElem.attr<std::string>(_Unicode(material)));
+  auto    sensorVis        =  desc.visAttributes(sensorElem.attr<std::string>(_Unicode(vis)));
+  auto    sensorSurf       =  surfMgr.opticalSurface(sensorElem.attr<std::string>(_Unicode(surface)));
+  double  sensorSide       =  sensorElem.attr<double>(_Unicode(side));
+  double  sensorGap        =  sensorElem.attr<double>(_Unicode(gap));
+  double  sensorThickness  =  sensorElem.attr<double>(_Unicode(thickness));
   // - sensor sphere
-  auto sensorSphElem = detElem.child(_Unicode(sensors)).child(_Unicode(sphere));
-  double sensorSphRadius = sensorSphElem.attr<double>(_Unicode(radius));
-  double sensorSphCenterX = sensorSphElem.attr<double>(_Unicode(centerx));
-  double sensorSphCenterY = sensorSphElem.attr<double>(_Unicode(centery));
-  double sensorSphCenterZ = sensorSphElem.attr<double>(_Unicode(centerz));
-  int sensorSphDebug = getAttrOrDefault<int>(sensorSphElem, _Unicode(debug), 0);
+  auto    sensorSphElem     =  detElem.child(_Unicode(sensors)).child(_Unicode(sphere));
+  double  sensorSphRadius   =  sensorSphElem.attr<double>(_Unicode(radius));
+  double  sensorSphCenterX  =  sensorSphElem.attr<double>(_Unicode(centerx));
+  double  sensorSphCenterZ  =  sensorSphElem.attr<double>(_Unicode(centerz));
   // - sensor sphere patch cuts
-  auto sensorSphPatchElem = detElem.child(_Unicode(sensors)).child(_Unicode(sphericalpatch));
-  double sensorSphPatchThetaMin = sensorSphPatchElem.attr<double>(_Unicode(thetamin));
-  double sensorSphPatchThetaMax = sensorSphPatchElem.attr<double>(_Unicode(thetamax));
-  double sensorSphPatchWidthFactor = sensorSphPatchElem.attr<double>(_Unicode(widthfactor));
-  double sensorSphPatchTaper = sensorSphPatchElem.attr<double>(_Unicode(taper));
+  auto    sensorSphPatchElem  =  detElem.child(_Unicode(sensors)).child(_Unicode(sphericalpatch));
+  double  sensorSphPatchPhiw  =  sensorSphPatchElem.attr<double>(_Unicode(phiw));
+  double  sensorSphPatchRmin  =  sensorSphPatchElem.attr<double>(_Unicode(rmin));
+  double  sensorSphPatchRmax  =  sensorSphPatchElem.attr<double>(_Unicode(rmax));
+  double  sensorSphPatchZmin  =  sensorSphPatchElem.attr<double>(_Unicode(zmin));
+  // - debugging switches
+  int   debug_optics_mode  =  detElem.attr<int>(_Unicode(debug_optics));
+  bool  debug_mirror       =  mirrorElem.attr<bool>(_Unicode(debug));
+  bool  debug_sensors      =  sensorSphElem.attr<bool>(_Unicode(debug));
+
+  // if debugging optics, override some settings
+  bool debug_optics = debug_optics_mode > 0;
+  if(debug_optics) {
+    printout(WARNING,"DRich_geo","DEBUGGING DRICH OPTICS");
+    switch(debug_optics_mode) {
+      case 1: vesselMat = aerogelMat = filterMat = sensorMat = gasvolMat = desc.material("VacuumOptical"); break;
+      case 2: vesselMat = aerogelMat = filterMat = sensorMat = desc.material("VacuumOptical"); break;
+      default: printout(FATAL,"DRich_geo","UNKNOWN debug_optics_mode"); return det;
+    };
+    aerogelVis = sensorVis = mirrorVis;
+    gasvolVis = vesselVis = desc.invisible();
+  };
+
 
-  
   // BUILD VESSEL ====================================================================
   /* - `vessel`: aluminum enclosure, the mother volume of the dRICh
    * - `gasvol`: gas volume, which fills `vessel`; all other volumes defined below
    *   are children of `gasvol`
-   * - the dRICh vessel geometry has two regions: the snout refers to the conic region 
+   * - the dRICh vessel geometry has two regions: the snout refers to the conic region
    *   in the front, housing the aerogel, while the tank refers to the cylindrical
    *   region, housing the rest of the detector components
    */
 
+  // derived attributes
+  double tankLength = vesselLength - snoutLength;
+  double vesselZmax = vesselZmin + vesselLength;
+
   // snout solids
   double boreDelta = vesselRmin1 - vesselRmin0;
+  double snoutDelta = vesselRmax1 - vesselRmax0;
   Cone vesselSnout(
       snoutLength/2.0,
       vesselRmin0,
@@ -139,14 +155,14 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
 
   // tank solids
   Cone vesselTank(
-      (vesselLength - snoutLength)/2.0,
+      tankLength/2.0,
       vesselSnout.rMin2(),
       vesselRmax2,
       vesselRmin1,
       vesselRmax2
       );
   Cone gasvolTank(
-      (vesselLength - snoutLength)/2.0 - windowThickness,
+      tankLength/2.0 - windowThickness,
       gasvolSnout.rMin2(),
       vesselRmax2 - wallThickness,
       vesselRmin1 + wallThickness,
@@ -154,126 +170,135 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
       );
 
   // snout + tank solids
-  UnionSolid vesselSolid(
+  UnionSolid vesselUnion(
       vesselTank,
       vesselSnout,
       Position(0., 0., -vesselLength/2.)
       );
-  UnionSolid gasvolSolid(
+  UnionSolid gasvolUnion(
       gasvolTank,
       gasvolSnout,
       Position(0., 0., -vesselLength/2. + windowThickness)
       );
 
+  //  extra solids for `debug_optics` only
+  Box vesselBox(1001,1001,1001);
+  Box gasvolBox(1000,1000,1000);
+
+  // choose vessel and gasvol solids (depending on `debug_optics_mode` (0=disabled))
+  Solid vesselSolid, gasvolSolid;
+  switch(debug_optics_mode) {
+    case 0:  vesselSolid=vesselUnion;  gasvolSolid=gasvolUnion;  break; // `!debug_optics`
+    case 1:  vesselSolid=vesselBox;    gasvolSolid=gasvolBox;    break;
+    case 2:  vesselSolid=vesselBox;    gasvolSolid=gasvolUnion;  break;
+  };
 
   // volumes
   Volume vesselVol(detName, vesselSolid, vesselMat);
-  Volume gasvolVol(detName + "_gas", gasvolSolid, gasvolMat);
+  Volume gasvolVol(detName+"_gas", gasvolSolid, gasvolMat);
   vesselVol.setVisAttributes(vesselVis);
   gasvolVol.setVisAttributes(gasvolVis);
 
-  // reference position
-  auto snoutFront = Position(0., 0., -(vesselLength + snoutLength)/2.);
+  // reference positions
+  // - the vessel is created such that the center of the cylindrical tank volume
+  //   coincides with the origin; this is called the "origin position" of the vessel
+  // - when the vessel (and its children volumes) is placed, it is translated in
+  //   the z-direction to be in the proper ATHENA-integration location
+  // - these reference positions are for the frontplane and backplane of the vessel,
+  //   with respect to the vessel origin position
+  auto originFront = Position(0., 0., -tankLength/2.0 - snoutLength );
+  auto originBack =  Position(0., 0., tankLength/2.0 );
+
+  // initialize sensor centroids (used for mirror parameterization below); this is
+  // the average (x,y,z) of the placed sensors, w.r.t. originFront
+  double sensorCentroidX = 0;
+  double sensorCentroidZ = 0;
+  int sensorCount = 0;
 
 
   // sensitive detector type
-  sens.setType("photoncounter");
+  sens.setType("tracker");
 
 
-  // SECTOR LOOP //////////////////////////////////
-  for(int isec=0; isec<nSectors; isec++) {
+  // BUILD RADIATOR ====================================================================
 
-    if(mirrorDebug*isec>0 || sensorSphDebug*isec>0) continue; // if debugging, draw only 1 sector
-    double sectorRotation = isec * 360/nSectors * degree; // sector rotation about z axis
+  // attributes
+  double airGap = 0.01*mm; // air gap between aerogel and filter (FIXME? actually it's currently a gas gap)
 
-    // BUILD RADIATOR ====================================================================
-
-    // derived attributes
-    auto radiatorPos = Position(0., 0., radiatorFrontplane) + snoutFront;
+  // solid and volume: create aerogel and filter
+  Cone aerogelSolid(
+      aerogelThickness/2,
+      radiatorRmin,
+      radiatorRmax,
+      radiatorRmin + boreDelta  * aerogelThickness / vesselLength,
+      radiatorRmax + snoutDelta * aerogelThickness / snoutLength
+      );
+  Cone filterSolid(
+      filterThickness/2,
+      radiatorRmin + boreDelta  * (aerogelThickness + airGap) / vesselLength,
+      radiatorRmax + snoutDelta * (aerogelThickness + airGap) / snoutLength,
+      radiatorRmin + boreDelta  * (aerogelThickness + airGap + filterThickness) / vesselLength,
+      radiatorRmax + snoutDelta * (aerogelThickness + airGap + filterThickness) / snoutLength
+      );
 
-    // solid and volume: create aerogel and filter sectors
-    Tube aerogelSolid(radiatorRmin, radiatorRmax, aerogelThickness/2, -radiatorPhiw/2.0, radiatorPhiw/2.0);
-    Tube filterSolid( radiatorRmin, radiatorRmax, filterThickness/2,  -radiatorPhiw/2.0, radiatorPhiw/2.0);
-    Volume aerogelVol("aerogel_v", aerogelSolid, aerogelMat);
-    Volume filterVol( "filter_v",  filterSolid,  filterMat);
-    aerogelVol.setVisAttributes(aerogelVis);
-    filterVol.setVisAttributes(filterVis);
+  Volume aerogelVol( detName+"_aerogel", aerogelSolid, aerogelMat );
+  Volume filterVol(  detName+"_filter",  filterSolid,  filterMat );
+  aerogelVol.setVisAttributes(aerogelVis);
+  filterVol.setVisAttributes(filterVis);
+
+  // aerogel placement and surface properties
+  // TODO [low-priority]: define skin properties for aerogel and filter
+  auto radiatorPos = Position(0., 0., radiatorFrontplane) + originFront;
+  auto aerogelPV = gasvolVol.placeVolume(aerogelVol,
+        Translation3D(radiatorPos.x(), radiatorPos.y(), radiatorPos.z()) // re-center to originFront
+      * RotationY(radiatorPitch) // change polar angle to specified pitch // FIXME: probably broken (not yet in use anyway)
+      );
+  DetElement aerogelDE(det, "aerogel_de", 0);
+  aerogelDE.setPlacement(aerogelPV);
+  //SkinSurface aerogelSkin(desc, aerogelDE, "mirror_optical_surface", aerogelSurf, aerogelVol);
+  //aerogelSkin.isValid();
 
-    // placement
-    auto aerogelPV = gasvolVol.placeVolume(aerogelVol,
-	  RotationZ(sectorRotation) // rotate about beam axis to sector
-	* Translation3D(radiatorPos.x(), radiatorPos.y(), radiatorPos.z()) // re-center to snoutFront
-	* RotationY(radiatorPitch) // change polar angle to specified pitch
-	);
+  // filter placement and surface properties
+  if(!debug_optics) {
     auto filterPV = gasvolVol.placeVolume(filterVol,
-  	  RotationZ(sectorRotation) // rotate about beam axis to sector
-	* Translation3D(radiatorPos.x(), radiatorPos.y(), radiatorPos.z()) // re-center to snoutFront
-	* RotationY(radiatorPitch) // change polar angle
-	* Translation3D(0., 0., (aerogelThickness+filterThickness)/2.) // move to aerogel backplane
-	);
-
-    // properties
-    // TODO [critical]: define skin properties for aerogel and filter
-    DetElement aerogelDE(det, Form("aerogel_de%d", isec), isec);
-    aerogelDE.setPlacement(aerogelPV);
-    //SkinSurface aerogelSkin(desc, aerogelDE, Form("mirror_optical_surface%d", isec), aerogelSurf, aerogelVol);
-    //aerogelSkin.isValid();
-    DetElement filterDE(det, Form("filter_de%d", isec), isec);
+          Translation3D(0., 0., airGap) // add an air gap
+        * Translation3D(radiatorPos.x(), radiatorPos.y(), radiatorPos.z()) // re-center to originFront
+        * RotationY(radiatorPitch) // change polar angle
+        * Translation3D(0., 0., (aerogelThickness+filterThickness)/2.) // move to aerogel backplane
+        );
+    DetElement filterDE(det, "filter_de", 0);
     filterDE.setPlacement(filterPV);
-    //SkinSurface filterSkin(desc, filterDE, Form("mirror_optical_surface%d", isec), filterSurf, filterVol);
+    //SkinSurface filterSkin(desc, filterDE, "mirror_optical_surface", filterSurf, filterVol);
     //filterSkin.isValid();
+  };
 
 
-    // BUILD MIRROR ====================================================================
+  // SECTOR LOOP //////////////////////////////////
+  for(int isec=0; isec<nSectors; isec++) {
 
-    // derived attributes
-    auto mirrorPos = Position(mirrorCenterX, 0., mirrorBackplane) + snoutFront;
-    double mirrorThetaRot = std::asin(mirrorCenterX/mirrorRadius);
-    double mirrorTheta1 = mirrorThetaRot - std::asin((mirrorCenterX-mirrorRmin)/mirrorRadius);
-    double mirrorTheta2 = mirrorThetaRot + std::asin((mirrorRmax-mirrorCenterX)/mirrorRadius);
+    // debugging filters, limiting the number of sectors
+    if( (debug_mirror||debug_sensors||debug_optics) && isec!=0) continue;
 
-    // if debugging, draw full sphere
-    if(mirrorDebug>0) { mirrorTheta1=0; mirrorTheta2=M_PI; mirrorPhiw=2*M_PI; };
-
-    // solid and volume: create sphere at origin, with specified angular limits
-    Sphere mirrorSolid(
-	mirrorRadius-mirrorThickness,
-	mirrorRadius,
-	mirrorTheta1,
-	mirrorTheta2,
-	-mirrorPhiw/2.0,
-	mirrorPhiw/2.0
-	);
-    Volume mirrorVol("mirror_v", mirrorSolid, mirrorMat);
-    mirrorVol.setVisAttributes(mirrorVis);
+    // sector rotation about z axis
+    double sectorRotation = isec * 360/nSectors * degree;
+    std::string secName = "sec" + std::to_string(isec);
 
-    // placement (note: transformations are in reverse order)
-    auto mirrorPV = gasvolVol.placeVolume(mirrorVol,
-	  RotationZ(sectorRotation) // rotate about beam axis to sector
-	* Translation3D(0,0,-mirrorRadius) // move longitudinally so it intersects snoutFront
-	* Translation3D(mirrorPos.x(), mirrorPos.y(), mirrorPos.z()) // re-center to snoutFront
-	* RotationY(-mirrorThetaRot) // rotate about origin
-	);
-
-    // properties
-    DetElement mirrorDE(det, Form("mirror_de%d", isec), isec);
-    mirrorDE.setPlacement(mirrorPV);
-    SkinSurface mirrorSkin(desc, mirrorDE, Form("mirror_optical_surface%d", isec), mirrorSurf, mirrorVol);
-    mirrorSkin.isValid();
 
 
     // BUILD SENSORS ====================================================================
 
     // if debugging sphere properties, restrict number of sensors drawn
-    if(sensorSphDebug>0) { sensorSide = 2*M_PI*sensorSphRadius / 64; };
+    if(debug_sensors) { sensorSide = 2*M_PI*sensorSphRadius / 64; };
 
     // solid and volume: single sensor module
     Box sensorSolid(sensorSide/2., sensorSide/2., sensorThickness/2.);
-    Volume sensorVol("sensor_v", sensorSolid, sensorMat);
+    Volume sensorVol(detName+"_sensor_"+secName, sensorSolid, sensorMat);
     sensorVol.setVisAttributes(sensorVis);
 
+    auto sensorSphPos = Position(sensorSphCenterX, 0., sensorSphCenterZ) + originFront;
+
     // sensitivity
-    sensorVol.setSensitiveDetector(sens);
+    if(!debug_optics) sensorVol.setSensitiveDetector(sens);
 
     // SENSOR MODULE LOOP ------------------------
     /* ALGORITHM: generate sphere of positions
@@ -299,13 +324,13 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
     int imod=0;
 
     // thetaGen loop: iterate less than "0.5 circumference / sensor size" times
-    double nTheta = M_PI*sensorSphRadius / (sensorSide+sensorGap); 
+    double nTheta = M_PI*sensorSphRadius / (sensorSide+sensorGap);
     for(int t=0; t<(int)(nTheta+0.5); t++) {
       double thetaGen = t/((double)nTheta) * M_PI;
 
       // phiGen loop: iterate less than "circumference at this latitude / sensor size" times
-      double nPhi = 2*M_PI * sensorSphRadius * std::sin(thetaGen) / (sensorSide+sensorGap); 
-      for(int p=0; p<(int)(nPhi+0.5); p++) { 
+      double nPhi = 2*M_PI * sensorSphRadius * std::sin(thetaGen) / (sensorSide+sensorGap);
+      for(int p=0; p<(int)(nPhi+0.5); p++) {
         double phiGen = p/((double)nPhi) * 2*M_PI - M_PI; // shift to [-pi,pi]
 
         // determine global phi and theta
@@ -321,35 +346,35 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
         double phi = std::atan2(y,x);
         double theta = std::acos(z/sensorSphRadius);
 
-        // cut spherical patch
-        // TODO [low priority]: instead of cutting a patch with complicated parameters,
-        //                      can we use something simpler such as Union or
-        //                      Intersection of solids?
-        // - applied on global coordinates
-        // - theta cuts are signed, since we will offset the patch along +x;
-        //   from the offset position, toward barrel is positive, toward beam is negative
-        double thetaSigned = (x<0?-1:1) * theta;
-        // - position of yz planes, associated to theta cuts
-        double xmin = sensorSphRadius * std::sin(sensorSphPatchThetaMin/2);
-        double xmax = sensorSphRadius * std::sin(sensorSphPatchThetaMax/2);
-        // - we want a phi wedge, but offset from the origin to allow more width, so
-        //   define phiCheck to account for the offset; the amount of the offset,
-        //   and hence the width, is controlled by `sensorSphPatchWidthFactor`
-        double phiCheck = std::atan2(y,(x+sensorSphCenterX)/sensorSphPatchWidthFactor);
-        // - apply cuts (only if not debugging)
-	// - NOTE: use `x<xmax` for straight cuts, or `theta<sensorSphPatchThetaMax` for
-	//   rounded cuts (which allows for more sensors)
-        if( ( std::fabs(phiCheck)<sensorSphPatchTaper && x>xmin && theta<sensorSphPatchThetaMax && z>0 )
-            || sensorSphDebug>0
-        ) {
+        // shift global coordinates so we can apply spherical patch cuts
+        double zCheck = z + sensorSphCenterZ;
+        double xCheck = x + sensorSphCenterX;
+        double yCheck = y;
+        double rCheck = std::hypot(xCheck,yCheck);
+        double phiCheck = std::atan2(yCheck,xCheck);
+
+        // patch cut
+        bool patchCut =
+          std::fabs(phiCheck) < sensorSphPatchPhiw
+          && zCheck > sensorSphPatchZmin
+          && rCheck > sensorSphPatchRmin
+          && rCheck < sensorSphPatchRmax;
+        if(debug_sensors) patchCut = std::fabs(phiCheck) < sensorSphPatchPhiw;
+        if(patchCut) {
+
+          // append sensor position to centroid calculation
+          if(isec==0) {
+            sensorCentroidX += xCheck;
+            sensorCentroidZ += zCheck;
+            sensorCount++;
+          };
 
           // placement (note: transformations are in reverse order)
           // - transformations operate on global coordinates; the corresponding
           //   generator coordinates are provided in the comments
           auto sensorPV = gasvolVol.placeVolume(sensorVol,
                 RotationZ(sectorRotation) // rotate about beam axis to sector
-              * Translation3D(sensorSphCenterX, sensorSphCenterY, sensorSphCenterZ) // move sphere to specified center
-              * Translation3D(snoutFront.x(), snoutFront.y(), snoutFront.z()) // move sphere to reference position
+              * Translation3D(sensorSphPos.x(), sensorSphPos.y(), sensorSphPos.z()) // move sphere to reference position
               * RotationX(phiGen) // rotate about `zGen`
               * RotationZ(thetaGen) // rotate about `yGen`
               * Translation3D(sensorSphRadius, 0., 0.) // push radially to spherical surface
@@ -357,15 +382,17 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
               * RotationZ(-M_PI/2) // correction for readout segmentation mapping
               );
 
-	  // generate LUT for module number -> sensor position, for readout mapping tests
-	  //if(isec==0) printf("%d %f %f\n",imod,sensorPV.position().x(),sensorPV.position().y());
+          // generate LUT for module number -> sensor position, for readout mapping tests
+          //if(isec==0) printf("%d %f %f\n",imod,sensorPV.position().x(),sensorPV.position().y());
 
           // properties
           sensorPV.addPhysVolID("sector", isec).addPhysVolID("module", imod);
-          DetElement sensorDE(det, Form("sensor_de%d_%d", isec, imod), 10000*isec+imod);
+          DetElement sensorDE(det, Form("sensor_de%d_%d", isec, imod), (imod<<3)|isec ); // id must match IRTAlgorithm usage
           sensorDE.setPlacement(sensorPV);
-          SkinSurface sensorSkin(desc, sensorDE, Form("sensor_optical_surface%d", isec), sensorSurf, sensorVol);
-          sensorSkin.isValid();
+          if(!debug_optics) {
+            SkinSurface sensorSkin(desc, sensorDE, Form("sensor_optical_surface%d", isec), sensorSurf, sensorVol);
+            sensorSkin.isValid();
+          };
 
           // increment sensor module number
           imod++;
@@ -373,9 +400,132 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
         }; // end patch cuts
       }; // end phiGen loop
     }; // end thetaGen loop
+
+    // calculate centroid sensor position
+    if(isec==0) {
+      sensorCentroidX /= sensorCount;
+      sensorCentroidZ /= sensorCount;
+    };
+
     // END SENSOR MODULE LOOP ------------------------
 
 
+    // BUILD MIRRORS ====================================================================
+
+    // derive spherical mirror parameters `(zM,xM,rM)`, for given image point
+    // coordinates `(zI,xI)` and `dO`, defined as the z-distance between the
+    // object and the mirror surface
+    // - all coordinates are specified w.r.t. the object point coordinates
+    // - this is point-to-point focusing, but it can be used to effectively steer
+    //   parallel-to-point focusing
+    double zM,xM,rM;
+    auto FocusMirror = [&zM,&xM,&rM](double zI, double xI, double dO) {
+      zM = dO*zI / (2*dO-zI);
+      xM = dO*xI / (2*dO-zI);
+      rM = dO - zM;
+    };
+
+    // attributes, re-defined w.r.t. IP, needed for mirror positioning
+    double zS = sensorSphCenterZ + vesselZmin; // sensor sphere attributes
+    double xS = sensorSphCenterX;
+    double rS = sensorSphRadius;
+    double B = vesselZmax - mirrorBackplane; // distance between IP and mirror back plane
+
+    // focus 1: set mirror to focus IP on center of sensor sphere `(zS,xS)`
+    /*double zF = zS;
+    double xF = xS;
+    FocusMirror(zF,xF,B);*/
+
+    // focus 2: move focal region along sensor sphere radius, according to `focusTuneLong`
+    // - specifically, along the radial line which passes through the approximate centroid
+    //   of the sensor region `(sensorCentroidZ,sensorCentroidX)`
+    // - `focusTuneLong` is the distance to move, given as a fraction of `sensorSphRadius`
+    // - `focusTuneLong==0` means `(zF,xF)==(zS,xS)`
+    // - `focusTuneLong==1` means `(zF,xF)` will be on the sensor sphere, near the centroid
+    /*
+    double zC = sensorCentroidZ + vesselZmin;
+    double xC = sensorCentroidX;
+    double slopeF = (xC-xS) / (zC-zS);
+    double thetaF = std::atan(std::fabs(slopeF));
+    double zF = zS + focusTuneLong * sensorSphRadius * std::cos(thetaF);
+    double xF = xS - focusTuneLong * sensorSphRadius * std::sin(thetaF);
+    //FocusMirror(zF,xF,B);
+
+    // focus 3: move along line perpendicular to focus 2's radial line,
+    // according to `focusTunePerp`, with the same numerical scale as `focusTuneLong`
+    zF += focusTunePerp * sensorSphRadius * std::cos(M_PI/2-thetaF);
+    xF += focusTunePerp * sensorSphRadius * std::sin(M_PI/2-thetaF);
+    FocusMirror(zF,xF,B);
+    */
+
+    // focus 4: use (z,x) coordinates for tune parameters
+    double zF = zS + focusTuneZ;
+    double xF = xS + focusTuneX;
+    FocusMirror(zF,xF,B);
+
+    // re-define mirror attributes to be w.r.t vessel front plane
+    double mirrorCenterZ = zM - vesselZmin;
+    double mirrorCenterX = xM;
+    double mirrorRadius = rM;
+
+    // spherical mirror patch cuts and rotation
+    double mirrorThetaRot = std::asin(mirrorCenterX/mirrorRadius);
+    double mirrorTheta1 = mirrorThetaRot - std::asin((mirrorCenterX-mirrorRmin)/mirrorRadius);
+    double mirrorTheta2 = mirrorThetaRot + std::asin((mirrorRmax-mirrorCenterX)/mirrorRadius);
+
+    // if debugging, draw full sphere
+    if(debug_mirror) { mirrorTheta1=0; mirrorTheta2=M_PI; /*mirrorPhiw=2*M_PI;*/ };
+
+    // solid : create sphere at origin, with specified angular limits;
+    // phi limits are increased to fill gaps (overlaps are cut away later)
+    Sphere mirrorSolid1(
+        mirrorRadius,
+        mirrorRadius + mirrorThickness,
+        mirrorTheta1,
+        mirrorTheta2,
+        -40*degree,
+        40*degree
+        );
+
+    /* CAUTION: if any of the relative placements or boolean operations below
+     * are changed, you MUST make sure this does not break access to the sphere
+     * primitive and positioning in Juggler `IRTAlgorithm`; cross check the
+     * mirror sphere attributes carefully!
+     */
+    /*
+    // PRINT MIRROR ATTRIBUTES (before any sector z-rotation)
+    printf("zM = %f\n",zM); // sphere centerZ, w.r.t. IP
+    printf("xM = %f\n",xM); // sphere centerX, w.r.t. IP
+    printf("rM = %f\n",rM); // sphere radius
+    */
+
+    // mirror placement transformation (note: transformations are in reverse order)
+    auto mirrorPos = Position(mirrorCenterX, 0., mirrorCenterZ) + originFront;
+    Transform3D mirrorPlacement(
+          Translation3D(mirrorPos.x(), mirrorPos.y(), mirrorPos.z()) // re-center to specified position
+        * RotationY(-mirrorThetaRot) // rotate about vertical axis, to be within vessel radial walls
+        );
+
+    // cut overlaps with other sectors using "pie slice" wedges, to the extent specified
+    // by `mirrorPhiw`
+    Tube pieSlice( 0.01*cm, vesselRmax2, tankLength/2.0, -mirrorPhiw/2.0, mirrorPhiw/2.0);
+    IntersectionSolid mirrorSolid2( pieSlice, mirrorSolid1, mirrorPlacement );
+
+    // mirror volume, attributes, and placement
+    Volume mirrorVol(detName+"_mirror_"+secName, mirrorSolid2, mirrorMat);
+    mirrorVol.setVisAttributes(mirrorVis);
+    auto mirrorPV2 = gasvolVol.placeVolume(mirrorVol,
+          RotationZ(sectorRotation) // rotate about beam axis to sector
+        * Translation3D(0,0,0)
+        );
+
+    // properties
+    DetElement mirrorDE(det, Form("mirror_de%d", isec), isec);
+    mirrorDE.setPlacement(mirrorPV2);
+    SkinSurface mirrorSkin(desc, mirrorDE, Form("mirror_optical_surface%d", isec), mirrorSurf, mirrorVol);
+    mirrorSkin.isValid();
+
+
   }; // END SECTOR LOOP //////////////////////////
 
 
@@ -387,7 +537,7 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
   // place mother volume (vessel)
   Volume motherVol = desc.pickMotherVolume(det);
   PlacedVolume vesselPV = motherVol.placeVolume(vesselVol,
-      Position(0, 0, vesselZ0) - snoutFront
+      Position(0, 0, vesselZmin) - originFront
       );
   vesselPV.addPhysVolID("system", detID);
   det.setPlacement(vesselPV);

src/FieldMapBrBz.cpp

 #include <DD4hep/DetFactoryHelper.h>
 #include <DD4hep/FieldTypes.h>
+#include <DD4hep/Printout.h>
 #include <XML/Utilities.h>
 
 #include <cstdlib>
@@ -12,6 +13,8 @@
 #include <tuple>
 namespace fs = std::filesystem;
 
+#include "FileLoaderHelper.h"
+
 using namespace dd4hep;
 
 
@@ -19,7 +22,7 @@ using namespace dd4hep;
 class FieldMapBrBz : public dd4hep::CartesianField::Object
 {
 public:
-    FieldMapBrBz(const std::string &field_type);
+    FieldMapBrBz(const std::string &field_type = "magnetic");
     void Configure(double rmin, double rmax, double rstep, double zmin, double zmax, double zstep);
     void LoadMap(const std::string &map_file, double scale);
     void GetIndices(double r, double z, int &ir, int &iz, double &dr, double &dz);
@@ -35,7 +38,7 @@ private:
 };
 
 // constructor
-FieldMapBrBz::FieldMapBrBz(const std::string &field_type = "magnetic")
+FieldMapBrBz::FieldMapBrBz(const std::string &field_type)
 {
     std::string ftype = field_type;
     for (auto &c : ftype) { c = tolower(c); }
@@ -185,20 +188,16 @@ static Ref_t create_field_map_brbz(Detector & /*lcdd*/, xml::Handle_t handle)
     std::string field_map_file = x_par.attr<std::string>(_Unicode(field_map));
     std::string field_map_url = x_par.attr<std::string>(_Unicode(url));
 
+    EnsureFileFromURLExists(field_map_url,field_map_file);
+
     double field_map_scale = x_par.attr<double>(_Unicode(scale));
 
     if( !fs::exists(fs::path(field_map_file))  ) {
-      auto ret = std::system(("mkdir -p fieldmaps && "
-                             "wget " +
-                             field_map_url + " -O " + field_map_file).c_str());
-
-      if (!fs::exists(fs::path(field_map_file))) {
-        std::cerr << "ERROR: file, " << field_map_file << ", does not exist\n";
-        std::quick_exit(1);
-      }
+        printout(ERROR, "FieldMapBrBz", "file " + field_map_file + " does not exist");
+        printout(ERROR, "FieldMapBrBz", "use a FileLoader plugin before the field element");
+        std::_Exit(EXIT_FAILURE);
     }
 
-
     auto map = new FieldMapBrBz(field_type);
     map->Configure(r_dim.rmin(), r_dim.rmax(), r_dim.step(), z_dim.zmin(), z_dim.zmax(), z_dim.step());
 

src/FileLoader.cpp

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

src/FileLoaderHelper.h

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

src/GaseousRICH_geo.cpp

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

src/HybridCalorimeter_geo.cpp

@@ -35,54 +35,69 @@ static Ref_t create_detector(Detector& desc, xml::Handle_t handle, SensitiveDete
     DetElement det(detName, detID);
     sens.setType("calorimeter");
 
-    auto glass_material = desc.material("PbGlass");
+    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 = desc.constantAsDouble("EcalEndcapN_rmin");
+    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 glass_shift_z = desc.constantAsDouble("GlassModule_z0");
-    double Thickness = desc.constantAsDouble("EcalEndcapN_thickness");
+    double thickness = desc.constantAsDouble("EcalEndcapN_thickness");
     double trans_radius = desc.constantAsDouble("EcalEndcapNCrystal_rmax");
-    double Glass_ShiftZ = desc.constantAsDouble("GlassModule_z0");
+    double Glass_z0 = desc.constantAsDouble("GlassModule_z0");
     double Glass_Width = desc.constantAsDouble("GlassModule_width");
-    double Glass_Thickness = desc.constantAsDouble("GlassModule_length");
+    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_thickness = desc.constantAsDouble("CrystalModule_length");
     double Crystal_Gap = desc.constantAsDouble("CrystalModule_wrap");
     double crystal_distance = desc.constantAsDouble("CrystalModule_distance");
-    double Crystal_shift_z = desc.constantAsDouble("CrystalModule_z0");
+    double Crystal_z0 = desc.constantAsDouble("CrystalModule_z0");
 
     // RIn and ROut will define outer tube embedding the calorimeter
-    // centers_rin/out define the maximum radius of module centers
+    // centers_rmin/out define the maximum radius of module centers
     // so that modules are not overlapping with mother tube volume
-    double hypotenuse = sqrt(0.5 * glass_distance * glass_distance);
-    double centers_rin = RIn + hypotenuse + 1*mm;
-    double centers_rout = ROut - hypotenuse - 1*mm;
+    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
-
-    // Assembly assembly(detName);
-
-    Tube outer_tube(RIn, ROut, SizeZ / 2.0, 0., 360.0 * deg);
-    Volume ecal_vol("negative_ecal", outer_tube, air_material);
+    // 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;
-    glass_shift_z = Thickness / 2. - Glass_Thickness / 2.;
 
     // Geometry of modules
-    Box glass_box("glass_box", Glass_Width * 0.5, Glass_Width * 0.5, Glass_Thickness * 0.5);
+    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);
+    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);
@@ -90,9 +105,11 @@ static Ref_t create_detector(Detector& desc, xml::Handle_t handle, SensitiveDete
     // GLASS
     double diameter = 2 * Glass_OuterR;
 
-    // How many towers do we have per row/columnt?
-    // Add a gap + diameter as if we have N towers, we have N-1 gaps;
-    int towersInRow = std::ceil((diameter + Glass_Gap) /  (Glass_Width + Glass_Gap));
+    // 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;
@@ -113,112 +130,45 @@ static Ref_t create_detector(Detector& desc, xml::Handle_t handle, SensitiveDete
 
     int moduleIndex = 0;
 
-//    fmt::print("\nCE EMCAL GLASS SQUARE START\n");
-//    fmt::print("Glass_Thickness = {} cm;\n", Glass_Thickness / cm);
-//    fmt::print("Glass_Width     = {} cm;\n", Glass_Width / cm);
-//    fmt::print("Glass_Gap       = {} cm;\n", Glass_Gap / cm);
-//    fmt::print("Glass_InnerR    = {} cm;\n", Glass_InnerR / cm);
-//    fmt::print("Glass_OuterR    = {} cm;\n", Glass_OuterR / cm);
-//    fmt::print("Glass_PosZ      = {} cm;\n", glass_shift_z / cm);
-//    fmt::print("Towers in Row/Col   = {} cm;\n", glass_shift_z / cm);
-//    fmt::print("Top left tower pos  = {:<10} {:<10} cm;\n", -leftTowerPos / cm, topTowerPos / cm);
-// fmt::print("#Towers info:\n");
-// fmt::print("#{:<5} {:<6} {:<3} {:<3} {:>10} {:>10}   {}\n", "idx",  "code", "col", "row", "x", "y", "name");
-
-
-    // We first do a "dry run", not really placing modules,
-    // but figuring out the ID scheme, number of modules, etc.
-    int glassModuleCount = 0;
-    int crystalModuleCount = 0;
-    int firstCrystRow = 1000000;   // The first row, where crystals are started
-    int firstCrystCol = 1000000;   // The fist column, where crystals are started
-    for(int rowIndex=0; rowIndex < towersInRow; rowIndex++) {
-      for(int colIndex=0; colIndex < towersInRow; colIndex++) {
-        double glass_x = leftTowerPos + colIndex * glass_distance;
-        double glass_y = topTowerPos + rowIndex * glass_distance;
-        double r = sqrt(glass_x * glass_x + glass_y * glass_y);
-
-        if (r < centers_rout && r > centers_rin) {
-          // we are placing something
-          if(r<trans_radius) {
-            // 4 crystal modules will be placed
-            crystalModuleCount+=4;
-
-            // Finding the first col and row where crystals start
-            // is the same algorithm as finding a minimum in array
-            if(colIndex<firstCrystCol) {
-              firstCrystCol = colIndex;
-            }
-            if(rowIndex<firstCrystRow) {
-              firstCrystRow = rowIndex;
-            }
-          }
-          else
-          {
-            // glass module will be places
-            glassModuleCount++;
-          }
-        }
-      }
-    }
-    // fmt::print("#Towers info:\n");
-    // fmt::print("#{:<5} {:<6} {:<3} {:<3} {:>10} {:>10}   {}\n", "idx",  "code", "col", "row", "x", "y", "name");
     int glass_module_index = 0;
     int cryst_module_index = 0;
     for(int rowIndex=0; rowIndex < towersInRow; rowIndex++) {
       for(int colIndex=0; colIndex < towersInRow; colIndex++) {
-        double glass_x = leftTowerPos + colIndex * (Glass_Width + Glass_Gap);
-        double glass_y = topTowerPos + rowIndex * (Glass_Width + Glass_Gap);
-        double r = sqrt(glass_x * glass_x + glass_y * glass_y);
-
-        if (r < centers_rout && r > centers_rin) {
-          int code = 1000 * rowIndex + colIndex;
-          std::string name = fmt::format("ce_EMCAL_glass_phys_{}", code);
-
-          if(r<trans_radius) {
-
-            // first crystal module
-            double crystal_x = glass_x - crystal_distance / 2;
-            double crystal_y = glass_y - crystal_distance / 2;
-            auto   placement = ecal_vol.placeVolume(crystal_module, Position(crystal_x, crystal_y, Crystal_shift_z));
-            placement.addPhysVolID("sector", 1);
-            placement.addPhysVolID("module", cryst_module_index++);
-
-            // second crystal module
-            crystal_x = glass_x + crystal_distance / 2;
-            crystal_y = glass_y - crystal_distance / 2;
-            placement = ecal_vol.placeVolume(crystal_module, Position(crystal_x, crystal_y, Crystal_shift_z));
-            placement.addPhysVolID("sector", 1);
-            placement.addPhysVolID("module", cryst_module_index++);
-
-            // third crystal module
-            crystal_x = glass_x - crystal_distance / 2;
-            crystal_y = glass_y + crystal_distance / 2;
-            placement = ecal_vol.placeVolume(crystal_module, Position(crystal_x, crystal_y, Crystal_shift_z));
-            placement.addPhysVolID("sector", 1);
-            placement.addPhysVolID("module", cryst_module_index++);
-
-            // forth crystal module
-            crystal_x = glass_x + crystal_distance / 2;
-            crystal_y = glass_y + crystal_distance / 2;
-            placement = ecal_vol.placeVolume(crystal_module, Position(crystal_x, crystal_y, Crystal_shift_z));
-            placement.addPhysVolID("sector", 1);
-            placement.addPhysVolID("module", cryst_module_index++);
-          }
-          else
-          {
-            // glass module
-            auto placement = ecal_vol.placeVolume(glass_module, Position(glass_x, glass_y, glass_shift_z));
-            placement.addPhysVolID("sector", 2);
-            placement.addPhysVolID("module", glass_module_index++);
+        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++);
+              }
+            }
           }
-
-
-          // fmt::print(" {:<5} {:<6} {:<3} {:<3} {:>10.4f} {:>10.4f}   {}\n", towerIndex, code, colIndex, rowIndex, x / cm, y / cm, name);
-          glass_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");
 

src/MRich_geo.cpp

@@ -29,21 +29,9 @@ static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDet
   string         det_name = x_det.nameStr();
   DetElement     sdet(det_name, x_det.id());
   Assembly       assembly(det_name);
-  sens.setType("photoncounter");
+  sens.setType("tracker");
   OpticalSurfaceManager surfMgr = description.surfaceManager();
 
-  // read module positions
-  std::vector<std::pair<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_pair(x_positions.scale() * x_position.x() * mm,
-                       x_positions.scale() * x_position.y() * mm));
-    }
-  }
-
   bool projective = getAttrOrDefault(x_det, _Unicode(projective), false);
   bool reflect    = x_det.reflect(true);
 
@@ -56,6 +44,7 @@ static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDet
 
   int                     n_sensor = 1;
 
+  // dimensions
   xml::Component dims   = x_det.dimensions();
   auto           rmin   = dims.rmin();
   auto           rmax   = dims.rmax();
@@ -63,6 +52,20 @@ static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDet
   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();
@@ -70,193 +73,231 @@ static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDet
   double     mod_height          = getAttrOrDefault(x_mod, _U(height), 130.0 * mm);
   double     mod_length          = getAttrOrDefault(x_mod, _U(length), 130.0 * mm);
 
-  // various components
-  xml_comp_t x_frame    = x_mod.child(_Unicode(frame));
-  xml_comp_t x_aerogel  = x_mod.child(_Unicode(aerogel));
-  xml_comp_t x_lens     = x_mod.child(_Unicode(lens));
-  xml_comp_t x_mirror   = x_mod.child(_Unicode(mirror));
-  xml_comp_t x_photodet = x_mod.child(_Unicode(photodet));
-
   // 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
-  double     frame_thickness  = getAttrOrDefault(x_frame, _U(thickness), 2.0 * mm);
+  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
-  xml_comp_t x_aerogel_frame = x_aerogel.child(_Unicode(frame));
-  double     aerogel_width   = getAttrOrDefault(x_aerogel, _U(width), 130.0 * mm);
-  double     aerogel_length  = getAttrOrDefault(x_aerogel, _U(length), 130.0 * mm);
-  double     foam_thickness  = getAttrOrDefault(x_aerogel_frame, _U(thickness), 2.0 * mm);
-  Material   foam_mat        = description.material(x_aerogel_frame.materialStr());
-  Material   aerogel_mat     = description.material(x_aerogel.materialStr());
-  auto       aerogel_vis     = getAttrOrDefault<std::string>(x_aerogel, _U(vis), std::string("InvisibleWithDaughters"));
-  auto       foam_vis        = getAttrOrDefault<std::string>(x_aerogel_frame, _U(vis), std::string("RedVis"));
-
-  // aerogel foam frame
-  Box foam_box(aerogel_width / 2.0 + foam_thickness, aerogel_width / 2.0 + foam_thickness, (aerogel_length + foam_thickness) / 2.0);
-  Box aerogel_sub_box(aerogel_width / 2.0, aerogel_width / 2.0, (aerogel_length + foam_thickness) / 2.0);
-  SubtractionSolid foam_frame_solid(foam_box, aerogel_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));
-  double foam_frame_zpos = -mod_length / 2.0 + frame_thickness + (aerogel_length + foam_thickness) / 2.0;
-  m_volume.placeVolume(foam_vol,Position(0,0,foam_frame_zpos));
-
-  // 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));
-  double aerogel_zpos = foam_frame_zpos + foam_thickness / 2.0;
-  pv = m_volume.placeVolume(aerogel_vol,Position(0,0,aerogel_zpos));
-  DetElement aerogel_de(mod_de, mod_name + std::string("_aerogel_de") + std::to_string(1), 1);
-  aerogel_de.setPlacement(pv);
-
-  auto aerogel_surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault(x_aerogel, _Unicode(surface), "MRICH_AerogelOpticalSurface"));
-  SkinSurface skin0(description, aerogel_de, Form("MRICH_aerogel_skin_surface_%d", 1), aerogel_surf, aerogel_vol);
-  skin0.isValid();
-
-  // Fresnel 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);//,Position(0,0,lens_zpos));
-
-  Solid  fresnel_lens_solid;
-
-  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);
-    //Volume   groove_vol(groove_solid, lens_mat, par->name.c_str(), 0, 0, 0);
-    //new G4PVPlacement(0, par->pos, Groove_log[i], par->name.c_str(), motherLV, false, 0, OverlapCheck());
-    //phi1 = phi1 + halfpi;  //g4 pre-defined: halfpi=pi/2
-    //Tube       sub_cylinder(r0, r1, 3*eff_diameter);
-    //IntersectionSolid groove_solid(lens_box,lens_sphere, Position(0,0,-eff_diameter/2.0 + lens_thickness/2.0+(t-lens_t)/2.0 ));
-    //IntersectionSolid lens_ring(groove_solid, sub_cylinder);
-    //if (i_groove == 0) {
-    //  fresnel_lens_solid = groove_solid;
-    //} else {
-    //  fresnel_lens_solid = UnionSolid(fresnel_lens_solid, groove_solid);
-    //}
-    //r0 = r1;
-    //if(i_groove > 3) { 
-    //  SubtractionSolid flat_lens(lens_box,Tube(0.0, r0, 3*eff_diameter));
-    //  fresnel_lens_solid = UnionSolid(fresnel_lens_solid, flat_lens);
-    //  break; // temporary
-    //}
-    i_groove++;
-    groove_rmin = (i_groove  )*groove_pitch;
-    groove_rmax = (i_groove+1)*groove_pitch;
+  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_solid = UnionSolid(fresnel_lens_solid, flat_lens);
-  //Volume   lens_vol(mod_name + "_lens", fresnel_lens_solid, lens_mat);
+  // 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));
-  double lens_zpos = aerogel_zpos +aerogel_length/ 2.0 + foam_thickness + lens_thickness/2.0;
-  pv = m_volume.placeVolume(lens_vol,Position(0,0,lens_zpos));
-  DetElement lens_de(mod_de, mod_name + std::string("_lens_de") + std::to_string(1), 1);
-  lens_de.setPlacement(pv);
+    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();
+  }
 
-  auto surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault(x_lens, _Unicode(surface), "MRICH_LensOpticalSurface"));
-  SkinSurface skin(description, lens_de, Form("MRichFresnelLens_skin_surface_%d", 1), surf, lens_vol);
-  skin.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
-  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);
-  double mirror_zpos = lens_zpos + lens_thickness/2.0 + foam_thickness + mirror_length/2.0;
-  pv = m_volume.placeVolume(mirror_vol,Position(0,0,mirror_zpos));
-  DetElement mirror_de(mod_de, mod_name + std::string("_mirror_de") + std::to_string(1), 1);
-  mirror_de.setPlacement(pv);
-
-  auto mirror_surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault(x_mirror, _Unicode(surface), "MRICH_MirrorOpticalSurface"));
-  SkinSurface skin1(description, mirror_de, Form("MRICH_mirror_skin_surface_%d", 1), mirror_surf, mirror_vol);
-  skin1.isValid();
+  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
-  xml_comp_t x_photodet_sensor  = x_photodet.child(_Unicode(sensor));
-  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);
-  double     sensor_thickness   = getAttrOrDefault(x_photodet_sensor, _U(thickness), 2.0 * mm);
-  Material   photodet_mat       = description.material(x_photodet.materialStr());
-  Material   sensor_mat         = description.material(x_photodet_sensor.materialStr());
-  int     sensor_nx   = getAttrOrDefault(x_photodet_sensor, _Unicode(nx), 2);
-  int     sensor_ny   = getAttrOrDefault(x_photodet_sensor, _Unicode(ny), 2);
-
-  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);
-  double window_zpos = mirror_zpos + mirror_length/2.0+photodet_thickness/2.0;
-  pv = m_volume.placeVolume(window_vol,Position(0,0,window_zpos));
-  DetElement   comp_de(mod_de, std::string("mod_sensor_de_") + std::to_string(1) ,  1);
-  comp_de.setPlacement(pv);
-  // sensitive
-  pv.addPhysVolID("sensor", n_sensor);
-  window_vol.setSensitiveDetector(sens);
-  sensitives[mod_name].push_back(pv);
-  ++n_sensor;
-
-  // photon detector electronics layers
-  double layer_zpos = window_zpos + photodet_thickness/2.0;
-  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);
-    layer_zpos += layer_thickness / 2.0;
-    pv = m_volume.placeVolume(layer_vol,Position(0,0,layer_zpos));
-    DetElement layer_de(mod_de, std::string("mod_layer_de_") + std::to_string(i_layer),  1);
-    layer_de.setPlacement(pv);
-    layer_zpos += layer_thickness / 2.0;
-    i_layer++;
+  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) {
@@ -277,11 +318,6 @@ static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDet
   module_assembly_delements[mod_name] = mod_de;
   // end module
 
-  // detector envelope
-  Tube   envShape(rmin, rmax, length / 2., 0., 2 * M_PI);
-  Volume envVol("MRICH_Envelope", envShape, air);
-  envVol.setVisAttributes(description.visAttributes(x_det.visStr()));
-
   // place modules in the sectors (disk)
   auto points = athena::geo::fillSquares({0., 0.}, mod_width, rmin, rmax);
 
@@ -291,45 +327,51 @@ static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDet
   // 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 = p.first;
-    double y = p.second;
-    double z = -zpos;
+    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}, {y,-x}, {-x,-y}, {-y,x}}) {
+    for (auto& p: decltype(positions){{x,y,z0}, {y,-x,z0}, {-x,-y,z0}, {-y,x,z0}}) {
 
       // get positions (intentional shadowing)
-      double x = p.first;
-      double y = p.second;
-
-      // get angles
-      double rotAngX = atan(y/z);
-      double rotAngY = -1.*atan(x/z);
-
-      /*
-      ROOT::Math::XYZVector x_location(p.x(), p.y(), zmin+std::signbit(zmin)*mod_length/2.0);
-      ROOT::Math::XYZVector z_dir(0, 0, 1);
-      ROOT::Math::XYZVector x_dir(1, 0, 0);
-      ROOT::Math::XYZVector rot_axis = x_location.Cross(z_dir);
-      double rot_angle = ROOT::Math::VectorUtil::Angle(z_dir,x_location);
-      ROOT::Math::AxisAngle proj_rot(rot_axis,-1.0*rot_angle);
-      ROOT::Math::AxisAngle grid_fix_rot(x_location,0.0*rot_angle);
-      auto new_x_dir = grid_fix_rot*x_dir;
-      // operations are inversely ordered
-      Transform3D tr = Translation3D(p.x(), p.y(), 0.) // move to position
-                      * RotationY(roty);              // facing z = 0.
-      */
+      double x = std::get<0>(p);
+      double y = std::get<1>(p);
+      double z0 = std::get<2>(p);
 
       Transform3D tr;
       if(projective) {
-        tr = Translation3D(x, y, 0)
-           * RotationX(rotAngX)
-           * RotationY(rotAngY);
+        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);
       }
@@ -346,6 +388,16 @@ static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDet
     }
   }
 
+  // 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) {
@@ -385,7 +437,7 @@ static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDet
 //    modVol.placeVolume(mVol, Position(0., 0., -0.1*mm));
 //
 //    modVol.setVisAttributes(description.visAttributes(mods.visStr()));
-//    sens.setType("photoncounter");
+//    sens.setType("tracker");
 //    modVol.setSensitiveDetector(sens);
 //
 //    // place modules in the sectors (disk)