GenericDetector/src/ProtoLayerCreator.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/GenericDetector/ProtoLayerCreator.hpp"
0010
0011 #include "Acts/Utilities/TransformHelpers.hpp"
0012
0013 using Acts::VectorHelpers::phi;
0014
0015 namespace ActsExamples::Generic {
0016
0017 std::vector<ProtoLayerSurfaces> ProtoLayerCreator::centralProtoLayers(
0018     const Acts::GeometryContext& gctx) const {
0019   // create the vector
0020   std::vector<ProtoLayerSurfaces> cpLayers;
0021
0022   // ----------------------- central layers -------------------------
0023   // the central layers
0024   std::size_t numcLayers = m_cfg.centralLayerRadii.size();
0025   if (numcLayers != 0u) {
0026     ACTS_DEBUG("Configured to build " << numcLayers
0027                                       << " active central layers.");
0028     cpLayers.reserve(numcLayers);
0029     // loop through
0030     for (std::size_t icl = 0; icl < numcLayers; ++icl) {
0031       // layer R/Z
0032       double layerR = m_cfg.centralLayerRadii.at(icl);
0033       // some screen output
0034       ACTS_DEBUG("Build layer " << icl << " with target radius = " << layerR);
0035
0036       // prepare the Surface vector
0037       std::vector<std::shared_ptr<Acts::Surface>> sVector;
0038       // assign the current envelope
0039       double layerEnvelopeCoverZ =
0040           !m_cfg.centralLayerEnvelopes.empty()
0041               ? m_cfg.centralLayerEnvelopes.at(icl).second
0042               : 0.;
0043       // module size & tilt
0044       double modulePhiTilt = m_cfg.centralModuleTiltPhi.at(icl);
0045       double moduleHalfX = m_cfg.centralModuleHalfX.at(icl);
0046       double moduleHalfY = m_cfg.centralModuleHalfY.at(icl);
0047       double moduleThickness = m_cfg.centralModuleThickness.at(icl);
0048       // create the shared module
0049       auto moduleBounds =
0050           std::make_shared<Acts::RectangleBounds>(moduleHalfX, moduleHalfY);
0051       std::size_t nCentralModules =
0052           m_cfg.centralModuleBinningSchema.at(icl).first *
0053           m_cfg.centralModuleBinningSchema.at(icl).second;
0054
0055       ACTS_DEBUG("- number of modules "
0056                  << nCentralModules << " ( from "
0057                  << m_cfg.centralModuleBinningSchema.at(icl).first << " x "
0058                  << m_cfg.centralModuleBinningSchema.at(icl).second << " )");
0059
0060       sVector.reserve(nCentralModules);
0061
0062       // prepartation :
0063       // create the Module material from input
0064       std::shared_ptr<const Acts::ISurfaceMaterial> moduleMaterialPtr = nullptr;
0065       if (!m_cfg.centralModuleMaterial.empty()) {
0066         // get the sensor material from configuration
0067         moduleMaterialPtr = m_cfg.centralModuleMaterial.at(icl);
0068       }
0069
0070       // confirm
0071       if (m_cfg.centralModulePositions.at(icl).size() != nCentralModules) {
0072         ACTS_WARNING("Mismatching module numbers, configuration error!");
0073         ACTS_WARNING("- Binning schema suggests : " << nCentralModules);
0074         ACTS_WARNING("- Positions provided are  : "
0075                      << m_cfg.centralModulePositions.at(icl).size());
0076       }
0077       // loop over the position, create the modules
0078       for (auto& moduleCenter : m_cfg.centralModulePositions.at(icl)) {
0079         // create the association transform
0080         double modulePhi = phi(moduleCenter);
0081         // the local z axis is the normal vector
0082         Acts::Vector3 moduleLocalZ(std::cos(modulePhi + modulePhiTilt),
0083                                    std::sin(modulePhi + modulePhiTilt), 0.);
0084         // the local y axis is the global z axis
0085         Acts::Vector3 moduleLocalY(0., 0., 1);
0086         // the local x axis the normal to local y,z
0087         Acts::Vector3 moduleLocalX(-std::sin(modulePhi + modulePhiTilt),
0088                                    std::cos(modulePhi + modulePhiTilt), 0.);
0089         // create the RotationMatrix
0090         Acts::RotationMatrix3 moduleRotation;
0091         moduleRotation.col(0) = moduleLocalX;
0092         moduleRotation.col(1) = moduleLocalY;
0093         moduleRotation.col(2) = moduleLocalZ;
0094         // get the moduleTransform
0095         Acts::Transform3 moduleTransform =
0096             Acts::makeTransform3(moduleRotation, moduleCenter);
0097         // stereo angle if necessary
0098         if (!m_cfg.centralModuleFrontsideStereo.empty() &&
0099             m_cfg.centralModuleFrontsideStereo.at(icl) != 0.) {
0100           // twist by the stereo angle
0101           double stereo = m_cfg.centralModuleFrontsideStereo.at(icl);
0102           moduleTransform *= Acts::AngleAxis3(-stereo, Acts::Vector3::UnitZ());
0103         }
0104         // create the module
0105         auto moduleElement = m_cfg.detectorElementFactory(
0106             moduleTransform, moduleBounds, moduleThickness, moduleMaterialPtr);
0107         // register the surface
0108         sVector.push_back(moduleElement->surface().getSharedPtr());
0109         // IF double modules exist
0110         // and the backside one (if configured to do so)
0111         if (!m_cfg.centralModuleBacksideGap.empty()) {
0112           // create the module identifier
0113
0114           Acts::Vector3 bsModuleCenter =
0115               moduleCenter +
0116               m_cfg.centralModuleBacksideGap.at(icl) * moduleLocalZ;
0117           Acts::Transform3 bsModuleTransform =
0118               Acts::makeTransform3(moduleRotation, bsModuleCenter);
0119           // apply the stereo
0120           if (!m_cfg.centralModuleBacksideStereo.empty()) {
0121             // twist by the stereo angle
0122             double stereoBackSide = m_cfg.centralModuleBacksideStereo.at(icl);
0123             bsModuleTransform *=
0124                 Acts::AngleAxis3(-stereoBackSide, Acts::Vector3::UnitZ());
0125           }
0126           // create the backseide moulde
0127           auto bsModuleElement =
0128               m_cfg.detectorElementFactory(bsModuleTransform, moduleBounds,
0129                                            moduleThickness, moduleMaterialPtr);
0130           // register the backside surface
0131           sVector.push_back(bsModuleElement->surface().getSharedPtr());
0132         }
0133       }
0134
0135       std::size_t phiBins = m_cfg.centralModuleBinningSchema.at(icl).first;
0136       phiBins *= m_cfg.centralLayerBinMultipliers.first;
0137       std::size_t zBins = m_cfg.centralModuleBinningSchema.at(icl).second;
0138       zBins *= m_cfg.centralLayerBinMultipliers.second;
0139       // create the surface array - it will also fill the accessible binmember
0140       // cache if available
0141       Acts::MutableProtoLayer pl(gctx, sVector);
0142       pl.envelope[Acts::AxisDirection::AxisR] = {m_cfg.approachSurfaceEnvelope,
0143                                                  m_cfg.approachSurfaceEnvelope};
0144       pl.envelope[Acts::AxisDirection::AxisZ] = {layerEnvelopeCoverZ,
0145                                                  layerEnvelopeCoverZ};
0146
0147       // Record the proto layer and the surfaces for the later layer building
0148       ProtoLayerSurfaces pls{std::move(pl), sVector, phiBins, zBins};
0149       cpLayers.push_back(std::move(pls));
0150     }
0151   }
0152   return cpLayers;
0153 }
0154
0155 std::vector<ProtoLayerSurfaces> ProtoLayerCreator::negativeProtoLayers(
0156     const Acts::GeometryContext& gctx) const {
0157   return createProtoLayers(gctx, -1);
0158 }
0159
0160 std::vector<ProtoLayerSurfaces> ProtoLayerCreator::positiveProtoLayers(
0161     const Acts::GeometryContext& gctx) const {
0162   return createProtoLayers(gctx, 1);
0163 }
0164
0165 ProtoLayerCreator::ProtoLayerCreator(const ProtoLayerCreator::Config& cfg,
0166                                      std::unique_ptr<const Acts::Logger> log)
0167     : m_cfg(cfg), m_logger(std::move(log)) {
0168   if (!m_cfg.detectorElementFactory) {
0169     throw std::invalid_argument("Detector element factory is not set");
0170   }
0171 }
0172
0173 std::vector<ProtoLayerSurfaces> ProtoLayerCreator::createProtoLayers(
0174     const Acts::GeometryContext& gctx, int side) const {
0175   // Count the current detector modules identifiers
0176   // the return layers
0177   std::vector<ProtoLayerSurfaces> epLayers;
0178   // -------------------------------- endcap type layers
0179   // pos/neg layers
0180   std::size_t numpnLayers = m_cfg.posnegLayerPositionsZ.size();
0181   if (numpnLayers != 0u) {
0182     ACTS_DEBUG("Configured to build 2 * "
0183                << numpnLayers << " passive positive/negative side layers.");
0184     epLayers.reserve(numpnLayers);
0185
0186     /// this is the loop over the layer positions
0187     for (std::size_t ipnl = 0; ipnl < numpnLayers; ++ipnl) {
0188       // some screen output
0189       ACTS_VERBOSE("- building layer "
0190                    << ipnl << " and " << numpnLayers + ipnl << " at z = "
0191                    << side * m_cfg.posnegLayerPositionsZ.at(ipnl));
0192       /// some preparation work
0193       // define the layer envelope
0194       double layerEnvelopeR = m_cfg.posnegLayerEnvelopeR.at(ipnl);
0195       // prepare for the r binning
0196       std::vector<std::shared_ptr<Acts::Surface>> esVector;
0197       // now fill the vectors
0198       std::size_t ipnR = 0;
0199       for (auto& discModulePositions : m_cfg.posnegModulePositions.at(ipnl)) {
0200         ACTS_VERBOSE("- building ring " << ipnR << " for this layer.");
0201         // now prepare all the shared stuff
0202         // (0) module specifications
0203         double moduleThickness = m_cfg.posnegModuleThickness.at(ipnl).at(ipnR);
0204         double moduleMinHalfX = m_cfg.posnegModuleMinHalfX.at(ipnl).at(ipnR);
0205         double moduleMaxHalfX = 0.;
0206         if (m_cfg.posnegModuleMaxHalfX.size() > ipnl &&
0207             m_cfg.posnegModuleMaxHalfX.at(ipnl).size() > ipnR) {
0208           moduleMaxHalfX = m_cfg.posnegModuleMaxHalfX.at(ipnl).at(ipnR);
0209         }
0210         double moduleHalfY = m_cfg.posnegModuleHalfY.at(ipnl).at(ipnR);
0211         // (1) module bounds
0212         // create the bounds
0213         std::shared_ptr<const Acts::PlanarBounds> moduleBounds;
0214         if (moduleMaxHalfX != 0. && moduleMinHalfX != moduleMaxHalfX) {
0215           moduleBounds = std::make_shared<Acts::TrapezoidBounds>(
0216               moduleMinHalfX, moduleMaxHalfX, moduleHalfY);
0217         } else {
0218           moduleBounds = std::make_shared<Acts::RectangleBounds>(moduleMinHalfX,
0219                                                                  moduleHalfY);
0220         }
0221         // (2)) module material
0222         // create the Module material from input
0223         std::shared_ptr<const Acts::ISurfaceMaterial> moduleMaterialPtr =
0224             nullptr;
0225         if (!m_cfg.posnegModuleMaterial.empty()) {
0226           // and create the shared pointer
0227           moduleMaterialPtr = m_cfg.posnegModuleMaterial.at(ipnl).at(ipnR);
0228         }
0229
0230         // low loop over the phi positions and build the stuff
0231         for (auto& ringModulePosition : discModulePositions) {
0232           // the module transform from the position
0233           double modulePhi = phi(ringModulePosition);
0234           // the center position of the modules
0235           Acts::Vector3 moduleCenter(ringModulePosition);
0236           moduleCenter.z() *= side;
0237           // the rotation matrix of the module
0238           Acts::Vector3 moduleLocalY(std::cos(modulePhi), std::sin(modulePhi),
0239                                      0.);
0240           // take different axis to have the same readout direction
0241           Acts::Vector3 moduleLocalZ(0., 0., side * 1.);
0242           Acts::Vector3 moduleLocalX = moduleLocalY.cross(moduleLocalZ);
0243           // local rotation matrices
0244           // create the RotationMatrix - negative side
0245           Acts::RotationMatrix3 moduleRotation;
0246           moduleRotation.col(0) = moduleLocalX;
0247           moduleRotation.col(1) = moduleLocalY;
0248           moduleRotation.col(2) = moduleLocalZ;
0249           // the transforms for the front module
0250           const Acts::Transform3 moduleTransform =
0251               Acts::makeTransform3(moduleRotation, moduleCenter);
0252
0253           // create the module
0254           auto moduleElement =
0255               m_cfg.detectorElementFactory(moduleTransform, moduleBounds,
0256                                            moduleThickness, moduleMaterialPtr);
0257           // register the surface
0258           esVector.push_back(moduleElement->surface().getSharedPtr());
0259           // now deal with the potential backside
0260           if (!m_cfg.posnegModuleBacksideGap.empty()) {
0261             // the new centers
0262             moduleCenter =
0263                 moduleCenter +
0264                 m_cfg.posnegModuleBacksideGap.at(ipnl).at(ipnR) * moduleLocalZ;
0265             // the backside transforms
0266             Acts::Transform3 bsModuleTransform =
0267                 Acts::makeTransform3(moduleRotation, moduleCenter);
0268             // apply the stereo
0269             if (!m_cfg.posnegModuleBacksideStereo.empty()) {
0270               // twist by the stereo angle
0271               double stereoBackSide =
0272                   m_cfg.posnegModuleBacksideStereo.at(ipnl).at(ipnR);
0273               bsModuleTransform *=
0274                   Acts::AngleAxis3(-stereoBackSide, Acts::Vector3::UnitZ());
0275             }
0276             // everything is set for the next module
0277             auto bsModuleElement = m_cfg.detectorElementFactory(
0278                 bsModuleTransform, moduleBounds, moduleThickness,
0279                 moduleMaterialPtr);
0280             // register the backside surface
0281             esVector.push_back(bsModuleElement->surface().getSharedPtr());
0282           }
0283         }
0284         // counter of rings
0285         ++ipnR;
0286       }
0287       // the binning
0288       std::size_t layerBinsR = m_cfg.posnegModulePhiBins.at(ipnl).size();
0289       // never multiply 1 single r-bin, does not make sense
0290       if (layerBinsR > 1) {
0291         // multiply with the given bin multiplier
0292         layerBinsR *= m_cfg.posnegLayerBinMultipliers.first;
0293       }
0294       std::size_t layerBinsPhi = 0;
0295       // take the minimum phi bins in that layer
0296       for (unsigned int phiBins : m_cfg.posnegModulePhiBins.at(ipnl)) {
0297         layerBinsPhi = phiBins < layerBinsPhi ? phiBins : layerBinsPhi;
0298         layerBinsPhi *= m_cfg.posnegLayerBinMultipliers.second;
0299       }
0300       // create the layers with the surface arrays
0301       Acts::MutableProtoLayer ple(gctx, esVector);
0302       ple.envelope[Acts::AxisDirection::AxisR] = {layerEnvelopeR,
0303                                                   layerEnvelopeR};
0304       ple.envelope[Acts::AxisDirection::AxisZ] = {
0305           m_cfg.approachSurfaceEnvelope, m_cfg.approachSurfaceEnvelope};
0306
0307       // push it into the layer vector
0308       ProtoLayerSurfaces ples{std::move(ple), esVector, layerBinsR,
0309                               layerBinsPhi};
0310       epLayers.push_back(std::move(ples));
0311     }
0312   }
0313   return epLayers;
0314 }
0315
0316 }  // namespace ActsExamples::Generic