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