diff --git a/compact/definitions.xml b/compact/definitions.xml
index 908de8a96403addcbb4cd8ff73ac60d9724fbd04..1452c29fdf3adff06a037064afea2ff367f0f6e3 100644
--- a/compact/definitions.xml
+++ b/compact/definitions.xml
@@ -356,16 +356,16 @@ Examples:
     <constant name="TrackerBarrel_zmax"         value="TrackerBarrel_length/2.0"/>
     <constant name="TrackerBarrelInside_length" value="VertexTrackingRegion_length"/>
     <constant name="TrackerBarrelInside_zmax"   value="TrackerBarrelInside_length/2.0"/>
-    <constant name="TrackerEndcapP_length"      value="675.0*mm"/>h
+    <constant name="TrackerEndcapP_length"      value="925.0*mm"/>h
     <constant name="TrackerEndcapN_length"      value="495.0*mm"/>h
 
     <documentation level="0">
 ### PID Detector Region Parameters
     </documentation>
 
-    <constant name="ForwardPID_length"            value="180.0*cm"/>
-    <constant name="ForwardPID_rmin1"             value="Beampipe_rmax + 80*mm"/>
-    <constant name="ForwardPID_rmin2"             value="19.0*cm"/>
+    <constant name="ForwardPID_length"            value="145.0*cm"/>
+    <constant name="ForwardPID_rmin1"             value="Beampipe_rmax + 90*mm"/>
+    <constant name="ForwardPID_rmin2"             value="18.0*cm"/>
 
     <constant name="BackwardPID_length"           value="40.0*cm"/>
     <constant name="BackwardPID_rmax"             value="TrackerBarrel_rmax"/>
@@ -424,7 +424,7 @@ Service gaps in FW direction (before endcapP ECAL) and BW direction (before endc
 
     <constant name="EcalBarrelEnvelope_thickness"      value="40.0*cm"/>
     <constant name="EcalBarrel_rmin"                value="CentralTracking_rmax + BarrelPIDThickness + BarrelExtraSpaceThickness"/>
-    <constant name="EcalBarrelForward_length"       value="4*cm"/>
+    <constant name="EcalBarrelForward_length"       value="-21*cm"/>
     <constant name="EcalBarrelForward_zmax"         value="ForwardPID_zmin + EcalBarrelForward_length"/>
     <constant name="EcalBarrelBackward_zmax"        value="BackwardPID_zmin + BackwardInnerEndcap_length + EcalEndcapN_length"/>
     <constant name="EcalBarrel_length"              value="EcalBarrelForward_zmax + EcalBarrelBackward_zmax"/>
diff --git a/compact/display.xml b/compact/display.xml
index 30d508d4d65b50a16156e924fa66b4ad7219ecdc..8c5edf5dbbfc6fa5b15b0617de7dedb3e86251e6 100644
--- a/compact/display.xml
+++ b/compact/display.xml
@@ -94,7 +94,7 @@
     <vis name="DRICH_aerogel_vis" ref="AnlTeal" showDaughters="true" visible="true" />
     <vis name="DRICH_filter_vis"  alpha="1.0" r="1.0" g="1.0" b="0.0" showDaughters="true" visible="true" />
     <vis name="DRICH_mirror_vis"  ref="AnlGray" showDaughters="true" visible="true" />
-    <vis name="DRICH_sensor_vis"  ref="AnlGreen" showDaughters="true" visible="true" />
+    <vis name="DRICH_sensor_vis"  ref="AnlBlue" showDaughters="true" visible="true" />
 
     <vis name="MRICH_aerogel_vis" ref="AnlTeal" showDaughters="true" visible="true" />
     <vis name="MRICH_frame_vis" ref="AnlGold" showDaughters="true" visible="true" />
diff --git a/compact/drich.xml b/compact/drich.xml
index ab7e3ff6995571745c0728485685577303cba092..dde6a2890ab106838ffeda97402b7720c8b8c128 100644
--- a/compact/drich.xml
+++ b/compact/drich.xml
@@ -1,194 +1,260 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <lccdd>
 
-  <define>
-    <!-- TODO [low priority]: some of these, viz. radii, could be parameterized
-	 with respect to other variables; for now they are hard coded in case
-	 other detectors' parameters are changing -->
-    <!-- parameters for re-scaling fun4all design to ATHENA -->
-    <constant name="DRICH_scale"              value="0.963"/> <!-- overall scale factor from fun4all to ATHENA -->
-    <constant name="DRICH_f4a_length"         value="161.0*cm"/> <!-- z-length of fun4all design -->
-    <!-- vessel (=snout+tank) geometry -->
-    <constant name="DRICH_zmin"               value="ForwardPID_zmin"/> <!-- vessel front -->
-    <constant name="DRICH_Length"             value="ForwardPID_length"/>  <!-- overall vessel length (including snout) -->
-    <constant name="DRICH_rmin0"              value="ForwardPID_rmin1"/>  <!-- bore radius at dRICh vessel frontplane -->
-    <constant name="DRICH_rmin1"              value="19.0*cm"/>  <!-- bore radius at dRICh vessel backplane -->
-    <constant name="DRICH_wall_thickness"     value="0.5*cm"/>  <!-- thickness of radial walls -->
-    <constant name="DRICH_window_thickness"   value="0.1*cm"/>  <!-- thickness of entrance and exit walls -->
-    <!-- tank geometry: cylinder, holding the majority of detector components -->
-    <constant name="DRICH_rmax2"              value="200*cm"/>  <!-- cylinder radius; 20 cm gap between dRICh and HCalBarrel -->
-    <!-- snout geometry: cone with front radius rmax0 and back radius of rmax1 -->
-    <constant name="DRICH_SnoutLength"        value="50.0*cm"/>
-    <constant name="DRICH_SnoutSlope"         value="DRICH_rmax2 / (DRICH_zmin + DRICH_Length)"/>
-    <constant name="DRICH_rmax0"              value="DRICH_SnoutSlope * DRICH_zmin"/>
-    <constant name="DRICH_rmax1"              value="DRICH_SnoutSlope * ( DRICH_zmin + DRICH_SnoutLength)"/>
-    <!-- additional parameters -->
-    <constant name="DRICH_aerogel_thickness"  value="4.0*cm"/>  <!-- aerogel thickness -->
-    <constant name="DRICH_sensor_size"        value="48.0*mm"/> <!-- sensor side length -->
-    <constant name="DRICH_sensor_thickness"   value="35.0*mm"/> <!-- sensor thickness -->
-    <constant name="DRICH_num_px"             value="16"/> <!-- number of pixels along one side of the sensor -->
-  </define>
-
-  <detectors>
-
-    <detector
-      id="ForwardRICH_ID"
-      name="DRICH"
-      type="athena_DRICH"
-      readout="DRICHHits"
-      gas="C2F6_DRICH"
-      material="Aluminum"
-      vis_vessel="DRICH_vessel_vis"
-      vis_gas="DRICH_gas_vis"
-      >
-
-      <!-- envelope dimensions (see above)
-	   - `wall_thickness`: thickness of radial walls
-	   - `window_thickness`: thickness of entrance and exit disks
-	   -->
-      <dimensions
-        z0="DRICH_zmin"
-        length="DRICH_Length"
-        snout_length="DRICH_SnoutLength"
-        rmin0="DRICH_rmin0"
-        rmin1="DRICH_rmin1"
-        rmax0="DRICH_rmax0"
-        rmax1="DRICH_rmax1"
-        rmax2="DRICH_rmax2"
-        nsectors="6"
-        wall_thickness="DRICH_wall_thickness"
-        window_thickness="DRICH_window_thickness"
-        />
-
-      <!-- radiator defined in a wedge of azimuthal space
-           - `phiw` is phi width of wedge
-           - `thickness` defined separately for aerogel and filter
-           - `frontplane` is the front of the aerogel, w.r.t. front plane of vessel envelope
-           - `pitch` controls the angle of the radiator (0=vertical)
-           - filter is applied to backplane of aerogel
-           -->
-      <radiator
-        rmin="DRICH_rmin0 + DRICH_wall_thickness + 2.0*cm"
-        rmax="DRICH_rmax0 - DRICH_wall_thickness - 2.0*cm"
-        phiw="56*degree"
-        frontplane="DRICH_window_thickness + 0.5*DRICH_aerogel_thickness"
-        pitch="0*degree"
-        >
-        <aerogel
-          material="Aerogel_DRICH"
-          vis="DRICH_aerogel_vis"
-          thickness="DRICH_aerogel_thickness"
-          />
-        <filter
-          material="Acrylic_DRICH"
-          vis="DRICH_filter_vis"
-          thickness="0.3*mm"
-          />
-      </radiator>
-
-
-      <!-- spherical mirror is part of a sphere
-           - `rmin` and `rmax` provide polar angle boundaries
-           - `phiw` is the azimuthal width
-           - `radius` is the radius of the sphere
-           - `centerx` is the transverse position of the center 
-             of the sphere, for the sector on the +x axis
-           - the back of the mirror will pass through `backplane`
-           - set `debug` to 1 so draw reference sphere instead, view with y-clipping
-           -->
-      <mirror
-        material="Acrylic_DRICH"
-        surface="MirrorSurface_DRICH"
-        vis="DRICH_mirror_vis"
-        backplane="DRICH_Length-2.0*cm"
-        thickness="0.2*cm"
-        radius="290*DRICH_scale*cm"
-        centerx="145*DRICH_scale*cm"
-        rmin="DRICH_rmin1 + DRICH_wall_thickness + 0.0*cm"
-        rmax="DRICH_rmax2 - DRICH_wall_thickness - 2.0*cm"
-        phiw="54*degree"
-        debug="0"
-        />
-
-      <sensors>
-        <!-- geometry for a single square sensor
-             - based on Hamamatsu H13700 MAPMT 
-               (https://www.hamamatsu.com/us/en/product/type/H13700/index.html)
-               - N.B. not ideal for a magnetic field, SiPM matrix would be better
-               - effective area: 48.5x48.5mm
-               - enclosure size: 52x52mm
-               - 16x16 channel matrix (see readout segmentation below)
-               - pixel size: 3x3mm
-             - `side` is the side length of the square
-             - `thickness` is the depth of the sensor
-             - `gap` provides room between the squares, to help
-               prevent them from overlapping
-             - the value of `side` will determine how many sensors there are,
-               since the sensor placement algorithm will try to place as many
-               as it can in the specified patch below
-             -->
-        <module
-          material="Silicon"
-          surface="SensorSurface_DRICH"
-          vis="DRICH_sensor_vis"
-          side="DRICH_sensor_size"
-          thickness="DRICH_sensor_thickness"
-          gap="0.5*(52-48)*mm + 2*mm"
-          />
-        <!-- sensors will be tiled on this sphere
-             - `center{x,y,z} is defined for sector on +x axis, defined w.r.t. snout frontplane,
-               and `radius` is the sphere radius; the first term of each of these comes from
-               the fun4all design
-             - these attributes were determined from a spherical fit to the
-               sensor placement in the fun4all port
-             - set `debug` to 1 so draw reference sphere instead, view with y-clipping
-             -->
-        <sphere
-          radius="159.76*DRICH_scale*cm"
-          centerx="144.91*DRICH_scale*cm"
-          centery="0*DRICH_scale*cm"
-          centerz="-197.25*DRICH_scale*cm + DRICH_Length - 0.5*DRICH_scale*DRICH_f4a_length"
-          debug="0"
-          />
-        <!-- sensors will be limited to a patch of the sphere
-             - `thetamin` and `thetamax` define pseudorapidity coverage
-             - `widthfactor` controls the azimuthal coverage, where lower=wider
-             - `taper` defines half the angle between the azimuthal boundaries
-             - the size of the sensor controls how many sensors are placed
-             -->
-        <sphericalpatch
-          thetamin="-10*degree"
-          thetamax="22*degree"
-          widthfactor="1.8"
-          taper="56*degree"
-          />
-      </sensors>
-
-    </detector>
-  </detectors>
-
-  <readouts>
-    <readout name="DRICHHits">
-    <!-- segmentation: square matrix of pixels
-         - note: for `grid_size`, divide sensor size by 1 less than the
-           number of pixels, to account for fenceposting
-         -->
-      <segmentation
-	 type="CartesianGridXY"
-	 grid_size_x="DRICH_sensor_size/(DRICH_num_px-1)"
-	 grid_size_y="DRICH_sensor_size/(DRICH_num_px-1)"
-	 offset_x="-DRICH_sensor_size/2.0"
-	 offset_y="-DRICH_sensor_size/2.0"
-	 />
-      <!-- cellID: 64bits
-	   - offset 0,  length 8: dRICh ID
-	   - offset 8,  length 3: sector number
-	   - offset 11, length 12: photosensor number
-	   - offset 23, length 16: x pixel
-	   - offset 39, length 16: y pixel
-           -->
-      <id>system:8,sector:3,module:12,x:23:16,y:16</id>
-    </readout>
-  </readouts>
+<define>
+<!-- vessel (=snout+tank) geometry -->
+<constant name="DRICH_zmin"               value="ForwardPID_zmin"/> <!-- vessel front -->
+<constant name="DRICH_Length"             value="ForwardPID_length"/>  <!-- overall vessel length (including snout) -->
+<constant name="DRICH_rmin0"              value="ForwardPID_rmin1"/>  <!-- bore radius at dRICh vessel frontplane -->
+<constant name="DRICH_rmin1"              value="ForwardPID_rmin2"/>  <!-- bore radius at dRICh vessel backplane -->
+<constant name="DRICH_wall_thickness"     value="0.5*cm"/>  <!-- thickness of radial walls -->
+<constant name="DRICH_window_thickness"   value="0.1*cm"/>  <!-- thickness of entrance and exit walls -->
+<!-- tank geometry: cylinder, holding the majority of detector components -->
+<constant name="DRICH_rmax2"              value="200*cm"/>  <!-- cylinder radius; 20 cm gap between dRICh and HCalBarrel -->
+<!-- snout geometry: cone with front radius rmax0 and back radius of rmax1 -->
+<constant name="DRICH_SnoutLength"        value="25.0*cm"/>
+<constant name="DRICH_SnoutSlope"         value="DRICH_rmax2 / (DRICH_zmin + DRICH_Length)"/>
+<constant name="DRICH_rmax0"              value="DRICH_SnoutSlope * DRICH_zmin"/>
+<constant name="DRICH_rmax1"              value="DRICH_SnoutSlope * ( DRICH_zmin + DRICH_SnoutLength)"/>
+<!-- additional parameters -->
+<constant name="DRICH_aerogel_thickness"  value="4.0*cm"/>  <!-- aerogel thickness -->
+<constant name="DRICH_sensor_size"        value="48.0*mm"/> <!-- sensor side length -->
+<constant name="DRICH_sensor_thickness"   value="35.0*mm"/> <!-- sensor thickness -->
+<constant name="DRICH_num_px"             value="16"/> <!-- number of pixels along one side of the sensor -->
+<!-- debugging switches -->
+<comment>
+- `DRICH_debug_optics`:  1 = all components become vacuum, except for mirrors; test opticalphotons from IP
+                         2 = all components become vacuum, except for mirrors and `gasvol`, test charged particles from IP
+                         0 = off
+- `DRICH_debug_mirror`:  1 = draw full mirror shape for single sector; 0 = off
+- `DRICH_debug_sensors`: 1 = draw full sensor sphere for a single sector; 0 = off
+</comment>
+<constant name="DRICH_debug_optics"  value="0"/>
+<constant name="DRICH_debug_mirror"  value="0"/>
+<constant name="DRICH_debug_sensors" value="0"/>
+</define>
+
+
+<detectors>
+
+
+<!-- /detectors/detector -->
+<documentation level="10">
+### dRICh: ***d***ual ***R***ing ***I***maging ***Ch***erenkov detector
+</documentation>
+<detector
+  id="ForwardRICH_ID"
+  name="DRICH"
+  type="athena_DRICH"
+  readout="DRICHHits"
+  gas="C2F6_DRICH"
+  material="Aluminum"
+  vis_vessel="DRICH_vessel_vis"
+  vis_gas="DRICH_gas_vis"
+  debug_optics="DRICH_debug_optics"
+  >
+
+
+<!-- /detectors/detector/dimensions -->
+<documentation level="10">
+#### Vessel
+- the dRICh vessel is composed of two parts:
+  - tank: cylindrical region containing most of the detector components
+  - snout: conical region at the front of the vessel, containing the aerogel
+- dimensions:
+  - `zmin`: z-position of vessel front plane
+  - `length`: overall z-length of the full vessel
+  - `snout_length`: length of cone-shaped snout region, housing aerogel
+  - `rmin0` and `rmin1`: bore radius at front plane and back plane, respectively
+  - `rmax0` and `rmax1`: outer radius of snout at front plane and snout-back (tank-front) plane, respectively
+  - `rmax2`: outer radius of tank, the main cylindrical vessel volume
+  - `nsectors`: number of azimuthal sectors
+  - `wall_thickness`: thickness of radial walls
+  - `window_thickness`: thickness of entrance and exit disks
+</documentation>
+<dimensions
+  zmin="DRICH_zmin"
+  length="DRICH_Length"
+  snout_length="DRICH_SnoutLength"
+  rmin0="DRICH_rmin0"
+  rmin1="DRICH_rmin1"
+  rmax0="DRICH_rmax0"
+  rmax1="DRICH_rmax1"
+  rmax2="DRICH_rmax2"
+  nsectors="6"
+  wall_thickness="DRICH_wall_thickness"
+  window_thickness="DRICH_window_thickness"
+  />
+
+
+<!-- /detectors/detector/radiator -->
+<documentation level="10">
+#### Radiator
+- radiator is defined in a wedge of azimuthal space, composed of aerogel and a
+  filter; the filter is applied to the back of the aerogel, so that it separates
+  the aerogel and gas radiators
+- dimensions:
+  - `phiw`: azimuthal width of wedge
+  - `thickness`: radiator thickness, defined separately for aerogel and filter
+  - `frontplane`: front of the aerogel, w.r.t. front plane of the vessel envelope
+  - `pitch`: controls the angle of the radiator (0=vertical)
+</documentation>
+<radiator
+  rmin="DRICH_rmin0 + DRICH_wall_thickness + 2.0*cm"
+  rmax="DRICH_rmax0 - DRICH_wall_thickness - 2.0*cm"
+  phiw="60*degree"
+  frontplane="DRICH_window_thickness + 0.5*DRICH_aerogel_thickness"
+  pitch="0*degree"
+  >
+  <aerogel
+    material="Aerogel_DRICH"
+    vis="DRICH_aerogel_vis"
+    thickness="DRICH_aerogel_thickness"
+    />
+  <filter
+    material="Acrylic_DRICH"
+    vis="DRICH_filter_vis"
+    thickness="0.3*mm"
+    />
+</radiator>
+
+
+<!-- /detectors/detector/mirror -->
+<documentation level="10">
+#### Spherical mirror
+- spherical mirrors are built from spherical patches, and positioned near the
+  vessel back plane, separately for each sector
+- dimensions:
+  - `backplane`: the position of the maximum z-plane intersected by the sphere,
+    w.r.t. the back plane of vessel envelope
+  - `rmin` and `rmax`: polar angle boundaries
+  - `phiw`: azimuthal width of one sector
+  - `thickness` is the radial thickness of the mirror; note that `backplane` is given for the 
+    reflective mirror surface, the inner radius of the sphere
+  - `focus_tune*` are tuning parameters for the focal plane:
+    - linearly retune focal point, using a real number; reference numbers are:
+      - 0.0: disabled, ideal focus at sensor sphere center
+      - 1.0: ideal focus re-tuned to be on sensor sphere
+      - `focus_tune_long` moves toward centroid sensor, while `focus_tune_perp` moves transverse
+        to this direction; movements are in the xz-plane only
+      - due to aberrations, ideal focus may not coincide with apparent focal region
+- other settings:
+  - `debug`: set to 1 so draw reference sphere instead, view with y-clipping
+</documentation>
+<mirror
+  material="Acrylic_DRICH"
+  surface="MirrorSurface_DRICH"
+  vis="DRICH_mirror_vis"
+  backplane="DRICH_window_thickness + 5.0*cm"
+  rmin="DRICH_rmin1 + DRICH_wall_thickness - 1.0*cm"
+  rmax="DRICH_rmax2 - DRICH_wall_thickness - 5.0*cm"
+  phiw="59.5*degree"
+  thickness="0.2*cm"
+  focus_tune_long="0.5"
+  focus_tune_perp="0.2"
+  debug="DRICH_debug_mirror"
+  />
+
+
+<!-- /detectors/detector/sensors -->
+<documentation level="10">
+#### Sensors
+</documentation>
+<sensors>
+
+
+<!-- /detectors/detector/sensors/module -->
+<documentation level="10">
+##### Sensor module
+- based on [Hamamatsu H13700 MAPMT](https://www.hamamatsu.com/us/en/product/type/H13700/index.html):
+  - not ideal for a magnetic field, SiPM matrix would be better
+  - effective area: 48.5x48.5 mm
+  - enclosure size: 52x52 mm
+  - 16x16 channel matrix (cf. readout segmentation below)
+  - pixel size: 3x3 mm
+- dimensions:
+  - `side`: side length of the square module
+  - `thickness`: thickness of the sensor module
+  - `gap`: provides room between the squares, to help prevent them from overlapping
+  - note: the value of `side` will determine how many sensors there are, since the
+    sensor placement algorithm will try to place as many as it can in the specified
+    spherical patch below
+</documentation>
+<module
+  material="Silicon"
+  surface="SensorSurface_DRICH"
+  vis="DRICH_sensor_vis"
+  side="DRICH_sensor_size"
+  thickness="DRICH_sensor_thickness"
+  gap="0.5*(52-48)*mm + 2*mm"
+  />
+
+
+<!-- /detectors/detector/sensors/{sphere,sphericalpatch} -->
+<documentation level="10">
+##### Sensor sphere
+- sensors will be placed on a sphere, using a "disco ball" tiling algorithm; each
+  sector has its own sensor sphere
+  - sphere dimensions:
+    - `centerx` and `centerz`: sphere center, defined w.r.t. vessel front plane,
+      for the sector on +x axis
+    - `radius`: radius of the sensor sphere
+  - other settings:
+    - `debug`: set to 1 so draw reference sphere instead, view with y-clipping
+- sensors will be limited to a patch of the sphere
+  - patch dimensions:
+    - `phiw`: defines half the angle between the azimuthal boundaries
+    - `rmin` and `rmax`: radial cut boundaries
+    - `zmin`: z-plane cut
+</documentation>
+<sphere
+  centerz="-112.0 * cm"
+  centerx="DRICH_rmax2 - 35.0*cm"
+  radius="160.0 * cm"
+  debug="DRICH_debug_sensors"
+  />
+<sphericalpatch
+  phiw="18*degree"
+  rmin="DRICH_rmax1 + 10.0*cm"
+  rmax="DRICH_rmax2 - 5.0*cm"
+  zmin="DRICH_SnoutLength + 5.0*cm"
+  />
+
+
+</sensors>
+</detector>
+</detectors>
+
+
+<documentation level="10">
+#### Readout
+- segmentation: square matrix of pixels
+  - `grid_size_x,y`: size of each sensor, but note we must divide sensor size
+    by 1 less than the number of pixels, to account for fenceposting
+  - `offset_x,y`: specified such that the `x` and `y` indicators are unsigned
+- indicators and `cellID` bits:
+
+  | indicator | offset | length |
+  |-----------|--------|--------|
+  | dRICh ID  | 0      | 8      |
+  | sector    | 8      | 3      |
+  | sensor    | 11     | 12     |
+  | x pixel   | 23     | 16     |
+  | y pixel   | 39     | 16     |
+
+</documentation>
+<readouts>
+  <readout name="DRICHHits">
+    <segmentation
+      type="CartesianGridXY"
+      grid_size_x="DRICH_sensor_size/(DRICH_num_px-1)"
+      grid_size_y="DRICH_sensor_size/(DRICH_num_px-1)"
+      offset_x="-DRICH_sensor_size/2.0"
+      offset_y="-DRICH_sensor_size/2.0"
+      />
+    <id>system:8,sector:3,module:12,x:23:16,y:16</id>
+  </readout>
+</readouts>
+
 
 </lccdd>
diff --git a/compact/optical_materials.xml b/compact/optical_materials.xml
index 66bd455a29f9abc461aefa46178c9e0fd8aad605..bcb78ee4d9ef72ba1456c56e634a1dcbbc11403c 100644
--- a/compact/optical_materials.xml
+++ b/compact/optical_materials.xml
@@ -675,6 +675,14 @@
       <property name="RINDEX" ref="RINDEX__Air"/>
       <property name="ABSLENGTH"               coldim="2" values="1*eV  200*m  5*eV  200*m"/>
     </material>
+    <material name="VacuumOptical">
+      <D type="density" unit="g/cm3" value="0.0000000001"/>
+      <fraction n="0.754" ref="N"/>
+      <fraction n="0.234" ref="O"/>
+      <fraction n="0.012" ref="Ar"/>
+      <property name="RINDEX" ref="RINDEX__Vacuum"/>
+      <property name="ABSLENGTH"               coldim="2" values="1*eV  2000*m  5*eV  2000*m"/>
+    </material>
     <material name="N2cherenkov">
       <D type="density" value="0.00125" unit="g/cm3"/>
       <composite n="1" ref="N"/>
diff --git a/compact/subsystem_views/drich_only.xml b/compact/subsystem_views/drich_only.xml
new file mode 100644
index 0000000000000000000000000000000000000000..f16c1460792f64811304f164d08a9a293c1bd2d7
--- /dev/null
+++ b/compact/subsystem_views/drich_only.xml
@@ -0,0 +1,161 @@
+<lccdd xmlns:compact="http://www.lcsim.org/schemas/compact/1.0"
+       xmlns:xs="http://www.w3.org/2001/XMLSchema"
+       xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/compact/1.0/compact.xsd">
+
+       <debug>
+         <type name="surface"       value="0"/>
+         <type name="material"      value="0"/>
+         <type name="readout"       value="0"/>
+         <type name="segmentation"  value="0"/>
+         <type name="limits"        value="0"/>
+         <type name="region"        value="0"/>
+         <type name="includes"      value="0"/>
+       </debug>
+
+  <documentation level="-1">
+  # Athena Detector 
+  - https://eicweb.phy.anl.gov/EIC/detectors/athena.git 
+  - https://eicweb.phy.anl.gov/EIC/detectors/ip6.git
+  </documentation>
+
+  <!-- Some information about detector  -->
+  <info name="Athena Detector" title="Athena Detector"
+        author="Athena Collaboration"
+	url="https://eicweb.phy.anl.gov/EIC/detectors/athena.git"
+	status="development"
+	version="v1 2021-03-16">
+  <comment> Athena </comment>
+  </info>
+  <define>
+  <documentation level="2">
+      ## Main Constant Definitions
+
+      The ip6 (or other ip) defines should be included first.
+      These files have only a define tags.
+  </documentation>
+    <include ref="ip6/ip6_defs.xml" /> 
+    <include ref="compact/definitions.xml" />
+  </define>
+
+  <includes>
+    <gdmlFile ref="compact/elements.xml"/>
+    <gdmlFile ref="compact/materials.xml"/>
+    <file     ref="compact/optical_materials.xml"/>
+  </includes>
+
+  <limits>
+    <limitset name="EICBeamlineLimits">
+      <limit name="step_length_max" particles="*" value="1.0" unit="mm" />
+      <limit name="track_length_max" particles="*" value="1.0" unit="mm" />
+      <limit name="time_max" particles="*" value="0.1" unit="ns" />
+      <limit name="ekin_min" particles="*" value="0.001" unit="MeV" />
+      <limit name="range_min" particles="*" value="0.1" unit="mm" />
+    </limitset>
+    <limitset name="cal_limits">
+      <limit name="step_length_max" particles="*" value="5.0" unit="mm"/>
+    </limitset>
+  </limits>
+
+  <display>
+  <include ref="compact/colors.xml" />
+  <!--include ref="compact/colors2.xml"/-->
+  <include ref="compact/display.xml" />
+  <!--include ref="compact/display_detailed.xml"/-->
+  </display>
+
+  <documentation level="0">
+    ## Detector Subsystems
+
+    ### IP Subsystems
+
+    The interaction point subsystems are included before the central detector subsystems.
+    This is becuase the IP subsystems, for exmaple the beampipe, will define paramters
+    which are subsquently used in the central detector construction -- e.g. the vertex tracker
+    uses the beampipe OD to help define its placement. 
+
+    The IP subsystems include the Far forward and backward regions. The list of subsystem includes:
+     - Interaction region beampipe 
+     - B0 tracker
+     - Off-momentum tracker
+     - Far forward roman pots
+     - Zero Degree Calorimeter
+     - Beam line magnets.
+     - and more...
+  </documentation>
+
+  <documentation level="10">
+  ### dRICh only
+  </documentation>
+  <include ref="compact/drich.xml" />
+
+  <fields>
+    <field name="B0PF_Magnet" type="MultipoleMagnet">
+      <position x="B0PF_XPosition" y="0" z="B0PF_CenterPosition"/>
+      <rotation x="0" y="B0PF_RotationAngle" z="0"/>
+      <shape type="Tube" rmin="0.0" rmax="B0PF_InnerRadius" dz="B0PF_Length*0.5"/>
+      <coefficient coefficient="B0PF_Bmax" skew="0.0*tesla"/>
+      <!--<coefficient coefficient="2.0*tesla/cm" skew="0.2*tesla/cm"/> -->
+    </field>
+    <field name="B0APF_Magnet" type="MultipoleMagnet">
+      <position x="B0APF_XPosition" y="0" z="B0APF_CenterPosition"/>
+      <rotation x="0" y="B0APF_RotationAngle" z="0"/>
+      <shape type="Tube" rmin="0.0" rmax="B0APF_InnerRadius" dz="B0APF_Length*0.5"/>
+      <coefficient coefficient="B0APF_Bmax" skew="0.0*tesla"/>
+      <!--<coefficient coefficient="2.0*tesla/cm" skew="0.2*tesla/cm"/> -->
+    </field>
+    <field name="Q1APF_Magnet" type="MultipoleMagnet">
+      <position x="Q1APF_XPosition" y="0" z="Q1APF_CenterPosition"/>
+      <rotation x="0" y="Q1APF_RotationAngle" z="0"/>
+      <shape type="Tube" rmin="0.0" rmax="Q1APF_InnerRadius" dz="Q1APF_Length*0.5"/>
+      <coefficient coefficient="Q1APF_Bmax" skew="0.0*tesla"/>
+      <coefficient coefficient="Q1APF_GradientMax" skew="0.0*tesla/cm"/>
+    </field>
+    <field name="Q1BPF_Magnet" type="MultipoleMagnet">
+      <position x="Q1BPF_XPosition" y="0" z="Q1BPF_CenterPosition"/>
+      <rotation x="0" y="Q1BPF_RotationAngle" z="0"/>
+      <shape type="Tube" rmin="0.0" rmax="Q1BPF_InnerRadius" dz="Q1BPF_Length*0.5"/>
+      <coefficient coefficient="Q1BPF_Bmax" skew="0.0*tesla"/>
+      <coefficient coefficient="Q1BPF_GradientMax" skew="0.0*tesla/cm"/>
+    </field>
+    <field name="Q2PF_Magnet" type="MultipoleMagnet">
+      <position x="Q2PF_XPosition" y="0" z="Q2PF_CenterPosition"/>
+      <rotation x="0" y="Q2PF_RotationAngle" z="pi/2.0"/>
+      <shape type="Tube" rmin="0.0" rmax="Q2PF_InnerRadius" dz="Q2PF_Length*0.5"/>
+      <coefficient coefficient="Q2PF_Bmax" skew="0.0*tesla"/>
+      <coefficient coefficient="Q2PF_GradientMax" skew="0.0*tesla/cm"/>
+    </field>
+    <field name="B1PF_Magnet" type="MultipoleMagnet">
+      <position x="B1PF_XPosition" y="0" z="B1PF_CenterPosition"/>
+      <rotation x="0" y="B1PF_RotationAngle" z="0"/>
+      <shape type="Tube" rmin="0.0" rmax="B1PF_InnerRadius" dz="B1PF_Length*0.5"/>
+      <coefficient coefficient="B1PF_Bmax" skew="0.0*tesla"/>
+      <coefficient coefficient="B1PF_GradientMax" skew="0.0*tesla/cm"/>
+    </field>
+    <field name="B1APF_Magnet" type="MultipoleMagnet">
+      <position x="B1APF_XPosition" y="0" z="B1APF_CenterPosition"/>
+      <rotation x="0" y="B1APF_RotationAngle" z="0"/>
+      <shape type="Tube" rmin="0.0" rmax="B1APF_InnerRadius" dz="B1APF_Length*0.5"/>
+      <coefficient coefficient="B1APF_Bmax" skew="0.0*tesla"/>
+      <coefficient coefficient="B1APF_GradientMax" skew="0.0*tesla/cm"/>
+    </field>
+    <field name="B2PF_Magnet" type="MultipoleMagnet">
+      <position x="B2PF_XPosition" y="0" z="B2PF_CenterPosition"/>
+      <rotation x="0" y="B2PF_RotationAngle" z="0"/>
+      <shape type="Tube" rmin="0.0" rmax="B2PF_InnerRadius" dz="B2PF_Length*0.5"/>
+      <coefficient coefficient="B2PF_Bmax" skew="0.0*tesla"/>
+      <coefficient coefficient="B2PF_GradientMax" skew="0.0*tesla/cm"/>
+    </field>
+  </fields>
+
+  <comment>
+      FB elements
+      -----------
+      None (TODO)
+
+      What is FB?
+  </comment>
+
+  <readouts>
+  </readouts>
+
+</lccdd>
diff --git a/src/DRich_geo.cpp b/src/DRich_geo.cpp
index 645d2965805b0555261ee33a41e8b5683980196a..e6225d7b950e5ebfdfce96496504c1d5ad615b21 100644
--- a/src/DRich_geo.cpp
+++ b/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,86 +34,102 @@ 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  focusTuneLong    =  mirrorElem.attr<double>(_Unicode(focus_tune_long));
+  double  focusTunePerp    =  mirrorElem.attr<double>(_Unicode(focus_tune_perp));
   // - 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;
   Cone vesselSnout(
@@ -139,14 +154,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,26 +169,50 @@ 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
@@ -183,97 +222,66 @@ static Ref_t createDetector(Detector& desc, xml::Handle_t handle, SensitiveDetec
   // SECTOR LOOP //////////////////////////////////
   for(int isec=0; isec<nSectors; isec++) {
 
-    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
+    // debugging filters, limiting the number of sectors
+    if( (debug_mirror||debug_sensors||debug_optics) && isec!=0) continue;
 
-    // BUILD RADIATOR ====================================================================
+    // sector rotation about z axis
+    double sectorRotation = isec * 360/nSectors * degree;
+    std::string secName = "sec" + std::to_string(isec);
 
-    // derived attributes
-    auto radiatorPos = Position(0., 0., radiatorFrontplane) + snoutFront;
+
+    // BUILD RADIATOR ====================================================================
 
     // 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);
+    Volume aerogelVol( detName+"_aerogel_"+secName, aerogelSolid, aerogelMat );
+    Volume filterVol(  detName+"_filter_"+secName,  filterSolid,  filterMat );
     aerogelVol.setVisAttributes(aerogelVis);
     filterVol.setVisAttributes(filterVis);
 
-    // placement
+    // 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,
-	  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
-	);
-    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
+          RotationZ(sectorRotation) // rotate about beam axis to sector
+        * Translation3D(radiatorPos.x(), radiatorPos.y(), radiatorPos.z()) // re-center to originFront
+        * RotationY(radiatorPitch) // change polar angle to specified pitch
+        );
     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);
-    filterDE.setPlacement(filterPV);
-    //SkinSurface filterSkin(desc, filterDE, Form("mirror_optical_surface%d", isec), filterSurf, filterVol);
-    //filterSkin.isValid();
-
 
-    // BUILD MIRROR ====================================================================
-
-    // 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);
-
-    // 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);
-
-    // 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();
+    // 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 originFront
+          * RotationY(radiatorPitch) // change polar angle
+          * Translation3D(0., 0., (aerogelThickness+filterThickness)/2.) // move to aerogel backplane
+          );
+      DetElement filterDE(det, Form("filter_de%d", isec), isec);
+      filterDE.setPlacement(filterPV);
+      //SkinSurface filterSkin(desc, filterDE, Form("mirror_optical_surface%d", isec), filterSurf, filterVol);
+      //filterSkin.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 +307,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 +329,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 +365,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);
           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 +383,121 @@ 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 ====================================================================
+
+    // 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
+
+    // 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
+    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;
+    };
+
+    // 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);
+
+    /*
+    printf("(zC,xC) = ( %.2f, %.2f )\n",zC,xC);
+    printf("zS = %f\n",zS);
+    printf("xS = %f\n",xS);
+    printf("B = %f\n",B);
+    printf("zM = %f\n",zM);
+    printf("xM = %f\n",xM);
+    printf("rM = %f\n",rM);
+    */
+
+    // 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
+        );
+
+    // 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 +509,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);