From dd54011b4a072dc386e77393aac1e2b4bffdc1e1 Mon Sep 17 00:00:00 2001
From: Chao Peng <cpeng@anl.gov>
Date: Sat, 18 Feb 2023 18:45:02 -0600
Subject: [PATCH 01/17] use grid instead of line as an intermediate layer of
 Volumes

---
 src/BarrelCalorimeterInterlayers_geo.cpp | 187 ++++++++++++-----------
 1 file changed, 101 insertions(+), 86 deletions(-)

diff --git a/src/BarrelCalorimeterInterlayers_geo.cpp b/src/BarrelCalorimeterInterlayers_geo.cpp
index 812fda6d5..bfe520703 100644
--- a/src/BarrelCalorimeterInterlayers_geo.cpp
+++ b/src/BarrelCalorimeterInterlayers_geo.cpp
@@ -18,6 +18,8 @@
 #include "Math/Point2D.h"
 #include "TGeoPolygon.h"
 #include "XML/Layering.h"
+#include <functional>
+#include <boost/range/adaptor/filtered.hpp>
 
 using namespace std;
 using namespace dd4hep;
@@ -25,10 +27,32 @@ using namespace dd4hep::detail;
 
 typedef ROOT::Math::XYPoint Point;
 // fiber placement helpers, defined below
-vector<vector<Point>> fiberPositions(double radius, double x_spacing, double z_spacing, double x, double z, double phi,
-                                     double spacing_tol = 1e-2);
+struct FiberGrid {
+  int ix = 0, iz = 0;
+  vector<Point> points;
+  Point mean_centroid = Point(0., 0.);
+  Assembly *assembly_ptr = nullptr;
+
+  // initialize with grid id and points
+  FiberGrid(int i, int j, const vector<Point> &pts) : ix(i), iz(j), points(pts) {
+    if (pts.empty()) {
+      return;
+    }
+
+    double mx = 0., my = 0.;
+    for (auto &p : pts) {
+      mx += p.x();
+      my += p.y();
+    }
+    mx /= static_cast<double>(pts.size());
+    my /= static_cast<double>(pts.size());
+    mean_centroid = Point(mx, my);
+  };
+};
+
+vector<Point> fiberPositions(double r, double sx, double sz, double trx, double trz, double phi, double stol = 1e-2);
 std::pair<int, int>   getNdivisions(double x, double z, double dx, double dz);
-vector<tuple<int, Point, Point, Point, Point>> gridPoints(int div_x, int div_z, double x, double z, double phi);
+vector<FiberGrid> gridPoints(int div_x, int div_z, double x, double z, double phi);
 
 // geometry helpers
 void buildFibers(Detector& desc, SensitiveDetector& sens, Volume& mother, xml_comp_t x_fiber,
@@ -192,71 +216,65 @@ void buildFibers(Detector& desc, SensitiveDetector& sens, Volume& s_vol, xml_com
   // SiPM chip (for GlueX 13mmx13mm: 4x4 grid 3mmx3mm with 3600 50×50 μm pixels each)
   // See, e.g., https://arxiv.org/abs/1801.03088 Fig. 2d
 
-  // Calculate number of divisions
-  auto grid_div = getNdivisions(s_trd_x1, s_thick, 2.0 * cm, 2.0 * cm);
-  // Calculate polygonal grid coordinates (vertices)
-  auto   grid_vtx = gridPoints(grid_div.first, grid_div.second, s_trd_x1, s_thick, hphi);
-  double f_radius_core = f_radius-f_cladding_thickness;
+  // fiber and its cladding
+  double f_radius_core = f_radius - f_cladding_thickness;
   Tube   f_tube_clad(f_radius_core, f_radius, s_length);
   Volume f_vol_clad("fiber_vol", f_tube_clad, desc.material(x_fiber.materialStr()));
   Tube   f_tube_core(0, f_radius_core, s_length);
   Volume f_vol_core("fiber_core_vol", f_tube_core, desc.material(x_fiber.materialStr()));
+  if (x_fiber.isSensitive()) {
+    f_vol_core.setSensitiveDetector(sens);
+  }
+  f_vol_core.setAttributes(desc, x_fiber.regionStr(), x_fiber.limitsStr(), x_fiber.visStr());
 
+
+  // Calculate number of divisions
+  auto        grid_div = getNdivisions(s_trd_x1, s_thick, 2.0 * cm, 2.0 * cm);
+  // Calculate polygonal grid coordinates (vertices)
+  auto        grids = gridPoints(grid_div.first, grid_div.second, s_trd_x1, s_thick, hphi);
   vector<int> f_id_count(grid_div.first * grid_div.second, 0);
   auto        f_pos = fiberPositions(f_radius, f_spacing_x, f_spacing_z, s_trd_x1, s_thick, hphi);
-  // std::cout << f_pos.size() << " lines, ~" << f_pos.front().size() << " fibers each line" << std::endl;
-  for (size_t il = 0; il < f_pos.size(); ++il) {
-    auto& line = f_pos[il];
-    if (line.empty()) {
-      continue;
-    }
-    double l_pos_y = line.front().y();
-    // use assembly as intermediate volume container to reduce number of daughter volumes
-    Assembly lfibers_clad(Form("fiber_clad_array_line_%lu", il));
-    Assembly lfibers_core(Form("fiber_core_array_line_%lu", il));
-    for (auto& p : line) {
-      int f_grid_id = -1;
-      int f_id      = -1;
-      // Check to which grid fiber belongs to
-      for (auto& poly_vtx : grid_vtx) {
-        if (p.y() != l_pos_y) {
-          std::cerr << Form("Expected the same y position from a same line: %.2f, but got %.2f", l_pos_y, p.y())
-                    << std::endl;
-          continue;
-        }
-        auto [grid_id, vtx_a, vtx_b, vtx_c, vtx_d] = poly_vtx;
-        double poly_x[4]                           = {vtx_a.x(), vtx_b.x(), vtx_c.x(), vtx_d.x()};
-        double poly_y[4]                           = {vtx_a.y(), vtx_b.y(), vtx_c.y(), vtx_d.y()};
-        double f_xy[2]                             = {p.x(), p.y()};
-
-        TGeoPolygon poly(4);
-        poly.SetXY(poly_x, poly_y);
-        poly.FinishPolygon();
-
-        if (poly.Contains(f_xy)) {
-          f_grid_id = grid_id;
-          f_id      = f_id_count[grid_id];
-          f_id_count[grid_id]++;
-        }
+  // a helper struct to speed up searching
+  struct Fiber {
+    Point pos;
+    bool assigned = false;
+    Fiber (const Point &p) : pos(p) {};
+  };
+  std::vector<Fiber> fibers(f_pos.begin(), f_pos.end());
+
+  // build assembly for each grid and put fibers in
+  for (auto &gr : grids) {
+    Assembly grid_vol(Form("fiber_grid_%i_%i", gr.ix, gr.iz));
+    // fiber is along y-axis of the layer volume, so grids are arranged on X-Z plane
+    Transform3D gr_tr(RotationZYX(0, 0, M_PI * 0.5), Position(gr.mean_centroid.x(), 0, gr.mean_centroid.y()));
+    auto grid_phv = s_vol.placeVolume(grid_vol, gr_tr);
+    grid_phv.addPhysVolID(f_id_grid, gr.ix + gr.iz * grid_div.first + 1);
+
+    // loop over all fibers that are not assigned to a grid
+    int f_id = 1;
+    for (auto &fi : fibers | boost::adaptors::filtered( [](const Fiber &f) {return not f.assigned;} )) {
+      // use TGeoPolygon to help check if fiber is inside a grid
+      TGeoPolygon poly(gr.points.size());
+      vector<double> vx, vy;
+      transform(gr.points.begin(), gr.points.end(), back_inserter(vx), mem_fn(&Point::x));
+      transform(gr.points.begin(), gr.points.end(), back_inserter(vy), mem_fn(&Point::y));
+      poly.SetXY(vx.data(), vy.data());
+      poly.FinishPolygon();
+
+      double f_xy[2] = {fi.pos.x(), fi.pos.y()};
+      if (not poly.Contains(f_xy)) {
+        continue;
       }
 
-      if (x_fiber.isSensitive()) {
-        f_vol_core.setSensitiveDetector(sens);
-      }
-      f_vol_core.setAttributes(desc, x_fiber.regionStr(), x_fiber.limitsStr(), x_fiber.visStr());
-
-      // Fiber placement
-      // Transform3D f_tr(RotationZYX(0,0,M_PI*0.5),Position(p.x(), 0, p.y()));
-      // PlacedVolume fiber_phv = s_vol.placeVolume(f_vol, Position(p.x(), 0., p.y()));
-      PlacedVolume core_phv = lfibers_core.placeVolume(f_vol_core, Position(p.x(), 0., 0.));
-      core_phv.addPhysVolID(f_id_grid, f_grid_id + 1).addPhysVolID(f_id_fiber, f_id + 1);
-      lfibers_clad.placeVolume(f_vol_clad, Position(p.x(), 0., 0.));
+      // place fiber in grid
+      auto p = fi.pos - gr.mean_centroid;
+      auto core_phv = grid_vol.placeVolume(f_vol_core, Position(p.x(), 0., p.y()));
+      core_phv.addPhysVolID(f_id_fiber, f_id);
+      grid_vol.placeVolume(f_vol_clad, Position(p.x(), 0., p.y()));
+      fi.assigned = true;
+      f_id ++;
     }
-    lfibers_core.ptr()->Voxelize("");
-    lfibers_clad.ptr()->Voxelize("");
-    Transform3D l_tr(RotationZYX(0, 0, M_PI * 0.5), Position(0., 0, l_pos_y));
-    s_vol.placeVolume(lfibers_core, l_tr);
-    s_vol.placeVolume(lfibers_clad, l_tr);
+    grid_vol.ptr()->Voxelize("");
   }
 }
 
@@ -350,36 +368,36 @@ void buildSupport(Detector& desc, Volume& mod_vol, xml_comp_t x_support,
 }
 
 // Fill fiber lattice into trapezoid starting from position (0,0) in x-z coordinate system
-vector<vector<Point>> fiberPositions(double radius, double x_spacing, double z_spacing, double x, double z, double phi,
-                                     double spacing_tol)
+vector<Point> fiberPositions(double r, double sx, double sz, double trx, double trz, double phi, double stol)
 {
-  // z_spacing - distance between fiber layers in z
-  // x_spacing - distance between fiber centers in x
-  // x - half-length of the shorter (bottom) base of the trapezoid
-  // z - height of the trapezoid
-  // phi - angle between z and trapezoid arm
+  // r      - fiber radius
+  // sx, sz - spacing between fibers in x, z
+  // trx    - half-length of the shorter (bottom) base of the trapezoid
+  // trz    - height of the trapezoid
+  // phi    - angle between z and trapezoid arm
+  // stol   - spacing tolerance
 
-  vector<vector<Point>> positions;
-  int z_layers = floor((z / 2 - radius - spacing_tol) / z_spacing); // number of layers that fit in z/2
+  vector<Point> positions;
+  int z_layers = floor((trz / 2 - r - stol) / sz); // number of layers that fits in half trapezoid-z
 
-  double z_pos = 0.;
-  double x_pos = 0.;
+  double px = 0., pz = 0.;
 
   for (int l = -z_layers; l < z_layers + 1; l++) {
     vector<Point> xline;
-    z_pos        = l * z_spacing;
-    double x_max = x + (z / 2. + z_pos) * tan(phi) - spacing_tol; // calculate max x at particular z_pos
-    (l % 2 == 0) ? x_pos = 0. : x_pos = x_spacing / 2;            // account for spacing/2 shift
-
-    while (x_pos < (x_max - radius)) {
-      xline.push_back(Point(x_pos, z_pos));
-      if (x_pos != 0.)
-        xline.push_back(Point(-x_pos, z_pos)); // using symmetry around x=0
-      x_pos += x_spacing;
+    pz           = l * sz;
+    double x_max = trx + (trz / 2. + pz) * tan(phi) - stol; // calculate max x at particular z_pos
+    (l % 2 == 0) ? px = 0. : px = sx / 2;                   // account for spacing/2 shift
+
+    while (px < (x_max - r)) {
+      xline.push_back(Point(px, pz));
+      if (px != 0.)
+        xline.push_back(Point(-px, pz)); // using symmetry around x=0
+      px += sx;
     }
+
     // Sort fiber IDs for a better organization
     sort(xline.begin(), xline.end(), [](const Point& p1, const Point& p2) { return p1.x() < p2.x(); });
-    positions.emplace_back(std::move(xline));
+    positions.insert(positions.end(), xline.begin(), xline.end());
   }
   return positions;
 }
@@ -418,14 +436,13 @@ std::pair<int, int> getNdivisions(double x, double z, double dx, double dz)
 }
 
 // Calculate dimensions of the polygonal grid in the cartesian coordinate system x-z
-vector<tuple<int, Point, Point, Point, Point>> gridPoints(int div_x, int div_z, double x, double z, double phi)
+vector<FiberGrid> gridPoints(int div_x, int div_z, double x, double z, double phi)
 {
   // x, z and phi defined as in vector<Point> fiberPositions
   // div_x, div_z - number of divisions in x and z
+  vector<FiberGrid> results;
   double dz = z / div_z;
 
-  std::vector<std::tuple<int, Point, Point, Point, Point>> points;
-
   for (int iz = 0; iz < div_z + 1; iz++) {
     for (int ix = 0; ix < div_x + 1; ix++) {
       double A_z = -z / 2 + iz * dz;
@@ -445,19 +462,17 @@ vector<tuple<int, Point, Point, Point, Point>> gridPoints(int div_x, int div_z,
       double C_x = B_x + dx_for_z_plus_1;
       double D_x = A_x + dx_for_z;
 
-      int id = ix + div_x * iz;
-
       auto A = Point(A_x, A_z);
       auto B = Point(B_x, B_z);
       auto C = Point(C_x, C_z);
       auto D = Point(D_x, D_z);
 
       // vertex points filled in the clock-wise direction
-      points.push_back(make_tuple(id, A, B, C, D));
+      results.emplace_back(FiberGrid(ix, iz, {A, B, C, D}));
     }
   }
 
-  return points;
+  return results;
 }
 
 DECLARE_DETELEMENT(epic_EcalBarrelInterlayers, create_detector)

From 006c43a2ea6cae1925ee60ffd2c25fe084c16547 Mon Sep 17 00:00:00 2001
From: Chao Peng <cpeng@anl.gov>
Date: Sun, 19 Feb 2023 00:13:34 -0600
Subject: [PATCH 02/17] remove c++20 feature

---
 src/BarrelCalorimeterInterlayers_geo.cpp | 16 +++++++++-------
 1 file changed, 9 insertions(+), 7 deletions(-)

diff --git a/src/BarrelCalorimeterInterlayers_geo.cpp b/src/BarrelCalorimeterInterlayers_geo.cpp
index bfe520703..633fd4e54 100644
--- a/src/BarrelCalorimeterInterlayers_geo.cpp
+++ b/src/BarrelCalorimeterInterlayers_geo.cpp
@@ -19,22 +19,20 @@
 #include "TGeoPolygon.h"
 #include "XML/Layering.h"
 #include <functional>
-#include <boost/range/adaptor/filtered.hpp>
 
 using namespace std;
 using namespace dd4hep;
-using namespace dd4hep::detail;
 
 typedef ROOT::Math::XYPoint Point;
 // fiber placement helpers, defined below
 struct FiberGrid {
-  int ix = 0, iz = 0;
+  int ix = 0, iy = 0;
   vector<Point> points;
   Point mean_centroid = Point(0., 0.);
   Assembly *assembly_ptr = nullptr;
 
   // initialize with grid id and points
-  FiberGrid(int i, int j, const vector<Point> &pts) : ix(i), iz(j), points(pts) {
+  FiberGrid(int i, int j, const vector<Point> &pts) : ix(i), iy(j), points(pts) {
     if (pts.empty()) {
       return;
     }
@@ -244,15 +242,19 @@ void buildFibers(Detector& desc, SensitiveDetector& sens, Volume& s_vol, xml_com
 
   // build assembly for each grid and put fibers in
   for (auto &gr : grids) {
-    Assembly grid_vol(Form("fiber_grid_%i_%i", gr.ix, gr.iz));
+    Assembly grid_vol(Form("fiber_grid_%i_%i", gr.ix, gr.iy));
     // fiber is along y-axis of the layer volume, so grids are arranged on X-Z plane
     Transform3D gr_tr(RotationZYX(0, 0, M_PI * 0.5), Position(gr.mean_centroid.x(), 0, gr.mean_centroid.y()));
     auto grid_phv = s_vol.placeVolume(grid_vol, gr_tr);
-    grid_phv.addPhysVolID(f_id_grid, gr.ix + gr.iz * grid_div.first + 1);
+    grid_phv.addPhysVolID(f_id_grid, gr.ix + gr.iy * grid_div.first + 1);
 
     // loop over all fibers that are not assigned to a grid
     int f_id = 1;
-    for (auto &fi : fibers | boost::adaptors::filtered( [](const Fiber &f) {return not f.assigned;} )) {
+    for (auto &fi : fibers) {
+      if (fi.assigned) {
+        continue;
+      }
+
       // use TGeoPolygon to help check if fiber is inside a grid
       TGeoPolygon poly(gr.points.size());
       vector<double> vx, vy;

From 03aa81ed08888daa40305935c2c822d633308750 Mon Sep 17 00:00:00 2001
From: Chao Peng <cpeng@anl.gov>
Date: Sun, 19 Feb 2023 10:27:10 -0600
Subject: [PATCH 03/17] initial commit for scfi grid visualization script

---
 .../subdetector_tests/draw_bemc_scfi_grids.py | 46 +++++++++++++++++++
 1 file changed, 46 insertions(+)
 create mode 100644 scripts/subdetector_tests/draw_bemc_scfi_grids.py

diff --git a/scripts/subdetector_tests/draw_bemc_scfi_grids.py b/scripts/subdetector_tests/draw_bemc_scfi_grids.py
new file mode 100644
index 000000000..05533d2b1
--- /dev/null
+++ b/scripts/subdetector_tests/draw_bemc_scfi_grids.py
@@ -0,0 +1,46 @@
+'''
+    A script to visualize the fibers and of one grid from BEMC ScFi part
+    Chao Peng (ANL)
+'''
+import os
+import DDG4
+import argparse
+from matplotlib import pyplot as plt
+import matplotlib.ticker as ticker
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+            '-c', '--compact',
+            dest='compact', required=True,
+            help='Top-level xml file of the detector description'
+            )
+    parser.add_argument(
+            '-o',
+            dest='outdir', default='.',
+            help='Output directory.'
+            )
+    parser.add_argument(
+            '--layer', type=int, default=1,
+            help='The BEMC ScFi layer number.'
+            )
+    parser.add_argument(
+            '--center-grid', dest='center_grid', type=int, default=1,
+            help='The grid number for centering in the plot.'
+            )
+    args = parser.parse_args()
+
+    # initialize dd4hep detector
+    kernel = DDG4.Kernel()
+    description = kernel.detectorDescription()
+    kernel.loadGeometry("file:{}".format(args.compact))
+
+    # dd4hep_decoder = description.readout(args.readout).idSpec().decoder()
+    # lindex = dd4hep_decoder.index('x')
+    # get_layer_id = np.vectorize(lambda cid: dd4hep_decoder.get(cid, lindex))
+    # df.loc[:, 'layerID'] = get_layer_id(df['cellID'].astype(int).values)
+
+    # always terminate dd4hep kernel
+    kernel.terminate()
+

From 719010b3b633c2426ff2204da46ef37207408b0a Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Sun, 19 Feb 2023 16:27:23 +0000
Subject: [PATCH 04/17] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 scripts/subdetector_tests/draw_bemc_scfi_grids.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/scripts/subdetector_tests/draw_bemc_scfi_grids.py b/scripts/subdetector_tests/draw_bemc_scfi_grids.py
index 05533d2b1..8afcbd53d 100644
--- a/scripts/subdetector_tests/draw_bemc_scfi_grids.py
+++ b/scripts/subdetector_tests/draw_bemc_scfi_grids.py
@@ -43,4 +43,3 @@
 
     # always terminate dd4hep kernel
     kernel.terminate()
-

From 3db523a3d8ea45852179334e82b705ca07c3513b Mon Sep 17 00:00:00 2001
From: Chao Peng <cpeng@anl.gov>
Date: Sun, 19 Feb 2023 17:31:58 -0600
Subject: [PATCH 05/17] update the script, now it can get geometry information
 of a fiber/grid

---
 .../subdetector_tests/draw_bemc_scfi_grids.py | 79 ++++++++++++++++---
 1 file changed, 68 insertions(+), 11 deletions(-)

diff --git a/scripts/subdetector_tests/draw_bemc_scfi_grids.py b/scripts/subdetector_tests/draw_bemc_scfi_grids.py
index 8afcbd53d..422a588a3 100644
--- a/scripts/subdetector_tests/draw_bemc_scfi_grids.py
+++ b/scripts/subdetector_tests/draw_bemc_scfi_grids.py
@@ -3,10 +3,22 @@
     Chao Peng (ANL)
 '''
 import os
-import DDG4
+import ROOT
+import dd4hep
+import DDRec
 import argparse
 from matplotlib import pyplot as plt
 import matplotlib.ticker as ticker
+from collections import OrderedDict
+
+
+def dict_to_cpp_vec(my_dict, dtype='int'):
+    # use ROOT to make sure the type is correct
+    vol_ids = ROOT.std.vector('pair<string, {}>'.format(dtype))()
+    for field, fid in my_dict.items():
+        vol_ids.push_back((field, fid))
+    return vol_ids
+
 
 
 if __name__ == '__main__':
@@ -14,16 +26,22 @@
     parser.add_argument(
             '-c', '--compact',
             dest='compact', required=True,
-            help='Top-level xml file of the detector description'
+            help='Top-level xml file of the detector description.'
             )
     parser.add_argument(
-            '-o',
+            '--dpath', default='EcalBarrelScFi/stave0/layer1/slice1',
+            dest='dpath',
+            help='Path to the scfi layer that contains the grids to draw.'
+            )
+    parser.add_argument(
+            '--outdir',
             dest='outdir', default='.',
             help='Output directory.'
             )
     parser.add_argument(
-            '--layer', type=int, default=1,
-            help='The BEMC ScFi layer number.'
+            '--readout',
+            default='EcalBarrelScFiHits',
+            help='Readout class for the EcalBarrelScFi.'
             )
     parser.add_argument(
             '--center-grid', dest='center_grid', type=int, default=1,
@@ -32,14 +50,53 @@
     args = parser.parse_args()
 
     # initialize dd4hep detector
-    kernel = DDG4.Kernel()
-    description = kernel.detectorDescription()
-    kernel.loadGeometry("file:{}".format(args.compact))
+    desc = dd4hep.Detector.getInstance()
+    desc.fromXML(args.compact)
+    converter = DDRec.CellIDPositionConverter(desc)
+
+    # search the DetElement start from the world
+    det_names = args.dpath.split('/')
+    det = desc.world()
+
+    for i, n in enumerate(det_names):
+        children = det.children()
+        try:
+            det = det.child(n)
+        except Exception:
+            print('Failed to find DetElement {} from {}'.format(n, '/'.join(['world'] + det_names[:i])))
+            print('Available children are listed below:')
+            for n, d in children:
+                print(' --- child: {}'.format(n))
+            exit(-1)
+
+    # build a volume manager so it triggers populating the volume IDs
+    vman = dd4hep.VolumeManager(det, desc.readout(args.readout))
+    idspec = desc.readout(args.readout).idSpec()
+
+    # build a dictionary for the readout ids
+    mvol_id = det.volumeID()
+    id_dict = OrderedDict()
+    for idstr in idspec.decoder().valueString(mvol_id).split(','):
+        field, fid = idstr.split(':')
+        id_dict[field] = int(fid)
+
+    # update
+    id_dict.update({'grid': 1})
+    vid = idspec.encode(dict_to_cpp_vec(id_dict))
+    print(idspec.decoder().valueString(vid))
+    vol = converter.findContext(vid).volumePlacement().volume()
+    print(vol.GetNdaughters())
+
+
+    # [(fid[0], int(fid[1])) for fid in idspec.decoder().valueString(mvol_id).split(',')]
+    # use ROOT to make sure the type is correct
+    # vol_ids = ROOT.std.vector('pair<string, int>')()
+    #     vol_ids.push_back((field, int(fid)))
+    # vol_ids[-2].second = 2
+    # print(vol_ids)
+    # print(idspec.encode(vol_ids))
 
     # dd4hep_decoder = description.readout(args.readout).idSpec().decoder()
     # lindex = dd4hep_decoder.index('x')
     # get_layer_id = np.vectorize(lambda cid: dd4hep_decoder.get(cid, lindex))
     # df.loc[:, 'layerID'] = get_layer_id(df['cellID'].astype(int).values)
-
-    # always terminate dd4hep kernel
-    kernel.terminate()

From 670fd5b25136942037ae20dcecd03ab9a7d2dc17 Mon Sep 17 00:00:00 2001
From: Chao Peng <cpeng@anl.gov>
Date: Sun, 19 Feb 2023 21:19:32 -0600
Subject: [PATCH 06/17] fix a typo in transform matrix

---
 src/BarrelCalorimeterInterlayers_geo.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/BarrelCalorimeterInterlayers_geo.cpp b/src/BarrelCalorimeterInterlayers_geo.cpp
index 633fd4e54..de8d49275 100644
--- a/src/BarrelCalorimeterInterlayers_geo.cpp
+++ b/src/BarrelCalorimeterInterlayers_geo.cpp
@@ -244,7 +244,7 @@ void buildFibers(Detector& desc, SensitiveDetector& sens, Volume& s_vol, xml_com
   for (auto &gr : grids) {
     Assembly grid_vol(Form("fiber_grid_%i_%i", gr.ix, gr.iy));
     // fiber is along y-axis of the layer volume, so grids are arranged on X-Z plane
-    Transform3D gr_tr(RotationZYX(0, 0, M_PI * 0.5), Position(gr.mean_centroid.x(), 0, gr.mean_centroid.y()));
+    Transform3D gr_tr(RotationZYX(0, M_PI*0.5, 0.), Position(gr.mean_centroid.x(), 0, gr.mean_centroid.y()));
     auto grid_phv = s_vol.placeVolume(grid_vol, gr_tr);
     grid_phv.addPhysVolID(f_id_grid, gr.ix + gr.iy * grid_div.first + 1);
 

From 93ac434b2c1fe3f246921585b74f0b42a53ab1a2 Mon Sep 17 00:00:00 2001
From: Chao Peng <cpeng@anl.gov>
Date: Sun, 19 Feb 2023 21:44:06 -0600
Subject: [PATCH 07/17] fix the rotation and position for fibers

---
 src/BarrelCalorimeterInterlayers_geo.cpp | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/BarrelCalorimeterInterlayers_geo.cpp b/src/BarrelCalorimeterInterlayers_geo.cpp
index de8d49275..210d88fe3 100644
--- a/src/BarrelCalorimeterInterlayers_geo.cpp
+++ b/src/BarrelCalorimeterInterlayers_geo.cpp
@@ -244,7 +244,7 @@ void buildFibers(Detector& desc, SensitiveDetector& sens, Volume& s_vol, xml_com
   for (auto &gr : grids) {
     Assembly grid_vol(Form("fiber_grid_%i_%i", gr.ix, gr.iy));
     // fiber is along y-axis of the layer volume, so grids are arranged on X-Z plane
-    Transform3D gr_tr(RotationZYX(0, M_PI*0.5, 0.), Position(gr.mean_centroid.x(), 0, gr.mean_centroid.y()));
+    Transform3D gr_tr(RotationZYX(0., 0., M_PI*0.5), Position(gr.mean_centroid.x(), 0., gr.mean_centroid.y()));
     auto grid_phv = s_vol.placeVolume(grid_vol, gr_tr);
     grid_phv.addPhysVolID(f_id_grid, gr.ix + gr.iy * grid_div.first + 1);
 
@@ -270,9 +270,9 @@ void buildFibers(Detector& desc, SensitiveDetector& sens, Volume& s_vol, xml_com
 
       // place fiber in grid
       auto p = fi.pos - gr.mean_centroid;
-      auto core_phv = grid_vol.placeVolume(f_vol_core, Position(p.x(), 0., p.y()));
+      auto core_phv = grid_vol.placeVolume(f_vol_core, Position(p.x(), p.y(), 0.));
       core_phv.addPhysVolID(f_id_fiber, f_id);
-      grid_vol.placeVolume(f_vol_clad, Position(p.x(), 0., p.y()));
+      grid_vol.placeVolume(f_vol_clad, Position(p.x(), p.y(), 0.));
       fi.assigned = true;
       f_id ++;
     }

From a2d33543533d211d63f35bd36f1204cf3501c78a Mon Sep 17 00:00:00 2001
From: Chao Peng <cpeng@anl.gov>
Date: Sun, 19 Feb 2023 21:44:19 -0600
Subject: [PATCH 08/17] finish geometry check script

---
 .../subdetector_tests/draw_bemc_scfi_grids.py | 109 ++++++++++++++----
 1 file changed, 88 insertions(+), 21 deletions(-)

diff --git a/scripts/subdetector_tests/draw_bemc_scfi_grids.py b/scripts/subdetector_tests/draw_bemc_scfi_grids.py
index 422a588a3..aeb0aebcb 100644
--- a/scripts/subdetector_tests/draw_bemc_scfi_grids.py
+++ b/scripts/subdetector_tests/draw_bemc_scfi_grids.py
@@ -7,11 +7,15 @@
 import dd4hep
 import DDRec
 import argparse
+import numpy as np
 from matplotlib import pyplot as plt
+from matplotlib.patches import Circle
+from matplotlib.collections import PatchCollection
 import matplotlib.ticker as ticker
 from collections import OrderedDict
 
 
+# helper function to do some type conversion for python wrapper of C++ function
 def dict_to_cpp_vec(my_dict, dtype='int'):
     # use ROOT to make sure the type is correct
     vol_ids = ROOT.std.vector('pair<string, {}>'.format(dtype))()
@@ -20,6 +24,51 @@ def dict_to_cpp_vec(my_dict, dtype='int'):
     return vol_ids
 
 
+# helper function to collect fibers under a grid
+def get_fibers(detelem, id_desc, id_conv, id_dict, grid=1):
+    gnode = detelem.volume().GetNode(grid-1)
+    global_trans = detelem.nominal().worldTransformation()
+    print(gnode.GetName())
+    # get grid
+    id_dict.update({'grid': grid, 'fiber': 0})
+    gid = id_desc.encode(dict_to_cpp_vec(id_dict))
+    # print(id_desc.decoder().valueString(vid))
+    vol = id_conv.findContext(gid).volumePlacement().volume()
+    grpos = id_conv.position(gid)
+    # print(vol.boundingBox().dimensions())
+
+    # use TGeoNode to get center positions
+    # here it can also use converter to get position but it's much slower (adds an additional lookup process)
+    fibers = []
+    for i in np.arange(vol.GetNdaughters()):
+        fnode = vol.GetNode(int(i))
+        # NOTE, this is defined in geometry plugin, fiber_core is the wanted sensitive detector
+        if 'fiber_core' not in fnode.GetName():
+            continue
+        fpos = np.array([0., 0., 0.])
+        gpos = np.array([0., 0., 0.])
+        pos = np.array([0., 0., 0.])
+        fnode.LocalToMaster(np.array([0., 0., 0.]), fpos)
+        gnode.LocalToMaster(fpos, gpos)
+        # the two method below are equivalent
+        global_trans.LocalToMaster(gpos, pos)
+        # detelem.nominal().localToWorld(gpos, pos)
+        """ a test with converter
+        if i < 50:
+            id_dict.update({'fiber': int(len(fibers) + 1)})
+            fid = id_desc.encode(dict_to_cpp_vec(id_dict))
+            tpos = id_conv.position(fid)
+            print(i,
+                  fnode.GetName(),
+                  np.asarray([tpos.X(), tpos.Y(), tpos.Z()]),
+                  pos,
+                  fpos,
+                  gpos)
+        """
+        fibers.append(pos)
+
+    return fibers, grpos
+
 
 if __name__ == '__main__':
     parser = argparse.ArgumentParser()
@@ -31,7 +80,7 @@ def dict_to_cpp_vec(my_dict, dtype='int'):
     parser.add_argument(
             '--dpath', default='EcalBarrelScFi/stave0/layer1/slice1',
             dest='dpath',
-            help='Path to the scfi layer that contains the grids to draw.'
+            help='Path to the closest DetElement (slice) that contains the grids and fibers.'
             )
     parser.add_argument(
             '--outdir',
@@ -44,7 +93,7 @@ def dict_to_cpp_vec(my_dict, dtype='int'):
             help='Readout class for the EcalBarrelScFi.'
             )
     parser.add_argument(
-            '--center-grid', dest='center_grid', type=int, default=1,
+            '--center-grid', dest='center_grid', type=int, default=2,
             help='The grid number for centering in the plot.'
             )
     args = parser.parse_args()
@@ -71,7 +120,10 @@ def dict_to_cpp_vec(my_dict, dtype='int'):
 
     # build a volume manager so it triggers populating the volume IDs
     vman = dd4hep.VolumeManager(det, desc.readout(args.readout))
-    idspec = desc.readout(args.readout).idSpec()
+
+    # manager instances
+    readout = desc.readout(args.readout)
+    idspec = readout.idSpec()
 
     # build a dictionary for the readout ids
     mvol_id = det.volumeID()
@@ -80,23 +132,38 @@ def dict_to_cpp_vec(my_dict, dtype='int'):
         field, fid = idstr.split(':')
         id_dict[field] = int(fid)
 
-    # update
-    id_dict.update({'grid': 1})
-    vid = idspec.encode(dict_to_cpp_vec(id_dict))
-    print(idspec.decoder().valueString(vid))
-    vol = converter.findContext(vid).volumePlacement().volume()
-    print(vol.GetNdaughters())
+    # get radius of the fiber
+    id_dict.update({'grid': args.center_grid, 'fiber': 1})
+    fid = idspec.encode(dict_to_cpp_vec(id_dict))
+    # need to do unit conversion from mm to cm (might from ROOT/DD4Hep difference)
+    fr = converter.cellDimensions(fid)[0]/2./10.
 
+    # plot fibers in the grid
+    fig, ax = plt.subplots(figsize=(12, 12), dpi=160)
+    # default color cycles
+    colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
+
+    # plot the center grid and adjacent grid
+    grids = [
+        (g, colors[i]) for i, g in enumerate([args.center_grid, args.center_grid - 1, args.center_grid + 1]) if g > 0
+        ]
+    for g, c in grids:
+        # get position of fibers
+        fibers, grid_pos = get_fibers(det, idspec, converter, id_dict, g)
+        patches = []
+        for fpos in fibers:
+            patches.append(Circle((fpos[0], fpos[1]), fr))
+        p = PatchCollection(patches, alpha=0.4, facecolors=(c,), edgecolors=('k',))
+        ax.plot([grid_pos.X()], [grid_pos.Y()], marker='P', mfc=c, mec='k', ms=9, label='grid {} center'.format(g))
+        ax.add_collection(p)
+        if g == args.center_grid:
+            ax.set_xlim(grid_pos.X()- 3., grid_pos.X() + 3.)
+            ax.set_ylim(grid_pos.Y() - 3., grid_pos.Y() + 3.)
+    ax.legend(fontsize=22)
+    ax.tick_params(labelsize=20, direction='in')
+    ax.set_xlabel('Global X (mm)', fontsize=22)
+    ax.set_ylabel('Global Y (mm)', fontsize=22)
+    ax.set_title('{} around grid {}'.format(args.dpath, args.center_grid), fontsize=22)
+
+    fig.savefig(os.path.join(args.outdir, 'grid_fibers.png'))
 
-    # [(fid[0], int(fid[1])) for fid in idspec.decoder().valueString(mvol_id).split(',')]
-    # use ROOT to make sure the type is correct
-    # vol_ids = ROOT.std.vector('pair<string, int>')()
-    #     vol_ids.push_back((field, int(fid)))
-    # vol_ids[-2].second = 2
-    # print(vol_ids)
-    # print(idspec.encode(vol_ids))
-
-    # dd4hep_decoder = description.readout(args.readout).idSpec().decoder()
-    # lindex = dd4hep_decoder.index('x')
-    # get_layer_id = np.vectorize(lambda cid: dd4hep_decoder.get(cid, lindex))
-    # df.loc[:, 'layerID'] = get_layer_id(df['cellID'].astype(int).values)

From 4441a4644ede12491d031b47f21dd49fb67b5614 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Mon, 20 Feb 2023 03:44:34 +0000
Subject: [PATCH 09/17] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 scripts/subdetector_tests/draw_bemc_scfi_grids.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/scripts/subdetector_tests/draw_bemc_scfi_grids.py b/scripts/subdetector_tests/draw_bemc_scfi_grids.py
index aeb0aebcb..d2707ed29 100644
--- a/scripts/subdetector_tests/draw_bemc_scfi_grids.py
+++ b/scripts/subdetector_tests/draw_bemc_scfi_grids.py
@@ -166,4 +166,3 @@ def get_fibers(detelem, id_desc, id_conv, id_dict, grid=1):
     ax.set_title('{} around grid {}'.format(args.dpath, args.center_grid), fontsize=22)
 
     fig.savefig(os.path.join(args.outdir, 'grid_fibers.png'))
-

From 500674ee7a707bf70db4af856032dda23bac6df2 Mon Sep 17 00:00:00 2001
From: Chao Peng <cpeng@anl.gov>
Date: Sun, 19 Feb 2023 21:52:27 -0600
Subject: [PATCH 10/17] follow up: a few minor changes [skip ci]

---
 scripts/subdetector_tests/draw_bemc_scfi_grids.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/scripts/subdetector_tests/draw_bemc_scfi_grids.py b/scripts/subdetector_tests/draw_bemc_scfi_grids.py
index d2707ed29..1a032c2e8 100644
--- a/scripts/subdetector_tests/draw_bemc_scfi_grids.py
+++ b/scripts/subdetector_tests/draw_bemc_scfi_grids.py
@@ -1,5 +1,6 @@
 '''
-    A script to visualize the fibers and of one grid from BEMC ScFi part
+    A script to visualize the fibers of some grids from BEMC ScFi part
+    02/19/2023
     Chao Peng (ANL)
 '''
 import os
@@ -28,7 +29,7 @@ def dict_to_cpp_vec(my_dict, dtype='int'):
 def get_fibers(detelem, id_desc, id_conv, id_dict, grid=1):
     gnode = detelem.volume().GetNode(grid-1)
     global_trans = detelem.nominal().worldTransformation()
-    print(gnode.GetName())
+    # print(gnode.GetName())
     # get grid
     id_dict.update({'grid': grid, 'fiber': 0})
     gid = id_desc.encode(dict_to_cpp_vec(id_dict))

From f29a3c376bf0f7f620c3912a148e9b87a2e91920 Mon Sep 17 00:00:00 2001
From: Chao Peng <cpeng@anl.gov>
Date: Sun, 19 Feb 2023 23:30:31 -0600
Subject: [PATCH 11/17] do not add empty grid if no actual fibers are in it

---
 src/BarrelCalorimeterInterlayers_geo.cpp | 26 +++++++++++++++++++-----
 1 file changed, 21 insertions(+), 5 deletions(-)

diff --git a/src/BarrelCalorimeterInterlayers_geo.cpp b/src/BarrelCalorimeterInterlayers_geo.cpp
index 210d88fe3..62dc93007 100644
--- a/src/BarrelCalorimeterInterlayers_geo.cpp
+++ b/src/BarrelCalorimeterInterlayers_geo.cpp
@@ -243,10 +243,6 @@ void buildFibers(Detector& desc, SensitiveDetector& sens, Volume& s_vol, xml_com
   // build assembly for each grid and put fibers in
   for (auto &gr : grids) {
     Assembly grid_vol(Form("fiber_grid_%i_%i", gr.ix, gr.iy));
-    // fiber is along y-axis of the layer volume, so grids are arranged on X-Z plane
-    Transform3D gr_tr(RotationZYX(0., 0., M_PI*0.5), Position(gr.mean_centroid.x(), 0., gr.mean_centroid.y()));
-    auto grid_phv = s_vol.placeVolume(grid_vol, gr_tr);
-    grid_phv.addPhysVolID(f_id_grid, gr.ix + gr.iy * grid_div.first + 1);
 
     // loop over all fibers that are not assigned to a grid
     int f_id = 1;
@@ -276,8 +272,28 @@ void buildFibers(Detector& desc, SensitiveDetector& sens, Volume& s_vol, xml_com
       fi.assigned = true;
       f_id ++;
     }
-    grid_vol.ptr()->Voxelize("");
+
+    // only add this if this grid has fibers
+    if (f_id > 1) {
+        // fiber is along y-axis of the layer volume, so grids are arranged on X-Z plane
+        Transform3D gr_tr(RotationZYX(0., 0., M_PI*0.5), Position(gr.mean_centroid.x(), 0., gr.mean_centroid.y()));
+        auto grid_phv = s_vol.placeVolume(grid_vol, gr_tr);
+        grid_phv.addPhysVolID(f_id_grid, gr.ix + gr.iy * grid_div.first + 1);
+        grid_vol.ptr()->Voxelize("");
+    }
+  }
+
+  /*
+  // sanity check
+  size_t missing_fibers = 0;
+  for (auto &fi : fibers) {
+    if (not fi.assigned) {
+      missing_fibers++;
+    }
   }
+  std::cout << "built " << fibers.size() << " fibers, "
+            << missing_fibers << " of them failed to find a grid" << std::endl;
+  */
 }
 
 // DAWN view seems to have some issue with overlapping solids even if they were unions

From 0bc0f4d2b3db1b88967a7ee3ca232fb2e7b7a18f Mon Sep 17 00:00:00 2001
From: Chao Peng <cpeng@anl.gov>
Date: Mon, 20 Feb 2023 12:30:40 -0600
Subject: [PATCH 12/17] add adjacent layer in the drawings

---
 .../subdetector_tests/draw_bemc_scfi_grids.py | 45 +++++++++----------
 1 file changed, 20 insertions(+), 25 deletions(-)

diff --git a/scripts/subdetector_tests/draw_bemc_scfi_grids.py b/scripts/subdetector_tests/draw_bemc_scfi_grids.py
index 1a032c2e8..0db4c0151 100644
--- a/scripts/subdetector_tests/draw_bemc_scfi_grids.py
+++ b/scripts/subdetector_tests/draw_bemc_scfi_grids.py
@@ -79,48 +79,43 @@ def get_fibers(detelem, id_desc, id_conv, id_dict, grid=1):
             help='Top-level xml file of the detector description.'
             )
     parser.add_argument(
-            '--dpath', default='EcalBarrelScFi/stave0/layer1/slice1',
-            dest='dpath',
-            help='Path to the closest DetElement (slice) that contains the grids and fibers.'
+            '--detector', default='EcalBarrelScFi',
+            dest='detector',
+            help='Detector name.'
             )
     parser.add_argument(
-            '--outdir',
-            dest='outdir', default='.',
-            help='Output directory.'
+            '--readout', default='EcalBarrelScFiHits',
+            help='Readout class for the detector.'
             )
     parser.add_argument(
-            '--readout',
-            default='EcalBarrelScFiHits',
-            help='Readout class for the EcalBarrelScFi.'
+            '--vpath', default='stave1/layer1/slice1/grid3',
+            hlelp='Path down to a grid volume to be centered at the plot'
             )
     parser.add_argument(
-            '--center-grid', dest='center_grid', type=int, default=2,
-            help='The grid number for centering in the plot.'
+            '--outdir',
+            dest='outdir', default='.',
+            help='Output directory.'
             )
     args = parser.parse_args()
 
     # initialize dd4hep detector
     desc = dd4hep.Detector.getInstance()
     desc.fromXML(args.compact)
-    converter = DDRec.CellIDPositionConverter(desc)
 
-    # search the DetElement start from the world
-    det_names = args.dpath.split('/')
+    # search the detector
     det = desc.world()
-
-    for i, n in enumerate(det_names):
-        children = det.children()
-        try:
-            det = det.child(n)
-        except Exception:
-            print('Failed to find DetElement {} from {}'.format(n, '/'.join(['world'] + det_names[:i])))
-            print('Available children are listed below:')
-            for n, d in children:
-                print(' --- child: {}'.format(n))
-            exit(-1)
+    try:
+        det = det.child(args.detector)
+    except Exception:
+        print('Failed to find detector {} from \"{}\"'.format(args.detector, args.compact)
+        print('Available detectors are listed below:')
+        for n, d in desc.world().children:
+            print(' --- detector: {}'.format(n))
+        exit(-1)
 
     # build a volume manager so it triggers populating the volume IDs
     vman = dd4hep.VolumeManager(det, desc.readout(args.readout))
+    converter = DDRec.CellIDPositionConverter(desc)
 
     # manager instances
     readout = desc.readout(args.readout)

From 67ac6d8758cfb012f2230102aa13a58a555e5cf5 Mon Sep 17 00:00:00 2001
From: Chao Peng <cpeng@anl.gov>
Date: Mon, 20 Feb 2023 17:55:28 -0600
Subject: [PATCH 13/17] follow up: finish revamping the visualization script

---
 .../subdetector_tests/draw_bemc_scfi_grids.py | 114 ++++++++++--------
 1 file changed, 61 insertions(+), 53 deletions(-)

diff --git a/scripts/subdetector_tests/draw_bemc_scfi_grids.py b/scripts/subdetector_tests/draw_bemc_scfi_grids.py
index 0db4c0151..052cc8a51 100644
--- a/scripts/subdetector_tests/draw_bemc_scfi_grids.py
+++ b/scripts/subdetector_tests/draw_bemc_scfi_grids.py
@@ -9,6 +9,7 @@
 import DDRec
 import argparse
 import numpy as np
+from pydoc import locate
 from matplotlib import pyplot as plt
 from matplotlib.patches import Circle
 from matplotlib.collections import PatchCollection
@@ -20,30 +21,41 @@
 def dict_to_cpp_vec(my_dict, dtype='int'):
     # use ROOT to make sure the type is correct
     vol_ids = ROOT.std.vector('pair<string, {}>'.format(dtype))()
+    dcast = locate(dtype)
     for field, fid in my_dict.items():
-        vol_ids.push_back((field, fid))
+        vol_ids.push_back((field, dcast(fid)))
     return vol_ids
 
 
 # helper function to collect fibers under a grid
-def get_fibers(detelem, id_desc, id_conv, id_dict, grid=1):
-    gnode = detelem.volume().GetNode(grid-1)
-    global_trans = detelem.nominal().worldTransformation()
-    # print(gnode.GetName())
-    # get grid
-    id_dict.update({'grid': grid, 'fiber': 0})
+def get_grid_fibers(det_elem, vol_man, id_conv, id_dict):
+    # locate the nearest DetElement
+    id_desc = vol_man.idSpec()
+    # get fiber radius
+    id_dict.update({'fiber': 1})
+    fid = id_desc.encode(dict_to_cpp_vec(id_dict))
+    # NOTE: for tube geometry, and it needs a cm to mm conversion
+    fr = id_conv.cellDimensions(fid)[0]/2./10.
+
+    # get the lowest level DetElement
+    sdet = id_conv.findDetElement(id_conv.position(fid))
+    gtrans = sdet.nominal().worldTransformation()
+
+    # get grid node (it's not a DetElement)
+    id_dict.update({'fiber': 0})
     gid = id_desc.encode(dict_to_cpp_vec(id_dict))
-    # print(id_desc.decoder().valueString(vid))
-    vol = id_conv.findContext(gid).volumePlacement().volume()
+    gnode = id_conv.findContext(gid).volumePlacement()
+    # print(id_desc.decoder().valueString(gid))
     grpos = id_conv.position(gid)
-    # print(vol.boundingBox().dimensions())
+    grpos = np.array([grpos.X(), grpos.Y(), grpos.Z()])
 
     # use TGeoNode to get center positions
-    # here it can also use converter to get position but it's much slower (adds an additional lookup process)
+    # here it can also use id_conv to do the same thing with cellIDs,
+    # but it's much slower (adds an additional lookup process)
     fibers = []
-    for i in np.arange(vol.GetNdaughters()):
-        fnode = vol.GetNode(int(i))
-        # NOTE, this is defined in geometry plugin, fiber_core is the wanted sensitive detector
+    for i in np.arange(gnode.GetNdaughters()):
+        fnode = gnode.GetDaughter(int(i))
+        # NOTE, this is defined in geometry plugin, fiber_core is the only wanted sensitive detector
         if 'fiber_core' not in fnode.GetName():
             continue
         fpos = np.array([0., 0., 0.])
@@ -52,7 +64,7 @@ def get_fibers(detelem, id_desc, id_conv, id_dict, grid=1):
         fnode.LocalToMaster(np.array([0., 0., 0.]), fpos)
         gnode.LocalToMaster(fpos, gpos)
         # the two method below are equivalent
-        global_trans.LocalToMaster(gpos, pos)
+        gtrans.LocalToMaster(gpos, pos)
         # detelem.nominal().localToWorld(gpos, pos)
         """ a test with converter
         if i < 50:
@@ -66,9 +78,9 @@ def get_fibers(detelem, id_desc, id_conv, id_dict, grid=1):
                   fpos,
                   gpos)
         """
-        fibers.append(pos)
+        fibers.append(np.hstack([pos, fr]))
 
-    return fibers, grpos
+    return np.array(fibers), grpos
 
 
 if __name__ == '__main__':
@@ -88,8 +100,8 @@ def get_fibers(detelem, id_desc, id_conv, id_dict, grid=1):
             help='Readout class for the detector.'
             )
     parser.add_argument(
-            '--vpath', default='stave1/layer1/slice1/grid3',
-            hlelp='Path down to a grid volume to be centered at the plot'
+            '--grid-path', default='module:1,layer:6,slice:1,grid:3',
+            help='Path down to a grid volume to be centered at the plot, with the format \"field1:i1,field2:i2,...\"',
             )
     parser.add_argument(
             '--outdir',
@@ -107,7 +119,7 @@ def get_fibers(detelem, id_desc, id_conv, id_dict, grid=1):
     try:
         det = det.child(args.detector)
     except Exception:
-        print('Failed to find detector {} from \"{}\"'.format(args.detector, args.compact)
+        print('Failed to find detector {} from \"{}\"'.format(args.detector, args.compact))
         print('Available detectors are listed below:')
         for n, d in desc.world().children:
             print(' --- detector: {}'.format(n))
@@ -117,48 +129,44 @@ def get_fibers(detelem, id_desc, id_conv, id_dict, grid=1):
     vman = dd4hep.VolumeManager(det, desc.readout(args.readout))
     converter = DDRec.CellIDPositionConverter(desc)
 
-    # manager instances
-    readout = desc.readout(args.readout)
-    idspec = readout.idSpec()
-
-    # build a dictionary for the readout ids
-    mvol_id = det.volumeID()
-    id_dict = OrderedDict()
-    for idstr in idspec.decoder().valueString(mvol_id).split(','):
-        field, fid = idstr.split(':')
-        id_dict[field] = int(fid)
-
-    # get radius of the fiber
-    id_dict.update({'grid': args.center_grid, 'fiber': 1})
-    fid = idspec.encode(dict_to_cpp_vec(id_dict))
-    # need to do unit conversion from mm to cm (might from ROOT/DD4Hep difference)
-    fr = converter.cellDimensions(fid)[0]/2./10.
+    fields = OrderedDict([['system', det.id()]] + [v.split(':') for v in args.grid_path.split(',')])
+    layer = int(fields.get('layer'))
+    grid = int(fields.get('grid'))
+    # add adjacent layers and grids, put the central one at the first
+    id_dicts = []
+    for dl in [0, 1, 2, -1, -2]:
+        for dg in [0, 1, 2, -1, -2]:
+            if layer + dl < 1 or grid + dg < 1:
+                continue
+            new_dict = fields.copy()
+            new_dict.update({'layer': layer + dl, 'grid': grid + dg})
+            id_dicts.append(new_dict)
 
     # plot fibers in the grid
     fig, ax = plt.subplots(figsize=(12, 12), dpi=160)
-    # default color cycles
+    # default color cycle
     colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
+    for i, ids in enumerate(id_dicts):
+        c = colors[i % len(colors)]
+        fibers, gr_pos = get_grid_fibers(det, vman, converter, ids)
+        if fibers is None:
+            print('ignored {} because the corresponding volume is not found.'.format(ids))
+            continue
 
-    # plot the center grid and adjacent grid
-    grids = [
-        (g, colors[i]) for i, g in enumerate([args.center_grid, args.center_grid - 1, args.center_grid + 1]) if g > 0
-        ]
-    for g, c in grids:
-        # get position of fibers
-        fibers, grid_pos = get_fibers(det, idspec, converter, id_dict, g)
         patches = []
-        for fpos in fibers:
-            patches.append(Circle((fpos[0], fpos[1]), fr))
+        for fi in fibers:
+            patches.append(Circle((fi[0], fi[1]), fi[3]))
         p = PatchCollection(patches, alpha=0.4, facecolors=(c,), edgecolors=('k',))
-        ax.plot([grid_pos.X()], [grid_pos.Y()], marker='P', mfc=c, mec='k', ms=9, label='grid {} center'.format(g))
+        ax.plot(gr_pos[0], gr_pos[1], marker='P', mfc=c, mec='k', ms=9, label='grid {}'.format(ids['grid']))
         ax.add_collection(p)
-        if g == args.center_grid:
-            ax.set_xlim(grid_pos.X()- 3., grid_pos.X() + 3.)
-            ax.set_ylim(grid_pos.Y() - 3., grid_pos.Y() + 3.)
-    ax.legend(fontsize=22)
+        # center at the first entry
+        if i == 0:
+            ax.set_xlim(gr_pos[0] - 4., gr_pos[0] + 4.)
+            ax.set_ylim(gr_pos[1] - 4., gr_pos[1] + 4.)
+
+    # ax.legend(fontsize=22)
     ax.tick_params(labelsize=20, direction='in')
     ax.set_xlabel('Global X (mm)', fontsize=22)
     ax.set_ylabel('Global Y (mm)', fontsize=22)
-    ax.set_title('{} around grid {}'.format(args.dpath, args.center_grid), fontsize=22)
-
+    ax.set_title('Centered at {}'.format('/'.join(['{}{}'.format(k, v) for k, v in fields.items()])), fontsize=22)
     fig.savefig(os.path.join(args.outdir, 'grid_fibers.png'))

From 597d5ad0cccb086489906bc021d66cc0bd86372c Mon Sep 17 00:00:00 2001
From: Chao Peng <cpeng@anl.gov>
Date: Mon, 20 Feb 2023 22:08:34 -0600
Subject: [PATCH 14/17] determine the fiber line shift with layer_number and
 fiber_line_number

---
 src/BarrelCalorimeterInterlayers_geo.cpp | 16 +++++++++-------
 1 file changed, 9 insertions(+), 7 deletions(-)

diff --git a/src/BarrelCalorimeterInterlayers_geo.cpp b/src/BarrelCalorimeterInterlayers_geo.cpp
index 62dc93007..8f8ed5fd1 100644
--- a/src/BarrelCalorimeterInterlayers_geo.cpp
+++ b/src/BarrelCalorimeterInterlayers_geo.cpp
@@ -48,12 +48,13 @@ struct FiberGrid {
   };
 };
 
-vector<Point> fiberPositions(double r, double sx, double sz, double trx, double trz, double phi, double stol = 1e-2);
+vector<Point> fiberPositions(double r, double sx, double sz, double trx, double trz, double phi,
+                             bool shift_first = false, double stol = 1e-2);
 std::pair<int, int>   getNdivisions(double x, double z, double dx, double dz);
 vector<FiberGrid> gridPoints(int div_x, int div_z, double x, double z, double phi);
 
 // geometry helpers
-void buildFibers(Detector& desc, SensitiveDetector& sens, Volume& mother, xml_comp_t x_fiber,
+void buildFibers(Detector& desc, SensitiveDetector& sens, Volume& mother, int layer_nunber, xml_comp_t x_fiber,
                  const std::tuple<double, double, double, double>& dimensions);
 void buildSupport(Detector& desc, Volume& mother, xml_comp_t x_support,
                   const std::tuple<double, double, double, double>& dimensions);
@@ -138,7 +139,7 @@ static Ref_t create_detector(Detector& desc, xml_h e, SensitiveDetector sens)
 
           // build fibers
           if (x_slice.hasChild(_Unicode(fiber))) {
-            buildFibers(desc, sens, s_vol, x_slice.child(_Unicode(fiber)), {s_trd_x1, s_thick, l_dim_y, hphi});
+            buildFibers(desc, sens, s_vol, l_num, x_slice.child(_Unicode(fiber)), {s_trd_x1, s_thick, l_dim_y, hphi});
           }
 
           if (x_slice.isSensitive()) {
@@ -196,7 +197,7 @@ static Ref_t create_detector(Detector& desc, xml_h e, SensitiveDetector sens)
   return sdet;
 }
 
-void buildFibers(Detector& desc, SensitiveDetector& sens, Volume& s_vol, xml_comp_t x_fiber,
+void buildFibers(Detector& desc, SensitiveDetector& sens, Volume& s_vol, int layer_number, xml_comp_t x_fiber,
                  const std::tuple<double, double, double, double>& dimensions)
 {
   auto [s_trd_x1, s_thick, s_length, hphi] = dimensions;
@@ -231,7 +232,7 @@ void buildFibers(Detector& desc, SensitiveDetector& sens, Volume& s_vol, xml_com
   // Calculate polygonal grid coordinates (vertices)
   auto        grids = gridPoints(grid_div.first, grid_div.second, s_trd_x1, s_thick, hphi);
   vector<int> f_id_count(grid_div.first * grid_div.second, 0);
-  auto        f_pos = fiberPositions(f_radius, f_spacing_x, f_spacing_z, s_trd_x1, s_thick, hphi);
+  auto        f_pos = fiberPositions(f_radius, f_spacing_x, f_spacing_z, s_trd_x1, s_thick, hphi, (layer_number % 2 == 0));
   // a helper struct to speed up searching
   struct Fiber {
     Point pos;
@@ -386,7 +387,7 @@ void buildSupport(Detector& desc, Volume& mod_vol, xml_comp_t x_support,
 }
 
 // Fill fiber lattice into trapezoid starting from position (0,0) in x-z coordinate system
-vector<Point> fiberPositions(double r, double sx, double sz, double trx, double trz, double phi, double stol)
+vector<Point> fiberPositions(double r, double sx, double sz, double trx, double trz, double phi, bool shift, double stol)
 {
   // r      - fiber radius
   // sx, sz - spacing between fibers in x, z
@@ -399,12 +400,13 @@ vector<Point> fiberPositions(double r, double sx, double sz, double trx, double
   int z_layers = floor((trz / 2 - r - stol) / sz); // number of layers that fits in half trapezoid-z
 
   double px = 0., pz = 0.;
+  int start_line = shift ? 1 : 0;
 
   for (int l = -z_layers; l < z_layers + 1; l++) {
     vector<Point> xline;
     pz           = l * sz;
     double x_max = trx + (trz / 2. + pz) * tan(phi) - stol; // calculate max x at particular z_pos
-    (l % 2 == 0) ? px = 0. : px = sx / 2;                   // account for spacing/2 shift
+    (abs(l) % 2 == start_line) ? px = 0. : px = sx / 2;          // account for spacing/2 shift
 
     while (px < (x_max - r)) {
       xline.push_back(Point(px, pz));

From 17fd9a367c33f56d8afc7a04621130aebf7af8be Mon Sep 17 00:00:00 2001
From: Chao Peng <cpeng@anl.gov>
Date: Mon, 20 Feb 2023 23:13:28 -0600
Subject: [PATCH 15/17] update the visualization script, add more arguments

---
 .../subdetector_tests/draw_bemc_scfi_grids.py | 67 ++++++++++++-------
 1 file changed, 44 insertions(+), 23 deletions(-)

diff --git a/scripts/subdetector_tests/draw_bemc_scfi_grids.py b/scripts/subdetector_tests/draw_bemc_scfi_grids.py
index 052cc8a51..600734817 100644
--- a/scripts/subdetector_tests/draw_bemc_scfi_grids.py
+++ b/scripts/subdetector_tests/draw_bemc_scfi_grids.py
@@ -1,5 +1,8 @@
 '''
     A script to visualize the fibers of some grids from BEMC ScFi part
+    use case:
+    python scripts/subdetector_tests/draw_bemc_scfi_grids.py -c epic_brycecanyon.xml
+
     02/19/2023
     Chao Peng (ANL)
 '''
@@ -31,23 +34,26 @@ def dict_to_cpp_vec(my_dict, dtype='int'):
 def get_grid_fibers(det_elem, vol_man, id_conv, id_dict):
     # locate the nearest DetElement
     id_desc = vol_man.idSpec()
-    # get fiber radius
-    id_dict.update({'fiber': 1})
-    fid = id_desc.encode(dict_to_cpp_vec(id_dict))
-    # NOTE: for tube geometry, and it needs a cm to mm conversion
-    fr = id_conv.cellDimensions(fid)[0]/2./10.
-
-    # get the lowest level DetElement
-    sdet = id_conv.findDetElement(id_conv.position(fid))
-    gtrans = sdet.nominal().worldTransformation()
-
-    # get grid node (it's not a DetElement)
-    id_dict.update({'fiber': 0})
-    gid = id_desc.encode(dict_to_cpp_vec(id_dict))
-    gnode = id_conv.findContext(gid).volumePlacement()
-    # print(id_desc.decoder().valueString(gid))
-    grpos = id_conv.position(gid)
-    grpos = np.array([grpos.X(), grpos.Y(), grpos.Z()])
+    try:
+        # get fiber radius
+        id_dict.update({'fiber': 1})
+        fid = id_desc.encode(dict_to_cpp_vec(id_dict))
+        # NOTE: for tube geometry, and it needs a cm to mm conversion
+        fr = id_conv.cellDimensions(fid)[0]/2./10.
+
+        # get the lowest level DetElement
+        sdet = id_conv.findDetElement(id_conv.position(fid))
+        gtrans = sdet.nominal().worldTransformation()
+
+        # get grid node (it's not a DetElement)
+        id_dict.update({'fiber': 0})
+        gid = id_desc.encode(dict_to_cpp_vec(id_dict))
+        gnode = id_conv.findContext(gid).volumePlacement()
+        # print(id_desc.decoder().valueString(gid))
+        grpos = id_conv.position(gid)
+        grpos = np.array([grpos.X(), grpos.Y(), grpos.Z()])
+    except Exception:
+        return None, None
 
     # use TGeoNode to get center positions
     # here it can also use id_conv to do the same thing with cellIDs,
@@ -108,6 +114,21 @@ def get_grid_fibers(det_elem, vol_man, id_conv, id_dict):
             dest='outdir', default='.',
             help='Output directory.'
             )
+    parser.add_argument(
+            '--adj-nlayers',
+            dest='nlayers', type=int, default=2,
+            help='number of adjacent layers to draw (+-n).'
+            )
+    parser.add_argument(
+            '--adj-ngrids',
+            dest='ngrids', type=int, default=2,
+            help='number of adjacent grids to draw (+-n).'
+            )
+    parser.add_argument(
+            '--window-size',
+            dest='wsize', type=float, default=4.,
+            help='Plot window size (mm).'
+            )
     args = parser.parse_args()
 
     # initialize dd4hep detector
@@ -132,10 +153,10 @@ def get_grid_fibers(det_elem, vol_man, id_conv, id_dict):
     fields = OrderedDict([['system', det.id()]] + [v.split(':') for v in args.grid_path.split(',')])
     layer = int(fields.get('layer'))
     grid = int(fields.get('grid'))
-    # add adjacent layers and grids, put the central one at the first
+    # add adjacent layers and grids, always put the central one (0, 0) first
     id_dicts = []
-    for dl in [0, 1, 2, -1, -2]:
-        for dg in [0, 1, 2, -1, -2]:
+    for dl in np.hstack([np.arange(0, args.nlayers + 1), np.arange(-1, -args.nlayers - 1, step=-1)]):
+        for dg in np.hstack([np.arange(0, args.ngrids + 1), np.arange(-1, -args.ngrids - 1, step=-1)]):
             if layer + dl < 1 or grid + dg < 1:
                 continue
             new_dict = fields.copy()
@@ -150,7 +171,7 @@ def get_grid_fibers(det_elem, vol_man, id_conv, id_dict):
         c = colors[i % len(colors)]
         fibers, gr_pos = get_grid_fibers(det, vman, converter, ids)
         if fibers is None:
-            print('ignored {} because the corresponding volume is not found.'.format(ids))
+            print('ignored {} because the volume might not exist.'.format(ids))
             continue
 
         patches = []
@@ -161,8 +182,8 @@ def get_grid_fibers(det_elem, vol_man, id_conv, id_dict):
         ax.add_collection(p)
         # center at the first entry
         if i == 0:
-            ax.set_xlim(gr_pos[0] - 4., gr_pos[0] + 4.)
-            ax.set_ylim(gr_pos[1] - 4., gr_pos[1] + 4.)
+            ax.set_xlim(gr_pos[0] - args.wsize, gr_pos[0] + args.wsize)
+            ax.set_ylim(gr_pos[1] - args.wsize, gr_pos[1] + args.wsize)
 
     # ax.legend(fontsize=22)
     ax.tick_params(labelsize=20, direction='in')

From 38c228c9f9f239721080466018646a979781dc43 Mon Sep 17 00:00:00 2001
From: Chao Peng <cpeng@anl.gov>
Date: Tue, 21 Feb 2023 15:49:49 -0600
Subject: [PATCH 16/17] add SPDX header

---
 scripts/subdetector_tests/draw_bemc_scfi_grids.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/scripts/subdetector_tests/draw_bemc_scfi_grids.py b/scripts/subdetector_tests/draw_bemc_scfi_grids.py
index 600734817..96791de64 100644
--- a/scripts/subdetector_tests/draw_bemc_scfi_grids.py
+++ b/scripts/subdetector_tests/draw_bemc_scfi_grids.py
@@ -1,10 +1,9 @@
+# SPDX-License-Identifier: LGPL-3.0-or-later
+# Copyright (C) 2022 Chao Peng
 '''
     A script to visualize the fibers of some grids from BEMC ScFi part
     use case:
     python scripts/subdetector_tests/draw_bemc_scfi_grids.py -c epic_brycecanyon.xml
-
-    02/19/2023
-    Chao Peng (ANL)
 '''
 import os
 import ROOT

From cddef5a01e000a595407b75f94bc6418e14e8107 Mon Sep 17 00:00:00 2001
From: Wouter Deconinck <wdconinc@gmail.com>
Date: Tue, 21 Feb 2023 15:56:01 -0600
Subject: [PATCH 17/17] it's 2023 already

---
 scripts/subdetector_tests/draw_bemc_scfi_grids.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/scripts/subdetector_tests/draw_bemc_scfi_grids.py b/scripts/subdetector_tests/draw_bemc_scfi_grids.py
index 96791de64..5bd43e729 100644
--- a/scripts/subdetector_tests/draw_bemc_scfi_grids.py
+++ b/scripts/subdetector_tests/draw_bemc_scfi_grids.py
@@ -1,5 +1,5 @@
 # SPDX-License-Identifier: LGPL-3.0-or-later
-# Copyright (C) 2022 Chao Peng
+# Copyright (C) 2023 Chao Peng
 '''
     A script to visualize the fibers of some grids from BEMC ScFi part
     use case: