epic/src/MPGDCylinderBarrelTracker_geo.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022-2024 Whitney Armstrong, Nivedith Ramasubramanian, Yann Bedfer
0003
0004 /** \addtogroup Trackers Trackers
0005  * \brief Type: **Barrel Cylinder MPGD with frames**.
0006  * \author Nivedith Ramasubramanian
0007  * Modified by M. Posik, Y. Bedfer
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 "DD4hepDetectorHelper.h"
0022
0023 using namespace std;
0024 using namespace dd4hep;
0025 using namespace dd4hep::rec;
0026 using namespace dd4hep::detail;
0027
0028 #include "Math/Vector2D.h"
0029 using ROOT::Math::XYVector;
0030
0031 /** Micromegas Barrel Tracker with space frame
0032  *
0033  * - Designed to process "mpgd_barrel.xml".
0034  *
0035  * - Derived from "BarrelTrackerWithFrame_geo.cpp".
0036  *
0037  * - "support" tag not addressed.
0038  *
0039  * - "frame" tag within the module element.
0040  *
0041  * - Several models of modules, with each a distinct radius of curvature
0042  *  but a single XML <module> and a single <layer>.
0043  *
0044  * - Two distinct versions of ".xml" are covered:
0045  *  I) Single Sensitive Volume (which name is then "GasGap"): "eicrecon" is to
0046  *   be executed while setting bit 0x1 of option "MPGD:SiFactoryPattern".
0047  * II) Multiple Sensitive Volume: Bit 0x1 of above-mentioned option not set.
0048  *
0049  * \code
0050  * \endcode
0051  *
0052  *
0053  * @author Yann Bedfer
0054  */
0055 static Ref_t create_MPGDCylinderBarrelTracker(Detector& description, xml_h e,
0056                                               SensitiveDetector sens) {
0057   xml_det_t x_det = e;
0058   Material air    = description.air();
0059   int det_id      = x_det.id();
0060   string det_name = x_det.nameStr();
0061   DetElement sdet(det_name, det_id);
0062
0063   vector<Volume> volumes;
0064   // Sensitive volumes and associated surfaces
0065   // - There can be either one or five.
0066   int sensitiveVolumeSet = 5; // 1: single volume, 5: 5 volumes, -1: error
0067   vector<PlacedVolume> sensitives[5];
0068   vector<VolPlane> volplane_surfaces[5];
0069
0070   PlacedVolume pv;
0071
0072   //#define DEBUG_MPGDCylinderBarrelTracker
0073 #ifdef DEBUG_MPGDCylinderBarrelTracker
0074   // TEMPORARILY INCREASE VERBOSITY level for debugging purposes
0075   PrintLevel priorPrintLevel = printLevel();
0076   setPrintLevel(DEBUG);
0077 #endif
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   Assembly assembly(det_name);
0092
0093   sens.setType("tracker");
0094
0095   // ********** MODULE
0096   // ***** ONE AND ONLY ONE MODULE
0097   xml_coll_t modules(x_det, _U(module));
0098   if (modules.size() != 1) {
0099     // Present detector constructor can only handle ONE <module> tag
0100     printout(ERROR, "MPGDCylinderBarrelTracker", "Number of modules = %u. Must be = 1",
0101              modules.size());
0102     throw runtime_error("Logics error in building modules.");
0103   }
0104   xml_comp_t x_mod   = modules;
0105   string m_nam       = x_mod.nameStr();
0106   double stave_width = dimensions.width(), det_length = dimensions.length();
0107   printout(DEBUG, "MPGDCylinderBarrelTracker", "Module \"%s\" width = %.2f cm length = %.2f cm",
0108            m_nam.c_str(), stave_width / cm, det_length / cm);
0109
0110   // ********** LAYER
0111   // ***** ONE AND ONLY ONE LAYER
0112   xml_coll_t li(x_det, _U(layer));
0113   if (li.size() != 1) {
0114     printout(ERROR, "MPGDCylinderBarrelTracker", "Number of layers = %d. Must be = 1",
0115              (int)li.size());
0116     throw runtime_error("Logics error in building modules.");
0117   }
0118   // ***** RETRIEVE PARAMETERS
0119   xml_comp_t x_layer = li;
0120   int lay_id         = x_layer.id();
0121   if (x_layer.moduleStr() != m_nam) {
0122     printout(ERROR, "MPGDCylinderBarrelTracker", "Layer \"%s\" does not match module \"%s\"",
0123              x_layer.moduleStr().c_str(), m_nam.c_str());
0124     throw runtime_error("Logics error in building layer.");
0125   }
0126   xml_comp_t x_barrel  = x_layer.child(_U(barrel_envelope));
0127   double barrel_length = x_barrel.z_length();
0128   double barrel_z0     = getAttrOrDefault(x_barrel, _U(z0), 0.);
0129   // ***** LAYOUTS
0130   xml_comp_t x_layout = x_layer.child(_U(rphi_layout));
0131   double phi0         = x_layout.phi0(); // Starting phi of first module.
0132   xml_comp_t z_layout = x_layer.child(_U(z_layout));
0133   double z_gap        = z_layout.gap();
0134   // ***** INVALID LAYOUT PROPERTIES
0135   // # of staves (along phi) and sections (along z) are now derived from stave
0136   // width and section length. Used to be specified directly. In order to remind
0137   // the user this is no longer the case, let's forbid the use of the
0138   // corresponding (and a few more) tags.
0139   const int nInvalids                     = 4;
0140   const xml::Tag_t unvalidTags[nInvalids] = {_U(phi_tilt), _U(nphi), _U(rc), _U(dr)};
0141   for (int uv = 0; uv < nInvalids; uv++) {
0142     if (x_barrel.hasChild(unvalidTags[uv])) {
0143       const string tag = _U(nphi);
0144       printout(ERROR, "MPGDCylinderBarrelTracker",
0145                "Layer \"%s\": Invalid property \"%s\" in \"rphi_layout\"", m_nam.c_str(),
0146                tag.c_str());
0147       throw runtime_error("Logics error in building modules.");
0148     }
0149   }
0150
0151   // ***** MODELS
0152   // Model = set of two radii of curvature used alternatively for staves
0153   //       + an offset to be applied radially.
0154   // - The offset may be null if the two radii are different enough that
0155   //  staves adjacent along phi do not overlap.
0156   // - There needs to be 2, as of 2024/03: Inner and outer.
0157   // StaveModel = radius of curvature (=> "nphi", given "rmin")
0158   //            + offset
0159   // section2Models[2][2]: Which model for [inner,outer][r1,r2]
0160   typedef struct {
0161     string name;
0162     double rmin;
0163     // "rsensor" is supposed to have been set equal to the radius specified in
0164     // the segmentation and is used here to double-check the consistency
0165     // between the stack of module components and that segmentation radius.
0166     double rsensor;
0167     int nphi;
0168     double offset;
0169     int io;
0170   } StaveModel;
0171   const int nSections = 4;
0172   StaveModel staveModels[nSections];
0173   int section2Models[2][2];
0174   xml_coll_t mi(x_mod, _U(model));
0175   if (mi.size() != 2) {
0176     printout(ERROR, "MPGDCylinderBarrelTracker", "Number of models = %d. Must be = 2",
0177              (int)mi.size());
0178     throw runtime_error("Logics error in building modules.");
0179   }
0180   int nStaveModels;
0181   for (nStaveModels = 0; mi; ++mi) {
0182     xml_comp_t x_model = mi;
0183     double rmin1 = x_model.rmin1(), rmin2 = x_model.rmin2();
0184     double rsensor = x_model.attr<double>(_Unicode(rsensor));
0185     // Determine "nphi" from stave width and radius
0186     int nphi  = int(2 * M_PI * rmin1 / stave_width + .5),
0187         nphi2 = int(2 * M_PI * rmin2 / stave_width + .5);
0188     if (nphi2 != nphi) {
0189       printout(ERROR, "MPGDCylinderBarrelTracker",
0190                "Model \"%s\": rmin1,2 = %.2f,%.2f cm are incompatible", x_model.nameStr().c_str(),
0191                rmin1, rmin2);
0192       throw runtime_error("Logics error in building modules.");
0193     }
0194     double offset          = x_model.offset();
0195     StaveModel& staveModel = staveModels[nStaveModels++];
0196     staveModel.rmin        = rmin1;
0197     staveModel.rsensor     = rsensor;
0198     staveModel.nphi        = nphi;
0199     staveModel.offset      = offset;
0200     int io;
0201     if (nStaveModels == 1) // 1st encountered "x_model" => ...
0202       io = 0;              // ...it's the inner one => io = 0.
0203     else
0204       io = 1;
0205     staveModel.io         = io;
0206     section2Models[io][0] = nStaveModels - 1;
0207     if (fabs(rmin2 - rmin1) > 1.e6) {
0208       StaveModel& staveMode2 = staveModels[nStaveModels++];
0209       staveMode2.rmin        = rmin2;
0210       staveModel.rsensor     = rsensor + rmin2 - rmin1;
0211       staveModel.nphi        = nphi;
0212       staveMode2.offset      = offset;
0213       staveMode2.io          = io;
0214       staveModel.name        = x_model.nameStr() + string("_1");
0215       staveMode2.name        = x_model.nameStr() + string("_2");
0216     } else
0217       staveModel.name = x_model.nameStr();
0218     section2Models[io][1] = nStaveModels - 1;
0219   }
0220   printout(DEBUG, "MPGDCylinderBarrelTracker", "%d Stave Models:", nStaveModels);
0221   for (int iSM = 0; iSM < nStaveModels; iSM++) {
0222     StaveModel& sM = staveModels[iSM];
0223     printout(DEBUG, "MPGDCylinderBarrelTracker",
0224              "Stave Model #%d \"%s\": rmin = %.2f cm offset = ±%.2f cm %s", iSM, sM.name.c_str(),
0225              sM.rmin, sM.offset / 2, sM.io ? "Outer" : "Inner");
0226   }
0227
0228   // ***** FRAMES
0229   // There must be two:
0230   // - Outward frame (wider, because supporting connectors)
0231   // - Otherwise frame.
0232   // Widest is taken as outward frame.
0233   typedef struct {
0234     string name;
0235     double width;
0236     string material;
0237     string vis;
0238   } Frame;
0239   Frame frames[2];
0240   xml_coll_t fi(x_mod, _U(frame));
0241   if (fi.size() != 2) {
0242     printout(ERROR, "MPGDCylinderBarrelTracker", "Number of frames = %d. Must be = 2",
0243              (int)fi.size());
0244     throw runtime_error("Logics error in building modules.");
0245   }
0246   printout(DEBUG, "MPGDCylinderBarrelTracker", "2 Frames:");
0247   int iFr;
0248   for (iFr = 0; fi; ++fi, iFr++) {
0249     xml_comp_t x_frame = fi;
0250     string name        = x_frame.nameStr();
0251     double width       = x_frame.width();
0252     string material = x_frame.materialStr(), vis = x_frame.visStr();
0253     Frame& frame   = frames[iFr];
0254     frame.name     = name;
0255     frame.width    = width;
0256     frame.material = material;
0257     frame.vis      = vis;
0258   }
0259   if (frames[0].width < frames[1].width) // Outward = widest must be first
0260     swap(frames[0], frames[1]);
0261   for (iFr = 0; iFr < 2; iFr++) {
0262     Frame& frame = frames[iFr];
0263     printout(DEBUG, "MPGDCylinderBarrelTracker",
0264              "Frame #%d \"%s\": width = %.2f cm material = \"%s\" vis = \"%s\"", iFr,
0265              frame.name.c_str(), frame.width, frame.material.c_str(), frame.vis.c_str());
0266   }
0267
0268   // ***** SERVICE
0269   xml_coll_t si(x_mod, _Unicode(service));
0270   if (si.size() != 1) {
0271     printout(ERROR, "MPGDCylinderBarrelTracker", "Number of services = %d. Must be = 1",
0272              (int)si.size());
0273     throw runtime_error("Logics error in building modules.");
0274   }
0275   xml_comp_t x_service     = si;
0276   double service_thickness = x_service.thickness();
0277   printout(DEBUG, "MPGDCylinderBarrelTracker",
0278            "1 Service \"%s\": thickness = %.4f cm, material = \"%s\"", x_service.nameStr().c_str(),
0279            service_thickness, x_service.materialStr().c_str());
0280
0281   // ***** TOTAL THICKNESS from components (used to build frames)
0282   double total_thickness = 0;
0283   xml_coll_t ci(x_mod, _U(module_component));
0284   for (ci.reset(), total_thickness = 0.0; ci; ++ci) {
0285     const xml_comp_t x_comp = ci;
0286     printout(DEBUG, "MPGDCylinderBarrelTracker", "\"%s\": \t total_thickness %.4f cm",
0287              x_comp.nameStr().c_str(), total_thickness / cm);
0288     total_thickness += x_comp.thickness();
0289   }
0290   printout(DEBUG, "MPGDCylinderBarrelTracker", " => total_thickness %.4f cm", total_thickness / cm);
0291
0292   // Now that we know total thickness, let's "printout" the characteristics
0293   // of the models (extrema, phi superposition).
0294   printout(DEBUG, "MPGDCylinderBarrelTracker", "2 Section Models:");
0295   double outerPhiSuperpos = 0; // Later used to define phi range of service to inner
0296   for (int io = 0; io < 2; io++) {
0297     StaveModel& sM0 = staveModels[section2Models[io][0]];
0298     StaveModel& sM1 = staveModels[section2Models[io][1]];
0299     XYVector vOffset(sM0.offset / 2, 0);
0300     double phiEdge0 = stave_width / 2 / sM0.rmin;
0301     XYVector vEdge0(sM0.rmin * cos(phiEdge0), sM0.rmin * sin(phiEdge0));
0302     vEdge0 -= vOffset;
0303     // Module0: Stave model has smaller curvature than circle centred on
0304     // beam axis. => Minimum is in the middle.
0305     double RMn = sM0.rmin - sM0.offset / 2, Mag0 = sqrt(vEdge0.Mag2());
0306     double phiEdge1 = stave_width / 2 / sM1.rmin;
0307     XYVector vEdge1(sM1.rmin * cos(phiEdge1), sM1.rmin * sin(phiEdge1));
0308     vEdge1 += vOffset;
0309     // Module1: Stave model has larger curvature than circle centred on
0310     // beam axis. => Maximum is in the middle.
0311     double RMx = sM1.rmin + sM1.offset / 2, Mag1 = sqrt(vEdge1.Mag2());
0312     double dPhi = acos(vEdge0.Dot(vEdge1) / Mag0 / Mag1);
0313     RMx += total_thickness;
0314     if (io == 1) {
0315       RMn -= service_thickness; // Outer section: account for services to inner section
0316       outerPhiSuperpos = dPhi;
0317     }
0318     printout(
0319         DEBUG, "MPGDCylinderBarrelTracker",
0320         "Sector Model #%d,\"%s\" = \"%s\"+\"%s\": phi overlap = %.1f deg, Extrema R = %.2f,%.2f",
0321         io, io ? "Outer" : "Inner", sM0.name.c_str(), sM1.name.c_str(), dPhi / M_PI * 180, RMn,
0322         RMx);
0323   }
0324
0325   // ********** LOOP OVER STAVE MODELS
0326   double total_length = 0; // Total length including frames
0327   for (int iSM = 0; iSM < nStaveModels; iSM++) {
0328     // "sensor_number" = Bit field in cellID identifying the sensitive surface.
0329     // It is used to set up the MultiSegmentation discrimination scheme
0330     // catering for the distinct radii of the inner/outer sectors.
0331     int sensor_number      = iSM;
0332     StaveModel& staveModel = staveModels[iSM];
0333     // phi range, when excluding frames
0334     double stave_rmin = staveModel.rmin;
0335     double dphi       = stave_width / stave_rmin / 2;
0336     double phi_start = -dphi, phi_end = dphi;
0337
0338     // ***** ASSEMBLY VOLUME: ONE PER STAVE MODEL
0339     Assembly m_vol(staveModel.name);
0340     volumes.push_back(m_vol);
0341     m_vol.setVisAttributes(description.visAttributes(x_mod.visStr()));
0342
0343     // Stave frames
0344     double zthickness, stave_rmax = stave_rmin + total_thickness;
0345     Frame &out_frame = frames[0], &frame = frames[1];
0346     total_length = det_length + out_frame.width + frame.width;
0347     // Outward frame
0348     zthickness = out_frame.width;
0349     Tube frame_tube_O(stave_rmin, stave_rmax, zthickness / 2, phi_start, phi_end);
0350     Volume frame_vol_O(out_frame.name, frame_tube_O, description.material(out_frame.material));
0351     m_vol.placeVolume(frame_vol_O, Position(0, 0, -(total_length - zthickness) / 2));
0352     // Inward frame
0353     zthickness = frame.width; // Update "zthickness"
0354     Tube frame_tube_I(stave_rmin, stave_rmax, zthickness / 2, phi_start, phi_end);
0355     Volume frame_vol_I(frame.name, frame_tube_I, description.material(frame.material));
0356     m_vol.placeVolume(frame_vol_I, Position(0, 0, (total_length - zthickness) / 2));
0357
0358     // Start/stop frames
0359     double frame_dphi =
0360         zthickness / stave_rmin; //converting the thickness of the frame to angular radians.
0361     Tube frame_tube_3(stave_rmin, stave_rmax, total_length / 2, phi_start - frame_dphi, phi_start);
0362     const string start_frame_nam("StartFrame");
0363     Volume frame_vol_3(start_frame_nam, frame_tube_3, description.material(frame.material));
0364     m_vol.placeVolume(frame_vol_3);
0365
0366     Tube frame_tube_4(stave_rmin, stave_rmax, total_length / 2, phi_end, phi_end + frame_dphi);
0367     const string stop_frame_nam("StopFrame");
0368     Volume frame_vol_4(stop_frame_nam, frame_tube_4, description.material(frame.material));
0369     m_vol.placeVolume(frame_vol_4);
0370
0371     frame_vol_O.setVisAttributes(description, out_frame.vis);
0372     frame_vol_I.setVisAttributes(description, frame.vis);
0373     frame_vol_3.setVisAttributes(description, frame.vis);
0374     frame_vol_4.setVisAttributes(description, frame.vis);
0375
0376     // ***** OUTER SECTORS: SERVICES TO INNER
0377     // Don't know what the purpose of the "(inner|outer)_thickness" arg.s
0378     // to the surface volume (see infra) and whether "inner_thickness" has to
0379     // include the extra thickness corresponding to the services.
0380     // Note: The phi range of the service volume is set so that it's smaller
0381     // than that of the rest of components. This, in order to avoid the
0382     // addition of a thickness of service to the already thick superposition
0383     // of 2 module thicknesses (or module+frame) on the edge.
0384     if (staveModel.io == 1) {
0385       // Superposition along Z
0386       double outerPos  = barrel_length / 2 - total_length / 2;
0387       double innerPos  = (total_length + z_gap) / 2;
0388       double zSuperpos = total_length - outerPos + innerPos;
0389       // Parameters
0390       double serv_thickness = service_thickness;
0391       double serv_rmin      = stave_rmin - service_thickness;
0392       double serv_length    = total_length - zSuperpos;
0393       // phi range: stay away from phi overlap and frame, add 1mm margin
0394       double dPhi       = outerPhiSuperpos / 2 + (frame.width + .1) / stave_rmin;
0395       double serv_start = phi_start + dPhi, serv_stop = phi_end - dPhi;
0396       Tube c_tube(serv_rmin, serv_rmin + serv_thickness, serv_length / 2, serv_start, serv_stop);
0397       Volume c_vol(x_service.nameStr(), c_tube, description.material(x_service.materialStr()));
0398       pv = m_vol.placeVolume(c_vol, Position(0, 0, -zSuperpos / 2));
0399       c_vol.setRegion(description, x_service.regionStr());
0400       c_vol.setLimitSet(description, x_service.limitsStr());
0401       c_vol.setVisAttributes(description, x_service.visStr());
0402     }
0403
0404     // ********** LOOP OVER COMPONENTS
0405     double comp_rmin        = stave_rmin;
0406     double thickness_so_far = 0;
0407     // Pattern of Multiple Sensitive Volumes
0408     // - The "GasGap" may be subdivided into SUBVOLUMES (this order to have one
0409     //  sensitive component per strip coordinate (and accessorily, some extras).
0410     int nSensitives = 0;
0411     for (xml_coll_t mci(x_mod, _U(module_component)); mci; ++mci) {
0412       xml_comp_t x_comp     = mci;
0413       const string c_nam    = x_comp.nameStr();
0414       double comp_thickness = x_comp.thickness();
0415       Tube c_tube(comp_rmin, comp_rmin + comp_thickness, det_length / 2, phi_start, phi_end);
0416       Volume c_vol(c_nam, c_tube, description.material(x_comp.materialStr()));
0417       pv = m_vol.placeVolume(c_vol, Position(0, 0, (out_frame.width - frame.width) / 2));
0418       c_vol.setRegion(description, x_comp.regionStr());
0419       c_vol.setLimitSet(description, x_comp.limitsStr());
0420       c_vol.setVisAttributes(description, x_comp.visStr());
0421       if (x_comp.isSensitive()) {
0422         // ***** SENSITIVE VOLUME
0423         if (nSensitives >= 5) {
0424           sensitiveVolumeSet = -1;
0425           break;
0426         }
0427         pv.addPhysVolID("sensor", sensor_number);
0428         // StripID. Single Sensitive Volume?
0429         if (c_nam == "GasGap") {
0430           if (nSensitives != 0) {
0431             sensitiveVolumeSet = -1;
0432             break;
0433           }
0434           sensitiveVolumeSet = 1;
0435         }
0436         int strip_id = x_comp.key();
0437         pv.addPhysVolID("strip", strip_id);
0438         c_vol.setSensitiveDetector(sens);
0439         sensitives[nSensitives].push_back(pv);
0440
0441         // -------- create a measurement plane for the tracking surface attached to the sensitive volume -----
0442         Vector3D u(-1., 0., 0.);
0443         Vector3D v(0., -1., 0.);
0444         Vector3D n(0., 0., 1.);
0445
0446         if (strip_id == 0) {
0447           // Consistency(+/-1um) check: segmentation consistent w/ stack of module components?
0448           double rXCheck = comp_rmin + comp_thickness / 2;
0449           if (fabs(staveModel.rsensor - rXCheck) > .0001 / cm) {
0450             printout(ERROR, "MPGDCylinderBarrelTracker",
0451                      "Sensitive Component \"%s\" of StaveModel #%d,\"%s\": rsensor(%.4f cm) != "
0452                      "radius @ sensitive surface(%.4f cm)",
0453                      c_nam.c_str(), iSM, staveModel.name.c_str(), staveModel.rsensor / cm,
0454                      rXCheck / cm);
0455             throw runtime_error("Logics error in building modules.");
0456           }
0457         }
0458
0459         // Compute the inner (i.e. thickness until mid-sensitive-volume) and
0460         //             outer (from mid-sensitive-volume to top)
0461         // thicknesses that need to be assigned to the tracking surface
0462         // depending on whether the support is above or below the sensor.
0463         double inner_thickness, outer_thickness;
0464         if (sensitiveVolumeSet == 1) {
0465           inner_thickness = thickness_so_far + comp_thickness / 2;
0466           outer_thickness = total_thickness - inner_thickness;
0467         } else if (nSensitives == 0) {
0468           inner_thickness = thickness_so_far + comp_thickness / 2;
0469           outer_thickness = comp_thickness / 2;
0470         } else if (nSensitives == 4) {
0471           inner_thickness = comp_thickness / 2;
0472           outer_thickness = total_thickness - thickness_so_far - comp_thickness / 2;
0473         } else {
0474           inner_thickness = outer_thickness = comp_thickness / 2;
0475         }
0476         printout(DEBUG, "MPGDCylinderBarrelTracker",
0477                  "Stave Model #%d,\"%s\": Sensitive surface @ R = %.4f (%.4f,%.4f) cm", iSM,
0478                  staveModel.name.c_str(), staveModel.rsensor / cm, inner_thickness / cm,
0479                  outer_thickness / cm);
0480
0481         SurfaceType type(SurfaceType::Sensitive);
0482         VolPlane surf(c_vol, type, inner_thickness, outer_thickness, u, v, n); //,o ) ;
0483         volplane_surfaces[nSensitives].push_back(surf);
0484         nSensitives++;
0485       }
0486       comp_rmin += comp_thickness;
0487       thickness_so_far += comp_thickness;
0488     } //end of module component loop
0489     if (sensitiveVolumeSet < 0 || sensitiveVolumeSet != nSensitives) {
0490       printout(ERROR, "MPGDCylinderBarrelTracker",
0491                "Invalid set of Sensitive Volumes: it's either one (named \"GasGap\") or 5");
0492       throw runtime_error("Logics error in building modules.");
0493     }
0494   } //end of stave model loop
0495
0496   // ********** LAYER
0497   // ***** ENVELOPE
0498   string lay_nam = det_name + _toString(x_layer.id(), "_layer%d");
0499   Tube lay_tub(x_barrel.inner_r(), x_barrel.outer_r(), barrel_length / 2);
0500   Volume lay_vol(lay_nam, lay_tub, air); // Create the layer envelope volume.
0501   Position lay_pos(0, 0, barrel_z0);
0502   lay_vol.setVisAttributes(description.visAttributes(x_layer.visStr()));
0503   printout(DEBUG, "MPGDCylinderBarrelTracker",
0504            "Layer \"%s\": rmin,max = %.2f,%.2f cm 1/2length = %.2f cm", lay_nam.c_str(),
0505            x_barrel.inner_r(), x_barrel.outer_r(), barrel_length / 2);
0506
0507   DetElement lay_elt(sdet, lay_nam, lay_id);
0508
0509   // the local coordinate systems of modules in dd4hep and acts differ
0510   // see http://acts.web.cern.ch/ACTS/latest/doc/group__DD4hepPlugins.html
0511   auto& layerParams =
0512       DD4hepDetectorHelper::ensureExtension<dd4hep::rec::VariantParameters>(lay_elt);
0513
0514   // ********** LOOP OVER THE SECTORS IN z
0515   // ***** SECTOR POSITIONS ALONG Z
0516   // These are the 4 central values in Z where the four sets of modules, called
0517   // sectors, will be placed.
0518   double modz_pos[nSections] = {-barrel_length / 2 + (total_length) / 2,
0519                                 -(total_length + z_gap) / 2, +(total_length + z_gap) / 2,
0520                                 +barrel_length / 2 - (total_length) / 2};
0521   int nModules               = 0;
0522   for (int iz = 0; iz < nSections; iz++) {
0523     int io                = (iz == 1 || iz == 2) ? 0 : 1;
0524     int iSMs[2]           = {section2Models[io][0], section2Models[io][1]};
0525     int iSM0              = iSMs[0];
0526     const StaveModel& sm0 = staveModels[iSM0];
0527     int nphi              = sm0.nphi;
0528     double offset         = sm0.offset;
0529     double phi_incr       = 2 * M_PI / nphi; // Phi increment for one module.
0530     // ***** LOOP OVER THE STAVES IN phi.
0531     int iphi;
0532     double phic;
0533     for (iphi = 0, phic = phi0; iphi < nphi; iphi++, phic += phi_incr, nModules++) {
0534       // Every other module...
0535       int jphi = iphi % 2, iV = iSMs[jphi];    // ...swap stave volume/sensitive
0536       double rc = (2 * jphi - 1) * offset / 2; // ...flip sign of offset
0537       if (iz >= 2)
0538         rc *= -1;    // Swap >0 and <0 offsets.
0539       double x1, y1; // Coordinates of the centre of curvature of
0540       x1                 = rc * std::cos(phic);
0541       y1                 = rc * std::sin(phic);
0542       string module_name = _toString(8 * iz + iphi, "module%02d");
0543       DetElement mod_elt(lay_elt, module_name, nModules);
0544       printout(DEBUG, "MPGDCylinderBarrelTracker",
0545                "System %d Layer \"%s\",id=%d Module \"%s\",id=%d: x,y,r: %7.4f,%7.4f,%7.4f cm",
0546                det_id, lay_nam.c_str(), lay_id, module_name.c_str(), nModules, x1 / cm, y1 / cm,
0547                sqrt(x1 * x1 + y1 * y1) / cm);
0548       RotationZYX rot;
0549       rot = RotationZYX(phic, 0, 0);
0550       if (iz >= 2) // Rotate so that outward-frame faces outwards.
0551         // Note: The product of rotations seems to have to be done in this
0552         // order (i.e. rotation by "phi" about Z times rotation by pi about Y.
0553         // This, for reason I don't fully understand. Anyway, a consequence of
0554         // that is that, in order to have iz=2,3 not out of phase w/ iz=0,1,
0555         // one has to swap >0 and <0 offsets (see instruction supra).
0556         rot *= RotationZYX(0, M_PI, 0);
0557       Transform3D tr(rot, Position(x1, y1, modz_pos[iz]));
0558       Volume& module_vol = volumes[iV];
0559       pv                 = lay_vol.placeVolume(module_vol, tr);
0560       pv.addPhysVolID("module", nModules);
0561       mod_elt.setPlacement(pv);
0562       for (int iSensitive = 0; iSensitive < sensitiveVolumeSet; iSensitive++) {
0563         // ***** SENSITIVE COMPONENTS
0564         PlacedVolume& sens_pv = sensitives[iSensitive][iV];
0565         int de_id             = nphi * nSections * iSensitive + nModules;
0566         DetElement comp_de(mod_elt,
0567                            std::string("de_") + sens_pv.volume().name() + _toString(de_id, "%02d"),
0568                            de_id);
0569         comp_de.setPlacement(sens_pv);
0570         auto& comp_de_params =
0571             DD4hepDetectorHelper::ensureExtension<dd4hep::rec::VariantParameters>(comp_de);
0572         comp_de_params.set<string>("axis_definitions", "XYZ");
0573         volSurfaceList(comp_de)->push_back(volplane_surfaces[iSensitive][iV]);
0574       }
0575     }
0576   }
0577
0578   for (xml_coll_t lmat(x_layer, _Unicode(layer_material)); lmat; ++lmat) {
0579     xml_comp_t x_layer_material = lmat;
0580     DD4hepDetectorHelper::xmlToProtoSurfaceMaterial(x_layer_material, layerParams,
0581                                                     "layer_material");
0582   }
0583
0584   // ***** CREATE THE PhysicalVolume FOR THE LAYER.
0585   pv = assembly.placeVolume(lay_vol, lay_pos); // Place layer in mother
0586   pv.addPhysVolID("layer", lay_id);            // Set the layer ID.
0587   lay_elt.setAttributes(description, lay_vol, x_layer.regionStr(), x_layer.limitsStr(),
0588                         x_layer.visStr());
0589   lay_elt.setPlacement(pv);
0590
0591   sdet.setAttributes(description, assembly, x_det.regionStr(), x_det.limitsStr(), x_det.visStr());
0592   assembly.setVisAttributes(description.invisible());
0593   pv = description.pickMotherVolume(sdet).placeVolume(assembly);
0594   pv.addPhysVolID("system", det_id); // Set the subdetector system ID.
0595   sdet.setPlacement(pv);
0596
0597 #ifdef DEBUG_MPGDCylinderBarrelTracker
0598   // Reset initial print level before exiting
0599   setPrintLevel(priorPrintLevel);
0600 #endif
0601
0602   return sdet;
0603 }
0604
0605 //@}
0606 // clang-format off
0607 DECLARE_DETELEMENT(epic_CylinderMPGDBarrel, create_MPGDCylinderBarrelTracker)