epic/src/BarrelTOFTracker_geo.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 - 2024, Whitney Armstrong, Chun Yuen Tsang
0003
0004 /** \addtogroup Trackers Trackers
0005  * \brief Type: **TOFBarrel**.
0006  * \author W. Armstrong
0007  * \modified by C.Y Tsang 3rd Aug, 2024
0008  *
0009  * \ingroup trackers
0010  *
0011  * @{
0012  */
0013 #include "DD4hep/DetFactoryHelper.h"
0014 #include "DD4hep/Printout.h"
0015 #include "DD4hep/Shapes.h"
0016 #include "DDRec/DetectorData.h"
0017 #include "DDRec/Surface.h"
0018 #include "XML/Layering.h"
0019 #include "XML/Utilities.h"
0020 #include <array>
0021 #include <iostream>
0022 #include "DD4hepDetectorHelper.h"
0023
0024 using namespace std;
0025 using namespace dd4hep;
0026 using namespace dd4hep::rec;
0027 using namespace dd4hep::detail;
0028
0029 /** Barrel Tracker with space frame.
0030  *
0031  * - Optional "support" tag within the detector element.
0032  *
0033  * The shapes are created using createShape which can be one of many basic geomtries.
0034  * See the examples Check_shape_*.xml in
0035  * [dd4hep's examples/ClientTests/compact](https://github.com/AIDASoft/DD4hep/tree/master/examples/ClientTests/compact)
0036  * directory.
0037  *
0038  *
0039  * - Optional "frame" tag within the module element.
0040  *
0041  * \ingroup trackers
0042  *
0043  * \code
0044  * \endcode
0045  *
0046  *
0047  * @author Whitney Armstrong
0048  */
0049 static Ref_t create_TOFBarrel(Detector& description, xml_h e, SensitiveDetector sens) {
0050   typedef vector<PlacedVolume> Placements;
0051   xml_det_t x_det = e;
0052   Material air    = description.air();
0053   int det_id      = x_det.id();
0054   string det_name = x_det.nameStr();
0055   DetElement sdet(det_name, det_id);
0056
0057   map<string, Volume> volumes;
0058   map<string, Placements> sensitives;
0059   map<string, std::vector<VolPlane>> volplane_surfaces;
0060   map<string, std::array<double, 2>> module_thicknesses;
0061
0062   PlacedVolume pv;
0063
0064   // Set detector type flag
0065   dd4hep::xml::setDetectorTypeFlag(x_det, sdet);
0066   auto& params = DD4hepDetectorHelper::ensureExtension<dd4hep::rec::VariantParameters>(sdet);
0067
0068   // Add the volume boundary material if configured
0069   for (xml_coll_t bmat(x_det, _Unicode(boundary_material)); bmat; ++bmat) {
0070     xml_comp_t x_boundary_material = bmat;
0071     DD4hepDetectorHelper::xmlToProtoSurfaceMaterial(x_boundary_material, params,
0072                                                     "boundary_material");
0073   }
0074
0075   // dd4hep::xml::Dimension dimensions(x_det.dimensions());
0076   // Tube topVolumeShape(dimensions.rmin(), dimensions.rmax(), dimensions.length() * 0.5);
0077   // Volume assembly(det_name,topVolumeShape,air);
0078   Assembly assembly(det_name);
0079
0080   sens.setType("tracker");
0081
0082   // Loop over the suports
0083   for (xml_coll_t su(x_det, _U(support)); su; ++su) {
0084     xml_comp_t x_support     = su;
0085     double support_thickness = getAttrOrDefault(x_support, _U(thickness), 2.0 * mm);
0086     double support_length    = getAttrOrDefault(x_support, _U(length), 2.0 * mm);
0087     double support_rmin      = getAttrOrDefault(x_support, _U(rmin), 2.0 * mm);
0088     double support_zstart    = getAttrOrDefault(x_support, _U(zstart), 2.0 * mm);
0089     std::string support_name =
0090         getAttrOrDefault<std::string>(x_support, _Unicode(name), "support_tube");
0091     std::string support_vis = getAttrOrDefault<std::string>(x_support, _Unicode(vis), "AnlRed");
0092     xml_dim_t pos(x_support.child(_U(position), false));
0093     xml_dim_t rot(x_support.child(_U(rotation), false));
0094     Solid support_solid;
0095     if (x_support.hasChild(_U(shape))) {
0096       xml_comp_t shape(x_support.child(_U(shape)));
0097       string shape_type = shape.typeStr();
0098       support_solid     = xml::createShape(description, shape_type, shape);
0099     } else {
0100       support_solid = Tube(support_rmin, support_rmin + support_thickness, support_length / 2);
0101     }
0102     Transform3D tr =
0103         Transform3D(Rotation3D(), Position(0, 0, (support_zstart + support_length / 2)));
0104     if (pos.ptr() && rot.ptr()) {
0105       Rotation3D rot3D(RotationZYX(rot.z(0), rot.y(0), rot.x(0)));
0106       Position pos3D(pos.x(0), pos.y(0), pos.z(0));
0107       tr = Transform3D(rot3D, pos3D);
0108     } else if (pos.ptr()) {
0109       tr = Transform3D(Rotation3D(), Position(pos.x(0), pos.y(0), pos.z(0)));
0110     } else if (rot.ptr()) {
0111       Rotation3D rot3D(RotationZYX(rot.z(0), rot.y(0), rot.x(0)));
0112       tr = Transform3D(rot3D, Position());
0113     }
0114     Material support_mat = description.material(x_support.materialStr());
0115     Volume support_vol(support_name, support_solid, support_mat);
0116     support_vol.setVisAttributes(description.visAttributes(support_vis));
0117     pv = assembly.placeVolume(support_vol, tr);
0118     // pv = assembly.placeVolume(support_vol, Position(0, 0, support_zstart + support_length / 2));
0119   }
0120
0121   // loop over the modules
0122   for (xml_coll_t mi(x_det, _U(module)); mi; ++mi) {
0123     xml_comp_t x_mod = mi;
0124     string m_nam     = x_mod.nameStr();
0125
0126     if (volumes.find(m_nam) != volumes.end()) {
0127       printout(ERROR, "TOFBarrel",
0128                string((string("Module with named ") + m_nam + string(" already exists."))).c_str());
0129       throw runtime_error("Logics error in building modules.");
0130     }
0131
0132     int ncomponents        = 0;
0133     int sensor_number      = 1;
0134     double total_thickness = 0;
0135
0136     // Compute module total thickness from components
0137     xml_coll_t ci(x_mod, _U(module_component));
0138     for (ci.reset(), total_thickness = 0.0; ci; ++ci) {
0139       if (!getAttrOrDefault<bool>(xml_comp_t(ci), _Unicode(keep_layer), false))
0140         total_thickness += xml_comp_t(ci).thickness();
0141     }
0142     // the module assembly volume
0143     Assembly m_vol(m_nam);
0144     volumes[m_nam] = m_vol;
0145     m_vol.setVisAttributes(description.visAttributes(x_mod.visStr()));
0146
0147     // Optional module frame.
0148     if (x_mod.hasChild(_U(frame))) {
0149       xml_comp_t m_frame = x_mod.child(_U(frame));
0150       // xmleles[m_nam]  = x_mod;
0151       double frame_thickness = m_frame.thickness();
0152       double frame_width     = m_frame.width();
0153       double frame_height    = getAttrOrDefault<double>(m_frame, _U(height), 5.0 * mm);
0154       double tanth           = frame_height / (frame_width / 2.0);
0155       double costh           = 1. / sqrt(1 + tanth * tanth);
0156       double frame_height2   = frame_height - frame_thickness - frame_thickness / costh;
0157       double frame_width2    = 2.0 * frame_height2 / tanth;
0158
0159       Trd1 moduleframe_part1(frame_width / 2, 0.001 * mm, m_frame.length() / 2, frame_height / 2);
0160       Trd1 moduleframe_part2(frame_width2 / 2, 0.001 * mm, m_frame.length() / 2 + 0.01 * mm,
0161                              frame_height2 / 2);
0162
0163       SubtractionSolid moduleframe(moduleframe_part1, moduleframe_part2,
0164                                    Position(0.0, frame_thickness, 0.0));
0165       Volume v_moduleframe(m_nam + "_vol", moduleframe,
0166                            description.material(m_frame.materialStr()));
0167       v_moduleframe.setVisAttributes(description, m_frame.visStr());
0168       m_vol.placeVolume(v_moduleframe,
0169                         Position(0.0, 0.0, frame_height / 2 + total_thickness / 2.0));
0170     }
0171
0172     double thickness_so_far = 0.0;
0173     double thickness_sum    = -total_thickness / 2.0;
0174     for (xml_coll_t mci(x_mod, _U(module_component)); mci; ++mci) {
0175       xml_comp_t x_comp = mci;
0176       xml_comp_t x_pos  = x_comp.position(false);
0177       xml_comp_t x_rot  = x_comp.rotation(false);
0178       auto make_box     = [&](double width, double length, double thickness, double pos_x = 0,
0179                           double pos_y = 0, double pos_z = 0, double rot_x = 0, double rot_y = 0,
0180                           double rot_z = 0, bool z_stacking = true) {
0181         // Utility variable for the relative z-offset based off the previous components
0182         const double zoff = thickness_sum + thickness / 2.0;
0183
0184         const string c_nam = _toString(ncomponents, "component%d");
0185         ++ncomponents;
0186         Box c_box(width / 2, length / 2, thickness / 2);
0187         Volume c_vol;
0188
0189         xml_coll_t ci_tube(x_comp, _Unicode(inner_tube));
0190         if (ci_tube) {
0191           double max_r = 0;
0192           for (; ci_tube; ++ci_tube) {
0193             // fill the hole with tube
0194             xml_comp_t ct = ci_tube;
0195             max_r         = std::max(max_r, ct.rmax());
0196             Tube c_tube(ct.rmin(), ct.rmax(), length / 2);
0197             Volume c_tubevol(c_nam + ct.nameStr(), c_tube, description.material(ct.materialStr()));
0198             if (ct.visStr() != "")
0199               c_tubevol.setVisAttributes(description, ct.visStr());
0200             m_vol.placeVolume(c_tubevol, Transform3D(RotationZYX(0, 0, -M_PI / 2),
0201                                                          Position(pos_x, pos_y, pos_z + zoff)));
0202           }
0203
0204           Tube c_fbox(0, max_r, length / 2 + 1);
0205           SubtractionSolid c_sbox(c_box, c_fbox,
0206                                       Transform3D(RotationZYX(0, 0, -M_PI / 2),
0207                                                   Position(0, 0, 0))); //pos_x, pos_y, pos_z + zoff)));
0208
0209           c_vol = Volume(c_nam, c_sbox, description.material(x_comp.materialStr()));
0210         } else
0211           c_vol = Volume(c_nam, c_box, description.material(x_comp.materialStr()));
0212
0213         Volume test;
0214         test = c_vol;
0215
0216         // center if off by half the box length if box length is cut in half
0217         Position c_pos(pos_x, pos_y, pos_z + zoff);
0218         RotationZYX c_rot(rot_z, rot_y, rot_x);
0219         pv = m_vol.placeVolume(c_vol, Transform3D(c_rot, c_pos));
0220
0221         c_vol.setRegion(description, x_comp.regionStr());
0222         c_vol.setLimitSet(description, x_comp.limitsStr());
0223         c_vol.setVisAttributes(description, x_comp.visStr());
0224         if (x_comp.isSensitive()) {
0225           pv.addPhysVolID("sensor", sensor_number++);
0226           c_vol.setSensitiveDetector(sens);
0227           sensitives[m_nam].push_back(pv);
0228           module_thicknesses[m_nam] = {thickness_so_far + thickness / 2.0,
0229                                        total_thickness - thickness_so_far - thickness / 2.0};
0230
0231           // -------- create a measurement plane for the tracking surface attched to the sensitive volume -----
0232           Vector3D u(-1., 0., 0.);
0233           Vector3D v(0., -1., 0.);
0234           Vector3D n(0., 0., 1.);
0235           //    Vector3D o( 0. , 0. , 0. ) ;
0236
0237           // compute the inner and outer thicknesses that need to be assigned to the tracking surface
0238           // depending on wether the support is above or below the sensor
0239           double inner_thickness = module_thicknesses[m_nam][0];
0240           double outer_thickness = module_thicknesses[m_nam][1];
0241
0242           SurfaceType type(SurfaceType::Sensitive);
0243
0244           // if( isStripDetector )
0245           //  type.setProperty( SurfaceType::Measurement1D , true ) ;
0246
0247           VolPlane surf(c_vol, type, inner_thickness, outer_thickness, u, v, n); //,o ) ;
0248           volplane_surfaces[m_nam].push_back(surf);
0249
0250           //--------------------------------------------
0251         }
0252         if (z_stacking) {
0253           thickness_sum += thickness;
0254           thickness_so_far += thickness;
0255           // apply relative offsets in z-position used to stack components side-by-side
0256           thickness_sum += pos_z;
0257           thickness_so_far += pos_z;
0258         }
0259       };
0260
0261       double pos_x = 0, pos_y = 0, pos_z = 0;
0262       double rot_x = 0, rot_y = 0, rot_z = 0;
0263       if (x_rot) {
0264         rot_x = x_rot.x(0);
0265         rot_y = x_rot.y(0);
0266         rot_z = x_rot.z(0);
0267       }
0268       if (x_pos) {
0269         pos_x = x_pos.x(0);
0270         pos_y = x_pos.y(0);
0271         pos_z = x_pos.z(0);
0272       }
0273       double width     = x_comp.width();
0274       double length    = x_comp.length();
0275       double thickness = x_comp.thickness();
0276       bool keep_layer  = getAttrOrDefault<bool>(x_comp, _Unicode(keep_layer), false);
0277
0278       if (x_comp.hasChild(_Unicode(GridSensors))) {
0279         auto x_comp_t = x_comp.child(_Unicode(GridSensors));
0280         // x-distance between centers of neighboring sensors
0281         double sensors_xdist = getAttrOrDefault<double>(x_comp_t, _Unicode(xdist), x_comp.width());
0282         // y-distance between centers of neighboring sensors
0283         double sensors_ydist = getAttrOrDefault<double>(x_comp_t, _Unicode(ydist), x_comp.length());
0284         // number of rows of sensors in a stave
0285         int nsensors_x = getAttrOrDefault<int>(x_comp_t, _Unicode(nx), 1);
0286         // number of column of sensors in a stave
0287         int nsensors_y = getAttrOrDefault<int>(x_comp_t, _Unicode(ny), 1);
0288         // x-location of the center of the leftmost sensor
0289         double start_x = getAttrOrDefault<double>(x_comp_t, _Unicode(start_x), 0);
0290         // y-location of the center of the uppermost sensor
0291         double start_y = getAttrOrDefault<double>(x_comp_t, _Unicode(start_y), 0);
0292         // z-locatino of the center of all sensors (All sensors appears at the same z-layer
0293         double start_z = getAttrOrDefault<double>(x_comp_t, _Unicode(start_z), 0);
0294         // central ring is located to the right of the ny_before_ring th sensor
0295         int ny_before_ring = getAttrOrDefault<int>(x_comp_t, _Unicode(ny_before_ring), 0);
0296         // Extra width caused by the ring
0297         // |<--sensors_ydist-->|<--sensors_ydist-->|<-----ring_extra_width------->|<--sensors_ydist-->|
0298         //   |<xcomp.width()>|   |<xcomp.width()>|                                   |<xcomp.width()>|
0299         //   ring_extra_width is the extra width between boundaries of the sensor boundaries (including dead space)
0300         double ring_extra_width = getAttrOrDefault<double>(x_comp_t, _Unicode(ring_extra_width), 0);
0301         auto half_length_str =
0302             getAttrOrDefault<std::string>(x_comp_t, _Unicode(half_length), "none");
0303
0304         double current_x = start_x;
0305         for (int nx = 0; nx < nsensors_x; ++nx) {
0306           double current_y = start_y;
0307           for (int ny = 0; ny < nsensors_y; ++ny) {
0308             double sensor_length     = length;
0309             double tmp_sensors_ydist = sensors_ydist;
0310             // when we draw half a sensor, the center has to be shifted by 0.25 times the length of a sensor
0311             // distance between centers to the next sensor also has to be reduced by 0.25 times the length of a sensor
0312             if ((half_length_str == "left" || half_length_str == "both") && ny == 0) {
0313               sensor_length = 0.5 * length;
0314               current_y += 0.25 * length;
0315               tmp_sensors_ydist -= 0.25 * length;
0316             }
0317             // same idea, but when you are drawing to the right, the right sensor center has to move in -y direction
0318             if ((half_length_str == "right" || half_length_str == "both") && ny == nsensors_y - 1) {
0319               sensor_length = 0.5 * length;
0320               current_y -= 0.25 * length;
0321             }
0322             make_box(width, sensor_length, thickness, current_x, current_y, start_z, rot_x, rot_y,
0323                      rot_z,
0324                      (((nx == nsensors_x - 1) && (ny == nsensors_y - 1))) &&
0325                          !keep_layer); // all sensors are located at the same z-layer
0326             // increment z-layers only at the end, after the last sensor is added
0327             current_y += tmp_sensors_ydist;
0328             if (ny + 1 == ny_before_ring)
0329               current_y += ring_extra_width;
0330           }
0331           current_x += sensors_xdist;
0332         }
0333       } else
0334         make_box(width, length, thickness, pos_x, pos_y, pos_z, rot_x, rot_y, rot_z, !keep_layer);
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 modules 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_TOFBarrel,       create_TOFBarrel)