From 579ca7b220535e34b08f7d7a1aa4906c7a061203 Mon Sep 17 00:00:00 2001
From: Wouter Deconinck <wouter.deconinck@umanitoba.ca>
Date: Sun, 12 Sep 2021 23:01:03 +0000
Subject: [PATCH] Resolve "mRICH alternative implementation"

---
 compact/mrich.xml     |   3 +-
 compact/mrich_alt.xml | 112 ++++++++++
 src/MRich_geo.cpp     | 481 ++++++++++++++++++++++--------------------
 3 files changed, 372 insertions(+), 224 deletions(-)
 create mode 100644 compact/mrich_alt.xml

diff --git a/compact/mrich.xml b/compact/mrich.xml
index 91fce35e..703e620d 100644
--- a/compact/mrich.xml
+++ b/compact/mrich.xml
@@ -77,7 +77,8 @@
               length="MRICHCarbonFrame_length">
         <frame vis="AnlGray" thickness="MRICHCarbonFrame_thickness" material="CarbonFiber"/>
         <aerogel vis="MRICH_aerogel_vis"
-                 length="MRICHAerogel_thickness" width="MRICHAerogel_width"
+                 length="MRICHAerogel_thickness"
+                 width="MRICHAerogel_width"
                  material="AerogelOptical">
           <frame vis="MRICH_frame_vis" thickness="MRICHFoam_thickness" material="PolystyreneFoam" />
         </aerogel>
diff --git a/compact/mrich_alt.xml b/compact/mrich_alt.xml
new file mode 100644
index 00000000..b12897eb
--- /dev/null
+++ b/compact/mrich_alt.xml
@@ -0,0 +1,112 @@
+<lccdd>
+  <comment> MRICH (alternative design) </comment>        
+  <define>
+    <constant name="MRICH_rmin"   value="10*cm"/>
+    <constant name="MRICH_rmax"   value="BackwardPID_rmax"/>
+    <constant name="MRICH_length" value="BackwardPID_length"/>
+    <constant name="MRICH_zmin"   value="BackwardPID_zmin"/>
+
+    <constant name="MRICHAerogel_thickness"           value="30.0*mm"/>
+    <constant name="MRICHAerogel_width"               value="126.5*mm"/>
+    <constant name="MRICHFoam_thickness"              value="2*mm"/>
+    <constant name="MRICHFresnelLens_thickness"       value="0.06*inch"/>
+    <constant name="MRICHAerogelLensGap_thickness"    value="2*mm"/>
+    <constant name="MRICHPhotoDet_thickness"          value="1.5*mm"/>
+    <constant name="MRICHPhotoDet_length"             value="48.5*mm"/>
+    <constant name="MRICHGlassWindow_width"           value="103.5*mm"/>
+    <constant name="MRICHGlassPhotoDet_thickness"     value="5.0*mm"/>
+    <constant name="MRICHPhotoDetMCPlate_thickness"   value="0.3*mm"/> <!-- FIXME: should be 1.2*mm with PyrexGlass25 -->
+    <constant name="MRICHPhotoDetAnode_thickness"     value="3.8*mm"/>
+    <constant name="MRICHPhotoDetPCB_thickness"       value="2.0*mm"/>
+    <constant name="MRICHPhotoDetCopper_thickness"    value="0.1*mm"/>
+    <constant name="MRICHPhotoDetKapton_thickness"    value="0.2*mm"/>
+    <constant name="MRICHRearExtraSpace_thickness"    value="0.3*mm"/>
+    <constant name="MRICHLensPhotoDet_length"         value="136.4*mm"/>
+
+    <constant name="MRICHFresnelLensEffectiveDiameter" value="6.0*inch"/>
+    <constant name="MRICHFresnelLensGroove_pitch"      value="inch/125"/>
+
+    <constant name="MRICHCarbonFrame_thickness"        value="1.0*mm"/>
+    <constant name="MRICHCarbonFrame_width"            value="MRICHAerogel_width+2.0*MRICHFoam_thickness + 2.0*MRICHCarbonFrame_thickness"/>
+
+    <constant name="MRICHModules_nx"                  value="floor((MRICH_rmax-MRICH_rmin)/MRICHCarbonFrame_width)"/>
+    <constant name="MRICHModules_ny"                  value="floor((MRICH_rmax-MRICH_rmin)/MRICHCarbonFrame_width)"/>
+
+    <constant name="MRICHCarbonFrame_length"
+      value="MRICHAerogel_thickness
+      + 2.0*MRICHCarbonFrame_thickness
+      + 2.0*MRICHFoam_thickness
+      + MRICHAerogelLensGap_thickness
+      + MRICHFresnelLens_thickness
+      + MRICHLensPhotoDet_length
+      + MRICHGlassPhotoDet_thickness
+      + 2.0*MRICHPhotoDetMCPlate_thickness
+      + MRICHPhotoDetAnode_thickness
+      + MRICHPhotoDetPCB_thickness
+      + MRICHPhotoDetCopper_thickness
+      + MRICHPhotoDetKapton_thickness
+      + MRICHRearExtraSpace_thickness "/>
+  </define>
+
+  <limits>
+  </limits>
+
+  <regions>
+  </regions>
+
+  <display>
+  </display>
+
+  <detectors>
+    <detector id="MRICH_ID" name="MRICH" type="athena_MRICH" 
+      readout="MRICHHits"
+      reflect="true"
+      projective="false"
+      vis="InvisibleWithDaughters"
+      material="Air">
+      <dimensions rmin="MRICH_rmin" rmax="MRICH_rmax" length="abs(MRICH_length)" zmin="MRICH_zmin"/>
+      <envelope thickness="MRICHCarbonFrame_thickness" material="CarbonFiber"/>
+      <module name="MRICH_module1"
+              vis="InvisibleWithDaughters"
+              width="MRICHCarbonFrame_width"
+              height="MRICHCarbonFrame_width"
+              length="MRICHCarbonFrame_length">
+        <aerogel vis="MRICH_aerogel_vis"
+                 length="MRICHAerogel_thickness"
+                 width="MRICHAerogel_width"
+                 material="AerogelOptical">
+          <frame vis="MRICH_frame_vis" thickness="MRICHFoam_thickness" material="PolystyreneFoam" />
+        </aerogel>
+        <lens vis="MRICH_lens_vis" thickness="MRICHFresnelLens_thickness"
+              pitch="MRICHFresnelLensGroove_pitch" focal_length="6.0*inch"
+              effective_diameter="MRICHFresnelLensEffectiveDiameter" 
+              width="MRICHAerogel_width"
+              material="AcrylicOptical"/>
+        <space thickness="MRICHLensPhotoDet_length"/>
+        <photodet width="MRICHGlassWindow_width" thickness="MRICHGlassPhotoDet_thickness"  material="PyrexGlassOptical">
+          <sensor nx="2" ny="2" thickness="MRICHPhotoDet_thickness" width="MRICHPhotoDet_length" material="SiliconOxide"/>
+          <layer thickness="MRICHPhotoDetMCPlate_thickness" material="PyrexGlass"/> <!-- FIXME: should be PyrexGlass25 with 1.2*mm thickness -->
+          <layer thickness="MRICHPhotoDetMCPlate_thickness" material="PyrexGlass"/> <!-- FIXME: should be PyrexGlass25 with 1.2*mm thickness -->
+          <layer thickness="MRICHPhotoDetAnode_thickness"   material="AluminumOxide"/>
+          <layer thickness="MRICHPhotoDetPCB_thickness"     material="Fr4"/>
+          <layer thickness="MRICHPhotoDetCopper_thickness" material="Copper"/>
+          <layer thickness="MRICHPhotoDetKapton_thickness" material="Kapton"/>
+        </photodet>
+      </module>
+    </detector>
+  </detectors>
+
+  <readouts>
+    <readout name="MRICHHits">
+      <segmentation type="CartesianGridXY" grid_size_x="3*mm" grid_size_y="3*mm" />
+      <id>system:8,module:14,sensor:8,x:32:-16,y:-16</id>
+    </readout>
+  </readouts>
+  
+  <!--Globals>
+    <Parameter Name="mrichInfo" Value="mrichmod/mrich_1_geoparams-0-0-4294967295-1527211159.xml"/>
+  </Globals-->
+
+
+
+</lccdd>
diff --git a/src/MRich_geo.cpp b/src/MRich_geo.cpp
index e00c2ece..c3c98668 100644
--- a/src/MRich_geo.cpp
+++ b/src/MRich_geo.cpp
@@ -32,19 +32,6 @@ static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDet
   sens.setType("photoncounter");
   OpticalSurfaceManager surfMgr = description.surfaceManager();
 
-  // read module positions
-  std::vector<std::tuple<double,double,double>> positions;
-  for (xml_coll_t x_positions_i(x_det, _Unicode(positions)); x_positions_i; ++x_positions_i) {
-    xml_comp_t x_positions = x_positions_i;
-    for (xml_coll_t x_position_i(x_positions, _U(position)); x_position_i; ++x_position_i) {
-      xml_comp_t x_position = x_position_i;
-      positions.push_back(
-        std::make_tuple(x_positions.scale() * x_position.x() * mm,
-                        x_positions.scale() * x_position.y() * mm,
-                        -x_positions.z0()));
-    }
-  }
-
   bool projective = getAttrOrDefault(x_det, _Unicode(projective), false);
   bool reflect    = x_det.reflect(true);
 
@@ -57,6 +44,7 @@ static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDet
 
   int                     n_sensor = 1;
 
+  // dimensions
   xml::Component dims   = x_det.dimensions();
   auto           rmin   = dims.rmin();
   auto           rmax   = dims.rmax();
@@ -64,6 +52,20 @@ static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDet
   auto           zmin   = dims.zmin();
   auto           zpos   = zmin + length / 2;
 
+  // envelope
+  Tube   envShape(rmin, rmax, length / 2., 0., 2 * M_PI);
+  Volume envVol("MRICH_Envelope", envShape, air);
+  envVol.setVisAttributes(description.visAttributes(x_det.visStr()));
+  if (x_det.hasChild(_Unicode(envelope))) {
+    xml_comp_t x_envelope = x_det.child(_Unicode(envelope));
+    double thickness  = x_envelope.thickness();
+    Material material = description.material(x_envelope.materialStr());
+    Tube   envInsideShape(rmin + thickness, rmax - thickness, length / 2. - thickness);
+    SubtractionSolid envShellShape(envShape, envInsideShape);
+    Volume envShell("MRICH_Envelope_Inside", envShellShape, material);
+    envVol.placeVolume(envShell);
+  }
+
   // expect only one module (for now)
   xml_comp_t x_mod = x_det.child(_U(module));
   string     mod_name            = x_mod.nameStr();
@@ -71,200 +73,231 @@ static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDet
   double     mod_height          = getAttrOrDefault(x_mod, _U(height), 130.0 * mm);
   double     mod_length          = getAttrOrDefault(x_mod, _U(length), 130.0 * mm);
 
-  // various components
-  xml_comp_t x_frame    = x_mod.child(_Unicode(frame));
-  xml_comp_t x_aerogel  = x_mod.child(_Unicode(aerogel));
-  xml_comp_t x_lens     = x_mod.child(_Unicode(lens));
-  xml_comp_t x_mirror   = x_mod.child(_Unicode(mirror));
-  xml_comp_t x_photodet = x_mod.child(_Unicode(photodet));
-
   // module
   Box    m_solid(mod_width / 2.0, mod_height / 2.0, mod_length / 2.0);
   Volume m_volume(mod_name, m_solid, air);
   m_volume.setVisAttributes(description.visAttributes(x_mod.visStr()));
   DetElement mod_de( mod_name + std::string("_mod_") + std::to_string(1), 1);
+  double z_placement = - mod_length / 2.0;
 
   // todo module frame
-  double     frame_thickness = getAttrOrDefault(x_frame, _U(thickness), 2.0 * mm);
-  Box        frame_inside(mod_width / 2.0 - frame_thickness, mod_height / 2.0 - frame_thickness, mod_length / 2.0 - frame_thickness);
-  SubtractionSolid frame_solid(m_solid, frame_inside);
-  Material   frame_mat       = description.material(x_frame.materialStr());
-  Volume     frame_vol(mod_name+"_frame", frame_solid, frame_mat);
-  auto       frame_vis       = getAttrOrDefault<std::string>(x_frame, _U(vis), std::string("GrayVis"));
-  frame_vol.setVisAttributes(description.visAttributes(frame_vis));
-  m_volume.placeVolume(frame_vol);
+  if (x_mod.hasChild(_Unicode(frame))) {
+    xml_comp_t x_frame    = x_mod.child(_Unicode(frame));
+    double     frame_thickness = getAttrOrDefault(x_frame, _U(thickness), 2.0 * mm);
+    Box        frame_inside(mod_width / 2.0 - frame_thickness, mod_height / 2.0 - frame_thickness, mod_length / 2.0 - frame_thickness);
+    SubtractionSolid frame_solid(m_solid, frame_inside);
+    Material   frame_mat       = description.material(x_frame.materialStr());
+    Volume     frame_vol(mod_name+"_frame", frame_solid, frame_mat);
+    auto       frame_vis       = getAttrOrDefault<std::string>(x_frame, _U(vis), std::string("GrayVis"));
+    frame_vol.setVisAttributes(description.visAttributes(frame_vis));
+    // update position
+    z_placement += frame_thickness / 2.0;
+    // place volume
+    m_volume.placeVolume(frame_vol);
+    // update position
+    z_placement += frame_thickness / 2.0;
+  }
 
   // aerogel box
-  xml_comp_t x_aerogel_frame = x_aerogel.child(_Unicode(frame));
-  double     aerogel_width   = getAttrOrDefault(x_aerogel, _U(width), 130.0 * mm);
-  double     aerogel_length  = getAttrOrDefault(x_aerogel, _U(length), 130.0 * mm);
-  double     foam_thickness  = getAttrOrDefault(x_aerogel_frame, _U(thickness), 2.0 * mm);
-  Material   foam_mat        = description.material(x_aerogel_frame.materialStr());
-  Material   aerogel_mat     = description.material(x_aerogel.materialStr());
-  auto       aerogel_vis     = getAttrOrDefault<std::string>(x_aerogel, _U(vis), std::string("InvisibleWithDaughters"));
-  auto       foam_vis        = getAttrOrDefault<std::string>(x_aerogel_frame, _U(vis), std::string("RedVis"));
-
-  // aerogel foam frame
-  Box foam_box(aerogel_width / 2.0 + foam_thickness, aerogel_width / 2.0 + foam_thickness, (aerogel_length + foam_thickness) / 2.0);
-  Box aerogel_sub_box(aerogel_width / 2.0, aerogel_width / 2.0, (aerogel_length + foam_thickness) / 2.0);
-  SubtractionSolid foam_frame_solid(foam_box, aerogel_sub_box, Position(0, 0, foam_thickness));
-  Volume           foam_vol(mod_name+"_aerogel_frame", foam_frame_solid, foam_mat);
-  foam_vol.setVisAttributes(description.visAttributes(foam_vis));
-  double foam_frame_zpos = -mod_length / 2.0 + frame_thickness + (aerogel_length + foam_thickness) / 2.0;
-  m_volume.placeVolume(foam_vol,Position(0,0,foam_frame_zpos));
-
-  // aerogel
-  Box              aerogel_box(aerogel_width / 2.0, aerogel_width / 2.0, (aerogel_length) / 2.0);
-  Volume           aerogel_vol(mod_name+"_aerogel", aerogel_box, aerogel_mat);
-  aerogel_vol.setVisAttributes(description.visAttributes(aerogel_vis));
-  double aerogel_zpos = foam_frame_zpos + foam_thickness / 2.0;
-  pv = m_volume.placeVolume(aerogel_vol,Position(0,0,aerogel_zpos));
-  DetElement aerogel_de(mod_de, mod_name + std::string("_aerogel_de") + std::to_string(1), 1);
-  aerogel_de.setPlacement(pv);
-
-  auto aerogel_surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault(x_aerogel, _Unicode(surface), "MRICH_AerogelOpticalSurface"));
-  SkinSurface skin0(description, aerogel_de, Form("MRICH_aerogel_skin_surface_%d", 1), aerogel_surf, aerogel_vol);
-  skin0.isValid();
-
-  // Fresnel Lens
-  //  - The lens has a constant groove pitch (delta r) as opposed to fixing the groove height.
-  //  - The lens area outside of the effective diamtere is flat.
-  //  - The grooves are not curved, rather they are polycone shaped, ie a flat approximating the curvature.
-  auto   lens_vis       = getAttrOrDefault<std::string>(x_lens, _U(vis), std::string("AnlBlue"));
-  double groove_pitch   = getAttrOrDefault(x_lens, _Unicode(pitch), 0.2 * mm);// 0.5 * mm);
-  double lens_f         = getAttrOrDefault(x_lens, _Unicode(focal_length), 6.0*2.54*cm);
-  double eff_diameter   = getAttrOrDefault(x_lens, _Unicode(effective_diameter), 152.4 * mm);
-  double lens_width     = getAttrOrDefault(x_lens, _Unicode(width), 6.7*2.54*cm);
-  double center_thickness = getAttrOrDefault(x_lens, _U(thickness), 0.068 * 2.54 * cm);//2.0 * mm);
-
-  double n_acrylic        = 1.49;
-  double lens_curvature   = 1.0 / (lens_f*(n_acrylic - 1.0)); //confirmed
-  double full_ring_rmax   = std::min(eff_diameter / 2.0, lens_width/2.0);
-
-  double N_grooves        = std::ceil((full_ring_rmax) / groove_pitch);
-  double groove_last_rmin = (N_grooves - 1) * groove_pitch;
-  double groove_last_rmax = N_grooves * groove_pitch;
-
-  auto   groove_sagitta = [&](double r) { return lens_curvature * std::pow(r, 2) / (1.0 + 1.0); };
-  double lens_thickness = groove_sagitta(groove_last_rmax) - groove_sagitta(groove_last_rmin) + center_thickness;
-
-  Material         lens_mat = description.material(x_lens.materialStr());
-  Box              lens_box(lens_width / 2.0, lens_width / 2.0, (center_thickness) / 2.0);
-  SubtractionSolid flat_lens(lens_box, Tube(0.0, full_ring_rmax, 2 * center_thickness));
-
-  Assembly lens_vol(mod_name + "_lens");
-  Volume   flatpart_lens_vol( "flatpart_lens", flat_lens, lens_mat);
-  lens_vol.placeVolume(flatpart_lens_vol);//,Position(0,0,lens_zpos));
-
-  Solid  fresnel_lens_solid;
-
-  int    i_groove           = 0;
-  double groove_rmax        = groove_pitch;
-  double groove_rmin        = 0;
-
-  while ( groove_rmax <= full_ring_rmax ) {
-    double   dZ = groove_sagitta(groove_rmax) - groove_sagitta(groove_rmin);
-    Polycone groove_solid(0, 2.0 * M_PI,
-                          {groove_rmin, groove_rmin, groove_rmin},
-                          {groove_rmax, groove_rmax, groove_rmin},
-                          {-lens_thickness/2.0, lens_thickness/2.0-dZ, lens_thickness/2.0});
-    Volume   lens_groove_vol("lens_groove_" + std::to_string(i_groove), groove_solid, lens_mat);
-    lens_vol.placeVolume(lens_groove_vol);
-    //Volume   groove_vol(groove_solid, lens_mat, par->name.c_str(), 0, 0, 0);
-    //new G4PVPlacement(0, par->pos, Groove_log[i], par->name.c_str(), motherLV, false, 0, OverlapCheck());
-    //phi1 = phi1 + halfpi;  //g4 pre-defined: halfpi=pi/2
-    //Tube       sub_cylinder(r0, r1, 3*eff_diameter);
-    //IntersectionSolid groove_solid(lens_box,lens_sphere, Position(0,0,-eff_diameter/2.0 + lens_thickness/2.0+(t-lens_t)/2.0 ));
-    //IntersectionSolid lens_ring(groove_solid, sub_cylinder);
-    //if (i_groove == 0) {
-    //  fresnel_lens_solid = groove_solid;
-    //} else {
-    //  fresnel_lens_solid = UnionSolid(fresnel_lens_solid, groove_solid);
-    //}
-    //r0 = r1;
-    //if(i_groove > 3) { 
-    //  SubtractionSolid flat_lens(lens_box,Tube(0.0, r0, 3*eff_diameter));
-    //  fresnel_lens_solid = UnionSolid(fresnel_lens_solid, flat_lens);
-    //  break; // temporary
-    //}
-    i_groove++;
-    groove_rmin = (i_groove  )*groove_pitch;
-    groove_rmax = (i_groove+1)*groove_pitch;
+  if (x_mod.hasChild(_Unicode(aerogel))) {
+    xml_comp_t x_aerogel  = x_mod.child(_Unicode(aerogel));
+    double     aerogel_width   = getAttrOrDefault(x_aerogel, _U(width), 130.0 * mm);
+    double     aerogel_length  = getAttrOrDefault(x_aerogel, _U(length), 130.0 * mm);
+    Material   aerogel_mat     = description.material(x_aerogel.materialStr());
+    auto       aerogel_vis     = getAttrOrDefault<std::string>(x_aerogel, _U(vis), std::string("InvisibleWithDaughters"));
+
+    xml_comp_t x_aerogel_frame = x_aerogel.child(_Unicode(frame));
+    double     foam_thickness  = getAttrOrDefault(x_aerogel_frame, _U(thickness), 2.0 * mm);
+    Material   foam_mat        = description.material(x_aerogel_frame.materialStr());
+    auto       foam_vis        = getAttrOrDefault<std::string>(x_aerogel_frame, _U(vis), std::string("RedVis"));
+
+    // foam frame
+    Box foam_box(aerogel_width / 2.0 + foam_thickness, aerogel_width / 2.0 + foam_thickness, (aerogel_length + foam_thickness) / 2.0);
+    Box foam_sub_box(aerogel_width / 2.0, aerogel_width / 2.0, (aerogel_length + foam_thickness) / 2.0);
+    SubtractionSolid foam_frame_solid(foam_box, foam_sub_box, Position(0, 0, foam_thickness));
+    Volume           foam_vol(mod_name+"_aerogel_frame", foam_frame_solid, foam_mat);
+    foam_vol.setVisAttributes(description.visAttributes(foam_vis));
+
+    // aerogel
+    Box              aerogel_box(aerogel_width / 2.0, aerogel_width / 2.0, (aerogel_length) / 2.0);
+    Volume           aerogel_vol(mod_name+"_aerogel", aerogel_box, aerogel_mat);
+    aerogel_vol.setVisAttributes(description.visAttributes(aerogel_vis));
+
+    // update position
+    z_placement += (aerogel_length + foam_thickness) / 2.0;
+    // place foam frame
+    pv = m_volume.placeVolume(foam_vol,Position(0,0,z_placement));
+    // place aerogel
+    z_placement += foam_thickness / 2.0;
+    pv = m_volume.placeVolume(aerogel_vol,Position(0,0,z_placement));
+    DetElement aerogel_de(mod_de, mod_name + std::string("_aerogel_de") + std::to_string(1), 1);
+    aerogel_de.setPlacement(pv);
+    // update position
+    z_placement += aerogel_length / 2.0;
+
+    // optical surfaces
+    auto aerogel_surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault(x_aerogel, _Unicode(surface), "MRICH_AerogelOpticalSurface"));
+    SkinSurface skin_surf(description, aerogel_de, Form("MRICH_aerogel_skin_surface_%d", 1), aerogel_surf, aerogel_vol);
+    skin_surf.isValid();
   }
 
-  //fresnel_lens_solid = UnionSolid(fresnel_lens_solid, flat_lens);
-  //Volume   lens_vol(mod_name + "_lens", fresnel_lens_solid, lens_mat);
+  // Fresnel Lens
+  if (x_mod.hasChild(_Unicode(lens))) {
+    xml_comp_t x_lens     = x_mod.child(_Unicode(lens));
+
+    //  - The lens has a constant groove pitch (delta r) as opposed to fixing the groove height.
+    //  - The lens area outside of the effective diamtere is flat.
+    //  - The grooves are not curved, rather they are polycone shaped, ie a flat approximating the curvature.
+    auto   lens_vis       = getAttrOrDefault<std::string>(x_lens, _U(vis), std::string("AnlBlue"));
+    double groove_pitch   = getAttrOrDefault(x_lens, _Unicode(pitch), 0.2 * mm);// 0.5 * mm);
+    double lens_f         = getAttrOrDefault(x_lens, _Unicode(focal_length), 6.0*2.54*cm);
+    double eff_diameter   = getAttrOrDefault(x_lens, _Unicode(effective_diameter), 152.4 * mm);
+    double lens_width     = getAttrOrDefault(x_lens, _Unicode(width), 6.7*2.54*cm);
+    double center_thickness = getAttrOrDefault(x_lens, _U(thickness), 0.068 * 2.54 * cm);//2.0 * mm);
+
+    double n_acrylic        = 1.49;
+    double lens_curvature   = 1.0 / (lens_f*(n_acrylic - 1.0)); //confirmed
+    double full_ring_rmax   = std::min(eff_diameter / 2.0, lens_width/2.0);
+
+    double N_grooves        = std::ceil((full_ring_rmax) / groove_pitch);
+    double groove_last_rmin = (N_grooves - 1) * groove_pitch;
+    double groove_last_rmax = N_grooves * groove_pitch;
+
+    auto   groove_sagitta = [&](double r) { return lens_curvature * std::pow(r, 2) / (1.0 + 1.0); };
+    double lens_thickness = groove_sagitta(groove_last_rmax) - groove_sagitta(groove_last_rmin) + center_thickness;
+
+    Material         lens_mat = description.material(x_lens.materialStr());
+    Box              lens_box(lens_width / 2.0, lens_width / 2.0, (center_thickness) / 2.0);
+    SubtractionSolid flat_lens(lens_box, Tube(0.0, full_ring_rmax, 2 * center_thickness));
+
+    Assembly lens_vol(mod_name + "_lens");
+    Volume   flatpart_lens_vol( "flatpart_lens", flat_lens, lens_mat);
+    lens_vol.placeVolume(flatpart_lens_vol);
+
+    int    i_groove           = 0;
+    double groove_rmax        = groove_pitch;
+    double groove_rmin        = 0;
+
+    while ( groove_rmax <= full_ring_rmax ) {
+      double   dZ = groove_sagitta(groove_rmax) - groove_sagitta(groove_rmin);
+      Polycone groove_solid(0, 2.0 * M_PI,
+                            {groove_rmin, groove_rmin, groove_rmin},
+                            {groove_rmax, groove_rmax, groove_rmin},
+                            {-lens_thickness/2.0, lens_thickness/2.0-dZ, lens_thickness/2.0});
+      Volume   lens_groove_vol("lens_groove_" + std::to_string(i_groove), groove_solid, lens_mat);
+      lens_vol.placeVolume(lens_groove_vol);
+
+      i_groove++;
+      groove_rmin = (i_groove  )*groove_pitch;
+      groove_rmax = (i_groove+1)*groove_pitch;
+    }
 
-  lens_vol.setVisAttributes(description.visAttributes(lens_vis));
-  double lens_zpos = aerogel_zpos +aerogel_length/ 2.0 + foam_thickness + lens_thickness/2.0;
-  pv = m_volume.placeVolume(lens_vol,Position(0,0,lens_zpos));
-  DetElement lens_de(mod_de, mod_name + std::string("_lens_de") + std::to_string(1), 1);
-  lens_de.setPlacement(pv);
+    lens_vol.setVisAttributes(description.visAttributes(lens_vis));
+
+    // update position
+    z_placement += lens_thickness/2.0;
+    // place volume
+    pv = m_volume.placeVolume(lens_vol,Position(0,0,z_placement));
+    DetElement lens_de(mod_de, mod_name + std::string("_lens_de") + std::to_string(1), 1);
+    lens_de.setPlacement(pv);
+    // update position
+    z_placement += lens_thickness/2.0;
+
+    // optical surfaces
+    auto lens_surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault(x_lens, _Unicode(surface), "MRICH_LensOpticalSurface"));
+    SkinSurface skin_surf(description, lens_de, Form("MRichFresnelLens_skin_surface_%d", 1), lens_surf, lens_vol);
+    skin_surf.isValid();
+  }
 
-  auto surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault(x_lens, _Unicode(surface), "MRICH_LensOpticalSurface"));
-  SkinSurface skin(description, lens_de, Form("MRichFresnelLens_skin_surface_%d", 1), surf, lens_vol);
-  skin.isValid();
+  // mirror
+  if (x_mod.hasChild(_Unicode(space))) {
+    xml_comp_t x_space = x_mod.child(_Unicode(space));
+    z_placement += getAttrOrDefault(x_space, _U(thickness), 0.0 * mm);
+  }
 
   // mirror
-  auto   mirror_vis      = getAttrOrDefault<std::string>(x_mirror, _U(vis), std::string("AnlGray"));
-  double mirror_x1      = getAttrOrDefault(x_mirror, _U(x1), 100.0 * mm);
-  double mirror_x2      = getAttrOrDefault(x_mirror, _U(x2), 80.0 * mm);
-  double mirror_length  = getAttrOrDefault(x_mirror, _U(length), 130.0 * mm);
-  double mirror_thickness  = getAttrOrDefault(x_mirror, _U(thickness), 2.0 * mm);
-  double outer_x1 = (mirror_x1+mirror_thickness)/2.0;
-  double outer_x2 = (mirror_x2+mirror_thickness)/2.0;
-  Trd2   outer_mirror_trd(outer_x1, outer_x2, outer_x1,  outer_x2, mirror_length/2.0);
-  Trd2   inner_mirror_trd(mirror_x1 / 2.0,  mirror_x2 / 2.0, mirror_x1 / 2.0,mirror_x2 / 2.0, mirror_length/2.0+0.1*mm);
-  SubtractionSolid mirror_solid(outer_mirror_trd, inner_mirror_trd);
-  Material   mirror_mat        = description.material(x_mirror.materialStr());
-  Volume           mirror_vol(mod_name+"_mirror", mirror_solid, mirror_mat);
-  double mirror_zpos = lens_zpos + lens_thickness/2.0 + foam_thickness + mirror_length/2.0;
-  pv = m_volume.placeVolume(mirror_vol,Position(0,0,mirror_zpos));
-  DetElement mirror_de(mod_de, mod_name + std::string("_mirror_de") + std::to_string(1), 1);
-  mirror_de.setPlacement(pv);
-
-  auto mirror_surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault(x_mirror, _Unicode(surface), "MRICH_MirrorOpticalSurface"));
-  SkinSurface skin1(description, mirror_de, Form("MRICH_mirror_skin_surface_%d", 1), mirror_surf, mirror_vol);
-  skin1.isValid();
+  if (x_mod.hasChild(_Unicode(mirror))) {
+    xml_comp_t x_mirror   = x_mod.child(_Unicode(mirror));
+    auto   mirror_vis     = getAttrOrDefault<std::string>(x_mirror, _U(vis), std::string("AnlGray"));
+    double mirror_x1      = getAttrOrDefault(x_mirror, _U(x1), 100.0 * mm);
+    double mirror_x2      = getAttrOrDefault(x_mirror, _U(x2), 80.0 * mm);
+    double mirror_length  = getAttrOrDefault(x_mirror, _U(length), 130.0 * mm);
+    double mirror_thickness  = getAttrOrDefault(x_mirror, _U(thickness), 2.0 * mm);
+    double outer_x1 = (mirror_x1+mirror_thickness)/2.0;
+    double outer_x2 = (mirror_x2+mirror_thickness)/2.0;
+    Trd2   outer_mirror_trd(outer_x1, outer_x2, outer_x1,  outer_x2, mirror_length/2.0);
+    Trd2   inner_mirror_trd(mirror_x1 / 2.0,  mirror_x2 / 2.0, mirror_x1 / 2.0,mirror_x2 / 2.0, mirror_length/2.0+0.1*mm);
+    SubtractionSolid mirror_solid(outer_mirror_trd, inner_mirror_trd);
+    Material mirror_mat        = description.material(x_mirror.materialStr());
+    Volume   mirror_vol(mod_name+"_mirror", mirror_solid, mirror_mat);
+
+    // update position
+    z_placement += mirror_length/2.0;
+    // place volume
+    pv = m_volume.placeVolume(mirror_vol,Position(0,0,z_placement));
+    DetElement mirror_de(mod_de, mod_name + std::string("_mirror_de") + std::to_string(1), 1);
+    mirror_de.setPlacement(pv);
+    // update position
+    z_placement += mirror_length/2.0;
+
+    // optical surfaces
+    auto mirror_surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault(x_mirror, _Unicode(surface), "MRICH_MirrorOpticalSurface"));
+    SkinSurface skin_surf(description, mirror_de, Form("MRICH_mirror_skin_surface_%d", 1), mirror_surf, mirror_vol);
+    skin_surf.isValid();
+  }
 
   // photon detector
-  xml_comp_t x_photodet_sensor  = x_photodet.child(_Unicode(sensor));
-  auto   photodet_vis      = getAttrOrDefault<std::string>(x_photodet, _U(vis), std::string("AnlRed"));
-  double     photodet_width     = getAttrOrDefault(x_photodet, _U(width), 130.0 * mm);
-  double     photodet_thickness = getAttrOrDefault(x_photodet, _U(thickness), 2.0 * mm);
-  double     sensor_thickness   = getAttrOrDefault(x_photodet_sensor, _U(thickness), 2.0 * mm);
-  Material   photodet_mat       = description.material(x_photodet.materialStr());
-  Material   sensor_mat         = description.material(x_photodet_sensor.materialStr());
-  int     sensor_nx   = getAttrOrDefault(x_photodet_sensor, _Unicode(nx), 2);
-  int     sensor_ny   = getAttrOrDefault(x_photodet_sensor, _Unicode(ny), 2);
-
-  Box    window_box(photodet_width/2.0,photodet_width/2.0,photodet_thickness/2.0);
-  Volume window_vol(mod_name+"_window", window_box, photodet_mat);
-  double window_zpos = mirror_zpos + mirror_length/2.0+photodet_thickness/2.0;
-  pv = m_volume.placeVolume(window_vol,Position(0,0,window_zpos));
-  DetElement   comp_de(mod_de, std::string("mod_sensor_de_") + std::to_string(1) ,  1);
-  comp_de.setPlacement(pv);
-  // sensitive
-  pv.addPhysVolID("sensor", n_sensor);
-  window_vol.setSensitiveDetector(sens);
-  sensitives[mod_name].push_back(pv);
-  ++n_sensor;
-
-  // photon detector electronics layers
-  double layer_zpos = window_zpos + photodet_thickness/2.0;
-  int i_layer = 1;
-  for (xml_coll_t li(x_photodet, _Unicode(layer)); li; ++li) {
-    xml_comp_t x_layer = li;
-    Material layer_mat = description.material(x_layer.materialStr());
-    double   layer_thickness = x_layer.thickness();
-    Box      layer_box(photodet_width/2.0,photodet_width/2.0,layer_thickness/2.0);
-    Volume   layer_vol(mod_name + "_layer_" + std::to_string(i_layer), layer_box, layer_mat);
-    layer_zpos += layer_thickness / 2.0;
-    pv = m_volume.placeVolume(layer_vol,Position(0,0,layer_zpos));
-    DetElement layer_de(mod_de, std::string("mod_layer_de_") + std::to_string(i_layer),  1);
-    layer_de.setPlacement(pv);
-    layer_zpos += layer_thickness / 2.0;
-    i_layer++;
+  if (x_mod.hasChild(_Unicode(photodet))) {
+    xml_comp_t x_photodet = x_mod.child(_Unicode(photodet));
+    auto       photodet_vis       = getAttrOrDefault<std::string>(x_photodet, _U(vis), std::string("AnlRed"));
+    double     photodet_width     = getAttrOrDefault(x_photodet, _U(width), 130.0 * mm);
+    double     photodet_thickness = getAttrOrDefault(x_photodet, _U(thickness), 2.0 * mm);
+    Material   photodet_mat       = description.material(x_photodet.materialStr());
+    Box        window_box(photodet_width/2.0,photodet_width/2.0,photodet_thickness/2.0);
+    Volume     window_vol(mod_name+"_window", window_box, photodet_mat);
+
+    // update position
+    z_placement += photodet_thickness/2.0;
+    // place volume
+    pv = m_volume.placeVolume(window_vol,Position(0,0,z_placement));
+    DetElement   comp_de(mod_de, mod_name + std::string("_sensor_de_") + std::to_string(1) ,  1);
+    comp_de.setPlacement(pv);
+    // update position
+    z_placement += photodet_thickness/2.0;
+
+    // sensitive
+    pv.addPhysVolID("sensor", n_sensor);
+    window_vol.setSensitiveDetector(sens);
+    sensitives[mod_name].push_back(pv);
+    ++n_sensor;
+
+    // sensor
+    xml_comp_t x_sensor  = x_photodet.child(_Unicode(sensor));
+    double     sensor_thickness   = getAttrOrDefault(x_sensor, _U(thickness), 2.0 * mm);
+    Material   sensor_mat         = description.material(x_sensor.materialStr());
+    int        sensor_nx          = getAttrOrDefault(x_sensor, _Unicode(nx), 2);
+    int        sensor_ny          = getAttrOrDefault(x_sensor, _Unicode(ny), 2);
+
+    // layers
+    int i_layer = 1;
+    for (xml_coll_t li(x_photodet, _Unicode(layer)); li; ++li) {
+      xml_comp_t x_layer = li;
+      Material layer_mat = description.material(x_layer.materialStr());
+      double   layer_thickness = x_layer.thickness();
+      Box      layer_box(photodet_width/2.0,photodet_width/2.0,layer_thickness/2.0);
+      Volume   layer_vol(mod_name + "_layer_" + std::to_string(i_layer), layer_box, layer_mat);
+
+      // update position
+      z_placement += layer_thickness / 2.0;
+      // place volume
+      pv = m_volume.placeVolume(layer_vol,Position(0,0,z_placement));
+      DetElement layer_de(mod_de, mod_name + std::string("_layer_de_") + std::to_string(i_layer),  1);
+      layer_de.setPlacement(pv);
+      // update position
+      z_placement += layer_thickness / 2.0;
+
+      i_layer++;
+    }
   }
 
   //for (size_t ic = 0; ic < sensVols.size(); ++ic) {
@@ -285,11 +318,6 @@ static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDet
   module_assembly_delements[mod_name] = mod_de;
   // end module
 
-  // detector envelope
-  Tube   envShape(rmin, rmax, length / 2., 0., 2 * M_PI);
-  Volume envVol("MRICH_Envelope", envShape, air);
-  envVol.setVisAttributes(description.visAttributes(x_det.visStr()));
-
   // place modules in the sectors (disk)
   auto points = athena::geo::fillSquares({0., 0.}, mod_width, rmin, rmax);
 
@@ -299,6 +327,29 @@ static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDet
   // determine module direction, always facing z = 0
   double roty = dims.zmin() < 0. ? -M_PI : 0 ;
 
+  // read module positions
+  std::vector<std::tuple<double,double,double>> positions;
+  for (xml_coll_t x_positions_i(x_det, _Unicode(positions)); x_positions_i; ++x_positions_i) {
+    xml_comp_t x_positions = x_positions_i;
+    for (xml_coll_t x_position_i(x_positions, _U(position)); x_position_i; ++x_position_i) {
+      xml_comp_t x_position = x_position_i;
+      positions.push_back(
+        std::make_tuple(x_positions.scale() * x_position.x() * mm,
+                        x_positions.scale() * x_position.y() * mm,
+                        -x_positions.z0()));
+    }
+  }
+  // if no positions, then autoplacement
+  if (positions.empty()) {
+    for (double x = mod_width / 2.0; x < rmax - mod_width / 2.0; x += mod_width) {
+      for (double y = mod_width / 2.0; y < rmax - mod_width / 2.0; y += mod_width) {
+        if (pow(x + mod_width / 2.0,2) + pow(y + mod_width / 2.0,2) > rmax*rmax) continue;
+        if (pow(x - mod_width / 2.0,2) + pow(y - mod_width / 2.0,2) < rmin*rmin) continue;
+        positions.push_back(std::make_tuple(x, y, 0));
+      }
+    }
+  }
+
   // place modules
   int i_mod = 1; // starts at 1
   for (auto& p: positions) {
@@ -316,29 +367,11 @@ static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDet
       double y = std::get<1>(p);
       double z0 = std::get<2>(p);
 
-      // get angles
-      double rotAngX = atan(y/z0);
-      double rotAngY = -1.*atan(x/z0);
-
-      /*
-      ROOT::Math::XYZVector x_location(p.x(), p.y(), zmin+std::signbit(zmin)*mod_length/2.0);
-      ROOT::Math::XYZVector z_dir(0, 0, 1);
-      ROOT::Math::XYZVector x_dir(1, 0, 0);
-      ROOT::Math::XYZVector rot_axis = x_location.Cross(z_dir);
-      double rot_angle = ROOT::Math::VectorUtil::Angle(z_dir,x_location);
-      ROOT::Math::AxisAngle proj_rot(rot_axis,-1.0*rot_angle);
-      ROOT::Math::AxisAngle grid_fix_rot(x_location,0.0*rot_angle);
-      auto new_x_dir = grid_fix_rot*x_dir;
-      // operations are inversely ordered
-      Transform3D tr = Translation3D(p.x(), p.y(), 0.) // move to position
-                      * RotationY(roty);              // facing z = 0.
-      */
-
       Transform3D tr;
       if(projective) {
-        tr = Translation3D(x, y, 0)
-           * RotationX(rotAngX)
-           * RotationY(rotAngY);
+        double rotAngX = atan(y/z0);
+        double rotAngY = -1.*atan(x/z0);
+        tr = Translation3D(x, y, 0) * RotationX(rotAngX) * RotationY(rotAngY);
       } else {
         tr = Translation3D(x, y, 0) * RotationX(0);
       }
@@ -355,13 +388,15 @@ static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDet
     }
   }
 
-  // place frame
-  xml_comp_t x_layer = x_det.child(_Unicode(layer));
-  double   layer_thickness = x_layer.thickness();
-  Material layer_mat = description.material(x_layer.materialStr());
-  Tube   frameShape(rmin, rmax, layer_thickness / 2., 0., 2 * M_PI);
-  Volume frameVol("MRICH_Frame", frameShape, layer_mat);
-  pv = envVol.placeVolume(frameVol, Position(0, 0, (length - layer_thickness) / 2.0));
+  // additional layers
+  if (x_det.hasChild(_Unicode(layer))) {
+    xml_comp_t x_layer = x_det.child(_Unicode(layer));
+    double   layer_thickness = x_layer.thickness();
+    Material layer_mat = description.material(x_layer.materialStr());
+    Tube   frameShape(rmin, rmax, layer_thickness / 2., 0., 2 * M_PI);
+    Volume frameVol("MRICH_Frame", frameShape, layer_mat);
+    pv = envVol.placeVolume(frameVol, Position(0, 0, (length - layer_thickness) / 2.0));
+  }
 
   // place envelope
   Volume motherVol = description.pickMotherVolume(sdet);
-- 
GitLab