epic/src/BarrelCalorimeterImaging_geo.cpp

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 "DD4hep/Printout.h"
0019 #include "DDRec/Surface.h"
0020 #include "Math/Point2D.h"
0021 #include "TGeoPolygon.h"
0022 #include "XML/Layering.h"
0023 #include <functional>
0024
0025 #include "DD4hepDetectorHelper.h"
0026
0027 using namespace dd4hep;
0028
0029 typedef ROOT::Math::XYPoint Point;
0030
0031 // geometry helpers
0032 static void buildSupport(Detector& desc, Volume& mother, xml_comp_t x_support,
0033                          const std::tuple<double, double, double, double>& dimensions);
0034
0035 // barrel ecal layers contained in an assembly
0036 static Ref_t create_detector(Detector& desc, xml_h e, SensitiveDetector sens) {
0037   xml_det_t x_detector      = e;
0038   int detector_id           = x_detector.id();
0039   std::string detector_name = x_detector.nameStr();
0040   double offset             = x_detector.attr<double>(_Unicode(offset));
0041   xml_comp_t x_dimensions   = x_detector.dimensions();
0042   int nsides                = x_dimensions.numsides();
0043   double inner_r            = x_dimensions.rmin();
0044   double dphi               = (2 * M_PI / nsides);
0045   double half_dphi          = dphi / 2;
0046
0047   DetElement sdet(detector_name, detector_id);
0048   Volume mother_volume = desc.pickMotherVolume(sdet);
0049
0050   Assembly detector_volume(detector_name);
0051   detector_volume.setAttributes(desc, x_detector.regionStr(), x_detector.limitsStr(),
0052                                 x_detector.visStr());
0053
0054   Transform3D detector_tr       = Translation3D(0, 0, offset) * RotationZ(half_dphi);
0055   PlacedVolume detector_physvol = mother_volume.placeVolume(detector_volume, detector_tr);
0056   sens.setType("calorimeter");
0057
0058   detector_physvol.addPhysVolID("system", detector_id);
0059   sdet.setPlacement(detector_physvol);
0060
0061   // build a single sector
0062   DetElement sector_element("sector0", detector_id);
0063   Assembly sector_volume("sector");
0064   if (x_detector.hasChild(_Unicode(sectors))) {
0065     xml_comp_t x_sectors = x_detector.child(_Unicode(sectors));
0066     sector_volume.setVisAttributes(desc.visAttributes(x_sectors.visStr()));
0067   }
0068
0069   // keep tracking of the total thickness
0070   double layer_pos_z = inner_r;
0071
0072   // Parameters for computing the layer X dimension:
0073   double tan_half_dphi = std::tan(half_dphi);
0074   double layer_dim_y   = x_dimensions.z() / 2.;
0075
0076   // Loop over the modules
0077   using dd4hep::rec::VolPlane;
0078   std::map<std::string, Volume> volumes;
0079   std::map<std::string, std::vector<PlacedVolume>> sensitives;
0080   std::map<std::string, std::vector<VolPlane>> volplane_surfaces;
0081   std::map<std::string, std::array<double, 2>> module_thicknesses;
0082   for (xml_coll_t i_module(x_detector, _U(module)); i_module; ++i_module) {
0083     xml_comp_t x_module     = i_module;
0084     std::string module_name = x_module.nameStr();
0085
0086     if (volumes.find(module_name) != volumes.end()) {
0087       printout(ERROR, "BarrelCalorimeterImaging", "Module with name %s already exists",
0088                module_name.c_str());
0089       throw std::runtime_error("Logics error in building modules.");
0090     }
0091
0092     // Compute module total thickness from components
0093     double total_thickness = 0;
0094     xml_coll_t ci(x_module, _U(module_component));
0095     for (ci.reset(), total_thickness = 0.0; ci; ++ci) {
0096       total_thickness += xml_comp_t(ci).thickness();
0097     }
0098
0099     // Create the module assembly volume
0100     Assembly module_volume(module_name);
0101     volumes[module_name] = module_volume;
0102     module_volume.setVisAttributes(desc.visAttributes(x_module.visStr()));
0103
0104     // Add components
0105     int sensor_number       = 1;
0106     int ncomponents         = 0;
0107     double thickness_so_far = 0.0;
0108     double thickness_sum    = -total_thickness / 2.0;
0109     for (xml_coll_t i_component(x_module, _U(module_component)); i_component;
0110          ++i_component, ++ncomponents) {
0111       xml_comp_t x_component           = i_component;
0112       xml_comp_t x_component_pos       = x_component.position(false);
0113       xml_comp_t x_component_rot       = x_component.rotation(false);
0114       const std::string component_name = _toString(ncomponents, "component%d");
0115       Box component_box(x_component.width() / 2, x_component.length() / 2,
0116                         x_component.thickness() / 2);
0117       Volume component_volume(component_name, component_box,
0118                               desc.material(x_component.materialStr()));
0119       component_volume.setAttributes(desc, x_component.regionStr(), x_component.limitsStr(),
0120                                      x_component.visStr());
0121
0122       // Loop over the slices for this component
0123       int slice_num      = 1;
0124       double slice_pos_z = -x_component.thickness() / 2;
0125       for (xml_coll_t i_slice(x_component, _U(slice)); i_slice; ++i_slice) {
0126         xml_comp_t x_slice     = i_slice;
0127         std::string slice_name = Form("slice%d", slice_num);
0128         double slice_dim_x     = x_component.width() / 2;
0129         double slice_dim_y     = x_component.length() / 2;
0130         double slice_dim_z     = x_slice.thickness() / 2;
0131         Box slice_shape(slice_dim_x, slice_dim_y, slice_dim_z);
0132         Volume slice_volume(slice_name, slice_shape, desc.material(x_slice.materialStr()));
0133         slice_volume.setAttributes(desc, x_slice.regionStr(), x_slice.limitsStr(),
0134                                    x_slice.visStr());
0135
0136         // Place slice
0137         PlacedVolume slice_physvol = component_volume.placeVolume(
0138             slice_volume, Position(0, 0, slice_pos_z + x_slice.thickness() / 2));
0139
0140         // Set sensitive
0141         if (x_slice.isSensitive()) {
0142           slice_physvol.addPhysVolID("slice", slice_num);
0143           slice_volume.setSensitiveDetector(sens);
0144           sensitives[module_name].push_back(slice_physvol);
0145         }
0146
0147         // Increment Z position of slice
0148         slice_pos_z += x_slice.thickness();
0149         ++slice_num;
0150       }
0151
0152       // Utility variable for the relative z-offset based off the previous components
0153       const double zoff = thickness_sum + x_component.thickness() / 2.0;
0154       PlacedVolume component_physvol;
0155       if (x_component_pos && x_component_rot) {
0156         Position component_pos(x_component_pos.x(0), x_component_pos.y(0),
0157                                x_component_pos.z(0) + zoff);
0158         RotationZYX component_rot(x_component_rot.z(0), x_component_rot.y(0), x_component_rot.x(0));
0159         component_physvol =
0160             module_volume.placeVolume(component_volume, Transform3D(component_rot, component_pos));
0161       } else if (x_component_rot) {
0162         Position component_pos(0, 0, zoff);
0163         RotationZYX component_rot(x_component_rot.z(0), x_component_rot.y(0), x_component_rot.x(0));
0164         component_physvol =
0165             module_volume.placeVolume(component_volume, Transform3D(component_rot, component_pos));
0166       } else if (x_component_pos) {
0167         Position component_pos(x_component_pos.x(0), x_component_pos.y(0),
0168                                x_component_pos.z(0) + zoff);
0169         component_physvol = module_volume.placeVolume(component_volume, component_pos);
0170       } else {
0171         Position component_pos(0, 0, zoff);
0172         component_physvol = module_volume.placeVolume(component_volume, component_pos);
0173       }
0174
0175       if (x_component.isSensitive()) {
0176         component_physvol.addPhysVolID("sensor", sensor_number++);
0177         component_volume.setSensitiveDetector(sens);
0178         sensitives[module_name].push_back(component_physvol);
0179         module_thicknesses[module_name] = {thickness_so_far + x_component.thickness() / 2.0,
0180                                            total_thickness - thickness_so_far -
0181                                                x_component.thickness() / 2.0};
0182       }
0183
0184       thickness_sum += x_component.thickness();
0185       thickness_so_far += x_component.thickness();
0186
0187       // apply relative offsets in z-position used to stack components side-by-side
0188       if (x_component_pos) {
0189         thickness_sum += x_component_pos.z(0);
0190         thickness_so_far += x_component_pos.z(0);
0191       }
0192     }
0193   }
0194   if (x_detector.hasChild(_Unicode(support))) {
0195     buildSupport(desc, sector_volume, x_detector.child(_Unicode(support)),
0196                  {inner_r, layer_pos_z, x_dimensions.z(), half_dphi});
0197   }
0198   //Loop over the sets of layer elements in the detector.
0199   int layer_num = 1;
0200   for (xml_coll_t i_layer(x_detector, _U(layer)); i_layer; ++i_layer) {
0201     xml_comp_t x_layer         = i_layer;
0202     int layer_repeat           = x_layer.repeat();
0203     double layer_thickness     = x_layer.thickness();
0204     bool layer_has_frame       = x_layer.hasChild(_Unicode(frame));
0205     double layer_space_between = getAttrOrDefault(x_layer, _Unicode(space_between), 0.);
0206     double layer_space_before  = getAttrOrDefault(x_layer, _Unicode(space_before), 0.);
0207     layer_pos_z += layer_space_before;
0208
0209     // Loop over number of repeats for this layer.
0210     for (int layer_j = 0; layer_j < layer_repeat; layer_j++) {
0211
0212       // Make an envelope for this layer
0213       std::string layer_name = Form("layer%d", layer_num);
0214       double layer_dim_x     = tan_half_dphi * layer_pos_z;
0215       auto layer_mat         = desc.air();
0216
0217       Position layer_pos(0, 0, layer_pos_z + layer_thickness / 2.);
0218       double layer_trd_x1 = layer_dim_x;
0219       double layer_trd_x2 = layer_dim_x + layer_thickness * tan_half_dphi;
0220       double layer_trd_y1 = layer_dim_y;
0221       double layer_trd_y2 = layer_trd_y1;
0222       double layer_trd_z  = layer_thickness / 2; // account for frame
0223       Trapezoid layer_shape(layer_trd_x1, layer_trd_x2, layer_trd_y1, layer_trd_y2, layer_trd_z);
0224       Volume layer_volume(layer_name, layer_shape, layer_mat);
0225       DetElement layer_element(sector_element, layer_name, detector_id);
0226
0227       // Set region, limitset, and vis of layer.
0228       layer_volume.setAttributes(desc, x_layer.regionStr(), x_layer.limitsStr(), x_layer.visStr());
0229
0230       // Loop over the staves for this layer, if the tray is enabled
0231       int stave_num = 1;
0232       for (xml_coll_t i_stave(x_layer, _U(stave)); i_stave; ++i_stave) {
0233         xml_comp_t x_stave = i_stave;
0234         // Check if we enabled this silicon tray (actual silicon detector in the slots) at the top level
0235         // allowing us to easily disable layers in the XML file without needing multiple
0236         // copies of the XML file
0237         if (getAttrOrDefault(x_stave, _Unicode(enable), 1) == 0) {
0238           // disabled
0239           continue;
0240         }
0241
0242         int stave_repeat   = x_stave.repeat();
0243         double stave_thick = x_stave.thickness();
0244         double stave_dim_x = x_stave.width() / 2.0;
0245         double stave_dim_y = x_stave.length() / 2.0;
0246         double stave_dim_z = stave_thick / 2.0;
0247         double stave_rot_y = getAttrOrDefault(x_stave, _Unicode(angle), 0.);
0248         double stave_off_z = getAttrOrDefault(x_stave, _Unicode(offset), 0.);
0249
0250         // Arrange staves symmetrically around center of layer
0251         double stave_pos_x = -layer_dim_x + stave_dim_x;
0252         double stave_pitch = -2.0 * stave_pos_x / (stave_repeat - 1);
0253
0254         // Make one stave
0255         std::string stave_name = Form("stave%d", stave_num);
0256         auto stave_material    = desc.air();
0257         Box stave_shape(stave_dim_x, stave_dim_y, stave_dim_z);
0258         Volume stave_volume(stave_name, stave_shape, stave_material);
0259         stave_volume.setAttributes(desc, x_stave.regionStr(), x_stave.limitsStr(),
0260                                    x_stave.visStr());
0261         DetElement stave_element(layer_element, stave_name, detector_id);
0262
0263         // Loop over the slices for this stave
0264         double slice_pos_z = -(stave_thick / 2.);
0265
0266         // Place in xy_layout
0267         if (x_stave.hasChild(_Unicode(xy_layout))) {
0268           auto module_str         = x_stave.moduleStr();
0269           auto& module_volume     = volumes[module_str];
0270           auto& module_sensitives = sensitives[module_str];
0271
0272           // Get layout grid pitch
0273           xml_comp_t x_xy_layout = x_stave.child(_Unicode(xy_layout));
0274           auto dx                = x_xy_layout.attr<double>(_Unicode(dx));
0275           auto dy                = x_xy_layout.attr<double>(_Unicode(dy));
0276
0277           // Default to filling
0278           auto nx = getAttrOrDefault<int>(x_xy_layout, _Unicode(nx), floor(2. * stave_dim_x / dx));
0279           auto ny = getAttrOrDefault<int>(x_xy_layout, _Unicode(ny), floor(2. * stave_dim_y / dy));
0280           printout(DEBUG, "BarrelCalorimeterImaging", "Stave %s layout with %d by %d modules",
0281                    stave_name.c_str(), nx, ny);
0282
0283           // Default to centered
0284           auto x0 = getAttrOrDefault<double>(x_xy_layout, _Unicode(x0), -(nx - 1) * dx / 2.);
0285           auto y0 = getAttrOrDefault<double>(x_xy_layout, _Unicode(x0), -(ny - 1) * dy / 2.);
0286           printout(DEBUG, "BarrelCalorimeterImaging", "Stave %s modules starting at x=%f, y=%f",
0287                    stave_name.c_str(), x0, y0);
0288
0289           // Place modules
0290           int i_module = 0;
0291           for (auto i_x = 0; i_x < nx; ++i_x) {
0292             for (auto i_y = 0; i_y < ny; ++i_y) {
0293
0294               // Create module
0295               std::string module_name = _toString(i_module, "module%d");
0296               DetElement module_element(stave_element, module_name, i_module);
0297
0298               // Place module
0299               auto x = x0 + i_x * dx;
0300               auto y = y0 + i_y * dy;
0301               Position module_pos(x, y, 0);
0302               PlacedVolume module_physvol = stave_volume.placeVolume(module_volume, module_pos);
0303               module_physvol.addPhysVolID("module", i_module);
0304               module_element.setPlacement(module_physvol);
0305
0306               // Add sensitive volumes
0307               for (auto& sensitive_physvol : module_sensitives) {
0308                 DetElement sensitive_element(module_element, sensitive_physvol.volume().name(),
0309                                              i_module);
0310                 sensitive_element.setPlacement(sensitive_physvol);
0311
0312                 auto& sensitive_element_params =
0313                     DD4hepDetectorHelper::ensureExtension<dd4hep::rec::VariantParameters>(
0314                         sensitive_element);
0315                 sensitive_element_params.set<std::string>("axis_definitions", "XYZ");
0316               }
0317
0318               i_module++;
0319             }
0320           }
0321           slice_pos_z += module_thicknesses[module_str][0] + module_thicknesses[module_str][1];
0322         }
0323
0324         // Loop over the slices for this stave
0325         int slice_num = 1;
0326         for (xml_coll_t i_slice(x_stave, _U(slice)); i_slice; ++i_slice) {
0327           xml_comp_t x_slice     = i_slice;
0328           std::string slice_name = Form("slice%d", slice_num);
0329           double slice_thick     = x_slice.thickness();
0330           double slice_dim_x     = stave_dim_x;
0331           double slice_dim_y     = stave_dim_y;
0332           double slice_dim_z     = slice_thick / 2.;
0333           Box slice_shape(slice_dim_x, slice_dim_y, slice_dim_z);
0334           Volume slice_volume(slice_name, slice_shape, desc.material(x_slice.materialStr()));
0335           slice_volume.setAttributes(desc, x_slice.regionStr(), x_slice.limitsStr(),
0336                                      x_slice.visStr());
0337           DetElement slice_element(stave_element, slice_name, detector_id);
0338
0339           // Set sensitive
0340           if (x_slice.isSensitive()) {
0341             slice_volume.setSensitiveDetector(sens);
0342           }
0343
0344           // Place slice
0345           PlacedVolume slice_physvol =
0346               stave_volume.placeVolume(slice_volume, Position(0, 0, slice_pos_z + slice_thick / 2));
0347           slice_physvol.addPhysVolID("slice", slice_num);
0348           slice_element.setPlacement(slice_physvol);
0349
0350           // Increment Z position of slice
0351           slice_pos_z += slice_thick;
0352           ++slice_num;
0353         }
0354
0355         // Loop over number of repeats for this stave
0356         for (int stave_j = 0; stave_j < stave_repeat; stave_j++) {
0357
0358           // Stave placement
0359           Position stave_pos(stave_pos_x, 0, layer_j % 2 == 0 ? +stave_off_z : -stave_off_z);
0360           RotationY stave_rot(layer_j % 2 == 0 ? +stave_rot_y : -stave_rot_y);
0361           Transform3D stave_tr(stave_rot, stave_pos);
0362           PlacedVolume stave_physvol = layer_volume.placeVolume(stave_volume, stave_tr);
0363           stave_physvol.addPhysVolID("stave", stave_num);
0364           stave_element.setPlacement(stave_physvol);
0365
0366           // Increment X position of stave
0367           stave_pos_x += stave_pitch;
0368           ++stave_num;
0369         }
0370       }
0371
0372       // Place frame if defined.
0373       if (layer_has_frame) {
0374         xml_comp_t x_frame     = x_layer.child(_Unicode(frame));
0375         double frame_height    = x_frame.height();
0376         double frame_thickness = x_frame.thickness();
0377         auto frame_material    = desc.material(x_frame.materialStr());
0378
0379         std::string frame1_name = Form("frame_inner%d", layer_num);
0380         double frame1_thick     = frame_thickness;
0381         double frame1_trd_x1 =
0382             layer_dim_x + (layer_thickness / 2 - frame_height / 2) * tan_half_dphi;
0383         double frame1_trd_x2 =
0384             layer_dim_x + (layer_thickness / 2 - frame_height / 2 + frame1_thick) * tan_half_dphi;
0385         double frame1_trd_y1 = layer_trd_y1;
0386         double frame1_trd_y2 = frame1_trd_y1;
0387         double frame1_trd_z  = frame1_thick / 2.;
0388         Trapezoid frame1_shape(frame1_trd_x1, frame1_trd_x2, frame1_trd_y1, frame1_trd_y2,
0389                                frame1_trd_z);
0390         Volume frame1_volume(frame1_name, frame1_shape, frame_material);
0391         layer_volume.placeVolume(frame1_volume,
0392                                  Position(0, 0, -frame_height / 2.0 + frame_thickness / 2.0));
0393
0394         std::string frame2_name = Form("frame_outer%d", layer_num);
0395         double frame2_thick     = frame_thickness;
0396         double frame2_trd_x1 =
0397             layer_dim_x + (layer_thickness / 2 - frame_height / 2) * tan_half_dphi;
0398         double frame2_trd_x2 =
0399             layer_dim_x + (layer_thickness / 2 - frame_height / 2 + frame2_thick) * tan_half_dphi;
0400         double frame2_trd_y1 = layer_trd_y1;
0401         double frame2_trd_y2 = frame2_trd_y1;
0402         double frame2_trd_z  = frame2_thick / 2.;
0403         Trapezoid frame2_shape(frame2_trd_x1, frame2_trd_x2, frame2_trd_y1, frame2_trd_y2,
0404                                frame2_trd_z);
0405         Volume frame2_volume(frame2_name, frame2_shape, frame_material);
0406         layer_volume.placeVolume(frame2_volume,
0407                                  Position(0, 0, +frame_height / 2.0 - frame_thickness / 2.0));
0408       }
0409
0410       // Place layer into sector
0411       PlacedVolume layer_physvol = sector_volume.placeVolume(layer_volume, layer_pos);
0412       layer_physvol.addPhysVolID("layer", layer_num);
0413       layer_element.setPlacement(layer_physvol);
0414
0415       // Increment to next layer Z position. Do not add space_between for the last layer
0416       layer_pos_z += layer_thickness;
0417       if (layer_j < layer_repeat - 1) {
0418         layer_pos_z += layer_space_between;
0419       }
0420       ++layer_num;
0421     }
0422   }
0423
0424   // Phi start for a sector.
0425   double phi = M_PI / nsides;
0426   // Create nsides sectors.
0427   for (int i = 0; i < nsides; i++, phi -= dphi) { // i is sector number
0428     // Compute the sector position
0429     Transform3D tr(RotationZYX(0, phi, M_PI * 0.5), Translation3D(0, 0, 0));
0430     PlacedVolume sector_physvol = detector_volume.placeVolume(sector_volume, tr);
0431     sector_physvol.addPhysVolID("sector", i + 1);
0432     DetElement sd = (i == 0) ? sector_element : sector_element.clone(Form("sector%d", i));
0433     sd.setPlacement(sector_physvol);
0434     sdet.add(sd);
0435   }
0436
0437   return sdet;
0438 }
0439
0440 // simple aluminum sheet cover
0441 // dimensions: (inner r, position in z, length, phi)
0442 static void buildSupport(Detector& desc, Volume& sector_volume, xml_comp_t x_support,
0443                          const std::tuple<double, double, double, double>& dimensions) {
0444   auto [inner_r, pos_z, sector_length, hphi] = dimensions;
0445   double support_thickness = getAttrOrDefault(x_support, _Unicode(thickness), 0.5 * cm);
0446   auto material            = desc.material(x_support.materialStr());
0447   double trd_y             = sector_length / 2.;
0448   double trd_x1_support    = std::tan(hphi) * pos_z;
0449   double trd_x2_support    = std::tan(hphi) * (pos_z + support_thickness);
0450
0451   Trapezoid s_shape(trd_x1_support, trd_x2_support, trd_y, trd_y, support_thickness / 2.);
0452   Volume s_vol("support_layer", s_shape, material);
0453   s_vol.setVisAttributes(desc.visAttributes(x_support.visStr()));
0454   sector_volume.placeVolume(s_vol, Position(0.0, 0.0, pos_z + support_thickness / 2.));
0455 }
0456 DECLARE_DETELEMENT(epic_EcalBarrelImaging, create_detector)