Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2024-06-16 07:06:58

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Chao Peng, Maria Zurek, Whitney Armstrong
0003 
0004 // Detector plugin to support a hybrid central barrel calorimeter
0005 // The detector consists of interlayers of Pb/ScFi (segmentation in global r, phi) and W/Si (segmentation in local x, y)
0006 // Assembly is used as the envelope so two different detectors can be interlayered with each other
0007 //
0008 //
0009 // 06/19/2021: Implementation of the Sci Fiber geometry. M. Żurek
0010 // 07/09/2021: Support interlayers between multiple detectors. C. Peng
0011 // 07/23/2021: Add assemblies as mother volumes of fibers to reduce the number of daughter volumes. C. Peng, M. Żurek
0012 //     Reference: TGeo performance issue with large number of daughter volumes
0013 //     https://indico.cern.ch/event/967418/contributions/4075358/attachments/2128099/3583278/201009_shKo_dd4hep.pdf
0014 // 07/24/2021: Changed support implementation to avoid too many uses of boolean geometries. DAWN view seems to have
0015 //     issue dealing with it. C. Peng
0016 
0017 #include "DD4hep/DetFactoryHelper.h"
0018 #include "Math/Point2D.h"
0019 #include "TGeoPolygon.h"
0020 #include "XML/Layering.h"
0021 #include <functional>
0022 
0023 using namespace dd4hep;
0024 
0025 typedef ROOT::Math::XYPoint Point;
0026 // fiber placement helpers, defined below
0027 struct FiberGrid {
0028   int ix = 0, iy = 0;
0029   std::vector<Point> points;
0030   Point mean_centroid    = Point(0., 0.);
0031   Assembly* assembly_ptr = nullptr;
0032 
0033   // initialize with grid id and points
0034   FiberGrid(int i, int j, const std::vector<Point>& pts) : ix(i), iy(j), points(pts) {
0035     if (pts.empty()) {
0036       return;
0037     }
0038 
0039     double mx = 0., my = 0.;
0040     for (auto& p : pts) {
0041       mx += p.x();
0042       my += p.y();
0043     }
0044     mx /= static_cast<double>(pts.size());
0045     my /= static_cast<double>(pts.size());
0046     mean_centroid = Point(mx, my);
0047   };
0048 };
0049 
0050 std::vector<Point> fiberPositions(double r, double sx, double sz, double trx, double trz,
0051                                   double phi, bool shift_first = false, double stol = 1e-2);
0052 std::vector<FiberGrid> gridPoints(int div_n_phi, double div_dr, double x, double z, double phi);
0053 
0054 // geometry helpers
0055 void buildFibers(Detector& desc, SensitiveDetector& sens, Volume& mother, int layer_nunber,
0056                  xml_comp_t x_fiber, const std::tuple<double, double, double, double>& dimensions);
0057 void buildSupport(Detector& desc, Volume& mother, xml_comp_t x_support,
0058                   const std::tuple<double, double, double, double>& dimensions);
0059 
0060 // barrel ecal layers contained in an assembly
0061 static Ref_t create_detector(Detector& desc, xml_h e, SensitiveDetector sens) {
0062   Layering layering(e);
0063   xml_det_t x_det      = e;
0064   Material air         = desc.air();
0065   int det_id           = x_det.id();
0066   double offset        = x_det.attr<double>(_Unicode(offset));
0067   xml_comp_t x_dim     = x_det.dimensions();
0068   int nsides           = x_dim.numsides();
0069   double inner_r       = x_dim.rmin();
0070   double dphi          = (2 * M_PI / nsides);
0071   double hphi          = dphi / 2;
0072   std::string det_name = x_det.nameStr();
0073 
0074   DetElement sdet(det_name, det_id);
0075   Volume motherVol = desc.pickMotherVolume(sdet);
0076 
0077   Assembly envelope(det_name);
0078   Transform3D tr_global = Translation3D(0, 0, offset) * RotationZ(hphi);
0079   PlacedVolume env_phv  = motherVol.placeVolume(envelope, tr_global);
0080   sens.setType("calorimeter");
0081 
0082   env_phv.addPhysVolID("system", det_id);
0083   sdet.setPlacement(env_phv);
0084 
0085   // build a single sector
0086   DetElement sector_det("sector0", det_id);
0087   Assembly mod_vol("sector");
0088 
0089   // keep tracking of the total thickness
0090   double l_pos_z = inner_r;
0091   { // =====  buildBarrelStave(desc, sens, module_volume) =====
0092     // Parameters for computing the layer X dimension:
0093     double tan_hphi = std::tan(hphi);
0094     double l_dim_y  = x_dim.z() / 2.;
0095 
0096     // Loop over the sets of layer elements in the detector.
0097     int l_num = 1;
0098     for (xml_coll_t li(x_det, _U(layer)); li; ++li) {
0099       xml_comp_t x_layer     = li;
0100       int repeat             = x_layer.repeat();
0101       double l_space_between = getAttrOrDefault(x_layer, _Unicode(space_between), 0.);
0102       double l_space_before  = getAttrOrDefault(x_layer, _Unicode(space_before), 0.);
0103       l_pos_z += l_space_before;
0104       // Loop over number of repeats for this layer.
0105       for (int j = 0; j < repeat; j++) {
0106         std::string l_name = Form("layer%d", l_num);
0107         double l_thickness = layering.layer(l_num - 1)->thickness(); // Layer's thickness.
0108         double l_dim_x     = tan_hphi * l_pos_z;
0109         l_pos_z += l_thickness;
0110 
0111         Position l_pos(0, 0, l_pos_z - l_thickness / 2.); // Position of the layer.
0112         double l_trd_x1 = l_dim_x;
0113         double l_trd_x2 = l_dim_x + l_thickness * tan_hphi;
0114         double l_trd_y1 = l_dim_y;
0115         double l_trd_y2 = l_trd_y1;
0116         double l_trd_z  = l_thickness / 2;
0117         Trapezoid l_shape(l_trd_x1, l_trd_x2, l_trd_y1, l_trd_y2, l_trd_z);
0118         Volume l_vol(l_name, l_shape, air);
0119         DetElement layer(sector_det, l_name, det_id);
0120 
0121         // Loop over the sublayers or slices for this layer.
0122         int s_num      = 1;
0123         double s_pos_z = -(l_thickness / 2.);
0124         for (xml_coll_t si(x_layer, _U(slice)); si; ++si) {
0125           xml_comp_t x_slice = si;
0126           std::string s_name = Form("slice%d", s_num);
0127           double s_thick     = x_slice.thickness();
0128           double s_trd_x1    = l_dim_x + (s_pos_z + l_thickness / 2) * tan_hphi;
0129           double s_trd_x2    = l_dim_x + (s_pos_z + l_thickness / 2 + s_thick) * tan_hphi;
0130           double s_trd_y1    = l_trd_y1;
0131           double s_trd_y2    = s_trd_y1;
0132           double s_trd_z     = s_thick / 2.;
0133           Trapezoid s_shape(s_trd_x1, s_trd_x2, s_trd_y1, s_trd_y2, s_trd_z);
0134           Volume s_vol(s_name, s_shape, desc.material(x_slice.materialStr()));
0135           DetElement slice(layer, s_name, det_id);
0136 
0137           // build fibers
0138           if (x_slice.hasChild(_Unicode(fiber))) {
0139             buildFibers(desc, sens, s_vol, l_num, x_slice.child(_Unicode(fiber)),
0140                         {s_trd_x1, s_thick, l_dim_y, hphi});
0141           }
0142 
0143           if (x_slice.isSensitive()) {
0144             s_vol.setSensitiveDetector(sens);
0145           }
0146           s_vol.setAttributes(desc, x_slice.regionStr(), x_slice.limitsStr(), x_slice.visStr());
0147 
0148           // Slice placement.
0149           PlacedVolume slice_phv = l_vol.placeVolume(s_vol, Position(0, 0, s_pos_z + s_thick / 2));
0150           slice_phv.addPhysVolID("slice", s_num);
0151           slice.setPlacement(slice_phv);
0152           // Increment Z position of slice.
0153           s_pos_z += s_thick;
0154           ++s_num;
0155         }
0156 
0157         // Set region, limitset, and vis of layer.
0158         l_vol.setAttributes(desc, x_layer.regionStr(), x_layer.limitsStr(), x_layer.visStr());
0159 
0160         PlacedVolume layer_phv = mod_vol.placeVolume(l_vol, l_pos);
0161         layer_phv.addPhysVolID("layer", l_num);
0162         layer.setPlacement(layer_phv);
0163         // Increment to next layer Z position. Do not add space_between for the last layer
0164         if (j < repeat - 1) {
0165           l_pos_z += l_space_between;
0166         }
0167         ++l_num;
0168       }
0169     }
0170   }
0171   // Phi start for a sector.
0172   double phi = M_PI / nsides;
0173   // Create nsides sectors.
0174   for (int i = 0; i < nsides; i++, phi -= dphi) { // i is module number
0175     // Compute the sector position
0176     Transform3D tr(RotationZYX(0, phi, M_PI * 0.5), Translation3D(0, 0, 0));
0177     PlacedVolume pv = envelope.placeVolume(mod_vol, tr);
0178     pv.addPhysVolID("sector", i + 1);
0179     DetElement sd = (i == 0) ? sector_det : sector_det.clone(Form("sector%d", i));
0180     sd.setPlacement(pv);
0181     sdet.add(sd);
0182   }
0183 
0184   // optional sector attributes
0185   if (x_det.hasChild(_Unicode(sectors))) {
0186     xml_comp_t x_sectors = x_det.child(_Unicode(sectors));
0187     mod_vol.setVisAttributes(desc.visAttributes(x_sectors.visStr()));
0188   }
0189   if (x_det.hasChild(_U(support))) {
0190     buildSupport(desc, mod_vol, x_det.child(_U(support)), {inner_r, l_pos_z, x_dim.z(), hphi});
0191   }
0192 
0193   // Set envelope volume attributes.
0194   envelope.setAttributes(desc, x_det.regionStr(), x_det.limitsStr(), x_det.visStr());
0195   return sdet;
0196 }
0197 
0198 void buildFibers(Detector& desc, SensitiveDetector& sens, Volume& s_vol, int layer_number,
0199                  xml_comp_t x_fiber, const std::tuple<double, double, double, double>& dimensions) {
0200   auto [s_trd_x1, s_thick, s_length, hphi] = dimensions;
0201   double f_radius                          = getAttrOrDefault(x_fiber, _U(radius), 0.1 * cm);
0202   double f_cladding_thickness = getAttrOrDefault(x_fiber, _Unicode(cladding_thickness), 0.0 * cm);
0203   double f_spacing_x          = getAttrOrDefault(x_fiber, _Unicode(spacing_x), 0.122 * cm);
0204   double f_spacing_z          = getAttrOrDefault(x_fiber, _Unicode(spacing_z), 0.134 * cm);
0205   int grid_n_phi              = getAttrOrDefault(x_fiber, _Unicode(grid_n_phi), 5);
0206   double grid_dr              = getAttrOrDefault(x_fiber, _Unicode(grid_dr), 2.0 * cm);
0207   std::string f_id_grid = getAttrOrDefault<std::string>(x_fiber, _Unicode(identifier_grid), "grid");
0208   std::string f_id_fiber =
0209       getAttrOrDefault<std::string>(x_fiber, _Unicode(identifier_fiber), "fiber");
0210 
0211   // Set up the readout grid for the fiber layers
0212   // Trapezoid is divided into segments with equal dz and equal number of divisions in x
0213   // Every segment is a polygon that can be attached later to the lightguide
0214   // The grid size is assumed to be ~2x2 cm (starting values). This is to be larger than
0215   // SiPM chip (for GlueX 13mmx13mm: 4x4 grid 3mmx3mm with 3600 50×50 μm pixels each)
0216   // See, e.g., https://arxiv.org/abs/1801.03088 Fig. 2d
0217 
0218   // fiber and its cladding
0219   double f_radius_core = f_radius - f_cladding_thickness;
0220   Tube f_tube_clad(0, f_radius, s_length);
0221   Volume f_vol_clad("fiber_vol", f_tube_clad, desc.material(x_fiber.materialStr()));
0222   Tube f_tube_core(0, f_radius_core, s_length);
0223   Volume f_vol_core("fiber_core_vol", f_tube_core, desc.material(x_fiber.materialStr()));
0224   if (x_fiber.isSensitive()) {
0225     f_vol_core.setSensitiveDetector(sens);
0226   }
0227   f_vol_core.setAttributes(desc, x_fiber.regionStr(), x_fiber.limitsStr(), x_fiber.visStr());
0228   f_vol_clad.placeVolume(f_vol_core);
0229 
0230   // calculate polygonal grid coordinates (vertices)
0231   auto grids = gridPoints(grid_n_phi, grid_dr, s_trd_x1, s_thick, hphi);
0232   std::vector<int> f_id_count(grids.size(), 0);
0233   // use layer_number % 2 to add correct shifts for the adjacent fibers at layer boundary
0234   auto f_pos = fiberPositions(f_radius, f_spacing_x, f_spacing_z, s_trd_x1, s_thick, hphi,
0235                               (layer_number % 2 == 0));
0236   // a helper struct to speed up searching
0237   struct Fiber {
0238     Point pos;
0239     bool assigned = false;
0240     Fiber(const Point& p) : pos(p) {};
0241   };
0242   std::vector<Fiber> fibers(f_pos.begin(), f_pos.end());
0243 
0244   // build assembly for each grid and put fibers in
0245   for (auto& gr : grids) {
0246     Assembly grid_vol(Form("fiber_grid_%i_%i", gr.ix, gr.iy));
0247 
0248     // loop over all fibers that are not assigned to a grid
0249     int f_id = 1;
0250     for (auto& fi : fibers) {
0251       if (fi.assigned) {
0252         continue;
0253       }
0254 
0255       // use TGeoPolygon to help check if fiber is inside a grid
0256       TGeoPolygon poly(gr.points.size());
0257       std::vector<double> vx, vy;
0258       std::transform(gr.points.begin(), gr.points.end(), back_inserter(vx), std::mem_fn(&Point::x));
0259       std::transform(gr.points.begin(), gr.points.end(), back_inserter(vy), std::mem_fn(&Point::y));
0260       poly.SetXY(vx.data(), vy.data());
0261       poly.FinishPolygon();
0262 
0263       double f_xy[2] = {fi.pos.x(), fi.pos.y()};
0264       if (not poly.Contains(f_xy)) {
0265         continue;
0266       }
0267 
0268       // place fiber in grid
0269       auto p        = fi.pos - gr.mean_centroid;
0270       auto clad_phv = grid_vol.placeVolume(f_vol_clad, Position(p.x(), p.y(), 0.));
0271       clad_phv.addPhysVolID(f_id_fiber, f_id);
0272       fi.assigned = true;
0273       f_id++;
0274     }
0275 
0276     // only add this if this grid has fibers
0277     if (f_id > 1) {
0278       // fiber is along y-axis of the layer volume, so grids are arranged on X-Z plane
0279       Transform3D gr_tr(RotationZYX(0., 0., M_PI * 0.5),
0280                         Position(gr.mean_centroid.x(), 0., gr.mean_centroid.y()));
0281       auto grid_phv = s_vol.placeVolume(grid_vol, gr_tr);
0282       grid_phv.addPhysVolID(f_id_grid, gr.ix + gr.iy * grid_n_phi + 1);
0283       grid_vol.ptr()->Voxelize("");
0284     }
0285   }
0286 
0287   /*
0288   // sanity check
0289   size_t missing_fibers = 0;
0290   for (auto &fi : fibers) {
0291     if (not fi.assigned) {
0292       missing_fibers++;
0293     }
0294   }
0295   std::cout << "built " << fibers.size() << " fibers, "
0296             << missing_fibers << " of them failed to find a grid" << std::endl;
0297   */
0298 }
0299 
0300 // simple aluminum sheet cover
0301 // dimensions: (inner r, position in z, length, phi)
0302 void buildSupport(Detector& desc, Volume& mod_vol, xml_comp_t x_support,
0303                   const std::tuple<double, double, double, double>& dimensions) {
0304   auto [inner_r, pos_z, sector_length, hphi] = dimensions;
0305   double support_thickness = getAttrOrDefault(x_support, _Unicode(thickness), 3. * cm);
0306   auto material            = desc.material(x_support.materialStr());
0307   double trd_y             = sector_length / 2.;
0308   double trd_x1_support    = std::tan(hphi) * pos_z;
0309   double trd_x2_support    = std::tan(hphi) * (pos_z + support_thickness);
0310 
0311   Trapezoid s_shape(trd_x1_support, trd_x2_support, trd_y, trd_y, support_thickness / 2.);
0312   Volume s_vol("support_layer", s_shape, material);
0313   s_vol.setVisAttributes(desc.visAttributes(x_support.visStr()));
0314   mod_vol.placeVolume(s_vol, Position(0.0, 0.0, pos_z + support_thickness / 2.));
0315 }
0316 
0317 // Fill fiber lattice into trapezoid starting from position (0,0) in x-z coordinate system
0318 std::vector<Point> fiberPositions(double r, double sx, double sz, double trx, double trz,
0319                                   double phi, bool shift, double stol) {
0320   // r      - fiber radius
0321   // sx, sz - spacing between fibers in x, z
0322   // trx    - half-length of the shorter (bottom) base of the trapezoid
0323   // trz    - height of the trapezoid
0324   // phi    - angle between z and trapezoid arm
0325   // stol   - spacing tolerance
0326 
0327   std::vector<Point> positions;
0328   int z_layers = floor((trz / 2 - r - stol) / sz); // number of layers that fits in half trapezoid-z
0329 
0330   double px = 0., pz = 0.;
0331   int start_line = shift ? 1 : 0;
0332 
0333   for (int l = -z_layers; l < z_layers + 1; l++) {
0334     std::vector<Point> xline;
0335     pz           = l * sz;
0336     double x_max = trx + (trz / 2. + pz) * tan(phi) - stol; // calculate max x at particular z_pos
0337     (abs(l) % 2 == start_line) ? px = 0. : px = sx / 2;     // account for spacing/2 shift
0338 
0339     while (px < (x_max - r)) {
0340       xline.push_back(Point(px, pz));
0341       if (px != 0.)
0342         xline.push_back(Point(-px, pz)); // using symmetry around x=0
0343       px += sx;
0344     }
0345 
0346     // Sort fiber IDs for a better organization
0347     sort(xline.begin(), xline.end(),
0348          [](const Point& p1, const Point& p2) { return p1.x() < p2.x(); });
0349     positions.insert(positions.end(), xline.begin(), xline.end());
0350   }
0351   return positions;
0352 }
0353 
0354 // Determine the number of divisions for the readout grid for the fiber layers
0355 // Calculate dimensions of the polygonal grid
0356 std::vector<FiberGrid> gridPoints(int div_n_phi, double div_dr, double trd_x1, double height,
0357                                   double phi) {
0358   /*
0359   // TODO: move this test to xml file
0360   double SiPMsize = 13.0 * mm;
0361   double grid_min = SiPMsize + 3.0 * mm;
0362 
0363   if (dz < grid_min) {
0364     dz = grid_min;
0365   }
0366 
0367   if (dx < grid_min) {
0368     dx = grid_min;
0369   }
0370   */
0371   // number of divisions
0372   int nph = div_n_phi;
0373   int nr  = floor(height / div_dr);
0374   if (nr == 0) {
0375     nr++;
0376   }
0377 
0378   // grid vertices
0379   std::vector<FiberGrid> results;
0380   double dr = height / nr;
0381 
0382   for (int ir = 0; ir <= nr; ir++) {
0383     for (int iph = 0; iph <= nph; iph++) {
0384       double A_y = -height / 2. + ir * dr;
0385       double B_y = -height / 2. + (ir + 1) * dr;
0386 
0387       double botl_dr = 2 * (trd_x1 + ir * dr * tan(phi));
0388       double topl_dr = 2 * (trd_x1 + (ir + 1) * dr * tan(phi));
0389 
0390       double botl_dph_dr = botl_dr / nph;
0391       double topl_dph_dr = topl_dr / nph;
0392 
0393       double A_x = -botl_dr / 2. + iph * botl_dph_dr;
0394       double B_x = -topl_dr / 2. + iph * topl_dph_dr;
0395 
0396       double C_y = B_y;
0397       double D_y = A_y;
0398       double C_x = B_x + topl_dph_dr;
0399       double D_x = A_x + botl_dph_dr;
0400 
0401       auto A = Point(A_x, A_y);
0402       auto B = Point(B_x, B_y);
0403       auto C = Point(C_x, C_y);
0404       auto D = Point(D_x, D_y);
0405 
0406       // vertex points filled in the clock-wise direction
0407       results.emplace_back(FiberGrid(iph, ir, {A, B, C, D}));
0408     }
0409   }
0410 
0411   return results;
0412 }
0413 
0414 DECLARE_DETELEMENT(epic_EcalBarrelScFi, create_detector)