epic/src/EndcapTOF_geo.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 - 2024, Nicolas Schmidt, Chun Yuen Tsang
0003
0004 /** \addtogroup Trackers Trackers
0005  * \brief Type: **Endcap Tracker with TOF**.
0006  *
0007  * \ingroup trackers
0008  *
0009  * @{
0010  */
0011 #include "DD4hep/DetFactoryHelper.h"
0012 #include "DD4hep/Printout.h"
0013 #include "DD4hep/Shapes.h"
0014 #include "DD4hepDetectorHelper.h"
0015 #include "DDRec/DetectorData.h"
0016 #include "DDRec/Surface.h"
0017 #include "XML/Layering.h"
0018 #include "XML/Utilities.h"
0019 #include <array>
0020 #include <map>
0021 #include <unordered_set>
0022
0023 using namespace std;
0024 using namespace dd4hep;
0025 using namespace dd4hep::rec;
0026 using namespace dd4hep::detail;
0027
0028 static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector sens) {
0029   typedef vector<PlacedVolume> Placements;
0030   xml_det_t x_det      = e;
0031   int det_id           = x_det.id();
0032   std::string det_name = x_det.nameStr();
0033   DetElement sdet(det_name, det_id);
0034   Material air = description.material("Air");
0035   PlacedVolume pv;
0036
0037   map<string, std::array<double, 2>> module_thicknesses;
0038
0039   // Set detector type flag
0040   dd4hep::xml::setDetectorTypeFlag(x_det, sdet);
0041   auto& params = DD4hepDetectorHelper::ensureExtension<dd4hep::rec::VariantParameters>(sdet);
0042   // Add the volume boundary material if configured
0043   for (xml_coll_t bmat(x_det, _Unicode(boundary_material)); bmat; ++bmat) {
0044     xml_comp_t x_boundary_material = bmat;
0045     DD4hepDetectorHelper::xmlToProtoSurfaceMaterial(x_boundary_material, params,
0046                                                     "boundary_material");
0047   }
0048
0049   Assembly assembly(det_name);
0050   assembly.setVisAttributes(description.invisible());
0051   sens.setType("tracker");
0052
0053   // dimensions of the modules (2x2 sensors)
0054   xml_comp_t x_modsz = x_det.child(_Unicode(modsize));
0055
0056   double module_x       = x_modsz.length();
0057   double module_y       = x_modsz.width();
0058   double module_overlap = getAttrOrDefault(x_modsz, _Unicode(overlap), 0.); // x_modsz.overlap();
0059   double module_spacing = getAttrOrDefault(x_modsz, _Unicode(spacing), 0.); // x_modsz.overlap();
0060   double board_gap      = getAttrOrDefault(x_modsz, _Unicode(board_gap), 0.);
0061
0062   map<string, Assembly> modules;
0063   map<string, Placements> sensitives;
0064   map<string, std::vector<VolPlane>> volplane_surfaces;
0065   map<string, double> mod_thickness;
0066
0067   for (xml_coll_t mi(x_det, _U(module)); mi; ++mi) {
0068     xml_comp_t x_mod       = mi;
0069     double total_thickness = 0;
0070     // Compute module total thickness from components
0071     for (xml_coll_t ci(x_mod, _U(module_component)); ci; ++ci) {
0072       xml_comp_t x_comp    = ci;
0073       bool keep_same_layer = getAttrOrDefault<bool>(x_comp, _Unicode(keep_layer), false);
0074       if (!keep_same_layer)
0075         total_thickness += x_comp.thickness();
0076     }
0077
0078     double thickness_so_far = 0.0;
0079     int sensitive_id        = 0;
0080     const std::string m_nam = x_mod.nameStr();
0081     mod_thickness[m_nam]    = total_thickness;
0082
0083     // the module assembly volume
0084     Assembly m_vol(m_nam);
0085     m_vol.setVisAttributes(description.visAttributes(x_mod.visStr()));
0086
0087     int ncomponents = 0;
0088     for (xml_coll_t mci(x_mod, _U(module_component)); mci; ++mci, ++ncomponents) {
0089       xml_comp_t x_comp  = mci;
0090       xml_comp_t x_pos   = x_comp.position(false);
0091       xml_comp_t x_rot   = x_comp.rotation(false);
0092       const string c_nam = x_comp.nameStr();
0093       bool cylindrical   = getAttrOrDefault<bool>(x_comp, _Unicode(cylindrical), false);
0094
0095       Volume c_vol;
0096       if (cylindrical) {
0097         Tube c_tub(x_comp.rmin(), x_comp.rmax(), x_comp.thickness() / 2.0, 0, 2 * M_PI);
0098         c_vol = Volume(c_nam, c_tub, description.material(x_comp.materialStr()));
0099       } else {
0100         Box c_box(x_comp.width() / 2, x_comp.length() / 2, x_comp.thickness() / 2);
0101         c_vol = Volume(c_nam, c_box, description.material(x_comp.materialStr()));
0102       }
0103       // Utility variable for the relative z-offset based off the previous components
0104       const double zoff = thickness_so_far + x_comp.thickness() / 2.0;
0105       if (x_pos && x_rot) {
0106         Position c_pos(x_pos.x(0), x_pos.y(0), x_pos.z(0) + zoff);
0107         RotationZYX c_rot(x_rot.z(0), x_rot.y(0), x_rot.x(0));
0108         pv = m_vol.placeVolume(c_vol, Transform3D(c_rot, c_pos));
0109       } else if (x_rot) {
0110         Position c_pos(0, 0, zoff);
0111         pv = m_vol.placeVolume(c_vol,
0112                                Transform3D(RotationZYX(x_rot.z(0), x_rot.y(0), x_rot.x(0)), c_pos));
0113       } else if (x_pos) {
0114         pv = m_vol.placeVolume(c_vol, Position(x_pos.x(0), x_pos.y(0), x_pos.z(0) + zoff));
0115       } else {
0116         pv = m_vol.placeVolume(c_vol, Position(0, 0, zoff));
0117       }
0118       c_vol.setRegion(description, x_comp.regionStr());
0119       c_vol.setLimitSet(description, x_comp.limitsStr());
0120       c_vol.setVisAttributes(description, x_comp.visStr());
0121       if (x_comp.isSensitive()) {
0122         pv.addPhysVolID("ids", sensitive_id);
0123         ++sensitive_id;
0124
0125         c_vol.setSensitiveDetector(sens);
0126         module_thicknesses[m_nam] = {thickness_so_far + x_comp.thickness() / 2.0,
0127                                      total_thickness - thickness_so_far - x_comp.thickness() / 2.0};
0128         sensitives[m_nam].push_back(pv);
0129
0130         // -------- create a measurement plane for the tracking surface attched to the sensitive volume -----
0131         Vector3D u(-1., 0., 0.);
0132         Vector3D v(0., -1., 0.);
0133         Vector3D n(0., 0., 1.);
0134
0135         // compute the inner and outer thicknesses that need to be assigned to the tracking surface
0136         // depending on wether the support is above or below the sensor
0137         double inner_thickness = module_thicknesses[m_nam][0];
0138         double outer_thickness = module_thicknesses[m_nam][1];
0139
0140         SurfaceType type(SurfaceType::Sensitive);
0141
0142         VolPlane surf(c_vol, type, inner_thickness, outer_thickness, u, v, n);
0143         volplane_surfaces[m_nam].push_back(surf);
0144
0145         //--------------------------------------------
0146       }
0147       bool keep_same_layer = getAttrOrDefault<bool>(x_comp, _Unicode(keep_layer), false);
0148       if (!keep_same_layer) {
0149         thickness_so_far += x_comp.thickness();
0150         // apply relative offsets in z-position used to stack components side-by-side
0151         if (x_pos) {
0152           thickness_so_far += x_pos.z(0);
0153         }
0154       }
0155     }
0156     modules[m_nam] = m_vol;
0157   }
0158
0159   int module       = 0;
0160   double current_z = std::numeric_limits<double>::lowest();
0161   for (xml_coll_t li(x_det, _U(layer)); li; ++li) {
0162     xml_comp_t x_layer = li;
0163     bool front         = x_layer.attr<bool>(_Unicode(front));
0164
0165     const std::string locStr = x_layer.nameStr();
0166
0167     // now build the envelope for the detector
0168     xml_comp_t envelope = x_layer.child(_Unicode(envelope), false);
0169     int lay_id          = x_layer.id();
0170     string lay_nam      = det_name + _toString(lay_id, "_layer%d");
0171     double phimin       = dd4hep::getAttrOrDefault<double>(envelope, _Unicode(phimin), 0.);
0172     double phimax       = dd4hep::getAttrOrDefault<double>(envelope, _Unicode(phimax), 2 * M_PI);
0173     double xoffset      = getAttrOrDefault<double>(envelope, _Unicode(xoffset), 0);
0174     bool zstack         = getAttrOrDefault<bool>(envelope, _Unicode(zstack), false);
0175
0176     // envelope thickness: use explicit thickness attribute if present, otherwise
0177     // compute as the max module thickness across all layouts in this layer
0178     double envelope_length = getAttrOrDefault<double>(envelope, _Unicode(thickness), 0.0);
0179     if (envelope_length == 0.0) {
0180       for (xml_coll_t llayout(x_layer, _Unicode(layout)); llayout; ++llayout) {
0181         xml_comp_t x_layout = llayout;
0182         string m_nam        = x_layout.moduleStr();
0183         envelope_length     = std::max(envelope_length, mod_thickness[m_nam]);
0184       }
0185     }
0186
0187     double zstart = 0;
0188     if (zstack) {
0189       if (current_z == std::numeric_limits<double>::lowest())
0190         throw std::runtime_error(
0191             "You cannot enable stack on the first layer. You need to place the first layer with "
0192             "'zstart', then you can stack the next layer.");
0193       else
0194         zstart = current_z;
0195     } else
0196       zstart = envelope.zstart();
0197
0198     Tube lay_tub(envelope.rmin(), envelope.rmax(), envelope_length / 2.0, phimin, phimax);
0199     Volume lay_vol(lay_nam, lay_tub, air); // Create the layer envelope volume.
0200     Position lay_pos(xoffset, 0, zstart + 0.5 * envelope_length);
0201     current_z = zstart + envelope_length;
0202     lay_vol.setVisAttributes(description.visAttributes(x_layer.visStr()));
0203
0204     DetElement lay_elt(sdet, lay_nam, lay_id);
0205     // Create the PhysicalVolume for the layer.
0206     pv = assembly.placeVolume(lay_vol, lay_pos); // Place layer in mother
0207     pv.addPhysVolID("layer", lay_id);            // Set the layer ID.
0208     lay_elt.setAttributes(description, lay_vol, x_layer.regionStr(), x_layer.limitsStr(),
0209                           x_layer.visStr());
0210     lay_elt.setPlacement(pv);
0211
0212     // the local coordinate systems of modules in dd4hep and acts differ
0213     // see http://acts.web.cern.ch/ACTS/latest/doc/group__DD4hepPlugins.html
0214     auto& layerParams =
0215         DD4hepDetectorHelper::ensureExtension<dd4hep::rec::VariantParameters>(lay_elt);
0216
0217     for (xml_coll_t lmat(x_layer, _Unicode(layer_material)); lmat; ++lmat) {
0218       xml_comp_t x_layer_material = lmat;
0219       DD4hepDetectorHelper::xmlToProtoSurfaceMaterial(x_layer_material, layerParams,
0220                                                       "layer_material");
0221     }
0222
0223     for (xml_coll_t llayout(x_layer, _Unicode(layout)); llayout; ++llayout, ++module) {
0224       xml_comp_t x_layout    = llayout;
0225       bool left              = x_layout.attr<bool>(_Unicode(left));
0226       string m_nam           = x_layout.moduleStr();
0227       Volume m_vol           = modules[m_nam];
0228       double total_thickness = mod_thickness[m_nam];
0229       Placements& sensVols   = sensitives[m_nam];
0230
0231       float ycoord = envelope.rmax() -
0232                      module_y / 2.; // y-center-coord of the top sensor. Start from the top row
0233       int iy       = 0;
0234
0235       for (xml_coll_t lrow(x_layout, _Unicode(row)); lrow; ++lrow) {
0236         xml_comp_t x_row = lrow;
0237         double deadspace = getAttrOrDefault<double>(x_row, _Unicode(deadspace), 0);
0238         if (deadspace > 0) {
0239           ycoord -= deadspace;
0240           continue;
0241         }
0242         double x_offset = getAttrOrDefault<double>(x_row, _Unicode(x_offset), 0);
0243         int nsensors    = getAttrOrDefault<int>(x_row, _Unicode(nsensors), 0);
0244
0245         // find the sensor id that corrsponds to the rightmost sensor in a board
0246         // we need to know where to apply additional spaces between neighboring board
0247         std::unordered_set<int> sensors_id_board_edge;
0248         int curr_ix = nsensors; // the first sensor to the right of center has ix of nsensors
0249         for (xml_coll_t lboard(x_row, _Unicode(board)); lboard; ++lboard) {
0250           xml_comp_t x_board = lboard;
0251           int nboard_sensors = getAttrOrDefault<int>(x_board, _Unicode(nsensors), 1);
0252           curr_ix += nboard_sensors;
0253           sensors_id_board_edge.insert(curr_ix);
0254           sensors_id_board_edge.insert(2 * nsensors - curr_ix -
0255                                        1); // reflected to sensor id on the left
0256         }
0257
0258         double accum_xoffset = x_offset;
0259
0260         for (int ix = (left ? nsensors - 1 : nsensors); (ix >= 0) && (ix < 2 * nsensors);
0261              ix     = ix + (left ? -1 : 1)) {
0262           // add board spacing
0263           if (sensors_id_board_edge.find(ix) != sensors_id_board_edge.end())
0264             accum_xoffset = accum_xoffset + board_gap;
0265
0266           // there is a hole in the middle, with radius = x_offset
0267           float xcoord = (ix - nsensors + 0.5) * (module_x + module_spacing) +
0268                          +(left ? -accum_xoffset : accum_xoffset);
0269           //! Note the module ordering is different for front and back side
0270
0271           double module_z = -0.5 * envelope_length;
0272           if (front)
0273             module_z = 0.5 * envelope_length - total_thickness;
0274
0275           // module built!
0276           Transform3D tr(RotationZYX(M_PI / 2, 0, 0), Position(xcoord, ycoord, module_z));
0277
0278           pv = lay_vol.placeVolume(m_vol, tr);
0279           pv.addPhysVolID("idx", ix);
0280           pv.addPhysVolID("idy", iy);
0281           pv.addPhysVolID("module", module);
0282
0283           string comp_nam = Form("%s_%s_ix%d_iy%d", lay_nam.c_str(), m_nam.c_str(), ix, iy);
0284           DetElement comp_elt(lay_elt, comp_nam, det_id);
0285           comp_elt.setPlacement(pv);
0286
0287           for (size_t ic = 0; ic < sensVols.size(); ++ic) {
0288             PlacedVolume sens_pv = sensVols[ic];
0289             DetElement sensor_elt(comp_elt, sens_pv.volume().name(), module);
0290             sensor_elt.setPlacement(sens_pv);
0291             auto& sensor_elt_params =
0292                 DD4hepDetectorHelper::ensureExtension<dd4hep::rec::VariantParameters>(sensor_elt);
0293             sensor_elt_params.set<string>("axis_definitions", "XYZ");
0294             volSurfaceList(sensor_elt)->push_back(volplane_surfaces[m_nam][ic]);
0295           }
0296         }
0297         ycoord -= (module_y - module_overlap);
0298         ++iy;
0299       }
0300     }
0301   }
0302   pv = description.pickMotherVolume(sdet).placeVolume(assembly, Position(0, 0, 0));
0303   pv.addPhysVolID("system", det_id);
0304   sdet.setPlacement(pv);
0305
0306   return sdet;
0307 }
0308
0309 DECLARE_DETELEMENT(epic_TOFEndcap, create_detector)