epic/src/BarrelTrackerWithCurves_geo.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2026 Whitney Armstrong, Aditya Yalavatti, Sam Henry
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 "DDRec/DetectorData.h"
0016 #include "DDRec/Surface.h"
0017 #include "XML/Layering.h"
0018 #include "XML/Utilities.h"
0019 #include <array>
0020 #include "DD4hepDetectorHelper.h"
0021
0022 /*
0023 This version is an adaptation of BarrelTrackerWithFrame_geo.cpp to handle the curved silicon surfaces of the outer barrels
0024
0025 A curved component (identified by the word "Curved" in the name) is a segment of a cylinder. This is constructed by a series of flat segments
0026 The number of segments is set in the xml file, together with the radius, anglular range, length, thickness, and the inner and outer thicknesses
0027 used to set the sensitive surface.
0028
0029 This approach is taken as it is not possibe to use the DD4HEP CylindricalGridPhiZ readout for a cylindrical surface where the axis is displaced
0030 from the beam line. Therefore the curved surface is modelled as a series of flat segments.
0031
0032 Unlike in BarrelTrackerWithFrame, the position of components (related to the centre of the stave) must now be given in the xml files.
0033
0034 17 April 2026. Sam Henry
0035
0036 */
0037
0038 using namespace std;
0039 using namespace dd4hep;
0040 using namespace dd4hep::rec;
0041 using namespace dd4hep::detail;
0042
0043 /** Barrel Tracker with space frame.
0044  *
0045  * - Optional "support" tag within the detector element.
0046  *
0047  * The shapes are created using createShape which can be one of many basic geomtries.
0048  * See the examples Check_shape_*.xml in
0049  * [dd4hep's examples/ClientTests/compact](https://github.com/AIDASoft/DD4hep/tree/master/examples/ClientTests/compact)
0050  * directory.
0051  *
0052  *
0053  * - Optional "frame" tag within the module element.
0054  *
0055  * \ingroup trackers
0056  *
0057  * \code
0058  * \endcode
0059  *
0060  *
0061  * @author Whitney Armstrong
0062  */
0063 static Ref_t create_BarrelTrackerWithCurves(Detector& description, xml_h e,
0064                                             SensitiveDetector sens) {
0065   typedef vector<PlacedVolume> Placements;
0066   xml_det_t x_det = e;
0067   Material air    = description.air();
0068   int det_id      = x_det.id();
0069   string det_name = x_det.nameStr();
0070   DetElement sdet(det_name, det_id);
0071
0072   map<string, Volume> volumes;
0073   map<string, Placements> sensitives;
0074   map<string, std::vector<VolPlane>> volplane_surfaces;
0075   map<string, std::array<double, 2>> module_thicknesses;
0076
0077   PlacedVolume pv;
0078
0079   // Set detector type flag
0080   dd4hep::xml::setDetectorTypeFlag(x_det, sdet);
0081   auto& params = DD4hepDetectorHelper::ensureExtension<dd4hep::rec::VariantParameters>(sdet);
0082
0083   // Add the volume boundary material if configured
0084   for (xml_coll_t bmat(x_det, _Unicode(boundary_material)); bmat; ++bmat) {
0085     xml_comp_t x_boundary_material = bmat;
0086     DD4hepDetectorHelper::xmlToProtoSurfaceMaterial(x_boundary_material, params,
0087                                                     "boundary_material");
0088   }
0089
0090   // dd4hep::xml::Dimension dimensions(x_det.dimensions());
0091   // Tube topVolumeShape(dimensions.rmin(), dimensions.rmax(), dimensions.length() * 0.5);
0092   // Volume assembly(det_name,topVolumeShape,air);
0093   Assembly assembly(det_name);
0094
0095   sens.setType("tracker");
0096
0097   // Loop over the suports
0098   for (xml_coll_t su(x_det, _U(support)); su; ++su) {
0099     xml_comp_t x_support     = su;
0100     double support_thickness = getAttrOrDefault(x_support, _U(thickness), 2.0 * mm);
0101     double support_length    = getAttrOrDefault(x_support, _U(length), 2.0 * mm);
0102     double support_rmin      = getAttrOrDefault(x_support, _U(rmin), 2.0 * mm);
0103     double support_zstart    = getAttrOrDefault(x_support, _U(zstart), 2.0 * mm);
0104     std::string support_name =
0105         getAttrOrDefault<std::string>(x_support, _Unicode(name), "support_tube");
0106     std::string support_vis = getAttrOrDefault<std::string>(x_support, _Unicode(vis), "AnlRed");
0107     xml_dim_t pos(x_support.child(_U(position), false));
0108     xml_dim_t rot(x_support.child(_U(rotation), false));
0109     Solid support_solid;
0110     if (x_support.hasChild(_U(shape))) {
0111       xml_comp_t shape(x_support.child(_U(shape)));
0112       string shape_type = shape.typeStr();
0113       support_solid     = xml::createShape(description, shape_type, shape);
0114     } else {
0115       support_solid = Tube(support_rmin, support_rmin + support_thickness, support_length / 2);
0116     }
0117     Transform3D tr =
0118         Transform3D(Rotation3D(), Position(0, 0, (support_zstart + support_length / 2)));
0119     if (pos.ptr() && rot.ptr()) {
0120       Rotation3D rot3D(RotationZYX(rot.z(0), rot.y(0), rot.x(0)));
0121       Position pos3D(pos.x(0), pos.y(0), pos.z(0));
0122       tr = Transform3D(rot3D, pos3D);
0123     } else if (pos.ptr()) {
0124       tr = Transform3D(Rotation3D(), Position(pos.x(0), pos.y(0), pos.z(0)));
0125     } else if (rot.ptr()) {
0126       Rotation3D rot3D(RotationZYX(rot.z(0), rot.y(0), rot.x(0)));
0127       tr = Transform3D(rot3D, Position());
0128     }
0129     Material support_mat = description.material(x_support.materialStr());
0130     Volume support_vol(support_name, support_solid, support_mat);
0131     support_vol.setVisAttributes(description.visAttributes(support_vis));
0132     pv = assembly.placeVolume(support_vol, tr);
0133     // pv = assembly.placeVolume(support_vol, Position(0, 0, support_zstart + support_length / 2));
0134   }
0135
0136   // loop over the modules
0137   for (xml_coll_t mi(x_det, _U(module)); mi; ++mi) {
0138     xml_comp_t x_mod = mi;
0139     string m_nam     = x_mod.nameStr();
0140
0141     if (volumes.find(m_nam) != volumes.end()) {
0142       printout(ERROR, "BarrelTrackerWithFrame",
0143                string((string("Module with named ") + m_nam + string(" already exists."))).c_str());
0144       throw runtime_error("Logics error in building modules.");
0145     }
0146
0147     int ncomponents        = 0;
0148     int sensor_number      = 1;
0149     double total_thickness = 0;
0150
0151     // Compute module total thickness from components
0152     // add pos z to allow several components being placed at the same z without double-counting the total thickness.
0153     xml_coll_t ci(x_mod, _U(module_component));
0154     for (ci.reset(), total_thickness = 0.0; ci; ++ci) {
0155       xml_comp_t x_pos = xml_comp_t(ci).position(false);
0156       double mod_z_off = 0;
0157       if (x_pos)
0158         mod_z_off = x_pos.z(0);
0159       total_thickness += xml_comp_t(ci).thickness() + mod_z_off;
0160     }
0161     // the module assembly volume
0162     Assembly m_vol(m_nam);
0163     volumes[m_nam] = m_vol;
0164     m_vol.setVisAttributes(description.visAttributes(x_mod.visStr()));
0165
0166     // Optional module frame.
0167     if (x_mod.hasChild(_U(frame))) {
0168       xml_comp_t m_frame = x_mod.child(_U(frame));
0169       // xmleles[m_nam]  = x_mod;
0170       double frame_thickness = m_frame.thickness();
0171       double frame_width     = m_frame.width();
0172       double frame_height    = getAttrOrDefault<double>(m_frame, _U(height), 5.0 * mm);
0173       double tanth           = frame_height / (frame_width / 2.0);
0174       double costh           = 1. / sqrt(1 + tanth * tanth);
0175       double frame_height2   = frame_height - frame_thickness - frame_thickness / costh;
0176       double frame_width2    = 2.0 * frame_height2 / tanth;
0177
0178       Trd1 moduleframe_part1(frame_width / 2, 0.001 * mm, m_frame.length() / 2, frame_height / 2);
0179       Trd1 moduleframe_part2(frame_width2 / 2, 0.001 * mm, m_frame.length() / 2 + 0.01 * mm,
0180                              frame_height2 / 2);
0181
0182       SubtractionSolid moduleframe(moduleframe_part1, moduleframe_part2,
0183                                    Position(0.0, frame_thickness, 0.0));
0184       Volume v_moduleframe(m_nam + "_vol", moduleframe,
0185                            description.material(m_frame.materialStr()));
0186       v_moduleframe.setVisAttributes(description, m_frame.visStr());
0187       m_vol.placeVolume(v_moduleframe,
0188                         Position(0.0, 0.0, frame_height / 2 + total_thickness / 2.0));
0189     }
0190
0191     double thickness_so_far = 0.0;
0192
0193     for (xml_coll_t mci(x_mod, _U(module_component)); mci; ++mci, ++ncomponents) {
0194       xml_comp_t x_comp  = mci;
0195       xml_comp_t x_pos   = x_comp.position(false);
0196       xml_comp_t x_rot   = x_comp.rotation(false);
0197       const string c_nam = _toString(ncomponents, "component%d");
0198
0199       Volume c_vol;
0200       if (x_comp.nameStr().find("Curved") != std::string::npos) {
0201         double c_rmin = x_comp.radius(); // radius of curvature in cm
0202         double c_phi0 =
0203             x_comp.phi0(); // start and stop angle of segment in rad - zero is centre of stave
0204         double c_phi1 = x_comp.phi1();
0205         double c_Nseg = x_comp.nsegments(); //  number of flat segments used to approximate cylinder
0206         double dphi   = (c_phi1 - c_phi0) / (double)c_Nseg;
0207         vector<double> Xp, Zp;
0208         double phiP = c_phi0;
0209
0210         for (int i = 0; i < c_Nseg + 1; ++i) {
0211           Xp.push_back(c_rmin * sin(phiP));
0212           Zp.push_back(c_rmin * cos(phiP));
0213           phiP += dphi;
0214         }
0215         phiP = c_phi0 + dphi / 2;
0216         for (int i = 0; i < c_Nseg; ++i) {
0217           double segwidth = sqrt((Xp[i + 1] - Xp[i]) * (Xp[i + 1] - Xp[i]) +
0218                                  (Zp[i + 1] - Zp[i]) * (Zp[i + 1] - Zp[i]));
0219           double midX     = 0.5 * (Xp[i] + Xp[i + 1]);
0220           double midZ     = 0.5 * (Zp[i] + Zp[i + 1]);
0221           Box seg_box((segwidth / 2) * ((c_rmin - x_comp.thickness() / 2) / c_rmin),
0222                       x_comp.length() / 2, x_comp.thickness() / 2);
0223           string seg_nam = c_nam + "." + _toString(i);
0224           Volume seg_vol(seg_nam, seg_box, description.material(x_comp.materialStr()));
0225
0226           if (x_pos && x_rot) {
0227             Position c_pos(midX + x_pos.x(0), x_pos.y(0), midZ + x_pos.z(0));
0228             RotationZYX c_rot(x_rot.z(0), phiP + x_rot.y(0), x_rot.x(0));
0229             pv = m_vol.placeVolume(seg_vol, Transform3D(c_rot, c_pos));
0230           } else if (x_rot) {
0231             Position c_pos(midX, 0, midZ);
0232             pv = m_vol.placeVolume(
0233                 seg_vol,
0234                 Transform3D(RotationZYX(x_rot.z(0), x_rot.y(0) + phiP, x_rot.x(0)), c_pos));
0235           } else if (x_pos) {
0236             RotationZYX c_rot(0, phiP, 0);
0237             Position c_pos(midX + x_pos.x(0), x_pos.y(0), x_pos.z(0) + midZ);
0238             pv = m_vol.placeVolume(seg_vol, Transform3D(c_rot, c_pos));
0239           } else {
0240             RotationZYX c_rot(0, phiP, 0);
0241             Position c_pos(midX, 0, midZ);
0242             pv = m_vol.placeVolume(seg_vol, Transform3D(c_rot, c_pos));
0243           }
0244           phiP += dphi;
0245           seg_vol.setRegion(description, x_comp.regionStr());
0246           seg_vol.setLimitSet(description, x_comp.limitsStr());
0247           seg_vol.setVisAttributes(description, x_comp.visStr());
0248           if (x_comp.isSensitive()) {
0249             pv.addPhysVolID("sensor", sensor_number++);
0250             seg_vol.setSensitiveDetector(sens);
0251             sensitives[m_nam].push_back(pv);
0252             //        module_thicknesses[m_nam] = {x_comp.innerthickness(),x_comp.outerthickness()};
0253
0254             // -------- create a measurement plane for the tracking surface attched to the sensitive volume -----
0255             Vector3D u(-1., 0., 0.);
0256             Vector3D v(0., -1., 0.);
0257             Vector3D n(0., 0., 1.);
0258             //    Vector3D o( 0. , 0. , 0. ) ;
0259
0260             // compute the inner and outer thicknesses that need to be assigned to the tracking surface
0261             // depending on wether the support is above or below the sensor
0262             double inner_thickness =
0263                 getAttrOrDefault<double>(x_comp, _Unicode(innerthickness), 0.5 * mm);
0264             double outer_thickness =
0265                 getAttrOrDefault<double>(x_comp, _Unicode(outerthickness), 0.5 * mm);
0266             SurfaceType type(SurfaceType::Sensitive);
0267
0268             VolPlane surf(seg_vol, type, inner_thickness, outer_thickness, u, v, n); //,o ) ;
0269             volplane_surfaces[m_nam].push_back(surf);
0270
0271             //--------------------------------------------
0272           }
0273         }
0274         // thickness of curved layer is height of curve
0275
0276         thickness_so_far += x_comp.thickness();
0277
0278       } else { // not a curved component
0279         Box c_box(x_comp.width() / 2, x_comp.length() / 2, x_comp.thickness() / 2);
0280         c_vol = Volume(c_nam, c_box, description.material(x_comp.materialStr()));
0281
0282         if (x_pos && x_rot) {
0283           Position c_pos(x_pos.x(0), x_pos.y(0), x_pos.z(0));
0284           RotationZYX c_rot(x_rot.z(0), x_rot.y(0), x_rot.x(0));
0285           pv = m_vol.placeVolume(c_vol, Transform3D(c_rot, c_pos));
0286         } else if (x_rot) {
0287           Position c_pos(0, 0, 0);
0288           pv = m_vol.placeVolume(
0289               c_vol, Transform3D(RotationZYX(x_rot.z(0), x_rot.y(0), x_rot.x(0)), c_pos));
0290         } else if (x_pos) {
0291           pv = m_vol.placeVolume(c_vol, Position(x_pos.x(0), x_pos.y(0), x_pos.z(0)));
0292         } else {
0293           pv = m_vol.placeVolume(c_vol, Position(0, 0, 0));
0294         }
0295         c_vol.setRegion(description, x_comp.regionStr());
0296         c_vol.setLimitSet(description, x_comp.limitsStr());
0297         c_vol.setVisAttributes(description, x_comp.visStr());
0298         if (x_comp.isSensitive()) {
0299           pv.addPhysVolID("sensor", sensor_number++);
0300           c_vol.setSensitiveDetector(sens);
0301           sensitives[m_nam].push_back(pv);
0302           module_thicknesses[m_nam] = {thickness_so_far + x_comp.thickness() / 2.0,
0303                                        total_thickness - thickness_so_far -
0304                                            x_comp.thickness() / 2.0};
0305
0306           // -------- create a measurement plane for the tracking surface attched to the sensitive volume -----
0307           Vector3D u(-1., 0., 0.);
0308           Vector3D v(0., -1., 0.);
0309           Vector3D n(0., 0., 1.);
0310           //    Vector3D o( 0. , 0. , 0. ) ;
0311
0312           // compute the inner and outer thicknesses that need to be assigned to the tracking surface
0313           // depending on wether the support is above or below the sensor
0314           double inner_thickness = module_thicknesses[m_nam][0];
0315           double outer_thickness = module_thicknesses[m_nam][1];
0316
0317           SurfaceType type(SurfaceType::Sensitive);
0318
0319           // if( isStripDetector )
0320           //  type.setProperty( SurfaceType::Measurement1D , true ) ;
0321
0322           VolPlane surf(c_vol, type, inner_thickness, outer_thickness, u, v, n); //,o ) ;
0323           volplane_surfaces[m_nam].push_back(surf);
0324
0325           //--------------------------------------------
0326         }
0327
0328         thickness_so_far += x_comp.thickness();
0329       }
0330       // apply relative offsets in z-position used to stack components side-by-side
0331       if (x_pos) {
0332
0333         thickness_so_far += x_pos.z(0);
0334       }
0335     }
0336   }
0337
0338   // now build the layers
0339   for (xml_coll_t li(x_det, _U(layer)); li; ++li) {
0340     xml_comp_t x_layer  = li;
0341     xml_comp_t x_barrel = x_layer.child(_U(barrel_envelope));
0342     xml_comp_t x_layout = x_layer.child(_U(rphi_layout));
0343     xml_comp_t z_layout = x_layer.child(_U(z_layout)); // Get the <z_layout> element.
0344     int lay_id          = x_layer.id();
0345     string m_nam        = x_layer.moduleStr();
0346     string lay_nam      = det_name + _toString(x_layer.id(), "_layer%d");
0347     Tube lay_tub(x_barrel.inner_r(), x_barrel.outer_r(), x_barrel.z_length() / 2.0);
0348     Volume lay_vol(lay_nam, lay_tub, air); // Create the layer envelope volume.
0349     Position lay_pos(0, 0, getAttrOrDefault(x_barrel, _U(z0), 0.));
0350     lay_vol.setVisAttributes(description.visAttributes(x_layer.visStr()));
0351
0352     double phi0     = x_layout.phi0();     // Starting phi of first module.
0353     double phi_tilt = x_layout.phi_tilt(); // Phi tilt of a module.
0354     double rc       = x_layout.rc();       // Radius of the module center.
0355     int nphi        = x_layout.nphi();     // Number of modfules in phi.
0356     double rphi_dr  = x_layout.dr();       // The delta radius of every other module.
0357     double phi_incr = (M_PI * 2) / nphi;   // Phi increment for one module.
0358     double phic     = phi0;                // Phi of the module center.
0359     double z0       = z_layout.z0();       // Z position of first module in phi.
0360     double nz       = z_layout.nz();       // Number of modules to place in z.
0361     double z_dr     = z_layout.dr();       // Radial displacement parameter, of every other module.
0362
0363     Volume module_env = volumes[m_nam];
0364     DetElement lay_elt(sdet, lay_nam, lay_id);
0365     Placements& sensVols = sensitives[m_nam];
0366
0367     // the local coordinate systems of modules in dd4hep and acts differ
0368     // see http://acts.web.cern.ch/ACTS/latest/doc/group__DD4hepPlugins.html
0369     auto& layerParams =
0370         DD4hepDetectorHelper::ensureExtension<dd4hep::rec::VariantParameters>(lay_elt);
0371
0372     for (xml_coll_t lmat(x_layer, _Unicode(layer_material)); lmat; ++lmat) {
0373       xml_comp_t x_layer_material = lmat;
0374       DD4hepDetectorHelper::xmlToProtoSurfaceMaterial(x_layer_material, layerParams,
0375                                                       "layer_material");
0376     }
0377
0378     // Z increment for module placement along Z axis.
0379     // Adjust for z0 at center of module rather than
0380     // the end of cylindrical envelope.
0381     double z_incr = nz > 1 ? (2.0 * z0) / (nz - 1) : 0.0;
0382     // Starting z for module placement along Z axis.
0383     double module_z = -z0;
0384     int module      = 1;
0385
0386     // Loop over the number of modules in phi.
0387     for (int ii = 0; ii < nphi; ii++) {
0388       double dx = z_dr * std::cos(phic + phi_tilt); // Delta x of module position.
0389       double dy = z_dr * std::sin(phic + phi_tilt); // Delta y of module position.
0390       double x  = rc * std::cos(phic);              // Basic x module position.
0391       double y  = rc * std::sin(phic);              // Basic y module position.
0392
0393       // Loop over the number of modules in z.
0394       for (int j = 0; j < nz; j++) {
0395         string module_name = _toString(module, "module%d");
0396         DetElement mod_elt(lay_elt, module_name, module);
0397
0398         Transform3D tr(RotationZYX(0, ((M_PI / 2) - phic - phi_tilt), -M_PI / 2),
0399                        Position(x, y, module_z));
0400
0401         pv = lay_vol.placeVolume(module_env, tr);
0402         pv.addPhysVolID("module", module);
0403         mod_elt.setPlacement(pv);
0404         for (size_t ic = 0; ic < sensVols.size(); ++ic) {
0405           PlacedVolume sens_pv = sensVols[ic];
0406           DetElement comp_de(mod_elt, std::string("de_") + sens_pv.volume().name(), module);
0407           comp_de.setPlacement(sens_pv);
0408
0409           auto& comp_de_params =
0410               DD4hepDetectorHelper::ensureExtension<dd4hep::rec::VariantParameters>(comp_de);
0411           comp_de_params.set<string>("axis_definitions", "XYZ");
0412           // comp_de.setAttributes(description, sens_pv.volume(), x_layer.regionStr(), x_layer.limitsStr(),
0413           //                       xml_det_t(xmleles[m_nam]).visStr());
0414           //
0415
0416           volSurfaceList(comp_de)->push_back(volplane_surfaces[m_nam][ic]);
0417         }
0418
0419         /// Increase counters etc.
0420         module++;
0421         // Adjust the x and y coordinates of the module.
0422         x += dx;
0423         y += dy;
0424         // Flip sign of x and y adjustments.
0425         dx *= -1;
0426         dy *= -1;
0427         // Add z increment to get next z placement pos.
0428         module_z += z_incr;
0429       }
0430       phic += phi_incr; // Increment the phi placement of module.
0431       rc += rphi_dr;    // Increment the center radius according to dr parameter.
0432       rphi_dr *= -1;    // Flip sign of dr parameter.
0433       module_z = -z0;   // Reset the Z placement parameter for module.
0434     }
0435     // Create the PhysicalVolume for the layer.
0436     pv = assembly.placeVolume(lay_vol, lay_pos); // Place layer in mother
0437     pv.addPhysVolID("layer", lay_id);            // Set the layer ID.
0438     lay_elt.setAttributes(description, lay_vol, x_layer.regionStr(), x_layer.limitsStr(),
0439                           x_layer.visStr());
0440     lay_elt.setPlacement(pv);
0441   }
0442   sdet.setAttributes(description, assembly, x_det.regionStr(), x_det.limitsStr(), x_det.visStr());
0443   assembly.setVisAttributes(description.invisible());
0444   pv = description.pickMotherVolume(sdet).placeVolume(assembly);
0445   pv.addPhysVolID("system", det_id); // Set the subdetector system ID.
0446   sdet.setPlacement(pv);
0447   return sdet;
0448 }
0449
0450 //@}
0451 // clang-format off
0452 DECLARE_DETELEMENT(epic_CurvedTrackerBarrel,   create_BarrelTrackerWithCurves)