scripts/plot_ce_ecal_placement.py

+'''
+    A script to calcualte placement of ecal endcap modules
+    lxml is not included in container, get it by simply typing 'pip install lxml'
+    Author: Chao Peng (ANL)
+    Date: 06/17/2021
+'''
+
+import os
+import numpy as np
+import argparse
+import DDG4
+from lxml import etree as ET
+from matplotlib import pyplot as plt
+from matplotlib.collections import PatchCollection
+from matplotlib.patches import Rectangle, Circle
+
+
+CRYSTAL_SIZE = (20., 20., 200.) # mm
+CRYSTAL_GAP = 0.5 # mm
+CRYSTAL_ALIGNMENT = [
+    (5, 21), (5, 21), (5, 21), (4, 22),
+    (3, 23), (0, 26), (0, 24), (0, 24),
+    (0, 24), (0, 24), (0, 24), (0, 24),
+    (0, 22), (0, 22), (0, 20), (0, 20),
+    (0, 18), (0, 18), (0, 16), (0, 16),
+    (0, 14), (0, 14), (0, 12), (0, 12),
+    (0,  6), (0,  6),
+]
+
+GLASS_SIZE = (40., 40., 400.) # mm
+GLASS_GAP = 1.0 # mm
+GLASS_ALIGNMENT = [
+    (13, 10), (13, 10), (13, 10), (12, 10),
+    (12, 10), (12, 10), (11, 11), (10, 11),
+    (9, 12),  (8, 13),  (7, 13),  (6, 14),
+    (3, 16),  (0, 18),  (0, 18),  (0, 16),
+    (0, 16),  (0, 14),  (0, 13),  (0, 11),
+    (0, 10),  (0, 7),   (0, 3),
+]
+
+# calculate positions of modules with a quad-alignment and module size
+def individual_placement(alignment, module_x, module_y, gap=0.):
+    placements = []
+    for row, (start, num) in enumerate(alignment):
+        for col in np.arange(start, start + num):
+            placements.append(((col + 0.5)*(module_y + gap), (row + 0.5)*(module_x + gap)))
+    placements = np.asarray(placements)
+    return np.vstack((placements,
+            np.vstack((placements.T[0]*-1., placements.T[1])).T,
+            np.vstack((placements.T[0], placements.T[1]*-1.)).T,
+            np.vstack((placements.T[0]*-1., placements.T[1]*-1.)).T))
+
+
+def draw_placement(axis, colors=('teal'), module_alignment=((CRYSTAL_SIZE, CRYSTAL_GAP, CRYSTAL_ALIGNMENT))):
+    xmin, ymin, xmax, ymax = 0., 0., 0., 0.
+    patches = []
+    numbers = []
+    for color, (mod_size, mod_gap, alignment) in zip(colors, module_alignment):
+        placements = individual_placement(alignment, *mod_size[:2], mod_gap)
+        boxes = [Rectangle((x - (mod_size[0] + mod_gap)/2., y - (mod_size[1] + mod_gap)/2.), mod_size[0], mod_size[1])
+                 for x, y in placements]
+        patches.append(Rectangle((0., 0.), *mod_size[:2], facecolor=color, alpha=0.5, edgecolor='k'))
+        numbers.append(len(placements))
+        pc = PatchCollection(boxes, facecolor=color, alpha=0.5, edgecolor='k')
+
+        xmin = min(xmin, placements.T[0].min() - 8.*(mod_size[0] + mod_gap))
+        ymin = min(ymin, placements.T[1].min() - 8.*(mod_size[1] + mod_gap))
+        xmax = max(xmax, placements.T[0].max() + 8.*(mod_size[0] + mod_gap))
+        ymax = max(ymax, placements.T[1].max() + 8.*(mod_size[1] + mod_gap))
+
+        # Add collection to axes
+        axis.add_collection(pc)
+    axis.set_xlim(xmin, xmax)
+    axis.set_ylim(ymin, ymax)
+    return axis, patches, numbers
+
+
+def compact_constants(path, names):
+    if not os.path.exists(path):
+        print('Cannot find compact file \"{}\".'.format(path))
+        return []
+    kernel = DDG4.Kernel()
+    description = kernel.detectorDescription()
+    kernel.loadGeometry("file:{}".format(path))
+    try:
+        vals = [description.constantAsDouble(n) for n in names]
+    except:
+        print('Fail to extract values from {}, return empty.'.format(names))
+        vals = []
+    kernel.terminate()
+    return vals
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser()
+    parser.add_argument('-s', '--save', default='ce_ecal_placement_test.xml',
+            help='path to save compact file.')
+    parser.add_argument('-c', '--compact', default='',
+            help='compact file to get contant to plot')
+    parser.add_argument('--radii-constants', dest='radii', default='EcalBarrel_rmin',
+            help='constant names in compact file to plot, seprate by \",\"')
+    parser.add_argument('--individual', dest='indiv', action='store_true',
+            help='individual block placements instead of line placements')
+    args = parser.parse_args()
+
+    data = ET.Element('lccdd')
+    defines = ET.SubElement(data, 'define')
+
+    # constants: name, value
+    CONSTANTS = [
+        ('CrystalModule_sx', '{:.2f}*mm'.format(CRYSTAL_SIZE[0])),
+        ('CrystalModule_sy', '{:.2f}*mm'.format(CRYSTAL_SIZE[1])),
+        ('CrystalModule_sz', '{:.2f}*mm'.format(CRYSTAL_SIZE[2])),
+        ('CrystalModule_wrap', '{:.2f}*mm'.format(CRYSTAL_GAP)),
+
+        ('GlassModule_sx', '{:.2f}*mm'.format(GLASS_SIZE[0])),
+        ('GlassModule_sy', '{:.2f}*mm'.format(GLASS_SIZE[1])),
+        ('GlassModule_sz', '{:.2f}*mm'.format(GLASS_SIZE[2])),
+        ('GlassModule_wrap', '{:.2f}*mm'.format(GLASS_GAP)),
+
+        ('CrystalModule_z0', '10.*cm'),
+        ('GlassModule_z0', '0.0*cm'),
+        ('EcalEndcapN_z0', '-EcalEndcapN_zmin-max(CrystalModule_sz,GlassModule_sz)/2.'),
+
+        ('CrystalModule_dx', 'CrystalModule_sx + CrystalModule_wrap'),
+        ('CrystalModule_dy', 'CrystalModule_sy + CrystalModule_wrap'),
+        ('GlassModule_dx', 'GlassModule_sx + GlassModule_wrap'),
+        ('GlassModule_dy', 'GlassModule_sy + GlassModule_wrap'),
+    ]
+
+# line-by-line alignment start pos, total number of blocks
+    for name, value in CONSTANTS:
+        constant = ET.SubElement(defines, 'constant')
+        constant.set('name', name)
+        constant.set('value', value)
+
+    # this value will be used multiple times, so define it here
+    readout_name = 'EcalEndcapNHits'
+
+    # detector and its dimension/position/rotation
+    dets = ET.SubElement(data, 'detectors')
+    cmt = ET.SubElement(dets, 'comment')
+    cmt.text = ' Backwards Endcap EM Calorimeter, placements generated by script '
+
+    det = ET.SubElement(dets, 'detector')
+    det.set('id', 'ECalEndcapN_ID')
+    det.set('name', 'EcalEndcapN')
+    det.set('type', 'HomogeneousCalorimeter')
+    det.set('readout', readout_name)
+
+    pos = ET.SubElement(det, 'position')
+    pos.set('x', '0')
+    pos.set('y', '0')
+    pos.set('z', 'EcalEndcapN_z0')
+
+    rot = ET.SubElement(det, 'rotation')
+    rot.set('x', '0')
+    rot.set('y', '0')
+    rot.set('z', '0')
+
+    # placements of modules
+    plm = ET.SubElement(det, 'placements')
+    pltype = 'individuals' if args.indiv else 'lines'
+
+    # crystal
+    crystal = ET.SubElement(plm, pltype)
+    crystal.set('sector', '1')
+    crystal_mod = ET.SubElement(crystal, 'module')
+    crystal_mod.set('sizex', 'CrystalModule_sx')
+    crystal_mod.set('sizey', 'CrystalModule_sy')
+    crystal_mod.set('sizez', 'CrystalModule_sz')
+    crystal_mod.set('material', 'PbWO4')
+    crystal_mod.set('vis', 'AnlTeal')
+    crystal_wrap = ET.SubElement(crystal, 'wrapper')
+    crystal_wrap.set('thickness', 'CrystalModule_wrap')
+    crystal_wrap.set('material', 'Epoxy')
+    crystal_wrap.set('vis', 'WhiteVis')
+    # crystal placements (for individuals)
+    if args.indiv:
+        for m, (x, y) in enumerate(individual_placement(CRYSTAL_ALIGNMENT, *CRYSTAL_SIZE[:2], CRYSTAL_GAP)):
+            module = ET.SubElement(crystal, 'placement')
+            module.set('x', '{:.3f}*mm'.format(x))
+            module.set('y', '{:.3f}*mm'.format(y))
+            module.set('z', 'CrystalModule_z0')
+            module.set('id', '{:d}'.format(m))
+    # crystal placements (for lines)
+    else:
+        crystal.set('mirrorx', 'true')
+        crystal.set('mirrory', 'true')
+        for row, (begin, nmods) in enumerate(CRYSTAL_ALIGNMENT):
+            line = ET.SubElement(crystal, 'line')
+            line.set('axis', 'x')
+            line.set('x', 'CrystalModule_dx/2.')
+            line.set('y', 'CrystalModule_dy*{:d}/2.'.format(row*2 + 1))
+            line.set('z', 'CrystalModule_z0')
+            line.set('begin', '{:d}'.format(begin))
+            line.set('nmods', '{:d}'.format(nmods))
+
+
+    # glass
+    glass = ET.SubElement(plm, pltype)
+    glass.set('sector', '2')
+    glass_mod = ET.SubElement(glass, 'module')
+    glass_mod.set('sizex', 'GlassModule_sx')
+    glass_mod.set('sizey', 'GlassModule_sy')
+    glass_mod.set('sizez', 'GlassModule_sz')
+    # TODO: change glass material
+    glass_mod.set('material', 'PbGlass')
+    glass_mod.set('vis', 'AnlBlue')
+    glass_wrap = ET.SubElement(glass, 'wrapper')
+    glass_wrap.set('thickness', 'GlassModule_wrap')
+    glass_wrap.set('material', 'Epoxy')
+    glass_wrap.set('vis', 'WhiteVis')
+    # crystal placements (for individuals)
+    if args.indiv:
+        for m, (x, y) in enumerate(individual_placement(GLASS_ALIGNMENT, *GLASS_SIZE[:2], GLASS_GAP)):
+            module = ET.SubElement(glass, 'placement')
+            module.set('x', '{:.3f}*mm'.format(x))
+            module.set('y', '{:.3f}*mm'.format(y))
+            module.set('z', 'GlassModule_z0')
+            module.set('id', '{:d}'.format(m))
+    # crystal placements (for lines)
+    else:
+        glass.set('mirrorx', 'true')
+        glass.set('mirrory', 'true')
+        for row, (begin, nmods) in enumerate(GLASS_ALIGNMENT):
+            line = ET.SubElement(glass, 'line')
+            line.set('axis', 'x')
+            line.set('x', 'GlassModule_dx/2.')
+            line.set('y', 'GlassModule_dy*{:d}/2.'.format(row*2 + 1))
+            line.set('z', 'GlassModule_z0')
+            line.set('begin', '{:d}'.format(begin))
+            line.set('nmods', '{:d}'.format(nmods))
+
+
+    # readout
+    readouts = ET.SubElement(data, 'readouts')
+    cmt = ET.SubElement(readouts, 'comment')
+    cmt.text = 'Effectively no segmentation, the segmentation is used to provide cell dimension info'
+    readout = ET.SubElement(readouts, 'readout')
+    readout.set('name', readout_name)
+    seg = ET.SubElement(readout, 'segmentation')
+    # need segmentation to provide cell dimension info
+    # seg.set('type', 'NoSegmentation')
+    seg.set('type', 'MultiSegmentation')
+    seg.set('key', 'sector')
+    crystal_seg = ET.SubElement(seg, 'segmentation')
+    crystal_seg.set('name', 'CrystalSeg')
+    crystal_seg.set('key_value', '1')
+    crystal_seg.set('type', 'CartesianGridXY')
+    crystal_seg.set('grid_size_x', 'CrystalModule_dx')
+    crystal_seg.set('grid_size_y', 'CrystalModule_dy')
+    glass_seg = ET.SubElement(seg, 'segmentation')
+    glass_seg.set('name', 'GlassSeg')
+    glass_seg.set('key_value', '2')
+    glass_seg.set('type', 'CartesianGridXY')
+    glass_seg.set('grid_size_x', 'GlassModule_dx')
+    glass_seg.set('grid_size_y', 'GlassModule_dy')
+    rid = ET.SubElement(readout, 'id')
+    rid.text = 'system:8,sector:4,module:20,x:32:-16,y:-16'
+
+
+    text = ET.tostring(data, pretty_print=True)
+    with open(args.save, 'wb') as f:
+        f.write(text)
+
+
+    fig, ax = plt.subplots(figsize=(12, 12), dpi=160)
+    ax, patches, nblocks = draw_placement(ax, ['teal', 'royalblue'],
+            [(CRYSTAL_SIZE, CRYSTAL_GAP, CRYSTAL_ALIGNMENT), (GLASS_SIZE, GLASS_GAP, GLASS_ALIGNMENT)])
+    ax.set_xlabel('x (mm)', fontsize=24)
+    ax.set_ylabel('y (mm)', fontsize=24)
+    ax.tick_params(direction='in', labelsize=22, which='both')
+    ax.set_axisbelow(True)
+    ax.grid(linestyle=':', which='both')
+    ax.legend(patches, ['{} {}'.format(num, name) for num, name in zip(nblocks, ['PbWO$_4$', 'SciGlass'])], fontsize=24)
+
+    if args.compact and args.radii:
+        names = [c.strip() for c in args.radii.split(',')]
+        radii = compact_constants(args.compact, names)
+        for name, radius in zip(names, radii):
+            ax.add_patch(Circle((0, 0), radius*10., facecolor='none', edgecolor='k', linewidth=2))
+            ax.annotate(name, xy=(radius*10/1.4, radius*10/1.4), fontsize=22)
+    fig.savefig('ce_ecal_placement.png')
+

scripts/run_detector_simulation.py

@@ -12,7 +12,7 @@ parser = argparse.ArgumentParser(
      This program should be run in the "compact" directory.
          ''')
 parser.add_argument("-v","--verbose", help="increase output verbosity", type=int, default=0)
-parser.add_argument("--compact", help="compact detector file",default="reference_detector.xml")
+parser.add_argument("--compact", help="compact detector file",default="athena.xml")
 parser.add_argument("--vis", help="vis true/false", action="store_true",default=False)
 parser.add_argument("--ui", help="ui setting tcsh or qt; default=qt", type=str,default="qt",dest="ui")
 parser.add_argument("-b","--batch", help="batch turns off vis/ui", action="store_true",default=False, dest="batch")
@@ -86,7 +86,7 @@ def run():
 
     outputfile = args.output
     if outputfile is None:
-        outputfile = 'data/reference_detector_' + time.strftime('%Y-%m-%d_%H-%M')
+        outputfile = 'data/athena_' + time.strftime('%Y-%m-%d_%H-%M')
     #rootoutput = geant4.setupROOTOutput('RootOutput', outputfile)
     #rootoutput.HandleMCTruth = True
 

scripts/test_ACTS.cxx

+R__LOAD_LIBRARY(libDDCore.so)
+R__LOAD_LIBRARY(libActsPluginDD4hep.so)
+R__LOAD_LIBRARY(libDDG4.so)
+R__LOAD_LIBRARY(libDDG4IO.so)
+#include "DD4hep/Detector.h"
+#include "DD4hep/DetElement.h"
+#include "DD4hep/Objects.h"
+#include "DD4hep/Detector.h"
+#include "DDG4/Geant4Data.h"
+#include "DDRec/CellIDPositionConverter.h"
+#include "DDRec/SurfaceManager.h"
+#include "DDRec/Surface.h"
+
+#include "TCanvas.h"
+#include "TChain.h"
+
+#include "Acts/Geometry/TrackingGeometry.hpp"
+#include "Acts/Geometry/TrackingVolume.hpp"
+#include "Acts/Plugins/DD4hep/ConvertDD4hepDetector.hpp"
+
+/** Example loading ACTs.
+ *
+ *
+ */
+void test_ACTS(const char* compact = "athena.xml"){
+
+  using namespace ROOT::Math;
+  // -------------------------
+  // Get the DD4hep instance
+  // Load the compact XML file
+  // Initialize the position converter tool
+  dd4hep::Detector& detector = dd4hep::Detector::getInstance();
+  detector.fromCompact(compact);
+  dd4hep::rec::CellIDPositionConverter cellid_converter(detector);
+
+  //std::unique_ptr<const Acts::TrackingGeometry> 
+  auto acts_tracking_geometry  = Acts::convertDD4hepDetector (detector.world(),Acts::Logging::Level::VERBOSE);
+  //acts_tracking_geometry  = Acts::convertDD4hepDetector (detector.world(),Acts::Logging::Level::INFO);
+
+  if(acts_tracking_geometry) {
+    std::cout << "success?\n";
+  }
+  //  if(acts_tracking_geometry->highestTrackingVolume()) {
+  //    std::cout << " volume name \n ";
+  //    std::cout << acts_tracking_geometry->highestTrackingVolume()->volumeName() << std::endl;
+  //  } else {
+  //    std::cout << "derp\n";
+  //  }
+  //}
+
+}

scripts/view1/generate_eps

@@ -44,28 +44,30 @@ done
 set -- "${POSITIONAL[@]}" # restore positional parameters
 
 
-# Side  view
-dawncut 1 0 0 1 ${INPUT_FILE} ${FILE_TAG}_temp0.prim 
-dawncut -1 0 0 1 ${FILE_TAG}_temp0.prim  ${FILE_TAG}.prim
-../../bin/dawn_tweak --mag 10 
-dawn -d ${FILE_TAG}.prim 
+# Side view (lines)
+dawncut -1 0 0 1 ${INPUT_FILE} ${FILE_TAG}.prim
+# dawncut 1 0 0 1 ${INPUT_FILE} ${FILE_TAG}_temp0.prim 
+# dawncut -1 0 0 1 ${FILE_TAG}_temp0.prim  ${FILE_TAG}.prim
+../../bin/dawn_tweak --mag 10 --draw 3 --theta 165 --phi 75 --light-theta 180 --light-phi 90
+
+dawn -d ${FILE_TAG}.prim
 ps2pdf ${FILE_TAG}.eps ${FILE_TAG}_full.pdf
 gs -o ${FILE_TAG}.pdf -sDEVICE=pdfwrite \
   -c "[/CropBox [51 250 550 590] /PAGES pdfmark" \
   -f ${FILE_TAG}_full.pdf
 
-pdftoppm ${FILE_TAG}.pdf ${FILE_TAG} -png -singlefile -cropbox
-
+pdftoppm ${FILE_TAG}.pdf ${FILE_TAG} -png -singlefile -cropbox -thinlinemode solid -aaVector yes -r 600
 
 # Top view
+#dawncut 0 -1 0 1 ${INPUT_FILE} ${FILE_TAG}.prim 
 dawncut 0 1 0 1 ${INPUT_FILE} ${FILE_TAG}_temp0.prim 
 dawncut 0 -1 0 1 ${FILE_TAG}_temp0.prim  ${FILE_TAG}.prim
-../../bin/dawn_tweak --theta 270 --mag 20
+../../bin/dawn_tweak --mag 10 --draw 1 --theta 90 --phi 90
 dawn -d ${FILE_TAG}.prim 
 ps2pdf ${FILE_TAG}.eps ${FILE_TAG}_top_full.pdf
 gs -o ${FILE_TAG}_top.pdf -sDEVICE=pdfwrite \
   -c "[/CropBox [51 250 550 590] /PAGES pdfmark" \
   -f ${FILE_TAG}_top_full.pdf
 
-pdftoppm ${FILE_TAG}_top.pdf ${FILE_TAG}_top -png -singlefile -cropbox
+pdftoppm ${FILE_TAG}_top.pdf ${FILE_TAG}_top -png -singlefile -cropbox -thinlinemode solid -aaVector yes
 

scripts/view11/.DAWN_1.history

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

scripts/view12/.DAWN_1.history

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

scripts/view12/generate_eps

@@ -18,7 +18,7 @@ tan ()
 
 
 function print_the_help {
-  echo "USAGE: $0 -i <PRIM_FILE>  "
+  echo "USAGE: $0 -i <PRIM_FILE> <slices ...> "
   echo "  OPTIONS: "
   echo "            -t,--tag           filename tag (default: view1)"
   exit 
@@ -48,12 +48,16 @@ do
       shift # past argument
       shift # past value
       ;;
-    *)    # unknown option
-      #POSITIONAL+=("$1") # save it in an array for later
+    -[a-zA-Z]*) # unknown option
+      POSITIONAL+=("$1") # save it in an array for later
       echo "unknown option $1"
       print_the_help
       shift # past argument
       ;;
+    *)     # positional options
+      POSITIONAL+=("$1") # save it in an array for later
+      shift # past argument
+      ;;
   esac
 done
 set -- "${POSITIONAL[@]}" # restore positional parameters
@@ -62,17 +66,6 @@ set -- "${POSITIONAL[@]}" # restore positional parameters
 
 echo "view12 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  2005 ${INPUT_FILE} ${FILE_TAG}b_temp0.prim 
-dawncut 0  0 -1 -2000 ${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
-
 original_file_tag="${FILE_TAG}"
 
 make_slice(){
@@ -90,9 +83,10 @@ make_slice(){
   rm "${FILE_TAG}_temp0.prim"
   rm "${FILE_TAG}.prim"
 }
-for zzz in $(seq 50 50 2000) ;
+
+for zzz in $@ ;
 do
-  make_slice ${zzz} & 
+  make_slice ${zzz}
 done
 
 wait

scripts/view13/generate_eps

@@ -9,7 +9,7 @@ function print_the_help {
   exit 
 }
 
-FILE_TAG="view2"
+FILE_TAG="view13"
 INPUT_FILE="g4_0000.prim"
 
 

scripts/view14/.DAWN_1.history

+0.0
+72.5
+189.5
+0
+0
+0
+1500
+55
+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/view14/generate_eps

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

scripts/view15/.DAWN_1.history

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

scripts/view15/generate_eps

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

scripts/view2/.DAWN_1.history

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

scripts/view20/.DAWN_1.history

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

scripts/view20/generate_eps

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

scripts/view5/.DAWN_1.history

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

scripts/view6/.DAWN_1.history

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

scripts/view6/generate_eps

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

src/BarrelBarDetectorWithSideFrame_geo.cpp

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

src/BarrelCalorimeterInterlayers_geo.cpp

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