Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-05-18 07:52:58 

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Sylvester Joosten
0003 
0004 /** \addtogroup Trackers Trackers
0005  * \brief Type: **BarrelTrackerWithFrame**.
0006  * \author W. Armstrong
0007  *
0008  * \ingroup trackers
0009  *
0010  * @{
0011  */
0012 #include <DD4hep/DetFactoryHelper.h>
0013 #include <DD4hep/Printout.h>
0014 #include <DD4hep/Shapes.h>
0015 #include <XML/Layering.h>
0016 #include <XML/Utilities.h>
0017 
0018 #include <cassert>
0019 
0020 using namespace std;
0021 using namespace dd4hep;
0022 
0023 namespace {
0024 std::pair<Volume, Transform3D> build_shape(const Detector& descr, const xml_det_t& x_det,
0025                                            const xml_comp_t& x_support, const xml_comp_t& x_child,
0026                                            const double offset = 0) {
0027   // Get Initial rotation/translation info
0028   xml_dim_t x_pos(x_child.child(_U(position), false));
0029   xml_dim_t x_rot(x_child.child(_U(rotation), false));
0030   Position pos3D{0, 0, 0};
0031   Rotation3D rot3D;
0032 
0033   if (x_rot) {
0034     rot3D = RotationZYX(x_rot.z(0), x_rot.y(0), x_rot.x(0));
0035   }
0036   if (x_pos) {
0037     pos3D = Position(x_pos.x(0), x_pos.y(0), x_pos.z(0));
0038   }
0039 
0040   // handle different known shapes and create solids
0041   Solid solid;
0042   const std::string type = x_support.attr<std::string>(_U(type));
0043   if (type == "Tube" || type == "Cylinder") {
0044     const double thickness = getAttrOrDefault(x_child, _U(thickness), x_support.thickness());
0045     const double length    = getAttrOrDefault(x_child, _U(length), x_support.length());
0046     const double rmin      = getAttrOrDefault(x_child, _U(rmin), x_support.rmin()) + offset;
0047     const double phimin    = getAttrOrDefault(
0048         x_child, _Unicode(phimin), getAttrOrDefault(x_support, _Unicode(phimin), 0.0 * deg));
0049     const double phimax = getAttrOrDefault(
0050         x_child, _Unicode(phimax), getAttrOrDefault(x_support, _Unicode(phimax), 360.0 * deg));
0051     solid = Tube(rmin, rmin + thickness, length / 2, phimin, phimax);
0052   } else if (type == "Box") {
0053     const double box_x = getAttrOrDefault(x_child, _U(x), x_support.x());
0054     const double box_y = getAttrOrDefault(x_child, _U(y), x_support.y());
0055     const double box_z = getAttrOrDefault(x_child, _U(z), x_support.z());
0056     solid              = Box(box_x / 2.0, box_y / 2.0, box_z / 2.0);
0057   }
0058   // A disk is a cylinder, constructed differently
0059   else if (type == "Disk") {
0060     const double thickness = getAttrOrDefault(x_child, _U(thickness), x_support.thickness());
0061     const double rmin      = getAttrOrDefault(x_child, _U(rmin), x_support.rmin());
0062     const double rmax      = getAttrOrDefault(x_child, _U(rmax), x_support.rmax());
0063     const double phimin    = getAttrOrDefault(
0064         x_child, _Unicode(phimin), getAttrOrDefault(x_support, _Unicode(phimin), 0.0 * deg));
0065     const double phimax = getAttrOrDefault(
0066         x_child, _Unicode(phimax), getAttrOrDefault(x_support, _Unicode(phimax), 360.0 * deg));
0067     pos3D = pos3D + Position(0, 0, -x_support.thickness() / 2 + thickness / 2 + offset);
0068     solid = Tube(rmin, rmax, thickness / 2, phimin, phimax);
0069   } else if (type == "Cone") {
0070     const double base_rmin1 = getAttrOrDefault(x_child, _U(rmin1), x_support.rmin1());
0071     const double base_rmin2 = getAttrOrDefault(x_child, _U(rmin2), x_support.rmin2());
0072     const double length     = getAttrOrDefault(x_child, _U(length), x_support.length());
0073     // Account for the fact that the distance between base_rmin1 and rmax2 is the projection
0074     // of the thickness on the transverse direction
0075     const double thickness = getAttrOrDefault(x_child, _U(thickness), x_support.thickness());
0076     const double transverse_thickness =
0077         thickness / cos(atan2(fabs(base_rmin2 - base_rmin1), length));
0078     // also account that the same is true for the offset
0079     const double transverse_offset = offset / cos(atan2(fabs(base_rmin2 - base_rmin1), length));
0080     const double rmin1             = base_rmin1 + transverse_offset;
0081     const double rmin2             = base_rmin2 + transverse_offset;
0082     const double rmax1             = rmin1 + transverse_thickness;
0083     const double rmax2             = rmin2 + transverse_thickness;
0084     // Allow for optional hard rmin/rmax cutoffs
0085     const double rmin = getAttrOrDefault(
0086         x_child, _U(rmin), getAttrOrDefault(x_support, _Unicode(rmin), min(rmin1, rmin2)));
0087     const double rmax = getAttrOrDefault(
0088         x_child, _U(rmax), getAttrOrDefault(x_support, _Unicode(rmax), max(rmax1, rmax2)));
0089     if (rmin > min(rmax1, rmax2)) {
0090       printout(ERROR, x_det.nameStr(),
0091                "%s: rmin (%f mm) must be smaller than the smallest rmax (%f %f mm)",
0092                x_support.nameStr().c_str(), rmin / mm, rmax1 / mm, rmax2 / mm);
0093       std::exit(1);
0094     }
0095     if (rmax < max(base_rmin1, base_rmin2)) {
0096       printout(ERROR, x_det.nameStr(),
0097                "%s: rmax (%f mm) must be larger than the largest rmin (%f %f mm)",
0098                x_support.nameStr().c_str(), rmax / mm, base_rmin1 / mm, base_rmin2 / mm);
0099       std::exit(1);
0100     }
0101     const double zmin  = -length / 2 + length * (rmin - rmin1) / (rmin2 - rmin1);
0102     const double zmax  = -length / 2 + length * (rmax - rmax1) / (rmax2 - rmax1);
0103     const auto rmin_at = [&](const double z) {
0104       return rmin1 + (z + length / 2) * (rmin2 - rmin1) / length;
0105     };
0106     const auto rmax_at = [&](const double z) {
0107       return rmax1 + (z + length / 2) * (rmax2 - rmax1) / length;
0108     };
0109     // Allow for optional phimin/phimax
0110     const double phimin = getAttrOrDefault<double>(
0111         x_child, _Unicode(phimin), getAttrOrDefault(x_support, _Unicode(phimin), 0.0 * deg));
0112     const double phimax = getAttrOrDefault<double>(
0113         x_child, _Unicode(phimax), getAttrOrDefault(x_support, _Unicode(phimax), 360.0 * deg));
0114     const double deltaphi = phimax - phimin;
0115     const double epsilon{TGeoShape::Tolerance()};
0116     if (fabs(zmin) >= length / 2 - epsilon && fabs(zmax) >= length / 2 - epsilon) {
0117       if (fabs(phimax - phimin - 360 * deg) < epsilon) {
0118         solid = Cone(length / 2, rmin1, rmax1, rmin2, rmax2);
0119       } else {
0120         solid = ConeSegment(length / 2, rmin1, rmax1, rmin2, rmax2, phimin, phimax);
0121       }
0122     } else {
0123       std::vector<double> v_rmin{max(rmin1, rmin), max(rmin2, rmin)},
0124           v_rmax{min(rmax1, rmax), min(rmax2, rmax)}, v_z{-length / 2, +length / 2};
0125       for (const auto& z :
0126            (zmin < zmax ? std::vector<double>{zmin, zmax} : std::vector<double>{zmax, zmin})) {
0127         if (-length / 2 + epsilon < z && z < -epsilon + length / 2) {
0128           v_rmin.insert(std::prev(v_rmin.end()), std::max(rmin, rmin_at(z)));
0129           v_rmax.insert(std::prev(v_rmax.end()), std::min(rmax, rmax_at(z)));
0130           v_z.insert(std::prev(v_z.end()), z);
0131         }
0132       }
0133       solid = Polycone(phimin, deltaphi, v_rmin, v_rmax, v_z);
0134     }
0135   } else if (type == "Disk") {
0136     const double thickness = getAttrOrDefault(x_child, _U(thickness), x_support.thickness());
0137     const double rmin      = getAttrOrDefault(x_child, _U(rmin), x_support.rmin());
0138     const double rmax      = getAttrOrDefault(x_child, _U(rmax), x_support.rmax());
0139     const double phimin    = getAttrOrDefault(
0140         x_child, _Unicode(phimin), getAttrOrDefault(x_support, _Unicode(phimin), 0.0 * deg));
0141     const double phimax = getAttrOrDefault(
0142         x_child, _Unicode(phimax), getAttrOrDefault(x_support, _Unicode(phimax), 360.0 * deg));
0143     pos3D = pos3D + Position(0, 0, -x_support.thickness() / 2 + thickness / 2 + offset);
0144     solid = Tube(rmin, rmax, thickness / 2, phimin, phimax);
0145   } else {
0146     printout(ERROR, x_det.nameStr(), "Unknown support type: %s", type.c_str());
0147     std::exit(1);
0148   }
0149   // Materials
0150   Material mat = descr.material(getAttrOrDefault<std::string>(x_child, _U(material), "Air"));
0151   // Create our volume
0152   Volume vol{getAttrOrDefault<std::string>(x_child, _U(name), "support_vol"), solid, mat};
0153 
0154   // Create full transformation
0155   Transform3D tr(rot3D, pos3D);
0156 
0157   // visualization?
0158   if (x_child.hasAttr(_U(vis))) {
0159     vol.setVisAttributes(descr.visAttributes(x_child.visStr()));
0160   }
0161   return {vol, tr};
0162 }
0163 
0164 // Function to create a subtraction of multiple shapes
0165 std::pair<Volume, Transform3D> build_subtraction(const Detector& descr, const xml_det_t& x_det,
0166                                                  const xml_comp_t& x_subtraction) {
0167 
0168   // Try generic multi-shape schema: <shape role="base"> + N x <shape role="hole">
0169   Solid baseSolid;
0170   bool baseSet           = false;
0171   bool haveGenericShapes = false;
0172   std::vector<std::pair<Solid, Transform3D>> holes;
0173 
0174   for (xml_coll_t s{x_subtraction, _U(shape)}; s; ++s) {
0175     haveGenericShapes = true;
0176     xml_comp_t x_s    = s;
0177 
0178     // Build solid for this child
0179     auto [v, tr_unused] = build_shape(descr, x_det, x_s, x_s);
0180     Solid sld           = v.solid();
0181 
0182     // Relative transform (default identity) for this child
0183     xml_dim_t x_pos(x_s.child(_U(position), false));
0184     xml_dim_t x_rot(x_s.child(_U(rotation), false));
0185     Position pos{0, 0, 0};
0186     Rotation3D rot;
0187     if (x_rot)
0188       rot = RotationZYX(x_rot.z(0), x_rot.y(0), x_rot.x(0));
0189     if (x_pos)
0190       pos = Position(x_pos.x(0), x_pos.y(0), x_pos.z(0));
0191     Transform3D tr_rel(rot, pos);
0192 
0193     const std::string role = getAttrOrDefault<std::string>(x_s, _Unicode(role), "hole");
0194     if (role == "base" || role == "primary") {
0195       baseSolid = sld;
0196       baseSet   = true;
0197     } else {
0198       holes.emplace_back(sld, tr_rel);
0199     }
0200   }
0201 
0202   // Backward compatibility: if no generic <shape> elements, use <shape1>/<shape2>
0203   if (!haveGenericShapes) {
0204     xml_comp_t x_shape1 = x_subtraction.child(_Unicode(shape1));
0205     xml_comp_t x_shape2 = x_subtraction.child(_Unicode(shape2));
0206     if (!x_shape1 || !x_shape2) {
0207       printout(ERROR, x_det.nameStr(),
0208                "Subtraction %s: expected shape1/shape2 or generic <shape> children.",
0209                x_subtraction.nameStr().c_str());
0210       std::exit(1);
0211     }
0212     auto [vol1, tr1_unused] = build_shape(descr, x_det, x_shape1, x_shape1);
0213     auto [vol2, tr2_unused] = build_shape(descr, x_det, x_shape2, x_shape2);
0214     baseSolid               = vol1.solid();
0215     baseSet                 = true;
0216 
0217     // Relative transform for shape2
0218     xml_dim_t x_pos2(x_shape2.child(_U(position), false));
0219     xml_dim_t x_rot2(x_shape2.child(_U(rotation), false));
0220     Position pos2{0, 0, 0};
0221     Rotation3D rot2;
0222     if (x_rot2)
0223       rot2 = RotationZYX(x_rot2.z(0), x_rot2.y(0), x_rot2.x(0));
0224     if (x_pos2)
0225       pos2 = Position(x_pos2.x(0), x_pos2.y(0), x_pos2.z(0));
0226     holes.emplace_back(vol2.solid(), Transform3D(rot2, pos2));
0227   }
0228 
0229   if (!baseSet) {
0230     printout(ERROR, x_det.nameStr(), "Subtraction %s: no base shape found (role=\"base\").",
0231              x_subtraction.nameStr().c_str());
0232     std::exit(1);
0233   }
0234 
0235   // Chain subtractions: base − hole1 − hole2 − ...
0236   Solid current = baseSolid;
0237   for (auto& h : holes) {
0238     current = SubtractionSolid(current, h.first, h.second);
0239   }
0240 
0241   // Create the resulting volume
0242   Material mat =
0243       descr.material(getAttrOrDefault<std::string>(x_subtraction, _Unicode(material), "Air"));
0244   std::string vol_name =
0245       getAttrOrDefault<std::string>(x_subtraction, _Unicode(name), "subtraction_vol");
0246   Volume vol{vol_name, current, mat};
0247 
0248   // Overall placement of the resulting volume
0249   xml_dim_t x_pos(x_subtraction.child(_U(position), false));
0250   xml_dim_t x_rot(x_subtraction.child(_U(rotation), false));
0251   Position pos3D{0, 0, 0};
0252   Rotation3D rot3D;
0253   if (x_rot)
0254     rot3D = RotationZYX(x_rot.z(0), x_rot.y(0), x_rot.x(0));
0255   if (x_pos)
0256     pos3D = Position(x_pos.x(0), x_pos.y(0), x_pos.z(0));
0257   Transform3D tr(rot3D, pos3D);
0258 
0259   // Set visualization if specified
0260   if (x_subtraction.hasAttr(_U(vis))) {
0261     vol.setVisAttributes(descr.visAttributes(x_subtraction.visStr()));
0262   }
0263 
0264   // Debug output
0265   printout(DEBUG, "SupportServiceMaterial", "Created subtraction volume: %s (holes=%d)",
0266            vol_name.c_str(), int(holes.size()));
0267 
0268   return {vol, tr};
0269 }
0270 
0271 std::pair<Volume, Transform3D> build_shape(const Detector& descr, const xml_det_t& x_det,
0272                                            const xml_comp_t& x_support, const double offset = 0) {
0273   return build_shape(descr, x_det, x_support, x_support, offset);
0274 }
0275 } // namespace
0276 
0277 /** Generic tracker support implementation, can consist of arbitrary shapes
0278  *
0279  * @author Sylvester Joosten
0280  */
0281 static Ref_t create_SupportServiceMaterial(Detector& description, xml_h e,
0282                                            [[maybe_unused]] SensitiveDetector sens) {
0283   const xml_det_t x_det = e;
0284   const int det_id      = x_det.id();
0285   const string det_name = x_det.nameStr();
0286 
0287   // global z-offset for the entire support assembly
0288   const double offset = getAttrOrDefault(x_det, _U(offset), 0.);
0289 
0290   DetElement det(det_name, det_id);
0291   Assembly assembly(det_name + "_assembly");
0292 
0293   // Loop over the supports
0294   for (xml_coll_t su{x_det, _U(support)}; su; ++su) {
0295     xml_comp_t x_sup   = su;
0296     const double rot_z = getAttrOrDefault(x_sup, _U(phi0), 0.);
0297     RotationZ rot(rot_z);
0298     Transform3D tr_rot(
0299         rot, Position(0, 0, 0)); // additional rotation of the module after position offset
0300 
0301     auto [vol, tr]           = build_shape(description, x_det, x_sup);
0302     [[maybe_unused]] auto pv = assembly.placeVolume(vol, tr_rot * tr);
0303     // Loop over support components, if any
0304     double cumulative_thickness = 0;
0305     for (xml_coll_t com{x_sup, _U(component)}; com; ++com) {
0306       xml_comp_t x_com = com;
0307       auto [cvol, ctr] = build_shape(description, x_det, x_sup, x_com, cumulative_thickness);
0308       vol.placeVolume(cvol, ctr);
0309       cumulative_thickness += x_com.thickness();
0310     }
0311   }
0312   // Loop over any subtraction volumes
0313   for (xml_coll_t sub{x_det, _Unicode(subtraction)}; sub; ++sub) {
0314     xml_comp_t x_sub = sub;
0315     auto [vol, tr]   = build_subtraction(description, x_det, x_sub);
0316     assembly.placeVolume(vol, tr);
0317   }
0318 
0319   // final placement
0320   Volume motherVol = description.pickMotherVolume(det);
0321   Position pos(0, 0, offset);
0322   PlacedVolume pv = motherVol.placeVolume(assembly, pos);
0323   pv.addPhysVolID("system", det.id());
0324   det.setPlacement(pv);
0325 
0326   return det;
0327 }
0328 
0329 // clang-format off
0330 DECLARE_DETELEMENT(epic_SupportServiceMaterial, create_SupportServiceMaterial)
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!