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  * - Two distinct versions of XML are covered:
0040  *  I) Multiple Sensitive Volume (fileName typically tagged "_2DStrip"),
0041  * II) Single Sensitive Volume (which single volume is then named "DriftGap").
0042  *   The XMLs have a distinctive setting of the <constant>
0043  *        "MPGDOuterBarrel_2DStrip",
0044  *   which "eicrecon" bases itself on to determine which version is in force.
0045  *   This can still be overridden by option "MPGD:SiFactoryPattern".
0046  *
0047  */
0048 static Ref_t create_BarrelPlanarMPGDTracker_geo(Detector& description, xml_h e,
0049                                                 SensitiveDetector sens) {
0050   xml_det_t x_det = e;
0051   // Material                                air      = description.air();
0052   int det_id      = x_det.id();
0053   string det_name = x_det.nameStr();
0054   DetElement sdet(det_name, det_id);
0055
0056   Volume* volume;
0057   // Sensitive volumes and associated surfaces
0058   // - There can be either one or five.
0059   int sensitiveVolumeSet = 5; // 1: single volume, 5: 5 volumes, -1: error
0060   vector<PlacedVolume> sensitives;
0061   vector<VolPlane> volplane_surfaces;
0062
0063   PlacedVolume pv;
0064
0065   //#define DEBUG_BarrelPlanarMPGDTracker
0066 #ifdef DEBUG_BarrelPlanarMPGDTracker
0067   // TEMPORARILY INCREASE VERBOSITY level for debugging purposes
0068   PrintLevel priorPrintLevel = printLevel();
0069   setPrintLevel(DEBUG);
0070 #endif
0071
0072   dd4hep::xml::Dimension dimensions(x_det.dimensions());
0073   xml_dim_t mpgd_pos = x_det.position();
0074   Assembly assembly(det_name);
0075
0076   double pcb_feb_ext = 0.0; //extension of PCB board to hold FEBs.
0077
0078   // Set detector type flag
0079   dd4hep::xml::setDetectorTypeFlag(x_det, sdet);
0080   auto& params = DD4hepDetectorHelper::ensureExtension<dd4hep::rec::VariantParameters>(sdet);
0081
0082   // Add the volume boundary material if configured
0083   for (xml_coll_t bmat(x_det, _Unicode(boundary_material)); bmat; ++bmat) {
0084     xml_comp_t x_boundary_material = bmat;
0085     DD4hepDetectorHelper::xmlToProtoSurfaceMaterial(x_boundary_material, params,
0086                                                     "boundary_material");
0087   }
0088
0089   sens.setType("tracker");
0090
0091   // ********** MODULE
0092   // ***** ONE AND ONLY ONE MODULE
0093   xml_coll_t modules(x_det, _U(module));
0094   if (modules.size() != 1) {
0095     // Present detector constructor can only handle ONE <module> tag
0096     printout(ERROR, "BarrelPlanarMPGDTracker_geo", "Number of modules = %u. Must be = 1",
0097              modules.size());
0098     throw runtime_error("Logics error in building modules.");
0099   }
0100   xml_comp_t x_mod = modules;
0101   string m_nam     = x_mod.nameStr();
0102
0103   int ncomponents        = 0;
0104   int sensor_number      = 0;
0105   double total_thickness = 0;
0106
0107   // Compute module total thickness from components
0108   xml_coll_t ci(x_mod, _U(module_component));
0109   for (ci.reset(), total_thickness = 0.0; ci; ++ci) {
0110     total_thickness += xml_comp_t(ci).thickness();
0111   }
0112   // the module assembly volume
0113   Assembly m_vol(m_nam);
0114   volume = &m_vol;
0115   m_vol.setVisAttributes(description, x_mod.visStr());
0116
0117   // Optional module frame.
0118   // frame is 4  bars around the perimeter of the rectangular module. The frame will eat the
0119   // overlapping module area
0120   //
0121   //  ___
0122   // |___| <-- example module cross section (x-y plane), frame is flush with the
0123   //           bottom of the module and protrudes on the top if needed
0124
0125   // Get frame width, as it impacts the main module for being built. We
0126   // construct the actual frame structure later (once we know the module width)
0127   double frame_width = 0;
0128   if (x_mod.hasChild(_U(frame))) {
0129     xml_comp_t m_frame = x_mod.child(_U(frame));
0130     frame_width        = m_frame.width();
0131   }
0132
0133   double thickness_so_far     = 0.0;
0134   double thickness_sum        = -total_thickness / 2.0;
0135   double max_component_width  = 0;
0136   double max_component_length = 0;
0137   double gas_thickness        = 0.0;
0138   // Pattern of Multiple Sensitive Volumes
0139   // - In order to have one sensitive component per strip coordinate (and
0140   //  accessorily, some extras), the "DriftGap" is subdivided into subVolumes.
0141   int nSensitives = 0;
0142   for (xml_coll_t mci(x_mod, _U(module_component)); mci; ++mci, ++ncomponents) {
0143     xml_comp_t x_comp     = mci;
0144     string c_nam          = x_comp.nameStr();
0145     double comp_thickness = x_comp.thickness();
0146     double box_width      = x_comp.width();
0147     double box_length     = x_comp.length();
0148     Box c_box;
0149     // Since MPGD frames are layed over the MPGD foils, the foil material is pressent under the frame as well.
0150     // The gas volumes are not present under the frames, so our frames must eat only the gas module areas
0151     //
0152     //   ------------------- MPGD foil
0153     //   --               -- Frame
0154     //   --  gas volume   -- Frame
0155     //   --               -- Frame
0156     //   ------------------- MPGD foil
0157
0158     // Look for gas modules to subtract frame thickness from
0159     // FIXME: these module names are hard coded for now. Should find
0160     // a way to set an attribute via the module tag to flag what components
0161     // need to have frame thickness subtracted.
0162     bool isDriftGap = c_nam == "DriftGap" ||
0163                       /* */ c_nam.find("ThinGap") != std::string::npos ||
0164                       /* */ c_nam.find("Radiator") != std::string::npos;
0165     if (isDriftGap || c_nam == "WindowGasGap") {
0166       box_width            = x_comp.width() - 2.0 * frame_width;
0167       box_length           = x_comp.length() - 2.0 * frame_width;
0168       max_component_width  = box_width;
0169       max_component_length = box_length;
0170       gas_thickness += comp_thickness;
0171       c_box = {box_width / 2, box_length / 2, comp_thickness / 2};
0172       printout(DEBUG, "BarrelPlanarMPGDTracker_geo", "gas: %s", c_nam.c_str());
0173       printout(DEBUG, "BarrelPlanarMPGDTracker_geo", "box_width: %f", box_width);
0174       printout(DEBUG, "BarrelPlanarMPGDTracker_geo", "box_length: %f", box_length);
0175       printout(DEBUG, "BarrelPlanarMPGDTracker_geo", "box_thickness: %f", comp_thickness);
0176     } else {
0177       c_box = {x_comp.width() / 2, x_comp.length() / 2, comp_thickness / 2};
0178       printout(DEBUG, "BarrelPlanarMPGDTracker_geo", "Not gas: %s", c_nam.c_str());
0179       printout(DEBUG, "BarrelPlanarMPGDTracker_geo", "box_comp_width: %f", x_comp.width());
0180       printout(DEBUG, "BarrelPlanarMPGDTracker_geo", "box_comp_length: %f", x_comp.length());
0181       printout(DEBUG, "BarrelPlanarMPGDTracker_geo", "box_comp_thickness: %f", comp_thickness);
0182     }
0183     Volume c_vol{c_nam, c_box, description.material(x_comp.materialStr())};
0184
0185     c_vol.setRegion(description, x_comp.regionStr());
0186     c_vol.setLimitSet(description, x_comp.limitsStr());
0187     c_vol.setVisAttributes(description, x_comp.visStr());
0188
0189     pcb_feb_ext = x_comp.offset();
0190
0191     pv = m_vol.placeVolume(c_vol,
0192                            Position(0, -pcb_feb_ext / 2.0, thickness_sum + comp_thickness / 2.0));
0193
0194     if (x_comp.isSensitive()) {
0195       // ***** SENSITIVE VOLUME
0196       if (nSensitives >= 5) {
0197         sensitiveVolumeSet = -1;
0198         break;
0199       }
0200       pv.addPhysVolID("sensor", sensor_number);
0201       // StripID. Single Sensitive Volume?
0202       if (c_nam == "DriftGap") {
0203         if (nSensitives != 0) {
0204           sensitiveVolumeSet = -1;
0205           break;
0206         }
0207         sensitiveVolumeSet = 1;
0208       }
0209       int strip_id = x_comp.key();
0210       pv.addPhysVolID("strip", strip_id);
0211       c_vol.setSensitiveDetector(sens);
0212       sensitives.push_back(pv);
0213
0214       // -------- create a measurement plane for the tracking surface attached to the sensitive volume -----
0215       Vector3D u(-1., 0., 0.);
0216       Vector3D v(0., -1., 0.);
0217       Vector3D n(0., 0., 1.);
0218
0219       // Compute the inner (i.e. thickness until mid-sensitive-volume) and
0220       //             outer (from mid-sensitive-volume to top)
0221       // thicknesses that need to be assigned to the tracking surface
0222       // depending on wether the support is above or below the sensor.
0223       double inner_thickness, outer_thickness;
0224       if (sensitiveVolumeSet == 1) {
0225         inner_thickness = thickness_so_far + comp_thickness / 2;
0226         outer_thickness = total_thickness - thickness_so_far - comp_thickness / 2;
0227       } else if (nSensitives == 0) {
0228         inner_thickness = thickness_so_far + comp_thickness / 2;
0229         outer_thickness = comp_thickness / 2;
0230       } else if (nSensitives == 4) {
0231         inner_thickness = comp_thickness / 2;
0232         outer_thickness = total_thickness - thickness_so_far - comp_thickness / 2;
0233       } else {
0234         inner_thickness = outer_thickness = comp_thickness / 2;
0235       }
0236       printout(DEBUG, "BarrelPlanarMPGDTracker_geo", "Sensitive surface @ R = %.4f (%.4f,%.4f) cm",
0237                (thickness_sum + comp_thickness / 2) / cm, inner_thickness / cm,
0238                outer_thickness / cm);
0239
0240       SurfaceType type(rec::SurfaceType::Sensitive);
0241
0242       VolPlane surf(c_vol, type, inner_thickness, outer_thickness, u, v, n); //,o ) ;
0243       volplane_surfaces.push_back(surf);
0244       nSensitives++;
0245     }
0246     thickness_sum += comp_thickness;
0247     thickness_so_far += comp_thickness;
0248   }
0249   if (sensitiveVolumeSet < 0 || sensitiveVolumeSet != nSensitives) {
0250     printout(ERROR, "BarrelPlanarMPGDTracker_geo",
0251              "Invalid set of Sensitive Volumes: it's either one (named \"DriftGap\") or 5");
0252     throw runtime_error("Logics error in building modules.");
0253   }
0254   // Now add-on the frame
0255   if (x_mod.hasChild(_U(frame))) {
0256     xml_comp_t m_frame     = x_mod.child(_U(frame));
0257     double frame_thickness = getAttrOrDefault<double>(m_frame, _U(thickness), total_thickness);
0258
0259     Box lframe_box{m_frame.width() / 2.0, (max_component_length + 2.0 * m_frame.width()) / 2.0,
0260                    frame_thickness / 2.0};
0261     Box rframe_box{m_frame.width() / 2.0, (max_component_length + 2.0 * m_frame.width()) / 2.0,
0262                    frame_thickness / 2.0};
0263     Box tframe_box{max_component_width / 2.0, m_frame.width() / 2.0, frame_thickness / 2.0};
0264     Box bframe_box{max_component_width / 2.0, m_frame.width() / 2.0, frame_thickness / 2.0};
0265
0266     // Keep track of frame with so we can adjust the module bars appropriately
0267
0268     Volume lframe_vol{"left_frame", lframe_box, description.material(m_frame.materialStr())};
0269     Volume rframe_vol{"right_frame", rframe_box, description.material(m_frame.materialStr())};
0270     Volume tframe_vol{"top_frame", tframe_box, description.material(m_frame.materialStr())};
0271     Volume bframe_vol{"bottom_frame", bframe_box, description.material(m_frame.materialStr())};
0272
0273     lframe_vol.setVisAttributes(description, m_frame.visStr());
0274     rframe_vol.setVisAttributes(description, m_frame.visStr());
0275     tframe_vol.setVisAttributes(description, m_frame.visStr());
0276     bframe_vol.setVisAttributes(description, m_frame.visStr());
0277
0278     printout(DEBUG, "BarrelPlanarMPGDTracker_geo", "frame_thickness: %f", frame_thickness);
0279     printout(DEBUG, "BarrelPlanarMPGDTracker_geo", "total_thickness: %f", total_thickness);
0280     printout(DEBUG, "BarrelPlanarMPGDTracker_geo", "frame_thickness + total_thickness: %f",
0281              frame_thickness + total_thickness);
0282
0283     m_vol.placeVolume(
0284         lframe_vol, Position(frame_width / 2.0 + max_component_width / 2, 0.0,
0285                              frame_thickness / 2.0 - total_thickness / 2.0 - gas_thickness / 2.0));
0286     m_vol.placeVolume(
0287         rframe_vol, Position(-frame_width / 2.0 - max_component_width / 2.0, 0.0,
0288                              frame_thickness / 2.0 - total_thickness / 2.0 - gas_thickness / 2.0));
0289     m_vol.placeVolume(
0290         tframe_vol, Position(0.0, frame_width / 2.0 + max_component_length / 2,
0291                              frame_thickness / 2.0 - total_thickness / 2.0 - gas_thickness / 2.0));
0292     m_vol.placeVolume(
0293         bframe_vol, Position(0.0, -frame_width / 2.0 - max_component_length / 2.0,
0294                              frame_thickness / 2.0 - total_thickness / 2.0 - gas_thickness / 2.0));
0295   }
0296
0297   // ********** LAYER
0298   // ***** ONE AND ONLY ONE LAYER
0299   xml_coll_t li(x_det, _U(layer));
0300   if (li.size() != 1) {
0301     printout(ERROR, "BarrelPlanarMPGDTracker_geo", "Number of layers = %d. Must be = 1",
0302              (int)li.size());
0303     throw runtime_error("Logics error in building modules.");
0304   }
0305   // build the layers the modules will be arranged around
0306   xml_comp_t x_layer            = li;
0307   xml_comp_t x_layout           = x_layer.child(_U(rphi_layout));
0308   xml_comp_t z_layout           = x_layer.child(_U(z_layout));
0309   int lay_id                    = x_layer.id();
0310   string lay_nam                = det_name + _toString(x_layer.id(), "_layer%d");
0311   xml_comp_t envelope_tolerance = x_layer.child(_Unicode(envelope_tolerance), false);
0312   double envelope_r_min         = 0;
0313   double envelope_r_max         = 0;
0314   double envelope_z_min         = 0;
0315   double envelope_z_max         = 0;
0316   if (envelope_tolerance) {
0317     envelope_r_min = getAttrOrDefault(envelope_tolerance, _Unicode(r_min), 0);
0318     envelope_r_max = getAttrOrDefault(envelope_tolerance, _Unicode(r_max), 0);
0319     envelope_z_min = getAttrOrDefault(envelope_tolerance, _Unicode(z_min), 0);
0320     envelope_z_max = getAttrOrDefault(envelope_tolerance, _Unicode(z_max), 0);
0321   }
0322
0323   double phi0     = x_layout.phi0();     // starting phi of first module
0324   double phi_tilt = x_layout.phi_tilt(); // Phi tilt of module
0325   double rc       = x_layout.rc();       // Radius of the module
0326   int nphi        = x_layout.nphi();     // Number of modules in phi
0327   double rphi_dr  = x_layout.dr();       // The delta radius of every other module
0328   double phi_incr = (2 * M_PI) / nphi;   // Phi increment for one module
0329   double phic     = phi0;                // Phi of the module
0330   int nz          = 2;                   // Number of modules placed in z
0331   double z_dr     = z_layout.dr();       // Radial offest of modules in z
0332   double z0       = z_layout.z0();       // Sets how much overlap in z the nz modules have
0333
0334   Assembly layer_assembly(lay_nam);
0335   Volume module_env = *volume;
0336   DetElement lay_elt(sdet, lay_nam, lay_id);
0337   auto& layerParams =
0338       DD4hepDetectorHelper::ensureExtension<dd4hep::rec::VariantParameters>(lay_elt);
0339
0340   pv = assembly.placeVolume(layer_assembly);
0341   pv.addPhysVolID("layer", lay_id);
0342   lay_elt.setPlacement(pv);
0343
0344   int module = 0;
0345   // loop over the modules in phi
0346   for (int ii = 0; ii < nphi; ii++) {
0347     double xc = rc * std::cos(phic);              // Basic x position of module
0348     double yc = rc * std::sin(phic);              // Basic y position of module
0349     double dx = z_dr * std::cos(phic + phi_tilt); // Delta x of module position
0350     double dy = z_dr * std::sin(phic + phi_tilt); // Delta y of module position
0351     // loop over the modules in z
0352     for (int j = 0; j < nz; j++) {
0353       string module_name = _toString(module, "module%02d");
0354       DetElement mod_elt(lay_elt, module_name, module);
0355       double mod_z       = 0.5 * dimensions.length();
0356       double z_placement = mod_z - 0.5 * pcb_feb_ext -
0357                            j * (nz * mod_z - pcb_feb_ext); // z location for module placement
0358       double z_offset =
0359           z_placement > 0 ? -z0 / 2.0
0360                           : z0 / 2.0; // determine the amount of overlap in z the z nz modules have
0361
0362       Transform3D tr(
0363           RotationZYX(0, ((M_PI / 2) - phic - phi_tilt), -M_PI / 2) * RotationZ(j * M_PI),
0364           Position(
0365               xc, yc,
0366               mpgd_pos.z() + z_placement +
0367                   z_offset)); //RotZYX rotates planes around azimuth, RotZ flips plane so pcb_feb_ext is facing endcaps
0368
0369       pv = layer_assembly.placeVolume(module_env, tr);
0370       pv.addPhysVolID("module", module);
0371       mod_elt.setPlacement(pv);
0372       for (int iSensitive = 0; iSensitive < sensitiveVolumeSet; iSensitive++) {
0373         // ***** SENSITIVE COMPONENTS
0374         PlacedVolume& sens_pv = sensitives[iSensitive];
0375         int de_id             = nphi * nz * iSensitive + module;
0376         DetElement comp_de(mod_elt,
0377                            std::string("de_") + sens_pv.volume().name() + _toString(de_id, "%02d"),
0378                            de_id);
0379         comp_de.setPlacement(sens_pv);
0380         auto& comp_de_params =
0381             DD4hepDetectorHelper::ensureExtension<dd4hep::rec::VariantParameters>(comp_de);
0382         comp_de_params.set<string>("axis_definitions", "XYZ");
0383         volSurfaceList(comp_de)->push_back(volplane_surfaces[iSensitive]);
0384       }
0385       // increas module counter
0386       module++;
0387       // adjust x and y coordinates
0388       xc += dx;
0389       yc += dy;
0390     }
0391     // increment counters
0392     phic += phi_incr;
0393     rc += rphi_dr;
0394   }
0395   layer_assembly->GetShape()->ComputeBBox();
0396   layerParams.set<double>("envelope_r_min", envelope_r_min);
0397   layerParams.set<double>("envelope_r_max", envelope_r_max);
0398   layerParams.set<double>("envelope_z_min", envelope_z_min);
0399   layerParams.set<double>("envelope_z_max", envelope_z_max);
0400
0401   for (xml_coll_t lmat(x_layer, _Unicode(layer_material)); lmat; ++lmat) {
0402     xml_comp_t x_layer_material = lmat;
0403     DD4hepDetectorHelper::xmlToProtoSurfaceMaterial(x_layer_material, layerParams,
0404                                                     "layer_material");
0405   }
0406   sdet.setAttributes(description, assembly, x_det.regionStr(), x_det.limitsStr(), x_det.visStr());
0407   assembly.setVisAttributes(description.invisible());
0408   pv = description.pickMotherVolume(sdet).placeVolume(assembly);
0409   pv.addPhysVolID("system", det_id); // Set the subdetector system ID
0410   sdet.setPlacement(pv);
0411
0412 #ifdef DEBUG_BarrelPlanarMPGDTracker
0413   // Reset initial print level before exiting
0414   setPrintLevel(priorPrintLevel);
0415 #endif
0416
0417   return sdet;
0418 }
0419
0420 //@}
0421 // clang-format off
0422 DECLARE_DETELEMENT(epic_OuterMPGDBarrel, create_BarrelPlanarMPGDTracker_geo)