epic/src/SupportServiceMaterial_geo.cpp

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   }
0053   // A disk is a cylinder, constructed differently
0054   else if (type == "Disk") {
0055     const double thickness = getAttrOrDefault(x_child, _U(thickness), x_support.thickness());
0056     const double rmin      = getAttrOrDefault(x_child, _U(rmin), x_support.rmin());
0057     const double rmax      = getAttrOrDefault(x_child, _U(rmax), x_support.rmax());
0058     const double phimin    = getAttrOrDefault(
0059         x_child, _Unicode(phimin), getAttrOrDefault(x_support, _Unicode(phimin), 0.0 * deg));
0060     const double phimax = getAttrOrDefault(
0061         x_child, _Unicode(phimax), getAttrOrDefault(x_support, _Unicode(phimax), 360.0 * deg));
0062     pos3D = pos3D + Position(0, 0, -x_support.thickness() / 2 + thickness / 2 + offset);
0063     solid = Tube(rmin, rmax, thickness / 2, phimin, phimax);
0064   } else if (type == "Cone") {
0065     const double base_rmin1 = getAttrOrDefault(x_child, _U(rmin1), x_support.rmin1());
0066     const double base_rmin2 = getAttrOrDefault(x_child, _U(rmin2), x_support.rmin2());
0067     const double length     = getAttrOrDefault(x_child, _U(length), x_support.length());
0068     // Account for the fact that the distance between base_rmin1 and rmax2 is the projection
0069     // of the thickness on the transverse direction
0070     const double thickness = getAttrOrDefault(x_child, _U(thickness), x_support.thickness());
0071     const double transverse_thickness =
0072         thickness / cos(atan2(fabs(base_rmin2 - base_rmin1), length));
0073     // also account that the same is true for the offset
0074     const double transverse_offset = offset / cos(atan2(fabs(base_rmin2 - base_rmin1), length));
0075     const double rmin1             = base_rmin1 + transverse_offset;
0076     const double rmin2             = base_rmin2 + transverse_offset;
0077     const double rmax1             = rmin1 + transverse_thickness;
0078     const double rmax2             = rmin2 + transverse_thickness;
0079     // Allow for optional hard rmin/rmax cutoffs
0080     const double rmin = getAttrOrDefault(
0081         x_child, _U(rmin), getAttrOrDefault(x_support, _Unicode(rmin), min(rmin1, rmin2)));
0082     const double rmax = getAttrOrDefault(
0083         x_child, _U(rmax), getAttrOrDefault(x_support, _Unicode(rmax), max(rmax1, rmax2)));
0084     if (rmin > min(rmax1, rmax2)) {
0085       printout(ERROR, x_det.nameStr(),
0086                "%s: rmin (%f mm) must be smaller than the smallest rmax (%f %f mm)",
0087                x_support.nameStr().c_str(), rmin / mm, rmax1 / mm, rmax2 / mm);
0088       std::exit(1);
0089     }
0090     if (rmax < max(base_rmin1, base_rmin2)) {
0091       printout(ERROR, x_det.nameStr(),
0092                "%s: rmax (%f mm) must be larger than the largest rmin (%f %f mm)",
0093                x_support.nameStr().c_str(), rmax / mm, base_rmin1 / mm, base_rmin2 / mm);
0094       std::exit(1);
0095     }
0096     const double zmin  = -length / 2 + length * (rmin - rmin1) / (rmin2 - rmin1);
0097     const double zmax  = -length / 2 + length * (rmax - rmax1) / (rmax2 - rmax1);
0098     const auto rmin_at = [&](const double z) {
0099       return rmin1 + (z + length / 2) * (rmin2 - rmin1) / length;
0100     };
0101     const auto rmax_at = [&](const double z) {
0102       return rmax1 + (z + length / 2) * (rmax2 - rmax1) / length;
0103     };
0104     // Allow for optional phimin/phimax
0105     const double phimin = getAttrOrDefault<double>(
0106         x_child, _Unicode(phimin), getAttrOrDefault(x_support, _Unicode(phimin), 0.0 * deg));
0107     const double phimax = getAttrOrDefault<double>(
0108         x_child, _Unicode(phimax), getAttrOrDefault(x_support, _Unicode(phimax), 360.0 * deg));
0109     const double deltaphi = phimax - phimin;
0110     const double epsilon{TGeoShape::Tolerance()};
0111     if (fabs(zmin) >= length / 2 - epsilon && fabs(zmax) >= length / 2 - epsilon) {
0112       if (fabs(phimax - phimin - 360 * deg) < epsilon) {
0113         solid = Cone(length / 2, rmin1, rmax1, rmin2, rmax2);
0114       } else {
0115         solid = ConeSegment(length / 2, rmin1, rmax1, rmin2, rmax2, phimin, phimax);
0116       }
0117     } else {
0118       std::vector<double> v_rmin{max(rmin1, rmin), max(rmin2, rmin)},
0119           v_rmax{min(rmax1, rmax), min(rmax2, rmax)}, v_z{-length / 2, +length / 2};
0120       for (const auto& z :
0121            (zmin < zmax ? std::vector<double>{zmin, zmax} : std::vector<double>{zmax, zmin})) {
0122         if (-length / 2 + epsilon < z && z < -epsilon + length / 2) {
0123           v_rmin.insert(std::prev(v_rmin.end()), std::max(rmin, rmin_at(z)));
0124           v_rmax.insert(std::prev(v_rmax.end()), std::min(rmax, rmax_at(z)));
0125           v_z.insert(std::prev(v_z.end()), z);
0126         }
0127       }
0128       solid = Polycone(phimin, deltaphi, v_rmin, v_rmax, v_z);
0129     }
0130   } else {
0131     printout(ERROR, x_det.nameStr(), "Unknown support type: %s", type.c_str());
0132     std::exit(1);
0133   }
0134   // Materials
0135   Material mat = descr.material(getAttrOrDefault<std::string>(x_child, _U(material), "Air"));
0136   // Create our volume
0137   Volume vol{getAttrOrDefault<std::string>(x_child, _U(name), "support_vol"), solid, mat};
0138
0139   // Create full transformation
0140   Transform3D tr(rot3D, pos3D);
0141
0142   // visualization?
0143   if (x_child.hasAttr(_U(vis))) {
0144     vol.setVisAttributes(descr.visAttributes(x_child.visStr()));
0145   }
0146   return {vol, tr};
0147 }
0148 std::pair<Volume, Transform3D> build_shape(const Detector& descr, const xml_det_t& x_det,
0149                                            const xml_comp_t& x_support, const double offset = 0) {
0150   return build_shape(descr, x_det, x_support, x_support, offset);
0151 }
0152 } // namespace
0153
0154 /** Generic tracker support implementation, can consist of arbitrary shapes
0155  *
0156  * @author Sylvester Joosten
0157  */
0158 static Ref_t create_SupportServiceMaterial(Detector& description, xml_h e,
0159                                            [[maybe_unused]] SensitiveDetector sens) {
0160   const xml_det_t x_det = e;
0161   const int det_id      = x_det.id();
0162   const string det_name = x_det.nameStr();
0163
0164   // global z-offset for the entire support assembly
0165   const double offset = getAttrOrDefault(x_det, _U(offset), 0.);
0166
0167   DetElement det(det_name, det_id);
0168   Assembly assembly(det_name + "_assembly");
0169
0170   // Loop over the supports
0171   for (xml_coll_t su{x_det, _U(support)}; su; ++su) {
0172     xml_comp_t x_sup         = su;
0173     auto [vol, tr]           = build_shape(description, x_det, x_sup);
0174     [[maybe_unused]] auto pv = assembly.placeVolume(vol, tr);
0175     // Loop over support components, if any
0176     double cumulative_thickness = 0;
0177     for (xml_coll_t com{x_sup, _U(component)}; com; ++com) {
0178       xml_comp_t x_com = com;
0179       auto [cvol, ctr] = build_shape(description, x_det, x_sup, x_com, cumulative_thickness);
0180       vol.placeVolume(cvol, ctr);
0181       cumulative_thickness += x_com.thickness();
0182     }
0183   }
0184
0185   // final placement
0186   Volume motherVol = description.pickMotherVolume(det);
0187   Position pos(0, 0, offset);
0188   PlacedVolume pv = motherVol.placeVolume(assembly, pos);
0189   pv.addPhysVolID("system", det.id());
0190   det.setPlacement(pv);
0191
0192   return det;
0193 }
0194
0195 // clang-format off
0196 DECLARE_DETELEMENT(epic_SupportServiceMaterial, create_SupportServiceMaterial)