MuonSpectrometerMockupDetector/src/MockupSectorBuilder.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/MockupSectorBuilder.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Detector/DetectorVolume.hpp"
0013 #include "Acts/Detector/PortalGenerators.hpp"
0014 #include "Acts/Geometry/CuboidVolumeBounds.hpp"
0015 #include "Acts/Geometry/CylinderVolumeBounds.hpp"
0016 #include "Acts/Geometry/GeometryContext.hpp"
0017 #include "Acts/Geometry/VolumeBounds.hpp"
0018 #include "Acts/Navigation/DetectorVolumeFinders.hpp"
0019 #include "Acts/Navigation/InternalNavigation.hpp"
0020 #include "Acts/Plugins/Geant4/Geant4Converters.hpp"
0021 #include "Acts/Plugins/Geant4/Geant4DetectorElement.hpp"
0022 #include "Acts/Plugins/Geant4/Geant4DetectorSurfaceFactory.hpp"
0023 #include "Acts/Plugins/Geant4/Geant4PhysicalVolumeSelectors.hpp"
0024 #include "Acts/Surfaces/StrawSurface.hpp"
0025 #include "Acts/Surfaces/Surface.hpp"
0026 #include "Acts/Surfaces/SurfaceBounds.hpp"
0027 #include "Acts/Utilities/Logger.hpp"
0028 #include "Acts/Visualization/GeometryView3D.hpp"
0029 #include "Acts/Visualization/ObjVisualization3D.hpp"
0030 #include "Acts/Visualization/ViewConfig.hpp"
0031 #include "ActsExamples/Geant4Detector/GdmlDetectorConstruction.hpp"
0032 #include "ActsExamples/Geant4Detector/Geant4Detector.hpp"
0033
0034 #include <algorithm>
0035 #include <array>
0036 #include <cmath>
0037 #include <cstddef>
0038 #include <limits>
0039 #include <numbers>
0040 #include <stdexcept>
0041 #include <string>
0042 #include <utility>
0043 #include <vector>
0044
0045 namespace ActsExamples {
0046
0047 MockupSectorBuilder::MockupSectorBuilder(
0048     const MockupSectorBuilder::Config& config) {
0049   mCfg = config;
0050   GdmlDetectorConstruction geo_gdml(mCfg.gdmlPath, {});
0051   g4World = geo_gdml.Construct();
0052 }
0053
0054 std::shared_ptr<Acts::Experimental::DetectorVolume>
0055 MockupSectorBuilder::buildChamber(
0056     const MockupSectorBuilder::ChamberConfig& chamberConfig) {
0057   if (g4World == nullptr) {
0058     throw std::invalid_argument("MockupSector: No g4World initialized");
0059   }
0060
0061   const Acts::GeometryContext gctx;
0062
0063   // Geant4Detector Config creator with the g4world from the gdml file
0064   auto g4WorldConfig = Geant4Detector::Config();
0065   g4WorldConfig.name = "Chamber";
0066   g4WorldConfig.g4World = g4World;
0067
0068   // Get the sensitive and passive surfaces and pass to the g4World Config
0069   auto g4Sensitive =
0070       std::make_shared<Acts::Geant4PhysicalVolumeSelectors::NameSelector>(
0071           chamberConfig.SensitiveNames);
0072   auto g4Passive =
0073       std::make_shared<Acts::Geant4PhysicalVolumeSelectors::NameSelector>(
0074           chamberConfig.PassiveNames);
0075
0076   auto g4SurfaceOptions = Acts::Geant4DetectorSurfaceFactory::Options();
0077   g4SurfaceOptions.sensitiveSurfaceSelector = g4Sensitive;
0078   g4SurfaceOptions.passiveSurfaceSelector = g4Passive;
0079   g4WorldConfig.g4SurfaceOptions = g4SurfaceOptions;
0080
0081   auto g4Detector = Geant4Detector(g4WorldConfig);
0082   // Trigger the build of the detector
0083   auto [surface, elements] = Geant4Detector::buildGeant4Volumes(
0084       g4WorldConfig,
0085       *Acts::getDefaultLogger("MockupSectorBuilder", Acts::Logging::INFO));
0086
0087   // The vector that holds the converted sensitive surfaces of the chamber
0088   std::vector<std::shared_ptr<Acts::Surface>> strawSurfaces;
0089
0090   std::array<std::pair<float, float>, 3> min_max;
0091   std::fill(min_max.begin(), min_max.end(),
0092             std::make_pair<float, float>(std::numeric_limits<float>::max(),
0093                                          -std::numeric_limits<float>::max()));
0094
0095   // Convert the physical volumes of the detector elements to straw surfaces
0096   for (const auto& detectorElement : elements) {
0097     auto context = Acts::GeometryContext();
0098     auto g4conv = Acts::Geant4PhysicalVolumeConverter();
0099
0100     g4conv.forcedType = Acts::Surface::SurfaceType::Straw;
0101     auto g4ConvSurf = g4conv.Geant4PhysicalVolumeConverter::surface(
0102         detectorElement->g4PhysicalVolume(),
0103         detectorElement->transform(context));
0104
0105     strawSurfaces.push_back(g4ConvSurf);
0106
0107     min_max[0].first = std::min(
0108         min_max[0].first, static_cast<float>(g4ConvSurf->center(context).x()));
0109     min_max[0].second = std::max(
0110         min_max[0].second, static_cast<float>(g4ConvSurf->center(context).x()));
0111
0112     min_max[1].first = std::min(
0113         min_max[1].first, static_cast<float>(g4ConvSurf->center(context).y()));
0114     min_max[1].second = std::max(
0115         min_max[1].second, static_cast<float>(g4ConvSurf->center(context).y()));
0116
0117     min_max[2].first = std::min(
0118         min_max[2].first, static_cast<float>(g4ConvSurf->center(context).z()));
0119     min_max[2].second = std::max(
0120         min_max[2].second, static_cast<float>(g4ConvSurf->center(context).z()));
0121   }
0122
0123   // Create the bounds of the detector volumes
0124   float radius = strawSurfaces.front()->bounds().values()[0];
0125
0126   Acts::Vector3 minValues = {min_max[0].first, min_max[1].first,
0127                              min_max[2].first};
0128   Acts::Vector3 maxValues = {min_max[0].second, min_max[1].second,
0129                              min_max[2].second};
0130
0131   double hx =
0132       strawSurfaces.front()->bounds().values()[1] + mCfg.toleranceOverlap;
0133   double hy = 0.5 * ((maxValues.y() + radius) - (minValues.y() - radius)) +
0134               mCfg.toleranceOverlap;
0135   double hz = 0.5 * ((maxValues.z() + radius) - (minValues.z() - radius)) +
0136               mCfg.toleranceOverlap;
0137
0138   auto detectorVolumeBounds =
0139       std::make_shared<Acts::CuboidVolumeBounds>(hx, hy, hz);
0140
0141   Acts::Vector3 chamber_position = {(maxValues.x() + minValues.x()) / 2,
0142                                     (maxValues.y() + minValues.y()) / 2,
0143                                     (maxValues.z() + minValues.z()) / 2};
0144
0145   // create the detector volume for the chamber
0146   auto detectorVolume = Acts::Experimental::DetectorVolumeFactory::construct(
0147       Acts::Experimental::defaultPortalAndSubPortalGenerator(), gctx,
0148       chamberConfig.name,
0149       Acts::Transform3(Acts::Translation3(chamber_position)),
0150       std::move(detectorVolumeBounds), strawSurfaces,
0151       std::vector<std::shared_ptr<Acts::Experimental::DetectorVolume>>{},
0152       Acts::Experimental::tryAllSubVolumes(),
0153       Acts::Experimental::tryAllPortalsAndSurfaces());
0154
0155   return detectorVolume;
0156 }
0157
0158 std::shared_ptr<Acts::Experimental::DetectorVolume>
0159 MockupSectorBuilder::buildSector(
0160     std::vector<std::shared_ptr<Acts::Experimental::DetectorVolume>>
0161         detVolumes) {
0162   if (mCfg.NumberOfSectors > maxNumberOfSectors) {
0163     throw std::invalid_argument("MockupSector:Number of max sectors exceeded");
0164   }
0165
0166   const Acts::GeometryContext gctx;
0167
0168   // sort the detector volumes by their radial distance (from
0169   // innermost---->outermost)
0170   std::ranges::sort(detVolumes, {},
0171                     [](const auto& detVol) { return detVol->center().y(); });
0172
0173   auto xA = detVolumes.back()->center().x() +
0174             detVolumes.back()->volumeBounds().values()[0];
0175   auto yA = detVolumes.back()->center().y() -
0176             detVolumes.back()->volumeBounds().values()[1];
0177   auto zA = detVolumes.back()->center().z();
0178
0179   auto xB = detVolumes.back()->center().x() -
0180             detVolumes.back()->volumeBounds().values()[0];
0181   auto yB = detVolumes.back()->center().y() -
0182             detVolumes.back()->volumeBounds().values()[1];
0183   auto zB = detVolumes.back()->center().z();
0184
0185   Acts::Vector3 pointA = {xA, yA, zA};
0186   Acts::Vector3 pointB = {xB, yB, zB};
0187
0188   // calculate the phi angles of the vectors
0189   auto phiA = Acts::VectorHelpers::phi(pointA);
0190   auto phiB = Acts::VectorHelpers::phi(pointB);
0191   double sectorAngle = std::numbers::pi;
0192
0193   double halfPhi = std::numbers::pi / mCfg.NumberOfSectors;
0194
0195   if (mCfg.NumberOfSectors == 1) {
0196     halfPhi = (phiB - phiA) / 2;
0197     sectorAngle = halfPhi;
0198   }
0199
0200   const int detVolumesSize = detVolumes.size();
0201
0202   std::vector<float> rmins(detVolumesSize);
0203   std::vector<float> rmaxs(detVolumesSize);
0204   std::vector<float> halfZ(detVolumesSize);
0205   std::vector<std::shared_ptr<Acts::CylinderVolumeBounds>>
0206       cylinderVolumesBounds(detVolumesSize);
0207
0208   for (int i = 0; i < detVolumesSize; i++) {
0209     const auto& detVol = detVolumes[i];
0210     rmins[i] = detVol->center().y() - detVol->volumeBounds().values()[1] -
0211                mCfg.toleranceOverlap;
0212     rmaxs[i] = std::sqrt(std::pow(detVol->volumeBounds().values()[0], 2) +
0213                          std::pow(detVol->center().y() +
0214                                       detVol->volumeBounds().values()[1],
0215                                   2)) +
0216                mCfg.toleranceOverlap;
0217     halfZ[i] = detVol->volumeBounds().values()[2];
0218
0219     cylinderVolumesBounds[i] = std::make_shared<Acts::CylinderVolumeBounds>(
0220         rmins[i], rmaxs[i], halfZ[i], sectorAngle);
0221   }
0222
0223   const Acts::Vector3 pos = {0., 0., 0.};
0224
0225   // the transform of the cylinder volume
0226   Acts::AngleAxis3 rotZ(std::numbers::pi / 2., Acts::Vector3(0., 0., 1));
0227   auto transform = Acts::Transform3(Acts::Translation3(pos));
0228   transform *= rotZ;
0229
0230   // create a vector for the shifted surfaces of each chamber
0231   std::vector<std::shared_ptr<Acts::Surface>> shiftedSurfaces = {};
0232
0233   // creare an array of vectors that holds all the chambers of each sector
0234   std::vector<std::vector<std::shared_ptr<Acts::Experimental::DetectorVolume>>>
0235       chambersOfSectors(detVolumesSize);
0236
0237   std::vector<std::shared_ptr<Acts::Experimental::DetectorVolume>>
0238       detectorCylinderVolumesOfSector = {};
0239
0240   for (int i = 0; i < mCfg.NumberOfSectors; i++) {
0241     Acts::AngleAxis3 rotation(2 * i * halfPhi, Acts::Vector3(0., 0., 1.));
0242
0243     for (int itr = 0; itr < detVolumesSize; itr++) {
0244       const auto& detVol = detVolumes[itr];
0245
0246       auto shift_vol =
0247           rotation * Acts::Transform3(Acts::Translation3(detVol->center()));
0248
0249       for (auto& detSurf : detVol->surfaces()) {
0250         auto shift_surf = Acts::Transform3::Identity() * rotation;
0251
0252         // create the shifted surfaces by creating copied surface objects
0253         auto strawSurfaceObject = Acts::Surface::makeShared<Acts::StrawSurface>(
0254             detSurf->transform(Acts::GeometryContext()),
0255             detSurf->bounds().values()[0], detSurf->bounds().values()[1]);
0256
0257         auto copiedTransformStrawSurface =
0258             Acts::Surface::makeShared<Acts::StrawSurface>(
0259                 Acts::GeometryContext(), *strawSurfaceObject, shift_surf);
0260
0261         shiftedSurfaces.push_back(copiedTransformStrawSurface);
0262       }
0263
0264       // create the bounds of the volumes of each chamber
0265       auto bounds = std::make_unique<Acts::CuboidVolumeBounds>(
0266           detVol->volumeBounds().values()[0],
0267           detVol->volumeBounds().values()[1],
0268           detVol->volumeBounds().values()[2]);
0269       // create the shifted chamber
0270       auto detectorVolumeSec =
0271           Acts::Experimental::DetectorVolumeFactory::construct(
0272               Acts::Experimental::defaultPortalAndSubPortalGenerator(), gctx,
0273               "detectorVolumeChamber_" + std::to_string(itr), shift_vol,
0274               std::move(bounds), shiftedSurfaces,
0275               std::vector<
0276                   std::shared_ptr<Acts::Experimental::DetectorVolume>>{},
0277               Acts::Experimental::tryAllSubVolumes(),
0278               Acts::Experimental::tryAllPortalsAndSurfaces());
0279
0280       chambersOfSectors[itr].push_back(detectorVolumeSec);
0281
0282       shiftedSurfaces.clear();
0283
0284     }  // end of detector volumes
0285
0286   }  // end of number of sectors
0287
0288   for (std::size_t i = 0; i < cylinderVolumesBounds.size(); ++i) {
0289     detectorCylinderVolumesOfSector.push_back(
0290         Acts::Experimental::DetectorVolumeFactory::construct(
0291             Acts::Experimental::defaultPortalAndSubPortalGenerator(), gctx,
0292             "cylinder_volume_" + std::to_string(i), transform,
0293             std::move(cylinderVolumesBounds[i]),
0294             std::vector<std::shared_ptr<Acts::Surface>>{}, chambersOfSectors[i],
0295             Acts::Experimental::tryAllSubVolumes(),
0296             Acts::Experimental::tryAllPortalsAndSurfaces()));
0297
0298   }  // end of cylinder volumes
0299
0300   auto cylinderVolumesBoundsOfMother =
0301       std::make_shared<Acts::CylinderVolumeBounds>(
0302           rmins.front(), rmaxs.back(),
0303           *std::max_element(halfZ.begin(), halfZ.end()), sectorAngle);
0304
0305   // creation of the mother volume
0306   auto detectorVolume = Acts::Experimental::DetectorVolumeFactory::construct(
0307       Acts::Experimental::defaultPortalAndSubPortalGenerator(), gctx,
0308       "detectorVolumeSector", transform,
0309       std::move(cylinderVolumesBoundsOfMother),
0310       std::vector<std::shared_ptr<Acts::Surface>>{},
0311       detectorCylinderVolumesOfSector, Acts::Experimental::tryAllSubVolumes(),
0312       Acts::Experimental::tryAllPortalsAndSurfaces());
0313
0314   return detectorVolume;
0315 }
0316
0317 void MockupSectorBuilder::drawSector(
0318     const std::shared_ptr<Acts::Experimental::DetectorVolume>&
0319         detectorVolumeSector,
0320     const std::string& nameObjFile) {
0321   Acts::ViewConfig sConfig = Acts::s_viewSensitive;
0322
0323   Acts::ObjVisualization3D objSector;
0324
0325   Acts::GeometryView3D::drawDetectorVolume(
0326       objSector, *detectorVolumeSector, Acts::GeometryContext(),
0327       Acts::Transform3::Identity(), sConfig);
0328
0329   objSector.write(nameObjFile);
0330 }
0331
0332 }  // namespace ActsExamples