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