Root/src/TGeoLayerBuilder.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 "ActsPlugins/Root/TGeoLayerBuilder.hpp"
0010
0011 #include "Acts/Geometry/Extent.hpp"
0012 #include "Acts/Geometry/LayerCreator.hpp"
0013 #include "Acts/Geometry/ProtoLayer.hpp"
0014 #include "Acts/Geometry/ProtoLayerHelper.hpp"
0015 #include "Acts/Utilities/Helpers.hpp"
0016 #include "ActsPlugins/Root/ITGeoDetectorElementSplitter.hpp"
0017 #include "ActsPlugins/Root/ITGeoIdentifierProvider.hpp"
0018 #include "ActsPlugins/Root/TGeoDetectorElement.hpp"
0019 #include "ActsPlugins/Root/TGeoParser.hpp"
0020 #include "ActsPlugins/Root/TGeoPrimitivesHelper.hpp"
0021
0022 #include <ostream>
0023 #include <stdexcept>
0024
0025 #include "TGeoManager.h"
0026 #include "TGeoMatrix.h"
0027
0028 namespace Acts {
0029 class ISurfaceMaterial;
0030 }  // namespace Acts
0031
0032 using namespace Acts;
0033
0034 ActsPlugins::TGeoLayerBuilder::TGeoLayerBuilder(
0035     const ActsPlugins::TGeoLayerBuilder::Config& config,
0036     std::unique_ptr<const Logger> logger)
0037     : m_logger(std::move(logger)) {
0038   setConfiguration(config);
0039 }
0040
0041 ActsPlugins::TGeoLayerBuilder::~TGeoLayerBuilder() = default;
0042
0043 void ActsPlugins::TGeoLayerBuilder::setConfiguration(
0044     const ActsPlugins::TGeoLayerBuilder::Config& config) {
0045   m_cfg = config;
0046 }
0047
0048 void ActsPlugins::TGeoLayerBuilder::setLogger(
0049     std::unique_ptr<const Logger> newLogger) {
0050   m_logger = std::move(newLogger);
0051 }
0052
0053 const LayerVector ActsPlugins::TGeoLayerBuilder::negativeLayers(
0054     const GeometryContext& gctx) const {
0055   // @todo Remove this hack once the m_elementStore mess is sorted out
0056   auto mutableThis = const_cast<TGeoLayerBuilder*>(this);
0057   LayerVector nVector;
0058   mutableThis->buildLayers(gctx, nVector, -1);
0059   return nVector;
0060 }
0061
0062 const LayerVector ActsPlugins::TGeoLayerBuilder::centralLayers(
0063     const GeometryContext& gctx) const {
0064   // @todo Remove this hack once the m_elementStore mess is sorted out
0065   auto mutableThis = const_cast<TGeoLayerBuilder*>(this);
0066   LayerVector cVector;
0067   mutableThis->buildLayers(gctx, cVector, 0);
0068   return cVector;
0069 }
0070
0071 const LayerVector ActsPlugins::TGeoLayerBuilder::positiveLayers(
0072     const GeometryContext& gctx) const {
0073   // @todo Remove this hack once the m_elementStore mess is sorted out
0074   auto mutableThis = const_cast<TGeoLayerBuilder*>(this);
0075   LayerVector pVector;
0076   mutableThis->buildLayers(gctx, pVector, 1);
0077   return pVector;
0078 }
0079
0080 void ActsPlugins::TGeoLayerBuilder::buildLayers(const GeometryContext& gctx,
0081                                                 LayerVector& layers, int type) {
0082   // Bail out if you have no gGeoManager
0083   if (gGeoManager == nullptr) {
0084     ACTS_WARNING("No gGeoManager found - bailing out.");
0085     return;
0086   }
0087
0088   using LayerSurfaceVector = std::vector<std::shared_ptr<const Surface>>;
0089   LayerSurfaceVector layerSurfaces;
0090
0091   std::vector<LayerConfig> layerConfigs = m_cfg.layerConfigurations[type + 1];
0092   std::string layerType = m_layerTypes[type + 1];
0093
0094   // Appropriate screen output
0095   std::string addonOutput = m_cfg.layerSplitToleranceR[type + 1] > 0.
0096                                 ? std::string(", splitting in r")
0097                                 : std::string("");
0098   addonOutput += m_cfg.layerSplitToleranceZ[type + 1] > 0.
0099                      ? std::string(", splitting in z")
0100                      : std::string("");
0101   addonOutput += std::string(".");
0102
0103   // Screen output of the configuration
0104   ACTS_DEBUG(layerType << " layers : found " << layerConfigs.size()
0105                        << " configuration(s)" + addonOutput);
0106
0107   // Helper function to fill the layer
0108   auto fillLayer = [&](const LayerSurfaceVector& lSurfaces,
0109                        const LayerConfig& lCfg, unsigned int pl_id = 0) {
0110     int nb0 = 0;
0111     int nt0 = 0;
0112     bool is_autobinning = ((lCfg.binning0.size() == 1) &&
0113                            (std::get<int>(lCfg.binning0.at(0)) <= 0));
0114     if (!is_autobinning && std::get<int>(lCfg.binning0.at(pl_id)) <= 0) {
0115       throw std::invalid_argument(
0116           "Incorrect binning configuration found for loc0 protolayer #" +
0117           std::to_string(pl_id) +
0118           ". Layer is autobinned: No mixed binning (manual and auto) for loc0 "
0119           "possible between layers in a single subvolume. Quitting");
0120     }
0121     if (is_autobinning) {
0122       // Set binning by hand if nb0 > 0 and nb1 > 0
0123       nb0 = std::get<int>(lCfg.binning0.at(0));
0124       // Read the binning type
0125       nt0 = std::get<BinningType>(lCfg.binning0.at(0));
0126     } else if (pl_id < lCfg.binning0.size()) {
0127       // Set binning by hand if nb0 > 0 and nb1 > 0
0128       nb0 = std::get<int>(lCfg.binning0.at(pl_id));
0129     }
0130
0131     int nb1 = 0;
0132     int nt1 = 0;
0133     is_autobinning = (lCfg.binning1.size() == 1) &&
0134                      (std::get<int>(lCfg.binning1.at(0)) <= 0);
0135     if (!is_autobinning && std::get<int>(lCfg.binning1.at(pl_id)) <= 0) {
0136       throw std::invalid_argument(
0137           "Incorrect binning configuration found for loc1 protolayer #" +
0138           std::to_string(pl_id) +
0139           ". Layer is autobinned: No mixed binning (manual and auto) for loc1 "
0140           "possible between layers in a single subvolume. Quitting");
0141     }
0142     if (is_autobinning) {
0143       // Set binning by hand if nb0 > 0 and nb1 > 0
0144       nb1 = std::get<int>(lCfg.binning1.at(0));
0145       // For a binning type
0146       nt1 = std::get<BinningType>(lCfg.binning1.at(0));
0147     } else if (pl_id < lCfg.binning1.size()) {
0148       // Set binning by hand if nb0 > 0 and nb1 > 0
0149       nb1 = std::get<int>(lCfg.binning1.at(pl_id));
0150     }
0151
0152     if (type == 0) {
0153       ProtoLayer pl(gctx, lSurfaces);
0154       ACTS_DEBUG("- creating CylinderLayer with "
0155                  << lSurfaces.size()
0156                  << " surfaces at r = " << pl.medium(AxisDirection::AxisR));
0157
0158       pl.envelope[AxisDirection::AxisR] = {lCfg.envelope.first,
0159                                            lCfg.envelope.second};
0160       pl.envelope[AxisDirection::AxisZ] = {lCfg.envelope.second,
0161                                            lCfg.envelope.second};
0162       if (nb0 >= 0 && nb1 >= 0) {
0163         layers.push_back(
0164             m_cfg.layerCreator->cylinderLayer(gctx, lSurfaces, nb0, nb1, pl));
0165       } else {
0166         layers.push_back(
0167             m_cfg.layerCreator->cylinderLayer(gctx, lSurfaces, nt0, nt1, pl));
0168       }
0169     } else {
0170       ProtoLayer pl(gctx, lSurfaces);
0171       ACTS_DEBUG("- creating DiscLayer with "
0172                  << lSurfaces.size()
0173                  << " surfaces at z = " << pl.medium(AxisDirection::AxisZ));
0174
0175       pl.envelope[AxisDirection::AxisR] = {lCfg.envelope.first,
0176                                            lCfg.envelope.second};
0177       pl.envelope[AxisDirection::AxisZ] = {lCfg.envelope.second,
0178                                            lCfg.envelope.second};
0179       if (nb0 >= 0 && nb1 >= 0) {
0180         layers.push_back(
0181             m_cfg.layerCreator->discLayer(gctx, lSurfaces, nb0, nb1, pl));
0182       } else {
0183         layers.push_back(
0184             m_cfg.layerCreator->discLayer(gctx, lSurfaces, nt0, nt1, pl));
0185       }
0186     }
0187   };
0188
0189   for (auto layerCfg : layerConfigs) {
0190     ACTS_DEBUG("- layer configuration found for layer " << layerCfg.volumeName
0191                                                         << " with sensors ");
0192     for (auto& sensor : layerCfg.sensorNames) {
0193       ACTS_DEBUG("  - sensor: " << sensor);
0194     }
0195     if (!layerCfg.parseRanges.empty()) {
0196       for (const auto& [axisDir, pRange] : layerCfg.parseRanges) {
0197         ACTS_DEBUG("- layer parsing restricted in "
0198                    << axisDirectionName(axisDir) << " to [" << pRange.first
0199                    << "/" << pRange.second << "].");
0200       }
0201     }
0202     if (!layerCfg.splitConfigs.empty()) {
0203       for (const auto& [axisDir, tConfig] : layerCfg.splitConfigs) {
0204         ACTS_DEBUG("- layer splitting attempt in " << axisDirectionName(axisDir)
0205                                                    << " with tolerance "
0206                                                    << tConfig << ".");
0207       }
0208     }
0209
0210     // Either pick the configured volume or take the top level volume
0211     // and retrieve its global transformation.
0212     TGeoVolume* tVolume =
0213         gGeoManager->FindVolumeFast(layerCfg.volumeName.c_str());
0214     TGeoHMatrix gmatrix = TGeoIdentity((layerCfg.volumeName + "ID").c_str());
0215
0216     if (tVolume == nullptr) {
0217       tVolume = gGeoManager->GetTopVolume();
0218       ACTS_DEBUG("- search volume is TGeo top volume");
0219     } else {
0220       ACTS_DEBUG("- setting search volume to " << tVolume->GetName());
0221
0222       auto node = TGeoParser::findNodeRecursive(gGeoManager->GetTopNode(),
0223                                                 tVolume->GetName());
0224
0225       if (node == nullptr) {
0226         std::string volname(tVolume->GetName());
0227         throw std::invalid_argument("Could not locate node for " + volname);
0228       }
0229
0230       gmatrix = *(node->GetMatrix());
0231     }
0232
0233     if (tVolume != nullptr) {
0234       TGeoParser::Options tgpOptions;
0235       tgpOptions.volumeNames = {layerCfg.volumeName};
0236       tgpOptions.targetNames = layerCfg.sensorNames;
0237       tgpOptions.parseRanges = layerCfg.parseRanges;
0238       tgpOptions.unit = m_cfg.unit;
0239       TGeoParser::State tgpState;
0240       tgpState.volume = tVolume;
0241
0242       ACTS_DEBUG("- applying  " << layerCfg.parseRanges.size()
0243                                 << " search restrictions.");
0244       for (const auto& [axisDir, pRange] : layerCfg.parseRanges) {
0245         ACTS_VERBOSE(" - range " << axisDirectionName(axisDir) << " within [ "
0246                                  << pRange.first << ", " << pRange.second
0247                                  << "]");
0248       }
0249
0250       TGeoParser::select(tgpState, tgpOptions, gmatrix);
0251
0252       ACTS_DEBUG("- number of selected nodes found : "
0253                  << tgpState.selectedNodes.size());
0254
0255       for (auto& snode : tgpState.selectedNodes) {
0256         auto identifier =
0257             m_cfg.identifierProvider != nullptr
0258                 ? m_cfg.identifierProvider->identify(gctx, *snode.node)
0259                 : TGeoDetectorElement::Identifier();
0260
0261         auto tgElement = m_cfg.detectorElementFactory(
0262             identifier, *snode.node, *snode.transform, layerCfg.localAxes,
0263             m_cfg.unit, nullptr);
0264
0265         std::vector<std::shared_ptr<const TGeoDetectorElement>> tgElements =
0266             (m_cfg.detectorElementSplitter == nullptr)
0267                 ? std::vector<
0268                       std::shared_ptr<const TGeoDetectorElement>>{tgElement}
0269                 : m_cfg.detectorElementSplitter->split(gctx, tgElement);
0270
0271         for (const auto& tge : tgElements) {
0272           m_elementStore.push_back(tge);
0273           layerSurfaces.push_back(tge->surface().getSharedPtr());
0274         }
0275       }
0276
0277       ACTS_DEBUG("- created TGeoDetectorElements : " << layerSurfaces.size());
0278
0279       if (m_cfg.protoLayerHelper != nullptr && !layerCfg.splitConfigs.empty()) {
0280         auto protoLayers = m_cfg.protoLayerHelper->protoLayers(
0281             gctx, unpackSmartPointers(layerSurfaces), layerCfg.splitConfigs);
0282         ACTS_DEBUG("- splitting into " << protoLayers.size() << " layers.");
0283
0284         // Number of options mismatch and has not been configured for
0285         // auto-binning
0286         const bool is_loc0_n_config =
0287             layerCfg.binning0.size() == protoLayers.size();
0288         const bool is_loc0_autobinning =
0289             (layerCfg.binning0.size() == 1) &&
0290             (std::get<int>(layerCfg.binning0.at(0)) <= 0);
0291         const bool is_loc1_n_config =
0292             layerCfg.binning1.size() == protoLayers.size();
0293         const bool is_loc1_autobinning =
0294             (layerCfg.binning1.size() == 1) &&
0295             (std::get<int>(layerCfg.binning1.at(0)) <= 0);
0296         if ((!is_loc0_n_config && !is_loc0_autobinning) ||
0297             (!is_loc1_n_config && !is_loc1_autobinning)) {
0298           throw std::invalid_argument(
0299               "Incorrect binning configuration found: Number of configurations "
0300               "does not match number of protolayers in subvolume " +
0301               layerCfg.volumeName + ". Quitting.");
0302         }
0303         unsigned int layer_id = 0;
0304         for (auto& pLayer : protoLayers) {
0305           layerSurfaces.clear();
0306
0307           for (const auto& lsurface : pLayer.surfaces()) {
0308             layerSurfaces.push_back(lsurface->getSharedPtr());
0309           }
0310           fillLayer(layerSurfaces, layerCfg, layer_id);
0311           layer_id++;
0312         }
0313       } else {
0314         fillLayer(layerSurfaces, layerCfg);
0315       }
0316     }
0317   }
0318   return;
0319 }
0320
0321 std::shared_ptr<ActsPlugins::TGeoDetectorElement>
0322 ActsPlugins::TGeoLayerBuilder::defaultElementFactory(
0323     const TGeoDetectorElement::Identifier& identifier, const TGeoNode& tGeoNode,
0324     const TGeoMatrix& tGeoMatrix, ActsPlugins::TGeoAxes axes, double scalor,
0325     std::shared_ptr<const ISurfaceMaterial> material) {
0326   return std::make_shared<TGeoDetectorElement>(
0327       identifier, tGeoNode, tGeoMatrix, axes, scalor, std::move(material));
0328 }