From 1e592bce12c8b4dee3bb28d114524f79fa8fa3ea Mon Sep 17 00:00:00 2001
From: Sylvester Joosten <sjoosten@anl.gov>
Date: Thu, 9 Sep 2021 02:45:18 +0000
Subject: [PATCH] Added scintillating glass, 55cm long bars

(cherry picked from commit 5755020d531eea569cb2c9bd4781e4eaa0998827)
---
 .gitlab-ci.yml                |   5 ++
 compact/ci_ecal_scfi.xml      |   6 +-
 compact/display_geoviewer.xml |   2 +-
 compact/hybrid_ecal.xml       |   2 +-
 compact/materials.xml         |   7 +++
 src/HybridCalorimeter_geo.cpp |  33 ++++++-----
 src/ScFiCalorimeter_geo.cpp   | 106 +++++++++++++++++-----------------
 7 files changed, 88 insertions(+), 73 deletions(-)

diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 5f0625d6..21ec1910 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -108,6 +108,8 @@ dump_geometry:
   script:
     - echo "dumping geometry"
     - mkdir -p geo
+    ## disable fibers in SciFi endcap ECAL for better performance
+    - sed -i '/<fiber/,+4d' ${DETECTOR_PATH}/compact/ci_ecal_scfi.xml
     ## subsystem views
     - |
       mv ${DETECTOR_PATH}/compact/display_geoviewer.xml ${DETECTOR_PATH}/compact/display.xml
@@ -121,7 +123,10 @@ dump_geometry:
     - dd_web_display --output geo/detector_geo_full.root ${DETECTOR_PATH}/athena.xml
     ## only central detector geo
     - sed -i '/forward_ion_beamline/d' ${DETECTOR_PATH}/athena.xml
+    - sed -i '/beampipe_hadron_B0/d' ${DETECTOR_PATH}/athena.xml
+    - sed -i '/B0_/d' ${DETECTOR_PATH}/athena.xml
     - sed -i '/far_forward/d' ${DETECTOR_PATH}/athena.xml
+    - sed -i '/roman_pots/d' ${DETECTOR_PATH}/athena.xml
     - dd_web_display --output geo/detector_geo.root ${DETECTOR_PATH}/athena.xml
     ## print some useful output
     - |
diff --git a/compact/ci_ecal_scfi.xml b/compact/ci_ecal_scfi.xml
index a5559643..9662825c 100644
--- a/compact/ci_ecal_scfi.xml
+++ b/compact/ci_ecal_scfi.xml
@@ -25,16 +25,16 @@
       ------------------------------------------
       Forward (Positive Z) Endcap EM Calorimeter
       ------------------------------------------
-      An EM calorimeter with shashlik hexagon modules
+      An EM calorimeter with ScFi modules
     </comment>
     <detector id="ECalEndcapP_ID"
       name="EcalEndcapP"
       type="ScFiCalorimeter"
-      vis="InvisibleWithDaughters"
+      vis="EcalEndcapVis"
       readout="EcalEndcapPHits">
       <position x="0" y="0" z="EcalEndcapP_zmin + EcalEndcapP_length/2."/>
       <dimensions rmin="EcalEndcapP_rmin" rmax="EcalEndcapP_rmax" length="EcalEndcapP_length"/>
-      <module sizex="25*mm" sizey="25*mm" sizez="170*mm" material="TungstenDens24" vis="EcalEndcapVis">
+      <module sizex="25*mm" sizey="25*mm" sizez="170*mm" material="TungstenDens24" vis="EcalEndcapLayerVis">
         <fiber material="Polystyrene"
           radius="EcalEndcapP_FiberRadius"
           offset="EcalEndcapP_FiberOffset"
diff --git a/compact/display_geoviewer.xml b/compact/display_geoviewer.xml
index c452c817..60a5ea60 100644
--- a/compact/display_geoviewer.xml
+++ b/compact/display_geoviewer.xml
@@ -33,7 +33,7 @@
     <vis name="EcalBarrelSupportVis"    ref="AnlOrange"/>
 
     <vis name="EcalVis"                 ref="AnlGold"   showDaughters="true"  visible="false"/>
-    <vis name="EcalEndcapVis"           ref="AnlGold"   showDaughters="true" visible="true"/>
+    <vis name="EcalEndcapVis"           ref="AnlGold"   showDaughters="false" visible="true"/>
     <vis name="EcalEndcapLayerVis"      ref="AnlGold"   showDaughters="false" visible="true"/>
 
     <vis name="EcalEndcapNModuleVis"    ref="AnlGold"   showDaughters="false" visible="true"/>
diff --git a/compact/hybrid_ecal.xml b/compact/hybrid_ecal.xml
index 9b38bb81..9115722f 100644
--- a/compact/hybrid_ecal.xml
+++ b/compact/hybrid_ecal.xml
@@ -27,7 +27,7 @@
     <constant name="CrystalModule_width" value="20.00*mm"/>
     <constant name="CrystalModule_length" value="200.00*mm"/>
     <constant name="CrystalModule_wrap" value="0.50*mm"/>
-    <constant name="CrystalModule_z0" value="10.*cm"/>
+    <constant name="CrystalModule_z0" value="0.0*cm"/>
 
     <constant name="GlassModule_width" value="2*CrystalModule_width"/>
     <constant name="GlassModule_length" value="55.00*cm"/>
diff --git a/compact/materials.xml b/compact/materials.xml
index a1ab3a2a..d21627f1 100644
--- a/compact/materials.xml
+++ b/compact/materials.xml
@@ -105,6 +105,13 @@
     <fraction n="0.0278" ref="SodiumOxide"/>
     <fraction n="0.0050" ref="ArsenicOxide"/>
   </material>
+  <material name="SciGlass">
+    <D type="density" value="4.22" unit="g / cm3"/>
+    <fraction n="0.3875" ref="Ba"/>
+    <fraction n="0.2146" ref="Gd"/>
+    <fraction n="0.1369" ref="Si"/>
+    <fraction n="0.2610" ref="O"/>
+  </material>
   <documentation level="3">
     #### Carbon fiber
       a level 3 doc
diff --git a/src/HybridCalorimeter_geo.cpp b/src/HybridCalorimeter_geo.cpp
index 20f4a38c..c7d61330 100644
--- a/src/HybridCalorimeter_geo.cpp
+++ b/src/HybridCalorimeter_geo.cpp
@@ -35,27 +35,26 @@ static Ref_t create_detector(Detector& desc, xml::Handle_t handle, SensitiveDete
     DetElement det(detName, detID);
     sens.setType("calorimeter");
 
-    auto glass_material = desc.material("PbGlass");
+    auto glass_material = desc.material("SciGlass");
     auto crystal_material = desc.material("PbWO4");
     auto air_material = desc.material("Air");
 
     double ROut = desc.constantAsDouble("EcalEndcapN_rmax");
     double RIn = desc.constantAsDouble("EcalEndcapN_rmin");
     double SizeZ = desc.constantAsDouble("EcalEndcapN_thickness");
-    double glass_shift_z = desc.constantAsDouble("GlassModule_z0");
-    double Thickness = desc.constantAsDouble("EcalEndcapN_thickness");
+    double thickness = desc.constantAsDouble("EcalEndcapN_thickness");
     double trans_radius = desc.constantAsDouble("EcalEndcapNCrystal_rmax");
-    double Glass_ShiftZ = desc.constantAsDouble("GlassModule_z0");
+    double Glass_z0 = desc.constantAsDouble("GlassModule_z0");
     double Glass_Width = desc.constantAsDouble("GlassModule_width");
-    double Glass_Thickness = desc.constantAsDouble("GlassModule_length");
+    double Glass_thickness = desc.constantAsDouble("GlassModule_length");
     double Glass_Gap = desc.constantAsDouble("GlassModule_wrap");
     double glass_distance = desc.constantAsDouble("GlassModule_distance");
 
     double Crystal_Width = desc.constantAsDouble("CrystalModule_width");
-    double Crystal_Thickness = desc.constantAsDouble("CrystalModule_length");
+    double Crystal_thickness = desc.constantAsDouble("CrystalModule_length");
     double Crystal_Gap = desc.constantAsDouble("CrystalModule_wrap");
     double crystal_distance = desc.constantAsDouble("CrystalModule_distance");
-    double Crystal_shift_z = desc.constantAsDouble("CrystalModule_z0");
+    double Crystal_z0 = desc.constantAsDouble("CrystalModule_z0");
 
     // RIn and ROut will define outer tube embedding the calorimeter
     // centers_rin/out define the maximum radius of module centers
@@ -64,6 +63,9 @@ static Ref_t create_detector(Detector& desc, xml::Handle_t handle, SensitiveDete
     double centers_rin = RIn + hypotenuse + 1*mm;
     double centers_rout = ROut - hypotenuse - 1*mm;
 
+    const double Crystal_offset = -0.5*(Crystal_thickness - thickness);
+    const double Glass_offset = -0.5*(Glass_thickness - thickness);
+
     // envelope
 
     // Assembly assembly(detName);
@@ -74,15 +76,14 @@ static Ref_t create_detector(Detector& desc, xml::Handle_t handle, SensitiveDete
 
     double Glass_OuterR = ROut - 1 * cm ;
     double Glass_InnerR = trans_radius;
-    glass_shift_z = Thickness / 2. - Glass_Thickness / 2.;
 
     // Geometry of modules
-    Box glass_box("glass_box", Glass_Width * 0.5, Glass_Width * 0.5, Glass_Thickness * 0.5);
+    Box glass_box("glass_box", Glass_Width * 0.5, Glass_Width * 0.5, Glass_thickness * 0.5);
     Volume glass_module("glass_module", glass_box, glass_material);
     glass_module.setVisAttributes(desc.visAttributes("EcalEndcapNModuleVis"));
     glass_module.setSensitiveDetector(sens);
     
-    Box crystal_box("crystal_box",  Crystal_Width* 0.5, Crystal_Width * 0.5, Crystal_Thickness * 0.5);
+    Box crystal_box("crystal_box",  Crystal_Width* 0.5, Crystal_Width * 0.5, Crystal_thickness * 0.5);
     Volume crystal_module("crystal_module", crystal_box, crystal_material);
     crystal_module.setVisAttributes(desc.visAttributes("EcalEndcapNModuleVis"));
     crystal_module.setSensitiveDetector(sens);
@@ -114,7 +115,7 @@ static Ref_t create_detector(Detector& desc, xml::Handle_t handle, SensitiveDete
     int moduleIndex = 0;
 
 //    fmt::print("\nCE EMCAL GLASS SQUARE START\n");
-//    fmt::print("Glass_Thickness = {} cm;\n", Glass_Thickness / cm);
+//    fmt::print("Glass_thickness = {} cm;\n", Glass_thickness / cm);
 //    fmt::print("Glass_Width     = {} cm;\n", Glass_Width / cm);
 //    fmt::print("Glass_Gap       = {} cm;\n", Glass_Gap / cm);
 //    fmt::print("Glass_InnerR    = {} cm;\n", Glass_InnerR / cm);
@@ -180,35 +181,35 @@ static Ref_t create_detector(Detector& desc, xml::Handle_t handle, SensitiveDete
             // first crystal module
             double crystal_x = glass_x - crystal_distance / 2;
             double crystal_y = glass_y - crystal_distance / 2;
-            auto   placement = ecal_vol.placeVolume(crystal_module, Position(crystal_x, crystal_y, Crystal_shift_z));
+            auto   placement = ecal_vol.placeVolume(crystal_module, Position(crystal_x, crystal_y, Crystal_z0 + Crystal_offset));
             placement.addPhysVolID("sector", 1);
             placement.addPhysVolID("module", cryst_module_index++);
 
             // second crystal module
             crystal_x = glass_x + crystal_distance / 2;
             crystal_y = glass_y - crystal_distance / 2;
-            placement = ecal_vol.placeVolume(crystal_module, Position(crystal_x, crystal_y, Crystal_shift_z));
+            placement = ecal_vol.placeVolume(crystal_module, Position(crystal_x, crystal_y, Crystal_z0 + Crystal_offset));
             placement.addPhysVolID("sector", 1);
             placement.addPhysVolID("module", cryst_module_index++);
 
             // third crystal module
             crystal_x = glass_x - crystal_distance / 2;
             crystal_y = glass_y + crystal_distance / 2;
-            placement = ecal_vol.placeVolume(crystal_module, Position(crystal_x, crystal_y, Crystal_shift_z));
+            placement = ecal_vol.placeVolume(crystal_module, Position(crystal_x, crystal_y, Crystal_z0 + Crystal_offset));
             placement.addPhysVolID("sector", 1);
             placement.addPhysVolID("module", cryst_module_index++);
 
             // forth crystal module
             crystal_x = glass_x + crystal_distance / 2;
             crystal_y = glass_y + crystal_distance / 2;
-            placement = ecal_vol.placeVolume(crystal_module, Position(crystal_x, crystal_y, Crystal_shift_z));
+            placement = ecal_vol.placeVolume(crystal_module, Position(crystal_x, crystal_y, Crystal_z0 + Crystal_offset));
             placement.addPhysVolID("sector", 1);
             placement.addPhysVolID("module", cryst_module_index++);
           }
           else
           {
             // glass module
-            auto placement = ecal_vol.placeVolume(glass_module, Position(glass_x, glass_y, glass_shift_z));
+            auto placement = ecal_vol.placeVolume(glass_module, Position(glass_x, glass_y, Glass_z0 + Glass_offset));
             placement.addPhysVolID("sector", 2);
             placement.addPhysVolID("module", glass_module_index++);
           }
diff --git a/src/ScFiCalorimeter_geo.cpp b/src/ScFiCalorimeter_geo.cpp
index 5fc14cef..5b170a00 100644
--- a/src/ScFiCalorimeter_geo.cpp
+++ b/src/ScFiCalorimeter_geo.cpp
@@ -113,58 +113,60 @@ std::tuple<Volume, Position> build_module(const Detector &desc, const xml::Compo
         modVol.setVisAttributes(desc.visAttributes(mod_x.attr<std::string>(_Unicode(vis))));
     }
 
-    auto fiber_x  = mod_x.child(_Unicode(fiber));
-    auto fr       = fiber_x.attr<double>(_Unicode(radius));
-    auto fsx      = fiber_x.attr<double>(_Unicode(spacex));
-    auto fsy      = fiber_x.attr<double>(_Unicode(spacey));
-    auto foff     = dd4hep::getAttrOrDefault<double>(fiber_x, _Unicode(offset), 0.5*mm);
-    auto fiberMat = desc.material(fiber_x.attr<std::string>(_Unicode(material)));
-    Tube fiberShape(0., fr, sz/2.);
-    Volume fiberVol("fiber_vol", fiberShape, fiberMat);
-    fiberVol.setSensitiveDetector(sens);
-
-    // Fibers are placed in a honeycomb with the radius = sqrt(3)/2. * hexagon side length
-    // So each fiber is fully contained in a regular hexagon, which are placed as
-    //           ______________________________________
-    //           |          ____        ____          |
-    // even:     |         /    \      /    \         |
-    //           |    ____/      \____/      \____    |
-    //           |   /    \      /    \      /    \   |
-    // odd:      |  /      \____/      \____/      \  |
-    //           |  \      /    \      /    \      /  |
-    //           |   \____/      \____/      \____/   |
-    // even:     |        \      /    \      /        |
-    //           |         \____/      \____/      ___|___
-    //           |____________________________________|___offset
-    //                                              | |
-    //                                              |offset
-    // the parameters space x and space y are used to add additional spaces between the hexagons
-    double fside  = 2. / std::sqrt(3.) * fr;
-    double fdistx = 2. * fside + fsx;
-    double fdisty = 2. * fr + fsy;
-
-    // maximum numbers of the fibers, help narrow the loop range
-    int nx = int(sx / (2.*fr)) + 1;
-    int ny = int(sy / (2.*fr)) + 1;
-
-    // std::cout << sx << ", " << sy << ", " << fr << ", " << nx << ", " << ny << std::endl;
-
-    // place the fibers
-    double y0 = (foff + fside);
-    int nfibers = 0;
-    for (int iy = 0; iy < ny; ++iy) {
-        double y = y0 + fdisty * iy;
-        // about to touch the boundary
-        if ((sy - y) < y0) { break; }
-        double x0 = (iy % 2) ? (foff + fside) : (foff + fside + fdistx / 2.);
-        for (int ix = 0; ix < nx; ++ix) {
-            double x = x0 + fdistx * ix;
-            // about to touch the boundary
-            if ((sx - x) < x0) { break; }
-            auto fiberPV = modVol.placeVolume(fiberVol, nfibers++, Position{x - sx/2., y - sy/2., 0});
-            std::cout << "(" << ix << ", " << iy << ", " << x - sx/2. << ", " << y - sy/2. << ", " << fr << "),\n";
-            fiberPV.addPhysVolID("fiber_x", ix + 1).addPhysVolID("fiber_y", iy + 1);
-        }
+    if (mod_x.hasChild("fiber")) {
+      auto fiber_x  = mod_x.child(_Unicode(fiber));
+      auto fr       = fiber_x.attr<double>(_Unicode(radius));
+      auto fsx      = fiber_x.attr<double>(_Unicode(spacex));
+      auto fsy      = fiber_x.attr<double>(_Unicode(spacey));
+      auto foff     = dd4hep::getAttrOrDefault<double>(fiber_x, _Unicode(offset), 0.5*mm);
+      auto fiberMat = desc.material(fiber_x.attr<std::string>(_Unicode(material)));
+      Tube fiberShape(0., fr, sz/2.);
+      Volume fiberVol("fiber_vol", fiberShape, fiberMat);
+      fiberVol.setSensitiveDetector(sens);
+
+      // Fibers are placed in a honeycomb with the radius = sqrt(3)/2. * hexagon side length
+      // So each fiber is fully contained in a regular hexagon, which are placed as
+      //           ______________________________________
+      //           |          ____        ____          |
+      // even:     |         /    \      /    \         |
+      //           |    ____/      \____/      \____    |
+      //           |   /    \      /    \      /    \   |
+      // odd:      |  /      \____/      \____/      \  |
+      //           |  \      /    \      /    \      /  |
+      //           |   \____/      \____/      \____/   |
+      // even:     |        \      /    \      /        |
+      //           |         \____/      \____/      ___|___
+      //           |____________________________________|___offset
+      //                                              | |
+      //                                              |offset
+      // the parameters space x and space y are used to add additional spaces between the hexagons
+      double fside  = 2. / std::sqrt(3.) * fr;
+      double fdistx = 2. * fside + fsx;
+      double fdisty = 2. * fr + fsy;
+
+      // maximum numbers of the fibers, help narrow the loop range
+      int nx = int(sx / (2.*fr)) + 1;
+      int ny = int(sy / (2.*fr)) + 1;
+
+      // std::cout << sx << ", " << sy << ", " << fr << ", " << nx << ", " << ny << std::endl;
+
+      // place the fibers
+      double y0 = (foff + fside);
+      int nfibers = 0;
+      for (int iy = 0; iy < ny; ++iy) {
+          double y = y0 + fdisty * iy;
+          // about to touch the boundary
+          if ((sy - y) < y0) { break; }
+          double x0 = (iy % 2) ? (foff + fside) : (foff + fside + fdistx / 2.);
+          for (int ix = 0; ix < nx; ++ix) {
+              double x = x0 + fdistx * ix;
+              // about to touch the boundary
+              if ((sx - x) < x0) { break; }
+              auto fiberPV = modVol.placeVolume(fiberVol, nfibers++, Position{x - sx/2., y - sy/2., 0});
+              //std::cout << "(" << ix << ", " << iy << ", " << x - sx/2. << ", " << y - sy/2. << ", " << fr << "),\n";
+              fiberPV.addPhysVolID("fiber_x", ix + 1).addPhysVolID("fiber_y", iy + 1);
+          }
+      }
     }
 
     return std::make_tuple(modVol, Position{sx, sy, sz});
-- 
GitLab