epic/src/BarrelPlanarMPGDTracker_geo.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Wouter Deconinck, Matt Posik
0003
0004 /** \addtogroup PID
0005  * \brief Type: **Barrel bar detector (rectangular geom.) with frames surrounding the perimeter.
0006  * \author S. Joosten
0007  * Modified by M. Posik
0008  *
0009  * \ingroup PID
0010  * \ingroup Tracking
0011  *
0012  *
0013  * \code
0014  * \endcode
0015  *
0016  * @{
0017  */
0018
0019 #include "DD4hep/DetFactoryHelper.h"
0020 #include "DD4hep/Printout.h"
0021 #include "DD4hep/Shapes.h"
0022 #include "DD4hepDetectorHelper.h"
0023 #include "DDRec/DetectorData.h"
0024 #include "DDRec/Surface.h"
0025 #include "XML/Layering.h"
0026 #include "XML/Utilities.h"
0027 #include <array>
0028
0029 using namespace std;
0030 using namespace dd4hep;
0031 using namespace dd4hep::rec;
0032
0033 /** Barrel Bar detector with optional frame
0034  *
0035  * - Optional "frame" tag within the module element (frame eats part of the module width and length
0036  *   and surrounds the rectangular perimeter)
0037  * - Detector is setup as a "tracker" so we can use the hits
0038  *
0039  */
0040 static Ref_t create_BarrelPlanarMPGDTracker_geo(Detector& description, xml_h e,
0041                                                 SensitiveDetector sens) {
0042   typedef vector<PlacedVolume> Placements;
0043   xml_det_t x_det = e;
0044   // Material                                air      = description.air();
0045   int det_id      = x_det.id();
0046   string det_name = x_det.nameStr();
0047   DetElement sdet(det_name, det_id);
0048   map<string, Volume> volumes;
0049   map<string, Placements> sensitives;
0050   map<string, std::vector<rec::VolPlane>> volplane_surfaces;
0051   PlacedVolume pv;
0052   dd4hep::xml::Dimension dimensions(x_det.dimensions());
0053   xml_dim_t mpgd_pos = x_det.position();
0054   Assembly assembly(det_name);
0055
0056   // Set detector type flag
0057   dd4hep::xml::setDetectorTypeFlag(x_det, sdet);
0058   auto& params = DD4hepDetectorHelper::ensureExtension<dd4hep::rec::VariantParameters>(sdet);
0059
0060   // Add the volume boundary material if configured
0061   for (xml_coll_t bmat(x_det, _Unicode(boundary_material)); bmat; ++bmat) {
0062     xml_comp_t x_boundary_material = bmat;
0063     DD4hepDetectorHelper::xmlToProtoSurfaceMaterial(x_boundary_material, params,
0064                                                     "boundary_material");
0065   }
0066
0067   map<string, std::array<double, 2>> module_thicknesses;
0068   sens.setType("tracker");
0069
0070   // loop over the modules
0071   for (xml_coll_t mi(x_det, _U(module)); mi; ++mi) {
0072     xml_comp_t x_mod = mi;
0073     string m_nam     = x_mod.nameStr();
0074
0075     if (volumes.find(m_nam) != volumes.end()) {
0076       printout(ERROR, "BarrelPlanarMPGDTracker_geo",
0077                string((string("Module with named ") + m_nam + string(" already exists."))).c_str());
0078       throw runtime_error("Logics error in building modules.");
0079     }
0080
0081     int ncomponents        = 0;
0082     int sensor_number      = 1;
0083     double total_thickness = 0;
0084
0085     // Compute module total thickness from components
0086     xml_coll_t ci(x_mod, _U(module_component));
0087     for (ci.reset(), total_thickness = 0.0; ci; ++ci) {
0088       total_thickness += xml_comp_t(ci).thickness();
0089     }
0090     // the module assembly volume
0091     Assembly m_vol(m_nam);
0092     volumes[m_nam] = m_vol;
0093     m_vol.setVisAttributes(description, x_mod.visStr());
0094
0095     // Optional module frame.
0096     // frame is 4  bars around the perimeter of the rectangular module. The frame will eat the
0097     // overlapping module area
0098     //
0099     //  ___
0100     // |___| <-- example module cross section (x-y plane), frame is flush with the
0101     //           bottom of the module and protrudes on the top if needed
0102
0103     // Get frame width, as it impacts the main module for being built. We
0104     // construct the actual frame structure later (once we know the module width)
0105     double frame_width = 0;
0106     if (x_mod.hasChild(_U(frame))) {
0107       xml_comp_t m_frame = x_mod.child(_U(frame));
0108       frame_width        = m_frame.width();
0109     }
0110
0111     double thickness_so_far     = 0.0;
0112     double thickness_sum        = -total_thickness / 2.0;
0113     double max_component_width  = 0;
0114     double max_component_length = 0;
0115     double gas_thickness        = 0.0;
0116     for (xml_coll_t mci(x_mod, _U(module_component)); mci; ++mci, ++ncomponents) {
0117       xml_comp_t x_comp = mci;
0118       string c_nam      = _toString(ncomponents, "component%d");
0119       string comp_name  = x_comp.nameStr();
0120
0121       double box_width  = x_comp.width();
0122       double box_length = x_comp.length();
0123       Box c_box;
0124       // Since MPGD frames are layed over the MPGD foils, the foil material is pressent under the frame as well.
0125       // The gas volumes are not present under the frames, so our frames must eat only the gas module areas
0126       //
0127       //   ------------------- MPGD foil
0128       //   --               -- Frame
0129       //   --  gas volume   -- Frame
0130       //   --               -- Frame
0131       //   ------------------- MPGD foil
0132
0133       // Look for gas modules to subtract frame thickness from
0134       // FIXME: these module names are hard coded for now. Should find
0135       // a way to set a arribut via the moduel tag to flag what components
0136       // need to have frame thickness subtracted.
0137       if ((comp_name == "DriftGap" || comp_name == "WindowGasGap")) {
0138         box_width            = x_comp.width() - 2.0 * frame_width;
0139         box_length           = x_comp.length() - 2.0 * frame_width;
0140         max_component_width  = box_width;
0141         max_component_length = box_length;
0142         gas_thickness += x_comp.thickness();
0143         c_box = {box_width / 2, box_length / 2, x_comp.thickness() / 2};
0144         printout(DEBUG, "BarrelPlanarMPGDTracker_geo", "gas: %s", comp_name.c_str());
0145         printout(DEBUG, "BarrelPlanarMPGDTracker_geo", "box_width: %f", box_width);
0146         printout(DEBUG, "BarrelPlanarMPGDTracker_geo", "box_length: %f", box_length);
0147         printout(DEBUG, "BarrelPlanarMPGDTracker_geo", "box_thickness: %f", x_comp.thickness());
0148       } else {
0149         c_box = {x_comp.width() / 2, x_comp.length() / 2, x_comp.thickness() / 2};
0150         printout(DEBUG, "BarrelPlanarMPGDTracker_geo", "Not gas: %s", comp_name.c_str());
0151         printout(DEBUG, "BarrelPlanarMPGDTracker_geo", "box_comp_width: %f", x_comp.width());
0152         printout(DEBUG, "BarrelPlanarMPGDTracker_geo", "box_comp_length: %f", x_comp.length());
0153         printout(DEBUG, "BarrelPlanarMPGDTracker_geo", "box_comp_thickness: %f",
0154                  x_comp.thickness());
0155       }
0156       Volume c_vol{c_nam, c_box, description.material(x_comp.materialStr())};
0157
0158       c_vol.setRegion(description, x_comp.regionStr());
0159       c_vol.setLimitSet(description, x_comp.limitsStr());
0160       c_vol.setVisAttributes(description, x_comp.visStr());
0161
0162       pv = m_vol.placeVolume(c_vol, Position(0, 0, thickness_sum + x_comp.thickness() / 2.0));
0163
0164       if (x_comp.isSensitive()) {
0165         pv.addPhysVolID("sensor", sensor_number++);
0166         c_vol.setSensitiveDetector(sens);
0167         sensitives[m_nam].push_back(pv);
0168         module_thicknesses[m_nam] = {thickness_so_far + x_comp.thickness() / 2.0,
0169                                      total_thickness - thickness_so_far - x_comp.thickness() / 2.0};
0170         // -------- create a measurement plane for the tracking surface attched to the sensitive volume -----
0171         Vector3D u(-1., 0., 0.);
0172         Vector3D v(0., -1., 0.);
0173         Vector3D n(0., 0., 1.);
0174
0175         // compute the inner and outer thicknesses that need to be assigned to the tracking surface
0176         // depending on wether the support is above or below the sensor
0177         double inner_thickness = module_thicknesses[m_nam][0];
0178         double outer_thickness = module_thicknesses[m_nam][1];
0179
0180         SurfaceType type(rec::SurfaceType::Sensitive);
0181
0182         VolPlane surf(c_vol, type, inner_thickness, outer_thickness, u, v, n); //,o ) ;
0183         volplane_surfaces[m_nam].push_back(surf);
0184       }
0185       thickness_sum += x_comp.thickness();
0186       thickness_so_far += x_comp.thickness();
0187     }
0188     // Now add-on the frame
0189     if (x_mod.hasChild(_U(frame))) {
0190       xml_comp_t m_frame     = x_mod.child(_U(frame));
0191       double frame_thickness = getAttrOrDefault<double>(m_frame, _U(thickness), total_thickness);
0192
0193       Box lframe_box{m_frame.width() / 2.0, (max_component_length + 2.0 * m_frame.width()) / 2.0,
0194                      frame_thickness / 2.0};
0195       Box rframe_box{m_frame.width() / 2.0, (max_component_length + 2.0 * m_frame.width()) / 2.0,
0196                      frame_thickness / 2.0};
0197       Box tframe_box{max_component_width / 2.0, m_frame.width() / 2.0, frame_thickness / 2.0};
0198       Box bframe_box{max_component_width / 2.0, m_frame.width() / 2.0, frame_thickness / 2.0};
0199
0200       // Keep track of frame with so we can adjust the module bars appropriately
0201
0202       Volume lframe_vol{"left_frame", lframe_box, description.material(m_frame.materialStr())};
0203       Volume rframe_vol{"right_frame", rframe_box, description.material(m_frame.materialStr())};
0204       Volume tframe_vol{"top_frame", tframe_box, description.material(m_frame.materialStr())};
0205       Volume bframe_vol{"bottom_frame", bframe_box, description.material(m_frame.materialStr())};
0206
0207       lframe_vol.setVisAttributes(description, m_frame.visStr());
0208       rframe_vol.setVisAttributes(description, m_frame.visStr());
0209       tframe_vol.setVisAttributes(description, m_frame.visStr());
0210       bframe_vol.setVisAttributes(description, m_frame.visStr());
0211
0212       printout(DEBUG, "BarrelPlanarMPGDTracker_geo", "frame_thickness: %f", frame_thickness);
0213       printout(DEBUG, "BarrelPlanarMPGDTracker_geo", "total_thickness: %f", total_thickness);
0214       printout(DEBUG, "BarrelPlanarMPGDTracker_geo", "frame_thickness + total_thickness: %f",
0215                frame_thickness + total_thickness);
0216
0217       m_vol.placeVolume(lframe_vol, Position(frame_width / 2.0 + max_component_width / 2, 0.0,
0218                                              frame_thickness / 2.0 - total_thickness / 2.0 -
0219                                                  gas_thickness / 2.0));
0220       m_vol.placeVolume(rframe_vol, Position(-frame_width / 2.0 - max_component_width / 2.0, 0.0,
0221                                              frame_thickness / 2.0 - total_thickness / 2.0 -
0222                                                  gas_thickness / 2.0));
0223       m_vol.placeVolume(tframe_vol, Position(0.0, frame_width / 2.0 + max_component_length / 2,
0224                                              frame_thickness / 2.0 - total_thickness / 2.0 -
0225                                                  gas_thickness / 2.0));
0226       m_vol.placeVolume(bframe_vol, Position(0.0, -frame_width / 2.0 - max_component_length / 2.0,
0227                                              frame_thickness / 2.0 - total_thickness / 2.0 -
0228                                                  gas_thickness / 2.0));
0229     }
0230   }
0231
0232   // build the layers the modules will be arranged around
0233   for (xml_coll_t li(x_det, _U(layer)); li; ++li) {
0234     xml_comp_t x_layer            = li;
0235     xml_comp_t x_layout           = x_layer.child(_U(rphi_layout));
0236     xml_comp_t z_layout           = x_layer.child(_U(z_layout));
0237     int lay_id                    = x_layer.id();
0238     string m_nam                  = x_layer.moduleStr();
0239     string lay_nam                = det_name + _toString(x_layer.id(), "_layer%d");
0240     xml_comp_t envelope_tolerance = x_layer.child(_Unicode(envelope_tolerance), false);
0241     double envelope_r_min         = 0;
0242     double envelope_r_max         = 0;
0243     double envelope_z_min         = 0;
0244     double envelope_z_max         = 0;
0245     if (envelope_tolerance) {
0246       envelope_r_min = getAttrOrDefault(envelope_tolerance, _Unicode(r_min), 0);
0247       envelope_r_max = getAttrOrDefault(envelope_tolerance, _Unicode(r_max), 0);
0248       envelope_z_min = getAttrOrDefault(envelope_tolerance, _Unicode(z_min), 0);
0249       envelope_z_max = getAttrOrDefault(envelope_tolerance, _Unicode(z_max), 0);
0250     }
0251
0252     double phi0     = x_layout.phi0();     // starting phi of first module
0253     double phi_tilt = x_layout.phi_tilt(); // Phi tilit of module
0254     double rc       = x_layout.rc();       // Radius of the module
0255     int nphi        = x_layout.nphi();     // Number of modules in phi
0256     double rphi_dr  = x_layout.dr();       // The delta radius of every other module
0257     double phi_incr = (2 * M_PI) / nphi;   // Phi increment for one module
0258     double phic     = phi0;                // Phi of the module
0259     double nz       = z_layout.nz();       // Number of modules placed in z
0260     double z_dr     = z_layout.dr();       // Radial offest of modules in z
0261     double z0       = z_layout.z0();       // Sets how much overlap in z the nz modules have
0262
0263     Assembly layer_assembly(lay_nam);
0264     Volume module_env = volumes[m_nam];
0265     DetElement lay_elt(sdet, lay_nam, lay_id);
0266     Placements& sensVols = sensitives[m_nam];
0267     auto& layerParams =
0268         DD4hepDetectorHelper::ensureExtension<dd4hep::rec::VariantParameters>(lay_elt);
0269
0270     pv = assembly.placeVolume(layer_assembly);
0271     pv.addPhysVolID("layer", lay_id);
0272     lay_elt.setPlacement(pv);
0273
0274     int module = 1;
0275     // loop over the modules in phi
0276     for (int ii = 0; ii < nphi; ii++) {
0277       double xc = rc * std::cos(phic);              // Basic x position of module
0278       double yc = rc * std::sin(phic);              // Basic y position of module
0279       double dx = z_dr * std::cos(phic + phi_tilt); // Deta x of module position
0280       double dy = z_dr * std::sin(phic + phi_tilt); // Deta y of module position
0281       // loop over the modules in z
0282       for (int j = 0; j < nz; j++) {
0283         string module_name = _toString(module, "module%d");
0284         DetElement mod_elt(lay_elt, module_name, module);
0285         double mod_z       = 0.5 * dimensions.length();
0286         double z_placement = mod_z - j * nz * mod_z; // z location for module placement
0287         double z_offset =
0288             z_placement > 0
0289                 ? -z0 / 2.0
0290                 : z0 / 2.0; // determine the amount of overlap in z the z nz modules have
0291
0292         Transform3D tr(RotationZYX(0.0, ((M_PI / 2) - phic - phi_tilt), -M_PI / 2),
0293                        Position(xc, yc, mpgd_pos.z() + z_placement + z_offset)); // in x-y plane,
0294         pv = layer_assembly.placeVolume(module_env, tr);
0295         pv.addPhysVolID("module", module);
0296         mod_elt.setPlacement(pv);
0297         for (size_t ic = 0; ic < sensVols.size(); ++ic) {
0298           PlacedVolume sens_pv = sensVols[ic];
0299           DetElement comp_de(mod_elt, std::string("de_") + sens_pv.volume().name(), module);
0300           comp_de.setPlacement(sens_pv);
0301         }
0302         // increas module counter
0303         module++;
0304         // adjust x and y coordinates
0305         xc += dx;
0306         yc += dy;
0307       }
0308       // increment counters
0309       phic += phi_incr;
0310       rc += rphi_dr;
0311     }
0312     layer_assembly->GetShape()->ComputeBBox();
0313     layerParams.set<double>("envelope_r_min", envelope_r_min);
0314     layerParams.set<double>("envelope_r_max", envelope_r_max);
0315     layerParams.set<double>("envelope_z_min", envelope_z_min);
0316     layerParams.set<double>("envelope_z_max", envelope_z_max);
0317
0318     for (xml_coll_t lmat(x_layer, _Unicode(layer_material)); lmat; ++lmat) {
0319       xml_comp_t x_layer_material = lmat;
0320       DD4hepDetectorHelper::xmlToProtoSurfaceMaterial(x_layer_material, layerParams,
0321                                                       "layer_material");
0322     }
0323   }
0324   sdet.setAttributes(description, assembly, x_det.regionStr(), x_det.limitsStr(), x_det.visStr());
0325   assembly.setVisAttributes(description.invisible());
0326   pv = description.pickMotherVolume(sdet).placeVolume(assembly);
0327   pv.addPhysVolID("system", det_id); // Set the subdetector system ID
0328   sdet.setPlacement(pv);
0329   return sdet;
0330 }
0331
0332 //@}
0333 // clang-format off
0334 DECLARE_DETELEMENT(epic_OuterMPGDBarrel, create_BarrelPlanarMPGDTracker_geo)