TGeo/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/TGeo/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/TGeo/ITGeoDetectorElementSplitter.hpp"
0016 #include "Acts/Plugins/TGeo/ITGeoIdentifierProvider.hpp"
0017 #include "Acts/Plugins/TGeo/TGeoDetectorElement.hpp"
0018 #include "Acts/Plugins/TGeo/TGeoParser.hpp"
0019 #include "Acts/Plugins/TGeo/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_cfg(), 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,
0108                        unsigned int pl_id = 0) -> void {
0109     int nb0 = 0, 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, nt1 = 0;
0130     is_autobinning = (lCfg.binning1.size() == 1) &&
0131                      (std::get<int>(lCfg.binning1.at(0)) <= 0);
0132     if (!is_autobinning && std::get<int>(lCfg.binning1.at(pl_id)) <= 0) {
0133       throw std::invalid_argument(
0134           "Incorrect binning configuration found for loc1 protolayer #" +
0135           std::to_string(pl_id) +
0136           ". Layer is autobinned: No mixed binning (manual and auto) for loc1 "
0137           "possible between layers in a single subvolume. Quitting");
0138     }
0139     if (is_autobinning) {
0140       // Set binning by hand if nb0 > 0 and nb1 > 0
0141       nb1 = std::get<int>(lCfg.binning1.at(0));
0142       // For a binning type
0143       nt1 = std::get<BinningType>(lCfg.binning1.at(0));
0144     } else if (pl_id < lCfg.binning1.size()) {
0145       // Set binning by hand if nb0 > 0 and nb1 > 0
0146       nb1 = std::get<int>(lCfg.binning1.at(pl_id));
0147     }
0148
0149     if (type == 0) {
0150       ProtoLayer pl(gctx, lSurfaces);
0151       ACTS_DEBUG("- creating CylinderLayer with "
0152                  << lSurfaces.size()
0153                  << " surfaces at r = " << pl.medium(AxisDirection::AxisR));
0154
0155       pl.envelope[Acts::AxisDirection::AxisR] = {lCfg.envelope.first,
0156                                                  lCfg.envelope.second};
0157       pl.envelope[Acts::AxisDirection::AxisZ] = {lCfg.envelope.second,
0158                                                  lCfg.envelope.second};
0159       if (nb0 >= 0 && nb1 >= 0) {
0160         layers.push_back(
0161             m_cfg.layerCreator->cylinderLayer(gctx, lSurfaces, nb0, nb1, pl));
0162       } else {
0163         layers.push_back(
0164             m_cfg.layerCreator->cylinderLayer(gctx, lSurfaces, nt0, nt1, pl));
0165       }
0166     } else {
0167       ProtoLayer pl(gctx, lSurfaces);
0168       ACTS_DEBUG("- creating DiscLayer with "
0169                  << lSurfaces.size()
0170                  << " surfaces at z = " << pl.medium(AxisDirection::AxisZ));
0171
0172       pl.envelope[Acts::AxisDirection::AxisR] = {lCfg.envelope.first,
0173                                                  lCfg.envelope.second};
0174       pl.envelope[Acts::AxisDirection::AxisZ] = {lCfg.envelope.second,
0175                                                  lCfg.envelope.second};
0176       if (nb0 >= 0 && nb1 >= 0) {
0177         layers.push_back(
0178             m_cfg.layerCreator->discLayer(gctx, lSurfaces, nb0, nb1, pl));
0179       } else {
0180         layers.push_back(
0181             m_cfg.layerCreator->discLayer(gctx, lSurfaces, nt0, nt1, pl));
0182       }
0183     }
0184   };
0185
0186   for (auto layerCfg : layerConfigs) {
0187     ACTS_DEBUG("- layer configuration found for layer " << layerCfg.volumeName
0188                                                         << " with sensors ");
0189     for (auto& sensor : layerCfg.sensorNames) {
0190       ACTS_DEBUG("  - sensor: " << sensor);
0191     }
0192     if (!layerCfg.parseRanges.empty()) {
0193       for (const auto& pRange : layerCfg.parseRanges) {
0194         ACTS_DEBUG("- layer parsing restricted in "
0195                    << axisDirectionName(pRange.first) << " to ["
0196                    << pRange.second.first << "/" << pRange.second.second
0197                    << "].");
0198       }
0199     }
0200     if (!layerCfg.splitConfigs.empty()) {
0201       for (const auto& sConfig : layerCfg.splitConfigs) {
0202         ACTS_DEBUG("- layer splitting attempt in "
0203                    << axisDirectionName(sConfig.first) << " with tolerance "
0204                    << sConfig.second << ".");
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& prange : layerCfg.parseRanges) {
0243         ACTS_VERBOSE(" - range " << axisDirectionName(prange.first)
0244                                  << " within [ " << prange.second.first << ", "
0245                                  << prange.second.second << "]");
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 =
0260             m_cfg.elementFactory(identifier, *snode.node, *snode.transform,
0261                                  layerCfg.localAxes, 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 }