MuonSpectrometerMockupDetector/src/MockupSectorBuilder.cpp

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