Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-10-17 08:04:12 

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 two 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   // Get the two shapes to subtract
0168   xml_comp_t x_shape1 = x_subtraction.child(_Unicode(shape1));
0169   xml_comp_t x_shape2 = x_subtraction.child(_Unicode(shape2));
0170 
0171   // Build the primary shape (without transformation)
0172   auto [vol1, tr1_unused] = build_shape(descr, x_det, x_shape1, x_shape1);
0173   Solid solid1            = vol1.solid();
0174 
0175   // Build the shape to subtract (without transformation)
0176   auto [vol2, tr2_unused] = build_shape(descr, x_det, x_shape2, x_shape2);
0177   Solid solid2            = vol2.solid();
0178 
0179   // Get relative transformation for shape2 (relative to shape1's center)
0180   xml_dim_t x_pos2(x_shape2.child(_U(position), false));
0181   xml_dim_t x_rot2(x_shape2.child(_U(rotation), false));
0182   Position pos2{0, 0, 0};
0183   Rotation3D rot2;
0184 
0185   if (x_rot2) {
0186     rot2 = RotationZYX(x_rot2.z(0), x_rot2.y(0), x_rot2.x(0));
0187   }
0188   if (x_pos2) {
0189     pos2 = Position(x_pos2.x(0), x_pos2.y(0), x_pos2.z(0));
0190   }
0191 
0192   // Create the subtraction solid
0193   Transform3D tr_relative(rot2, pos2);
0194   SubtractionSolid subtracted_solid(solid1, solid2, tr_relative);
0195 
0196   // Create volume with the subtracted solid
0197   Material mat = descr.material(getAttrOrDefault<std::string>(x_subtraction, _U(material), "Air"));
0198   std::string vol_name = getAttrOrDefault<std::string>(x_subtraction, _U(name), "subtraction_vol");
0199   Volume vol{vol_name, subtracted_solid, mat};
0200 
0201   // Get overall position/rotation for the subtracted volume
0202   xml_dim_t x_pos(x_subtraction.child(_U(position), false));
0203   xml_dim_t x_rot(x_subtraction.child(_U(rotation), false));
0204   Position pos3D{0, 0, 0};
0205   Rotation3D rot3D;
0206 
0207   if (x_rot) {
0208     rot3D = RotationZYX(x_rot.z(0), x_rot.y(0), x_rot.x(0));
0209   }
0210   if (x_pos) {
0211     pos3D = Position(x_pos.x(0), x_pos.y(0), x_pos.z(0));
0212   }
0213 
0214   Transform3D tr(rot3D, pos3D);
0215 
0216   // Set visualization if specified
0217   if (x_subtraction.hasAttr(_U(vis))) {
0218     vol.setVisAttributes(descr.visAttributes(x_subtraction.visStr()));
0219   }
0220 
0221   // Debug output
0222   printout(DEBUG, "SupportServiceMaterial", "Created subtraction volume: %s", vol_name.c_str());
0223 
0224   return {vol, tr};
0225 }
0226 
0227 std::pair<Volume, Transform3D> build_shape(const Detector& descr, const xml_det_t& x_det,
0228                                            const xml_comp_t& x_support, const double offset = 0) {
0229   return build_shape(descr, x_det, x_support, x_support, offset);
0230 }
0231 } // namespace
0232 
0233 /** Generic tracker support implementation, can consist of arbitrary shapes
0234  *
0235  * @author Sylvester Joosten
0236  */
0237 static Ref_t create_SupportServiceMaterial(Detector& description, xml_h e,
0238                                            [[maybe_unused]] SensitiveDetector sens) {
0239   const xml_det_t x_det = e;
0240   const int det_id      = x_det.id();
0241   const string det_name = x_det.nameStr();
0242 
0243   // global z-offset for the entire support assembly
0244   const double offset = getAttrOrDefault(x_det, _U(offset), 0.);
0245 
0246   DetElement det(det_name, det_id);
0247   Assembly assembly(det_name + "_assembly");
0248 
0249   // Loop over the supports
0250   for (xml_coll_t su{x_det, _U(support)}; su; ++su) {
0251     xml_comp_t x_sup   = su;
0252     const double rot_z = getAttrOrDefault(x_sup, _U(phi0), 0.);
0253     RotationZ rot(rot_z);
0254     Transform3D tr_rot(
0255         rot, Position(0, 0, 0)); // additional rotation of the module after position offset
0256 
0257     auto [vol, tr]           = build_shape(description, x_det, x_sup);
0258     [[maybe_unused]] auto pv = assembly.placeVolume(vol, tr_rot * tr);
0259     // Loop over support components, if any
0260     double cumulative_thickness = 0;
0261     for (xml_coll_t com{x_sup, _U(component)}; com; ++com) {
0262       xml_comp_t x_com = com;
0263       auto [cvol, ctr] = build_shape(description, x_det, x_sup, x_com, cumulative_thickness);
0264       vol.placeVolume(cvol, ctr);
0265       cumulative_thickness += x_com.thickness();
0266     }
0267   }
0268   // Loop over any subtraction volumes
0269   for (xml_coll_t sub{x_det, _Unicode(subtraction)}; sub; ++sub) {
0270     xml_comp_t x_sub = sub;
0271     auto [vol, tr]   = build_subtraction(description, x_det, x_sub);
0272     assembly.placeVolume(vol, tr);
0273   }
0274 
0275   // final placement
0276   Volume motherVol = description.pickMotherVolume(det);
0277   Position pos(0, 0, offset);
0278   PlacedVolume pv = motherVol.placeVolume(assembly, pos);
0279   pv.addPhysVolID("system", det.id());
0280   det.setPlacement(pv);
0281 
0282   return det;
0283 }
0284 
0285 // clang-format off
0286 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!