MuonSpectrometerMockupDetector/src/GeoMuonMockupExperiment.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsExamples/MuonSpectrometerMockupDetector/GeoMuonMockupExperiment.hpp"
0010
0011 #include "Acts/Utilities/Helpers.hpp"
0012 #include "Acts/Utilities/MathHelpers.hpp"
0013 #include "Acts/Utilities/StringHelpers.hpp"
0014
0015 #include <format>
0016 #include <iostream>
0017
0018 #include "GeoGenericFunctions/Variable.h"
0019 #include "GeoModelHelpers/MaterialManager.h"
0020 #include "GeoModelHelpers/defineWorld.h"
0021 #include "GeoModelHelpers/printVolume.h"
0022 #include "GeoModelIOHelpers/GMIO.h"
0023 #include "GeoModelKernel/GeoBox.h"
0024 #include "GeoModelKernel/GeoSerialTransformer.h"
0025 #include "GeoModelKernel/GeoTransform.h"
0026 #include "GeoModelKernel/GeoTrd.h"
0027 #include "GeoModelKernel/GeoTube.h"
0028 #include "GeoModelKernel/GeoXF.h"
0029 #include "GeoModelKernel/throwExcept.h"
0030
0031 using namespace Acts;
0032 namespace {
0033 constexpr double rot90deg = 90. * GeoModelKernelUnits::deg;
0034 }
0035 namespace ActsExamples {
0036
0037 std::string to_string(GeoMuonMockupExperiment::MuonLayer layer) {
0038   switch (layer) {
0039     using enum GeoMuonMockupExperiment::MuonLayer;
0040     case Inner:
0041       return "Inner";
0042     case Middle:
0043       return "Middle";
0044     case Outer:
0045       return "Outer";
0046     case nLayers:
0047       return "nLayers";
0048   }
0049   return "UNKNOWN";
0050 }
0051
0052 using FpvLink = GeoMuonMockupExperiment::FpvLink;
0053
0054 GeoMuonMockupExperiment::GeoMuonMockupExperiment(
0055     const Config& cfg, std::unique_ptr<const Acts::Logger> logger)
0056     : m_cfg{cfg}, m_logger{std::move(logger)} {}
0057 ActsPlugins::GeoModelTree GeoMuonMockupExperiment::constructMS() {
0058   const double worldR = m_cfg.barrelRadii[2] + 1. * GeoModelKernelUnits::m;
0059
0060   const double barrelZ =
0061       (m_cfg.nEtaStations + 1) * (m_chamberLength + m_cfg.stationDistInZ) +
0062       0.5 * GeoModelKernelUnits::m;
0063   const double worldZ =
0064       barrelZ + m_cfg.bigWheelDistZ + 2. * m_stationHeightEndcap;
0065
0066   PVLink world = createGeoWorld(worldR, worldR, worldZ);
0067
0068   setupMaterials();
0069
0070   m_publisher->setName("Muon");
0071
0072   const double dX = 2. * (m_cfg.barrelRadii[2] - 0.5 * m_stationHeightBarrel) *
0073                     std::sin(0.5 * m_sectorSize);
0074
0075   double outerRadius = Acts::fastHypot(
0076       m_cfg.barrelRadii[2] + 0.5 * m_stationHeightBarrel, 0.5 * dX);
0077
0078   auto barrelCylinder = make_intrusive<GeoTube>(
0079       (m_cfg.barrelRadii[0] - 0.5 * m_stationHeightBarrel), outerRadius,
0080       (m_cfg.nEtaStations + 1) * (m_chamberLength + m_cfg.stationDistInZ));
0081   auto barrelLogVol = make_intrusive<GeoLogVol>(
0082       "BarrelEnvelope", barrelCylinder,
0083       MaterialManager::getManager()->getMaterial("std::air"));
0084
0085   auto barrelEnvelope = make_intrusive<GeoPhysVol>(barrelLogVol);
0086
0087   auto toyBox = make_intrusive<GeoBox>(10. * GeoModelKernelUnits::cm,
0088                                        10. * GeoModelKernelUnits::cm,
0089                                        10. * GeoModelKernelUnits::cm);
0090   auto muonEnvelope = make_intrusive<GeoPhysVol>(make_intrusive<GeoLogVol>(
0091       "MuonEnvelope", cacheShape(toyBox),
0092       MaterialManager::getManager()->getMaterial("special::Ether")));
0093
0094   /// @brief Axis transformation from ACTS chamber frame to Geomodel chamber frame.
0095   ///        This will be used for barrel boxes (cuboids) only. In the ACTS
0096   ///        chamber frame the Z axis is along the thickness of the chamber, Y
0097   ///        along the chamber length (bending direction) and X along the
0098   ///        chamber width (tube direction)
0099   const GeoTrf::Transform3D ActsToGeomodelChamberFrame{
0100       GeoTrf::GeoRotation{rot90deg, rot90deg, 0.}};
0101   for (MuonLayer layer :
0102        {MuonLayer::Inner, MuonLayer::Middle, MuonLayer::Outer}) {
0103     for (unsigned sector = 1; sector <= m_cfg.nSectors; ++sector) {
0104       for (unsigned etaIdx = 1; etaIdx <= m_cfg.nEtaStations; ++etaIdx) {
0105         const double z_displacement =
0106             0.25 * m_chamberLength +
0107             etaIdx * (m_chamberLength + m_cfg.stationDistInZ);
0108         const double radius = m_cfg.barrelRadii[toUnderlying(layer)];
0109         barrelEnvelope->add(makeTransform(
0110             GeoTrf::TranslateZ3D(z_displacement) *
0111             GeoTrf::RotateZ3D(sector * m_sectorSize) *
0112             GeoTrf::TranslateX3D(radius) * ActsToGeomodelChamberFrame));
0113         barrelEnvelope->add(assembleBarrelStation(layer, sector, etaIdx));
0114         ///
0115         barrelEnvelope->add(
0116             makeTransform(GeoTrf::TranslateZ3D(-z_displacement) *
0117                           GeoTrf::RotateZ3D(sector * m_sectorSize) *
0118                           GeoTrf::TranslateX3D(radius) *
0119                           GeoTrf::RotateX3D(180. * GeoModelKernelUnits::deg) *
0120                           ActsToGeomodelChamberFrame));
0121         barrelEnvelope->add(
0122             assembleBarrelStation(layer, sector, -static_cast<int>(etaIdx)));
0123       }
0124     }
0125   }
0126   muonEnvelope->add(barrelEnvelope);
0127   /// Construct the endcaps
0128   if (m_cfg.buildEndcaps) {
0129     const double midWheelZ = barrelZ + 0.5 * m_stationHeightEndcap;
0130     const double outWheelZ =
0131         midWheelZ + m_stationHeightEndcap + m_cfg.bigWheelDistZ;
0132     using enum ActsExamples::GeoMuonMockupExperiment::MuonLayer;
0133     assembleBigWheel(muonEnvelope, Middle, midWheelZ);
0134     assembleBigWheel(muonEnvelope, Middle, -midWheelZ);
0135     assembleBigWheel(muonEnvelope, Outer, outWheelZ);
0136     assembleBigWheel(muonEnvelope, Outer, -outWheelZ);
0137     const double innerWheelZ = 0.8 * barrelZ;
0138     const double innerWheelR =
0139         0.95 * m_cfg.barrelRadii[toUnderlying(MuonLayer::Inner)];
0140     assembleSmallWheel(muonEnvelope, innerWheelR, innerWheelZ);
0141     assembleSmallWheel(muonEnvelope, innerWheelR, -innerWheelZ);
0142   }
0143   const unsigned nChambers =
0144       2 * m_cfg.nSectors * m_cfg.nEtaStations *
0145       static_cast<unsigned>(MuonLayer::nLayers);  // barrel part
0146   const unsigned nMultiLayers = 2 * nChambers;
0147   const unsigned nTubes = nMultiLayers * m_cfg.nTubeLayers * m_cfg.nTubes;
0148   const unsigned nRpc = 2 * nChambers * m_cfg.nRpcAlongZ * m_cfg.nRpcAlongPhi;
0149   ACTS_INFO("Constructed a muon system with "
0150             << nChambers << " muon stations containing in total "
0151             << nMultiLayers << " Mdt multilayers & " << nRpc
0152             << " Rpc chambers. Total: " << (nMultiLayers + nRpc));
0153   ACTS_INFO("Each multilayer contains "
0154             << m_cfg.nTubeLayers << " tube-layers with " << m_cfg.nTubes
0155             << " tubes each giving in total " << nTubes << " placed tubes.");
0156
0157   world->add(nameTag(m_publisher->getName()));
0158   world->add(muonEnvelope);
0159
0160   ACTS_VERBOSE("Printout of the  entire world \n " << printVolume(world));
0161
0162   clearSharedCaches();
0163   ActsPlugins::GeoModelTree outTree{};
0164   outTree.worldVolume = world;
0165
0166   using VolumeMap_t = ActsPlugins::GeoModelTree::VolumePublisher::VolumeMap_t;
0167   VolumeMap_t publishedVol{};
0168   for (const auto& [fpV, pubKey] : m_publisher->getPublishedFPV()) {
0169     try {
0170       const auto key = std::any_cast<std::string>(pubKey);
0171       if (!publishedVol
0172                .insert(std::make_pair(key, static_cast<GeoFullPhysVol*>(fpV)))
0173                .second) {
0174         throw std::invalid_argument("GeoMuonMockupExperiment() - Key " + key +
0175                                     " is no longer unique");
0176       }
0177     } catch (const std::bad_any_cast& e) {
0178       throw std::domain_error(
0179           "GeoMuonMockupExperiment() - Failed to cast the key to string " +
0180           std::string{e.what()});
0181     }
0182   }
0183   outTree.publisher->publishVolumes(m_publisher->getName(),
0184                                     std::move(publishedVol));
0185
0186   if (m_cfg.dumpTree) {
0187     // open the DB connection
0188     GMDBManager db{m_cfg.dbName};
0189     // check the DB connection
0190     if (!db.checkIsDBOpen()) {
0191       throw std::runtime_error(
0192           "GeoMuonMockupExperiment::constructMS() - It was not possible to "
0193           "open the DB correctly!");
0194     }
0195     // init the GeoModel node action
0196     GeoModelIO::WriteGeoModel writeGeoDB{db};
0197     world->exec(&writeGeoDB);  // visit all GeoModel nodes
0198     writeGeoDB.saveToDB(m_publisher.get());
0199   }
0200
0201   m_publisher.reset();
0202   return outTree;
0203 }
0204 PVLink GeoMuonMockupExperiment::assembleEndcapStation(const double lowR,
0205                                                       const MuonLayer layer,
0206                                                       const unsigned sector,
0207                                                       const int etaIdx) {
0208   const double angularScale = 2. * std::sin(0.5 * m_sectorSize);
0209
0210   const double lowTubeLength = angularScale * lowR;
0211   const double upperTubeLength = angularScale * (lowR + m_chamberLength);
0212
0213   auto envelopeTrd = make_intrusive<GeoTrd>(
0214       0.5 * m_stationHeightEndcap, 0.5 * m_stationHeightEndcap,
0215       0.5 * lowTubeLength, 0.5 * upperTubeLength, 0.5 * m_chamberLength);
0216   auto logVol = make_intrusive<GeoLogVol>(
0217       "MuonEndcapStation", cacheShape(envelopeTrd),
0218       MaterialManager::getManager()->getMaterial("std::air"));
0219   auto envelopeVol = make_intrusive<GeoPhysVol>(cacheVolume(logVol));
0220   double currentX = -envelopeTrd->getXHalfLength1() + 0.5 * m_tgcChamberHeight;
0221   if (layer == MuonLayer::Middle) {
0222     envelopeVol->add(makeTransform(GeoTrf::TranslateX3D(currentX)));
0223     publishFPV(envelopeVol, assembleTgcChamber(lowTubeLength, upperTubeLength),
0224                std::format("TGC_T1E_{:}_{:}", etaIdx, sector));
0225   }
0226   currentX +=
0227       0.5 * m_tgcChamberHeight + s_tgcMdtSeparation + 0.5 * m_multiLayerHeight;
0228
0229   envelopeVol->add(makeTransform(GeoTrf::TranslateX3D(currentX)));
0230   publishFPV(envelopeVol,
0231              assembleMultilayerEndcap(1, lowTubeLength, upperTubeLength),
0232              std::format("{}_EMDT_{}_{}_1", to_string(layer), etaIdx, sector));
0233
0234   currentX += 0.5 * m_multiLayerHeight + m_cfg.multiLayerSeparation +
0235               0.5 * m_multiLayerHeight;
0236   envelopeVol->add(makeTransform(GeoTrf::TranslateX3D(currentX)));
0237   publishFPV(envelopeVol,
0238              assembleMultilayerEndcap(2, lowTubeLength, upperTubeLength),
0239              std::format("{}_EMDT_{}_{}_2", to_string(layer), etaIdx, sector));
0240   currentX += s_tgcMdtSeparation + 0.5 * m_tgcChamberHeight;
0241   if (layer == MuonLayer::Middle) {
0242     envelopeVol->add(makeTransform(GeoTrf::TranslateX3D(currentX)));
0243     publishFPV(envelopeVol, assembleTgcChamber(lowTubeLength, upperTubeLength),
0244                std::format("TGC_T2E_{:}_{:}", etaIdx, sector));
0245     currentX += 0.5 * m_tgcChamberHeight + s_tgcChamberSeparation +
0246                 0.5 * m_tgcChamberHeight;
0247     envelopeVol->add(makeTransform(GeoTrf::TranslateX3D(currentX)));
0248     publishFPV(envelopeVol, assembleTgcChamber(lowTubeLength, upperTubeLength),
0249                std::format("TGC_T3E_{:}_{:}", etaIdx, sector));
0250   }
0251   return envelopeVol;
0252 }
0253 FpvLink GeoMuonMockupExperiment::assembleTgcChamber(const double bottomWidth,
0254                                                     const double topWidth) {
0255   auto envelopeTrd = make_intrusive<GeoTrd>(
0256       0.5 * m_tgcChamberHeight, 0.5 * m_tgcChamberHeight, 0.5 * bottomWidth,
0257       0.5 * topWidth, 0.48 * m_chamberLength);
0258   auto logVol = make_intrusive<GeoLogVol>(
0259       "TgcEnvelope", cacheShape(envelopeTrd),
0260       MaterialManager::getManager()->getMaterial("std::G10"));
0261   auto tgcPhysVol = make_intrusive<GeoFullPhysVol>(cacheVolume(logVol));
0262   double currentX =
0263       -envelopeTrd->getXHalfLength1() + 0.5 * s_tgcGasSingletSeparation;
0264   for (unsigned gap = 1; gap <= m_cfg.nTgcGasGaps; ++gap) {
0265     currentX += 0.5 * s_tgcGasHeight;
0266
0267     auto gasTrd =
0268         make_intrusive<GeoTrd>(0.5 * s_tgcGasHeight, 0.5 * s_tgcGasHeight,
0269                                0.95 * envelopeTrd->getYHalfLength1(),
0270                                0.95 * envelopeTrd->getYHalfLength2(),
0271                                0.95 * envelopeTrd->getZHalfLength());
0272     auto gasLogVol = cacheVolume(make_intrusive<GeoLogVol>(
0273         "TgcGasGap", cacheShape(gasTrd),
0274         MaterialManager::getManager()->getMaterial("std::ArCO2")));
0275
0276     tgcPhysVol->add(geoId(gap));
0277     tgcPhysVol->add(makeTransform(GeoTrf::TranslateX3D(currentX)));
0278     tgcPhysVol->add(cacheVolume(make_intrusive<GeoPhysVol>(gasLogVol)));
0279
0280     currentX += 0.5 * s_tgcGasHeight + s_tgcGasSingletSeparation;
0281   }
0282
0283   return tgcPhysVol;
0284 }
0285
0286 void GeoMuonMockupExperiment::assembleBigWheel(const PVLink& envelopeVol,
0287                                                const MuonLayer layer,
0288                                                const double wheelZ) {
0289   envelopeVol->add(makeTransform(GeoTrf::TranslateZ3D(wheelZ)));
0290   const double lowR = m_cfg.endCapWheelLowR;
0291   const double effR = (m_chamberLength + m_cfg.stationDistInR);
0292
0293   const unsigned nEta =
0294       1 +
0295       static_cast<unsigned>(
0296           (m_cfg.barrelRadii[toUnderlying(MuonLayer::Outer)] - lowR) / effR);
0297   const double highR = lowR + nEta * effR;
0298   auto envelopeShape =
0299       make_intrusive<GeoTube>(lowR, highR, 0.5 * m_stationHeightEndcap);
0300   auto envelopeLogVol = make_intrusive<GeoLogVol>(
0301       "EndcapEnvelope", cacheShape(envelopeShape),
0302       MaterialManager::getManager()->getMaterial("std::air"));
0303
0304   auto wheelEnvelope = make_intrusive<GeoPhysVol>(cacheVolume(envelopeLogVol));
0305
0306   for (unsigned stationEta = 0; stationEta < nEta; ++stationEta) {
0307     const double radius = lowR + stationEta * effR;
0308     for (unsigned sector = 1; sector <= m_cfg.nSectors; ++sector) {
0309       if (wheelZ > 0) {
0310         wheelEnvelope->add(
0311             makeTransform(GeoTrf::RotateZ3D(sector * m_sectorSize) *
0312                           GeoTrf::TranslateX3D(radius + 0.5 * m_chamberLength) *
0313                           GeoTrf::RotateY3D(90. * GeoModelKernelUnits::deg) *
0314                           GeoTrf::RotateZ3D(180. * GeoModelKernelUnits::deg)));
0315       } else {
0316         wheelEnvelope->add(
0317             makeTransform(GeoTrf::RotateZ3D(sector * m_sectorSize) *
0318                           GeoTrf::TranslateX3D(radius + 0.5 * m_chamberLength) *
0319                           GeoTrf::RotateY3D(90. * GeoModelKernelUnits::deg)));
0320       }
0321       const int castEta =
0322           copySign(1, wheelZ) * static_cast<int>(stationEta + 1);
0323       wheelEnvelope->add(assembleEndcapStation(radius, layer, sector, castEta));
0324     }
0325   }
0326   envelopeVol->add(wheelEnvelope);
0327 }
0328
0329 void GeoMuonMockupExperiment::setupMaterials() {
0330   auto* matMan = MaterialManager::getManager();
0331
0332   matMan->setMaterialNamespace("std");
0333   ACTS_DEBUG("Create the chemical elements.");
0334   /// Table taken from
0335   /// https://gitlab.cern.ch/atlas/geomodelatlas/GeoModelData/-/blob/master/Materials/elements.xml
0336   matMan->addElement("Carbon", "C", 6, 12.0112);
0337   matMan->addElement("Aluminium", "Al", 13, 26.9815);
0338   matMan->addElement("Iron", "Fe", 26, 55.847);
0339   matMan->addElement("Copper", "Cu", 29, 63.54);
0340   matMan->addElement("Nitrogen", "N", 7.0, 14.0031);
0341   matMan->addElement("Oxygen", "O", 8.0, 15.9949);
0342   matMan->addElement("Argon", "Ar", 18.0, 39.9624);
0343   matMan->addElement("Hydrogen", "H", 1.0, 1.00782503081372);
0344   matMan->addElement("Chlorine", "Cl", 17.0, 35.453);
0345   matMan->addElement("Fluorine", "F", 9., 18.9984);
0346   using MatComposition_t = std::vector<std::pair<std::string, double>>;
0347
0348   auto appendMaterial = [matMan, this](const std::string& matName,
0349                                        const MatComposition_t& composition,
0350                                        const double density) {
0351     ACTS_DEBUG("Create new material " << matName << " with density "
0352                                       << density);
0353     matMan->addMaterial(matName, density);
0354     for (const auto& [compName, fraction] : composition) {
0355       ACTS_DEBUG("Try to add new material component "
0356                  << compName << " contributing by " << fraction
0357                  << " parts to the total material");
0358       matMan->addMatComponent(compName, fraction);
0359     }
0360     ACTS_DEBUG("Material defined.");
0361     matMan->lockMaterial();
0362   };
0363   appendMaterial("air",
0364                  {{"Nitrogen", 0.7494},
0365                   {"Oxygen", 0.2369},
0366                   {"Argon", 0.0129},
0367                   {"Hydrogen", 0.0008}},
0368                  0.001290);
0369   appendMaterial("Aluminium", {{"Aluminium", 1.}}, 2.7);
0370   appendMaterial("Copper", {{"Copper", 1.}}, 8.96);
0371   appendMaterial("CO2", {{"Carbon", 1.}, {"Oxygen", 2.}}, 0.00184);
0372   appendMaterial("ArCO2", {{"Argon", .93}, {"std::CO2", .07}}, .0054);
0373
0374   appendMaterial("Forex",
0375                  {{"Carbon", 0.3843626433635827},
0376                   {"Hydrogen", 0.0483830941493594},
0377                   {"Chlorine", 0.5672542624870579}},
0378                  0.7);
0379   appendMaterial("G10", {{"Carbon", 2}, {"Hydrogen", 2}, {"Fluorine", 4}}, 2);
0380
0381   if (logger().level() == Acts::Logging::Level::DEBUG) {
0382     matMan->printAll();
0383   }
0384 }
0385 void GeoMuonMockupExperiment::publishFPV(const PVLink& envelopeVol,
0386                                          const FpvLink& publishMe,
0387                                          const std::string& pubName) {
0388   m_publisher->publishNode(static_cast<GeoVFullPhysVol*>(publishMe.get()),
0389                            pubName);
0390   ACTS_DEBUG("Publish new volume " << pubName);
0391   envelopeVol->add(nameTag(pubName));
0392   envelopeVol->add(publishMe);
0393 }
0394 PVLink GeoMuonMockupExperiment::assembleBarrelStation(const MuonLayer layer,
0395                                                       const unsigned sector,
0396                                                       const int etaIdx) {
0397   /// We are working now in the ACTS chamber frame, where Z is along the
0398   /// thickness of the chamber, Y along the chamber length (bending direction)
0399   /// and X along the chamber width (tube direction)
0400   const double envelopeWidth =
0401       2. *
0402       (m_cfg.barrelRadii[toUnderlying(layer)] - 0.5 * m_stationHeightBarrel) *
0403       std::sin(0.5 * m_sectorSize);
0404
0405   auto box = make_intrusive<GeoBox>(
0406       0.5 * envelopeWidth,
0407       0.5 * m_chamberLength + 0.1 * GeoModelKernelUnits::mm,
0408       0.5 * m_stationHeightBarrel);
0409   auto logVol = make_intrusive<GeoLogVol>(
0410       "MuonBarrelLogVol", cacheShape(box),
0411       MaterialManager::getManager()->getMaterial("std::air"));
0412   auto envelopeVol = make_intrusive<GeoPhysVol>(cacheVolume(logVol));
0413
0414   /// add the rpc at doubletR = 1
0415   auto placeRpc = [&](const double currentZ, unsigned dRIdx) {
0416     const double stepdY = m_chamberLength / m_cfg.nRpcAlongZ;
0417     const double stepdX = envelopeWidth / m_cfg.nRpcAlongPhi;
0418
0419     for (unsigned dY = 0; dY < m_cfg.nRpcAlongZ; ++dY) {
0420       for (unsigned dX = 0; dX < m_cfg.nRpcAlongPhi; ++dX) {
0421         envelopeVol->add(makeTransform(GeoTrf::Translate3D(
0422             -0.5 * envelopeWidth + stepdX * (dX + 0.5),
0423             -0.5 * m_chamberLength + stepdY * (dY + 0.5), currentZ)));
0424         publishFPV(envelopeVol, assembleRpcChamber(envelopeWidth),
0425                    std::format("{:}_RPC_{:}_{:}_{:}_{:}_{:}", to_string(layer),
0426                                etaIdx, sector, dRIdx, dX, dY));
0427       }
0428     }
0429   };
0430
0431   double currentZ = -box->getZHalfLength();
0432   currentZ += 0.5 * m_rpcChamberHeight;
0433   placeRpc(currentZ, 1);
0434   currentZ +=
0435       0.5 * m_rpcChamberHeight + s_rpcMdtDistance + 0.5 * m_multiLayerHeight;
0436   envelopeVol->add(makeTransform(GeoTrf::TranslateZ3D(currentZ)));
0437   publishFPV(envelopeVol, assembleMultilayerBarrel(1, envelopeWidth),
0438              std::format("{}_BMDT_{}_{}_1", to_string(layer), etaIdx, sector));
0439   currentZ += 0.5 * m_multiLayerHeight + m_cfg.multiLayerSeparation +
0440               0.5 * m_multiLayerHeight;
0441   envelopeVol->add(makeTransform(GeoTrf::TranslateZ3D(currentZ)));
0442   publishFPV(envelopeVol, assembleMultilayerBarrel(2, envelopeWidth),
0443              std::format("{}_BMDT_{}_{}_2", to_string(layer), etaIdx, sector));
0444   currentZ +=
0445       0.5 * m_rpcChamberHeight + s_rpcMdtDistance + 0.5 * m_multiLayerHeight;
0446   placeRpc(currentZ, 2);
0447   return envelopeVol;
0448 }
0449
0450 PVLink GeoMuonMockupExperiment::assembleTube(const double tubeLength) {
0451   auto* matMan = MaterialManager::getManager();
0452
0453   auto outerTube =
0454       make_intrusive<GeoTube>(0., m_outerTubeRadius, 0.5 * tubeLength);
0455   auto outerTubeLogVol = make_intrusive<GeoLogVol>(
0456       "MdtDriftWall", outerTube, matMan->getMaterial("std::Aluminium"));
0457   auto outerTubeVol = make_intrusive<GeoPhysVol>(cacheVolume(outerTubeLogVol));
0458   /// Place the drift gas inside the outer tube
0459   auto innerTube =
0460       make_intrusive<GeoTube>(0., m_cfg.innerTubeRadius, 0.5 * tubeLength);
0461   auto innerTubeLogVol = make_intrusive<GeoLogVol>(
0462       "MDTDriftGas", innerTube, matMan->getMaterial("std::ArCO2"));
0463   outerTubeVol->add(make_intrusive<GeoPhysVol>(cacheVolume(innerTubeLogVol)));
0464   return cacheVolume(outerTubeVol);
0465 }
0466
0467 PVLink GeoMuonMockupExperiment::buildTubes(const double lowerTubeLength,
0468                                            const double upperTubeLength) {
0469   auto* matMan = MaterialManager::getManager();
0470   GeoShapePtr envShape{};
0471   if (std::abs(lowerTubeLength - upperTubeLength) <
0472       std::numeric_limits<double>::epsilon()) {
0473     envShape = make_intrusive<GeoBox>(
0474         0.5 * m_tubeLayersHeight, 0.5 * m_chamberLength, 0.5 * lowerTubeLength);
0475   } else {
0476     envShape = make_intrusive<GeoTrd>(
0477         0.5 * m_tubeLayersHeight, 0.5 * m_tubeLayersHeight,
0478         0.5 * lowerTubeLength, 0.5 * upperTubeLength, 0.5 * m_chamberLength);
0479   }
0480   auto tubeLogVol = make_intrusive<GeoLogVol>(
0481       "MdtTubeEnvelope", cacheShape(envShape), matMan->getMaterial("std::air"));
0482
0483   auto envelopeVol = make_intrusive<GeoPhysVol>(cacheVolume(tubeLogVol));
0484   /// Place the tubes inside the envelope
0485   const Acts::Vector3 posStag{
0486       m_tubePitch * std::sin(60. * GeoModelKernelUnits::deg),
0487       m_tubePitch * std::cos(60. * GeoModelKernelUnits::deg), 0.};
0488   const Acts::Vector3 negStag{
0489       m_tubePitch * std::sin(60. * GeoModelKernelUnits::deg),
0490       -m_tubePitch * std::cos(60. * GeoModelKernelUnits::deg), 0.};
0491
0492   Acts::Vector3 firstTubePos{-0.5 * m_tubeLayersHeight + 0.5 * m_tubePitch,
0493                              -0.5 * m_chamberLength + 0.5 * m_tubePitch, 0.};
0494
0495   //// Put the tube into a separate container
0496   auto toyBox = make_intrusive<GeoBox>(10. * GeoModelKernelUnits::cm,
0497                                        10. * GeoModelKernelUnits::cm,
0498                                        10. * GeoModelKernelUnits::cm);
0499   auto toyBoxLogVol = cacheVolume(
0500       make_intrusive<GeoLogVol>("TubeLayerLog", cacheShape(toyBox),
0501                                 matMan->getMaterial("special::Ether")));
0502
0503   const double dTube = (upperTubeLength - lowerTubeLength) / (m_cfg.nTubes - 1);
0504
0505   GeoGenfun::Variable K;
0506   GeoGenfun::GENFUNCTION F = K * m_tubePitch;
0507   GeoXF::TRANSFUNCTION T = GeoXF::Pow(GeoTrf::TranslateY3D(1.0), F);
0508
0509   auto barrelTubeLayer = make_intrusive<GeoSerialTransformer>(
0510       assembleTube(lowerTubeLength - 1. * GeoModelKernelUnits::cm), &T,
0511       m_cfg.nTubes);
0512
0513   for (unsigned tL = 0; tL < m_cfg.nTubeLayers; ++tL) {
0514     auto layerVol = make_intrusive<GeoPhysVol>(toyBoxLogVol);
0515     /// For endcap chambers the tube need to be placed individually
0516     if (envShape->typeID() == GeoTrd::getClassTypeID()) {
0517       envelopeVol->add(makeTransform(GeoTrf::RotateX3D(rot90deg)));
0518       for (unsigned t = 0; t < m_cfg.nTubes; ++t) {
0519         layerVol->add(makeTransform(GeoTrf::TranslateY3D(t * m_tubePitch)));
0520         layerVol->add(assembleTube(lowerTubeLength + dTube * t));
0521       }
0522     } else {
0523       /// Simple serial transformer for the barrel
0524       layerVol->add(serialId(tL));
0525       layerVol->add(barrelTubeLayer);
0526     }
0527     envelopeVol->add(makeTransform(GeoTrf::Translate3D(firstTubePos)));
0528     envelopeVol->add(cacheVolume(layerVol));
0529     firstTubePos = firstTubePos + ((tL % 2) != 1 ? posStag : negStag);
0530   }
0531   return cacheVolume(envelopeVol);
0532 }
0533
0534 PVLink GeoMuonMockupExperiment::buildAbsorber(const double thickness,
0535                                               const double widthS,
0536                                               const double widthL,
0537                                               const double length) {
0538   GeoShapePtr shape{};
0539   if (std::abs(widthS - widthL) < std::numeric_limits<double>::epsilon()) {
0540     // For geoBoxes we use the ACTS chamber frame
0541     shape = make_intrusive<GeoBox>(0.5 * widthS, 0.5 * length, 0.5 * thickness);
0542   } else {
0543     // For geoTrd we use the GeoModel chamber frame
0544     shape = make_intrusive<GeoTrd>(0.5 * thickness, 0.5 * thickness,
0545                                    0.5 * widthS, 0.5 * widthL, 0.5 * length);
0546   }
0547   auto logVol = cacheVolume(make_intrusive<GeoLogVol>(
0548       "PassiveMat", cacheShape(shape),
0549       MaterialManager::getManager()->getMaterial("std::Forex")));
0550   return cacheVolume(make_intrusive<GeoPhysVol>(logVol));
0551 }
0552
0553 FpvLink GeoMuonMockupExperiment::assembleRpcChamber(const double chamberWidth) {
0554   auto* matMan = MaterialManager::getManager();
0555   constexpr double margin = 0.98;
0556   auto rpcBox = make_intrusive<GeoBox>(
0557       0.5 * (chamberWidth / m_cfg.nRpcAlongPhi) * margin,
0558       0.5 * (m_chamberLength / m_cfg.nRpcAlongZ) * margin,
0559       0.5 * m_rpcChamberHeight);
0560   auto envLogVol = cacheVolume(make_intrusive<GeoLogVol>(
0561       "RpcChamber", cacheShape(rpcBox), matMan->getMaterial("std::Copper")));
0562   auto rpcEnvelope = make_intrusive<GeoFullPhysVol>(envLogVol);
0563   ///
0564   double currentZ = -rpcBox->getZHalfLength() + 0.5 * s_rpcGasSingletSeparation;
0565
0566   for (unsigned gap = 1; gap <= m_cfg.nRpcGasGaps; ++gap) {
0567     currentZ += 0.5 * s_rpcGasHeight;
0568
0569     auto gasBox =
0570         make_intrusive<GeoBox>(rpcBox->getXHalfLength(),
0571                                rpcBox->getYHalfLength(), 0.5 * s_rpcGasHeight);
0572     auto gasLogVol = cacheVolume(make_intrusive<GeoLogVol>(
0573         "RpcGasGap", cacheShape(gasBox), matMan->getMaterial("std::ArCO2")));
0574
0575     rpcEnvelope->add(geoId(gap));
0576     rpcEnvelope->add(makeTransform(GeoTrf::TranslateZ3D(currentZ)));
0577     rpcEnvelope->add(cacheVolume(make_intrusive<GeoPhysVol>(gasLogVol)));
0578
0579     currentZ += 0.5 * s_rpcGasHeight;
0580     if (gap == m_cfg.nRpcGasGaps) {
0581       break;
0582     }
0583     currentZ += 0.5 * s_rpcGasSingletSeparation;
0584     rpcEnvelope->add(makeTransform(GeoTrf::TranslateZ3D(currentZ)));
0585     rpcEnvelope->add(buildAbsorber(
0586         s_rpcGasSingletSeparation, 2. * rpcBox->getXHalfLength(),
0587         2. * rpcBox->getXHalfLength(), 2. * rpcBox->getYHalfLength()));
0588     currentZ += 0.5 * s_rpcGasSingletSeparation;
0589   }
0590   return rpcEnvelope;
0591 };
0592 FpvLink GeoMuonMockupExperiment::assembleMultilayerEndcap(
0593     const unsigned ml, const double lowerTubeLength,
0594     const double upperTubeLength) {
0595   auto* matMan = MaterialManager::getManager();
0596   auto envelopeTrd = make_intrusive<GeoTrd>(
0597       0.5 * m_multiLayerHeight, 0.5 * m_multiLayerHeight,
0598       0.5 * lowerTubeLength + 0.05 * GeoModelKernelUnits::mm,
0599       0.5 * upperTubeLength + 0.05 * GeoModelKernelUnits::mm,
0600       0.5 * m_chamberLength);
0601
0602   auto envelopeLogVol =
0603       make_intrusive<GeoLogVol>("MultilayerEnv", cacheShape(envelopeTrd),
0604                                 matMan->getMaterial("std::air"));
0605
0606   auto envelopeVol =
0607       make_intrusive<GeoFullPhysVol>(cacheVolume(envelopeLogVol));
0608
0609   double currentX = -envelopeTrd->getXHalfLength1();
0610   if (ml == 1 && m_cfg.mdtFoamThickness > 0.) {
0611     currentX += 0.5 * m_cfg.mdtFoamThickness;
0612     envelopeVol->add(makeTransform(GeoTrf::TranslateX3D(currentX)));
0613     envelopeVol->add(buildAbsorber(m_cfg.mdtFoamThickness, lowerTubeLength,
0614                                    upperTubeLength, m_chamberLength));
0615     currentX += 0.5 * m_cfg.mdtFoamThickness + s_mdtFoamTubeDistance;
0616   }
0617   currentX += 0.5 * m_tubeLayersHeight;
0618   envelopeVol->add(makeTransform(GeoTrf::TranslateX3D(currentX)));
0619   envelopeVol->add(buildTubes(lowerTubeLength, upperTubeLength));
0620   if (ml == 2 && m_cfg.mdtFoamThickness > 0.) {
0621     currentX += 0.5 * m_tubeLayersHeight + 0.5 * m_cfg.mdtFoamThickness +
0622                 s_mdtFoamTubeDistance;
0623     envelopeVol->add(makeTransform(GeoTrf::TranslateX3D(currentX)));
0624     envelopeVol->add(buildAbsorber(m_cfg.mdtFoamThickness, lowerTubeLength,
0625                                    upperTubeLength, m_chamberLength));
0626   }
0627   return envelopeVol;
0628 }
0629 FpvLink GeoMuonMockupExperiment::assembleMultilayerBarrel(
0630     const unsigned ml, const double tubeLength) {
0631   auto* matMan = MaterialManager::getManager();
0632
0633   const double envelopeWidth = tubeLength;
0634
0635   auto envelopeBox = make_intrusive<GeoBox>(
0636       0.5 * envelopeWidth, 0.5 * m_chamberLength, 0.5 * m_multiLayerHeight);
0637   auto envelopeLogVol =
0638       make_intrusive<GeoLogVol>("MultilayerEnv", cacheShape(envelopeBox),
0639                                 matMan->getMaterial("std::air"));
0640
0641   auto envelopeVol =
0642       make_intrusive<GeoFullPhysVol>(cacheVolume(envelopeLogVol));
0643
0644   double currentZ = -envelopeBox->getZHalfLength();
0645   if (ml == 1 && m_cfg.mdtFoamThickness > 0.) {
0646     currentZ += 0.5 * m_cfg.mdtFoamThickness;
0647     envelopeVol->add(makeTransform(GeoTrf::TranslateZ3D(currentZ)));
0648     envelopeVol->add(buildAbsorber(m_cfg.mdtFoamThickness, envelopeWidth,
0649                                    envelopeWidth, m_chamberLength));
0650     currentZ += 0.5 * m_cfg.mdtFoamThickness + s_mdtFoamTubeDistance;
0651   }
0652   PVLink tubeVol = buildTubes(envelopeWidth, envelopeWidth);
0653   currentZ += 0.5 * m_tubeLayersHeight;
0654   envelopeVol->add(makeTransform(GeoTrf::TranslateZ3D(currentZ) *
0655                                  GeoTrf::RotateY3D(-rot90deg)));
0656   envelopeVol->add(tubeVol);
0657   if (ml == 2 && m_cfg.mdtFoamThickness > 0.) {
0658     currentZ += 0.5 * m_tubeLayersHeight + 0.5 * m_cfg.mdtFoamThickness +
0659                 s_mdtFoamTubeDistance;
0660     envelopeVol->add(makeTransform(GeoTrf::TranslateZ3D(currentZ)));
0661     envelopeVol->add(buildAbsorber(m_cfg.mdtFoamThickness, envelopeWidth,
0662                                    envelopeWidth, m_chamberLength));
0663   }
0664   return envelopeVol;
0665 }
0666 PVLink GeoMuonMockupExperiment::assembleSmallWheelSector(const double wedgeL,
0667                                                          const int etaIdx,
0668                                                          const int sector) {
0669   const double angularScale = 2. * std::sin(0.5 * m_sectorSize);
0670
0671   auto envelopeTrd = make_intrusive<GeoTrd>(
0672       0.5 * m_innerWheelHeight, 0.5 * m_innerWheelHeight,
0673       0.5 * angularScale * m_innerWheelHeight,
0674       0.5 * angularScale * (m_innerWheelHeight + wedgeL), 0.5 * wedgeL);
0675   auto envelopeLog = make_intrusive<GeoLogVol>(
0676       "InnerWheelWedgeLog", cacheShape(envelopeTrd),
0677       MaterialManager::getManager()->getMaterial("std::air"));
0678
0679   auto envelopeWedge = make_intrusive<GeoPhysVol>(cacheVolume(envelopeLog));
0680
0681   double currentX = -envelopeTrd->getXHalfLength1() + 0.5 * m_innerWheelWedgeH;
0682   for (unsigned ml = 1; ml <= m_cfg.nInnerMultiplets; ++ml) {
0683     envelopeWedge->add(makeTransform(GeoTrf::TranslateX3D(currentX)));
0684
0685     auto wedgeTrd = make_intrusive<GeoTrd>(
0686         0.5 * m_innerWheelWedgeH, 0.5 * m_innerWheelWedgeH,
0687         envelopeTrd->getYHalfLength1(), envelopeTrd->getYHalfLength2(),
0688         envelopeTrd->getZHalfLength() - 1. * GeoModelKernelUnits::mm);
0689
0690     auto wedgeLogVol = make_intrusive<GeoLogVol>(
0691         "SmallWheelMultiplet", cacheShape(wedgeTrd),
0692         MaterialManager::getManager()->getMaterial("std::G10"));
0693
0694     auto wedgePhysVol =
0695         make_intrusive<GeoFullPhysVol>(cacheVolume(wedgeLogVol));
0696     publishFPV(envelopeWedge, wedgePhysVol,
0697                std::format("SmallWheel_{}_{}_{}", etaIdx, sector, ml));
0698     currentX += 0.5 * m_innerWheelWedgeH + s_swMultipletSeparation;
0699
0700     double wedgeX = -wedgeTrd->getXHalfLength1() +
0701                     0.5 * (s_swlGasGapHeight + s_swGasGapSeparation);
0702     auto gasShape = make_intrusive<GeoTrd>(
0703         0.5 * s_swlGasGapHeight, 0.5 * s_swlGasGapHeight,
0704         wedgeTrd->getYHalfLength1() - 1. * GeoModelKernelUnits::mm,
0705         wedgeTrd->getYHalfLength2() - 1. * GeoModelKernelUnits::mm,
0706         wedgeTrd->getZHalfLength() - 1. * GeoModelKernelUnits::mm);
0707     auto gasLogVol = make_intrusive<GeoLogVol>(
0708         "SmallWheelGasGap", cacheShape(gasShape),
0709         MaterialManager::getManager()->getMaterial("std::ArCO2"));
0710
0711     auto gasGapVol =
0712         cacheVolume(make_intrusive<GeoPhysVol>(cacheVolume(gasLogVol)));
0713
0714     for (unsigned int gasGap = 1; gasGap <= m_cfg.nInnerGasGapsPerMl;
0715          ++gasGap) {
0716       wedgePhysVol->add(makeTransform(GeoTrf::TranslateX3D(wedgeX)));
0717       wedgePhysVol->add(geoId(gasGap));
0718       wedgePhysVol->add(gasGapVol);
0719       wedgeX += 0.5 * (s_swlGasGapHeight + s_swGasGapSeparation);
0720     }
0721   }
0722   return envelopeWedge;
0723 }
0724 void GeoMuonMockupExperiment::assembleSmallWheel(const PVLink& envelope,
0725                                                  const double outerR,
0726                                                  const double wheelZ) {
0727   if (m_cfg.nInnerMultiplets == 0u) {
0728     ACTS_DEBUG("Small wheel will not be assembled");
0729     return;
0730   }
0731
0732   auto envelopeShape = make_intrusive<GeoTube>(m_cfg.endCapWheelLowR, outerR,
0733                                                0.5 * m_innerWheelHeight);
0734   auto envelopeVol = make_intrusive<GeoLogVol>(
0735       "InnerWheelEnvelope", cacheShape(envelopeShape),
0736       MaterialManager::getManager()->getMaterial("std::air"));
0737   auto envelopeWheel = make_intrusive<GeoPhysVol>(cacheVolume(envelopeVol));
0738
0739   const double wedgeL = envelopeShape->getRMax() - envelopeShape->getRMin();
0740
0741   for (unsigned sector = 1; sector <= m_cfg.nSectors; ++sector) {
0742     envelopeWheel->add(
0743         makeTransform(GeoTrf::RotateZ3D(sector * m_sectorSize) *
0744                       GeoTrf::TranslateX3D(0.5 * (envelopeShape->getRMax() +
0745                                                   envelopeShape->getRMin())) *
0746                       GeoTrf::RotateY3D(90. * GeoModelKernelUnits::deg)));
0747
0748     envelopeWheel->add(
0749         assembleSmallWheelSector(wedgeL, copySign(1, wheelZ), sector));
0750   }
0751   envelope->add(makeTransform(GeoTrf::TranslateZ3D(wheelZ)));
0752   envelope->add(envelopeWheel);
0753 }
0754 }  // namespace ActsExamples